马尔可夫决策过程(Markov decision process,MDP)是强化学习的重要概念。要学好强化学习,我们首先要掌握马尔可夫决策过程的基础知识。前两章所说的强化学习中的环境一般就是一个马尔可夫决策过程。与多臂老虎机问题不同,马尔可夫决策过程包含状态信息以及状态之间的转移机制。如果要用强化学习去解决一个实际问题,第一步要做的事情就是把这个实际问题抽象为一个马尔可夫决策过程,也就是明确马尔可夫决策过程的各个组成要素。本章将从马尔可夫过程出发,一步一步地进行介绍,最后引出马尔可夫决策过程。
随机过程(stochastic process)是概率论的“动力学”部分。概率论的研究对象是静态的随机现象,而随机过程的研究对象是随时间演变的随机现象(例如天气随时间的变化、城市交通随时间的变化)。在随机过程中,随机现象在某时刻
当且仅当某时刻的状态只取决于上一时刻的状态时,一个随机过程被称为具有马尔可夫性质(Markov property),用公式表示为
马尔可夫过程(Markov process)指具有马尔可夫性质的随机过程,也被称为马尔可夫链(Markov chain)。我们通常用元组
矩阵
图 3-1 是一个具有 6 个状态的马尔可夫过程的简单例子。其中每个绿色圆圈表示一个状态,每个状态都有一定概率(包括概率为 0)转移到其他状态,其中
我们可以写出这个马尔可夫过程的状态转移矩阵:
其中第
给定一个马尔可夫过程,我们就可以从某个状态出发,根据它的状态转移矩阵生成一个状态序列(episode),这个步骤也被叫做采样(sampling)。例如,从
在马尔可夫过程的基础上加入奖励函数
在一个马尔可夫奖励过程中,从第
其中,
比如选取
接下来我们用代码表示图 3-2 中的马尔可夫奖励过程,并且定义计算回报的函数。
import numpy as npnp.random.seed(0)# 定义状态转移概率矩阵PP = [[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 = 0for i in reversed(range(start_index, len(chain))):G = gamma * G + rewards[chain[i] - 1]return G# 一个状态序列,s1-s2-s3-s6chain = [1, 2, 3, 6]start_index = 0G = compute_return(start_index, chain, gamma)print("根据本序列计算得到回报为:%s。" % G)
根据本序列计算得到回报为:-2.5。
在马尔可夫奖励过程中,一个状态的期望回报(即从这个状态出发的未来累积奖励的期望)被称为这个状态的价值(value)。所有状态的价值就组成了价值函数(value function),价值函数的输入为某个状态,输出为这个状态的价值。我们将价值函数写成
在上式的最后一个等号中,一方面,即时奖励的期望正是奖励函数的输出,即
上式就是马尔可夫奖励过程中非常有名的贝尔曼方程(Bellman equation),对每一个状态都成立。若一个马尔可夫奖励过程一共有
我们可以直接根据矩阵运算求解,得到以下解析解:
以上解析解的计算复杂度是
接下来编写代码来实现求解价值函数的解析解方法,并据此计算该马尔可夫奖励过程中所有状态的价值。
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 valueV = compute(P, rewards, gamma, 6)print("MRP中每个状态价值分别为\n", V)
MRP中每个状态价值分别为[[-2.01950168][-2.21451846][ 1.16142785][10.53809283][ 3.58728554][ 0. ]]
根据以上代码,求解得到各个状态的价值
我们现在用贝尔曼方程来进行简单的验证。例如,对于状态
可以发现左右两边的值几乎是相等的,说明我们求解得到的价值函数是满足状态为
3.2 节和 3.3 节讨论到的马尔可夫过程和马尔可夫奖励过程都是自发改变的随机过程;而如果有一个外界的“刺激”来共同改变这个随机过程,就有了马尔可夫决策过程(Markov decision process,MDP)。我们将这个来自外界的刺激称为智能体(agent)的动作,在马尔可夫奖励过程(MRP)的基础上加入动作,就得到了马尔可夫决策过程(MDP)。马尔可夫决策过程由元组
我们发现 MDP 与 MRP 非常相像,主要区别为 MDP 中的状态转移函数和奖励函数都比 MRP 多了动作
不同于马尔可夫奖励过程,在马尔可夫决策过程中,通常存在一个智能体来执行动作。例如,一艘小船在大海中随着水流自由飘荡的过程就是一个马尔可夫奖励过程,它如果凭借运气漂到了一个目的地,就能获得比较大的奖励;如果有个水手在控制着这条船往哪个方向前进,就可以主动选择前往目的地获得比较大的奖励。马尔可夫决策过程是一个与时间相关的不断进行的过程,在智能体和环境 MDP 之间存在一个不断交互的过程。一般而言,它们之间的交互是如图 3-3 循环过程:智能体根据当前状态
智能体的策略(Policy)通常用字母
我们用
不同于 MRP,在 MDP 中,由于动作的存在,我们额外定义一个动作价值函数(action-value function)。我们用
状态价值函数和动作价值函数之间的关系:在使用策略
使用策略
在贝尔曼方程中加上“期望”二字是为了与接下来的贝尔曼最优方程进行区分。我们通过简单推导就可以分别得到两个价值函数的贝尔曼期望方程(Bellman Expectation Equation):
价值函数和贝尔曼方程是强化学习非常重要的组成部分,之后的一些强化学习算法都是据此推导出来的,读者需要明确掌握!
图 3-4 是一个马尔可夫决策过程的简单例子,其中每个深色圆圈代表一个状态,一共有从
每个浅色小圆圈旁的数字代表在某个状态下采取某个动作能获得的奖励。虚线箭头代表采取动作后可能转移到的状态,箭头边上的数字代表转移概率,如果没有数字则表示转移概率为 1。例如,在
接下来我们编写代码来表示图 3-4 中的马尔可夫决策过程,并定义两个策略。第一个策略是一个完全随机策略,即在每个状态下,智能体会以同样的概率选取它可能采取的动作。例如,在
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,}# 策略2Pi_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 下,一个策略
同理,我们计算采取动作的概率与使
于是,我们构建得到了一个 MRP:
我们接下来就编写代码来实现该方法,计算用随机策略(也就是代码中的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. ]]
知道了状态价值函数
这个 MRP 解析解的方法在状态动作集合比较大的时候不是很适用,那有没有其他的方法呢?第 4 章将介绍用动态规划算法来计算得到价值函数。3.5 节将介绍用蒙特卡洛方法来近似估计这个价值函数,用蒙特卡洛方法的好处在于我们不需要知道 MDP 的状态转移函数和奖励函数,它可以得到一个近似值,并且采样数越多越准确。
蒙特卡洛方法(Monte-Carlo methods)也被称为统计模拟方法,是一种基于概率统计的数值计算方法。运用蒙特卡洛方法时,我们通常使用重复随机抽样,然后运用概率统计方法来从抽样结果中归纳出我们想求的目标的数值估计。一个简单的例子是用蒙特卡洛方法来计算圆的面积。例如,在图 3-5 所示的正方形内部随机产生若干个点,细数落在圆中点的个数,圆的面积与正方形面积之比就等于圆中点的个数与正方形中点的个数之比。如果我们随机产生的点的个数越多,计算得到圆的面积就越接近于真实的圆的面积。
我们现在介绍如何用蒙特卡洛方法来估计一个策略在一个马尔可夫决策过程中的状态价值函数。回忆一下,一个状态的价值是它的期望回报,那么一个很直观的想法就是用策略在 MDP 上采样很多条序列,计算从这个状态出发的回报再求其期望就可以了,公式如下:
在一条序列中,可能没有出现过这个状态,可能只出现过一次这个状态,也可能出现过很多次这个状态。我们介绍的蒙特卡洛价值估计方法会在该状态每一次出现时计算它的回报。还有一种选择是一条序列只计算一次回报,也就是这条序列第一次出现该状态时计算后面的累积奖励,而后面再次出现该状态时,该状态就被忽略了。假设我们现在用策略
(1) 使用策略
(2) 对每一条序列中的每一时间步
(3) 每一个状态的价值被估计为回报的平均值
根据大数定律,当
这种增量式更新期望的方法已经在第 2 章中展示过。
接下来我们用代码定义一个采样函数。采样函数需要遵守状态转移矩阵和相应的策略,每次将(s,a,r,s_next)
元组放入序列中,直到到达终止序列。然后我们通过该函数,用随机策略在图 3-4 的 MDP 中随机采样几条序列。
def sample(MDP, Pi, timestep_max, number):''' 采样函数,策略Pi,限制最长时间步timestep_max,总共采样序列数number '''S, A, P, R, gamma = MDPepisodes = []for _ in range(number):episode = []timestep = 0s = S[np.random.randint(4)] # 随机选择一个除s5以外的状态s作为起点# 当前状态为终止状态或者时间步太长时,一次采样结束while s != "s5" and timestep <= timestep_max:timestep += 1rand, temp = np.random.rand(), 0# 在状态s下根据策略选择动作for a_opt in A:temp += Pi.get(join(s, a_opt), 0)if temp > rand:a = a_optr = R.get(join(s, a), 0)breakrand, temp = np.random.rand(), 0# 根据状态转移概率得到下一个状态s_nextfor s_opt in S:temp += P.get(join(join(s, a), s_opt), 0)if temp > rand:s_next = s_optbreakepisode.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 = 0for i in range(len(episode) - 1, -1, -1): #一个序列从后往前计算(s, a, r, s_next) = episode[i]G = r + gamma * GN[s] = N[s] + 1V[s] = V[s] + (G - V[s]) / N[s]timestep_max = 20# 采样1000次,可以自行修改episodes = sample(MDP, Pi_1, timestep_max, 1000)gamma = 0.5V = {"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.4 节提到,不同策略的价值函数是不一样的。这是因为对于同一个 MDP,不同策略会访问到的状态的概率分布是不同的。想象一下,图 3-4 的 MDP 中现在有一个策略,它的动作执行会使得智能体尽快到达终止状态
首先我们定义 MDP 的初始状态分布为
其中,
此外,我们还可以定义策略的占用度量(occupancy measure):
它表示动作状态对
进一步得出如下两个定理。定理 1:智能体分别以策略
定理 2:给定一合法占用度量
注意:以上提到的“合法”占用度量是指存在一个策略使智能体与 MDP 交互产生的状态动作对被访问到的概率。
接下来我们编写代码来近似估计占用度量。这里我们采用近似估计,即设置一个较大的采样轨迹长度的最大值,然后采样很多次,用状态动作对出现的频率估计实际概率。
def occupancy(episodes, s, a, timestep_max, gamma):''' 计算状态动作对(s,a)出现的频率,以此来估算策略的占用度量 '''rho = 0total_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] += 1if s == s_opt and a == a_opt:occur_times[i] += 1for i in reversed(range(timestep_max)):if total_times[i]:rho += gamma**i * occur_times[i] / total_times[i]return (1 - gamma) * rhogamma = 0.5timestep_max = 1000episodes_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
通过以上结果可以发现,不同策略对于同一个状态动作对的占用度量是不一样的。
强化学习的目标通常是找到一个策略,使得智能体从初始状态出发能获得最多的期望回报。我们首先定义策略之间的偏序关系:当且仅当对于任意的状态
最优策略都有相同的状态价值函数,我们称之为最优状态价值函数,表示为:
同理,我们定义最优动作价值函数:
为了使
这与在普通策略下的状态价值函数和动作价值函数之间的关系是一样的。另一方面,最优状态价值是选择此时使最优动作价值最大的那一个动作时的状态价值:
根据
第 4 章将介绍如何用动态规划算法得到最优策略。
本章从零开始介绍了马尔可夫决策过程的基础概念知识,并讲解了如何通过求解贝尔曼方程得到状态价值的解析解以及如何用蒙特卡洛方法估计各个状态的价值。马尔可夫决策过程是强化学习中的基础概念,强化学习中的环境就是一个马尔可夫决策过程。我们接下来将要介绍的强化学习算法通常都是在求解马尔可夫决策过程中的最优策略。
[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.