失眠网 > 【多目标进化优化】NSGAII 算法原理与代码实现

【多目标进化优化】NSGAII 算法原理与代码实现

时间：2022-05-04 13:49:47

🎉博主相信：有足够的积累，并且一直在路上，就有无限的可能！！！

👨‍🎓个人主页：青年有志的博客
💯Gitee 源码地址：/futurelqh/Multi-objective-evolutionary-optimization

理论分析

\quad\quad Srinivas 和 Deb 于 1993 年提出了 NSGA (non-dominated sorting in genetic algorithm) (Srinivas et al, 1994)。NSGA 主要有三个方面的不足：一是没有最优个体 (elitist) 保留机制，有关研究表明（如文献Zitzler et al, 2000； Rudolph, 2001），最优个体保留机制一方面可以提高 MOEA 的性能，同时也能防止优秀解的丢失；二是共享参数问题，在进化过程中，主要是采用共享参数 σshare\sigma_{share}σshare 来维持解群体的分布性（基于共享机制的小生境技术），但共享参数的大小不容易确定，参数的动态修改和调整更是一件困难的工作；三是构造 Pareto 最优解集（通常是构造进化群体的非支配集）的时间复杂度高，为 O(rN3)O(rN^3)O(rN3) （这里 rrr 为目标数，NNN 为进化群体的规模），因为每一代进化都需要构造非支配集，这样一来，当进化群体规模较大时，算法执行的时间开销就很大。为此，Deb 等于 2000 年在 NSGA 的基础上，提出了 NSGAII (Deb et al, 2000)。

\quad\quad 在 NSGA-II 中，将进化群体按支配关系分为若干层，第一层为进化群体的非支配个体集合，第二层为在进化群体中去掉第一层个体后所求得的非支配个体集合，第三层为在进化群体中去掉第一层和第二层个体后所求得的非支配个体集合，依此类推。选择操作首先考虑第一层非支配集，按照某种策略从第一层中选取个体；然后再考虑在第二层非支配个体集合中选择个体，依此类推，直至满足新进化群体的大小要求。

\quad\quad 下面从非支配集的构造方法、维持解集分布性的策略等方面进行比较详细的讨论。

\quad\quad 1. 非支配集的构造方法

\quad\quad 设群体 PopPopPop 的规模大小为 NNN，将群体 PopPopPop 按某种策略进行分类排序为 mmm 个子集 P1、P2、...、PmP_1、P_2、...、P_mP1、P2、...、Pm，且满足下列性质：

\quad\quad ① ∪p∈{P1,P2,...,Pm}P=Pop\cup_{p \in \{P_1, P_2, ..., P_m\}} P = Pop∪p∈{P1,P2,...,Pm}P=Pop（所有子集的并集为 PopPopPop）

\quad\quad ② ∀i,j∈{1,2,...,m}\forall i, j \in \{1, 2, ..., m\}∀i,j∈{1,2,...,m} 且 i≠j，Pi∩Pj=∅i ≠ j，P_i \cap P_j = \emptyseti=j，Pi∩Pj=∅（所有子集不存在相同的个体）

\quad\quad ③ P1≻P1≻⋅⋅⋅≻PmP_1 \succ P_1 \succ ··· \succ P_mP1≻P1≻⋅⋅⋅≻Pm，即 Pk+1P_{k+1}Pk+1 中的个体直接受 PkP_kPk 中个体的支配（k=1,2,...,m−1）（k=1, 2, ..., m-1）（k=1,2,...,m−1）。

\quad\quad 对群体 PopPopPop 进行分类排序的目的是为了将其划分成若干个满足上述三个性质的互不相交的子群体。对个体分类排序的依据为 ParetoParetoPareto 支配关系。

\quad\quad 设两个向量 {np}\{n_p\}{np} 和 {sp}\{s_p\}{sp}，其中p∈Pop，npp \in Pop，n_pp∈Pop，np 记录支配个体 ppp 的个体数，{sp}\{s_p\}{sp} 记录被个体 ppp 支配的个体的集合，即有

np=∣{q∣q≻pp,q∈Pop}∣n_p = |\{q|q \succ p \quad p, q \in Pop\}|np=∣{q∣q≻pp,q∈Pop}∣

sp=∣{q∣p≻qp,q∈Pop}∣s_p = |\{q|p \succ q \quad p, q \in Pop\}|sp=∣{q∣p≻qp,q∈Pop}∣

\quad\quad 首先通过一个二重循环计算每个个体的 npn_pnp 和 sps_psp ，则 P1={q∣nq=0，q∈Pop}P_1 = \{q~|~n_q = 0，q \in Pop\}P1={q∣nq=0，q∈Pop}。然后依次按方法 Pk={所有个体q∣nq−k+1=0}P_k = \{所有个体~ q~ | ~n_q - k + 1 = 0\}Pk={所有个体q∣nq−k+1=0} 求 P2、P3⋅⋅⋅P_2、P_3 ···P2、P3⋅⋅⋅。

