粒子群算法matlab代码（注释很详细哦，图像也美美哒，任意维度）

2022.07.09更新：
更简洁的实现请点击这里，有时候链接会挂掉，在这句话上方↑数学建模专栏里有。

原博文：

答疑：
我把答疑放在最前面了。很多小伙伴真的是零基础，既不懂matlab也不懂粒子群算法。因此我把最常问我的问题总结一下。

---

Q1：为什么我把直接把代码复制过去不能运行？
A1：因为这个代码没有目标函数，你需要设置求解的函数才能计算吧。下文有推荐的测试函数，copy过来可以直接用。

---

Q2：为什么设置了目标函数还是报错？
A2：设置目标函数后，需要针对目标更改参数。最起码你得把维度数和求解范围设置对（pso1_mian.m的第16行和第18行）。不然就很可能不能运行。

---

Q3：为什么每次运行的结果不一样？
A3：进化算法都会有这种现象，尤其是在目标函数很复杂的情况下，可能每次的结果差异都很大。调整参数也许能缓解这种情况，但是需要你对进化算法有很深的理解。

---

Q4：为什么我运行后没有函数图像，只有收敛过程图像？
A4：只有2维和3维的目标函数能画图，高维是没法画图的。

以下是正文：
如果你还不懂粒子群的原理，可以参考我这个博文：初探粒子群优化算法

整个程序分为5个脚本

• pso1_mian.m：主程序，在此脚本内设置参数。
• pso1_im.m：画出函数图像（仅1维和2维）
• pso1_in.m:初始化
• pso1_in2.m:迭代寻优并输出结果
• 另外还有一个目标函数，单独为一个脚本。

• 推荐的测试函数—>这里

• 支持任意维度目标函数的寻优

• 支持动态调整参数（惯性因子随迭代减小，$c_1,c_2$的比值随迭代增大）

• 可视化粒子运动位置

---

先上运行结果图

---

下面是源码

1.pso1_mian.m
这里的目标函数用函数句柄的形式调用（第15行）

%% 粒子群算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pso1_im.m：画出函数图像（仅1维和2维）。
% pso1_in.m:初始化
% pso1_in2.m:迭代寻优并输出结果
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 可支持N维函数
% 可使用动态参数算法
% 多维目标函数的自变量输入格式为行向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc,clear;close all;

% file = 'E:\日常工作\学习笔记\matlab代码\粒子群优化算法\pso1\数据shubert1.xlsx';
funct = @langer;                          %目标函数
fun_range = [0,3.08;0,2.1];    %每个维度的求解范围，行代表维度，第一列最小值，第二列最大值
density = 0.01;                         %图像中取点的密度
d = 2;                                  %维度
popsize = 10;                           %种群大小
maxgen = 100;                           %最大迭代次数
limit_v = 10;                          %每个维度的最大运动速度
maxormin = 0;                           %求最大值还是最小值。最大值为1，最小值为0
dynamic = 1;                            %是否采用动态参数
W = [0.7,1.4];                          %惯性因子。动态参数时为行向量，第一位是最小值，第二位是最大值
C2 = [0.5,2.5];
C1 = [0.5,2.5];                         %学习因子。动态参数时为行向量，第一位是最小值，第二位是最大值


Z = pso1_im(funct,fun_range,density,d);
[x,v,fitness,pbest,gbest,gbest_fitness,gbest_index,pbest_fitness] =...
   pso1_in(funct,popsize,d,fun_range,limit_v,maxormin);
[gbest_fitness,gbest_index,record] =...
    pso1_in2(funct,d,popsize,maxgen,dynamic,W,C1,C2,fun_range,limit_v,density...
    ,x,v,fitness,pbest,gbest,gbest_fitness,gbest_index,pbest_fitness,maxormin);

% canshu = {'惯性因子',W,'种群规模',popsize,'最大迭代次数',maxgen,'速度限制',limit_v,...
%     '最大值还是最小值',maxormin};
% xlswrite(file,record,'record');
% xlswrite(file,canshu,'参数');

2. pso1_im.m

 function z = pso1_im(funct,fun_range,density,D)
%PSO1_IM 显示目标函数的图像（1维或2维）
%
%INPUT:
%   funct为目标函数
%   fun_range为求解区域矩阵，行代表维度，第一列最小值，第二列最大值
%   density为取值密度，为一个数字
%   D为空间维度，为正整数
if D == 2
    figure(1)
    [x1,x2] = meshgrid(fun_range(1,1):density:fun_range(1,2),fun_range(2,1):density:fun_range(2,2));
    for i = 1 : size(x1,1)
        for j = 1 : size(x1,2)
            z(i,j) = funct([x1(i,j),x2(i,j)]);
        end
    end
    mesh(x1,x2,z);title('目标函数图形');
    hold on;
elseif D == 1
    figure(1)
    x1 = fun_range(1):density:fun_range(2);
    for i = 1 : length(x1)
        z(i) = funct(x1(i));
    end
    plot(x1,z);
    hold on;
else 
    print('高维函数不能画出图像');
end
end

3.pso1_in.m

