NSGA-II算法在MATLAB中求解单目标和多目标函数极值

NSGA-II算法在MATLAB中求解单目标和多目标函数极值

---

一、算法原理与适用场景

NSGA-II（非支配排序遗传算法II） 是解决多目标优化问题的经典算法，通过非支配排序和拥挤度计算维护种群多样性。虽然其设计初衷是多目标优化，但通过调整也可处理单目标问题。以下是其核心特点：

• 多目标优化：同时优化多个冲突目标，生成Pareto最优解集
• 单目标优化：通过添加辅助目标（如约束违反度）转化为多目标问题
• 关键优势：收敛性强、解集分布均匀、适合高维问题

---

二、MATLAB实现框架

1. 基础结构（多目标优化）

function nsga2_demo()
    % 参数设置
    nPop = 100;      % 种群大小
    maxGen = 50;     % 最大迭代次数
    nObj = 2;        % 目标函数维度
    varSize = ; % 决策变量维度
    lb = ; % 下限
    ub = ; % 上限
    
    % 初始化种群
    pop = repmat(struct(), nPop, 1);
    for i = 1:nPop
        pop(i).pos = lb + (ub-lb).*rand(varSize);
        [pop(i).f, pop(i).g] = objFunc(pop(i).pos);
    end
    
    % 主循环
    for gen = 1:maxGen
        % 选择操作（锦标赛选择）
        parents = tournamentSelect(pop, nPop/2);
        
        % 交叉变异
        offspring = geneticOperators(parents, varSize, lb, ub);
        
        % 合并种群
        combined = [pop; offspring];
        
        % 非支配排序
        [fronts, ~] = nonDominatedSort(combined);
        
        % 拥挤度计算
        crowding = calcCrowdingDistance(combined, fronts);
        
        % 选择下一代
        pop = selectNextGen(combined, fronts, crowding, nPop);
    end
    
    % 绘制Pareto前沿
    plotPareto(fronts{1}.f);
end

2. 关键函数实现

(1) 目标函数定义

function [f, g] = objFunc(x)
    % 示例：ZDT1测试函数（多目标）
    f1 = x(1);
    g = 1 + 9*sum(x(2:end))/(length(x)-1);
    f2 = g*(1 - sqrt(f1/g));
    f = [f1, f2];
end

(2) 非支配排序

function [fronts, ranks] = nonDominatedSort(pop)
    n = numel(pop);
    fronts = {};
    ranks = zeros(n,1);
    dominated = zeros(n,1);
    dominatesList = cell(n,1);
    
    for i = 1:n
        for j = 1:n
            if i ~= j
                if dominates(pop(i).f, pop(j).f)
                    dominatesList{i} = [dominatesList{i}, j];
                elseif dominates(pop(j).f, pop(i).f)
                    dominated(i) = dominated(i) + 1;
                end
            end
        end
        if dominated(i) == 0
            fronts{1} = [fronts{1}, i];
            ranks(i) = 1;
        end
    end
    
    k = 1;
    while ~isempty(fronts{k})
        nextFront = [];
        for i = fronts{k}
            for j = dominatesList{i}
                dominated(j) = dominated(j) - 1;
                if dominated(j) == 0
                    nextFront = [nextFront, j];
                end
            end
        end
        k = k + 1;
        fronts{k} = nextFront;
    end
end

function isDominate = dominates(a, b)
    isDominate = all(a <= b) && any(a < b);
end

(3) 拥挤度计算

function crowding = calcCrowdingDistance(combined, fronts)
    crowding = zeros(size(combined));
    for f = fronts
        frontIdx = f{1};
        n = numel(frontIdx);
        if n == 1
            crowding(frontIdx) = Inf;
            continue;
        end
        for m = 1:size(combined(1).f, 2)
            [~, sortIdx] = sort([combined(frontIdx).f(:,m)]);
            crowding(frontIdx(sortIdx(1))) = Inf;
            crowding(frontIdx(sortIdx(end))) = Inf;
            for i = 2:n-1
                crowding(frontIdx(sortIdx(i))) = crowding(frontIdx(sortIdx(i))) + ...
                    (combined(frontIdx(sortIdx(i+1))).f(:,m) - combined(frontIdx(sortIdx(i-1))).f(:,m)) / ...
                    (max(combined.f(:,m)) - min(combined.f(:,m)));
            end
        end
    end
end

---

三、单目标优化改造方案

虽然NSGA-II专为多目标设计，但可通过以下方式处理单目标问题：

1. 伪多目标转化

   添加虚拟目标函数（如约束违反度）：

   function [f, g] = singleObjWrapper(x)
       f1 = objective(x);  % 原目标函数
       g1 = 0;             % 虚拟目标（无约束时）
       f = [f1, g1];
   end
   

2. 参数调整建议 种群大小：增大至200-500 交叉概率：提高至0.9-0.95 变异概率：降低至0.01-0.05

---

四、工程应用案例

案例1：液压系统参数优化

% 目标函数：压力波动最小化 + 能耗最小化
function [f, g] = hydraulicObj(x)
    f1 = std(x(1:3));    % 压力波动
    f2 = sum(x(4:6));    % 能耗
    f = [f1, f2];
end

% 约束条件：工作压力范围
function [c, ceq] = constraints(x)
    c = [x(1) - 200;    % 最小压力
         300 - x(1)];   % 最大压力
    ceq = [];
end

案例2：电力系统经济调度

% 目标函数：发电成本最小化 + 排放最小化
function [f, g] = powerObj(x)
    f1 = sum(x(1:3).*costMatrix);  % 发电成本
    f2 = sum(x(4:6).*emissionMatrix); % 排放
    f = [f1, f2];
end

% 约束条件：功率平衡
function [c, ceq] = constraints(x)
    c = sum(x(1:3)) - 1000;  % 总发电量=1000MW
    ceq = [];
end

---

五、优化

1. 并行计算加速

   parfor i = 1:nPop
       pop(i).f = objFunc(pop(i).pos);
   end
   

2. 自适应参数调整

   function pc = adaptiveCrossover(gen, maxGen)
       pc = 0.8 + 0.1*(gen/maxGen);  % 线性递增
   end
   

3. 混合优化策略

   % NSGA-II + 模拟退火
   if rand < 0.1
       pop = localSearch(pop);  % 局部扰动
   end
   

---

六、结果分析

% 绘制Pareto前沿
figure;
hold on;
colors = hsv(size(fronts{1}.f,1));
for i = 1:size(fronts{1}.f,1)
    plot(fronts{1}.f(i,1), fronts{1}.f(i,2), 'o', 'Color', colors(i,:), 'MarkerSize',10);
end
xlabel('Objective 1');
ylabel('Objective 2');
title('Pareto Front');
grid on;

---

七、参考

1. MATLAB官方工具箱 Global Optimization Toolbox Reinforcement Learning Toolbox
2. 参考代码 NSGA2求解单目标或者多目标函数的极值 www.youwenfan.com/contentcsi/63794.html
3. 进阶文献 Deb, K., et al. "A fast and elitist multiobjective genetic algorithm: NSGA-II." IEEE TEVC 2002. 《多目标优化算法与应用》（胡包钢）