MATLAB实现粒子群算法的进阶讲解（多维+约束条件）

我们在之前的博客中，对粒子群算法的实现进行了讲解，主要讲解了粒子群算法的产生，实现步骤，并且通过代码来实现了当输入变量为n维向量时的粒子群算法。
许多网友对之前的代码有些疑惑，并且提到了几个问题：
1，对于之前的代码，觉得有些看不懂？
2，如果输入变量不再是n维向量，而是n*n维的矩阵，甚至n*n*n…维的矩阵，该如何编写代码？
3，如果目标函数存在约束条件，该如何编写代码？
本节就基于之前博客的内容，出现的问题，以及这两天我的一些学习，来对上述问题进行一个解答。

上一讲博客链接在此:(https://blog.csdn.net/DBLLLLLLLL/article/details/82965722)

文章目录

1.自变量为2维向量的粒子群算法
2.自变量为n维向量的粒子群算法
3.自变量为n维向量且存在约束条件的粒子群算法
4.自变量为n*n维矩阵且存在约束条件的粒子群算法

1.自变量为2维向量的粒子群算法

首先还是从最简单的开始，关于自变量为2维的粒子群算法，其代码的实现过程可以参考上一讲的博客，这里我对于该代码进行了一些改进。主要改进就是增加了一个自适应变异机制，即每次算法计算周期，都有一定的概率随机改变某一个个体的数值，从而提高了整个粒子群的寻优能力。关于这个能力，感兴趣的人可以尝试将代码中的自适应变异的代码段注释掉，观察结果，可以有更直观的感受。

%% 清空环境
clc;clear all;close all;
%% 目标函数
%  下面这个函数叫做Schaffer函数，最小值在（0,0）处，为0
fun= @(a,b)(0.5+(sin(sqrt(a.^2+b.^2)).^2-0.5)./((1+0.001.*(a.^2+b.^2)).^2));
%  绘图函数，可舍弃
figure(1);
[x0_1, x0_2]=meshgrid(-5:0.1:5);
y0=fun(x0_1,x0_2);
mesh(x0_1, x0_2, y0);
hold on;
%% 设置种群参数
%   需要自行配置
sizepop = 500;                         % 初始种群个数
dim = 2;                           % 空间维数
ger = 300;                       % 最大迭代次数
xlimit = [ -5,5 ; -5 ,5 ];        % 设置位置参数限制(矩阵的形式可以多维)
vlimit = [ -1.5,1.5 ; -1.5,1.5]; % 设置速度限制
c_1 = 0.8;                       % 惯性权重
c_2 = 0.5;                       % 自我学习因子
c_3 = 0.5;                       % 群体学习因子
%% 生成初始种群
%  首先随机生成初始种群位置
%  然后随机生成初始种群速度
%  然后初始化个体历史最佳位置，以及个体历史最佳适应度
%  然后初始化群体历史最佳位置，以及群体历史最佳适应度
%  绘图（可舍弃）for i = 1:dimpop_x(i,:) = xlimit(i, 1) + (xlimit(i, 2) - xlimit(i, 1)) * rand(1, sizepop);%初始种群的位置
end
pop_v = rand(dim, sizepop);                   % 初始种群的速度
gbest = pop_x;                                % 每个个体的历史最佳位置
fitness_gbest = fun(pop_x(1,:),pop_x(2,:));   % 每个个体的历史最佳适应度
zbest = pop_x(:,1);                          % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1);             % 种群的历史最佳适应度
for j=1:sizepopif fitness_gbest(j) < fitness_zbest       % 如果求最小值，则为<; 如果求最大值，则为>; zbest = pop_x(:,j);fitness_zbest=fitness_gbest(j);end
endplot3(gbest(1,:),gbest(2,:),fun(gbest(1,:),gbest(2,:)), 'ro');title('初始状态图');
hold on;
figure(2);
mesh(x0_1, x0_2, y0);
hold on;
plot3(gbest(1,:),gbest(2,:),fun(gbest(1,:),gbest(2,:)), 'ro');
hold on;%% 粒子群迭代
%    更新速度并对速度进行边界处理
%    更新位置并对位置进行边界处理
%    进行自适应变异
%    计算新种群各个个体位置的适应度
%    新适应度与个体历史最佳适应度做比较
%    个体历史最佳适应度与种群历史最佳适应度做比较
%    再次循环或结束iter = 1;                        %迭代次数
times = 1;                       %用于显示，可舍弃
record = zeros(ger, 1);          % 记录器
while iter <= ger%    更新速度并对速度进行边界处理 pop_v =  c_1 * pop_v  + c_2 * rand *(gbest - pop_x) + c_3 * rand *(repmat(zbest, 1, sizepop) - pop_x);% 速度更新for i=1:dim for j=1:sizepopif  pop_v(i,j)>vlimit(i,2)pop_v(i,j)=vlimit(i,2);endif  pop_v(i,j) < vlimit(i,1)pop_v(i,j)=vlimit(i,1);endendend %    更新位置并对位置进行边界处理pop_x = pop_x + pop_v;% 位置更新% 边界位置处理for i=1:dim for j=1:sizepopif  pop_x(i,j)>xlimit(i,2)pop_x(i,j)=xlimit(i,2);endif  pop_x(i,j) < xlimit(i,1)pop_x(i,j)=xlimit(i,1);endendend%    自适应变异for j=1:sizepopif rand > 0.85i=ceil(dim*rand);pop_x(i,j)=xlimit(i, 1) + (xlimit(i, 2) - xlimit(i, 1)) * rand;endend%    计算新种群各个个体位置的适应度fitness_pop = fun(pop_x(1,:),pop_x(2,:)) ; % 当前所有个体的适应度%    新适应度与个体历史最佳适应度做比较for j = 1:sizepop      if fitness_pop(j) < fitness_gbest(j)       % 如果求最小值，则为<; 如果求最大值，则为>; gbest(:,j) = pop_x(:,j);               % 更新个体历史最佳位置            fitness_gbest(j) = fitness_pop(j);     % 更新个体历史最佳适应度end end%    个体历史最佳适应度与种群历史最佳适应度做比较for j = 1:sizepop  if fitness_gbest(j) < fitness_zbest        % 如果求最小值，则为<; 如果求最大值，则为>; zbest = gbest(:,j);                    % 更新群体历史最佳位置  fitness_zbest=fitness_gbest(j);        % 更新群体历史最佳适应度  endendrecord(iter) = fitness_zbest;%最大值记录if times >= 10        %显示函数 可以舍去cla;mesh(x0_1, x0_2, y0);plot3(pop_x(1,:),pop_x(2,:),fun(pop_x(1,:),pop_x(2,:)), 'ro');title('状态位置变化');pause(0.5);times=0;enditer = iter+1;times=times+1;        %显示函数 可以舍去
end
%% 迭代结果输出figure(3);plot(record);title('收敛过程')
figure(4);
mesh(x0_1, x0_2, y0);
hold on;
plot3(pop_x(1,:),pop_x(2,:),fun(pop_x(1,:),pop_x(2,:)), 'ro');title('最终状态图');disp(['最优值：',num2str(fitness_zbest)]);
disp('变量取值：');
disp(zbest);

