多目标粒子群算法可以出pareto图
1）函数CreateEmptyParticle：构建粒子结构体
function particle=CreateEmptyParticle(n) if nargin<1n=1;endempty_particle.Position=[];empty_particle.Velocity=[];empty_particle.Cost=[];empty_particle.Dominated=false;empty_particle.Best.Position=[];empty_particle.Best.Cost=[];empty_particle.GridIndex=[];empty_particle.GridSubIndex=[];%%函数B = repmat(A,m,n):将矩阵A复制m×n块，即B由m×n块A平铺而成。particle=repmat(empty_particle,n,1);end2）初始化粒子群
for i=1:nPop
particle(i).Velocity=0;
%%unifrnd函数生成区间[a,b]上连续型均匀分布的m行n列的随机数矩阵。particle(i).Position=unifrnd(VarMin,VarMax,VarSize);particle(i).Cost=CostFunction(particle(i).Position);particle(i).Best.Position=particle(i).Position;particle(i).Best.Cost=particle(i).Cost;end3）初始化rep
particle=DetermineDomination(particle);
rep=GetNonDominatedParticles(particle);函数DetermineDomination：判断粒子群的非支配性
function pop=DetermineDomination(pop)npop=numel(pop); %% numel元素的数目 for i=1:npoppop(i).Dominated=false;for j=1:i-1if ~pop(j).Dominatedif Dominates(pop(i),pop(j))pop(j).Dominated=true;elseif Dominates(pop(j),pop(i))pop(i).Dominated=true;break;endendendend
endfunction dom=Dominates(x,y)%%函数isstruct(x) Determine if input is a MATLAB structure arrayif isstruct(x)x=x.Cost;endif isstruct(y)y=y.Cost;enddom=all(x<=y) && any(x<y);end函数GetNonDominatedParticles(particle)为获取群中非支配粒子
function nd_pop=GetNonDominatedParticles(pop)ND=~[pop.Dominated]; nd_pop=pop(ND);
end4）构建立方体，并确定rep中粒子的位置。
rep_costs=GetCosts(rep);
G=CreateHypercubes(rep_costs,nGrid,alpha);
for i=1:numel(rep)[rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);
End函数说明：
function costs=GetCosts(pop)nobj=numel(pop(1).Cost);costs=reshape([pop.Cost],nobj,[]);%%reshape转化矩阵的格式endfunction G=CreateHypercubes(costs,ngrid,alpha)nobj=size(costs,1); empty_grid.Lower=[];empty_grid.Upper=[];G=repmat(empty_grid,nobj,1);for j=1:nobj min_cj=min(costs(j,:));max_cj=max(costs(j,:)); dcj=alpha*(max_cj-min_cj); min_cj=min_cj-dcj;max_cj=max_cj+dcj;%%linespace(x1,x2,N),以x1为起始元素，x2为结尾元素，生成等间距的列向量。gx=linspace(min_cj,max_cj,ngrid-1); G(j).Lower=[-inf gx];G(j).Upper=[gx inf]; end
endfunction [Index SubIndex]=GetGridIndex(particle,G)c=particle.Cost; nobj=numel(c);
ngrid=numel(G(1).Upper); %%sub2ind（size,I,j）函数表示将矩阵的下标（I,j）转换为matlab的矩阵索引。 str=['sub2ind(' mat2str(ones(1,nobj)*ngrid)];SubIndex=zeros(1,nobj);for j=1:nobj U=G(j).Upper; i=find(c(j)<U,1,'first'); SubIndex(j)=i; str=[str ',' num2str(i)];end str=[str ');'];%% eval('str')的作用是将字符串按matlab中的命令执行，也就是相当于在matlab主窗口中输入 str 再运行。（str代表特定的命令字符串）Index=eval(str);
End5）MOPSO主循环程序
for it=1:MaxItfor i=1:nPop %%a) selec the leaderrep_h=SelectLeader(rep,beta);%%b) compute the new Velocity and positionparticle(i).Velocity=w*particle(i).Velocity ...+c1*rand*(particle(i).Best.Position - particle(i).Position) ...+c2*rand*(rep_h.Position - particle(i).Position);%%c) matain the new Velocity in the search spaceparticle(i).Velocity=min(max(particle(i).Velocity,-VelMax),+VelMax);particle(i).Position=particle(i).Position + particle(i).Velocity;%%flag为logical（1*n）行向量flag=(particle(i).Position<VarMin | particle(i).Position>VarMax);particle(i).Velocity(flag)=-particle(i).Velocity(flag); particle(i).Position=min(max(particle(i).Position,VarMin),VarMax);particle(i).Cost=CostFunction(particle(i).Position);%%d) update the PBEST(i)if Dominates(particle(i),particle(i).Best)particle(i).Best.Position=particle(i).Position;particle(i).Best.Cost=particle(i).Cost; elseif ~Dominates(particle(i).Best,particle(i))if rand<0.5 %%随机选择particle(i).Best.Position=particle(i).Position;particle(i).Best.Cost=particle(i).Cost;endendend%%f) update the the REPparticle=DetermineDomination(particle);nd_particle=GetNonDominatedParticles(particle); rep=[rep;nd_particle]; rep=DetermineDomination(rep);rep=GetNonDominatedParticles(rep); for i=1:numel(rep)[rep(i).GridIndex rep(i).GridSubIndex]=GetGridIndex(rep(i),G);end%7)if the rep is full,delete the particle in highly populated regionsif numel(rep)>nRepEXTRA=numel(rep)-nRep;rep=DeleteFromRep(rep,EXTRA,gamma); rep_costs=GetCosts(rep);G=CreateHypercubes(rep_costs,nGrid,alpha); end disp(['Iteration ' num2str(it) ': Number of Repository Particles = ' num2str(numel(rep))]); w=w*wdamp; %%惯性因子线性下降end函数说明：
function rep_h=SelectLeader(rep,beta)if nargin<2beta=1;end%%1）get the ocuupied cells information（index and count）
[occ_cell_index occ_cell_member_count]=GetOccupiedCells(rep); %%2）select the hepercube using the RouletteWheelSelectionp=occ_cell_member_count.^(-beta);p=p/sum(p); selected_cell_index=occ_cell_index(RouletteWheelSelection(p)); GridIndices=[rep.GridIndex]; selected_cell_members=find(GridIndices==selected_cell_index);%%3）randomly select the cell_membern=numel(selected_cell_members); selected_memebr_index=randint(1,1,[1 n]); h=selected_cell_members(selected_memebr_index); rep_h=rep(h);end函数GetOccupied：get the ocuupied cells information
function [occ_cell_index occ_cell_member_count]=GetOccupiedCells(pop)GridIndices=[pop.GridIndex];
%%unique函数取集合的不重复元素。
occ_cell_index=unique(GridIndices); occ_cell_member_count=zeros(size(occ_cell_index));m=numel(occ_cell_index);
for k=1:m%%sum(A>0)，A>0的结果是得到一个逻辑矩阵，大小跟原来的A一致，sum求和即A中每列大于零的元素个数occ_cell_member_count(k)=sum(GridIndices==occ_cell_index(k));end
end函数RouletteWheelSelection轮盘赌方法选择
function i=RouletteWheelSelection(p)r=rand;c=cumsum(p);%% cumsum(A)求取每列累积和i=find(r<=c,1,'first');
end函数DeleteFromRep 裁剪REP集合
function rep=DeleteFromRep(rep,EXTRA,gamma)if nargin<3gamma=1;endfor k=1:EXTRA[occ_cell_index occ_cell_member_count]=GetOccupiedCells(rep);p=occ_cell_member_count.^gamma;p=p/sum(p);selected_cell_index=occ_cell_index(RouletteWheelSelection(p));GridIndices=[rep.GridIndex];selected_cell_members=find(GridIndices==selected_cell_index);n=numel(selected_cell_members);selected_memebr_index=randint(1,1,[1 n]);j=selected_cell_members(selected_memebr_index); rep=[rep(1:j-1); rep(j+1:end)];endend6）结果显示
%% Results
costs=GetCosts(particle);
rep_costs=GetCosts(rep);figure;
plot(costs(1,:),costs(2,:),'b.');
hold on;
plot(rep_costs(1,:),rep_costs(2,:),'rx');
legend('Main Population','Repository');完整源码 多目标粒子群算法可以出pareto图