粒子群优化算法的实现方式_matlab粒子群优化算法

2022-11-10 15:43:14 浏览数 (1)

文章目录

  • 1 算法基本概念
  • 2 算法的MATLAB实现
    • 2.1 算法的基本程序
    • 2.2 适应度函数
      • 示例
    • 2.3 免疫粒子群算法的MATLAB应用
  • 3 粒子群算法的权重控制
    • 3.1 线性递减法
    • 3.2 自适应法
      • 3.2.1 根据全局最优点距离进行调整
      • 3.2.2 依据早熟收敛程度和适应值进行调整权重
  • 4 混合粒子群算法
  • 参考文献

1 算法基本概念

粒子群优化算法属于进化算法的一种，通过追随当前搜索到的最优值来寻找全局最优。粒子群算法也称粒子群优化算法（Particle Swarm Optimization，PSO），PSO有几个关键概念：粒子、优化函数、适值（Fitness Value）、飞行方向、飞行距离。

粒子群优化算法实现容易、精度高、收敛快，在解决实际问题中展示了其优越性。粒子群算法通用性较好，适合处理多种类型的目标函数和约束，并且容易与传统的优化方法结合，从而改进自身的局限性，更高效地解决问题。因此，将粒子群算法应用于解决多目标优化问题上具有很大的优势。

2 算法的MATLAB实现

基本粒子群算法使用固定长度的二进制符号串来表示群体中的个体，其等位基因是由二值符号集 { 0 , 1 } {0,1} { 0,1} 所组成。初始群体中各个个体的基因值可用均匀分布的随机数来生成。

2.1 算法的基本程序

基本粒子群（PSO）算法描述如下：

代码语言：javascript复制
begin
	Initalize; %包括初始化粒子群数,粒子初始速度和位置
		[x,xd] = judge(x,pop_size); %调用judge函数,初始化第一次值

	fornum=2:最大迭代次数
		wk=wmax-num*(wmax-wmin)/max_gen; %计算惯性权重
		r1= ; r2=  %随机产生加速权重
		PSO算法
	迭代求vk,xk;
	While 判断 vk 是否满足条件
		再次重新生成加速权重系数r1;r2
		PSO算法
		再次迭代求vk,xk数值
	end
	调用[x,xd] = judge(x,pop_size); 重新计算出目标函数值
	判断并重新生成pj数值;
	判断并重新生成pjd数值
	if 迭代前数值 > 迭代后的数值
		累加迭代次数值
	end
	输出随机数种子、进度、最优迭代次数、每个函数的数值和目标函数的数值
	用ASCII保存粒子位移的数值
	用ASCII保存粒子速度的数值
end

在MATLAB中，编程实现的基本粒子群算法基本函数为PSO，其调用格式如下：

代码语言：javascript复制
[xm, dv] = PSO(fitness, N, c1, c2, w, M, D)

其中，fitness为待优化的目标函数（适应度函数），N是粒子数目，c1是学习因子 1 1 1，c2是学习因子 2 2 2，w是惯性权重，M是最大迭代次数，D是自变量个数，xm是目标函数取最小值时自变量，fv是目标函数最小值。

使用MATLAB实现基本粒子群（PSO）算法代码如下：

代码语言：javascript复制
function[xm,fv]=PSO(fitness, N, c1, c2, w, M, D)
%%%%%%%给定初始化条件%%%%%%%
% c1 学习因子1
% c2 学习因子2
% w 惯性权重
% M 最大迭代次数
% D 搜索空间维数
% N 初始化群体个体数目
%%%%%%%初始化种群个体（限定位置和速度范围）%%%%%%%
format long;
for i = 1: N
for j = 1: D
x(i,j) = randn;	% 随机初始化位置
v(i,j) = randn; % 随机初始化速度
end
end
%%%%%%%先计算各个粒子的适应度，并初始化Pi和Pg%%%%%%%
for i = 1: N
p(i) = fitness(x(i,:));
y(i,:) = x(i,:);
end
pg = x(N,:); % pg为全局最优
for i = 1:(N-1)
if fitness(x(i,:)) < fitness(pg)
pg = x(i,:);
end
end
%%%%%%%进入主循环，按照公式依次迭代，知道满足精度要求%%%%%%%
for t = 1:M
for i = 1:N % 更新速度和位移
v(i,:) = w * v(i,:)   c1 * rand * (y(i,:)-x(i,:))   c2 * rand * (pg - x(i,:));
x(i,:) = x(i,:)   v(i,:);
if fitness(x(i,:)) < p(i)
p(i) = fitness(x(i,:));
y(i,:) = x(i,:);
end
if p(i) < fitness(pg)
pg = y(i,:);
end
end
Pbest(t) = fitness(pg);
end
%%%%%%%最后给出计算结果%%%%%%%
disp('***********************')
disp('目标函数取最小值时自变量：')
xm = pg' disp('目标函数最小值为：') fv = fitness(pg) disp('***********************')

