想必大家都知道优化算法,这种智能优化算法咱们就不多说,现在我们不讲原理,直接展示怎么实现怎么用的一个过程,最终目的在与我们会改和调用这个算法。
clc;clear;close all;
f= @(x) - (x - 10) .^ 2 + x .* sin(x) .* cos(2 * x) - x .* sin(2 * x) ;
首先清屏,防止出现其他问题,这是一个匿名函数,这是就是我们要求就目标函数,现在,你一眼看过去,就会能看出来,他是带sinX和cosX,毋庸置疑他是一个非线性函数,一般智能算法都是解决非线性函数。
figure(1);
fplot(f, [0 20], 'b-'); % 画出初始图像
d = 1; % 空间维数(该例子是1维)
N = 15; % 初始种群个数
A_limit = [0, 20]; % 设置位置限制
Z_limit = [-1, 1]; % 设置速度限制
因为是非线性函数,所以他是不规则的,所有我们必须要在一定范围找最优解 ,这这里我们限定的大小是0到20内移动。
x = A_limit(1) + (A_limit(2) - A_limit(1)) * rand(N, d); %初始每个粒子的位置
v = rand(N, d); % 初始每个粒子的速度
pbest = x; % 初始化每个个体的历史最佳位置
gbest = zeros(1, d); % 初始化种群的历史最佳位置
best = zeros(N, 1); % 初始化每个个体的历史最佳适应度 为 0
great = -inf; % 初始化种群历史最佳适应度 为 负无穷
iter = 50; % 最大迭代次数
w = 0.8; % 惯性权重
c1 = 0.5; % 自我学习因子
c2 = 0.2; % 群体学习因子
这里是一些初始设定,这些基本都是通用的,如果不了解,你直接用,不需要改。
hold on;
plot(x, f(x), 'ro');
title('初始状态图');
rd = zeros(iter, 1); % 记录器(用于记录 fg_best 的变化过程)
figure(2);
i = 1;
while i <= iter
fx = f(x) ; % 计算个体当前适应度
for j = 1:N
if best(j) < fx(j)
best(j) = fx(j); % 更新个体历史最佳适应度
pbest(j) = x(j); % 更新个体历史最佳位置
end
end
if great < max(best)
[great,ind_max] = max(best); % 更新群体历史最佳适应度
gbest = pbest(ind_max); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
% 注:将上面的两个式子换成下面两个是不可以的
% [gbest, ind_max] = max(pbest); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
% fg_best = fp_best(ind_max); % 更新群体历史最佳适应度
end
v = v * w + c1 * rand() * (pbest - x) + c2 * rand() * (repmat(gbest, N, 1) - x); % 速度更新
% 注: repmat(A,r1,r2):可将矩阵 扩充 为每个单位为A的r1*r2的矩阵
% 边界速度处理
v(v > Z_limit(2)) = Z_limit(2);
v(v < Z_limit(1)) = Z_limit(1);
x = x + v;% 位置更新
% 边界位置处理
x(x > A_limit(2)) = A_limit(2);
x(x < A_limit(1)) = A_limit(1);
rd(i) = great;%最大值记录
% 画动态展示图
zuo_x = 0 : 0.01 : 20;
plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
title('状态位置变化')
pause(0.1)
i = i + 1;
% if mod(i,10) == 0 % 显示进度
% i
% end
end
figure(3);
plot(rd);
title('收敛过程');
zuo_x = 0 : 0.01 : 20;
figure(4);
plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
title('最终状态图')
disp(['最佳适应度:',num2str(great)]);
disp(['最佳粒子的位置x:',num2str(gbest)]);
以上是为了展示动画过程,把每个值都记录下来,还有就是寻找每次的最佳值的过程,用循环做。
最后输出结果图。
结果图1:初始位置画出来的图
结果2:寻找过程画出来的图。
结果三:最后结果出来的图
结果四:整个寻找过程的图
这是最后找到的结果:最优解为:11.2572。
下面是全部的代码:
clc;clear;close all;
f= @(x) - (x - 10) .^ 2 + x.* cos(4 * x) - x .* sin(2 * x) ; % 适应度函数表达式(求这个函数的最大值)
figure(1);
fplot(f, [0 20], 'b-'); % 画出初始图像
d = 1; % 空间维数(该例子是1维)
N = 15; % 初始种群个数
A_limit = [0, 20]; % 设置位置限制
Z_limit = [-1, 1]; % 设置速度限制
x = A_limit(1) + (A_limit(2) - A_limit(1)) * rand(N, d); %初始每个粒子的位置
v = rand(N, d); % 初始每个粒子的速度
pbest = x; % 初始化每个个体的历史最佳位置
gbest = zeros(1, d); % 初始化种群的历史最佳位置
best = zeros(N, 1); % 初始化每个个体的历史最佳适应度 为 0
great = -inf; % 初始化种群历史最佳适应度 为 负无穷
iter = 50; % 最大迭代次数
w = 0.8; % 惯性权重
c1 = 0.5; % 自我学习因子
c2 = 0.2; % 群体学习因子
hold on;
plot(x, f(x), 'ro');
title('初始状态图');
rd = zeros(iter, 1); % 记录器(用于记录 fg_best 的变化过程)
figure(2);
i = 1;
while i <= iter
fx = f(x) ; % 计算个体当前适应度
for j = 1:N
if best(j) < fx(j)
best(j) = fx(j); % 更新个体历史最佳适应度
pbest(j) = x(j); % 更新个体历史最佳位置
end
end
if great < max(best)
[great,ind_max] = max(best); % 更新群体历史最佳适应度
gbest = pbest(ind_max); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
% 注:将上面的两个式子换成下面两个是不可以的
% [gbest, ind_max] = max(pbest); % 更新群体历史最佳位置,其中 ind_max 是最大值所在的下标
% fg_best = fp_best(ind_max); % 更新群体历史最佳适应度
end
v = v * w + c1 * rand() * (pbest - x) + c2 * rand() * (repmat(gbest, N, 1) - x); % 速度更新
% 注: repmat(A,r1,r2):可将矩阵 扩充 为每个单位为A的r1*r2的矩阵
% 边界速度处理
v(v > Z_limit(2)) = Z_limit(2);
v(v < Z_limit(1)) = Z_limit(1);
x = x + v;% 位置更新
% 边界位置处理
x(x > A_limit(2)) = A_limit(2);
x(x < A_limit(1)) = A_limit(1);
rd(i) = great;%最大值记录
% 画动态展示图
zuo_x = 0 : 0.01 : 20;
plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
title('状态位置变化')
pause(0.1)
i = i + 1;
% if mod(i,10) == 0 % 显示进度
% i
% end
end
figure(3);
plot(rd);
title('收敛过程');
zuo_x = 0 : 0.01 : 20;
figure(4);
plot(zuo_x, f(zuo_x), 'b-', x, f(x), 'ro');
title('最终状态图')
disp(['最佳适应度:',num2str(great)]);
disp(['最佳粒子的位置x:',num2str(gbest)]);
不多,唯手熟尔!