粒子群算法 笔记 数学建模

引入:

如何找到全局最大值：如果只是贪心的话，容易被局部最大解锁定

方法有：盲目搜索，启发式搜索

盲目搜索：枚举法和蒙特卡洛模拟，但是样例太多花费巨量时间

所以启发式算法就来了，通过经验和规律减少尝试的次数去解决这个问题

启发式算法的总结：

1）有助于找到最优解，但是不一定能找到最优解

2）有助于加速找到最优解

粒子群算法：

由肯尼迪和埃伯哈特于 1995 年提出， 是一种受鸟群、鱼群等群体生物觅食行为启发的智能优化算法，在多领域广泛应用。以下从原理、算法流程、应用场景、优缺点展开介绍：

基本原理：

将优化问题的潜在解视为搜索空间中的粒子，每个粒子都有位置和速度两个属性。粒子在解空间中以一定速度飞行，通过不断调整自身位置来搜索最优解。在飞行过程中，粒子会参考自身历史最优位置（pbest）和群体历史最优位置（gbest）来更新自己的速度和位置。例如，在一个二维平面上寻找函数最小值问题中，每个粒子代表一个可能的解点，粒子根据自身曾经到达过的最优位置以及整个群体目前找到的最优位置来决定下一步往哪个方向移动以及移动多远。

这只鸟第d步所在的位置 = 第d‐1步所在的位置 + 第d‐1步的速度 * 运动的时间

x(d) = x(d‐1) + v(d‐1)  * t    (每一步运动的时间t一般取1)

这只鸟第d步的速度 =上一步自身的速度惯性 + 自我认知部分 +社会认知部分

v(d) = w*v(d‐1)  + c1*r1*(pbest(d)‐x(d))  + c2*r2*(gbest(d)‐x(d))   （三个部分的和）

粒子群算法中的基本概念

Ø粒子：优化问题的候选解

Ø位置：候选解所在的位置

Ø速度：候选解移动的速度

Ø适应度：评价粒子优劣的值，一般设置为目标函数值

Ø个体最佳位置：单个粒子迄今为止找到的最佳位置

Ø群体最佳位置：所有粒子迄今为止找到的最佳位置

算法流程

初始化：随机生成一群粒子，确定每个粒子的初始位置和速度。这些初始值在问题的解空间内随机分布，为后续搜索奠定基础。

适应度计算：根据优化问题的目标函数，计算每个粒子的适应度值，以此评估粒子当前位置的优劣。例如在求解函数最大值问题中，将粒子位置代入函数计算得到的函数值就是该粒子的适应度。

更新个体最优和全局最优：每个粒子将自身当前适应度与其历史最优位置的适应度进行比较，若当前更好，则更新自身历史最优位置（pbest）。同时，所有粒子的适应度相互比较，找出整个群体当前的最优位置（gbest）。

速度和位置更新：依据特定公式，结合粒子自身的速度、当前位置与个体最优位置、全局最优位置的距离，更新粒子的速度和位置。公式中通常包含惯性权重、学习因子等参数，用于平衡粒子的全局搜索和局部搜索能力。例如，惯性权重较大时，粒子倾向于保持原有飞行方向，进行全局搜索；学习因子较大时，粒子更关注自身经验和群体经验，加强局部搜索。

终止判断：检查是否满足终止条件，如达到最大迭代次数、适应度值收敛到一定精度等。若满足，则输出全局最优解；否则，返回适应度计算步骤继续迭代。

例题：求一元函数的最大值

first先画出一开始1函数图

x = -3:0.01:3;

y = 11*sin(x) + 7*cos(5*x);

figure(1)

plot(x,y,'b-')

title('y = 11*sin(x) + 7*cos(5*x)')

hold on % 不关闭图形，继续在上面画图

然后设置粒子的参数

n = 10; % 粒子数量

narvs = 1; % 变量个数

c1 = 2; % 每个粒子的个体学习因子，也称为个体加速常数

c2 = 2; % 每个粒子的社会学习因子，也称为社会加速常数

w = 0.9; % 惯性权重

K = 50; % 迭代的次数

vmax = 1.2; % 粒子的最大速度

x_lb = -3; % x的下界

x_ub = 3; % 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]）

计算粒子的适应度(也就是对应的y值)

fit = zeros(n,1); % 初始化这n个粒子的适应度全为0

for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度

fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度（这里写成x(i,:)主要是为了和以后遇到的多元函数互通）

end

找到适应度最大的下标

pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）

ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标

并记录最大值

gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

模拟绘制出过程

h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数,80是我设置的散点显示的大小（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）

接下来迭代k次来更新位置和速度

伪代码:

迭代k次

更新每个例子的信息

调整速度

调整位置

检测速度是否超出范围

检测位置是否超出定义域

更新当前的适应度

检测自学习是否需要更新

检测社会学习是否要更新

fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度

for d = 1:K % 开始迭代，一共迭代K次

