馬可夫決策過程與動態規劃

馬可夫決策過程 Markov Decision Process

強化學習的問題用馬可夫決策過程（MDP）描述，由五個部分組成：

• 狀態集合 \(\mathcal S\)、動作集合 \(\mathcal A\)
• 轉移機率 \(P(s' \mid s, a)\)：在狀態 s 做動作 a 之後，到達 s′ 的機率
• 報酬 \(r(s, a, s')\)
• 折扣因子 \(\gamma \in [0, 1]\)：未來的報酬打幾折

「馬可夫」指下一個狀態只取決於現在的狀態與動作，與更早的歷史無關，這正是狀態空間模型中「狀態」的意義。代理人的策略 \(\pi(a \mid s)\) 決定在每個狀態做什麼動作，目標是讓報酬總和（return）的期望值最大：

價值函數與 Bellman 方程式 Value Functions

狀態 s 的價值是從 s 出發、照策略 π 走下去的報酬總和期望值；動作價值 Q 則是先做動作 a、之後再照 π 走：

因為 \(G_t = r_{t+1} + \gamma G_{t+1}\)，價值滿足遞迴關係，也就是 Bellman 方程式：

最佳策略的價值 \(V^*\) 滿足 Bellman 最佳方程式：每一步都選讓「立即報酬＋折扣後的未來價值」最大的動作：

動態規劃：值迭代與策略迭代 Dynamic Programming

模型（P 與 r）已知時，可以直接解 Bellman 方程式：

• 值迭代：從任意的 V 開始，反覆把 Bellman 最佳方程式的右邊當成新的 V。γ < 1 時這是壓縮映射，一定收斂到 \(V^*\)。
• 策略迭代：輪流做「策略評估」（解出目前策略的 \(V^\pi\)）與「策略改善」（對 \(V^\pi\) 取貪婪動作），通常幾輪就收斂。

範例：4×3 網格世界（Russell 與 Norvig《人工智慧：現代方法》的經典例子）。代理人從左下角出發，(3, 2) 是 +1 的出口、(3, 1) 是 −1 的出口，(1, 1) 是牆。動作會「滑」：往意圖方向的機率 0.8，往左右兩側各 0.1，撞牆留在原地。每走一步報酬 −0.04，γ = 1：

# 值迭代：4×3 網格世界（Russell & Norvig 的經典例子）
# 每走一步報酬 r，(3,2) 終點 +1、(3,1) 終點 −1，(1,1) 是牆；往意圖方向走的機率 0.8，往左右兩側各 0.1
import numpy as np

ACTS = {"↑": (0, 1), "→": (1, 0), "↓": (0, -1), "←": (-1, 0)}
WALL, TERM = (1, 1), {(3, 2): 1.0, (3, 1): -1.0}
STATES = [(x, y) for y in range(3) for x in range(4) if (x, y) != WALL]

def move(s, a):
    n = (s[0] + ACTS[a][0], s[1] + ACTS[a][1])
    return n if 0 <= n[0] < 4 and 0 <= n[1] < 3 and n != WALL else s      # 撞牆就留在原地

def q(V, s, a, r, gamma=1.0):
    names = list(ACTS); i = names.index(a)
    outcomes = [(a, 0.8), (names[(i - 1) % 4], 0.1), (names[(i + 1) % 4], 0.1)]
    return sum(p * (r + gamma * V[move(s, b)]) for b, p in outcomes)

def value_iteration(r):
    V = {s: TERM.get(s, 0.0) for s in STATES}
    for it in range(1, 1000):
        newV = {s: V[s] if s in TERM else max(q(V, s, a, r) for a in ACTS) for s in STATES}   # Bellman 最佳方程式
        delta = max(abs(newV[s] - V[s]) for s in STATES)
        V = newV
        if delta < 1e-10:
            break
    pi = {s: max(ACTS, key=lambda a: q(V, s, a, r)) for s in STATES if s not in TERM}
    return V, pi, it

for r in (-0.04, -2.0):
    V, pi, it = value_iteration(r)
    print(f"每步報酬 r = {r}：{it} 次迭代收斂")
    for y in (2, 1, 0):
        vals = "  ".join("  牆  " if (x, y) == WALL else f"{V[(x, y)]:6.3f}" for x in range(4))
        arrows = " ".join("■" if (x, y) == WALL else ("+1" if (x, y) == (3, 2) else "-1" if (x, y) == (3, 1) else pi[(x, y)]) for x in range(4))
        print(f"  {vals}    {arrows}")

⬇ 下載 value_iteration.py（Python 3＋numpy）

實際執行的輸出：

每步報酬 r = -0.04：40 次迭代收斂
   0.812   0.868   0.918   1.000    → → → +1
   0.762    牆     0.660  -1.000    ↑ ■ ↑ -1
   0.705   0.655   0.611   0.388    ↑ ← ← ←
每步報酬 r = -2.0：35 次迭代收斂
  -7.043  -4.230  -1.730   1.000    → → → +1
  -9.543    牆    -3.570  -1.000    ↑ ■ → -1
  -10.815  -8.474  -5.974  -3.775    → → → ↑

每步報酬 −0.04 時，算出的價值與教科書的數值完全相同（0.812、0.868、0.918……0.388，verify_rl.py 核對）。注意 (2, 0) 往左走、(2, 1) 往上走：寧可繞遠路，也不冒 10% 的機率滑進 −1。若每步要付 −2，留在場上太貴，(2, 1) 乾脆直接跳進 −1 結束。策略迭代只要 5 輪就得到同樣的策略。值迭代的收斂與 γ 有關：γ = 1 要 40 次，γ = 0.9 只要 31 次（下圖左）。

Bellman 方程式與 LQR Bellman Equation and LQR

把同樣的想法用在連續狀態的線性系統 \(\mathbf x_{k+1} = \mathbf A\mathbf x_k + \mathbf B\mathbf u_k\)，成本（負報酬）為 \(\sum_k (\mathbf x_k^{\mathsf T}\mathbf Q\mathbf x_k + \mathbf u_k^{\mathsf T}\mathbf R\mathbf u_k)\)。Bellman 方程式成為

猜 \(V(\mathbf x) = \mathbf x^{\mathsf T}\mathbf P\mathbf x\)，代入後對 \(\mathbf u\) 微分令為 0，得到最佳回授 \(\mathbf u = -\mathbf K\mathbf x\)，而 \(\mathbf P\) 必須滿足

這正是離散時間 LQR 的代數 Riccati 方程式。從 \(\mathbf P = 0\) 反覆套用右邊（也就是值迭代），就會收斂到 Riccati 方程式的解。以模糊控制章的 MCK 系統為例（取樣 0.1 秒）：

# 用 Bellman 方程式（動態規劃）解 LQR：MCK 系統（m = 1、c = 1、k = 2，同模糊控制章）離散化後，
# 從 V = 0 開始反覆套用 Bellman 方程式，就是 Riccati 迭代，結果和 scipy 解離散代數 Riccati 方程式相同
import numpy as np
from scipy.linalg import expm, solve_discrete_are

dt = 0.1
Ac, Bc = np.array([[0., 1.], [-2., -1.]]), np.array([[0.], [1.]])
M = expm(np.block([[Ac, Bc], [np.zeros((1, 3))]]) * dt)                 # 零階保持離散化
A, B = M[:2, :2], M[:2, 2:]
Q, R = np.diag([10., 1.]), np.array([[0.1]])                           # 成本 Σ (xᵀQx + uᵀRu)

# V(x) = xᵀPx 代入 V(x) = min_u [xᵀQx + uᵀRu + V(Ax + Bu)]，最小化 u 後得到 P 的更新式
P = np.zeros((2, 2))
for k in range(1, 201):
    K = np.linalg.solve(R + B.T @ P @ B, B.T @ P @ A)                   # 這一步的最佳回授 u = −Kx
    P_new = Q + A.T @ P @ (A - B @ K)
    if np.abs(P_new - P).max() < 1e-9 * np.abs(P_new).max():
        P = P_new
        break
    P = P_new
print(f"Bellman 迭代 {k} 次收斂，P =\n{P.round(4)}")
P_are = solve_discrete_are(A, B, Q, R)
print("與 solve_discrete_are 的最大差：", f"{np.abs(P - P_are).max():.1e}")
K = np.linalg.solve(R + B.T @ P @ B, B.T @ P @ A)
print("最佳回授增益 K =", K.round(4), "，閉迴路特徵值大小", np.abs(np.linalg.eigvals(A - B @ K)).round(4))

x = np.array([1.0, 0.0])                                               # 從位移 1、速度 0 出發
for t in range(30):
    x = (A - B @ K) @ x
x_open = np.linalg.matrix_power(A, 30) @ np.array([1.0, 0.0])
print(f"3 秒後 x = {x.round(4)}（不加控制時 x = {x_open.round(4)}）")

⬇ 下載 lqr_bellman.py（Python 3＋numpy＋scipy）

實際執行的輸出：

Bellman 迭代 42 次收斂，P =
[[56.7106  8.2748]
 [ 8.2748  4.8288]]
與 solve_discrete_are 的最大差： 5.9e-08
最佳回授增益 K = [[6.2624 3.5862]] ，閉迴路特徵值大小 [0.7707 0.7707]
3 秒後 x = [-0.0001  0.0016]（不加控制時 x = [-0.2131  0.2483]）

所以 LQR 可以看成強化學習在「模型已知、線性、二次成本」時的特例，答案是解析的。當模型未知、系統非線性，或報酬不是二次式時，就需要下一頁不靠模型的學習方法。