【电力系统经济调度】多元宇宙算法求解电力系统多目标优化问题（Matlab实现）【电气期刊论文复现】

clc;
clear;
close all;
tStart=tic;
% global costdata emissiondata B B0 B00 Pd VarMin VarMax nVar
global data B B0 B00 Pd VarMin VarMax nVar
Pd=2000;
data=xlsread('IEEE10.xls');
B1=xlsread('B10.xls');
B=B1(1:10,1:10);
B0=[0 0 0 0 0 0 0 0 0 0];
B00=0;
%  B=B1(1:3,1:3);
%  B0=[0 0 0 0 0 0 0 0 0 0];
%  B00=0;
%%
Max_time=600; %迭代次数
N=100;
ArchiveMaxSize=100;
% max_iter=Max_time;
 nVar=10;             % 机组个数
VarSize=[1 nVar];   % 决策变量矩阵的大小
VarMin=data(:,2);          %机组出力下限
VarMax= data(:,3);          % 机组出力上限
fobj=@(x) IEEE3aobj(x);
dim=nVar;
lb=VarMin';
ub=VarMax';
obj_no=2;
Best_universe=zeros(1,dim);
Best_universe_Inflation_rate=inf*ones(1,obj_no);
Archive_X=zeros(ArchiveMaxSize,dim);
Archive_F=ones(ArchiveMaxSize,obj_no)*inf;
Archive_member_no=0;
WEP_Max=1;
WEP_Min=0.2;
for i=1:N
   Universes(i,:)=lcheck3; 
end
Time=1;
while Time<Max_time+1
    WEP=WEP_Min+Time*((WEP_Max-WEP_Min)/Max_time);
    TDR=1-((Time)^(1/6)/(Max_time)^(1/6));
    for i=1:size(Universes,1)
        %边界检查（如果宇宙超出边界，则将它们带回搜索空间内）
        Flag4ub=Universes(i,:)>ub;
        Flag4lb=Universes(i,:)<lb;
        Universes(i,:)=(Universes(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        Universes(i,:)=lbcoff3bus(Universes(i,:));
        %计算宇宙的通货膨胀率(适合度)
        Inflation_rates(i,:)=fobj(Universes(i,:));
        %精英主义
        if dominates(Inflation_rates(i,:),Best_universe_Inflation_rate)
            Best_universe_Inflation_rate=Inflation_rates(i,:);
            Best_universe=Universes(i,:);
        end
    end
    [sorted_Inflation_rates,sorted_indexes]=sort(Inflation_rates);
    for newindex=1:N
        Sorted_universes(newindex,:)=Universes(sorted_indexes(newindex),:);
    end
    %原始MVO论文中的标准化通货膨胀率
    normalized_sorted_Inflation_rates=normr(sorted_Inflation_rates);
    Universes(1,:)= Sorted_universes(1,:);
%     Universes(1,:)=lchecktf1(Universes(1,:));
    [Archive_X, Archive_F, Archive_member_no]=UpdateArchive(Archive_X, Archive_F, Universes, Inflation_rates, Archive_member_no);
    if Archive_member_no>ArchiveMaxSize
        Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
        [Archive_X, Archive_F, Archive_mem_ranks, Archive_member_no]=HandleFullArchive(Archive_X, Archive_F, Archive_member_no, Archive_mem_ranks, ArchiveMaxSize);
    else
        Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
    end
    Archive_mem_ranks=RankingProcess(Archive_F, ArchiveMaxSize, obj_no);
    % 提高复盖率
    index=RouletteWheelSelection(1./Archive_mem_ranks);
    if index==-1
        index=1;
    end
   Best_universe_Inflation_rate=Archive_F(index,:);
   Best_universe=Archive_X(index,:); 
    %更新宇宙的位置
    for i=2:size(Universes,1)%从2开始，因为第1位是精英
        Back_hole_index=i;
        for j=1:size(Universes,2)
            r1=rand();
            if r1<normalized_sorted_Inflation_rates(i)
                White_hole_index=RouletteWheelSelection(-sorted_Inflation_rates);% 对于最大化问题，排序的通货膨胀率应该写成排序的通货膨胀率
                if White_hole_index==-1
                    White_hole_index=1;
                end
                %Eq. (3.1) 
                Universes(Back_hole_index,j)=Sorted_universes(White_hole_index,j);
%                 Universes(Back_hole_index,j)=lchecktf1(Universes(Back_hole_index,j));
            end
            if (size(lb',1)==1)
                %如果边界都是一样的，那么原MVO论文中的公式(3.2)就会出现
                r2=rand();
                if r2<WEP
                    r3=rand();
                    if r3<0.5
                        Universes(i,j)=Best_universe(1,j)+TDR*((ub-lb)*rand+lb);
                    end
                    if r3>0.5
                        Universes(i,j)=Best_universe(1,j)-TDR*((ub-lb)*rand+lb);
                    end
                end
            end
            if (size(lb',1)~=1)
            %公式( 3.2 )在原始MVO论文中，如果对每个变量的上下界不同
                r2=rand();
                if r2<WEP
                    r3=rand();
                    if r3<0.5
                        Universes(i,j)=Best_universe(1,j)+TDR*((ub(j)-lb(j))*rand+lb(j));
                    end
                    if r3>0.5
                        Universes(i,j)=Best_universe(1,j)-TDR*((ub(j)-lb(j))*rand+lb(j));
                    end
                end
            end
        end
        Universes(i,:)=lbcoff3bus(Universes(i,:));
    end
    display(['At the iteration ', num2str(Time), ' there are ', num2str(Archive_member_no), ' non-dominated solutions in the archive']);
    Time=Time+1;
%
end
plot(Archive_F(:,1),Archive_F(:,2),'Ro','LineWidth',2,...
        'MarkerEdgeColor','r',...
        'MarkerFaceColor','r',...
        'MarkerSize',2);
xlabel('污染排放量')
ylabel('煤耗量')
title('Pareto最前沿')
% Universes
Archive_F(:,1)
Archive_F(:,2)
Best_universe