剖析强化学习 - 第一部分

分手后的思念是犯贱 2022-05-28 12:09 264阅读 0赞

**作者：Massimiliano Patacchiola**

**前言** \[本文是**对强化学习****的介绍**，适合已经有一些机器学习背景，并且懂一些数学和Python的读者。当我研究一种新算法时，我总是希望了解底层机制，从这个意义上讲，使用一种编程语言从头开始实现算法对理解算法是有帮助的。我在这篇文章中采用了这种方法，虽然需要花更长时间阅读但值得这样。我不是以英语为母语的人，所以如果你发现一些难以理解的错误句子，请在下面的评论中说明，我会解决它\]

当我开始学习强化学习时，我没有找到任何从基础上解释什么是强化学习的好的在线资源。大多数（非常好）博客都关注现代方法（DeepReinforcement Learning），都会介绍[Bellman方程][Bellman]但没有令人满意的解释。我把注意力转向这本书中学习，**Russel和Norvig**的《**人工智能：一种现代的方法》**。

![20180405104223093][]

本文基于第二版的**第****17章**，可以被认为是本章的延伸阅读。我将使用与作者相同的数学符号，通过这种方式，您可以使用这本书来覆盖本文一些缺失的部分，反之亦然。你可以在[我的github版本库中][github]找到这篇文章中使用的完整代码，并附带pdf版本的帖子。在下一节我将介绍**马尔可夫链**，如果你已经知道这个概念，可以跳过下一节...

## Andrey Markov ##

[Andrey Markov][]是一位研究随机过程的俄罗斯数学家，他对遵循一系列相关事件的系统特别感兴趣。1906年，马尔可夫对他称之为**链（****Chain）的**离散过程的结果产生了兴趣。一个**马尔可夫链**具有一组**状态** S=\{ s0,s1,...,sm\}以及可以从一个状态连续移动到另一个状态的**过程**，每一次移动都是一个**步骤（****Step）**，并基于**转移模型** T。您应该努力记住粗体字，因为我们将在文章的其余部分广泛使用它们。总结马尔可夫链的定义如下：

1.一组可能的状态：S=\{ s0,s1,...,sm\}

2.初始状态：s0

3.转移模型：T(s,s′)

马尔科夫链中有一个我没有提及的重要特性，马尔可夫链是基于**马尔可夫属性**，**当前状态的未来变化条件上独立于过去的状态**，就是过程中当前的状态只依赖于它在t-1时的状态。一个例子可以简化马尔可夫链的理解，假设我们有一个只有两个状态s0和s1，其中s0是初始状态。该过程在s0保持90％的时间，剩下的10％的时间它移动到s1，当过程处于状态s1时它将保持50％的时间。根据这些数据，我们可以创建一个如下的**转移矩阵**T ：

![20180405104809884][]

转移矩阵总是一个方形矩阵，并且由于我们处理的是概率分布，因此所有条目都在0和1之间，单个行总和为1。**我们可以用图形表示马尔可夫链**。在下面的表示中，链的每个状态都是一个节点，转移概率用边表示，最高的概率具有最厚边：

![20180405104831197][]

到目前为止，我们没有提及**时间**，但我们必须这样做，因为马尔可夫链是随时间发展的动态过程。假设我们猜测这个过程在3个步骤之后和50个步骤之后的状态是什么，我们该怎么做？我们注意到这是个具有有限数量的状态并且时间均匀的链，这意味着转移矩阵不随时间变化。考虑到这些假设，**我们可以将k步转移概率计算为转移矩阵的k次幂**，让我们在Numpy中计算：

import numpy as np
    
    #Declaring the Transition Matrix T
    T = np.array([[0.90, 0.10],
                  [0.50, 0.50]])
    
    #Obtaining T after 3 steps
    T_3 = np.linalg.matrix_power(T, 3)
    #Obtaining T after 50 steps
    T_50 = np.linalg.matrix_power(T, 50)
    #Obtaining T after 100 steps
    T_100 = np.linalg.matrix_power(T, 100)
    
    #Printing the matrices
    print("T: " + str(T))
    print("T_3: " + str(T_3))
    print("T_50: " + str(T_50))
    print("T_100: " + str(T_100))

T: [[ 0.9  0.1]
        [ 0.5  0.5]]
    
    T_3: [[ 0.844  0.156]
          [ 0.78   0.22 ]]
    
    T_50: [[ 0.83333333  0.16666667]
           [ 0.83333333  0.16666667]]
    
    T_100: [[ 0.83333333  0.16666667]
            [ 0.83333333  0.16666667]]

