PSO算法代码详解
前期文章介绍了PSO算法的基本原理，来帮助同学们去了解PSO算法的流程。因为篇幅的问题呢，就没有在前期文章附上代码详解，为了弥补遗憾，我把代码解释作为独立的一篇博客发表，方便不同需求的同学们来交流学习。 希望大家可以踊跃评论，一起交流学习，也可以给我指出改进的点和错误的点。
前期文章链接：进化算法之粒子群算法介绍附代码——PSO
一元函数的优化
fun函数，要重新建立文件
fun函数，要重新建立文件
function y = fun(x)
%函数用于计算粒子适应度值
%x input 输入粒子
%y output 粒子适应度值
y = sin(10*pi*x) ./ x;
主文件代码
%% I. 清空环境
clc
clear all //基本的，必有的操作。背下来
%% II. 绘制目标函数曲线图
x = 1:0.01:2;
/* 因为是函数，y值在是x的范围内求解最大值的，这里就是在设置x的范围，x的范围是1到2，Matlab中范围表示为a:b。
大家都知道极限这个概念，所以说，x从1到2不能是从1到无穷，在从无穷到2吧，所以我们要设置x从1到2的步伐，即每一步为0.01*/
y = sin(10*pi*x) ./ x; //Matlab中除法的特色 ./ ①
/* 开始绘制函数图像了吧 ②*/
figure //figure表示新建一个图窗口，以免后续的绘图语句覆盖原图
plot(x, y) //因为是一元函数，所以plot点设置，x和y两个轴就可以了，默认是视线
hold on //hold on表示在原图的基础上绘制新的图像
那么这里我们就可以看到在Matlab中生成的函数图像，即
%% III. 参数初始化
//这里大家会有疑问，为什么参数里面没有重要惯性权重W呢，其实这种操作，程序会自动把w的值赋为1
c1 = 1.49445;
c2 = 1.49445; //c1和c2是学习因子，他的值为多少，是由我们来设定的
maxgen = 50; % 进化次数 // 也就是迭代次数
sizepop = 10; %种群规模 // 就是说由几个粒子组成的一个粒子群，这里就是10个粒子
Vmax = 0.5; %速度的范围，超过则用边界值。
Vmin = -0.5; // 同样的速度的最大值Vmax和速度最小值Vmin都是由我们来设定的
popmax = 2; %个体的变化范围
popmin = 1; // 因为在绘制目标函数曲线图部分，我们是设置了x的范围，所以要严格按照该部分设定值，不能随心所欲了。
%% IV. 产生初始粒子和速度
for i = 1:sizepop //进入for循环；循环是从1 到 sizepop的值（10），也就是说通过for循环，产生了10个粒子。
% 随机产生一个种群
pop(i,:) = (rands(1) + 1) / 2 + 1; %初始种群，rands产生(-1,1)，调整到(1,2)
/*这里的初始化，就是要如何实现在x的范围，即1-2，rands(1)表示为生成1行1列均匀分布在（-1~1）之间的伪随机数，
自定义的公式通过计算的，pop的值实现了x的范围。pop也是自定义的，表示的粒子的位置。重要的是(i,:)的含义 ③
*/
V(i,:) = 0.5 * rands(1); %初始化速度 //同理
% 计算适应度
fitness(i) = fun(pop(i,:));
/* 通过fun()函数来寻找，来计算适应度的值。
这里的值pop的10行1列，每一行的值通过fun函数计算出适应度的值，产生了一个1行n列的值*/
end
这里我们调试一个代码，就会发现，已经生成了10行1列的矩阵，那么每一行的值就可以代表每一个粒子在一维空间内的pop值，和v值。
- pop值
- V值
- fitness值
③ Matlab中rand，randn，rands和randi函数使用
这里要注意的是，此案例中rands(1)括号里面仅有一个值，是因为我们是在一维空间里面求解，一维空间在矩阵中的表示就是一个值单独占到1行。n个值占n行。如果在n维的空间里面，就要表示为rands(1,n)，即表示为1行n列，即一个粒子，在矩阵中，在一行中产生了20个位置。
X(:,i)与X(i,:)的区别
%% V. 个体极值和群体极值
[bestfitness bestindex] = max(fitness);
/*这里就是寻找fitness的值的最大值，赋值给bestfitness，bestindex就是说在bestfitness在fitness1行10列中的位置。*/
zbest = pop(bestindex,:); %全局最佳
/*通过函数，对比每一个值，可以找到，全局最佳的值，
根据bestindex在fitness的位置可以找到粒子在pop中的位置，
比如bestindex=6，那么pop中第6行的值，就是全局最佳值。*/
gbest = pop; %个体最佳
/*pop的值赋值给gbest。即个人最佳值，因为刚开始每一个粒子都是初始化的，
那么10个粒子的最佳位置，就是初始化pop值。*/
fitnessgbest = fitness; %个体最佳适应度值 //赋值
fitnesszbest = bestfitness; %全局最佳适应度值 //赋值
-
bestfitness和bestindex
-
zbest
-
gbest
%% VI. 迭代寻优
for i = 1:maxgen //需要嵌套循环，外部循环——迭代需要10次，从第1次到第10次的循环
for j = 1:sizepop //内部循环，每一个粒子都需要更新一次
% 速度更新
V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:)); //速度更新公式 w为1
V(j,find(V(j,:)>Vmax)) = Vmax;
V(j,find(V(j,:)<Vmin)) = Vmin;
/*更新后的速度，不能超出定义的速度范围，如果超出范围，
就把定义速度的最大（小）值赋值给更新后速度*/
% 种群更新
pop(j,:) = pop(j,:) + V(j,:); //位置更新公式
pop(j,find(pop(j,:)>popmax)) = popmax;
pop(j,find(pop(j,:)<popmin)) = popmin;
/*更新后的位置，不能超出定义的位置范围，如果超出范围，
就把定义位置的最大（小）值赋值给更新后位置*/
% 适应度值更新
fitness(j) = fun(pop(j,:)); //更新每一个粒子的适应度的值。
end
for j = 1:sizepop //内部循环
% 个体最优更新
if fitness(j) > fitnessgbest(j) //如果更新后的fitness大于初始化的个体最优值
gbest(j,:) = pop(j,:); //符合条件。更新gbest的值
fitnessgbest(j) = fitness(j); //将fitnessgbest 个体最优值更新
end
% 群体最优更新
if fitness(j) > fitnesszbest //如果更新后的fitness大于初始化的群体最优值
zbest = pop(j,:); //符合条件。更新zbest的值
fitnesszbest = fitness(j); //将fitnesszbest 群体最优值更新
end
end
yy(i) = fitnesszbest; //将群体最优值赋值给一个点
end
%% VII. 输出结果并绘图
[fitnesszbest zbest] //输出群体最佳适应度值(函数值)和 群体最佳值(粒子的值)
plot(zbest, fitnesszbest,'r*')
/*规定点的位置，用红色*号表示，*/
figure //重新绘图
plot(yy) //表示yy点
title('最优个体适应度','fontsize',12); //标题 字体大小
xlabel('进化代数','fontsize',12);//横坐标 字体大小
ylabel('适应度','fontsize',12);//纵坐标 字体大小
-
输出值
-
点在函数中的位置
-
点的具体位置
二元函数的最大值
代码一已经对PSO算法应该如何书写，做了详细介绍了，大致就是那个样子的，基本上就是套用模板。对一些参数进行修改。所以我这里附上代码，就不做详细介绍了。
fun函数，要重新建立文件
function y = fun(x)
%函数用于计算粒子适应度值
%x input 输入粒子
%y output 粒子适应度值
y = x(1).^2 + x(2).^2 - 10*cos(2*pi*x(1)) - 10*cos(2*pi*x(2)) + 20;
主文件代码
%% 题目2： 利用粒子群算法计算二元函数的最大值
%% I. 清空环境
clc
clear
%% II. 绘制目标函数曲线
figure
[x, y] = meshgrid(-5: 0.1: 5, -5: 0.1: 5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
mesh(x,y,z)
hold on
%% III. 参数初始化
c1 = 1.49445;
c2 = 1.49445;
maxgen = 1000; % 进化次数
sizepop = 100; % 种群规模
dimension = 2; % 这里因为是二元函数的求解，即二维，故列数为2
% 速度的边界
Vmax = 1;
Vmin = -1;
% 种群的边界
popmax = 5;
popmin = -5;
% 用于计算惯性权重，经验值
ws = 0.9;
we = 0.4;
% 给矩阵预分配内存
pop = zeros(sizepop, dimension);
V = zeros(sizepop, dimension);
fitness = zeros(sizepop, 1);
yy = zeros(maxgen);
w = zeros(maxgen);
%% IV. 产生初始粒子和速度
for i = 1: sizepop
% 随机产生一个种群
pop(i, :) = 5 * rands(1, 2); % 初始种群
V(i, :) = rands(1, 2); % 初始化速度
% 计算适应度
fitness(i) = fun(pop(i, :));
end
%% V. 初始化Personal best和Global best
[bestfitness, bestindex] = max(fitness);
gbest = pop(bestindex, :); % Global best
pbest = pop; % 个体最佳
fitnesspbest = fitness; % 个体最佳适应度值
fitnessgbest = bestfitness; % 全局最佳适应度值
%% VI. 迭代寻优
for i = 1: maxgen
w(i) = ws - (ws - we) * (i / maxgen); % 让惯性权重随着迭代次数而动态改变，控制搜索精度
for j = 1: sizepop
% 速度更新
V(j, :) = w(i)*V(j, :) + c1*rand*(pbest(j, :) - pop(j, :)) + c2*rand*(gbest - pop(j, :));
for k = 1: dimension
if V(j, k) > Vmax
V(j, k) = Vmax;
end
if V(j, k) < Vmin
V(j, k) = Vmin;
end
end
% 种群更新(位置更新)
pop(j, :) = pop(j, :) + V(j, :);
for k = 1: dimension
if pop(j, k) > popmax
pop(j, k) = popmax;
end
if pop(j, k) < popmin
pop(j, k) = popmin;
end
end
% 适应度值更新
fitness(j) = fun(pop(j, :));
end
for j = 1: sizepop
% 个体最优更新
if fitness(j) > fitnesspbest(j)
pbest(j, :) = pop(j, :);
fitnesspbest(j) = fitness(j);
end
%群体最优更新
if fitness(j) > fitnessgbest
gbest = pop(j, :);
fitnessgbest = fitness(j);
end
end
yy(i) = fitnessgbest; % 记录每次迭代完毕的群体最优解
end
%% VII.输出结果
[fitnessgbest, gbest]
plot3(gbest(1), gbest(2), fitnessgbest, 'bo','linewidth', 1.5)
figure
plot(yy)
title('最优个体适应度', 'fontsize', 12);
xlabel('进化代数', 'fontsize', 12);
ylabel('适应度', 'fontsize', 12);
这里给大家推荐一篇文章，可以参考了学习。
3228
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
