粒子群优化算法(PSO)-MATLAB代码

关于粒子群优化算法(PSO)的介绍与一种C++实现可以参考链接: PSO介绍及其一种C++实现 ,这里不再赘述。

本片博文目的在于提供并简要介绍一种粒子群优化算法(PSO)的MATLAB代码实现。

本文提供的MATLAB代码中,PSO算法本身被封装成一个函数,优化目标函数的句柄作为PSO的输入参数,从而成为了一个较高独立性的函数模块。

以下为pso算法对应的函数输入输出说明。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%函数名:zpso
%函数描述:粒子群优化算法函数
%输入参数:
%   objfunction:最大化目标函数句柄,目标函数是向量到实数的一个多元函数(定义域为超立方)
%   option:规定了pso算法的控制参数,利用ZPSO_parameters函数生成
%   evolutiontimes:粒子群进化代数限制
%输出参数:
%   bestPosition:搜索得到的目标函数最大值点
%   bestFitness:搜索得到的目标函数最大值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bestPosition,bestFitness] = zpso(objfunction,option,evolutiontimes)

zpso函数的输入输入在注释中已经表述得较清楚,不再赘述,需要注意的是objfunction是一个以列向量为输入,标量为输出的目标函数。
以下介绍pso算法控制参数结构体——option 的生成。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%函数名:zpso_parameters
%函数描述:生成PSO算法所需要的参数结构体
%输入参数:
%   dimension:候选解维度
%   minVal:候选解下界,列向量形式
%   maxVal:候选解上界,列向量形式
%   particalCount:粒子数量
%   maxSpeed:粒子最大运动速度
%   globalGuideCoe:全局最优引导加速度因子
%   localGuideCoe:个体最优引导加速度因子
%   disturbanceRate:粒子速度扰动概率
%   maxDisturbanceSpeed:最大扰动速度
%输出参数:
%   op:PSO算法所需参数结构体
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function op = zpso_parameters(dimension,minVal,maxVal,particalCount,maxSpeed,...
                                                              globalGuideCoe,localGuideCoe,...
                                                              disturbanceRate,maxDisturbanceSpeed)

需要说明的是 disturbanceRatemaxDisturbanceSpeed 参数。为了提高算法局部搜索能力,避免加速度因子的选择要求过于严苛,本文实现的pso算法中,粒子运动过程中有一定概率获得一个随机速度。 disturbanceRate 描述了粒子获得随机速度的概率,当其为0,算法退化为最基本的pso算法。而 maxDisturbanceSpeed 则限定了随机速度的最大值。

PSO算法的运行过程如下所示。其中红点代表当前粒子群搜索得到的全局最优解。黑色箭头代表了粒子位置及其运动速度。待优化函数为:

粒子群优化人工神经网络 粒子群优化算法代码_PSO

运行效果:

粒子群优化人工神经网络 粒子群优化算法代码_PSO_02

以下附上与PSO算法相关的MATLAB代码。

<main.m> ——用于测试pso算法:

pso_op = zpso_parameters(2,[-10;-10],[10;10],200,4,1,1,0.1,2);
[bestPosition,bestFitness] = zpso(@objfunction,pso_op,100);

<objfunction.m>——用于定义优化目标函数

function y = objfunction(x)
    x = x - [0.8;-0.4];
    y = sin(sqrt(x'*x))/(sqrt(x'*x)+eps);
end

<zpso_parameters.m>——用于生成pso算法控制参数结构体

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%函数名:zpso_parameters
%函数描述:生成PSO算法所需要的参数结构体
%输入参数:
%   dimension:候选解维度
%   minVal:候选解下界,列向量形式
%   maxVal:候选解上界,列向量形式
%   particalCount:粒子数量
%   maxSpeed:粒子最大运动速度
%   globalGuideCoe:全局最优引导加速度因子
%   localGuideCoe:个体最优引导加速度因子
%   disturbanceRate:粒子速度扰动概率
%   maxDisturbanceSpeed:最大扰动速度
%输出参数:
%   op:PSO算法所需参数结构体
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function op = zpso_parameters(dimension,minVal,maxVal,particalCount,maxSpeed,...
                                                              globalGuideCoe,localGuideCoe,...
                                                              disturbanceRate,maxDisturbanceSpeed)
    %粒子群维度
    op.dimension = dimension;
    %粒子位置下界
    op.minVal = minVal;
    %粒子位置上界
    op.maxVal = maxVal;
    %粒子个数
    op.particalCount = particalCount;
    %粒子最大速度
    op.maxSpeed = maxSpeed;
    %全局最优加速度因子
    op.globalGuideCoe = globalGuideCoe;
    %个体最优加速度因子
    op.localGuideCoe = localGuideCoe;
    %粒子速度扰动概率
    op.disturbanceRate = disturbanceRate;
    %粒子最大扰动速度
    op.maxDisturbanceSpeed = maxDisturbanceSpeed;