现在我们定义代表k = 0时系统状态的**初始分布**，系统由两个状态组成，我们可以将初始分布建模为具有两个元素的向量，向量的第一个元素表示停留在状态s0的概率，第二个元素是停留在状态s1的概率。我们假设从s0开始，向量v 代表下列形式的初始分布：

v=（1，0）

我们可以计算出**在****k次迭代之后的特定状态的概率：**初始分布和转移矩阵相乘：v⋅Tk。让我们在Numpy中实现这一点：

import numpy as np
    
    #Declaring the initial distribution
    v = np.array([[1.0, 0.0]])
    #Declaring the Transition Matrix T
    T = np.array([[0.90, 0.10],
                  [0.50, 0.50]])
    
    #Obtaining T after 3 steps
    T_3 = np.linalg.matrix_power(T, 3)
    #Obtaining T after 50 steps
    T_50 = np.linalg.matrix_power(T, 50)
    #Obtaining T after 100 steps
    T_100 = np.linalg.matrix_power(T, 100)
    
    #Printing the initial distribution
    print("v: " + str(v))
    print("v_1: " + str(np.dot(v,T)))
    print("v_3: " + str(np.dot(v,T_3)))
    print("v_50: " + str(np.dot(v,T_50)))
    print("v_100: " + str(np.dot(v,T_100)))

v: [[ 1.  0.]]
    
    v_1: [[ 0.9  0.1]]
    
    v_3: [[ 0.844  0.156]]
    
    v_50: [[ 0.83333333  0.16666667]]
    
    v_100: [[ 0.83333333  0.16666667]]

**发生了什么？**该过程从s0开始经过一次迭代后，我们有90%确定它仍然处于这个状态。这很容易理解，我们的转移模型说明，该过程有90％的概率可以留在s0状态，没有什么新意。展望k = 3时的状态分布，我们注意到有一些不同，未来有更多的可能性进入不同的状态。如果我们想找到经过三次迭代后状态处于s0的概率，我们应该把所有可能导致s0的分支概率加起来。一张图片胜过千言万语：

![20180405105059248][]

在k = 3时处于s0的概率由（0.729 + 0.045 +0.045 + 0.025）给出，它等于0.844，我们得到了相同的结果。现在假设在开始时我们对于过程的开始状态有一些不确定性，让我们定义另一个起始向量，如下所示：

v =（0.5，0.5）

就是说，有50％的概率，我们可以从s0开始。再次运行Python脚本，我们在1次，3次，50次和100次迭代后打印结果：

v: [[ 0.5, 0.5]]
    
    v_1: [[ 0.7  0.3]]
    
    v_3: [[ 0.812  0.188]]
    
    v_50: [[ 0.83333333  0.16666667]]
    
    v_100: [[ 0.83333333  0.16666667]]

这一次在k = 3时进入s0的概率较低（0.812），但在长时间执行后，我们有相同的结果（0.8333333）。**长时间运行发生了什么？**经过50和100次迭代的结果是相同的，不管开始是什么概率始终`v_50`=`v_100`。该链**收敛****于平衡状态，**意味着随着时间的推移，它忘记了开始的分配。但我们必须小心，收敛并不总能得到保证，马尔可夫链的动态性可能非常复杂，特别是可能存在**瞬态和循环状态**。对于我们的讨论范围，这些就已经足够了。我建议你看看[setosa.io][]博客，因为他们有一个用于马尔可夫链可视化的交互式页面。

## 马尔科夫决策过程（MDP） ##

在强化学习中，经常使用与马尔可夫链有关的概念，下面我们说说**马尔可夫决策过程（****MDP）**。MDP是马尔可夫链的重新释义，包括**Agent**和**决策**过程。MDP由这些组件定义：

1.  可能的状态集合: S=\{ s0,s1,...,sm\}

2.  初始状态: s0

3.  可能的动作（action）集合: A=\{ a0,a1,...,an\}

4.  转移模型: T(s,a,s′)

5.  奖励（Reward）函数: R(s)

正如您所看到的，我们正在引入一些关于马尔可夫链的新元素，特别是依赖于当前状态、下一个状态和Agent动作的转移模型。转移模型返回在状态s完成动作a到达状态s′的概率。但是给定s和a的该模型有条件地独立于以前的所有状态和行为（马尔可夫属性）。此外还有**奖励函数** R(s)，每当Agent从一个状态移动到另一个状态时，它都会返回一个实数值（注意：将奖励函数定义为仅依赖于s可能会引起混淆，Russel和Norvig在书中使用这种表示法来简化描述，但并没有改变问题的本质）。由于有奖励功能，我们可以说**一些状态更期望某些状态，**因为当Agent在这些状态移动时，它会得到更高的奖励。相反，有些**状态恰恰相反**，因为当Agent移动到这些状态时，会收到负的奖励。

