【强化学习】03 ——马尔可夫决策过程

在此推荐另一篇文章【自动驾驶决策规划】POMDP之Introduction

1. 马尔科夫决策过程(Markov Decision Process，MDP)

• 提供了一套在结果部分随机、部分在决策者的控制下的决策过程建模的数学框架。
• MDP形式化地描述了一种强化学习的环境
  • 环境完全可观测
  • 当前状态可以完全表征过程(马尔科夫性质)
• 几乎所有的RL问题都可以转换到MDP上进行解决
  • 最优控制主要处理连续MDP
  • 部分不可观的问题可以转为MDP
  • 多臂老虎机是单状态的MDP

1.1. 马尔科夫性质

“The future is independent of the past given the present”

未来状态的概率分布只与当前状态有关，而与过去状态无关。

定义：

• 状态$S_t$是Markov的，当且仅当$$\mathbb{P}\left[S_{t+1}\mid S_t\right]=\mathbb{P}\left[S_{t+1}\mid S_1,...,S_t\right]$$

性质：

• 状态捕获历史记录中的所有相关信息
• 一旦知道了这个状态，可以抛弃历史
• 即当前状态是对未来的充分统计

1.2. 状态转移矩阵

${\mathbf{P}}_{s{s}^{\prime }}$为从状态$s$转移到状态$s^{\prime }$的概率，又称一步状态转移概率。$\mathbf{P}$为一步状态转移矩阵。

P [ S t + 1 ∣ S t ] = P [ S t + 1 ∣ S 1 , … , S t ] P s s ′ = P [ S t + 1 = s ′ ∣ S t = s ] P = [ P 11 P 12 … P 1 n P 21 P 22 … P 2 n … … P n 1 P n 2 … P n n ] P[St+1|St]=P[St+1|S1,…,St]Pss′=P[St+1=s′|St=s]P=[P11P12…P1nP21P22…P2n……Pn1Pn2…Pnn]

P[St+1​∣St​]=P[St+1​∣S1​,…,St​]Pss′​=P[St+1​=s′∣St​=s]P= ​P11​P21​……Pn1​​P12​P22​Pn2​​………​P1n​P2n​Pnn​​ ​​

矩阵有以下性质：

1. 非负性性质,$P_{ij}\ge 0$
2. 行元素和为1,$\sum P_{ij}=1,i=1,2,...,n$

以上图为例，$S_1$转移到自身和$S2$的概率分别为0.1，0.9；$S_2$转移到自身和$S1$的概率分别为0.2，0.8。状态转移矩阵可表示为： P = [ 0.1 0.9 0.8 0.2 ] \boldsymbol{P}=[0.10.90.80.2]

P=[0.10.8​0.90.2​]

1.3. 马尔可夫过程

马尔可夫过程（Markov process）指具有马尔可夫性质的随机过程，也被称为马尔可夫链（Markov chain）。马尔可夫链可由二元组$(S,P)$描述。状态转移概率不随时间发生变化。

1.3.1. 一个简单的例子

$Class1$为起始状态，$Sleep$为终止状态，之中的数字代表了各状态之间的转移概率，从每个状态出发转移到其他状态的概率总和为 1。

我们可以写出这个马尔可夫过程的状态转移矩阵：

给定一个马尔可夫过程，我们就可以从某个状态出发，根据它的状态转移矩阵生成一个状态序列（episode），这个步骤也被叫做采样（sampling）。例如，我们可能得到一下采样结果：

• C1 C2 C3 Pass Sleep
• C1 FB FB C1 C2 Sleep
• C1 C2 C3 Pub C2 C3 Pass Sleep
• C1 FB FB C1 C2 C3 Pub C1 FB FB
• FB C1 C2 C3 Pub C2 Sleep

2. 马尔可夫奖励过程

在马尔可夫过程的基础上加入奖励函数$r$ 和折扣因子$\gamma$，就可以得到马尔可夫奖励过程（Markov reward process）。一个马尔可夫奖励过程由$⟨\mathcal{S},\mathcal{P},\mathcal{R},\gamma ⟩$构成，各个组成元素的含义如下所示。

• $S$：状态集合， S = { s 1 , s 2 , . . . , s n } S=\{s_1, s_2, ..., s_n\} S={s1​,s2​,...,sn​}
• $P(s^{\prime }|s)$：状态转移概率，表示在当前状态 $s$ 下，采取动作 $a$ 后转移到状态 $s^{\prime }$ 的概率，$s^{\prime }\in S$，$s,s^{\prime }\in S$。${\mathcal{P}}_{s{s}^{\prime }}=\mathbb{P}\left[S_{t+1}=s^{\prime }\mid S_t=s\right]$
• $R(s,s^{\prime })$：即时奖励函数，表示在当前状态 $s$ 下，采取动作 $a$ 后转移到状态 $s^{\prime }$ 所获得的即时奖励，$s,s^{\prime }\in S$。${\mathcal{R}}_s=\mathbb{E}\left[R_{t+1}\mid S_t=s\right]$
• $\gamma \in (0,1)$：折扣因子，表示当前奖励和未来奖励的权重比例。
  • 如果 $\gamma \to 1$,则代表更相信现在的状态，即“短视myopic”；
  • 如果 $\gamma \to 0$,则代表更相信未来的状态，即“远见far-sighted”；

