探寻“群魔乱舞”的优化秘密！[智能优化算法]——粒子群算法

最新推荐文章于 2024-10-02 23:40:03 发布

涵娃子

最新推荐文章于 2024-10-02 23:40:03 发布

阅读量828

点赞数 15

文章标签：算法 matlab 程序人生经验分享

本文链接：https://blog.csdn.net/yzhh1124/article/details/141564068

版权

/写在前面/

这篇文章和我的上一篇文章一样都是总结自大二下学期的智能优化算法课程。整个代码由我自己编写。期间，老师也给我提供了很多思考的方向。如果同学们有其他的思路或者想法，欢迎在评论区交流哦！

个人能力有限，以下内容仅供参考！

一.算法简介

1.群智能（Swarm Intelligence，SI）

群智能是一种基于群体行为模式和集体智慧的智能系统。它是由多个相对简单的个体组成的群体，通过个体间的相互作用和局部规则，展现出复杂而智能的群体行为。这些个体可以是实体，如机器人、计算机或其他物理系统，也可以是抽象概念，如算法或数据模型。

群智能的主要特点是集体智慧，即群体作为一个整体能够展现出超越个体能力的智能表现。个体间的相互作用和局部规则使得群体能够自我组织、学习和适应环境，解决复杂的问题。这种智能系统通过模拟自然界中的群体行为，如鸟群飞行、蚂蚁搬运食物等，来启发人们设计和构建智能算法。

图1 鸟群觅食

图2 鱼群觅食

群智能的应用领域非常广泛，包括优化问题求解、预测模型、协同工作、社交网络分析、智能控制等。例如，在优化问题求解中，通过模拟群体的竞争和合作行为，群体算法可以找到全局最优解；在协同工作中，个体之间的信息交流和合作可以帮助完成复杂的任务；在社交网络中，可以通过模拟社交网络中的信息传播机制进行舆情分析或舆情控制等任务。

2.算法历史和发展

PSO算法最初是由学者Kennedy和Eberhart在1995年提出的。当时，他们观察到鸟群捕食行为的模式，并从中获得了灵感，试图将这种模式应用于计算机上的优化问题求解。在提出之初，PSO算法主要应用于连续空间的优化问题求解，并逐渐在多个领域展现出较好的效果。

随着研究的发展和应用需求的扩展，PSO算法开始与其他技术相结合，如与神经网络结合应用于分类和预测问题、与模糊逻辑结合用于解决不确定问题等。此外，PSO算法也逐渐扩展到离散优化问题求解领域，如旅行商问题、作业车间调度问题等。这些扩展和与其他技术的结合进一步增强了PSO算法的适应性和解决问题的能力。

除了算法的改进和发展外，PSO算法的应用范围也在不断扩大。它已被广泛应用于机器学习、数据挖掘、图像处理、工程优化等领域。随着人工智能和大数据技术的快速发展，PSO算法也在不断更新和完善，如自适应粒子群优化算法的出现进一步提高了算法的收敛速度和全局搜索能力。

图3 Kennedy和Eberhart

3.算法适用范围

PSO算法应用范围较广，可用于约束优化问题、多目标优化问题、离散空间组合优化问题以及动态跟踪优化问题的求解。同时它也适用于模式识别和图像处理。

在课堂上，我们主要运用PSO算法解决了高维函数优化问题。例如Rosenbrock函数优化问题，在这类问题上，PSO算法展示出了强大的求解能力。

二.算法思想

1.思想概述

既然PSO算法是从自然界获取的灵感并发明出来的，那么我们就以鸟群觅食作为切入点开始理解PSO。

已知有一群鸟在一块区域当中寻找食物。在这块区域中只有一块食物，并且所有的鸟都不知道食物在哪里，但是它们能感受到当前位置距离食物还有多远。

此时，若它们搜寻目前距离食物最近的鸟的周围区域，并根据自己的飞行经验去判断食物所在。那么它们能够较快地找到食物。

这里便体现出了PSO的特点，即个人与社会的相互作用。

2.数学表达

