Hands on Reinforcement Learning Basic Chapter

1 初探强化学习

1.1 简介

亲爱的读者，欢迎来到强化学习的世界。初探强化学习，你是否充满了好奇和期待呢？我们想说，首先感谢你的选择，学习本书不仅能够帮助你理解强化学习的算法原理，提高代码实践能力，更能让你了解自己是否喜欢决策智能这个方向，从而更好地决策未来是否从事人工智能方面的研究和实践工作。人生中充满选择，每次选择就是一次决策，我们正是从一次次决策中，把自己带领到人生的下一段旅程中。在回忆往事时，我们会对生命中某些时刻的决策印象深刻：“还好我当时选择了读博，我在那几年找到了自己的兴趣所在，现在我能做自己喜欢的工作！”“唉，当初我要是去那家公司实习就好了，在那里做的技术研究现在带来了巨大的社会价值。”通过这些反思，我们或许能领悟一些道理，变得更加睿智和成熟，以更积极的精神来迎接未来的选择和成长。

在机器学习领域，有一类重要的任务和人生选择很相似，即序贯决策（sequential decision making）任务。决策和预测任务不同，决策往往会带来“后果”，因此决策者需要为未来负责，在未来的时间点做出进一步的决策。实现序贯决策的机器学习方法就是本书讨论的主题—强化学习（reinforcement learning）。预测仅仅产生一个针对输入数据的信号，并期望它和未来可观测到的信号一致，这不会使未来情况发生任何改变。

本章主要讨论强化学习的基本概念和思维方式。希望通过本章的讨论，读者能了解强化学习在解决什么任务，其基本的数学刻画是什么样的，学习的目标是什么，以及它和预测型的有监督学习方法有什么根本性的区别。而关于如何设计强化学习算法，我们会在接下来的章节里面细细讨论。

1.2 什么是强化学习

广泛地讲，强化学习是机器通过与环境交互来实现目标的一种计算方法。机器和环境的一轮交互是指，机器在环境的一个状态下做一个动作决策，把这个动作作用到环境当中，这个环境发生相应的改变并且将相应的奖励反馈和下一轮状态传回机器。这种交互是迭代进行的，机器的目标是最大化在多轮交互过程中获得的累积奖励的期望。强化学习用智能体（agent）这个概念来表示做决策的机器。相比于有监督学习中的“模型”，强化学习中的“智能体”强调机器不但可以感知周围的环境信息，还可以通过做决策来直接改变这个环境，而不只是给出一些预测信号。

智能体和环境之间具体的交互方式如图所示。在每一轮交互中，智能体感知到环境目前所处的状态，经过自身的计算给出本轮的动作，将其作用到环境中；环境得到智能体的动作后，产生相应的即时奖励信号并发生相应的状态转移。智能体则在下一轮交互中感知到新的环境状态，依次类推。

这里，智能体有3种关键要素，即感知、决策和奖励。

• 感知。智能体在某种程度上感知环境的状态，从而知道自己所处的现状。例如，下围棋的智能体感知当前的棋盘情况；无人车感知周围道路的车辆、行人和红绿灯等情况；机器狗通过摄像头感知面前的图像，通过脚底的力学传感器来感知地面的摩擦功率和倾斜度等情况。

• 智能体根据当前的状态计算出达到目标需要采取的动作的过程叫作决策。例如，针对当前的棋盘决定下一颗落子的位置；针对当前的路况，无人车计算出方向盘的角度和刹车、油门的力度；针对当前收集到的视觉和力觉信号，机器狗给出4条腿的齿轮的角速度。策略是智能体最终体现出的智能形式，是不同智能体之间的核心区别。

• 奖励。环境根据状态和智能体采取的动作，产生一个标量信号作为奖励反馈。这个标量信号衡量智能体这一轮动作的好坏。例如，围棋博弈是否胜利；无人车是否安全、平稳且快速地行驶；机器狗是否在前进而没有摔倒。最大化累积奖励期望是智能体提升策略的目标，也是衡量智能体策略好坏的关键指标。

从以上分析可以看出，面向决策任务的强化学习和面向预测任务的有监督学习在形式上是有不少区别的。首先，决策任务往往涉及多轮交互，即序贯决策；而预测任务总是单轮的独立任务。如果决策也是单轮的，那么它可以转化为“判别最优动作”的预测任务。其次，因为决策任务是多轮的，智能体就需要在每轮做决策时考虑未来环境相应的改变，所以当前轮带来最大奖励反馈的动作，在长期来看并不一定是最优的。

1.3 强化学习的环境

我们从1.2节可以看到，强化学习的智能体是在和一个动态环境的交互中完成序贯决策的。我们说一个环境是动态的，意思就是它会随着某些因素的变化而不断演变，这在数学和物理中往往用随机过程来刻画。其实，生活中几乎所有的系统都在进行演变，例如一座城市的交通、一片湖中的生态、一场足球比赛、一个星系等。对于一个随机过程，其最关键的要素就是状态以及状态转移的条件概率分布。这就好比一个微粒在水中的布朗运动可以由它的起始位置以及下一刻的位置相对当前位置的条件概率分布来刻画。

如果在环境这样一个自身演变的随机过程中加入一个外来的干扰因素，即智能体的动作，那么环境的下一刻状态的概率分布将由当前状态和智能体的动作来共同决定，用最简单的数学公式表示则是

$$\begin{aligned} \text{下一个状态} &\sim P(\cdot | \text{当前状态, 智能体的动作})\\ s_{t+1} &\sim P(\cdot | s_t, a_t) \end{aligned}$$

根据上式可知，智能体决策的动作作用到环境中，使得环境发生相应的状态改变，而智能体接下来则需要在新的状态下进一步给出决策。

由此我们看到，与面向决策任务的智能体进行交互的环境是一个动态的随机过程，其未来状态的分布由当前状态和智能体决策的动作来共同决定，并且每一轮状态转移都伴随着两方面的随机性：一是智能体决策的动作的随机性，二是环境基于当前状态和智能体动作来采样下一刻状态的随机性。通过对环境的动态随机过程的刻画，我们能清楚地感受到，在动态随机过程中学习和在一个固定的数据分布下学习是非常不同的。

1.4 强化学习的目标

在上述动态环境下，智能体和环境每次进行交互时，环境会产生相应的奖励信号，其往往由实数标量来表示。这个奖励信号一般是诠释当前状态或动作的好坏的及时反馈信号，好比在玩游戏的过程中某一个操作获得的分数值。整个交互过程的每一轮获得的奖励信号可以进行累加，形成智能体的整体回报（return），好比一盘游戏最后的分数值。根据环境的动态性我们可以知道，即使环境和智能体策略不变，智能体的初始状态也不变，智能体和环境交互产生的结果也很可能是不同的，对应获得的回报也会不同。因此，在强化学习中，我们关注回报的期望，并将其定义为价值（value），这就是强化学习中智能体学习的优化目标。

价值的计算有些复杂，因为需要对交互过程中每一轮智能体采取动作的概率分布和环境相应的状态转移的概率分布做积分运算。强化学习和有监督学习的学习目标其实是一致的，即在某个数据分布下优化一个分数值的期望。不过，经过后面的分析我们会发现，强化学习和有监督学习的优化途径是不同的。

1.5 强化学习中的数据

接下来我们从数据层面谈谈有监督学习和强化学习的区别。

有监督学习的任务建立在从给定的数据分布中采样得到的训练数据集上，通过优化在训练数据集中设定的目标函数（如最小化预测误差）来找到模型的最优参数。这里，训练数据集背后的数据分布是完全不变的。

在强化学习中，数据是在智能体与环境交互的过程中得到的。如果智能体不采取某个决策动作，那么该动作对应的数据就永远无法被观测到，所以当前智能体的训练数据来自之前智能体的决策结果。因此，智能体的策略不同，与环境交互所产生的数据分布就不同，如图所示。

具体而言，强化学习中有一个关于数据分布的概念，叫作占用度量（occupancy measure），其具体的数学定义和性质会在第3章讨论，在这里我们只做简要的陈述：归一化的占用度量用于衡量在一个智能体决策与一个动态环境的交互过程中，采样到一个具体的状态动作对（state-action pair）的概率分布。

占用度量有一个很重要的性质：给定两个策略及其与一个动态环境交互得到的两个占用度量，那么当且仅当这两个占用度量相同时，这两个策略相同。也就是说，如果一个智能体的策略有所改变，那么它和环境交互得到的占用度量也会相应改变。

根据占用度量这一重要的性质，我们可以领悟到强化学习本质的思维方式。

• 强化学习的策略在训练中会不断更新，其对应的数据分布（即占用度量）也会相应地改变。因此，强化学习的一大难点就在于，智能体看到的数据分布是随着智能体的学习而不断发生改变的。

• 由于奖励建立在状态动作对之上，一个策略对应的价值其实就是一个占用度量下对应的奖励的期望，因此寻找最优策略对应着寻找最优占用度量。

1.6 强化学习的独特性

通过前面5节的讲解，读者现在应该已经对强化学习的基本数学概念有了一定的了解。这里我们回过头来再看看一般的有监督学习和强化学习的区别。

对于一般的有监督学习任务，我们的目标是找到一个最优的模型函数，使其在训练数据集上最小化一个给定的损失函数。在训练数据独立同分布的假设下，这个优化目标表示最小化模型在整个数据分布上的泛化误差（generalization error），用简要的公式可以概括为：

$$\text{最优模型} = \underset{\text{模型}}{\arg\min} \mathbb{E}_{\text{特征,标签} \sim \text{数据分布}} [\text{损失函数(标签, 模型(特征))}]$$

相比之下，强化学习任务的最终优化目标是最大化智能体策略在和动态环境交互过程中的价值。根据1.5节的分析，策略的价值可以等价转换成奖励函数在策略的占用度量上的期望，即：

$$\text{最优策略} = \underset{\text{策略}}{\arg\min} \mathbb{E}_{\text{状态,动作} \sim \text{策略的占用度量}} [\text{奖励函数(状态, 动作)}]$$

观察以上两个优化公式，我们可以回顾1.4节，总结出两者的相似点和不同点。

• 有监督学习和强化学习的优化目标相似，即都是在优化某个数据分布下的一个分数值的期望。

• 二者优化的途径是不同的，有监督学习直接通过优化模型对于数据特征的输出来优化目标，即修改目标函数而数据分布不变；强化学习则通过改变策略来调整智能体和环境交互数据的分布，进而优化目标，即修改数据分布而目标函数不变。

1.7 小结

本章通过简短的篇幅，大致介绍了强化学习的样貌，梳理了强化学习和有监督学习在范式以及思维方式上的相似点和不同点。在大多数情况下，强化学习任务往往比一般的有监督学习任务更难，因为一旦策略有所改变，其交互产生的数据分布也会随之改变，并且这样的改变是高度复杂、不可追踪的，往往不能用显式的数学公式刻画。这就好像一个混沌系统，我们无法得到其中一个初始设置对应的最终状态分布，而一般的有监督学习任务并没有这样的混沌效应。

好了，接下来该是我们躬身入局，通过理论学习和代码实践来学习强化学习的时候了。你准备好了吗？我们开始吧！

2 多臂老虎机

2.1 简介

我们在第 1 章中了解到，强化学习关注智能体和环境交互过程中的学习，这是一种试错型学习（trial-and-error learning）范式。在正式学习强化学习之前，我们需要先了解多臂老虎机问题，它可以被看作简化版的强化学习问题。与强化学习不同，多臂老虎机不存在状态信息，只有动作和奖励，算是最简单的“和环境交互中的学习”的一种形式。多臂老虎机中的探索与利用（exploration vs. exploitation）问题一直以来都是一个特别经典的问题，理解它能够帮助我们学习强化学习。

2.2 问题介绍

2.2.1 问题定义

在多臂老虎机（multi-armed bandit，MAB）问题（见图）中，有一个拥有 $K$ 根拉杆的老虎机，拉动每一根拉杆都对应一个关于奖励的概率分布 $\mathcal{R}$ 。我们每次拉动其中一根拉杆，就可以从该拉杆对应的奖励概率分布中获得一个奖励 $r$ 。我们在各根拉杆的奖励概率分布未知的情况下，从头开始尝试，目标是在操作 $T$ 次拉杆后获得尽可能高的累积奖励。由于奖励的概率分布是未知的，因此我们需要在“探索拉杆的获奖概率”和“根据经验选择获奖最多的拉杆”中进行权衡。“采用怎样的操作策略才能使获得的累积奖励最高”便是多臂老虎机问题。

2.2.2 形式化描述