2.1. 回报

定义：
在一个马尔可夫奖励过程中，从第$t$时刻状态$S_t$开始，直到终止状态时，所有奖励的衰减之和称为回报（Return）$G_t$，公式如下：

$$G_t=R_{t+1}+\gamma R_{t+2}+...=\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}$$

仍然沿用上面的例子：
例如，进入状态$Class2$可以得到奖励$-2$，表明我们不希望进入，进入$Pass$可以获得最高的奖励$10$，但是进入$Sleep$之后奖励为零，并且此时序列也终止了。一个计算例子：$Class1\to Class2\to Class3\to Pass\to Sleep$：$-2+0.5\times (-2)+(0.5{)}^2\times (-2)+(0.5{)}^3\times 10=-2.25$

用代码进行表示：

import numpy as np

# Define the transition Matrix
P = [
    [0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0],
    [0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.2],
    [0.0, 0.0, 0.0, 0.6, 0.4, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
    [0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0],
    [0.1, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
]
P = np.array(P)
RewardVector = [-2, -2, -2, 10, 1, -1, 0]

# 给定一条序列,计算从某个索引（起始状态）开始到序列最后（终止状态）得到的回报
def ComputeSequenceReward(Start_idx, Sequence, RewardVector, gamma=0.5):
    TotalReward = 0.0
    for i in reversed(range(Start_idx, len(Sequence))):
        TotalReward = gamma * TotalReward + RewardVector[Sequence[i] - 1]
    return TotalReward

def test01():
    chain = [1, 2, 3, 4, 7]
    start_index = 0
    print("根据本序列计算得到回报为：%s。"% ComputeSequenceReward(start_index, chain, RewardVector, gamma=0.5))
    
if __name__ == "__main__":
    test01()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

根据本序列计算得到回报为：-2.25。
1

$\gamma =0，0.9，1$时的例子(图中数字为各节点的价值)：

2.2. 价值函数

在马尔可夫奖励过程中，一个状态的期望回报（即从这个状态出发的未来累积奖励的期望）被称为这个状态的价值（value）。所有状态的价值就组成了价值函数（value function），价值函数的输入为某个状态，输出为这个状态的价值。我们将价值函数写成$v(s)=\mathbb{E}\left[G_t\mid S_t=s\right]$，展开为 V ( s ) = E [ G t ∣ S t = s ] = E [ R t + 1 + γ R t + 2 + γ 2 R t + 3 + . . . ∣ S t = s ] = E [ R t + 1 + γ ( R t + 2 + γ R t + 3 + . . . ) ∣ S t = s ] = E [ R t + 1 + γ G t + 1 ∣ S t = s ] = E [ R t + 1 + γ V ( S t + 1 ) ∣ S t = s ] V(s)=E[Gt∣St=s]=E[Rt+1+γRt+2+γ2Rt+3+...∣St=s]=E[Rt+1+γ(Rt+2+γRt+3+...)∣St=s]=E[Rt+1+γGt+1∣St=s]=E[Rt+1+γV(St+1)∣St=s]

V(s)​=E[Gt​∣St​=s]=E[Rt+1​+γRt+2​+γ2Rt+3​+...∣St​=s]=E[Rt+1​+γ(Rt+2​+γRt+3​+...)∣St​=s]=E[Rt+1​+γGt+1​∣St​=s]=E[Rt+1​+γV(St+1​)∣St​=s]​

在上式的最后一个等号中,一方面,即时奖励的期望正是奖励函数的输出,即$\mathbb{E}[R_{t+1}|S_t=s]={\mathcal{R}}_s]$;另一
方面,等式中剩余部分$\mathbb{E}[\gamma V(S_{t+1})|S_t=s]$可以根据从状态$s$出发的转移概率得到,即可以得到$$V(s)={\mathcal{R}}_s+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}V(s^{\prime })$$

上式就是马尔可夫奖励过程中非常有名的贝尔曼方程（Bellman equation），对每一个状态都成立。我们可以将其写成矩阵形式：
$$\mathcal{V}=\mathcal{R}+\gamma \mathcal{P}\mathcal{V}$$ [ V ( s 1 ) ⋮ V ( s n ) ] = [ R 1 ⋮ R n ] + γ [ P 11 … P 1 n ⋮ P 11 … P n n ] [ V ( s 1 ) ⋮ V ( s n ) ] [V(s1)⋮V(sn)]

=[R1⋮Rn]+\gamma[P11…P1n⋮P11…Pnn][V(s1)⋮V(sn)] ​V(s1​)⋮V(sn​)​ ​= ​R1​⋮Rn​​ ​+γ ​P11​⋮P11​​……​P1n​Pnn​​ ​ ​V(s1​)⋮V(sn​)​ ​
我们可以直接根据矩阵运算求解，得到以下解析解：$$\begin{array}{rl}\mathcal{V}& =\mathcal{R}+\gamma \mathcal{P}\mathcal{V}\\ (I-\gamma \mathcal{P})\mathcal{V}& =\mathcal{R}\\ \mathcal{V}& =(I-\gamma \mathcal{P}{)}^{-1}\mathcal{R}\end{array}$$

