粒子群优化算法原理_耀世娱乐-耀世注册登录官方入口

粒子群优化（Particle Swarm Optimization,PSO)算法是1995年由美国学者Kennedy等人提出的，该算法是模拟鸟类觅食等群体智能行为的智能优化算法。在自然界中，鸟群在觅食的时候，一般存在个体和群体协同的行为。有时鸟群分散觅食，有时鸟群也全体觅食。在每次觅食的过程中，都会存在一些搜索能力强的鸟，这些搜索能力强的鸟，会给其他鸟传递信息，带领其他鸟到食物源位置。

在粒子群优化算法中，目标空间中的每个解都可以用一只鸟（粒子）表示，问题中的需求解就是鸟群所要寻找的食物源。在寻找最优解的过程中，每个粒子都存在个体行为和群体行为。每个粒子都会学习同伴的飞行经验和借鉴自己的飞行经验去寻找最优解。每个粒子都会向两个值学习，一个值是个体的历史最优值 $p_{best}$ ;另一个值是群体的历史最优值（全局最优值） $g_{best}$ 。粒子会根据这两个值来调整自身的速度和位置，而每个位置的优劣都是根据适应度值来确定的。适应度函数是优化的目标函数。

在一个D维的目标搜索空间中，由N个粒子组成一个粒子群，其中每个粒子都是一个D维向量，其空间位置可以表示为

$x_i=(x_{i1},x_{i2},\\cdots,x_{iD}),i=1,2,\\cdots,N$

粒子的空间位置是目标优化问题中的一个解，将其代入适应度函数可以计算出适应度值，根据适应度值的大小衡量粒子的优劣。

第i个粒子的飞行速度也是一个D维向量，记为

$v_i=(v_{i1},v_{i2},\\cdots,v_{iD}),i=1,2,\\cdots,N$

粒子的位置和速度均值都在给定的范围内随机生成。

第i个粒子经历过的具有最优适应度值的位置称为个体历史最优位置，记为

$p_{best i}=(p_{besti1},p_{besti2},\\cdots,p_{bestiD}),i=1,2,\\cdots,N$

整个粒子群经历过的最优位置称为全局历史最优位置，记为

$g_{best i}=(g_{besti1},g_{besti2},\\cdots,g_{bestiD}),i=1,2,\\cdots,N$

粒子群的位置更新操作可用速度更新和位置更新表示。

速度更新为

$v_{ij}(t+1)=v_{ij}(t)+c_1r_1(p_{bestij}(t)-x_{ij}(t))+c_2r_2(g_{best j}-x_{ij}(t))$

位置更新为

$x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)$

其中，下标j表示粒子的第j维，下标i表示第i个粒子，t表示当前迭代次数，c1与c2均加速常量，通常在区间(0,2)内取值，r1与r2为两个相互独立的取值范围在[0,1]的随机数。从上述方程可以看出，c1与c2将粒子向个体学习和向群体学习联合起来，使得粒子能够借鉴体自身的搜索经验和群体的搜索经验。

粒子群优化算法的流程如下：

步骤1:初始化粒子群参数c1与c2，设置位置边界范围与速度边界范围，初始化粒子群群，初始化粒子群速度。

步骤2:根据适应度函数计算适应度值，记录历史最优值 $p_{best}$ 与全局最优值 $g_{best}$ 。

步骤3:利用速度更新公式 $v_{ij}(t+1)=v_{ij}(t)+c_1r_1(p_{bestij}(t)-x_{ij}(t))+c_2r_2(g_{best j}-x_{ij}(t))$ 对粒子群的速度进行更新，并对越界的速度进行约束。

步骤4:利用位置更新公式 $x_{ij}(t+1)=x_{ij}(t)+v_{ij}(t+1)$ 对粒子群的位置进行更新，并对越界的位置进行约束。

步骤5:根据适应度函数计算适应度值。

步骤6:对于每个粒子，将其适应度值与它的历史最优适应度值相比较，若更好，则将作为历史最优值 $p_{best}$ 。

