多目标优化（一）简单的 NSGA-Ⅱ

多目标优化（一）简单的 NSGA-Ⅱ

写在前面：

1、接触到的第一个多目标优化算法，对于2~3个目标真挺有效的；
2、本文主要介绍自己理解的 NSGA-II Matlab 程序实现，在学习中迷糊的地方都会点出来；
3、本文程序基于 [link] https://blog.csdn.net/qq_40434430/article/details/82876572 改编（觉得太罗嗦了）；
4、欢迎指正。

算法概要：

1、基础 Pareto 概念请转：
2、算法目标：我们知道遗传算法会构建一个带有 NP 个个体的种群，单目标输出的就是迭代后种群中最好的 个体。但在多目标中，输出的是迭代后整个种群作为 NP 个 Pareto 最优解，构成一个 Pareto 实际前沿面。（单目标是要找最优个体，多目标是要找最接近理论前沿面的实际前沿面）
3、快速非支配排序：将当代种群根据各个目标函数值进行 Pareto 等级分层（用到非支配解的概念）
4、拥挤度计算：不同层的解分别计算拥挤度（就是解之间间隔越大，拥挤度越大，是我们更中意的解）
5、选择、交叉、变异：GA基本操作，NSGA-II论文并没有限定选用哪种方法，我觉得可以随意选用
6、精英保留：就是将更新的候选种群和原种群合并，重新 Pareto 分层和计算拥挤度，然后按照 Pareto 层级和拥挤度选择前 NP 个个体进入下一循环。

主程序：

clear; close all; clc
global NP iteration f_num x_num pc pm yita1 yita2 k % 定义全局变量，写子程序就不用那么麻烦了
% NP：种群大小
% iteration：最大迭代次数
% f_num：目标个数
% x_num：解得维数
% pc、pm：交叉、变异概率
% yita1、yita2、k：GA更新方式公式用到的参数

%% 测试ZDT1函数，参数如下：
f_num = 2; % 目标函数个数
x_num = 30; % 决策变量个数
x_min = zeros(1,x_num); % 决策变量的最小值
x_max = ones(1,x_num); % 决策变量的最大值

NP = 100;
iteration = 500;
k = ceil(NP / 2);
yita1 = 20; yita2 = 20;
pc = 0.9; pm = 1 / x_num;
    
% 初始化种群
chromo_ori = initialize(); % 参数都是全局变量
[F0, chromo] = non_domination_sort(chromo_ori);  % 非支配排序 种群个数 NP
chromo = crowding_distance_sort(F0, chromo);  % 计算拥挤度 种群个数 NP

% 开始迭代
for i = 1:iteration

chromo_parent = chromo(:,1:x_num+f_num+2); % 定义父代种群
% 选择 交叉 变异
pick_sel = selection(chromo_parent); % 这里输出的是选中的个体索引（我用来观察都选那些个体了）
chromo_cro = crossover(pick_sel, chromo_parent); 
chromo_off = mutation(chromo_cro); % GA更新得到的种群定义为子代种群
    
% 种群合并
chromo_co = [chromo_parent(:,1:x_num+f_num); chromo_off(:,1:x_num+f_num)]; % 种群个数 2*NP
[F_i, chromo_co] = non_domination_sort(chromo_co); % 非支配排序
chromo_co = crowding_distance_sort(F_i, chromo_co);% 计算拥挤度
% 精英保留
chromo = elitism(chromo_co); % 按 Pareto 等级和拥挤度从2*NP个个体中选出NP个个体进入下一代
end
 
%% 可视化 （这个不重要）
chromo = sortrows(chromo, x_num+1);
plot(chromo(:, x_num+1), chromo(:,x_num+2),'*');

end

子程序原理：

1、initialize：在约束范围内随机生成 NP * x_num+f_num 的矩阵，每一行是一个解，后面几列是各目标函数值。

function chromo = initialize()
    
    global NP f_num x_num x_min x_max % 需要用到的参数
    
    chromo = zeros(NP, x_num+f_num); 
    for i = 1:NP
        for j = 1:x_num
            chromo(i, j) = x_min(j) + (x_max(j) - x_min(j)) * rand;
        end
        chromo(i, x_num+1:x_num+f_num) = object_fun(chromo(i,1:x_num)); % 计算各目标函数值放在后面几列
    end
    
