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

马尔可夫决策过程始终贯穿强化学习，要学好强化学习，必须掌握马尔可夫决策过程的基础知识。与多臂老虎机不同，马尔可夫决策过程包含状态信息以及状态转移机制。

马尔可夫过程

当且仅当某时刻的状态只取决于上个时刻的状态时，一个随机过程被称为具有马尔可夫性质。
$$P(S_{t+1}|S_t)=P(S_{t+1}|S_1,S_2,...,S_t)$$

而具有马尔可夫性质的随机过程就称为马尔可夫过程，我们通常用一个元组表示 $<S,P>$。$S$ 是状态集合 { s 1 , s 2 , s 3 , . . . } \{s_1, s_2, s_3, ...\} {s1​,s2​,s3​,...}，$P$ 是状态转移矩阵，表明状态之间相互转移的概率。
$$\mathcal{P}=\left[\begin{array}{ccc}P(s_1|s_1)& \cdots & P(s_n|s_1)\\ ⋮& \ddots & ⋮\\ P(s_1|s_n)& \cdots & P(s_n|s_n)\end{array}\right]$$

例如对于上面这个随机过程我们就可以表示为：

S = [0, 1, 2, 3, 4, 5]
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]]
P = np.array(P)

马尔可夫奖励过程

在马尔可夫过程的基础上加入奖励函数 $r$ 和折扣因子 $\gamma$，就可以得到马尔可夫奖励过程，我们通常用 $(\mathcal{S},\mathcal{P},r,\gamma )$ 表示。

在一个马尔可夫奖励过程中， 从 t 时刻状态 $S_t$ 开始，到终止状态时，衰减奖励之和记作回报 $G_t$，即
$$G_t=R_t+\gamma G_{t+1}$$

比如我们选择 $s_5$ 为起点，设置 $\gamma =0.5$，选择的路径是 $s_5\to s_2\to s_3\to s_4\to s_6$。则在这个过程中 $s_5$ 的回报$G_1=1+0.5\ast -2+0.25\ast -2+0.125\ast 10+0.0625\ast 0=0.75$

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]]
P = np.array(P)

reward = [-1, -2, -2, 10, 1, 0]
gamma = 0.5

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

chain = [5, 2, 3, 4, 6]
start_index = 0
G = compute_return(start_index, chain, gamma)
print("根据本序列计算得到回报为：%s。" % G)

在马尔可夫奖励过程中，一个状态的期望回报被称为这个状态的价值。所有状态的价值就组成了价值函数 $V(s)$。
$$V(s)=r(s)+\gamma \sum _{s^{\prime }\in S}p(s^{\prime }|s)V(s^{\prime })$$
其中，$r(s)$ 是即时奖励，$p(s^{\prime }|s)$ 为状态转移概率。

上式就是马尔可夫奖励过程中非常有名的贝尔曼方程。相信大家已经看出来了，这个方程具有解析解。
$$\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}$$

def compute(P, reward, gamma, status_num):
    rewards = np.array(reward).reshape((-1, 1))
    value = np.dot(np.linalg.inv(np.eye(status_num, status_num)-gamma * P), rewards)
    return value

v = compute(P, reward, gamma, 6)
print("MRP中每个状态价值分别为\n", v)

马尔可夫决策过程

马尔可夫过程和马尔可夫奖励过程都是自发改变的随机过程；而如果有一个外界的“刺激”来共同改变这个随机过程，就有了马尔可夫决策过程。

我们将这个来自外界的刺激称为智能体的动作，在马尔可夫奖励过程（MRP）的基础上加入动作，就得到了马尔可夫决策过程（MDP），记作 $(\mathcal{S},\mathcal{A},\mathcal{P},r,\gamma )$。

在马尔可夫决策过程中，通常存在一个智能体来执行动作。智能体根据当前状态选择动作；MDP 根据奖励函数和状态转移函数得到和并反馈给智能体。智能体的目标是最大化得到的累计奖励。智能体根据当前状态从动作的集合中选择一个动作的函数，被称为策略。

智能体的策略通常用 $\pi$ 来表示。$\pi (a|s)$ 表示在输入状态 $s$ 下采取动作 $a$ 的概率。和 MRP 相似，MDP也定义了类似的价值函数。

