15-2利用粒子群算法实现二元函数最大值的优化（附matlab程序）

1.简述

      

算法的基本原理

（1）PSO的基础：信息的社会共享

（2）粒子的两个属性：速度和位置（算法的两个核心要素）

       速度表示粒子下一步迭代时移动的方向和距离，位置是所求解问题的一个解。

粒子群算法的流程图

粒子群解释

① 第一项：惯性部分

由惯性权重和粒子自身速度构成，表示粒子对先前自身运动状态的信任。

② 第二项：认知部分

表示粒子本身的思考，即粒子自己经验的部分，可理解为粒子当前位置与自身历史最优位置之间的距离和方向。

③ 第三项：社会部分

表示粒子之间的信息共享与合作，即来源于群体中其他优秀粒子的经验，可理解为粒子当前位置与群体历史最优位置之间的距离和方向。

算法参数的详细解释

（1）粒子群规模：

一个正整数，推荐取值范围：[20,1000]，简单问题一般取20~40，较难或特定类别的问题可以取100~200。较小的种群规模容易陷入局部最优；较大的种群规模可以提高收敛性，更快找到全局最优解，但是相应地每次迭代的计算量也会增大；当种群规模增大至一定水平时，再增大将不再有显著的作用。

（2）粒子维度： 

粒子搜索的空间维数即为自变量的个数。

（3）迭代次数：

推荐取值范围：[50,100]，典型取值：60、70、100；这需要在优化的过程中根据实际情况进行调整，迭代次数太小的话解不稳定，太大的话非常耗时，没有必要。

（4）惯性权重：

1998年，Yuhui Shi和Russell Eberhart对基本粒子群算法引入了惯性权重(inertia weight)�，并提出动态调整惯性权重以平衡收敛的全局性和收敛速度，该算法被称为标准PSO算法（本文所介绍）[2]。

• 参数意义

惯性权重  表示上一代粒子的速度对当代粒子的速度的影响，或者说粒子对当前自身运动状态的信任程度，粒子依据自身的速度进行惯性运动。惯性权重使粒子保持运动的惯性和搜索扩展空间的趋势。值越大，探索新区域的能力越强，全局寻优能力越强，但是局部寻优能力越弱。反之，全局寻优能力越弱，局部寻优能力强。较大的  值有利于全局搜索，跳出局部极值，不至于陷入局部最优；而较小的 有利于局部搜索，让算法快速收敛到最优解。当问题空间较大时，为了在搜索速度和搜索精度之间达到平衡，通常做法是使算法在前期有较高的全局搜索能力以得到合适的种子，而在后期有较高的局部搜索能力以提高收敛精度，所以 不宜为一个固定的常数[3]。

当 =1 时，退化成基本粒子群算法，当 =0 时，失去对粒子本身经验的思考。推荐取值范围：[0.4,2]，典型取值：0.9、1.2、1.5、1.8

• 改善惯性权重 

在解决实际优化问题时，往往希望先采用全局搜索，使搜索空间快速收敛于某一区域，然后采用局部精细搜索以获得高精度的解。因此提出了自适应调整的策略，即随着迭代的进行，线性地减小值。这里提供一个简单常用的方法——线性变化策略：随着迭代次数的增加，惯性权重不断减小，从而使得粒子群算法在初期具有较强的全局收敛能力，在后期具有较强的局部收敛能力。

算法的一些重要概念和技巧

（1）适应值（fitness values）

即优化目标函数的值，用来评价粒子位置的好坏程度，决定是否更新粒子个体的历史最优位置和群体的历史最优位置，保证粒子朝着最优解的方向搜索。

（2）位置限制

限制粒子搜索的空间，即自变量的取值范围，对于无约束问题此处可以省略。

（3）速度限制

为了平衡算法的探索能力与开发能力，需要设定一个合理的速度范围，限制粒子的最大速度 ，即粒子下一步迭代可以移动的最大距离。

• 值较大时，粒子飞行速度快，探索能力强，但粒子容易飞过最优解；
• 值较小时，飞行速度慢，开发能力强，但是收敛速度慢，且容易陷入局部最优解；
• 一般设为粒子变化范围的10％~20％，可根据实际情况调试，但不能大于粒子(解)的变化范围。

（4）优化停止准则

一般有两种：

① 最大迭代步数

② 可接受的满意解：上一次迭代后最优解的适应值与本次迭代后最优解的适应值之差小于某个值后停止优化

（5）初始化

粒子群算法优化的结果受很多因素的影响，其中受粒子初始值的影响比较大，而且较难调控。如果粒子初始值是随机初始化的，在不改变任何参数的情况下，多次优化的结果不一定都收敛到一个全局或局部最优解，也可能会得到一个无效解。所以粒子初始化是一个十分重要的步骤，它关系到整个优化过程中优化收敛的速度与方向。如果粒子的初始化范围选择得好的话可以大大缩短优化的收敛时间，也不易于陷入局部最优解。我们需要根据具体的问题进行分析，如果根据我们的经验判断出最优解一定在某个范围内，则就在这个范围内初始化粒子。如果无法确定，则以粒子的取值边界作为初始化范围。

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;   %种群规模

Vmax = 1;
Vmin = -1;
popmax = 5;    %跟范围设置一致
popmin = -5;

%% IV. 产生初始粒子和速度
for i = 1:sizepop
    % 随机产生一个种群
    pop(i,:) = 5*rands(1,2);    %初始种群，保持范围一致，由于是二元函数，所以随机产生的是两组
    V(i,:) = rands(1,2);  %初始化速度，两组
    % 计算适应度
    fitness(i) = fun(pop(i,:));   %染色体的适应度
end

%% V. 个体极值和群体极值
[bestfitness bestindex] = max(fitness);
zbest = pop(bestindex,:);   %全局最佳
gbest = pop;    %个体最佳
fitnessgbest = fitness;   %个体最佳适应度值
fitnesszbest = bestfitness;   %全局最佳适应度值

%% VI. 迭代寻优
for i = 1:maxgen
    
    for j = 1:sizepop
        % 速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        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)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        
        % 群体最优更新
        if fitness(j) > fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end 
    yy(i) = fitnesszbest;            
end
%% VII.输出结果
[fitnesszbest, zbest]
plot3(zbest(1), zbest(2), fitnesszbest,'bo','linewidth',1.5)

figure
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

3.运行结果