差分进化算法解决多目标优化问题--内附matlab代码

Zephus

一 多目标优化

所谓的多目标优化是尝试同时最小化K个独立的目标函数。其目标是：
求





由此可知只有当有一个单一的向量同时取得K个目标函数的最小值，及可以说明其是多目标的解。故此解向量x满足条件：






往往这k个目标是相互冲突的，想要找到一个最优解满足k目标这是比较困难的，所以需要解决这一问题，排除彼此冲突的情况。最优解是一个折中的方案，取决于哪个目标更重要。如果每个目标函数能被指定一个权重来衡量其相对重要性，那么最优解就更加明显。
二 差分进化算法解决多目标问题

在了解多目标之后，现在就利用差分进化算法(DE)实战应用在其上。这里采用基准测试集（ZDT测试函数，来源于文献：Comparison of Multiobjective Evolutionary Algorithms: Empirical Results；SCH,FON,POL,KUR,来源于文献：A Fast and Elitist Multiobjective Genetic Algorithm:
NSGA-II；DTLZ，来源于文献：Scalable Multi-Objective Optimization Test Problems；MOP，来源于文献：Decomposition of a Multiobjective Optimization Problem into a Number of Simple Multiobjective Subproblems）

测试函数文件：

%%%%%%%%%%%%%%%%%%多目标测试函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=muti_fun(x,type)
global D f_num
if type==1  %SCH----维度1---边界范围[-10^3,10^3]
    y(1)=x.^2;  
    y(2)=(x-2).^2;
elseif type==2  %FON-------维度3----边界范围[-4,4]
    L=3;
    sumFON1=0;
    sumFON2=0;
    for i=1:L
        sumFON1=sumFON1+(x(i)-1/sqrt(3)).^2;
        sumFON2=sumFON2+(x(i)+1/sqrt(3)).^2;
    end
    y(1)=1-exp(-sumFON1);
    y(2)=1-exp(-sumFON2);
elseif type==3   %POL-----维度2-----边界范围[-pi,pi]
   A1=0.5*sin(1)-2*cos(1)+sin(2)-1.5*cos(2);
   A2=1.5*sin(1)-cos(1)+sin(2)-0.5*cos(2);
   B1=0.5*sin(x(1))-2*cos(x(1))+sin(x(2))-1.5*cos(x(2));
   B2=1.5*sin(x(1))-cos(x(1))+2*sin(x(2))-0.5*cos(x(2));
   y(1)=1+(A1-B1).^2+(A2-B2).^2;
   y(2)=(x(1)+3).^2+(x(2)+1).^2;
elseif type==4  %KUR-----------维度3-----边界范围[-5,5]
    sumKUR1=0;
    sumKUR2=0;
    for i=1:D-1
        sumKUR1=sumKUR1+(-10*exp(-0.2*sqrt(x(i).^2+x(i+1).^2)));
        sumKUR2=sumKUR2+abs(x(i)).^(0.8)+5*sin(x(i).^3);
    end
    y(1)=sumKUR1;
    y(2)=sumKUR2;
elseif type==5  %ZDT1----维度30----边界范围[0,1]
    y(1)=x(1);
    y(2)=1-sqrt(y(1)/g(x));
elseif type==6  %ZDT2----维度30----边界范围[0,1]
    y(1)=x(1);
    y(2)=1-(y(1)/g1(x)).^2;
elseif type==7  %ZDT3----维度30----边界范围[0,1]
    y(1)=x(1);
    y(2)=1-sqrt(y(1)/g1(x))-(y(1)/g1(x))*sin(10*pi*y(1));
elseif type==8  %ZDT4----维度10----边界范围x1=[0,1],x(2,...,10)=[-5,5]------PF distribution of objection space have problems
    y(1)=x(1);
   % y(2)=g2(x)*(1-(y(1)/g2(x)).^2);
    y(2)=1-sqrt(y(1)/g2(x));
elseif type==9  %ZDT6----维度10----边界范围[0,1]
    y(1)=1-exp(-4*x(1))*sin(6*pi*x(1)).^6;
    y(2)=1-(y(1)/g3(x)).^2;
elseif type==10  %MOP1----维度10----边界范围[0,1]
    y(1)=(1+g4(x))*x(1);
    y(2)=(1+g4(x))*(1-sqrt(x(1)));
elseif type==11 %MOP2----维度10----边界范围[0,1]
    y(1)=(1+g5(x))*x(1);
    y(2)=(1+g5(x))*(1-x(1).^2);
elseif type==12 %MOP3----维度10----边界范围[0,1]
    y(1)=(1+g6(x))*cos((pi*x(1))/2);
    y(2)=(1+g6(x))*sin((pi*x(1))/2);
elseif type==13  %MOP4----维度10----边界范围[0,1]
    y(1)=(1+g7(x))*x(1);
    y(2)=(1+g7(x))*(1-x(1).^2*cos(2*pi*x(1)).^2);