\quad\quad 构造分类子集的具体过程如算法 6.2 所示，其中，P1P_1P1 为非支配集。

\quad\quad 算法 6.2 由两部分组成，第一部分用于计算 nin_ini 和 sis_isi，并求得 P1P_1P1，所需要的时间为 (rN2)(rN^2)(rN2) ，这里 rrr 为目标数，NNN 为进化群体规模；第二部分用于求 P2、P3、⋅⋅⋅、PmP_2、P_3、···、P_mP2、P3、⋅⋅⋅、Pm，最坏情况下一个规模为 NNN 的进化群体有 NNN 层边界集 (Pareto front)，即 m=Nm=Nm=N , 此时其时间复杂度为 O(N2)O(N^2)O(N2)。由此可得，算法 6.2 的总时间复杂度为 O(rN2)十O(N2)O(rN^2)十O(N^2)O(rN2)十O(N2)，即为 O(rN2)O(rN^2)O(rN2)。

\quad\quad Deb 等于 2002 年又提出了一种新的构造非支配集的方法 (Debet al, 2002)，其时间复杂度仍为 O(N2)O(N^2)O(N2)

\quad\quad 2. 保持解群体分布性和多样性的方法

\quad\quad 为了保持解群体的分布性和多样性，Deb 等在文献 (Debet al, 2000) 中，首先通过计算进化群体中每个个体的聚集距离，然后依据个体所处的层次及其聚集距离，定义一个偏序集 (partial order set)，构造新群体时依次在偏序集中选择个体。

\quad\quad 在产生新群体时，通常将优秀且聚集密度比较小的个体保留并参与下一代进化。聚集密度小的个体其聚集距离反而大，一个个体的聚集距离可以通过计算与其相邻的两个个体在每个子目标上的距离差之和来求取。如图 6.6 所示，设有两个子目标 f1f_1f1 和 f2f_2f2，个体 iii 的聚集距离是图中虚线四边形的长与宽之和。设 P[i]distanceP[i]_{distance}P[i]distance 为个体 iii 的聚集距离，P[i].mP[i].mP[i].m 为个体 iii 在子目标 mmm 上的函数值，则图 6.6 中个体 iii 的聚集距离为

P[i]distance=(P[i+1].f1−P[i−1].f1)+(P[i+1].f2−P[i−1].f2)P[i]_{distance} = (P[i + 1].f_1 - P[i-1].f_1) + (P[i + 1].f_2 - P[i-1].f_2)P[i]distance=(P[i+1].f1−P[i−1].f1)+(P[i+1].f2−P[i−1].f2)

\quad\quad 一般情况下，当有 rrr 个子目标时个体 iii 的聚集距离为

P[i]distance=∑k=1r(P[i+1].fk−P[i−1].fk)P[i]_{distance} = \sum_{k=1}^{r}( P[i + 1].f_k - P[i-1].f_k)P[i]distance=k=1∑r(P[i+1].fk−P[i−1].fk)

\quad\quad 为了计算每个个体的聚集距离，需要对群体 PPP 按每个子目标函数值进行排序，当釆用最好的排序方法时(如快速排序、堆排序等)，若群体规模为 NNN，最坏情况下对 rrr 个子目标分别进行排序的时间复杂度为 O(rNlogN)O(rNlogN)O(rNlogN)

\quad\quad 计算个体聚集距离的方法如算法 6.3 所示(TamakietaL 1995)

\quad\quad 算法 6.3 在最坏情况下对 rrr 个子目标分别进行排序的时间为 O(rNlogN)O(rNlogN)O(rNlogN)，计算每个个体的聚集距离的时间为 O(rN)，O(rN)，O(rN)，因此算法的时间复杂度为 O(rNlogN)O(rNlogN)O(rNlogN)。

\quad\quad 3. Deb 的 NSGA-II

\quad\quad 在具体讨论 NSGA-II 之前，先讨论建立在进化群体上的一类偏序关系，因为 NSGA-II 在构造新群体时，将依据这种偏序关系进行选择操作。定义进化群体的偏序关系时，主要考虑下列两个因素：

\quad\quad ① 个体 iii 的分类序号 irank，irank=ki_{rank}，i_{rank}=kirank，irank=k 当且仅当 i∈Pki \in P_ki∈Pk \quad（根据支配关系给每个个体编号： iranki_{rank}irank = 支配 iii 的个体数 + 1，i∈Pki \in P_ki∈Pk：表示 iii 是属于第 k 层非支配集）

\quad\quad ② 个体 iii 的聚集距离 P[i]distanceP[i]_{distance}P[i]distance \quad\quad 详细介绍见 3.2

\quad\quad 得到偏序关系的定义如下。

\quad\quad定义 6.1\quad 设个体 iii 和 jjj 为进化群体中的任意个体，个体之间的偏序(partial order) 关系 ≻n\succ_n≻n 为

