无模型自适应预测控制与迭代学习控制仿真：MATLAB对比验证平台

1. 项目概述

1.1 为什么写这个仿真程序

先交代一下背景。我在做过程控制相关的研究时，经常要验证各种新型控制算法，但每次都要从零手写仿真环境，真是够折腾的。后来索性整理了一套基于 MATLAB 的无模型自适应预测控制（MFAPC）和迭代学习控制（MFAILC）的数值验证仿真程序，专门用来对比这两类算法在非线性、强耦合、时变系统中的表现。

这套程序解决的核心痛点是：没有精确数学模型时，控制器该怎么设计？传统的 PID 或者基于模型的控制（如 MPC）都需要先建模型，但很多工业过程（比如化工反应釜温度控制、电机转速控制、无人机姿态控制）本身就很难建准模型，或者模型参数随工况漂移。无模型自适应控制（MFAC）的思路是"在线估计系统动态的伪偏导数，再用这个估计值设计控制器"，完全不依赖被控对象的结构信息，只依赖输入输出数据。预测控制则是在这个基础上引入多步预测和滚动优化，提高控制品质。而迭代学习控制（ILC）则是针对重复运行的批次过程，利用历史批次数据修正当前批次的控制信号，让跟踪误差逐批减小。

我把两类算法放在同一个仿真框架里，方便横向对比：同样一个非线性被控对象，用 MFAPC 和 MFAILC 分别跑，看跟踪效果、抗扰能力、超调量、控制能量消耗等指标。程序支持用户自定义被控对象模型、参考轨迹、噪声类型和控制器参数，可以说是一个"开箱即用"的验证平台。

1.2 这套程序适合谁用

如果你是以下角色，这套程序能帮你省掉大量重复造轮子的时间：

• 控制理论与控制工程专业的研究生，正在做无模型自适应控制相关课题，需要快速出仿真对比数据；
• 刚接触 MFAC 理论、想通过仿真加深理解的本科生，需要可复现的代码和清晰的参数注释；
• 企业 R&D 工程师，要评估 MFAPC 或 MFAILC 在具体产线上的可行性，先用仿真离线验证控制效果；
• 喜欢自己折腾算法、写代码验证新 idea 的工科爱好者。

程序的设计目标不是追求理论上的极致创新，而是提供一个结构清晰、参数可调、注释完整、结果可复现的仿真实验平台。再复杂的算法，没有一组能跑的仿真代码配合，说服力也会打折。

1.3 核心思路一句话总结

在被控对象模型未知（或仅知道输入输出数据）的前提下，MFAPC 通过对系统动态的在线估计 + 多步预测 + 滚动优化来实现高精度跟踪；MFAILC 则通过对批次误差的学习修正来实现重复轨迹的高精度跟踪。两者都只需要输入输出数据，但工作方式和适用场景明显不同。这套程序让你在同一个被控对象上直接对比它们的差异。

2. 顶层设计与关键模块拆解

2.1 程序总架构：把"控制器"和"被控对象"彻底解耦

写仿真程序最容易犯的错是"算法和对象混在一起，测一个变量就得动全局"。我的做法是把整体框架分成三个独立模块：

第一个模块是被控对象模块。你可以把任意非线性离散系统写成以下形式：

[ y(k+1) = f(y(k), y(k-1), ..., u(k), u(k-1), ...) + d(k) ]

其中 ( d(k) ) 是外部扰动。程序里默认给出三个典型对象：

• 非线性单入单出（SISO）系统：( y(k+1) = \frac{y(k)}{1+y(k)^2} + u(k)^3 + d(k) )；
• 带纯滞后的系统：( y(k+1) = 0.6y(k) - 0.2y(k-1) + 0.5u(k-2) + 0.3u(k-3) )；
• 时变参数系统：系数随 ( k ) 缓慢变化。

第二个模块是控制器模块。MFAPC 和 MFAILC 都实现为独立的函数，输入是历史输入输出数据、参考轨迹和控制器参数，输出是当前时刻的控制量。这样你可以换对象不换控制器，也可以固定对象去调不同算法参数。

第三个模块是仿真主脚本（main）。负责初始化、循环迭代、数据存储、结果绘图。整个过程无需修改算法实现，只需要改配置参数块。

模块解耦有一个实实在在的好处：调试时能单独验证"这个控制量到底有没有被正确计算"，也能快速对比不同算法在同一轨迹、同一噪声水平下的表现。如果你有新的控制算法，只要按照接口规范写一个u = myController(history, ref, params)，就能无缝接入现有框架，不用改其他任何代码。

2.2 MFAPC 控制器：从模型到实现的关键步骤

MFAPC（无模型自适应预测控制）本质上是在无模型自适应控制（MFAC）的基础上引入了预测控制的"多步预测 + 滚动优化"思想。纯 MFAC 用的是当前时刻的估计值来计算当前控制量，模式是"看到现在、控制现在"；而 MFAPC 会利用当前已知数据预测未来 ( N ) 步的输出轨迹，再在这 ( N ) 步窗口内求解优化问题，让预测输出尽量贴合参考轨迹，是一种"看到未来、控制现在"的模式。