elseif type==14   %MOP5----维度10----边界范围[0,1]
    y(1)=(1+g8(x))*x(1);
    y(2)=(1+g8(x))*(1-sqrt(x(1)));
elseif type==15 %MOP6----三目标----维度10----边界范围[0,1]
    y(1)=(1+g9(x))*x(1)*x(2);
    y(2)=(1+g9(x))*x(1)*(1-x(2));
    y(3)=(1+g9(x))*(1-x(1));
elseif type==16  %MOP7----三目标----维度10----边界范围[0,1]
    y(1)=(1+g10(x))*cos((x(1)*pi)/2)*cos((x(2)*pi)/2);
    y(2)=(1+g10(x))*cos((x(1)*pi)/2)*sin((x(2)*pi)/2);
    y(3)=(1+g10(x))*sin((x(1)*pi)/2);
elseif type==17   %DTLZ1----三目标----维度2----边界范围[0,1]
    for i=1:f_num
        y(1)=(1/2)*x(1)*x(2)*(1+g11(x));
        y(2)=(1/2)*x(1)*(1-x(2))*(1+g11(x));
        y(3)=(1/2)*(1-x(1))*(1+g11(x));
    end
   % y(1)=1/2*x(1)*(1+g11(x));
   % y(2)=1/2*(1-x(1))*(1+g11(x));
elseif type==18    %DTLZ2----三目标----维度2----边界范围[0,1]
    y(1)=(1+g12(x))*cos((pi/2)*x(1))*cos((pi/2)*x(2));
    y(2)=(1+g12(x))*sin((pi/2)*x(1))*sin((pi/2)*x(2));
    y(3)=(1+g12(x))*sin((pi/2)*x(1));
elseif type==19   %DTLZ3----
   y(1)=(1+g11(x))*cos((pi/2)*x(1));
   y(2)=(1+g11(x))*sin((pi/2)*x(1));
elseif type==20    %DTLZ4----
    a=100;
    y(1)=(1+g12(x))*cos((pi/2)*x(1).^a);
    y(2)=(1+g12(x))*sin((pi/2)*x(1).^a);
elseif type==21     %DTLZ7----
    y(1)=x(1);
    y(2)=(1+g13(x))*h(x);
end
end
   
% function f=fun1(x,a,b)
%     if x>a && x<=b
%     f=exp(-2*((x-0.1)/0.8).^2*log(2))*abs(sin(5*pi*x)).^0.5;
%      else
%     f=exp(-2*((x-0.1)/0.8).^2*log(2))*abs(sin(5*pi*x)).^6;
%     end
% end
%
% function f=fun2(x,e)
% f=sin(5*pi*(x.^e-0.05)).^6;
% end

function f=g(x)
global D
f=1+9*(sum(x(2:D))/(D-1));
end

function f=g1(x)
global D

f=1+9*(sum(x(2:D))/(D-1));
end

function f=g2(x)
global D
D=10;
f=0;
for i=2:D
f=f+x(i).^2-10*cos(4*pi*x(i));
end
f=f+1+10*(D-1);
end

function f=g3(x)
global D

f=1+9*(sum(x(:,2:D))/(D-1)).^0.25;
end

function f=g4(x) %MOP1
global D
f=0;
for i=2:D
    t(i)=x(i)-sin(0.5*pi*x(1));
    f=-0.9*t(i).^2+abs(t(i)).^0.6+f;
end
f=f*2*sin(pi*x(1));
end

function f=g5(x) %MOP2
global D
f=0;
for i=2:D
    t(i)=x(i)-sin(0.5*pi*x(1));
    f=(abs(t(i))/(1+exp(5*abs(t(i)))))+f;
end
f=f*10*sin(pi*x(1));
end

function f=g6(x) %MOP3
global D
f=0;
for i=2:D
    t(i)=x(i)-sin(0.5*pi*x(1));
    f=(abs(t(i))/(1+exp(5*abs(t(i)))))+f;
end
f=f*10*sin(pi*x(1)*(1/2));
end

function f=g7(x) %MOP4
global D
f=0;
for i=2:D
    t(i)=x(i)-sin(0.5*pi*x(1));
    f=(abs(t(i))/(1+exp(5*abs(t(i)))))+f;
end
f=f*10*sin(pi*x(1));
end

function f=g8(x) %MOP5
global D
f=0;
for i=2:D
    t(i)=x(i)-sin(0.5*pi*x(1));
    f=-0.9*t(i).^2+abs(t(i)).^0.6+f;
end
f=f*2*abs(cos(pi*x(1)));
end

function f=g9(x) %MOP6
global D
f=0;
for i=2:D
    t(i)=x(i)-x(1)*x(2);
    f=-0.9*t(i).^2+abs(t(i)).^0.6+f;
end
f=f*2*sin(pi*x(1));
end

function f=g10(x) %MOP7
global D
f=0;
for i=3:D
    t(i)=x(i)-x(1)*x(2);
    f=-0.9*t(i).^2+abs(t(i)).^0.6+f;
end
f=f*2*sin(pi*x(1));
end