在 MDP 中基于策略 $\pi$ 的状态价值函数 $V^{\pi }(s)$ 定义为从状态 $s$ 出发遵循策略 $\pi$ 能获得的期望回报。
$$V^{\pi }(s)={\mathbb{E}}_{\pi }\left[G_t\mid S_t=s\right]$$
在 MDP 中，由于动作的存在，我们额外定义一个动作价值函数 $Q\pi (s,a)$。
$$V^{\pi }(s,a)={\mathbb{E}}_{\pi }\left[G_t\mid S_t=s,A_t=a\right]$$

状态的价值等于在该状态下基于策略采取所有动作的概率与相应的价值相乘再求和的结果。
$$V^{\pi }(s)=\sum _{a\in A}\pi (a|s)V^{\pi }(s,a)$$

同时，和MRP类似，状态 $s$ 下采取动作 $a$ 的价值为
$$V^{\pi }(s,a)=r(s,a)+\gamma \sum _{s^{\prime }\in S}p(s^{\prime }|s,a)V^{\pi }(s^{\prime })$$

蒙特卡洛方法

蒙特卡洛方法（Monte-Carlo methods）也被称为统计模拟方法，是一种基于概率统计的数值计算方法。运用蒙特卡洛方法时，我们通常使用重复随机抽样，然后运用概率统计方法来从抽样结果中归纳出我们想求的目标的数值估计。一个简单的例子是用蒙特卡洛方法来计算圆的面积。例如，如图所示的正方形内部随机产生若干个点，细数落在圆中点的个数，圆的面积与正方形面积之比就等于圆中点的个数与正方形中点的个数之比。

对于MDP问题应用以上算法：

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,
}


# 把输入的两个字符串通过“-”连接,便于使用上述定义的P、R变量
def join(str1, str2):
    return str1 + '-' + str2
gamma = 0.5
# 转化后的MRP的状态转移矩阵
P_from_mdp_to_mrp = [
    [0.5, 0.5, 0.0, 0.0, 0.0],
    [0.5, 0.0, 0.5, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.5, 0.5],
    [0.0, 0.1, 0.2, 0.2, 0.5],
    [0.0, 0.0, 0.0, 0.0, 1.0],
]
P_from_mdp_to_mrp = np.array(P_from_mdp_to_mrp)
R_from_mdp_to_mrp = [-0.5, -1.5, -1.0, 5.5, 0]

V = compute(P_from_mdp_to_mrp, R_from_mdp_to_mrp, gamma, 5)
print("MDP中每个状态价值分别为\n", V)
def sample(MDP, Pi, timestep_max, number):
    ''' 采样函数,策略Pi,限制最长时间步timestep_max,总共采样序列数number '''
    S, A, P, R, gamma = MDP
    episodes = []
    for _ in range(number):
        episode = []
        timestep = 0
        s = S[np.random.randint(4)]  # 随机选择一个除s5以外的状态s作为起点
        # 当前状态为终止状态或者时间步太长时,一次采样结束
        while s != "s5" 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
            # 根据状态转移概率得到下一个状态s_next
            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,a,r,s_next）元组放入序列中
            s = s_next  # s_next变成当前状态,开始接下来的循环
        episodes.append(episode)
    return episodes


# 采样5次,每个序列最长不超过20步
episodes = sample(MDP, Pi_1, 20, 5)
print('第一条序列\n', episodes[0])
print('第二条序列\n', episodes[1])
print('第五条序列\n', episodes[4])
# 对所有采样序列计算所有状态的价值
def MC(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] = N[s] + 1
            V[s] = V[s] + (G - V[s]) / N[s]


timestep_max = 20
# 采样1000次,可以自行修改
episodes = sample(MDP, Pi_1, timestep_max, 1000)
gamma = 0.5
V = {"s1": 0, "s2": 0, "s3": 0, "s4": 0, "s5": 0}
N = {"s1": 0, "s2": 0, "s3": 0, "s4": 0, "s5": 0}
MC(episodes, V, N, gamma)
print("使用蒙特卡洛方法计算MDP的状态价值为\n", V)