【智能优化算法】基于粒子群结合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

            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代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。