最近把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'PAP是离散代数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]单输入二阶系统Qceye(2)状态权重R1输入权重K0[0.8, 0.5]必须稳定闭环T400大于参数个数的20倍以上sigma0.05可调范围0.01到0.2max_iter30一般十几步就够tol1e-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，增多初始状态，检查TK发散或振荡单轮估计失真折半更新，增加数据，检查符号K收敛但偏K*连续离散混用或u符号反统一符号约定，核对K提取公式收敛慢sigma过小或T过长导致信息冗余适当增大sigma，缩短T并重置初值写到最后，分享一点复现过程中的个人感受。TAC这类长文之所以难啃，往往不是数学有多高深，而是理论证明、算法流程和实验代码之间存在大量隐性约定，你不亲手把代码跑起来，很难发现“哦，原来这里的数据要做这个预处理”。我复现这篇LQR数据驱动策略优化最大的收获，不是多背了几条公式，而是建立了一种直觉：在模型不可用的场景里，贝尔曼方程加参数化函数加最小二乘，这套老组合依旧是可信赖的骨架。如果你的下一步是读更多数据驱动控制论文，我强烈建议先把这篇的Matlab demo复现透，再往鲁棒、约束、输出反馈方向扩展。踩过这些坑之后，你会发现很多新论文本质上都是在同一套骨架上加工程细节。祝复现顺利。
最近把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'PAP是离散代数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]单输入二阶系统Qceye(2)状态权重R1输入权重K0[0.8, 0.5]必须稳定闭环T400大于参数个数的20倍以上sigma0.05可调范围0.01到0.2max_iter30一般十几步就够tol1e-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，增多初始状态，检查TK发散或振荡单轮估计失真折半更新，增加数据，检查符号K收敛但偏K*连续离散混用或u符号反统一符号约定，核对K提取公式收敛慢sigma过小或T过长导致信息冗余适当增大sigma，缩短T并重置初值写到最后，分享一点复现过程中的个人感受。TAC这类长文之所以难啃，往往不是数学有多高深，而是理论证明、算法流程和实验代码之间存在大量隐性约定，你不亲手把代码跑起来，很难发现“哦，原来这里的数据要做这个预处理”。我复现这篇LQR数据驱动策略优化最大的收获，不是多背了几条公式，而是建立了一种直觉：在模型不可用的场景里，贝尔曼方程加参数化函数加最小二乘，这套老组合依旧是可信赖的骨架。如果你的下一步是读更多数据驱动控制论文，我强烈建议先把这篇的Matlab demo复现透，再往鲁棒、约束、输出反馈方向扩展。踩过这些坑之后，你会发现很多新论文本质上都是在同一套骨架上加工程细节。祝复现顺利。