end

2、object_fun：计算各目标函数值，f 返回的是一个行向量

function f = object_fun(x)
        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

3、non_domination_sort：思想如下，但程序操作不是这样做的

① 输入一个种群，输出一个带 Pareto 分层后的种群，层级标注放在种群最后一列
② 找到非支配个体，Pareto 层级标注为1，去掉层级为1的个体（少了这些，有些个体按规则就是非支配个体了呀）
③ 找到非支配个体，Pareto 层级标注为2，去掉层级为2的个体
④ ……
⑤ 直到整个种群 Pareto 层级标注完

function [F, chromo] = non_domination_sort(chromo)
    
    global f_num x_num
    NP = size(chromo, 1); % 在合并种群后种群大小实际 2倍的，所以这里不用全局变量
   
    pareto_rank = 1;  % 初始化 pareto 等级为 1
    F(pareto_rank).ss = []; % 这是一个结构体，可以看成是 F 是一个向量，向量元素是一个向量（之前都是用元组实现的，学习了）
    
    for i = 1:NP
        
        % p 也是一个结构体（用用就知道是怎么回事了）
        p(i).n = 0; % 记录个体 i 被多少个个体支配了
        p(i).s = []; % 记录个体 i 支配了哪些个体
        
        for j = 1:NP
            gov = 0; equal = 0; gov_ed = 0;
            for k = 1:f_num 
                if chromo(i, x_num+k) < chromo(j, x_num+k) % 个体 i 与其他个体在第 k 个目标上的适应度值比较
                    gov = gov + 1; 
                elseif chromo(i, x_num+k) == chromo(j, x_num+k)
                    equal = equal + 1; 
                else
                    gov_ed = gov_ed + 1; 
                end
            end
            
            if gov == 0 && equal ~= f_num % 表示 i 被个体 j 支配了，计数 +1
                p(i).n = p(i).n + 1; 
            elseif gov_ed == 0 && equal ~= f_num % 表示个体 i 支配了个体 j，将个体记录下来
                p(i).s = [p(i).s j]; 
            end
        end
        
        if p(i).n == 0 % 表示没有人能支配 i 呀，那 i 就是非支配个体
            chromo(i, f_num+x_num+1) = 1; % 在种群函数值后标记 Pareto 层级
            F(1).ss = [F(1).ss i]; % 统计层级为 1 的个体
        end
        
    end % Pareto 层级为 1 的统计完啦
 
    while ~isempty(F(pareto_rank).ss) % F(pareto_rank).ss 非空，执行
        temp = [];
        len = length(F(pareto_rank).ss); % 该层级有几个个体
        for i = 1:len
            if ~isempty(p(F(pareto_rank).ss(i)).s) % 个体 i 有支配个体，执行
                ind = F(pareto_rank).ss(i); % 第 i 各非支配解的索引（对应p(i).n=0）
                len_i = length(p(ind).s); % 非支配解 ind 支配了几个个体        
                for j = 1:len_i  
                    ind_i = p(ind).s(j); % 第 j 个被个体 ind 支配的个体索引
                    p(ind_i).n = p(ind_i).n - 1; % 个体 ind_i 被支配数-1
                    if p(ind_i).n == 0 % 变成了 非支配解？
                        chromo(ind_i, f_num+x_num+1) = pareto_rank + 1;
                        temp = [temp ind_i]; % 统计 pareto_rank+1 等级的解
                    end
                end
            end
        end
        pareto_rank = pareto_rank + 1; % 接着统计下一层级
        F(pareto_rank).ss = temp; % 记录 第 Pareto_rank 层级包含了哪些个体
    end
end

4、crowding_distance_sort：计算各层级个体拥挤度

① 取出各层级所有个体，按第 i 个目标值进行排序
② 将目标值最小和最大的个体的间隔设置为 无穷大
③ 按公式计算其他个体的间隔
④ 所有目标的间隔相加起来就是个体的拥挤度
⑤ 拥挤度理解请参考：

