platEMO：一款强大的多目标优化工具（MATLAB）的扩展操作——操作实例
孙承
前言
后台有很多人催更我关于platEMO的操作实例，实在是由于科研期间琐事太多了，因此文章更新就停滞了（其实说白了就是自己太懒了），根据大家问我最多的几个问题，这次博客主要分为两部分：增加问题和增加算法。一些架构和官方的介绍可以参考我之前的两篇博客。这篇博客直接说怎么添加。除了matlab以外，还有一个python有关进化算法的库叫做geatpy，也非常不错，他的网址：http://geatpy.com，里面有详细的教程，感兴趣的朋友可以去学习一下。
增加问题
每一个多目标优化问题在platEMO中由一个.m文件表示，该文件应该放在文件夹Problems中。即我们首先需要在Problems这个文件夹新建一个我们问题的文件夹。我们定义一个多目标优化问题一般需要考虑他的问题个数Global.M，变量的个数Global.D，变量的上界Global.upper，变量的下界Global.lower，以及变量的编码方式Global.encoding，最后就是目标函数的形式(包括约束条件)了。以下我们以带约束的DTLZ1()函数为例说明问题是如何增加的。
classdef C1_DTLZ1 < PROBLEM
% <problem> <DTLZ variant>%如果你想在GUI界面去调用你的问题，这里<>必须起一个你问题自己的名字
% Constrained DTLZ1%如果你想在GUI界面去调用你的问题，这里<>必须起一个你问题自己的名字
methods
%% Initialization
function obj = C1_DTLZ1()
if isempty(obj.Global.M)
obj.Global.M = 3;%定义你优化问题的目标个数
end
if isempty(obj.Global.D)
obj.Global.D = obj.Global.M + 4;%定义你优化问题的变量个数，在DTLZ1中为目标个数加4
end
obj.Global.lower = zeros(1,obj.Global.D);定义你的优化问题的下界
obj.Global.upper = ones(1,obj.Global.D);定义你优化问题的上界
obj.Global.encoding = 'real';%定义你优化问题的编码方式，例如实数编码码'real',platEMO提供三种编码方式：'binary'二进制编码；'permutation'排列编码以及其他实数编码，根据你的优化问题的实际搜索空间可以自行选择
end
%% Calculate objective values计算目标函数值
function PopObj = CalObj(obj,PopDec)
[N,D] = size(PopDec);%首先我们需要得到种群大小N和自变量个数D
M = obj.Global.M;%得到目标个数
g = 100*(D-M+1+sum((PopDec(:,M:end)-0.5).^2-cos(20.*pi.*(PopDec(:,M:end)-0.5)),2));
PopObj = 0.5*repmat(1+g,1,M).*fliplr(cumprod([ones(N,1),PopDec(:,1:M-1)],2)).*[ones(N,1),1-PopDec(:,M-1:-1:1)];%给出自己问题的表达式，如果是黑盒优化问题，这里需要将黑盒问题抽象成包含输入和输出的接口形式，注意输入的矩阵PopDec每一行代表种群中的一条染色体，黑盒问题需要可以接受这样的输入形式，如果黑盒优化的函数只能一个个染色体计算，那么需要用for循环算出每一个染色体在不同目标上对应的目标函数值，再合成矩阵PopObj。
end
%% Calculate constraint violations计算目标函数的约束值，没有约束的话，可以不写这个接口函数
function PopCon = CalCon(obj,PopDec)
PopObj = obj.CalObj(PopDec);
PopCon = PopObj(:,end)/0.6 + sum(PopObj(:,1:end-1)/0.5,2) - 1;%值得注意的是，约束必须要写成g(x)>=0的形式哦，这里计算的g(x)的值
end
%% Sample reference points on Pareto front%这里是优化问题的最优PF，在实际问题中我们是得不到最优的PF的，因此实际问题时直接省略就可以的，这里C1_DTLZ1函数是学者们构造出来用来测试算法性能的，因此我们可以得到PF的理论解。
function P = PF(obj,N)
P = UniformPoint(N,obj.Global.M)/2;
end
end
end
首先由于C1_DTLZ1是PROBLEM的一个子类，每个操作都是一个重载方法。以上每条代码都给出来详细的注释，参考这些就可以加入自己的方法啦。一下有两点值得注意的事项：
- 决策变量在三种情况下可能是不可行的。首先，对于连续MOPs，它可能大于全局上界（GLOBAL.upper）或小于下界（GLOBAL.lower）。在这种情况下，它将被INDIVIDUAL的类设置为边界值，而MOP类不需要处理这种情况。其次，它可能没有满足约束(一个正的约束违背表明这个约束没有被满足)，在这种情况下，约束违背应该通过重载方法PROBLEM.CalCon()来计算。最后，对于组合MOPs，它可能是一个非法字符，在这种情况下，它应该通过重载方法PROBLEM.CalDec()来修复，例如在MOKP.m中：
2. function PopDec = CalDec(obj,PopDec)
3. C = sum(obj.W,2)/2;
4. [~,rank] = sort(max(obj.P./obj.W));
5. for i = 1 : size(PopDec,1)
6. while any(obj.W*PopDec(i,:)’>C)
7. k = find(PopDec(i,rank),1);
8. PopDec(i,rank(k)) = 0;
9. end
10. end
11. end
- 问题类还可以接受用户的参数设置，比如我们目标函数中有一个参数可能需要经常改变，我们可以把他定义为该问题类的参数，例如：WFG1中的参数K
classdef WFG1 < PROBLEM
% <problem> <WFG>
% Benchmark MOP proposed by Walking Fish Group
% K --- --- The position parameter, which should be a multiple of M-1%如果你想在GUI界面显示，这里必须写上
properties(Access = private)
K; % Position parameter
end
methods
%% Initialization
function obj = WFG1()
if isempty(obj.Global.M)
obj.Global.M = 3;
end
if isempty(obj.Global.D)
obj.Global.D = obj.Global.M + 9;
end
obj.K = obj.Global.ParameterSet(obj.Global.M-1);
obj.Global.lower = zeros(1,obj.Global.D);
obj.Global.upper = 2 : 2 : 2*obj.Global.D;
obj.Global.encoding = 'real';
end
%% Calculate objective values
function PopObj = CalObj(obj,PopDec)
[N,D] = size(PopDec);
M = obj.Global.M;
K = obj.K;
L = D - K;
D = 1;
S = 2 : 2 : 2*M;
A = ones(1,M-1);
z01 = PopDec./repmat(2:2:size(PopDec,2)*2,N,1);
t1 = zeros(N,K+L);
t1(:,1:K) = z01(:,1:K);
t1(:,K+1:end) = s_linear(z01(:,K+1:end),0.35);
t2 = zeros(N,K+L);
t2(:,1:K) = t1(:,1:K);
t2(:,K+1:end) = b_flat(t1(:,K+1:end),0.8,0.75,0.85);
t3 = zeros(N,K+L);
t3 = b_poly(t2,0.02);
t4 = zeros(N,M);
for i = 1 : M-1
t4(:,i) = r_sum(t3(:,(i-1)*K/(M-1)+1:i*K/(M-1)),2*((i-1)*K/(M-1)+1):2:2*i*K/(M-1));
end
t4(:,M) = r_sum(t3(:,K+1:K+L),2*(K+1):2:2*(K+L));
x = zeros(N,M);
for i = 1 : M-1
x(:,i) = max(t4(:,M),A(i)).*(t4(:,i)-0.5)+0.5;
end
x(:,M) = t4(:,M);
h = convex(x);
h(:,M) = mixed(x);
PopObj = repmat(D*x(:,M),1,M) + repmat(S,N,1).*h;
end
%% Sample reference points on Pareto front
function P = PF(obj,N)
M = obj.Global.M;
P = UniformPoint(N,M);
c = ones(size(P,1),M);
for i = 1 : size(P,1)
for j = 2 : M
temp = P(i,j)/P(i,1)*prod(1-c(i,M-j+2:M-1));
c(i,M-j+1) = (temp^2-temp+sqrt(2*temp))/(temp^2+1);
end
end
x = acos(c)*2/pi;
temp = (1-sin(pi/2*x(:,2))).*P(:,M)./P(:,M-1);
a = 0 : 0.0001 : 1;
E = abs(temp*(1-cos(pi/2*a))-1+repmat(a+cos(10*pi*a+pi/2)/10/pi,size(x,1),1));
[~,rank] = sort(E,2);
for i = 1 : size(x,1)
x(i,1) = a(min(rank(i,1:10)));
end
P = convex(x);
P(:,M) = mixed(x);
P = repmat(2:2:2*M,size(P,1),1).*P;
end
end
end
function Output = s_linear(y,A)
Output = abs(y-A)./abs(floor(A-y)+A);
end
function Output = b_flat(y,A,B,C)
Output = A+min(0,floor(y-B))*A.*(B-y)/B-min(0,floor(C-y))*(1-A).*(y-C)/(1-C);
Output = roundn(Output,-6);
end
function Output = b_poly(y,a)
Output = y.^a;
end
function Output = r_sum(y,w)
Output = sum(y.*repmat(w,size(y,1),1),2)./sum(w);
end
function Output = convex(x)
Output = fliplr(cumprod([ones(size(x,1),1),1-cos(x(:,1:end-1)*pi/2)],2)).*[ones(size(x,1),1),1-sin(x(:,end-1:-1:1)*pi/2)];
end
function Output = mixed(x)
Output = 1-x(:,1)-cos(10*pi*x(:,1)+pi/2)/10/pi;
end
这样定义之后，我们就可以直接在外部通过下面命令接受用户设定的参数了：
1. obj.K = obj.Global.ParameterSet(obj.Global.M-1);
增加算法
接下来，我们主要讲如何增加自己设计的多目标进化算法(MOEA)了。MOEA函数由PlatEMO中的.m文件表示，该文件应该放在Algorithm文件夹中。即我们首先需要在Algorithm这个文件夹新建一个我们算法的文件夹。我们知道一个进化算法最重要的就是五个部分，初始化，函数评价，交配选择，交配产解，环境选择，接下来我们以NSGAII为例,说明如何在PlatEMO中增加一个多目标进化算法。
首先我们给出NSGAII的主流程代码NSGAII.m：
function NSGAII(Global)
% <algorithm> <N>
% Nondominated sorting genetic algorithm II%如果你想在GUI界面去调用你的算法，这里<>必须起一个你算法的名字
%% Generate random population
Population = Global.Initialization();%这里是对种群进行初始化，Global类中已经完成了，不需要我们自己编写，如果你想设计特殊的初始化方式，可以在Population的基础上进行自己的操作即可。值得注意的是，在platEMO中， 个体被初始化或者利用遗传算子产生之后（INDIVIDUAL目标被实例化）就会被自动评价，这是不需要我们人为操作的。
[~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Global.N);%这里是对种群中的个体进行选择，详情见后面的分析。
%% Optimization
while Global.NotTermination(Population)%这里是对算法停止条件的判断，里面包含了函数评估次数的检查，以及一些数据保存和画图等操作。如果需要修改算法的停止条件，需要自己修改Global，个人建议不要修改，因为如果对架构不了解有可能改错。函数评估次数是目前多目标优化算法中最常用的停止条件了，毕竟MOEA的收敛性并不好最为停止条件去度量。
MatingPool = TournamentSelection(2,Global.N,FrontNo,-CrowdDis);%这里是进行交配选择，即给出交配的父代，platEMO中给出了常用的选择方式,比如二进制锦标赛和轮盘赌选择，他们在Public这个文件夹中，我们可以直接调用，需要自己设计选择方式的话，按照这些文件修改即可。
Offspring = GA(Population(MatingPool));%这里是进行子代产生，paltEMO提供了很多子代产生方式，包含：DE,FEP,GA,PSO等，我们可以直接调用，这些函数位于Operater中。对于函数GA()，如果输入是INDIVIDUAL对象的数组，那么输出也是INDIVIDUAL对象的数组。如果输入是决策变量矩阵，那么输出也是决策变量矩阵。
[Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Global.N);%这里是进行自然选择，在MOEA中，如何选择解以兼顾收敛性和多样性是学者们研究的重点之一，这里给出的是NSGAII中最经典的方法.
end
end
EnvironmentalSelection.m
function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N)
% The environmental selection of NSGA-II
%% Non-dominated sorting
[FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);%这里进行非支配排序，这是MOEA中最常用的一种方法
Next = FrontNo < MaxFNo;
%% Calculate the crowding distance of each solution
CrowdDis = CrowdingDistance(Population.objs,FrontNo);%在非支配排序的基础上NSGAII提出来拥挤度距离以考虑了解的多样性。
%% Select the solutions in the last front based on their crowding distances
Last = find(FrontNo==MaxFNo);
[~,Rank] = sort(CrowdDis(Last),'descend');
Next(Last(Rank(1:N-sum(Next)))) = true;
%% Population for next generation
Population = Population(Next);
FrontNo = FrontNo(Next);
CrowdDis = CrowdDis(Next);
end
除了上述提到的基于支配关系的方法，基于分解的方法（eg.MOEA/D）以及基于指标的方法(eg.ARMOEA)也非常常见，这些方法需要给出针对性的参数设置以及参考点设置等等，大家参考这类方法的代表性方法的伪代码即可。
[1]: Ye Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, “PlatEMO: A MATLAB
platform for evolutionary multi‐objective optimization [educational
forum],” IEEE Computational Intelligence Magazine, 2017, 12(4): 73‐87.
类似资料：