for i = 1:n % 依次更新第i个粒子的速度与位置

v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)); % 更新第i个粒子的速度

% 如果粒子的速度超过了最大速度限制，就对其进行调整

for j = 1: narvs

if v(i,j) < -vmax(j)

v(i,j) = -vmax(j);

elseif v(i,j) > vmax(j)

v(i,j) = vmax(j);

end

end

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_fun1(x(i,:)); % 重新计算第i个粒子的适应度

if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度

pbest(i,:) = x(i,:); % 那就更新第i个粒子迄今为止找到的最佳位置

end

if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度

gbest = pbest(i,:); % 那就更新所有粒子迄今为止找到的最佳位置

end

end

fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度

pause(0.1) % 暂停0.1s

h.XData = x; % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）

h.YData = fit; % 更新散点图句柄的y轴的数据（此时粒子的位置在图上发生了变化）

end

figure(2)

plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图

xlabel('迭代次数');

disp('最佳的位置是：'); disp(gbest)

disp('此时最优值是：'); disp(Obj_fun1(gbest))

例题2:求二元函数的最小值

x1 = -15:1:15;

x2 = -15:1:15;

[x1,x2]=meshgrid(x1,x2);

y = x1.^2+x2.^2-x1.*x2-10*x1-4*x2+60;

figure(1);

mesh(x1,x2,y);

title('y = 11*sin(x) + 7*cos(5*x)')

hold on % 不关闭图形，继续在上面画图

n = 50; % 粒子数量

narvs = 2; % 变量个数

c1 = 2; % 每个粒子的个体学习因子，也称为个体加速常数

c2 = 2; % 每个粒子的社会学习因子，也称为社会加速常数

w = 0.9; % 惯性权重

K = 50; % 迭代的次数

vmax =[4,4]; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x_lb = [-15,-15]; % x%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x_ub = [15,15]; % 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);

fit = zeros(n,1); % 初始化这n个粒子的适应度全为0

for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度

fit(i) = only(x(i,:)); %&&&&%%%%%%%%%

end

pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）

ind = find(fit == min(fit), 1); % 找到适应度最大的那个粒子的下标

gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）

h = scatter3(x(:,1),x(:,2),fit,100,'*r'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度

for d = 1:K % 开始迭代，一共迭代K次

for i = 1:n % 依次更新第i个粒子的速度与位置

v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)); % 更新第i个粒子的速度

% 如果粒子的速度超过了最大速度限制，就对其进行调整

for j = 1: narvs

if v(i,j) < -vmax(j)

v(i,j) = -vmax(j);

elseif v(i,j) > vmax(j)

v(i,j) = vmax(j);

end

end

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) = only(x(i,:)); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if fit(i) < only(pbest(i,:)) % 555555555555555555455555555555555555%%%%%%

pbest(i,:) = x(i,:); % 那就更新第i个粒子迄今为止找到的最佳位置

end

if fit(i) < only(gbest) % %%%%%%%%%%%%%%%%%%%%%%%%%%%

gbest = pbest(i,:); % 那就更新所有粒子迄今为止找到的最佳位置

end

end

fitnessbest(d) = only(gbest); % %%%%%%%%%%%%%%%%%%%

pause(0.1) % 暂停0.1s

h.XData = x(:,1); % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）

h.YData = x(:,2); % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）

h.ZData = fit; % 更新散点图句柄的y轴的数据（此时粒子的位置在图上发生了变化）

end

figure(2)

plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图

xlabel('迭代次数');

disp('最佳的位置是：'); disp(gbest)

disp('此时最优值是：'); disp(only(gbest))

改进::

-线性递减惯性权重

惯性权重w体现的是粒子继承先前的速度的能力，Shi,Y最先将惯性权重w引入到粒子群算法中，并分析指出一个较大的惯性 权值有利于全局搜索，而一个较小的权值则更利于局部搜索。为了更好地平衡算法的全局搜索以及局部搜索能力，,公式如下

参数的设定

w_max = 0.9; % 最大惯性权重，通常取0.9   w_min = 0.4; % 最小惯性权重，通常取0.4

枚举粒子的时候加入

w = w_start - (w_start - w_end) * d / K;

-非线性递减惯性权重

三种递减惯性权重

-自适应惯性权重

(假如我们现在在求最小值问题)

适应度越小,距离我们要求的最小值越近,就要局部搜,不然步幅太大,容易错过,当我们的适应度越大,代表我们距离最小值越远,越要全局搜索

(假如我们现在在求最大值问题)

参数的设定

w_max = 0.9; % 最大惯性权重，通常取0.9   w_min = 0.4; % 最小惯性权重，通常取0.4

枚举粒子都时候加入

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

-随机惯性权重:

参数的设定:

w_max = 0.9; % 最大惯性权重，通常取0.9

w_min = 0.4; % 最小惯性权重，通常取0.4

