使用粒子群与NSGA2智能算法解决多目标优化问题的实例教程及Matlab实现代码
1 内容介绍
为解决高度复杂的热电联合经济排放调度问题，本研究提出了一种将非支配排序遗传算法II和多目标粒子群优化算法相结合的协同混合元启发式算法，以经济地运行电力系统并减少环境污染的影响。 .在迭代过程中，根据排名，人口被分成两半。探索是通过非支配排序遗传算法II使用人口的上半部分进行的。通过增加个人学习系数、降低全局学习系数和使用自适应变异算子来修改多目标粒子群优化以有效利用下半部分人口。为了满足线性、非线性约束，并确保人口始终位于热电联产厂的可行运行区域内，开发了一种有效的约束处理机制。所提出的具有有效约束处理机制的混合算法通过有效的信息交换增强了搜索能力。将该算法应用于标准测试功能和测试系统，同时考虑火力发电厂的阀点效应、传输功率损耗、机组边界和热电联产机组的可行运行区域。混合算法可以得到一个分布良好且多样化的帕累托最优解，并且比现有的一些算法更快地收敛到实际的帕累托最优前沿。统计分析表明，所提出的混合算法是解决这一复杂而重要问题的可行替代方案。
2 仿真代码
% Hybrid NSGAII-MOPSO algorithm developed by Dr.Arunachalam Sundaram
clc;
clear;
close all;
global lb ub
%% Problem Definition
CostFunction=@(x)tf1(x); % Cost Function
nVar=2; % Number of Decision Variables
VarSize=[1 nVar]; % Size of Decision Variables Matrix
VarMin=[0 -10^20]; % Lower Bound of Variables
VarMax=[10^20 7]; % Upper Bound of Variables
lb=VarMin;
ub=VarMax;
% Number of Objective Functions
nObj=2;
%% NSGA-II Parameters
MaxIt=120; % Maximum Number of Iterations
nPop=200; % Population Size
pCrossover=0.7; % Crossover Percentage
nCrossover=2*round(pCrossover*nPop/2); % Number of Parnets (Offsprings)
nCrossover=round(nCrossover/2);
pMutation=0.4; % Mutation Percentage
nMutation=round(pMutation*nPop); % Number of Mutants
nMutation=round(nMutation/2);
mug=0.02; % Mutation Rate
sigma=0.1*(VarMax-VarMin); % Mutation Step Size
%% MOPSO Parameters
% MaxIt=2; % Maximum Number of Iterations
% nPop=length(popul1); % Population Size
nRep=50; % Repository Size
w=0.5; % Inertia Weight
wdamp=0.99; % Intertia Weight Damping Rate
c1=4; % Personal Learning Coefficient
c2=1; % 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_individual.Position=[];
empty_individual.Cost=[];
empty_individual.Rank=[];
empty_individual.DominationSet=[];
empty_individual.DominatedCount=[];
empty_individual.CrowdingDistance=[];
pop=repmat(empty_individual,nPop,1);
%% 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=[];
popul=repmat(empty_particle,nPop/2,1);
%%
for i=1:nPop
pop(i).Position=inittf1;
pop(i).Cost=CostFunction(pop(i).Position);
end
% Non-Dominated Sorting
[pop, F]=NonDominatedSorting(pop);
% Calculate Crowding Distance
pop=CalcCrowdingDistance(pop,F);
% Sort Population
[pop, F]=SortPopulation(pop);
% Truncate and divide into two halves
popul1=pop((nPop/2)+1:nPop);
pop=pop(1:nPop/2);
%%
for it=1:MaxIt
% Crossover
popc=repmat(empty_individual,round(nCrossover/2),2);
for k=1:round(nCrossover/2)
i1=randi([1 nPop/2]);
p1=pop(i1);
i2=randi([1 nPop/2]);
p2=pop(i2);
[popc(k,1).Position, popc(k,2).Position]=Crossover(p1.Position,p2.Position);
popc(k,1).Position=lchecktf1(popc(k,1).Position);
popc(k,2).Position=lchecktf1(popc(k,2).Position);
popc(k,1).Cost=CostFunction(popc(k,1).Position);
popc(k,2).Cost=CostFunction(popc(k,2).Position);
end
popc=popc(:);
% Mutation
popm=repmat(empty_individual,nMutation,1);
for k=1:nMutation
i=randi([1 nPop/2]);
p=pop(i);
popm(k).Position=Mutate(p.Position,mug,sigma);
popm(k).Position=lchecktf1(popm(k).Position);
popm(k).Cost=CostFunction(popm(k).Position);
end
% Merge
pop=[pop
popc
popm]; %#ok
% Non-Dominated Sorting
[pop, F]=NonDominatedSorting(pop);
% Calculate Crowding Distance
pop=CalcCrowdingDistance(pop,F);
% Sort Population
pop=SortPopulation(pop);
% Truncate
pop=pop(1:nPop/2);
% Non-Dominated Sorting
[pop, F]=NonDominatedSorting(pop);
% Calculate Crowding Distance
pop=CalcCrowdingDistance(pop,F);
% Sort Population
[pop, F]=SortPopulation(pop);
% Store F1
F1=pop(F{1});
%%
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Number of F1 Members = ' num2str(numel(F1))]);
% Plot F1 Costs
figure(1)
subplot(2,1,1)
PlotCosts(pop);
title('Plot of the population in the search space')
xlabel('F1');
ylabel('F2')
subplot(2,1,2)
PlotCosts(F1);
title('Pareto Front for test function I')
xlabel('F1');
ylabel('F2');
f1(it)=numel(F1);
pause(0.01);
% end
%% MOPSO
for i=1:nPop/2
popul(i).Position=popul1(i).Position;
popul(i).Velocity=zeros(VarSize);
popul(i).Cost=popul1(i).Cost;
% Update Personal Best
popul(i).Best.Position=popul(i).Position;
popul(i).Best.Cost=popul(i).Cost;
end
% Determine Domination
popul=DetermineDomination(popul);
rep=popul(~[popul.IsDominated]);
Grid=CreateGrid(rep,nGrid,alpha);
for i=1:numel(rep)
rep(i)=FindGridIndex(rep(i),Grid);
end
%% MOPSO Main Loop
% for it=1:MaxIt
for i=1:nPop/2
leader=SelectLeader(rep,beta);
popul(i).Velocity = w*popul(i).Velocity ...
+c1*rand(VarSize).*(popul(i).Best.Position-popul(i).Position) ...
+c2*rand(VarSize).*(leader.Position-popul(i).Position);
for j=1:nVar
popul(i).Position(j) = popul(i).Position(j) + popul(i).Velocity(j);
end
popul(i).Position=lchecktf1(popul(i).Position);
popul(i).Cost = CostFunction(popul(i).Position);
% Apply Mutation
pm=(1-(it-1)/(MaxIt-1))^(1/mu);
pm1(it)=pm;
if rand<pm
NewSol.Position=Mutatetf1(popul(i).Position,pm,VarMin,VarMax);
NewSol.Position=lchecktf1(NewSol.Position);
NewSol.Cost=CostFunction(NewSol.Position);
if Dominates(NewSol,popul(i))
popul(i).Position=NewSol.Position;
popul(i).Cost=NewSol.Cost;
elseif Dominates(popul(i),NewSol)
% Do Nothing
else
if rand<0.5
popul(i).Position=NewSol.Position;
popul(i).Cost=NewSol.Cost;
end
end
end
if Dominates(popul(i),popul(i).Best)
popul(i).Best.Position=popul(i).Position;
popul(i).Best.Cost=popul(i).Cost;
elseif Dominates(popul(i).Best,popul(i))
% Do Nothing
else
if rand<0.5
popul(i).Best.Position=popul(i).Position;
popul(i).Best.Cost=popul(i).Cost;
end
end
end
% Add Non-Dominated Particles to REPOSITORY
rep=[rep
popul(~[popul.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(popul,rep);
% pause(0.01);
%
% % Show Iteration Information
% disp(['Iteration ' num2str(it) ': Number of Rep Members = ' num2str(numel(rep))]);
%
% % Damping Inertia Weight
rep1(it)=numel(rep);
w=w*wdamp;
pop1=repmat(empty_individual,nPop/2,1);
for i=1:nPop/2
if i<=length(rep)
pop1(i).Position=rep(i).Position;
pop1(i).Cost=rep(i).Cost;
else
pop1(i).Position=popul(i-length(rep)).Position;
pop1(i).Cost=popul(i-length(rep)).Cost;
end
end
pop2=[pop
pop1
];
% Non-Dominated Sorting
[pop2, F]=NonDominatedSorting(pop2);
% Calculate Crowding Distance
pop2=CalcCrowdingDistance(pop2,F);
% Sort Population
[pop2, F]=SortPopulation(pop2);
% Truncate and divide into two halves
% popul1=pop((nPop/2)+1:nPop)
pop=pop2(1:nPop/2);
% end
end
figure(2)
plot(f1)
title('The nondominated solutions stored in external repository')
xlabel('Iterations');
ylabel('Number of solutions in the repository')
hold on
% figure(3)
plot(rep1,'g')
figure(3)
plot(pm1)
title('The behaviour of the mutation operator when mu=0.1');
xlabel('Iterations');
ylabel('Population in Percentage');
3 运行结果
4 参考文献
[1] Sundaram A . Combined Heat and Power Economic Emission Dispatch using Hybrid NSGA II-MOPSO Algorithm incorporating an Effective Constraint Handling Mechanism[J]. IEEE Access, 2020:1-1.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。