参考清风数学建模资料粒子群章节
改进粒子群算法
改进一 ：自适应惯性权重
（1）假设现在求最小值问题，那么：
注意！！！
代码段：
f_i = fit(i); % 取出第i个粒子的适应度
f_avg = sum(fit)/n; % 计算此时适应度的平均值
f_min = min(fit); % 计算此时适应度的最小值
if f_i <= f_avg
if f_avg ~= f_min % 如果分母为0，我们就干脆令w=w_min
w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
else
w = w_min;
end
else
w = w_max;
end
假设现在求最大值问题，那么：
代码段：
f_i = fit(i); % 取出第i个粒子的适应度
f_avg = sum(fit)/n; % 计算此时适应度的平均值
f_max = max(fit); % 计算此时适应度的最大值
if f_i >= f_avg
if f_avg ~= f_max % 如果分母为0，我们就干脆令w=w_min
w = w_min + (w_max - w_min)*(f_max - f_i)/(f_max - f_avg);
else
w = w_min;
end
else
w = w_max;
end
说明：
一个较大的惯性权值有利于全局搜索
而一个较小的权值则更利于局部搜索
假设现在一共五个粒子ABCDE，此时它们的适应度分别为1, 2, 3, 4, 5
取最大惯性权重为0.9，最小惯性权重为0.4
那么，这五个粒子的惯性权重应该为：0.4, 0.65, 0.9, 0.9, 0.9
适应度越小，说明距离最优解越近，此时更需要局部搜索
适应度越大，说明距离最优解越远，此时更需要全局搜索
改进二：压缩(收缩)因子法
Clerc构造了引入收缩因子的 PSO模型，采用了压缩因子，这种调整方法通过合适选取参 数，可确保PSO算法的收敛性，并可取消对速度的边界限制。
代码：
w = w_min + (w_max - w_min)*rand(1) + sigma*randn(1);
例子：
求下面这个函数的最大值。
运行结果：
最佳的位置是：
11.6255 5.7250
此时最优值是：
38.8503
演示图：
迭代图：
例题代码：
main：
%% 综合应用：粒子群算法PSO:
clear; clc
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 1000; % 粒子数量
narvs = 2; % 变量个数
c1 = 2.05; % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2.05; % 每个粒子的社会学习因子，也称为社会加速常数
C = c1+c2;
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w_max = 0.9; % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
K = 100; % 迭代的次数
vmax = [1.51 0.17]; % 粒子的最大速度
x_lb = [-3 4.1]; % x的下界
x_ub = [12.1 5.8]; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1); % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度（这里我们设置为[-vmax,vmax]）
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun_2(x(i,:)); % 调用Obj_fun_2函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置
f_i = fit(i); % 取出第i个粒子的适应度
f_avg = sum(fit)/n; % 计算此时适应度的平均值
f_max = max(fit); % 计算此时适应度的最大值
if f_i >= f_avg
if f_avg ~= f_max % 如果分母为0，我们就干脆令w=w_min
w = w_min + (w_max - w_min)*(f_max - f_i)/(f_max - f_avg);
else
w = w_min;
end
else
w = w_max;
end
v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:))); % 更新第i个粒子的速度
x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
% 如果粒子的位置超出了定义域，就对其进行调整
for j = 1: narvs
if x(i,j) < x_lb(j)
x(i,j) = x_lb(j);
elseif x(i,j) > x_ub(j)
x(i,j) = x_ub(j);
end
end
fit(i) = Obj_fun_2(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun_2(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 那就更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun_2(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 那就更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun_2(gbest); % 更新第d次迭代得到的最佳的适应度
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun_2(gbest))子程序：
function y = Obj_fun_2(x)
y = 21.5 + x(1)*sin(4*pi*x(1)) + x(2)*sin(20*pi*x(2));
end