·        **问题：**Agent必须最大化获得的奖励，避免哪些返回负值的状态，而是选择返回正值的状态。

·        **解决方案：**找到一个**策略** π(s) 它返回最高奖励的动作。

Agent可以尝试不同的策略，但只有其中一个可以被认为是**最优策略**，用π\*表示，它会产生最高的预期效用(utility)。现在开始介绍我将在所有帖子中使用的例子，这个例子的灵感来自Russel和Norving在他们书中第17.1章中提出的简单环境，假设我们有一台必须到达充电站的**清洁机器人**，我们的简单世界是一个4×3矩阵，其中起点s0在（1,1）位置，（4,3）位置是充电站，（4,2）位置是危险楼梯，以及（2,2）位置是障碍物。**机器人必须找到到达充电站的最佳方式**（奖励\+1）**并避免跌落楼梯**（奖励\-1）。每当机器人做出决定时，就有可能受到随机因素（例如地面滑溜，一只恶猫刺痛机器人）的干扰，这会使机器人在20％的时间内偏离原始路径。机器人向前走时将在10％的情况下向左，10％的情况向右的状态，如果机器人撞到墙壁或障碍物，它会弹回到原来的位置。这个世界的主要特点如下：

·        离散的时间和空间

·        完全可观察

·        无限的视野

·        已知的转移模型

环境是**完全可观察的**，这意味着**机器人总是知道它处于哪个状态**。应该进一步解释术语无限时域，**无限时域**意味着，**没有一个固定的时间限制**。如果Agent有相同的两种终止状态的策略，它将永远持续下去，这个假设并不意味着每个剧集(episode )都需要Agent走过一系列无限状态，当两个终止状态中的一个到达时，该剧集结束。下面说明了这个世界和转移模型的表示，小心Russell和Norvig使用的状态索引，这可能会让人困惑，他们从左下角开始按列和行来命名世界的每个状态。

![2018040510520954][]

机器人的目标是找到**到达充电站的最佳方式**，这意味着什么才是**最好的方式**呢？根据机器人在每个中间状态下获得的奖励类型，我们可以有不同的最优策略π\*。假设我们正在设计机器人的固件，根据电池级别，我们在每个时间步骤给予不同的奖励，两个终止状态的奖励保持不变（充电器= + 1，楼梯= -1），（2,2）处的障碍不是有效的状态，因此没有与其相关的奖励。鉴于这些假设，我们可以有四种不同的情况：

1.R(s)≤−1.6284R(s)≤−1.6284 **极低的电池**

2.−0.4278≤R(s)≤−0.085−0.4278≤R(s)≤−0.085 **较低的电池**

3.−0.0221≤R(s)≤0−0.0221≤R(s)≤0 **稍低的电池**

4.R(s)>0R(s)>0 **完全充满**

对于这些条件中的每一个，我们都可以尝试猜测Agent将选择哪种策略。在**极低的电池**情况下，Agent会受到如此高的惩罚，以至于只想尽快让伤痛停止，生活是如此痛苦，以至于坠下楼梯是一个不错的选择。在**非常低的电池**情况下，Agent将选择到充电站的最短路径，它不关心是否跌倒。在**电池电量稍低的**情况下，机器人完全不承担风险，并躲避楼梯避免撞墙的代价。最后，在**充满电的**情况下，Agent避免了两个终止状态并保持在稳定状态，在每个时间步骤都获得积极回报。

![20180405105325496][]

到目前为止，我们知道哪些策略可以在具有确定奖励的特定环境中出现，但我仍然没有谈到：Agent如何选择最佳策略？

## Bellman方程 ##

上一节结束时提出了一个问题：Agent**如何选择最佳策略？**为了回答这个问题，我将介绍Bellman方程。首先，必须找到一种方法来比较两种策略，我们可以使用每个状态给出的奖励来获得**状态序列效用的**度量。我们定义状态**历史** h的效用如：

Uh=R(s0)+γR(s1)+γ2R(s2)+...+γnR(sn)

前面的公式定义了状态序列的**折扣奖励**，其中的 γ∈\[0,1\]被称为**折扣因子**。折扣因子描述了Agent关于当前奖励对于未来奖励的偏好，折扣因子为1.0会导致先前公式退化为**累加奖励**。折扣奖励不是我们估算效用的唯一方式，但它是减少问题的方法。例如，在无限状态序列的情况下，折扣奖励会产生**确定的效用**（使用无限序列的奖励汇总），而且我们还可以使用每个时间步骤获得的平均奖励来比较无限序列。**如何比较单一状态的效用**？效用U(s)可以定义为：

![20180405105407816][]

