（零基础可以看懂）强化学习中的动态规划（贝尔曼方程）（含代码）-《强化学习系列专栏第1篇》

介绍

    首先简单介绍下强化学习模型是干什么用的，强化学习模型其实就是一种做事情的策略或者方法，可以理解成一种做事的步骤。强化学习模型会基于当前的情形（确切地说，在强化学习模型里面叫做环境），做出可以获得最大收益的行为。举个例子来说，强化学习模型可以应用在自动驾驶里面，因为自动驾驶的系统可以等价于在当前的情形下，我应该做什么动作的这种系统。那么自然，就可以使用强化学习模型。

    强化学习模型是这样的一种模型：模型在学习阶段的时候，并没有像机器学习或深度学习中的算法一样，模型被告知应该往哪样的“标签”方向去训练。取而代之的是，模型在每一次迭代的时候，会被告知这种行为会获得正的收益，还是负的收益。然后根据这种收益的反馈，去更新模型。如果是使用了神经网络的模型，那么这种反馈就会去求解梯度，去更新神经网络中的参数。

    下面我们开始正式介绍强化学习，强化学习中一共有5个要素，分别是智能体（Agent）、状态（State）、行为（Action）、奖励（Reward）和策略（Policy）。智能体指的是在环境中存在的个体，比如自动驾驶中的汽车，玩Atari（类似小霸王那种游戏的统称）游戏中的人物或者角色等。状态指的是当前环境中的情形，比如自动驾驶中，传感器所识别到的图像，图像中各种东西的组合最终形成了这张图片，这个汽车所处的环境就可以叫做当前汽车所处的状态。行为指的是智能体根据策略所作出的行为。比如自动驾驶中，系统发出指令，告诉四个轮子往前旋转1弧度。奖励就是设计模型的设计者所给予模型的一种反馈，可以将其类比于监督学习中的标签。模型最终会往奖励最大的方向去更新。策略指的是给予当前的状态，智能体应该做出什么样的行为的一种决策。

    到此，我们初步的已经介绍了什么是强化学习、强化学习所使用的场景的举例、以及强化学习的基本概念，接下去，我们就介绍一种利用动态规划去求解强化学习中的模型的参数的例子。

动态规划求解模型参数

    动态规划求解方法是一种具有严格的数学推导的一种模型，它需要一个完整的、精确的环境模型。
想象这样的一个游戏，假设有一个4*4的方格，游戏开始时，你随机的从一个方格开始，然后每走一步，你将会扣1分（游戏开始时你的分数是0分），当你走到左上角的方格或者右下角的方格的时候，游戏停止。该游戏如下图所示，

    那么问题来了，如何走，你的最后的分数最高？
    首先，我们设定$t_0$时刻，所处的状态（格子）为$s_0$，然后我们走了一步，记为执行了$a_0$操作，此时，获得了奖励$r_1$，状态变成了$s_1$，然后我们又走了一步，记为执行了$a_1$操作，此时，获得了奖励$r_2$，状态变成了$s_2$。以此类推，我们可以得到一系列的结果，也就是
$$s_0,a_0,r_1,s_1,a_1,r_2,......,s_t,a_t,r_{t+1},s_{t+1},......$$
等序列。
    这个时候要注意，由于我们举的例子比较简单，因此，当执行了$a_0$操作后，状态从$s_0$转成$s_1$的时候，只有一种选项，而有些更复杂的情况下，当执行了某个操作后，状态是有可能从$s_t$转到$s_{t+1}$的时候有多种可能性，也就是说，在执行某个操作后，新的状态是某种概率分布。同理，在执行某个操作后，获得的奖励也是某种概率分布。
    这里要提一下马尔科夫性，我们从上面的例子中可以看到，当我当前状态确定后，我下一步将要走到哪儿去，只和我当前的状态有关，和之前的状态是无关的。因此它是符合马尔科夫性的。因此，以上的这种过程也叫做马尔科夫决策过程。
    基于以上的马尔科夫决策过程，从状态0开始，我们可以用公式写出总的奖励是