核心公式有两部分。

第一部分是伪偏导数（Pseudo Partial Derivative, PPD）估计。对于 SISO 系统，引入一个时变参数 ( \phi_c(k) )，使它满足：

[ \Delta y(k+1) = \phi_c(k) \cdot \Delta u(k) ]

这个 ( \phi_c(k) ) 不是被控对象的真实梯度，而是"数据驱动意义下"的等价梯度。它由当前和历史输入输出数据在线估计得到：

[ \hat{\phi}_c(k) = \hat{\phi}_c(k-1) + \frac{\eta \Delta u(k-1)}{\mu + |\Delta u(k-1)|^2} \left[ \Delta y(k) - \hat{\phi}_c(k-1)\Delta u(k-1) \right] ]

其中 ( \eta ) 是步长因子（通常取 0.5~1），( \mu ) 是权重因子（防止分母为零，并限制 PPD 变化速度，通常取 0.001~0.01）。还要对 PPD 做重置机制：如果 ( |\hat{\phi}_c(k)| \leq \varepsilon ) 或 ( |\Delta u(k-1)| ) 过小，就将 PPD 重置为初始值，防止估计发散。

第二部分是预测控制律。设预测步数为 ( N_p )，控制步数为 ( N_c )（( N_c \leq N_p )）。在当前时刻 ( k )，根据 PPD 估计值可以预测未来输出：

[ \hat{y}(k+i) = y(k) + \sum_{j=1}^{i} \hat{\phi}_c(k) \Delta u(k+j-1) \quad (i=1,2,...,N_p) ]

优化目标为：

[ J(k) = \sum_{i=1}^{N_p} \left[ y^*(k+i) - \hat{y}(k+i) \right]^2 + \lambda \sum_{j=0}^{N_c-1} \Delta u^2(k+j) ]

其中 ( y^*(k+i) ) 是未来参考轨迹，( \lambda ) 是控制量变化惩罚因子。这里要注意：因为在预测模型里 PPD 用的是当前时刻的估计值，并把未来窗口内的 PPD 视为恒定，所以当系统时变剧烈时，预测模型会有偏差。解决办法是让 PPD 估计器跑得足够快（适当加大 ( \eta )），同时控制步数 ( N_c ) 不要取太大。

上述优化问题是一个无约束的二次规划问题，可以解析求解。在单入单出且 ( N_c=1 ) 时，控制律退化为：

[ \Delta u(k) = \frac{\hat{\phi}_c(k) \left[ y^*(k+1) - y(k) \right]}{\lambda + \hat{\phi}_c(k)^2} ]

这正是 MFAC 的标准形式。所以 MFAPC 是 MFAC 的推广，当预测步数为 1 时就回到 MFAC。这也是我在程序里同时实现两者的原因——只需要改一个参数Np，你就能看到预测控制带来的平滑效果。

2.3 MFAILC 控制器：批次维度上的学习修正

MFAILC（无模型自适应迭代学习控制）处理的场景和 MFAPC 完全不同。它假设被控对象在一条固定的参考轨迹上重复运行多次（比如焊接机器人每次焊接同一焊缝、注塑机每次成型同一产品）。每个批次称为一次"迭代"，控制目标是让跟踪误差随着迭代次数增加而减小。

在批次 ( i )，时刻 ( k )，系统动态为（以 SISO 为例）：

[ y_i(k+1) = f(y_i(k), ..., u_i(k), ...) + d_i(k) ]

MFAILC 的核心思想是：不需要知道 ( f ) 的准确形式，而是定义一个"伪偏导数和输入差分的乘积"来等价描述系统动态。具体做法是对系统进行沿迭代轴的线性化：

[ y_{i+1}(k+1) - y_i(k+1) = \phi(k) \cdot \left[ u_{i+1}(k) - u_i(k) \right] ]

注意这里的下标是批次索引 ( i )，时刻是 ( k )。也就是说，沿着批次轴方向，输出变化与输入变化成线性关系，( \phi(k) ) 就是该时刻的伪偏导数，它与批次无关（或缓慢变化）。

控制律的设计目标是：让 ( y_{i+1}(k+1) ) 更接近参考轨迹 ( y_d(k+1) )。令输出误差：

[ e_i(k+1) = y_d(k+1) - y_i(k+1) ]

控制目标是 ( e_{i+1}(k+1) ) 尽量小。由此得到学习律：

[ u_{i+1}(k) = u_i(k) + \rho(k) \cdot \hat{\phi}(k) \cdot e_i(k+1) ]

