Matlab新手快速上手2（粒子群算法）

爱打代码的小高

已于 2024-04-21 22:50:03 修改

阅读量1k

点赞数 28

分类专栏： Matlab智能优化算法文章标签：算法 matlab 启发式算法

于 2024-04-21 22:49:26 首次发布

本文链接：https://blog.csdn.net/gr1423428723/article/details/137966782

版权

Matlab智能优化算法专栏收录该内容

7 篇文章 2 订阅

订阅专栏

本文根据一个较为简单的粒子群算法框架详细分析粒子群算法的实现过程，对matlab新手友好，源码在文末给出。

粒子群算法简介

粒子群算法（Particle Swarm Optimization，PSO）是一种群体智能优化算法，灵感来源于鸟群或鱼群等生物群体的行为。在PSO中，每个潜在解被表示为粒子，这些粒子在解空间中移动，并且其位置和速度根据个体经验和群体经验进行更新。

PSO 算法的核心思想是通过模拟群体行为来搜索解空间中的最优解。每个粒子代表了解空间中的一个潜在解，并且根据其当前位置和速度来调整其移动方向。粒子的速度和位置的更新依赖于两个重要的信息：个体最优（局部最优）和全局最优。

PSO 算法的基本步骤如下：

初始化粒子群：随机生成一组粒子，并随机初始化它们的位置和速度。
评估适应度：对每个粒子根据其当前位置计算适应度值。
更新个体最优位置：对于每个粒子，根据其当前位置和个体历史最优位置，更新其个体最优位置。
更新全局最优位置：根据整个粒子群的个体最优位置，更新全局最优位置。
更新粒子速度和位置：根据个体历史最优位置和全局历史最优位置，以及一些随机因素，更新每个粒子的速度和位置。
重复迭代：重复步骤 3~5 直到达到停止条件（例如达到最大迭代次数或找到满意的解）。

PSO 算法的优点包括简单易实现、不需要太多的参数调整、能够处理高维问题以及对于非线性、非凸优化问题有较好的收敛性。然而，它也有一些缺点，例如对于高度多模态函数可能会陷入局部最优、收敛速度慢等。

参数与子函数定义：

function PSO
%---------------------------------- 共性参数 -------------------------------
NP=20;D=2;                    % 种群规模,变量个数
selfw=2.0;globalw=2.0;        % 自身因子,全局因子
w=0.5;maxgen=1000;       % 惯性因子,限定步数
%---------------------------------- 个性参数 -------------------------------
MinX=-65.536;MaxX=65.536; %变量范围
%------------------------------- 粒子位置初始化 ----------------------------
X=MinX+(MaxX-MinX)*rand(NP,D);
F=fun(X);
selfX=X;selfF=F;dX=zeros(NP,D);
%--------------------------------- 适应度统计 ------------------------------
[Bestf,Indexf]=sort(F); globalfi=Bestf(NP);
globalBestX=X(Indexf(NP),:);
function F=fun(X)
a=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;...
   -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];
F=0.002*ones(size(X,1),1);
for j=1:1:25
    F=F+1./(j+(X(:,1)-a(1,j)).^6+(X(:,2)-a(2,j)).^6);%计算适应度
end

如代码中，省略号的意思是将代码延续到下一行，使代码更易读。无实际意义

参数定义

X矩阵表示NP行2列，元素为0-1的随机值的矩阵，每一行表示一个粒子，两个参数分别为粒子的x坐标和y坐标。

F矩阵与计算适应度：

ones函数在上一个遗传算法的文章中已经介绍，F=0.002*ones(size(X,1),1);这个代码生成了一个列数为1，与X矩阵行数相等的矩阵，矩阵元素全为0.002。

下面的for循环就是计算适应度的过程，如代码所示，这个函数根据距离权重也就是j的值（1-25）计算了每个粒子到25个参考点之间的距离之和，得出与最优距离的适应度（0-1）越接近1表示适应度越好，最终目的也就是使用算法找到最结果最趋近1的解，也就是第一个点（-32，-32）的解，为什么是第一个点，带入点的坐标计算就得知，若距离每个点都不相近，那么结果就是一个很大的值除1，就是趋近于0的值，当粒子趋近于第一个点时，结果就趋近1/1也就是1。

zeros函数：

