简介面向离散控制与最优控制研究者的Matlab代码包聚焦Hamilton-Jacobi-Bellman方程的牛顿策略迭代求解并采用Tensor Train格式缓解高维状态空间的计算与存储压力。包内共15个文件以13个.m函数与脚本为主体另有README说明和许可证文件压缩包约24KB体量轻便、结构紧凑便于按求解器、数值测试与辅助函数逐层阅读。代码提供一维和二维Allen-Cahn方程、Fokker-Planck方程等数值实验入口可通过umax参数开启控制约束并配有交互式参数输入与默认值提示适合复现实验、对照算法流程或改造为自身问题。读者可从中观察牛顿策略迭代如何配合TT线性求解与交叉近似理解谱离散和Legendre基下微分矩阵的构造并借助测试脚本验证精度与效率。已有1200人学习对具备Matlab基础、希望理解TT近似与策略迭代实现细节的研究生和工程人员较友好也能作为高维HJB数值求解的参考实现。1. 一张 550 GB 的表HJB 方程为什么在六维就“算不动”把一个六维离散控制问题铺在每维 64 点的规则网格上值函数就是一张 64^6 ≈ 6.9×10^10 个元素的表双精度存下来要 550 GB——光是分配内存这一步就失败了。同一道题换成 Tensor Train 格式占不到 1 MB配上牛顿策略迭代跑六轮残差就压到 1e-10 量级。差距不在算力而在一个被默认接受的前提我们一直假设 HJB 方程的解必须“查表”。可这个解本身是光滑的又被某种可分离的动力学和代价约束着它在张量展开下的秩远小于维度。这套系统开源的 MATLAB 实现把这条路径走通了而想搞懂它为什么 work得先把 Bellman 不动点、牛顿法和张量分解三者的接缝看清楚。2. 从 Bellman 不动点到牛顿步策略迭代凭什么只用 6 次2.1 离散控制问题的三件套动力学、代价、折扣先把问题摆正。离散时间、无限时域、折扣型最优控制的标准形式是三样东西状态转移 x_{k1} f(x_k, u_k)单步代价 g(x_k, u_k)折扣因子 γ ∈ (0,1)。控制量取自一个有限集合 U {u¹, …, u^m}也就是所谓“离散控制”——不是连续优化出来的而是从 m 个候选里挑一个。目标是最小化累计折扣代价 Σ γ^k g(x_k, u_k)对应的值函数满足 Bellman 最优方程V*(x) min_{u∈U} { g(x,u) γ·V*(f(x,u)) }这个式子里有两个让数值方法难受的地方。第一min 是非线性的而且不可微——它是一族仿射函数的逐点下确界。第二f(x,u) 一般不在网格点上所以 V*(f(x,u)) 必须靠插值取值。要在计算机上真解它只能把状态空间离散成乘积网格 X X₁×…×X_d每维 n_i 个节点然后在这个网格上定义值函数数组 V(i₁,…,i_d)总规模 N Π n_i。这就是 550 GB 的来源也是所有后续设计要绕开的障碍。顺带说一句网格不能随便省。策略改进那一步是逐点取 min如果网格太粗最优策略会在相邻网格点之间来回切换得到一个物理上不存在的“抖振”解。所以离散化误差和后面的张量截断误差必须匹配着调这个坑在第 5 章会专门讲。2.2 策略评估是一个线性方程不是一次迭代Bellman 方程里的 min 让它变成非线性不动点问题但只要把策略 π: X → U 固定住整个方程立刻退化成线性的V_π g_π γ·P_π·V_π ⟺ (I − γP_π)·V_π g_π这里 g_π 是“按策略 π 取的那个动作对应的单步代价”P_π 是策略诱导的转移算子。P_π 的每一行对应一个网格点非零元落在 f(x, π(x)) 周围的插值节点上——如果做 d 线性插值每行最多 2^d 个非零元。这个线性系统有两个非常好的性质。其一是解唯一P_π 是广义随机矩阵谱半径不超过 1所以 γP_π 的谱半径 ≤ γ 1(I − γP_π) 必然可逆。其二是条件数可控||(I − γP_π)^{-1}||_∞ ≤ 1/(1−γ)γ 0.95 时大概 20 量级。这一点极其关键——很多张量线性求解器DMRG 式扫描、交替最小二乘在病态问题上会崩而策略评估天生良态这几乎是“牛顿 TT”这个组合能成立的前提条件。策略改进则是另一回事π(x) argmin_{u∈U} { g(x,u) γ·V_π(f(x,u)) }。对每个网格点算 m 个候选值取最小纯逐点操作代价 O(m·N)。在稠密表示下这一步和策略评估一样贵在 TT 表示下它反而成了最麻烦的一步因为它不保秩。2.3 牛顿视角二次收敛从哪里来又在哪里失效把策略迭代看成牛顿法这个视角不是类比是严格成立的。定义算子 F(V) V − T(V)其中 T 是 Bellman 最优算子求解 F(V) 0。T 不可微min 的锅但它是半光滑的取它的广义雅可比恰好就是某个策略对应的 (I − γP_π)。于是牛顿步 V_{k1} V_k − J^{-1}F(V_k) 展开后正是“策略评估 策略改进”这两步。在线性二次型问题里这就是经典的 Kleinman 迭代收敛性有闭式结论二次收敛。二次收敛有多值钱对比一下就明白。值迭代就是 Picard 迭代 V ← T(V)压缩因子 γ每轮误差乘 γ。要把残差从 1 压到 1e-10γ 0.99 时需要 log(1e-10)/log(0.99) ≈ 2300 轮γ 0.999 时需要两万多轮。策略迭代一般 5 到 8 轮就到机器精度。代价是每轮要从一次矩阵向量乘升级成一次线性求解成本比大约 10 到 30 倍——即使这样算总账仍然是两个数量级的净胜。但二次收敛有个前提每一步的线性系统要解得足够准。牛顿法的收敛速率对线性求解精度极度敏感。经典结论是若内层求解的相对残差为 η_k则外层满足 ||F_{k1}|| ≤ (C·η_k)·||F_k||——η_k 是常数退化成线性收敛η_k 随 ||F_k|| 一起趋于零才保住二次收敛。在 TT 格式里内层精度不由浮点决定而由截断容差决定这就是整个实现里最容易翻车的地方第 4 章会专门处理。3. Tensor Train把 n^d 的表拆成 d 个三阶核3.1 TT 秩到底是什么从矩阵展开看Tensor Train 分解把一个 d 阶张量写成 d 个三阶核的连乘V(i₁,…,i_d) G₁(i₁)·G₂(i₂)···G_d(i_d)其中 G_k 的尺寸是 r_{k−1} × n_k × r_k边界上 r₀ r_d 1保证连乘出来是个标量。存储量从 Π n_k 降到 Σ r_{k−1}·n_k·r_k。用前面那组数验一下d 6、n 64、所有内部秩取 20存储是 2×(1×64×20) 4×(20×64×20) ≈ 1.05×10^5 个双精度数约 840 KB。相对 550 GB压缩比接近 65 万倍。TT 秩不是人为设定的参数它有明确的代数含义r_k 等于张量按前 k 个维度展开成矩阵 V_{(1..k),(k1..d)} 的秩。这句话解释了很多现象。完全可分离的函数比如 V v₁(i₁)·v₂(i₂)···v_d(i_d)所有展开秩都是 1TT 表示退化成 d 个向量零误差。s 个可分离函数之和秩不超过 s。而一般的光滑函数展开秩随维度指数增长——但增长得“慢”在中等维度上往往只有几十。值函数为什么恰好落在“低秩”这一侧我一般这么理解折扣最优控制的值函数是高维空间里的一个光滑曲面它的非可分离性来自状态之间的耦合强度。弱耦合的链式系统第 k 维只跟第 k±1 维有交互展开秩与系统长度无关是个常数全耦合的强非线性系统秩会随维度指数爆炸这时候 TT 也救不了。所以 TT 不是万能压缩它是一个建模假设——假设这个问题的解其交互阶数是有限的。3.2 TT 上的四类基本运算与复杂度账实际写代码时真正用到的运算只有四类但要清楚各自对秩做了什么。第一类是加法。两个 TT 张量相加逐核做直和拼接秩直接相加r_k ← r_k^A r_k^B。所以加法必须紧跟一次取整否则秩会一路膨胀。第二类是取整rounding。用 TT-SVD 或者正交化之后的奇异值截断把秩从 R 压回 r相对误差控制在给定容差 ε。代价量级 O(d·n·r²·R)是全流程里最常被调用、也最贵的一步。所有“秩失控”的问题最终都靠它兜底。第三类是 TT-矩阵也叫 MPO乘 TT-向量。算子写成 A(i₁..i_d; j₁..j_d) A₁(i₁,j₁)···A_d(i_d,j_d)做一次乘法后秩变成两个秩的乘积 r_A·r_V然后必须取整。这一步是策略评估内层迭代的核心开销。第四类是逐点非线性运算比如取 min、取绝对值。这一类没有代数上的封闭形式——把两个小秩 TT 张量逐点取 min结果的秩既不是和也不是积理论上可以很大。这是第 4 章要单独处理的问题。把账算清楚一次 TT 矩阵向量乘加取整复杂度在 O(d·n·r_A²·r_V²) 到 O(d·n·r³) 之间关键结论是它对 d 是线性、对 n 是线性或线性对数而稠密表示对 d 是指数。这就是整套方法能跑起来的全部理由。3.3 转移算子为什么能写成小秩 TT-矩阵值函数能压成 TT只是成功了一半。如果 (I − γP_π) 这个算子本身必须是稠密的 N×N 矩阵那压缩值函数毫无意义——算子还是存不下。所以真正支撑整个方案的第二根支柱是转移算子本身也有小秩结构。先看没有策略耦合的理想情况。如果动力学的第 k 个分量只依赖 x_k 和 u即 f(x,u) (f₁(x₁,u), …, f_d(x_d,u))那么对固定的动作 ud 线性插值算子在每一维上是一个稀疏的双对角矩阵 P_k整体算子精确等于 Kronecker 积 P₁ ⊗ … ⊗ P_d。Kronecker 积在 TT 里的秩恒为 1。写出来是一个很朴素的单维插值函数function [lo, w] interp_weights(grid, y) % grid: 单维单调网格升序 % y : 查询点落在 [grid(1), grid(end)] 内 % lo : 所在单元的左端点下标 % w : 线性插值权重满足 y (1-w)*grid(lo) w*grid(lo1) n numel(grid); lo discretize(y, grid); lo min(max(lo, 1), n-1); % 夹到边界防止越界索引 w (y - grid(lo)) ./ (grid(lo1) - grid(lo)); % 归一化局部坐标 end这个函数看着平平无奇但它是可分离性的落脚点。d 维插值时角点 (i₁,…,i_d) 的权重恰好是 Π_k weight_k即每维权重相乘——这正是 Kronecker 积的定义。所以只要动力学满足“第 k 维只依赖第 k 维”转移算子就是秩 1 的 TT-矩阵无论 d 多大存储都是 O(d·n)。再看策略耦合。策略 π(x) 依赖全部坐标不同网格点选不同动作算子变成至多 m 个“动作掩码 × Kronecker 积”之和秩上界是 m。如果 f_k 依赖的是一个小邻域块比如第 k−1、k、k1 维这块的联合维度是 B那么秩上界是 B 的函数通常写成 m 乘以块大小的若干次幂。注意这里藏着整套方法最硬的假设。动力学耦合得越厉害P_π 的 TT 秩越高高到某个点就吃光全部内存优势。看到有人在十维以上跑通这套代码先别急着佩服算法去看看他的动力学是不是链式的。4. 牛顿步的 TT 落地线性求解器、残差与秩控制4.1 交替最小二乘把全局线性系统拆成一串小系统策略评估要解 (I − γP_π)V g_π。在 TT 格式下全局系统有 N 个未知量但未知量之间有 TT 结构约束于是标准做法是交替最小二乘ALS固定除第 k 个核之外的所有核此时残差范数 ||A·V − b||² 关于第 k 个核的元素是一个二次型令其导数为零就得到一个局部线性系统M_k · vec(G_k) rhs_k其中 M_k 是把 A 的核与左右环境张量收缩后得到的“有效算子”规模是 (r_{k−1}·n_k·r_k) × (r_{k−1}·n_k·r_k)。这个规模与维度 d 无关只与局部秩和单维网格数有关r 20、n 64 时约 25600 阶用共轭梯度解足够r 8 时只有 4096 阶直接分解都行。从左到右扫一遍、再从右到左扫一遍构成一次 sweep。关键实现细节是正交化。扫描前要把当前张量“左正交化”或“右正交化”让未参与优化的那部分核满足正交条件否则 M_k 的条件数会随着扫描迅速恶化解出来的核数值上互相抵消。这就是下面这段骨架里 tt_orthogonalize 的作用也是很多手写实现最先忽略、后来靠血泪经验补上的一步。function [V, res] tt_dmrg_solve(A, b, V, opts) % A : TT-矩阵即 (I - gamma*P_pi)按 3.3 的可分离结构构造 % b : TT-向量即阶段代价 g_pi % V : 初值直接用上一轮策略评估的解做热启动 % opts: .tol 内层相对残差, .maxRank 秩上限, .maxSweep 扫描上限 d numel(V); res inf; for it 1:opts.maxSweep for k 1:d-1 [L, R] tt_orthogonalize(V, k); % 左右正交化稳定局部系统 theta local_solve(A, b, L, R, k, opts); % 解局部线性系统 [V{k}, V{k1}] tt_split(theta, opts.tol, opts.maxRank); % SVD 拆回 end res tt_residual(A, V, b) / tt_norm(b); % 残差也在 TT 里算 if res opts.tol, break; end endlocal_solve 是全部开销的所在它要用 L 和 R 把 A 的中间核收缩成有效算子规模只与局部秩有关tt_split 里的 SVD 同时完成两件事——把合并核拆回两个核以及把超过 opts.maxRank 的奇异值砍掉。这一步决定了秩能不能长也决定了截断误差有多大。4.2 单站点会卡死双站点才能长秩上面这段骨架是“双站点”版本因为 local_solve 合并了第 k 和第 k1 个核。这个选择不是风格问题是能不能收敛的问题。单站点变体每次只优化一个核其余全固定。此时 V 的所有 TT 秩是写死的常量优化在秩固定的流形上进行。如果真解的秩比你设的上限高扫描会在流形上的一个局部极小处停住——残差降一两个数量级之后就不动了加扫描次数也没用。这就是观察到的“扫了 200 轮只降一个数量级”。双站点把相邻两个核合并成一个 (r_{k−1}) × (n_k·n_{k1}) × (r_{k1}) 的大核解完局部系统再 SVD 拆开。拆开时奇异值的个数可以比原来多——秩在一次扫描中最多能涨到 min(r_{k−1}·n_k, n_{k1}·r_{k1})。也就是说秩是“自适应长出来的”不是一开始猜的。再进一步是 AMEn 式的残差增强当扫描停滞时往局部基里补几个方向——通常取当前残差张量在局部的一小块或者几步 Krylov 方向。补完之后再做一次双站点优化等价于给秩加了一点“有依据的”增量比盲目加秩收敛快得多。这套代码里如果只有纯双站点而没做增强在 γ 接近 1 的时候会明显吃力。4.3 非精确牛顿把截断容差挂到残差上回到 2.3 那个结论外层收敛速率由内层精度 η_k 决定。问题在于 TT 求解器的精度不是一个连续可调的旋钮它由截断容差 tol 控制而 tol 直接决定了每次取整丢多少信息、秩能压到多低、一次策略评估要扫多少轮。tol 设太松外层退化成线性收敛设太紧内层扫描次数暴涨白白烧时间。工程上的解法是 Eisenstat-Walker 型的自适应 forcing term让内层容差随当前外层残差一起收紧。% 非精确牛顿外层策略迭代内层 TT 线性求解 for k 1:maxNewton r_k tt_norm(V - bellman_operator(V, policy)); % 当前 Bellman 残差 eta min(0.1, sqrt(r_k)); % forcing term opts.tol max(eta * r_k, 1e-12); % 内层只解到这个精度 V tt_dmrg_solve(A_pi, g_pi, V, opts); % 策略评估热启动 policy policy_improve(V); % 策略改进逐点取 min end几个参数值得说清楚。eta 的上界取 0.1 是经验值超过 0.1外层基本退化成线性收敛低于 1e-3内层要多扫好几轮而外层收益有限。sqrt(r_k) 的形式让 eta 在远离解时保持较大省扫描接近解时快速收紧保住二次收敛。opts.tol 的下限 1e-12 是浮点地板再往下压没有意义只会让双站点 SVD 反复切掉噪声奇异值。热启动这一行也很要紧——把上一轮的 V 直接当内层初值扫描次数通常能砍掉一半以上因为相邻两次策略评估的解差别不大。不同 forcing term 策略的实际表现大致如下forcing term 策略外层迭代数内层总扫描数到达 1e-10 残差固定 η 1e-2962可以但偏慢固定 η 1e-86190可以内层过度求解η min(0.1, sqrt(r_k))674推荐无策略评估纯值迭代2000—不可行表里的数字是量级示意具体会随问题和秩上限浮动但相对关系是稳定的自适应 forcing term 同时拿到了二次收敛的迭代数和接近最优的内层开销。5. 避坑排查TT 版策略迭代最容易翻车的五个地方5.1 前三次迭代飞快第四次开始原地踏步现象很典型残差序列是 1e-1、1e-3、1e-6然后卡在 8e-7 一动不动扫描次数加到 500 也没用。原因几乎一定是截断容差写死了。如果 opts.tol 固定成 1e-8而当前残差已经在 1e-6 量级那么每次取整引入的误差就与牛顿方向的量级相当牛顿步走的方向被噪声主导——本质上是拿一个精度只有 1e-8 的解去算一个需要 1e-12 精度的修正量。解决办法就是 4.3 那套tol 必须跟着残差走tol_k ≈ 0.1×r_k并且下限设在 1e-12 附近。改完之后通常能看到残差重新出现超线性下降最后一次迭代从 1e-6 直接跳到 1e-11。5.2 单站点扫描 200 轮残差只降一个数量级这是 4.2 描述过的秩固定陷阱最直接的表现。判断方法很简单打印每一轮扫描后各核的实际秩如果从头到尾一个数都没变那就是被秩上限锁死了。要区分两种情况。真解的秩确实低于上限那就是局部极小问题加残差增强AMEn能出来真解的秩高于上限那就是上限设置错误得放宽 maxRank 再看。后一种情况有个实用判据把 maxRank 从 r 提到 1.5r如果残差明显继续下降说明瓶颈是秩如果几乎不动说明瓶颈在别处多半是 tol 或者算子本身构造错了。5.3 策略改进一步秩从 12 涨到 60现象是内存和时间在第 2、3 轮突然爆掉profile 一跑发现时间全花在取整上。原因是策略改进那一步做了逐点 min。min 是非线性操作没有任何保秩的性质。m 个候选值张量每个都是“代价 γ·插值后的 V”逐点比较结果的秩上界可以随 m 和插值带宽一起涨。更麻烦的是策略在网格上切换的位置会形成“界面”界面越不规则秩越高。可行的做法有两个方向。一是绝不在全网格上显式计算 min用 TT-cross 式的自适应采样只在一批选定的网格点上算 min然后用这些样本拟合一个新的 TT 张量拟合时直接带上 maxRank 上限。二是把秩控制前置候选值先各自取整到较低精度再比较虽然引入了额外误差但秩可控。前一种更准后一种更快实际代码里常见的是两者结合——策略改进轮用采样拟合策略评估轮用严格取整。5.4 折扣因子提到 0.99耗时涨 20 倍γ 从 0.9 提到 0.99物理上只是“看得更远”计算上却是双重打击。第一重是算子条件数 1/(1−γ) 从 10 变成 100内层扫描次数线性上升第二重更隐蔽——值函数变得更“平坦”最优策略的切换界面变宽TT 秩需求上升。先试热启动几乎零成本能砍掉一半扫描。再不行就做 γ 延拓先用 γ 0.9 解出一个 V把它当 γ 0.95 的初值逐级往上推。每一步都是良态问题最后一步虽然难但初值已经非常接近真解。这套做法在折扣因子接近 1 的问题上几乎是标配代码里加十几行就能落地。5.5 网格加密后精度不升反降把每维网格数从 32 加到 64理论离散误差应该降 4 倍结果残差反而变大策略在相邻点之间抖。原因通常是两个误差源失配。网格加密让插值误差 O(h²) 降下来但 TT 截断容差没跟着调截断误差成了主导更糟的是策略抖振——网格变细之后最优策略的切换边界被分辨出来边界附近的网格点对 V 的微小扰动极度敏感一扰动就换动作换了动作又扰动 V。这不是数值 bug是离散化本身的病态。处理办法是让 tol 与 h² 匹配h 减半tol 降 4 倍并在策略改进里加一个“迟滞”只有当候选动作的代价比当前动作低出某个阈值 δ 时才切换δ 取 h² 量级。代价是最优性上损失一点点收益是迭代终于能收敛。6. 设计取舍什么时候该用牛顿TT什么时候不该把几条路线摊开对比取舍就清楚了。方案可处理维度每轮成本收敛速度精度保证主要风险稠密网格 策略迭代d ≤ 4O(N) 内存二次到浮点精度内存墙值迭代 TTd ≤ 10一次 TT 矩阵向量乘线性速率 γ到截断容差γ 接近 1 时不可用牛顿策略迭代 TTd ≤ 10一次 TT 线性求解二次到截断容差秩爆炸、容差失配神经网络拟合值函数任意 d一次反向传播不确定无无证书、调参玄学这张表里最容易被误读的是最后一行。神经网络方法能上高维但它没有残差证书——你无法验证解对不对只能看训练损失。而 TT 方法的价值恰恰在于残差可以在 TT 格式里严格算出来tt_residual收敛有判据精度有边界。对于控制律要上真实系统的场景这个差别比速度重要得多。那 TT 方法的边界在哪我认为就一条值函数的 TT 秩是否随维度增长而失控。这条边界不能靠理论预估只能实测。可靠的验证方法是做秩收敛检查——用 maxRank 20 和 maxRank 40 各跑一遍比较两个解在随机采样点上的最大偏差。如果偏差远小于容差说明秩已经饱和结果可信如果偏差和容差同量级那这个问题的解根本不在低秩流形上再怎么调参也没用。第二个检查是直接验 HJB 残差在随机采样的状态点上算 V(x) − min_u{g γV(f(x,u))}这个量绕过了所有 TT 内部表示是对物理问题的直接体检。还有一个不太被提及的取舍这套方法把复杂度从“维度”转移到了“交互阶数”。稠密网格的代价是 n^dTT 的代价大致是 d·n·r³r 由动力学的耦合结构决定。换句话说你不是免费拿到了高维而是把“维度灾难”换成了“耦合灾难”。对链式弱耦合系统这是巨大的胜利对全耦合的强非线性系统则毫无优势。选型的时候先看动力学怎么写再看维度有多高——顺序反了会白忙一场。从那以后我每次拿到一个新的控制问题都强制先走一遍这个顺序先把动力学按耦合块画出来估一下交互阶数再决定用不用 TT跑通之后一定做那两个验证——双秩对比和采样点上的 HJB 残差两个都过了才敢往下走。这套资源里的 MATLAB 实现把策略评估、双站点扫描、非精确牛顿的容差调度都写全了对照着源码把这几处接缝看一遍比读十篇综述都管用。希望帮到你。本文还有配套的精品资源点击获取