其中 ( \hat{\phi}(k) ) 是通过历史批次数据估计出来的，( \rho(k) ) 是步长因子，用于控制学习收敛速度。这是一个非常简洁的形式，类似于"沿迭代轴的 P 控制器"，但增益不是人为设定的，而是来自数据驱动估计。

伪偏导数 ( \phi(k) ) 的估计可以通过前两个批次的输入输出数据直接算：

[ \hat{\phi}(k) = \frac{\Delta y_{i-1}(k+1)}{\Delta u_{i-1}(k)} ]

其中 ( \Delta y_{i-1}(k+1)=y_i(k+1)-y_{i-1}(k+1) )，( \Delta u_{i-1}(k)=u_i(k)-u_{i-1}(k) )。如果分母接近零，就保持上一批估计结果不变，防止数值爆炸。

这里的关键判断是：MFAPC 面向的是"连续运行、每次都不同"的场景，而 MFAILC 面向的是"批次重复、每次轨迹相同"的场景。如果你把 MFAPC 用在重复轨道上，它也能工作，因为可以参考轨迹在线跟踪；但如果你把 MFAILC 用在新轨迹上，它就失效了，因为它依赖历史批次的信息。程序里我引入了"轨迹重复次数"这个参数：设成 1，就是单次轨迹跟踪，适合 MFAPC；设成比如 20，就是 20 批次重复，适合 MFAILC。

2.4 为什么不用深度学习或传统系统辨识

也许有人会问：既然模型未知，为什么要绕一圈用 PPD 估计，而不是用神经网络直接拟合系统？或者先做系统辨识再上 MPC？

我的体会是，这是一个"成本和收益"的权衡问题。

神经网络需要大量训练数据，在线实时部署时还需要不断重训练，对于工业控制场合，样本少、安全性要求高、算力受限，神经网络并不总是合适。系统辨识要设计激励信号（如 PRBS），辨识过程本身会干扰生产过程，而且线性模型辨识得到的结果往往是非线性系统的一个局部近似，模型失配后控制品质会明显下降。

MFAC 家族的思路属于"在线等价动态估计"，计算复杂度非常低。PPD 是一个标量（SISO 情况），估计它只需要加减乘除和一次除法，实时性远超神经网络；PDD 估计不存在"局部最优"问题，它就是一个梯度类估计，参数整定也相对直观。这是我选择这两类算法做验证平台的原因——它们最适合工程现场。

3. 实操过程与核心环节实现

3.1 环境准备与代码结构

开发环境我建议统一用 MATLAB R2020b 及以上版本，不需要额外工具箱，纯手写脚本实现。如果偏好开源环境，可以移植到 GNU Octave，绝大多数函数都兼容。以下是我的项目目录结构：

MFAPC_MFAILC_Sim/ ├── main.m # 主仿真脚本，选择对象和控制器 ├── plant/ │ ├── plant_nonlinear.m # 非线性对象 │ ├── plant_delay.m # 纯滞后对象 │ ├── plant_timevar.m # 时变对象 │ └── plant_setup.m # 初始化对象参数 ├── controllers/ │ ├── mfapc_controller.m # MFAPC 控制器 │ ├── mfailc_controller.m# MFAILC 控制器 │ └── common_utils.m # 公共函数：轨迹生成、PPD估计 ├── config/ │ └── config_params.m # 所有可调参数集中管理 ├── results/ │ ├── plot_results.m # 绘图脚本 │ └── compare_algo.m # 性能指标对比 └── README.md

我不建议把参数散落在各个函数里，那样调参时找半天。所有参数集中在一个config_params.m结构体里，例如：

% config_params.m params.sim_time = 500; % 仿真步数 params.Ts = 0.01; % 采样周期（秒） params.ref_type = 'sine'; % 参考轨迹类型：'step','sine','square','multi' params.ref_amp = 1.0; params.ref_freq = 0.1; % 正弦频率（Hz） params.noise_amp = 0.01; % 输出噪声幅值 params.Np = 5; % 预测步数（MFAPC） params.Nc = 2; % 控制步数（MFAPC） params.lambda = 0.01; % 控制增量惩罚 params.eta = 0.8; % PPD估计步长 params.mu = 0.01; % PPD估计分母权重 params.rho = 0.4; % MFAILC学习步长 params.batch_num = 20; % MFAILC迭代批次数 params.batch_len = 500; % 每批次采样点数

这种做法的好处是一眼能看到全部关键参数，改起来也方便，不用担心改了A处漏了B处。

3.2 MFAPC 核心代码实现与参数选择

MFAPC 控制器函数的核心结构如下（简化版）：