让我们记住效用是根据策略π定义的，为了简洁起见我上面没有提到。一旦我们有了效用，我们如何为下一个状态选择最佳的动作？利用**最大期望效用**原则，即理性的Agent应该选择一种最大化其预期效用的动作。我们离Bellman公式更近了一步，一个状态s的效用和与它相邻的状态 s′的效用相关，即：

![20180405105437656][]

我们刚刚推导了**Bellman方程**！使用Bellman方程，Agent可以估计最佳的动作并找出最佳策略。我们试着剖析这个等式,**首先**，术语R(s)是我们必须在等式中加上的东西，我们知道效用必须考虑在状态s给予该状态的奖励；**其次**，请注意，该公式使用转移模型T乘以下一个状态s的效用。如果你认为这是有道理的，那么发生概率较低的状态（如我们简化世界中左侧和右侧移动的概率为10％）在总和中将具有最低的权重。

我们将在简化的4x3环境中使用我们的清洁机器人，以经验方式测试Bellman方程。在这个例子中，每个非终止状态的奖励是 R(s)=−0.04。我们可以想当然每个状态的效用值，**因为目前你不需要知道我们如何得到这些值，可以想象它们是神奇的**。以同样神奇的方式，我们为世界获得了最佳策略（重新检查我们从Bellman方程中获得的结果是否合理）。这张图片非常重要，牢记在心。

![20180405105456796][]

在我们的例子中，我们假设**机器人从状态（1,1）开始**，使用Bellman方程，我们必须**找到** UP、LEFT、DOWN和RIGHT中**效用最高的动作**。我们没有最优策略，但我们有每个状态的转移模型和效用值。你必须记住我们环境的两个主要规则，首先，如果机器人碰到墙，它会回到之前的状态；其次，根据转移模型，选定的动作仅以80％的概率执行。我不想说明简陋的数字，而是想向您展示可能结果的图形说明：

![20180405105521490][]

对于每个可能的结果，我都说明了效用和转移模型给出的概率，这对应于Bellman方程的第一部分。下一步是计算**效用和转移概率之间的乘积，然后总结每个动作的值**。

![20180405105627749][]

我们发现，对于状态（1,1）**，UP的行为具有最高的价值**。这符合我们神奇般获得的最佳策略。Bellman方程的这一部分返回最大化后续状态的期望效用的动作，这是最优策略应该做的事情：

![20180405105653885][]

现在我们有了所有的元素，并且我们可以将值带入到Bellman方程得到状态（1,1）的效用：

U(s11)=−0.04+1.0×0.7456=0.7056

这就是Bellman方程的工作原理！下面说明在我们的模拟世界中使用的等式的**Python实现**，我们将使用前面部分的相同术语。我们的世界有4x3 = 12个可能的状态，起始矢量包含12个值，转移矩阵是一个巨大的12x12x4矩阵（12个起始状态，12个下一个状态，4个动作），其中大部分值为零（我们只能从一个状态移动到其邻近状态）。我使用脚本生成转移矩阵，并将其保存为Numpy矩阵（[可以在此处下载][Link 1]）。在脚本中我定义了这个函数`return_state_utility()`，它是Bellman方程的一个实现。使用这个函数，我们将打印状态（1,1）的效用并检查它是否与我们之前找到的一样：

import numpy as np
    
    def return_state_utility(v, T, u, reward, gamma):
        """Return the state utility. @param v the state vector @param T transition matrix @param u utility vector @param reward for that state @param gamma discount factor @return the utility of the state """
        action_array = np.zeros(4)
        for action in range(0, 4):
            action_array[action] = np.sum(np.multiply(u, np.dot(v, T[:,:,action])))
        return reward + gamma * np.max(action_array)
    
    def main():
        #Starting state vector
        #The agent starts from (1, 1)
        v = np.array([[0.0, 0.0, 0.0, 0.0, 
                       0.0, 0.0, 0.0, 0.0, 
                       1.0, 0.0, 0.0, 0.0]])
    
        #Transition matrix loaded from file
        #(It is too big to write here)
        T = np.load("T.npy")
    
        #Utility vector
        u = np.array([[0.812, 0.868, 0.918,   1.0,
                       0.762,   0.0, 0.660,  -1.0,
                       0.705, 0.655, 0.611, 0.388]])
    
        #Defining the reward for state (1,1)
        reward = -0.04
        #Assuming that the discount factor is equal to 1.0
        gamma = 1.0
    
        #Use the Bellman equation to find the utility of state (1,1)
        utility_11 = return_state_utility(v, T, u, reward, gamma)
        print("Utility of state (1,1): " + str(utility_11))
    
    if __name__ == "__main__":
        main()

Utility of state (1,1): 0.7056

