精华内容
下载资源
问答
  • 多目标优化算法NSGAⅢMATLAB代码及注释
  • 多目标优化算法NSGAⅡ MATLAB代码及详细注释
  • Matlab编写多目标优化算法NSGA-Ⅱ
  • 可运行的NSGA-2优化算法,可自由设置目标函数及约束函数,作者几乎在每行代码后面加入了中文注解,以便使用者更好的理解算法原理。
  • 基于matlab的遗传算法的非支配排序遗传算法,多目标优化算法 可以用于多目标优化变量回归,求解最优值
  • 多目标优化免费NSGA-II代码+详细解释(详见文章)该函数基于求解目标最优解的进化算法,即目标的帕累托前沿。最初只输入种群大小和回采标准,或算法自动停止的总代数。您将被要求输入目标函数的数量、决策变量的...
  • 本资源可以用于目标函数以及个变量,例如三目标三变量
  • 1)本程序主要针对测试函数集ZDT1进行的NSGA-Ⅱ算法的编写; 2)本程序有详细的备注解释; 3)本文件里包含论文《非支配排序遗传算法(NSGA)的研究与应用》.pdf,用来指导学习NSGA-Ⅱ算法
  • 收敛性极佳的经典多目标优化算法NSGA2源码。用c++实现的。
  • 本代码资源是关于NSGA-2的Python实现,是原始论文A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II的代码复现结果,有整个NSGA-2的实现流程包括初始化种群,基因生成,染色体交叉变异等
  • NSGA2优化算法Matlab求解多目标优化问题,遗传算法优化+帕累托排序,有效地解决了多目标优化问题,算例可行有效。
  • 资源整理不易,欢迎下载交流学习! NSGA2优化算法Matlab求解多目标优化问题,遗传算法优化+帕累托排序,有效地解决了多目标优化问题,算例可行有效。
  • 考虑自平衡能力的并网型微电网多目标容量优化设计nsga2算法matlab实现
  • 多目标优化算法(一)NSGA-Ⅱ(NSGA2)

    万次阅读 多人点赞 2018-09-28 20:39:42
    多目标优化算法(一)NSGA-Ⅱ 0.前言 这个算法是本人接触科研学习实现的第一个算法,因此想在这里和大家分享一下心得。...Pareto支配关系:对于最小化多目标优化问题,对于n个目标分量 fi(x),i=1...nf_i(...

    多目标优化算法(一)NSGA-Ⅱ(NSGA2)

    注:没有想到这篇博客竟然有很多人查看,这是我去年写的算法,里面难免会有一些错误,大家可以看看评论区,这里的matlab代码写的不太好,是以C语言思维写的,基本上没有并行,初学者建议可以看看platEMO上的源代码,提前培养好写代码的习惯!

    0. 前言

    这个算法是本人接触科研学习实现的第一个算法,因此想在这里和大家分享一下心得。

    1. 算法简介

    NSGA-Ⅱ算法,即带有精英保留策略的快速非支配多目标优化算法,是一种基于Pareto最优解的多目标优化算法。

    1.1 Pareto支配关系以及Pareto等级

    Pareto支配关系:对于最小化多目标优化问题,对于n个目标分量 f i ( x ) , i = 1... n f_i(x), i=1...n fi(x),i=1...n,任意给定两个决策变量 X a X_a Xa X b X_b Xb,如果有以下两个条件成立,则称 X a X_a Xa支配 X b X_b Xb
    1.对于 ∀ i ∈ 1 , 2 , . . . , n \forall i \in {1,2,...,n} i1,2,...,n,都有 f i ( X a ) ≤ f i ( X b ) f_i(X_a) \leq f_i(X_b) fi(Xa)fi(Xb)成立。
    2. ∃ i ∈ 1 , 2 , . . . , n \exists i \in {1,2,...,n} i1,2,...,n,使得 f i ( X a ) < f i ( X b ) f_i(X_a) <f_i(X_b) fi(Xa)<fi(Xb)成立。
    如果对于一个决策变量,不存在其他决策变量能够支配他,那么就称该决策变量为非支配解。
    Pareto等级:在一组解中,非支配解Pareto等级定义为1,将非支配解从解的集合中删除,剩下解的Pareto等级定义为2,依次类推,可以得到该解集合中所有解的Pareto等级。示意图如图1所示。
    在这里插入图片描述

    1.2 快速非支配排序

    假设种群大小为P,该算法需要计算每个个体p的被支配个数 n p n_p np和该个体支配的解的集合 S p S_p Sp这两个参数。遍历整个种群,该参数的计算复杂度为 O ( m N 2 ) O(mN^2) O(mN2)。该算法的伪代码如下:
    1.计算出种群中每个个体的两个参数 n p n_p np S p S_p Sp
    2.将种群中参数 n p = 0 n_p=0 np=0的个体放入集合 F 1 F_1 F1中。
    3.for 个体 i ∈ F 1 i \in F_1 iF1:
            for 个体 l ∈ S i l \in S_i lSi
                    n l = n l − 1 n_l=n_l-1 nl=nl1;(该步骤即消除Pareto等级1对其余个体的支配,相当于将Pareto等级1的个体从集合中删除)
                    if  n l = 0 n_l=0 nl=0
                           将个体l加入集合 F 2 F_2 F2中。
                    end
            end
    end
    4.上面得到Pareto等级2的个体的集合 F 2 F_2 F2,对集合 F 2 F_2 F2中的个体继续重复步骤3,依次类推直到种群等级被全部划分。
    matlab 代码如下:

    function [F,chromo] = non_domination_sort( pop,chromo,f_num,x_num )
    %non_domination_sort 初始种群的非支配排序和计算拥挤度
    %初始化pareto等级为1
    pareto_rank=1;
    F(pareto_rank).ss=[];%pareto等级为pareto_rank的集合
    p=[];%每个个体p的集合
    for i=1:pop
        %%%计算出种群中每个个体p的被支配个数n和该个体支配的解的集合s
        p(i).n=0;%被支配个体数目n
        p(i).s=[];%支配解的集合s
        for j=1:pop
            less=0;%y'的目标函数值小于个体的目标函数值数目
            equal=0;%y'的目标函数值等于个体的目标函数值数目
            greater=0;%y'的目标函数值大于个体的目标函数值数目
            for k=1:f_num
                if(chromo(i,x_num+k)<chromo(j,x_num+k))
                    less=less+1;
                elseif(chromo(i,x_num+k)==chromo(j,x_num+k))
                    equal=equal+1;
                else
                    greater=greater+1;
                end
            end
            if(less==0 && equal~=f_num)%if(less==0 && greater~=0)
                p(i).n=p(i).n+1;%被支配个体数目n+1
            elseif(greater==0 && equal~=f_num)%elseif(greater==0 && less~=0)
                p(i).s=[p(i).s j];
            end
        end
        %%%将种群中参数为n的个体放入集合F(1)if(p(i).n==0)
            chromo(i,f_num+x_num+1)=1;%储存个体的等级信息
            F(pareto_rank).ss=[F(pareto_rank).ss i];
        end
    end
    %%%求pareto等级为pareto_rank+1的个体
    while ~isempty(F(pareto_rank).ss)
        temp=[];
        for i=1:length(F(pareto_rank).ss)
            if ~isempty(p(F(pareto_rank).ss(i)).s)
                for j=1:length(p(F(pareto_rank).ss(i)).s)
                    p(p(F(pareto_rank).ss(i)).s(j)).n=p(p(F(pareto_rank).ss(i)).s(j)).n - 1;%nl=nl-1
                    if p(p(F(pareto_rank).ss(i)).s(j)).n==0
                        chromo(p(F(pareto_rank).ss(i)).s(j),f_num+x_num+1)=pareto_rank+1;%储存个体的等级信息
                        temp=[temp p(F(pareto_rank).ss(i)).s(j)];
                    end
                end
            end
        end
        pareto_rank=pareto_rank+1;
        F(pareto_rank).ss=temp;
    end
    

    1.3 拥挤度

    为了使得到的解在目标空间中更加均匀,这里引入了拥挤度 n d n_d nd,其算法如下:

    1. 令参数 n d = 0 , n ∈ 1 … N n_d=0,n∈1…N nd=0n1N
    2. for 每个目标函数 f m f_m fm
              ① 根据该目标函数对该等级的个体进行排序,记 f m m a x f_m^{max} fmmax为个体目标函数值 f m f_m fm的最大值, f m m i n f_m^{min} fmmin为个体目标函数值 f m f_m fm的最小值;
              ② 对于排序后两个边界的拥挤度 1 d 1_d 1d N d N_d Nd置为∞;
              ③ 计算 n d = n d + ( f m ( i + 1 ) − f m ( i − 1 ) ) / ( f m m a x − f m m i n ) n_d=n_d+(f_m (i+1)-f_m (i-1))/(f_m^{max}-f_m^{min}) nd=nd+(fm(i+1)fm(i1))/(fmmaxfmmin),其中 f m ( i + 1 ) f_m (i+1) fm(i+1)是该个体排序后后一位的目标函数值。

        end

          从二目标优化问题来看,就像是该个体在目标空间所能生成的最大的矩形(该矩形不能触碰目标空间其他的点)的边长之和。拥挤度示意图如图2所示。
    拥挤度示意图
    matlab代码如下:

    function chromo = crowding_distance_sort( F,chromo,f_num,x_num )
    %计算拥挤度
    %%%按照pareto等级对种群中的个体进行排序
    [~,index]=sort(chromo(:,f_num+x_num+1));
    [~,mm1]=size(chromo);
    temp=zeros(length(index),mm1);
    for i=1:length(index)%=pop
        temp(i,:)=chromo(index(i),:);%按照pareto等级排序后种群
    end
    
    %%%对于每个等级的个体开始计算拥挤度
    current_index = 0;
    for pareto_rank=1:(length(F)-1)%计算F的循环时多了一次空,所以减掉
        %%拥挤度初始化为0
        nd=[];
        nd(:,1)=zeros(length(F(pareto_rank).ss),1);
        %y=[];%储存当前处理的等级的个体
        [~,mm2]=size(temp);
        y=zeros(length(F(pareto_rank).ss),mm2);%储存当前处理的等级的个体
        previous_index=current_index + 1;
        for i=1:length(F(pareto_rank).ss)
            y(i,:)=temp(current_index + i,:);
        end
        current_index=current_index + i;
        %%对于每一个目标函数fm
        for i=1:f_num
            %%根据该目标函数值对该等级的个体进行排序
            [~,index_objective]=sort(y(:,x_num+i));
            %objective_sort=[];%通过目标函数排序后的个体
            [~,mm3]=size(y);
            objective_sort=zeros(length(index_objective),mm3);%通过目标函数排序后的个体
            for j=1:length(index_objective)
                objective_sort(j,:)=y(index_objective(j),:);
            end
            %%记fmax为最大值,fmin为最小值
            fmin=objective_sort(1,x_num+i);
            fmax=objective_sort(length(index_objective),x_num+i);
            %%对排序后的两个边界拥挤度设为1d和nd设为无穷
            y(index_objective(1),f_num+x_num+1+i)=Inf;
            y(index_objective(length(index_objective)),f_num+x_num+1+i)=Inf;
            %%计算nd=nd+(fm(i+1)-fm(i-1))/(fmax-fmin)
            for j=2:(length(index_objective)-1)
                pre_f=objective_sort(j-1,x_num+i);
                next_f=objective_sort(j+1,x_num+i);
                if (fmax-fmin==0)
                    y(index_objective(j),f_num+x_num+1+i)=Inf;
                else
                    y(index_objective(j),f_num+x_num+1+i)=(next_f-pre_f)/(fmax-fmin);
                end
            end
        end
        %多个目标函数拥挤度求和
        for i=1:f_num
            nd(:,1)=nd(:,1)+y(:,f_num+x_num+1+i);
        end
        %2列保存拥挤度,其他的覆盖掉
        y(:,f_num+x_num+2)=nd;
        y=y(:,1:(f_num+x_num+2));
        temp_two(previous_index:current_index,:) = y;
    end
    chromo=temp_two;
    end
    

    1.4 精英保留策略

    1.首先将父代种群 C i C_i Ci和子代种群 D i D_i Di合成种群 R i R_i Ri
    2.根据以下规则从种群 R i R_i Ri生成新的父代种群 C i + 1 C_{i+1} Ci+1
          ①根据Pareto等级从低到高的顺序,将整层种群放入父代种群 C i + 1 C_{i+1} Ci+1,直到某一层该层个体不能全部放入父代种群 C i + 1 C_{i+1} Ci+1
          ②将该层个体根据拥挤度从大到小排列,依次放入父代种群 C i + 1 C_{i+1} Ci+1中,直到父代种群 C i + 1 C_{i+1} Ci+1填满。
    matlab代码如下:

    function chromo = elitism( pop,combine_chromo2,f_num,x_num )
    %精英保留策略
    [pop2,temp]=size(combine_chromo2);
    chromo_rank=zeros(pop2,temp);
    chromo=zeros(pop,(f_num+x_num+2));
    %根据pareto等级从高到低进行排序
    [~,index_rank]=sort(combine_chromo2(:,f_num+x_num+1));
    for i=1:pop2
        chromo_rank(i,:)=combine_chromo2(index_rank(i),:);
    end
    %求出最高的pareto等级
    max_rank=max(combine_chromo2(:,f_num+x_num+1));
    %根据排序后的顺序,将等级相同的种群整个放入父代种群中,直到某一层不能全部放下为止
    prev_index=0;
    for i=1:max_rank
        %寻找当前等级i个体里的最大索引
        current_index=max(find(chromo_rank(:,(f_num+x_num+1))==i));
        %不能放下为止
        if(current_index>pop)
            %剩余群体数
            remain_pop=pop-prev_index;
            %等级为i的个体
            temp=chromo_rank(((prev_index+1):current_index),:);
            %根据拥挤度从大到小排序
            [~,index_crowd]=sort(temp(:,(f_num+x_num+2)),'descend');
            %填满父代
            for j=1:remain_pop
                chromo(prev_index+j,:)=temp(index_crowd(j),:);
            end
            return;
        elseif(current_index<pop)
            % 将所有等级为i的个体直接放入父代种群
            chromo(((prev_index+1):current_index),:)=chromo_rank(((prev_index+1):current_index),:);
        else
            % 将所有等级为i的个体直接放入父代种群
            chromo(((prev_index+1):current_index),:)=chromo_rank(((prev_index+1):current_index),:);
            return;
        end
        prev_index = current_index;
    end
    end
    

    1.5 实数编码的交叉操作(SBX)

    模拟二进制交叉:
    x 1 j ( t ) = 0.5 × [ ( 1 + γ j ) x 1 j ( t ) + ( 1 − γ j ) x 2 j ( t ) ] x _{1j}(t)=0.5×[(1+γ_j ) x_{1j}(t)+(1-γ_j ) x_{2j} (t)] x1j(t)=0.5×[(1+γj)x1j(t)+(1γj)x2j(t)]
    x 2 j ( t ) = 0.5 × [ ( 1 − γ j ) x 1 j ( t ) + ( 1 + γ j ) x 2 j ( t ) ] x _{2j} (t)=0.5×[(1-γ_j ) x_{1j}(t)+(1+γ_j ) x_{2j}(t)] x2j(t)=0.5×[(1γj)x1j(t)+(1+γj)x2j(t)]
    其中
    γ j = { ( 2 u j ) 1 / ( η + 1 )                          u j < 0.5 ( 1 / ( 2 ( 1 − u j ) ) 1 / ( η + 1 )          e l s e γ_j=\left\{\begin{matrix}(2u_j)^{1/(η+1)}\;\;\;\;\;\;\;\;\;\;\;\;u_j<0.5 \\ (1/(2(1-u_j))^{1/(η+1)}\;\;\;\;else \end{matrix}\right. γj={(2uj)1/(η+1)uj<0.5(1/(2(1uj))1/(η+1)else

    1.6 多项式变异(polynomial mutation)

    多项式变异:
    x 1 j ( t ) = x 1 j ( t ) + ∆ j x _{1j} (t)=x_{1j} (t)+∆_j x1j(t)=x1j(t)+j
    其中
    ∆ j = { ( 2 u j ) 1 / ( η + 1 ) − 1                          u j < 0.5 ( 1 − ( 2 ( 1 − u j ) ) 1 / ( η + 1 )          e l s e ∆_j=\left\{\begin{matrix}(2u_j)^{1/(η+1)}-1\;\;\;\;\;\;\;\;\;\;\;\;u_j<0.5 \\ (1-(2(1-u_j))^{1/(η+1)}\;\;\;\;else \end{matrix}\right. j={(2uj)1/(η+1)1uj<0.5(1(2(1uj))1/(η+1)else
    并且0≤ u j u_j uj≤1。
    matlab代码如下:

    function chromo_offspring = cross_mutation( chromo_parent,f_num,x_num,x_min,x_max,pc,pm,yita1,yita2,fun )
    %模拟二进制交叉与多项式变异
    [pop,~]=size(chromo_parent);
    suoyin=1;
    for i=1:pop
       %%%模拟二进制交叉
       %初始化子代种群
       %随机选取两个父代个体
       parent_1=round(pop*rand(1));
       if (parent_1<1)
           parent_1=1;
       end
       parent_2=round(pop*rand(1));
       if (parent_2<1)
           parent_2=1;
       end
       %确定两个父代个体不是同一个
       while isequal(chromo_parent(parent_1,:),chromo_parent(parent_2,:))
           parent_2=round(pop*rand(1));
           if(parent_2<1)
               parent_2=1;
           end
       end
       chromo_parent_1=chromo_parent(parent_1,:);
       chromo_parent_2=chromo_parent(parent_2,:);
       off_1=chromo_parent_1;
       off_2=chromo_parent_1;
       if(rand(1)<pc)
           %进行模拟二进制交叉
           u1=zeros(1,x_num);
           gama=zeros(1,x_num);
           for j=1:x_num
               u1(j)=rand(1);
               if u1(j)<0.5
                   gama(j)=(2*u1(j))^(1/(yita1+1));
               else
                   gama(j)=(1/(2*(1-u1(j))))^(1/(yita1+1));
               end
               off_1(j)=0.5*((1+gama(j))*chromo_parent_1(j)+(1-gama(j))*chromo_parent_2(j));
               off_2(j)=0.5*((1-gama(j))*chromo_parent_1(j)+(1+gama(j))*chromo_parent_2(j));
               %使子代在定义域内
               if(off_1(j)>x_max(j))
                   off_1(j)=x_max(j);
               elseif(off_1(j)<x_min(j))
                   off_1(j)=x_min(j);
               end
               if(off_2(j)>x_max(j))
                   off_2(j)=x_max(j);
               elseif(off_2(j)<x_min(j))
                   off_2(j)=x_min(j);
               end
           end
           %计算子代个体的目标函数值
           off_1(1,(x_num+1):(x_num+f_num))=object_fun(off_1,f_num,x_num,fun);
           off_2(1,(x_num+1):(x_num+f_num))=object_fun(off_2,f_num,x_num,fun);
       end
       %%%多项式变异
       if(rand(1)<pm)
           u2=zeros(1,x_num);
           delta=zeros(1,x_num);
           for j=1:x_num
               u2(j)=rand(1);
               if(u2(j)<0.5)
                   delta(j)=(2*u2(j))^(1/(yita2+1))-1;
               else
                   delta(j)=1-(2*(1-u2(j)))^(1/(yita2+1));
               end
               off_1(j)=off_1(j)+delta(j);
               %使子代在定义域内
               if(off_1(j)>x_max(j))
                   off_1(j)=x_max(j);
               elseif(off_1(j)<x_min(j))
                   off_1(j)=x_min(j);
               end
           end
           %计算子代个体的目标函数值
           off_1(1,(x_num+1):(x_num+f_num))=object_fun(off_1,f_num,x_num,fun);
       end
       if(rand(1)<pm)
           u2=zeros(1,x_num);
           delta=zeros(1,x_num);
           for j=1:x_num
               u2(j)=rand(1);
               if(u2(j)<0.5)
                   delta(j)=(2*u2(j))^(1/(yita2+1))-1;
               else
                   delta(j)=1-(2*(1-u2(j)))^(1/(yita2+1));
               end
               off_2(j)=off_2(j)+delta(j);
               %使子代在定义域内
               if(off_2(j)>x_max(j))
                   off_2(j)=x_max(j);
               elseif(off_2(j)<x_min(j))
                   off_2(j)=x_min(j);
               end
           end
           %计算子代个体的目标函数值
           off_2(1,(x_num+1):(x_num+f_num))=object_fun(off_2,f_num,x_num,fun);
       end
       off(suoyin,:)=off_1;
       off(suoyin+1,:)=off_2;
       suoyin=suoyin+2;
    end
    chromo_offspring=off;
    end
    

    1.7 竞标赛选择(tournament selection)

    锦标赛法是选择操作的一种常用方法,二进制竞标赛用的最多。
    假设种群规模为N,该法的步骤为:
    1.从这N个个体中随机选择k(k<n)个个体,k的取值小,效率就高(节省运行时间),但不宜太小,一般取为n/2(取整)。(二进制即取2)
    2.根据每个个体的适应度值,选择其中适应度值最好的个体进入下一代种群。
    3.重复1-2步,至得到新的N个个体。
    二进制选择的matlab代码如下:

    function chromo_parent = tournament_selection( chromo )
    %二进制竞标赛
    [pop, suoyin] = size(chromo);
    touranment=2;
    a=round(pop/2);
    chromo_candidate=zeros(touranment,1);
    chromo_rank=zeros(touranment,1);
    chromo_distance=zeros(touranment,1);
    chromo_parent=zeros(a,suoyin);
    % 获取等级的索引
    rank = suoyin - 1;
    % 获得拥挤度的索引
    distance = suoyin;
    for i=1:a
      for j=1:touranment
          chromo_candidate(j)=round(pop*rand(1));%随机产生候选者
          if(chromo_candidate(j)==0)%索引从1开始
              chromo_candidate(j)=1;
          end
          if(j>1)
              while (~isempty(find(chromo_candidate(1:j-1)==chromo_candidate(j))))
                  chromo_candidate(j)=round(pop*rand(1));
                  if(chromo_candidate(j)==0)%索引从1开始
                      chromo_candidate(j)=1;
                  end
              end
          end
      end
      for j=1:touranment
          chromo_rank(j)=chromo(chromo_candidate(j),rank);
          chromo_distance(j)=chromo(chromo_candidate(j),distance);
      end
      %取出低等级的个体索引
      minchromo_candidate=find(chromo_rank==min(chromo_rank));
      %多个索引按拥挤度排序
      if (length(minchromo_candidate)~=1)
          maxchromo_candidate=find(chromo_distance(minchromo_candidate)==max(chromo_distance(minchromo_candidate)));
          if(length(maxchromo_candidate)~=1)
              maxchromo_candidate = maxchromo_candidate(1);
          end
          chromo_parent(i,:)=chromo(chromo_candidate(minchromo_candidate(maxchromo_candidate)),:);
      else
          chromo_parent(i,:)=chromo(chromo_candidate(minchromo_candidate(1)),:);
      end
    end
    end
    

    2. 算法实现框图

    NSGA-Ⅱ算法的基本思想为:首先,随机产生规模为N的初始种群,非支配排序后通过遗传算法的选择、交叉、变异三个基本操作得到第一代子代种群;其次,从第二代开始,将父代种群与子代种群合并,进行快速非支配排序,同时对每个非支配层中的个体进行拥挤度计算,根据非支配关系以及个体的拥挤度选取合适的个体组成新的父代种群;最后,通过遗传算法的基本操作产生新的子代种群:依此类推,直到满足程序结束的条件。该算法的流程图如图3所示。
     NSGA-Ⅱ算法流程图
    初始化matlab代码如下:

    function chromo = initialize( pop,f_num,x_num,x_min,x_max,fun )
    %   种群初始化
    for i=1:pop
       for j=1:x_num
           chromo(i,j)=x_min(j)+(x_max(j)-x_min(j))*rand(1);
       end
       chromo(i,x_num+1:x_num+f_num) = object_fun(chromo(i,:),f_num,x_num,fun);
    end
    

    主程序matlab代码如下:

    %---------------------------------------------------------------------
    %程序功能:实现nsga2算法,测试函数为ZDT1,ZDT2,ZDT3,ZDT4,ZDT6,DTLZ1,DTLZ2
    %说明:遗传算子为二进制竞赛选择,模拟二进制交叉和多项式变异
    %      需要设置的参数:pop,gen,pc,pm,yita1,yita2
    %作者:(晓风)
    %最初建立时间:2018.9.3
    %最近修改时间:2018.9.20
    %----------------------------------------------------------
    clear all
    clc
    tic;
    %参数设置
    fun='DTLZ1';
    funfun;%函数选择
    pop=300;%种群大小100
    gen=250;%250进化代数
    pc=1;%交叉概率
    pm=1/x_num;%变异概率
    % yita1=1;%模拟二进制交叉参数
    % yita2=10;%多项式变异参数
    yita1=2;%模拟二进制交叉参数10
    yita2=5;%多项式变异参数50
    
    %%初始化种群
    chromo=initialize(pop,f_num,x_num,x_min,x_max,fun);
    %%初始种群的非支配排序
    [F1,chromo_non]=non_domination_sort(pop,chromo,f_num,x_num);%F为pareto等级为pareto_rank的集合,%p为每个个体p的集合(包括每个个体p的被支配个数n和该个体支配的解的集合s),chromo最后一列加入个体的等级
    %%计算拥挤度进行排序
    chromo=crowding_distance_sort(F1,chromo_non,f_num,x_num);
    i=1;
    %%循环开始
    for i=1:gen
       %%二进制竞赛选择(k取了pop/2,所以选两次)
       chromo_parent_1 = tournament_selection(chromo);
       chromo_parent_2 = tournament_selection(chromo);
       chromo_parent=[chromo_parent_1;chromo_parent_2];
       %%模拟二进制交叉与多项式变异
       chromo_offspring=cross_mutation(chromo_parent,f_num,x_num,x_min,x_max,pc,pm,yita1,yita2,fun );
       %%精英保留策略
       %将父代和子代合并
       [pop_parent,~]=size(chromo);
       [pop_offspring,~]=size(chromo_offspring);
       combine_chromo(1:pop_parent,1:(f_num+x_num))=chromo(:,1:(f_num+x_num));
       combine_chromo((pop_parent+1):(pop_parent+pop_offspring),1:(f_num+x_num))=chromo_offspring(:,1:(f_num+x_num));
       %快速非支配排序
       [pop2,~]=size(combine_chromo);
       [F2,combine_chromo1]=non_domination_sort(pop2,combine_chromo,f_num,x_num);
       %计算拥挤度进行排序
       combine_chromo2=crowding_distance_sort(F2,combine_chromo1,f_num,x_num);
       %精英保留产生下一代种群
       chromo=elitism(pop,combine_chromo2,f_num,x_num);
       if mod(i,10) == 0
           fprintf('%d gen has completed!\n',i);
       end
    end
    toc;
    aaa=toc;
    
    hold on
    if(f_num==2)
       plot(chromo(:,x_num+1),chromo(:,x_num+2),'r*');
    end
    if(f_num==3)
       plot3(chromo(:,x_num+1),chromo(:,x_num+2),chromo(:,x_num+2),'r*');
    end
    %%--------------------delta度量--------------------------------
    % [~,index_f1]=sort(chromo(:,x_num+1));
    % d=zeros(length(chromo)-1,1);
    % for i=1:(length(chromo)-1)
    %     d(i)=((chromo(index_f1(i),x_num+1)-chromo(index_f1(i+1),x_num+1))^2+(chromo(index_f1(i),x_num+2)-chromo(index_f1(i+1),x_num+2))^2)^0.5;
    % end
    % d_mean1=mean(d);
    % d_mean=d_mean1*ones(length(chromo)-1,1);
    % d_sum=sum(abs(d-d_mean));
    % delta=(d(1)+d(length(chromo)-1)+d_sum)/(d(1)+d(length(chromo)-1)+(length(chromo)-1)*d_mean1);
    % disp(delta);
    %mean(a)
    %(std(a))^2
    %--------------------Coverage(C-metric)---------------------
    A=PP;B=chromo(:,(x_num+1):(x_num+f_num));%%%%%%%%%%%%%%%%%%%
    [temp_A,~]=size(A);
    [temp_B,~]=size(B);
    number=0;
    for i=1:temp_B
       nn=0;
       for j=1:temp_A
           less=0;%当前个体的目标函数值小于多少个体的数目
           equal=0;%当前个体的目标函数值等于多少个体的数目
           for k=1:f_num
               if(B(i,k)<A(j,k))
                   less=less+1;
               elseif(B(i,k)==A(j,k))
                   equal=equal+1;
               end
           end
           if(less==0 && equal~=f_num)
               nn=nn+1;%被支配个体数目n+1
           end
       end
       if(nn~=0)
           number=number+1;
       end
    end
    C_AB=number/temp_B;
    disp("C_AB:");
    disp(C_AB);
    %-----Distance from Representatives in the PF(D-metric)-----
    A=chromo(:,(x_num+1):(x_num+f_num));P=PP;%%%%%%与论文公式保持相同,注意A和上述不一样
    [temp_A,~]=size(A);
    [temp_P,~]=size(P);
    min_d=0;
    for v=1:temp_P
       d_va=(A-repmat(P(v,:),temp_A,1)).^2;
       min_d=min_d+min(sqrt(sum(d_va,2)));
    end
    D_AP=(min_d/temp_P);
    disp("D_AP:");
    disp(D_AP);
    filepath=pwd;          
    cd('E:\GA\MOEA D\MOEAD_王超\nsga2对比算子修改\DTLZ1');
    save C_AB4.txt C_AB -ASCII
    save D_AP4.txt D_AP -ASCII
    save solution4.txt chromo -ASCII
    save toc4.txt aaa -ASCII
    cd(filepath);
    

    3. 实验结果

    本文使用实数编码的模拟二进制交叉和多项式变异,交叉概率 p c = 0.9 p_c=0.9 pc=0.9,变异概率 p m = 1 / n p_m=1/n pm=1/n η c = 20 η_c=20 ηc=20 η m = 20 η_m=20 ηm=20。以下为选取的5个非凸非均匀的多目标函数的运行结果如图4到图8所示。
    函数范围的matlab代码如下:

    %测试函数
    %--------------------ZDT1--------------------
    if strcmp(fun,'ZDT1')
      f_num=2;%目标函数个数
      x_num=30;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('ZDT1.txt');%导入正确的前端解
      plot(ZDT1(:,1),ZDT1(:,2),'g*');
      PP=ZDT1;
    end
    %--------------------ZDT2--------------------
    if strcmp(fun,'ZDT2')
      f_num=2;%目标函数个数
      x_num=30;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('ZDT2.txt');%导入正确的前端解
      plot(ZDT2(:,1),ZDT2(:,2),'g*');
      PP=ZDT2;
    end
    %--------------------ZDT3--------------------
    if strcmp(fun,'ZDT3')
      f_num=2;%目标函数个数
      x_num=30;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('ZDT3.txt');%导入正确的前端解
      plot(ZDT3(:,1),ZDT3(:,2),'g*');
      PP=ZDT3;
    end
    %--------------------ZDT4--------------------
    if strcmp(fun,'ZDT4')
      f_num=2;%目标函数个数
      x_num=10;%决策变量个数
      x_min=[0,-5,-5,-5,-5,-5,-5,-5,-5,-5];%决策变量的最小值
      x_max=[1,5,5,5,5,5,5,5,5,5];%决策变量的最大值
      load('ZDT4.txt');%导入正确的前端解
      plot(ZDT4(:,1),ZDT4(:,2),'g*');
      PP=ZDT4;
    end
    %--------------------ZDT6--------------------
    if strcmp(fun,'ZDT6')
      f_num=2;%目标函数个数
      x_num=10;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('ZDT6.txt');%导入正确的前端解
      plot(ZDT6(:,1),ZDT6(:,2),'g*');
      PP=ZDT6;
    end
    %--------------------------------------------
    %--------------------DTLZ1--------------------
    if strcmp(fun,'DTLZ1')
      f_num=3;%目标函数个数
      x_num=10;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('DTLZ1.txt');%导入正确的前端解
    %     plot3(DTLZ1(:,1),DTLZ1(:,2),DTLZ1(:,3),'g*');
      PP=DTLZ1;
    end
    %--------------------------------------------
    %--------------------DTLZ2--------------------
    if strcmp(fun,'DTLZ2')
      f_num=3;%目标函数个数
      x_num=10;%决策变量个数
      x_min=zeros(1,x_num);%决策变量的最小值
      x_max=ones(1,x_num);%决策变量的最大值
      load('DTLZ2.txt');%导入正确的前端解
    %     plot3(DTLZ2(:,1),DTLZ2(:,2),DTLZ2(:,3),'g*');
      PP=DTLZ2;
    end
    %--------------------------------------------
    

    函数的matlab代码如下:

    function f = object_fun( x,f_num,x_num,fun )
    % --------------------ZDT1--------------------
    if strcmp(fun,'ZDT1')
      f=[];
      f(1)=x(1);
      sum=0;
      for i=2:x_num
          sum = sum+x(i);
      end
      g=1+9*(sum/(x_num-1));
      f(2)=g*(1-(f(1)/g)^0.5);
    end
    % --------------------ZDT2--------------------
    if strcmp(fun,'ZDT2')
      f=[];
      f(1)=x(1);
      sum=0;
      for i=2:x_num
          sum = sum+x(i);
      end
      g=1+9*(sum/(x_num-1));
      f(2)=g*(1-(f(1)/g)^2);
    end
    % --------------------ZDT3--------------------
    if strcmp(fun,'ZDT3')
      f=[];
      f(1)=x(1);
      sum=0;
      for i=2:x_num
          sum = sum+x(i);
      end
      g=1+9*(sum/(x_num-1));
      f(2)=g*(1-(f(1)/g)^0.5-(f(1)/g)*sin(10*pi*f(1)));
    end
    % --------------------ZDT4--------------------
    if strcmp(fun,'ZDT4')
      f=[];
      f(1)=x(1);
      sum=0;
      for i=2:x_num
          sum = sum+(x(i)^2-10*cos(4*pi*x(i)));
      end
      g=1+9*10+sum;
      f(2)=g*(1-(f(1)/g)^0.5);
    end
    % --------------------ZDT6--------------------
    if strcmp(fun,'ZDT6')
      f=[];
      f(1)=1-(exp(-4*x(1)))*((sin(6*pi*x(1)))^6);
      sum=0;
      for i=2:x_num
          sum = sum+x(i);
      end
      g=1+9*((sum/(x_num-1))^0.25);
      f(2)=g*(1-(f(1)/g)^2);
    end
    % --------------------------------------------
    % --------------------DTLZ1--------------------
    if strcmp(fun,'DTLZ1')
      f=[];
      sum=0;
      for i=3:x_num
          sum = sum+((x(i)-0.5)^2-cos(20*pi*(x(i)-0.5)));
      end
      g=100*(x_num-2)+100*sum;
      f(1)=(1+g)*x(1)*x(2);
      f(2)=(1+g)*x(1)*(1-x(2));
      f(3)=(1+g)*(1-x(1));
    end
    % --------------------------------------------
    % --------------------DTLZ2--------------------
    if strcmp(fun,'DTLZ2')
      f=[];
      sum=0;
      for i=3:x_num
          sum = sum+(x(i))^2;
      end
      g=sum;
      f(1)=(1+g)*cos(x(1)*pi*0.5)*cos(x(2)*pi*0.5);
      f(2)=(1+g)*cos(x(1)*pi*0.5)*sin(x(2)*pi*0.5);
      f(3)=(1+g)*sin(x(1)*pi*0.5);
    end
    % --------------------------------------------
    end
    

    3.1 ZDT1

    f 1 = x 1 f_1=x_1 f1=x1
    g = 1 + 9 ( ( ∑ i = 2 n x i ) / ( n − 1 ) ) g=1+9((∑_{i=2}^nx_i )/(n-1)) g=1+9((i=2nxi)/(n1))
    f 2 = g ( 1 − ( f 1 / g ) 0.5 ) . n = 30 , x i ∈ [ 0 , 1 ] f_2=g(1-(f_1/g)^{0.5} ).n=30,x_i∈[0,1] f2=g(1(f1/g)0.5).n=30,xi[0,1]
    选取种群大小为500,迭代次数500次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图4(上)所示。
    选取种群大小为100,迭代次数2000次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图4(下)所示。
    在这里插入图片描述
    在这里插入图片描述
    图4 ZDT1 pareto最优解对比图(绿色为理论值,红色为实验值)

    3.2 ZDT2

    f 1 = x 1 f_1=x_1 f1=x1
    g = 1 + 9 ( ( ∑ i = 2 n x i ) / ( n − 1 ) ) g=1+9((∑_{i=2}^n x_i )/(n-1)) g=1+9((i=2nxi)/(n1))
    f 2 = g ( 1 − ( f 1 / g ) 2 ) . n = 30 , x i ∈ [ 0 , 1 ] f_2=g(1-(f_1/g)^2 ).n=30,x_i∈[0,1] f2=g(1(f1/g)2).n=30,xi[0,1]
    选取种群大小为500,迭代次数500次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图5(上)所示。
    选取种群大小为100,迭代次数2000次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图5(下)所示。
    在这里插入图片描述
    在这里插入图片描述
    图5 ZDT2 pareto最优解对比图(绿色为理论值,红色为实验值)

    3.3 ZDT3

    f 1 = x 1 f_1=x_1 f1=x1
    g = 1 + 9 ( ( ∑ i = 2 n x i ) / ( n − 1 ) ) g=1+9((∑_{i=2}^nx_i )/(n-1)) g=1+9((i=2nxi)/(n1))
    f 2 = g ( 1 − ( f 1 / g ) 0.5 − ( f 1 / g ) s i n ⁡ ( 10 π f 1 ) ) . n = 30 , x i ∈ [ 0 , 1 ] f_2=g(1-(f_1/g)^{0.5}-(f_1/g)sin⁡(10πf_1)).n=30,x_i∈[0,1] f2=g(1(f1/g)0.5(f1/g)sin(10πf1)).n=30,xi[0,1]
    选取种群大小为136,迭代次数500次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图6示。
    选取种群大小为136,迭代次数2000次,pc=0.9,pm=1/30,ηc=20,ηm=20,得到结果如图6示。
    在这里插入图片描述
    在这里插入图片描述
    图6 ZDT3 pareto最优解对比图(绿色为理论值,红色为实验值)

    3.4 ZDT4

    f 1 = x 1 f_1=x_1 f1=x1
    g = 1 + 9 ∗ 10 + ∑ i = 2 n ( ( x i ) 2 − 10 c o s ⁡ ( 4 π x i ) ) g=1+9*10+∑_{i=2}^n( (x_i )^2-10cos⁡(4πx_i)) g=1+910+i=2n((xi)210cos(4πxi))
    f 2 = g ( 1 − ( f 1 / g ) 0.5 ) . n = 10 , x 1 ∈ [ 0 , 1 ] ; x ( 2 … 6 ) ∈ [ − 5 , 5 ] f_2=g(1-(f_1/g)^{0.5} ).n=10,x_1∈[0,1];x_(2…6)∈[-5,5] f2=g(1(f1/g)0.5).n=10,x1[0,1];x(26)[5,5]
    选取种群大小为100,迭代次数500次,pc=0.9,pm=0.1,ηc=20,ηm=10,得到结果如图7(上)所示。
    选取种群大小为100,迭代次数2000次,pc=0.9,pm=0.1,ηc=20,ηm=10,得到结果如图7(下)所示。
    在这里插入图片描述
    在这里插入图片描述
    图7 ZDT4 pareto最优解对比图(绿色为理论值,红色为实验值)

    3.5 ZDT6

    f 1 = 1 − e − 4 x 1 s i n 6 ( 6 π x 1 ) f_1=1-e^{-4x_1} sin^6 (6πx_1) f1=1e4x1sin6(6πx1)
    g = 1 + 9 ( ( ∑ i = 2 n x i ) / ( n − 1 ) ) 0.25 g=1+9((∑_{i=2}^nx_i )/(n-1))^{0.25} g=1+9((i=2nxi)/(n1))0.25
    f 2 = g ( 1 − ( f 1 / g ) 2 ) . n = 10 , x i ∈ [ 0 , 1 ] f_2=g(1-(f_1/g)^2 ).n=10,x_i∈[0,1] f2=g(1(f1/g)2).n=10,xi[0,1]
    选取种群大小为100,迭代次数500次,pc=0.9,pm=0.1,ηc=20,ηm=20,得到结果如图8(上)所示。
    选取种群大小为100,迭代次数2000次,pc=0.9,pm=0.1,ηc=20,ηm=20,得到结果如图8(下)所示。
    在这里插入图片描述
    在这里插入图片描述
    图8 ZDT6 pareto最优解对比图(绿色为理论值,红色为实验值)
    从结果中可以看出,除ZDT4以外,找到的解几乎全部是pareto前端上的点,并且解在目标空间上分布十分均匀,该算法对于非凸非均匀的多目标函数最优解的寻找还是十分有效的。由于ZDT4具有许多局部pareto最优前沿,为了摆脱这些局部前沿,需要对决策向量进行大的改变,因此选取ηm=10,但仍离真正的pareto前沿还有很大的差距。

    3.6 对比

    表1 γ的均值和方差(500代)

    项目ZDT1ZDT2ZDT3ZDT4ZDT6
    均值0.05320.09720.171220.23040.3284
    方差0.00210.00510.00480.63100.0109

    表2 Δ的均值和方差(500代)

    项目ZDT1ZDT2ZDT3ZDT4ZDT6
    均值0.45750.46150.66400.52030.9538
    方差0.00210.00150.00300.00130.0055

    表3 γ的均值和方差(2000代)

    项目ZDT1ZDT2ZDT3ZDT4ZDT6
    均值0.00110.00270.022512.03410.2092
    方差0.00020.00020.00050.33200.0021

    表4 Δ的均值和方差(2000代)

    项目ZDT1ZDT2ZDT3ZDT4ZDT6
    均值0.44780.43370.65860.51680.4529
    方差0.00010.00030.00060.00180.0001

    根据上述实验数据,进行五组实验,计算distance metric γ,diversity metric Δ的均值和方差,迭代500次后得到结果如表格1-2所示。迭代2000次后得到的结果如表格3-4所示。
          表1显示了NSGA-Ⅱ(实数编码500代)在5个问题上运行获得的distance metric γ的均值和方差,实验结果表明,除了ZDT4和ZDT6以外的问题均取得了良好的收敛效果,除ZDT4以外,其他问题5次运行结果的差异也都很小。
          表2显示了NSGA-Ⅱ(实数编码500代)在5个问题上运行获得的diversity metric Δ的均值和方差,实验结果表明,NSGA-Ⅱ在上述问题上均具有较强的传播多样性,并且5次运行结果的差异都很小。
          像500代的时候一样保留所有其他参数,但将代数增加至2000代,得到表3的distance metric γ的均值和方差以及表4的diversity metric Δ的均值和方差。将表3表4与表1表2对比可知,5个问题的distance metric γ随着代数的增加都有改善,其中ZDT3和ZDT4的改善尤为显著,除此之外,因为在500代时问题的结果已经有了良好的多样性分布了,所以代数的增加并没有对diversity metric Δ产生巨大的影响。
    [1]:Deb K, Pratap A, Agarwal S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2):182-197.

    展开全文
  • 【MATLAB】多目标优化算法 NSGA-II (gamultiobj) 的使用

    万次阅读 多人点赞 2020-08-27 09:57:09
    【MATLAB】多目标优化 简单写一写 一个例子 建立模型如下 min⁡F(X)=(−Z(x),f(x)) \min \text{F}\left( \text{X} \right) =\left( -\text{Z}\left( \text{x} \right) ,\text{f}\left( \text{x} \right) \right) ...

    程序输出

    【MATLAB】多目标优化算法NSGA-II(gamultiobj)的使用精解

    原始博文因为写的比较潦草,评论中有疑问的网友较多,所以重新写了一下 2021-4-24

    增加了一些说明与参考文献,修改了几处笔误 2021-5-20

    对于多目标优化(multiobjective optimization)算法NSGA-II实现的细节与原理不在此说明。感兴趣的读者可另行查阅

    gamultiobj的使用范式

    编写程序

    清除所有变量(非必须,但注意变量不能和下面所用的冲突)

    clear
    
    • 需求解模型的参数设置部分:(模型导入)
    %% 模型设置
    % 适应度函数的函数句柄
    fitnessfcn=@Fun;
    % 变量个数
    nvars=4;
    % 约束条件形式1:下限与上限(若无取空数组[])
    % lb<= X <= ub
    lb=[0,0,0,0];
    ub=[];
    
    % 约束条件形式2:线性不等式约束(若无取空数组[])
    % A*X <= b 
    A = [0    0 1 1
        -1/3  0 0 0
         0 -1/2 0 0
         0    0 0 0];
    
    b = [48 ; 30 ; 30 ; 0];
    
    % 约束条件形式3:线性等式约束(若无取空数组[])
    % Aeq*X == beq
    Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
    beq=[120;0;0;0];
    

    目标函数: (这一段需放在脚本最后或单独放在一个文件里)

    function y=Fun(x)
    	% y是目标函数向量。有几个目标函数y就有多少个维度(数组y的长度)
    	% 因为gamultiobj是以目标函数分量取极小值为目标,
    	% 因此有些取极大值的目标函数注意取相反数
    	y(1)=-(x(1)*100/3 + x(3)*90/3  + x(2)*80/2+x(4)*70/2);
    	y(2)=x(3)+x(4);
    end
    
    • gamultiobj求解器设置部分:
    %% 求解器设置
    % 最优个体系数paretoFraction
    % 种群大小populationsize
    % 最大进化代数generations
    % 停止代数stallGenLimit
    % 适应度函数偏差TolFun
    % 函数gaplotpareto:绘制Pareto前沿 
    options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
    
    • gamultiobj求解与结果输出部分:
    %% 主求解
    [x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
    
    %% 结果提取
    % 因为gamultiobj是以目标函数分量取极小值为目标,
    % 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
    plot(-fval(:,1),fval(:,2),'pr')
    xlabel('f_1(x)')
    ylabel('f_2(x)')
    title('Pareto front')
    grid on
    
    

    RUN

    求解时间受求解器设置影响,可能会较长,请耐心等待

    求解过程中会实时显示当前种群的情况:

    程序输出

    如果已经达到满意,也可点击stop按钮提前结束求解

    最后的求解结果,即Pareto最优解集储存在[x,fval]中,fvalx对应的目标函数值。fval大致构成了一条空间曲线——Pareto前沿。若各个解较为均匀分布,说明该图包含了大部分最优解情况,全局性优,适用性强。在满足Pareto最优的条件下,是没有办法在不让某一优化目标受损的情况下,令另一方目标获得更优的。所以这些解均为最优,对最优解的具体选择可以根据实际情况。

    程序输出

    程序输出

    例子1

    Image

    表1 工厂产品生产规格表

    产品生产时间(h/公斤)利润(元/公斤)加班时利润(元/公斤)
    A310090
    B28070

    设工厂每周生产产品A、B的常规生产时长为 x 1 x_1 x1 x 2 x_2 x2(h),加班生产时长为 x 3 x_3 x3 x 4 x_4 x4 (h)。令 x = ( x 1 , x 2 , x 3 , x 4 ) x=\left( x_1,x_2,x_3,x_4 \right) x=(x1,x2,x3,x4) 。设每周的利润函数为 Z ( x ) Z(x) Z(x),加班时长函数为 f ( x ) f(x) f(x)

    则目标函数为:
    Z ( x ) = x 1 3 × 100 + x 3 3 × 90 + x 2 2 × 80 + x 4 2 × 70 Z\left( x \right) =\frac{x_1}{3}\times 100+\frac{x_3}{3}\times 90+\frac{x_2}{2}\times 80+\frac{x_4}{2}\times 70 Z(x)=3x1×100+3x3×90+2x2×80+2x4×70

    f ( x ) = x 3 + x 4 f\left( x \right) =x_3+x_4 f(x)=x3+x4

    约束条件为:
    s . t . { x 1 + x 2 = 120 x 3 + x 4 ⩽ 48 x 1 3 ⩾ 30 x 2 2 ⩾ 30 x 1 , x 2 , x 3 , x 4 ⩾ 0 \mathrm{s}.\mathrm{t}.\left\{ \begin{array}{c} \begin{array}{c} \mathrm{x}_1+\mathrm{x}_2=120\\ \mathrm{x}_3+\mathrm{x}_4\leqslant 48\\ \end{array}\\ \frac{\mathrm{x}_1}{3}\geqslant 30\\ \frac{\mathrm{x}_2}{2}\geqslant 30\\ \mathrm{x}_1,\mathrm{x}_2,\mathrm{x}_3,\mathrm{x}_4\geqslant 0\\ \end{array} \right. s.t.x1+x2=120x3+x4483x1302x230x1,x2,x3,x40

    因此,数学模型可以归纳为:

    m i n    F ( X ) = ( − Z ( x ) , f ( x ) ) min\,\,F\left( X \right) =\left( -Z\left( x \right) ,f\left( x \right) \right) minF(X)=(Z(x),f(x))

    s . t . { x 1 + x 2 = 120 x 3 + x 4 ⩽ 48 x 1 3 ⩾ 30 x 2 2 ⩾ 30 x 1 , x 2 , x 3 , x 4 ⩾ 0 \mathrm{s}.\mathrm{t}.\left\{ \begin{array}{c} \begin{array}{c} \mathrm{x}_1+\mathrm{x}_2=120\\ \mathrm{x}_3+\mathrm{x}_4\leqslant 48\\ \end{array}\\ \frac{\mathrm{x}_1}{3}\geqslant 30\\ \frac{\mathrm{x}_2}{2}\geqslant 30\\ \mathrm{x}_1,\mathrm{x}_2,\mathrm{x}_3,\mathrm{x}_4\geqslant 0\\ \end{array} \right. s.t.x1+x2=120x3+x4483x1302x230x1,x2,x3,x40

    MATLAB求解如下:

    clear
    clc
    fitnessfcn=@Fun;
    % 变量个数
    nvars=4;
    % lb<= X <= ub
    lb=[0,0,0,0];
    ub=[270 240 460 130];
    % A*X <= b 
    A = [[13 13.5 14 11.5]
        -[13 13.5 14 11.5]
          0 0 -1 0
         [0.015 0.02 0.018 0.011]];
    
    b = [300*48 ; -300*40 ; -150 ; 20];
    
    % Aeq*X = beq
    Aeq=[];beq=[];
    %最优个体系数paretoFraction
    %种群大小populationsize
    %最大进化代数generations
    %停止代数stallGenLimit
    %适应度函数偏差TolFun
    %函数gaplotpareto:绘制Pareto前沿 
    options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
    
    [x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
    
    plot(-fval(:,1),fval(:,2),'pr')
    xlabel('f_1(x)')
    ylabel('f_2(x)')
    title('Pareto front')
    grid on
    
    
    function y=Fun(x)
    	b = [270 240 460 130];
    	c = [300 300 600 200];
    	t = [190 210 148 100];
    	s = [200 230 160 114];
    	a = [0.015 0.02 0.018 0.011];
    	d = [13 13.5 14 11.5];
    	y(1)=-sum(x.*(s-t));
    	y(2)=sum(a.*x);
    end
    

    得到结果:

    Image

    x1 x2 x3 x4如下:

    119.231391258967	0.769608488165712	0	0
    
    119.231391258967	0.769608488165712	0	0
    
    0.000499510291359209	120.000482867813	13.1757491344817	34.8252467588411
    
    71.3391090591218	48.6614868846125	3.11344686170493	9.36001224382937
    
    …… …… …… ……
    
    27.0549008917871	92.9459248170927	9.47711506311976	25.3969178809142
    
    8.65187243477257	111.349067997015	13.3680683558073	31.7052095195761
    

    例子2

    Image

    i = 1 , 2 , 3 , 4 \mathrm{i}=1,2,3,4 i=1,2,3,4分别表示A、B、C、D四种产品, x i x_i xi表示第i种产品的产量(kg)。设最大产量为 b i b_i bi,销售量为 c i c_i ci,成本为 t i t_i ti,售价为 s i s_i si,能耗为 a i a_i ai,生产时间为 d i d_i di。设该问题的利润函数为 Z ( x ) Z\left( x \right) Z(x),能耗函数为 f ( x ) f\left( x \right) f(x)
    则利润函数为:

    Z ( X ) = ∑ i = 1 4 x i ( s i − t i ) Z\left( X \right) =\sum_{i=1}^4{x_i\left( s_i-t_i \right)} Z(X)=i=14xi(siti)

    能耗函数为:

    f ( x ) = ∑ i = 1 4 a i x i f\left( x \right) =\sum_{i=1}^4{a_ix_i} f(x)=i=14aixi

    • 建立模型如下

    m i n F ( X ) = ( − Z ( x ) , f ( x ) ) minF\left( X \right) =\left( -Z\left( x \right) ,f\left( x \right) \right) minF(X)=(Z(x),f(x))

    s . t . { x i ⩽ b i ( i = 1 , 2 , 3 , 4 ) 300 × 40 ⩽ ∑ i = 1 4 x i d i ⩽ 300 × 48 150 ⩽ x 3 ∑ i = 1 4 a i x i ⩽ 20 x i ⩾ 0 ( i = 1 , 2 , 3 , 4 ) \mathrm{s}.\mathrm{t}.\left\{ \begin{array}{c} \begin{array}{c} x_i\leqslant b_i\left( i=1,2,3,4 \right)\\ 300\times 40\leqslant \sum_{i=1}^4{x_id_i\leqslant 300\times 48}\\ 150\leqslant x_3\\ \end{array}\\ \sum_{i=1}^4{a_ix_i}\leqslant 20\\ x_i\geqslant 0\left( i=1,2,3,4 \right)\\ \end{array} \right. s.t.xibi(i=1,2,3,4)300×40i=14xidi300×48150x3i=14aixi20xi0(i=1,2,3,4)

    • MATLAB求解如下
    % 清除所有变量(非必须)
    clear
    
    %% 模型设置
    % 获取目标函数的函数句柄
    fitnessfcn=@Fun;
    % 变量个数
    nvars=4;
    % 约束条件形式1:(若无取空数组[])
    % lb<= X <= ub
    lb=[0,0,0,0];
    ub=[];
    
    % 约束条件形式2:(若无取空数组[])
    % A*X <= b 
    A = [0    0 1 1
        -1/3  0 0 0
         0 -1/2 0 0
         0    0 0 0];
    
    b = [48 ; 30 ; 30 ; 0];
    
    % 约束条件形式3:(若无取空数组[])
    % Aeq*X = beq
    Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
    beq=[120;0;0;0];
    
    %% 求解器设置
    % 最优个体系数paretoFraction
    % 种群大小populationsize
    % 最大进化代数generations
    % 停止代数stallGenLimit
    % 适应度函数偏差TolFun
    % 函数gaplotpareto:绘制Pareto前沿 
    options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
    
    %% 主求解
    [x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
    
    %% 结果提取
    % 因为gamultiobj是以目标函数分量取极小值为目标,
    % 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
    plot(-fval(:,1),fval(:,2),'pr')
    xlabel('f_1(x)')
    ylabel('f_2(x)')
    title('Pareto front')
    grid on
    
    
    function y=Fun(x)
    % y是目标函数向量。有几个目标函数y就有多少个维度(数组y的长度)
    % 因为gamultiobj是以目标函数分量取极小值为目标,
    % 因此有些取极大值的目标函数注意取相反数
    y(1)=-(x(1)*100/3 + x(3)*90/3  + x(2)*80/2+x(4)*70/2);
    y(2)=x(3)+x(4);
    end
    

    求解结果为:

    Image

    x1 x2 x3 x4如下:

    257.911499184609	147.920309053797	368.392357989384	129.803238959239
    
    204.527452370415	215.831670926692	376.135110383218	129.541339137201
    
    251.942563886570	239.988410149935	456.713154231118	129.635721179650
    
    …… …… …… ……
    
    245.261897051381	238.784443203755	429.044675830294	129.548264747046
    
    216.531507460989	201.737471951873	368.856283121162	129.726418090506
    

    参考文献:

    K. Deb, S. Agrawal, A. Pratap, T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 2002 6(2): 182-197.

    展开全文
  • python拓展包之pymoo使用方法:多目标优化一、pymoo的安装二、多目标优化的一般模式三、pymoo处理多目标优化问题的格式python中pymoo的使用步骤一级目录 一、pymoo的安装 pip安装 pip install -U pymoo 二、多目标...

    一、pymoo的安装

    pip安装

    pip install -U pymoo
    

    二、多目标优化的一般模式

    一般来说,多目标优化具有几个受不等式和等式约束的目标函数。其目标是找到一组满足所有约束条件并尽可能好地处理其所有目标值的解决方案。问题的一般形式的定义为:
    在这里插入图片描述
    目标函数f(x)
    不等式约束g(x)
    等式约束和h(x)
    变量x的上下约束

    在下面,我们说明了一个具有两个约束条件的双目标优化问题。

    在这里插入图片描述
    它的可视化图像如下:

    在这里插入图片描述
    可以看到,它的最优解在(图中橙色区域)在这里插入图片描述
    下面我们在python中使用pymoo来实现上述问题公式。

    三、pymoo处理多目标优化问题的格式

    一、pymoo支能够处理目标函数最小化的优化问题

    然而,在不失去一般性的情况下,一个应该最大化的目标可以乘以−1并最小化。
    在上述优化问题中,我们需要将目标函数
    max f2(x)
    转化为
    min -f2(x)
    此外,所有的约束函数都需要表述为一个小于等于的约束。
    我们需要将约束条件
    g2(x)>=0
    转化为
    -g2(x)<=0
    我们建议将约束标准化,使它们都同等重视。将约束条件的每一个因式的常数部分相乘得到一个约束常量,
    对于g1(x),约束常量为 2 ⋅ (−0.1) ⋅(−0.9) = 0.18
    对于g2(x),约束常量为 20 ⋅ (−0.4) ⋅ (−0.6) = 4.8
    最后,将使用pymoo进行优化的优化问题定义为:

    pymoo的目标函数及约束

    四、python中pymoo的使用

    接下来,用Python实现该优化问题。pymoo中的每个优化问题都必须从Problem类中继承。
    首先,通过调用super()函数,可以初始化问题属性。

    from pymoo.model.problem import Problem
    class MyProblem(Problem):
        def __init__(self,cost_matrix: list, ):
            self.cost_matrix = cost_matrix
            super().__init__(n_var=2,   # 变量数
                             n_obj=2,   # 目标数
                             n_constr=2,    # 约束数
                             xl=np.array([-2, -2]),     # 变量下界
                             xu=np.array([2, 2]),   # 变量上界
                             )
    
        def _evaluate(self, x, out, *args, **kwargs):
    
            # 定义目标函数
            f1 = x[:, 0]**2 + x[:, 1]**2    # x1放在x的第0列,x2放在x的第一列
            f2 = (x[:, 0] - 1)**2 + x[:, 1]**2
            # 定义约束条件
            g1 = 2*(x[:, 0] - 0.1) * (x[:, 0] - 0.9) / 0.18
            g2 = -20*(x[:, 0] - 0.4) * (x[:, 0] - 0.6) / 4.8
            # todo
            out["F"] = np.column_stack([f1, f2])
            out["G"] = np.column_stack([g1, g2])
    
    

    super初始化中
    变量数(n_var)
    目标(n_obj)
    约束(n_constr)
    此外,下(xl)和上变量边界(xu)作为NumPy数组提供。


    _evaluate函数中
    定义目标函数以键“F”添加到字典中
    定义约束条件以键“G”添加到字典中

    五、选择优化算法

    在pymoo中,需要创建一个算法对象来进行优化。对于每种算法,都有API文档,通过提供不同的参数,可以以即插即用的方式定制算法。一般来说,选择一个合适的优化问题的算法本身就是一个挑战。当事先知道问题特征时,我们建议通过定制操作符使用这些特征。然而,在我们的情况下,优化问题相当简单,但应该考虑有两个目标和两个约束的方面。因此,我们决定使用NSGA-II[12]及其默认配置,并进行少量修改。我们选择了40人的人口规模,但我们没有产生相同数量的后代,而是每一代只创造10个后代。这是NSGA-II的一个稳态变体,对于比较简单的优化问题,可以提高收敛性,存在局部帕雷托前沿。此外,我们启用了一个重复的检查,以确保交配产生的后代,以及与现有种群关于其可变向量的不同。
    使用所提供的参数调用NSGA2的构造函数并返回一个初始化的算法对象。

    from pymoo.algorithms.nsga2 import NSGA2
    from pymoo.factory import get_sampling, get_crossover, get_mutation
    from pymoo.optimize import minimize
    from example import MyProblem
    
    # 定义遗传算法
    algorithm = NSGA2(
        pop_size=40,
        n_offsprings=10,
        sampling=get_sampling("real_random"),
        crossover=get_crossover("real_sbx", prob=0.9, eta=15),
        mutation=get_mutation("real_pm", eta=20),
        eliminate_duplicates=True
    )
    

    接下来,我们使用初始化的算法对象来优化所定义的问题。因此,将具有实例问题和算法的最小化函数作为参数。此外,我们提供了运行40代算法的终止标准,这将导致40个+40×10=440函数评估。此外,我们定义了一个随机种子,以确保重现性,并使冗长标志能够看到每一代的打印输出。该方法返回一个结果对象,其中包含该算法找到的非支配解集。

    res = minimize(MyProblem(),
                   algorithm,
                   ('n_gen', 40),
                   seed=1,
                   verbose=True
                   )
    

    打印输出

    D:\Soft\Anaconda3\envs\datadeal\python.exe D:/WorkSpace/pymooto/alg_example.py
    =====================================================================================
    n_gen |  n_eval |   cv (min)   |   cv (avg)   |  n_nds  |     eps      |  indicator  
    =====================================================================================
        1 |      40 |  0.00000E+00 |  2.36399E+01 |       1 |            - |            -
        2 |      50 |  0.00000E+00 |  1.15486E+01 |       1 |  0.00000E+00 |            f
        3 |      60 |  0.00000E+00 |  5.277918607 |       1 |  0.00000E+00 |            f
        4 |      70 |  0.00000E+00 |  2.406068542 |       2 |  1.000000000 |        ideal
        5 |      80 |  0.00000E+00 |  0.908316880 |       3 |  0.869706146 |        ideal
        6 |      90 |  0.00000E+00 |  0.264746300 |       3 |  0.00000E+00 |            f
        7 |     100 |  0.00000E+00 |  0.054063822 |       4 |  0.023775686 |        ideal
        8 |     110 |  0.00000E+00 |  0.003060876 |       5 |  0.127815454 |        ideal
        9 |     120 |  0.00000E+00 |  0.00000E+00 |       6 |  0.004633441 |        ideal
       10 |     130 |  0.00000E+00 |  0.00000E+00 |       6 |  0.062649485 |        nadir
       11 |     140 |  0.00000E+00 |  0.00000E+00 |       7 |  0.026769546 |            f
       12 |     150 |  0.00000E+00 |  0.00000E+00 |       7 |  0.047566009 |            f
       13 |     160 |  0.00000E+00 |  0.00000E+00 |       7 |  0.000678729 |            f
       14 |     170 |  0.00000E+00 |  0.00000E+00 |       9 |  0.043888006 |        ideal
       15 |     180 |  0.00000E+00 |  0.00000E+00 |      10 |  0.013506283 |            f
       16 |     190 |  0.00000E+00 |  0.00000E+00 |      11 |  0.052719639 |        ideal
       17 |     200 |  0.00000E+00 |  0.00000E+00 |      14 |  0.008531768 |            f
       18 |     210 |  0.00000E+00 |  0.00000E+00 |      16 |  0.011165361 |            f
       19 |     220 |  0.00000E+00 |  0.00000E+00 |      18 |  0.004040787 |            f
       20 |     230 |  0.00000E+00 |  0.00000E+00 |      20 |  0.001296614 |            f
       21 |     240 |  0.00000E+00 |  0.00000E+00 |      22 |  0.000911500 |            f
       22 |     250 |  0.00000E+00 |  0.00000E+00 |      24 |  0.003461153 |        nadir
       23 |     260 |  0.00000E+00 |  0.00000E+00 |      27 |  0.004261571 |        ideal
       24 |     270 |  0.00000E+00 |  0.00000E+00 |      28 |  0.029927769 |        nadir
       25 |     280 |  0.00000E+00 |  0.00000E+00 |      28 |  0.001150567 |            f
       26 |     290 |  0.00000E+00 |  0.00000E+00 |      30 |  0.002181257 |            f
       27 |     300 |  0.00000E+00 |  0.00000E+00 |      33 |  0.006887739 |        nadir
       28 |     310 |  0.00000E+00 |  0.00000E+00 |      36 |  0.001237981 |            f
       29 |     320 |  0.00000E+00 |  0.00000E+00 |      37 |  0.000102543 |            f
       30 |     330 |  0.00000E+00 |  0.00000E+00 |      38 |  0.000747546 |            f
       31 |     340 |  0.00000E+00 |  0.00000E+00 |      40 |  0.005157126 |        nadir
       32 |     350 |  0.00000E+00 |  0.00000E+00 |      40 |  0.00000E+00 |            f
       33 |     360 |  0.00000E+00 |  0.00000E+00 |      40 |  0.002321697 |            f
       34 |     370 |  0.00000E+00 |  0.00000E+00 |      40 |  0.016674348 |        nadir
       35 |     380 |  0.00000E+00 |  0.00000E+00 |      40 |  0.000760584 |            f
       36 |     390 |  0.00000E+00 |  0.00000E+00 |      40 |  0.016011922 |        nadir
       37 |     400 |  0.00000E+00 |  0.00000E+00 |      40 |  0.000333129 |            f
       38 |     410 |  0.00000E+00 |  0.00000E+00 |      40 |  0.000980877 |            f
       39 |     420 |  0.00000E+00 |  0.00000E+00 |      40 |  0.001264690 |            f
       40 |     430 |  0.00000E+00 |  0.00000E+00 |      40 |  0.001571863 |            f
    

    关于pymoo的中文使用教程很少,本文内容大都翻译自pymoo的官方使用说明。

    参考资料

    pymoo: Multi-objective Optimization in Python

    展开全文
  • 考虑自平衡能力的并网型微电网多目标容量优化设计nsga2算法matlab实现
  • 遗传算法代码,外加个人理解希望大吉你可以多多交流

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,014
精华内容 405
关键字:

多目标优化问题nsga