function chromo = crowding_distance_sort(F, chromo)

    global f_num x_num
    
    temp = sortrows(chromo, f_num+x_num+1); % 按 pareto_rank 排序后的种群
    chromo = []; 
    
    for pareto_rank = 1:(length(F)-1) 
        
        len = length(F(pareto_rank).ss); % 该 pareto_rank 层包含了多少各个体
        y = temp(1:len, :); % 取出该 pareto_rank 层的所有个体。（temp是排过序的）
        temp(1:len, :) = []; % 清除 pareto_rank 层个体（便于后面取下一层级个体）
        
        for i = 1:f_num
            
            [y_sort, index] = sortrows(y, x_num+i); % 按第 i 个函数值对取出的个体排序
            f_min = y_sort(1, x_num+i); % y_sort 是第 Pareto_rank 层按第 i 个目标排序后的种群
            f_max = y_sort(len, x_num+i);
            y(index(1), x_num+f_num+i+1) = inf; % 将目标值最小的个体间隔设为 无穷
            y(index(len), x_num+f_num+i+1) = inf; % 将目标值最大的个体间隔设为 无穷 
            for j = 2:len-1 % 按公式计算其他个体的间隔
                y(index(j), x_num+f_num+i+1) = (y_sort(j+1, x_num+i) - y_sort(j-1, x_num+i)) / (f_max - f_min);
                if isnan(y(j, x_num+f_num+i+1))
                    y(j, x_num+f_num+i+1) = inf;
                end
            end
        end
        chromo = [chromo; y]; % 列数变成了 x_num + f_num + f_num
    end

    % 多个目标函数拥挤度求和
    f_y = sum(chromo(:, (x_num+f_num+2):(x_num+f_num*2+1)), 2); % 间隔相加起来就是拥挤度了
    chromo = [chromo(:, 1:x_num+f_num+1) f_y]; % 清除间隔所在的列数，添加拥挤度。列数变成了 x_num + f_num + 1
    
end

5、selection：

① 网上大部分参考选用的方法是：锦标赛选择（但程序和方法对不上，所以这里我是按锦标赛选择方法重新编写了程序）
② 要从 NP 个父代种群中选出 NP 个个体（较好的个体希望被多选几次）
③ 重复 NP 次，每次从随机挑出的 k 个个体中选最好的。（k 太大的话，每次都包含最好个体，那选出的个体不都是同一个了；k 太小的话，不包含最好的个体，最好的个体没被选中也不合理。所以一般 k=NP/2）
④ 个体好不好首先看 Pareto 层级，层级一样的话就看拥挤度，拥挤度又一样？随机选一个吧。

function pick = selection(chromo)

    global NP f_num x_num k
    pick = []; % 记录选中的个体
    
    for i = 1:NP
        
      	index = randperm(NP, k)'; % 随机选出 k 个个体进行比较
        rank = chromo(index, x_num+f_num+1); % 先比 Pareto 层级
        dis = chromo(index, x_num+f_num+2); % 再比 拥挤度
        idx_rank = find(rank == min(rank));
        idx_dis = find(dis(idx_rank) == max(dis(idx_rank)));
        chose = randperm(length(idx_dis),1); % 比完还是不只一个？随机选一个
        idx = index(idx_rank(idx_dis(chose))); % 找到选中的个体在父代中的索引
        pick = [pick; idx]; % 记录索引比记录个体更好观察呀。（其实都是一样的）
        
    end
    
end

6、crossover：参考的程序将交叉和变异放一起，极度不合理；所以这里我是将交叉和变异分开；而且这里是先判断是否交叉，然后每维向量都交叉（也挺奇怪的，我认为交叉方法是可以自选的，想怎么交叉变异看自己喜好）

① 从selection中得到的种群选出两个不一样的（这种群里有太多一样的了，要限制以下）
② 判断是否要交叉（概率还是挺大的，所以一般都是交叉了），要的话每维都按公式交叉
③ 不管交叉了没都将这两个个体保存起来
④ 重复 NP 次，就保存了 2*NP 个个体了，那就再随机算出 NP 个吧。（所以要不干脆重复 NP/2 次就好了，妥）