i≻njif(irank<jrank)or(irank=jrank)and(P[i]distance>P[j]distance)i \succ_n j \quad if(i_{rank} < j_{rank})~ or~ (i_{rank} = j_{rank}) ~and~ (P[i]_{distance} > P[j]_{distance})i≻njif(irank<jrank)or(irank=jrank)and(P[i]distance>P[j]distance)

\quad\quad 定义 6.1 表明，当两个个体属于不同的分类排序子集时，优先考虑序号 iranki_{rank}irank 小的个体; 当序号 iranki_{rank}irank 相同时，则优先考虑聚集距离大或者说聚集密度小的个体。聚集距离，聚集密度详细介绍见 3.

\quad\quad 在 NSGA-II 中，开始时随机产生一个初始群体 P0P_0P0，在此基础上采用二元锦标赛选择 (binary tournament selection) 、交叉(crossover)和变异操作(mutation) 产生一个新群体 Q0Q_0Q0，P0P_0P0 和 Q0Q_0Q0 的群体规模均为 NNN。将 PtP_tPt 和 QtQ_tQt 并入到 RtR_tRt 中（初始时 t=0t = 0t=0）,对 RtR_tRt 进行分类排序，然后根据需要计算某个分类排序子集中所有个体的聚集距离，并按照定义 6.1 建立偏序集。接下来从偏序集依次选取个体进入 Pt+1P_{t+1}Pt+1，直至 Pt+1P_{t+1}Pt+1 的规模为 NNN。具体过程如算法 6.4 所示。

\quad\quad 算法 6.4 中，通过 F = nondominated-sort(RtR_tRt) 产生了若干个分类子集 F=(F1,F2,⋅⋅⋅)F=(F_1, F_2, ···)F=(F1,F2,⋅⋅⋅)，但被选入新群体的只有一部分。如图 6.7 所示，分类子集 F1F_1F1 和 F2F_2F2 中的所有个体均被选入了新群体 Pt+1P_{t +1}Pt+1，但分类子集 F3F_3F3 中只有一部分个体被选入新群体 Pt+1P_{t + 1}Pt+1。一般地，若 ∣F1∣+∣F2∣+⋅⋅⋅+∣Fi−1∣≤N|F_1| + |F_2| + ··· + |F_{i-1}| ≤ N∣F1∣+∣F2∣+⋅⋅⋅+∣Fi−1∣≤N 且 ∣F1∣+∣F2∣+⋅⋅⋅+∣Fi∣>N|F_1| + |F_2| + ··· + |F_{i}| > N∣F1∣+∣F2∣+⋅⋅⋅+∣Fi∣>N，则称 FiF_iFi 为临界层分类子集，图 6.7 中的 F3F_3F3 为临界层分类子集。

\quad\quad 算法 6-4 的时间开销主要由三部分组成(其中 rrr 为目标数)：

\quad\quad ① 构造分类子集(non-dominated sort) ：O(r(2N)2)O(r(2N)^2)O(r(2N)2)。

\quad\quad ② 计算聚集距离(crowding distance assignment)： O(r(2N)log(2N))O(r(2N )log(2N))O(r(2N)log(2N))。

\quad\quad ③ 构造偏序集(sorting on >〃)： O(2NIog(2N))O(2NIog(2N))O(2NIog(2N))。

\quad\quad 由此可得算法 6.4 的总时间复杂度为 O(rN2)O(rN^2)O(rN2)，其中主要的时间开销花费在构造边界集上，因此一个快速的构造分类子集（或构造非支配集）的方法有利于提高 MOEA 的效率。

算法运行效果

当选用DTLZ1作为测试函数，目标数设置为 3 时，MATLAB 命令：NSGAII('DTLZ1', 3)

当选用ZDT1作为测试函数时， MATLAB 命令：NSGAII('ZDT1', 2)

利用综合指标（HV、IGD）测试算法在测试函数上的收敛性与分布性，在迭代完成后加上下面两段代码即可：

zdt1.mat：为测试函数 zdt1 的真是 Pareto 解集，在文章开头给出的 Gitee 地址中有

load('zdt1.mat'); % 真实 Pareto 解集HVResult = HV(FunctionValue, zdt1)IGDResult = IGD(FunctionValue, zdt1)

代码实现

main.m

function main()% UserInput.Probleme = problem;% UserInput.ObjectNum = dimension;[UserInput] = GetUserInput()switch (UserInput.Probleme)case 1AlgoName = 'DTLZ1';case 2AlgoName = 'DTLZ2';case 3AlgoName = 'DTLZ3';case 4AlgoName = 'DTLZ4';endNSGAII(AlgoName, UserInput.ObjectNum);end

NSGAII.m