sigma = 0.3; % 正态分布的随机扰动项的标准差（一般取0.2-0.5之间）

枚举粒子都时候加入:

w = w_min + (w_max - w_min)*rand(1) + sigma*randn(1);

-收缩因子法:

个体学习因子c1和社会(群体)学习因子c2决定了粒子本身经验信息和其他粒子的经验信息对粒子运行轨迹的影响，其反映了粒子群之间的信息交流。设置c1较大的值，会使粒子

过多地在自身的局部范围内搜索，而较大的c2的值，则又会促使粒子过早收敛到局部最优值。为了有效地控制粒子的飞行速度，使算法达到全局搜索与局部搜索两者间的有效平衡，Clerc构造了引入收缩因子的PSO模型，采用了压缩因子，这种调整方法通过合适选取参

数，可确保PSO算法的收敛性，并可取消对速度的边界限制。

参数设定

c1 = 2.05; % 每个粒子的个体学习因子，也称为个体加速常数

c2 = 2.05; % 每个粒子的社会学习因子，也称为社会加速常数

C = c1+c2;

fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子

w = 0.9; % 惯性权重

更新粒子速度

v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:))); % 更新第i个粒子的速度

-非对称学习因子

设计的思路:一开始强调各干各的,最后朝着全局最优解靠拢

参数设定:

c1_ini = 2.5; % 个体学习因子的初始值

c1_fin = 0.5; % 个体学习因子的最终值

c2_ini = 1; % 社会学习因子的初始值

c2_fin = 2.25; % 社会学习因子的最终值

每次循环更新因子

c1 = c1_ini + (c1_fin - c1_ini)*d/K;

c2 = c2_ini + (c2_fin - c2_ini)*d/K;

每次更新粒子的速度

v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)); % 更新第i个粒子的速度

-自动退出迭代循环

一个记录最大值,一个记录最大值计数器

如果多次在最大值极小附近变化,那么计数器+1

反之如果大于最大值一定量,就计数器清空

%% 改进：自动退出迭代循环

Count = 0; % 计数器初始化为0

max_Count = 30; % 最大计数值初始化为30

tolerance = 1e-6; % 函数变化量容忍度，取10^(-6)

每次循环完一系列的粒子,检测最大值的delta

delta_fit = abs(Obj_fun3(gbest) - Obj_fun3(tem)); % 计算相邻两次迭代适应度的变化量

if delta_fit < tolerance % 判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1

Count = Count + 1;

else

Count = 0; % 否则计数器清0

end 

if Count > max_Count % 如果计数器的值达到了最大计数值

break; % 跳出循环

end

matlab自带的粒子群函数(只能求最小值)

narvs = 2; % 变量个数

x_lb = [-15 -15]; % x的下界(长度等于变量的个数，每个变量对应一个下界约束)

x_ub = [15 15]; % x的上界

[x,fval,exitflag,output] = particleswarm(@Obj_fun2, narvs, x_lb, x_ub)

粒子群算法的应用:

1求解方程组

将方程组移动到右边,变成一个==0的表达式f,我们尝试去求得f^2(x)或者|f(x)|的最小值,这个时候的x的取值就是解

function F = my_fun(x)

F(1) = abs(x(1)+x(2))-abs(x(3));

F(2) = x(1) * x(2) * x(3) + 18;

F(3) = x(1)^2*x(2)+3*x(3);

end

clear; clc

narvs = 3;

% 使用粒子群算法，不需要指定初始值，只需要给定一个搜索的范围

lb = -10*ones(1,3); ub = 10*ones(1,3);

options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',20000,'SwarmSize',100);

[x, fval] = particleswarm(@Obj_fun,narvs,lb,ub,options)

2,拟合函数

同样的将多个点的值求平方和,通过多次的搜索确定一个最接近0的值

如果求解有约束的,那么可以是fmincon

如果是没有约束的,可以是fminsearch和fminunc

k0 = [0 0]; %初始值，注意哦，这个初始值的选取确实挺难，需要尝试。。。

lb = []; ub = [];

% lb = [-1 -1]; ub = [2 2];

[k, fval] = fmincon(@Obj_fun,k0,[],[],[],[],lb,ub)

%% 求解无约束非线性函数的最小值的两个函数

[k, fval] = fminsearch(@Obj_fun,k0) % Nelder-Mead单纯形法求解最小值，适用于解决不可导或求导复杂的函数优化问题

[k, fval] = fminunc(@Obj_fun,k0) % 拟牛顿法求解无约束最小值，适用于解决求导容易的函数优化问题

多元函数拟合

使用matlab的lsqcurvefit(@fun,k0,)

注意这里的fun函数,对于fmincon返回的是残差平方和,这个函数返回的是表达式

load data_x_y

k0 = [0 0];

lb = []; ub = [];

[k, fval] = lsqcurvefit(@fit_fun,k0,x,y,lb,ub)