关于NSGA-II的介绍可以参考这篇文章

一、伪代码

初始化种群 P0
评估种群目标函数值
非支配排序和拥挤度计算

while (未达到最大代数) do
    基于非支配排序和拥挤度选择父代
    交叉和变异生成子代 Qt
    评估子代目标函数值
    
    合并种群 Rt = Pt ∪ Qt
    非支配排序和拥挤度计算
    选择前 N 个个体作为下一代种群 Pt+1
end while

输出Pareto前沿

二、完整代码

% NSGA-II 多目标优化示例
% 流程

% 步骤〇、初始化种群
% 评估种群目标函数值
% 非支配排序和拥挤度计算
% 
% while (未达到最大代数) do
%     基于非支配排序和拥挤度选择父代
%     交叉和变异生成子代 Qt
%     评估子代目标函数值
% 
%     合并种群 Rt = Pt ∪ Qt
%     非支配排序和拥挤度计算
%     选择前 N 个个体作为下一代种群 Pt+1
% end while
% 
% 输出Pareto前沿



clc;clear;
% 定义参数
pop_size = 50; % 种群大小
max_gen = 100; % 最大代数
num_obj = 2; % 目标函数数量
num_var = 2; % 决策变量数量
lb = [-5, -5]; % 下边界
ub = [5, 5]; % 上边界

% 初始化种群
pop = initializePopulation(pop_size, num_var, lb, ub); % 生成初始种群
pop_obj = evaluatePopulation(pop, num_obj); % 评估种群目标值

% 主循环
for gen = 1:max_gen % 遍历每一代

    % 步骤一、采用二元锦标赛选择选出优质父代
    mating_pool = tournamentSelection(pop, pop_obj, pop_size);

    % 步骤二、生成子代（随机选择2个父代，交叉，变异为2个子代），并求适应度
    offspring = generateOffspring(mating_pool, lb, ub); 
    offspring_obj = evaluatePopulation(offspring, num_obj); 
    
    % 步骤三、合并父代和子代的 个体和适应度
    combined_pop = [pop; offspring]; 
    combined_obj = [pop_obj; offspring_obj];
    
    % 步骤四、进行快速非支配排序和拥挤度计算
    [fronts, crowding_distance] = fastNonDominatedSorting(combined_obj);
    
    % 步骤五、根据快速非支配排序和拥挤度结果，选择pop_size个下一代
    [pop, pop_obj] = selectNextGeneration(combined_pop, combined_obj, fronts, crowding_distance, pop_size); % 选择下一代
    
end

% 绘制Pareto前沿
figure; % 创建新图形窗口
plot(pop_obj(:, 1), pop_obj(:, 2), 'ro'); % 绘制Pareto前沿
xlabel('目标1'); % 设置x轴标签
ylabel('目标2'); % 设置y轴标签
title('Pareto 前沿'); % 设置图形标题
grid on; % 打开网格

% 输出最优解
disp('最优解：'); % 显示最优解
disp(pop); % 输出种群
disp('对应的目标值：'); % 显示对应的目标值
disp(pop_obj); % 输出目标值

% 初始化种群函数
function pop = initializePopulation(pop_size, num_var, lb, ub)
    pop = repmat(lb, pop_size, 1) + repmat(ub - lb, pop_size, 1) .* rand(pop_size, num_var); % 生成大小为pop_size的种群，每个个体有num_var个变量，变量值在lb和ub之间
end

% 评估种群函数
function obj = evaluatePopulation(pop, num_obj)
    obj = zeros(size(pop, 1), num_obj); % 初始化目标值矩阵
    for i = 1:size(pop, 1) % 对每个个体计算目标值
        obj(i, :) = myMultiObjectiveFcn(pop(i, :)); % 计算个体的目标值
    end
end

% 定义多目标函数
function f = myMultiObjectiveFcn(x)
    f1 = x(1)^2 + x(2)^2; % 目标函数1
    f2 = (x(1) - 1)^2 + x(2)^2; % 目标函数2
    f = [f1, f2]; % 返回目标值
end