具体而言，在一个 $D$ 维空间中，有 $N$ 个粒子在进行搜寻。每个粒子都有一个对应的速度 $v_{i}=(v_{i1},v_{i2},v_{i3},...v_{iD})$ 和位置 $x_{i}=(x_{i1},x_{i2},x_{i3},...x_{iD})$ ，在我看来，这里的速度更像是一个带有方向的步长。此处的速度和位置都有相应的范围。第 $d(1\leqslant d\leq D)$ 维的位置变化范围： $[X_{min},X_{max}]$ ，速度变化范围： $[-V_{max},V_{max}]$ 。

在搜寻目标时，每个粒子会综合考虑自己的最佳位置 $pbest_{i}=(p_{i1},p_{i2},p_{i3},...p_{iD})$ 和种群的最佳位置 $gbest_{i}=(g_{i1},g_{i2},g_{i3},...g_{iD})$ 去调整搜寻速度。

每个粒子都根据同一个适应度函数来评估自己当前位置的好坏。

粒子 $i$ 的第 $d$ 维速度更新公式：

$v_{id}^{k} = w*v_{id}^{k-1}+c_{1}*r_{1}*(pbest_{id}-x_{id}^{k-1}) + c_{2}*r_{2}*(gbest_{id}-x_{id}^{k-1})$

$v_{id}^{k}$ ：第 $k$ 次迭代粒子 $i$ 搜寻速度的第d维分量。

$x_{id}^{k}$ ：第 $k$ 次迭代粒子 $i$ 搜寻位置的第d维分量。

$c_{1},c_{2}$ ：加速度常数，调节学习最大步长。

$r_{1},r_{2}$ ：两个随机数，取值范围[0,1]，以增加搜索随机性。

$w$ ：惯性权重，非负数。

从速度更新公式可以看出，粒子的速度综合考虑了上一次的速度、目前位置相对于个体最佳的距离和目前位置相对于种群最佳的距离。该公式也体现出了PSO算法中个体与集体的相互作用。

图4 速度更新公式

粒子 $i$ 的第 $d$ 维位置更新公式：

$x_{id}^{k} = x_{id}^{k-1} +v _{id}^{k-1}$

3.算法流程

图5 算法流程图

三.案例分析（PID控制器优化设计）

1.题目要求

在Simulink环境下建立如下模型，利用PSO算法对该PID控制器进行参数优化，即找到一组合适的参数 $K_{p} K_{i} K_{d}$ 使误差性能指标达到最优。

图6 Simulink模型

2.解题思路

该题我们需要确定3个参数，即 $K_{p} K_{i} K_{d}$ ，因此PSO算法解的维数便是3。在MATLAB脚本文件中我们需要进行PSO算法，根据误差性能指标更新粒子的速度、位置、gbest、pbest。而在Simulink中，我们需要根据PSO算法计算的PID参数进行仿真得到误差性能指标。

图7 解题过程示意图

3.解题难点

因为该题PSO算法部分的适应度函数是Simulink模型中ITSE性能指标，而Simulink模型中PID控制器的PID参数又源于PSO计算之后的解 $x_{i} = (x_{i1},x_{i2},x_{i3})$ 其中 $(x_{i1},x_{i2},x_{i3})$ 分别对应 $(k_{p},k_{i},k_{d})$ ，那么MATLAB中的代码（.m脚本文件）如何与Simulink模型进行数据传输呢？

解答：

首先解决Simulink模型向MATLAB中的代码（.m脚本文件）单向传输数据。

在Simulink模型部分，我使用To Workspace模块输出误差性能指标。

To Workspace 模块将输入信号数据写入到工作区。在仿真期间，模块将数据写入到内部缓冲区。暂停仿真或仿真完成后，该数据将写入到工作区。在仿真暂停或停止之前，数据不可用。

图8 To Workspace 模块

在MATLAB代码（.m脚本文件)中我们只需要用两行代码即可获取Simulink模型传输过来的数据。