以上解析解的计算复杂度是$O(n^3)$,其中$n$是状态个数,因此这种方法只适用很小的马尔可夫奖励过程。求解较大规模的马尔可夫奖励过程中的价值函数时，可以使用动态规划（dynamic programming）算法、蒙特卡洛方法（Monte-Carlo method）(见3.6)和时序差分（temporal difference）

接下来编写代码来实现求解价值函数的解析解方法，并据此计算该马尔可夫奖励过程中所有状态的价值。

# Exploit Bellman equation to compute value of all states
def ComputeValue(RewardVector, Statesize, TransitionMatrix=P, gamma=0.5):
    RewardVector = np.array(RewardVector).reshape(-1, 1)
    Value = np.dot(np.linalg.inv(np.eye(Statesize, Statesize) - gamma * TransitionMatrix),
                   RewardVector)
    return Value

print("MRP中每个状态价值分别为\n", ComputeValue(RewardVector, 7))
1
2
3
4
5
6
7
8

MRP中每个状态价值分别为
 [[-2.90815722]
 [-1.55006913]
 [ 1.12482718]
 [10.        ]
 [ 0.62413589]
 [-2.08255975]
 [ 0.        ]]
1
2
3
4
5
6
7
8

3. 马尔科夫决策过程

3.1. MDP五元组

MDP（Markov Decision Process）是用于描述智能体与环境交互的数学模型。MDP 可以表示为五元组 $(S,A,P,R,\gamma )$：

• $S$：状态集合， S = { s 1 , s 2 , . . . , s n } S=\{s_1, s_2, ..., s_n\} S={s1​,s2​,...,sn​}，包含了车道、环境、世界模型等信息。
• $A$：动作集合， A = { a 1 , a 2 , . . . , a m } A=\{a_1, a_2, ..., a_m\} A={a1​,a2​,...,am​}，车辆的决策空间，包含换道、跟随、超车等等。
• $P(s^{\prime }|s,a)$：状态转移概率函数，表示在当前状态 $s$ 下，采取动作 $a$ 后转移到状态 $s^{\prime }$ 的概率，$s^{\prime }\in S$，$s,s^{\prime }\in S$，$a\in A$。${\mathcal{P}}_{{\mathbf{ss}}^{\prime }}^a=\mathbb{P}\left[S_{t+1}=s^{\prime }\mid S_t=s,A_t=a\right]$
• $R(s,a,s^{\prime })$：即时奖励函数，表示在当前状态 $s$ 下，采取动作 $a$ 后转移到状态 $s^{\prime }$ 所获得的即时奖励，$s,s^{\prime }\in S$，$a\in A$。${\mathcal{R}}_{\mathfrak{s}}^a=\mathbb{E}\left[R_{t+1}\mid S_t=s,A_t=a\right]$
• $\gamma \in (0,1)$：折扣因子，表示当前奖励和未来奖励的权重比例。

注意: 在上面 MDP 的定义中，我们不再使用类似 MRP 定义中的状态转移矩阵方式，而是直接表示成了状态转移函数。

• 这样做一是因为此时状态转移与动作也有关，变成了一个三维数组，而不再是一个矩阵（二维数组）；
• 二是因为状态转移函数更具有一般意义，例如，如果状态集合不是有限的，就无法用数组表示，但仍然可以用状态转移函数表示。

不同于马尔可夫奖励过程，在马尔可夫决策过程中，通常存在一个智能体来执行动作。马尔可夫决策过程是一个与时间相关的不断进行的过程，在智能体和环境 之间存在一个不断交互的过程。一般而言，它们之间的交互是如下图的循环过程：智能体根据当前状态$S_t$选择动作$A_t$；对于状态$S_t$和动作$A_t$，环境根据奖励函数和状态转移函数得到$S_{t+1}$和$R_{t+1}$并反馈给智能体。智能体的目标是最大化得到的累计奖励。智能体根据当前状态从动作的集合中选择一个动作的函数，被称为策略（policy）。

3.2. 策略

智能体的策略（Policy）通常用字母$\pi$表示。在输入状态$s$情况下采取动作$a$的概率的策略如下式所示。$$\pi (a|s)=\mathbb{P}\left[A_t=a\mid S_t=s\right]$$

• 策略只需要与当前状态有关，不需要考虑历史状态。
• 策略是固定的（与时间无关）$A_t\sim \pi (\cdot |S_t),\forall t>0$

3.3. 价值函数

3.3.1. 状态价值函数

我们用$V^{\pi }(s)$表示在 MDP 中基于策略$\pi$的状态价值函数（state-value function），定义为从状态$s$出发遵循策略$\pi$能获得的期望回报，数学表达为：
$$V^{\pi }(s)={\mathbb{E}}_{\pi }[G_t|S_t=s]$$

3.3.2. 动作价值函数

不同于 MRP，在 MDP 中，由于动作的存在，我们额外定义一个动作价值函数（action-value function）。我们用$Q^{\pi }(s,a)$表示在 MDP 遵循策略$\pi$时，对当前状态$s$执行动作$a$得到的期望回报：
$$Q^{\pi }(s,a)={\mathbb{E}}_{\pi }[G_t|S_t=s,A_t=a]$$

3.4. 贝尔曼期望方程

在贝尔曼方程中加上“期望”二字是为了与接下来的贝尔曼最优方程进行区分。