很好，我们获得了完全相同的价值！到现在为止，**我们认为效用价值神奇般地出现**。我们不想使用魔法，而是想**找到一种算法来获得这些值**。有一个问题，对于n个可能的状态有n个Bellman方程，每个方程包含n个未知数。使用任何线性代数包将有可能解决这些方程，问题是因为max的操作导致它们不是线性的。该怎么办？我们可以使用值迭代算法...

## 值迭代算法 ##

Bellman方程是用于求解MDP的值迭代算法的核心。**我们的目标是找到每个状态的效用（也称为价值）**。正如前面所说的，不能使用线性代数库，我们需要一种迭代方法，从任意初始效用值（通常为零）开始，然后使用Bellman方程计算一个状态的效用，并将其分配给状态，这个迭代过程被称为**Bellman更新**。无限地应用Bellman更新可以**确保达到平衡，**一旦我们达到了平衡状态，我们就可以找到我们正在寻找的效用值，并且我们可以用它们来估计每个状态的最佳动作。如何理解算法何时达到平衡？我们需要一个**停止标准**，考虑到两次连续迭代之间的效用，我们可以在没有状态变化的情况下停止算法。

![20180405110330479][]

这个结果是**收缩属性的结果**，我会跳过它，因为它在书中的第17.2章已经很好地解释了。好的，是时候用Python实现算法了。我将重新使用`return_state_utility()`函数来更新效用向量`u`。

import numpy as np
    
    def return_state_utility(v, T, u, reward, gamma):
        """Return the state utility. @param v the value vector @param T transition matrix @param u utility vector @param reward for that state @param gamma discount factor @return the utility of the state """
        action_array = np.zeros(4)
        for action in range(0, 4):
            action_array[action] = np.sum(np.multiply(u, np.dot(v, T[:,:,action])))
        return reward + gamma * np.max(action_array)
    
    def main():
        #Change as you want
        tot_states = 12
        gamma = 0.999 #Discount factor
        iteration = 0 #Iteration counter
        epsilon = 0.01 #Stopping criteria small value
    
        #List containing the data for each iteation
        graph_list = list()
    
        #Transition matrix loaded from file (It is too big to write here)
        T = np.load("T.npy")
    
        #Reward vector
        r = np.array([-0.04, -0.04, -0.04,  +1.0,
                      -0.04,   0.0, -0.04,  -1.0,
                      -0.04, -0.04, -0.04, -0.04])    
    
        #Utility vectors
        u = np.array([0.0, 0.0, 0.0,  0.0,
                       0.0, 0.0, 0.0,  0.0,
                       0.0, 0.0, 0.0,  0.0])
        u1 = np.array([0.0, 0.0, 0.0,  0.0,
                        0.0, 0.0, 0.0,  0.0,
                        0.0, 0.0, 0.0,  0.0])
    
        while True:
            delta = 0
            u = u1.copy()
            iteration += 1
            graph_list.append(u)
            for s in range(tot_states):
                reward = r[s]
                v = np.zeros((1,tot_states))
                v[0,s] = 1.0
                u1[s] = return_state_utility(v, T, u, reward, gamma)
                delta = max(delta, np.abs(u1[s] - u[s])) #Stopping criteria 
            if delta < epsilon * (1 - gamma) / gamma:
                    print("=================== FINAL RESULT ==================")
                    print("Iterations: " + str(iteration))
                    print("Delta: " + str(delta))
                    print("Gamma: " + str(gamma))
                    print("Epsilon: " + str(epsilon))
                    print("===================================================")
                    print(u[0:4])
                    print(u[4:8])
                    print(u[8:12])
                    print("===================================================")
                    break
    
    if __name__ == "__main__":
        main()

在收敛期间看看每个效用的稳定是很有趣的，我使用matplotlib绘制**每个状态25次迭代**的**效用值**。

![20180405110412178][]

使用相同的代码，针对折扣因子`gamma`使用不同值执行**不同的模拟**，当折扣因子接近1.0时，我们对效用的预测变得更加精确，**在** `gamma = 1.0` **的极限情况下，算法将永远不会结束，因为我们永远不会达到停止规则**。

=================== FINAL RESULT ==================
    Iterations: 9
    Delta: 0.000304045
    Gamma: 0.5
    Epsilon: 0.001
    ===================================================
    [[ 0.00854086  0.12551955  0.38243452  1.        ]]
    [[-0.04081336  0.          0.06628399 -1.        ]]
    [[-0.06241921 -0.05337728 -0.01991461 -0.07463402]]
    ===================================================