function NSGAII(Problem,M)% Function: NSGAII(Problem,M)%% Description: 执行 NSGAII 算法的所有逻辑%%% Syntax:% %% Parameters:% Problem：选择测试函数，例如：'DTLZ1',注意输出字符% M ：目标函数个数%% Return:% 无%%Young99% Revision:1.0 Data: -12-06%*************************************************************************clc;format compact;% 空格紧凑tic; % 记录运行时间% 参数设定Generations = 700; % 迭代次数if M == 2 % 目标数量N = 100; % 种群大小else M == 3N = 105;end% 初始化种群，返回初始化种群和决策变量上下限[Population,Boundary] = Objective(0, Problem, M, N);FunctionValue = Objective(1,Problem,M,Population);% 计算目标函数值% 进行非支配排序FrontValue = NonDominateSort(FunctionValue,0);CrowdDistance = CrowdDistances(FunctionValue,FrontValue); % 计算聚集距离%开始迭代for Gene = 1 : Generations%产生子代。MatingPool = Mating(Population,FrontValue,CrowdDistance); % 交配池选择。2 的锦标赛选择方式Offspring = NSGA_Gen(MatingPool,Boundary,N); % 交叉,变异，越界处理并生成新的种群%% 精英保留策略Population = [Population; Offspring]; % 种群合并FunctionValue = Objective(1,Problem,M,Population); % 计算目标函数值[FrontValue,MaxFront] = NonDominateSort(FunctionValue,1); % 进行非支配排序CrowdDistance = CrowdDistances(FunctionValue,FrontValue); % 计算聚集距离% 选出非支配的个体Next = zeros(1,N);NoN = numel(FrontValue,FrontValue<MaxFront);Next(1:NoN) = find(FrontValue<MaxFront);% 选出最后一个面的个体Last = find(FrontValue==MaxFront);[~,Rank] = sort(CrowdDistance(Last),'descend');Next(NoN+1:N) = Last(Rank(1:N-NoN));% 下一代种群Population = Population(Next,:);FrontValue = FrontValue(Next);CrowdDistance = CrowdDistance(Next);FunctionValue = Objective(1,Problem,M,Population);cla;%cla 从当前坐标区删除包含可见句柄的所有图形对象%绘图for i = 1 : MaxFrontFrontCurrent = find(FrontValue==i);DrawGraph(FunctionValue(FrontCurrent,:));hold on;switch Problem %case 'DTLZ1'if M == 2pareto_x = linspace(0,0.5);pareto_y = 0.5 - pareto_x;plot(pareto_x, pareto_y, 'r');elseif M == 3[pareto_x,pareto_y] = meshgrid(linspace(0,0.5));pareto_z = 0.5 - pareto_x - pareto_y;axis([0,1,0,1,0,1]);mesh(pareto_x, pareto_y, pareto_z);endotherwiseif M == 2pareto_x = linspace(0,1);pareto_y = sqrt(1-pareto_x.^2);plot(pareto_x, pareto_y, 'r');elseif M == 3[pareto_x,pareto_y,pareto_z] =sphere(50);axis([0,1,0,1,0,1]);mesh(1*pareto_x,1*pareto_y,1*pareto_z);endendpause(0.01);endclc;endend

Objective.m：进行种群初始化

function [Output,Boundary] = Objective(Operation,Problem,M,Input)% Function: [Output,Boundary] = Objective(Operation,Problem,M,Input)%% Description: Initial population and compute objective function value%%% Syntax:% %% Parameters:% Operation: if Operation is 0, Initial population. if Operation is 1,% comput funcitn value% Problem：Test function name. such as 'DTLZ1'% M: function numbers% Input: Population numbers%% Return:% Output: if Operation is 0, Output is Initial result, Operation is 1,% Output is functionvalue%% Boundary: function boundary% %Young99%Revision: 1.0 Data: -12-06%*************************************************************************% input:N 种群规模 persistent K; % 定义持久变量% 持久变量是声明它们的函数的局部变量；% 但其值保留在对该函数的各次调用所使用的内存中。Boundary = NaN; switch Operation % 选择操作模式 0/1% 0：初始化种群case 0k = find(~isstrprop(Problem,'digit'),1,'last'); % 判断输入的 Problem 有几个英文字母%“~” 对0/1数组取反K = [5 10 10 10 10];K = K(str2double(Problem(k+1:end))); %'DTZL2'K=10%str2double是一种函数，其功能是把字符串转换数值，D = M + K-1; % D：变量数量, M：目标函数个数 K： K = D - M + 1MaxValue = ones(1,D);MinValue = zeros(1,D);Population = rand(Input,D); % 确定种群大小：( N * D)Population = Population.*repmat(MaxValue,Input,1)+(1-Population).*repmat(MinValue, Input, 1);% 产生新的初始种群Output = Population;Boundary = [MaxValue;MinValue];% 1：计算目标函数值；这里只包含 DTLZ1~DTZL4 函数问题case 1Population = Input; % 已经初始化完成的种群FunctionValue = zeros(size(Population,1),M);switch Problemcase 'DTLZ1'g = 100*(K+sum((Population(:,M:end)-0.5).^2-cos(20.*pi.*(Population(:,M:end)-0.5)),2));for i = 1 : M % 计算第 i 维目标函数值FunctionValue(:,i) = 0.5.*prod(Population(:,1:M-i),2).*(1+g);if i > 1FunctionValue(:,i) = FunctionValue(:,i).*(1-Population(:,M-i+1));endendcase 'DTLZ2'g = sum((Population(:,M:end)-0.5).^2,2);for i = 1 : MFunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);if i > 1FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));endendcase 'DTLZ3'g = 100*(K+sum((Population(:,M:end)-0.5).^2-cos(20.*pi.*(Population(:,M:end)-0.5)),2));for i = 1 : MFunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);if i > 1FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));endendcase 'DTLZ4'Population(:,1:M-1) = Population(:,1:M-1).^100;g = sum((Population(:,M:end)-0.5).^2,2);for i = 1 : MFunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);if i > 1FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));endendcase 'DTLZ5'g = sum((Population(:,M:end)-0.5).^2,2);for i = 1 : Mtheta =(pi./(4*(1+g))).*(1+2.*g.*Population(:,1:M-i));FunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*theta),2);if i > 1FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*(pi./(4*(1+g))).*(1+2.*g.*Population(:,M-i+1)));endend case 'ZDT1'n = size(Population, 2);f1 = Population(:,1);g = 1+9/(n-1)*sum(Population(:, 2:end), 2);h = 1-sqrt(f1 ./ g);f2 = g .* h;FunctionValue = [f1 f2];endOutput = FunctionValue;endend

