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

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

👨‍🎓个人主页： 青年有志的博客

💯 Gitee 源码地址： https://gitee.com/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 的性能，同时也能防止优秀解的丢失；二是共享参数问题，在进化过程中，主要是采用共享参数 ${\sigma }_{share}$ 来维持解群体的分布性（基于共享机制的小生境技术），但共享参数的大小不容易确定， 参数的动态修改和调整更是一件困难的工作；三是构造 Pareto 最优解集（通常是构造进化群体的非支配集）的时间复杂度高，为 $O(r{N}^3)$ （这里 $r$ 为目标数，$N$ 为进化群体的规模），因为每一代进化都需要构造非支配集，这样一来，当进化群体规模较大时，算法执行的时间开销就很大。为此，Deb 等于 2000 年在 NSGA 的基础上，提出了 NSGAII (Deb et al, 2000)。

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

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

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

\quad\quad 设群体 $Pop$ 的规模大小为 $N$，将群体 $Pop$ 按某种策略进行分类排序为 $m$ 个子集 $P_1、P_2、...、P_m$， 且满足下列性质：

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

\quad\quad ② ∀ i , j ∈ { 1 , 2 , . . . , m } \forall i, j \in \{1, 2, ..., m\} ∀i,j∈{1,2,...,m} 且 $i\ne j，P_i\cap P_j=\varnothing$ （所有子集不存在相同的个体）

\quad\quad ③ $P_1\succ P_1\succ \cdot \cdot \cdot \succ P_m$，即 $P_{k+1}$ 中的个体直接受 $P_k$ 中个体的支配 $（k=1,2,...,m-1）$。

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

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

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

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

\quad\quad 首先通过一个二重循环计算每个个体的 $n_p$ 和 $s_p$ ，则 P 1 = { q   ∣   n q = 0 ， q ∈ P o p } P_1 = \{q~|~n_q = 0，q \in Pop\} P1​={q ∣ nq​=0，q∈Pop}。 然后依次按方法 P k = { 所有个体  q   ∣   n q − k + 1 = 0 } P_k = \{所有个体~ q~ | ~n_q - k + 1 = 0\} Pk​={所有个体 q ∣ nq​−k+1=0} 求 $P_2、P_3\cdot \cdot \cdot$。

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

\quad\quad 算法 6.2 由两部分组成，第一部分用于计算 $n_i$ 和 $s_i$，并求得 $P_1$， 所需要的时间为 $(r{N}^2)$ ，这里 $r$ 为目标数，$N$ 为进化群体规模；第二部分用于求 $P_2、P_3、\cdot \cdot \cdot 、P_m$， 最坏情况下一个规模为 $N$ 的进化群体有 $N$ 层边界集 (Pareto front)，即 $m=N$ , 此时其时间复杂度为 $O(N^2)$。由此可得，算法 6.2 的总时间复杂度为 $O(r{N}^2)十O(N^2)$，即为 $O(r{N}^2)$。

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

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

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

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

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

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

$$P[i{]}_{distance}=\sum _{k=1}^r(P[i+1].f_k-P[i-1].f_k)$$

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

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

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

\quad\quad 3. Deb 的 NSGA-II

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

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

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

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

\quad\quad 定义 6.1 \quad 设个体 $i$ 和 $j$ 为进化群体中的任意个体，个体之间的偏序(partial order) 关系 ${\succ }_n$ 为

$$i{\succ }_{n}jif(i_{rank}<j_{rank})or(i_{rank}=j_{rank})and(P[i{]}_{distance}>P[j{]}_{distance})$$

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

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

\quad\quad 算法 6.4 中，通过 F = nondominated-sort($R_t$) 产生了若干个分类子集 $F=(F_1,F_2,\cdot \cdot \cdot )$，但被选入新群体的只有一部分。如图 6.7 所示，分类子集 $F_1$ 和 $F_2$ 中的所有个体均被选入了新群体 $P_{t+1}$，但分类子集 $F_3$ 中只有一部分个体被选入新群体 $P_{t+1}$。 一般地，若 $|F_1|+|F_2|+\cdot \cdot \cdot +|F_{i-1}|\le N$ 且 $|F_1|+|F_2|+\cdot \cdot \cdot +|F_i|>N$，则称 $F_i$ 为临界层分类子集，图 6.7 中的 $F_3$ 为临界层分类子集。

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

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

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

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

\quad\quad 由此可得算法 6.4 的总时间复杂度为 $O(r{N}^2)$，其中主要的时间开销花费在构造边界集上，因此一个快速的构造分类子集（或构造非支配集）的方法有利于提高 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)
1
2
3

代码实现

main.m

function main()

% UserInput.Probleme = problem;
% UserInput.ObjectNum = dimension;

[UserInput] = GetUserInput()
switch (UserInput.Probleme)
	case 1
		AlgoName = 'DTLZ1';
	case 2
		AlgoName = 'DTLZ2';
	case 3
		AlgoName = 'DTLZ3';
	case 4
		AlgoName = 'DTLZ4';
end

		NSGAII(AlgoName, UserInput.ObjectNum);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