步骤7:对于每个粒子，比较其适应度值和群体所经历的最优位置的适应度值，若更好，将其作为全局最优值 $p_{best}$

步骤8:判断是否达到结束条件（达到最大达代次数）,若达到，则输出最优位置，否则复步骤3~8。

数rand()是MATLAB自带的随机数生成函数，会生成区间[0,1]内的随机数。

>> rand()

生成多个随机数：

rand(3,4)

若要生成指定范围内的随机数，可以利用如下公式：

$r=lb+(ub-lb)*rand()$

(4-0).*rand(1,5)+0 % 0到4内的随机数5个
ans=列 1 至 2

         0.428760196877498          3.81908193731818

  列 3 至 4

          2.38316963892275         0.828332505066357

  列 5

           3.0342077213501

%% 粒子群初始化函数
function X = initialization(pop,ub,lb,dim)
    %pop:为种群数量
    %dim:每个粒子群的维度
    %ub:为每个维度的变量上边界，维度为[1,dim];
    %lb:为每个维度的变量下边界，维度为[1,dim];
    %X:为输出的种群，维度[pop,dim];
    X = zeros(pop,dim); %为X事先分配空间
    for i = 1:pop
       for j = 1:dim
           X(i,j) = (ub(j) - lb(j))*rand() + lb(j);  %生成[lb,ub]之间的随机数
       end
    end
end

假设种群数量为10，每个粒子群维度为5，每个维度的边界均为[-5,5],利用粒子群初始化函数初始化种群。

pop = 10;
dim =5;
ub=[5,5,5,5,5];
lb=[-5,-5,-5,-5,-5];
x=initialization(pop,ub.lb,dim)

优化问题的目标函数，根据不同应用设计相应的适应度函数。

%% 适应度函数
function fitness = fun(x)
%x为输入一个粒子，维度为[1,dim];
%fitness 为输出的适应度值
    fitness = sum(x.^2);
end

边界检查和约束函数

防止变量超过规定的范围，一般当变量大于边界时，直接将其设置为上边界，当变量小于下边界时，直接将其设置为下边界。

%% 边界检查函数
function [X]= BoundaryCheck(x,ub,lb,dim)
    %dim为数据的维度大小
    %x为输入数据，维度为[1,dim];
    %ub为数据上边界，维度为[1,dim]
    %lb为数据下边界，维度为[1,dim]
    for i = 1:dim
        if x(i)>ub(i)
           x(i) = ub(i); 
        end
        if x(i)<lb(i)
            x(i) = lb(i);
        end
    end
    X = x;
end