多臂老虎机问题可以表示为一个元组 $<\mathcal{A,R}>$ ，其中：

• $\mathcal{A}$ 为动作集合，其中一个动作表示拉动一个拉杆。若多臂老虎机一共有 $K$ 根拉杆，那动作空间就是集合$\{a_1,a_2,\cdots,a_K\}$，我们用 $a_t \in \mathcal{A}$ 表示任意一个动作；
• $\mathcal{R}$ 为奖励概率分布，拉动每一根拉杆的动作 $a$ 都对应一个奖励概率分布 $\mathcal{R}(r|a)$ ，不同拉杆的奖励分布通常是不同的。

假设每个时间步只能拉动一个拉杆，多臂老虎机的目标为最大化一段时间步 $T$ 内累积的奖励: $\max \sum\limits_{t=1}^T r_t$，$r_t \sim \mathcal{R}(\cdot|a_t)$，其中 $a_t$ 表示在第 $t$ 时间步拉动某一拉杆的动作，$r_t$ 表示动作 $a_t$ 获得的奖励。

2.2.3 累积懊悔

对于每一个动作 $a_t$，我们定义其期望奖励为 $Q(a) = \mathbb{E}_{r\sim R(\cdot|a_t)}(r)$。于是，至少存在一根拉杆，它的期望奖励不小于拉动其他任意一根拉杆，我们将该最优期望奖励表示为 $Q^* = \max\limits_{a\in \mathcal{A}}Q(a)$。为了更加直观、方便地观察拉动一根拉杆的期望奖励离最优拉杆期望奖励的差距，我们引入懊悔（regret）概念。懊悔定义为拉动当前拉杆的动作与最优拉杆的期望奖励差，即 $R(a) = Q^* - Q(a)$。累积懊悔（cumulative regret）即操作 $T$ 次拉杆后累积的懊悔总量，对于一次完整的 $T$ 步决策$\{a_1,a_2,\cdots,a_T\}$，累积懊悔为 $\sigma_R=\sum_{t=1}^TR(a_t)$ 。MAB问题的目标为最大化累积奖励，等价于最小化累积懊悔。

2.2.4 估计期望奖励

为了知道拉动哪一根拉杆能获得更高的奖励，我们需要估计拉动这根拉杆的期望奖励。由于只拉动一次拉杆获得的奖励存在随机性，所以需要多次拉动一根拉杆，然后计算得到的多次奖励的期望，其算法流程如下所示:

• 对于 $\forall a \in \mathcal{A}$，初始化计数器 $N(a)=0$ 和期望奖励估值 $\hat{Q}(a)=0$

• for $t=1 \rightarrow T$ do

  • 选取某根拉杆，该动作记为 $a_t$

  • 得到奖励 $r_t$

  • 更新计数器: $N(a_t) = N(a_t) + 1$

  • 更新期望奖励估值：$\hat{Q}(a_t) = \hat{Q}(a_t) + \dfrac{1}{N(a_t)} [r_t - \hat{Q}(a_t)] $

• end for

以上 for 循环中的第四步如此更新估值，是因为这样可以进行增量式的期望更新，公式如下。

$$\begin{aligned} Q_k &= \dfrac{1}{k} \sum_{i=1}^k r_i \\ &= \dfrac{1}{k} \Big(r_k + \sum_{i=1}^{k-1} r_i\Big) \\ &= \dfrac{1}{k} \Big(r_k + (k-1)Q_{k-1}\Big) \\ &= \dfrac{1}{k} \Big(r_k + kQ_{k-1} - Q_{k-1}\Big) \\ &= Q_{k-1} + \dfrac{1}{k} (r_k - Q_{k-1}) \\ \end{aligned}$$

如果将所有数求和再除以次数，其缺点是每次更新的时间复杂度和空间复杂度均为 $O(n)$ 。而采用增量式更新，时间复杂度和空间复杂度均为 $O(1)$ 。

下面我们编写代码来实现一个拉杆数为 $T=10$ 的多臂老虎机。其中拉动每根拉杆的奖励服从伯努利分布（Bernoulli distribution），即每次拉下拉杆有 $p$ 的概率获得的奖励为 $1$，有 $1-p$ 的概率获得的奖励为 $0$。奖励为 $1$ 代表获奖，奖励为 $0$ 代表没有获奖。

## 导入需要使用的库,其中numpy是支持数组和矩阵运算的科学计算库,而matplotlib是绘图库
import numpy as np
import matplotlib.pyplot as plt


class BernoulliBandit:
    """ 伯努利多臂老虎机,输入K表示拉杆个数 """
    def __init__(self, K):
        self.probs = np.random.uniform(size=K)  ## 随机生成K个0～1的数,作为拉动每根拉杆的获奖
        ## 概率
        self.best_idx = np.argmax(self.probs)  ## 获奖概率最大的拉杆
        self.best_prob = self.probs[self.best_idx]  ## 最大的获奖概率
        self.K = K

    def step(self, k):
        ## 当玩家选择了k号拉杆后,根据拉动该老虎机的k号拉杆获得奖励的概率返回1（获奖）或0（未
        ## 获奖）
        if np.random.rand() < self.probs[k]:
            return 1
        else:
            return 0

np.random.seed(1)  ## 设定随机种子,使实验具有可重复性
K = 10
bandit_10_arm = BernoulliBandit(K)
print("随机生成了一个%d臂伯努利老虎机" % K)
print("获奖概率最大的拉杆为%d号,其获奖概率为%.4f" %
      (bandit_10_arm.best_idx, bandit_10_arm.best_prob))

随机生成了一个10臂伯努利老虎机
获奖概率最大的拉杆为1号,其获奖概率为0.7203

接下来我们用一个 Solver 基础类来实现上述的多臂老虎机的求解方案。根据前文的算法流程，我们需要实现下列函数功能：根据策略选择动作、根据动作获取奖励、更新期望奖励估值、更新累积懊悔和计数。在下面的 MAB 算法基本框架中，我们将根据策略选择动作、根据动作获取奖励和更新期望奖励估值放在 run_one_step() 函数中，由每个继承 Solver 类的策略具体实现。而更新累积懊悔和计数则直接放在主循环 run() 中。

class Solver:
    """ 多臂老虎机算法基本框架 """
    def __init__(self, bandit):
        self.bandit = bandit
        self.counts = np.zeros(self.bandit.K)  ## 每根拉杆的尝试次数
        self.regret = 0.  ## 当前步的累积懊悔
        self.actions = []  ## 维护一个列表,记录每一步的动作
        self.regrets = []  ## 维护一个列表,记录每一步的累积懊悔

    def update_regret(self, k):
        ## 计算累积懊悔并保存,k为本次动作选择的拉杆的编号
        self.regret += self.bandit.best_prob - self.bandit.probs[k]
        self.regrets.append(self.regret)

    def run_one_step(self):
        ## 返回当前动作选择哪一根拉杆,由每个具体的策略实现
        raise NotImplementedError

    def run(self, num_steps):
        ## 运行一定次数,num_steps为总运行次数
        for _ in range(num_steps):
            k = self.run_one_step()
            self.counts[k] += 1
            self.actions.append(k)
            self.update_regret(k)

2.3 探索与利用的平衡

在上述算法框架中，还没有一个策略告诉我们应该采取哪个动作，即拉动哪根拉杆，所以接下来我们将学习如何设计一个策略。例如，一个最简单的策略就是一直采取第一个动作，但这就非常依赖运气的好坏。如果运气绝佳，可能拉动的刚好是能获得最大期望奖励的拉杆，即最优拉杆；但如果运气很糟糕，获得的就有可能是最小的期望奖励。在多臂老虎机问题中，一个经典的问题就是探索与利用的平衡问题。探索（exploration）是指尝试拉动更多可能的拉杆，这根拉杆不一定会获得最大的奖励，但这种方案能够摸清楚所有拉杆的获奖情况。例如，对于一个 $10$ 臂老虎机，我们要把所有的拉杆都拉动一下才知道哪根拉杆可能获得最大的奖励。利用（exploitation）是指拉动已知期望奖励最大的那根拉杆，由于已知的信息仅仅来自有限次的交互观测，所以当前的最优拉杆不一定是全局最优的。例如，对于一个 $10$ 臂老虎机，我们只拉动过其中 $3$ 根拉杆，接下来就一直拉动这 $3$ 根拉杆中期望奖励最大的那根拉杆，但很有可能期望奖励最大的拉杆在剩下的 $7$ 根当中，即使我们对 $10$ 根拉杆各自都尝试了 $20$ 次，发现 $5$ 号拉杆的经验期望奖励是最高的，但仍然存在着微小的概率，即另一根 $6$ 号拉杆的真实期望奖励是比 $5$ 号拉杆更高的。

于是在多臂老虎机问题中，设计策略时就需要平衡探索和利用的次数，使得累积奖励最大化。一个比较常用的思路是在开始时做比较多的探索，在对每根拉杆都有比较准确的估计后，再进行利用。目前已有一些比较经典的算法来解决这个问题，例如 ϵ-贪婪算法、上置信界算法 和 汤普森采样算法 等，我们接下来将分别介绍这几种算法。

2.4 ϵ-贪心算法

完全贪婪算法即在每一时刻采取期望奖励估值最大的动作（拉动拉杆），这就是纯粹的利用，而没有探索，所以我们通常需要对完全贪婪算法进行一些修改，其中比较经典的一种方法为 ϵ-贪婪（ϵ-Greedy）算法。

ϵ-贪婪算法在完全贪婪算法的基础上添加了噪声，每次以概率 $1-\epsilon$ 选择以往经验中期望奖励估值最大的那根拉杆（利用），以概率 $\epsilon$ 随机选择一根拉杆（探索），公式如下：

$$a_t = \begin{cases} \underset{a \in \mathcal{A}}{\arg\max} \hat{Q}(x), & \text{采样概率}: 1-\epsilon\\ \text{从} \mathcal{A} \text{中随机选择}, & \text{采样概率}: \epsilon \end{cases}$$

随着探索次数的不断增加，我们对各个动作的奖励估计得越来越准，此时我们就没必要继续花大力气进行探索。所以在 ϵ-贪婪 算法的具体实现中，我们可以令 ϵ 随时间衰减，即探索的概率将会不断降低。但是请注意， ϵ 不会在有限的步数内衰减至 0，因为基于有限步数观测的完全贪婪算法仍然是一个局部信息的贪婪算法，永远距离最优解有一个固定的差距。

我们接下来编写代码来实现一个 ϵ-贪婪 算法，并用它去解决之前生成的 $10$ 臂老虎机的问题。设置 $\epsilon=0.01$，以及 $T=5000$。

class EpsilonGreedy(Solver):
    """ epsilon贪婪算法,继承Solver类 """
    def __init__(self, bandit, epsilon=0.01, init_prob=1.0):
        super(EpsilonGreedy, self).__init__(bandit)
        self.epsilon = epsilon
        #初始化拉动所有拉杆的期望奖励估值
        self.estimates = np.array([init_prob] * self.bandit.K)

    def run_one_step(self):
        if np.random.random() < self.epsilon:
            k = np.random.randint(0, self.bandit.K)  ## 随机选择一根拉杆
        else:
            k = np.argmax(self.estimates)  ## 选择期望奖励估值最大的拉杆
        r = self.bandit.step(k)  ## 得到本次动作的奖励
        self.estimates[k] += 1. / (self.counts[k] + 1) * (r - self.estimates[k])
        return k

为了更加直观地展示，可以把每一时间步的累积函数绘制出来。于是我们定义了以下绘图函数，方便之后调用。

def plot_results(solvers, solver_names):
    """生成累积懊悔随时间变化的图像。输入solvers是一个列表,列表中的每个元素是一种特定的策略。
    而solver_names也是一个列表,存储每个策略的名称"""
    for idx, solver in enumerate(solvers):
        time_list = range(len(solver.regrets))
        plt.plot(time_list, solver.regrets, label=solver_names[idx])
    plt.xlabel('Time steps')
    plt.ylabel('Cumulative regrets')
    plt.title('%d-armed bandit' % solvers[0].bandit.K)
    plt.legend()
    plt.show()


np.random.seed(1)
epsilon_greedy_solver = EpsilonGreedy(bandit_10_arm, epsilon=0.01)
epsilon_greedy_solver.run(5000)
print('epsilon-贪婪算法的累积懊悔为：', epsilon_greedy_solver.regret)
plot_results([epsilon_greedy_solver], ["EpsilonGreedy"])

psilon-贪婪算法的累积懊悔为：25.526630933945313