NonDominateSort.m：构造非支配集

function [FrontValue,MaxFront] = NonDominateSort(FunctionValue,Operation)% Function: [FrontValue,MaxFront] = NonDominateSort(FunctionValue,Operation)%% Description: 进行非支配排序%%% Syntax:% %% Parameters:% FunctionValue, 待排序的种群(目标空间)，的目标函数值% Operation,可指定仅排序第一个面,排序前一半个体,或是排序所有的个体, 默认为排序所有的个体%% Return:% FrontValue, 排序后的每个个体所在的前沿面编号, 未排序的个体前沿面编号为inf% MaxFront, 排序的最大前面编号%%%Young99%Revision: 1.0 Data: -12-06%*************************************************************************if Operation == 1Kind = 2; elseKind = 1; end[N,M] = size(FunctionValue);MaxFront = 0;cz = zeros(1,N); %% 记录个体是否已被分配编号FrontValue = zeros(1,N)+inf; % 每个个体的前沿面编号初始化为无穷：inf[FunctionValue,Rank] = sortrows(FunctionValue);% sortrows：由小到大以某行的方式进行排序（默认为第一行），返回排序结果和检索到的数据(按相关度排序)在原始数据中的索引%开始迭代判断每个个体的前沿面,采用改进的 deductive sort，Deb 非支配排序算法% 代码思路：最外层的 while 表示，% 当 kind = 1，,且所有个体均编号了，退出循环，% 或 kind = 2，一般个体已经被选出编号了% 或 kind = 3，并且只选出第一层非支配集% while 内嵌套了三层 for 循环，最外层 for 循环固定一个个体，% 内层 for 循环将后面其他所有个体与此个体进行两两比较，将被次个体支配的个体标号 d(j) = 1,% 在后面循环中直接跳过（~d(i)），下次进入 if 语句的第一个个体，即为与个体 i 同层级的个体（与 i 是非支配关系）% 从而在 for 遍历完所有个体后，可选出当前层级的所有个体，并执行 d = cz，在后续筛选个体时，跳过上一次已经筛选出的个体% MaxFront 用于对当前层级的个体进行编号，即第一轮为支配集 P1，第二轮个体为 P2while (Kind==1 && sum(cz)<N) || (Kind==2 && sum(cz)<N/2) || (Kind==3 && MaxFront<1)MaxFront = MaxFront + 1;d = cz;% 最外层两个 for 循环用于个体之间两两比较，最内层的 for 循环用于，选出的两个个体进行每个函数值的比较for i = 1 : Nif ~d(i)for j = i + 1 : Nif ~d(j)k = 1; % 作为标记for m = 2 : Mif FunctionValue(i,m) > FunctionValue(j,m) %比较函数值，判断个体的支配关系k = 0;break;endendif k == 1d(j) = 1; % 将所有被 i 支配的都标记为 1endendendFrontValue(Rank(i)) = MaxFront;cz(i) = 1; % 将当前的个体标记为 1，表示已经确定所在的非支配层级endendendend%% 非支配排序伪代码% def fast_nondominated_sort( P ):% F = [ ]% for p in P:% Sp = [ ] #记录被个体p支配的个体% np = 0 #支配个体p的个数% for q in P:% if p > q:#如果p支配q，把q添加到Sp列表中% Sp.append( q ) #被个体p支配的个体% else if p < q: #如果p被q支配，则把np加1% np += 1 #支配个体p的个数% if np == 0:% p_rank = 1 #如果该个体的np为0，则该个体为Pareto第一级% F1.append( p )% F.append( F1 )% i = 0% while F[i]:% Q = [ ]% for p in F[i]:% for q in Sp: #对所有在Sp集合中的个体进行排序% nq -= 1% if nq == 0:#如果该个体的支配个数为0，则该个体是非支配个体% q_rank = i+2 #该个体Pareto级别为当前最高级别加1。此时i初始值为0，所以要加2% Q.append( q )% F.append( Q )% i += 1

