强化学习——马尔可夫决策过程（MDP）【附 python 代码】

Python资料网推荐强化学习——马尔可夫决策过程（MDP）【附 python 代码】这篇文章给大家，欢迎收藏Python资料网享受知识的乐趣

一、马尔可夫过程

过程	介绍
随机过程	在某时刻 t 的状态 S t S_t St​ 通常取决于 t 时刻之前的状态。状态 S t + 1 S_{t+1} St+1​ 的概率表示为： P ( S t + 1 ∣ S 1 , . . . , S t ) P(S_{t+1}|S_1,...,S_t) P(St+1​∣S1​,...,St​)
马尔可夫过程	某时刻 t 的状态只取决于上一个时刻 t-1 的状态。状态 S t + 1 S_{t+1} St+1​ 的概率表示为： P ( S t + 1 ∣ S t ) = P ( S t + 1 ∣ S 1 , . . . , S t ) P(S_{t+1}|S_t)=P(S_{t+1}|S_1,...,S_t) P(St+1​∣St​)=P(St+1​∣S1​,...,St​)

  马尔可夫过程也被称为马尔可夫链，通常用元组 $<S,P>$ 来描述，其中 S 是有限数量的状态集合，P 是状态转移矩阵。假设有 n 个状态，则 S = { s 1 , s 2 , . . . , s n } S=\{s_1,s_2,...,s_n\} S={s1​,s2​,...,sn​} ,
P = [ P ( s 1 ∣ s 1 ) ⋯ P ( s n ∣ s 1 ) ⋮ ⋱ ⋮ P ( s 1 ∣ s n ) ⋯ P ( s n ∣ s n ) ] P=\begin{bmatrix} P(s_1|s_1) & \cdots & P(s_n|s_1) \\ \vdots & \ddots & \vdots\\ P(s_1|s_n) & \cdots & P(s_n|s_n) \end{bmatrix} P= ​P(s1​∣s1​)⋮P(s1​∣sn​)​⋯⋱⋯​P(sn​∣s1​)⋮P(sn​∣sn​)​ ​
  矩阵 P 中第 i 行第 j 列元素 P ( s j ∣ s i ) = P ( S t + 1 = s j ∣ S t = s i ) P(s_j|s_i)=P(S_{t+1}=s_j|S_t=s_i) P(sj​∣si​)=P(St+1​=sj​∣St​=si​) 表示从状态 s i s_i si​ 转移到状态 s j s_j sj​ 的概率，称 P ( s ′ ∣ s ) P(s'|s) P(s′∣s) 为状态转移函数。从某个状态出发，到达其他状态的概率和必须为 1 。即状态转移矩阵 P 的每一行和为 1 。

示例：

S = { S i , 1 ≤ i ≤ 6 }    状态转移矩阵    P = [ 0.9 0.1 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0.6 0 0.4 0 0 0 0 0.3 0.7 0 0.2 0.3 0.5 0 0 0 0 0 0 0 1 ] S=\{S_i,1\le i \le 6\} \\ \; \\ 状态转移矩阵\;P= \begin{bmatrix} 0.9 & 0.1 & 0 & 0 & 0 & 0 \\ 0.5 & 0 & 0.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.6 & 0 & 0.4 \\ 0 & 0 & 0 & 0 & 0.3 & 0.7 \\ 0 & 0.2 & 0.3 & 0.5 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} S={Si​,1≤i≤6}状态转移矩阵P= ​0.90.50000​0.10000.20​00.5000.30​000.600.50​0000.300​000.40.701​ ​

二、马尔可夫奖励过程

2.1 定义

  马尔可夫奖励过程由 $<S,P,r,\gamma >$ 构成，其中：

• S 是有限状态的集合
• P 是状态转移矩阵
• r 是奖励函数，$r(s)$ 是指转移到该状态时可以获得的奖励期望
• $\gamma$ 是折扣因子，取值范围为：$[0,1]$ 。引入折扣因子是因为远期利益具有一定的不确定性，有时希望能尽快获得有些奖励，所以需要对远期利益打一些折扣。接近 1 则更关注长期的累积奖励，接近 0 则更关注短期奖励。

2.2 回报

  回报是指从状态 S t S_t St​ 开始，一直到终止状态，所有奖励的衰减之和。具体公式如下：
G t = ∑ k = 0 ∞ γ k R t + k G_t=\sum^{\infty}_{k=0}\gamma^kR_{t+k} Gt​=k=0∑∞​γkRt+k​
  其中 R t R_t Rt​ 是指在时刻 t 获得的奖励。

示例：

【状态旁的数字表示进入该状态获得的奖励】

新建项目 MDP，项目结构如下：

在 MDP 目录下新建文件 MRP.py ，文件代码如下：

import numpy as np

np.random.seed(0)

P = [
    [0.9, 0.1, 0.0, 0.0, 0.0, 0.0],
    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4],
    [0.0, 0.0, 0.0, 0.0, 0.3, 0.7],
    [0.0, 0.2, 0.3, 0.5, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
]
P = np.array(P)
rewards = [-1, -2, -2, 10, 1, 0]
gamma = 0.5


def compute_return(start, chain, gamma):
    G = 0
    for i in reversed(range(start, len(chain))):
        G = gamma * G + rewards[chain[i] - 1]
    return G


chain = [1, 2, 3, 6]
start = 0
G = compute_return(start, chain, gamma)
print('计算回报为：%s' % G)

运行 MRP.py 文件，运行结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\MRP.py 
计算回报为：-2.5

进程已结束，退出代码为 0

2.3 价值函数

