最近把TAC(控制领域公认的顶级期刊)上一篇关于LQR直接自适应学习的数据驱动策略优化长文完整复现了一遍,Matlab版本跑通之后,第一反应是:原来标题看上去那么硬核的东西,核心逻辑拆开之后竟然这么干净。论文要解决的问题一句话就能讲清楚:系统模型未知时,不通过模型辨识,直接利用在线数据学习线性二次调节器(LQR)的最优反馈控制律。这个问题的实际价值很大——实验室里我们习惯用精确模型设计控制器,但工程现场很多被控对象要么拿不到模型,要么模型参数一直在漂,这时候数据驱动策略优化就是一条非常现实的出路。如果你在做最优控制、自适应控制相关方向,或者正在被“模型失配”折磨,这篇复现笔记应该能帮你节省大量试错时间。我会把问题动机、数学原理、Matlab实现、参数调节和踩坑过程一五一十列出来。
1. 这篇TAC长文到底解决了什么问题
1.1 传统LQR设计的“模型依赖”病根在哪里
LQR本身是控制理论里最老牌的问题:线性系统,二次型代价,最优控制器是状态反馈u=-Kx,K直接由代价权重和系统矩阵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_{k+1}=Ax_k+Bu_k,代价是Σ(x'Qx+u'Ru),最优反馈为u=-Kx,其中
K=(R+B'PB)^(-1)B'PA
P是离散代数Riccati方程P=Q+A'PA-A'PB(R+B'PB)^(-1)B'PA的唯一正定解。这套理论非常成熟,模型已知时一条dare命令就能出结果。但我复现时始终提醒自己:K的表达式只是“答案的形态”,真正关键的是“如何在不使用A、B的情况下逼近这个形态”。如果从头到尾都把A、B当作已知量,那这个“复现”就只复现了实验结果,没有复现出方法本身的价值。
2.2 把最优控制问题改写成Q函数参数估计
数据驱动策略优化的支点是动作值函数Q(x,u)。在LQR框架下,对于给定策略K,Q(x,u)是关于[x;u]的二次型,可以写成z'Hz,其中z=[x;u],H是一个(n+m)×(n+m)对称矩阵。为什么这个转换重要?因为最优反馈增益可以直接从H里读出来。把H分块成
[[H_xx, H_xu], [H_ux, H_uu]]
对Q关于u求梯度并令其为零,得到u=-H_uu^(-1)H_ux x,也就是
K=H_uu^(-1)H_ux
这意味着系统模型A、B完全不出现在新策略的表达式里,只要把H估计准,策略改进就是一次矩阵运算。这个性质是整个数据驱动设计能够成立的支点。
2.3 策略评估与策略改进:最小二乘如何驱动迭代
接下来要解决的是H怎么估计。对当前策略K,采集数据时实际执行u_k=-Kx_k+e_k,其中e_k是有界零均值探索噪声。考虑贝尔曼方程
Q(z_k)=l(x_k,u_k)+Q(z_{k+1},-Kx_{k+1})
这里右边第二项的输入必须是策略K在x_{k+1}处选定的动作-Kx_{k+1},而不是实际执行的动作。为什么要这样?因为动作值函数的定义是“在某状态执行某动作后,后续一直按照当前策略走”的代价。如果直接使用下一个时刻带噪声的实际输入,等于把探索噪声也当成了策略的一部分,理论框架就乱了。把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*,误差后面两三位小数快速往下掉。这里硬编码了n=2、m=1的结构,方便讲原理;扩展到更高阶系统时,只需把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, T+1); 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(:,T+1) = 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(:,k+1); 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个独立参数对应n=2、m=1时H矩阵的对称结构,原理上n维状态加m维输入需要(n+m)(n+m+1)/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] | 单输入二阶系统 |
| Qc | eye(2) | 状态权重 |
| R | 1 | 输入权重 |
| K0 | [0.8, 0.5] | 必须稳定闭环 |
| T | 400 | 大于参数个数的20倍以上 |
| sigma | 0.05 | 可调范围0.01到0.2 |
| max_iter | 30 | 一般十几步就够 |
| tol | 1e-6 | 残差阈值 |
这套参数可以反复用,因为A选成了开环临界不稳定系统,数据驱动学习的效果一目了然。如果选一个本身就很稳定的系统,可能初始策略已经不错,学习收益不明显,反而不容易看出算法有没有真正学对。
4.2 探索噪声幅度怎么选:激励与扰动之间找平衡
噪声幅值的选择是数据驱动LQR里最微妙的手工活。太小,回归矩阵条件数爆炸,最小二乘解被噪声主导;太大,轨迹偏离工作点,非线性效应显现,甚至输入饱和,而且代价中u'Ru权重会放大噪声影响。我的经验是:先看R的数值和输入量级,让sigma的平方约等于输入信号均方值的1%到5%。在这个示例里R=1,输入合理范围大约在正负1之间,sigma取0.05到0.1就足够。判断是否合适的硬指标是cond(Phi'*Phi),一般应该小于1e10。如果条件数已经大于1e12,别急着改代码,先加噪声试试。
4.3 数据量、信息矩阵条件数与迭代轮次的配合
LSPI对数据量要求不高,但也不是够估计参数就行。回归矩阵的行数是T,参数个数是6,T=50可能也能跑,但估计方差大,策略会抖动。我推荐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传进学习模块。我见过不少第一次接触这个方向的同学,随便初始化K=zeros(1,2),结果第一轮数据轨迹就冲出去了,后面再怎么调都没用。
5.2 条件数爆炸:激励不足的典型症状
症状是fit_q_function能跑但theta绝对值大且符号不断跳变,cond(Phi'*Phi)>1e12。原因基本两个:噪声太小,或者初始状态落在某些方向上使z_k长期线性相关。排查时先打印条件数,把它降到1e8以下;如果加噪声仍然不降,考虑增加独立初始状态数。另一个细节:数据里u_k=-Kx_k+e_k,若K太大导致状态快速收敛,后半段数据几乎没有信息量,信息矩阵有效秩不够。此时可缩短每次采集的步数,多段重置状态,而不是单纯拉长T。
5.3 策略更新后系统发散:加一层“策略保护”
我在迭代过程中遇到过某一轮K_new让闭环发散的情况,原因是该轮数据量偏小或噪声序列偶然相关性较强,H估计失真。保护办法很简单:在下一轮采集前检查K_new的稳定裕度。仿真中有模型可以直接验极点和稳定裕度;实际系统中可以用一段非常短的小信号试验测试闭环响应,或者简单做“折半更新”,K_new = K_old + 0.5*(K_new - K_old),降低一步走太远的风险。这个思路对应论文里的步长或正则化处理,虽然理论收敛性会略打折扣,但工程上很管用。
5.4 符号与公式陷阱:连续系统公式不能直接套
复现时最容易翻车的是符号约定。离散还是连续、状态反馈是u=-Kx还是u=+Kx、代价里有没有1/2系数、K=H_uu^(-1)H_ux还是加负号,每个细节错了结果都会偏。我建议动手前先把自己用的符号约定写在注释里,再对照论文伪代码逐行翻译。如果学习后的K与K*差一个符号或者差一个转置,十有八九是这里的问题。这里多说一句:TAC长文里的证明通常用连续时间表述,而实际仿真大多用离散模型,两条线的Riccati方程形式长得像但细节不同,切忌照搬公式。
5.5 常见问题速查表
| 现象 | 可能原因 | 处理方法 |
|---|---|---|
| 数据轨迹发散 | K0不稳定 | 更换初始策略,先稳住系统 |
| 条件数大于1e12 | 激励不足 | 增大sigma,增多初始状态,检查T |
| K发散或振荡 | 单轮估计失真 | 折半更新,增加数据,检查符号 |
| K收敛但偏K* | 连续离散混用或u符号反 | 统一符号约定,核对K提取公式 |
| 收敛慢 | sigma过小或T过长导致信息冗余 | 适当增大sigma,缩短T并重置初值 |
写到最后,分享一点复现过程中的个人感受。TAC这类长文之所以难啃,往往不是数学有多高深,而是理论证明、算法流程和实验代码之间存在大量隐性约定,你不亲手把代码跑起来,很难发现“哦,原来这里的数据要做这个预处理”。我复现这篇LQR数据驱动策略优化最大的收获,不是多背了几条公式,而是建立了一种直觉:在模型不可用的场景里,贝尔曼方程加参数化函数加最小二乘,这套老组合依旧是可信赖的骨架。如果你的下一步是读更多数据驱动控制论文,我强烈建议先把这篇的Matlab demo复现透,再往鲁棒、约束、输出反馈方向扩展。踩过这些坑之后,你会发现很多新论文本质上都是在同一套骨架上加工程细节。祝复现顺利。