function [u, phi_hat] = mfapc_controller(history, ref_hist, params) % history: 结构体，包含 y(1:end), u(1:end-1) % ref_hist: 参考轨迹历史数组 % 返回当前控制量 u 和 PPD 估计值 persistent phi_prev; % 注意：在循环中调用时用 persistent 或外部传递 if isempty(phi_prev) phi_prev = params.phi_init; end y = history.y(end); y_prev = history.y(end-1); u_prev = history.u(end-1); Delta_u = u_prev - history.u(end-2); Delta_y = y - y_prev; % 1. 伪偏导数估计 denom = params.mu + Delta_u^2; phi_new = phi_prev + params.eta * Delta_u / denom * (Delta_y - phi_prev * Delta_u); % 重置机制 if abs(phi_new) < params.phi_min || abs(Delta_u) < 1e-5 phi_new = params.phi_init; end phi_hat = phi_new; % 2. 预测模型与控制律（以 Nc=1 为例） % 实际程序中这里需要根据 Np 展开预测窗口，这里展示标量形式 e = ref_hist(end) - y; Delta_u_now = phi_hat * e / (params.lambda + phi_hat^2); u = u_prev + Delta_u_now; phi_prev = phi_new; end

在完整实现中，当 ( Np>1 ) 时，需要构造一个向量化的预测方程。我建议用递归数组方式构建预测输出，而不是用符号工具箱。以下是我实现的predict_output.m的核心片段：

function y_pred = predict_output(y_current, phi_hat, delta_u_seq) % delta_u_seq: 未来控制增量序列，长度为 Np % y_pred: 未来输出预测 N = length(delta_u_seq); y_pred = zeros(1,N); y_pred(1) = y_current + phi_hat * delta_u_seq(1); for i = 2:N y_pred(i) = y_pred(i-1) + phi_hat * delta_u_seq(i); end end

求解最优控制序列时，可以用二次规划函数quadprog，但因为是简单的稀疏二次函数，我更喜欢直接解析解。对于 ( Np=5, Nc=2 )，可以写出 5 个预测方程，代入 2 个控制增量，求偏导并令其为零，得到一个二元一次方程组，直接矩阵求逆即可。这样做速度快，也不依赖优化工具箱。

参数选择经验（以下参数基于默认被控对象的调试结果，供参考）：

第一次调试时先设 ( Np=1, Nc=1 )，这样算法退化成 MFAC，参数好调；调到跟踪稳定后，再逐步增加 ( Np )，观察控制量是否变得更平滑，超调是否减小。这样由简入繁，能少走很多弯路。

3.3 MFAILC 代码实现与批次循环逻辑

MFAILC 的仿真结构跟 MFAPC 完全不同。它不是在一个 time loop 里跑完所有步数，而是外层迭代批次，内层迭代时间：

% main_MFAILC.m params = config_params(); % 初始化存储 u_log = zeros(params.batch_num, params.batch_len); y_log = zeros(params.batch_num, params.batch_len); e_log = zeros(params.batch_num, params.batch_len); phi_log = zeros(params.batch_len, 1); % 注意：phi与批次无关，但要随批次更新 % 初始控制量和输出设为 0 u_log(1, :) = 0; for k = 1:params.batch_len-1 y_log(1, k+1) = plant_nonlinear(y_log(1, k), u_log(1, k)); end e_log(1, :) = ref - y_log(1, :); for trial = 2:params.batch_num % 当前批次初始状态 u_log(trial, 1) = u_log(trial-1, 1); y_log(trial, 1) = y_log(trial-1, 1); for k = 1:params.batch_len-1 % 1. 先用上一批的数据更新 phi(k) if k > 1 && abs(u_log(trial-1, k) - u_log(trial-1, k-1)) > 1e-8 phi_hat = (y_log(trial-1, k+1) - y_log(trial-1, k)) / ... (u_log(trial-1, k) - u_log(trial-1, k-1)); else phi_hat = phi_log(k); end % 引入投影修正，限制 phi 在合理区间 if abs(phi_hat) > params.phi_max phi_hat = sign(phi_hat) * params.phi_max; end phi_log(k) = phi_hat; % 2. 计算本批次 k 时刻的控制量（沿批次轴学习） if trial == 2 % 第一批学习前，先用上一批的控制量加一个小的学习增量 u_log(trial, k) = u_log(trial-1, k) + params.rho * e_log(trial-1, k+1); else u_log(trial, k) = u_log(trial-1, k) + params.rho * phi_log(k) * e_log(trial-1, k+1); end % 3. 把控制量输入对象得到当前输出 y_log(trial, k+1) = plant_nonlinear(y_log(trial, k), u_log(trial, k)); e_log(trial, k+1) = ref(k+1) - y_log(trial, k+1); end end

这段代码的关键点在于( \phi(k) ) 不是实时估计的，而是通过前一批次的数据事后估计的。所以第一轮迭代（trial=2）学习律只用误差，不用 ( \phi )，相当于先跑一个粗调；从 trial=3 开始，用前两批的数据估计出 ( \phi(k) )，再加速学习。这非常类似人类学骑自行车——先随机摸索，积累了几次经验后开始有意识调整转向幅度。

