应用多目标遗传算法（NSGA-II）优化PID参数

一、 算法原理与核心思想

在工业控制中，PID控制器的三个参数（$K_p$, $K_i$, $K_d$）直接决定了系统的性能。然而，传统的Ziegler-Nichols整定法或单目标优化往往顾此失彼：响应快了就超调大，超调小了就响应慢。

多目标遗传算法（Multi-Objective Genetic Algorithm, MOGA），特别是NSGA-II (Non-dominated Sorting Genetic Algorithm II)，是解决这一矛盾的终极武器。它不再寻找单一的“最佳”解，而是寻找一组Pareto最优解集。

核心隐喻：
想象你要买一辆车，你有两个目标：速度快和油耗低。显然，这两个目标是冲突的。MOGA就像一个聪明的购车顾问，它会给你展示一系列选择：

• 方案A：速度极快，但油耗高（激进型PID）。
• 方案B：速度一般，油耗极低（保守型PID）。
• 方案C：速度与油耗平衡（中庸型PID）。

这组方案就是Pareto前沿（Pareto Front）。作为决策者，你可以根据当前的生产需求（比如节能模式还是高效模式）从这组解中灵活选择，而不是被一个固定的参数绑架。

---

二、 问题建模：定义战场

要对PID进行多目标优化，首先需要明确我们要优化的战场（目标函数）和士兵（决策变量）。

2.1 决策变量（士兵）

PID的三个参数，它们是有物理边界的：
$X = [K_p, K_i, K_d]$
约束条件：
$K_{p\_min} \le K_p \le K_{p\_max}$
$K_{i\_min} \le K_i \le K_{i\_max}$
$K_{d\_min} \le K_d \le K_{d\_max}$

2.2 目标函数（战场）

我们需要同时最小化以下几个冲突的指标：

1. 超调量 ($f_1$)：越小越好，防止系统震荡。
   $f_1 = \text{Overshoot}(\%)$
2. 调节时间 ($f_2$)：越小越好，系统越快达到稳定。
   $f_2 = \text{Settling Time (s)}$
3. 稳态误差 ($f_3$)：越小越好，控制精度高。
   $f_3 = \text{Steady State Error}$
4. 控制输入能量 ($f_4$)：越小越好，保护执行机构（如阀门、电机），防止剧烈动作。
   $f_4 = \int u(t)^2 dt$

最终目标：找到一组 $[K_p, K_i, K_d]$，使得 $F(X) = [f_1, f_2, f_3, f_4]$ 达到 Pareto 最优。

---

三、 NSGA-II 算法执行流程

初始化种群: 随机生成N组PID参数 [Kp, Ki, Kd]
  |
  v
[快速非支配排序]: 将种群分层 (Rank 1, Rank 2...)
  |- Rank 1: 最优层 (不被任何其他解支配)
  |- Rank 2: 次优层
  |      ...
  |
  v
[拥挤度计算]: 在同一层中，计算解之间的距离
  |- 目的: 保持种群多样性，防止解都挤在一起
  |
  v
[锦标赛选择]: 随机选两个个体PK
  |- 胜者: 等级(Rank)更高 或 拥挤度更大
  |
  v
[交叉与变异]: 模拟生物进化
  |- 交叉: 两个父代PID参数交换基因 (如 Kp1 + Kp2 -> 新Kp)
  |- 变异: 随机微调参数 (如 Kp = Kp * 1.05)
  |
  v
[生成新一代种群]
  |
  v
[终止条件?] --> 是 --> 输出 Pareto 最优解集
  |       |
  |      否
  |       v
  +-------+

---

四、 MATLAB 核心实现代码

MATLAB 提供了强大的 Global Optimization Toolbox，我们可以直接调用 gamultiobj 函数来实现 NSGA-II。

4.1 定义被控对象（以二阶系统为例）

假设我们要控制的系统传递函数为：
$G(s) = \frac{1}{s^2 + 2s + 1}$

4.2 定义多目标评价函数 (evaluate_pid.m)

这是最核心的部分，它接收一个PID参数，运行Simulink仿真或数值计算，返回四个目标值。

