1 算法介绍
2 部分代码
clc; clear; close all; %% 问题定义,这里可以删除 TestProblem=1; % Set to 1, 2, or 3 switch TestProblem case 1 CostFunction=@(x) MyCost1(x); nVar=50; VarMin=0; VarMax=1; case 2 CostFunction=@(x) MyCost2(x); nVar=3; VarMin=-5; VarMax=5; case 3 CostFunction=@(x) MyCost3(x); nVar=2; VarMin=0; VarMax=1; end VarSize=[1 nVar]; VelMax=(VarMax-VarMin)/10; %% MOPSO 设置 nPop=500; % Population Size nRep=100; % Repository Size MaxIt=500; % Maximum Number of Iterations phi1=2.05;%学习因子 phi2=2.05; phi=phi1+phi2; chi=2/(phi-2+sqrt(phi^2-4*phi)); w=chi; % Inertia Weight wdamp=1; % Inertia Weight Damping Ratio c1=chi*phi1; % Personal Learning Coefficient c2=chi*phi2; % Global Learning Coefficient alpha=0.1; % Grid Inflation Parameter nGrid=10; % Number of Grids per each Dimension beta=4; % Leader Selection Pressure Parameter gamma=2; % Extra (to be deleted) Repository Member Selection Pressure %% 初始化 particle=CreateEmptyParticle(nPop); for i=1:nPop particle(i).Velocity=0; particle(i).Position=unifrnd(VarMin,VarMax,VarSize); %unifrnd在[VarMin,VarMax]区间模拟VarSize的数组 % 初始化了自变量X的位置 particle(i).Cost=CostFunction(particle(i).Position); % 求解目标函数 particle(i).Best.Position=particle(i).Position; particle(i).Best.Cost=particle(i).Cost; end particle=DetermineDomination(particle); rep=GetNonDominatedParticles(particle); 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 %% MOPSO 主循环 for it=1:MaxIt for i=1:nPop rep_h=SelectLeader(rep,beta); particle(i).Velocity=w*particle(i).Velocity ... +c1*rand*(particle(i).Best.Position - particle(i).Position) ... +c2*rand*(rep_h.Position - particle(i).Position); particle(i).Velocity=min(max(particle(i).Velocity,-VelMax),+VelMax); particle(i).Position=particle(i).Position + particle(i).Velocity; flag=(particle(i).Position<VarMin | particle(i).Position>VarMax); %非劣检查:flag 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); 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 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 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 %% 结果 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');3 仿真结果
4 参考文献
[1]王越, 吕光宏. 改进的粒子群求解多目标优化算法[J]. 计算机技术与发展, 2014, 000(002):42-45.
博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。
部分理论引用网络文献,若有侵权联系博主删除。