有一个非常隐蔽但常见的坑：当输出轨迹含有噪声时，直接用差分方式估计 ( \phi(k) ) 会被噪声严重污染。我建议先对 ( y ) 做移动平均滤波，再计算差分，或者改用递推最小二乘估计 ( \phi(k) )。在程序中我增加了一个开关use_filter，开启时会先执行：

y_smooth = movmean(y_log(trial-1, :), 5);

然后用平滑后的数据进行差分。实测表明，在 0.01 噪声幅值下，不滤波时 MFAILC 到第 10 批才开始收敛，滤波后第 4 批就基本收敛了。

3.4 被控对象设计：如何让对比更公平

为了公平对比 MFAPC 和 MFAILC，我给两个控制器设置了相同的被控对象和参考轨迹。但在被控对象上做了点文章：

function y_next = plant_nonlinear(y, u) % 非线性 + 扰动 y_next = 0.8 * y * sin(y) + 0.5 * u^3 + 0.2 * u; end

这个对象故意设计了两个特点：一是 ( y ) 有非线性自激励（( y\sin y )），如果只靠线性控制增益，容易在小范围振荡；二是 ( u^3 ) 引入了强非线性，控制量越大，需要的控制灵敏度就越高，很考验算法的自适应能力。

另外我还在对象中加入了时变增益：

function y_next = plant_timevar(y, u, k) a = 0.6 + 0.2 * sin(2 * pi * k / 200); y_next = a * y - 0.1 * y^2 + u + 0.05 * u^2; end

这样可以让 PPD 随时间变化，考验 MFAPC 的 PPD 估计器能否跟得上。而 MFAILC 在这种情况下还能通过迭代修正弥补，两者差异会非常明显——MFAPC 跟踪误差会随参数波动而起伏，而 MFAILC 在批次维度上会持续收敛。

3.5 性能指标与对比可视化

只凭曲线"看起来差不多"是不够的，我设置了几个量化指标，主脚本运行完会直接在命令行输出：

=================== 性能对比 =================== 算法 RMSE MaxAbsErr TV(u) 批次收敛率 MFAPC 0.042 0.117 0.538 - MFAILC 0.019 0.046 0.271 第6批收敛 ================================================

其中：

• RMSE（均方根误差）：( \sqrt{\frac{1}{N}\sum_{k=1}^{N} e(k)^2} )；
• MaxAbsErr（最大绝对误差）；
• TV(u)（控制器总变化）：( \sum |u(k)-u(k-1)| )，体现控制动作的平滑程度；
• 批次收敛率：MFAILC 中达到基准误差（比如首次批次的 10% 以内）所需的迭代次数。

绘图部分，我分了三个子图：

• 子图一：参考轨迹和实际输出的跟踪对比；
• 子图二：控制量变化过程；
• 子图三：PPD 估计值变化过程。

MFAPC 和 MFAILC 的结果用不同颜色叠加在同一坐标系里，方便直接比较。我在plot_results.m里还加了一个"动态轨迹生成视频"的选项（save_movie=true），可以把每批次输出画成动画，直观地看到 MFAILC 的误差随着批次收缩，这个视觉冲击力远超静态图。

4. 常见问题与排查技巧实录

4.1 仿真发散怎么办：PPD 重置机制和参数调整顺序

第一类典型问题是"程序跑起来就是发散"，通常表现为输出曲线瞬间到几千甚至无穷大。我排查这类问题有个固定顺序：

先查 PPD。MFAPC 的发散几乎都源于 ( \hat{\phi}(k) ) 估计不准或不受控。最有效的防御机制是"强制区间限制"：

phi_hat = max(params.phi_min, min(params.phi_max, phi_hat));

我一般把 ( \phi_{min}=0.01 )，( \phi_{max}=5 )。如果系统是逆符号（比如加热系统中增大控制量反而降低输出），还要允许符号变化，也就是取 ( \phi_{min}=-5, \phi_{max}=-0.01 \cup [0.01,5] )，这需要根据被控对象特性预设。如果一个系统在运行中符号会反转，PPD 估计会特别困难，需要缩小 ( \eta )，同时让重置机制更灵敏。

其次是控制量惩罚 ( \lambda )。反馈回路中，如果 ( \lambda ) 太小，控制增量会非常大，导致被控对象输入饱和。这就引出一个很实用的经验：先调 PPD 重置与区间，再调 lambda，最后调 step length。如果一上来就调 ( \eta )，很容易陷入"加了 ( \eta ) 效果好一点，但加了又发散"的恶性循环。

4.2 MFAILC 在特定批次突然发散

MFAILC 的一个常见病是：前面几批收敛得很好，到某一批突然发散。原因通常是学习律中的 ( \hat{\phi}(k) ) 估计到了病态值。

想象一下，当 ( \Delta u_{i-1}(k) ) 非常接近零时（比如系统已经收敛得很好，控制量不再变化），直接用差分公式算 ( \phi ) 会得到极大值。这时候乘以误差再放大，控制量就会突然跳变。解决办法是加一个阈值判断：