状态价值函数可以分解为当前时刻的奖励加上后继状态的经过折扣后的状态价值
$$V_{\pi }(s)={\mathbb{E}}_{\pi }\left[R_{t+1}+\gamma V_{\pi }(S_{t+1})\mid S_t=s\right]$$

动作价值函数同样也可以分解：
$$Q_{\pi }(s,a)={\mathbb{E}}_{\pi }\left[R_{t+1}+\gamma Q_{\pi }(S_{t+1},A_{t+1})\mid S_t=s,A_t=a\right]$$

状态价值函数和动作价值函数之间的关系：在使用策略$\pi$中，状态$s$的价值等于在该状态下基于策略$\pi$采取所有动作$a_i$的概率与相应的价值相乘再求和的结果：
$$V_{\pi }(s)=\sum _{a\in \mathcal{A}}\pi (a|s)Q_{\pi }(s,a)$$

使用策略$\pi$时，状态$s$下采取动作$a$的价值等于即时奖励加上经过衰减后的所有可能的下一个状态的状态转移概率与相应的价值的乘积：

$$Q_{\pi }(s,a)={\mathcal{R}}_s^a+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}^a{V}_{\pi }(s^{\prime })$$

将两者叠加，于是可以得到

$$V_{\pi }(s)=\sum _{a\in \mathcal{A}}\pi (a|s)\left({\mathcal{R}}_s^a+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}^a{V}_{\pi }(s^{\prime })\right)$$

以及

$$Q_{\pi }(s,a)={\mathcal{R}}_s^a+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}^a\sum _{a^{\prime }\in \mathcal{A}}\pi (a^{\prime }|s^{\prime })Q_{\pi }(s^{\prime },a^{\prime })$$
以上就是通过简单推导分别得到两个价值函数的贝尔曼期望方程。

同样用一个例子进行解释：图中基于策略$\pi$采取所有动作$a_i$的概率为：$\pi (a|s)=0.5$。计算$Class3$的状态价值需要计算其下一步的状态价值之和，因此可以得到图中式子所示：$7.4=0.5\ast (1+0.2\ast -1.3+0.4\ast 2.7+0.4\ast 7.4)+0.5\ast 10$

同样可以利用MRP进行求解。我们可以将策略的动作选择进行边缘化（marginalization)，就可以得到没有动作的 MRP 了。具体来说，对于某一个状态，我们根据策略所有动作的概率进行加权，得到的奖励和就可以认为是一个 MRP 在该状态下的奖励，即：$${\mathcal{R}}^{\pi }={\mathcal{R}}_s^{{}^{\prime }}=\sum _{a\in \mathcal{A}}\pi (a|s){\mathcal{R}}_s^a$$
同理，我们计算采取动作的概率与使$s$转移到$s^{\prime }$的概率的乘积，再将这些乘积相加，其和就是一个 MRP 的状态$s$从转移至$s^{\prime }$的概率：$${\mathcal{P}}^{\pi }={\mathcal{P}}_{s{s}^{\prime }}^{{}^{\prime }}=\sum _{a\in \mathcal{A}}\pi (a|s){\mathcal{P}}_{s{s}^{\prime }}^a$$

如此，便可以继续利用之前MRP的步骤了：$$\begin{gathered} V_{\pi }={\mathcal{R}}^{\pi }+\gamma {\mathcal{P}}^{\pi }{V}_{\pi } \\ {V}_{\pi }=(1-\gamma {\mathcal{P}}^{\pi }{)}^{-1}{\mathcal{R}}^{\pi } \end{gathered}$$

def test02():
    # Define the transition Matrix
    # C1 C2  Pass FB Sleep
    P_TransformMDP2MRP = [
        [0.0, 0.5, 0.0, 0.5, 0.0],
        [0.0, 0.0, 0.5, 0.0, 0.5],
        [0.1, 0.2, 0.2, 0.0, 0.5],
        [0.5, 0.0, 0.0, 0.5, 0.0],
        [0.0, 0.0, 0.0, 0.0, 1.0]
    ]
    P_TransformMDP2MRP = np.array(P_TransformMDP2MRP)
    R_TransformMDP2MRP = [-1.5, -1, 5.5, -0.5, 0.0]
    print("MDP中每个状态价值分别为\n", ComputeValue(R_TransformMDP2MRP, 5, P_TransformMDP2MRP, gamma=1))
1
2
3
4
5
6
7
8
9
10
11
12
13

MDP中每个状态价值分别为
 [[-1.30769231]
 [ 2.69230769]
 [ 7.38461538]
 [-2.30769231]
 [ 0.        ]]
1
2
3
4
5
6

3.5. 最优策略

强化学习的目标通常是找到一个策略，使得智能体从初始状态出发能获得最多的期望回报。我们首先定义策略之间的偏序关系：当且仅当对于任意的状态$s$都有$V^{\pi }(s)\ge V^{{\pi }^{\prime }}(s)$ ，记$\pi >{\pi }^{{}^{\prime }}$。于是在有限状态和动作集合的 MDP 中，至少存在一个策略比其他所有策略都好或者至少存在一个策略不差于其他所有策略，这个策略就是最优策略（optimal policy）。最优策略可能有很多个，我们都将其表示为${\pi }^{\ast }(s)$。

最优策略都有相同的状态价值函数，我们称之为最优状态价值函数，表示为：
$$V^{\ast }(s)=\underset{\pi }{\max }{V}^{\pi }(s),\forall s\in \mathcal{S}$$