% 选择父代（二元锦标赛选择算法）
% 生成一些适应度较高的个体。具体来说，通过在当前种群中进行比较和选择，
% 将那些在非支配排序和拥挤度方面表现较好的个体挑选出来，形成一个新的交配池。
function mating_pool = tournamentSelection(pop, obj, pool_size)
    mating_pool = zeros(pool_size, size(pop, 2)); 
    for i = 1:pool_size
        a = randi(size(pop, 1));
        b = randi(size(pop, 1));
        if dominates(obj(a, :), obj(b, :)) || rand < 0.5
            winner = a;
        else
            winner = b;
        end
        mating_pool(i, :) = pop(winner, :);
    end
end


% 生成后代函数
function offspring = generateOffspring(pop, lb, ub)
    pop_size = size(pop, 1); % 获取种群大小
    num_var = size(pop, 2); % 获取变量数量
    offspring = zeros(pop_size, num_var); % 初始化后代矩阵
    for i = 1:2:pop_size % 生成后代
        p1 = pop(randi([1 pop_size]), :); % 随机选择父母个体1
        p2 = pop(randi([1 pop_size]), :); % 随机选择父母个体2
        [c1, c2] = crossover(p1, p2); % 交叉操作生成两个子代
        c1 = mutate(c1, lb, ub); % 变异操作子代1
        c2 = mutate(c2, lb, ub); % 变异操作子代2
        offspring(i, :) = c1; % 将子代1加入后代种群
        offspring(i+1, :) = c2; % 将子代2加入后代种群
    end
end

% 交叉操作函数
function [c1, c2] = crossover(p1, p2)
    alpha = rand(size(p1)); % 生成随机权重alpha
    c1 = alpha .* p1 + (1 - alpha) .* p2; % 生成子代1
    c2 = alpha .* p2 + (1 - alpha) .* p1; % 生成子代2
end

% 变异操作函数
function mutated = mutate(ind, lb, ub)
    mutation_rate = 0.1; % 变异率
    for i = 1:length(ind) % 对个体的每个变量进行变异
        if rand < mutation_rate % 判断是否进行变异
            ind(i) = lb(i) + (ub(i) - lb(i)) * rand; % 生成新的变量值
        end
    end
    mutated = ind; % 返回变异后的个体
end

% 非支配排序函数
function [fronts, crowding_distance] = fastNonDominatedSorting(obj)
    % obj是父代+子代合并的100个个体
    num_obj = size(obj, 2); % 获取目标函数数量
    pop_size = size(obj, 1); % 获取种群大小
    fronts = {}; % 初始化 fronts
    dom_count = zeros(pop_size, 1); % 初始化每个个体被支配的数量
    dom_set = cell(pop_size, 1); % 初始化支配集合
    fronts{1} = []; % 初始化 fronts{1}
    
    % 步骤一、生成第一个前沿
    for i = 1:pop_size % 计算支配关系
        for j = 1:pop_size % 遍历所有个体
            if dominates(obj(i, :), obj(j, :)) % 判断个体i是否支配个体j
                dom_set{i} = [dom_set{i}, j]; % 将j加入i的支配集合
            elseif dominates(obj(j, :), obj(i, :)) % 判断个体j是否支配个体i
                dom_count(i) = dom_count(i) + 1; % 增加i被支配的数量
            end
        end
        if dom_count(i) == 0 % 如果个体i不被任何个体支配
            fronts{1} = [fronts{1}, i]; % 将其放入第一前沿
        end
    end
    
    % 步骤二、生成剩余前沿
    k = 1; % 初始化前沿索引
    while ~isempty(fronts{k}) 
        next_front = []; % 初始化下一前沿
        for i = fronts{k} % 遍历当前前沿的每个个体
            for j = dom_set{i} % 遍历个体i的支配集合
                dom_count(j) = dom_count(j) - 1; % 减少j被支配的数量
                if dom_count(j) == 0 % 如果j不再被任何个体支配
                    next_front = [next_front, j]; % 将j加入下一前沿
                end
            end
        end
        k = k + 1; % 增加前沿索引
        fronts{k} = next_front; % 设置下一前沿
    end
    
    % 步骤三、计算每个个体在每个目标函数下的拥挤度（与旁边解的距离）
    % 【区别】这是NSGA-II新增的，NSGA没有这一步
    crowding_distance = zeros(pop_size, 1); % 初始化拥挤度
    for k = 1:length(fronts) % 计算每个前沿的拥挤度
        front = fronts{k}; % 获取当前前沿
        if length(front) > 2 % 如果前沿包含多个个体
            for m = 1:num_obj % 遍历每个目标函数
                [~, sorted_idx] = sort(obj(front, m)); % 按目标函数值排序
                crowding_distance(front(sorted_idx(1))) = inf; % 设置边界个体的拥挤度为无穷大
                crowding_distance(front(sorted_idx(end))) = inf; % 设置边界个体的拥挤度为无穷大
                for i = 2:length(front) - 1 % 计算内部个体的拥挤度
                    crowding_distance(front(sorted_idx(i))) = crowding_distance(front(sorted_idx(i))) + ...
                        (obj(front(sorted_idx(i+1)), m) - obj(front(sorted_idx(i-1)), m)) / ...
                        (max(obj(:, m)) - min(obj(:, m))); % 计算拥挤度
                end
            end
        else
            crowding_distance(front) = inf; % 如果前沿只有一个个体，设置其拥挤度为无穷大
        end
    end