2.自变量为n维向量的粒子群算法

根据自变量为2位的粒子群算法的代码，我们可以很容易就类比出来自变量为n维的粒子群算法的代码。这里很多人对于这个类比过程不太擅长，那我把已经写好的5维的代码发出来供大家参考。

%% 清空环境
clc;clear all;close all;
%% 目标函数
%    下列函数为Rastrigin函数，是一个多峰值的函数，在0向量处取得全局最小值0
fun= @(X)(sum((X.^2-10*cos(2*pi*X)+10),1));%% 设置种群参数
%   需要自行配置
sizepop = 500;                    % 初始种群个数
dim = 5;                          % 空间维数
ger = 200;                       % 最大迭代次数
xlimit_max = 5.12*ones(dim,1);    % 设置位置参数限制(矩阵的形式可以多维)
xlimit_min = -5.12*ones(dim,1);
vlimit_max = 1*ones(dim,1);       % 设置速度限制
vlimit_min = -1*ones(dim,1);
c_1 = 0.8;                        % 惯性权重
c_2 = 0.5;                        % 自我学习因子
c_3 = 0.5;                        % 群体学习因子 %% 生成初始种群
%  首先随机生成初始种群位置
%  然后随机生成初始种群速度
%  然后初始化个体历史最佳位置，以及个体历史最佳适应度
%  然后初始化群体历史最佳位置，以及群体历史最佳适应度
for i=1:dimfor j=1:sizepoppop_x(i,j) = xlimit_min(i)+(xlimit_max(i) - xlimit_min(i))*rand;pop_v(i,j) = vlimit_min(i)+(vlimit_max(i) - vlimit_min(i))*rand;  % 初始种群的速度end
end
gbest = pop_x;                                % 每个个体的历史最佳位置
fitness_gbest = fun(pop_x);                   % 每个个体的历史最佳适应度
zbest = pop_x(:,1);                           % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1);             % 种群的历史最佳适应度
for j=1:sizepopif fitness_gbest(j) < fitness_zbest       % 如果求最小值，则为<; 如果求最大值，则为>; zbest = pop_x(:,j);fitness_zbest=fitness_gbest(j);end
end%% 粒子群迭代
%    更新速度并对速度进行边界处理
%    更新位置并对位置进行边界处理
%    进行自适应变异
%    计算新种群各个个体位置的适应度
%    新适应度与个体历史最佳适应度做比较
%    个体历史最佳适应度与种群历史最佳适应度做比较
%    再次循环或结束iter = 1;                        %迭代次数
record = zeros(ger, 1);          % 记录器
while iter <= gerfor j=1:sizepop%    更新速度并对速度进行边界处理 pop_v(:,j)= c_1 * pop_v(:,j) + c_2*rand*(gbest(:,j)-pop_x(:,j))+c_3*rand*(zbest-pop_x(:,j));% 速度更新for i=1:dimif  pop_v(i,j) > vlimit_max(i)pop_v(i,j) = vlimit_max(i);endif  pop_v(i,j) < vlimit_min(i)pop_v(i,j) = vlimit_min(i);endend%    更新位置并对位置进行边界处理pop_x(:,j) = pop_x(:,j) + pop_v(:,j);% 位置更新for i=1:dimif  pop_x(i,j) > xlimit_max(i)pop_x(i,j) = xlimit_max(i);endif  pop_x(i,j) < xlimit_min(i)pop_x(i,j) = xlimit_min(i);endend%    进行自适应变异if rand > 0.85i=ceil(dim*rand);pop_x(i,j)=xlimit_min(i) + (xlimit_max(i) - xlimit_min(i)) * rand;end%    计算新种群各个个体位置的适应度fitness_pop(j) = fun(pop_x(:,j));                      % 当前个体的适应度%    新适应度与个体历史最佳适应度做比较if fitness_pop(j) < fitness_gbest(j)       % 如果求最小值，则为<; 如果求最大值，则为>; gbest(:,j) = pop_x(:,j);               % 更新个体历史最佳位置            fitness_gbest(j) = fitness_pop(j);     % 更新个体历史最佳适应度end   %    个体历史最佳适应度与种群历史最佳适应度做比较if fitness_gbest(j) < fitness_zbest        % 如果求最小值，则为<; 如果求最大值，则为>; zbest = gbest(:,j);                    % 更新群体历史最佳位置  fitness_zbest=fitness_gbest(j);        % 更新群体历史最佳适应度  end    endrecord(iter) = fitness_zbest;%最大值记录iter = iter+1;end
%% 迭代结果输出plot(record);title('收敛过程')
disp(['最优值：',num2str(fitness_zbest)]);
disp('变量取值：');
disp(zbest);