$$G_0=R_1+R_2+R_3+...$$
    上面的$R_1$代表着的意思是它是个随机变量，它的值还是没确定下来，$r_1$与其不同的是，$r_1$代表着一个确定的，已经获得的值。同时，我们设定从状态s开始，智能体按照策略$\pi$，进行决策操作后，所获得的回报是
$$G_t=R_{t+1}+R_{t+2}+......$$
    当时刻为无穷大时，我们为了让该值收敛，我们给后面加上一个折扣系数$\gamma$，于是公式就变成了
$$\begin{array}{rl}{G}_t& =R_{t+1}+\gamma \ast R_{t+2}+{\gamma }^2\ast R_{t+3}+{\gamma }^3\ast R_{t+4}+...{\gamma }^{T-(t+1)}\ast R_T\\ & =\sum _{k=t+1}^T{\gamma }^{k-(t+1)}\ast R_k\end{array}（1）$$
    为什么会存在一个T时刻？因为$\gamma$是小于1的，随着t趋于正无穷大，某一项的奖励值将收敛于0（当然，R不是发散的）。
    由于上面的$G_t$没办法衡量大小，于是，我们定义，状态s的状态价值函数是$G_t$的期望，也就是
$$\begin{array}{rl}{v}_{\pi }(s)& =E_{\pi }[G_t|S_t=s]\\ & =E_{\pi }[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s]（2）\end{array}$$
    这里要注意式子（1）和式子（2）k的区别。这个状态价值函数衡量了智能体在当前pi策略下，在给定状态下有多好，这里的有多好是以未来所获得的奖励的期望来衡量的。
    下面我们进行一个推导，来求解出$v_{\pi }(s)$和$v_{\pi }(s_{t+1})$之间的关系。我们先了解下各个符号代表什么意思，

①$\pi (a|s)$：代表在某种策略$\pi$的情况下，在状态$s$情况下，执行动作a的概率，举个例子来说，假如我们处于上面4*4格子中的10号格子，那么我们依据策略$\pi$，我们向上走的概率是0.6，向右走的概率是0.1，向下走的概率是0.1，向左走的概率是0.2。这个概率的数值就是$\pi (a|s)$
②$p(s^{\prime }|S_t=s,A_t=a)$：意思是在状态$s$情况下，执行了动作a后，状态从$s$变成（转移）到$s^{\prime }$的概率。
③$p(r|S_t=s,A_t=a,S_{t+1}=s^{\prime })$：意思是在状态$s$情况下，执行了动作a后，状态从$s$变成（转移）到$s^{\prime }$的后，获得奖励r的概率。
④$\gamma$：是折扣系数，意思是我要考虑多少的将来的情况到现在中去，也就是折算多少的未来到当前。

$$\begin{array}{rl}{v}_{\pi }(s)& =E_{\pi }[G_t|S_t=s]\\ & =E_{\pi }[R_{t+1}+\gamma G_{t+1}|S_t=s]\\ & =E_{\pi }[R_{t+1}|S_t=s]+E_{\pi }[\gamma G_{t+1}|S_t=s]\\ & =\sum _a\pi (a|s)\ast E_{\pi }[R_{t+1}|S_t=s,A_t=a]+\sum _a\pi (a|s)\ast E_{\pi }[\gamma G_{t+1}|S_t=s,A_t=a]\\ & =\sum _a\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|S_t=s,A_t=a)\ast E_{\pi }[R_{t+1}|S_t=s,A_t=a,S_{t+1}=s^{\prime }]\\ & +\sum _a\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|S_t=s,A_t=a)\ast E_{\pi }[\gamma G_{t+1}|S_t=s,A_t=a,S_{t+1}=s^{\prime }]\\ & =\sum _a\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|S_t=s,A_t=a)\ast \sum _{r}p(r|S_t=s,A_t=a,S_{t+1}=s^{\prime })\ast E_{\pi }[R_{t+1}|S_t=s,A_t=a,S_{t+1}=s^{\prime },R_{t+1}=r]\\ & +\sum _a\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|S_t=s,A_t=a)\sum _{r}p(r|S_t=s,A_t=a,S_{t+1}=s^{\prime })\ast E_{\pi }[\gamma G_{t+1}|S_t=s,A_t=a,S_{t+1}=s^{\prime },R_{t+1}=r]\\ & =\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast E_{\pi }[r|S_t=s,A_t=a,S_{t+1}=s^{\prime },R_{t+1}=r]\\ & +\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast \gamma E_{\pi }[G_{t+1}|S_t=s,A_t=a,S_{t+1}=s^{\prime },R_{t+1}=r]\\ & =\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast r+\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast \gamma E_{\pi }[G_{t+1}|S_{t+1}=s^{\prime }]\\ & =\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast r+\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)\ast \gamma v_{\pi }(s^{\prime })\\ & =\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)[r+\gamma v_{\pi }(s^{\prime })]\end{array}$$
    也就是说