end

<zpso.m>——用于实现pso算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%函数名:zpso
%函数描述:粒子群优化算法函数
%输入参数:
%   objfunction:最大化目标函数句柄,目标函数是向量到实数的一个多元函数(定义域为超立方)
%   option:pso算法参数,利用ZPSO_parameters函数生成
%   evolutiontimes:粒子群进化代数限制
%输出参数:
%   bestPosition:搜索得到的目标函数最大值点
%   bestFitness:搜索得到的目标函数最大值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bestPosition,bestFitness] = zpso(objfunction,option,evolutiontimes)

%% 获取pso算法参数
    dimension = option.dimension;
    minVal = option.minVal;
    maxVal = option.maxVal;
    particalCount = option.particalCount;
    globalGuideCoe = option.globalGuideCoe;
    localGuideCoe = option.localGuideCoe;
    maxSpeed = option.maxSpeed;
    disturbanceRate = option.disturbanceRate;
    maxDisturbanceSpeed = option.maxDisturbanceSpeed;

%% 初始化粒子群
    %初始化粒子位置,在[minVal,maxVal]超立方空间内
    pso_position = minVal*ones(1,particalCount);
    pso_position = pso_position+rand(dimension,particalCount).*((maxVal-minVal)*ones(1,particalCount));
    %初始化粒子速度
    pso_velocity = rand(dimension,particalCount); %生成随机速度
    temp = ones(dimension,1) * sqrt(diag(pso_velocity' * pso_velocity))'; %求取随机速度模长
    pso_velocity = pso_velocity./temp; %随机速度模长归一化
    velocityMod = maxSpeed*ones(dimension,1)*rand(1,particalCount); %初始化速度大小随机
    pso_velocity = velocityMod.*pso_velocity;
    %初始化粒子适应度、最佳适应度与个体最优解位置
    pso_fitness = zeros(1,particalCount);
    for j = 1:particalCount
        pso_fitness(j) = objfunction(pso_position(:,j));
    end
    pso_bestposition = pso_position;
    pso_bestfitness = pso_fitness;
    %更新全局最优适应度与最优解位置
    [pso_globalbestfitness,temp] = max(pso_bestfitness);
    pso_globalbestposition = pso_position(:,temp);
    %% 粒子群进化
    for generationCount = 1:evolutiontimes
        %全局最优引导
        pso_velocity = pso_velocity+globalGuideCoe*ones(dimension,1)*rand(1,particalCount)...
                                    .*(pso_globalbestposition*ones(1,particalCount) - pso_position);
        %局部最优引导
        pso_velocity = pso_velocity+localGuideCoe*ones(dimension,1)*rand(1,particalCount)...
                                    .*(pso_bestposition - pso_position);
        %速度扰动
        disturbanceVelocity = rand(dimension,particalCount);
        temp = ones(dimension,1)*sqrt(diag(disturbanceVelocity'*disturbanceVelocity))';
        disturbanceVelocity = disturbanceVelocity./temp;
        temp = maxDisturbanceSpeed*ones(dimension,1)*rand(1,particalCount);
        temp = temp .* (ones(dimension,1)*(rand(1,particalCount)<disturbanceRate));
        disturbanceVelocity = temp.*disturbanceVelocity;
        pso_velocity = pso_velocity+disturbanceVelocity;
        %速度受限
        temp = ones(dimension,1)*sqrt(diag(pso_velocity'*pso_velocity))';
        pso_velocity = pso_velocity./temp;
        temp = min(temp,maxSpeed);
        pso_velocity = temp.*pso_velocity;
        %位置更新
        pso_position = pso_position + pso_velocity;
        %位置受限
        pso_position = max(pso_position,minVal*ones(1,particalCount));
        pso_position = min(pso_position,maxVal*ones(1,particalCount));
        %适应度更新
        for j = 1:particalCount
            pso_fitness(j) = objfunction(pso_position(:,j));
        end
        %最佳适度更新
        index = pso_fitness>pso_bestfitness;
        pso_bestfitness = max(pso_fitness,pso_bestfitness);
        pso_bestposition(:,index) = pso_position(:,index);
        %全局最优更新
        [bestfitness,index] = max(pso_bestfitness);
        index = index(1);
        if(bestfitness > pso_globalbestfitness)
            pso_globalbestfitness = bestfitness;
            pso_globalbestposition = pso_bestposition(:,index);
        end
    end
    %% 输出搜索结果
    bestPosition = pso_globalbestposition;
    bestFitness = pso_globalbestfitness;
end

受作者水平所限,matlab代码可能不尽如人意,仅用于记录、学习与分享。