3.自变量为n维向量且存在约束条件的粒子群算法

解决了输入为n维向量的粒子群算法的代码，接下来我们来看，当存在约束条件时，如何构造这个算法。这里采用的思想是，对每一个粒子进行判断，判断该粒子的取值是否满足约束条件，如果满足，则进行正常的适应度计算，如果不满足，则将该粒子的适应度直接赋值为一个很大的数（求最小值时赋值成大数，求最大值时应该赋值成小数）。具体代码实现如下：

%% 清空环境
clc;clear all;close all;
%% 目标函数
%    下列函数为：
%    y=0.072*x1+0.063*x2+0.057*x3+0.05*x4+0.032*x5+0.0442*x6+0.0675*x7+7
%    约束条件为：
%    0.072*x1+0.063*x2+0.057*x3+0.05*x4+0.032*x5+0.0442*x6+0.0675*x7<=264.4
%    128*x1+78.1*x2+64.1*x3+43*x4+58.1*x5+36.9*x6+50.5*x7<=67919
a=[0.072,0.063,0.057,0.05,0.032,0.0442,0.0675];
b=[0.072,0.063,0.057,0.005,0.032,0.0442,0.0675];
c=[128,78.1,64.1,43,58.1,36.9,50.5];
fun= @(X)(a*X+7);
cons1= @(X)(b*X<=264.4);
cons2= @(X)(c*X>=-20);
%% 设置种群参数
%   需要自行配置
sizepop = 500;                         % 初始种群个数
dim = 7;                           % 空间维数
ger = 500;                       % 最大迭代次数
xlimit_max = 50*ones(dim,1);     % 设置位置参数限制(矩阵的形式可以多维)
xlimit_min = -50*ones(dim,1);
vlimit_max = 1*ones(dim,1);      % 设置速度限制
vlimit_min = -1*ones(dim,1);
c_1 = 0.8;                       % 惯性权重
c_2 = 0.5;                       % 自我学习因子
c_3 = 0.5;                       % 群体学习因子 %% 生成初始种群
%  首先随机生成初始种群位置
%  然后随机生成初始种群速度
%  然后初始化个体历史最佳位置，以及个体历史最佳适应度
%  然后初始化群体历史最佳位置，以及群体历史最佳适应度
for i=1:dimfor j=1:sizepoppop_x(i,j) = xlimit_min(i)+(xlimit_max(i) - xlimit_min(i))*rand;  % 初始种群的位置pop_v(i,j) = vlimit_min(i)+(vlimit_max(i) - vlimit_min(i))*rand;  % 初始种群的速度end
end
gbest = pop_x;                                % 每个个体的历史最佳位置
for j=1:sizepopif cons1(pop_x(:,j))if cons2(pop_x(:,j))fitness_gbest(j) = fun(pop_x(:,j));                      % 每个个体的历史最佳适应度elsefitness_gbest(j) = 10^10; endelsefitness_gbest(j) = 10^10; end
end
zbest = pop_x(:,1);                           % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1);             % 种群的历史最佳适应度
for j=1:sizepopif fitness_gbest(j) < fitness_zbest       % 如果求最小值，则为<; 如果求最大值，则为>; zbest = pop_x(:,j);fitness_zbest=fitness_gbest(j);end
end%% 粒子群迭代
%    更新速度并对速度进行边界处理
%    更新位置并对位置进行边界处理
%    进行自适应变异
%    进行约束条件判断并计算新种群各个个体位置的适应度
%    新适应度与个体历史最佳适应度做比较
%    个体历史最佳适应度与种群历史最佳适应度做比较
%    再次循环或结束iter = 1;                        %迭代次数
record = zeros(ger, 1);          % 记录器
while iter <= gerfor j=1:sizepop%    更新速度并对速度进行边界处理 pop_v(:,j)= c_1 * pop_v(:,j) + c_2*rand*(gbest(:,j)-pop_x(:,j))+c_3*rand*(zbest-pop_x(:,j));% 速度更新for i=1:dimif  pop_v(i,j) > vlimit_max(i)pop_v(i,j) = vlimit_max(i);endif  pop_v(i,j) < vlimit_min(i)pop_v(i,j) = vlimit_min(i);endend%    更新位置并对位置进行边界处理pop_x(:,j) = pop_x(:,j) + pop_v(:,j);% 位置更新for i=1:dimif  pop_x(i,j) > xlimit_max(i)pop_x(i,j) = xlimit_max(i);endif  pop_x(i,j) < xlimit_min(i)pop_x(i,j) = xlimit_min(i);endend%    进行自适应变异if rand > 0.85i=ceil(dim*rand);pop_x(i,j)=xlimit_min(i) + (xlimit_max(i) - xlimit_min(i)) * rand;end%    进行约束条件判断并计算新种群各个个体位置的适应度if cons1(pop_x(:,j))if cons2(pop_x(:,j))fitness_pop(j) = fun(pop_x(:,j));                      % 当前个体的适应度elsefitness_pop(j) = 10^10; endelsefitness_pop(j) = 10^10; end%    新适应度与个体历史最佳适应度做比较if fitness_pop(j) < fitness_gbest(j)       % 如果求最小值，则为<; 如果求最大值，则为>; gbest(:,j) = pop_x(:,j);               % 更新个体历史最佳位置            fitness_gbest(j) = fitness_pop(j);     % 更新个体历史最佳适应度end   %    个体历史最佳适应度与种群历史最佳适应度做比较if fitness_gbest(j) < fitness_zbest        % 如果求最小值，则为<; 如果求最大值，则为>; zbest = gbest(:,j);                    % 更新群体历史最佳位置  fitness_zbest=fitness_gbest(j);        % 更新群体历史最佳适应度  end    endrecord(iter) = fitness_zbest;%最大值记录iter = iter+1;end
%% 迭代结果输出plot(record);title('收敛过程')
disp(['最优值：',num2str(fitness_zbest)]);
disp('变量取值：');
disp(zbest);