$$v_{\pi }(s)=\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)[r+\gamma v_{\pi }(s^{\prime })]（3）$$
    上面这个式子说明了，$t$时刻的状态价值函数$v_{\pi }(s)$与$t+1$时刻的状态价值函数$v_{\pi }(s^{\prime })$的关系。式（3）也被称作贝尔曼方程（Bellman Equation）.

    找任务中找出一个最优的策略，就相当于是找出一种策略，能够在长期中，使智能体所获得的收益尽可能的大。如果策略A比策略B好，那么策略A应该在每一个状态上的状态价值函数的值应该$比B大或者相等$。因为但凡有一个状态A的价值函数的值比B小，那就说明在该状态上开始执行动作，策略A长远的收益一定比B小，同时，计算其他状态的状态价值函数的时候，有可能会考虑到该状态的价值函数，这会有可能导致其他状态的价值函数的值也比B小。如果这样的话，就说明策略A并没有比B好（大家深思下这句话，有点绕）。
这里，为了找出贝尔曼最优方程，我们定义动作价值函数，我们把在策略$\pi$下，在状态$s$时采取动作$a$的价值记为$q_{\pi }(s,a)$，它相当于是在状态$s$，在策略$\pi$下，采取动作$a$后，所有可能的回报的期望值。它被定义为，
$$q_{\pi }(s,a)=E_{\pi }[G_t|S_t=s,A_t=a]=E_{\pi }[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a]（4）$$
    我们可以很容量的推导出$v_{\pi }(s)$和$q_{\pi }(s,a)$的关系，为
$$v_{\pi }(s)=\sum _a\pi (a|s)q_{\pi }(s,a)$$
    我们设${\pi }^{\ast }$代表最优策略（最优策略有时候不止一个，可能有多个，这里我们统一用${\pi }^{\ast }$代表），则$v_{\ast }(s)$代表最优策略下的状态价值函数，$q_{\ast }(s,a)$代表最优策略下的动作价值函数，它们分别是，
$$v_{\ast }(s)=\underset{\pi }{\max }{v}_{\pi }(s)$$

$$q_{\ast }(s,a)=\underset{\pi }{\max }{q}_{\pi }(s,a)$$
    用$v_{\ast }(s)$来表示$q_{\ast }(s,a)$，就是，
$$\begin{array}{rl}{q}_{\ast }(s,a)& =\underset{\pi }{\max }{q}_{\pi }(s,a)\\ & =\underset{\pi }{\max }{E}_{\pi }[R_{t+1}+\sum _{k=1}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a]\\ & =E_{{\pi }^{\ast }}[R_{t+1}+\sum _{k=1}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a]\\ & =E_{{\pi }^{\ast }}[R_{t+1}|S_t=s,A_t=a]+E_{{\pi }^{\ast }}[\sum _{k=1}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a]\\ & =E_{{\pi }^{\ast }}[R_{t+1}|S_t=s,A_t=a]+\gamma v_{\ast }(S_{t+1})\\ & =E_{{\pi }^{\ast }}[R_{t+1}+\gamma v_{\ast }(S_{t+1})|S_t=s,A_t=a]\end{array}$$
    由于$R_{t+1}$已经被“剥离”出来了，因此$R_{t+2}$时刻的值与t时刻的$S_t$的值和$A_t$值无关。原因是马尔科夫性，$R_{t+2}$的值只受$S_{t+1}$和$A_{t+1}$影响。
    $v_{\ast }(s)$说明了在状态$s$下，按照最优策略去行动，所获得的回报的期望。
    $q_{\ast }(s,a))$说明了在状态$s$下，执行了动作a后，注意，这里的动作a是一个确切的动作，而不是动作的随机变量A，按照最优策略去行动，所获得的回报的期望。
    我们可以可以很容易想到，在最优策略的情况下，每个状态的价值一定等于这个状态下最优动作的期望回报。假如不等于最优动作的期望回报，那么回溯回来的时候，就会发现该状态的价值比其他策略下该状态的价值低。由此，我们可以写出以下公式，