2.2 适应度函数

适应度表示个体 x x x 对环境的适应程度，分为针对被优化目标函数优化行适应度和 针对约束函数的约束型适应度。粒子群算法使用的适应度函数多样， G W GW GW 函数和 R A RA RA 函数是两类经典的适应度函数。

  • 优化型适应度 F o b j ( X ) = f ( x ) F_{obj}(X)=f(x) Fobj​(X)=f(x)
  • 约束型适应度 F i ( X ) = { 0 , g i ( X ) ≤ 0 g i ( X ) , g i ( X ) ≥ 0 F_{i}(X)=left{ begin{matrix} 0,qquad g_i(X)leq0\ \ g_i(X),quad g_i(X)geq0\ end{matrix} right. Fi​(X)=⎩⎨⎧​0,gi​(X)≤0gi​(X),gi​(X)≥0​

示例

利用基本粒子群算法求解函数 f ( x ) = ∑ i = 1 10 ( x i 2 2 x i − 3 ) f(x)=sum_{i=1}^{10}(x_i^2 2x_i-3) f(x)=∑i=110​(xi2​ 2xi​−3) 最小值。

解析

利用PSO求解最小值，需要确认迭代次数对结果的影响。设定题中函数最小点均为 0 0 0，粒子群规模为 40 40 40，惯性权重为 0.6 0.6 0.6，学习因子1为 1.2 1.2 1.2，学习因子2为 2.2 2.2 2.2，迭代次数为 100 100 100 和 300 300 300。

基本粒子群PSO算法代码见上。

目标函数代码如下：

代码语言：javascript复制
function F = fitness(x)
F = 0;
for i = 1: 10
F = F   x(i)^2   2 * x(i) - 3;
end

求解函数最小值代码如下：

  • 迭代次数为100
代码语言：javascript复制
clear all
clc
x  = zeros(1,10);
[xm,fv] = PSO(@fitness,40,1.2,2.2,0.6,100,10);	% 迭代次数为100
% 取自变量
xm;
% 取函数最小值
fv;

结果如下：

代码语言：javascript复制
***********************
目标函数取最小值时自变量：
xm =
-0.991104610834485
-0.997666388064225
-0.993499168365228
-0.995209292046358
-0.997319510025391
-0.997398011693752
-0.997812988297172
-1.003889705997851
-0.999468835940385
-1.006062399124978
目标函数最小值为：
fv =
-39.999779311591290
***********************
  • 迭代次数为300
代码语言：javascript复制
clear all
clc
x  = zeros(1,10);
[xm,fv] = PSO(@fitness,40,1.2,2.2,0.6,300,10);	% 迭代次数为300
% 取自变量
xm;
% 取函数最小值
fv;

结果如下：

代码语言：javascript复制
***********************
目标函数取最小值时自变量：
xm =
-1.005827335897396
-1.003735840310715
-1.001088656877167
-1.001724852313095
-1.003933895880758
-1.002330669862129
-1.001909424295409
-0.996959460575888
-0.993648871189404
-1.001008751197640
目标函数最小值为：
fv =
-39.999872772608128
***********************

PSO算法是一种随机算法，同样的参数也会算出不同结果，且迭代次数越大，获得解的精度不一定越高。在粒子群算法中，要想获得精度高的解，关键各个参数之间的合理搭配。

2.3 免疫粒子群算法的MATLAB应用

使用基于模拟退火的混合粒子群算法求解 f ( x ) = c o s x 1 2 − x 2 2 − 3 [ 2 ( x 1 2 x 2 2 ) ] 2 0.8 f(x)=frac{cossqrt{x_1^2-x_2^2}-3}{[2 (x_1^2 x_2^2)]^2} 0.8 f(x)=[2 (x12​ x22​)]2cosx12​−x22​ ​−3​ 0.8 最小值，其中 − 10 ≤ x i ≤ 10 -10 leq x_i leq 10 −10≤xi​≤10 ，粒子数为 60 60 60，学习因子均为 1.2 1.2 1.2，退火常数为 0.8 0.8 0.8，迭代次数为 800 800 800。

解析

免疫粒子群算法代码如下：

代码语言：javascript复制
function [x,y,Result]=PSO_immu(func,N,c1,c2,w,MaxDT,D,eps,DS,replaceP,minD,Psum)
format long;
%%%%%%给定初始化条件%%%%%%%%%%%%%%%%%%%%%%%%%%%
c1=1.2;             %学习因子1
c2=1.2;             %学习因子2
w=0.8;            %惯性权重
MaxDT=800;        %最大迭代次数
D=2;              %搜索空间维数（未知数个数）
N=60;            %初始化群体个体数目
eps=10^(-10);     %设置精度(在已知最小值时候用)
DS=8;             %每隔DS次循环就检查最优个体是否变优
replaceP=0.5;     %粒子的概率大于replaceP将被免疫替换
minD=1e-10;       %粒子间的最小距离
Psum=0;           %个体最佳的和
range=100;
count = 0;
%%%%%%初始化种群的个体%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N
for j=1:D
x(i,j)=-range 2*range*rand;  %随机初始化位置
v(i,j)=randn;  %随机初始化速度
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%先计算各个粒子的适应度，并初始化Pi和Pg%%%%%%%%%%%%%%%%%%% 
for i=1:N    
p(i)=feval(func,x(i,:));
y(i,:)=x(i,:);
end
pg=x(1,:);             %Pg为全局最优
for i=2:N
if feval(func,x(i,:))<feval(func,pg)    
pg=x(i,:);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%主循环，按照公式依次迭代，直到满足精度要求%%%%%%% 
for t=1:MaxDT
for i=1:N
v(i,:)=w*v(i,:) c1*rand*(y(i,:)-x(i,:)) c2*rand*(pg-x(i,:));
x(i,:)=x(i,:) v(i,:);
if feval(func,x(i,:))<p(i) 
p(i)=feval(func,x(i,:)); 
y(i,:)=x(i,:);
end
if p(i)<feval(func,pg)  
pg=y(i,:);
subplot(1,2,1);            
bar(pg,0.25); 
axis([0 3 -40 40 ]) ;
title (['Iteration ', num2str(t)]); pause (0.1);
subplot(1,2,2); 
plot(pg(1,1),pg(1,2),'rs','MarkerFaceColor','r', 'MarkerSize',8)
hold on;
plot(x(:,1),x(:,2),'k.');
set(gca,'Color','g')
hold off;
grid on;
axis([-100 100 -100 100 ]) ;
title(['Global Min = ',num2str(p(i))]);
xlabel(['Min_x= ',num2str(pg(1,1)),' Min_y= ',num2str(pg(1,2))]);
end
end
Pbest(t)=feval(func,pg) ;   
%     if Foxhole(pg,D)<eps                %如果结果满足精度要求则跳出循环
%         break;
%     end
%%%%%开始进行免疫%%%%%%%%%%%%%%%%%
if t>DS
if mod(t,DS)==0 && (Pbest(t-DS 1)-Pbest(t))<1e-020    %如果连续DS代数，群体中的最优没有明显变优，则进行免疫.
%在函数测试的过程中发现，经过一定代数的更新，个体最优不完全相等，但变化非常非常小，
for i=1:N                            %先计算出个体最优的和
Psum=Psum p(i);
end
for i=1:N                            %免疫程序              
for j=1:N                        %计算每个个体与个体i的距离
distance(j)=abs(p(j)-p(i));
end
num=0;   
for j=1:N                        %计算与第i个个体距离小于minD的个数
if distance(j)<minD
num=num 1;
end
end
PF(i)=p(N-i 1)/Psum;             %计算适应度概率
PD(i)=num/N;                     %计算个体浓度
a=rand;                          %随机生成计算替换概率的因子
PR(i)=a*PF(i) (1-a)*PD(i);       %计算替换概率
end
for i=1:N
if PR(i)>replaceP
x(i,:)=-range 2*range*rand(1,D);
count=count 1;
end
end
end
end    
end
%%%%%%%最后给出计算结果%%%%%%%%%%%%%%%%%%%%
x=pg(1,1);
y=pg(1,2);
Result=feval(func,pg);
%%%%%%%%%%算法结束%%%%%%%%%%%%%%%%%%
function probabolity(N,i)
PF=p(N-i)/Psum;%适应度概率
disp(PF);
for jj=1:N
distance(jj)=abs(P(jj)-P(i));
end
num=0;
for ii=1:N
if distance(ii)<minD
num=num 1;
end
end
PD=num/N;            %个体浓度
PR=a*PF (1-a)*PD;     %替换概率

目标函数代码如下：

代码语言：javascript复制
function y = imF(x)
y = (cos(x(1)^2-x(2)^2)-3)/((2 (x(1)^2 x(2)^2))^2) 0.8;
end

目标函数最小值计算代码如下：

代码语言：javascript复制
clear all
clc
x=zeros(1,10);
[x1,x2,f] = PSO_im(@imF,60,2,2,0.8,800,5,0.0000001,10,0.6,0.0000000000000000001,0);
% 得到出计算结果
disp('*************************************************');
disp('目标函数取最小值时的自变量：');
x1
x2
disp('目标函数的最小值为：')
f
disp('**************************************************');

结果如下所示：

代码语言：javascript复制
*************************************************
目标函数取最小值时的自变量：
x1 =
-9.576147568073508e-09
x2 =
4.596695783685527e-09
目标函数的最小值为：
f =
0.300000000000000
**************************************************

3 粒子群算法的权重控制

惯性权重控制前一变化量对当前变化量的影响。 w w w 较大，全局搜索能力较强； w w w 较小，局部搜索能力较强。常见的PSO算法有自适应权重法、随机权重法、线性递减权重法等。

3.1 线性递减法

针对PSO算法容易早熟及后期容易在全局最优解附近产生振荡的现象，提出了线性递减权重法。即惯性权重依照线性从大到小递减，其变化公式为 w = w m a x − t ∗ ( w m a x − w m i n ) t m a x w = w_{max}-frac{t * (w_{max}-w_{min})}{t_{max}} w=wmax​−tmax​t∗(wmax​−wmin​)​ 其中， w m a x w_{max} wmax​ 表示惯性权重最大值， w m i n w_{min} wmin​ 表示惯性权重最小值， t t t 表示当前迭代步数。

3.2 自适应法

3.2.1 根据全局最优点距离进行调整

目前大多采用非线性动态惯性权重系数公式，如下： w = { w m i n − ( w m a x − w m i n ) ∗ ( f − f m i n ) f a v g − f m i n , f ≤ f a v g w m a x , f > f a v g w=left{ begin{matrix} w_{min}-frac{(w_{max}-w_{min})*(f-f_{min})}{f_{avg}-f_{min}},quad f leq f_{avg}\ \ w_{max},quad f>f_{avg}\ end{matrix} right. w=⎩⎨⎧​wmin​−favg​−fmin​(wmax​−wmin​)∗(f−fmin​)​,f≤favg​wmax​,f>favg​​ 其中， f f f 表示粒子实时的目标函数值， f a v g f_{avg} favg​ 和 f m i n f_{min} fmin​ 分别表示当前所有粒子的平均值和最小目标值。

从上面公式可以看出,惯性权重随着粒子目标函数值的改变而改变。当粒子目标值分散时，减小惯性权重；粒子目标值一致时，增加惯性权重。

3.2.2 依据早熟收敛程度和适应值进行调整权重

根据群里的早熟收敛程度和个体适应值，可以确定惯性权重的变化。

设定粒子 p i p_i pi​ 的适应值为 f i f_i fi​，最优粒子适应度为 f m i n f_{min} fmin​，则粒子群的平均适应值是 f a v g = 1 n ∑ i = 1 n f i f_{avg}=frac{1}{n} sum_{i=1}^{n}f_i favg​=n1​∑i=1n​fi​，将优于平均适应值的粒子适应值求平均（记为 f a v g ′ f_{avg}^{‘} favg′​），定义 Δ = ∣ f m − f a v g ′ ∣ Delta=|f_m-f_{avg}^{‘}| Δ=∣fm​−favg′​∣。

依据 f i f_i fi​、 f m f_m fm​、 f a v g f_{avg} favg​ 将群体分为 3 3 3 个子群，分别进行不同的自适应操作。其惯性权重的调整如下：

（1）如果 f i f_i fi​ 优于 f a v g ′ f_{avg}^{‘} favg′​，那么 w = w − ( w − w m i n ) ∗ ∣ f i − f a v g ′ f m − f a v g ′ ∣ w=w-(w-w_{min})*|frac{f_i-f_{avg}^{‘}}{f_m-f_{avg}^{‘}}| w=w−(w−wmin​)∗∣fm​−favg′​fi​−favg′​​∣

（2）如果 f i f_i fi​ 优于 f a v g ′ f_{avg}^{‘} favg′​，且次于 f m f_m fm​，则惯性权重不变。

（3）如果 f i f_i fi​ 次于 f a v g ′ f_{avg}^{‘} favg′​，那么 w = 1.5 − 1 1 k 1 ∗ e x p ( − k 2 ∗ Δ ) w=1.5-frac{1}{1 k_1*exp(-k_2*Delta)} w=1.5−1 k1​∗exp(−k2​∗Δ)1​。

其中， k 1 k_1 k1​、 k 2 k_2 k2​ 为控制参数， k 1 k_1 k1​ 用来控制 w w w 的上限， k 2 k_2 k2​ 主要用来控制 w = 1.5 − 1 1 k 1 ∗ e x p ( − k 2 ∗ Δ ) w=1.5-frac{1}{1 k_1*exp(-k_2*Delta)} w=1.5−1 k1​∗exp(−k2​∗Δ)1​ 的调节能力。

当算法停止时，如果粒子的分布分散，则 Δ Delta Δ 比较大， w w w 变小，此时算法局部搜索能力加强，从而使得群体趋于收敛；若粒子的分布聚集，则 Δ Delta Δ 比较小， w w w 变大，使粒子具有较强的探查能力，从而有效地跳出局部最优。

4 混合粒子群算法

混合策略PSO就是将其他进化算法或传统优化算法或其他技术应用到PSO中，用于提高局部开发能力、增强收敛速度与精度，或者提高粒子多样性、增强粒子地全局探索能力。包括基于模拟退火的混合粒子群算法、基于杂交的混合粒子群算法等。下面以基于的混合粒子群算法为例。

基于的混合粒子群算法是借鉴遗传算法中杂交的概念，在每次迭代中，根据杂交率选取指定数量的粒子放入杂交池内，池内的粒子随机两两杂交，产生同样数目的子代粒子( n n n)，并用子代粒子替代父代粒子( m m m)。子代位置由父代位置进行交叉得到 n x = i ∗ m x ( 1 ) ( 1 − i ) ∗ m x ( 2 ) nx=i*mx(1) (1-i)*mx(2) nx=i∗mx(1) (1−i)∗mx(2) 其中， m x mx mx 表示父代粒子的位置， n x nx nx 表示子代粒子的位置， i i i 是 0 0 0 到 1 1 1 之间的随机数。子代的速度由下式计算： n v = m v ( 1 ) m v ( 2 ) ∣ m v ( 1 ) m v ( 2 ) ∣ ∣ m v ∣ nv=frac{mv(1) mv(2)}{|mv(1) mv(2)|}|mv| nv=∣mv(1) mv(2)∣mv(1) mv(2)​∣mv∣ 其中， m v mv mv 表示父代粒子的速度， n v nv nv 表示子代粒子的速度。

参考文献

[1] MATLAB优化算法/科学与工程计算技术丛书

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/183169.html原文链接：https://javaforall.cn

编程算法 matlab

0 人点赞