• 价值：一个状态的期望回报。【就是给回报取个别名叫价值】
• 价值函数：输入为某个状态，输出为该状态的价值。
  • 公式为： V ( s ) = E [ G t ∣ S t = s ] V(s)=E[G_t|S_t=s] V(s)=E[Gt​∣St​=s] ，等价于 V ( s ) = E [ R t + γ V ( S t + 1 ) ∣ S t = s ] V(s)=E[R_t+\gamma V(S_{t+1})|S_t=s] V(s)=E[Rt​+γV(St+1​)∣St​=s] 。【即：当前状态的价值 = 当前状态获得的奖励 + 折扣因子 $\ast$ 下一个状态的价值。下一个状态不是确定的状态，而是取决于转移矩阵 P 】
  • 一方面， E [ R t ∣ S t = s ] = r ( s ) E[R_t|S_t=s]=r(s) E[Rt​∣St​=s]=r(s) ；另一方面， E [ γ V ( S t + 1 ) ∣ S t = s ] E[\gamma V(S_{t+1})|S_t=s] E[γV(St+1​)∣St​=s] 可以从状态 s 出发的转移概率得到，即： E [ γ V ( S t + 1 ) ∣ S t = s ] = γ ∑ s ′ ∈ S P ( s ′ ∣ s ) V ( s ′ ) E[\gamma V(S_{t+1})|S_t=s] = \gamma \sum_{s'\in S}P(s'|s)V(s') E[γV(St+1​)∣St​=s]=γs′∈S∑​P(s′∣s)V(s′)所以
    V ( s ) = r ( s ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s ) V ( s ′ ) V(s)=r(s)+\gamma \sum_{s'\in S}P(s'|s)V(s') V(s)=r(s)+γs′∈S∑​P(s′∣s)V(s′)这就是有名的贝尔曼方程，对每一个状态都成立。
  • 马尔可夫奖励过程有 n 个状态，即 S = { s 1 , s 2 , ⋯   , s n } S=\{s_1,s_2,\cdots,s_n\} S={s1​,s2​,⋯,sn​} ，所有状态的价值表示成一个列向量 V = [ V ( s 1 ) , V ( s 2 ) , ⋯   , V ( s n ) ] T V=[V(s_1),V(s_2),\cdots,V(s_n)]^T V=[V(s1​),V(s2​),⋯,V(sn​)]T ，同理，奖励函数列向量 R = [ r ( s 1 ) , r ( s 2 ) , ⋯   , r ( s n ) ] T R=[r(s_1),r(s_2),\cdots,r(s_n)]^T R=[r(s1​),r(s2​),⋯,r(sn​)]T 。
    于是，贝尔曼方程写成矩阵的形式：$V=R+\gamma PV$ ，解得： V = ( E − γ P ) − 1 R V=(E-\gamma P)^{-1}R V=(E−γP)−1R ，其中 E 是单位矩阵。
    【解析解复杂度为 O ( n 3 ) O(n^3) O(n3) ，只适合求解小规模的马尔可夫奖励过程，大规模的可以使用 DP算法，蒙特卡洛方法，时序差分算法】

示例：

修改文件 MRP.py ，文件代码如下：

import numpy as np

np.random.seed(0)

P = [
    [0.9, 0.1, 0.0, 0.0, 0.0, 0.0],
    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4],
    [0.0, 0.0, 0.0, 0.0, 0.3, 0.7],
    [0.0, 0.2, 0.3, 0.5, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
]
P = np.array(P)
rewards = [-1, -2, -2, 10, 1, 0]
gamma = 0.5


def compute_return(start, chain, gamma):
    G = 0
    for i in reversed(range(start, len(chain))):
        G = gamma * G + rewards[chain[i] - 1]
    return G


def compute(P, rewards, gamma, states_num):
    rewards = np.array(rewards).reshape((-1, 1))  # 将 rewards 改写为列向量
    value = np.dot(np.linalg.inv(np.eye(states_num, states_num) - gamma * P), rewards)
    return value


def compute_return_sample():
    chain = [1, 2, 3, 6]
    start = 0
    G = compute_return(start, chain, gamma)
    print('计算回报为：%s' % G)


def compute_value_sample():
    V = compute(P, rewards, gamma, 6)
    print("马尔可夫奖励过程中每个状态的价值分别为：\n", V)


compute_value_sample()

运行结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\MRP.py 
马尔可夫奖励过程中每个状态的价值分别为：
 [[-2.01950168]
 [-2.21451846]
 [ 1.16142785]
 [10.53809283]
 [ 3.58728554]
 [ 0.        ]]

进程已结束，退出代码为 0

三、马尔可夫决策过程

  在马尔可夫奖励过程加上动作，就得到了马尔可夫决策过程（MDP）。MDP 由元组 $<S,A,P,r,\gamma >$ 构成，其中：

• S 是状态的集合
• A 是动作的集合
• $\gamma$ 是折扣因子
• $r(s,a)$ 是奖励函数，此时奖励同时取决于状态 s 和动作 a
• P ( s ′ ∣ s , a ) P(s'|s,a) P(s′∣s,a) 是状态转移函数，表示在状态 s 执行动作 a 之后到达状态 s ′ s' s′ 的概率

3.1 基本概念

• 策略 π ( a ∣ s ) = P ( A t = a ∣ S t = s ) \pi(a|s)=P(A_t=a|S_t=s) π(a∣s)=P(At​=a∣St​=s) ，表示在输入状态为 s 的情况下采取动作 a 的概率。
  • 确定性策略：在每个状态下只输出一个确定性的动作，即只有该动作的概率为 1 ，其他动作的概率为 0 。
  • 随机性策略：在每个状态下输出的是关于动作的概率分布，根据该分布采样得到一个动作。
• 状态价值函数 V π ( s ) = E π [ G t ∣ S t = s ] V^{\pi}(s)=E_{\pi}[G_t|S_t=s] Vπ(s)=Eπ​[Gt​∣St​=s] ，表示从状态 s 出发遵循策略 $\pi$ 能获得的期望回报
• 动作价值函数 Q π ( s , a ) = E π [ G t ∣ S t = s , A t = a ] Q^{\pi}(s,a)=E_{\pi}[G_t|S_t=s,A_t=a] Qπ(s,a)=Eπ​[Gt​∣St​=s,At​=a] ，表示遵循策略 $\pi$ 对当前状态 s 执行动作 a 得到的期望回报。
• 关系：
  • 状态 s 的价值等于在该状态下基于策略 $\pi$ 采取所有动作的概率与相应价值的乘积和，即： V π ( s ) = ∑ a ∈ A π ( a ∣ s )    Q π ( s , a ) V^{\pi}(s)=\sum_{a\in A}\pi(a|s)\;Q^{\pi}(s,a) Vπ(s)=a∈A∑​π(a∣s)Qπ(s,a)
  • 在状态 s 下采取动作 a 的价值等于即时奖励加上经过衰减的所有可能的下一个状态的转移概率与相应价值的乘积，即： Q π ( s , a ) = r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a )    V π ( s ′ ) Q^{\pi}(s,a)=r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)\;V^{\pi}(s') Qπ(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)

3.2 贝尔曼期望方程

• 状态价值函数：
  V π ( s ) = ∑ a ∈ A π ( a ∣ s )    Q π ( s , a ) = ∑ a ∈ A π ( a ∣ s ) [ r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a )    V π ( s ′ ) ] \begin{align*} V^{\pi}(s)&=\sum_{a\in A}\pi(a|s)\;Q^{\pi}(s,a) \\ & =\sum_{a\in A}\pi(a|s)[r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)\;V^{\pi}(s')] \end{align*} Vπ(s)​=a∈A∑​π(a∣s)Qπ(s,a)=a∈A∑​π(a∣s)[r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)]​
• 动作价值函数：
  Q π ( s , a ) = r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a )    V π ( s ′ ) = r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a )    ∑ a ∈ A π ( a ′ ∣ s ′ )    Q π ( s ′ , a ′ ) \begin{align*} Q^{\pi}(s,a) & =r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)\;V^{\pi}(s')\\ & =r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)\;\sum_{a\in A}\pi(a'|s')\;Q^{\pi}(s',a') \end{align*} Qπ(s,a)​=r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)=r(s,a)+γs′∈S∑​P(s′∣s,a)a∈A∑​π(a′∣s′)Qπ(s′,a′)​

【 V π ( s ) V^{\pi}(s) Vπ(s) 和 Q π ( s , a ) Q^{\pi}(s,a) Qπ(s,a) 互相代入就可以得到】

示例：

【白色圆圈表示状态，蓝色圆圈表示动作，动作旁边的数字表示奖励，虚线旁边没有数字表示概率为1】

在 MDP 目录下新建 MDP.py 文件，文件代码如下：

S = ['S1', 'S2', 'S3', 'S4', 'S5']

A = ['保持S1', '前往S1', '前往S2', '前往S3', '前往S4', '前往S5', '概率前往']

P = {'S1-保持S1-S1': 1.0, 'S1-前往S2-S2': 1.0,
     'S2-前往S1-S1': 1.0, 'S2-前往S3-S3': 1.0,
     'S3-前往S4-S4': 1.0, 'S3-前往S5-S5': 1.0,
     'S4-前往S5-S5': 1.0, 'S4-概率前往-S2': 0.2, 'S4-概率前往-S3': 0.4, 'S4-概率前往-S4': 0.4}

R = {'S1-保持S1': -1, 'S1-前往S2': 0,
     'S2-前往S1': -1, 'S2-前往S3': -2,
     'S3-前往S4': -2, 'S3-前往S5': 0,
     'S4-前往S5': 10, 'S4-概率前往': 1}

gamma = 0.5

MDP = (S, A, P, R, gamma)

# 策略 1 ：随机策略
Pi_1 = {'S1-保持S1': 0.5, 'S1-前往S2': 0.5,
        'S2-前往S1': 0.5, 'S2-前往S3': 0.5,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.5, 'S4-概率前往': 0.5}

# 策略 2 ：给定策略
Pi_2 = {'S1-保持S1': 0.6, 'S1-前往S2': 0.4,
        'S2-前往S1': 0.3, 'S2-前往S3': 0.7,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.1, 'S4-概率前往': 0.9}


def join(s1, s2):
    return s1 + '-' + s2

将 MDP 转化为 MRP ，公式如下：

1. r ′ ( s ) = ∑ a ∈ A π ( a ∣ s ) r ( s , a ) r'(s)=\sum_{a\in A}\pi (a|s)r(s,a) r′(s)=∑a∈A​π(a∣s)r(s,a)
2. p ′ ( s ′ ∣ s ) = ∑ a ∈ A π ( a ∣ s ) P ( s ′ ∣ s , a ) p'(s'|s)=\sum_{a\in A}\pi (a|s)P(s'|s,a) p′(s′∣s)=∑a∈A​π(a∣s)P(s′∣s,a)

示例：

修改 MDP.py 文件，文件代码如下：

import numpy
import numpy as np

S = ['S1', 'S2', 'S3', 'S4', 'S5']

A = ['保持S1', '前往S1', '前往S2', '前往S3', '前往S4', '前往S5', '概率前往']

P = {'S1-保持S1-S1': 1.0, 'S1-前往S2-S2': 1.0,
     'S2-前往S1-S1': 1.0, 'S2-前往S3-S3': 1.0,
     'S3-前往S4-S4': 1.0, 'S3-前往S5-S5': 1.0,
     'S4-前往S5-S5': 1.0, 'S4-概率前往-S2': 0.2, 'S4-概率前往-S3': 0.4, 'S4-概率前往-S4': 0.4}

R = {'S1-保持S1': -1, 'S1-前往S2': 0,
     'S2-前往S1': -1, 'S2-前往S3': -2,
     'S3-前往S4': -2, 'S3-前往S5': 0,
     'S4-前往S5': 10, 'S4-概率前往': 1}

gamma = 0.5

MDP = (S, A, P, R, gamma)

# 策略 1 ：随机策略
Pi_1 = {'S1-保持S1': 0.5, 'S1-前往S2': 0.5,
        'S2-前往S1': 0.5, 'S2-前往S3': 0.5,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.5, 'S4-概率前往': 0.5}

# 策略 2 ：给定策略
Pi_2 = {'S1-保持S1': 0.6, 'S1-前往S2': 0.4,
        'S2-前往S1': 0.3, 'S2-前往S3': 0.7,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.1, 'S4-概率前往': 0.9}


def join(s1, s2):
    return s1 + '-' + s2


# 从字典中取出对应状态的参数数组
# 例如：从策略1中取出状态 S1 的动作概率分布，即[0.5,0.5]
def get_state_parameter(x, s):
    p = []
    for i in x.keys():
        if i.split('-')[0] == s:
            p.append(x[i])
    return p


# 计算a,b两个向量的内积
def compute_sum(a, b):
    a = np.array(a)
    b = np.array(b)
    return np.inner(a, b)


# 将 MDP 的奖励函数转为 MRP 的奖励函数
# 思想：对于每个状态，根据策略将所有动作的概率进行加权，得到的奖励和就可以被认为是在 MRP 中该状态的奖励
def R_MDP_to_MRP(R, Pi, S):
    MRP_R = []
    for i in range(len(S)):
        MRP_R.append(compute_sum(get_state_parameter(Pi, S[i]), get_state_parameter(R, S[i])))
    return MRP_R


# 计算 MDP 的转移概率
def compute_P(P, Pi):
    for i in P.keys():
        P[i] *= Pi[join(i.split('-')[0], i.split('-')[1])]
    return P


# 根据转移概率创建 MRP 的转移矩阵
def set_MRP_P(P, S):
    MRP_P = np.zeros((len(S), len(S)))
    for i in P.keys():
        start_index = S.index(i.split('-')[0])
        end_index = S.index(i.split('-')[2])
        MRP_P[start_index][end_index] = P[i]
    MRP_P[len(S) - 1][len(S) - 1] = 1.0  # 终止状态设置转移概率
    return MRP_P


# 将 MDP 的转移函数转为 MRP 的转移矩阵
# 思想：对于每个状态转移到其他状态，计算策略的转移概率与状态转移函数的转移概率乘积和作为 MRP 的转移概率
def P_MDP_to_MRP(P, Pi, S):
    return set_MRP_P(compute_P(P, Pi), S)

使用 MRP 方法计算每个状态的价值：

示例：

新建 Main.py 文件，文件代码如下：

from MDP import *
from MRP import compute

R=R_MDP_to_MRP(R, Pi_1, S)
P=P_MDP_to_MRP(P, Pi_1, S)

print('使用策略 1，将 MDP 转化为 MRP')
print('转化后的 MRP 奖励函数：\n', R)
print('\n转化后的 MRP 状态转移矩阵：\n', P)

V=compute(P,R,gamma,len(S))
print("\nMDP 中每个状态价值分别为\n",V)

运行 Main.py 文件，运行结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\Main.py 
使用策略 1，将 MDP 转化为 MRP
转化后的 MRP 奖励函数：
 [np.float64(-0.5), np.float64(-1.5), np.float64(-1.0), np.float64(5.5), np.float64(0.0)]

转化后的 MRP 状态转移矩阵：
 [[0.5 0.5 0.  0.  0. ]
 [0.5 0.  0.5 0.  0. ]
 [0.  0.  0.  0.5 0.5]
 [0.  0.1 0.2 0.2 0.5]
 [0.  0.  0.  0.  1. ]]

MDP 中每个状态价值分别为
 [[-1.22555411]
 [-1.67666232]
 [ 0.51890482]
 [ 6.0756193 ]
 [ 0.        ]]

进程已结束，退出代码为 0

计算在状态 S 下采取动作 A 的价值：

修改 MDP.py 文件，文件代码如下：

import numpy
import numpy as np

S = ['S1', 'S2', 'S3', 'S4', 'S5']

A = ['保持S1', '前往S1', '前往S2', '前往S3', '前往S4', '前往S5', '概率前往']

P = {'S1-保持S1-S1': 1.0, 'S1-前往S2-S2': 1.0,
     'S2-前往S1-S1': 1.0, 'S2-前往S3-S3': 1.0,
     'S3-前往S4-S4': 1.0, 'S3-前往S5-S5': 1.0,
     'S4-前往S5-S5': 1.0, 'S4-概率前往-S2': 0.2, 'S4-概率前往-S3': 0.4, 'S4-概率前往-S4': 0.4}

R = {'S1-保持S1': -1, 'S1-前往S2': 0,
     'S2-前往S1': -1, 'S2-前往S3': -2,
     'S3-前往S4': -2, 'S3-前往S5': 0,
     'S4-前往S5': 10, 'S4-概率前往': 1}

gamma = 0.5

MDP = (S, A, P, R, gamma)

# 策略 1 ：随机策略
Pi_1 = {'S1-保持S1': 0.5, 'S1-前往S2': 0.5,
        'S2-前往S1': 0.5, 'S2-前往S3': 0.5,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.5, 'S4-概率前往': 0.5}

# 策略 2 ：给定策略
Pi_2 = {'S1-保持S1': 0.6, 'S1-前往S2': 0.4,
        'S2-前往S1': 0.3, 'S2-前往S3': 0.7,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.1, 'S4-概率前往': 0.9}


def join(s1, s2):
    return s1 + '-' + s2


# 从字典中取出对应状态的参数数组
# 例如：从策略1中取出状态 S1 的动作概率分布，即[0.5,0.5]
def get_state_parameter(x, s):
    p = []
    for i in x.keys():
        if i.split('-')[0] == s:
            p.append(x[i])
    return p


# 计算a,b两个向量的内积
def compute_sum(a, b):
    a = np.array(a)
    b = np.array(b)
    return np.inner(a, b)


# 将 MDP 的奖励函数转为 MRP 的奖励函数
# 思想：对于每个状态，根据策略将所有动作的概率进行加权，得到的奖励和就可以被认为是在 MRP 中该状态的奖励
def R_MDP_to_MRP(R, Pi, S):
    MRP_R = []
    for i in range(len(S)):
        MRP_R.append(compute_sum(get_state_parameter(Pi, S[i]), get_state_parameter(R, S[i])))
    return MRP_R


# 计算 MDP 的转移概率
def compute_P(P, Pi):
    P1 = P.copy()
    for i in P1.keys():
        P1[i] *= Pi[join(i.split('-')[0], i.split('-')[1])]
    return P1


# 根据转移概率创建 MRP 的转移矩阵
def set_MRP_P(P, S):
    MRP_P = np.zeros((len(S), len(S)))
    for i in P.keys():
        start_index = S.index(i.split('-')[0])
        end_index = S.index(i.split('-')[2])
        MRP_P[start_index][end_index] = P[i]
    MRP_P[len(S) - 1][len(S) - 1] = 1.0  # 终止状态设置转移概率
    return MRP_P


# 将 MDP 的转移函数转为 MRP 的转移矩阵
# 思想：对于每个状态转移到其他状态，计算策略的转移概率与状态转移函数的转移概率乘积和作为 MRP 的转移概率
def P_MDP_to_MRP(P, Pi, S):
    return set_MRP_P(compute_P(P, Pi), S)


# MDP = (0--S, 1--A, 2--P, 3--R, 4--gamma)
# 计算在状态 S 下采取动作 A 的价值 Q 
def compute_Q(s, a, MDP, V):
    r = MDP[3][join(s, a)]
    sum_PV = 0
    for i in range(len(MDP[0])):
        p = join(join(s, a), MDP[0][i])
        if p in MDP[2].keys():
            sum_PV += MDP[2][p] * V[i]

    return r + MDP[4] * sum_PV

修改 Main.py 文件，文件代码如下：

from MDP import *
from MRP import compute

R = R_MDP_to_MRP(R, Pi_1, S)
P = P_MDP_to_MRP(P, Pi_1, S)

print('使用策略 1，将 MDP 转化为 MRP')
print('转化后的 MRP 奖励函数：\n', R)
print('\n转化后的 MRP 状态转移矩阵：\n', P)

V = compute(P, R, gamma, len(S))
print("\nMDP 中每个状态价值分别为\n", V)

print("\n在状态为 S4 时采取动作 概率前往 的价值为：", compute_Q(S[3], A[6], MDP, V))

运行 Main.py 文件，运行结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\Main.py 
使用策略 1，将 MDP 转化为 MRP
转化后的 MRP 奖励函数：
 [np.float64(-0.5), np.float64(-1.5), np.float64(-1.0), np.float64(5.5), np.float64(0.0)]

转化后的 MRP 状态转移矩阵：
 [[0.5 0.5 0.  0.  0. ]
 [0.5 0.  0.5 0.  0. ]
 [0.  0.  0.  0.5 0.5]
 [0.  0.1 0.2 0.2 0.5]
 [0.  0.  0.  0.  1. ]]

MDP 中每个状态价值分别为
 [[-1.22555411]
 [-1.67666232]
 [ 0.51890482]
 [ 6.0756193 ]
 [ 0.        ]]

在状态为 S4 时采取动作 概率前往 的价值为： [2.15123859]

进程已结束，退出代码为 0

---

   MRP 解析解的方法在状态动作集合比较大的时候比较不适用。

四、蒙特卡洛方法

用蒙特卡洛方法估计 MDP 中的状态价值

  用策略在 MDP 上采样很多序列，计算从这个状态出发的回报再求其期望即可，公式如下： V π ( s ) = E π [ G t ∣ S t = s ] ≈ 1 N ∑ i = 1 N G t ( i ) V^{\pi}(s)=E_{\pi}[G_t|S_t=s]\approx\frac1N\sum^N_{i=1}G^{(i)}_t Vπ(s)=Eπ​[Gt​∣St​=s]≈N1​i=1∑N​Gt(i)​

具体过程如下：

1. 使用策略 $\pi$ 采样若干条序列：
   s 0 ( i ) ⟶ a 0 ( i ) r 0 ( i ) ,    s 1 ( i ) ⟶ a 1 ( i ) r 1 ( i ) ,    ⋯   ,    s T − 1 ( i ) ⟶ a T − 1 ( i ) r T − 1 ( i ) ,    s T ( i ) s^{(i)}_0\stackrel{a^{(i)}_0}{\longrightarrow}r^{(i)}_0,\;s^{(i)}_1\stackrel{a^{(i)}_1}{\longrightarrow}r^{(i)}_1,\;\cdots,\;s^{(i)}_{T-1}\stackrel{a^{(i)}_{T-1}}{\longrightarrow}r^{(i)}_{T-1},\;s^{(i)}_{T} s0(i)​⟶a0(i)​​r0(i)​,s1(i)​⟶a1(i)​​r1(i)​,⋯,sT−1(i)​⟶aT−1(i)​​rT−1(i)​,sT(i)​
2. 对每一条序列中的每一时间步 t 的状态 s 进行以下操作：
   • 更新状态 s 的计数器 $N(s)+=1$
   • 更新状态 s 的总回报 $M(s)+=G$
3. 每一个状态的价值被估计为回报的期望 $V(s)=M(s)/N(s)$
   实际中，采用增量式更新会更好，即：
   • $N(s)+=1$
   • V ( s ) + = 1 N ( s ) ( G − V ( s ) ) V(s)+=\frac{1}{N(s)}(G-V(s)) V(s)+=N(s)1​(G−V(s))

示例：

新建 MonteCarlo.py 文件，文件代码如下：

import numpy as np

from MDP import join


class MonteCarlo:
    # 采样序列
    @staticmethod
    def sample(MDP, Pi, timestep_max, number):
        S, A, P, R, gamma = MDP
        episodes = []
        for _ in range(number):
            episode = []
            timestep = 0
            s = S[np.random.randint(len(S) - 1)]  # 随机选择除一个终止状态外的状态作为起点
            # 一次采样 【到达终止状态或者到达最大时间步】
            while s != S[len(S) - 1] and timestep <= timestep_max:
                timestep += 1
                rand, temp = np.random.rand(), 0
                # 在状态 s 下根据策略选择动作
                for a_opt in A:
                    temp += Pi.get(join(s, a_opt), 0)
                    if temp > rand:
                        a = a_opt
                        r = R.get(join(s, a), 0)
                        break
                rand, temp = np.random.rand(), 0
                # 根据状态转移函数得到下一个状态
                for s_opt in S:
                    temp += P.get(join(join(s, a), s_opt), 0)
                    if temp > rand:
                        s_next = s_opt
                        break
                episode.append((s, a, r, s_next))
                s = s_next
            episodes.append(episode)
        return episodes

    # 计算价值
    @staticmethod
    def compute(episodes, V, N, gamma):
        for episode in episodes:
            G = 0
            # 从后往前计算
            for i in range(len(episode) - 1, -1, -1):
                (s, a, r, s_next) = episode[i]
                G = r + gamma * G
                N[s] += 1
                V[s] += (G - V[s]) / N[s]

修改 Main.py 文件，文件代码如下：

from MDP import *
from MRP import compute
from MonteCarlo import MonteCarlo

S = ['S1', 'S2', 'S3', 'S4', 'S5']

A = ['保持S1', '前往S1', '前往S2', '前往S3', '前往S4', '前往S5', '概率前往']

P = {'S1-保持S1-S1': 1.0, 'S1-前往S2-S2': 1.0,
     'S2-前往S1-S1': 1.0, 'S2-前往S3-S3': 1.0,
     'S3-前往S4-S4': 1.0, 'S3-前往S5-S5': 1.0,
     'S4-前往S5-S5': 1.0, 'S4-概率前往-S2': 0.2, 'S4-概率前往-S3': 0.4, 'S4-概率前往-S4': 0.4}

R = {'S1-保持S1': -1, 'S1-前往S2': 0,
     'S2-前往S1': -1, 'S2-前往S3': -2,
     'S3-前往S4': -2, 'S3-前往S5': 0,
     'S4-前往S5': 10, 'S4-概率前往': 1}

gamma = 0.5

MDP = (S, A, P, R, gamma)

# 策略 1 ：随机策略
Pi_1 = {'S1-保持S1': 0.5, 'S1-前往S2': 0.5,
        'S2-前往S1': 0.5, 'S2-前往S3': 0.5,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.5, 'S4-概率前往': 0.5}

# 策略 2 ：给定策略
Pi_2 = {'S1-保持S1': 0.6, 'S1-前往S2': 0.4,
        'S2-前往S1': 0.3, 'S2-前往S3': 0.7,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.1, 'S4-概率前往': 0.9}


def MDP_to_MRP():
    (S, A, P, R, gamma) = MDP
    R = R_MDP_to_MRP(R, Pi_1, S)
    P = P_MDP_to_MRP(P, Pi_1, S)

    print('使用策略 1，将 MDP 转化为 MRP')
    print('转化后的 MRP 奖励函数：\n', R)
    print('\n转化后的 MRP 状态转移矩阵：\n', P)

    V = compute(P, R, gamma, len(S))
    print("\nMDP 中每个状态价值分别为\n", V)
    print("\n在状态为 S4 时采取动作 概率前往 的价值为：", compute_Q(S[3], A[6], MDP, V))


# MDP_to_MRP(MDP, Pi_1)

def MC():
    timestep_max = 20
    num = 1000
    V = {'S1': 0, 'S2': 0, 'S3': 0, 'S4': 0, 'S5': 0}
    N = {'S1': 0, 'S2': 0, 'S3': 0, 'S4': 0, 'S5': 0}
    episodes = MonteCarlo.sample(MDP, Pi_1, timestep_max, num)
    print("采样前 5 条序列为：\n")
    for i in range(5):
        if i >= len(episodes):
            break
        print("序列 %d：%s" % (i + 1, episodes[i]))
    MonteCarlo.compute(episodes, V, N, gamma)
    print('\n使用蒙特卡洛方法计算 MDP 的状态价值为\n', V)


MC()

运行 Main.py 文件，运行结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\Main.py 
采样前 5 条序列为：

序列 1：[('S1', '前往S2', 0, 'S2'), ('S2', '前往S3', -2, 'S3'), ('S3', '前往S5', 0, 'S5')]
序列 2：[('S4', '概率前往', 1, 'S4'), ('S4', '前往S5', 10, 'S5')]
序列 3：[('S4', '前往S5', 10, 'S5')]
序列 4：[('S2', '前往S1', -1, 'S1'), ('S1', '保持S1', -1, 'S1'), ('S1', '前往S2', 0, 'S2'), ('S2', '前往S3', -2, 'S3'), ('S3', '前往S4', -2, 'S4'), ('S4', '前往S5', 10, 'S5')]
序列 5：[('S2', '前往S3', -2, 'S3'), ('S3', '前往S4', -2, 'S4'), ('S4', '前往S5', 10, 'S5')]

使用蒙特卡洛方法计算 MDP 的状态价值为
 {'S1': -1.2250047295780488, 'S2': -1.6935084269467038, 'S3': 0.48519918492538405, 'S4': 5.97601313133763, 'S5': 0}

进程已结束，退出代码为 0

五、占用度量

  定义 MDP 的初始状态分布为 v 0 ( s ) v_0(s) v0​(s) ，用 P t π ( s ) P^{\pi}_t(s) Ptπ​(s) 表示采取策略 $\pi$ 使得智能体在 t 时刻状态为 s 的概率，所以 P 0 π ( s ) = V 0 ( s ) P^{\pi}_0(s)=V_0(s) P0π​(s)=V0​(s)，然后定义一个策略的状态访问分布：

V π ( s ) = ( 1 − γ ) ∑ t = 0 ∞ γ t P t π ( s ) V^{\pi}(s)=(1-\gamma)\sum^{\infty}_{t=0}\gamma^tP^{\pi}_t(s) Vπ(s)=(1−γ)t=0∑∞​γtPtπ​(s)

  其中，$1-\gamma$ 是归一化因子。状态访问概率表示在一个策略下智能体和 MDP 交互会访问到的状态的分布。其有如下性质：
V π ( s ′ ) = ( 1 − γ ) V 0 ( s ′ ) + γ ∫ P ( s ′ ∣ s , a ) π ( a ∣ s ) V π ( s ) d s d a V^{\pi}(s')=(1-\gamma)V_0(s')+\gamma \int P(s'|s,a)\pi(a|s)V^{\pi}(s)dsda Vπ(s′)=(1−γ)V0​(s′)+γ∫P(s′∣s,a)π(a∣s)Vπ(s)dsda

【解释说明： ∑ t = 0 ∞ γ t      是等比数列求和，即： ∑ t = 0 ∞ γ t = lim ⁡ n → ∞ 1 − r n 1 − r = 1 1 − r ，所以 ( 1 − γ ) ∑ t = 0 ∞ γ t = 1 \sum^{\infty}_{t=0}\gamma^t\;\;是等比数列求和，即：\sum^{\infty}_{t=0}\gamma^t=\lim_{n\rightarrow\infty}\frac{1-r^n}{1-r}=\frac{1}{1-r}，所以(1-\gamma)\sum^{\infty}_{t=0}\gamma^t=1 t=0∑∞​γt是等比数列求和，即：t=0∑∞​γt=n→∞lim​1−r1−rn​=1−r1​，所以(1−γ)t=0∑∞​γt=1
注：这只是解释说明为什么 $1-\gamma$ 是归一化因子，不可认为
V π ( s ) = ( 1 − γ ) ∑ t = 0 ∞ γ t P t π ( s ) = 1 ⋅ P t π ( s ) = P t π ( s ) V^{\pi}(s)=(1-\gamma)\sum^{\infty}_{t=0}\gamma^tP^{\pi}_t(s)=1\cdot P^{\pi}_t(s)=P^{\pi}_t(s) Vπ(s)=(1−γ)t=0∑∞​γtPtπ​(s)=1⋅Ptπ​(s)=Ptπ​(s)】

  此外，定义策略的占用度量为：

ρ π ( s , a ) = V π ( s ) π ( a ∣ s ) = ( 1 − γ ) ∑ t = 0 ∞ γ t P t π ( s ) π ( a ∣ s ) \begin{align*} \rho^{\pi}(s,a) & =V^{\pi}(s)\pi(a|s) \\ & =(1-\gamma)\sum^{\infty}_{t=0}\gamma^tP^{\pi}_t(s)\pi(a|s) \end{align*} ρπ(s,a)​=Vπ(s)π(a∣s)=(1−γ)t=0∑∞​γtPtπ​(s)π(a∣s)​

  它表示状态动作对 $(s,a)$ 被访问到的概率。它有两个定理：

1. ρ π 1 = ρ π 2 ⇔ π 1 = π 2 \rho^{\pi_1}=\rho^{\pi_2}\Leftrightarrow \pi_1=\pi_2 ρπ1​=ρπ2​⇔π1​=π2​
2. 给定合法的占用度量 $\rho$ ，可生成该占用度量的唯一策略 π ρ \pi_{\rho} πρ​：
   π ρ = ρ ( s , a ) ∑ a ′ ρ ( s , a ′ ) \pi_{\rho}=\frac{\rho(s,a)}{\sum_{a'}\rho(s,a')} πρ​=∑a′​ρ(s,a′)ρ(s,a)​

示例：

MDP.py 文件新增函数 occupancy：

# 计算状态动作（s, a）出现的频率，以此估计策略的占用度量
def occupancy(episodes, s, a, timestep_max, gamma):
    rho = 0
    total_times = np.zeros(timestep_max)
    occur_times = np.zeros(timestep_max)
    for episode in episodes:
        for i in range(len(episode)):
            (s_opt, a_opt, r, s_next) = episode[i]
            total_times[i] += 1
            if s == s_opt and a == a_opt:
                occur_times[i] += 1
    for i in reversed(range(timestep_max)):
        if total_times[i]:
            rho += gamma ** i * occur_times[i] / total_times[i]
    return (1 - gamma) * rho

修改 Main.py 文件，文件代码如下：

from MDP import *
from MRP import compute
from MonteCarlo import MonteCarlo

S = ['S1', 'S2', 'S3', 'S4', 'S5']

A = ['保持S1', '前往S1', '前往S2', '前往S3', '前往S4', '前往S5', '概率前往']

P = {'S1-保持S1-S1': 1.0, 'S1-前往S2-S2': 1.0,
     'S2-前往S1-S1': 1.0, 'S2-前往S3-S3': 1.0,
     'S3-前往S4-S4': 1.0, 'S3-前往S5-S5': 1.0,
     'S4-前往S5-S5': 1.0, 'S4-概率前往-S2': 0.2, 'S4-概率前往-S3': 0.4, 'S4-概率前往-S4': 0.4}

R = {'S1-保持S1': -1, 'S1-前往S2': 0,
     'S2-前往S1': -1, 'S2-前往S3': -2,
     'S3-前往S4': -2, 'S3-前往S5': 0,
     'S4-前往S5': 10, 'S4-概率前往': 1}

gamma = 0.5

MDP = (S, A, P, R, gamma)

# 策略 1 ：随机策略
Pi_1 = {'S1-保持S1': 0.5, 'S1-前往S2': 0.5,
        'S2-前往S1': 0.5, 'S2-前往S3': 0.5,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.5, 'S4-概率前往': 0.5}

# 策略 2 ：给定策略
Pi_2 = {'S1-保持S1': 0.6, 'S1-前往S2': 0.4,
        'S2-前往S1': 0.3, 'S2-前往S3': 0.7,
        'S3-前往S4': 0.5, 'S3-前往S5': 0.5,
        'S4-前往S5': 0.1, 'S4-概率前往': 0.9}


def MDP_to_MRP():
    (S, A, P, R, gamma) = MDP
    R = R_MDP_to_MRP(R, Pi_1, S)
    P = P_MDP_to_MRP(P, Pi_1, S)

    print('使用策略 1，将 MDP 转化为 MRP')
    print('转化后的 MRP 奖励函数：\n', R)
    print('\n转化后的 MRP 状态转移矩阵：\n', P)

    V = compute(P, R, gamma, len(S))
    print("\nMDP 中每个状态价值分别为\n", V)
    print("\n在状态为 S4 时采取动作 概率前往 的价值为：", compute_Q(S[3], A[6], MDP, V))


# MDP_to_MRP(MDP, Pi_1)

def MC():
    timestep_max = 20
    num = 1000
    V = {'S1': 0, 'S2': 0, 'S3': 0, 'S4': 0, 'S5': 0}
    N = {'S1': 0, 'S2': 0, 'S3': 0, 'S4': 0, 'S5': 0}
    episodes = MonteCarlo.sample(MDP, Pi_1, timestep_max, num)
    print("采样前 5 条序列为：\n")
    for i in range(5):
        if i >= len(episodes):
            break
        print("序列 %d：%s" % (i + 1, episodes[i]))
    MonteCarlo.compute(episodes, V, N, gamma)
    print('\n使用蒙特卡洛方法计算 MDP 的状态价值为\n', V)


# MC()


def occupancy_instance():
    gamma = 0.5
    timestep_max = 1000
    num = 1000
    s = S[3]
    a = A[6]
    episodes_1 = MonteCarlo.sample(MDP, Pi_1, timestep_max, num)
    episodes_2 = MonteCarlo.sample(MDP, Pi_2, timestep_max, num)
    rho_1 = occupancy(episodes_1, s, a, timestep_max, gamma)
    rho_2 = occupancy(episodes_2, s, a, timestep_max, gamma)
    print('策略1对状态动作对 (%s, %s) 的占用度量为：%f' % (s, a, rho_1))
    print('策略2对状态动作对 (%s, %s) 的占用度量为：%f' % (s, a, rho_2))


occupancy_instance()

运行 Main.py 文件，结果如下：

D:\RL\MDP\.venv\Scripts\python.exe D:\RL\MDP\Main.py 
策略1对状态动作对 (S4, 概率前往) 的占用度量为：0.109338
策略2对状态动作对 (S4, 概率前往) 的占用度量为：0.226629

进程已结束，退出代码为 0

六、最优策略

  强化学习的目标通常是找到一个策略，使得智能体从初始状态出发能获得最多的期望回报。

• 最优策略 π ∗ ( s ) \pi^*(s) π∗(s)：在 MDP 中，至少存在一个不差于其他所有策略的策略。
  当且仅当对任意的状态 s 都有 V π ( s ) ≥ V π ′ ( s ) V^{\pi}(s)\ge V^{\pi'}(s) Vπ(s)≥Vπ′(s) ，记 π ⪰ π ′ \pi \succeq \pi' π⪰π′ ，MDP 至少存在一个最优策略
• 最优状态价值函数 V ∗ ( s ) V^*(s) V∗(s)： V ∗ ( s ) = max ⁡ π V π ( s ) ,      ∀ s ∈ S V^*(s)=\max\limits_{\pi} V^{\pi}(s),\;\; \forall s\in S V∗(s)=πmax​Vπ(s),∀s∈S
• 最优动作价值函数 Q ∗ ( s , a ) Q^*(s,a) Q∗(s,a)： Q ∗ ( s , a ) = max ⁡ π Q π ( s , a ) ,      ∀ s ∈ S , a ∈ A Q^*(s,a)=\max\limits_{\pi}Q^{\pi}(s,a),\;\;\forall s \in S,a \in A Q∗(s,a)=πmax​Qπ(s,a),∀s∈S,a∈A
• 最优状态价值函数与最优动作价值函数的关系：
  Q ∗ ( s , a ) = r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a ) V ∗ ( s ′ )    V ∗ ( s ) = max ⁡ a ∈ A Q ∗ ( s , a ) Q^*(s,a)=r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)V^*(s') \\ \;\\ V^*(s)=\max\limits_{a\in A}Q^*(s,a) Q∗(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)V∗(s′)V∗(s)=a∈Amax​Q∗(s,a)
• 贝尔曼最优方程：
  Q ∗ ( s , a ) = r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a ) max ⁡ a ′ ∈ A Q ∗ ( s ′ , a ′ )    V ∗ ( s ) = max ⁡ a ∈ A { r ( s , a ) + γ ∑ s ′ ∈ S P ( s ′ ∣ s , a ) V ∗ ( s ′ ) } Q^*(s,a)=r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)\max\limits_{a'\in A}Q^*(s',a') \\ \;\\ V^*(s)=\max\limits_{a\in A}\{r(s,a)+\gamma \sum_{s'\in S}P(s'|s,a)V^*(s')\} Q∗(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)a′∈Amax​Q∗(s′,a′)V∗(s)=a∈Amax​{r(s,a)+γs′∈S∑​P(s′∣s,a)V∗(s′)}
• 得到最优策略的方法：
  • 第4章 动态规划算法
  • 第5章 时序差分算法