function objectives = evaluate_pid(x)
    % x(1) = Kp, x(2) = Ki, x(3) = Kd
    Kp = x(1);
    Ki = x(2);
    Kd = x(3);
    
    % --- 1. 构建PID控制系统 ---
    % 使用数值积分（欧拉法）模拟闭环响应
    dt = 0.01;          % 仿真步长
    T = 10;             % 仿真总时间 10秒
    time = 0:dt:T;
    N = length(time);
    
    % 参考输入（阶跃信号）
    r = ones(1, N) * 1.0; 
    
    % 初始化状态
    y = zeros(1, N);    % 系统输出
    u = zeros(1, N);    % 控制输入
    error = zeros(1, N);
    integral = 0;
    prev_error = 0;
    
    % --- 2. 时域仿真 ---
    for i = 1:N-1
        error(i) = r(i) - y(i);
        integral = integral + error(i) * dt;
        
        % PID 控制律
        derivative = (error(i) - prev_error) / dt;
        u(i) = Kp * error(i) + Ki * integral + Kd * derivative;
        
        % 限制控制量饱和
        u(i) = max(min(u(i), 10), -10); 
        
        % 模拟被控对象 (G(s) = 1/(s^2+2s+1))
        % 使用简单的欧拉法求解微分方程
        % d^2y/dt^2 + 2*dy/dt + y = u
        if i == 1
            y_dot = 0;
            y_ddot = 0;
        else
            y_dot = (y(i) - y(i-1)) / dt;
            y_ddot = (y_dot - prev_y_dot) / dt;
        end
        
        y_ddot = u(i) - 2*y_dot - y(i);
        y_dot = y_dot + y_ddot * dt;
        y(i+1) = y(i) + y_dot * dt;
        
        prev_error = error(i);
        prev_y_dot = y_dot;
    end
    
    % --- 3. 计算性能指标 (目标函数) ---
    % 1. 超调量
    overshoot = max(0, max(y) - 1.0) * 100;
    
    % 2. 调节时间 (2%误差带)
    settling_idx = find(abs(y - 1.0) < 0.02, 1, 'last');
    settling_time = settling_idx * dt;
    if isempty(settling_time)
        settling_time = T; % 未稳定
    end
    
    % 3. 稳态误差
    steady_state_error = abs(y(end) - 1.0);
    
    % 4. 控制能量 (ITSE)
    control_effort = sum(u.^2) * dt;
    
    % 返回目标向量
    objectives = [overshoot, settling_time, steady_state_error, control_effort];
end

4.3 主优化脚本 (main_nsga2.m)

clc; clear; close all;

% --- 1. 设置优化参数 ---
nvars = 3; % 变量个数 (Kp, Ki, Kd)

% 变量边界 (根据实际情况调整)
lb = [0.1, 0.001, 0.01];  % 下限
ub = [100, 10, 10];       % 上限

% --- 2. 设置NSGA-II选项 ---
options = optimoptions('gamultiobj', ...
    'PopulationSize', 100, ...      % 种群规模
    'MaxGenerations', 50, ...       % 最大迭代代数
    'CrossoverFraction', 0.8, ...   % 交叉概率
    'MutationFcn', {@mutationadaptfeasible}, ...
    'PlotFcn', {'gaplotpareto'}, ... % 绘制Pareto前沿
    'Display', 'iter');

% --- 3. 运行多目标遗传算法 ---
fprintf('开始多目标PID参数优化...\n');
[x_pareto, f_pareto] = gamultiobj(@evaluate_pid, nvars, [], [], [], [], lb, ub, options);

% --- 4. 分析结果 ---
% x_pareto: 每一行是一组最优PID参数
% f_pareto: 对应的四个目标值

fprintf('\n优化完成！共找到 %d 组Pareto最优解。\n', size(x_pareto, 1));

% 显示第一组解
best_idx = 1; % 选择第一组解
Kp_opt = x_pareto(best_idx, 1);
Ki_opt = x_pareto(best_idx, 2);
Kd_opt = x_pareto(best_idx, 3);

fprintf('\n推荐的PID参数：\n');
fprintf('Kp = %.4f\n', Kp_opt);
fprintf('Ki = %.4f\n', Ki_opt);
fprintf('Kd = %.4f\n', Kd_opt);

% 绘制Pareto前沿
figure;
scatter3(f_pareto(:,1), f_pareto(:,2), f_pareto(:,3), 'filled');
xlabel('超调量 (%)');
ylabel('调节时间 (s)');
zlabel('稳态误差');
title('PID参数优化的Pareto前沿');
grid on;

参考代码 应用多目标遗传算法优化对PID参数进行多目标优化 www.youwenfan.com/contentcsu/60041.html

五、 结果分析与决策

运行上述代码后，你会得到一个三维（或四维）的散点图，这就是 Pareto前沿。

如何选择最终参数？

1. 保守场景（安全第一）：如果你的系统是锅炉或高压反应釜，选择 超调量 $f_1$ 最小 的解。哪怕响应慢一点，也不能炸。
2. 高效场景（速度第一）：如果是数控机床进给，选择 调节时间 $f_2$ 最小 的解。允许少量超调，换取最快到位。
3. 节能场景（成本第一）：如果是暖通空调，选择 控制能量 $f_4$ 最小 的解。牺牲一点响应速度，换来电费的大幅降低。

实例对比：

方案	$K_p$	$K_i$	$K_d$	超调量	调节时间	评价
传统Z-N	12.5	0.8	0.1	25%	1.2s	震荡大
MOGA-A	8.2	0.2	0.5	0%	2.5s	平稳
MOGA-B	15.1	1.5	0.01	15%	0.8s	快速