=================== FINAL RESULT ==================
    Iterations: 16
    Delta: 0.000104779638547
    Gamma: 0.9
    Epsilon: 0.001
    ===================================================
    [[ 0.50939438  0.64958568  0.79536209  1.        ]]
    [[ 0.39844322  0.          0.48644002 -1.        ]]
    [[ 0.29628832  0.253867    0.34475423  0.12987275]]
    ===================================================

=================== FINAL RESULT ==================
    Iterations: 29
    Delta: 9.97973302774e-07
    Gamma: 0.999
    Epsilon: 0.001
    ===================================================
    [[ 0.80796344  0.86539911  0.91653199  1.        ]]
    [[ 0.75696623  0.          0.65836281 -1.        ]]
    [[ 0.69968285  0.64882069  0.6047189   0.38150244]]
    ===================================================

还有另一种算法可以让我们找到效用向量，同时还有一个最优策略，那就是策略迭代算法。

## 策略迭代算法 ##

利用值迭代算法，我们可以估计每个状态的效用。我们仍然想找到估算最佳策略的另一种方法，在本节中，我将向您展示如何**使用策略迭代算法来寻找**最大化期望奖励**的最优策略**，没有其它策略会产生比最优策略π\*更多的奖励。**策略迭代保证收敛，当前策略及其效用函数是最优策略和最优效用函数**。首先，我们定义一个策略π它为每个状态分配一个动作，我们可以为此策略分配随机动作，这并不重要。使用`return_state_utility()`函数（Bellman方程），我们可以计算策略的预期效用。有一个好消息,我们并不需要**Bellman方程**的**完整版本，**即：

![20180405110647313][]

由于我们有一个策略和策略中每个状态相关联的动作，我们可以去掉max运算符并使用**Bellman等式**的**简化版本**：

![20180405110457115][]

当评估策略时，我们可以改进它，该策略**的改进**是算法的第二个和最后一个步骤。我们的环境具有有限数量的状态及有限数量的策略，每次迭代都会产生更好的策略。我实现了一个函数`return_policy_evaluation()`，它包含Bellman方程的简化版本。此外，我们需要函数`return_expected_action()`返回基于当前值`u`和`T`具有最高效用的动作。为了检查发生了什么，我还创建了一个print函数，它将策略向量`p`中包含的每个动作映射到一个符号并将其打印在终端上。

import numpy as np
    
    def return_policy_evaluation(p, u, r, T, gamma):
      """Return the policy utility. @param p policy vector @param u utility vector @param r reward vector @param T transition matrix @param gamma discount factor @return the utility vector u """
      for s in range(12):
        if not np.isnan(p[s]):
          v = np.zeros((1,12))
          v[0,s] = 1.0
          action = int(p[s])
          u[s] = r[s] + gamma * np.sum(np.multiply(u, np.dot(v, T[:,:,action])))
      return u
    
    def return_expected_action(u, T, v):
        """Return the expected action. It returns an action based on the expected utility of doing a in state s, according to T and u. This action is the one that maximize the expected utility. @param u utility vector @param T transition matrix @param v starting vector @return expected action (int) """
        actions_array = np.zeros(4)
        for action in range(4):
           #Expected utility of doing a in state s, according to T and u.
           actions_array[action] = np.sum(np.multiply(u, np.dot(v, T[:,:,action])))
        return np.argmax(actions_array)
    
    def print_policy(p, shape):
        """Printing utility. Print the policy actions using symbols: ^, v, <, > up, down, left, right * terminal states # obstacles """
        counter = 0
        policy_string = ""
        for row in range(shape[0]):
            for col in range(shape[1]):
                if(p[counter] == -1): policy_string += " * "            
                elif(p[counter] == 0): policy_string += " ^ "
                elif(p[counter] == 1): policy_string += " < "
                elif(p[counter] == 2): policy_string += " v "           
                elif(p[counter] == 3): policy_string += " > "
                elif(np.isnan(p[counter])): policy_string += " # "
                counter += 1
            policy_string += '\n'
        print(policy_string)

现在我将在main函数的循环中使用这些函数，这是一个策略迭代算法的实现。我声明了一个新的向量`p`，其中包含每个状态的动作，该算法的**停止条件**是两次连续迭代的效用向量之间的差，当步骤中对效用的改进不会发生变化（或发生非常小的变化）时，算法终止。

def main():