Simulink_Progress = sim("Simulink_Code1");%双引号内填写自己的Simulink模型文件名
Fitness(i) = Simulink_Progress.simout.Data(end);

此处sim() 使用现有模型配置参数对指定模型进行仿真。

图9 To Workspace模块参数

在这里我们的To Workspace模块变量名为simout，因此我们可以通过Simulink_Progress.simout.Data对To Workspace模块的数值进行输出。

在命令行窗口展示如下：

>> Simulink_Progress = sim("Simulink_Code1")

Simulink_Progress = 

  Simulink.SimulationOutput:

                 simout: [1x1 timeseries] 
                   tout: [126x1 double] 

     SimulationMetadata: [1x1 Simulink.SimulationMetadata] 
           ErrorMessage: [0x0 char] 

>> Simulink_Progress.simout
  timeseries

  常见属性:
            Name: ''
            Time: [126x1 double]
        TimeInfo: [1x1 tsdata.timemetadata]
            Data: [126x1 double]
        DataInfo: [1x1 tsdata.datametadata]

  更多属性, 方法

>> Simulink_Progress.simout.Data

ans =

                   0
                   0
                   0
                   0
                   0
                   0
                   0
   0.000000000000014
   0.000019531571930
   0.000039063516017
   0.000058595840123
   0.000078128545472

接着，我们解决MATLAB中的代码（.m脚本文件）向Simulink模型单向传输数据。

同样，我仅用三行代码即可解决这个问题。

mublockhandle = getSimulinkBlockHandle('Simulink_Code1/PID Controller',true);
%从模块路径中获取模块句柄，这里获取Simulink_Code1文件中的PID Controller这个模块的句柄。
P_Val = num2str(Kp);
%将要传输的数据转换为字符串。
set_param(mublockhandle,'P',P_Val);
%修改PID Controller模块中的参数P的值。

在MATLAB命令行中输入以下代码可以输出Simulink模型中所有使用的块的名称。

blockpaths = find_system('Simulink_Code1','Type','Block')%Simulink_Code1改为自己的Simulink文件名称。

4.部分代码

（1）设置相关参数

参数设置部分主要完成速度、位置的范围、学习因子、惯性因子和迭代次数的设置，以及个体最佳和种群最佳的空间创建。

%% 设置参数
D = 3;%维数
N = 10;%种群规模 z
X_Scale = [0,40];%X的范围
XI_Scale = [0,40];
V_Scale = [-2,2];%V的范围
W = 1;
C1 = 2;
C2 = 2;
P_best_F = zeros(N,1);
P_best_X = zeros(N,D);
G_best_F = 0;
G_best_X = zeros(1,D);
LOOP_Count = 50; %50
His_G_best_F = zeros(1,LOOP_Count);
His_G_best_X = zeros(LOOP_Count,3);
His_Kpid = zeros(LOOP_Count,3);
His_Fitness = zeros(N,LOOP_Count);
mublockhandle = getSimulinkBlockHandle('Simulink_Code1/PID Controller',true);

（2）生成初始解并计算适应度

在给定范围内，生成速度和位置的初始解，给PID模块赋值后，启用Simulink仿真，得到适应度并返回到工作区。

%% 生成初始解
[X,V] = Generate_Solution (D,N,X_Scale,V_Scale,XI_Scale);
%% 计算初始解的适应度
Fitness = zeros(N,1);
for i=1:N
    Kp = X(i,1);
    Ki = X(i,2);
    Kd = X(i,3);
    P_Val = num2str(Kp);
    I_Val = num2str(Ki);
    D_Val = num2str(Kd);
    set_param(mublockhandle,'P',P_Val);
    set_param(mublockhandle,'I',I_Val);
    set_param(mublockhandle,'D',D_Val);
    
    Simulink_Progress = sim("Simulink_Code1");
    Fitness(i) = Simulink_Progress.simout.Data(end);
end
% Fitness = Calculate_Fitness(X,D,N)
P_best_F = Fitness;
P_best_X = X;
[a,b]=find(P_best_F==min(P_best_F));
G_best_F = min(Fitness);
G_best_X = P_best_X(a,:);