通过上面的实验可以发现，在经历了开始的一小段时间后，ϵ-贪婪算法的累积懊悔几乎是线性增长的。这是 $\epsilon=0.01$ 时的结果，因为一旦做出了随机拉杆的探索，那么产生的懊悔值是固定的。其他不同的 ϵ 取值又会带来怎样的变化呢？我们继续使用该 $10$ 臂老虎机，我们尝试不同的参数 $\{10^{-4}, 0.01, 0.1, 0.25, 0.5\}$，查看相应的实验结果。

np.random.seed(0)
epsilons = [1e-4, 0.01, 0.1, 0.25, 0.5]
epsilon_greedy_solver_list = [
    EpsilonGreedy(bandit_10_arm, epsilon=e) for e in epsilons
]
epsilon_greedy_solver_names = ["epsilon={}".format(e) for e in epsilons]
for solver in epsilon_greedy_solver_list:
    solver.run(5000)

plot_results(epsilon_greedy_solver_list, epsilon_greedy_solver_names)

通过实验结果可以发现，基本上无论 ϵ 取值多少，累积懊悔都是线性增长的。在这个例子中，随着 ϵ 的增大，累积懊悔增长的速率也会增大。 接下来我们尝试 ϵ 值随时间衰减的 ϵ-贪婪 算法，采取的具体衰减形式为反比例衰减，公式为 $\epsilon_t = \dfrac{1}{t}$

class DecayingEpsilonGreedy(Solver):
    """ epsilon值随时间衰减的epsilon-贪婪算法,继承Solver类 """
    def __init__(self, bandit, init_prob=1.0):
        super(DecayingEpsilonGreedy, self).__init__(bandit)
        self.estimates = np.array([init_prob] * self.bandit.K)
        self.total_count = 0

    def run_one_step(self):
        self.total_count += 1
        if np.random.random() < 1 / self.total_count:  ## epsilon值随时间衰减
            k = np.random.randint(0, self.bandit.K)
        else:
            k = np.argmax(self.estimates)

        r = self.bandit.step(k)
        self.estimates[k] += 1. / (self.counts[k] + 1) * (r - self.estimates[k])

        return k


np.random.seed(1)
decaying_epsilon_greedy_solver = DecayingEpsilonGreedy(bandit_10_arm)
decaying_epsilon_greedy_solver.run(5000)
print('epsilon值衰减的贪婪算法的累积懊悔为：', decaying_epsilon_greedy_solver.regret)
plot_results([decaying_epsilon_greedy_solver], ["DecayingEpsilonGreedy"])

epsilon值衰减的贪婪算法的累积懊悔为：10.114334931260183

从实验结果图中可以发现，随时间做反比例衰减的 ϵ-贪婪 算法能够使累积懊悔与时间步的关系变成次线性（sublinear）的，这明显优于固定 ϵ 值的 ϵ-贪婪算法。

2.5 上置信界算法

设想这样一种情况：对于一台双臂老虎机，其中第一根拉杆只被拉动过一次，得到的奖励为 $0$；第二根拉杆被拉动过很多次，我们对它的奖励分布已经有了大致的把握。这时你会怎么做？或许你会进一步尝试拉动第一根拉杆，从而更加确定其奖励分布。这种思路主要是基于不确定性，因为此时第一根拉杆只被拉动过一次，它的不确定性很高。一根拉杆的不确定性越大，它就越具有探索的价值，因为探索之后我们可能发现它的期望奖励很大。我们在此引入不确定性度量 $U(a)$，它会随着一个动作被尝试次数的增加而减小。我们可以使用一种基于不确定性的策略来综合考虑现有的期望奖励估值和不确定性，其核心问题是如何估计不确定性。

上置信界（upper confidence bound，UCB）算法是一种经典的基于不确定性的策略算法，它的思想用到了一个非常著名的数学原理：霍夫丁不等式（Hoeffding’s inequality）。在霍夫丁不等式中，令 $X_1,\cdots,X_n$ 为 $n$ 个独立同分布的随机变量，取值范围为 $[0,1]$ ，其经验期望为 $\bar{x}_n = \dfrac{1}{n}\sum\limits_{j=1}^n x_j$，则有

$$\mathbb{P}\{ \mathbb{E} (x) \ge \bar{x}_n + u \} \le e^{-2nu^2}$$

现在我们将霍夫丁不等式运用于多臂老虎机问题中。将 $\hat{Q}_t(a)$ 代入 $\bar{x}_t$，不等式中的参数 $u=\hat{U}_t(a)$ 代表不确定性度量。给定一个概率 $p=e^{-2N_t(a)U_t(a)^2}$，根据上述不等式，$Q_t(a) < \hat{Q}_t(a) + \hat{U}_t(a)$ 至少以概率 $1-p$ 成立。当 $p$ 很小时，$Q_t(a) < \hat{Q}_t(a) + \hat{U}_t(a)$ 就以很大概率成立，$\hat{Q}_t(a) + \hat{U}_t(a)$ 便是期望奖励上界。此时，上置信界算法便选取期望奖励上界最大的动作，即 $a = \underset{a\in\mathcal{A}}{\arg\max}[\hat{Q}(a) + \hat{U}(a)]$。那其中 $\hat{U}_t(a)$ 具体是什么呢？根据等式 $e^{-2N_t(a)U_t(a)^2}$，解之即得 $\hat{U}_t(a) = \sqrt{\dfrac{-\log p}{2N_t(a)}}$。因此，设定一个概率 $p$ 后，就可以计算相应的不确定性度量 $\hat{U}_t(a)$ 了。更直观地说，UCB 算法在每次选择拉杆前，先估计每根拉杆的期望奖励的上界，使得拉动每根拉杆的期望奖励只有一个较小的概率 $p$ 超过这个上界，接着选出期望奖励上界最大的拉杆，从而选择最有可能获得最大期望奖励的拉杆。

我们编写代码来实现 UCB 算法，并且仍然使用ϵ-贪心算法中定义的$10$臂老虎机来观察实验结果。在具体的实现过程中，设置 $p=\dfrac{1}{t}$，并且在分母中为拉动每根拉杆的次数加上常数 1，以免出现分母为 0 的情形，即此时 $\hat{U}_t(a) = \sqrt{\dfrac{\log t}{2(N_t(a) + 1)}}$ 。同时，我们设定一个系数 $c$ 来控制不确定性的比重，此时 $a = \underset{a\in\mathcal{A}}{\arg\max} \hat{Q}(a) + c \cdot \hat{U}(a)$。

class UCB(Solver):
    """ UCB算法,继承Solver类 """
    def __init__(self, bandit, coef, init_prob=1.0):
        super(UCB, self).__init__(bandit)
        self.total_count = 0
        self.estimates = np.array([init_prob] * self.bandit.K)
        self.coef = coef

    def run_one_step(self):
        self.total_count += 1
        ucb = self.estimates + self.coef * np.sqrt(
            np.log(self.total_count) / (2 * (self.counts + 1)))  ## 计算上置信界
        k = np.argmax(ucb)  ## 选出上置信界最大的拉杆
        r = self.bandit.step(k)
        self.estimates[k] += 1. / (self.counts[k] + 1) * (r - self.estimates[k])
        return k


np.random.seed(1)
coef = 1  ## 控制不确定性比重的系数
UCB_solver = UCB(bandit_10_arm, coef)
UCB_solver.run(5000)
print('上置信界算法的累积懊悔为：', UCB_solver.regret)
plot_results([UCB_solver], ["UCB"])

2.6 汤普森采样算法

MAB 中还有一种经典算法——汤普森采样（Thompson sampling），先假设拉动每根拉杆的奖励服从一个特定的概率分布，然后根据拉动每根拉杆的期望奖励来进行选择。但是由于计算所有拉杆的期望奖励的代价比较高，汤普森采样算法使用采样的方式，即根据当前每个动作 $a$ 的奖励概率分布进行一轮采样，得到一组各根拉杆的奖励样本，再选择样本中奖励最大的动作。可以看出，汤普森采样是一种计算所有拉杆的最高奖励概率的蒙特卡洛采样方法。

了解了汤普森采样算法的基本思路后，我们需要解决另一个问题：怎样得到当前每个动作 $a$ 的奖励概率分布并且在过程中进行更新？在实际情况中，我们通常用 Beta 分布对当前每个动作的奖励概率分布进行建模。具体来说，若某拉杆被选择了 $k$ 次，其中 $m_1$ 次奖励为 $1$，$m_2$ 次奖励为 $0$，则该拉杆的奖励服从参数为 $(m_1 + 1, m_2 + 1)$ 的 Beta 分布。下图是汤普森采样的一个示例。

我们编写代码来实现汤普森采样算法，并且仍然使用先前定义的$10$臂老虎机来观察实验结果。

class ThompsonSampling(Solver):
    """ 汤普森采样算法,继承Solver类 """
    def __init__(self, bandit):
        super(ThompsonSampling, self).__init__(bandit)
        self._a = np.ones(self.bandit.K)  ## 列表,表示每根拉杆奖励为1的次数
        self._b = np.ones(self.bandit.K)  ## 列表,表示每根拉杆奖励为0的次数

    def run_one_step(self):
        samples = np.random.beta(self._a, self._b)  ## 按照Beta分布采样一组奖励样本
        k = np.argmax(samples)  ## 选出采样奖励最大的拉杆
        r = self.bandit.step(k)

        self._a[k] += r  ## 更新Beta分布的第一个参数
        self._b[k] += (1 - r)  ## 更新Beta分布的第二个参数
        return k


np.random.seed(1)
thompson_sampling_solver = ThompsonSampling(bandit_10_arm)
thompson_sampling_solver.run(5000)
print('汤普森采样算法的累积懊悔为：', thompson_sampling_solver.regret)
plot_results([thompson_sampling_solver], ["ThompsonSampling"])

汤普森采样算法的累积懊悔为：57.19161964443925

通过实验我们可以得到以下结论：ϵ-贪心算法的累积懊悔是随时间线性增长的，而另外3种算法（ϵ-衰减贪婪算法、上置信界算法、汤普森采样算法）的累积懊悔都是随时间次线性增长的（具体为对数形式增长）。

2.7 总结

探索与利用是与环境做交互学习的重要问题，是强化学习试错法中的必备技术，而多臂老虎机问题是研究探索与利用技术理论的最佳环境。了解多臂老虎机的探索与利用问题，对接下来我们学习强化学习环境探索有很重要的帮助。对于多臂老虎机各种算法的累积懊悔理论分析，有兴趣的同学可以自行查阅相关资料。ϵ-贪婪算法、上置信界算法和汤普森采样算法在多臂老虎机问题中十分常用，其中上置信界算法和汤普森采样方法均能保证对数的渐进最优累积懊悔。

多臂老虎机问题与强化学习的一大区别在于其与环境的交互并不会改变环境，即多臂老虎机的每次交互的结果和以往的动作无关，所以可看作无状态的强化学习（stateless reinforcement learning）。第 3 章将开始在有状态的环境下讨论强化学习，即马尔可夫决策过程。

2.8 参考文献

[1] AUER P, CESA-BIANCHI N,FISCHER P. Finite-time analysis of the multiarmed bandit problem[J]. Machine learning,2002, 47 (2) : 235-256.

[2] AUER P. Using confidence bounds for exploitation-exploration trade-offs[J]. Journal of Machine Learning Research,2002, 3(3): 397-422.

[3] GITTINS J, GLAZEBROOK K,WEBER R. Multi-armed bandit allocation indices[M]. 2nd ed. America: John Wiley & Sons, 2011.

[4] CHAPELLE O, LI L. An empirical evaluation of thompson sampling [J]. Advances in neural information processing systems, 2011, 24: 2249-2257.

3 马尔可夫决策过程

3.1 简介

马尔可夫决策过程（Markov decision process，MDP）是强化学习的重要概念。要学好强化学习，我们首先要掌握马尔可夫决策过程的基础知识。前两章所说的强化学习中的环境一般就是一个马尔可夫决策过程。与多臂老虎机问题不同，马尔可夫决策过程包含状态信息以及状态之间的转移机制。如果要用强化学习去解决一个实际问题，第一步要做的事情就是把这个实际问题抽象为一个马尔可夫决策过程，也就是明确马尔可夫决策过程的各个组成要素。本章将从马尔可夫过程出发，一步一步地进行介绍，最后引出马尔可夫决策过程。

3.2 马尔可夫过程

3.2.1 随机过程