同理，我们定义最优动作价值函数:
$$Q^{\ast }(s,a)=\underset{\pi }{\max }{Q}^{\pi }(s,a),\forall s\in \mathcal{S},a\in \mathcal{A}$$

最优状态价值的例子

最优动作价值的例子

通过最大化$Q^{\ast }(s,a)$,我们可以找到最优策略：
π ∗ ( a ∣ s ) = { 1   i f   a = a r g m a x   q ∗ ( s , a ) 0 o t h e r w i s e \left.\pi_*(a|s)=\left\{1 if a=argmax q∗(s,a)0otherwise

\right.\right. π∗​(a∣s)={10​ if a=argmax q∗​(s,a)otherwise​

• 对于所有的MDP，都有一个确定的最优策略
• 只要找到$Q^{\ast }(s,a)$，就能得到最优策略

下图红色弧线箭头部分代表了最优策略的选择

3.5.1. 贝尔曼最优方程

同样通过递归的方式推导方程：

最优状态价值是选择此时使最优动作价值最大的那一个动作时的状态价值：

$$V_{\ast }(s)=\underset{a}{\max }{Q}_{\ast }(s,a)$$

最优状态价值函数和最优动作价值函数之间的关系：

$$Q_{\ast }(s,a)={\mathcal{R}}_s^a+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}^a{V}_{\ast }(s^{\prime })$$
之后叠加，可得

V ∗ ( s ) = max ⁡ a ∈ A { R s a + γ ∑ s ′ ∈ S P s s ′ a V ∗ ( s ′ ) } V_*(s)=\max_{a\in\mathcal{A}}\{\mathcal{R}_s^a+\gamma\sum_{s^{\prime}\in\mathcal{S}}\mathcal{P}_{ss^{\prime}}^aV_*(s^{\prime})\} V∗​(s)=a∈Amax​{Rsa​+γs′∈S∑​Pss′a​V∗​(s′)}

$$Q_{\ast }(s,a)={\mathcal{R}}_s^a+\gamma \sum _{s^{\prime }\in \mathcal{S}}{\mathcal{P}}_{s{s}^{\prime }}^a\underset{a^{\prime }\in \mathcal{A}}{\max }{Q}_{\ast }(s^{\prime },a^{\prime })$$
如此，我们便得到了贝尔曼最优方程（Bellman optimality equation）。

3.5.2. 最优策略求解

• 贝尔曼最优方程是非线性的
• 通常来说，没有闭式解
• 一些迭代方法：
  • Value Iteration
  • Policy Iteration
  • Q-learning
  • Sarsa

3.5.3. 扩展MDP

• Infinite and continuous MDPs
• Partially observable MDPs(POMDP)
• Undiscounted, average reward MDPs（Ergodic Markov Process）

3.6 蒙特卡洛方法（Monte-Carlo method）

求解较大规模的马尔可夫奖励过程中的价值函数时，可以使用动态规划（dynamic programming）算法、蒙特卡洛方法（Monte-Carlo method）和时序差分（temporal difference）。本小节介绍蒙特卡洛方法（Monte-Carlo method）。

蒙特卡洛方法（Monte-Carlo methods）也被称为统计模拟方法，是一种基于概率统计的数值计算方法。运用蒙特卡洛方法时，我们通常使用重复随机抽样，然后运用概率统计方法来从抽样结果中归纳出我们想求的目标的数值估计。一个简单的例子是用蒙特卡洛方法来计算圆的面积。例如，在下图所示的正方形内部随机产生若干个点，细数落在圆中点的个数，圆的面积与正方形面积之比就等于圆中点的个数与正方形中点的个数之比。如果我们随机产生的点的个数越多，计算得到圆的面积就越接近于真实的圆的面积。
我们现在介绍如何用蒙特卡洛方法来估计一个策略在一个马尔可夫决策过程中的状态价值函数。回忆一下，一个状态的价值是它的期望回报，那么一个很直观的想法就是用策略在 MDP 上采样很多条序列，计算从这个状态出发的回报再求其期望就可以了，公式如下：$$V^{\pi }(s)={\mathbb{E}}_{\pi }[G_t|S_t=s]\approx \frac{1}{N}\sum _{i=1}^N{G}_t^{(i)}$$

在一条序列中，可能没有出现过这个状态，可能只出现过一次这个状态，也可能出现过很多次这个状态。我们介绍的蒙特卡洛价值估计方法会在该状态每一次出现时计算它的回报。还有一种选择是一条序列只计算一次回报，也就是这条序列第一次出现该状态时计算后面的累积奖励，而后面再次出现该状态时，该状态就被忽略了。

假设我们现在用策略$\pi$从状态$s$开始采样序列，据此来计算状态价值。我们为每一个状态维护一个计数器和总回报，计算状态价值的具体过程如下所示。

1. 使用策略$\pi$采样若干条序列：$$s_0^{(i)}\stackrel{a_0^{(i)}}{\to }{r}_0^{(i)},s_1^{(i)}\stackrel{a_1^{(i)}}{\to }{r}_1^{(i)},s_2^{(i)}\stackrel{a_2^{(i)}}{\to }\cdots \stackrel{a_{T-1}^{(i)}}{\to }{r}_{T-1}^{(i)},s_T^{(i)}$$
2. 对每一条序列中的每一时间步$t$的状态$s$进行以下操作：

• 更新状态$s$的计数器$N(s)=N(s)+1$
• 更新状态$s$的总回报$M(s)=M(s)+G_t$;
• 或者采用增量式更新策略：$V(s)\leftarrow V(s)+\frac{1}{N(s)}(G-V(S))$

3. 每一个状态的价值被估计为回报的平均值 $V(s)=M(s)/N(s)$。根据大数定律，当$N(s)\to \infty$，有$V(s)\to V_{\pi }(s)$ 。

代码采样部分结果

[('C2', 'Study', -2, 'Pass'), ('Pass', 'Pub', 1, 'Pass'), ('Pass', 'Pub', 1, 'Pass'), ('Pass', 'Pub', 1, 'Pass'), ('Pass', 'Study', 10, 'Sleep')]
[('C1', 'Study', -2, 'C2'), ('C2', 'Sleep', 0, 'Sleep')]
[('FB', 'Quit', 0, 'C1'), ('C1', 'Facebook', -1, 'FB'), ('FB', 'Quit', 0, 'C1'), ('C1', 'Facebook', -1, 'FB'), ('FB', 'Quit', 0, 'C1'), ('C1', 'Study', -2, 'C2'), ('C2', 'Study', -2, 'Pass'), ('Pass', 'Pub', 1, 'Pass'), ('Pass', 'Pub', 1, 'Pass')]
[('C1', 'Study', -2, 'C2'), ('C2', 'Sleep', 0, 'Sleep')]
[('C2', 'Study', -2, 'Pass'), ('Pass', 'Pub', 1, 'C2'), ('C2', 'Sleep', 0, 'Sleep')]
1
2
3
4
5

可以看到最后结果还是比较接近的。

使用蒙特卡洛方法计算MDP的状态价值为
 {'C1': -1.6584167352261565, 'C2': 0.5744913689985154, 'Pass': 6.330419227770518, 'FB': -1.1820907116805823, 'Sleep': 0}
MDP中每个状态价值分别为
 [[-1.67666232]
 [ 0.51890482]
 [ 6.0756193 ]
 [-1.22555411]
 [ 0.        ]]
1
2
3
4
5
6
7
8

3.7. 占有度量

不同策略的价值函数是不一样的。这是因为对于同一个 MDP，不同策略会访问到的状态的概率分布是不同的。我们需要理解不同策略会使智能体访问到不同概率分布的状态这个事实，这会影响到策略的价值函数。

首先我们定义 MDP 的初始状态分布为${\nu }_0(s)$，在有些资料中，初始状态分布会被定义进 MDP 的组成元素中。我们用$P_t^{\pi }(s)$表示采取策略$\pi$使得智能体在时刻$t$状态为$s$的概率，所以我们有$P_0^{\pi }(s)={\nu }_0(s)$，然后就可以定义一个策略的状态访问分布（state visitation distribution）：

$${\nu }^{\pi }(s)=(1-\gamma )\sum _{t=0}^{\infty }{\gamma }^t{P}_t^{\pi }(s)$$
其中，$1-\gamma$是用来使得概率加和为 1 的归一化因子。状态访问概率表示一个策略和 MDP 交互会访问到的状态的分布。需要注意的是，理论上在计算该分布时需要交互到无穷步之后，但实际上智能体和 MDP 的交互在一个序列中是有限的。不过我们仍然可以用以上公式来表达状态访问概率的思想，状态访问概率有如下性质：$${\nu }^{\pi }(s^{\prime })=(1-\gamma ){\nu }_0(s^{\prime })+\gamma \int P(s^{\prime }|s,a)\pi (a|s){\nu }^{\pi }(s)dsda$$
此外，我们还可以定义策略的占用度量（occupancy measure）：$${\rho }^{\pi }(s,a)=(1-\gamma )\sum _{t=0}^{\infty }{\gamma }^t{P}_t^{\pi }(s)\pi (a|s)$$
它表示动作状态对$(s,a)$被访问到的概率。二者之间存在如下关系：$${\rho }^{\pi }(s,a)={\nu }^{\pi }(s)\pi (a|s)$$
进一步得出如下两个定理。定理 1：智能体分别以策略${\pi }_1$和${\pi }_2$和同一个 MDP 交互得到的占用度量${\rho }^{{\pi }_1}$和${\rho }^{{\pi }_2}$满足$${\rho }^{{\pi }_1}={\rho }^{{\pi }_2}⟺{\pi }_1={\pi }_2$$
定理 2：给定一合法占用度量，可生成该占用度量的唯一策略是$${\pi }_{\rho }=\frac{\rho (s,a)}{{\sum }_{a^{\prime }}\rho (s,a^{\prime })}$$
注意：以上提到的“合法”占用度量是指存在一个策略使智能体与 MDP 交互产生的状态动作对被访问到的概率。

# Occupancy
def test04():
    # 策略2
    Policy_2 = {
        "C1-Study": 0.6,
        "C1-Facebook": 0.4,
        "FB-Facebook": 0.3,
        "FB-Quit": 0.7,
        "C2-Study": 0.5,
        "C2-Sleep": 0.5,
        "Pass-Study": 0.1,
        "Pass-Pub": 0.9,
    }
    MAXTimeStep = 8
    MDP, Policy_1 = Set_MDPParameterAndPolicy()
    Sequences1 = MonteCarloSampling(MDP, Policy_1, MAXTimeStep, SamplingNum=1000)
    Sequences2 = MonteCarloSampling(MDP, Policy_2, MAXTimeStep, SamplingNum=1000)
    # Sequences1 = sample(MDP, Policy_1, MAXTimeStep, 1000)
    # Sequences2 = sample(MDP, Policy_2, MAXTimeStep, 1000)
    rho1 = ComputeOccupancy("Pass", "Pub", Sequences1, MAXTimeStep, MDP)
    rho2 = ComputeOccupancy("Pass", "Pub", Sequences2, MAXTimeStep, MDP)
    print(rho1, rho2)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

最后得到的占有度量有所不同

0.058 0.1145
1

代码

import numpy as np

# 给定一条序列,计算从某个索引（起始状态）开始到序列最后（终止状态）得到的回报
def ComputeSequenceReward(Start_idx, Sequence, RewardVector, gamma=0.5):
    TotalReward = 0.0
    for i in reversed(range(Start_idx, len(Sequence))):
        TotalReward = gamma * TotalReward + RewardVector[Sequence[i] - 1]
    return TotalReward

# Exploit Bellman equation to compute value of all states
def ComputeValue(RewardVector, Statesize, TransitionMatrix, gamma=0.5):
    RewardVector = np.array(RewardVector).reshape(-1, 1)
    try:
        Value = np.dot(np.linalg.inv(np.eye(Statesize, Statesize) - gamma * TransitionMatrix),
                   RewardVector)
    except:
        print("-------------状态转移矩阵为奇异矩阵，存在求解误差-------------")
        TransitionMatrix[Statesize - 1][Statesize - 1] += 1e-7
        I = np.eye(Statesize, Statesize)
        Value = np.dot(np.linalg.inv(I - gamma * TransitionMatrix),
                   RewardVector)
    return Value

def Set_MDPParameterAndPolicy():
    # 状态集合
    S = ["C1", "C2", "Pass", "FB", "Sleep"]
    # 动作集合
    A = ["Facebook", "Study", "Sleep", "Pub", "Quit"]
    # 状态转移函数
    P = {
        "C1-Study-C2": 1.0,
        "C1-Facebook-FB": 1.0,
        "FB-Facebook-FB": 1.0,
        "FB-Quit-C1": 1.0,
        "C2-Study-Pass": 1.0,
        "C2-Sleep-Sleep": 1.0,
        "Pass-Study-Sleep": 1.0,
        "Pass-Pub-C1": 0.2,
        "Pass-Pub-C2": 0.4,
        "Pass-Pub-Pass": 0.4,
    }
    # 奖励函数
    R = {
        "C1-Study": -2,
        "C1-Facebook": -1,
        "FB-Facebook": -1,
        "FB-Quit": 0,
        "C2-Study": -2,
        "C2-Sleep": 0,
        "Pass-Study": 10,
        "Pass-Pub": 1,
    }
    # 折扣因子
    gamma = 0.5
    MDP = (S, A, P, R, gamma)

    # 策略1,随机策略
    Pi_1 = {
        "C1-Study": 0.5,
        "C1-Facebook": 0.5,
        "FB-Facebook": 0.5,
        "FB-Quit": 0.5,
        "C2-Study": 0.5,
        "C2-Sleep": 0.5,
        "Pass-Study": 0.5,
        "Pass-Pub": 0.5,
    }
    # 策略2
    Pi_2 = {
        "C1-Study": 0.7,
        "C1-Facebook": 0.3,
        "FB-Facebook": 0.3,
        "FB-Quit": 0.7,
        "C2-Study": 0.5,
        "C2-Sleep": 0.5,
        "Pass-Study": 0.2,
        "Pass-Pub": 0.8,
    }
    return MDP, Pi_1

# 把输入的两个字符串通过“-”连接,便于使用上述定义的P、R变量
def join(str1, str2):
    return str1 + '-' + str2

def MonteCarloSampling(MDP, Policy, MAXTimeStep, SamplingNum):
    ''' 采样函数,策略Pi,限制最长时间步MaxTimeStep,总共采样序列数SamplingNum '''
    S, A, P, R, gamma = MDP
    StateNum = len(S)
    Sequences = []
    for _ in range(SamplingNum):
        Sequence = []
        TimeStep = 0
        # 随机选择一个除Sleep以外的状态s作为起点
        s = S[np.random.randint(StateNum - 1)]
        # 当前状态为终止状态或者时间步太长时,一次采样结束
        while s != "Sleep" and TimeStep <= MAXTimeStep:
            TimeStep += 1
            rand, temp = np.random.rand(), 0
            # 在状态s下根据策略选择动作
            for a_ in A:
                temp += Policy.get(join(s, a_), 0.0)
                if temp >= rand:
                    a = a_
                    r = R.get(join(s, a_), 0.0)
                    break
            rand, temp = np.random.rand(), 0
            # 根据状态转移概率得到下一个状态s_next
            for s_ in S:
                temp += P.get(join(join(s, a), s_), 0.0)
                if temp >= rand:
                    s_next = s_
                    break
            # 把（s,a,r,s_next）元组放入序列中
            Sequence.append((s, a, r, s_next))
            # s_next变成当前状态,开始接下来的循环
            s = s_next
        Sequences.append(Sequence)
    return Sequences

# 对所有采样序列计算所有状态的价值
def MonteCarloComputeValue(Sequences, MDP):
    gamma = MDP[4]
    V = {"C1": 0, "C2": 0, "Pass": 0, "FB": 0, "Sleep": 0}
    N = {"C1": 0, "C2": 0, "Pass": 0, "FB": 0, "Sleep": 0}
    for Sequence in Sequences:
        G = 0
        # 一个序列从后往前计算
        for i in reversed(range(len(Sequence))):
            s, r = Sequence[i][0], Sequence[i][2]
            G = r + gamma * G
            N[s] = N[s] + 1
            V[s] = V[s] + (G - V[s]) / N[s]
    return V

def ComputeOccupancy(s, a, Sequences, MAXTimeStep, MDP):
    ''' 计算状态动作对（s,a）出现的频率,以此来估算策略的占用度量 '''
    gamma = MDP[4]
    rho = 0
    total_times = np.zeros(MAXTimeStep)  # 记录每个时间步t各被经历过几次
    occur_times = np.zeros(MAXTimeStep)  # 记录(s_t,a_t)=(s,a)的次数
    for Sequence in Sequences:
        for i in range(len(Sequence)):
            try:
                s_, a_ = Sequence[i][0], Sequence[i][1]
                total_times[i] += 1
                if s_ == s and a_ == a:
                    occur_times[i] += 1
            except IndexError:
                continue
    for i in reversed(range(MAXTimeStep)):
        if total_times[i]:
            # 用频率来估算策略的占用度量
            rho = gamma ** i * occur_times[i] / total_times[i]
    return (1 - gamma) * rho

def SampleTEXT():
    MDP, Policy = Set_MDPParameterAndPolicy()
    Sequences = MonteCarloSampling(MDP, Policy, MAXTimeStep=8, SamplingNum=5)
    for Sequence in Sequences:
        print(Sequence)

def MonteCarloTEXT():
    MDP, Policy = Set_MDPParameterAndPolicy()
    Sequences = MonteCarloSampling(MDP, Policy, MAXTimeStep=8, SamplingNum=5000)
    V = MonteCarloComputeValue(Sequences, MDP)
    print("使用蒙特卡洛方法计算MDP的状态价值为\n", V)

def test01():
    # Define the transition Matrix
    # C1 C2 C3 Pass Pub FB Sleep
    P = [
        [0.0, 0.5, 0.0, 0.0, 0.0, 0.5, 0.0],
        [0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.2],
        [0.0, 0.0, 0.0, 0.6, 0.4, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
        [0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0],
        [0.1, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
    ]
    P = np.array(P)
    RewardVector = [-2, -2, -2, 10, 1, -1, 0]
    chain = [1, 6, 6, 1, 2, 7]
    start_index = 0
    print("根据本序列计算得到回报为：%s。"% ComputeSequenceReward(start_index, chain, RewardVector, gamma=0.5))
    print("MRP中每个状态价值分别为\n", ComputeValue(RewardVector, 7, P))

# MDP2MRP
def test02():
    # Define the transition Matrix
    # C1 C2  Pass FB Sleep
    P_TransformMDP2MRP = [
        [0.0, 0.5, 0.0, 0.5, 0.0],
        [0.0, 0.0, 0.5, 0.0, 0.5],
        [0.1, 0.2, 0.2, 0.0, 0.5],
        [0.5, 0.0, 0.0, 0.5, 0.0],
        [0.0, 0.0, 0.0, 0.0, 1.0]
    ]
    P_TransformMDP2MRP = np.array(P_TransformMDP2MRP)
    R_TransformMDP2MRP = [-1.5, -1, 5.5, -0.5, 0.0]
    print("MDP中每个状态价值分别为\n", ComputeValue(R_TransformMDP2MRP, 5, P_TransformMDP2MRP, gamma=0.5))

# MonteCarlo
def test03():
    # SampleTEXT(
    MonteCarloTEXT()
    test02()

# Occupancy
def test04():
    # 策略2
    Policy_2 = {
        "C1-Study": 0.6,
        "C1-Facebook": 0.4,
        "FB-Facebook": 0.3,
        "FB-Quit": 0.7,
        "C2-Study": 0.5,
        "C2-Sleep": 0.5,
        "Pass-Study": 0.1,
        "Pass-Pub": 0.9,
    }
    MAXTimeStep = 8
    MDP, Policy_1 = Set_MDPParameterAndPolicy()
    Sequences1 = MonteCarloSampling(MDP, Policy_1, MAXTimeStep, SamplingNum=1000)
    Sequences2 = MonteCarloSampling(MDP, Policy_2, MAXTimeStep, SamplingNum=1000)
    rho1 = ComputeOccupancy("Pass", "Pub", Sequences1, MAXTimeStep, MDP)
    rho2 = ComputeOccupancy("Pass", "Pub", Sequences2, MAXTimeStep, MDP)
    print(rho1, rho2)

if __name__ == "__main__":
    test04()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

参考

[1] 伯禹AI
[2] https://www.deepmind.com/learning-resources/introduction-to-reinforcement-learning-with-david-silver
[3] 动手学强化学习
[4] Reinforcement Learning