""" ============================================================================= 第 2 讲 · 代码 1：马尔可夫决策过程（MDP）的基本表示 ============================================================================= 对应讲义：lectures/lecture02_mdp_planning.md 第 1、2、3 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 第 1 讲的 MRP 只有 (S, P, R, γ)：火星车只能被动地随机游走。 这一讲我们把"动作"加进来，得到 MDP (S, A, P, R, γ)：火星车可以**选择**往哪走。 在写任何算法（策略评估、策略迭代、价值迭代）之前，得先回答一个很"工程"的问题： **一个 MDP 在代码里长什么样？一个策略在代码里长什么样？** 本文件回答这个问题，并完成本讲最关键的一个观察： **MDP + 一个固定的策略 = 一个 MRP** 所以第 1 讲学过的 MRP 价值算法，可以直接拿来评估 MDP 中的策略。 ----------------------------------------------------------------------------- 【本文件的数据约定（第 2 讲所有代码通用）】 ----------------------------------------------------------------------------- N = 状态数，A = 动作数 P : 形状 (N, A, N) 的三维数组，P[s, a, s2] = 在 s 执行 a 后到达 s2 的概率 R : 形状 (N, A) 的二维数组，R[s, a] = 在 s 执行 a 获得的期望即时奖励 策略有两种写法： 确定性策略 policy : 长度为 N 的整数数组，policy[s] = 在 s 选择的动作 随机策略 pi : 形状 (N, A) 的数组， pi[s, a] = 在 s 选择 a 的概率 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 01_mdp_basics.py 【下标约定】代码中的状态 i <==> 讲义中的 s_{i+1} """ import numpy as np np.set_printoptions(precision=3, suppress=True) # 动作编号与显示用的名字、箭头 TRY_LEFT = 0 TRY_RIGHT = 1 ACTION_NAMES = ["TryLeft", "TryRight"] ARROWS = ["←", "→"] # ============================================================================= # 第 1 部分：构建 MDP # ============================================================================= def build_mars_rover_mdp(success_prob=0.8): """ 构建讲义 1.3 节的火星车 MDP。 【模型】 - 7 个状态（格子），2 个动作（TryLeft=0, TryRight=1）。 - 执行动作后：以 success_prob 的概率朝想去的方向移动一格， 否则原地不动（打滑）。撞墙也原地不动。 - 奖励只取决于所在的格子（与动作无关）：s1 为 +1，s7 为 +10，其余为 0。 所以 R[s, TryLeft] == R[s, TryRight]。 【参数】 success_prob : float，范围 [0, 1]。 1.0 → 确定性世界（CS234 课件中的版本）； 0.0 → 火星车永远动不了（任何策略都一样，一个有趣的退化情况）。 【返回】 (P, R)：P 形状 (7, 2, 7)，R 形状 (7, 2)。 【Python 小知识】三维数组 np.zeros((7, 2, 7)) 创建一个"7 × 2 × 7"的数组，可以想象成： 7 张表（每个状态一张），每张表有 2 行（每个动作一行），每行 7 个数（到每个状态的概率）。 P[s, a] 取出的是长度为 7 的一维数组——"在 s 执行 a 后，下一状态的概率分布"。 P[:, a] 取出的是 7×7 的矩阵——"所有状态都执行动作 a 时"的转移矩阵。 冒号 : 表示"这一维全都要"。 """ if not (0.0 <= success_prob <= 1.0): raise ValueError(f"success_prob 必须在 [0, 1] 之间，但收到了 {success_prob}") n_states, n_actions = 7, 2 P = np.zeros((n_states, n_actions, n_states)) for s in range(n_states): for a in range(n_actions): direction = -1 if a == TRY_LEFT else +1 # 目标格子；撞墙时 min/max 会把它"夹"回边界，即原地不动 target = min(max(s + direction, 0), n_states - 1) # 注意这里用 += 而不是 =：撞墙时 target == s， # "成功"和"失败"两部分概率都要加到同一个格子上，总和才是 1。 # 如果写成 =，后一句会覆盖前一句，撞墙时这一行的概率和就只剩 0.2 了——一个很隐蔽的 bug。 P[s, a, target] += success_prob P[s, a, s] += 1.0 - success_prob reward_per_state = np.array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0]) # 【Python 小知识】np.tile(x, (1, 2)) 把 x 在第二个方向上复制 2 份。 # 这里先用 reward_per_state[:, None] 把形状 (7,) 变成 (7, 1)（None 表示"插入一个新维度"）， # 再复制成 (7, 2)：两个动作的奖励相同。 R = np.tile(reward_per_state[:, None], (1, n_actions)) return P, R def check_mdp(P, R, atol=1e-8): """ 检查 (P, R) 是否构成一个合法的 MDP。不合法就抛出 ValueError。 【检查内容】 1. P 是三维数组，形状 (N, A, N)；R 是二维数组，形状 (N, A)。 2. P 所有元素 >= 0。 3. 对**每一对** (s, a)，P[s, a, :] 之和等于 1。 （MRP 里是"每一行"和为 1；MDP 里是"每个状态的每个动作"和为 1。） 4. P、R 中没有 NaN 或无穷大。 【参数】 atol : float，浮点误差容忍度。 【边界】 - A = 0（没有动作）或 N = 0（没有状态）都视为非法。 """ P = np.asarray(P, dtype=float) R = np.asarray(R, dtype=float) if P.ndim != 3 or P.shape[0] != P.shape[2]: raise ValueError(f"P 的形状必须是 (N, A, N)，但收到了 {P.shape}") n, n_actions, _ = P.shape if n == 0 or n_actions == 0: raise ValueError("状态数和动作数都必须 >= 1") if R.shape != (n, n_actions): raise ValueError(f"R 的形状必须是 {(n, n_actions)}，但收到了 {R.shape}") # np.isfinite 检查"是有限的正常数字"（不是 NaN、不是 ±inf） if not (np.isfinite(P).all() and np.isfinite(R).all()): raise ValueError("P 或 R 中含有 NaN 或无穷大") if (P < -atol).any(): raise ValueError("P 中有负的概率") # P.sum(axis=2)：沿最后一维（下一状态）求和，结果形状 (N, A)，每个元素应当是 1 sums = P.sum(axis=2) if not np.allclose(sums, 1.0, atol=atol): # np.argwhere 返回所有"不满足条件"的 (s, a) 下标对 bad = np.argwhere(~np.isclose(sums, 1.0, atol=atol)) raise ValueError(f"以下 (状态, 动作) 的转移概率之和不等于 1：{bad.tolist()}") # ============================================================================= # 第 2 部分：策略的表示 # ============================================================================= def deterministic_to_stochastic(policy, n_actions): """ 把确定性策略（每个状态一个动作）转换成随机策略矩阵（每个状态一个概率分布）。 【思路】 确定性策略是随机策略的特例：选中的动作概率为 1，其余为 0（"独热编码 one-hot"）。 统一用 (N, A) 矩阵表示策略后，后面的算法只需要写一个版本，就能同时处理两种策略。 例：policy = [1, 1, 0]，n_actions = 2 → pi = [[0, 1], [0, 1], [1, 0]] 【参数】 policy : 长度为 N 的整数数组，每个元素在 0 到 n_actions-1 之间。 n_actions : int，动作数。 【返回】 numpy 数组，形状 (N, n_actions)。 【边界】 - policy 里出现越界的动作编号 → ValueError。 """ policy = np.asarray(policy, dtype=int) if ((policy < 0) | (policy >= n_actions)).any(): # | 是"逐元素或"。注意 numpy 数组之间要用 | 和 &，而不是 Python 的 or 和 and。 raise ValueError(f"policy 中的动作必须在 0 到 {n_actions - 1} 之间，但收到了 {policy}") n = len(policy) pi = np.zeros((n, n_actions)) # 【Python 小知识】花式索引：pi[行下标数组, 列下标数组] = 1 # np.arange(n) = [0, 1, ..., n-1]，policy = [a_0, a_1, ...]， # 这一行等价于 for s in range(n): pi[s, policy[s]] = 1.0 pi[np.arange(n), policy] = 1.0 return pi def check_stochastic_policy(pi, n_states, n_actions, atol=1e-8): """ 检查随机策略矩阵是否合法：形状 (N, A)、非负、每行和为 1。不合法就抛出 ValueError。 """ pi = np.asarray(pi, dtype=float) if pi.shape != (n_states, n_actions): raise ValueError(f"策略矩阵形状必须是 {(n_states, n_actions)}，但收到了 {pi.shape}") if (pi < -atol).any() or not np.allclose(pi.sum(axis=1), 1.0, atol=atol): raise ValueError("策略矩阵每一行必须是概率分布（非负且和为 1）") def format_policy(policy): """把确定性策略 [0, 1, 1] 显示成 "s1:← s2:→ s3:→" 这样的字符串。""" return " ".join(f"s{s + 1}:{ARROWS[a]}" for s, a in enumerate(policy)) def count_deterministic_policies(n_states, n_actions): """ 一共有多少个确定性（平稳）策略？答案是 A^N。 【思路】 每个状态独立地从 A 个动作里挑一个，N 个状态的选择互不影响， 由乘法原理，总数 = A × A × ... × A（N 个）= A^N。 【边界】 Python 的整数没有上限，所以即使 A^N 是天文数字也能精确算出来（不会溢出）。 这一点和 C/Java 等语言不同。 """ return n_actions ** n_states # ============================================================================= # 第 3 部分：MDP + 策略 = MRP # ============================================================================= def mdp_to_mrp_loop(P, R, pi): """ 【教学版】用最直白的 for 循环，把"MDP + 策略"变成一个 MRP。 【思路】（讲义 3 节） 固定策略 π 之后，"选哪个动作"不再是决策，而是一个已知的随机过程。 把动作的随机性"平均掉"： R_π(s) = Σ_a π(a|s) · R(s, a) 在 s 的期望即时奖励 P_π(s'|s) = Σ_a π(a|s) · P(s'|s, a) 从 s 到 s' 的总概率 直觉：到达 s' 有好几条"路"（通过不同的动作）， 每条路的概率 = 选这个动作的概率 × 这个动作把你带到 s' 的概率， 把所有路的概率加起来（全概率公式）。 【参数】 P : (N, A, N)；R : (N, A)；pi : (N, A) 随机策略矩阵。 【返回】 (P_pi, R_pi)：P_pi 形状 (N, N)，R_pi 形状 (N,)。 """ n, n_actions, _ = P.shape P_pi = np.zeros((n, n)) R_pi = np.zeros(n) for s in range(n): for a in range(n_actions): R_pi[s] += pi[s, a] * R[s, a] for s_next in range(n): P_pi[s, s_next] += pi[s, a] * P[s, a, s_next] return P_pi, R_pi def mdp_to_mrp(P, R, pi): """ 【实用版】与 mdp_to_mrp_loop 完全相同的计算，但用 numpy 一行完成，快得多。 【Python 小知识】np.einsum（爱因斯坦求和） einsum 用一个字符串描述"哪些下标相乘、哪些下标求和"： np.einsum("sa,sat->st", pi, P) 读法：第一个数组的下标叫 s,a；第二个叫 s,a,t；结果保留 s,t。 **在箭头右边消失的下标（这里是 a）会被求和**。所以它计算的是： 结果[s, t] = Σ_a pi[s, a] * P[s, a, t] 正是 P_π 的定义。同理 "sa,sa->s" 计算 Σ_a pi[s,a] * R[s,a]。 einsum 第一次看会有点懵，但它是把数学公式翻译成代码最直接的方式： 公式里有哪些下标、对哪个求和，照着写就行。 如果不放心，就像本文件这样，先写一个 for 循环版本，再对比两者结果是否一致。 """ P_pi = np.einsum("sa,sat->st", pi, P) R_pi = np.einsum("sa,sa->s", pi, R) return P_pi, R_pi # ============================================================================= # 第 4 部分：在 MDP 中按策略采样一条轨迹 # ============================================================================= def sample_episode(P, R, pi, start_state, horizon, rng): """ 按策略 π 与 MDP 交互 horizon 步，返回 [(s_0, a_0, r_0), (s_1, a_1, r_1), ...]。 【思路】每一步做两次"掷骰子"： 1. 智能体掷骰子选动作： a ~ π(·|s) 2. 世界掷骰子决定下一状态：s' ~ P(·|s, a) 对比第 1 讲的 MRP 采样只有第 2 步——这就是"加入动作"在代码层面的全部区别。 【参数】 pi : (N, A) 随机策略矩阵（确定性策略请先用 deterministic_to_stochastic 转换）。 start_state : int，起点。 horizon : int，步数，必须 >= 1。 rng : numpy 随机数生成器。 【返回】 list，每个元素是 (状态, 动作, 奖励) 三元组。 """ if horizon < 1: raise ValueError(f"horizon 必须 >= 1，但收到了 {horizon}") n, n_actions, _ = P.shape s = start_state episode = [] for _ in range(horizon): a = int(rng.choice(n_actions, p=pi[s])) # 智能体的随机性 r = float(R[s, a]) episode.append((s, a, r)) s = int(rng.choice(n, p=P[s, a])) # 世界的随机性 return episode # ============================================================================= # 第 5 部分：演示 # ============================================================================= def main(): P, R = build_mars_rover_mdp(success_prob=0.8) n, n_actions, _ = P.shape # --------------------------------------------------------------------- print("=" * 78) print("演示 1：火星车 MDP 的转移概率——每个动作各有一张 7×7 的转移矩阵") print("=" * 78) check_mdp(P, R) print(f" P 的形状 = {P.shape}（状态, 动作, 下一状态），R 的形状 = {R.shape}（状态, 动作）") print(" ✓ 通过合法性检查\n") for a in range(n_actions): print(f" 动作 {ACTION_NAMES[a]}：P[:, {a}, :] =") print(" " + " ".join(f"s{j + 1}" for j in range(n))) for s in range(n): print(f" s{s + 1} " + " ".join(f"{x:.1f}" for x in P[s, a])) print() print(" 注意 s1 执行 TryLeft、s7 执行 TryRight 时撞墙，所以原地不动的概率是 1.0。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 2：三种策略及其矩阵表示") print("=" * 78) always_left = np.zeros(n, dtype=int) # 全是 0 = TryLeft always_right = np.ones(n, dtype=int) # 全是 1 = TryRight random_pi = np.full((n, n_actions), 1.0 / n_actions) # np.full(形状, 值)：用同一个值填满 print(f" 一直向左：{format_policy(always_left)}") print(f" 一直向右：{format_policy(always_right)}") print(f" 随机策略：每个状态都以 0.5 的概率选择每个动作") print(f"\n 一直向右转成随机策略矩阵（前 3 行）：\n{deterministic_to_stochastic(always_right, n_actions)[:3]}") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 3：MDP + 策略 = MRP（计算 P_π 和 R_π）") print("=" * 78) pi_right = deterministic_to_stochastic(always_right, n_actions) P_pi, R_pi = mdp_to_mrp(P, R, pi_right) print(" 策略\"一直向右\"下的 P_π（它就是一个普通的马尔可夫链转移矩阵）：") print(P_pi) print(f" R_π = {R_pi}") P_rand, R_rand = mdp_to_mrp(P, R, random_pi) print("\n 随机策略下 P_π 的第 4 行（从 s4 出发）：", P_rand[3]) print(" 解读：0.5×(0.8 向左) = 0.4 到 s3；0.5×(0.8 向右) = 0.4 到 s5；两个动作各有 0.2 打滑，") print(" 0.5×0.2 + 0.5×0.2 = 0.2 留在 s4。——这恰好就是第 1 讲的火星车马尔可夫链！") # 验证 einsum 版本和 for 循环版本完全一致 P_loop, R_loop = mdp_to_mrp_loop(P, R, random_pi) print(f"\n for 循环版与 einsum 版的最大差异：P {np.abs(P_loop - P_rand).max():.1e}，" f"R {np.abs(R_loop - R_rand).max():.1e}（一致 ✓）") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 4：按策略采样一条轨迹（随机策略，从 s4 出发，8 步）") print("=" * 78) rng = np.random.default_rng(0) for s, a, r in sample_episode(P, R, random_pi, start_state=3, horizon=8, rng=rng): print(f" 在 s{s + 1} 选择 {ACTION_NAMES[a]:<8s} 获得奖励 {r:+.0f}") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 5：确定性策略一共有多少个？A^N 增长得有多快？") print("=" * 78) for (ns, na, desc) in [(7, 2, "火星车"), (16, 4, "4×4 网格世界"), (100, 4, "10×10 网格世界"), (1000, 10, "中等规模问题")]: count = count_deterministic_policies(ns, na) # len(str(count)) 是这个整数的位数 shown = f"{count}" if count < 10**6 else f"约 10^{len(str(count)) - 1}" # 名字放在行尾：中文在终端里占两格宽，放行首用 {:<12s} 对齐会歪 print(f" N={ns:5d}, A={na:3d} → {shown:>12s} 个确定性策略 （{desc}）") print(" → 除了极小的问题，把所有策略挨个试一遍是不可能的。所以我们需要策略迭代、价值迭代。") # --------------------------------------------------------------------- print("\n" + "=" * 78) print("演示 6：边界情况") print("=" * 78) bad_P = P.copy() bad_P[2, 1, 3] = 0.5 # 让 (s3, TryRight) 的概率和变成 0.2 + 0.5 = 0.7 try: check_mdp(bad_P, R) except ValueError as e: print(f" [概率和≠1] ValueError: {e}") try: deterministic_to_stochastic([0, 1, 2], n_actions=2) except ValueError as e: print(f" [非法动作] ValueError: {e}") P0, R0 = build_mars_rover_mdp(success_prob=0.0) P0_pi, _ = mdp_to_mrp(P0, R0, pi_right) print(f" [success_prob=0] P_π 是否为单位矩阵（永远原地不动）：{np.allclose(P0_pi, np.eye(n))}") if __name__ == "__main__": main()