随机过程（stochastic process）是概率论的“动力学”部分。概率论的研究对象是静态的随机现象，而随机过程的研究对象是随时间演变的随机现象（例如天气随时间的变化、城市交通随时间的变化）。在随机过程中，随机现象在某时刻 $t$ 的取值是一个向量随机变量，用 $S_t$ 表示，所有可能的状态组成状态集合 $\mathcal{S}$ 。随机现象便是状态的变化过程。在某时刻 $t$ 的状态 $S_t$ 通常取决于 $t$ 时刻之前的状态。我们将已知历史信息 $(S_1,\cdots,S_t)$ 时下一个时刻状态为 $S_{t+1}$ 的概率表示成 $P(S_{t+1}|S_1,\cdots,S_t)$ 。

3.2.2 马尔可夫性质

当且仅当某时刻的状态只取决于上一时刻的状态时，一个随机过程被称为具有马尔可夫性质（Markov property），用公式表示为 $P(S_{t+1}|S_t) = P(S_{t+1}|S_1,\cdots,S_t)$。也就是说，当前状态是未来的充分统计量，即下一个状态只取决于当前状态，而不会受到过去状态的影响。需要明确的是，具有马尔可夫性并不代表这个随机过程就和历史完全没有关系。因为虽然 $t+1$ 时刻的状态只与 $t$ 时刻的状态有关，但是 $t$ 时刻的状态其实包含了 $t-1$ 时刻的状态的信息，通过这种链式的关系，历史的信息被传递到了现在。马尔可夫性可以大大简化运算，因为只要当前状态可知，所有的历史信息都不再需要了，利用当前状态信息就可以决定未来。

3.2.3 马尔可夫过程

马尔可夫过程（Markov process）指具有马尔可夫性质的随机过程，也被称为马尔可夫链（Markov chain）。我们通常用元组 $(\mathcal{S,P})$ 描述一个马尔可夫过程，其中 $\mathcal{S}$ 是有限数量的状态集合， $\mathcal{P}$ 是状态转移矩阵（state transition matrix）。假设一共有 $n$ 个状态，此时 $\mathcal{S} = \{s_1, s_2, \cdots, s_n\}$ 。状态转移矩阵 $\mathcal{P}$ 定义了所有状态对之间的转移概率，即

$$\mathcal{P} = \begin{bmatrix} P(s_1|s_1) & \cdots & P(s_n|s_1)\\ \vdots & \ddots & \vdots\\ P(s_1|s_n) & \cdots & P(s_n|s_n)\\ \end{bmatrix}$$

矩阵 $\mathcal{P}$ 中第 $i$ 行第 $j$ 列元素 $P(s_j | s_i)$ 表示从状态 $s_i$ 转移到状态 $s_j$ 的概率，我们称 $P(s'|s)$ 为状态转移函数。从某个状态出发，到达其他状态的概率和必须为 $1$，即状态转移矩阵的每一行的和为 $1$。

下图是一个具有 $6$ 个状态的马尔可夫过程的简单例子。其中每个绿色圆圈表示一个状态，每个状态都有一定概率（包括概率为 $0$）转移到其他状态，其中 $s_6$ 通常被称为终止状态（terminal state），因为它不会再转移到其他状态，可以理解为它永远以概率 $1$ 转移到自己。状态之间的虚线箭头表示状态的转移，箭头旁的数字表示该状态转移发生的概率。从每个状态出发转移到其他状态的概率总和为 $1$。例如，$s_1$ 有 $90\%$ 概率保持不变，有 $10\%$ 概率转移到 $s_2$ ，而在 $s_2$ 又有 $50\%$ 概率回到 $s_1$ ，有 $50\%$ 概率转移到 $s_3$ 。

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

$$\mathcal{P} = \begin{bmatrix} 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 \\ \end{bmatrix}$$

其中第 $i$ 行 $j$ 列的值 $\mathcal{P}_{i,j}$ 则代表从状态 $s_i$ 转移到 $s_j$ 的概率。

给定一个马尔可夫过程，我们就可以从某个状态出发，根据它的状态转移矩阵生成一个状态序列（episode），这个步骤也被叫做采样（sampling）。例如，从 $s_1$ 出发，可以生成序列 $s_1 \rightarrow s_2 \rightarrow s_3 \rightarrow s_6$ 或序列 $s_1 \rightarrow s_1 \rightarrow s_2 \rightarrow s_3 \rightarrow s_4 \rightarrow s_5 \rightarrow s_3 \rightarrow s_6$ 等。生成这些序列的概率和状态转移矩阵有关。

3.3 马尔可夫奖励过程

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

• $\mathcal{S}$ 是有限状态的集合。

• $\mathcal{P}$ 是状态转移矩阵。

• $\mathcal{r}$ 是奖励函数，某个状态 $s$ 的奖励 $r(s)$ 指转移到该状态时可以获得奖励的期望。

• $\gamma$ 是折扣因子（discount factor），$\gamma$ 的取值范围为 $[0,1)$。引入折扣因子的理由为远期利益具有一定不确定性，有时我们更希望能够尽快获得一些奖励，所以我们需要对远期利益打一些折扣。接近 $1$ 的 $\gamma$ 更关注长期的累计奖励，接近 $0$ 的 $\gamma$ 更考虑短期奖励。

3.3.1 回报

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

$$G_t = R_t + \gamma R_{t+1} + \gamma^2 R_{t+2} + \cdots = \sum_{k=1}^\infin \gamma^k R_{t+k}$$

其中，$R_t$ 表示在时刻 $t$ 获得的奖励。继续沿用之前的马尔可夫过程的例子，并在其基础上添加奖励函数，构建成一个马尔可夫奖励过程。例如，进入状态 $s_2$ 可以得到奖励 $-2$ ，表明我们不希望进入 $s_2$ ，进入 $s_4$ 可以获得最高的奖励 $10$，但是进入 $s_6$ 之后奖励为零，并且此时序列也终止了。

比如选取 $s_1$ 为起始状态，设置 $\gamma = 0.5$ ，采样到一条状态序列为 $s_1 \rightarrow s_2 \rightarrow s_3 \rightarrow s_6$ ，就可以计算 $s_1$ 的回报 $G_1$，得到 $G_1 = -1 + 0.5 \times (-2) + 0.5^2 \times (-2) = -2.5$。

接下来我们用代码表示上图中的马尔可夫奖励过程，并且定义计算回报的函数。

import numpy as np
np.random.seed(0)
## 定义状态转移概率矩阵P
P = [
    [0.9, 0.1, 0.0, 0.0, 0.0, 0.0],
    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4],
    [0.0, 0.0, 0.0, 0.0, 0.3, 0.7],
    [0.0, 0.2, 0.3, 0.5, 0.0, 0.0],
    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
]
P = np.array(P)

