多目标粒子群算法可以出pareto图
1）函数CreateEmptyParticle：构建粒子结构体
function particle=CreateEmptyParticle(n)
if nargin<1
n=1;
end
empty_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);
end
2）初始化粒子群
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;
end
3）初始化rep
particle=DetermineDomination(particle);
rep=GetNonDominatedParticles(particle);函数DetermineDomination：判断粒子群的非支配性
function pop=DetermineDomination(pop)
npop=numel(pop); %% numel元素的数目
for i=1:npop
pop(i).Dominated=false;
for j=1:i-1
if ~pop(j).Dominated
if Dominates(pop(i),pop(j))
pop(j).Dominated=true;
elseif Dominates(pop(j),pop(i))
pop(i).Dominated=true;
break;
end
end
end
end
end
function dom=Dominates(x,y)
%%函数isstruct(x) Determine if input is a MATLAB structure array
if isstruct(x)
x=x.Cost;
end
if isstruct(y)
y=y.Cost;
end
dom=all(x<=y) && any(x<y);
end函数GetNonDominatedParticles(particle)为获取群中非支配粒子
function nd_pop=GetNonDominatedParticles(pop)
ND=~[pop.Dominated];
nd_pop=pop(ND);
end
4）构建立方体，并确定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转化矩阵的格式
end
function 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
end
function [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);
End
5）MOPSO主循环程序
for it=1:MaxIt
for i=1:nPop
%%a) selec the leader
rep_h=SelectLeader(rep,beta);
%%b) compute the new Velocity and position
particle(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 space
particle(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;
end
end
end
%%f) update the the REP
particle=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 regions
if numel(rep)>nRep
EXTRA=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<2
beta=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 RouletteWheelSelection
p=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_member
n=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<3
gamma=1;
end
for 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)];
end
end6）结果显示
%% 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');完整源码
浙公网安备 33010602011771号