简介：这是一套基于MATLAB实现的马尔科夫决策过程编程资源，面向运筹学、强化学习初学者及需要求解动态规划问题的工程师，讲解了从状态空间、动作集合到转移概率和奖励函数的完整建模框架。压缩包内共二十七个文件，其中二十四个为脚本文件，三个为文本说明，整体大小约十七KB，体量精简便于快速部署实验。脚本实现涵盖价值迭代、策略迭代、Q值函数计算、网格世界环境构建等核心模块，并附有Sutton与Russell的经典算例演示；文本说明则用于介绍基本概念与来源信息，方便读者结合代码对照学习。目前已有1705人学习下载。通过运行这些程序，读者可以系统掌握状态转移矩阵与奖励函数的搭建方法，深入理解价值函数迭代与策略改进的数学推导，并可将该建模能力迁移到机器人路径规划、资源调度等实际决策场景中。
1. 当贝尔曼方程撞上MATLAB：MDP编程实现的第一步是建模不是算法
我第一次动手写马尔科夫决策过程（MDP）的代码时，以为难点在值迭代那个循环上。真正上手才发现，公式里一笔带过的状态转移，到了MATLAB里必须变成高维数组，维度对不齐，循环根本跑不起来。这份资源把MDP从理论公式完整翻译成可运行的MATLAB代码，从状态空间定义、转移概率存储、值迭代求解，到策略迭代与Q学习，每一段都有带注释的代码和参数说明。适合已有MDP概念、想亲手跑通一遍代码的初学者，也适合需要快速验证算法效果的研究者。网格世界是贯穿全文的示例，所有代码复制到MATLAB里就能直接改参数跑。
2. 从公式到张量：MDP四元素在MATLAB中的建模方式
2.1 状态空间与动作空间的离散化设计
MDP的五元组 (S, A, P, R, γ) 里，S和A都是抽象集合，写代码的第一步是给它们分配整数索引。常见的做法是把状态編成连续正整数，比如一个5×5的网格世界，状态数就是25。如果按行优先编号，任意一个网格位置 (r, c) 都能和标量索引互相转换：
% ---- 构建5x5网格世界的状态索引 ---- n_rows = 5; n_cols = 5; num_states = n_rows * n_cols; % 25个状态 % 行优先编号：s=1 对应 (1,1)，s=25 对应 (5,5) row_of = @(s) floor((s - 1) / n_cols) + 1; col_of = @(s) mod(s - 1, n_cols) + 1; state_of = @(r, c) (r - 1) * n_cols + c;这段代码里三个匿名函数是核心工具。row_of把一维状态索引还原成行号，col_of还原成列号，state_of把行列坐标编码回索引。用匿名函数而不是普通函数的好处是不占额外文件，网格规模调整时只需要改n_rows和n_cols，后面所有代码自动适配。
动作空间的离散化更直接。网格世界里我常用4个动作：上、下、左、右，编号1到4。这样做的好处是动作编号可以直接用作转移概率元胞数组的索引：
num_actions = 4; % 1:上 2:下 3:左 4:右 % 也可以用一个向量存动作名，方便调试时打印 action_labels = {'up', 'down', 'left', 'right'};动作编号本身没有任何数学含义，纯粹是给编程用。只要保证后面填充P矩阵时，P{1}永远对应动作1，P{2}永远对应动作2，就不会混乱。我在实际项目中习惯把动作编号写进注释，或者用action_labels这类变量辅助可视化，避免过几天回来看代码还要重新猜。
2.2 转移概率与奖励：元胞数组和稀疏矩阵的取舍
转移概率 P(s,a,s') 是一个三维张量。初学者最容易犯的错误是直接开一个num_states * num_actions * num_states的全矩阵，25个状态时也就25*4*25=2500个元素，看似不痛不痒，但一旦状态规模到几千，内存立刻爆掉。我一般用MATLAB的元胞数组嵌套稀疏矩阵，外层维度是动作，内层每个元素是num_states × num_states的稀疏转移矩阵：
% ---- 转移概率存储：元胞数组 + 稀疏矩阵 ---- P = cell(num_actions, 1); for a = 1:num_actions P{a} = sparse(num_states, num_states); end % 填充动作2（向下移动）的确定性转移 for s = 1:num_states r = row_of(s); c = col_of(s); if r < n_rows s_next = state_of(r + 1, c); P{2}(s, s_next) = 1; % 向下移动一格，概率为1 else P{2}(s, s) = 1; % 已在底边界，原地停留 end end这段代码是典型的确定性转移写法：P{2}(s, s_next)=1 表示从状态s执行动作2后，以概率1落到s_next。边界条件的处理也在这里体现——当Agent在最低一行尝试向下走时，我没有让它穿出网格，而是让它停留在原状态。这是MDP实现里最常见的边界约定，既能避免非法状态索引，又不破坏概率归一化。
为什么不直接用三维数组P(s, a, s_next)？两个原因。第一，稀疏矩阵只存储非零元素，对网格世界这种大多数转移只涉及相邻状态的问题，稀疏表示的内存占用是线性规模的而不是立方规模。第二，后续值迭代里我们需要反复计算P{a} * V，稀疏矩阵乘稀疏向量在MATLAB里经过高度优化，比全矩阵的squeeze(P(:,a,:)) * V快一个量级。
奖励矩阵的存储相对简单。当奖励只依赖状态和动作而不依赖下一状态时，用num_states × num_actions的普通矩阵即可：
% ---- 奖励定义 ---- R = zeros(num_states, num_actions); % 每走一步的步代价 for a = 1:num_actions R(:, a) = -0.1; end % 陷阱状态：进入即获得较大的负奖励 trap = state_of(2, 3); R(trap, :) = -5; % 目标状态：到达后终止，奖励为0 target = state_of(5, 5); R(target, :) = 0;把奖励设计成每个动作一样的负值，会让Agent倾向于尽快到达终点而不是在原地兜圈子。陷阱状态设置 -5 的负奖励，则是在教它绕开危险区域。实际项目中如果奖励只依赖于状态而不依赖于动作，也可以直接用num_states × 1的列向量，然后在代码里通过repmat扩展。不过我在项目里更偏好矩阵形式，因为一旦后续想设计“向某个方向走有额外惩罚”这类动作相关奖励，矩阵形式不需要改动代码结构。
到这里，MDP的四个核心元素在MATLAB里全部有了明确的物理载体：S是整数索引，A是整数编号，P是元胞数组内嵌稀疏矩阵，R是二维矩阵。狄拉克符号写的抽象公式已经落地成可计算的数据结构。
3. 值迭代：手写贝尔曼方程收敛循环的核心代码
3.1 值迭代的数学拆解与收敛性判定
值迭代的更新公式是：
V_{k+1}(s) = max_a [ R(s,a) + γ ∑_{s'} P(s'|s,a) V_k(s') ]
这个公式的含义是：状态s的价值等于在最优动作下，立即奖励加上折扣后的未来期望价值。迭代开始前给V一个任意初始值，每轮都对所有状态做一次更新，直到V不再变化，就认为收敛到了最优价值函数。
收敛性判定在代码里体现为相邻两轮V的最大绝对差小于阈值。MATLAB里我习惯用无穷范数：
diff = norm(V_new - V, Inf);用无穷范数而不是二范数，是因为我们关心的是最坏情况下单个状态的价值变化。如果最坏的那个状态都稳定了，那整体自然就稳了。
3.2 完整实现：代码注释、参数设置和运行结果
值迭代的主循环代码不长，但每一步都对应贝尔曼方程的一个运算：
% ---- 值迭代主循环 ---- gamma = 0.9; % 折扣因子 V = zeros(num_states, 1); Q = zeros(num_states, num_actions); max_iter = 200; tol = 1e-6; for iter = 1:max_iter % 对每个动作计算 Q(s,a) = R(s,a) + gamma * sum(P(s'|s,a) * V(s')) for a = 1:num_actions Q(:, a) = R(:, a) + gamma * (P{a} * V); end % 贝尔曼最优更新：取所有动作中的最大值 V_new = max(Q, [], 2); % 收敛判定 diff = norm(V_new - V, Inf); V = V_new; if diff < tol fprintf('值迭代收敛：%d 轮，残差 %.2e\n', iter, diff); break; end end这里最需要理解的是P{a} * V这一步。P{a}是 num_states × num_states 的稀疏矩阵，V是 num_states × 1 的列向量。按矩阵乘法规则，结果的第s行正好是 ∑_{s'} P(s'|s,a) V(s')，也就是从状态s出发执行动作a后，下一状态价值按概率加权求和。这一步做完，Q(:,a)的每一行都是在说“如果我在此状态选择动作a，能拿到多少期望回报”。之后max(Q, [], 2)对每一行取最大值，完成max_a操作。
参数设置的逻辑参考如下：gamma = 0.9 是常用的中间值，太大（比如0.99）会拖慢收敛，太小（比如0.5）会让Agent变得短视。tol = 1e-6 是判定阈值，比这个更小意义不大，因为浮点运算本身有精度极限。max_iter 设置成200，是给最坏情况一个逃生门——如果gamma很接近1或者网格很大，200轮可能才勉强够用，但配合tol已经足够正常工作。
初始值 V = zeros 是最常见的做法。严格从理论讲，初始值只要是有界向量，最终都会收敛到同一个最优价值函数，所以不必担心初始值带来的偏差。唯一需要注意的是，如果状态空间中有终止状态，应该让终止状态的价值保持为0，否则终止状态会通过转移概率把错误的价值传播给前驱状态。
收敛之后提取策略也有一个标准动作：
% ---- 从Q矩阵提取最优策略 ---- [~, policy] = max(Q, [], 2); % 打印成网格形状，便于直观检查 policy_grid = reshape(policy, [n_rows, n_cols]); disp('策略图（1:上 2:下 3:左 4:右）：'); disp(policy_grid);这个policy向量的每个元素就是该状态下的最优动作编号。把它 reshape 成网格形状后，一眼就能看出Agent在哪些位置会选择哪条路径。如果发现某些状态的动作明显不合理，先检查奖励设置和转移矩阵，而不要先怀疑迭代公式。
4. 策略迭代与Q学习：另外两条路线的MATLAB对照
4.1 策略迭代：线性方程组求解与策略稳定判定
值迭代每次只更新一步V，轮数往往较多。策略迭代走的是另一条路：既然最终要的是最优策略，不如先固定策略，把这个策略下的V解出来，再做策略改进，直到策略不再变化。
策略评估需要求解的线性方程组是：
V = R_π + γ P_π V
其中 R_π 和 P_π 是由当前策略 π 确定的联合奖励向量和联合转移矩阵。在MATLAB里用\运算符直接解这个方程：
% ---- 策略迭代主循环 ---- policy = ones(num_states, 1); % 初始策略：所有状态选动作1 for iter = 1:100 % 策略评估：组装该策略下的 R_pi 和 P_pi R_pi = zeros(num_states, 1); P_pi = sparse(num_states, num_states); for s = 1:num_states a = policy(s); R_pi(s) = R(s, a); P_pi(s, :) = P{a}(s, :); end % 解线性方程组 V = R_pi + gamma * P_pi * V V = (speye(num_states) - gamma * P_pi) \ R_pi; % 策略改进：用当前V重新计算Q，选新策略 Q = zeros(num_states, num_actions); for a = 1:num_actions Q(:, a) = R(:, a) + gamma * (P{a} * V); end [~, new_policy] = max(Q, [], 2); % 策略稳定则收敛 if isequal(new_policy, policy) fprintf('策略迭代收敛：%d 轮\n', iter); break; end policy = new_policy; end这里speye(num_states)是稀疏单位矩阵，减去gamma * P_pi后做反斜杠求解，得到的是当前策略下的精确价值函数，而不是像值迭代那样逐步逼近。策略改进一步则利用这个精确V重新计算每个状态的最优动作。
策略迭代和值迭代的实际差异体现在轮数上。网格世界这种小规模问题，策略迭代通常5到10轮就收敛，值迭代可能要几十轮甚至上百轮。但策略迭代每一轮都要解一次线性方程组，单个轮次的耗时远高于值迭代的一次矩阵乘法。我自己的习惯是：状态数在几百以内且转移矩阵稀疏度较高时，用策略迭代；状态数大、稀疏度一般时，值迭代的轮次优势更明显。
4.2 Q学习：无模型场景下的实践变体
前面两种方法都假设转移概率 P 和奖励 R 已知，属于动态规划求解类别。如果实际环境给不出 P，只有交互接口，Q学习才是出路。Q学习的更新公式不再累加所有可能转移，而是基于单次采样：
Q(s,a) ← Q(s,a) + α [ r + γ max_{a'} Q(s',a') - Q(s,a) ]
MATLAB实现里，为了运行需要先模拟一个环境，我这里直接复用前面定义的P和R采样，但请注意真正的Q学习场景里这两者是黑匣子：
% ---- Q学习：无模型场景下的采样式更新 ---- alpha = 0.1; % 学习率 epsilon = 0.1; % 探索率 Q = zeros(num_states, num_actions); for episode = 1:500 % 每个回合随机重置初始状态，避免样本分布偏斜 s = randi(num_states); for step = 1:200 % epsilon-greedy 动作选择 if rand < epsilon a = randi(num_actions); % 探索 else [~, a] = max(Q(s, :)); % 利用 end % 环境模拟：按 P{a}(s,:) 采样下一状态 trans_probs = full(P{a}(s, :)); s_next = find(cumsum(trans_probs) >= rand, 1); r = R(s, a); % Q更新 Q(s, a) = Q(s, a) + alpha * (r + gamma * max(Q(s_next, :)) - Q(s, a)); s = s_next; end end这段代码里最关键的是cumsum那一行：它在按离散分布采样——如果P{a}(s,:)是 [0.2, 0.3, 0.5]，cumsum变成 [0.2, 0.5, 1.0]，然后任取一个 [0,1] 均匀随机数，落进哪个区间就选中哪个状态。
Q学习参数的意义与值迭代完全不同。alpha = 0.1决定每次更新的步长，太小则收敛慢，太大则Q值在目标附近来回震荡。epsilon = 0.1控制探索概率，平衡“尝试新动作”和“利用当前最优动作”。我在实际项目中还遇到过需要动态调整alpha的场景：早期用0.2加速探索，后期降到0.01精调Q值。这种学习率衰减机制在MATLAB里实现并不复杂，但对收敛稳定性帮助很大。
5. 避坑指南：MATLAB实现MDP的高频翻车点
5.1 转移概率矩阵行未归一化导致价值函数膨胀
现象：值迭代每轮V都在增大，残差始终降不下去，最后出现NaN。
原因：手写转移矩阵时，某些行的非零元素之和大于1或小于1。大于1时，每一轮迭代都在人为放大未来期望回报，V自然膨胀；小于1时，V又会异常偏低。尤其注意从环境采样得到的转移概率，直接把频次写进矩阵而没有归一化。
解决：填充完P之后强制做一次行归一化，用sum(P{a}, 2)得到每行之和，逐行除掉：
for a = 1:num_actions row_sums = sum(P{a}, 2); % 防止除以0：全零行保持全零 row_sums(row_sums == 0) = 1; P{a} = P{a} ./ row_sums; end这个检查我每次都写，已经成为条件反射。哪怕理论上确定性转移每行只有1，也跑一遍归一化，成本极低但能挡掉后续几乎所有诡异现象。
5.2 折扣因子取1，值迭代不收敛
现象：迭代到max_iter也没有触发diff < tol，V仍缓慢变化。
原因：γ=1意味着未来回报和当前回报权重完全一样。在没有终止状态或者终止状态没有被特殊处理时，Agent可以无限期获得步代价，总回报发散，价值函数不存在有限固定点。
解决：将gamma设为小于1的值，0.9到0.99之间都能正常工作。如果业务需求必须用γ=1，需要给MDP添加吸收终止状态，并且确保终止状态的转移和奖励定义能让价值有界。
5.3 元胞数组访问符写反：P(2) 与 P{2} 的维度错乱
现象：代码提示维度不一致，或者Q矩阵变成元胞数组无法参与乘法运算。
原因：P(2)从元胞数组中取出的是“一个包含索引为2的子元胞”，尺寸仍是1×1的cell，而不是里面的矩阵。P{2}才取出真正的稀疏转移矩阵。
解决：把P(a)全部改成P{a}，同时在写循环时用a遍历1:num_actions而不是1:numel(P)，再从花括号取矩阵。检查代码时重点看所有索引P的地方，一秒钟就能定位。
这个问题看起来低级，但我在调试别人的代码时多次见过。MATLAB的花括号和圆括号语义和其它语言不同，写得快了确实容易混。
5.4 初始状态固定，Q学习只学到了局部策略
现象：收敛后策略在部分状态明显不合理，比如某个方向明明能更快到达终点却频繁绕路。
原因：Q学习循环里s = randi(num_states)如果被写成了s = 1，那么每个episode都从同一个状态出发，远离起点的那部分状态很难被访问，样本分布严重不均，经验不足的Q值更新只能收敛到局部最优。
解决：每个episode开始都重置随机初始状态。更深一层的做法是把epsilon调大一点，比如0.2，让Agent有更多概率随机探索。使用rng shuffle固定随机种子前要注意复现需求：需要可复现时用固定种子，否则让RNG自由初始化。
5.5 步代价和终止奖励量级失衡，策略变得保守
现象：Agent几乎原地不动，或者绕开一切非终止状态，策略图全是同一个动作。
原因：步代价设成 -0.1，终点奖励却只有0，那么Agent每多走一步都会累计负收益，到达终点也只能止损。理性策略当然是原地不动。
解决：把步代价改成 -0.01，或者把终点奖励提高到 1 以上，让“到达终点的净收益”明显优于“原地徘徊”。奖励量级通常要差一个数量级以上。设计阶段先算一下最优路径的步数和每一步的代价，再做一次头脑验证，能省很多调参时间。
6. 让迭代跑得更快：贝尔曼残差检查与一套我的调试习惯
算法跑通之后，真正要养成的好习惯是验证结果而不只是看收敛日志。我常用的验证手段有两个，都很简单但非常有效。
第一个是计算贝尔曼残差。值迭代收敛后，把求得的V代回贝尔曼最优方程右侧，重新算一次Q，再取最大值得到V_hat，然后计算max(abs(V_hat - V))。残差如果在1e-6量级，说明V确实满足了贝尔曼最优性，迭代不是恰好卡在某个局部波动：
% ---- 收敛后的贝尔曼残差验证 ---- Q_check = zeros(num_states, num_actions); for a = 1:num_actions Q_check(:, a) = R(:, a) + gamma * P{a} * V; end V_hat = max(Q_check, [], 2); residual = max(abs(V_hat - V)); fprintf('贝尔曼残差：%.3e\n', residual);这个步骤相当于给整个MDP求解过程留了一张后悔药。如果残差远大于1e-6，别急着调参，先回去检查P是不是归一化了，gamma是不是小于1，奖励矩阵的维度是不是对得上。
第二个习惯是策略不一致性检查。值迭代收敛后提取一次策略，再随机用别的初始值重新跑一遍，对比两次策略是否一致。如果一致，说明策略搜索空间里没有第二个歧义解；如果不一致，则要检查是不是存在多个等价最优策略，或者收敛判定过松导致提前停止。
从那以后，我每次拿到别人的MDP代码或者写完一段新的MDP逻辑，都强制自己依次走完这几步：先查行归一化，再确认gamma范围，跑完值迭代后做贝尔曼残差验证，最后对比初始策略重跑的一致性。这套流程帮我过滤掉了绝大部分看似诡异、实则低级的错误。希望帮到你。
本文还有配套的精品资源，点击获取