if abs(du_prev) < 1e-4 phi_hat = phi_log(k); % 保持上一批估计值 else phi_hat = dy_prev / du_prev; end

这是"保持-重置"策略。另外，为了抑制振荡，我给学习律加了一个简单的低通滤波：

u_log(trial, k) = 0.7 * (u_log(trial-1, k) + params.rho * phi_log(k) * e_log(trial-1, k+1)) + ... 0.3 * u_log(trial, k-1);

注意这个低通是在"批次方向"上混合还是"时间方向"上混合？我这里是时间方向的平滑，可以让控制量更温和，但会稍微牺牲收敛速度。如果发散现象严重，我建议先这样做，再把rho调小到 0.2 左右。

4.3 两个算法的参数整定对比

我发现很多初学者会把 MFAPC 的参数照搬到 MFAILC，这正是踩坑的关键。两者的参数作用维度完全不同：

我在 GUI 调试面板上特意把两个控制器的参数区隔开来，避免混淆。有一次我一个师弟把 MFAPC 的 ( \lambda=0.01 ) 当成了 MFAILC 的 ( \rho )，结果第一轮学习幅度很小，到第十批才开始收敛，他以为是收敛慢，其实是参数理解错了。

4.4 测量噪声导致 PPD 估计抖动

如果在输出端加入了高斯噪声（幅度 0.01），PPD 估计值会像心电图一样上下乱跳。即使控制效果还行，但绘图时很不好看，更重要的是它会让控制量变躁。

我的解决方案有三层：

1. 在数据进控制器前做低通滤波（如filter函数或移动平均），但注意会增加相位延迟，对需要快速响应的系统不友好；
2. 在 PPD 估计式中增加衰减因子，让估计值“惯性”大一点，不是每次都用最新误差更新。具体做法是把eta变成时变的，从 0.9 逐步衰减到 0.2；
3. 在控制律中加入死区：当误差小于阈值时，不更新控制量，防止微小的噪声反复扰动执行器。

三层结合，实测下来噪声抑制效果非常好，RMSE 在噪声环境下能改善 30% 以上。

4.5 仿真速度优化

纯 MATLAB 仿真如果批次多、步长密，循环嵌套会非常慢。我之前跑 100 批、每批 1000 步、做了 20 次蒙特卡洛实验，一共花了 25 分钟。后来把内层循环向量化了，时间降到 4 分钟。

向量化的主要思路是：在被控对象可用向量化表达式时（例如线性部分），直接对整个数组用数组操作代替for；对非线性对象，可以先用parfor（Parallel Computing Toolbox）并行跑不同批次，因为每个批次的计算互相独立。我推荐优先用parfor，改动成本最低。另外，不要在循环内部用persistent声明变量并频繁读写工作区，优先用预分配的数组存储，避免动态数组增长。

4.6 结果不可复现问题

如果程序中用了随机数生成噪声，每次跑的结果可能不同。为了让实验可复现，我固定随机种子：

rng(20240415); % 固定随机种子

或者在配置参数中加一个params.seed。这个细节看似不起眼，但在写论文时非常重要，审稿人让你提供"具体某一次实验的数据"，你必须能重新生成一模一样的结果。不要觉得随机种子无关紧要，它真的救过我的论文。

5. 实用扩展与场景适配

5.1 从 SISO 到 MIMO 的扩展思路

我最初程序只支持 SISO，但后来很多场景（如双输入双输出温度湿度系统）需要 MIMO。扩展的难点在于 PPD 从标量变成了矩阵。理论上的 MIMO-MFAC 要用到"伪分块雅可比矩阵"，实现比较复杂。我在程序里采用了一种更简单的解耦方案：假设输入输出配对已知，把 MIMO 系统分解为多个 SISO 回路分别控制，只在被控对象层面加入耦合项。这种工程近似做法虽然不能从理论上保证最优，但很多实际系统在弱耦合条件下效果足够好。

如果你想深入研究 MIMO-MFAC 理论，可以改造 PPD 估计模块，把标量phi_hat改成矩阵Phi_hat，学习律变成矩阵形式：

[ \Delta u(k) = \left[ \lambda I + \Phi(k)^T \Phi(k) \right]^{-1} \Phi(k)^T e(k) ]

这是标准的最小二乘形式，很多论文里都这么写。但要注意矩阵求逆的数值稳定性，建议加正则化项。

5.2 与 PID 的对比实验设计

为了体现无模型自适应控制的优势，我还在框架里预留了增量式 PID 控制器作为 baseline：

function u = pid_controller(history, ref, params) e = ref - history.y(end); de = e - history.e_prev; u = history.u_prev + params.Kp*(e - history.e_prev) + ... params.Ki*history.e_int + params.Kd*de; end

建议你在对比实验中设置这样几个场景：

