大家好，又见面了，我是你们的朋友全栈君。
我们在之前的博客中，对粒子群算法的实现进行了讲解，主要讲解了粒子群算法的产生，实现步骤，并且通过代码来实现了当输入变量为n维向量时的粒子群算法。
许多网友对之前的代码有些疑惑，并且提到了几个问题：
1，对于之前的代码，觉得有些看不懂？
2，如果输入变量不再是n维向量，而是n*n维的矩阵，甚至n*n*n…维的矩阵，该如何编写代码？
3，如果目标函数存在约束条件，该如何编写代码？
本节就基于之前博客的内容，出现的问题，以及这两天我的一些学习，来对上述问题进行一个解答。
上一讲博客链接在此:(https://blog.csdn.net/DBLLLLLLLL/article/details/82965722)
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:dim
pop_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:sizepop
if fitness_gbest(j) < fitness_zbest % 如果求最小值，则为<; 如果求最大值，则为>;
zbest = pop_x(:,j);
fitness_zbest=fitness_gbest(j);
end
end
plot3(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:sizepop
if pop_v(i,j)>vlimit(i,2)
pop_v(i,j)=vlimit(i,2);
end
if pop_v(i,j) < vlimit(i,1)
pop_v(i,j)=vlimit(i,1);
end
end
end
% 更新位置并对位置进行边界处理
pop_x = pop_x + pop_v;% 位置更新
% 边界位置处理
for i=1:dim
for j=1:sizepop
if pop_x(i,j)>xlimit(i,2)
pop_x(i,j)=xlimit(i,2);
end
if pop_x(i,j) < xlimit(i,1)
pop_x(i,j)=xlimit(i,1);
end
end
end
% 自适应变异
for j=1:sizepop
if rand > 0.85
i=ceil(dim*rand);
pop_x(i,j)=xlimit(i, 1) + (xlimit(i, 2) - xlimit(i, 1)) * rand;
end
end
% 计算新种群各个个体位置的适应度
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); % 更新群体历史最佳适应度
end
end
record(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;
end
iter = 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:dim
for j=1:sizepop
pop_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:sizepop
if fitness_gbest(j) < fitness_zbest % 如果求最小值，则为<; 如果求最大值，则为>;
zbest = pop_x(:,j);
fitness_zbest=fitness_gbest(j);
end
end
%% 粒子群迭代
% 更新速度并对速度进行边界处理
% 更新位置并对位置进行边界处理
% 进行自适应变异
% 计算新种群各个个体位置的适应度
% 新适应度与个体历史最佳适应度做比较
% 个体历史最佳适应度与种群历史最佳适应度做比较
% 再次循环或结束
iter = 1; %迭代次数
record = zeros(ger, 1); % 记录器
while iter <= ger
for 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:dim
if pop_v(i,j) > vlimit_max(i)
pop_v(i,j) = vlimit_max(i);
end
if pop_v(i,j) < vlimit_min(i)
pop_v(i,j) = vlimit_min(i);
end
end
% 更新位置并对位置进行边界处理
pop_x(:,j) = pop_x(:,j) + pop_v(:,j);% 位置更新
for i=1:dim
if pop_x(i,j) > xlimit_max(i)
pop_x(i,j) = xlimit_max(i);
end
if pop_x(i,j) < xlimit_min(i)
pop_x(i,j) = xlimit_min(i);
end
end
% 进行自适应变异
if rand > 0.85
i=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
end
record(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:dim
for j=1:sizepop
pop_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:sizepop
if cons1(pop_x(:,j))
if cons2(pop_x(:,j))
fitness_gbest(j) = fun(pop_x(:,j)); % 每个个体的历史最佳适应度
else
fitness_gbest(j) = 10^10;
end
else
fitness_gbest(j) = 10^10;
end
end
zbest = pop_x(:,1); % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1); % 种群的历史最佳适应度
for j=1:sizepop
if fitness_gbest(j) < fitness_zbest % 如果求最小值，则为<; 如果求最大值，则为>;
zbest = pop_x(:,j);
fitness_zbest=fitness_gbest(j);
end
end
%% 粒子群迭代
% 更新速度并对速度进行边界处理
% 更新位置并对位置进行边界处理
% 进行自适应变异
% 进行约束条件判断并计算新种群各个个体位置的适应度
% 新适应度与个体历史最佳适应度做比较
% 个体历史最佳适应度与种群历史最佳适应度做比较
% 再次循环或结束
iter = 1; %迭代次数
record = zeros(ger, 1); % 记录器
while iter <= ger
for 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:dim
if pop_v(i,j) > vlimit_max(i)
pop_v(i,j) = vlimit_max(i);
end
if pop_v(i,j) < vlimit_min(i)
pop_v(i,j) = vlimit_min(i);
end
end
% 更新位置并对位置进行边界处理
pop_x(:,j) = pop_x(:,j) + pop_v(:,j);% 位置更新
for i=1:dim
if pop_x(i,j) > xlimit_max(i)
pop_x(i,j) = xlimit_max(i);
end
if pop_x(i,j) < xlimit_min(i)
pop_x(i,j) = xlimit_min(i);
end
end
% 进行自适应变异
if rand > 0.85
i=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)); % 当前个体的适应度
else
fitness_pop(j) = 10^10;
end
else
fitness_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
end
record(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:dim1
for i_2=1:dim2
for j=1:sizepop
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; % 初始种群的位置
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; % 初始种群的速度
end
end
end
gbest = pop_x; % 每个个体的历史最佳位置
for j=1:sizepop % 每个个体的历史最佳适应度
if cons1(pop_x(:,:,j)) % 约束条件
fitness_gbest(j)=fun(pop_x(:,:,j));
else
fitness_gbest(j) = 10^10;
end
end
zbest = pop_x(:,:,1); % 种群的历史最佳位置
fitness_zbest = fitness_gbest(1); % 种群的历史最佳适应度
for j=1:sizepop
if fitness_gbest(j) < fitness_zbest % 如果求最小值，则为<; 如果求最大值，则为>;
zbest = pop_x(:,:,j);
fitness_zbest=fitness_gbest(j);
end
end
%% 粒子群迭代
% 更新速度并对速度进行边界处理
% 更新位置并对位置进行边界处理
% 进行自适应变异
% 进行约束条件判断并计算新种群各个个体位置的适应度
% 新适应度与个体历史最佳适应度做比较
% 个体历史最佳适应度与种群历史最佳适应度做比较
% 再次循环或结束
iter = 1; %迭代次数
record = zeros(ger, 1); % 记录器
while iter <= ger
for 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:dim1
for i_2=1:dim2
if 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);
end
if 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
end
end
% 更新位置并对位置进行边界处理
pop_x(:,:,j) = pop_x(:,:,j) + pop_v(:,:,j);% 位置更新
for i_1=1:dim1
for i_2=1:dim2
if 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);
end
if 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);
end
end
end
% 进行自适应变异
if rand > 0.85
i_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)); % 当前个体的适应度
else
fitness_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
end
record(iter) = fitness_zbest;%最大值记录
iter = iter+1;
end
%% 迭代结果输出
plot(record);title('收敛过程')
disp(['最优值：',num2str(fitness_zbest)]);
disp('变量取值：');
disp(zbest);
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/133327.html原文链接：https://javaforall.net
