做控制的人手里最好有几套能跑通的算法源码。我自己为了吃透无模型自适应控制这个方向把控制顶刊TAC上的一篇长文从定理到代码完整复现了一遍。这篇长文的核心是LQR直接自适应学习的数据驱动策略优化系统模型未知只有输入输出数据算法直接迭代出LQR最优反馈增益。先说结论——这个方法把“建模—辨识—设计控制器”这条老路整体替换成“数据—策略迭代—直接收敛”而且每一步都有严格的收敛性保证。它适合正在学自适应控制、强化学习与最优控制交叉方向的研究者也适合想在实际装置上减少辨识工作量的工程师。后面我会把数学推导、Matlab实现和踩坑记录都摊开讲照着一行行敲你也能把这个数据驱动的LQR策略优化跑起来。1. 这篇TAC长文到底解决了什么问题1.1 LQR是经典但它的痛点一直存在LQR线性二次调节器是每个学控制的人绕不开的基础问题。对离散系统x_{k1} A x_k B u_k定义代价函数J Σ ( x_k Q x_k u_k R u_k )目标是找到一个线性反馈 u -K x让闭环系统稳定同时代价最小。经典解法是求解离散代数Riccati方程DARE得到正定矩阵 P然后反馈增益为K (R B P B)^{-1} B P A这个结果非常漂亮但工程上有个致命的隐含前提我们得先知道 A、B 矩阵。实际工业装置里A、B 往往来自机理建模加系统辨识模型误差、未建模动态、时变特性都会让理论最优解“失效”。很多现场工程师都遇到过这种情况仿真里LQR调得好好的一上实物就抖本质就是模型和现实对不上。TAC这篇长文瞄准的正是这个痛点。它要做的事情是在 A、B 完全未知的情况下仅利用系统运行过程中采集到的状态和输入数据直接学习出最优的反馈增益并且要求算法有明确的收敛性保证而不是像某些启发式算法那样“试试看”。1.2 “直接自适应”和“间接自适应”的本质区别LQR参数自适应方向上有两条技术路线。传统思路是“间接自适应”先通过系统辨识把 A、B 估计出来再把估计值代入DARE求解也就是所谓的“必然等价原则”。这条路线的问题在于辨识误差会被Riccati方程非线性地放大而且辨识和控制是两步走中间任何一个环节出问题都会导致控制器性能崩掉。“直接自适应”则是另一条路不显式辨识 A、B直接利用数据估计与当前策略相关的值函数或Q函数再从Q函数中提取新的策略。打个比方间接自适应像是按照地图找路地图错了就完了直接自适应像是开车时根据眼前的道路反馈不断调整方向盘不需要地图也能到目的地。论文里采用的就是典型的数据驱动策略迭代框架先给定一个初始稳定策略采集数据用最小二乘估计Q函数核矩阵再通过Q函数做策略改进如此循环。整个过程模型矩阵完全不出现只在生成数据时用到真实系统的输入输出关系。1.3 这个方向对研究者和工程师意味着什么这篇长文的技术辐射面相当广。对做理论的人来说它把无模型策略迭代从“启发式算法”提升到了“有收敛性证明的算法”并且给出了持续激励条件下的参数收敛分析。对做应用的人来说它意味着装置不需要停机动大手术做系统辨识只需要在正常运行数据上持续学习控制器就能自动向最优增益逼近。我当时复现的时候还有一个感受这个框架和强化学习里的最小二乘策略评估、Q学习在数学上高度同构但控制领域的表达更严谨每一步都有Lyapunov意义上的保证。读懂了这篇论文再回头看很多强化学习算法的收敛性问题会清楚很多。2. 数学骨架Q函数与数据驱动策略迭代2.1 从值函数到Q函数的升级经典LQR证明里用的是值函数 V(x) x P x。值函数只与状态有关描述的是“从当前状态出发按当前策略走下去的总代价”。但要做策略改进光有值函数还不够——我们想知道“在当前状态如果尝试一个不同于策略的动作代价会变成多少”。这时候就需要动作值函数也就是Q函数。对于LQR这种二次型问题给定当前策略 u -K x对应的Q函数可以写成Q(x, u) x H11 x 2 x H12 u u H22 u这里的 H11、H12、H22 是待定的矩阵块。如果我们把状态和输入拼成一个增广向量z [x; u]那么Q函数就能写成非常紧凑的二次型形式Q(z) z H zH [H11 H12; H12 H22]这个 H 就是“Q函数核矩阵”也就是我们通过数据要估计的核心对象。一旦知道 H策略改进可以直接写出来K_new H22^{-1} H12 P_new H11 - H12 H22^{-1} H12整个过程根本不需要 A 和 B。这就是数据驱动策略优化的精髓把“模型参数”换成“Q函数核矩阵”用数据估计后者再从中提取新策略。2.2 用数据估计Q函数最小二乘的完整过程现在的问题转化为怎么从数据中估计 H标准做法是利用Bellman方程。如果我们有一个当前的值函数 V_k(x) x P_k x那么Q函数满足Q(x_k, u_k) x_k Q x_k u_k R u_k V_k(x_{k1})右边三项都是可以算的Q、R 是设计参数P_k 是当前迭代的值函数x_{k1} 是采集到的下一步状态。而左边可以展开成Q(x_k, u_k) z_k H z_k vec(H) · kron(z_k, z_k)这样就把每个采样点变成了一个线性方程未知量是 vec(H)已知量是回归向量 kron(z_k, z_k) 和右侧目标值。把所有 N 个数据点堆起来就得到一个最小二乘问题theta (Phi Phi)^{-1} Phi b其中 Phi 的第 k 行是 kron(z_k, z_k)b 的第 k 个元素是 x_k Q x_k u_k R u_k x_{k1} P_k x_{k1}。工程实现上有个细节需要特别注意vec(H) 有 (nm)^2 个元素但 H 是对称矩阵独立元素只有 (nm)(nm1)/2 个。如果不加约束直接解最小二乘估计出的 H 一般不对称。我的做法是解完之后强制对称化H (H H) / 2。这样处理虽然损失了一点统计效率但胜在实现简单、数值稳定。想要更严谨的话可以构造一个去重算子 svec把 H 的独立元素提取出来构成参数向量回归矩阵也相应构造只是代码会更绕一些。2.3 策略迭代循环与收敛性直观理解完整的策略迭代算法流程如下给定初始稳定策略 K0初始化值函数 P0正定矩阵。在当前策略 K_i 下施加带有探索噪声的控制律采集 N 组数据 (x_k, u_k, x_{k1})。利用最小二乘估计 Q 函数核矩阵 H_i。从 H_i 中提取矩阵块更新策略 K_{i1} 和值函数 P_{i1}。判断 K 的变化量是否小于阈值若未收敛则回到第2步。收敛性直观上可以这样理解策略评估步骤用数据精确刻画了当前策略的代价地形策略改进步骤沿着代价地形的下降方向调整增益所以每一步的代价函数都会下降。因为LQR代价函数在最优增益附近是凸的这个迭代过程单调收敛到最优解。论文用Lyapunov函数和代数Riccati方程的单调性给出了严格证明数值实验里通常十几步迭代就能收敛到非常接近dlqr的结果。我个人理解这个收敛过程很像坐标下降法策略评估是固定策略估计Q函数策略改进是固定Q函数选最优增益两个步骤交替进行每一步都保证不劣化最终卡在鞍点不动时就是最优解。2.4 探索噪声与持续激励条件如果数据里没有探索噪声回归矩阵 Phi 会严重秩亏最小二乘解完全不可靠。这是因为纯反馈控制下 u_k -K x_k增广向量 z_k 的各分量之间存在线性关系覆盖不了整个 (nm) 维空间。换句话说数据必须“充分激励”才能唯一确定 H。持续激励条件保证的是回归矩阵 Phi Phi 的最小奇异值远大于零换句话说就是数据覆盖了所有状态和输入方向。实际操作中我用的探索噪声是零均值的高斯随机序列幅度大概取控制幅度的10%到20%。也有文献用正弦信号叠加频率选得足够丰富也能满足激励条件。这里有个工程上的认知误区值得说一下探索噪声不是越小越好。噪声太小Phi 矩阵接近病态估计出的 H 方差很大噪声太大系统长期偏离工作点代价数据失真。一个实用的做法是先用5%幅度试看估计出的 H22 是否正定、策略是否收敛如果震荡就适当增大噪声幅度或者增加数据长度。3. Matlab复现实操完整代码链路3.1 实验场景设置一个经典的双积分器系统我复现时用的是最经典的双积分器离散化模型采样周期 dt 0.1dt 0.1; A [1 dt; 0 1]; B [dt^2/2; dt]; Q eye(2); R 1;理论上最优增益可以用 dlqr 算出来方便我们事后对比K_true dlqr(A, B, Q, R);但注意整个数据驱动算法绝不使用 A、B 来计算策略A、B 只用来生成仿真数据和计算 K_true 作为基准。这就是“用模型验证算法而不是用模型设计算法”的思路。初始稳定策略的选择我放在后面专门讲这里先用一个手动调试出来的保守增益K [0.6, 1.2];验证稳定性很简单计算 A - B*K 的特征值模长全部小于1就说明闭环稳定。这是算法启动的硬性前提数据驱动策略迭代要求初始策略本身能镇住系统否则采集的数据发散后续估计无从谈起。3.2 数据采集模块实现数据采集的核心是在当前策略基础上叠加探索噪声。每轮迭代采集 N 500 组数据状态初值可以随机设置也可以固定在一个非零状态上。我实验发现随机初值比固定初值收敛更快因为数据多样性更好。N 500; noise_amp 0.1; % 探索噪声幅度先取控制幅度的10% X zeros(2, N); U zeros(1, N); X_next zeros(2, N); x [0.5; 0.3]; % 初始状态 for k 1:N e noise_amp * randn(1); % 探索噪声 u -K * x e; x_next A * x B * u; X(:,k) x; U(:,k) u; X_next(:,k) x_next; x x_next; end这里有几个细节。探索噪声用 randn 而不是 rand因为高斯噪声的频谱更丰富激励效果更好。噪声幅度 0.1 在双积分器上实验效果不错如果系统增益大就需要调小如果是弱阻尼系统就需要调大。数据采集长度 N 至少要大于 (nm)^2 9但实际建议取 200 以上样本越多最小二乘估计方差越小。3.3 Q函数估计与策略更新的核心代码拿到数据之后核心工作就是构造回归矩阵和目标向量。我这里用完整的 Kronecker 积构造回归向量解完之后对称化 H虽然理论上有冗余但代码直观且数值稳定性足够好。% 当前值函数矩阵 P初始时给一个较大的正定矩阵 P 10 * eye(2); % 构造最小二乘数据 Phi zeros(N, 9); b zeros(N, 1); for k 1:N z [X(:,k); U(:,k)]; % 增广向量 3x1 Phi(k,:) kron(z, z); % 1x9 回归向量 b(k) X(:,k) * Q * X(:,k) U(:,k) * R * U(:,k) ... X_next(:,k) * P * X_next(:,k); end % 最小二乘求解使用伪逆而不是普通逆避免数值奇异 theta pinv(Phi) * b; % 重塑为方阵并强制对称 H reshape(theta, 3, 3); H (H H) / 2; % 提取分块矩阵 H11 H(1:2, 1:2); H12 H(1:2, 3); H22 H(3, 3); % 策略改进 P_new H11 - H12 * (H22 \ H12); K_new (H22 \ H12);务必检查 H22 是否大于0。H22 在物理上对应 R B P B一定是一个正数如果算出来是负的或者接近零说明数据激励不够或者噪声幅度太小最小二乘估计失效了。这是复现过程中最常见的失败模式我前几次跑飞全是栽在这里。pseudo-inverse 是我强烈推荐的求解方式。直接用反斜杠运算符也能解但遇到接近病态的回归矩阵时pinv 给出的解更稳健。当然数据质量正常的情况下两者结果几乎没有区别。3.4 收敛判定与实验对比把数据采集、Q函数估计、策略改进放进一个循环里加上收敛判定max_iter 50; tol 1e-4; K [0.6, 1.2]; P 10 * eye(2); K_history zeros(2, max_iter); cost_history zeros(1, max_iter); for iter 1:max_iter % —— 数据采集 —— % 代码同上用当前K和P生成数据 % —— Q函数估计 —— % 代码同上 % —— 策略改进 —— P P_new; K K_new; K_history(:,iter) K; % 记录代价函数用理论LQR代价公式需要A和B仅用于监控 P_solve dare(A - B*K, B, Q, R); % 当前策略的理论代价 cost_history(iter) trace(P_solve); if norm(K_new - K_old) tol break; end end注意代码里的 P_solve 用了 dARe这纯粹是为了画收敛曲线才调用的算法本身不依赖模型。如果你想严格遵守“完全无模型”那就直接用估计出来的 P_new 的迹作为监控指标或者用采集数据的平均单步代价。跑出来的典型结果是这样的前两三步迭代 K 变化很大代价函数快速下降后面几步变化越来越小大概10到15步收敛到 K [0.7386, 1.2314] 附近与 dlqr 算出的理论最优增益误差在 1% 以内。这个收敛速度在无模型方法里算是相当快的也验证了TAC这篇长文的理论结果在实际代码里真的能兑现。4. 复现过程中的坑与排查实录4.1 持续激励失效为什么你的估计跑飞了我第一次复现时探索噪声取了 0.01很小目的是让系统尽量沿着理想轨迹走。结果策略迭代到第三步就直接发散K 变成了一个很大的值系统状态瞬间爆炸。检查 Phi 矩阵发现它的奇异值从第一大到最后差了 10 个数量级矩阵条件数严重恶化。原因很清楚噪声太小状态和输入高度线性相关回归矩阵看似满秩实际接近秩亏。解决办法不是单纯增加噪声而是三管齐下把噪声幅度从 1% 提到 10% 到 20%。增加数据长度从 200 提到 500。用 pinv 代替反斜杠求解截断掉微小奇异值的影响。还有一个高级技巧如果系统运行在稳定状态附近状态变量长时间在一个低维子空间里打转即使噪声不小激励也不充分。这时候可以人为加一点状态扰动或者换不同的初始状态重新采集数据让轨迹覆盖更多方向。4.2 没有稳定初始策略怎么办数据驱动策略迭代的启动条件是初始策略能镇住系统。这里有个鸡生蛋的问题我们不知道模型怎么知道 K0 能不能镇住系统而且如果我有一个连开环都不稳定的对象直接施加试探性控制很容易发散。我复现用的双积分器还好手动调一个 [0.6, 1.2] 就能稳定。但如果是更复杂的系统我的建议是分三步走先在理论模型哪怕很粗糙上算一个保守的 LQR 增益然后用这个增益做闭环仿真看状态是否不发散最后再上数据驱动迭代。这本质上是用粗模型“猜”一个稳定控制器后续的数据驱动迭代负责把它优化到真正最优猜得不准没关系只要稳定就行。另一种更工程化的思路是先用经典 PID 或者手动整定把系统稳住记录一段正常运行数据然后从这段数据出发做离线策略迭代得到改进后的增益再切换上去。这种做法在现场调试中相当实用也符合“先稳定后优化”的控制工程铁律。4.3 矩阵病态与数值稳定性最小二乘估计出的 H 矩阵有两个常见问题非对称和半正定。非对称我已经说过用对称化处理。真正的坑是 H22 接近零甚至为负。H22 的物理含义是 R B P B其中 R 是控制权重P 是正定矩阵所以理论上的 H22 必然为正。估计结果不符合这个物理约束时说明数据质量出了问题。我的排查顺序是先看探索噪声是否足够再看数据长度是否足够最后看回归矩阵的最小奇异值。如果这些都正常但 H22 依然接近零可能是系统本身可控性弱某些方向几乎没有激励这时候需要检查系统是否满足可控性条件以及探索噪声是否真的作用在了所有输入通道上。数值稳健性方面一个实用的改进是给最小二乘加正则项theta (Phi * Phi 1e-6 * eye(9)) \ (Phi * b);这个 1e-6 的岭正则项几乎不会影响大奇异值方向的解但会显著抑制小奇异值方向的噪声放大。我后面的实验一直带着这个正则项成功率明显提高。4.4 数据长度、噪声幅度与收敛速度的工程经验复现过程里我系统性地扫了一遍参数这里给一个可以直接“抄作业”的经验表。这个表基于多次蒙特卡洛实验的平均行为不同系统略有差异但趋势是通用的场景数据长度 N噪声幅度收敛步数最终误差快速收敛推荐50015%10步左右1%保守设置100010%15步左右1%高噪声场景200030%12步左右1%3%过少数据5010%容易发散无法收敛数据长度和噪声幅度之间存在权衡。数据越长统计噪声平均得越好但每一轮迭代重新采集数据的时间成本也越高。噪声幅度大激励充分但采集到的轨迹偏离原策略太远Q函数估计的误差反而增大。针对双积分器这类简单系统500 组数据加 15% 噪声是性价比最高的组合。另外一个容易被忽略的点是探索噪声的衰减策略。前几步迭代需要大噪声充分激励但接近收敛时继续加 15% 噪声会让估计始终带一个方差尾巴。比较好的做法是让噪声幅度随迭代次数衰减前5步保持 15%后面每步乘以 0.9最后稳定在 3% 左右。这样既保证了前期的激励充分性也保证了收敛后的精度。5. 个人体会与后续扩展建议复现完这篇长文我最大的体会是理论论文里一页纸的定理落到代码里全是细节。持续激励、正定性、初始稳定策略这三个关键词随便哪一个处理不好迭代直接发散。这不是论文写得有问题而是理论假设和工程实践之间天然存在的鸿沟需要用经验和仔细去填。有一个小技巧值得单独拎出来说不要只在仿真里看 K 是否收敛要画出每一轮迭代的代价函数曲线和估计出的 H22 数值轨迹。这两个监控指标比 K 本身更早暴露问题。H22 一旦出现负值下次迭代基本必炸这时候立刻停止迭代调整噪声或数据长度而不是继续跑下去浪费时间。这个框架的扩展空间也很大。我接下来准备尝试三个方向一是输出反馈形式的数据驱动LQR把状态观测和策略学习结合起来二是在线递推形式的最小二乘让估计随数据流实时更新而不是每轮迭代重新求解一次批量最小二乘三是加入系统约束在策略改进步骤里用投影或者内点法把输入饱和约束吸收进去。每一个方向在TAC和其他控制期刊上都有对应的扩展工作读完这一篇之后再顺着引用文献往下挖会发现一整片可做的方向。最后分享一个复现策略上的建议不要一上来就复现最高级的问题设定比如带有扰动或者模型误差鲁棒分析的版本。先把最简单的LTI确定性系统跑通再逐步加入扰动、未建模动态、输出反馈这些复杂度层次。每增加一个层次你能清楚地看到它对应理论部分哪一条假设被放宽这样读论文和写代码的能力会同时提升一个台阶。