作用：创建一个全为0的矩阵。
用法：Z = zeros(m,n)；表示创建了一个mxn的元素全为0的矩阵。
举例：Z = zreos(3,4);这样就创建了一个3行4列的元素全为0的矩阵

sort函数：

[Bestf,Indexf]=sort(F);这个sort函数是给Bestf排序的，其中排序后的结果存储在 Bestf中，排序后元素在原向量中的索引存储在 Indexf中。例如：

假设有以下向量 F= [10, 30, 20, 50, 40]运行[Bestf,Indexf] = sort(F)
Bestf 存储了排序后的结果：Bestf = [10, 20, 30, 40, 50]。
Indexf 存储了排序后元素在原向量中的索引：Indexf = [1, 3, 2, 5, 4]。

变量定义：

globalfi =Bestf(NP);表示记录粒子与所有目标点的最佳的权重距离之和，最佳值记录在globalfi 中。

globalBestX=X(Indexf(NP),:);表示索引这个最佳粒子并记录在globalBestX中。

a矩阵：

a矩阵存储了15个的参考点坐标，第i个参考点的坐标就是 $(x_{i},y_{i})$ = $(a(1,i),a(2,i))$ ,例如第一个参考点坐标就是（32，32）

size函数：

作用：函数用于获取数组的尺寸。
一个参数：例如：A = [1, 2, 3; 4, 5, 6]; s = size(A);这时输出[2,3]表示这是2行3列的矩阵
两个参数：s = size(A, 1);这样会返回A的行数2，s = size(A, 2);这样会返回A的列数3

算法开始

for gen=1:1:maxgen
    time(gen)=gen;
    %---------------------------- 粒子位置移动 -----------------------------
    for i=1:1:NP
        for j=1:1:D
            dX(i,j)=w*dX(i,j)+selfw*rand*(selfX(i,j)-X(i,j))+...
                globalw*rand*(globalBestX(1,j)-X(i,j));
            X(i,j)=X(i,j)+dX(i,j);    %移动后的位置
            if X(i,j)>MaxX X(i,j)=MaxX;end
            if X(i,j)<MinX X(i,j)=MinX;end
        end
    end
    F=fun(X);
    %----------------------------- 适应度统计 ------------------------------
    [Bestf,Indexf]=sort(F); Bestfi=Bestf(NP);
    BestX=X(Indexf(NP),:);
    %---------------------------- 更新自身最优 -----------------------------
    for i=1:1:NP
        if F(i)>=selfF(i)
            selfF(i)=F(i);selfX(i,:)=X(i,:);
        end
    end
    %---------------------------- 更新全局最优 -----------------------------
    if Bestfi>=globalfi
        globalfi=Bestfi;globalBestX=BestX;
    end
    %----------------------------- 记录结果 --------------------------------
    BestJ(gen)=globalfi;
    if mod(gen,10)==0
        disp(sprintf('当前代数：%d；当前结果：%f',gen,globalfi));
    end
    plot(time,BestJ,'r');axis([1,maxgen,0,1.1]);
    xlabel('迭代步数');ylabel('优化结果');drawnow;pause(0.1);
    if globalfi>1.0 break;end
end
disp(sprintf('迭代步数：%d；优化结果：%f',gen,globalfi));
%---------------------------- 子函数1：目标函数 ----------------------------
disp(F);
disp(BestX);

首先for循环规定迭代次数为1000次，也就是粒子最多行进1000次，time存储for循环的的第几次迭代。

粒子移动：

接下来是粒子移动，这是整个代码的核心部分，首先两个for循环更新所有粒子的x和y坐标，dX用于存储每个粒子的坐标变化包括x和y两个值也就是两个变化方向，可以理解为粒子的变化向量，也就是粒子的行进方向和速度。

惯性因子：

w表示粒子的惯性因子，w*dX(i,j)表示粒子沿上一次的速度与方向行进，w的值为0.5表示速度减慢为一半，后面的算式表示根据自身最优与全局最优改变方向并增加速度，详细实现在下面给出：。

自身最优：

    for i=1:1:NP
        if F(i)>=selfF(i)
            selfF(i)=F(i);selfX(i,:)=X(i,:);
        end
    end