CrowdDistances.m：计算拥挤距离

function CrowdDistance = CrowdDistances(FunctionValue,FrontValue)% Function: CrowdDistance = CrowdDistances(FunctionValue,FrontValue)%% Description: 分前沿面计算种群中每个个体的拥挤距离%%% Syntax:% %% Parameters:% FunctionValue：每个个体的目标函数值% FrontValue：每个个体的最大前沿面编号%% Return:% CrowdDistance：每个个体的拥挤距离%% Young99% Revision: Data: %*************************************************************************[N,M] = size(FunctionValue);CrowdDistance = zeros(1,N);Fronts = setdiff(unique(FrontValue),inf);for f = 1 : length(Fronts)Front = find(FrontValue==Fronts(f));Fmax = max(FunctionValue(Front,:),[],1);Fmin = min(FunctionValue(Front,:),[],1);for i = 1 : M[~,Rank] = sortrows(FunctionValue(Front,i));CrowdDistance(Front(Rank(1))) = inf;CrowdDistance(Front(Rank(end))) = inf;for j = 2 : length(Front)-1CrowdDistance(Front(Rank(j))) = CrowdDistance(Front(Rank(j)))+(FunctionValue(Front(Rank(j+1)),i)-FunctionValue(Front(Rank(j-1)),i)).../(Fmax(i)-Fmin(i)); % 除以最大值 - 最小值，即对其标准化处理endendendend

Mating.m：二元锦标赛选择个体

function MatingPool = Mating(Population,FrontValue,CrowdDistance)% Function:MatingPool = Mating(Population,FrontValue,CrowdDistance)%% Description: 通过二元锦标赛选择方式，选择后代个体% 模拟自然界中一个普遍存在的规律，即优秀个体有更大的概率能够繁衍后代。%% Syntax:% %% Parameters:% Population:种群% FrontValue：每个种群对应的最大前沿面% CrowdDistance：每个种群的拥挤距离%% Return:% MatingPool：经过交配池留下的种群%% Young99% Revision:1.0 Data: -12-06%*************************************************************************[N,D] = size(Population);% 二元联赛选择MatingPool = zeros(N,D);Rank = randperm(N); % randperm 将一列序号随机打乱,序号必须是整数Pointer = 1;for i = 1 : 2 : N%选择父母k = zeros(1,2);for j = 1 : 2if Pointer >= NRank = randperm(N);Pointer = 1;endp = Rank(Pointer);q = Rank(Pointer+1); % 规模为2的锦标赛选择方式if FrontValue(p) < FrontValue(q)k(j) = p;elseif FrontValue(p) > FrontValue(q)k(j) = q;elseif CrowdDistance(p) > CrowdDistance(q)k(j) = p;elsek(j) = q;endPointer = Pointer + 2;endMatingPool(i,:) = Population(k(1),:);MatingPool(i+1,:) = Population(k(2),:);endend

NSGA_Gen.m：模拟二进制交叉、多项式变异