rewards = [-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 + rewards[chain[i] - 1]
    return G


## 一个状态序列,s1-s2-s3-s6
chain = [1, 2, 3, 6]
start_index = 0
G = compute_return(start_index, chain, gamma)
print("根据本序列计算得到回报为：%s。" % G)

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

3.3.2 价值函数

在马尔可夫奖励过程中，一个状态的期望回报（即从这个状态出发的未来累积奖励的期望）被称为这个状态的价值（value）。所有状态的价值就组成了价值函数（value function），价值函数的输入为某个状态，输出为这个状态的价值。我们将价值函数写成 $V(s)=\mathbb{E}[G_t|S_t=s]$ ，展开为

$$\begin{aligned} V(s) &= \mathbb{E}[G_t|S_t=s]\\ &= \mathbb{E}[R_t + \gamma R_{t+1} + \gamma^2 R_{t+2} + \cdots|S_t=s]\\ &= \mathbb{E}[R_t + \gamma (R_{t+1} + \gamma R_{t+2} + \cdots)|S_t=s]\\ &= \mathbb{E}[R_t + \gamma G_{t+1}|S_t=s]\\ &= \mathbb{E}[R_t + \gamma V(S_{t+1})|S_t=s]\\ \end{aligned}$$

在上式的最后一个等号中，一方面，即时奖励的期望正是奖励函数的输出，即 $\mathbb{E}[R_t|S_t=s]=r(s)$；另一方面，等式中剩余部分 $\mathbb{E}[\gamma V(S_{t+1})|S_t=s]$ 可以根据从状态 $s$ 出发的转移概率得到，即可以得到

$$V(s) = r(s) + \gamma \sum_{s' \in \mathcal{S}} p(s'|s) V(s')$$

上式就是马尔可夫奖励过程中非常有名的贝尔曼方程（Bellman equation），对每一个状态都成立。若一个马尔可夫奖励过程一共有 $n$ 个状态，即 $\mathcal{S} = \{s_1, s_2, \cdots, s_n\}$，我们将所有状态的价值表示成一个列向量 $\mathcal{V} = [V(s_1),V(s_2),\cdots,V(s_n)]^T$，同理，将奖励函数写成一个列向量 $\mathcal{R} = [r(s_1),r(s_2),\cdots,r(s_n),]^T$。于是我们可以将贝尔曼方程写成矩阵的形式：

$$\begin{aligned} \mathcal{V}&=\mathcal{R+\gamma PV}\\ \begin{bmatrix*} V(s_1)\\ V(s_2)\\ \cdots \\ V(s_n) \end{bmatrix*} &= \begin{bmatrix*} r(s_1)\\ r(s_2)\\ \cdots \\ r(s_n) \end{bmatrix*} + \gamma \begin{bmatrix} P(s_1|s_1) & P(s_2|s_1) & \cdots & P(s_n|s_1)\\ P(s_1|s_2) & P(s_2|s_2) & \cdots & P(s_n|s_2)\\ \cdots & \cdots & \cdots & \cdots \\ P(s_1|s_n) & P(s_2|s_n) & \cdots & P(s_n|s_n)\\ \end{bmatrix} \begin{bmatrix*} V(s_1)\\ V(s_2)\\ \cdots \\ V(s_n) \end{bmatrix*} \end{aligned}$$

我们可以直接根据矩阵运算求解，得到以下解析解：

$$\begin{aligned} \mathcal{V}&=\mathcal{R+\gamma PV}\\ (I - \gamma \mathcal{P})\mathcal{V}&=\mathcal{R}\\ \mathcal{V}&=(I - \gamma \mathcal{P})^{-1}\mathcal{R}\\ \end{aligned}$$

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

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

def compute(P, rewards, gamma, states_num):
    ''' 利用贝尔曼方程的矩阵形式计算解析解,states_num是MRP的状态数 '''
    rewards = np.array(rewards).reshape((-1, 1))  #将rewards写成列向量形式
    value = np.dot(np.linalg.inv(np.eye(states_num, states_num) - gamma * P),
                   rewards)
    return value


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

MRP中每个状态价值分别为
 [[-2.01950168]
 [-2.21451846]
 [ 1.16142785]
 [10.53809283]
 [ 3.58728554]
 [ 0.        ]]

根据以上代码，求解得到各个状态的价值 $V(s)$，具体如下：

$$\begin{bmatrix*} V(s_1)\\ V(s_2)\\ V(s_3)\\ V(s_4)\\ V(s_5) \\ V(s_6) \end{bmatrix*} = \begin{bmatrix*} -2.02\\-2.21\\1.16\\10.54\\3.59\\0 \end{bmatrix*}$$

我们现在用贝尔曼方程来进行简单的验证。例如，对于状态 $s_4$ 来说，当 $\gamma=0.5$ 时，有

$$\begin{aligned} V(s_4) &= r(s_4) + \gamma \sum_{s'\in\mathcal{S}}P(s'|s_4)V(s')\\ 10.54 &= 10 + 0.5 \times (0.7 \times 0 + 0.3 \times 3.59)\\ \end{aligned}$$

可以发现左右两边的值几乎是相等的，说明我们求解得到的价值函数是满足状态为 $s_4$ 时的贝尔曼方程。读者可以自行验证在其他状态时贝尔曼方程是否也成立。若贝尔曼方程对于所有状态都成立，就可以说明我们求解得到的价值函数是正确的。除了使用动态规划算法，马尔可夫奖励过程中的价值函数也可以通过蒙特卡洛方法估计得到，我们将在 3.5 节中介绍该方法。

3.4 马尔可夫决策过程

3.2 节和 3.3 节讨论到的马尔可夫过程和马尔可夫奖励过程都是自发改变的随机过程；而如果有一个外界的“刺激”来共同改变这个随机过程，就有了马尔可夫决策过程（Markov decision process，MDP）。我们将这个来自外界的刺激称为智能体（agent）的动作，在马尔可夫奖励过程（MRP）的基础上加入动作，就得到了马尔可夫决策过程（MDP）。马尔可夫决策过程由元组 $(\mathcal{S,A,P,r,\gamma})$ 构成，其中：

• $\mathcal{S}$ 是状态的集合;

• $\mathcal{A}$ 是动作的集合;

• $\gamma$ 是折扣因子;

• $\mathcal{r(s,a)}$ 是奖励函数，此时奖励可以同时取决于状态 $s$ 和动作 $a$，在奖励函数只取决于状态 $s$ 时，则退化为 $r(s)$；

• $\mathcal{P(s'|s,a)}$ 是状态转移函数，表示在状态 $s$ 执行动作 $a$ 之后到达状态 $s'$ 的概率。

我们发现 MDP 与 MRP 非常相像，主要区别为 MDP 中的状态转移函数和奖励函数都比 MRP 多了动作 $a$ 作为自变量。注意，在上面 MDP 的定义中，我们不再使用类似 MRP 定义中的状态转移矩阵方式，而是直接表示成了状态转移函数。这样做一是因为此时状态转移与动作也有关，变成了一个三维数组，而不再是一个矩阵（二维数组）；二是因为状态转移函数更具有一般意义，例如，如果状态集合不是有限的，就无法用数组表示，但仍然可以用状态转移函数表示。我们在之后的课程学习中会遇到连续状态的 MDP 环境，那时状态集合都不是有限的。现在我们主要关注于离散状态的 MDP 环境，此时状态集合是有限的。

不同于马尔可夫奖励过程，在马尔可夫决策过程中，通常存在一个智能体来执行动作。例如，一艘小船在大海中随着水流自由飘荡的过程就是一个马尔可夫奖励过程，它如果凭借运气漂到了一个目的地，就能获得比较大的奖励；如果有个水手在控制着这条船往哪个方向前进，就可以主动选择前往目的地获得比较大的奖励。马尔可夫决策过程是一个与时间相关的不断进行的过程，在智能体和环境 MDP 之间存在一个不断交互的过程。一般而言，它们之间的交互是如下图的循环过程：智能体根据当前状态 $S_t$ 选择动作 $A_t$；对于状态 $S_t$ 和动作 $A_t$，MDP 根据奖励函数和状态转移函数得到 $S_{t+1}$ 和 $R_t$ 并反馈给智能体。智能体的目标是最大化得到的累计奖励。智能体根据当前状态从动作的集合 $\mathcal{A}$ 中选择一个动作的函数，被称为策略。

3.4.1 策略

智能体的策略（Policy）通常用字母 $\pi$ 表示。策略 $\pi(a|s) = P(A_t=a|S_t=s)$ 是一个函数，表示在输入状态情况下采取动作的概率。当一个策略是确定性策略（deterministic policy）时，它在每个状态时只输出一个确定性的动作，即只有该动作的概率为 $1$，其他动作的概率为 $0$；当一个策略是随机性策略（stochastic policy）时，它在每个状态时输出的是关于动作的概率分布，然后根据该分布进行采样就可以得到一个动作。在 MDP 中，由于马尔可夫性质的存在，策略只需要与当前状态有关，不需要考虑历史状态。回顾一下在 MRP 中的价值函数，在 MDP 中也同样可以定义类似的价值函数。但此时的价值函数与策略有关，这意味着对于两个不同的策略来说，它们在同一个状态下的价值也很可能是不同的。这很好理解，因为不同的策略会采取不同的动作，从而之后会遇到不同的状态，以及获得不同的奖励，所以它们的累积奖励的期望也就不同，即状态价值不同。

3.4.2 状态价值函数

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

$$V^{\pi}(s) = \mathbb{E}_{\pi}[G_t | S_t = s]$$

3.4.3 动作价值函数

不同于 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]$$

状态价值函数和动作价值函数之间的关系：在使用策略 $\pi$ 中，状态 $s$ 的价值等于在该状态下基于策略 $\pi$ 采取所有动作的概率与相应的价值相乘再求和的结果：

$$V^{\pi}(s) = \sum_{a\in \mathcal{A}} \pi(a|s)Q^\pi(s,a)$$

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

$$Q^\pi (s,a) = r(s,a) + \gamma \sum_{s' \in \mathcal{S}} P(s'|s,a) V^\pi (s')$$

3.4.4 贝尔曼期望方程

在贝尔曼方程中加上“期望”二字是为了与接下来的贝尔曼最优方程进行区分。我们通过简单推导就可以分别得到两个价值函数的贝尔曼期望方程（Bellman Expectation Equation）：

$$\begin{aligned} V^{\pi}(s) &= \mathbb{E}_{\pi}[G_t | S_t = s]\\ &= \sum_{a\in\mathcal{A}} \pi(a|s) Q^\pi(s,a)\\ &= \sum_{a\in\mathcal{A}} \pi(a|s) \Big(r(s,a) + \gamma \sum_{s'\in\mathcal{S}}p(s'|s,a)V^{\pi}(s') \Big)\\ Q^\pi (s,a) &= \mathbb{E}_{\pi} [G_t | S_t = s, A_t = a]\\ &= \mathbb{E}_{\pi} [R_t + \gamma Q^\pi(S_{t+1}, A_{t+1}) | S_t = s, A_t = a]\\ &= r(s,a) + \gamma \sum_{s' \in \mathcal{S}} p(s'|s,a) \sum_{a'\in\mathcal{A}} \pi(a'|s')Q^\pi(s',a') \end{aligned}$$

价值函数和贝尔曼方程是强化学习非常重要的组成部分，之后的一些强化学习算法都是据此推导出来的，读者需要明确掌握！

下图是一个马尔可夫决策过程的简单例子，其中每个深色圆圈代表一个状态，一共有从 $s_1 \sim s_5$ 这 $5$ 个状态。黑色实线箭头代表可以采取的动作，浅色小圆圈代表动作，需要注意，并非在每个状态都能采取所有动作，例如在状态 $s_1$，智能体只能采取“保持$s_1$”和“前往$s_2$”这两个动作，无法采取其他动作。

每个浅色小圆圈旁的数字代表在某个状态下采取某个动作能获得的奖励。虚线箭头代表采取动作后可能转移到的状态，箭头边上的数字代表转移概率，如果没有数字则表示转移概率为 $1$。例如，在$s_2$下， 如果采取动作“前往$s_3$”，就能得到奖励$-2$，并且转移到$s_3$；在$s_4$下，如果采取“概率前往”，就能得到奖励 $1$，并且会分别以概率 $0.2, 0.4, 0.4$ 转移到$s_2, s_3, s_4$。

接下来我们编写代码来表示上图中的马尔可夫决策过程，并定义两个策略。第一个策略是一个完全随机策略，即在每个状态下，智能体会以同样的概率选取它可能采取的动作。例如，在$s_1$下，智能体会以 $0.5$ 和 $0.5$ 的概率选取动作“保持$s_1$”和“前往$s_2$”。第二个策略是一个提前设定的一个策略。

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

接下来我们想要计算该 MDP 下，一个策略 $\pi$ 的状态价值函数。我们现在有的工具是 MRP 的解析解方法。于是，一个很自然的想法是：给定一个 MDP 和一个策略，我们是否可以将其转化为一个 MRP？答案是肯定的。我们可以将策略的动作选择进行边缘化（marginalization)，就可以得到没有动作的 MRP 了。具体来说，对于某一个状态，我们根据策略所有动作的概率进行加权，得到的奖励和就可以认为是一个 MRP 在该状态下的奖励，即：

$$r'(s) = \sum_{a \in \mathcal{A}} \pi(a|s) r(s, a)$$

同理，我们计算采取动作的概率与使 $s$ 转移到的 $s'$ 概率的乘积，再将这些乘积相加，其和就是一个 MRP 的状态从 $s$ 转移至 $s'$ 的概率：

$$P'(s'|s) = \sum_{a \in \mathcal{A}} \pi(a|s) P(s'|s, a)$$

于是，我们构建得到了一个 MRP: $(\mathcal{S, P', r', \gamma})$。根据价值函数的定义可以发现，转化前的 MDP 的状态价值函数和转化后的 MRP 的价值函数是一样的。于是我们可以用 MRP 中计算价值函数的解析解来计算这个 MDP 中该策略的状态价值函数。

P_from_mdp_to_mrp = np.zeros((5, 5))
R_from_mdp_to_mrp = np.zeros(5)

for s in range(5):
    pi_sum = 0
    for s_next in range(5):
        for action in A:
            P_from_mdp_to_mrp[s][s_next] += Pi_1.get(join(S[s], action), 0) * \
                                            P.get(join(join(S[s], action), S[s_next]), 0)
            pi_sum += Pi_1.get(join(S[s], action), 0)
    if pi_sum == 0:
        P_from_mdp_to_mrp[s][s] = 1.0

for s in range(5):
    for action in A:
        R_from_mdp_to_mrp[s] += Pi_1.get(join(S[s], action), 0) * \
                             R.get(join(S[s], action), 0)

print(P_from_mdp_to_mrp, '\n', R_from_mdp_to_mrp)

[[0.5 0.5 0.  0.  0. ]
 [0.5 0.  0.5 0.  0. ]
 [0.  0.  0.  0.5 0.5]
 [0.  0.1 0.2 0.2 0.5]
 [0.  0.  0.  0.  1. ]] 
 [-0.5 -1.5 -1.   5.5  0. ]

我们接下来就编写代码来实现该方法，计算用随机策略（也就是代码中的Pi_1）时的状态价值函数。为了简单起见，我们直接给出转化后的 MRP 的状态转移矩阵和奖励函数，感兴趣的读者可以自行验证。

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)

MDP中每个状态价值分别为
 [[-1.22555411]
 [-1.67666232]
 [ 0.51890482]
 [ 6.0756193 ]
 [ 0.        ]]

知道了状态价值函数 $V^{\pi}(s)$ 后，我们可以计算动作价值函数 $Q^{\pi}(s, a)$。例如（$s_4$，概率前往）的动作价值为 $2.152$，根据以下公式可以计算得到：

$$Q^\pi(s,a) = r(s,a) + \gamma \sum_{s'\in\mathcal{S}}P(s'|s,a)V_{\pi}(s') = 1 + 0.5 \times [0.2 \times (-1.68) + 0.4 \times 0.52 + 0.4 \times 6.08] = 2.152$$

这个 MRP 解析解的方法在状态动作集合比较大的时候不是很适用，那有没有其他的方法呢？第 4 章将介绍用动态规划算法来计算得到价值函数。3.5 节将介绍用蒙特卡洛方法来近似估计这个价值函数，用蒙特卡洛方法的好处在于我们不需要知道 MDP 的状态转移函数和奖励函数，它可以得到一个近似值，并且采样数越多越准确。

3.5 蒙特卡洛方法

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

我们现在介绍如何用蒙特卡洛方法来估计一个策略在一个马尔可夫决策过程中的状态价值函数。回忆一下，一个状态的价值是它的期望回报，那么一个很直观的想法就是用策略在 MDP 上采样很多条序列，计算从这个状态出发的回报再求其期望就可以了，公式如下：

$$V^{\pi}(s) = \mathbb{E}_{\pi} [G_t | S_t=s] \approx \dfrac{1}{N} \sum_{i=1}^N G_t^{(i)}$$

在一条序列中，可能没有出现过这个状态，可能只出现过一次这个状态，也可能出现过很多次这个状态。我们介绍的蒙特卡洛价值估计方法会在该状态每一次出现时计算它的回报。还有一种选择是一条序列只计算一次回报，也就是这条序列第一次出现该状态时计算后面的累积奖励，而后面再次出现该状态时，该状态就被忽略了。假设我们现在用策略 $\pi$ 从状态 $s$ 开始采样序列，据此来计算状态价值。我们为每一个状态维护一个计数器和总回报，计算状态价值的具体过程如下所示。

(1) 使用策略 $\pi$ 采样若干条序列：

$$s_0^{(i)} \overset{a_0^{(i)}}{\longrightarrow} r_0^{(i)}, s_1^{(i)} \overset{a_1^{(i)}}{\longrightarrow} r_1^{(i)}, s_2^{(i)} \overset{a_2^{(i)}}{\longrightarrow} \cdots \overset{a_{T-1}^{(i)}}{\longrightarrow} r_{T-1}^{(i)}, s_T^{(i)}$$

(2) 对每一条序列中的每一时间步的状态进行以下操作：

• 更新状态的计数器 $N(s) \leftarrow N(s) + 1$；

• 更新状态的总回报 $M(s) \leftarrow M(s) + G_t$；

(3) 每一个状态的价值被估计为回报的平均值 $V(s) = \dfrac{M(s)}{N(s)}$

根据大数定律，当 $N(s) \rightarrow \infin$，有 $V(s) \rightarrow V^{\pi}(s)$。计算回报的期望时，除了可以把所有的回报加起来除以次数，还有一种增量更新的方法。对于每个状态 $s$ 和对应回报 $G$ ，进行如下计算：

• $N(s) \leftarrow N(s) + 1$；

• $V(s) \leftarrow V(s) + \dfrac{1}{N(s)}(G - V(s))$；

这种增量式更新期望的方法已经在第 2 章中展示过。

接下来我们用代码定义一个采样函数。采样函数需要遵守状态转移矩阵和相应的策略，每次将 (s,a,r,s_next) 元组放入序列中，直到到达终止序列。然后我们通过该函数，用随机策略在MDP例子中随机采样几条序列。

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])

第一条序列
 [('s1', '前往s2', 0, 's2'), ('s2', '前往s3', -2, 's3'), ('s3', '前往s5', 0, 's5')]
第二条序列
 [('s4', '概率前往', 1, 's4'), ('s4', '前往s5', 10, 's5')]
第五条序列
 [('s2', '前往s3', -2, 's3'), ('s3', '前往s4', -2, 's4'), ('s4', '前往s5', 10, 's5')]

## 对所有采样序列计算所有状态的价值
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)

使用蒙特卡洛方法计算MDP的状态价值为
 {'s1': -1.228923788722258, 's2': -1.6955696284402704, 's3': 0.4823809701532294, 's4': 5.967514743019431, 's5': 0}

可以看到用蒙特卡洛方法估计得到的状态价值和我们用 MRP 解析解得到的状态价值是很接近的。这得益于我们采样了比较多的序列，感兴趣的读者可以尝试修改采样次数，然后观察蒙特卡洛方法的结果。

3.6 占用度量

3.4 节提到，不同策略的价值函数是不一样的。这是因为对于同一个 MDP，不同策略会访问到的状态的概率分布是不同的。想象一下，图 3-4 的 MDP 中现在有一个策略，它的动作执行会使得智能体尽快到达终止状态$s_5$，于是当智能体处于状态$s_3$时，不会采取“前往$s_4$”的动作，而只会以 $1$ 的概率采取“前往$s_5$”的动作，所以智能体也不会获得在状态$s_4$下采取“前往$s_5$”可以得到的很大的奖励 $10$。可想而知，根据贝尔曼方程，这个策略在状态$s_3$的概率会比较小，究其原因是因为它没法到达状态$s_4$。因此我们需要理解不同策略会使智能体访问到不同概率分布的状态这个事实，这会影响到策略的价值函数。

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

$$\nu^\pi(s) = (1-\gamma) \sum_{t=0}^\infin \gamma^t P_t^\pi(s)$$

其中，$1-\gamma$是用来使得概率加和为 $1$ 的归一化因子。状态访问概率表示一个策略和 MDP 交互会访问到的状态的分布。需要注意的是，理论上在计算该分布时需要交互到无穷步之后，但实际上智能体和 MDP 的交互在一个序列中是有限的。不过我们仍然可以用以上公式来表达状态访问概率的思想，状态访问概率有如下性质：

$$\nu^\pi(s') = (1-\gamma)\nu_0(s') + \gamma \iint P(s'|s,a) \pi(a|s) \nu^\pi(s) \mathbf{dsda}$$

此外，我们还可以定义策略的占用度量（occupancy measure）：

$$\rho^\pi (s,a) = (1 - \gamma) \sum_{t=0}^\infin \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} \Leftrightarrow \pi_1 = \pi_2$$

定理 2：给定一合法占用度量$\rho$，可生成该占用度量的唯一策略是

$$\pi_\rho = \dfrac{\rho(s,a)}{\sum\limits_{a' \in \mathcal{A}}\rho(s,a')}$$

注意：以上提到的“合法”占用度量是指存在一个策略使智能体与 MDP 交互产生的状态动作对被访问到的概率。

接下来我们编写代码来近似估计占用度量。这里我们采用近似估计，即设置一个较大的采样轨迹长度的最大值，然后采样很多次，用状态动作对出现的频率估计实际概率。

def occupancy(episodes, s, a, timestep_max, gamma):
    ''' 计算状态动作对（s,a）出现的频率,以此来估算策略的占用度量 '''
    rho = 0
    total_times = np.zeros(timestep_max)  # 记录每个时间步t各被经历过几次
    occur_times = np.zeros(timestep_max)  # 记录(s_t,a_t)=(s,a)的次数
    for episode in episodes:
        for i in range(len(episode)):
            (s_opt, a_opt, r, s_next) = episode[i]
            total_times[i] += 1
            if s == s_opt and a == a_opt:
                occur_times[i] += 1
    for i in reversed(range(timestep_max)):
        if total_times[i]:
            rho += gamma**i * occur_times[i] / total_times[i]
    return (1 - gamma) * rho


gamma = 0.5
timestep_max = 1000

episodes_1 = sample(MDP, Pi_1, timestep_max, 1000)
episodes_2 = sample(MDP, Pi_2, timestep_max, 1000)
rho_1 = occupancy(episodes_1, "s4", "概率前往", timestep_max, gamma)
rho_2 = occupancy(episodes_2, "s4", "概率前往", timestep_max, gamma)
print(rho_1, rho_2)

0.112567796310472 0.23199480615618912

通过以上结果可以发现，不同策略对于同一个状态动作对的占用度量是不一样的。

3.7 最优策略

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

最优策略都有相同的状态价值函数，我们称之为最优状态价值函数，表示为：

$$V^*(s) = \max_\pi V^\pi(s), \quad \forall s \in \mathcal{S}$$

同理，我们定义最优动作价值函数:

$$Q^*(s,a) = \max_\pi Q^\pi(s,a), \quad \forall s\in\mathcal{S}, a \in \mathcal{A}$$

为了使$Q^\pi(s,a)$最大，我们需要在当前的状态动作对 $(s,a)$ 之后都执行最优策略。于是我们得到了最优状态价值函数和最优动作价值函数之间的关系：

$$Q^*(s,a) = r(s,a) + \gamma \sum_{s' \in \mathcal{S}} P(s'|s,a)V^*(s')$$

这与在普通策略下的状态价值函数和动作价值函数之间的关系是一样的。另一方面，最优状态价值是选择此时使最优动作价值最大的那一个动作时的状态价值：

$$V^*(s) = \max_{a\in\mathcal{A}}Q^*(s,a)$$

3.7.1 贝尔曼最优方程

根据$V^*(s)$和$Q^*(s,a)$的关系，我们可以得到贝尔曼最优方程（Bellman optimality equation）：

$$\begin{aligned} V^*(s) &= \max_{a\in\mathcal{A}} \Big\{ r(s,a) + \gamma \sum_{s' \in \mathcal{S}} p(s'|s,a) V^*(s') \Big\}\\ Q^*(s,a) &= r(s,a) + \gamma \sum_{s'\in\mathcal{S}} p(s'|s,a) \max_{a'\in\mathcal{A}}Q^*(s',a') \end{aligned}$$

第 4 章将介绍如何用动态规划算法得到最优策略。

3.8 总结

本章从零开始介绍了马尔可夫决策过程的基础概念知识，并讲解了如何通过求解贝尔曼方程得到状态价值的解析解以及如何用蒙特卡洛方法估计各个状态的价值。马尔可夫决策过程是强化学习中的基础概念，强化学习中的环境就是一个马尔可夫决策过程。我们接下来将要介绍的强化学习算法通常都是在求解马尔可夫决策过程中的最优策略。

3.9 参考文献

[1] SUTTON R S, BARTO A G. Reinforcement learning: an introduction [M]. Cambridge:MIT press, 2018.

[2] OTTERLO M V, WIERING M. Reinforcement learning and markov decision processes [M]. Berlin, Heidelberg: Springer, 2012: 3-42.

4 动态规划算法

4.1 简介

动态规划（dynamic programming）是程序设计算法中非常重要的内容，能够高效解决一些经典问题，例如背包问题和最短路径规划。动态规划的基本思想是将待求解问题分解成若干个子问题，先求解子问题，然后从这些子问题的解得到目标问题的解。动态规划会保存已解决的子问题的答案，在求解目标问题的过程中，需要这些子问题答案时就可以直接利用，避免重复计算。本章介绍如何用动态规划的思想来求解在马尔可夫决策过程中的最优策略。

基于动态规划的强化学习算法主要有两种：一是策略迭代（policy iteration），二是价值迭代（value iteration）。其中，策略迭代由两部分组成：策略评估（policy evaluation）和策略提升（policy improvement）。具体来说，策略迭代中的策略评估使用贝尔曼期望方程来得到一个策略的状态价值函数，这是一个动态规划的过程；而价值迭代直接使用贝尔曼最优方程来进行动态规划，得到最终的最优状态价值。

不同于 3.5 节介绍的蒙特卡洛方法和第 5 章将要介绍的时序差分算法，基于动态规划的这两种强化学习算法要求事先知道环境的状态转移函数和奖励函数，也就是需要知道整个马尔可夫决策过程。在这样一个白盒环境中，不需要通过智能体和环境的大量交互来学习，可以直接用动态规划求解状态价值函数。但是，现实中的白盒环境很少，这也是动态规划算法的局限之处，我们无法将其运用到很多实际场景中。另外，策略迭代和价值迭代通常只适用于有限马尔可夫决策过程，即状态空间和动作空间是离散且有限的。

4.2 悬崖漫步环境

本节使用策略迭代和价值迭代来求解悬崖漫步（Cliff Walking）这个环境中的最优策略。接下来先简单介绍一下该环境。

悬崖漫步是一个非常经典的强化学习环境，它要求一个智能体从起点出发，避开悬崖行走，最终到达目标位置。如图 4-1 所示，有一个 4×12 的网格世界，每一个网格表示一个状态。智能体的起点是左下角的状态，目标是右下角的状态，智能体在每一个状态都可以采取 4 种动作：上、下、左、右。如果智能体采取动作后触碰到边界墙壁则状态不发生改变，否则就会相应到达下一个状态。环境中有一段悬崖，智能体掉入悬崖或到达目标状态都会结束动作并回到起点，也就是说掉入悬崖或者达到目标状态是终止状态。智能体每走一步的奖励是 −1，掉入悬崖的奖励是 −100。

图4-1 Cliff Walking环境示意图

接下来一起来看一看 Cliff Walking 环境的代码吧。

import copy


class CliffWalkingEnv:
    """ 悬崖漫步环境"""
    def __init__(self, ncol=12, nrow=4):
        self.ncol = ncol  ## 定义网格世界的列
        self.nrow = nrow  ## 定义网格世界的行
        ## 转移矩阵P[state][action] = [(p, next_state, reward, done)]包含下一个状态和奖励
        self.P = self.createP()

    def createP(self):
        ## 初始化
        P = [[[] for j in range(4)] for i in range(self.nrow * self.ncol)]
        ## 4种动作, change[0]:上,change[1]:下, change[2]:左, change[3]:右。坐标系原点(0,0)
        ## 定义在左上角
        change = [[0, -1], [0, 1], [-1, 0], [1, 0]]
        for i in range(self.nrow):
            for j in range(self.ncol):
                for a in range(4):
                    ## 位置在悬崖或者目标状态,因为无法继续交互,任何动作奖励都为0
                    if i == self.nrow - 1 and j > 0:
                        P[i * self.ncol + j][a] = [(1, i * self.ncol + j, 0,
                                                    True)]
                        continue
                    ## 其他位置
                    next_x = min(self.ncol - 1, max(0, j + change[a][0]))
                    next_y = min(self.nrow - 1, max(0, i + change[a][1]))
                    next_state = next_y * self.ncol + next_x
                    reward = -1
                    done = False
                    ## 下一个位置在悬崖或者终点
                    if next_y == self.nrow - 1 and next_x > 0:
                        done = True
                        if next_x != self.ncol - 1:  ## 下一个位置在悬崖
                            reward = -100
                    P[i * self.ncol + j][a] = [(1, next_state, reward, done)]
        return P    

4.3 策略迭代算法

策略迭代是策略评估和策略提升不断循环交替，直至最后得到最优策略的过程。本节分别对这两个过程进行详细介绍。

4.3.1 策略评估

策略评估这一过程用来计算一个策略的状态价值函数。回顾一下之前学习的贝尔曼期望方程：

$$V^\pi(s) = \sum_{a \in \mathcal{A}} \pi(a|s) \Big( r(s,a) + \gamma \sum_{s'\in\mathcal{S}} p(s'|s,a)V^\pi(s') \Big)$$

其中，$\pi(a|s)$是策略$\pi$在状态$s$下采取动作$a$的概率。可以看到，当知道奖励函数和状态转移函数时，我们可以根据下一个状态的价值来计算当前状态的价值。因此，根据动态规划的思想，可以把计算下一个可能状态的价值当成一个子问题，把计算当前状态的价值看作当前问题。在得知子问题的解后，就可以求解当前问题。更一般的，考虑所有的状态，就变成了用上一轮的状态价值函数来计算当前这一轮的状态价值函数，即

$$V^{k+1}(s) = \sum_{a \in \mathcal{A}} \pi(a|s) \Big( r(s,a) + \gamma \sum_{s'\in\mathcal{S}} p(s'|s,a)V^{k}(s') \Big)$$

我们可以选定任意初始值$V_0$。根据贝尔曼期望方程，可以得知$V^k=V^\pi$是以上更新公式的一个不动点（fixed point）。事实上，可以证明当$k \rightarrow \infin$时，序列$\{V^k\}$会收敛到$V^\pi$，所以可以据此来计算得到一个策略的状态价值函数。可以看到，由于需要不断做贝尔曼期望方程迭代，策略评估其实会耗费很大的计算代价。在实际的实现过程中，如果某一轮 $\underset{s\in\mathcal{S}}{\max}|V^{k+1}(s) - V^k(s)|$ 的值非常小，可以提前结束策略评估。这样做可以提升效率，并且得到的价值也非常接近真实的价值。

4.3.2 策略提升

使用策略评估计算得到当前策略的状态价值函数之后，我们可以据此来改进该策略。假设此时对于策略$\pi$，我们已经知道其价值$V^\pi$，也就是知道了在策略$\pi$下从每一个状态$s$出发最终得到的期望回报。我们要如何改变策略来获得在状态$s$下更高的期望回报呢？假设智能体在状态$s$下采取动作$a$，之后的动作依旧遵循策略$\pi$，此时得到的期望回报其实就是动作价值$Q^\pi(s,a)$。如果我们有$Q^\pi(s,a) > V^\pi(s)$，则说明在状态$s$下采取动作$a$会比原来的策略$\pi(a|s)$得到更高的期望回报。以上假设只是针对一个状态，现在假设存在一个确定性策略$\pi'$，在任意一个状态$s$下，都满足

$$Q^{\pi}(s,\pi'(s)) \ge V^\pi(s)$$

于是在任意状态$s$下，我们有

$$V^{\pi'}(s) \ge V^\pi(s)$$

这便是策略提升定理（policy improvement theorem）。于是我们可以直接贪心地在每一个状态选择动作价值最大的动作，也就是

$$\pi'(s) = \underset{a}{\arg\max} Q^\pi (s,a) = \underset{a}{\arg\max} \{r(s,a) + \gamma \sum_{s'}P(s'|s,a)V^\pi(s')\}$$

我们发现构造的贪心策略$\pi'$满足策略提升定理的条件，所以策略$\pi'$能够比策略$\pi$更好或者至少与其一样好。这个根据贪心法选取动作从而得到新的策略的过程称为策略提升。当策略提升之后得到的策略$\pi'$和之前的策略$\pi$一样时，说明策略迭代达到了收敛，此时$\pi$和$\pi'$就是最优策略。

策略提升定理的证明可以通过以下推导过程得证，使用上述提升公式得到的新策略$\pi'$在每个状态的价值不低于原策略$\pi$在该状态的价值。

$$\begin{aligned} V^\pi(s) &\le Q^\pi(s, \pi'(s))\\ &= \mathbb{E}_{\pi'} [R_t + \gamma V^\pi(S_{t+1})|S_t = s]\\ &\le \mathbb{E}_{\pi'} [R_t + \gamma Q^\pi(S_{t+1}, \pi'(S_{t+1}))|S_t = s]\\ &= \mathbb{E}_{\pi'} [R_t + \gamma R_{t+1} + \gamma^2V^\pi(S_{t+2})|S_t = s]\\ &\le \mathbb{E}_{\pi'} [R_t + \gamma R_{t+1} + \gamma^2R_{t+2} + \gamma^3V^\pi(S_{t+3})|S_t = s]\\ & \cdots\\ &\le \mathbb{E}_{\pi'} [R_t + \gamma R_{t+1} + \gamma^2R_{t+2} + \gamma^3R_{t+3} + \cdots|S_t = s]\\ &=V^{\pi'}(s) \end{aligned}$$

可以看到，推导过程中的每一个时间步都用到局部动作价值优势$V^\pi(S_{t+1}) \le Q^\pi(S_{t+1}, \pi'(S_{t+1}))$，累积到无穷步或者终止状态时，我们就得到了整个策略价值提升的不等式。

4.3.3 策略迭代算法

总体来说，策略迭代算法的过程如下：对当前的策略进行策略评估，得到其状态价值函数，然后根据该状态价值函数进行策略提升以得到一个更好的新策略，接着继续评估新策略、提升策略……直至最后收敛到最优策略（收敛性证明参见 4.7 节）：

$$\pi^0 \overset{\text{策略评估}}{\longrightarrow} V^{\pi^0} \overset{\text{策略提升}}{\longrightarrow} \pi^1 \overset{\text{策略评估}}{\longrightarrow} V^{\pi^1} \overset{\text{策略提升}}{\longrightarrow} \pi^2 \overset{\text{策略评估}}{\longrightarrow} \cdots \overset{\text{策略提升}}{\longrightarrow} \pi^*$$

结合策略评估和策略提升，我们得到以下策略迭代算法：

• 随机初始化策略$\pi(s)$和价值函数$V(s)$

• while $\Delta > \theta$ do: (策略评估循环)

  • $\Delta \leftarrow 0$

  • 对于每一个状态 $s\in\mathcal{S}$:

    • $v \leftarrow V(S)$

    • $V(s) \leftarrow r(s, \pi(s)) + \gamma \sum_{s'} P(s'|s, \pi(s))V(s')$

    • $\Delta \leftarrow \max(\Delta, |v - V(s)|)$

• end while

• $\pi_{\text{old}} \leftarrow \pi$

• 对于每一个状态 $s\in\mathcal{S}$：

  • $\pi(s) \leftarrow \underset{a\in\mathcal{A}}{\arg\max} \{r(s,a) + \gamma\underset{s'}{\sum} p(s'|s,a)V(s')\}$

• 若$\pi_{\text{old}} = \pi$，则停止算法并返回$V$和$\pi$; 否则转到策略评估循环

我们现在来看一下策略迭代算法的代码实现过程。

class PolicyIteration:
    """ 策略迭代算法 """
    def __init__(self, env, theta, gamma):
        self.env = env
        self.v = [0] * self.env.ncol * self.env.nrow  ## 初始化价值为0
        self.pi = [[0.25, 0.25, 0.25, 0.25]
                   for i in range(self.env.ncol * self.env.nrow)]  ## 初始化为均匀随机策略
        self.theta = theta  ## 策略评估收敛阈值
        self.gamma = gamma  ## 折扣因子

    def policy_evaluation(self):  ## 策略评估
        cnt = 1  ## 计数器
        while 1:
            max_diff = 0
            new_v = [0] * self.env.ncol * self.env.nrow
            for s in range(self.env.ncol * self.env.nrow):
                qsa_list = []  ## 开始计算状态s下的所有Q(s,a)价值
                for a in range(4):
                    qsa = 0
                    for res in self.env.P[s][a]:
                        p, next_state, r, done = res
                        qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))
                        ## 本章环境比较特殊,奖励和下一个状态有关,所以需要和状态转移概率相乘
                    qsa_list.append(self.pi[s][a] * qsa)
                new_v[s] = sum(qsa_list)  ## 状态价值函数和动作价值函数之间的关系
                max_diff = max(max_diff, abs(new_v[s] - self.v[s]))
            self.v = new_v
            if max_diff < self.theta: break  ## 满足收敛条件,退出评估迭代
            cnt += 1
        print("策略评估进行%d轮后完成" % cnt)

    def policy_improvement(self):  ## 策略提升
        for s in range(self.env.nrow * self.env.ncol):
            qsa_list = []
            for a in range(4):
                qsa = 0
                for res in self.env.P[s][a]:
                    p, next_state, r, done = res
                    qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))
                qsa_list.append(qsa)
            maxq = max(qsa_list)
            cntq = qsa_list.count(maxq)  ## 计算有几个动作得到了最大的Q值
            ## 让这些动作均分概率
            self.pi[s] = [1 / cntq if q == maxq else 0 for q in qsa_list]
        print("策略提升完成")
        return self.pi

    def policy_iteration(self):  ## 策略迭代
        while 1:
            self.policy_evaluation()
            old_pi = copy.deepcopy(self.pi)  ## 将列表进行深拷贝,方便接下来进行比较
            new_pi = self.policy_improvement()
            if old_pi == new_pi: break

现在我们已经写好了环境代码和策略迭代代码。为了更好地展现最终的策略，接下来增加一个打印策略的函数，用于打印当前策略在每个状态下的价值以及智能体会采取的动作。对于打印出来的动作，我们用^o<o表示等概率采取向左和向上两种动作，ooo>表示在当前状态只采取向右动作。

def print_agent(agent, action_meaning, disaster=[], end=[]):
    print("状态价值：")
    for i in range(agent.env.nrow):
        for j in range(agent.env.ncol):
            ## 为了输出美观,保持输出6个字符
            print('%6.6s' % ('%.3f' % agent.v[i * agent.env.ncol + j]), end=' ')
        print()

    print("策略：")
    for i in range(agent.env.nrow):
        for j in range(agent.env.ncol):
            ## 一些特殊的状态,例如悬崖漫步中的悬崖
            if (i * agent.env.ncol + j) in disaster:
                print('****', end=' ')
            elif (i * agent.env.ncol + j) in end:  ## 目标状态
                print('EEEE', end=' ')
            else:
                a = agent.pi[i * agent.env.ncol + j]
                pi_str = ''
                for k in range(len(action_meaning)):
                    pi_str += action_meaning[k] if a[k] > 0 else 'o'
                print(pi_str, end=' ')
        print()


env = CliffWalkingEnv()
action_meaning = ['^', 'v', '<', '>']
theta = 0.001
gamma = 0.9
agent = PolicyIteration(env, theta, gamma)
agent.policy_iteration()
print_agent(agent, action_meaning, list(range(37, 47)), [47])

策略评估进行60轮后完成
策略提升完成
策略评估进行72轮后完成
策略提升完成
策略评估进行44轮后完成
策略提升完成
策略评估进行12轮后完成
策略提升完成
策略评估进行1轮后完成
策略提升完成
状态价值：
-7.712 -7.458 -7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710
-7.458 -7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710 -1.900
-7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710 -1.900 -1.000
-7.458  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
策略：
ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovoo
ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovoo
ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ovoo
^ooo **** **** **** **** **** **** **** **** **** **** EEEE

经过 5 次策略评估和策略提升的循环迭代，策略收敛了，此时将获得的策略打印出来。用贝尔曼最优方程去检验其中每一个状态的价值，可以发现最终输出的策略的确是最优策略。

4.4 价值迭代算法

从上面的代码运行结果中我们能发现，策略迭代中的策略评估需要进行很多轮才能收敛得到某一策略的状态函数，这需要很大的计算量，尤其是在状态和动作空间比较大的情况下。我们是否必须要完全等到策略评估完成后再进行策略提升呢？试想一下，可能出现这样的情况：虽然状态价值函数还没有收敛，但是不论接下来怎么更新状态价值，策略提升得到的都是同一个策略。如果只在策略评估中进行一轮价值更新，然后直接根据更新后的价值进行策略提升，这样是否可以呢？答案是肯定的，这其实就是本节将要讲解的价值迭代算法，它可以被认为是一种策略评估只进行了一轮更新的策略迭代算法。需要注意的是，价值迭代中不存在显式的策略，我们只维护一个状态价值函数。

确切来说，价值迭代可以看成一种动态规划过程，它利用的是贝尔曼最优方程：

$$V^*(s) = \max_{a\in\mathcal{A}} \{ r(s,a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a)V^*(s') \}$$

将其写成迭代更新的方式为

$$V^{k+1}(s) = \max_{a\in\mathcal{A}} \{ r(s,a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a)V^k(s') \}$$

价值迭代便是按照以上更新方式进行的。等到$V^{k+1}$和$V^k$相同时，它就是贝尔曼最优方程的不动点，此时对应着最优状态价值函数$V^*$。然后我们利用$\pi(s) = \underset{a\in\mathcal{A}}{\arg\max} \{r(s,a) + \gamma\underset{s'}{\sum} p(s'|s,a)V^*(s')\}$，从中恢复出最优策略即可。

价值迭代算法流程如下：

• 随机初始化$V(s)$

• while $\Delta > \theta$ do:

  • $\Delta \leftarrow 0$

  • 对于每一个状态 $s\in\mathcal{S}$:

    • $v \leftarrow V(S)$

    • $V(s) \leftarrow \max_a \{r(s, a) + \gamma \sum_{s'} P(s'|s, a)V(s')\}$

    • $\Delta \leftarrow \max(\Delta, |v - V(s)|)$

• end while

• 返回一个确定性策略 $\pi(s) = \underset{a\in\mathcal{A}}{\arg\max} \{r(s,a) + \gamma\underset{s'}{\sum} p(s'|s,a)V(s')\}$

我们现在来编写价值迭代的代码。

class ValueIteration:
    """ 价值迭代算法 """
    def __init__(self, env, theta, gamma):
        self.env = env
        self.v = [0] * self.env.ncol * self.env.nrow  ## 初始化价值为0
        self.theta = theta  ## 价值收敛阈值
        self.gamma = gamma
        ## 价值迭代结束后得到的策略
        self.pi = [None for i in range(self.env.ncol * self.env.nrow)]

    def value_iteration(self):
        cnt = 0
        while 1:
            max_diff = 0
            new_v = [0] * self.env.ncol * self.env.nrow
            for s in range(self.env.ncol * self.env.nrow):
                qsa_list = []  ## 开始计算状态s下的所有Q(s,a)价值
                for a in range(4):
                    qsa = 0
                    for res in self.env.P[s][a]:
                        p, next_state, r, done = res
                        qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))
                    qsa_list.append(qsa)  ## 这一行和下一行代码是价值迭代和策略迭代的主要区别
                new_v[s] = max(qsa_list)
                max_diff = max(max_diff, abs(new_v[s] - self.v[s]))
            self.v = new_v
            if max_diff < self.theta: break  ## 满足收敛条件,退出评估迭代
            cnt += 1
        print("价值迭代一共进行%d轮" % cnt)
        self.get_policy()

    def get_policy(self):  ## 根据价值函数导出一个贪婪策略
        for s in range(self.env.nrow * self.env.ncol):
            qsa_list = []
            for a in range(4):
                qsa = 0
                for res in self.env.P[s][a]:
                    p, next_state, r, done = res
                    qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))
                qsa_list.append(qsa)
            maxq = max(qsa_list)
            cntq = qsa_list.count(maxq)  ## 计算有几个动作得到了最大的Q值
            ## 让这些动作均分概率
            self.pi[s] = [1 / cntq if q == maxq else 0 for q in qsa_list]


env = CliffWalkingEnv()
action_meaning = ['^', 'v', '<', '>']
theta = 0.001
gamma = 0.9
agent = ValueIteration(env, theta, gamma)
agent.value_iteration()
print_agent(agent, action_meaning, list(range(37, 47)), [47])

价值迭代一共进行14轮
状态价值：
-7.712 -7.458 -7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710
-7.458 -7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710 -1.900
-7.176 -6.862 -6.513 -6.126 -5.695 -5.217 -4.686 -4.095 -3.439 -2.710 -1.900 -1.000
-7.458  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
策略：
ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovoo
ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovo> ovoo
ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ooo> ovoo
^ooo **** **** **** **** **** **** **** **** **** **** EEEE

可以看到，解决同样的训练任务，价值迭代总共进行了数十轮，而策略迭代中的策略评估总共进行了数百轮，价值迭代中的循环次数远少于策略迭代。

4.5 冰湖环境

除了悬崖漫步环境，本章还准备了另一个环境——冰湖（Frozen Lake）。冰湖环境的状态空间和动作空间是有限的，我们在该环境中也尝试一下策略迭代算法和价值迭代算法，以便更好地理解这两个算法。

冰湖是 OpenAI Gym 库中的一个环境。OpenAI Gym 库中包含了很多有名的环境，例如 Atari 和 MuJoCo，并且支持我们定制自己的环境。在之后的章节中，我们还会使用到更多来自 OpenAI Gym 库的环境。如图 4-2 所示，冰湖环境和悬崖漫步环境相似，也是一个网格世界，大小为$4\times 4$。每一个方格是一个状态，智能体起点状态在左上角，目标状态在右下角，中间还有若干冰洞$H$。在每一个状态都可以采取上、下、左、右 4 个动作。由于智能体在冰面行走，因此每次行走都有一定的概率滑行到附近的其它状态，并且到达冰洞或目标状态时行走会提前结束。每一步行走的奖励是 0，到达目标的奖励是 1。

图4-2 Frozen Lake环境示意图

我们先创建 OpenAI Gym 中的 FrozenLake-v0 环境，并简单查看环境信息，然后找出冰洞和目标状态。

import gym
env = gym.make("FrozenLake-v0")  ## 创建环境
env = env.unwrapped  ## 解封装才能访问状态转移矩阵P
env.render()  ## 环境渲染,通常是弹窗显示或打印出可视化的环境

holes = set()
ends = set()
for s in env.P:
    for a in env.P[s]:
        for s_ in env.P[s][a]:
            if s_[2] == 1.0:  ## 获得奖励为1,代表是目标
                ends.add(s_[1])
            if s_[3] == True:
                holes.add(s_[1])
holes = holes - ends
print("冰洞的索引:", holes)
print("目标的索引:", ends)

for a in env.P[14]:  ## 查看目标左边一格的状态转移信息
    print(env.P[14][a])

SFFF
FHFH
FFFH
HFFG
冰洞的索引: {11, 12, 5, 7}
目标的索引: {15}
[(0.3333333333333333, 10, 0.0, False), (0.3333333333333333, 13, 0.0, False), (0.3333333333333333, 14, 0.0, False)]
[(0.3333333333333333, 13, 0.0, False), (0.3333333333333333, 14, 0.0, False), (0.3333333333333333, 15, 1.0, True)]
[(0.3333333333333333, 14, 0.0, False), (0.3333333333333333, 15, 1.0, True), (0.3333333333333333, 10, 0.0, False)]
[(0.3333333333333333, 15, 1.0, True), (0.3333333333333333, 10, 0.0, False), (0.3333333333333333, 13, 0.0, False)]

首先，我们发现冰洞的索引是$\{5,7,11,12\}$（集合 set 的索引是无序的），起点状态（索引为 0）在左上角，和悬崖漫步环境一样。其次，根据第 15 个状态（即目标左边一格，数组下标索引为 14）的信息，我们可以看到每个动作都会等概率“滑行”到 3 种可能的结果，这一点和悬崖漫步环境是不一样的。我们接下来先在冰湖环境中尝试一下策略迭代算法。

## 这个动作意义是Gym库针对冰湖环境事先规定好的
action_meaning = ['<', 'v', '>', '^']
theta = 1e-5
gamma = 0.9
agent = PolicyIteration(env, theta, gamma)
agent.policy_iteration()
print_agent(agent, action_meaning, [5, 7, 11, 12], [15])

策略评估进行25轮后完成
策略提升完成
策略评估进行58轮后完成
策略提升完成
状态价值：
 0.069  0.061  0.074  0.056
 0.092  0.000  0.112  0.000
 0.145  0.247  0.300  0.000
 0.000  0.380  0.639  0.000
策略：
<ooo ooo^ <ooo ooo^
<ooo **** <o>o ****
ooo^ ovoo <ooo ****
**** oo>o ovoo EEEE

这个最优策略很看上去比较反直觉，其原因是这是一个智能体会随机滑向其他状态的冰冻湖面。例如，在目标左边一格的状态，采取向右的动作时，它有可能会滑到目标左上角的位置，从该位置再次到达目标会更加困难，所以此时采取向下的动作是更为保险的，并且有一定概率能够滑到目标。我们再来尝试一下价值迭代算法。

action_meaning = ['<', 'v', '>', '^']
theta = 1e-5
gamma = 0.9
agent = ValueIteration(env, theta, gamma)
agent.value_iteration()
print_agent(agent, action_meaning, [5, 7, 11, 12], [15])

价值迭代一共进行60轮
状态价值：
 0.069  0.061  0.074  0.056
 0.092  0.000  0.112  0.000
 0.145  0.247  0.300  0.000
 0.000  0.380  0.639  0.000
策略：
<ooo ooo^ <ooo ooo^
<ooo **** <o>o ****
ooo^ ovoo <ooo ****
**** oo>o ovoo EEEE

可以发现价值迭代算法的结果和策略迭代算法的结果完全一致，这也互相验证了各自的结果。

4.6 小结

本章讲解了强化学习中两个经典的动态规划算法：策略迭代算法和价值迭代算法，它们都能用于求解最优价值和最优策略。动态规划的主要思想是利用贝尔曼方程对所有状态进行更新。需要注意的是，在利用贝尔曼方程进行状态更新时，我们会用到马尔可夫决策过程中的奖励函数和状态转移函数。如果智能体无法事先得知奖励函数和状态转移函数，就只能通过和环境进行交互来采样（状态-动作-奖励-下一状态）这样的数据，我们将在之后的章节中讲解如何求解这种情况下的最优策略。

4.7 扩展阅读：收敛性证明

4.7.1 策略迭代

策略迭代的过程如下：

$$\pi^0 \overset{\text{策略评估}}{\longrightarrow} V^{\pi^0} \overset{\text{策略提升}}{\longrightarrow} \pi^1 \overset{\text{策略评估}}{\longrightarrow} V^{\pi^1} \overset{\text{策略提升}}{\longrightarrow} \pi^2 \overset{\text{策略评估}}{\longrightarrow} \cdots \overset{\text{策略提升}}{\longrightarrow} \pi^*$$

根据策略提升定理，我们知道更新后的策略的价值函数满足单调性，即 $V^{\pi^{k+1}}\ge V^{\pi^{k}}$。所以只要所有可能的策略个数是有限的，策略迭代就能收敛到最优策略。假设 MDP 的状态空间大小为$|\mathcal S|$，动作空间大小为$|\mathcal A|$，此时所有可能的策略个数为$|\mathcal A|^{|\mathcal S|}$，是有限个，所以策略迭代能够在有限步找到其中的最优策略。

还有另一种类似的证明思路。在有限马尔可夫决策过程中，如果$\gamma<1$，那么很显然存在一个上界$C=\dfrac{R_{\text{max}}}{1-\gamma}$（这里的$R_{\text{max}}$为最大单步奖励值），使得对于任意策略$\pi$和状态$s$，其价值$V^\pi(s)<C$。因此，对于每个状态$s$，我们可以将策略迭代得到的价值写成数列$\{V^{\pi^k}\}_{k=1,\cdots,\infin}$。根据实数列的单调有界收敛定理，该数列一定收敛，也即是策略迭代算法一定收敛。

4.7.2 价值迭代

价值迭代的更新公式为：

$$V^{k+1}(s) = \max_{a\in\mathcal{A}} \{r(s,a) + \gamma \sum_{s' \in \mathcal{S}} P(s'|s,a) V^k(s')\}$$

我们将其定义为一个贝尔曼最优算子$\mathcal{T}$:

$$V^{k+1}(s) = \mathcal{T}V^k(s) = \max_{a\in\mathcal{A}} \{r(s,a) + \gamma \sum_{s' \in \mathcal{S}} P(s'|s,a) V^k(s')\}$$

然后我们引入压缩算子（contraction operator）：若$\mathcal O$是一个算子，如果满足$||\mathcal OV - \mathcal OV'||_q \le ||V - V'||_q$条件，则我们称$\mathcal O$是一个压缩算子。其中$||x||_q$表示x的$L_q$范数，包括我们将会用到的无穷范数$||x||_\infin = \max_i |x_i|$。

我们接下来证明当$\gamma < 1$时，贝尔曼最优算子$\mathcal{T}$是一个$\gamma$-压缩算子。

$$\begin{aligned} ||\mathcal{T}V - \mathcal{T}V'||_\infin &= \max_{s\in\mathcal{S}} \Bigg| \max_{a\in\mathcal{A}} \Big\{ r(s,a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a)V(s')\Big\} - \max_{a'\in\mathcal{A}} \Big\{ r(s,a') + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a')V'(s')\Big\} \Bigg|\\ &\le \max_{s,a} \Bigg| r(s,a) + \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a)V(s') - r(s,a) - \gamma \sum_{s'\in\mathcal{S}} P(s'|s,a)V'(s') \Bigg|\\ &= \gamma \max_{s,a}\Bigg|\sum_{s'\in\mathcal{S}} P(s'|s,a)\big(V(s') - V'(s')\big)\Bigg|\\ &\le \gamma \max_{s,a}\sum_{s'\in\mathcal{S}} P(s'|s,a) \max_{s'}\Big|V(s') - V'(s')\Big|\\ &=\gamma||V - V'||_{\infin} \end{aligned}$$