简单的粒子群优化算法的代码讲解及展示

       想必大家都知道优化算法，这种智能优化算法咱们就不多说，现在我们不讲原理，直接展示怎么实现怎么用的一个过程，最终目的在与我们会改和调用这个算法。

  

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)]);  

 不多，唯手熟尔！