（3）更新速度

在速度更新部分综合考虑了自我认识和社会经验部分。

function [New_V] = Change_Speed(V,D,N,V_Scale,W,C1,C2,P_best_X,G_best_X,X)
New_V = rand(N,D);
r1 = rand();
r2 = rand();
for i=1:N
    for j=1:D
        New_V(i,j) = W*V(i,j) + r1*C1*(P_best_X(i,j) - X(i,j)) + r2*C2*(G_best_X(j) - X(i,j));
        if  New_V(i,j) >= V_Scale(2)
            New_V(i,j) = V_Scale(2);
        else if New_V(i,j) <= V_Scale(1)
                New_V(i,j) = V_Scale(1);
            end
        end
    end
end
end

（4）更新位置

function [New_X] = Change_Location(X,V,D,N,X_Scale,XI_Scale)
New_X = zeros(N,D);
for i=1:N
    for j=1:D
        New_X(i,j) = X(i,j) + V(i,j);
        if j==2
            if  New_X(i,j) >= XI_Scale(2)
                New_X(i,j) = XI_Scale(2);
            else if New_X(i,j) <= XI_Scale(1)
                    New_X(i,j) = XI_Scale(1);
                end
            end
        else
            if  New_X(i,j) >= X_Scale(2)
                New_X(i,j) = X_Scale(2);
            else if New_X(i,j) <= X_Scale(1)
                    New_X(i,j) = X_Scale(1);
                end
            end 
        end
    end
end
end

（5）更新gbest和pbest

function [New_G_best_F,New_G_best_X] = Change_G_best(G_best_F,G_best_X,P_best_F,P_best_X)
Min_P_best_F = min(P_best_F);
[r,~] = find(P_best_F == Min_P_best_F);
if Min_P_best_F < G_best_F
    New_G_best_F = Min_P_best_F;
    New_G_best_X = P_best_X(r,:);
else
    New_G_best_F = G_best_F;
    New_G_best_X = G_best_X;
end
end

function [New_P_best_F,New_P_best_X] = Change_P_best(P_best_F,P_best_X,Fitness,D,N,X)
New_P_best_F = zeros(N,1);
New_P_best_X = zeros(N,D);
for i=1:N
    if Fitness(i) < P_best_F(i)
        New_P_best_F(i) = Fitness(i);
        New_P_best_X(i,:) = X(i,:);
    else
        New_P_best_F(i) = P_best_F(i);
        New_P_best_X(i,:) = P_best_X(i,:);
    end
end
end

（6）迭代部分

for n=1:LOOP_Count
    %% 计算适应度
    for i=1:N
        Kp = X(i,1);
        Ki = X(i,2);
        Kd = X(i,3);
        P_Val = num2str(Kp);
        I_Val = num2str(Ki);
        D_Val = num2str(Kd);
        set_param(mublockhandle,'P',P_Val);
        set_param(mublockhandle,'I',I_Val);
        set_param(mublockhandle,'D',D_Val);
        
        Simulink_Progress = sim("Simulink_Code1");
        Fitness(i) = Simulink_Progress.simout.Data(end);
    end
    %% 更新P_best_F、P_best_X
    [P_best_F,P_best_X] = Change_P_best(P_best_F,P_best_X,Fitness,D,N,X);
    %% 更新G_best_F、G_best_X
    [G_best_F,G_best_X] = Change_G_best(G_best_F,G_best_X,P_best_F,P_best_X);
    %% 更新速度
    V = Change_Speed(V,D,N,V_Scale,W,C1,C2,P_best_X,G_best_X,X);
    %% 更新位置
    X = Change_Location(X,V,D,N,X_Scale,XI_Scale); 
    %% 记录G_best_F
    His_G_best_F(1,n) = G_best_F;
    %% 记录G_best_X
    for i=1:D
        His_G_best_X(n,i) = G_best_X(1,i);
    end
    %% 记录每代个体适应度
    His_Fitness(:,n) = Fitness(:,1);
end