gamma = 0.999
        epsilon = 0.0001
        iteration = 0
        T = np.load("T.npy")
        #Generate the first policy randomly
        # NaN=Nothing, -1=Terminal, 0=Up, 1=Left, 2=Down, 3=Right
        p = np.random.randint(0, 4, size=(12)).astype(np.float32)
        p[5] = np.NaN
        p[3] = p[7] = -1
        #Utility vectors
        u = np.array([0.0, 0.0, 0.0,  0.0,
                      0.0, 0.0, 0.0,  0.0,
                      0.0, 0.0, 0.0,  0.0])
        #Reward vector
        r = np.array([-0.04, -0.04, -0.04,  +1.0,
                      -0.04,   0.0, -0.04,  -1.0,
                      -0.04, -0.04, -0.04, -0.04])
    
        while True:
            iteration += 1
            #1- Policy evaluation
            u_0 = u.copy()
            u = return_policy_evaluation(p, u, r, T, gamma)
            #Stopping criteria
            delta = np.absolute(u - u_0).max()
            if delta < epsilon * (1 - gamma) / gamma: break
            for s in range(12):
                if not np.isnan(p[s]) and not p[s]==-1:
                    v = np.zeros((1,12))
                    v[0,s] = 1.0
                    #2- Policy improvement
                    a = return_expected_action(u, T, v)         
                    if a != p[s]: p[s] = a
            print_policy(p, shape=(3,4))
    
        print("=================== FINAL RESULT ==================")
        print("Iterations: " + str(iteration))
        print("Delta: " + str(delta))
        print("Gamma: " + str(gamma))
        print("Epsilon: " + str(epsilon))
        print("===================================================")
        print(u[0:4])
        print(u[4:8])
        print(u[8:12])
        print("===================================================")
        print_policy(p, shape=(3,4))
        print("===================================================")
    
    if __name__ == "__main__":
        main()

使用`gamma=0.999`和`epsilon=0.0001`运行脚本，在22次迭代后收敛，我们得到下面的结果:

=================== FINAL RESULT ==================
    Iterations: 22
    Delta: 9.03617490833e-08
    Gamma: 0.999
    Epsilon: 0.0001
    ===================================================
    [ 0.80796344  0.86539911  0.91653199  1.        ]
    [ 0.75696624  0.          0.65836281 -1.        ]
    [ 0.69968295  0.64882105  0.60471972  0.38150427]
    ===================================================
     >   >   >   *  
     ^   #   ^   *  
     ^   <   <   <  
    ===================================================

**算法返回的最终策略等于最优策略**。此外，使用简化的Bellman方程，该算法设法找到效用向量的最优值。如果我们看看策略的演变，我们会注意到一些有趣的现象，开始时，策略是随机生成的，经过四次迭代后，算法找到次最优策略并坚持到第10次迭代，直到找到最优策略。从迭代10到迭代22，算法根本不改变策略。一个**次优的策略在无模型强化学习中可能是一个问题**，因为贪婪的Agent可能停留在次优策略，就目前对我们来说它不是一个问题。

![20180405112448717][]

策略迭代和值迭代，哪个更好？如果您有很多动作，或者您从合理的策略开始，请选择策略迭代；如果您有很少的动作，并且转移是非循环的，那么选择值迭代。如果你想从两个中获得最好的结果，那么请看看[修改过的策略迭代算法][Link 2]。

## 使用线性代数的策略评估 ##

我说过，从Bellman方程中消除max操作使得我们的世界变得更加简单，因为我们可以使用任何线性代数包来计算效用。在最后一节中，我想说明**如何使用线性代数方法达到相同的结果**。在Bellman方程中，有一个带n个变量和n个约束的线性系统，我们需要处理矩阵和向量。给定一个策略`p`、与状态`s`相关的动作、奖励向量`r`、转移矩阵`T`和折扣因子`gamma`，我们可以用一行代码来估计效用：

u\[s\]**=** np**.**linalg**.**solve(np**.**identity(12)**-** gamma**\***T\[:,:,p\[s\]\],r)\[s\]

我使用了Numpy方法`np.linalg.solve`，它将系数矩阵`A`和一个数组`b`作为参数输入，该方法返回给系统`A x = b`的解。对于参数矩阵`A`我传递单位矩阵`I`和`gamma * T`之间的差，对于参数数组`b`我传递奖励矢量`r`。**为什么用**`I -gamma*T`**作为第一个参数传递** ？我们可以从简化的Bellman方程推导出这个值：

![20180405112519196][]

事实上，我们可以在Numpy中实现最后一个等式得到`u`：

u[s]=np.dot(np.linalg.inv(np.identity(12)-gamma*T[:,:,p[s]]),r)[s]

不过，我更喜欢使用`np.linalg.solve`，因为它更具可读性。

## 结论 ##

在第一部分中，我总结了强化学习的基本思想。例如，我使用了具有预定义转移模型的有限环境。**如果我们没有转移模型，会发生什么？**在接下来的部分中，我将介绍无模型强化学习，它用一组新的有趣工具回答这个问题。你可以在我的github存储库中找到完整的[代码][github]。

## 索引 ##

1. **[\[第一篇][Link 3]\]**马尔科夫决策过程，贝尔曼方程，值迭代和策略迭代算法。