NSGAII.m

function NSGAII(Problem,M)
% Function: NSGAII(Problem,M)
%
% Description: 执行 NSGAII 算法的所有逻辑
%
%
% Syntax:
%   
%
% Parameters:
%   Problem：选择测试函数，例如：'DTLZ1',注意输出字符
%   M ： 目标函数个数
%
% Return:
%   无
%
%                Young99
%         Revision:1.0    Data: 2022-12-06
%*************************************************************************

clc;
format compact;% 空格紧凑
tic; % 记录运行时间


% 参数设定
Generations = 700;   % 迭代次数
if M == 2  % 目标数量
    N = 100; % 种群大小
else M == 3
    N = 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

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: 2022-12-06
%*************************************************************************


% input:N 种群规模 
    persistent K; % 定义持久变量
    % 持久变量是声明它们的函数的局部变量；
    % 但其值保留在对该函数的各次调用所使用的内存中。
    
    Boundary = NaN; 
    switch Operation  % 选择操作模式 0/1
        
        % 0：初始化种群
        case 0
            k = 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 + 1
            MaxValue   = 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 1
            Population    = Input;  % 已经初始化完成的种群
            FunctionValue = zeros(size(Population,1),M);
            switch Problem
                case '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 > 1
                            FunctionValue(:,i) = FunctionValue(:,i).*(1-Population(:,M-i+1));
                        end
                    end
                case 'DTLZ2'
                    g = sum((Population(:,M:end)-0.5).^2,2);
                    for i = 1 : M
                        FunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);
                        if i > 1
                            FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));
                        end
                    end
                case 'DTLZ3'
                    g = 100*(K+sum((Population(:,M:end)-0.5).^2-cos(20.*pi.*(Population(:,M:end)-0.5)),2));
                    for i = 1 : M
                        FunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);
                        if i > 1
                            FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));
                        end
                    end
                case 'DTLZ4'
                    Population(:,1:M-1) = Population(:,1:M-1).^100;
                    g = sum((Population(:,M:end)-0.5).^2,2);
                    for i = 1 : M
                        FunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*Population(:,1:M-i)),2);
                        if i > 1
                            FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*Population(:,M-i+1));
                        end
                    end
                case 'DTLZ5'
                    g = sum((Population(:,M:end)-0.5).^2,2);
                    
                    for i = 1 : M
                        theta =(pi./(4*(1+g))).*(1+2.*g.*Population(:,1:M-i));
                        FunctionValue(:,i) = (1+g).*prod(cos(0.5.*pi.*theta),2);
                        if i > 1
                            FunctionValue(:,i) = FunctionValue(:,i).*sin(0.5.*pi.*(pi./(4*(1+g))).*(1+2.*g.*Population(:,M-i+1)));
                        end
                    end                  
            
                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];
                end
            Output = FunctionValue;
         end
            
    end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

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: 2022-12-06
%*************************************************************************

    if Operation == 1
        Kind = 2; 
    else
        Kind = 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，第二轮个体为 P2


    while (Kind==1 && sum(cz) FunctionValue(j,m)  %比较函数值，判断个体的支配关系
                                k = 0;
                                break;
                            end
                        end
                        if k == 1
                            d(j) = 1; % 将所有被 i 支配的都标记为 1
                        end
                    end
                end
                FrontValue(Rank(i)) = MaxFront;
                cz(i) = 1; % 将当前的个体标记为 1，表示已经确定所在的非支配层级
            end
        end
    end
end


%% 非支配排序伪代码
% 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

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)-1
                CrowdDistance(Front(Rank(j))) = CrowdDistance(Front(Rank(j)))+(FunctionValue(Front(Rank(j+1)),i)-FunctionValue(Front(Rank(j-1)),i))...
                    /(Fmax(i)-Fmin(i)); % 除以最大值 - 最小值，即对其标准化处理
            end
        end
    end
end


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

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: 2022-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 : 2
            if Pointer >= N
                Rank = randperm(N);
                Pointer = 1;
            end

            p = 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;
            else
                k(j) = q;
            end
            Pointer = Pointer + 2;
        end
        MatingPool(i,:) = Population(k(1),:);
        MatingPool(i+1,:) = Population(k(2),:);
    end
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

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 > N
        MaxOffspring = N;
    end
    
            % 遗传操作参数
            ProC = 1;       % 交叉概率
            ProM = 1/D;     % 变异概率
            DisC = 20;   	% 交叉参数
            DisM = 20;   	% 变异参数
 
            % 模拟二进制交叉
            Offspring = zeros(N,D);
            for i = 1 : 2 : N
                beta = 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;
            end
            Offspring = Offspring(1:MaxOffspring,:); % 用于选取自带个体数 MaxOffspring

            % 多项式变异
            if MaxOffspring == 1
                MaxValue = Boundary(1,:);
                MinValue = Boundary(2,:);
            else
                MaxValue = repmat(Boundary(1,:),MaxOffspring,1);
                MinValue = repmat(Boundary(2,:),MaxOffspring,1);
            end
            k    = 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