赞
踩
先简单介绍一下粒子群优化算法(Particle Swarm Optimization),后边会介绍一些改进的粒子群算法。
受到鸟群觅食行为的启发(鸟群觅食,通过信息共享使种群找到最优的觅食点),由社会心理学家JamesKennedy和电气工程师RussellEberhart于1995年提出,用于解决科学工程领域的非线性,非凸性,组合优化问题;在函数优化,图像处理也有广泛的应用。
粒子群优化算法是一种基于数值的优化算法,粒子群优化算法的基础是“信息共享”。具有收敛速度快,参数少,算法简单易于实现的优点(对于高维度比遗传算法更快)。存在早熟收敛,维数灾难,陷入局部最优的问题。
场景,一群鸟随机搜寻食物,一片森林只有一块地方有食物,所有鸟都不知道食物在哪里,只知道自己和其他鸟距离食物的位置,寻找食物的策略就是所有鸟寻找当前距离食物最近的鸟。
粒子群优化算法的相关发展方向:
用粒子来描述鸟类个体,每个粒子在搜索空间中搜索,飞行过程为搜索过程,飞行速度可以根据粒子历史最优位置和种群历史最优位置动态调整,粒子当前位置就是对应优化问题的一个候选解 ,称为个体极值,粒子群中的最优的个体极值称为全局最优解。不断地迭代,粒子更新速度和我位置,最终满足条件的全局最优解即为所求目标。
鸟群觅食 | 粒子群算法 |
鸟 | 粒子 |
森林 | 求解空间 |
食物的量 | 目标函数值(适应值) |
每只鸟所在位置 | 空间中的一个解(粒子位置) |
食物量最多的位置 | 全局最优解 |
MATLAB伪代码
- FOR each particle(粒子) i
- FOR each dimension(维度) d
- 初始化粒子位置 Xid 为允许的随机数
- 初始化粒子速度 Vid 为允许的随机数
- END FOR
- END FOR
-
- Iteration k=1 %迭代次数
- DO
- % 找到每个粒子的最优适应值
- FOR each particle i
- 计算适应值
- IF 这个适应值比历史的最优适应值好
- SET 当前的适应值为最优的适应值
- END IF
- END FOR
- 挑选出最优适应值的粒子就是当前迭代次数的群体最优值
- FOR each particle i
- FOR each dimension d
- 计算更新粒子速度 % 速度更新公式
- 计算更新位置 % 位置更新公式
- END FOR
- END FOR
- k=k+1
- While 当满足最大迭代次数或者两次迭代之间的适应值误差满足目标时停止迭代
(1)初始化:设置最大迭代次数,目标函数的自变量个数,粒子的最大速度,位置信息为整个搜索空间,我们在速度区间和搜索空间上随机初始化速度和位置,设置粒子群规模为N,每个粒子随机初始化一个飞翔速度。
假设在D维空间中:
其他参数:
(2)定义适应度函数
即需要优化的目标函数
(3)更新速度公式
说明:
速度方向更新示意图
公式定义参考:
(4)更新位置公式
(5)终止条件
(6)算法参数解释
粒子群规模N:推荐取值范围 [20-1000] ,简单问题取20-40;较难问题取100-200。较小的种群规模容易陷入局部最优,较大的会提高收敛性(迭代计算量大),种群规模达到一定水平,再增大不会再有显著作用。
粒子维度D:粒子搜索的空间维度数(自变量的个数)
迭代次数K:推荐取值 [50-100] ,典型取值60,70,100,迭代次数太小解不稳定,太大浪费资源,没必要。
惯性因子(惯性权重):w 值较大,全局寻优能力强。参数意义:使粒子保持运动惯性,使其有搜索扩展空间的趋势, w 越大,探索新的区域的能力越强。也表示粒子对当前自身运动状态的信任,依据自身的速度进行惯性运动。 较大的 w 有利于跳出局部极值,而较小的 w 有利于算法收敛。大惯性权重有利于全局搜索,不至于陷入局部最优,小惯性权重有利于在局部搜索,可以快速收敛得到最优解。推荐取值 [0.4-2],典型取值0.9, 1.2, 1.5,1.8
改善w:解决实际优化问题时,往往希望先采用全局搜索,使搜索空间快速收敛于某一区域,然后采用局部精细化搜索获得高精度的解。提出了自适应调整的策略,即随着迭代次数,线性的减小 w 的值。常用方法:线性变化策略:随着迭代次数的增加,惯性权重 w 不断减,从而使得粒子群算法在初期具有较强的全局收敛能力,在后期具有较强的局部收敛能力。
学习因子C1,C2:也称为加速系数或加速因子
(7)粒子群算法的一些技巧
适应值:即优化目标函数的值,用来评价粒子位置的好坏程度,决定是否更新粒子个体的历史最优位置和群体的历史最优位置,保证粒子朝着最优解的方向搜索。
位置限制:限制粒子搜索的空间,即自变量的取值范围,对于无约束问题此处可以省略。
速度限制:为了平衡算法探索能力与开发能力,需要设定一个合理的速度范围,限制粒子的最大速度。粒子飞行速度快,探索能力强,但粒子容易飞过最优解。粒子飞行速度快。飞行速度慢,开发能力强,但是收敛速度慢,且容易陷入局部最优解。一般设为每维变量变化范围的10%~20%。
初始化:如果粒子的初始化范围选择得好的话可以缩短优化的收敛时间,需要根据具体的问题进行分析。如果根据经验判断出最优解一定在某个范围内,则就在这个范围内初始化粒子;如果无法确定,则以粒子的取值边界作为初始化范围。
1.一维
- clc;
- clear;
- close all;
-
- f= @(x) - (x - 10) .^ 2 + x .* sin(x) .* cos(2 * x) - 5 * x .* sin(3 * x) ; % 适应度函数表达式(求这个函数的最大值)
- figure(1);
- fplot(f, [0 20], 'b-'); % 画出初始图像
-
- d = 1; % 空间维数(该例子是1维)
- N = 15; % 初始种群个数
- x_limit = [0, 20]; % 设置位置限制
- v_limit = [-1, 1]; % 设置速度限制
-
- x = x_limit(1) + (x_limit(2) - x_limit(1)) * rand(N, d); %初始每个粒子的位置
- v = rand(N, d); % 初始每个粒子的速度
-
- pbest = x; % 初始化每个个体的历史最佳位置
- gbest = zeros(1, d); % 初始化种群的历史最佳位置 ,创建全0数组
-
- fp_best = zeros(N, 1); % 初始化每个个体的历史最佳适应度 为 0
- fg_best = -inf; % 初始化种群历史最佳适应度 为 负无穷
-
- iter = 50; % 最大迭代次数
-
- c1 = 1.6; % 自我学习因子
- c2 = 1.8; % 群体学习因子
- w=0.4; % 惯性学习因子
- hold on;
- plot(x, f(x), 'ro');
- title('初始状态图');
-
- record = zeros(iter, 1); % 记录器(用于记录 fg_best 的变化过程)
- figure(2);
- i = 1;
- while i <= iter
- fx = f(x) ; % 计算个体当前适应度
- % 个体寻优
- for j = 1:N
- if fp_best(j) < fx(j) % 第j个个体的历史(最优)适应度小于当前的适应度则更新
- fp_best(j) = fx(j); % 更新个体历史最佳适应度
- pbest(j) = x(j); % 更新个体历史最佳位置
- end
- end
- % 群体寻优
- if fg_best < max(fp_best) % 第i次迭代,种群历史最优适应值,小于当前迭代次数的群体最优适应值则更新
- [fg_best,ind_max] = max(fp_best); % 更新群体历史最佳适应度
- gbest = pbest(ind_max); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
- % 注:将上面的两个式子换成下面两个是不可以的
- % [gbest, ind_max] = max(pbest); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
- % fg_best = fp_best(ind_max); % 更新群体历史最佳适应度
- end
-
- v = v * w + c1 * rand() * (pbest - x) + c2 * rand() * (repmat(gbest, N, 1) - x); % 速度更新
- % 注: repmat(A,r1,r2):可将矩阵 扩充 为每个单位为A的r1*r2的矩阵
-
- % 边界速度处理
- v(v > v_limit(2)) = v_limit(2);
- v(v < v_limit(1)) = v_limit(1);
-
- x = x + v;% 位置更新
-
- % 边界位置处理
- x(x > x_limit(2)) = x_limit(2);
- x(x < x_limit(1)) = x_limit(1);
-
- record(i) = fg_best;%最大值记录
-
- % 画动态展示图
- zuo_x = 0 : 0.01 : 20;
- plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
- title('状态位置变化')
- pause(0.1)
- i = i + 1;
- % if mod(i,10) == 0 % 显示进度
- % i
- % end
- end
- figure(3);
- plot(record);
- title('收敛过程');
- zuo_x = 0 : 0.01 : 20;
-
- figure(4);
- plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
- title('最终状态图')
-
- disp(['最佳适应度:',num2str(fg_best)]);
- disp(['最佳粒子的位置x:',num2str(gbest)]);
2.二维
- % 主调用函数,以下新建一个文件,命名为f_xy
- % function z = f_xy(x, y)
- % % 适应度函数
- % z = - x.^2 - y.^2 + 10*cos(2*pi*x) + 10*cos(2*pi*y) + 100;
- % end
-
- clc;
- clear;
- close all;
-
- [ZuoBiao_x, ZuoBiao_y] = meshgrid(-10:0.1:10,-5:0.1:5); % 产生二维坐标
- ZuoBiao_z = f_xy(ZuoBiao_x, ZuoBiao_y);
- figure(1);
- s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z); % 画网格曲面图
- s.FaceColor = 'flat'; % 修改网格图的属性
-
- % 初始化种群
- N = 100; % 初始种群个数
- D = 2; % 空间维数
- iter = 50; % 迭代次数
- x_limit = [-10, 10; -5, 5]; % 位置限制
- v_limit = [ -2, 2; -1, 1]; % 速度限制
-
- x = zeros(N, D);
- for i = 1:D
- x(:,i) = x_limit(i, 1) + (x_limit(i, 2) - x_limit(i, 1)) * rand(N, 1);%初始种群的位置
- end
- v(:,1) = rands(N, 1) * 1; % 初始种群x方向的速度
- v(:,2) = rands(N, 1) * 2; % 初始种群y方向的速度
-
- p_best = x; % (初始化)每个个体的历史最佳位置
- f_best = zeros(1, D); % (初始化)种群的历史最佳位置
-
- fp_best = zeros(N, 1) - inf; % (初始化)每个个体的历史最佳适应度为负无穷
- fg_best = -inf; % (初始化)种群历史最佳适应度为负无穷
-
- w = 0.8; % 惯性权重
- c1 = 0.5; % 自我学习因子
- c2 = 0.5; % 群体学习因子
-
- figure(2);
- s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z); % 画网格曲面图
- s.FaceColor = 'flat'; % 修改网格图的属性
- hold on;
- plot3(x(:,1), x(:,2),f_xy(x(:,1), x(:,2)), 'ro');
- title('初始状态图');
-
- figure(3);
- i = 1;
- record = zeros(iter, 1); % 记录器
- while i <= iter
- fx = f_xy(x(:,1), x(:,2)); % 个体当前适应度
- for j = 1:N
- if fp_best(j) < fx(j) % 记录最大值
- fp_best(j) = fx(j); % 更新个体历史最佳适应度
- p_best(j,:) = x(j,:); % 更新个体历史最佳位置
- end
- end
- if fg_best < max(fp_best)
- [fg_best, ind_max] = max(fp_best); % 更新群体历史最佳适应度
- f_best = p_best(ind_max, :); % 更新群体历史最佳位置
- end
-
- v = v * w + c1 * rand() * (p_best - x) + c2 * rand() * (repmat(f_best, N, 1) - x); % 速度更新
-
- % 速度处理
- for t=1:N
- for k=1:D
- if v(t,k) > v_limit(k,2) % 超速处理
- v(t,k) = v_limit(k,2);
- elseif v(t,k) < v_limit(k,1) % 慢速处理
- v(t,k) = v_limit(k,1);
- end
- end
- end
-
- x = x + v; % 位置更新
-
- % 边界处理
- for t=1:N
- for k=1:D
- if x(t,k) > x_limit(k,2) % 超过边界上限
- x(t,k) = x_limit(k,2);
- elseif x(t,k) < x_limit(k,1) % 超过边界下限
- x(t,k) = x_limit(k,1);
- end
- end
- end
-
- record(i) = fg_best; % 最大值记录
- i = i + 1;
- if mod(i,10) == 0
- i % 收敛进度输出
- end
- pause(0.1)
-
- end
-
- figure(4)
- plot(record);
- title('收敛过程');
-
- % 画最终状态图
- figure(5);
- [ZuoBiao_x, ZuoBiao_y] = meshgrid(-10:0.1:10,-5:0.1:5); % 产生二维坐标
- s = mesh(ZuoBiao_x, ZuoBiao_y, ZuoBiao_z); % 画网格曲面图
- s.FaceColor = 'flat'; % 修改网格图的属性
- hold on;
- scatter3(x(:,1), x(:,2), f_xy(x(:,1), x(:,2)), 'ro');
- title('最终状态图')
-
- disp(['最大值:',num2str(fg_best)]);
- disp(['变量取值:',num2str(f_best)]);
-
- % 下面一段用来输出 粒子群最佳适应度的排行榜
- [socre,ind] = sort(fp_best, 'descend') % 降序排序
- DaAn = zeros(3,2);
- for i=1:N
- row = ind(i);
- DaAn(i,:) = p_best(row,:);
- end
改进思路:
(1)修改惯性权重
A:
B:
后续看到文献会补充!!!
小白一个,参照多篇,自我总结,侵删(私聊我),勿喷!!!
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。