%%--------------粒子群函数----------------------%%
%% 输入：
%   pop:种群数量
%   dim:单个粒子的维度
%   ub:粒子上边界信息，维度为[1,dim];
%   lb:粒子下边界信息，维度为[1,dim];
%   fobj:为适应度函数接口
%   vmax: 为速度的上边界信息，维度为[1,dim];
%   vmin: 为速度的下边界信息，维度为[1,dim];
%   maxIter: 算法的最大迭代次数，用于控制算法的停止。
%% 输出：
%   Best_Pos：为粒子群找到的最优位置
%   Best_fitness: 最优位置对应的适应度值
%   IterCure:  用于记录每次迭代的最佳适应度，即后续用来绘制迭代曲线。
function [Best_Pos,Best_fitness,IterCurve]= pso(pop,dim,ub,lb,fobj,vmax,vmin,maxIter)
    %%设置c1,c2参数
    c1 = 2.0;
    c2 = 2.0;
    %% 初始化种群速度
    V = initialization(pop,vmax,vmin,dim);
    %% 初始化种群位置
    X = initialization(pop,ub,lb,dim);

    %% 计算适应度值
    fitness = zeros(1,pop);
    for i = 1:pop
       fitness(i) = fobj(X(i,:));
    end
    %% 将初始种群作为历史最优
    pBest = X;
    pBestFitness = fitness;
    %% 记录初始全局最优解,默认优化最小值。
    %寻找适应度最小的位置
    [~,index] = min(fitness);
    %记录适应度值和位置
    gBestFitness = fitness(index);
    gBest = X(index,:);

    Xnew = X; %新位置
    fitnessNew = fitness;%新位置适应度值

    IterCurve = zeros(1,maxIter);
    %% 开始迭代
    for t = 1:maxIter
       %对每个粒子进行更新 
       for i = 1:pop
          %速度更新
          r1 = rand(1,dim);
          r2 = rand(1,dim);
          V(i,:) = V(i,:) + c1.*r1.*(pBest(i,:) - X(i,:)) + c2.*r2.*(gBest - X(i,:));
          %速度边界检查及约束
          V(i,:) = BoundaryCheck(V(i,:),vmax,vmin,dim);
          %位置更新
          Xnew(i,:) = X(i,:) + V(i,:);
          %位置边界检查及约束
          Xnew(i,:) = BoundaryCheck(Xnew(i,:),ub,lb,dim);
          %计算新位置适应度值
          fitnessNew(i) = fobj(Xnew(i,:));
          %更新历史最优值
          if fitnessNew(i) < pBestFitness(i)
              pBest(i,:) = Xnew(i,:);
              pBestFitness(i) = fitnessNew(i);    
          end
          %更新全局最优值
          if fitnessNew(i)<gBestFitness
              gBestFitness = fitnessNew(i);
              gBest = Xnew(i,:);
          end   
       end
       X = Xnew;
       fitness = fitnessNew;
       %% 记录当前迭代最优值和最优适应度值
       %记录最优解
       Best_Pos = gBest;
       %记录最优解的适应度值
       Best_fitness = gBestFitness;
       %记录当前迭代的最优解适应度值
       IterCurve(t) = gBestFitness;   
    end
end

求解函数的极值， $f(x_1,x_2)=x_1^2+x_2^2,-10\\leq x_1,x_2\\leq10$

绘制函数的曲面图。

x1=-10:0.01:10;
x2=-10:0.01:10;
for i=1:size(x1,2)
    for j =1:size(x2,2)
        x1(i,j)=x1(i);
        x2(i,j)=x2(j);
        f(i,j)=x1(i)^2+x2(j)^2;
    end
end
surfc(x1,x2,f,'LineStyle','none');%绘制搜索曲面

适应度函数

%% 适应度函数
function fitness = fun(x)
%x为输入一个粒子，维度为[1,dim];
%fitness 为输出的适应度值
    fitness = x(1)^2 + x(2)^2;
end

设置种群规模pop=50，最大迭代次数maxIter=100,dim=2,ub=[10,10],下边界lb=[-10,-10]，速度上边界v_max=[2,2],下边界v_min=[-2,-2]。

%% 粒子群算法求解x1^2 + x2^2的最小值
clc;clear all;close all;
%粒子群参数设定
pop = 50;%种群数量
dim = 2;%变量维度
ub = [10,10];%粒子上边界信息
lb = [-10,-10];%粒子下边界信息
vmax = [2,2];%粒子的速度上边界
vmin = [-2,-2];%粒子的速度下边界
maxIter = 100;%最大迭代次数
fobj = @(x) fun(x);%设置适应度函数为fun(x); 
%粒子群求解问题
[Best_Pos,Best_fitness,IterCurve] = pso(pop,dim,ub,lb,fobj,vmax,vmin,maxIter);
%绘制迭代曲线
figure
plot(IterCurve,'r-','linewidth',1.5);
grid on;%网格开
title('粒子群迭代曲线')
xlabel('迭代次数')
ylabel('适应度值')

disp(['求解得到的x1，x2为',num2str(Best_Pos(1)),'   ',num2str(Best_Pos(2))]);
disp(['最优解对应的函数值为：',num2str(Best_fitness)]);

运行结果：

求解得到的x1，x2为0.0016402 -0.0056427

最优解对应的函数值为：3.4531e-05

为了更加直观地了解粒子群在每代的分布、前后迭代、粒子群的位置变化，以及每代中最优粒子的位置，需要将粒子群优化算法的中间结果绘制出来。为了达到此目的，我们需要记录每代粒子群的位置（History Position),同时记录每代最优粒子的位置（History Best),然后通过MATLAB绘图函数，将图像绘制出来。

