我做控制仿真这几年有一个特别深的体会很多算法论文写得天花乱坠但真到自己动手做数值验证的时候光是把伪偏导数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(k1) f(y(k), y(k-1), ..., y(k-ny), u(k), u(k-1), ..., u(k-nu))在满足一定条件偏导数连续、广义Lipschitz等的情况下可以写成Δy(k1) φ(k) · Δu(k)其中φ(k)就是伪偏导数是一个时变标量。Δy(k1) y(k1) - 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(k1) y(k) φ̂(k)·Δu(k) y(k2) y(k1) φ̂(k)·Δu(k1) y(k) φ̂(k)·(Δu(k) Δu(k1)) ... y(kNu) y(k) φ̂(k)·Σ(i1到Nu) Δu(ki-1) ...这里Nu是控制时域N是预测时域通常Nu ≤ N。预测采用的是冻结PPD的策略即认为在当前时刻估计的φ̂(k)在预测时域内保持不变。这是MFAPC和显式MPC的一个重要区别——MFAPC不需要预测PPD的未来变化直接用当前估计值就行这大幅降低了计算复杂度。优化准则函数选取J Σ(i1到N) λ_i·(y(ki) - y*(ki))² ρ·Σ(j1到Nu) (Δu(kj-1))²第一项是预测输出跟参考轨迹的误差惩罚第二项是控制增量惩罚。λ_i是误差加权系数ρ是控制量加权系数。ρ的作用是限制控制量的剧烈变化ρ越大控制越温柔但响应也越慢。把这个二次型优化问题求解出来得到控制增量序列的最优解取第一个分量作用于系统然后在下一时刻重新估计PPD重新求解。这就是滚动优化Receding Horizon策略。具体的矩阵形式推导不在这里铺开了实现的时候用二次规划或者直接解析求解都能做因为目标函数是二次的约束如果不加的话可以解析求解速度快得多。2.3 MFAILC控制律设计MFAILC针对的是批次过程。设第i次运行的输入输出序列分别为u_i(k)y_i(k)k 0, 1, ..., N。期望轨迹为y*(k)。在每次运行中同样利用紧格式动态线性化Δy_i(k1) φ_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}(k1)这里e_{i-1}(k1) y*(k1) - y_{i-1}(k1)是上一批次在k1时刻的跟踪误差。ρ是学习增益。φ̂_i(k)是当前批次的PPD估计值。PPD的估计也会跨批次进行。一种常见做法是φ̂_i(k) φ̂_{i-1}(k) η·Δu_{i-1}(k)/(μ Δu_{i-1}(k)²) · (Δy_{i-1}(k1) - φ̂_{i-1}(k)·Δu_{i-1}(k))这样做的好处是PPD估计的历史信息可以跨批次传递随着批次增加PPD估计越来越准控制性能也逐步提升。MFAILC对初值的敏感度比MFAPC更明显。第一批次的输入u_0(k)怎么给、PPD初值φ̂_0(k)怎么设直接影响收敛速度。通常的做法是把第一批次的输入设为一个简单的基准输入比如期望轨迹的静态前馈值或者零输入PPD初值设为一个小常数。3. 仿真程序的结构与关键实现3.1 被控对象选取与仿真配置这套仿真程序选择了一个典型的非线性被控对象来做验证。工业过程控制里最常用来做数据驱动算法验证的例子之一就是y(k1) y(k) / (1 y(k)²) u(k)³这个对象够非线性但又不至于复杂到让人无法分析——y(k)/(1y(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那么从Nu1步到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(k1) 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(k1) * e_prev(k1); 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(k1) y_i(k1) - y_{i-1}(k1)。用错了方向算法的收敛性直接就没了。3.4 结果分析与可视化仿真跑完以后需要从几个维度评估算法性能跟踪误差MFAPC看稳态阶段的均方根误差RMSEMFAILC看每批次的最大绝对误差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的假设在长时域内可能严重失真导致预测输出跟实际输出偏离过大控制效果反而恶化。对这个测试对象来说N8是个甜点值N超过15以后性能明显下降。控制时域Nu的选择跟被控对象的相对阶和系统惯性有关。系统惯性越大Nu越应该大一些让控制器有足够的自由度来安排控制增量的变化。但Nu太大会让计算量增加而且对惩罚项ρ的调节压力增大。通常取值Nu N/2左右是个不错的起点比如N8时Nu4。这里有一个具体验算过的例子。我的程序中预测矩阵A的构造方式是第j行j1到N第l列l1到Nu的元素如果j l则是φ̂(k)否则为0。原因很简单y(k1)受Δu(k)影响y(k2)受Δu(k)和Δu(k1)共同影响以此类推。如果没有这个累积效应预测模型就是错的控制效果一定崩溃。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框架的理解会完全不一样。