function Offspring = NSGA_Gen(MatingPool,Boundary,MaxOffspring)% Function: Offspring = NSGA_Gen(MatingPool,Boundary,MaxOffspring)%% Description: 通过遗传算法利用模拟二进制交叉,多项式变异生成新的种群 %%% Syntax:% %% Parameters:% MatingPool, 交配池, 其中每第 i 个和第 i+1 个个体交叉产生两个子代, i 为奇数% Boundary,决策空间, 其第一行为空间中每维的上界, 第二行为下界% MaxOffspring, 返回的子代数目, 若缺省则返回所有产生的子代, 即和交配池的大小相同% % Return:% Offspring, 产生的子代新种群%% Young99% Revision: Data: %*************************************************************************[N,D] = size(MatingPool);if nargin < 3 || MaxOffspring < 1 || MaxOffspring > NMaxOffspring = N;end% 遗传操作参数ProC = 1; % 交叉概率ProM = 1/D;% 变异概率DisC = 20; % 交叉参数DisM = 20; % 变异参数% 模拟二进制交叉Offspring = zeros(N,D);for i = 1 : 2 : Nbeta = zeros(1,D);miu = rand(1,D);beta(miu<=0.5) = (2*miu(miu<=0.5)).^(1/(DisC+1));beta(miu>0.5) = (2-2*miu(miu>0.5)).^(-1/(DisC+1));beta = beta.*(-1).^randi([0,1],1,D);beta(rand(1,D)>ProC) = 1;Offspring(i,:) = (MatingPool(i,:)+MatingPool(i+1,:))/2+beta.*(MatingPool(i,:)-MatingPool(i+1,:))/2;Offspring(i+1,:) = (MatingPool(i,:)+MatingPool(i+1,:))/2-beta.*(MatingPool(i,:)-MatingPool(i+1,:))/2;endOffspring = Offspring(1:MaxOffspring,:); % 用于选取自带个体数 MaxOffspring% 多项式变异if MaxOffspring == 1MaxValue = Boundary(1,:);MinValue = Boundary(2,:);elseMaxValue = repmat(Boundary(1,:),MaxOffspring,1);MinValue = repmat(Boundary(2,:),MaxOffspring,1);endk = rand(MaxOffspring,D); % 随机选定要变异的基因位miu = rand(MaxOffspring,D); % 采用多项式变异Temp = k<=ProM & miu<0.5; % 要变异的基因位Offspring(Temp) = Offspring(Temp)+(MaxValue(Temp)-MinValue(Temp)).*((2.*miu(Temp)+(1-2.*miu(Temp)).*(1-(Offspring(Temp)-MinValue(Temp))./(MaxValue(Temp)-MinValue(Temp))).^(DisM+1)).^(1/(DisM+1))-1);Temp = k<=ProM & miu>=0.5; Offspring(Temp) = Offspring(Temp)+(MaxValue(Temp)-MinValue(Temp)).*(1-(2.*(1-miu(Temp))+2.*(miu(Temp)-0.5).*(1-(MaxValue(Temp)-Offspring(Temp))./(MaxValue(Temp)-MinValue(Temp))).^(DisM+1)).^(1/(DisM+1)));%越界处理Offspring(Offspring>MaxValue) = MaxValue(Offspring>MaxValue); %子代越上界处理Offspring(Offspring<MinValue) = MinValue(Offspring<MinValue); %子代越下界处理end

DrawGraph.m：绘图

function Handle = DrawGraph(FigureData,FigureFormat)% Function:Handle = DrawGraph(FigureData,FigureFormat)%% Description: 绘制出指定数据的图形%%% Syntax:% %% Parameters:% FigureData, 待绘制的数据矩阵(即点的坐标集合)% FigureFormat, 图形的格式, 可省略%% Return:% Handle, 图形的句柄%% Young99% Revision: 1.0 Data: -12-06%*************************************************************************A = 1;switch size(FigureData,2)case 2hold on;box onset(gca, 'Fontname', 'Times New Roman', 'Fontsize', 13)if nargin < 2 || ~ischar(FigureFormat)Handle = plot(FigureData(:,1),FigureData(:,2),'ok','MarkerSize',6,'Marker','o','Markerfacecolor',[0 0 0]+0.4,'Markeredgecolor',[0 0 0]+0.7);elseHandle = plot(FigureData(:,1),FigureData(:,2),FigureFormat);endaxis tightxlabel('f_1')ylabel('f_2')case 3hold on;box onset(gca, 'Fontname', 'Times New Roman', 'Fontsize', 13)if nargin < 2 || ~ischar(FigureFormat)Handle = plot3(FigureData(:,1),FigureData(:,2),FigureData(:,3),'ok','MarkerSize',8,'Marker','o','Markerfacecolor',[0 0 0]+0.4,'Markeredgecolor',[0 0 0]+0.7);elseHandle = plot3(FigureData(:,1),FigureData(:,2),FigureData(:,3),FigureFormat);endaxis tightxlabel('f_1')ylabel('f_2')zlabel('f_3')view(135,30)otherwisehold on;box on;if A == 1Num = size(FigureData,1);for i = 1 : Numplot(FigureData(i,:),'bo-');endHandle = NaN;elseif A == 2M = size(FigureData,2);SumM = sum([1:M-1]);RowValue = floor(sqrt(SumM));while mod(SumM,RowValue) ~= 0RowValue = RowValue - 1;endLineValue = SumM / RowValue;Num = 1;for i = 1 : Mfor j= i+1 : Msubplot(RowValue,LineValue,Num);plot(FigureData(:,i),FigureData(:,j),'b*');Num = Num + 1;endendHandle = NaN;endendend

HV.m