end

% 判断支配关系函数
function is_dominated = dominates(p, q)
    % all(p <= q)：如果p在所有目标都不劣于q
    % any(p < q)：如果p在至少一个目标上强于q（确保两个解不是完全相同）
    is_dominated = all(p <= q) && any(p < q); % 判断个体p是否支配个体q
end

% 选择下一代函数
function [next_pop, next_obj] = selectNextGeneration(pop, obj, fronts, crowding_distance, pop_size)
    next_pop = []; % 初始化下一代种群
    next_obj = []; % 初始化下一代目标值
    i = 1; % 初始化前沿索引

    % 步骤一、优先选择 前沿靠前的非支配解
    while size(next_pop, 1) + length(fronts{i}) <= pop_size % 按前沿逐个添加个体，直到达到种群大小
        next_pop = [next_pop; pop(fronts{i}, :)]; % 添加当前前沿的个体
        next_obj = [next_obj; obj(fronts{i}, :)]; % 添加当前前沿的目标值
        i = i + 1; % 增加前沿索引
    end


    % 步骤二、如果没有选满，剩余的位置通过拥挤度排序来填补
    % 具体来说，假如pop_size=100，
    %   情况一、前两个前沿分别有40和60个个体，那么正好填满pop_size，不需要拥挤度排序
    %   情况二、前两个前沿分别有40和70个个体，总数110超过pop_size，需要采取算法从70个中选出60个
    % 【区别】关于怎么选这60个，NSGA是直接取前60个；NSGA-II是根据拥挤度排序，优先选拥挤度高的，确保解的多样性
    if ~isempty(fronts{i}) 

        % 以下是NSGA-II
        [~, idx] = sort(crowding_distance(fronts{i}), 'descend'); % 按拥挤度排序，优先选距离远的
        remaining = pop_size - size(next_pop, 1); % 计算剩余数量
        next_pop = [next_pop; pop(fronts{i}(idx(1:remaining)), :)]; % 添加拥挤度最高的个体
        next_obj = [next_obj; obj(fronts{i}(idx(1:remaining)), :)]; % 添加拥挤度最高的目标值

        % % 以下是NSGA
        % % 1）NSGA 没有拥挤度的概念，也就没有拥挤度排序（精英策略）
        % % 2）NSGA直接根据前沿里个体顺序取前几个。NSGA-II优先选择拥挤度高的个体(比较分散的)，确保解的多样性；
        % remaining = pop_size - size(next_pop, 1); % 计算剩余需要添加的个体数量
        % next_pop = [next_pop; pop(fronts{i}(1:remaining), :)]; % 添加当前前沿中的前remaining个个体
        % next_obj = [next_obj; obj(fronts{i}(1:remaining), :)]; % 添加当前前沿中的前remaining个个体的目标值
    end