• 场景一：对象参数固定，三个控制器（PID、MFAPC、MFAILC）都整定到最佳状态，比较在阶跃响应下的超调量和调节时间；
• 场景二：被控对象增益在实验中途突变（比如在第 200 步把增益从 1.0 调到 1.5），观察谁受的影响小；
• 场景三：加入输出噪声，观察谁的控制量更平滑。

通常会发现：PID 在场景一里表现很好，但场景二里会有明显退化；MFAPC 在场景二和三里表现稳定；MFAILC 在重复轨迹场景三（如果是批次重复）表现最优。这样一个三脚猫的对比，已经足够说明"无模型自适应控制"的价值所在。

5.3 拓展到实际工业现场的关键点

仿真只是第一步，实际落地还要注意几个问题：

输入约束：真实执行器有幅值饱和与速率限制。仿真程序里可以简单加限幅：

u = max(u_min, min(u_max, u));

但直接裁剪会导致 PPD 估计失真，因为模型假定的是未限幅的控制序列。解决办法是记录实际实施的控制量（限幅后）作为历史数据，而控制量计算使用限幅前的值时要谨慎。更优的方案是把幅值约束纳入优化问题，变成一个带约束的二次规划，但这对实时性有考验。

采样时间匹配：MFAPC 的控制律是离散形式，采样周期Ts必须大于控制器计算时间。如果控制器计算耗时为 1 ms，而采样周期设为 0.1 ms，实际执行时会丢步或抖动。我建议采样周期至少是控制器计算耗时的 5 倍以上。

故障处理：实际系统的传感器可能断线，输出数据会跳变到 0 或不可信。仿真程序里我加了一个fault_inject()函数，可以注入暂时性传感器故障，测试算法的鲁棒性。MFAC 对传感器故障是比较脆弱的，因为 PPD 估计依赖输出反馈，一旦输出跳变，PPD 会瞬间失真，需要额外的故障诊断与容错机制。

5.4 与其他控制算法的混合与切换

在实际系统中，我经常使用一种"双模控制"策略：当误差较大时，用 MFAPC 跑快速跟踪；当误差较小，等于进入稳态微调阶段时，切换成定增益 PID 省去不必要的计算量。或者反过来：当输出接近目标时，PID 容易产生振荡，这时切换到 MFAPC，利用它的在线估计特性保持稳定。程序里我预留了switch_controller()接口，你可以设定误差阈值来自动切换。

这个思路其实源自工程中常见的"启停分离"思想。你可以用仿真程序预先测试切换阈值，观察切换瞬间是否有冲击。我这里有一个经验：切换瞬间，控制器输出的初始值要设置为前一个控制器最后的输出值，否则会有一个阶跃跳跃。程序里加了这个初始化，能够显著减少切换震荡。

6. 现场调试实录与经验总结

6.1 一次调参的完整心路历程

拿 MFAPC 调试默认非线性对象来说，我的第一版参数是：( \eta=0.5, \mu=0.01, \lambda=0.01, Np=3, Nc=1 )。初次运行，发现输出经过约 50 步后能跟踪上正弦参考，但每个波峰处都有明显的 10% 超调，控制量曲线不光滑，且 PPD 估计值波动很大。

我的第一个怀疑点是：( \lambda ) 太小，控制增量惩罚不足，导致控制动作太"猛"。于是把 ( \lambda ) 从 0.01 调到 0.1，超调确实降到了 4%，但跟踪速度变慢，前 30 步误差明显增大。接着我把 ( \eta ) 从 0.5 提到 0.9，PPD 估计加快，跟踪延迟缩小，但控制量有点抖动。然后我把 ( Np ) 从 3 增到 5，这时控制量明显变平滑，超调维持在 4%，但 RMSE 几乎不变。

最后微调 ( \mu ) 从 0.01 到 0.005，PPD 追踪更快，但整体变化不大。最终参数是：( \eta=0.9, \mu=0.005, \lambda=0.05, Np=5, Nc=2 )。结果：RMSE 0.031，最大超调 2.8%，控制量平滑度比初版提升了 40%。这个调参过程大概花了 40 分钟，大部分时间耗在理解"哪个参数影响哪个指标"上。如果你没耐心手动调，可以写一个简单的网格搜索脚本，但建议你先手动调一遍找到感觉，再交给脚本去精细化搜索。

6.2 MFAILC 的学习曲线与批次收敛规律

观察 MFAILC 的批次收敛曲线有一个很典型的模式：从第一批到第三批，误差下降非常快（50% 以上）；从第三批到第六批，下降变缓；第六批以后，误差几乎不变，进入平台期。这个平台期的误差水平取决于系统噪声下限，不能无限制降低。

如果你发现批次误差平台期太高，比如在 0.1 左右下不去，首先检查是不是 PPD 估计噪声太大。可以关闭噪声（noise_amp=0）再测试，如果平台期明显降低，就是噪声导致的；如果完全没有变化，则需要检查参考轨迹是否真的可被系统实现——比如参考轨迹的带宽高于对象的物理允许范围，任何学习算法都追不上。

