简介一份面向机器学习与强化学习初学者的MDP策略迭代MATLAB实践资源以马尔科夫决策过程为背景重点演示策略评估与策略改进两个核心步骤的实现过程。压缩包共93个文件以m源码、html文档和dat数据为主另有txt说明、ps/eps图像以及readme等辅助材料MDPtoolbox工具包提供了策略迭代、值迭代、有限期决策等函数及配套文档页面便于对照学习。整个包约2.32MB结构清晰已有312人学习下载。通过运行策略迭代相关脚本读者可以理解状态值函数计算、策略更新与收敛判断的完整流程并能基于自定义转移矩阵和奖励函数进一步扩展适合作为强化学习入门与动态规划算法实验的参考。1. 从 MDP.zip 说到策略迭代为什么手写一遍 MATLAB 实现比调库更靠谱看到 “MDP.zip” 这个压缩包名字第一反应就是里面装了一套基于马尔可夫决策过程MDP的策略迭代Policy IterationMATLAB 实现。策略迭代是强化学习里最经典的求解方法和值迭代相比它在状态空间不大的离散问题上往往能靠更少的决策轮次收敛到最优策略。很多人觉得 MATLAB 有现成工具箱、有 Python 库没必要手写问题在于 MDP 建模、转移矩阵排布和收敛判定这些细节都会在调库时被藏起来。这里就把这些边界讲清楚从一个能跑通的 5×5 网格世界切进去看完你能自己写出策略迭代也能拿同一套思路去排查别人的代码。适合刚接触马尔可夫决策过程、需要在 MATLAB 里做对比实验的从业者也适合想把老代码从值迭代迁到策略迭代的人。2. 策略迭代的核心逻辑为什么它总能收敛到最优策略2.1 MDP 五元组如何映射成 MATLAB 矩阵策略迭代解决的是这样一个问题在马尔可夫决策过程中找一个让累计折扣回报最大的策略。一个 MDP 有五个部分状态空间 S、动作空间 A、转移概率 P、奖励函数 R、折扣因子 gamma。在 MATLAB 里最直接的映射方式是状态用 1 到 nState 的整数编号动作用 1 到 nAction 的整数编号转移概率用三维数组 P 来存P(s, s, a) 表示在状态 s 执行动作 a 之后转移到 s 的概率。奖励用 nState 行 nAction 列的矩阵 R 来存R(s, a) 表示状态 s 执行动作 a 后拿到的即时收益。状态编号从 1 开始不是个人习惯而是 MATLAB 的硬性约束。任何从伪代码迁移过来的状态 0 都会在索引上翻车这点后面第 5 章还会细说。下面这段代码构建了一个 5×5 网格世界的转移概率和奖励我认为它是很多 MDP 实验的公共底座。% 5x5 grid, nState25, nAction4 % 动作编号1上, 2下, 3左, 4右 nState 25; nAction 4; P zeros(nState, nState, nAction); % P(s, s, a) R zeros(nState, nAction); % R(s, a) for s 1:nState [row, col] ind2sub([5 5], s); for a 1:nAction nextRow row; nextCol col; if a 1, nextRow row - 1; end % 上 if a 2, nextRow row 1; end % 下 if a 3, nextCol col - 1; end % 左 if a 4, nextCol col 1; end % 右 if nextRow 1 || nextRow 5 || nextCol 1 || nextCol 5 nextRow row; nextCol col; % 撞墙就留在原地 end sNext sub2ind([5 5], nextRow, nextCol); P(sNext, s, a) 1; % 确定性转移概率集中在 sNext R(s, a) -1; % 每走一步代价 -1 end end % 把最后一个格子设为终止状态 P(:, nState, :) zeros(nState, 1, nAction); for a 1:nAction P(nState, nState, a) 1; % 终止状态自环 R(nState, a) 0; % 到达终止不再扣分 end这里把 P 设计成 P(s, s, a) 是为了后面策略评估时用squeeze(P(:, s, a))取到“从某状态出发的转移向量”第二维正好是当前状态不容易出现转置错位。终止状态被设成在所有动作下都自环、奖励为 0这样它不会影响其它状态的价值回灌逻辑也不会让状态本身积累不存在的收益。2.2 策略评估固定策略下价值函数的线性方程形式策略迭代第一步是评估当前策略的价值。给定一个确定性策略 pi(s)每个状态下只选一个动作价值函数满足贝尔曼方程V_pi R_pi gamma * P_pi * V_pi。这是一个线性方程组MATLAB 里可以直接用左除来解也就是写成(I - gamma * P_pi) \ r_pi。左除比 inv 更稳也更符合实际工程的习惯。function V policy_eval(P, R, gamma, pi) % pi: nState x 1每个元素取值为 1 到 nAction nState size(P, 2); P_pi zeros(nState, nState); r_pi zeros(nState, 1); for s 1:nState a pi(s); P_pi(s, :) squeeze(P(:, s, a)); % 从 s 执行 a 后转移到各状态的概率 r_pi(s) R(s, a); end V (eye(nState) - gamma * P_pi) \ r_pi; end这个函数里最容易错的一行是squeeze(P(:, s, a))。P 是三维数组取出的P(:, s, a)在 MATLAB 中会被保留成 nState×1×1 的数组squeeze 去掉两个单例维度后变成列向量转置后才是行向量这样才能直接赋值给 P_pi 的某一行。如果你丢掉转置结果会变成行写入失败或者形状不匹配。参数说明gamma 在 0.9 附近时矩阵I - gamma * P_pi是对角占优的左除很稳定gamma 一旦到 0.99矩阵条件数会变大解出来的价值表面可能仍然正确但离真实价值会有一个被放大的误差项。对于 nState 只有 25 的小网格直接左除是最省事的不需要迭代。2.3 策略改进从价值函数到更优的贪婪策略得到当前策略的价值函数 V 后策略改进要做的是对每个状态比较所有动作的 Q 值然后选最大。Q(s, a) R(s, a) gamma * sum_{s} P(s, s, a) * V(s)。写成 MATLAB 就是下面的函数。function pi_new policy_improve(P, R, gamma, V) nState size(P, 2); nAction size(P, 3); Q zeros(nState, nAction); for a 1:nAction Q(:, a) R(:, a) gamma * squeeze(P(:, :, a)) * V; end [~, pi_new] max(Q, [], 2); end这里squeeze(P(:, :, a))得到一个 nState×nState 的矩阵它的第一维是 s第二维是 s。矩阵乘 V 相当于对每个 s 做了一个加权求和结果正好是gamma * E[V(s)]再累加即时奖励就是 Q 值。max(Q, [], 2)是沿着第二维求最大值返回每个状态的最优动作索引。策略改进有一个很隐蔽的性质如果新策略和旧策略不同那么它带来的价值函数一定严格优于旧策略如果完全相同说明当前策略已经满足贝尔曼最优条件。因此策略迭代不需要像值迭代那样设一个价值阈值来判定停止只要盯着策略是否变化就行。这个收敛判定是整个算法里最省心的地方。2.4 为什么有限 MDP 上策略迭代不会无限循环确定性策略的总数是 nAction 的 nState 次方是一个有限集合。策略迭代从某个确定性策略出发每一轮要么严格改进策略要么不动。因为改进操作保证价值单调上升所以算法不可能陷入不同策略之间来回切换的循环一旦出现重复策略只可能是已经达到最优。这种“停止条件天然存在”的优势在工程里很关键你不需要猜阈值。和值迭代对比时有一个常见误区值迭代必须用相邻两轮价值差的绝对值判断收敛阈值设大了会早停得到次优策略设小了会多跑很多轮。策略迭代则直接把收敛条件放在策略本身上。我整理了一张对比表方便你决定在什么场景选哪种算法。对比项策略迭代值迭代核心更新评估当前策略 - 贪心改进每轮直接应用贝尔曼最优算子停止条件策略不再变化相邻价值差小于阈值单轮计算成本高需要解线性方程组低只需一次 max 和矩阵乘对于 gamma 接近 1可能出现病态矩阵但策略通常稳定收敛速度变慢阈值选择变难典型迭代轮数很少网格世界常见 3 到 5 轮可能几十轮以上小规模状态空间我基本选策略迭代因为轮数少、停止条件干净当状态上千后策略迭代每轮线性方程组成本按 O(nState^3) 增长这时候反而值得考虑值迭代或者后面第 6 章提到的修正策略迭代。这些边界不是书本上的玄学是实际跑出来需要权衡的成本。3. 在 MATLAB 里写一个最小可复现的策略迭代工程五元组、函数与可视化3.1 状态、动作、奖励的全量初始化有了原理接下来把它编译成能直接跑的脚本。我的习惯是先把地图、状态编号、动作编号固定下来再写算法主循环这样任何一步出错都能用whos和画出转移矩阵来定位。下面就是一个完整脚本的主干它使用上一章的 5×5 网格世界初始策略选全向上gamma 先取 0.95。clear; clc; gamma 0.95; nState 25; nAction 4; % 构建 P 和 R逻辑与第二章相同 P zeros(nState, nState, nAction); R zeros(nState, nAction); for s 1:nState [r, c] ind2sub([5 5], s); for a 1:nAction nr r; nc c; if a 1, nr r - 1; end if a 2, nr r 1; end if a 3, nc c - 1; end if a 4, nc c 1; end if nr 1 || nr 5, nr r; end if nc 1 || nc 5, nc c; end sNext sub2ind([5 5], nr, nc); P(sNext, s, a) 1; R(s, a) -1; end end for a 1:nAction P(25, 25, a) 1; % 终止状态自环 R(25, a) 0; end pi ones(nState, 1); % 初始策略全部向上 for iter 1:50 V policy_eval(P, R, gamma, pi); pi_new policy_improve(P, R, gamma, V); if isequal(pi_new, pi) fprintf(converged at iter%d, optimal strategy found\n, iter); pi pi_new; break; end pi pi_new; end这段代码把主循环压缩到最简。iter 1:50是一个保险网防止因为终止状态或者概率构造问题导致策略在两个动作间抖动时无限循环。实际调试时如果发现 50 次还没停我会先去查转移矩阵是否归一化而不是盲目调大 iter 上限。3.2 策略评估的两种写法直接解 vs 迭代更新策略评估可以有两种实现各有各的适用面。直接解法适合状态数小于几千的离线分析迭代法适合状态多到连(eye - gamma*P)\r都解不动的情况。直接解法在 2.2 已经给过这里补充迭代法的实现。function V eval_policy_iter(P, R, gamma, pi, tol, maxIter) nState size(P, 2); V zeros(nState, 1); for k 1:maxIter Vnew zeros(nState, 1); for s 1:nState a pi(s); % 注意这里 squeeze(P(:,s,a)) 是行向量点乘 V 是加权平均 Vnew(s) R(s, a) gamma * squeeze(P(:, s, a)) * V; end if norm(Vnew - V, inf) tol V Vnew; break; end V Vnew; end end迭代法里 tol 和 maxIter 是真正的止血点。我一般把 tol 设为 1e-6maxIter 设为 10000gamma 到 0.99 时迭代次数会明显上升但不会因为没设上限而卡死。如果你把 tol 设成 1e-12在 gamma 接近 1 时往往要跑到接近上限才满足这是没必要的因为策略改进真正依赖的是 V 的相对排序而不是绝对精确到最后几位小数。3.3 主循环与策略可视化用图验证收敛结果调试策略迭代最直观的办法是把最终的价值函数和策略方向画出来。下面的代码用 imagesc 显示价值再在格子中心用文字标注动作编号。figure; imagesc(reshape(V, 5, 5)); colorbar; title(State Value after Policy Iteration); axis xy; % 让行坐标从上到下符合网格直觉 for s 1:nState [r, c] ind2sub([5 5], s); % text 的坐标参数是 (x,y)所以第一参数必须是列 c第二参数是行 r text(c, r, num2str(pi(s)), HorizontalAlignment, center, ... Color, w, FontWeight, bold, FontSize, 14); end一句话经验text(c, r, ...)里 c 和 r 的位置非常容易写反我在这上面翻过车。imagesc 默认 y 轴方向是自下而上加上axis xy之后第 1 行才会显示在图像顶部和ind2sub得到的行列语义保持一致。这样一张图能同时验证最优动作是否指向目标、价值是否从目标状态向外递减比盯着 25 个数字排查效率高得多。4. 调参实战gamma、评估精度与停止条件怎么选才算稳4.1 折扣因子 gamma既要任务语义也要矩阵健康gamma 决定智能体有多看重远期收益。在网格世界里gamma 越接近 1智能体越愿意绕远路去拿远处更大的奖励gamma 偏小时它会更倾向于立刻结束。但 gamma 直接影响策略评估中I - gamma * P_pi的条件数gamma 越接近 1矩阵越接近奇异左除可能给出数值很大的中间结果不过最终对策略排序的影响不一定致命。我的参数选择策略是先让逻辑跑通再调任务语义。第一次跑通时用 gamma0.9确认策略合理后再按任务需求调到 0.95 或 0.99。如果一上来就 0.99评估矩阵的条件数大了一旦转移概率有轻微归一化误差价值函数可能直接爆出 NaN到那时候很难判断是算法错还是参数错。可以随时在命令行看一眼条件数% P_pi 需要按照当前策略构造这里只展示命令 cond(eye(nState) - gamma * P_pi)条件数量级超过 1e8 就要怀疑是不是 gamma 太接近 1 或者转移矩阵有误。gamma1 在无终止的 MDP 上会让价值函数不收敛除非改用平均奖励准则否则不要把它设成 1。gamma 范围适用任务风险提示0.8 - 0.95短视任务、奖励密集矩阵健康收敛快0.96 - 0.99需要远期规划病态风险上升迭代法收敛变慢1.0仅在有保证终止的任务无终止时价值发散不建议4.2 策略评估的精度不是越严越好如果你采用迭代法评估tol 的选择会影响策略改进的质量。设得太小比如 1e-12加上 gamma0.98评估这一轮可能要几千次迭代设得太大V 里有较明显的误差下一轮改进可能选出实际次优的动作。我的启发式做法是让 tol 比奖励量级低两到三个数量级奖励是 1tol 用 1e-6奖励是 0.1tol 用 1e-8。这样既不会浪费算力也不会让误差压过真实差异。还有一种“故意欠评估”的做法对应第 6 章要展开的修正策略迭代。现在需要记住的是策略迭代的每一轮策略评估不必完全收敛但误差必须被控制在同一个数量级。实际调试时我常用这行代码检查当前策略是否已经稳定V_old eval_policy_iter(P, R, gamma, pi, 1e-8, 5000); [~, pi_new] max(compute_Q(P, R, gamma, V_old), [], 2); if isequal(pi_new, pi) fprintf(strategy stable at iter %d\n, iter); end这里的 compute_Q 可以理解为 2.3 中 policy_improve 里的 Q 计算部分。如果策略稳定但价值差仍然很大说明只是评估没有收敛完全收敛条件本身没有错。4.3 停止条件看策略还是看价值差这是策略迭代代码里最常见的“黑匣子”问题。有些人在主循环里模仿值迭代用norm(V_new - V, inf) tol作为停止条件。策略迭代这样做会让程序在最优策略已经出现后继续空转因为价值函数的绝对变化还需要更多轮次才能降到阈值以下更危险的是如果价值差阈值设得比策略改进的边际价值差更大程序会在策略还没有达到最优时提前退出得到一份看着合理其实次优的策略。正确做法是主循环只比较策略是否变化价值评估只作为策略改进的输入。我见过的最稳妥判定是同时使用两个条件策略完全一致且价值函数差值已经很小。不过在纯策略迭代里策略一致已经蕴含了最优解价值差条件只是额外保险。4.4 策略更新顺序同步改进和异步改进不能乱用标准策略迭代是全状态同步改进先用当前 V 对所有状态计算 Q然后统一生成 pi_new下一步再重新评估。如果你在同一个 for 循环里更新一个状态的动作后立刻用它修改后的 V 再更新下一个状态就变成了异步策略改进。异步在某些特定 MDP 上可能不收敛尤其当状态之间有强耦合时价值函数会被反复抹掉。这不是理论上的吹毛求疵我在一个线性链结构的模拟项目X上遇到过策略反复跳变的情况正是因为把同步写成了异步。实现上如何保证同步不要写成for s1:nState, QR(s,:)...; [~,a]max(Q); pi(s)a; end这样逐个更新而是先循环动作把整个 Q 矩阵算完再一次max(Q, [], 2)。2.3 里的 policy_improve 就是标准同步写法建议直接作为模板。5. 策略迭代 MATLAB 实现的五个典型翻车现场与排查清单5.1 转移矩阵概率行不归一价值函数悄悄膨胀现象价值函数里出现超过合理范围几个数量级的数值比如网格世界最优值应该在目标附近收敛到某个有限值结果出现了一百多。策略在角落里反复横跳评估结果对初始策略极其敏感。原因构建转移矩阵时只设置了有效转移的概率却把“撞墙停在原地”的转移漏掉。转移概率行加起来小于或大于 1贝尔曼算子不再是一个收缩映射价值函数就会失真。这个问题在确定性转移中尤其隐蔽因为你不会看到任何报错。解决每次构建完 P 后对所有 (s, a) 组合做归一化检查。终止状态允许概率和为 1其它状态必须有sum(P(:, s, a)) 1。检查代码放在构建函数末尾几乎是零成本。for s 1:nState for a 1:nAction if s terminal % terminal 需要提前定义为终止状态的编号 continue end if abs(sum(P(:, s, a)) - 1) 1e-12 warning(un-normalized transition at state %d action %d, s, a); end end end5.2 MATLAB 索引从 1 开始状态 0 混入导致矩阵越界现象主循环有时候能跑有时候报Index in position 1 exceeds array bounds。更诡异的是报错经常只在某些动作下出现换一个初始策略就不报。原因很多伪代码用 0-based 状态编号重写成 MATLAB 时把 0 直接当索引用了。MATLAB 的线性状态编号不支持 0强行使用会让数组维度错位。还有一个同源问题算格子坐标时直接用row1或col-1在边界处得到 0然后sub2ind传进了 0 或负坐标。解决统一从 1 到 nState 编号所有坐标计算后用sub2ind转线性索引并且先夹紧边界。动作编号也按 1 到 4 对应上下左右。我通常在脚本开头贴一行注释状态、动作都是 1-based保证自己在连续写两三小时后不会犯这个低端错误。5.3 max 返回第一个最优动作策略在等价最优间反复震荡现象策略迭代不收敛但价值函数已经非常接近。每一轮 pi_new 和 pi 只有一个或几个状态的动作不同价值差也小可主循环始终不停。原因max 函数在并列最大值时返回第一个索引。由于浮点噪声两个等价动作的 Q 值可能一个比另一个多 1e-15导致改进结果在两次评估之间抖动。如果收敛判定只比较策略你就会看到策略来回变。解决不要用策略完全一致作为唯一停止条件。我常用的做法是先检查策略差异如果差异为零直接停如果策略差异不是零但价值差小于一个小阈值再额外看差异状态上的 Q 值差距若差距小于浮点噪声量级则强制停止。另一种工程做法是对 Q 值做微小扰动但这不推荐因为它会掩盖真实策略差异。if isequal(pi_new, pi) || norm(V_new - V_old, inf) 1e-10 break; end上面的写法等于在策略比较之外加一道保险能覆盖大多数浮点抖动场景。5.4 终止状态参与价值回灌目标附近价值被污染现象目标状态本身的价值不是 0而是随着迭代变化目标附近状态的价值出现凹凸不平甚至有几个状态的价值高于目标状态。原因策略评估时把终止状态当作普通状态参与更新而且它的转移又被设成一个任意的动作导致终止状态通过自环回灌价值。如果R(term, a)不为 0每次更新都在累积收益如果 P 里 term 没有自环而指向别处还会把价值广播到错误状态。解决在构建 MDP 时就把终止状态“焊死”。我这里说的焊死是指终止状态下所有动作转移到自身奖励为 0并且迭代评估时强制该状态价值保持为 0。直接解法里则把 V(term) 在方程中固定为 0不要让它出现在未知数向量里。这样周围状态可以通过转移进入终止状态承担 0 价值但终止状态本身不会有额外回灌。5.5 迭代评估不收敛gamma 太大、阈值太小、maxIter 不够现象策略迭代跑了很久评估函数的 maxIter 耗尽V 返回一个中间结果后续策略改进变得很怪最终出现 NaN。原因gamma0.99 时贝尔曼算子的收缩因子接近 1每次迭代价值变化很小要达到 1e-12 需要几千甚至上万次更新。如果 maxIter 设成 1000评估还没收敛就被掐断改进用的 V 并不是当前策略的真实价值策略就会偏离。解决评估迭代的 tol 放宽到 1e-6maxIter 至少 10000或者干脆在状态数不大时改用左解。如果确实需要在 gamma 接近 1 时用迭代法可以考虑下面第 6 章的修正策略迭代用有限步评估替代收敛评估但要明白它牺牲的是评估精度换来的是单轮速度。6. 从玩具代码到能用的策略迭代一个能减少等待时间的固定步数评估技巧策略迭代在状态数变大后最痛的一步就是评估。每一轮都解一个 nState 阶线性系统nState 到几千后左除成本按立方增长。一个务实技巧是用“修正策略迭代”评估阶段不收敛只固定跑 K 步价值迭代然后马上做策略改进。K 一般取 5 到 10我习惯从 8 开始。很多任务在这个设置下和完整评估得到的最优策略几乎一样但每一轮评估时间能少一个数量级。K 8; for iter 1:30 V eval_policy_iter(P, R, gamma, pi, 1e-12, K); % maxIterK等于固定步数 pi_new policy_improve(P, R, gamma, V); if isequal(pi_new, pi) break; end pi pi_new; end注意这里的 eval_policy_iter 最后两个参数被我用 1e-12 和 K 占住了所以它不会因为 tol 太严而多跑也不会在 maxIter 处停下来而是严格执行 K 步。这样评估结果有误差但策略改进对误差的敏感度通常不高尤其当奖励稀疏时早期策略的区别主要来自“有没有通向目标的路径”而不是“路径上的微小差异”。我一般会在固定步数评估收敛后再用完整左解评估一次最终策略拿到一个精确的 V 来做策略可视化或和值迭代结果对拍。对拍的方法是把两种算法求出的价值函数做差norm(V_pi - V_vi, inf)小于 1e-6 就说明两个实现互相印证。这比单看策略图像可靠因为策略在离散空间里可能看起来相同但价值差异能暴露建模错误。我的习惯是任何 MDP 求解器都先在 3×3 手工可算的地图上对拍手算出期望的最优路径再跑随机地图。这个习惯帮我挡住了至少三次转移矩阵方向写反的问题。把固定步数评估、双重停止条件和 3×3 对拍合在一起策略迭代就不再是一个需要反复猜测参数的黑匣子。希望帮到你。本文还有配套的精品资源点击获取