end

三、几个简单问题

1、为什么要合并父代和子代？对比GA

• 选择机制：

  • GA：传统遗传算法通过选择操作从当前代选择个体，然后通过交叉和变异生成下一代。每一代仅保留通过选择操作选出的个体。
  • NSGA-II：通过合并父代和子代，进行快速非支配排序和拥挤度计算，选择出下一代种群。保证最优解和多样性同时得到保留。

• 精英保留策略：

  • GA：通常不包含显式的精英保留策略，可能会导致最优解在进化过程中丢失。
  • NSGA-II：通过合并父代和子代，确保每一代的最优解能够被保留，并在后续迭代中继续优化。

• 多目标优化：

  • GA：主要用于单目标优化，选择机制较为简单，通常使用适应度值进行选择。
  • NSGA-II：专为多目标优化设计，通过非支配排序和拥挤度计算处理多个目标函数，确保优化结果在多个目标函数上同时优化。

2、什么是拥挤度？为什么要计算拥挤度？

目的是选稀疏的个体。这是因为选稀疏的个体有助于保持种群的多样性，防止过早收敛到局部最优解，并确保帕累托前沿的均匀分布。

1）计算公式

2）拥挤度高低区别

• 高拥挤度：

  • 如果一个个体 $i$ 的拥挤度 $d(i)$ 高，则意味着它与其相邻个体之间的距离较大，即它位于一个相对稀疏的区域。
  • 在NSGA-II的选择过程中，这些高拥挤度的个体会被优先选择，因为它们代表了目标空间中稀疏区域，有助于保持种群的多样性。

• 低拥挤度：

  • 如果一个个体 $i$ 的拥挤度 $d(i)$ 低，则意味着它与其相邻个体之间的距离较小，即它位于一个相对稠密的区域。
  • 这些低拥挤度的个体在选择过程中可能会被舍弃，以保证种群的多样性和覆盖面。

3、为什么通过二元锦标赛选出父代？

• 简单高效：

  • 实现简单，只需要随机选择和比较两个个体即可，不需要复杂的计算。
  • 选择过程快速，适用于大规模种群。

• 平衡选择压力：

  • 二元锦标赛选择能够平衡选择压力，避免过强或过弱的选择压力。
  • 选择压力适中，既能保证优秀个体有较高的选择概率，又能保持种群多样性。

• 适应度比例选择：

  • 通过比较两个个体的适应度（或非支配排序等级和拥挤度），可以根据个体的相对优劣进行选择。
  • 优秀个体在多次锦标赛中获胜的概率更高，从而增加了它们的选择概率。

• 保持多样性：

  • 由于随机选择和比较的过程，二元锦标赛选择能够在一定程度上保持种群的多样性。
  • 防止种群过早收敛到局部最优解。

4、二元锦标赛会选出重复个体？

二元锦标赛选择可能会选择到重复的个体，这是进化算法中的一个正常现象。

在某些情况下，重复的个体可以帮助保持优良的基因在种群中的比例，特别是当这些个体在适应度上表现较好时。

为了减少重复个体对种群多样性的负面影响，可以在生成子代时引入更多的变异操作，确保生成的子代具有多样性。

5、选择下一代是否会用不到拥挤度？

存在这种功能，假如pop_size=100

• 情况一、前两个前沿分别有40和60个个体，那么正好填满pop_size，不需要拥挤度排序
• 情况二、前两个前沿分别有40和70个个体，总数110超过pop_size，第一个前沿40个会加入下一代，第二个前沿70个需要通过拥挤度排序，选出拥挤度高的

6、怎么理解 is_dominated = all(p <= q) && any(p < q)

• all(p <= q)：如果p在所有目标都不劣于q，即 小于或等于
• any(p < q)：如果p在至少一个目标上强于q（确保两个解不是完全相等，即相同）

同时满足这2点，才能说p支配q

...待补充