这个表示更新自身最优的过程，selfF表示粒子每次的与最优解的适应度（0-1），越接近1表示越接近最优解，selfX表示记录所有粒子的自身最优解的位置，因此selfw*rand*(selfX(i,j)-X(i,j))就是得出当前粒子坐标与最优坐标的差，并乘上自身因子与随机值selfw*rand，表示速度。

全局最优：

    if Bestfi>=globalfi
        globalfi=Bestfi;globalBestX=BestX;
    end

与上面自身最优同理，根据最好粒子的适应度，更新全局最优解的坐标，其实也就是适应度最好的粒子，同理globalw*rand*(globalBestX(1,j)-X(i,j))也就是向全局最优解方向移动一定的值。

总结：根据惯性因子、自身最优、全局最优实现粒子的更新，惯性因子让粒子保持上一次的移动方向，自身最优与全局最优都是让粒子改变方向，根据权重向着这两个方向偏移一定的角度和速度，最终实现粒子的寻优过程。

源代码：

%---------------------------------- 程序说明 -------------------------------

%                           该程序实现了普通粒子群算法

%---------------------------------- 程序正文 -------------------------------
function PSO
%---------------------------------- 共性参数 -------------------------------
NP=20;D=2;                    % 种群规模,变量个数
selfw=2.0;globalw=2.0;        % 自身因子,全局因子
w=0.5;maxgen=1000;       % 惯性因子,限定步数
%---------------------------------- 个性参数 -------------------------------
MinX=-65.536;MaxX=65.536; %变量范围
%------------------------------- 粒子位置初始化 ----------------------------
X=MinX+(MaxX-MinX)*rand(NP,D);
F=fun(X);
selfX=X;selfF=F;dX=zeros(NP,D);
%--------------------------------- 适应度统计 ------------------------------
[Bestf,Indexf]=sort(F); globalfi=Bestf(NP);
globalBestX=X(Indexf(NP),:);
%------------------------------- 程序主循环开始 ----------------------------
for gen=1:1:maxgen
    time(gen)=gen;
    %---------------------------- 粒子位置移动 -----------------------------
    for i=1:1:NP
        for j=1:1:D
            dX(i,j)=w*dX(i,j)+selfw*rand*(selfX(i,j)-X(i,j))+...
                globalw*rand*(globalBestX(1,j)-X(i,j));
            X(i,j)=X(i,j)+dX(i,j);    %移动后的位置
            if X(i,j)>MaxX X(i,j)=MaxX;end
            if X(i,j)<MinX X(i,j)=MinX;end
        end
    end
    F=fun(X);
    %----------------------------- 适应度统计 ------------------------------
    [Bestf,Indexf]=sort(F); Bestfi=Bestf(NP);
    BestX=X(Indexf(NP),:);
    %---------------------------- 更新自身最优 -----------------------------
    for i=1:1:NP
        if F(i)>=selfF(i)
            selfF(i)=F(i);selfX(i,:)=X(i,:);
        end
    end
    %---------------------------- 更新全局最优 -----------------------------
    if Bestfi>=globalfi
        globalfi=Bestfi;globalBestX=BestX;
    end
    %----------------------------- 记录结果 --------------------------------
    BestJ(gen)=globalfi;
    if mod(gen,10)==0
        disp(sprintf('当前代数：%d；当前结果：%f',gen,globalfi));
    end
    plot(time,BestJ,'r');axis([1,maxgen,0,1.1]);
    xlabel('迭代步数');ylabel('优化结果');drawnow;pause(0.1);
    if globalfi>1.0 break;end
end
disp(sprintf('迭代步数：%d；优化结果：%f',gen,globalfi));
%---------------------------- 子函数1：目标函数 ----------------------------
disp(F);
disp(BestX);
function F=fun(X)
a=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;...
   -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];
F=0.002*ones(size(X,1),1);
for j=1:1:25
    F=F+1./(j+(X(:,1)-a(1,j)).^6+(X(:,2)-a(2,j)).^6);%计算适应度
end

爱打代码的小高

关注

28
点赞
踩
30

收藏

觉得还不错? 一键收藏
0
评论
Matlab新手快速上手2（粒子群算法）

粒子群算法（Particle Swarm Optimization，PSO）是一种群体智能优化算法，灵感来源于鸟群或鱼群等生物群体的行为。本文根据一个较为简单的粒子群算法框架详细分析粒子群算法的实现过程，一步步分析代码实现，对matlab新手友好
复制链接

扫一扫