$$\begin{array}{rl}{v}_{\ast }(s)& =\underset{a\in A(s)}{\max }{q}_{{\pi }^{\ast }}(s,a)\\ & =\underset{a}{\max }{E}_{{\pi }^{\ast }}[\sum _{k=0}^{\infty }{\gamma }^k{R}_{t+k+1}|S_t=s,A_t=a]\\ & =\underset{a}{\max }{E}_{{\pi }^{\ast }}[R_{t+1}+\gamma v_{\ast }(S_{t+1})|S_t=s,A_t=a]（5）\\ & =\underset{a}{\max }\sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+\gamma v_{\ast }(s^{\prime })]（6）\end{array}$$
式子（5）和式子（6）是$v_{\ast }$的贝尔曼方程的两种形式。$q_{\ast }$的贝尔曼方程如下，
$$\begin{array}{rl}{q}_{\ast }(s,a)& =E_{{\pi }^{\ast }}[R_{t+1}+\gamma v_{\ast }(S_{t+1})|S_t=s,A_t=a]（7）\\ & =E_{{\pi }^{\ast }}[R_{t+1}+\underset{a^{\prime }}{\max }{q}_{\ast }(S_{t+1},a^{\prime })|S_t=s,A_t=a]（8）\\ & =\sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+\gamma \underset{a^{\prime }}{\max }{q}_{\ast }(s^{\prime },a^{\prime })]（9）\end{array}$$
    式子（7）可以转换成式子（8）的原因我们之前说过，在最优策略的情况下，每个状态的价值一定等于这个状态下最优动作的期望回报，也就是在每个状态下，一定会执行有最优回报的动作，这就是贪心策略，它在局部贪心的时候也能够达到全局最优，因为最优的状态价值函数以及最优的动作价值函数已经保证了这一点。
    对于有限的马尔科夫决策过程来说（MDP，也就是说，状态有限，动作有限，奖励有限），式子（6）和（9）是有唯一解的。每个状态都会有一个式子（6）和（9）。通过解非线性方程组的方法，是可以求解的。但是一旦状态比较多，动作比较多，奖励比较多的时候，直接通过解方程组来求解，是很困难的，原因是需要很大的算力（计算复杂度极高）。那么如何近似的求解呢？

近似方法求解最优状态贝尔曼方程和最优动作贝尔曼方程

    我们继续以这个4*4格子例子为例

    然后我们仔细观察式子（10）
$$v_{\pi }(s)=\sum _a\pi (a|s)\sum _{s^{\prime },r}p(s^{\prime },r|S_t=s,A_t=a)[r+\gamma v_{\pi }(s^{\prime })]（10）$$
    上面这个式子说明了下一个状态$s^{\prime }$和当前状态$s$的迭代关系，有些人可能已经想到了，使用迭代法来近似求解。迭代法求解的一些充分条件这里就不说了，主要就是式子中的$\gamma$小于1。具体的大家可以参考Sutton的《强化学习》（第2版）。
    但是使用这个方法一定要注意！使用这个方法的前提是我们已经对环境（也就是状态）拥有完备的知识，也就是说，环境（状态）变化的规律我们是知道的，比如我们从上面的格子6向下走，就是走到了格子10，这样的环境知识我们是需要知道的。否则的话就不能使用该方法，如果环境对智能体是未知的话，我们有其他方法来求解最优策略，比如加入蒙特卡洛方法等等。大家可以想象，如果环境的状态有非常多的话，那么通过这种方法来求解计算复杂度将会非常高。因此，环境复杂的话，我们会更偏向于使用其他方法如加入蒙特卡洛采样等方法来求解最优策略。
    下面我直接将代码放上来。