function f=g11(x)  %DTLZ1 DTLZ3
global D
f=0;
for i=1:D
    f=(x(i)-0.5).^2-cos(20*pi*(x(i)-0.5))+f;
end
f=(f+sqrt(sum(x.^2)))*100;
end

function f=g12(x)   %DTLZ2 DTLZ4
global D
f=0;
for i=1:D
    f=f+(x(i)-0.5).^2;
end
end

function f=g13(x) %DTLZ7
f=1+(9/sqrt(sum(x.^2)))*sum(x);
end
function f=h(x)
global f_num
f=f_num-(x(1)/(1+g13(x)))*(1+sin(3*pi*x(1)));
end测试主m文件

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DE解决多目标优化问题-----传统的差分进化算法
%时间：2021.10.27
%开发者：gouhuipeng
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
Np=100;  %种群数
D=30;     %个体数（问题的维度）
F=0.5;   %突变因子（变异因子）
Cr=0.3;  %交叉因子
%xmax(1)=1;   %问题上界
%xmin(1)=0;   %问题下届
for i=1:D
    xmin(:,i)=0;
    xmax(:,i)=1;
end
type=5;   %测试问题的类型
K=500;    %迭代次数

for j=1:D
    x(:,j)=xmin(:,j)+(xmax(:,j)-xmin(:,j))*rand(Np,1);  %产生初始值------第一代----对于这里的初始值可不可以进行给定一个最优初始值
end
best=x(1,:);  %全局最优个体 ---之后不断更新--令其等于我们的初始值的第一个值
%end
%plot(x(:,1),x(:,2),'.');
  for i=2:Np
     if(muti_fun(x(i,:),type)<muti_fun(best(1,:),type))   %编写函数有问题，下去思考
         best=x(i,:);                     %选择最优个体
     end
      end         
muti_y=muti_fun(best,type);   %保证个体全局最优

%%进入循环直到满足精度要求或者迭代次数达到
  for n=1:K            %做迭代
     time(n)=n;
     %F=0.9-(n*(0.9-0.5)./K);   %采用线性递减的方式做
     %第二步 变异
      for i=1:Np
           for j=1:D %检查是否越界
         r1=1;r2=1;r3=1;%使得个体满足变异条件.此处与Java有点不一样，他是从1开始
          while(r1==r2||r1==r3||r2==r3||r1==i||r2==i||r3==i)
            r1=ceil(Np*rand(1));          %保持其中的r1,r2,r3互异，这样做的目的是为了防止种群的单一性
            r2=ceil(Np*rand(1));
            r3=ceil(Np*rand(1));
            end
          v(i,j)=x(r1,j)+F*(x(r2,j)-x(r3,j));  
         
         %做一个防止越界
              if v(i,j)<xmin(:,j)
                  v(i,j)=xmin(:,j);
              elseif x(i,j) >xmax(:,j)
                  x(i,j)=xmax(:,j);
              end
           end
         
          %交叉
        for j=1:D
        temper=rand;  %随机产生一个数，用于进行多点交叉
        if(temper<Cr)
            u(i,j)=v(i,j);
        else
            u(i,j)=x(i,j);
        end
        end
        mutifun_Uresult(i,:)=muti_fun(u(i,:),type);
        mutifun_Xresult(i,:)=muti_fun(x(i,:),type);
        %选择----------修改选择方法
       % if(muti_fun(u(i,:),type)<=muti_fun(x(i,:),type))
       %     x(i,:)=u(i,:);
       % end
       % plot(i,x,'bo');
        if  mutifun_Uresult(i,:)<= mutifun_Xresult(i,:)
            x(i,:)=u(i,:);
            
        end
        fun_muti_x(i,:)=muti_fun(x(i,:),type);
        % plot( fun_muti_x(i,1), fun_muti_x(i,2),'b--');
        if fun_muti_x(i,:)< muti_y
            %fi=fun_DE1(x(i,:));
            best=x(i,:);
         
        end
           Best_f(n,:)=muti_fun(best,type);  
      end   
  end
  %ZDT1理论前沿面
  x1=0:0.01:1;
  f2=1-sqrt(x1);
  plot(x1,f2,'k-');title('ZDT1')
  %ZDT2理论前沿面
% x1=0:0.01:1;
% f2=1-x1.^2;
% plot(x1,f2,'k-');
   %ZDT3理论前沿面
    %ZDT4理论前沿面
     %ZDT6理论前沿面
  hold on
  plot( fun_muti_x(:,1), fun_muti_x(:,2),'k+');
  xlabel('f-1');ylabel('f-2');zlabel('f-3');
% legend('optimal PF','DE')
  %optvalue=Best_f;
  %fprintf('最优解结果为%f,%f,%f,%f,%f,%f,%f,%f,%f,%f',best);
  disp(['最小函数值为:', num2str(Best_f(n,:))]);三 输出结果




对于想学习matlab和差分进化算法的，推荐如下几本书，这里建议大家做一下机器学习和进化算法的结合。
欢迎关注知乎，留言咨询！