一、获取代码方式
获取代码方式1：
完整代码已上传我的资源：【优化算法】多目标粒子群优化算法（MOPSO）【含Matlab源码 033期】
获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。
备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；
二、粒子群算法简介
1 粒子群算法的概念
粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解．
PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。
2 粒子群算法分析
2.1基本思想
粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。下面的动图很形象地展示了PSO算法的过程：
2 更新规则
PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。
公式(1)的第一部分称为【记忆项】，表示上次速度大小和方向的影响；公式(1)的第二部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；公式(1)的第三部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。以上面两个公式为基础，形成了PSO的标准形式。
公式(2)和 公式(3)被视为标准PSO算法。
3 PSO算法的流程和伪代码
三、部分源代码
%
%
%
mopso;
clc;
clear;
close all;
%% Problem Definition
CostFunction=@(x) ZDT(x); % Cost Function
nVar=5; % Number of Decision Variables
VarSize=[1 nVar]; % Size of Decision Variables Matrix
VarMin=0; % Lower Bound of Variables
VarMax=1; % Upper Bound of Variables
%% MOPSO Parameters
MaxIt=200; % Maximum Number of Iterations
nPop=200; % Population Size
nRep=100; % Repository Size 库的大小
w=0.5; % Inertia Weight
wdamp=0.99; % Intertia Weight Damping Rate 惯性重量阻尼率
c1=1; % Personal Learning Coefficient
c2=2; % Global Learning Coefficient
nGrid=7; % Number of Grids per Dimension 每个维度的网格数
alpha=0.1; % Inflation Rate 通货膨胀率
beta=2; % Leader Selection Pressure 领导人选择压力
gamma=2; % Deletion Selection Pressure 删除选择压力
mu=0.1; % Mutation Rate 突变率
%% Initialization 初始化
empty_particle.Position=[];
empty_particle.Velocity=[];
empty_particle.Cost=[];
empty_particle.Best.Position=[];
empty_particle.Best.Cost=[];
empty_particle.IsDominated=[];
empty_particle.GridIndex=[];
empty_particle.GridSubIndex=[]; %网格的子索引
pop=repmat(empty_particle,nPop,1);
for i=1:nPop
pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
pop(i).Velocity=zeros(VarSize);
pop(i).Cost=CostFunction(pop(i).Position);
% Update Personal Best 更新个人最好
pop(i).Best.Position=pop(i).Position;
pop(i).Best.Cost=pop(i).Cost;
end
% Determine Domination 确定统治
pop=DetermineDomination(pop);
rep=pop(~[pop.IsDominated]);
Grid=CreateGrid(rep,nGrid,alpha);
for i=1:numel(rep)
rep(i)=FindGridIndex(rep(i),Grid);
end
%% MOPSO Main Loop MOPSO主循环
for it=1:MaxIt
for i=1:nPop
leader=SelectLeader(rep,beta);
pop(i).Velocity = w*pop(i).Velocity ...
+c1*rand(VarSize).*(pop(i).Best.Position-pop(i).Position) ...
+c2*rand(VarSize).*(leader.Position-pop(i).Position);
pop(i).Position = pop(i).Position + pop(i).Velocity;
pop(i).Position = max(pop(i).Position, VarMin);
pop(i).Position = min(pop(i).Position, VarMax);
pop(i).Cost = CostFunction(pop(i).Position);
% Apply Mutation 应用突变
pm=(1-(it-1)/(MaxIt-1))^(1/mu);
if rand<pm
NewSol.Position=Mutate(pop(i).Position,pm,VarMin,VarMax);
NewSol.Cost=CostFunction(NewSol.Position);
if Dominates(NewSol,pop(i))
pop(i).Position=NewSol.Position;
pop(i).Cost=NewSol.Cost;
elseif Dominates(pop(i),NewSol)
% Do Nothing
else
if rand<0.5
pop(i).Position=NewSol.Position;
pop(i).Cost=NewSol.Cost;
end
end
end
if Dominates(pop(i),pop(i).Best)
pop(i).Best.Position=pop(i).Position;
pop(i).Best.Cost=pop(i).Cost;
elseif Dominates(pop(i).Best,pop(i))
% Do Nothing
else
if rand<0.5
pop(i).Best.Position=pop(i).Position;
pop(i).Best.Cost=pop(i).Cost;
end
end
end
% Add Non-Dominated Particles to REPOSITORY 向存储库添加非支配粒子
rep=[rep
pop(~[pop.IsDominated])]; %#ok
% Determine Domination of New Resository Members 决定新成员的支配地位
rep=DetermineDomination(rep);
% Keep only Non-Dminated Memebrs in the Repository 在存储库中只保留非支配成员
rep=rep(~[rep.IsDominated]);
% Update Grid 更新网格
Grid=CreateGrid(rep,nGrid,alpha);
% Update Grid Indices 更新网格索引
for i=1:numel(rep)
rep(i)=FindGridIndex(rep(i),Grid);
end
% Check if Repository is Full 检查存储库是否已满
if numel(rep)>nRep
Extra=numel(rep)-nRep;
for e=1:Extra
rep=DeleteOneRepMemebr(rep,gamma);
end
end
% Plot Costs 绘制成本图
figure(1);
PlotCosts(pop,rep);
pause(0.01);
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Number of Rep Members = ' num2str(numel(rep))]);
% Damping Inertia Weight 显示迭代信息
w=w*wdamp;
end
%% Resluts四、运行结果
五、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
