我做控制仿真这几年，有一个特别深的体会：很多算法论文写得天花乱坠，但真到自己动手做数值验证的时候，光是把伪偏导数（PPD）的初值调好、把学习增益选对，就够让人折腾一整天。
MFAPC和MFAILC这两个名字，研究数据驱动控制的人应该都不陌生。MFAPC是无模型自适应预测控制，MFAILC是无模型自适应迭代学习控制，两者都出自侯忠生教授提出的无模型自适应控制（MFAC）框架。这套仿真程序就是围绕这两种算法做的数值验证工具，用来在典型非线性被控对象上检验控制效果、对比算法差异、调试控制器参数。
这篇文章我把这套程序的来龙去脉、算法推导、仿真实现、参数整定和踩坑经验全部摊开讲，既有理论层面的逻辑拆解，也有能直接照着跑的代码思路。想搞懂MFAPC和MFAILC到底怎么落地、怎么验证的，不管是写论文需要对比仿真结果，还是做工程预研想看看数据驱动控制的实际表现，这篇文章都能给你省下不少时间。
1. 为什么要写这套仿真程序
1.1 MFAC框架解决的是什么问题
先说个最简单的场景。假设你要控制一个电加热炉的温度，炉子的热惯性、环境温度波动、电网电压变化都会影响输出。你尝试用PID控制，发现负载一变参数就不对了；想用模型预测控制（MPC），但你需要先建立一个足够精确的传热模型，还得在线辨识——这个过程本身就耗费大量精力。
MFAC的思路完全不一样。它不依赖被控对象的数学模型，而是利用被控对象的输入输出数据，在每一个工作点上构造一个等价的动态线性化模型，然后基于这个线性化模型设计控制器。这个思路的核心是伪偏导数（Pseudo Partial Derivative, PPD）的概念——把非线性系统在局部工作点用一个带时变参数的线性模型来逼近，而这个时变参数就是PPD。
PPD不需要精确建模，只需要根据实时输入输出数据在线估计。这就像一个不用看地图也能开车的人，他不需要知道每条路的具体走向，只需要根据当前看到的道路状况，不断修正自己的方向盘角度，照样能把车开到目的地。
1.2 MFAPC和MFAILC各自的定位
MFAPC把MFAC和预测控制的思路结合起来。传统的MFAC控制器，控制律推导时只考虑当前一步的跟踪误差，属于一种"贪心"策略。而MFAPC引入了预测时域的概念，在当前时刻把未来N步的输出预测出来，然后在一个预测时域内优化控制增量序列，只取第一步作用于被控对象，滚动优化。这个思路跟传统MPC类似，但区别在于MPC需要显式的模型，而MFAPC用的是PPD在线估计出来的等价线性化模型。
MFAILC则是针对重复运行过程的。很多工业场景是批次式的，比如注塑成型、间歇反应、机器人重复搬运轨迹跟踪，每个批次从起点跑到终点，过程高度重复。迭代学习控制（ILC）的核心思想就是利用上一次运行产生的误差信号来修正当前批次的控制输入，随着批次增加，跟踪误差逐步收敛。MFAILC把这个思想和MFAC的数据驱动框架结合起来，不需要被控对象的模型，直接利用每次运行的输入输出数据更新PPD估计和控制律。
这两类算法放在一起做数值验证，正好覆盖了两类典型需求：一类是连续运行过程的实时控制，另一类是重复批次过程的逐次改进。仿真程序把它们放到同一个被控对象上跑，可以从收敛速度、稳态精度、抗扰动能力、参数敏感性等方面做横向对比。
2. 核心算法推导与实现逻辑
2.1 紧格式动态线性化与PPD估计
MFAC系列算法的基础是紧格式动态线性化（Compact Form Dynamic Linearization, CFDL）。对于单输入单输出的非线性离散时间系统：
y(k+1) = f(y(k), y(k-1), ..., y(k-ny), u(k), u(k-1), ..., u(k-nu))
在满足一定条件（偏导数连续、广义Lipschitz等）的情况下，可以写成：
Δy(k+1) = φ(k) · Δu(k)
其中φ(k)就是伪偏导数，是一个时变标量。Δy(k+1) = y(k+1) - y(k)，Δu(k) = u(k) - u(k-1)。
这个式子的意义非常直观：在当前工作点附近，系统输出的变化量近似等于PPD乘以输入的变化量。PPD体现了"这个工作点附近的等效增益"，它随工作点变化而变化，通过在线估计来更新。
PPD的估计采用带惩罚项的准则函数：
J = (Δy(k) - φ̂(k)·Δu(k-1))² + μ·(φ̂(k) - φ̂(k-1))²
第一项让估计误差尽可能小，第二项让PPD估计值的变化不要太剧烈。μ是惩罚因子，μ越大，PPD估计越平滑，对噪声的鲁棒性越强；μ太小，PPD估计值可能剧烈跳动，导致控制量抖动。
对φ̂(k)求极值，用梯度法得到PPD估计算法：
φ̂(k) = φ̂(k-1) + (η·Δu(k-1))/(μ + Δu(k-1)²) · (Δy(k) - φ̂(k-1)·Δu(k-1))
这里η是PPD估计的步长因子，取值通常在(0, 1]之间。
有个必须处理的工程细节：当Δu(k-1)接近零的时候，PPD估计会退化。所以复位机制是必不可少的——当|φ̂(k)| ≤ ε 或者 |Δu(k-1)| ≤ ε 的时候，把φ̂(k)重置为一个预设定的初值或者上一个有效估计值。这个细节在仿真中非常关键，很多跑飞的现象都是这里没有处理好。
2.2 MFAPC控制律推导
MFAPC的核心是在预测时域内进行滚动优化。在k时刻，利用当前PPD估计值φ̂(k)，可以预测未来N步的输出：
y(k+1) = y(k) + φ̂(k)·Δu(k) y(k+2) = y(k+1) + φ̂(k)·Δu(k+1) = y(k) + φ̂(k)·(Δu(k) + Δu(k+1)) ... y(k+Nu) = y(k) + φ̂(k)·Σ(i=1到Nu) Δu(k+i-1) ...
这里Nu是控制时域，N是预测时域（通常Nu ≤ N）。预测采用的是"冻结"PPD的策略，即认为在当前时刻估计的φ̂(k)在预测时域内保持不变。这是MFAPC和显式MPC的一个重要区别——MFAPC不需要预测PPD的未来变化，直接用当前估计值就行，这大幅降低了计算复杂度。
优化准则函数选取：
J = Σ(i=1到N) λ_i·(y(k+i) - y*(k+i))² + ρ·Σ(j=1到Nu) (Δu(k+j-1))²
第一项是预测输出跟参考轨迹的误差惩罚，第二项是控制增量惩罚。λ_i是误差加权系数，ρ是控制量加权系数。ρ的作用是限制控制量的剧烈变化，ρ越大控制越"温柔"，但响应也越慢。
把这个二次型优化问题求解出来，得到控制增量序列的最优解，取第一个分量作用于系统，然后在下一时刻重新估计PPD，重新求解。这就是滚动优化（Receding Horizon）策略。具体的矩阵形式推导不在这里铺开了，实现的时候用二次规划或者直接解析求解都能做，因为目标函数是二次的，约束如果不加的话可以解析求解，速度快得多。
2.3 MFAILC控制律设计
MFAILC针对的是批次过程。设第i次运行的输入输出序列分别为u_i(k)，y_i(k)，k = 0, 1, ..., N。期望轨迹为y*(k)。
在每次运行中，同样利用紧格式动态线性化：
Δy_i(k+1) = φ_i(k)·Δu_i(k)
注意这里的φ_i(k)是第i次运行在第k时刻的PPD，它既随时间k变化，也随批次i变化。
MFAILC的控制律设计思路是：当前批次的输入=上一批次的输入+修正项。修正项由上一批次的跟踪误差驱动：
u_i(k) = u_{i-1}(k) + ρ·φ̂_i(k)·e_{i-1}(k+1)
这里e_{i-1}(k+1) = y*(k+1) - y_{i-1}(k+1)是上一批次在k+1时刻的跟踪误差。ρ是学习增益。φ̂_i(k)是当前批次的PPD估计值。
PPD的估计也会跨批次进行。一种常见做法是：
φ̂_i(k) = φ̂_{i-1}(k) + η·Δu_{i-1}(k)/(μ + Δu_{i-1}(k)²) · (Δy_{i-1}(k+1) - φ̂_{i-1}(k)·Δu_{i-1}(k))
这样做的好处是，PPD估计的历史信息可以跨批次传递，随着批次增加，PPD估计越来越准，控制性能也逐步提升。
MFAILC对初值的敏感度比MFAPC更明显。第一批次的输入u_0(k)怎么给、PPD初值φ̂_0(k)怎么设，直接影响收敛速度。通常的做法是把第一批次的输入设为一个简单的基准输入（比如期望轨迹的静态前馈值或者零输入），PPD初值设为一个小常数。
3. 仿真程序的结构与关键实现
3.1 被控对象选取与仿真配置
这套仿真程序选择了一个典型的非线性被控对象来做验证。工业过程控制里最常用来做数据驱动算法验证的例子之一就是：
y(k+1) = y(k) / (1 + y(k)²) + u(k)³
这个对象够"非线性"，但又不至于复杂到让人无法分析——y(k)/(1+y(k)²)项让系统增益随输出水平变化，u(k)³项让系统对控制输入的响应呈非线性放大。用这个对象跑MFAPC和MFAILC，能比较充分地展示算法的非线性适应能力。
仿真配置方面需要明确以下几组参数：
- 仿真总时长/总批次：连续运行仿真跑1000步，迭代学习仿真跑50个批次，每个批次200步
- 参考轨迹：连续运行用阶跃信号+正弦信号的组合；批次运行用一条固定的期望轨迹y*(k)
- 采样周期：设为单位1，离散时间系统天然满足
- 扰动设置：在输出端加入幅值0.01的白噪声，模拟测量噪声
程序中把系统模型、控制器、参数配置分成独立的模块，方便替换被控对象和算法。我用MATLAB写的这套程序，整体结构类似下面这样：
main_MFAPC.m % MFAPC主仿真脚本 main_MFAILC.m % MFAILC主仿真脚本 plant_nonlinear.m % 被控对象模型函数 controller_MFAPC.m % MFAPC控制器函数 controller_MFAILC.m % MFAILC控制器函数 ppd_estimator.m % PPD估计通用函数 plot_results.m % 结果绘图脚本分模块的好处是显而易见的：想换被控对象，只需要改plant函数；想对比不同参数的作用，只需要在配置区修改参数而不动核心算法代码。
3.2 MFAPC的仿真流程
MFAPC仿真的主循环每一步做四件事：估计PPD、预测输出、计算控制增量、施加控制并采集数据。
关键代码逻辑如下：
% 参数配置 alpha = 1; % PPD重置值 eta = 0.5; % PPD估计步长 mu = 1; % PPD估计惩罚因子 rho = 0.5; % 控制增量加权 N = 8; % 预测时域 Nu = 4; % 控制时域 lambda = ones(1, N); % 误差加权 % PPD估计与重置 phi_hat = alpha; delta_u = u(k-1) - u(k-2); delta_y = y(k) - y(k-1); if abs(delta_u) < 1e-6 phi_hat = alpha; % 输入变化太小，重置 else phi_hat = phi_hat + ... eta * delta_u / (mu + delta_u^2) * (delta_y - phi_hat * delta_u); end if abs(phi_hat) < 1e-4 phi_hat = alpha; % PPD退化，重置 end % 构建预测矩阵并解析求解控制增量 % A矩阵由phi_hat构成，维度Nu x N % 控制增量 = inv(A'*diag(lambda)*A + rho*I) * A'*diag(lambda)*(y* - y)这里有个关键点必须说明：预测矩阵的构造。由于预测用的是冻结PPD策略，未来N步的输出预测可以用当前y(k)加上φ̂(k)乘以控制增量的累积和来表示。如果Nu < N，那么从Nu+1步到N步的控制增量视为零，预测值只取决于前Nu个控制增量。这个细节决定了矩阵的维度和求解方式，写程序的时候很容易在这出错——把控制时域和预测时域搞混，导致矩阵维度不匹配。
我实测下来，预测时域N取8、控制时域Nu取4，对这个对象比较合适。N太小，预测优势体现不出来；N太大，冻结PPD的假设失真，性能反而变差。ρ取0.5左右能在快速性和平滑性之间取得较好的平衡。ρ太小会出现控制量高频抖动的情况，ρ太大则系统的上升时间明显变长。
3.3 MFAILC的仿真流程
MFAILC仿真采用批次循环加时域循环的双层结构。外层循环控制批次，内层循环控制每个批次内部的时域推进。
关键流程如下：
% 初始化 u = zeros(N_steps, 1); % 当前批次输入 u_prev = zeros(N_steps, 1); % 上一批次输入 phi_hat_matrix = ones(N_steps, 1) * 0.5; % PPD初值 for i = 1:N_batches % 运行当前批次 y = zeros(N_steps, 1); y(1) = 0; for k = 1:N_steps-1 y(k+1) = y(k) / (1 + y(k)^2) + u(k)^3; end % 计算跟踪误差 e = y_ref - y; % 更新PPD估计（跨批次） for k = 2:N_steps du = u(k-1) - u_prev(k-1); % 这里用的是相邻批次输入的差 dy = y(k) - y_prev(k); phi_hat_matrix(k) = phi_hat_matrix(k-1) + ... eta * du / (mu + du^2) * (dy - phi_hat_matrix(k-1) * du); end % 更新控制输入 for k = 1:N_steps-1 u(k) = u_prev(k) + rho * phi_hat_matrix(k+1) * e_prev(k+1); end % 保存当前批次数据，更新迭代变量 u_prev = u; y_prev = y; e_prev = e; end注意这里更新控制律用到的误差是上一批次的误差e_prev，而PPD估计用到了当前批次和上一批次的输入输出差。这种"批次间前馈"结构是ILC的本质：本次的控制输入不直接依赖当前批次的实时误差，而是利用历史批次的误差积累经验来修正。这也是ILC和实时反馈控制最大的区别——它更像是在"越做越好"，而不是"边做边防"。
实际操作中还有一个细节很容易踩坑：PPD更新里用到的du到底是时间方向上相邻控制量的差，还是批次方向上相邻控制输入的差？我最初写程序的时候混用了这两种差分，导致PPD估计完全发散。正确做法是：在MFAILC里，紧格式动态线性化是沿着批次方向定义的，即Δu_i(k) = u_i(k) - u_{i-1}(k)，Δy_i(k+1) = y_i(k+1) - y_{i-1}(k+1)。用错了方向，算法的收敛性直接就没了。
3.4 结果分析与可视化
仿真跑完以后，需要从几个维度评估算法性能：
- 跟踪误差：MFAPC看稳态阶段的均方根误差（RMSE），MFAILC看每批次的最大绝对误差（MAE）随批次的变化
- 收敛速度：MFAILC要看误差收敛到稳定水平需要的批次数量，MFAPC看上升时间和超调量
- 控制量行为：观察控制量是否平滑，是否有频繁饱和或振荡的迹象
画图方面，我习惯用三张图来呈现MFAPC的结果：输出跟踪曲线（纵轴是y和y*）、控制输入曲线（纵轴是u）、PPD估计值曲线（纵轴是φ̂）。PPD估计曲线特别值得看，它能直观反映算法对系统增益变化的"感知"——如果PPD在系统输出变化剧烈的位置有大幅波动，说明估计在正常工作；如果PPD一直不变或者剧烈震荡，那参数一定有问题。
MFAILC的图则更关注批次维度的演化：可以画热力图展示不同批次的输出轨迹，也可以画误差随批次变化的收敛曲线。我一般画两个图：一个是最终批次和第一批次的输出轨迹对比，另一个是各批次最大绝对误差的对数坐标曲线。后一张图能清晰看出收敛趋势——正常情况应该是误差随批次呈近似指数衰减，如果误差曲线出现平台或者发散，说明学习增益ρ或者PPD参数需要调整。
4. 参数整定与调优实战
4.1 PPD估计参数的影响规律
PPD估计器里的η和μ是一对需要配合调整的参数。
η是步长因子，决定PPD估计的更新速度。η偏大，PPD跟踪能力强但容易受噪声影响，估计值可能出现高频抖动；η偏小，PPD估计平滑但可能跟不上实际的增益变化。实测下来，η在0.3到0.7之间是大多数对象的合理区间。注意η跟控制量增量Δu的幅值有耦合——如果控制量本身数值很大，η就要适当调小，保证η·Δu/(μ+Δu²)这个增益不至于过大。
μ是惩罚因子，直观效果是"阻尼"——μ越大，PPD变化越慢。μ取值如果太小，当Δu接近零的时候，η·Δu/(μ+Δu²)会变得非常大，PPD估计会出现尖峰。这就是为什么我前面强调μ不能太小，一般取1左右，具体要看控制量的量级。如果控制量是千量级，μ就得按1000的量级来取。
我自己的调参顺序是先定μ，保证PPD估计不发散；然后调η，让PPD能跟上系统增益变化；最后才动控制端的ρ和预测时域。很多新手一上来就四个参数一起调，出了问题根本分不清是哪个参数引起的，这是调参的大忌。
4.2 MFAPC的预测时域与控制时域匹配
N和Nu的选择有很强的经验性。
预测时域N决定了算法"看得多远"。N越大，控制决策越有前瞻性，对延迟系统的好处越明显；但N过大时，冻结PPD的假设在长时域内可能严重失真，导致预测输出跟实际输出偏离过大，控制效果反而恶化。对这个测试对象来说，N=8是个甜点值，N超过15以后性能明显下降。
控制时域Nu的选择跟被控对象的相对阶和系统惯性有关。系统惯性越大，Nu越应该大一些，让控制器有足够的自由度来安排控制增量的变化。但Nu太大会让计算量增加，而且对惩罚项ρ的调节压力增大。通常取值Nu = N/2左右是个不错的起点，比如N=8时Nu=4。
这里有一个具体验算过的例子。我的程序中，预测矩阵A的构造方式是：第j行（j=1到N）第l列（l=1到Nu）的元素，如果j >= l则是φ̂(k)，否则为0。原因很简单：y(k+1)受Δu(k)影响，y(k+2)受Δu(k)和Δu(k+1)共同影响，以此类推。如果没有这个累积效应，预测模型就是错的，控制效果一定崩溃。
4.3 MFAILC的学习增益与收敛性
MFAILC里最关键的就是学习增益ρ。ρ的理论取值范围需要满足收敛条件，在紧格式动态线性化框架下，通常要求|1 - ρ·φ̂(k)| < 1，也就是0 < ρ·φ̂(k) < 2。
但实际调试中发现，这个界限只是一个必要条件，不是充分条件。当ρ·φ̂(k)接近2的时候，批次间的误差会出现震荡收敛——前期下降很快，但后期出现波浪状的波动，很难收敛到很高的精度。把ρ·φ̂(k)控制在0.3到0.8之间，收敛最平滑。
还有个心得是：MFAILC的PPD初始值对前几个批次的影响很大。如果φ̂_0(k)跟实际增益偏差过大，前几批次的误差可能反而不降反升，看起来"发散了"。但不要急着判死刑，多做几个批次再看。我遇到过的情况是前3批误差上升，第5批才开始下降，第15批左右收敛到满意的水平。这说明ILC本身需要一定的"学习期"，评估收敛性要拉长批次看趋势。
5. 常见问题与排查技巧实录
5.1 程序跑飞与数值发散
最常见的现象就是y(k)迅速涨到NaN或者Inf。排查顺序我建议是固定的：
先查PPD的复位机制。看Δu(k-1)是否出现过零值——如果控制量在一段时间内保持不变，Δu就是零，那么PPD估计公式的分母就是μ，虽然不至于除零，但PPD估计会漂移。必须在估计之前检查|Δu|是否小于阈值，小于就强制复位。
再查PPD是否出现负值或者异常大值。对于这个对象，PPD应该是正值（系统在输入增加时输出增加），如果PPD估计出负值，控制方向就反了，系统必然发散。所以复位条件里一定要加|φ̂(k)| < ε的判断，用alpha重新赋值。
最后查控制增量是否过大。MFAPC求解出来的Δu如果乘上φ̂之后预测输出变化远超实际允许范围，说明ρ太小或者λ矩阵没设好。可以加一个控制增量限幅：|Δu(k)| ≤ Δu_max。虽然理论上MFAC框架不强制要求限幅，但工程上加上限幅能显著提高鲁棒性，代价是可能牺牲一点理论上的完美性。
5.2 MFAILC误差不收敛怎么办
MFAILC跑完50个批次，误差还纹丝不动，甚至越来越大。这是很多人初次跑MFAILC时最崩溃的时刻。
第一个排查点还是PPD的差分方向。百分之八十的"不收敛"都是因为程序里把批次方向和时间方向的差分搞混了。检查一下你更新PPD的时候用的Δu到底是u_i(k) - u_{i-1}(k)（正确）还是u_i(k) - u_i(k-1)（错误）。后者把迭代学习变成了实时差分，概念上就错了。
第二个排查点是第一批次的输入。如果第一批次的输入跟期望轨迹所需的输入量级差太多，学习过程需要很长时间才能追上。建议先跑一个开环仿真，看看给定一个合理输入时对象输出大概是什么量级，然后把第一批量输入设置在这个合理值附近，别让系统从完全错误的起点开始学习。
第三个排查点是ρ是否过于保守。我见过有人把ρ设成0.01，跑50批次的误差曲线跟水平线似的。学习增益太小，每批次只能修正一点点误差，需要极多批次才能看到效果。把ρ提到0.3左右，通常几批次就能看到明显下降。
5.3 仿真速度优化技巧
MATLAB跑这个仿真，如果批次多、时域长，循环的写法会极大影响速度。一个非常实用的优化是：把内层时域循环尽量向量化。
我最初写的MFAILC程序用了三层嵌套循环——批次循环、时域循环、PPD更新循环，跑50个批次每个批次200步居然花了十几秒。改成部分向量化之后，同样的仿真只需要几秒。关键是把PPD更新和输入更新用矩阵运算替代逐点循环。如果不太擅长向量化，也可以把整个仿真放到parfor并行循环里——不同批次虽然存在递推关系，但从第二批到第N批的PPD初始化不依赖上一批的完整结果，可以并行计算各自批次内的响应，最后再汇总误差收敛曲线。
不过说实话，对于单次仿真跑几秒钟这种规模，优化不优化都没太大关系。真正要优化的是批量调参的场景——比如你想扫描ρ的20个候选值，每个值跑50批次，那就要考虑并行或者把扫描任务拆到多个工作日脚本里跑。这里可以用一个朴素的技巧：先跑粗扫描，确定大致区间，再在区间内细扫。别一上来就全体网格搜索，时间全浪费在完全没希望的参数组合上了。
6. 扩展思路与后续方向
这套仿真程序的价值不只是在单机上跑出两张收敛曲线。沿着现在这个框架，还能扩展出不少有价值的方向。
首先是抗扰动和鲁棒性验证。目前仿真是在理想化的设定下跑的，实际系统里还有未建模动态、时变参数、输入饱和等约束。可以在程序里把被控对象换成更苛刻的版本——比如加上时变增益系数、加入输入端限幅、在模型参数上叠加随机漂移，看看MFAPC和MFAILC的鲁棒性到底怎么样。我自己跑下来发现，MFAPC对参数时变的适应能力明显强于固定参数的MPC，这跟PPD在线更新的机制直接相关。
其次是推广到多输入多输出系统。现在的实现是单输入单输出，扩展的思路是引入分块PPD矩阵，把标量PPD变成矩阵PPD，控制律推导也从标量代数变成矩阵运算。MIMO版本的MFAPC用到了矩阵求逆，计算复杂度会明显上升，但对这类系统的鲁棒性提升很有研究价值。
再一个方向是跟我做过的其他控制算法做对比研究。把MFAPC和MFAILC的结果与标准MPC、传统ILC放在一起比较，量化在建模成本、在线计算量、控制性能三个维度上的差异。这不仅是论文需要的对比实验，也能帮自己判断：在什么场景下值得用无模型方法，什么场景下建个简单模型用传统方法反而更省事。
最后提一个工程化的建议：这套程序的模块化结构很适合改造成一个通用的数据驱动控制算法库。把PPD估计器封装成独立的类或函数，把控制器策略做成可切换的模式，就能快速测试新的改进算法——比如多步预测时域变化策略、变学习增益MFAILC、带遗忘因子的PPD估计等。封装好之后，再加新算法只是加一个函数的事，而不是从头重写整套仿真。
我在实际使用中的体会是，这类仿真程序的价值不在于代码本身有多精巧，而在于它给了你一个反复试验的沙盒。参数怎么调、算法怎么改、理论界和工程实践的差距在哪里，跑上几十组仿真就全清楚了。数据驱动控制看起来门槛高，但真正动手把MFAPC和MFAILC跑通了之后，你对整个MFAC框架的理解会完全不一样。