function [Score,PopObj] = HV(PopObj,PF)% <metric> <max>% Hypervolume%------------------------------- Reference --------------------------------% E. Zitzler and L. Thiele, Multiobjective evolutionary algorithms: A% comparative case study and the strength Pareto approach, IEEE% Transactions on Evolutionary Computation, 1999, 3(4): 257-271.%------------------------------- Copyright --------------------------------% Copyright (c) - BIMK Group. You are free to use the PlatEMO for% research purposes. All publications which use this platform or any code% in the platform should acknowledge the use of "PlatEMO" and reference "Ye% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform% for evolutionary multi-objective optimization [educational forum], IEEE% Computational Intelligence Magazine, , 12(4): 73-87".%--------------------------------------------------------------------------% % 这段代码是用来计算超体积指标的。它包含了一个 HV 函数，用于计算超体积指标，并输出其值和输入种群的结果。% % HV 函数的输入参数包括种群结果 PopObj 和理论最优结果 PF。它首先会将种群结果进行归一化，然后根据超体积指标的定义，计算种群结果与理论最优结果的相似度。% % 最后，HV 函数会输出计算得到的超体积指标值和归一化后的种群结果。% Normalize the population according to the reference point set[N,M] = size(PopObj);fmin = min(min(PopObj,[],1),zeros(1,M));fmax = max(PF,[],1);PopObj = (PopObj-repmat(fmin,N,1))./repmat((fmax-fmin)*1.1,N,1);PopObj(any(PopObj>1,2),:) = [];% The reference point is set to (1,1,...)RefPoint = ones(1,M);if isempty(PopObj)Score = 0;elseif M < 4% Calculate the exact HV valuepl = sortrows(PopObj);S = {1,pl};for k = 1 : M-1S_ = {};for i = 1 : size(S,1)Stemp = Slice(cell2mat(S(i,2)),k,RefPoint);for j = 1 : size(Stemp,1)temp(1) = {cell2mat(Stemp(j,1))*cell2mat(S(i,1))};temp(2) = Stemp(j,2);S_= Add(temp,S_);endendS = S_;endScore = 0;for i = 1 : size(S,1)p= Head(cell2mat(S(i,2)));Score = Score + cell2mat(S(i,1))*abs(p(M)-RefPoint(M));endelse% Estimate the HV value by Monte Carlo estimationSampleNum = 1000000;MaxValue = RefPoint;MinValue = min(PopObj,[],1);Samples = unifrnd(repmat(MinValue,SampleNum,1),repmat(MaxValue,SampleNum,1));if gpuDeviceCount > 0% GPU accelerationSamples = gpuArray(single(Samples));PopObj = gpuArray(single(PopObj));endfor i = 1 : size(PopObj,1)drawnow();domi = true(size(Samples,1),1);m = 1;while m <= M && any(domi)domi = domi & PopObj(i,m) <= Samples(:,m);m = m + 1;endSamples(domi,:) = [];endScore = prod(MaxValue-MinValue)*(1-size(Samples,1)/SampleNum);endendfunction S = Slice(pl,k,RefPoint)p = Head(pl);pl = Tail(pl);ql = [];S = {};while ~isempty(pl)ql = Insert(p,k+1,ql);p_ = Head(pl);cell_(1,1) = {abs(p(k)-p_(k))};cell_(1,2) = {ql};S = Add(cell_,S);p = p_;pl = Tail(pl);endql = Insert(p,k+1,ql);cell_(1,1) = {abs(p(k)-RefPoint(k))};cell_(1,2) = {ql};S = Add(cell_,S);endfunction ql = Insert(p,k,pl)flag1 = 0;flag2 = 0;ql = [];hp = Head(pl);while ~isempty(pl) && hp(k) < p(k)ql = [ql;hp];pl = Tail(pl);hp = Head(pl);endql = [ql;p];m = length(p);while ~isempty(pl)q = Head(pl);for i = k : mif p(i) < q(i)flag1 = 1;elseif p(i) > q(i)flag2 = 1;endendendif ~(flag1 == 1 && flag2 == 0)ql = [ql;Head(pl)];endpl = Tail(pl);end endfunction p = Head(pl)if isempty(pl)p = [];elsep = pl(1,:);endendfunction ql = Tail(pl)if size(pl,1) < 2ql = [];elseql = pl(2:end,:);endendfunction S_ = Add(cell_,S)n = size(S,1);m = 0;for k = 1 : nif isequal(cell_(1,2),S(k,2))S(k,1) = {cell2mat(S(k,1))+cell2mat(cell_(1,1))};m = 1;break;endendif m == 0S(n+1,:) = cell_(1,:);endS_ = S;end

IGD.m

function Score = IGD(PopObj,PF)% PopObej：算法求得的非支配解集% PF：真实 Pareto 解集% <metric> <min>% Inverted generational distance%------------------------------- Reference --------------------------------% C. A. Coello Coello and N. C. Cortes, Solving multiobjective optimization% problems using an artificial immune system, Genetic Programming and% Evolvable Machines, , 6(2): 163-190.%------------------------------- Copyright --------------------------------% Copyright (c) - BIMK Group. You are free to use the PlatEMO for% research purposes. All publications which use this platform or any code% in the platform should acknowledge the use of "PlatEMO" and reference "Ye% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform% for evolutionary multi-objective optimization [educational forum], IEEE% Computational Intelligence Magazine, , 12(4): 73-87".%--------------------------------------------------------------------------% 通过 pdist 计算 PF 与 PopObj 两两欧式距离，通过 min 取每一行的最小值作为列向量% 最后再取所有距离的平均值，越小说明分布性和收敛性越Distance = min(pdist2(PF,PopObj),[],2); Score = mean(Distance);end