pyswarm粒子群算法python解决约束规划 粒子群算法怎么优化

文章目录

• 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\}$

2.1 算法的基本程序

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

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，其调用格式如下：

[xm, dv] = PSO(fitness, N, c1, c2, w, M, D)

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

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

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$ 对环境的适应程度，分为针对被优化目标函数优化行适应度和 针对约束函数的约束型适应度。粒子群算法使用的适应度函数多样， $GW$ 函数和 $RA$

• 优化型适应度
  $F_{obj}(X)=f(x)$
• 约束型适应度
  $F_{i}(X)=\left\{ \begin{matrix} 0,\qquad g_i(X)\leq0\\ \\ g_i(X),\quad g_i(X)\geq0\\ \end{matrix} \right.$

示例

利用基本粒子群算法求解函数 $f(x)=\sum_{i=1}^{10}(x_i^2+2x_i-3)$

解析

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

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

目标函数代码如下：

function F = fitness(x)
F = 0;
for i = 1: 10
    F = F + x(i)^2 + 2 * x(i) - 3;
end

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

• 迭代次数为100

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

结果如下：

***********************
目标函数取最小值时自变量：

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

clear all
clc
x  = zeros(1,10);
[xm,fv] = PSO(@fitness,40,1.2,2.2,0.6,300,10);	% 迭代次数为300
% 取自变量
xm;
% 取函数最小值
fv;

结果如下：

***********************
目标函数取最小值时自变量：

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)=\frac{cos\sqrt{x_1^2-x_2^2}-3}{[2+(x_1^2+x_2^2)]^2}+0.8$ 最小值，其中 $-10 \leq x_i \leq 10$ ，粒子数为 $60$，学习因子均为 $1.2$，退火常数为 $0.8$，迭代次数为 $800$。

解析

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

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;     %替换概率

目标函数代码如下：

function y = imF(x)
	y = (cos(x(1)^2-x(2)^2)-3)/((2+(x(1)^2+x(2)^2))^2)+0.8;
end

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

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('**************************************************');

结果如下所示：

*************************************************
目标函数取最小值时的自变量：

x1 =

    -9.576147568073508e-09


x2 =

     4.596695783685527e-09

目标函数的最小值为：

f =

   0.300000000000000

**************************************************

3 粒子群算法的权重控制

惯性权重控制前一变化量对当前变化量的影响。$w$ 较大，全局搜索能力较强；$w$

3.1 线性递减法

针对PSO算法容易早熟及后期容易在全局最优解附近产生振荡的现象，提出了线性递减权重法。即惯性权重依照线性从大到小递减，其变化公式为
$w = w_{max}-\frac{t * (w_{max}-w_{min})}{t_{max}}$
其中，$w_{max}$ 表示惯性权重最大值，$w_{min}$ 表示惯性权重最小值，$t$

3.2 自适应法

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

目前大多采用非线性动态惯性权重系数公式，如下：
$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.$
其中，$f$ 表示粒子实时的目标函数值，$f_{avg}$ 和 $f_{min}$

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

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

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

设定粒子 $p_i$ 的适应值为 $f_i$，最优粒子适应度为 $f_{min}$，则粒子群的平均适应值是 $f_{avg}=\frac{1}{n} \sum_{i=1}^{n}f_i$，将优于平均适应值的粒子适应值求平均（记为 $f_{avg}^{$），定义 $\Delta=|f_m-f_{avg}^{$。

依据 $f_i$、$f_m$、$f_{avg}$ 将群体分为 $3$

（1）如果 $f_i$ 优于 $f_{avg}^{$，那么 $w=w-(w-w_{min})*|\frac{f_i-f_{avg}^{$

（2）如果 $f_i$ 优于 $f_{avg}^{$，且次于 $f_m$，则惯性权重不变。

（3）如果 $f_i$ 次于 $f_{avg}^{$，那么 $w=1.5-\frac{1}{1+k_1*exp(-k_2*\Delta)}$。

其中，$k_1$、$k_2$ 为控制参数，$k_1$ 用来控制 $w$ 的上限，$k_2$ 主要用来控制 $w=1.5-\frac{1}{1+k_1*exp(-k_2*\Delta)}$

当算法停止时，如果粒子的分布分散，则 $\Delta$ 比较大，$w$ 变小，此时算法局部搜索能力加强，从而使得群体趋于收敛；若粒子的分布聚集，则 $\Delta$ 比较小，$w$

4 混合粒子群算法

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

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

参考文献

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