function [x,v,fitness,pbest,gbest,gbest_fitness,gbest_index,pbest_fitness] = pso1_in (funct,popsize,d,fun_range,limit_v,maxormin)
%PSO1_IN 进行初始化
%
%INPUT：
%   size：种群规模
%   maxgen：最大迭代次数
%   d：维数
%   fun_range：运动范围矩阵，行数为维数，第一列最小值，第二列最大值
%   limit_v:速度范围，为一个数字
%   maxormin：0为求最小值，1为求最大值

%位置和速度初始化
for i = 1 : popsize
    for j = 1 : d
        x(i,j) = (fun_range(j,2) - fun_range(j,1)) *  rand + fun_range(j,1); % 位置初始化
        v(i,j) = rand * limit_v;% 速度初始化
    end
end
for i = 1 : popsize
    fitness(i,:) = funct(x(i,:)); % 计算初始适应度
end

pbest = x;              % 个体最优值位置
pbest_fitness = fitness;% 个体最优值
if maxormin == 1
    [gbest_fitness,gbest_index] = max(fitness);
elseif maxormin == 0
    [gbest_fitness,gbest_index] = min(fitness);
else print('maxormin的值非法')
        return

end
gbest = x(gbest_index,:);

end

4. pso1_in2.m

function [gbest_fitness,gbest_index,record] =...
    pso1_in2(funct,d,popsize,maxgen,dynamic,W,C1,C2,fun_range,limit_v,density...
    ,x,v,fitness,pbest,gbest,gbest_fitness,gbest_index,pbest_fitness,maxormin)
%PSO1_IN2 开始迭代计算
%   
%INPUT:
%   size：种群规模
%   maxgen：最大迭代次数
%   d：维数
%   fun_range：运动范围矩阵，行数为维数，第一列最小值，第二列最大值
%   limit_v:速度范围，为一个数字
%   maxormin：0为求最小值，1为求最大值
times = 1;

for i = 1 : maxgen      % i表示迭代次数
    if dynamic == 1
        w = W(2) - (W(2) - W(1)) * i/maxgen; % 动态调整学习因子
        c1 = (C1(1) - C1(2)) * i / maxgen + C1(2);
        c2 = (C2(2) - C2(1)) * i / maxgen + C2(1);
    elseif dynamic == 0
        w = W;
        c1 = C1;
        c2 = C2;
    end
    for j = 1 : popsize    % j表示粒子索引
        % 更新速度
        v(j,:) = w * v(j,:) + c1 * rand * (pbest(j,:) - x(j,:)) + c2 * rand * (gbest - x(j,:));
        % 限速
    for g=1:d           % g表示维度
        for k=1:popsize    % k表示粒子索引
        if  v(k,g) > limit_v
            v(k,g) = limit_v;
        end
        if  v(k,g) < (-1 * limit_v)
            v(k,g) = (-1 * limit_v);
        end
        end
    end     
        % 更新位置
        x(j,:) = x(j,:) + v(j,:);
        %限制运动范围
    for g=1:d           % g表示维度
        for k=1:popsize    % k表示粒子索引
        if  x(k,g) > fun_range(g,2)
            x(k,g) = fun_range(g,2);
        end
        if  x(k,g) < fun_range(g,1)
            x(k,g) = fun_range(g,1);
        end
        end
    end   
        % 更新适应度
        for u = 1 : popsize
            fitness(u,:) = funct(x(u,:)); % 计算初始适应度
        end
    % 更新全局最优
    if maxormin == 0
        if min(fitness) < gbest_fitness
            [gbest_fitness,gbest_index] = min(fitness);
            gbest = x(gbest_index,:);
        end
    elseif maxormin == 1
        if max(fitness) > gbest_fitness
            [gbest_fitness,gbest_index] = max(fitness);
            gbest = x(gbest_index,:);
        end
    end
    
    % 更新局部最优
    for k = 1 : popsize
      if maxormin == 0
        if fitness(k,:) < pbest_fitness(k,:)
            pbest_fitness(k,:) = fitness(k,:);
            pbest(k,:) = x(k,:);
        end
      elseif maxormin == 1 
        if fitness(k,:) > pbest_fitness(k,:)
            pbest_fitness(k,:) = fitness(k,:);
            pbest(k,:) = x(k,:);
        end
      end
    end
        record(i,:) = gbest_fitness; % record记录每代的全局最优
    end
    
    if d == 2
        if times >= 10
            cla;
            [x1,x2] = meshgrid(fun_range(1,1):density:fun_range(1,2),fun_range(2,1):density:fun_range(2,2));
            for i = 1 : size(x1,1)
                for j = 1 : size(x1,2)
                    z(i,j) = funct([x1(i,j),x2(i,j)]);
                end
            end
            mesh(x1, x2, z,'EdgeColor','none','FaceColor','flat','FaceAlpha','0.4');
            plot3(x(:,1),x(:,2),fitness, 'ro');title('状态位置变化');
            pause(0.5);
            times=0;
        end
    end
    times = times + 1;
end
disp( ['目标函数最小值:' num2str(gbest_fitness)]);
disp( ['最小值对应的（x,y）:' num2str(x(gbest_index,:))]);
figure
plot(record,'LineWidth',2);title('收敛过程')
xlabel('迭代次数');
ylabel('得到的目标函数值');
end

写代码不容易，要拿去用的小伙伴来个三连吧！