绘制每代信息的主函数

%% 粒子群算法求解x1^2 + x2^2的最小值
clc;clear all;close all;
%粒子群参数设定
pop = 50;%种群数量
dim = 2;%变量维度
ub = [10,10];%粒子上边界信息
lb = [-10,-10];%粒子下边界信息
vmax = [2,2];%粒子的速度上边界
vmin = [-2,-2];%粒子的速度下边界
maxIter = 50;%最大迭代次数
fobj = @(x) fun(x);%设置适应度函数为fun(x); 
%粒子群求解问题
[Best_Pos,Best_fitness,IterCurve,HistoryPosition,HistoryBest] = pso1(pop,dim,ub,lb,fobj,vmax,vmin,maxIter);
%绘制迭代曲线
figure
plot(IterCurve,'r-','linewidth',1.5);
grid on;%网格开
title('粒子群迭代曲线')
xlabel('迭代次数')
ylabel('适应度值')

disp(['求解得到的x1，x2为',num2str(Best_Pos(1)),'   ',num2str(Best_Pos(2))]);
disp(['最优解对应的函数值为：',num2str(Best_fitness)]);

%% 绘制每一代粒子群分布
for i = 1:maxIter
   Position = HistoryPosition{i};%获取当前代位置
   BestPosition = HistoryBest{i};%获取当前代最佳位置
   figure(3)
   plot(Position(:,1),Position(:,2),'*','linewidth',3); 
   hold on;
    plot(BestPosition(1),BestPosition(2),'ro','linewidth',3); 
   grid on;
   axis([-10 10,-10,10])
   legend('粒子','最佳粒子');
   title(['第',num2str(i),'次迭代']);
   hold off
end

记录每代位置的函数

%%--------------粒子群函数----------------------%%
%% 输入：
%   pop:种群数量
%   dim:单个粒子的维度
%   ub:粒子上边界信息，维度为[1,dim];
%   lb:粒子下边界信息，维度为[1,dim];
%   fobj:为适应度函数接口
%   vmax: 为速度的上边界信息，维度为[1,dim];
%   vmin: 为速度的下边界信息，维度为[1,dim];
%   maxIter: 算法的最大迭代次数，用于控制算法的停止。
%% 输出：
%   Best_Pos：为粒子群找到的最优位置
%   Best_fitness: 最优位置对应的适应度值
%   IterCure:  用于记录每次迭代的最佳适应度，即后续用来绘制迭代曲线。
%   HistoryPosition: 用于记录每代，粒子群的位置
%   HistoryBest: 用于记录每代，粒子群的最佳位置
function [Best_Pos,Best_fitness,IterCurve,HistoryPosition,HistoryBest]= pso1(pop,dim,ub,lb,fobj,vmax,vmin,maxIter)
    %%设置c1,c2参数
    c1 = 2.0;
    c2 = 2.0;
    %% 初始化种群速度
    V = initialization(pop,vmax,vmin,dim);
    %% 初始化种群位置
    X = initialization(pop,ub,lb,dim);

    %% 计算适应度值
    fitness = zeros(1,pop);
    for i = 1:pop
       fitness(i) = fobj(X(i,:));
    end
    %% 将初始种群作为历史最优
    pBest = X;
    pBestFitness = fitness;
    %% 记录初始全局最优解,默认优化最小值。
    %寻找适应度最小的位置
    [~,index] = min(fitness);
    %记录适应度值和位置
    gBestFitness = fitness(index);
    gBest = X(index,:);

    Xnew = X; %新位置
    fitnessNew = fitness;%新位置适应度值

    IterCurve = zeros(1,maxIter);
    %% 开始迭代
    for t = 1:maxIter
       %对每个粒子进行更新 
       for i = 1:pop
          %速度更新
          r1 = rand(1,dim);
          r2 = rand(1,dim);
          V(i,:) = V(i,:) + c1.*r1.*(pBest(i,:) - X(i,:)) + c2.*r2.*(gBest - X(i,:));
          %速度边界检查及约束
          V(i,:) = BoundaryCheck(V(i,:),vmax,vmin,dim);
          %位置更新
          Xnew(i,:) = X(i,:) + V(i,:);
          %位置边界检查及约束
          Xnew(i,:) = BoundaryCheck(Xnew(i,:),ub,lb,dim);
          %计算新位置适应度值
          fitnessNew(i) = fobj(Xnew(i,:));
          %更新历史最优值
          if fitnessNew(i) < pBestFitness(i)
              pBest(i,:) = Xnew(i,:);
              pBestFitness(i) = fitnessNew(i);    
          end
          %更新全局最优值
          if fitnessNew(i)<gBestFitness
              gBestFitness = fitnessNew(i);
              gBest = Xnew(i,:);
          end   
       end
       X = Xnew;
       fitness = fitnessNew;
       %% 记录当前迭代最优值和最优适应度值
       %记录最优解
       Best_Pos = gBest;
       %记录最优解的适应度值
       Best_fitness = gBestFitness;
       %记录当前迭代的最优解适应度值
       IterCurve(t) = gBestFitness;  
       HistoryBest{t} = Best_Pos;
       %记录当前代粒子群的位置
       HistoryPosition{t} = X;
       
    end