function chromo = crossover(pick_sel, chromo)

    global NP f_num x_num x_min x_max pc yita1
    temp = chromo(pick_sel, 1:x_num); % 取出选中的个体（放这步在这里是为了更方便用索引选出两个不同的个体）
    
    off_s = []; % 存放子代个体
    for i = 1:NP
        
        % 选出两个不同的个体交叉
        pick = randperm(NP, 2);
        while pick_sel(pick(1)) == pick_sel(pick(2))
            pick = randperm(NP, 2);
        end
        
        off_1 = temp(pick(1), 1:x_num); off_2 = temp(pick(2), 1:x_num); % 选出了不一样的个体啦
        if rand < pc 
        	% 每一维都按公式交叉以下。。
            for j = 1:x_num
                u = rand;
                if u < 0.5
                    gama = (2 * u)^(1 / (yita1 + 1));
                else
                    gama = (1 / (2 * (1 - u)))^(1 / (yita1+1));
                end
                off_1(j) = 0.5 * ((1 + gama) * off_1(j) + (1 - gama) * off_2(j));
                off_2(j) = 0.5 * ((1 - gama) * off_1(j) + (1 + gama) * off_2(j));
            end
            % 超界了就修改下
            off_1 = min(off_1, x_max); off_1 = max(off_1, x_min);
            off_2 = min(off_2, x_max); off_2 = max(off_2, x_min);
            
        end
        
        off_s = [off_s; off_1; off_2]; % 存放起来
    end
    
    % 从 off_s 中随机选出 NP 个个体
    pick = randperm(2*NP, NP)'; % 从 2*NP 中随机选出 NP 个
    chromo = off_s(pick, :);
    
end

7、mutation：

function chromo = mutation(chromo_cro)

    global NP x_num f_num x_min x_max pm yita2 
    chromo = chromo_cro(:,1:x_num);  % 我忘记 chromo_cro 列数是多少了，反正我只用到 x_num 列
    
    for i = 1:NP
        
        if rand < pm
        	% 每维都要按公式进行变异啊（变异概率是比较小的，没几个会变异）
            for j = 1:x_num
                u = rand;
                if u < 0.5
                    delta = (2 * u)^(1 / (yita2+1)) - 1;
                else
                    delta = 1 - (2 * (1 - u))^(1 / (yita2+1));
                end
                chromo(i,j) = chromo(i,j) + delta;
            end
        end
        % 个体超界了就修改下
        chromo(i,1:x_num) = min(chromo(i,1:x_num), x_max);
        chromo(i,1:x_num) = max(chromo(i,1:x_num), x_min);
        chromo(i,x_num+1:x_num+f_num) = object_fun(chromo(i,1:x_num)); % 顺便计算下各目标函数值吧
        
    end
    
end

8、 elitism：

① 变异后的种群就是子代种群了，和父代种群合并以下吧（变成 2NP 个个体了）
② 将这 2NP 个个体 Pareto 分层和计算拥挤度吧，
③ 先将 Pareto 等级为1的选中，再选 Pareto 等级为2的，……
④ 选 NP 个就够啦，比如 Pareto 等级为1的个体有 NP-1 个， Pareto等级为2的再选一个最好的就好了，按拥挤度选呗，拥挤度最大的不只一个呀，那就随机选一个呗。

function chromo = elitism(chromo_co)
    
    global NP f_num x_num
    
    chromo_co_sort = sortrows(chromo_co, x_num+f_num+1); % 按 pareto_rank 排序
    chromo = [];
    pareto_rank = 1;
    
    while 1 == 1
        
        % 取出 pareto_rank 整层
        temp = chromo_co_sort(chromo_co_sort(:,x_num+f_num+1)==pareto_rank,:);
        if (size(temp, 1) + size(chromo, 1)) < NP % 没放满就一直放
            chromo = [chromo; temp];
        else
            remain = NP - size(chromo, 1);
            temp = sortrows(temp, -(x_num+f_num+2)); % 按拥挤度排序
            temp = temp(1:remain,:);
            chromo = [chromo; temp];
            break; % 放满了就收工
        end
        
        pareto_rank = pareto_rank + 1;
        
    end
    
end

结果比较一下子

实现：每个程序无脑复制保存，打开主程序，运行，这张图就出来了：
红点是理论前沿面，蓝色的是算法求出来的。参数调一下，这两条线很容易就重合了
再给一下测试函数理论前沿面在哪找吧：https://sop.tik.ee.ethz.ch/download/supplementary/testproblems/