离散时间非线性系统自适应动态规划：值迭代与策略迭代的Matlab复现

最近把一篇关于离散时间非线性系统自适应动态规划ADP的经典文章复现做完了趁热打铁记录一下。这篇博客针对的是标题里那类“非线性值迭代自适应动态规划”研究目标是用Matlab实现离散时间非线性系统的策略迭代/值迭代算法。如果你正在做最优控制方向或者被论文里的HJB方程、递增迭代序列、容许控制这些概念弄得焦头烂额这篇笔记应该能帮你少走不少弯路。我会从算法原理出发把值迭代和策略迭代的数学逻辑讲清楚然后给出一套可以直接抄作业的Matlab实现方案最后分享几个我实测踩过的坑。1. 项目背景与核心问题拆解1.1 非线性最优控制为什么难HJB方程这堵墙先聊点背景。学过线性系统最优控制的都知道LQR给定一个线性系统Q和R一选Riccati方程一解最优状态反馈增益就出来了整个过程非常标准。但现实世界的系统几乎都是非线性的倒立摆、机械臂、电力系统、飞行器哪个都不是简单的dx/dt Ax Bu。对非线性系统做最优控制理论上的解法是把Hamilton-Jacobi-Bellman方程解出来得到一个关于状态的值函数。问题是HJB方程是一个非线性偏微分方程除了极少数特殊系统根本不存在解析解。动态规划倒是能从数值上逼近最优解但它需要对连续状态空间逐点求值每增加一个状态维度计算量就指数膨胀这就是所谓的“维数灾难”。ADP的出现就是冲着这个问题去的。它不再逐点求值而是用一个参数化的函数结构——多项式、RBF基函数或者神经网络——去逼近值函数把偏微分方程的求解问题转化为参数估计问题。这是整个ADP家族的核心思想也是这篇文章复现的基础。1.2 值迭代和策略迭代到底在迭代什么标题里同时出现了“值迭代”和“策略迭代”乍看有点乱实际是同一类算法的两个分支。值迭代的出发点是随便给一个初始值函数然后反复套用Bellman方程更新它每轮更新都隐式地改进一次策略。它的优点是无需初始稳定策略条件放宽很多。缺点是收敛速度通常比较慢尤其当初始值函数与真实最优值差得很远时前几十轮迭代基本都在“热身”。策略迭代则是从一条明确的策略出发交替执行策略评估和策略改进。策略评估是在固定策略下求解Bellman方程得到对应的值函数策略改进是根据新值函数更新控制律。它的收敛速度通常更快但代价是需要一个初始的稳定策略——也就是能镇定系统的控制律这一步往往不是白给的。下面这个表格直观对比了两者的差异复现之前最好心里有个数。-值迭代策略迭代初始条件任意非负值函数一个稳定的初始策略迭代方式直接更新值函数策略评估 策略改进交替收敛速度较慢较快实现难度较简单稍复杂典型风险逼近误差导致不收敛初始策略不稳定导致发散说句实在话论文里这两种方法经常放在一起讲不是要你二选一而是帮你理解同一个最优控制问题的两条求解路径。复现时先把值迭代跑通再改造成策略迭代理解起来会顺很多。2. 算法原理两大迭代机制的数学逻辑2.1 系统建模与性能指标回到基础。考虑一类离散时间非线性系统状态方程写成$$x_{k1} f(x_k) g(x_k)u_k$$其中$x_k \in \mathbb{R}^n$是状态$u_k \in \mathbb{R}^m$是控制输入。这里要求$g(x)$已知这在很多机械系统里是成立的因为执行器往往以线性方式注入控制量。性能指标用无限时域折扣或非折扣代价函数$$J(x_0, u) \sum_{k0}^{\infty} (x_k^T Q x_k u_k^T R u_k)$$其中$Q$和$R$是对称正定矩阵。最优控制问题的目标就是找到控制序列 $u_k h(x_k)$让这个代价最小。根据Bellman最优性原理最优值函数满足$$V^(x_k) \min_{u_k} { x_k^T Q x_k u_k^T R u_k V^(x_{k1}) }$$这就是离散时间的HJB方程。它把无限时域的最优性问题转化成一个函数方程。但$V^(x)$未知而且嵌套在$V^(x_{k1})$里没法直接解。ADP的思路是既然直接解不行那就一层一层迭代逼近它。2.2 值迭代更新方程与收敛逻辑值迭代的更新规则很直接。设第$i$轮的值函数为$V_i(x)$首先根据当前值函数求一个“贪心”控制$$h_i(x_k) \arg\min_{u_k} { x_k^T Q x_k u_k^T R u_k V_i(x_{k1}) }$$然后用这个控制去更新值函数$$V_{i1}(x_k) x_k^T Q x_k h_i^T(x_k) R h_i(x_k) V_i(x_{k1})$$从任意非负初始函数$V_0$出发在比较温和的条件下这个序列会单调收敛到最优值函数$V^*$。这就是为什么我前面说“不需要初始稳定策略”。实际仿真中初始$V_0$常直接设为零或者用一个随机的小权重向量。这样做的风险是前期迭代中策略质量很差状态轨迹可能发散但由于只是值函数在更新数值上并不会崩溃只要采样区域覆盖合理就能收回来。需要留意的是理论上的收敛性建立在“精确更新”的前提下。一旦引入函数近似每次更新的投影误差就会累积。所以复现时控制逼近误差很关键这也是后面调参的重点。2.3 策略迭代的评估-改进循环策略迭代走的是另一条路。假设当前策略为$h_i$第一步是策略评估把策略固定下来求解Bellman方程$$V_i(x_k) x_k^T Q x_k h_i^T(x_k) R h_i(x_k) V_i(x_{k1})$$这是一个关于$V_i$的函数方程比HJB方程简单因为它去掉了对$u$的极小化操作。第二步是策略改进利用新值函数更新控制$$h_{i1}(x_k) \arg\min_{u_k} { x_k^T Q x_k u_k^T R u_k V_i(x_{k1}) }$$两个步骤交替直到策略不再变化。这个流程的关键条件初始策略$h_0$必须是容许控制也就是能让系统在该策略下稳定同时使代价函数有限。如何获得这样一个初始策略工程上最省事的办法就是用LQR先算一个线性反馈作为起点只要系统的非线性项不是太强这个反馈通常能镇定系统。策略迭代与值迭代在“策略改进”这一步用的是同一个极小化操作核心区别在于值函数如何更新。前者需要先评估当前策略后者直接递推值函数。3. Matlab仿真环境搭建与模型参数确定3.1 复现需要准备的工具与环境复现这套算法不需要多新的环境。我用的版本是Matlab R2022b实测在R2018b上也能跑通。核心依赖只有几个基本矩阵运算、优化工具箱的fminbnd函数用于策略改进中求极小值、控制系统工具箱的dlqr函数用于生成初始线性策略。如果你连控制系统工具箱都没有手写Riccati方程迭代也完全可行后面我会提一句。神经网络工具箱不是必须的。虽然不少论文用BP神经网络做值函数逼近但我在复现时采用了“线性参数化基函数 最小二乘”的方案也就是把值函数写成$$V(x) W^T \phi(x)$$其中$\phi(x)$是选定的基函数向量$W$是待定权重。这个方案的好处是权重更新可以直接用最小二乘求解绕开了BP训练的收敛震荡问题。很多学生一上来就上工具箱被trainlm的学习率、动量项折腾半天其实完全没必要。对复现论文而言线性参数化是最稳的起点。3.2 仿真系统模型设计我采用的仿真对象是一个近线性系统配非线性扰动这种模型在很多ADP论文中出现过既能体现非线性特征又方便用LQR生成初始稳定策略。具体状态方程$$x_{k1} A x_k B u_k \Delta(x_k)$$其中$$A \begin{bmatrix} 0.8 0.1 \ 0.05 0.9 \end{bmatrix}, \quad B \begin{bmatrix} 0 \ 1 \end{bmatrix}, \quad \Delta(x) \begin{bmatrix} 0.1 \sin(x_1) \ 0.05 x_2 \sin^2(x_1) \end{bmatrix}$$性能指标中取$QI_{2\times 2}$$R1$。系统有两个状态单控制输入。$\Delta(x)$是带有界扰动项满足Lipschitz条件保证系统在非线性扰动下仍然有良好行为。采样范围取$x_1, x_2 \in [-2, 2]$网格步长0.2一共401个状态点这个规模在Matlab单次循环中大概几十毫秒能算完调试时非常友好。为什么选这个模型关键因素是$B$向量使得系统关于控制输入是线性的这样在策略改进步骤中对每个状态点求最优$u$的问题就退化成一个单变量极值问题fminbnd处理起来很稳定。如果系统在$u$上也是非线性就需要更复杂的优化手段复现难度会陡增。4. 核心代码实现与关键环节拆解4.1 基函数选取与值函数参数化线性参数化方案里基函数的选择直接决定逼近精度。我测试了两类基函数多项式基函数和高斯RBF基函数。多项式形式简单比如二次型加交叉项但在状态空间跨度大时容易产生数值不稳定。RBF形式更稳表达式为$$\phi_j(x) \exp \left( -\frac{|x - c_j|^2}{\sigma^2} \right)$$其中$c_j$是中心点$\sigma$是宽度参数。我在$[-2, 2]^2$区域均匀布了9个中心点3x3网格$\sigma$取1.2。这个配置下即使状态跑到采样边界之外RBF的输出也不会像多项式那样极速膨胀安全性好得多。权重个数上我用了9个基函数加上一个常数项总共10个参数。值函数写成矩阵形式$V_W(x) W^T \phi(x)$这个参数规模对当前二维系统来说完全够用。如果你复现的是更高维系统经验是基函数个数要大于状态维数的5到10倍中心点覆盖采样区域宽度不能太小否则容易造成“孤岛效应”。4.2 值迭代MATLAB核心实现先给出值迭代的主循环。初始化$V_0 0$对应$W$全零然后在每个采样状态点$X(:, j)$上求优化控制再更新值函数。% 采样状态点 [X1, X2] meshgrid(-2:0.2:2, -2:0.2:2); X [X1(:); X2(:)]; N size(X, 2); % RBF基函数初始化 centers [-1.5 -0.5 0.5 1.5 -1.5 -0.5 0.5 1.5 -0.5; -1.5 -1.5 -1.5 -1.5 0.5 0.5 0.5 0.5 1.5]; % 简化的9个中心 sigma 1.2; phi (x) [1; exp(-sum((x - centers).^2, 1) / sigma^2)]; % 值迭代参数 tol 1e-5; maxIter 500; errList zeros(maxIter, 1); W zeros(10, 1); for iter 1:maxIter V_new zeros(1, N); u_new zeros(1, N); for j 1:N x X(:, j); % 策略改进对当前状态求最优控制 fun (u) x*Q*x u*R*u W * phi(f(x) g(x)*u); u_opt fminbnd(fun, -5, 5); u_new(j) u_opt; V_new(j) fun(u_opt); end % 最小二乘更新权重W * PhiPhi ≈ V_new PhiMat zeros(10, N); for j 1:N PhiMat(:, j) phi(X(:, j)); end W_new (PhiMat * PhiMat) \ (PhiMat * V_new); % 收敛判据 err max(abs(W_new - W)); errList(iter) err; W W_new; if err tol fprintf(值迭代收敛于第 %d 轮误差 %e\n, iter, err); break; end end这里有两个关键细节需要展开说。第一fminbnd优化时边界选$[-5, 5]$这个边界要在仿真初期就测试好。边界太大浪费计算太小可能漏掉最优控制点。4.2中第二处关键是权重更新那行$W (\Phi \Phi^T)^{-1} \Phi V_{new}^T$这是典型的最小二乘正规方程。在实际运行中如果$\Phi \Phi^T$接近奇异需要用pinv代替反斜杠。我测试的9个RBF中心在网格上位置分布均匀矩阵条件数在可接受范围内直接用反斜杠没问题。这段代码运行完W就是值函数网络参数。接着用这个值函数做闭环仿真验证控制效果。4.3 策略迭代MATLAB核心实现策略迭代的代码结构有些不同核心在于策略评估和策略改进两阶段。初始策略我直接用dlqr基于线性部分$A, B$求$K_0$保证系统在初始策略下稳定。策略评估阶段需要解关于值函数权重的线性方程组。% 初始线性策略 K dlqr(A, B, Q, R); h_fun (x, K) -K * x; % 策略迭代参数 tol 1e-5; maxIterPI 100; W_pi zeros(10, 1); for iter 1:maxIterPI % ---- 策略评估固定当前策略求 W 满足 % V(x) xQx uRu V(f(x)g(x)u) A_mat zeros(N, 10); b_vec zeros(N, 1); for j 1:N x X(:, j); u h_fun(x, K); x_next f(x) g(x) * u; A_mat(j, :) phi(x) - phi(x_next); b_vec(j) x*Q*x u*R*u; end W_pi (A_mat * A_mat) \ (A_mat * b_vec); % ---- 策略改进更新控制增益这里用神经网络的局部线性化近似 K_new update_gain_from_value(X, W_pi, phi, f, g, Q, R, K); gain_change norm(K_new - K); K K_new; fprintf(策略迭代第 %d 轮增益变化 %e\n, iter, gain_change); if gain_change 1e-6 break; end end策略改进部分的实现我采用的思路是在每个采样点计算出贪心控制$u^(x) \arg\min_u { x^TQx u^TRu W^T \phi(f(x)g(x)u) }$得到一组$(x, u^)$数据然后用最小二乘拟合一个线性反馈$u -Kx$。这样做的好处是保持策略的线性结构策略迭代过程非常稳定。如果你希望策略提高表达能力可以用RBF网络去拟合任意非线性策略但要注意避免过拟合。4.4 主流程与闭环仿真验证学习完成后需要做闭环仿真验证。从初始状态$x_0 [1.5; -1]$出发用学到的控制律跑50个时间步记录状态和控制轨迹。% 闭环仿真 x0 [1.5; -1]; Tmax 50; x_traj zeros(2, Tmax1); u_traj zeros(1, Tmax); x_traj(:, 1) x0; for k 1:Tmax x x_traj(:, k); % 使用值迭代学到的策略min_u或策略迭代的线性反馈 fun (u) x*Q*x u*R*u W * phi(f(x) g(x)*u); u_k fminbnd(fun, -5, 5); u_traj(k) u_k; x_traj(:, k1) f(x) g(x) * u_k; end % 绘图 figure; subplot(2, 1, 1); plot(0:Tmax, x_traj(1, :), b-, LineWidth, 1.5); hold on; plot(0:Tmax, x_traj(2, :), r--, LineWidth, 1.5); grid on; legend(x1, x2); title(闭环状态轨迹); subplot(2, 1, 2); plot(0:Tmax-1, u_traj, k-, LineWidth, 1.5); grid on; title(控制轨迹);把闭环轨迹画出来这是最能直接反映算法好坏的图。5. 仿真结果分析与收敛性验证5.1 判断算法是否真正收敛了复现时最容易犯的错误是迭代误差降下去了就宣布收敛结果闭环仿真一塌糊涂。判断收敛不能只看单一指标至少要看三个信号。第一个信号是相邻迭代的权重差。程序中我用$err \max |W_{new} - W|$作为停止条件。这个值反映的是参数层面的稳定程度权重变化小说明值函数逼近基本饱和。第二个信号是值函数曲面形态。把学到的值函数画成三维曲面最优点应该是一个光滑的“漏斗”形状底部在原点附近没有任何异常的尖峰或者凹陷。如果曲面出现大范围震荡基本可以判定RBF宽度或中心点设置有问题。第三个信号是闭环状态轨迹。稳定收敛的控制策略应当让状态快速收敛到原点控制量不过大、不震荡。这三个信号同时满足才能说算法真的收住了。我在实验中发现有时权重误差已经很小但状态轨迹仍然有缓慢漂移原因是最优值函数在原点邻域内的曲率没有被充分刻画。解决办法是增加原点附近的采样点密度或者在基函数里加二次型硬编码项。5.2 三张必出的仿真图与判读方法复现论文时审稿人或者导师最关心的图一般有三类。第一类是收敛曲线图横轴迭代轮次纵轴权重误差或值函数误差这条曲线应当单调下降然后趋于平缓。第二类是最终的闭环状态轨迹图反映控制器对系统的镇定能力。第三类是控制输入图检查控制量是否在合理幅度内。如果做对比实验建议把值迭代和策略迭代的收敛曲线叠在一张图里。策略迭代通常会在更少的轮次内收敛但每一轮的计算量更大因为策略评估需要解一次线性方程组。值迭代虽然轮次多但单轮计算简单整体时间未必差很多。把两种算法的运行时间和收敛轮次记录下来对比会非常有说服力。5.3 参数敏感性哪些参数决定成败参数选择对ADP算法的影响是决定性的。我整理了一个快速参考表所有结论都来自实测。- 参数影响我的推荐值RBF中心分布覆盖采样区域中心3x3或4x4均匀网格$\sigma$宽度过小震荡过大模糊0.8~1.5按状态范围缩放采样范围过小学到局部解±2倍期望工作区域采样密度过低逼近误差大每维20个点左右$Q$影响状态收敛速度对角元素1~5$R$影响控制幅度0.1~1fminbnd边界过小截断最优解±5左右逐步放宽停止阈值过小浪费计算1e-5~1e-6这里面最容易出问题的参数是$\sigma$。$\sigma$太小时每个RBF中心只能影响极小区域基函数之间没有重叠值函数曲面会出现一个个“鼓包”$\sigma$太大时所有基函数几乎一样逼近能力直线下降。一个实用的调参方法是先固定中心点把$\sigma$从0.5开始逐渐增大观察拟合误差的变化找到误差平台期的起始点。6. 常见问题与调参经验6.1 值迭代收敛慢或不收敛怎么办值迭代收敛慢往往是三个原因初始值函数离最优解太远、采样范围过大导致无效点占比高、fminbnd的搜索边界太窄。我的经验是把采样区域分阶段考虑先用较粗网格跑通流程再用细网格精调。不要一上来就是密密麻麻的状态点几百个点循环一次都要半天排错效率太低了。另一个常见问题是不收敛。如果你看到误差曲线在某个数值附近反复震荡首先检查值函数更新里是不是忘了包含代价项$x^TQx$。这个错误非常隐蔽因为程序不会报错但迭代就是降不下去。其次检查最小二乘的矩阵条件数条件数飙升时用pinv代替反斜杠问题通常会立刻缓解。6.2 策略迭代初始策略不稳定的问题策略迭代对初始策略的稳定性极其敏感。我最初偷懒用全零控制作为初始策略结果策略评估阶段解出的值函数直接爆炸数值上出现NaN。后来老老实实先用LQR求线性反馈K0仿真才顺利跑通。如果你的系统非线性太强LQR无法镇定系统可以改用状态相关Riccati方程SDRE生成初始策略。SDRE逐点求解线性矩阵方程非线性项被当作“冻结”状态进入矩阵基本能处理大部分工程系统。再不行可以先人为给一个阻尼项比如$u -Kx - \alpha x^3$让系统先稳定下来再做迭代。6.3 网络逼近过拟合与欠拟合的表现线性参数化RBF网络也会过拟合或者欠拟合。欠拟合的表现是值函数曲面过于平滑细节特征丢失闭环轨迹收敛很慢。过拟合的表现是值函数在采样点上完全贴合但采样点之间剧烈跳变闭环仿真时控制量来回振荡。判断方法是随机抽出一部分状态点做验证集比较验证误差和训练误差。缓解过拟合的办法很简单增加采样点密度、适当增大$\sigma$、减少RBF中心数量。缓解欠拟合则相反。很多时候你只需要调整中心点的布局方式把网格中心改成“中间密、四周疏”的非均匀分布就能在不增加参数数量的前提下提升逼近效果。6.4 Matlab实现路上的三个坑最后说几个Matlab特有的坑每一个我都亲自踩过。第一个坑是关于fminbnd的这个函数默认使用黄金分割搜索要求目标函数在边界内是单峰的。我们的代价函数通常不是严格凸的如果最优解落在边界上fminbnd返回的结果可能不可靠。处理方式是把边界放宽或者在搜索后检查解是否在边界上如果是扩大边界重新搜索。第二个坑是反斜杠运算与伪逆的取舍。矩阵$\Phi\Phi^T$条件数过大时反斜杠结果数值不稳定。我建议在复现阶段一律改成pinv虽然稍微慢一点但换来的是稳定的收敛行为。第三个坑是要特别注意避免变量名覆盖内置函数。我曾经把存储状态网格的变量命名为data结果无意中覆盖了某工具箱的数据类引用程序到处报错。后来养成习惯所有自定义变量都用带前缀的名字比如xGrid、WCritic、KPolicy从此再没犯过这种低级错误。最后分享一个我个人的习惯。复现这类文章时我会刻意把值迭代和策略迭代放在同一个脚本里共享同一套基函数和采样点。这样做的好处是对比实验时只需要切换一个开关就能同时看到两种算法的收敛曲线和闭环效果。这套代码跑通之后换任何系统模型只需要改f和g函数基函数自动按新状态的取值范围重新生成中心。做ADP方向的研究生这套代码基本能应付绝大多数论文复现任务。