最近把TAC控制领域公认的顶级期刊上一篇关于LQR直接自适应学习的数据驱动策略优化长文完整复现了一遍Matlab版本跑通之后第一反应是原来标题看上去那么硬核的东西核心逻辑拆开之后竟然这么干净。论文要解决的问题一句话就能讲清楚系统模型未知时不通过模型辨识直接利用在线数据学习线性二次调节器LQR的最优反馈控制律。这个问题的实际价值很大——实验室里我们习惯用精确模型设计控制器但工程现场很多被控对象要么拿不到模型要么模型参数一直在漂这时候数据驱动策略优化就是一条非常现实的出路。如果你在做最优控制、自适应控制相关方向或者正在被“模型失配”折磨这篇复现笔记应该能帮你节省大量试错时间。我会把问题动机、数学原理、Matlab实现、参数调节和踩坑过程一五一十列出来。1. 这篇TAC长文到底解决了什么问题1.1 传统LQR设计的“模型依赖”病根在哪里LQR本身是控制理论里最老牌的问题线性系统二次型代价最优控制器是状态反馈u-KxK直接由代价权重和系统矩阵A、B通过代数Riccati方程解出来。看起来非常简单但整个设计链条的第一步就卡在“要拿到A和B”。实际系统里状态空间模型通常是机理建模加实验辨识拼出来的辨识结果和真实对象之间总存在误差更麻烦的是很多对象的动态特性会随工况变化比如温度、负载、老化都会让A和B慢慢偏离名义值。用失配模型算出来的K轻则性能下降重则闭环不稳定。所以“模型未知”不是理论家的洁癖而是工程现场每天都要面对的现实。1.2 直接自适应学习为什么比“先辨识再设计”更直接面对模型缺失传统做法是“间接自适应控制”先在线辨识模型再把辨识结果代入Riccati方程。这个思路没问题但存在一个结构性的弱点——辨识误差会在后续设计步骤中被放大而且为了满足辨识条件往往需要设计专门的激励信号这会让控制输入偏离正常操作。数据驱动策略优化走的是另外一条路不显式估计A、B而是直接把状态的测量值、实际施加的输入和下一步状态这三个量拿到手里从贝尔曼方程出发用最小二乘去拟合当前策略下的动作值函数再从动作值函数里提取更好的反馈增益。这很像人学骑车你不会先建模再解方程而是在一次次“试错加修正”中直接学会怎么把方向、油门配合好。1.3 我为什么选择这篇长文作为复现对象我挑这篇TAC长文来复现原因有三。第一研究对象是最基本的LQR没有非线性、没有复杂约束数学结构干净方便用几十行代码对应到论文里的每个公式非常适合作为数据驱动控制方向的入门复现项目。第二这类长文通常会把理论证明、算法流程、仿真实验打包得很完整算法伪代码往往是自洽的不像某些会议短文那样关键步骤要靠猜。第三LQR的数据驱动策略优化是后续很多工作的基石——模型参考、鲁棒控制、输出反馈、以及带约束的MPC式学习很多都能从这个起点延伸出去。把这一篇跑通后面读同类文章会轻松得多。我复现时给自己定的验收标准是学习出来的反馈增益与用真实模型解Riccati方程得到的K*基本一致并且迭代残差单调下降。2. 核心数学原理从Riccati方程到数据驱动策略优化2.1 LQR最优反馈与代数Riccati方程的关系离散系统x_{k1}Ax_kBu_k代价是Σ(xQxuRu)最优反馈为u-Kx其中K(RBPB)^(-1)BPAP是离散代数Riccati方程PQAPA-APB(RBPB)^(-1)BPA的唯一正定解。这套理论非常成熟模型已知时一条dare命令就能出结果。但我复现时始终提醒自己K的表达式只是“答案的形态”真正关键的是“如何在不使用A、B的情况下逼近这个形态”。如果从头到尾都把A、B当作已知量那这个“复现”就只复现了实验结果没有复现出方法本身的价值。2.2 把最优控制问题改写成Q函数参数估计数据驱动策略优化的支点是动作值函数Q(x,u)。在LQR框架下对于给定策略KQ(x,u)是关于[x;u]的二次型可以写成zHz其中z[x;u]H是一个(nm)×(nm)对称矩阵。为什么这个转换重要因为最优反馈增益可以直接从H里读出来。把H分块成[[H_xx, H_xu], [H_ux, H_uu]]对Q关于u求梯度并令其为零得到u-H_uu^(-1)H_ux x也就是KH_uu^(-1)H_ux这意味着系统模型A、B完全不出现在新策略的表达式里只要把H估计准策略改进就是一次矩阵运算。这个性质是整个数据驱动设计能够成立的支点。2.3 策略评估与策略改进最小二乘如何驱动迭代接下来要解决的是H怎么估计。对当前策略K采集数据时实际执行u_k-Kx_ke_k其中e_k是有界零均值探索噪声。考虑贝尔曼方程Q(z_k)l(x_k,u_k)Q(z_{k1},-Kx_{k1})这里右边第二项的输入必须是策略K在x_{k1}处选定的动作-Kx_{k1}而不是实际执行的动作。为什么要这样因为动作值函数的定义是“在某状态执行某动作后后续一直按照当前策略走”的代价。如果直接使用下一个时刻带噪声的实际输入等于把探索噪声也当成了策略的一部分理论框架就乱了。把Q的二次型形式代入方程左右两边关于H都是线性的于是可以通过线性最小二乘求解H。这个过程就是策略评估得到H之后再做策略改进得到新的K然后重新采集数据循环往复。这就是LSPI风格的数据驱动策略迭代也是我复现这篇长文时采用的主算法骨架。2.4 探索噪声、持久激励与收敛性直觉关于探索噪声e_k许多初学者会觉得“加噪声会把数据搞脏”其实恰恰相反。线性最小二乘要解一个满秩法方程需要数据在状态-输入空间充分“各方向都有信息”这就是持久激励条件。没有探索噪声u_k完全由x_k决定z_k落在低维子空间里法方程必然奇异H估计不出来。e_k的作用是迫使系统在被控的同时产生足够的动态信息。理论分析会给出PE条件的定量要求工程上更关心如何验证——你可以检查回归矩阵Phi的条件数条件数过大就是激励不足的明确信号。收敛性方面LQR二次型结构决定了策略迭代往往非常快十几轮就能逼近最优解。直觉上可以把它类比为用牛顿法解方程每轮策略评估相当于在当前点做局部线性化策略改进相当于沿着梯度走一步二次型目标函数让这个过程没有太多“地形障碍”。所以复现时你经常会看到误差曲线呈指数下降而不是缓慢爬行。3. Matlab复现代码架构与关键实现3.1 复现环境与文件组织我用的是Matlab R2021a不需要额外安装深度学习或强化学习工具箱纯脚本级代码。唯一可能用到的函数是dare控制系统工具箱它是用来计算真实K*做对比的如果你没有工具箱也可以用固定点迭代解Riccati方程或者直接在收敛后拿dare结果当参考。文件结构很简单main_lqr_dd.m是主脚本collect_data_on_policy.m负责数据采集fit_q_function.m做Q函数估计unpack_H.m辅助恢复H矩阵。模块化不是为了好看而是为了后面换系统、换采样策略时只改一个文件。3.2 主流程策略迭代循环怎么搭下面给出一份可以直接运行的主脚本。注意A、B只出现在数据采集和最终验证部分算法的学习模块完全没有用到它们。clear; clc; rng(42); % 仿真用模型仅在数据生成阶段和最终验证时使用 A [1.05, 0.1; 0, 1.0]; B [0; 1]; Qc eye(2); R 1; % 初始稳定策略 K [0.8, 0.5]; T 400; % 单次采样步数 sigma 0.05; % 探索噪声标准差 max_iter 30; tol 1e-6; for iter 1:max_iter [X, U] collect_data_on_policy(A, B, K, T, sigma); theta fit_q_function(X, U, K, Qc, R); H unpack_H(theta); K_new H(3,1:2) / H(3,3); err norm(K_new - K, inf); fprintf(iter %2d: K[%.4f, %.4f], err%.3e\n, iter, K_new, err); K K_new; if err tol break; end end % 用真实模型计算最优解仅供验证 P dare(A, B, Qc, R); K_star (R B*P*B) \ (B*P*A); disp(Learned K vs K*:); disp([K; K_star]);这段代码对应的是策略迭代的完整闭环采集数据、估计Q函数、提取新策略、判断收敛。实际运行时你会看到K从[0.8, 0.5]逐步逼近K*误差后面两三位小数快速往下掉。这里硬编码了n2、m1的结构方便讲原理扩展到更高阶系统时只需把H的索引逻辑改成通用的半向量化构造。3.3 数据采集模块生成“能学”的样本数据采集看似平淡却是整个复现中最容易翻车的地方之一。函数如下function [X, U] collect_data_on_policy(A, B, K, T, sigma) n size(A,1); m size(B,2); X zeros(n, T1); U zeros(m, T); x [0.8; -0.4] 0.1*randn(n,1); for k 1:T e sigma * randn(m,1); u -K*x e; X(:,k) x; U(:,k) u; x A*x B*u; end X(:,T1) x; end初值随机化是为了增加数据多样性。现场实验中做不到每一次都重置状态的话可以分成多段数据在不同工况、不同初始条件下采集然后拼接成同一批训练数据。不要小看这个细节采样轨迹的质量直接决定最小二乘能否成功。如果系统状态很快收敛到零点后面的样本基本是重复的信息量集中在初始阶段。解决办法是适当调小K或者增加初始随机扰动或者做多段不同初值的实验。3.4 Q函数估计与策略更新最小二乘核心代码核心估计函数如下它把贝尔曼方程直接翻译成线性回归function theta fit_q_function(X, U, K, Qc, R) T size(U,2); Phi zeros(T, 6); Gamma zeros(T, 1); for k 1:T xk X(:,k); uk U(:,k); xp X(:,k1); zk [xk; uk]; % 当前状态-实际输入 zp [xp; -K*xp]; % 下一状态-策略动作 Phi(k,:) phi_vec(zk) - phi_vec(zp); Gamma(k) xk*Qc*xk uk*R*uk; end theta Phi \ Gamma; end function p phi_vec(z) x1 z(1); x2 z(2); u z(3); p [x1^2; 2*x1*x2; 2*x1*u; x2^2; 2*x2*u; u^2]; end function H unpack_H(theta) H [theta(1), theta(2), theta(3); theta(2), theta(4), theta(5); theta(3), theta(5), theta(6)]; end6个独立参数对应n2、m1时H矩阵的对称结构原理上n维状态加m维输入需要(nm)(nm1)/2个参数。单步代价里用的是目标代价矩阵Qc和R这是LQR特有的便利。回归矩阵Phi的每一行是从当前样本特征中减去下一时刻策略特征对应贝尔曼方程。最小二乘用反斜杠运算符简单直接数据量不够时它会警告矩阵病态这个警告要重视。很多复现者会在这里走弯路Gamma里的u_k用的是实际执行值包含噪声而zp里的u_policy用的是纯策略输出不含噪声两者不能搞混。实际执行值决定系统下一状态纯策略输出是未来代价的“意图值”。混用后K会偏而且偏差方向很隐蔽往往不会被仿真稳定性发现却会让学习结果离真解差很远。3.5 收敛性验证拿什么证明学出来了验证阶段我会额外记录K_history和err_history并绘制收敛曲线figure; subplot(2,1,1); plot(K_history, LineWidth, 1.5); legend(K1, K2); yline(K_star(1), --); yline(K_star(2), --); grid on; title(Gain convergence); subplot(2,1,2); semilogy(err_history); grid on; title(Iteration residual);理想结果有两层一是K曲线收敛到K*所在位置二是残差曲线在semilogy坐标下接近直线下降。再用学习到的K做一次闭环仿真状态轨迹应该明显收敛初始策略下可能还有振荡学习后则平稳快速回到零点。我还会顺手检查一下闭环极点的位置保证学习后的系统稳定裕度不是靠运气拿到的。4. 参数调节与实操经验让代码稳定收敛4.1 一套可以直接抄的数值示例下面这组参数在我这里跑得非常稳迭代约10到15次就能收敛到K*附近。建议先用它建立基线再逐步调整参数推荐值备注A[1.05, 0.1; 0, 1.0]开环有一个单位圆外的特征值B[0; 1]单输入二阶系统Qceye(2)状态权重R1输入权重K0[0.8, 0.5]必须稳定闭环T400大于参数个数的20倍以上sigma0.05可调范围0.01到0.2max_iter30一般十几步就够tol1e-6残差阈值这套参数可以反复用因为A选成了开环临界不稳定系统数据驱动学习的效果一目了然。如果选一个本身就很稳定的系统可能初始策略已经不错学习收益不明显反而不容易看出算法有没有真正学对。4.2 探索噪声幅度怎么选激励与扰动之间找平衡噪声幅值的选择是数据驱动LQR里最微妙的手工活。太小回归矩阵条件数爆炸最小二乘解被噪声主导太大轨迹偏离工作点非线性效应显现甚至输入饱和而且代价中uRu权重会放大噪声影响。我的经验是先看R的数值和输入量级让sigma的平方约等于输入信号均方值的1%到5%。在这个示例里R1输入合理范围大约在正负1之间sigma取0.05到0.1就足够。判断是否合适的硬指标是cond(Phi*Phi)一般应该小于1e10。如果条件数已经大于1e12别急着改代码先加噪声试试。4.3 数据量、信息矩阵条件数与迭代轮次的配合LSPI对数据量要求不高但也不是够估计参数就行。回归矩阵的行数是T参数个数是6T50可能也能跑但估计方差大策略会抖动。我推荐T大于20倍参数个数也就是120以上考虑到系统动态时间尺度更需要覆盖若干倍上升时间。此外单组初值的样本往往在空间分布上不均匀建议多段数据拼接或者每轮采集前重置状态到一个随机点。如果数据采集来自连续过程无法重置就用递推最小二乘实时更新在线估计这部分偏扩展有兴趣可以后续再深入。轮次方面如果残差下降很慢先看是不是噪声太小或者K0离K*太远而不是盲目增加迭代次数。4.4 判断学习是否成功的几个硬指标怎么判断复现是成功的只看K收敛不够我给自己定了四个指标第一K_history相邻误差单调下降最终低于1e-6第二学习得到的K与模型已知时的K*误差在1e-3以内第三用学习到的K做闭环仿真所有闭环极点位于单位圆内状态不发散第四闭环累计代价低于初始策略下的代价。这四条都满足才算真正跑通了。如果K收敛但误差一直停在某个量级下不去大概率是符号约定或公式提取有问题而不是随机噪声造成的。5. 踩坑实录复现中常见的五个问题与排查方法5.1 初始稳定策略怎么找别让第一轮就发散数据驱动策略迭代不是随便从零开始的。K0必须保证闭环稳定否则采集的数据自带发散轨迹最小二乘出来全是噪音。实际中没有A、B怎么办可以用简化模型设计一个粗K0或直接现场手动整定一个能稳住系统的增益。仿真复现时可以先用模型验证eig(A-B*K0)1但算法流程里绝不能把A、B传进学习模块。我见过不少第一次接触这个方向的同学随便初始化Kzeros(1,2)结果第一轮数据轨迹就冲出去了后面再怎么调都没用。5.2 条件数爆炸激励不足的典型症状症状是fit_q_function能跑但theta绝对值大且符号不断跳变cond(Phi*Phi)1e12。原因基本两个噪声太小或者初始状态落在某些方向上使z_k长期线性相关。排查时先打印条件数把它降到1e8以下如果加噪声仍然不降考虑增加独立初始状态数。另一个细节数据里u_k-Kx_ke_k若K太大导致状态快速收敛后半段数据几乎没有信息量信息矩阵有效秩不够。此时可缩短每次采集的步数多段重置状态而不是单纯拉长T。5.3 策略更新后系统发散加一层“策略保护”我在迭代过程中遇到过某一轮K_new让闭环发散的情况原因是该轮数据量偏小或噪声序列偶然相关性较强H估计失真。保护办法很简单在下一轮采集前检查K_new的稳定裕度。仿真中有模型可以直接验极点和稳定裕度实际系统中可以用一段非常短的小信号试验测试闭环响应或者简单做“折半更新”K_new K_old 0.5*(K_new - K_old)降低一步走太远的风险。这个思路对应论文里的步长或正则化处理虽然理论收敛性会略打折扣但工程上很管用。5.4 符号与公式陷阱连续系统公式不能直接套复现时最容易翻车的是符号约定。离散还是连续、状态反馈是u-Kx还是uKx、代价里有没有1/2系数、KH_uu^(-1)H_ux还是加负号每个细节错了结果都会偏。我建议动手前先把自己用的符号约定写在注释里再对照论文伪代码逐行翻译。如果学习后的K与K*差一个符号或者差一个转置十有八九是这里的问题。这里多说一句TAC长文里的证明通常用连续时间表述而实际仿真大多用离散模型两条线的Riccati方程形式长得像但细节不同切忌照搬公式。5.5 常见问题速查表现象可能原因处理方法数据轨迹发散K0不稳定更换初始策略先稳住系统条件数大于1e12激励不足增大sigma增多初始状态检查TK发散或振荡单轮估计失真折半更新增加数据检查符号K收敛但偏K*连续离散混用或u符号反统一符号约定核对K提取公式收敛慢sigma过小或T过长导致信息冗余适当增大sigma缩短T并重置初值写到最后分享一点复现过程中的个人感受。TAC这类长文之所以难啃往往不是数学有多高深而是理论证明、算法流程和实验代码之间存在大量隐性约定你不亲手把代码跑起来很难发现“哦原来这里的数据要做这个预处理”。我复现这篇LQR数据驱动策略优化最大的收获不是多背了几条公式而是建立了一种直觉在模型不可用的场景里贝尔曼方程加参数化函数加最小二乘这套老组合依旧是可信赖的骨架。如果你的下一步是读更多数据驱动控制论文我强烈建议先把这篇的Matlab demo复现透再往鲁棒、约束、输出反馈方向扩展。踩过这些坑之后你会发现很多新论文本质上都是在同一套骨架上加工程细节。祝复现顺利。