4.自变量为n*n维矩阵且存在约束条件的粒子群算法

有些网友说，自己所要解决的情况，输入变量不仅仅是一个n维的向量，而是一个n*n维的矩阵，甚至n*n*n*…维的空间矩阵。想要采用粒子群算法，该如何实现呢？
其实每一个n*n维矩阵，可以看成是n维向量，每个向量的元素又是一个向量，这样不断嵌套得到的。因此总的来说，整个代码的框架并没有太大的变化，仅仅是将原来的计算多嵌套几层就好。
代码如下：

%% 清空环境
clc;clear all;close all;
%% 目标函数
%    该函数的变量是一个3*4的矩阵形式的变量
%    cons1=@(X)(1);表示无约束条件
%    cons1=@(X)(sum(sum(X))>=0);可以尝试添加约束条件
a=[2.5,4.5,-2];
b=[2;4.2;4;-4];
c=[1,2,3,4;-4,-3,-2,-1;1,3,2,4];
fun= @(X)(sum(sum(c.*X.^2))+a*X*b);
cons1=@(X)(sum(sum(X))>=0);
%% 设置种群参数
%   需要自行配置
%   下文中所有出现dim2，或者：的地方，都需要根据实际的矩阵维度来进行更改
sizepop = 500;                         % 初始种群个数
dim1 = 3;                           % 空间维数
dim2 = 4;
ger = 500;                       % 最大迭代次数
xlimit_max = 50*ones(dim1,dim2,1);     % 设置位置参数限制(矩阵的形式可以多维)
xlimit_min = -50*ones(dim1,dim2,1);
vlimit_max = 1*ones(dim1,dim2,1);      % 设置速度限制
vlimit_min = -1*ones(dim1,dim2,1);
c_1 = 0.8;                       % 惯性权重
c_2 = 0.5;                       % 自我学习因子
c_3 = 0.5;                       % 群体学习因子 %% 生成初始种群
%  首先随机生成初始种群位置
%  然后随机生成初始种群速度
%  然后初始化个体历史最佳位置，以及个体历史最佳适应度
%  然后初始化群体历史最佳位置，以及群体历史最佳适应度
for i_1=1:dim1for i_2=1:dim2      for j=1:sizepoppop_x(i_1,i_2,j) = xlimit_min(i_1,i_2)+(xlimit_max(i_1,i_2) - xlimit_min(i_1,i_2))*rand;  % 初始种群的位置pop_v(i_1,i_2,j) = vlimit_min(i_1,i_2)+(vlimit_max(i_1,i_2) - vlimit_min(i_1,i_2))*rand;  % 初始种群的速度endend
end
gbest = pop_x;                                % 每个个体的历史最佳位置
for j=1:sizepop                               % 每个个体的历史最佳适应度if cons1(pop_x(:,:,j))                      % 约束条件       fitness_gbest(j)=fun(pop_x(:,:,j));                    elsefitness_gbest(j) = 10^10; end
end
zbest = pop_x(:,:,1);                           % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1);             % 种群的历史最佳适应度
for j=1:sizepopif fitness_gbest(j) < fitness_zbest       % 如果求最小值，则为<; 如果求最大值，则为>; zbest = pop_x(:,:,j);fitness_zbest=fitness_gbest(j);end
end%% 粒子群迭代
%    更新速度并对速度进行边界处理
%    更新位置并对位置进行边界处理
%    进行自适应变异
%    进行约束条件判断并计算新种群各个个体位置的适应度
%    新适应度与个体历史最佳适应度做比较
%    个体历史最佳适应度与种群历史最佳适应度做比较
%    再次循环或结束iter = 1;                        %迭代次数
record = zeros(ger, 1);          % 记录器
while iter <= gerfor j=1:sizepop%    更新速度并对速度进行边界处理 pop_v(:,:,j)= c_1 * pop_v(:,:,j) + c_2*rand*(gbest(:,:,j)-pop_x(:,:,j))+c_3*rand*(zbest-pop_x(:,:,j));% 速度更新for i_1=1:dim1for i_2=1:dim2if  pop_v(i_1,i_2,j) > vlimit_max(i_1,i_2)pop_v(i_1,i_2,j) = vlimit_max(i_1,i_2);endif  pop_v(i_1,i_2,j) < vlimit_min(i_1,i_2)pop_v(i_1,i_2,j) = vlimit_min(i_1,i_2);end                endend%    更新位置并对位置进行边界处理pop_x(:,:,j) = pop_x(:,:,j) + pop_v(:,:,j);% 位置更新for i_1=1:dim1for i_2=1:dim2if  pop_x(i_1,i_2,j) > xlimit_max(i_1,i_2)pop_x(i_1,i_2,j) = xlimit_max(i_1,i_2);endif  pop_x(i_1,i_2,j) < xlimit_min(i_1,i_2)pop_x(i_1,i_2,j) = xlimit_min(i_1,i_2);endendend%    进行自适应变异if rand > 0.85i_1=ceil(dim1*rand);i_2=ceil(dim2*rand);pop_x(i_1,i_2,j)=xlimit_min(i_1,i_2) + (xlimit_max(i_1,i_2) - xlimit_min(i_1,i_2)) * rand;end%    进行约束条件判断并计算新种群各个个体位置的适应度if cons1(pop_x(:,:,j))                                       % 约束条件fitness_pop(j) = fun(pop_x(:,:,j));                      % 当前个体的适应度elsefitness_pop(j) = 10^10; end%    新适应度与个体历史最佳适应度做比较if fitness_pop(j) < fitness_gbest(j)       % 如果求最小值，则为<; 如果求最大值，则为>; gbest(:,:,j) = pop_x(:,:,j);               % 更新个体历史最佳位置            fitness_gbest(j) = fitness_pop(j);     % 更新个体历史最佳适应度end   %    个体历史最佳适应度与种群历史最佳适应度做比较if fitness_gbest(j) < fitness_zbest        % 如果求最小值，则为<; 如果求最大值，则为>; zbest = gbest(:,:,j);                    % 更新群体历史最佳位置  fitness_zbest=fitness_gbest(j);        % 更新群体历史最佳适应度  end    endrecord(iter) = fitness_zbest;%最大值记录iter = iter+1;end
%% 迭代结果输出plot(record);title('收敛过程')
disp(['最优值：',num2str(fitness_zbest)]);
disp('变量取值：');
disp(zbest);