另一个规律是，步长 ( \rho ) 越大，初期收敛越快，但后期更容易在平台期附近震荡。我一般采取"递减步长"策略：

rho_k = params.rho * (1 - 0.5 * trial / params.batch_num);

前期大步长快速逼近，后期小步长精细收敛。实测能在保证收敛的同时把平台期震荡降低 70%。

6.3 容易出现却容易被忽视的逻辑错误

以下这几个 bug，我在自己写代码和帮别人改代码时见过无数次，哪怕是很熟练的程序员也容易掉坑：

第一个是索引越界。在 MFAILC 中，( k+1 ) 时刻参考轨迹与输出对应，循环里经常出现访问倒数第二个元素后再访问最后一个，一不小心就数组越界。建议先在脚本第一行加assert(size(ref, 2) == batch_len)，防止隐性错误。

第二个是历史数据初始化错误。MFAPC 在 ( k=1 ) 时没有 ( u(0) )、( y(0) )，需要自定义初始值。我一般设 ( u(0)=0, y(0)=0 )，但有些系统初始状态不是零，比如输出被强制到某个工作点，这会跟 PPD 估计打架。所以在plant_setup.m里特意把初始状态作为参数提供，你要根据实际系统去改。

第三个是控制量和输出命名搞混。因为程序里同时有u_log和y_log，矩阵索引又都是(trial, k)，很容易把u_log(trial, k)写成y_log(trial, k)。一旦写错，结果看起来会非常离谱。我建议在变量命名中用前缀区分：cmd_u表示控制量，sys_y表示系统输出。

第四个是偏移量错误。由于实际被控对象的输出响应有惯性，控制量是在 ( k ) 时刻作用，而输出在 ( k+1 ) 时刻才能看到。如果把 ( k ) 时刻的控制量与 ( k ) 时刻的输出对齐，看起来会有一步延迟，但控制律不会错。很多人误以为这是算法问题，实际上是数据对齐问题。

6.4 蒙特卡洛实验与鲁棒性验证

在论文或者项目报告中，光有确定性实验不够，通常要做蒙特卡洛实验。程序里我提供了monte_carlo.m，输入参数范围，自动生成多组随机参数（比如噪声幅度在 0~0.05 之间随机，被控对象参数在 ±20% 范围内随机），每组跑 50 次，统计 RMSE 的均值、标准差、P95 值。

这样做有一个很实际的收益：能够量化说明"MFAPC 比 PID 在抗参数失配方面平均提升 XX%，在 95% 置信区间内……"。审稿人喜欢看到这种有统计意义的数据，工程决策者也更认可这种鲁棒性证据。

有一点要提醒：蒙特卡洛实验必须确保所有算法都在同一组随机参数下运行，否则公平性存疑。我的做法是先生成一个随机种子池，每个实验场景存一份种子编号，不同算法跑相同的种子顺序。

7. 最后再聊几句心里话

在我的项目中做这个仿真程序的过程，有几点体会很想分享给后来者。

第一个体会是，无模型自适应控制这一类算法的核心价值不在于理论多漂亮，而在于"它给你一种不依赖模型但有理论保障的控制设计路径"。当你在现场被一个复杂的非线性对象搞得焦头烂额时，MFAPC 和 MFAILC 能成为你工具箱里几把比较趁手的工具。仿真程序最大的作用是逼你把算法细节理解到"能用代码写出来"的程度，这个过程中你会看透很多公式之外的工程细节。

第二个体会是，代码组织能力和控制理论功底同样重要。我见过很多研究者的控制算法理论上完美，但仿真代码逻辑混乱，无法复现，最后自己都忘了当初怎么调的参数。我写这套程序的初衷之一，就是希望即便三年后再打开文件夹，我还能在三分钟内重新跑通全部实验并看懂每个参数的含义。为了实现这个目标，我在代码里写了大量中文注释，并且把参数集中管理。这是一个费工夫但长期受益的习惯，建议你也试试。

第三个体会是，仿真验证不是终点，而是起点。跑通仿真后一定要问自己：如果被控对象变化了，算法还 work 吗？如果噪声大了，算法会崩溃吗？如果换一个工程对象，要怎么调参数？顺着这个思路去扩展你的实验，你收获的远比搞定一个仿真要多。

我后续的计划是把这套程序扩展成一个支持拖拽式 GUI 的验证平台，让不会写代码的工程师也能直观地调节参数、观察曲线、导出报告。目前已经完成了基础框架，等测试稳定后再专门写一篇文章分享。在这之前，建议你先下载代码，自己动手把默认对象跑一遍，改几个参数看看曲线怎么变化，这个直接感受过程比看我文字描述有用得多。

如果调试过程中碰到具体问题，按照文章里的排查顺序基本都能解决。实在不行，返回去看你的 PPD 估计函数，八成问题出在那里。