简介这是一套面向机器学习与强化学习初学者的马尔科夫决策过程MDP实践资料聚焦策略迭代算法的 MATLAB 实现与原理讲解可帮助读者在离散状态和动作空间中快速上手 MDP 建模、策略评估与策略改进流程。压缩包共 93 个文件大小约 2.32MB包含 33 个 MATLAB 函数脚本、20 个 HTML 算法说明文档、18 个 dat 数据文件、txt 说明与 README以及少量图例、EPS 图和 LaTeX 排版文件代码与文档配套清晰。以策略迭代、有限期、价值迭代等函数为核心并附有 race 等示例场景便于对照运行和验证结果文档部分则系统说明各函数的输入输出与算法逻辑。目前已有 312 人学习浏览。通过阅读源码和文档学习者既能理解策略评估与策略改进的交替收敛机制也能直接修改示例模型进行实验适合作为 MDP 入门、课程作业或强化学习基础研究的参考工具。1. 策略迭代不是论文专属MDPtoolbox 把策略迭代变成 MATLAB 的一行调用初次接触马尔可夫决策过程MDP的人对“策略迭代”四个字通常又熟悉又陌生熟悉是因为它是机器学习、强化学习教材里的常客陌生是因为公式推到一半就不知道怎么写代码。MDPtoolbox 提供了一个直接的答案——把策略迭代相关的核心算法封装成 mdp_policy_iteration、mdp_eval_policy、mdp_finite_horizon 一整套 MATLAB 函数输入状态转移矩阵、奖励矩阵和折扣因子就可以拿回最优策略和值函数。它的价值不在于替代理论学习而在于把“模型定义 → 代数求解 → 策略输出”这条链路打通。正在做强化学习课程设计或者在调度、库存、设备维护这类离散决策场景里建模的人用这套工具能少走大量弯路。2. 马尔可夫决策过程与策略迭代先把“评估—改进”循环讲透2.1 MDP 的建模要素状态、动作、转移矩阵和奖励矩阵怎么对应到代码MDP 在教材里往往写成五元组状态集合 S、动作集合 A、状态转移概率 P、奖励函数 R、折扣因子 γ。但真正写 MATLAB 代码时只需要准备两个矩阵 P 和 R状态和动作的数量都隐含在矩阵维度里。P 的形状是 S×S×A读法是“当前状态 s 执行动作 a 后转移到下一状态 s 的概率”P(s, s, a) 就是这个概率值。R 的形状通常是 S×A代表状态 s 下执行动作 a 的即时奖励也有部分实现允许 S×S×A 的奖励即把“转移到哪个状态”也考虑进奖励但 MDPtoolbox 的常规用法是前者。这里有一个新手最容易搞反的地方P 的维度顺序。MATLAB 的索引是按列优先存储的P(s, s, a) 里前两个下标是状态第三个下标是动作。我在第一次定义 P 时习惯性地把它写成了 P(s, a, s)结果算出来的策略全部错位后来发现是维度顺序的问题。因此定义矩阵之前先给变量写下注释% P: S×S×A, P(s, s, a) Pr(next_states | states, actiona)。这个注释值十行代码。折扣因子 γ 在 MDPtoolbox 里是独立参数它的含义是未来奖励折算到当前时刻的比率。γ 越接近 1算法越看重长期收益越接近 0模型越短视只看眼前一两步的奖励。策略迭代对 γ 的要求是 0 ≤ γ 1如果 γ 恰好等于 1策略评估对应的线性方程组可能变成奇异矩阵值函数解不出来。这一点在 5.2 节还会展开。实际使用中我一般把 γ 当作一个可调超参数而不是固定值库存补货问题用 0.95设备维护问题用 0.9单步决策问题直接用 0因为不需要跨周期累积。2.2 策略评估贝尔曼方程在 mdp_eval_policy 里做了什么策略迭代分两步走第一步是策略评估Policy Evaluation计算当前策略 π 下每个状态的值函数 Vπ(s)。定义式是贝尔曼方程Vπ(s) Σ_a π(a|s) Σ_s P(s|s,a) [ R(s,a) γ Vπ(s) ]这个公式初看有点绕拆开看就三件事当前状态 s 下按策略 π 选动作 a执行 a 之后按概率 P(s|s,a) 跳到下一个状态跳到 s 之后拿到即时奖励 R(s,a)再加上折算后的未来回报 γVπ(s)。把所有状态都列出来它就是一组线性方程。策略固定之后动作选择的不确定性被消除可以把转移矩阵聚合成一个 S×S 的矩阵 Pπ把奖励聚合成向量 Rπ于是方程变成(I - γ Pπ) Vπ Rπ这组方程的解就是当前策略的值函数。工具箱里的 mdp_policy_iteration 在评估阶段做的事本质上就是解这个线性方程组。标准策略迭代每次评估都从零开始解一遍而修改版本 mdp_policy_iteration_modified 则允许你直接指定初始策略这一点在第 4 章再细说。如果自己写一个最简单的策略评估MATLAB 里几十行就能跑通function V eval_policy_naive(P, R, discount, policy, tol) % 极简策略评估迭代法解贝尔曼方程 % 输入 P: S×S×AR: S×Apolicy: S×1tol: 收敛容差 S size(P, 1); V zeros(S, 1); diff inf; while diff tol V_new zeros(S, 1); for s 1:S a policy(s); V_new(s) R(s, a) discount * (squeeze(P(s, :, a)) * V); end diff max(abs(V_new - V)); % 前后两轮值函数的最大变化量 V V_new; end end这个函数对应 P(s,:,a) 表示状态 s 执行动作 a 后的一行转移概率squeeze 把它从 S×1×1 压成 S×1 向量乘以 V 就是对下一状态值函数求期望。实际 MDPtoolbox 内部不是这么循环的它直接用稀疏矩阵解线性方程速度会快一个量级但语义完全一致。容差 tol 取 1e-8 就够如果取到 1e-12 反而可能在边界上费额外的迭代次数收益很小。2.3 策略改进贪心更新与收敛判据策略评估算出了 Vπ第二步策略改进Policy Improvement就做一件事对每个状态 s贪心地选一个让动作值函数最大的动作作为新策略 π 在这个状态下的输出。动作值函数的定义是Qπ(s,a) R(s,a) γ Σ_s P(s|s,a) Vπ(s)对所有动作计算完 Q 值之后取最大者。如果新策略 π 和旧策略 π 在每个状态上都一样说明策略已经收敛到最优迭代结束只要有一个状态的动作变了就用 π 替换 π回到策略评估继续循环。这个“评估—改进—再评估—再改进”的循环是策略迭代的核心结构。理论上有限状态和有限动作空间下策略迭代经过有限轮就一定收敛不会像值迭代那样需要预设收敛容差。这里有个容易忽略的细节策略改进阶段取最大值时如果没有特别处理并列情况工具箱通常返回序号最小的那个动作。也就是说如果状态 s 下动作 1 和动作 2 的 Q 值完全相同返回哪个完全取决于遍历顺序。这在数值上没有任何问题但如果你手工推导的结果和工具箱输出差一个动作先检查是不是并列 Q 值造成的而不是怀疑算法写错了。这个细节在文档里一般不会写但实际对上答案时非常关键。3. 跑通第一个策略迭代解压配置、定义矩阵、调用求解3.1 解压后先看 README目录结构决定后面少踩哪些坑这份压缩包解压出来是一个完整的工具包不是单个脚本里面有 mdp_LP.m、mdp_policy_iteration.m、mdp_value_iteration.m、mdp_finite_horizon.m、mdp_check.m 这些核心函数还有一批 html 文档和 examples 示例目录。拿到手第一件事不是急着运行而是把它放到一个稳定的目录比如D:/projects/mdptoolbox然后把它加进 MATLAB 搜索路径。addpath(D:/projects/mdptoolbox); savepath;addpath 的作用是让当前 MATLAB 会话能够找到这个目录下的 m 文件savepath 把路径配置保存到用户路径文件这样下次启动 MATLAB 时就不需要重新添加。如果只在某一个脚本里临时使用也可以不加 savepath每次运行前在脚本开头写 addpath 就行。验证路径是否生效用 which 命令which mdp_policy_iteration这条命令会返回 mdp_policy_iteration.m 的完整路径。如果返回空或者提示未找到说明路径没有配置成功检查目录名是否写错或者改成相对路径再试一遍。3.2 定义 P 和 R4 状态 2 动作的设备维护小模型直接跑文档里的例子之前最好先自己定义一个小模型跑一遍完整链路。这里用 4 个状态、2 个动作的设备维护问题状态 0 到 3 表示设备磨损程度递增动作 1 是“维护”动作 2 是“不维护”。维护会花钱但不维护会让状态越来越差。S 4; % 状态数 A 2; % 动作数 P zeros(S, S, A); % P 维度顺序是 当前状态 × 下一状态 × 动作 % 动作 1维护设备大概率回到状态 1小概率保持原状态或恶化 P(1, 1, 1) 0.9; P(1, 2, 1) 0.1; P(2, 1, 1) 0.8; P(2, 3, 1) 0.2; P(3, 1, 1) 0.7; P(3, 4, 1) 0.3; P(4, 1, 1) 0.6; P(4, 4, 1) 0.4; % 动作 2不维护状态向更差的方向游走 P(1, 2, 2) 1.0; P(2, 3, 2) 1.0; P(3, 4, 2) 1.0; P(4, 4, 2) 1.0; % 奖励矩阵S×A行对应状态列对应动作 R zeros(S, A); R(:, 1) [-1; -2; -3; -4]; % 维护成本状态越差维护费用越高 R(:, 2) [ 5; 3; 0; -5]; % 不维护短期收益高但状态差时损失大定义 P 的核心要求是对每个动作 aP(:,:,a) 的每一行之和必须等于 1也就是从任意状态出发转移到所有状态的完整概率分布总和为 1。上面代码里动作 2 的第 1 行只写了 P(1,2,2)1.0这一行只剩一个非零项行和正好等于 1没问题。奖励矩阵的符号也值得注意正奖励是吸引算法去选择这个动作负奖励是促使算法回避。状态 1 不维护给 5状态 4 不维护给 -5算法自然会倾向于在早期不维护、在晚期选择维护。3.3 策略迭代与值迭代对比参数含义和结果解读矩阵就绪之后调用策略迭代只需要四行代码discount 0.95; max_iter 1000; [policy, V, cpu_time] mdp_policy_iteration(P, R, discount, max_iter); disp(policy); disp(V);参数说明discount 是折扣因子按 2.1 节的原则设置max_iter 是策略迭代的最大轮数工具箱会在达到 max_iter 之前提前收敛所以设成 1000 足够设太大没有好处。返回值 policy 是一个 S×1 的向量第 s 个元素表示状态 s 下应该执行的动作编号V 是对应的状态值函数cpu_time 是求解耗时。第一次跑如果看到 policy 全部是 1 或者全部是 2不用意外那是奖励设置导致的结果不是 bug。验证策略是否合理的一个简单方法把 V 打印出来看趋势随着状态变差值函数应该递减如果状态 4 的 V 反而比状态 1 高说明奖励或转移概率出了问题。对比值迭代[V_vi, policy_vi, cpu_time_vi] mdp_value_iteration(P, R, discount);值迭代不需要传 max_iter工具箱内部会根据收敛容差自行判断何时停止。策略迭代和值迭代的一句话区别策略迭代每次都做完整的策略评估然后一次性改进策略值迭代则是在评估一步之后立刻做改进交替进行。对状态规模小、动作数少的模型差别不大状态数上千时策略迭代的迭代次数明显更少但每次迭代要做线性方程组求解单次开销大。实际项目中我通常会先用值迭代跑一个结果作为对照再和策略迭代的输出对比两个结果一致才放心。4. 拆开工具箱看实现校验、变体、有限期和调试输出4.1 mdp_check 与 mdp_check_square_stochastic求解前的闸门压缩包里有 mdp_check.m 和 mdp_check_square_stochastic.m 两个校验函数名字类似职责不同。mdp_check_square_stochastic 只检查状态转移矩阵的合法性P 必须是 S×S×A 的三维数组每个动作对应的 S×S 子矩阵所有元素非负且每一行之和都等于 1。mdp_check 是更全面的检查除了转移矩阵还会校验奖励矩阵的形状、数值范围和维度一致性。我把这两个函数当作一切求解操作的入口闸门。无论 P 和 R 是自己写的还是从外部数据转换来的先跑一遍mdp_check(P, R); mdp_check_square_stochastic(P);如果矩阵有问题这两个函数会直接报错并给出错误类型的提示省去后面策略迭代跑出奇怪结果的排查时间。常见的一个坑奖励矩阵中意外出现 NaN 时mdp_check 能查出非有限数值但如果出现的是 Inf早期版本的工具箱不一定能拦得住Inf 会顺着值函数计算扩散到策略结果里。因此我养成了在调用 mdp_check 之后额外加一句话的习惯assert(all(isfinite(R(:))), Reward matrix contains Inf or NaN);这一行 assert 的成本几乎为零但能避免最隐蔽的一类数值污染问题。4.2 mdp_policy_iteration_modified初始策略偏置的作用mdp_policy_iteration_modified 是标准政策迭代的变体它允许传入一个初始策略。标准版本默认从全 0 策略开始不需要额外参数修改版本则要先给一个初始策略然后从这个策略出发做评估和改进。用途在哪儿如果业务上已经有一个启发式策略比如“状态 1 和 2 不做动作状态 3 和 4 维护”你把它作为初始策略传进去观察算法在它基础上做出的改进能够很直观地看出原策略在哪些状态上不是最优的。另外对某些状态空间巨大的问题从一个好的初始策略开始能少走几轮策略评估的弯路。init_policy [1; 1; 2; 2]; % 业务给定的启发式策略长度必须等于 S [policy_m, V_m] mdp_policy_iteration_modified(P, R, discount, init_policy);这里 init_policy 的每个元素是动作编号取值必须在 1 到 A 之间否则 mdp_check 会给初始策略的合法性做校验。使用这个函数时特别注意在离散 MDP 中最优策略是确定性策略所以初始策略本身是否随机并不影响最终收敛结果但它会影响收敛过程中的中间步骤数。如果传一个随机初始策略收敛路径可能变长但不影响最终输出的一致性。4.3 mdp_finite_horizon有限期模型的输出结构和 horizon 设置mdp_finite_horizon.m 解决的是有限期 MDP也就是决策只在固定步数内有效的问题。典型的例子补货周期还剩 5 天、设备还能再运营 3 个批次。这类问题的最优策略不是平稳的——同一个状态在“还有 5 天”和“还有 1 天”时的最优动作可能完全不同越接近结束期策略越偏向短期收益。horizon 5; [V_fh, policy_fh] mdp_finite_horizon(P, R, discount, horizon);这里有一个特别容易踩坑的地方policy_fh 的输出形状是 S×horizon不是 S×1。每一列对应一个时间步也就是说 policy_fh(:, t) 是“距离结束还有 t 步时”每个状态的最优动作。如果你不看文档直接把它当普通策略向量用MATLAB 会给出维度不匹配的错误或者更麻烦自动广播导致结果完全错位。我在第一次使用这个函数时就是没注意输出的维度把一整块策略矩阵直接送进后续的计算逻辑浪费了小半天时间。horizon 的参数选择没有绝对标准但有一个实用原则horizon 越大前几个时间步的策略就越接近无限期最优策略。如果你取 horizon100前 50 步的策略基本和策略迭代算出来的无限期平稳策略一致后 50 步才开始出现“末期效应”。所以如果你手头的问题其实没有明确的结束点就不要为了套 mdp_finite_horizon 随意设一个大的 horizon直接使用 mdp_policy_iteration 更合理。4.4 mdp_verbose 与 mdp_silent把迭代过程打开看调试策略迭代不收敛时最难判断的是问题出在矩阵定义还是算法参数上。MDPtoolbox 提供了两个辅助函数mdp_verbose 和 mdp_silent用来控制求解器是否打印每一轮迭代的内部信息。开启 verbose 之后策略迭代会输出每轮评估过后值函数的更新幅度、策略是否发生改变以及当前轮数。这比盲改参数要直观得多。mdp_verbose(1); % 打开详细输出 [policy, V] mdp_policy_iteration(P, R, discount); mdp_silent(); % 恢复正常静默模式我在调试阶段基本都会开 verbose分散注意力也是准确的印证尤其是在需要给导师或客户展示迭代过程时它能直观证明算法确实是一轮一轮收敛的。提交结果的时候再恢复 silent避免不必要的刷屏输出。还有一个辅助函数 mdp_span计算值函数在状态维度上的跨度最大值减最小值可以用来判断值迭代是否停止。5. 避坑跑策略迭代最容易翻车的 5 个场景5.1 现象矩阵维度不一致一调用就报错这是初次使用工具箱时几乎必现的错误。现象是调用 mdp_policy_iteration 或 mdp_check 时报错提示索引超出矩阵维度或维度不匹配错误位置往往指向求解函数内部而不是你自己的脚本容易让人误以为是工具箱的 bug。原因P 被定义成了 S×S 而不是 S×S×A或者 R 被定义成了 S×S 而不是 S×A。前者的常见诱因是只写了一个动作的转移矩阵忘了补第三个维度后者常见于把奖励直接写成矩阵形式却遗漏了动作维。解决在调用任何求解函数之前先加一个形状自检assert(ndims(P) 3 size(P, 3) A, P must be S×S×A); assert(size(R, 1) S size(R, 2) A, R must be S×A);这两行断言能拦截掉九成以上的维度问题。之后再用 mdp_check 做一次完整校验解决概率接近百分之百。5.2 现象策略迭代不收敛V 值无限膨胀现象是运行时间很长返回的 V 向量里出现数量级极大的值甚至 Infpolicy 则在一个状态上反复切换。原因最常见的是 discount 设置成了 1 或者非常接近 1同时状态转移矩阵存在闭环。折扣因子为 1 时(I - γPπ) 可能变成奇异矩阵线性方程组解不出来迭代就会发散。另一个可能原因是 P 矩阵中有一行概率和不为 1导致每个周期都在放大值函数。解决把 discount 降到 0.99 以下尽量使用 0.95 或 0.9验证方法是在求解前计算cond(eye(S) - discount * squeeze(P(:,:,1)))如果条件数巨大说明 P 或 discount 有问题。这也与 2.2 节提到的线性解有关。5.3 现象策略和手算结果对不上差一个动作现象是工具箱返回的 policy 和教材或手推结果在一个状态上有差异反复检查模型设置也找不到逻辑错误。原因大概率是并列 Q 值。两个动作在某个状态下的动作值函数差异小到浮点精度以下工具箱按索引顺序取了较小的动作编号。4.2 节提到的初始策略索引顺序也与此相关。解决把 Q 值矩阵算出来手动比对两个动作之间的差异。在 MATLAB 中计算R_plus_gamma zeros(S, A); for a 1:A R_plus_gamma(:, a) R(:, a) discount * squeeze(P(:, :, a)) * V; end % 查看状态 2 的两个动作值差异 disp(R_plus_gamma(2, :) - max(R_plus_gamma(2, :)));如果差异小于 1e-8说明数值上等价策略选哪一个都合理不必纠结。5.4 现象mdp_finite_horizon 的结果与预期方向相悖现象是同样一组 P 和 R用有限期求解得到的策略和无限期策略迭代结果差别很大甚至在某个状态上选择了明显“更差”的动作。原因这不是 bug而是有限期问题的本质属性。越靠近 horizon 最后一步策略越倾向于最大化即时奖励长期代价在剩余步数不足时被忽略。比如库存问题中剩下最后一天时选择不补货以避免超储是正确决策而并非错误。解决先问自己的问题是否有明确的结束期限。有明确期限就用 mdp_finite_horizon并理解输出是 S×horizon 的逐期策略矩阵每一期策略需要按时间戳索引没有明确期限就使用 mdp_policy_iteration 求解无限期模型两种工具不能混用。5.5 现象文档引用的函数当前版本找不到现象是按压缩包里的 html 文档调用某个函数比如 mdp_spanMATLAB 提示未定义函数或变量。原因工具包包含多个模块和历年版本的文档文档页面和 m 文件目录可能来自不同版本。mdp_span、mdp_computePR 等函数在部分版本中存在在另一部分版本中缺失或改名。文档索引与实际文件不一定同步。解决先用 which 或 dir 确认函数是否存在which mdp_span dir(*.m)如果确认当前目录没有这个文件就把文档里对应的公式自己实现一遍或者改用已有的同类函数。这类问题在开源工具包里很常见不用怀疑是自己配置错了也无需纠结版本号以实际文件清单为准。6. 把策略迭代用到实际决策库存补货模型与验证习惯6.1 把业务规则编码成 MDP 矩阵策略迭代能落地的典型场景之一是库存补货。假设一个零售商库存量就是状态 s取值范围 0 到 4补货数量是动作 a可取 0 或 2每天有随机需求需求数量服从一个已知的概率分布。于是状态转移矩阵 P 的每个元素就变成了当前库存 s 执行补货动作 a 之后需求量 d 消耗掉一部分库存下一状态 s 由 s a - d 决定并截断到边界 0 到 4。奖励矩阵 R 可以设计为销售收入减去补货成本、持货成本和缺货损失。这个建模过程中最花时间的是把业务参数翻译成 P 和 R而一旦翻译完成求解只是 mdp_policy_iteration 一行调用的事。6.2 最优策略的残差验证和模拟对比拿到最优策略后我从不直接采信输出结果每次都强制自己走一遍验证流程。第一步是把 policy 代回贝尔曼方程计算残差V_hat zeros(S, 1); for s 1:S a policy(s); V_hat(s) R(s, a) discount * (squeeze(P(s, :, a)) * V); end residual max(abs(V_hat - V)); disp(residual);残差应该等于 0如果大于 1e-6 就说明 V 和 policy 并不自洽需要回溯检查求解过程。第二步是做一次蒙特卡洛仿真用最优策略执行 100 步统计累计回报再跟一个简单启发式策略比如“库存低于 1 就补货”做对比。如果最优策略竟然比启发式还差那几乎可以断定 P 的某一个值写错了尤其要检查转移概率的行归一化是否在边界状态被破坏。从那以后我每次拿到新的求解结果都会强制过一遍“残差验证 模拟对比”这两步宁可多写十行检查代码也不愿意面对一个数值上正确、业务上不合理的最优策略。希望你拿到这套工具后也能先跑通校验再追逐算法输出的精确数字这样策略迭代才能真正帮你做决策。本文还有配套的精品资源点击获取