2.  [\[第二篇\]][Link 4]蒙特卡罗概念，蒙特卡洛方法，预测与控制，广义策略迭代，Q函数。

3.  [\[第三篇\]][Link 5]时间差分概念，动物学习，TD(0),TD(λ)和资格痕迹，SARSA，Q-learning。

4.  [\[第四篇\]][Link 6]Actor-Critic方法背后的神经生物学，计算Actor-Critic方法，Actor-only和Critic-only方法。

5.  [\[第五篇\]][Link 7]进化算法介绍，强化学习中的遗传算法，遗传算法的策略选择。

6.  [\[第六篇\]][Link 8]强化学习应用，多臂老虎机，山地车，倒立摆，无人机着陆，难题。

7.  [\[第七篇][Link 9]\]函数逼近概念，线性逼近器，应用，高阶逼近器。

8. **[\[第八篇\]][Link 10]** 非线性函数逼近，感知器，多层感知器，应用，政策梯度。

## 资源 ##

The [dissecting-reinforcement-learning][github] repository.

The [setosa blog][setosa.io] containing a good-lookingsimulator for Markov chains.

Official [github repository][] for the book *“ArtificialIntelligence: a Modern Approach”*.

## 参考 ##

Bellman,R. (1957). A Markovian decision process (No. P-1066). RAND CORP SANTA MONICACA.

Russell,S. J., Norvig, P., Canny, J. F., Malik, J. M., & Edwards, D. D. (2003).Artificial intelligence: a modern approach (Vol. 2). Upper Saddle River:Prentice hall.

[Bellman]: https://en.wikipedia.org/wiki/Bellman_equation
[20180405104223093]: /images/20220528/b82b49d0ea9b47238dca154aea38636f.png
[github]: https://github.com/mpatacchiola/dissecting-reinforcement-learning
[Andrey Markov]: https://en.wikipedia.org/wiki/Andrey_Markov
[20180405104809884]: /images/20220528/cc5d087c776f4751b8505c9d16497e40.png
[20180405104831197]: /images/20220528/f77d4c8c0c3b42f3b899df6b50a0f05d.png
[20180405105059248]: /images/20220528/b5c4ed0614ed4a3088449be969ebb961.png
[setosa.io]: http://setosa.io/
[2018040510520954]: /images/20220528/a2a71e9f55a24f138d3f3af7b8b81f64.png
[20180405105325496]: /images/20220528/8a5fb4c13c2b42b1958ac5eefa43b13b.png
[20180405105407816]: /images/20220528/f2304bf310724fa0b4eca2b6701cbb1b.png
[20180405105437656]: /images/20220528/52d550d232104829a294de37015fdf0b.png
[20180405105456796]: /images/20220528/cc8ce4819808452690e0a17c7b5e69d6.png
[20180405105521490]: /images/20220528/9cdf26efe4d940a585d0a82a627cd056.png
[20180405105627749]: /images/20220528/c3649f238be445c2954e1a430a7eb7dd.png
[20180405105653885]: /images/20220528/aa6bf420335a4fa9a43448643a27f993.png
[Link 1]: https://github.com/mpatacchiola/dissecting-reinforcement-learning/raw/master/src/1/T.npy
[20180405110330479]: /images/20220528/2ab360735d7b42788ce5f93f01c00b2a.png
[20180405110412178]: /images/20220528/b3d97a2629de4858a29634527d5fc286.png
[20180405110647313]: /images/20220528/141ad6a3f0f24aaf9e826907e8060d84.png
[20180405110457115]: /images/20220528/57205aa40a70433fbbcd26f7d2b958d5.png
[20180405112448717]: /images/20220528/8630524e01204175af8bf051b6c64a2e.png
[Link 2]: http://pubsonline.informs.org/doi/abs/10.1287/mnsc.24.11.1127
[20180405112519196]: /images/20220528/6847f8978eec4a0ea42469b1c186dff8.png
[Link 3]: https://blog.csdn.net/wilbertzhou/article/details/79822237
[Link 4]: https://blog.csdn.net/wilbertzhou/article/details/79845415
[Link 5]: https://blog.csdn.net/wilbertzhou/article/details/79953950
[Link 6]: https://blog.csdn.net/wilbertzhou/article/details/80041495
[Link 7]: https://blog.csdn.net/wilbertzhou/article/details/80113538
[Link 8]: https://blog.csdn.net/wilbertzhou/article/details/80187257
[Link 9]: https://blog.csdn.net/wilbertzhou/article/details/80260205
[Link 10]: https://blog.csdn.net/wilbertzhou/article/details/86743877
[github repository]: https://github.com/aimacode