end

$\\min f(x)=0.6224x_1x_3x_4+1.7781x_2x_3^2+3.1661x_1^2x_4+19.84x_1^2x_3$

约束条件：

$g_1(x)=-x_1+0.0193x_3\\leq0$

$g_2(x)=-x_2+0.00954x_3\\leq0$

$g_3(x)=-\\pi x_3^2-4\\pi x_3^3/3+129600 \\leq0$

$g_4(x)=x_4-240\\leq0$

$0\\leq x_1,x_2,x_3,x_4\\leq100$

% 有约束的适应度函数
function fitness = fun(x)
    x1 = x(1); %Ts
    x2 = x(2); %Th
    x3 = x(3); %R
    x4 = x(4); %L

    %% 约束条件判断
    g1 = -x1+0.0193*x3;
    g2=-x2+0.00954*x3;
    g3=-pi*x3^2-4*pi*x3^3/3+1296000;
    g4=x4-240;
    if(g1<=0&&g2<=0&&g3<=0&&g4<=0)%如果满足约束条件则计算适应度值
       fitness = 0.6224*x1*x3*x4 + 1.7781*x2*x3^2 + 3.1661*x1^2*x4 + 19.84*x1^2*x3;
    else%否则适应度值无效
        fitness = inf;
    end

end

%% 基于粒子群算法的压力容器设计
clc;clear all;close all;
%粒子群参数设定
pop = 50;%种群数量
dim = 4;%变量维度
ub =[100,100,100,100];%粒子上边界信息
lb = [0,0,10,10];%粒子下边界信息
vmax = [2,2,2,2];%粒子的速度上边界
vmin = [-2,-2,-2,-2];%粒子的速度下边界
maxIter = 500;%最大迭代次数
fobj = @(x) fun1(x);%设置适应度函数为fun1(x); 
%粒子群求解问题
[Best_Pos,Best_fitness,IterCurve] = pso(pop,dim,ub,lb,fobj,vmax,vmin,maxIter);
%绘制迭代曲线
figure
plot(IterCurve,'r-','linewidth',1.5);
grid on;%网格开
title('粒子群迭代曲线')
xlabel('迭代次数')
ylabel('适应度值')

disp(['求解得到的x1,x2,x3,x4为:',num2str(Best_Pos(1)),'   ',num2str(Best_Pos(2)),' ',num2str(Best_Pos(3)),' ',num2str(Best_Pos(4))]);
disp(['最优解对应的函数值为：',num2str(Best_fitness)]);

运行结果：

求解得到的x1,x2,x3,x4为:1.3363 0.68178 68.6181 28.1582

最优解对应的函数值为：9905.3211

智能优化算法及其MATLAB实现，陈克伟，电子工业出版社，2022年1月