简介：这是一套面向机器学习与强化学习初学者的马尔科夫决策过程（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_state=s' | state=s, action=a)。这个注释值十行代码。
折扣因子 γ 在 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×A，R: S×A，policy: S×1，tol: 收敛容差 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 和 R：4 状态 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 的核心要求是：对每个动作 a，P(:,:,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 越大，前几个时间步的策略就越接近无限期最优策略。如果你取 horizon=100，前 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 向量里出现数量级极大的值，甚至 Inf，policy 则在一个状态上反复切换。
原因：最常见的是 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_span），MATLAB 提示未定义函数或变量。
原因：工具包包含多个模块和历年版本的文档，文档页面和 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 的某一个值写错了，尤其要检查转移概率的行归一化是否在边界状态被破坏。从那以后我每次拿到新的求解结果，都会强制过一遍“残差验证 + 模拟对比”这两步，宁可多写十行检查代码，也不愿意面对一个数值上正确、业务上不合理的最优策略。希望你拿到这套工具后，也能先跑通校验，再追逐算法输出的精确数字，这样策略迭代才能真正帮你做决策。
本文还有配套的精品资源，点击获取
