简介面向自动控制、机器学习领域的研究者与工程师鲁棒自适应动态规划仿真代码源自一篇2012年期刊论文用于解决含不确定性和非线性特性系统的在线控制与最优决策问题。压缩包共3个文件包含两个Matlab源文件与一个Markdown说明文档整体仅6KB轻量易得其中Matlab源文件分别对应原始与优化实现md文档对算法流程和仿真设置作了注释。已有1323人浏览学习适合学生、科研人员快速掌握自适应动态规划原理或对照论文直接复现与二次开发。代码覆盖系统建模、不确定性扰动建模、价值函数近似、策略迭代等核心环节支持离线训练与在线调节并附优化版本各模块注释清晰运行逻辑分明可帮助读者深入理解鲁棒动态规划的迭代机制也能便捷移植至无人机、机器人等非线性控制任务中。1. 鲁棒自适应动态规划仿真代码从“跑不起来”到“改一改就能用”的完整拆解这套《鲁棒自适应动态规划仿真代码.zip》第一眼看到的人多半会愣一下里面就俩.m文件加一个README看起来“简陋”得不像一篇Automatica 2012论文的官方开源。但真正跑过鲁棒自适应动态规划Robust Adaptive Dynamic Programming, RDP的人会明白论文级别仿真代码的最大价值恰恰在这——没有工业项目里那堆配置文件的干扰核心算法裸露在外适合做两件事一是把H∞鲁棒控制的迭代逻辑彻底看懂二是拿自己的系统模型替换进去几分钟内得到一个能收敛的仿真基线。适合的人群很明确做鲁棒控制、自适应动态规划方向的研究生以及想把ADP类算法落地到实际被控对象上的工程师。本文会从理论骨架、代码逐段拆解、参数设置一路写到最常见的翻车点和验证技巧。2. 理论先行RDP 怎么把 H∞ 鲁棒控制变成一场零和博弈2.1 标准动态规划的痛点维数灾和模型依赖经典动态规划的核心是Bellman最优性原理——把多阶段决策问题拆成一层层子问题从末端往前递推。这套理论本身非常优雅但到了控制工程里马上撞上两堵墙。第一堵墙是维数灾状态变量一多值函数在一个高维网格上的存储和计算代价指数膨胀离散网格根本没法用。第二堵墙更隐蔽标准DP要求你精确知道系统的数学模型包括状态转移方程的全部参数。代入到连续时间线性时不变系统里问题变成求解Algebraic Riccati EquationARE这已经算幸运的了因为至少还有解析求解路径。但如果系统带不确定性、外部扰动或者模型参数本身就是黑匣子标准DP的“精确模型”前提瞬间崩塌。这也解释了为什么ADP自适应动态规划在这类问题上会被反复提起ADP用函数近似器来表达值函数用在线数据迭代更新绕开了“必须有精确模型”这个前提。2.2 零和博弈视角把扰动当成对面的对手RDP处理的场景本质上是一个带外部扰动的连续时间线性系统。系统方程可以写成ẋ A x B₁ w B₂ u这个式子里的w就是扰动输入可能是外界干扰、建模误差、参数摄动。工程师的直觉是“把扰动压下去”但RDP换了个视角把扰动w看作一个试图让控制性能变差的“对手”控制器u是另一个玩家两者在玩一场零和博弈。控制器的目标是让某个性能指标最小化扰动则想把同一指标最大化。这样一来H∞控制的那种“最坏情况优化”思想就能自然嵌进动态规划框架里。原来标准的性能指标J ∫(xᵀQx uᵀRu)dt现在要改造成带博弈的形式J ∫(xᵀQx uᵀRu − γ²wᵀw)dt。注意这里γ是一个设计参数叫扰动衰减水平物理含义是系统对外部扰动的抑制能力。γ越小要求控制器把扰动压得越狠但如果γ小于系统本身能达到的极限博弈方程就无解了体现为代数黎卡提方程不存在正定解。2.3 策略迭代骨架评估与改进不停交替RDP在实现层面用的仍然是策略迭代Policy Iteration的经典骨架但每一轮里嵌入了鲁棒修正项。整个流程可以拆成三步第一步给定一个稳定的初始反馈增益K₀保证闭环系统A − B₂K₀是稳定的。这一步极其重要后面避坑章节会展开讲它为什么是最高频翻车点。第二步是策略评估。对当前增益Kᵢ求解修正后的Lyapunov方程Pᵢ(A − B₂Kᵢ) (A − B₂Kᵢ)ᵀPᵢ Q KᵢᵀRKᵢ γ⁻²PᵢB₁B₁ᵀPᵢ 0这个方程比标准Lyapunov方程多了最后一项γ⁻²PᵢB₁B₁ᵀPᵢ它的来源就是上一节说的博弈项。直观理解是在评估当前策略时同时把最坏情况扰动的影响算进去得到的是一个“保守”的值函数Pᵢ。第三步是策略改进Kᵢ₊₁ R⁻¹B₂ᵀPᵢ。这一步和线性二次型调节器LQR的增益更新公式形式一致但P的来源已经被鲁棒修正过了。三个步骤循环往复直到P的Frobenius范数变化小于预设阈值。这套迭代写起来不算复杂但真正落到代码里你会发现有几个细节直接决定收敛与否初值的稳定性、γ的取值、以及每一轮策略评估时用的是care函数还是自己写迭代。这些正是接下来拆代码时要重点盯的位置。3. 拆解 Jiang2012Automatica.m从状态方程到策略迭代的每一行3.1 系统模型与不确定性权矩阵的入口打开Jiang2012Automatica.m最前面的一段代码定义的就是被控对象。原实现里系统矩阵A、B、B₁、性能权重Q、R和γ值都集中在这个区域。为了讲清楚我重构了一个可读性更好的版本方便逐段对照。% 系统模型定义连续时间线性时不变系统 A [-0.5 1.0; -0.5 0.0]; % 系统矩阵二阶不稳定对象 B2 [0; 1]; % 控制输入矩阵 B1 [1 0; 0 1]; % 扰动输入矩阵 % 性能权重与扰动衰减水平 Q eye(2); % 状态权重对角线为1表示两个状态同等重要 R 0.1; % 控制权重越小说明控制器越“敢用”控制量 gamma 1.5; % 扰动衰减水平小于某个临界值才存在可行解这几行是整个仿真里最需要花心思理解的部分。B1矩阵代表外部扰动的入口通道在RDP问题里它直接影响修正项γ⁻²P B₁B₁ᵀP的形态。如果B1取成单位阵意味着两个状态通道都受到同等强度的扰动这对控制器来说是最恶劣的情形之一。Q和R的比例决定了控制性能和能耗之间的权衡R取0.1意味着控制器不太在乎控制量大小更容易把状态压到零要是把R调大到10控制器就会变得“惜力”状态收敛会明显变慢。gamma的取值则需要先试算后面避坑章节会说怎么快速判断一个gamma是否可行。3.2 策略评估与改进的完整循环核心的迭代循环在代码中段。这里的实现方式是把策略评估转换成一次care函数调用而不是手动做Lyapunov迭代。很多第一次接触这套代码的人会困惑为什么能用care——关键就是策略评估方程里那个γ⁻²P B₁B₁ᵀP项它是P的二次项正好落在代数黎卡提方程的标准形式里。% 初始稳定增益先用LQR凑一个“安全”的起点 K0 -lqr(A, B2, Q, R); P care(A - B2*K0, B1, Q K0*R*K0, -gamma^-2 * eye(size(B1, 2))); K K0; % 策略迭代主循环 for iter 1:200 Acl A - B2 * K; % 当前闭环系统矩阵 P_new care(Acl, B1, Q K*R*K, -gamma^-2 * eye(size(B1, 2))); K_new R \ B2 * P_new; % 策略改进更新增益 % 收敛判定值函数矩阵的变化量 if norm(P_new - P, fro) 1e-8 P P_new; K K_new; fprintf(在第 %d 步收敛\n, iter); break; end P P_new; K K_new; end这里需要解释几个关键位置。K0用lqr求解而不是随便给一个数是为了保证初始策略就稳定这个习惯应当成为条件反射。care函数的四个参数分别是闭环系统矩阵、扰动输入矩阵、二次型权重矩阵、以及“控制权重”矩阵——最后一个取负定矩阵−γ⁻²I数学上对应零和博弈中扰动玩家的“代价”这是RDP与普通LQR在代码上最直观的差异。for循环里P_new的计算顺序也值得注意先用当前K构造权重Q KᵀRK再解care得到P_new然后用P_new更新K。这对应理论部分的评估-改进交替。收敛判定用Frobenius范数而不是逐元素比较为的是对矩阵所有元素同时敏感。迭代次数上限取200实际一般十几步就能收敛设上限是为了防死循环。3.3 收敛结果的可视化与残差核对循环跑完后代码通常还会有一段后处理画出状态轨迹和值函数收敛曲线。更重要的一步是自己额外做一个残差校验。% 收敛后校验把P代回鲁棒Bellman方程看残差 P_final P; Acl_final A - B2 * K; residual norm(P_final * Acl_final Acl_final * P_final ... Q K*R*K gamma^-2 * P_final * B1 * B1 * P_final, fro); disp([Bellman残差 , num2str(residual)]); % 状态轨迹仿真 tspan [0 10]; x0 [1; -0.5]; [t, x] ode45((t, x) (Acl_final * x), tspan, x0); plot(t, x); grid on;这段代码里ode45只用了闭环状态转移矩阵因为仿真时我们观察的是确定性部分的状态响应扰动w被当作“对手”在策略评估时已经隐含处理过了不再显式加入。残差校验这一步我建议每一次改参数后都跑一遍它可能比收敛曲线更能暴露问题——如果残差在1e-6量级说明P确实是方程的解如果残差大到1e-2就要怀疑care调用时参数顺序填错了。值得再强调一次care的第三个参数是状态权重第四个参数是“扰动代价”这两个位置极易填反填反的直接后果是P对不上方程迭代猛一看是收敛的实则得到的是另一个问题的解。4. 跑通这包代码MATLAB 环境、参数与收敛判断4.1 解压、目录结构与环境检查拿到压缩包后先把文件解压到一个工作目录比如D:\RDP_Project纯英文路径MATLAB里中文路径偶尔会在绘图保存时出诡异问题。然后清空MATLAB工作区把目录加入Path确认能正确读取三个文件Jiang2012Automatica.m、Jiang2012Automatica_optimized_version.m、README.md。两个.m文件的关系需要注意原版更贴近论文推导过程的逐句复刻优化版则做了一系列针对MATLAB运行效率的改写。初次学习先跑原版因为它和论文的公式编号几乎一一对应确认理解后再跑优化版观察两者输出是否一致——这也是一种交叉验证。启动运行前命令行检查必要的工具箱。care函数位于Control System Toolboxlqr也在其中如果缺失后面两步根本走不动。检查方法很简单ver(control) which care which lqr如果返回空或者报“未定义”说明工具箱没装全去Add-On Explorer里补装再回来。ode45在MATLAB基础模块里不需要额外工具箱。4.2 运行主程序与预期输出直接在编辑器里运行Jiang2012Automatica.m正常情况下命令行窗口会输出“在第 X 步收敛”并弹出一张状态响应曲线图。我第一次跑这套程序时把收敛迭代步数和图中的曲线形态都记录了下来这些数据既用于确认仿真成功也便于之后改参数做对比。为了确认结果可信建议做一组对照实验分别运行原版和优化版把两者输出的最终增益K打印出来比较。如果两位小数位上完全一致说明两个版本逻辑一致如果出现可见差异优先怀疑优化版是否在某种条件下修改了收敛阈值。这是很常见的差异来源因为优化版往往会人为放宽迭代次数或调整判定阈值来换取速度。4.3 核心参数速查表在反复改参数的过程中我整理了一张参数速查表每个参数的影响范围都做了标注方便后续排查问题。参数默认值作用调参方向Qeye(2)状态权重决定状态收敛速度与稳态精度的优先级增大Q让状态更快回零R0.1控制权重限制控制量大小增大R使控制更平滑收敛变慢gamma1.5扰动衰减水平越小鲁棒性越强逐渐减小到临界值附近观察是否发散K0lqr求解初始稳定增益必须保证闭环稳定不手填用lqr或care初始化收敛阈值1e-8判定策略迭代是否终止过小会卡死过大则精度不足迭代上限200防止死循环的保护性参数不收敛时先查前面四个参数这张表里最值得反复品味的是gamma和K0两行。gamma的临界值直接取决于B1的维度与数值——如果B1是单位阵系统要同时抵抗两个通道的扰动临界gamma会明显偏大反过来如果B1退化为单个向量可行域会宽不少。K0则纯粹是稳定性的问题后面专门讲。5. 避坑记录鲁棒 ADP 仿真里最常翻车的六个位置5.1 初值K0不稳定仿真直接发散现象策略迭代第一轮就报错或者状态轨迹迅速发散到Inf控制量巨大到像在打摆子。原因K0必须保证闭环系统A − B₂K0是稳定矩阵。手填一个K0很容易碰上不稳定的情况特别是系统本来就有开环右半平面极点的时候。解决用lqr先解一次初始增益代码如下。lqr即使参数选得很随意返回的增益也能保证闭环稳定这是最省心的初始化方式K0 -lqr(A, B2, Q, R); % 任何Q、R正定组合下都闭环稳定5.2 care函数报错或返回复数矩阵现象care调用报“solution does not exist”警告或返回的P含有虚部。原因gamma取得太小修正黎卡提方程的正定解不存在本质上是系统在当前gamma下达不到所需的扰动衰减水平。解决把gamma增大到2或者3重新试确认能出正定P之后再逐步减小gamma逼近临界值。临界值本身有工程意义——它就是这个控制器结构下能达到的最强鲁棒性值得记录下来。5.3 收敛阈值过小导致迭代卡死现象循环一直不触发breakiter一直奔着200去。原因阈值1e-12高估了care数值解的精度相邻两次P的差异在某个迭代步之后不会继续单调下降而是进入数值噪声区间。解决改成1e-8或者1e-9。对绝大多数控制仿真来说这个精度已经远高于工程需求。如果非得追求更高精度先检查残差指标是不是已经到1e-7量级——到了就没必要继续跑。5.4 与论文结果对不上符号和转置的隐性错误现象自己改写的代码输出K和论文表格里的K差了正负号或者转置位置不对导致矩阵维度报错。原因策略改进公式K R⁻¹B₂ᵀP有些论文定义u −Kx有些定义u Kx符号体系不同。还有MATLAB中care的公式约定和教科书里常见的ARE写法有差异需要小心对照。解决以残差校验为准不要以“看起来像”为准。算出K之后代进P A_cl A_clᵀP Q KᵀRK γ⁻²P B₁B₁ᵀP看残差是不是接近零。残差说话论文符号只是参考。5.5 zip伪加密导致的解压失败现象从网盘下载后双击解压报错“文件损坏”或“密码错误”但压缩包明明没有密码。原因这类压缩包有时会被平台或中转工具加上伪加密标记CRC校验被改掉导致解压软件误判。伪加密不是真加密不需要密码是解压软件被文件头里的标记骗了。解决换用7-Zip打开它能识别伪加密并正常解压。如果7-Zip也报错用虚拟机里的老版本WinRAR再试一次。总之第一反应不应该是重新下载——先换工具大概率一步解决。5.6 优化版与原版输出不一致现象两个版本的最终K或P存在微小但可见的差异。原因优化版可能修改了策略评估的求解路径比如用低秩近似或放宽了收敛阈值也可能对状态方程做了坐标变换。解决先判断差异量级。如果相对误差在1e-6以下视作一致继续推进如果差异在1e-2量级把两个版本的P矩阵逐项打印出来看差异集中在哪些元素上重点检查是不是gamma或B1的赋值被优化版改掉了。6. 进阶把 RDP 移植到自己的系统并验证 H∞ 指标6.1 三步替换法从示例对象到自己的被控对象把示例系统换成自己的模型只需要动三处系统矩阵A、输入矩阵B2、扰动通道B1。替换之后照例全流程跑一遍。% 以风洞系统模型为例展示如何替换 A [0 1; -3 -2]; % 换成你自己的系统矩阵 B2 [0; 1.5]; % 控制通道 B1 [0.1 0; 0.1 0.2]; % 扰动通道反映实际干扰注入路径 Q diag([2, 1]); % 状态权重 R 0.05; gamma 2.0;替换之后第一件事不是直接跑策略迭代而是先算一下开环极点eig(A)确认系统本身的稳定性和振荡模态。然后检查能控性rank(ctrb(A, B2))是否等于状态维数——如果不可控后面的LQR初始化和策略迭代全部没有意义。这两个检查花不了几十行代码但能省掉后续大量的无用调试时间。6.2 验证鲁棒性闭环H∞范数难道真的低于gamma代码跑完只能说明“策略迭代收敛了”不能说明“控制器达成了预期的H∞性能”。真正的验证手段是计算闭环系统到扰动的H∞范数看是否严格小于gamma。可以用下面的代码验证% 构造闭环系统对扰动的传递函数 sys_eval ss(Acl_final, B1, eye(2), zeros(2, size(B1, 2))); hinf_gain norm(sys_eval, inf); fprintf(闭环H∞范数 %.4f, gamma %.4f\n, hinf_gain, gamma);这里norm(sys, inf)要在MATLAB R2018a以上版本才有直接调用方式。范数结果如果比gamma小说明控制器确实实现了声称的扰动抑制水平如果比gamma大说明策略迭代收敛到了“错误的解”检查care的第四参是否真的写成了−γ⁻²I。6.3 我固化下来的调试习惯从那以后我每拿到一套新的ADP仿真代码都强制自己先走一遍固定流程先用LQR求出初始稳定增益再检查care四参的正负号位置最后跑完必做残差校验和H∞范数校验。这套流程帮我挡下了无数次“看起来收敛但实际无效”的隐性错误尤其是刚换新模型的时候稳得一批。另外提醒一句gamma的临界值不是算出来的是试出来的。从大往小试每次衰减10%直到收敛失败那一步往回退一格这就是当前系统结构下的最低gamma。这个过程有点枯燥却比任何理论分析都直观。希望这套整理过的方法和踩坑记录能帮到你。本文还有配套的精品资源点击获取
鲁棒自适应动态规划仿真代码拆解与调试指南
简介面向自动控制、机器学习领域的研究者与工程师鲁棒自适应动态规划仿真代码源自一篇2012年期刊论文用于解决含不确定性和非线性特性系统的在线控制与最优决策问题。压缩包共3个文件包含两个Matlab源文件与一个Markdown说明文档整体仅6KB轻量易得其中Matlab源文件分别对应原始与优化实现md文档对算法流程和仿真设置作了注释。已有1323人浏览学习适合学生、科研人员快速掌握自适应动态规划原理或对照论文直接复现与二次开发。代码覆盖系统建模、不确定性扰动建模、价值函数近似、策略迭代等核心环节支持离线训练与在线调节并附优化版本各模块注释清晰运行逻辑分明可帮助读者深入理解鲁棒动态规划的迭代机制也能便捷移植至无人机、机器人等非线性控制任务中。1. 鲁棒自适应动态规划仿真代码从“跑不起来”到“改一改就能用”的完整拆解这套《鲁棒自适应动态规划仿真代码.zip》第一眼看到的人多半会愣一下里面就俩.m文件加一个README看起来“简陋”得不像一篇Automatica 2012论文的官方开源。但真正跑过鲁棒自适应动态规划Robust Adaptive Dynamic Programming, RDP的人会明白论文级别仿真代码的最大价值恰恰在这——没有工业项目里那堆配置文件的干扰核心算法裸露在外适合做两件事一是把H∞鲁棒控制的迭代逻辑彻底看懂二是拿自己的系统模型替换进去几分钟内得到一个能收敛的仿真基线。适合的人群很明确做鲁棒控制、自适应动态规划方向的研究生以及想把ADP类算法落地到实际被控对象上的工程师。本文会从理论骨架、代码逐段拆解、参数设置一路写到最常见的翻车点和验证技巧。2. 理论先行RDP 怎么把 H∞ 鲁棒控制变成一场零和博弈2.1 标准动态规划的痛点维数灾和模型依赖经典动态规划的核心是Bellman最优性原理——把多阶段决策问题拆成一层层子问题从末端往前递推。这套理论本身非常优雅但到了控制工程里马上撞上两堵墙。第一堵墙是维数灾状态变量一多值函数在一个高维网格上的存储和计算代价指数膨胀离散网格根本没法用。第二堵墙更隐蔽标准DP要求你精确知道系统的数学模型包括状态转移方程的全部参数。代入到连续时间线性时不变系统里问题变成求解Algebraic Riccati EquationARE这已经算幸运的了因为至少还有解析求解路径。但如果系统带不确定性、外部扰动或者模型参数本身就是黑匣子标准DP的“精确模型”前提瞬间崩塌。这也解释了为什么ADP自适应动态规划在这类问题上会被反复提起ADP用函数近似器来表达值函数用在线数据迭代更新绕开了“必须有精确模型”这个前提。2.2 零和博弈视角把扰动当成对面的对手RDP处理的场景本质上是一个带外部扰动的连续时间线性系统。系统方程可以写成ẋ A x B₁ w B₂ u这个式子里的w就是扰动输入可能是外界干扰、建模误差、参数摄动。工程师的直觉是“把扰动压下去”但RDP换了个视角把扰动w看作一个试图让控制性能变差的“对手”控制器u是另一个玩家两者在玩一场零和博弈。控制器的目标是让某个性能指标最小化扰动则想把同一指标最大化。这样一来H∞控制的那种“最坏情况优化”思想就能自然嵌进动态规划框架里。原来标准的性能指标J ∫(xᵀQx uᵀRu)dt现在要改造成带博弈的形式J ∫(xᵀQx uᵀRu − γ²wᵀw)dt。注意这里γ是一个设计参数叫扰动衰减水平物理含义是系统对外部扰动的抑制能力。γ越小要求控制器把扰动压得越狠但如果γ小于系统本身能达到的极限博弈方程就无解了体现为代数黎卡提方程不存在正定解。2.3 策略迭代骨架评估与改进不停交替RDP在实现层面用的仍然是策略迭代Policy Iteration的经典骨架但每一轮里嵌入了鲁棒修正项。整个流程可以拆成三步第一步给定一个稳定的初始反馈增益K₀保证闭环系统A − B₂K₀是稳定的。这一步极其重要后面避坑章节会展开讲它为什么是最高频翻车点。第二步是策略评估。对当前增益Kᵢ求解修正后的Lyapunov方程Pᵢ(A − B₂Kᵢ) (A − B₂Kᵢ)ᵀPᵢ Q KᵢᵀRKᵢ γ⁻²PᵢB₁B₁ᵀPᵢ 0这个方程比标准Lyapunov方程多了最后一项γ⁻²PᵢB₁B₁ᵀPᵢ它的来源就是上一节说的博弈项。直观理解是在评估当前策略时同时把最坏情况扰动的影响算进去得到的是一个“保守”的值函数Pᵢ。第三步是策略改进Kᵢ₊₁ R⁻¹B₂ᵀPᵢ。这一步和线性二次型调节器LQR的增益更新公式形式一致但P的来源已经被鲁棒修正过了。三个步骤循环往复直到P的Frobenius范数变化小于预设阈值。这套迭代写起来不算复杂但真正落到代码里你会发现有几个细节直接决定收敛与否初值的稳定性、γ的取值、以及每一轮策略评估时用的是care函数还是自己写迭代。这些正是接下来拆代码时要重点盯的位置。3. 拆解 Jiang2012Automatica.m从状态方程到策略迭代的每一行3.1 系统模型与不确定性权矩阵的入口打开Jiang2012Automatica.m最前面的一段代码定义的就是被控对象。原实现里系统矩阵A、B、B₁、性能权重Q、R和γ值都集中在这个区域。为了讲清楚我重构了一个可读性更好的版本方便逐段对照。% 系统模型定义连续时间线性时不变系统 A [-0.5 1.0; -0.5 0.0]; % 系统矩阵二阶不稳定对象 B2 [0; 1]; % 控制输入矩阵 B1 [1 0; 0 1]; % 扰动输入矩阵 % 性能权重与扰动衰减水平 Q eye(2); % 状态权重对角线为1表示两个状态同等重要 R 0.1; % 控制权重越小说明控制器越“敢用”控制量 gamma 1.5; % 扰动衰减水平小于某个临界值才存在可行解这几行是整个仿真里最需要花心思理解的部分。B1矩阵代表外部扰动的入口通道在RDP问题里它直接影响修正项γ⁻²P B₁B₁ᵀP的形态。如果B1取成单位阵意味着两个状态通道都受到同等强度的扰动这对控制器来说是最恶劣的情形之一。Q和R的比例决定了控制性能和能耗之间的权衡R取0.1意味着控制器不太在乎控制量大小更容易把状态压到零要是把R调大到10控制器就会变得“惜力”状态收敛会明显变慢。gamma的取值则需要先试算后面避坑章节会说怎么快速判断一个gamma是否可行。3.2 策略评估与改进的完整循环核心的迭代循环在代码中段。这里的实现方式是把策略评估转换成一次care函数调用而不是手动做Lyapunov迭代。很多第一次接触这套代码的人会困惑为什么能用care——关键就是策略评估方程里那个γ⁻²P B₁B₁ᵀP项它是P的二次项正好落在代数黎卡提方程的标准形式里。% 初始稳定增益先用LQR凑一个“安全”的起点 K0 -lqr(A, B2, Q, R); P care(A - B2*K0, B1, Q K0*R*K0, -gamma^-2 * eye(size(B1, 2))); K K0; % 策略迭代主循环 for iter 1:200 Acl A - B2 * K; % 当前闭环系统矩阵 P_new care(Acl, B1, Q K*R*K, -gamma^-2 * eye(size(B1, 2))); K_new R \ B2 * P_new; % 策略改进更新增益 % 收敛判定值函数矩阵的变化量 if norm(P_new - P, fro) 1e-8 P P_new; K K_new; fprintf(在第 %d 步收敛\n, iter); break; end P P_new; K K_new; end这里需要解释几个关键位置。K0用lqr求解而不是随便给一个数是为了保证初始策略就稳定这个习惯应当成为条件反射。care函数的四个参数分别是闭环系统矩阵、扰动输入矩阵、二次型权重矩阵、以及“控制权重”矩阵——最后一个取负定矩阵−γ⁻²I数学上对应零和博弈中扰动玩家的“代价”这是RDP与普通LQR在代码上最直观的差异。for循环里P_new的计算顺序也值得注意先用当前K构造权重Q KᵀRK再解care得到P_new然后用P_new更新K。这对应理论部分的评估-改进交替。收敛判定用Frobenius范数而不是逐元素比较为的是对矩阵所有元素同时敏感。迭代次数上限取200实际一般十几步就能收敛设上限是为了防死循环。3.3 收敛结果的可视化与残差核对循环跑完后代码通常还会有一段后处理画出状态轨迹和值函数收敛曲线。更重要的一步是自己额外做一个残差校验。% 收敛后校验把P代回鲁棒Bellman方程看残差 P_final P; Acl_final A - B2 * K; residual norm(P_final * Acl_final Acl_final * P_final ... Q K*R*K gamma^-2 * P_final * B1 * B1 * P_final, fro); disp([Bellman残差 , num2str(residual)]); % 状态轨迹仿真 tspan [0 10]; x0 [1; -0.5]; [t, x] ode45((t, x) (Acl_final * x), tspan, x0); plot(t, x); grid on;这段代码里ode45只用了闭环状态转移矩阵因为仿真时我们观察的是确定性部分的状态响应扰动w被当作“对手”在策略评估时已经隐含处理过了不再显式加入。残差校验这一步我建议每一次改参数后都跑一遍它可能比收敛曲线更能暴露问题——如果残差在1e-6量级说明P确实是方程的解如果残差大到1e-2就要怀疑care调用时参数顺序填错了。值得再强调一次care的第三个参数是状态权重第四个参数是“扰动代价”这两个位置极易填反填反的直接后果是P对不上方程迭代猛一看是收敛的实则得到的是另一个问题的解。4. 跑通这包代码MATLAB 环境、参数与收敛判断4.1 解压、目录结构与环境检查拿到压缩包后先把文件解压到一个工作目录比如D:\RDP_Project纯英文路径MATLAB里中文路径偶尔会在绘图保存时出诡异问题。然后清空MATLAB工作区把目录加入Path确认能正确读取三个文件Jiang2012Automatica.m、Jiang2012Automatica_optimized_version.m、README.md。两个.m文件的关系需要注意原版更贴近论文推导过程的逐句复刻优化版则做了一系列针对MATLAB运行效率的改写。初次学习先跑原版因为它和论文的公式编号几乎一一对应确认理解后再跑优化版观察两者输出是否一致——这也是一种交叉验证。启动运行前命令行检查必要的工具箱。care函数位于Control System Toolboxlqr也在其中如果缺失后面两步根本走不动。检查方法很简单ver(control) which care which lqr如果返回空或者报“未定义”说明工具箱没装全去Add-On Explorer里补装再回来。ode45在MATLAB基础模块里不需要额外工具箱。4.2 运行主程序与预期输出直接在编辑器里运行Jiang2012Automatica.m正常情况下命令行窗口会输出“在第 X 步收敛”并弹出一张状态响应曲线图。我第一次跑这套程序时把收敛迭代步数和图中的曲线形态都记录了下来这些数据既用于确认仿真成功也便于之后改参数做对比。为了确认结果可信建议做一组对照实验分别运行原版和优化版把两者输出的最终增益K打印出来比较。如果两位小数位上完全一致说明两个版本逻辑一致如果出现可见差异优先怀疑优化版是否在某种条件下修改了收敛阈值。这是很常见的差异来源因为优化版往往会人为放宽迭代次数或调整判定阈值来换取速度。4.3 核心参数速查表在反复改参数的过程中我整理了一张参数速查表每个参数的影响范围都做了标注方便后续排查问题。参数默认值作用调参方向Qeye(2)状态权重决定状态收敛速度与稳态精度的优先级增大Q让状态更快回零R0.1控制权重限制控制量大小增大R使控制更平滑收敛变慢gamma1.5扰动衰减水平越小鲁棒性越强逐渐减小到临界值附近观察是否发散K0lqr求解初始稳定增益必须保证闭环稳定不手填用lqr或care初始化收敛阈值1e-8判定策略迭代是否终止过小会卡死过大则精度不足迭代上限200防止死循环的保护性参数不收敛时先查前面四个参数这张表里最值得反复品味的是gamma和K0两行。gamma的临界值直接取决于B1的维度与数值——如果B1是单位阵系统要同时抵抗两个通道的扰动临界gamma会明显偏大反过来如果B1退化为单个向量可行域会宽不少。K0则纯粹是稳定性的问题后面专门讲。5. 避坑记录鲁棒 ADP 仿真里最常翻车的六个位置5.1 初值K0不稳定仿真直接发散现象策略迭代第一轮就报错或者状态轨迹迅速发散到Inf控制量巨大到像在打摆子。原因K0必须保证闭环系统A − B₂K0是稳定矩阵。手填一个K0很容易碰上不稳定的情况特别是系统本来就有开环右半平面极点的时候。解决用lqr先解一次初始增益代码如下。lqr即使参数选得很随意返回的增益也能保证闭环稳定这是最省心的初始化方式K0 -lqr(A, B2, Q, R); % 任何Q、R正定组合下都闭环稳定5.2 care函数报错或返回复数矩阵现象care调用报“solution does not exist”警告或返回的P含有虚部。原因gamma取得太小修正黎卡提方程的正定解不存在本质上是系统在当前gamma下达不到所需的扰动衰减水平。解决把gamma增大到2或者3重新试确认能出正定P之后再逐步减小gamma逼近临界值。临界值本身有工程意义——它就是这个控制器结构下能达到的最强鲁棒性值得记录下来。5.3 收敛阈值过小导致迭代卡死现象循环一直不触发breakiter一直奔着200去。原因阈值1e-12高估了care数值解的精度相邻两次P的差异在某个迭代步之后不会继续单调下降而是进入数值噪声区间。解决改成1e-8或者1e-9。对绝大多数控制仿真来说这个精度已经远高于工程需求。如果非得追求更高精度先检查残差指标是不是已经到1e-7量级——到了就没必要继续跑。5.4 与论文结果对不上符号和转置的隐性错误现象自己改写的代码输出K和论文表格里的K差了正负号或者转置位置不对导致矩阵维度报错。原因策略改进公式K R⁻¹B₂ᵀP有些论文定义u −Kx有些定义u Kx符号体系不同。还有MATLAB中care的公式约定和教科书里常见的ARE写法有差异需要小心对照。解决以残差校验为准不要以“看起来像”为准。算出K之后代进P A_cl A_clᵀP Q KᵀRK γ⁻²P B₁B₁ᵀP看残差是不是接近零。残差说话论文符号只是参考。5.5 zip伪加密导致的解压失败现象从网盘下载后双击解压报错“文件损坏”或“密码错误”但压缩包明明没有密码。原因这类压缩包有时会被平台或中转工具加上伪加密标记CRC校验被改掉导致解压软件误判。伪加密不是真加密不需要密码是解压软件被文件头里的标记骗了。解决换用7-Zip打开它能识别伪加密并正常解压。如果7-Zip也报错用虚拟机里的老版本WinRAR再试一次。总之第一反应不应该是重新下载——先换工具大概率一步解决。5.6 优化版与原版输出不一致现象两个版本的最终K或P存在微小但可见的差异。原因优化版可能修改了策略评估的求解路径比如用低秩近似或放宽了收敛阈值也可能对状态方程做了坐标变换。解决先判断差异量级。如果相对误差在1e-6以下视作一致继续推进如果差异在1e-2量级把两个版本的P矩阵逐项打印出来看差异集中在哪些元素上重点检查是不是gamma或B1的赋值被优化版改掉了。6. 进阶把 RDP 移植到自己的系统并验证 H∞ 指标6.1 三步替换法从示例对象到自己的被控对象把示例系统换成自己的模型只需要动三处系统矩阵A、输入矩阵B2、扰动通道B1。替换之后照例全流程跑一遍。% 以风洞系统模型为例展示如何替换 A [0 1; -3 -2]; % 换成你自己的系统矩阵 B2 [0; 1.5]; % 控制通道 B1 [0.1 0; 0.1 0.2]; % 扰动通道反映实际干扰注入路径 Q diag([2, 1]); % 状态权重 R 0.05; gamma 2.0;替换之后第一件事不是直接跑策略迭代而是先算一下开环极点eig(A)确认系统本身的稳定性和振荡模态。然后检查能控性rank(ctrb(A, B2))是否等于状态维数——如果不可控后面的LQR初始化和策略迭代全部没有意义。这两个检查花不了几十行代码但能省掉后续大量的无用调试时间。6.2 验证鲁棒性闭环H∞范数难道真的低于gamma代码跑完只能说明“策略迭代收敛了”不能说明“控制器达成了预期的H∞性能”。真正的验证手段是计算闭环系统到扰动的H∞范数看是否严格小于gamma。可以用下面的代码验证% 构造闭环系统对扰动的传递函数 sys_eval ss(Acl_final, B1, eye(2), zeros(2, size(B1, 2))); hinf_gain norm(sys_eval, inf); fprintf(闭环H∞范数 %.4f, gamma %.4f\n, hinf_gain, gamma);这里norm(sys, inf)要在MATLAB R2018a以上版本才有直接调用方式。范数结果如果比gamma小说明控制器确实实现了声称的扰动抑制水平如果比gamma大说明策略迭代收敛到了“错误的解”检查care的第四参是否真的写成了−γ⁻²I。6.3 我固化下来的调试习惯从那以后我每拿到一套新的ADP仿真代码都强制自己先走一遍固定流程先用LQR求出初始稳定增益再检查care四参的正负号位置最后跑完必做残差校验和H∞范数校验。这套流程帮我挡下了无数次“看起来收敛但实际无效”的隐性错误尤其是刚换新模型的时候稳得一批。另外提醒一句gamma的临界值不是算出来的是试出来的。从大往小试每次衰减10%直到收敛失败那一步往回退一格这就是当前系统结构下的最低gamma。这个过程有点枯燥却比任何理论分析都直观。希望这套整理过的方法和踩坑记录能帮到你。本文还有配套的精品资源点击获取