$代码说明$

REWARD_SETTING 变量：指的是每走一步，要被扣多少分

    本质上，还是通过迭代法逼近动态规划中的最优解。当计算完成后，V会长这个样子，它说明了每个状态的价值。那么很明显，我们的策略就出来了，当我们处在6号格子的时候，6号格子的价值是-17.99915625。接着我们要决定往哪走，然后根据往状态价值大的方向走，于是我们往2号或者5号走（它们价值一样，那就随便选一个）。走到2号或者5号后，我们再往1号走，就走到终点了。

# encoding=utf-8
'''
Author: Haitaifantuan
Create Date: 2020-09-07 23:10:17
Author Email: 47970915@qq.com
Description: Should you have any question, do not hesitate to contact me via E-mail.
'''
import numpy as np
import copy


TOTAL_ROWS = 4
TOTAL_COLUMN = 4
TOTAL_ACTIONS_NUM = 4  # 0代表上，1代表右，2代表下，3代表左
STOP_CRITERION = 1e-4
GAMMA = 1
REWARD_SETTING = -1
ACTION_DICT = {0: '上', 1: '右', 2: '下', 3: '左'}
FOUR_ACTION_PROBABILITY = {'上': 0.25, '右': 0.25, '下': 0.25, '左': 0.25}  # 分别是走上、下、左、右的概率。
IDX_CHANGE_DICT = {'上': (-1, 0), '右': (0, 1), '下': (1, 0), '左': (0, -1)}  # 左边这个是行的索引的改变，右边这个是列的索引的改变

def get_current_reward_and_next_state(current_state, action):
    '''
    根据当前的状态，以及行为，计算当前行为的奖励以及下一个状态
    '''
    # 先判断是否到了终点，如果是终点，不管执行什么操作
    # 奖励都是0，并且都会回到终点
    row_idx, column_idx = current_state
    if (row_idx == 0 and column_idx == 0):
        return 0, (0, 0)

    if (row_idx == 3 and column_idx == 3):
        return 0, (3, 3)

    # 否则的话就计算下下一步的state和reward
    next_row_idx = row_idx + IDX_CHANGE_DICT[action][0]
    next_column_idx = column_idx + IDX_CHANGE_DICT[action][1]

    # 再判断是否在边缘，如果是的话，那就回到该位置。
    if next_row_idx < 0 or next_row_idx > TOTAL_ROWS - 1 or next_column_idx < 0 or next_column_idx > TOTAL_COLUMN - 1:
        return REWARD_SETTING, (row_idx, column_idx)
    else:
        return REWARD_SETTING, (next_row_idx, next_column_idx)


# 初始化状态价值函数V
V = np.zeros((TOTAL_ROWS, TOTAL_COLUMN))

# 开始迭代更新状态价值函数
iteration = 0
flag = True
while flag:
    delta = 0
    old_V = copy.deepcopy(V)
    # 遍历每一个状态，对其进行更新
    for row_idx in range(TOTAL_ROWS):
        for column_idx in range(TOTAL_COLUMN):
            new_final_value = 0
            # 根据sutton的《强化学习》第72页公式4.5进行更新
            for each_action in range(TOTAL_ACTIONS_NUM):
                action = ACTION_DICT[each_action]
                action_proba = FOUR_ACTION_PROBABILITY[action]
                current_action_reward, next_state = get_current_reward_and_next_state((row_idx, column_idx), action)
                new_final_value = new_final_value + action_proba * (1 * (current_action_reward +
                                                                         GAMMA * V[next_state[0]][next_state[1]]))

            V[row_idx][column_idx] = new_final_value


    delta = max(delta, abs(old_V - V).max())

    if delta < STOP_CRITERION:
        flag = False

    iteration += 1

print(V)
print(iteration)

代码复现、详细讲解及我的Github地址

完整代码地址：https://github.com/haitaifantuan/reinforcement_leanring