强化学习：计算马尔可夫奖励过程状态价值函数代码示例

参考书目：王琦，杨毅远，江季，Easy RL：强化学习教程，人民邮电出版社，https://github.com/datawhalechina/easy-rl, 2022.

状态价值函数的计算是强化学习中很重要的一步。在马尔科夫过程中，并不是哪一个状态得到的奖励越大，这个状态就很有价值，相比之下有些状态虽然没有获得环境的奖励，但该状态却又很高的状态价值，我们来看下面这个例子。

根据参考书中第二章马尔可夫过程的例子（如下图所示），我们给出三种方法计算该过程的状态价值函数V(S)及其代码示例。

首先我们在python中定义它的奖励向量R以及状态转移矩阵P等参量，代码如下：

import numpy as np

gama = 0.8  # 折扣因子
state_space = range(7)  # 状态空间
R = np.mat([5, 0, 0, 0, 0, 0, 10])  # 奖励矩阵
P = np.mat([[0.6, 0.4, 0, 0, 0, 0, 0],
            [0.4, 0.2, 0.4, 0, 0, 0, 0],
            [0, 0.4, 0.2, 0.4, 0, 0, 0],
            [0, 0, 0.4, 0.2, 0.4, 0, 0],
            [0, 0, 0, 0.4, 0.2, 0.4, 0],
            [0, 0, 0, 0, 0.4, 0.2, 0.4],
            [0, 0, 0, 0, 0, 0.4, 0.6]])  # 状态转移矩阵
I = np.eye(7)  # 单位矩阵

然后开始介绍三种计算价值函数的方法。

第一种方法：解析法

由于该过程的状态比较少，只有7种，因此状态转移矩阵P是个7*7的方阵，比较容易求逆，

书中给出了解析法求状态价值函数的公式：

我们直接根据该公式就能很快计算出状态价值函数矩阵，代码如下：

# 贝尔曼方程的解析解（状态价值函数）
def analytical_solution():
    V = np.dot(np.linalg.inv(I - gama * P), R.T)
    print('解析法计算出的状态价值函数矩阵：', np.array(V.T)[0])

结果如下：

解析法计算出的状态价值函数矩阵： [13.82495062  6.84054476  4.13147937  4.00458859  6.38056568 12.74439632
 27.07347466]

可以看出，并不是只有s1和s7才这两个状态有价值，其他状态也有价值，而且我们看到越接近s7的状态，其价值也越大。我们把s1看作是在南京市区，s7看作是在上海市区，价值就是当地的房价，很显然，越靠近上海市区，房价越高，靠近南京市区房价也比较高，但不比上海来得高，s6这个状态就好比在上海比较偏的郊区，这里房价虽然不及市区，但快赶上南京市区了，s6状态虽然没有环境的奖励（奖励为0），但因其靠近s7，它的状态价值就很高。相比之下，我们看s3和s4两个状态，由于离两个市区都比较远，他们的状态价值就相对比较低。

第二种方法：蒙特卡罗法

书中给出了蒙特卡罗法，伪代码如下所示：

 相应的python代码如下：

# 计算马尔可夫奖励过程价值的蒙特卡罗法
def Monte_Carlo_method():
    N = 100  # 随机游走的回合数
    H = 100  # 每回合时间长度
    V = []  # 初始化状态价值矩阵
    for initial_state in state_space:
        Gt = 0
        for i in range(N):
            g = 0
            current_state = initial_state
            for t in range(H):
                g += (gama ** t) * R[0, current_state]
                next_state = np.random.choice(state_space, p=np.array(P[current_state])[0])
                current_state = next_state
            Gt += g
        V.append(Gt / N)
    print('蒙特卡罗法计算出来的状态价值函数矩阵：', V)

蒙特卡罗法计算得比较慢，得到的结果如下：

蒙特卡罗法计算出来的状态价值函数矩阵： [13.309069553663592, 6.386829247553726, 3.880879633756227, 3.2405430246716445, 6.714685548749097, 13.40892479703931, 26.20065500434438]

可以看到和直接用解析法得到的结果还是有一些误差，毕竟这种迭代算法求的是近似解，而且回合数N越高，精度越高。

第三种方法：动态规划法

书中给出的伪代码如下：

 对应的python代码如下：
 

# 计算马尔可夫奖励过程价值的动态规划算法
def dynamic_planning():
    V_ = np.array([0., 0., 0., 0., 0., 0., 0.])
    V = np.array([1., 1., 1., 1., 1., 1., 1.])
    while np.linalg.norm(V - V_) > 0.0001:
        V = V_.copy()
        for state in state_space:
            future_value = 0
            for state_ in state_space:
                future_value += P[state, state_]*V[state_]
            V_[state] = R[0, state] + gama * future_value
    print('动态规划法计算出来的状态价值函数矩阵：', V_)

计算出来的结果如下：

动态规划法计算出来的状态价值函数矩阵： [13.82482916  6.84042312  4.13135743  4.00446626  6.38044296 12.74427328
 27.07335145]

这个结果和第一种方法还是十分接近的，而且计算速度也很快，书中提到，对于比较多的状态，由于矩阵求逆运算比较复杂，用这种动态规划算法是比较靠谱的，第一种解析法虽然精确，但不适用与维数多的马尔可夫过程。