动态规划9：马尔可夫决策过程

马尔可夫决策过程(MDP)

MDP可以看作是一个随机版本的动态规划，本期我们暂停一下“确定”版本的动态规划，来讨论一下MDP的基本设定，以及动态规划方法。

MDP的基本设定

马尔可夫决策过程(MDP)可以视为是马尔可夫过程的一个推广，首先我们要给定一个状态空间$\mathcal{S}$，本文总假定这个状态空间是有限的，用整数表示 S = { 0 , 1 , ⋯   , n − 1 } \mathcal{S}=\{0,1,\cdots,n-1\} S={0,1,⋯,n−1}。我们的决策主体，后面称为代理人(agent)，在$T$个阶段都作出决策，决策会产生状态转移，但与前面的股票买卖问题不同，这里的状态转移是随机的，取决于当前状态(state)和采取的行动(action)，每次状态转移都会产生回报，代理人的目标是指定一个策略(policy)，按照这个策略进行行动产生最大的累积回报。

基于这个想法，我们给出MDP的一些基本设定：

状态空间： 代理人可能处于的所有状态的集合，记为$\mathcal{S}$，我们这里假定状态空间是个有限集合，状态用自然数表示， S = { 0 , 1 , ⋯   , n − 1 } \mathcal{S}=\{0,1,\cdots,n-1\} S={0,1,⋯,n−1}

行为空间： 每个状态上我们都有能采取的行动集合，称为行为空间，行为空间是状态空间的函数，不同状态有不同的行为空间，状态$s\in \mathcal{S}$的行为空间我们记为$\mathcal{A}(s)$，行为空间可以是离散的，也可以是连续的，我们后面假设所有的行为空间都是有限的，用整数表示。

转移概率： 代理人处于某个状态$s_t\in \mathcal{S}$中采取行动$a\in \mathcal{A}(s)$，会产生状态转移，但是状态转移是随机的，转移到$s_{t+1}\in \mathcal{S}$有一定的概率，这一系列条件概率分布$
p(s_{t+1}|s_t,a_t)$称为转移概率，既然这是概率分布，就要满足
$$\begin{array}{rlr}& 0\le p(s|s_t,a_t)\le 1& s,s_t\in \mathcal{S},a_t\in \mathcal{A}(s_t)\\ & & \\ & \sum _{s\in \mathcal{S}}p(s|s_t,a_t)=1& s_t\in \mathcal{S},a_t\in \mathcal{A}(s_t)\end{array}$$

单步回报： 代理人某个状态$s_t\in \mathcal{S}$中采取行动$a\in \mathcal{A}(s)$，并转移到状态$s_{t+1}$时，会产生回报$R(s_t,a_t,s_{t+1})$，这个回报或确定，或随机，如果随机我们取其期望，这样我们总假定这个回报是确定的，与当前状态，采取行动以及下一步状态有关。

策略： 代理人要求解一个最优策略，使得累积回报最大化，在MDP中，策略是一系列状态的行为值函数或一系列状态空间上的概率分布；如果是前者，在有$T<\infty$的MDP中，策略表述为$$\pi =({\pi }_0,{\pi }_1,\cdots ,{\pi }_{T-1})$$
其中，${\pi }_i$是状态空间到行为空间的一个映射。那么按照这个策略，代理人会如何采取行动呢？

1. 在0时刻，代理人处于状态$s_0$，按照策略采取行动$a_0={\pi }_0(s_0)$，转移到状态$s_1$，产生回报$R(s_0,a_0,s_1)$
2. 在1时刻，代理人处于状态$s_1$，按照策略采取行动$a_1={\pi }_1(s_1)$，转移到状态$s_2$，产生回报$R(s_1,a_1,s_2)$
3. 以此类推，产生序列
   { s 0 , a 1 , s 1 , a 2 , ⋯   , s T − 1 , a T − 1 , s T } \{s_0,a_1,s_1,a_2,\cdots,s_{T-1},a_{T-1},s_T\} {s0​,a1​,s1​,a2​,⋯,sT−1​,aT−1​,sT​}
   用图表示就是：
   
   如果是随机策略，则$\pi$表示一个取决于$s_i$的行为空间上的概率分布，按照这个概率分布$\pi (s_i)$随机选择一个行为$a_i$，产生状态转移和回报，也会产生上图的一条轨迹。

轨迹： 代理人有一个确定的策略，给定一个初始状态，代理人按照策略行动，就会产生一条上图表示的状态和行为交替出现的轨迹，当然由于状态转移的随机性，以上轨迹也是随机的，一条轨迹通常用$\tau$表示

累计回报： 轨迹中每一步状态转移都会产生单步回报$R(s_t,a_t,s_{t+1})$，在MDP问题中，我们通常更看重当期的回报，所以未来的回报要通过一个折现率$\gamma (0<\gamma <1)$折现到当期，得到累积回报
$$\begin{array}{rl}R(\tau )=& \sum _{t=0}^{T-1}{\gamma }^{t}R(s_t,a_t,s_{t+1})\\ =& R(s_0,a_0,s_1)+\gamma R(s_1,a_1,s_2)+\cdots +{\gamma }^{T-1}R(s_{T-1},a_{T-1},s_T)\end{array}$$由于轨迹是随机的，所以累计回报也是随机的，轨迹取决于我们的策略$\pi$，在给定策略$\pi$时，不同轨迹出现的概率分布就能确定下来，就有累积期望回报$$E^{\pi }(R(\tau ))$$我们的目标是选一个策略$\pi$，使得累积期望回报最大化$$\underset{\pi }{\max }\left[E^{\pi }(R(\tau ))\right]$$

举例说明

我们给出一个例子来阐述上一节给出的MDP的设定（后面的例子纯属本人杜撰）。假设有一个商人，在三座小岛之间旅行贸易来赚钱，我们给这三座小岛编号为 { 0 , 1 , 2 } \{0,1,2\} {0,1,2}。这三个岛每个时刻都有两艘小船，但是商人不知道这两艘小船的目的地，只知道目的地的概率分布（有可能是观光船，环岛一周后回到原地）。每个小岛都有一个特产，商人到达一个新的小岛上，就要采购一个这个小岛的特产，然后到新的小岛上卖出去。

设0,1,2号小岛产出的特产分别是A,B,C，买卖都是按照市场价进行，市场价见下表

小岛	0	1	2
A	1	3	4
B	5	2	6
C	8	6	3

我们的状态空间就定义为小岛的集合，$s_t=1$就表示$t$时刻商人在小岛$1$上进行采购。行为空间就定义为商人坐哪艘船离开小岛。如果当前时刻商人在小岛$0$，采购当地的特产$A$，花费1块金币，再坐上某一艘小船离开小岛$0$；

1. 如果小船目的地是$0$，商人到达目的地后卖出$A$，得到1块金币，这样的利润是0
2. 如果小船目的地是$1$，商人到达目的地后卖出$A$，得到3块金币，这样的利润是2
3. 如果小船目的地是$2$，商人到达目的地后卖出$A$，得到4块金币，这样的利润是3

按照这种计算方法，可以得到回报函数矩阵

	0	1	2
0	0	2	3
1	3	0	4
2	5	3	0

这样，$R(1,a,2)=R(1,2)=4$，我们就计算了所有的回报函数值，如果我们排成一个三维矩阵，第一维是当前状态，第二维是行动，第三维是下一时刻的状态，回报函数矩阵如下

st=0	0	1	2
a1	0	2	3
a2	0	2	3

st=1	0	1	2
a3	3	0	4
a4	3	0	4

st=2	0	1	2
a5	5	3	0
a6	5	3	0

再给出转移概率矩阵，我们还是用一个三维矩阵来表示，第一维是当前状态，第二维是行动，第三维是下一时刻的状态，转移概率矩阵如下

st=0	0	1	2
a1	0.2	0.3	0.5
a2	0.3	0.3	0.4

st=1	0	1	2
a3	0.1	0.2	0.7
a4	0.2	0.1	0.7

st=2	0	1	2
a5	0.2	0.4	0.4
a6	0.5	0.3	0.2

以上设定可以用一个带权的有向图表示

边上画横线的数字是转移概率，不画横线的数字是回报。

状态价值函数与状态-行为价值函数

状态价值函数用来度量某个时刻某个状态的价值，那么很显然，价值函数我们就应该定义成该时刻的累积回报，即$$V_i(s)=E\left[\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})\Vert s_i=s\right]$$
但我们上面已经讨论过了，只有给定了某个具体的策略之后不同轨迹出现的概率才能计算出来，所以$V$还应该是$\pi$的函数，我们记为$$V_i^{\pi }(s)=E^{\pi }\left[\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})\Vert s_i=s\right]$$
那么状态价值函数该如何计算呢？我们假定策略是确定性策略

1. 如果$i=T-1$，那么
   $$\begin{array}{rl}{V}_{T-1}^{\pi }(s)=& E\left[R(s_{T-1},a_{T-1},s_T)|s_{T-1}=s\right]\\ =& E[R(s_{T-1},{\pi }_{T-1}(s_{T-1}),s_T)|s_{T-1}=s]\\ =& \sum _{s_T\in \mathcal{S}}p(s_T|s,{\pi }_{T-1}(s))R(s,{\pi }_{T-1}(s),s_T)\end{array}$$

2. 如果$i<T-1$，那么
   $$\begin{array}{rl}{V}_i^{\pi }(s)=& E^{\pi }\left[\sum _{t=i,s_i=s}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})|s_i=s\right]\\ =& E^{\pi }\left[R(s_i,a_i,s_{i+1})+\gamma \sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_i=s\right]\\ =& E^{\pi }\left[R(s_i,a_i,s_{i+1})|s_i=s\right]+\gamma E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_i=s\right]\end{array}$$
   上式的第一项很容易求解$$E^{\pi }\left[R(s_i,a_i,s_{i+1})\Vert s_i=s\right]=\sum _{s_{i+1}\in \mathcal{S}}p(s_{i+1}\Vert s,{\pi }_i(s))R(s,{\pi }_i(s),s_{i+1})$$
   我们再求解第二项：
   $$\begin{array}{rl}& E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_i=s\right]\\ =& \sum _{s_{i+1},a_{i+1},\cdots ,s_T}\left[P^{\pi }(s_{i+1},a_{i+1},\cdots ,s_T|s_i=s)\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})\right]\end{array}$$
   我们再来看上面的的条件概率
   $$\begin{array}{rl}& P^{\pi }(s_{i+1}=s_{i+1}^{\prime },a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_i=s)\\ =& \frac{P^{\pi }(s_i=s,s_{i+1}=s_{i+1}^{\prime },a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime })}{P^{\pi }(s_i=s)}\\ =& \frac{P^{\pi }(s_i=s,s_{i+1}=s_{i+1}^{\prime },a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime })}{P^{\pi }(s_i=s,s_{i+1}=s_{i+1}^{\prime })}\frac{P^{\pi }(s_i=s,s_{i+1}=s_{i+1}^{\prime })}{P^{\pi }(s_i=s)}\\ =& P^{\pi }(a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_i=s,s_{i+1}=s_{i+1}^{\prime })P^{\pi }(s_{i+1}=s_{i+1}^{\prime }|s_i=s)\end{array}$$
   对于上面的分解式的第二项，在确定性策略的条件下，给定$s_i$，那么$a_i$也由策略${\pi }_i$给定
   $$\begin{array}{rl}& P^{\pi }(s_{i+1}=s_{i+1}^{\prime },s_i=s)=P^{\pi }(s_i=s,a_i={\pi }_i(s),s_{i+1}=s_{i+1}^{\prime })\\ & P^{\pi }(s_i=s)=P^{\pi }(s_i=s,a_i={\pi }_i(s))\end{array}$$
   于是
   $$P^{\pi }(s_{i+1}=s_{i+1}^{\prime }|s_i=s)=P^{\pi }(s_{i+1}=s_{i+1}^{\prime }|s_i=s,a_i={\pi }_i(s))$$
   我们假定转移概率是时齐的，这个概率就是我们设定中的转移概率。再看第一项，我们假定整个过程具有马尔可夫性，也就是下一时刻的状态只与当前状态和行为有关，与之前的状态和行为都无关；$a_{i+1}$只与$s_{i+1}$有关，由策略${\pi }_{i+1}$给出，而$s_{i+2},\cdots ,s_{T-1}$都与$s_i$无关（马尔可夫性），于是
   $$\begin{array}{rl}& P^{\pi }(a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_i=s,s_{i+1}=s_{i+1}^{\prime })\\ =& P^{\pi }(a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_{i+1}=s_{i+1}^{\prime })\end{array}$$
   因此
   $$\begin{array}{rl}& P^{\pi }(s_{i+1}=s_{i+1}^{\prime },a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_i=s)\\ =& p(s_{i+1}^{\prime }|s,{\pi }_i(s))P^{\pi }(a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_{i+1}=s_{i+1}^{\prime })\end{array}$$
   于是
   $$\begin{array}{rl}& E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_i=s\right]\\ =& \sum _{s_{i+1},a_{i+1},\cdots ,s_T}\left[P^{\pi }(s_{i+1},a_{i+1},\cdots ,s_T|s_i=s)\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})\right]\\ =& \sum _{s_{i+1}=s^{\prime }\in \mathcal{S}}p(s^{\prime }|s,{\pi }_i(s))\sum _{a_{i+1},\cdots ,s_T}{P}^{\pi }(a_{i+1}=a_{i+1}^{\prime },\cdots ,s_T=s_T^{\prime }|s_{i+1}=s_{i+1}^{\prime })\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})\\ =& \sum _{s_{i+1}=s^{\prime }\in \mathcal{S}}p(s^{\prime }|s,{\pi }_i(s))E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_{i+1}=s_{i+1}^{\prime }\right]\\ =& \sum _{s_{i+1}=s^{\prime }\in \mathcal{S}}p(s^{\prime }|s,{\pi }_i(s))V_{i+1}^{\pi }(s^{\prime })\\ =& E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=s\right]\end{array}$$
   所以我们得到以下等式
   $$V_i^{\pi }(s)=E^{\pi }[R(s,{\pi }_i(s),s_{i+1})|s_i=s]+\gamma E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=s\right]$$
   这个等式我们称为Bellman等式，在随机策略的情况下也成立，是MDP的一个最重要的公式之一。这意味着求解MDP可以通过递归的方式进行，不难看出Bellman方程和动态规划的状态转移方程非常像，因此可以自底向上地求解出所有的状态价值的函数的值，本文在后面会给出求解所有状态价值函数值的C++代码。

我们可以从另外一个角度看Bellman方程，给定策略$\pi =({\pi }_0,\cdots ,\pi )$，我们以第二节给的海岛商人的例子来说明，假设海岛商人要在三道小岛之间转移五轮。

比如我们现在处于1时刻，在小岛1上，按照策略$\pi$我们要乘坐小船$a_3$，那么有$0.1$的概率开到小岛0上，转移的过程获得当期收益3，然后获得小岛0的价值$V_2^{\pi }(0)$，这个价值是未来的，所以要折现一期折现到现在，因此这个路线的价值为$$3+\gamma V_2^{\pi }(0)$$
有$0.2$的概率还是回到小岛1，当期收益为0，小岛1的价值为$V_2^{\pi }(1)$，这个收益是未来的，因此要折现到现在，价值为$$0+\gamma V_2^{\pi }(1)$$
有$0.7$的概率开到小岛2上，当期收益为4，小岛2的价值为$V_2^{\pi }(2)$，这个收益是未来的，要折算到现在$$4+\gamma V_2^{\pi }(2)$$
把这三种情况按照概率进行加权，就得到当期状态的价值函数
$$\begin{array}{rl}{V}_1^{\pi }(1)=& 0.1(3+\gamma V_2^{\pi }(0))+0.2(0+\gamma V_2^{\pi }(1))+0.7(4+V_2^{\pi }(2))\\ & \\ =& (0.1\times 3+0.2\times 0+0.7\times 4)+\gamma (0.1V_2^{\pi }(0)+0.2V_2^{\pi }(1)+0.7V_2^{\pi }(2))\end{array}$$

上面分解式第一项就是Bellman方程的第一项，第二项就是Bellman方程的第二项。第一项是按照策略采取行动获得当期收益的期望值，第二项是按照策略采取行动获得的未来价值的期望值再按照折现率$\gamma$进行折现，就得到当期的价值函数，这是Bellman方程比较直观的理解方式。这个例子也给出了给定确定性策略情况下计算所有状态价值函数的方法，我们来给出相应的C++代码。

首先，状态转移概率矩阵$P$和回报函数矩阵$R$用一个三维矩阵，第一个维度表示$s_t$，第二个维度表示$a_t$，第三个维度表示$s_{t+1}$，如$P[i][j][k]$表示在当前状态为$i$，采取行动$j$的条件下，转移到状态$k$的概率；价值函数我们用一个二维矩阵表示，$value\_function[i][j]$表示$V_i^{\pi }(j)$；策略我们也用一个二维矩阵表示，$policy[i][j]$表示$i$时刻处于状态$j$要采取的行为；我们定义一个MDP类，来储存这些参数，并给出构造函数

struct MDP_exception{
    int exception_id;
    string exception_info;
    MDP_exception(int id,const string& info):exception_id(id),exception_info(info){}
};
class MDP{
    int states,T;//state:状态的数量,T:阶段的数量
    double gamma;//折现因子
    vector<vector<vector<double>>> P;//转移概率矩阵
    vector<vector<vector<double>>> R;//回报矩阵
    vector<vector<int>> policy;//策略矩阵
    vector<vector<double>> value_function;//价值函数
    public:
        /*构造函数*/
        MDP(int nums_states=0,
        double discount_factor=0.5,
        int T=5,
        const vector<vector<vector<double>>>& transition_probability=vector<vector<vector<double>>>(),
        const vector<vector<vector<double>>>& rewards=vector<vector<vector<double>>>()){
            /*
            设定：
            nums_states：状态的数量，状态空间用[0,1,...,nums_states-1]表示
            T：阶段数，必须大于或等于1，默认为5
            initial_states：代理人所处的初始状态，如果不是空的MDP，即状态数为0，那么就必须在0和nums_states-1之间
            discount_factor：折现因子，储存在MDP内是discount_factor，折现因子大于0，并且小于1，如果违反会抛出异常
            transition_probability：转移概率矩阵，在构造函数中会对转移概率矩阵进行合法性检查，否则会抛出异常，转移概率矩阵合法的条件如下：
                1. 转移概率矩阵的第一个维度要等于状态数，即transition_probability.size()==nums_states
                2. 转移概率矩阵每一个矩阵的每一行之和都要等于1，且每一个元都要在0到1之间
                3. 转移概率矩阵的每一列的维度要等于状态数，即transition_probability[i][j].size()==nums_states
                除此以外，每个转移概率矩阵的行数没有限制，表示每个状态可以采取的行动的数量
            rewards：回报矩阵的维度必须和转移概率矩阵的维度一模一样
            */
            if(discount_factor<=0||discount_factor>=1) throw MDP_exception(1,"Discount factor <=0 || >=1");//检查折现因子的合法性
            else if(T<=0) throw MDP_exception(2,"nums_of_stages <=0");//检查T的合法性
            else if(nums_states<0) throw MDP_exception(3,"illegal states number");
            else if(transition_probability.size()!=nums_states) throw MDP_exception(5,"transition_probability.size() != nums_states");
            else if(rewards.size()!=nums_states) throw MDP_exception(6,"rewards.size()!=nums_states");
            else{
                //对转移概率矩阵和回报矩阵进行逐行检查
                for(int i=0;i<nums_states;i++){
                    //检查行数是否相等
                    if(transition_probability[i].size()!=rewards[i].size()) throw MDP_exception(7,"shape of transition probability != shape of rewards");
                    int actions = transition_probability[i].size();
                    if(actions<=0) throw MDP_exception(8,"actions<=0");
                    for(int j=0;j<actions;j++){
                        if(transition_probability[i][j].size()!=nums_states) throw MDP_exception(8,"transition_probability[i][j].size()!=nums_states");
                        if(rewards[i][j].size()!=nums_states) throw MDP_exception(7,"shape of transition probability != shape of rewards");
                        double sum=0;
                        for(int k=0;k<nums_states;k++){
                            if(transition_probability[i][j][k]<0||transition_probability[i][j][k]>1) throw MDP_exception(9,"transition_probability[i][j][k]<0||transition_probability[i][j][k]>1");
                            sum+=transition_probability[i][j][k];
                        }
                        if(sum-1<-10e-8||sum-1>10e-8) throw MDP_exception(10,"sum(transition_probability[i][j])!=1");//由于浮点数相等经常会出错误，所以判断是否和为1要设定精度
                    }
                }
                //通过以上合法性检查之后MDP对象就是合法的，可以进行赋值
                states=nums_states;
                this->T=T;
                gamma=discount_factor;
                P=transition_probability;
                R=rewards;
                //创建初始策略矩阵
                policy=vector<vector<int>>(T,vector<int>(states,0));
                //创建初始价值函数矩阵
                value_function=vector<vector<double>>(T,vector<double>(states,0));
                policy_evaluation();
            }
        }
        /*其他成员函数暂不展示*/
};

上面的policy_evalution函数就是用来求出给定策略的价值函数，接下来我们按照Bellman等式求出所有状态价值函数

void policy_evaluation(){
        //DP
        for(int state=0;state<states;state++){
            double& v=value_function[T-1][state];
            int a=policy[T-1][state];
            v=0;
            //Expectation
            for(int j=0;j<states;j++) v+=P[state][a][j]*R[state][a][j];
        }
        for(int t=T-2;t>=0;t--){
            for(int s=0;s<states;s++){
                double& v=value_function[t][s];
                int a=policy[t][s];
                v=0;
                //Bellman Equation
                for(int j=0;j<states;j++) v+=P[s][a][j]*(R[s][a][j]+gamma*value_function[t+1][j]);
             }
        }
    }

不难看出上面的算法的时间复杂度是$O(T{S}^2)$，其中$T$是阶段数，$S$是状态数。

状态-行为价值函数(state-action value function) 度量的是在当前状态为$s_t$，采取行动$a_t$，后续行动按照策略$\pi$进行的价值。也就是说，状态-行为价值函数允许代理人在当前时刻偏离策略来采取行动。我们把这个函数记为
$$Q_i^{\pi }(s,a)=E^{\pi }\left[\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})|s_t=s,a_t=a\right]$$
同样地我们可以对状态-行为价值函数也进行相同的推导：
$$\begin{array}{rlr}{Q}_{T-1}^{\pi }(s,a)=& E^{\pi }\left[R(s,a,s_T)|s_t=s,a_t=a\right]=\sum _{s_{t+1}\in \mathcal{S}}\left[p(s_T|s,a)R(s,a,s_T)\right]& \\ Q_i^{\pi }(s,a)=& E^{\pi }\left[R(s,a,s_{i+1}|s_t=s,a_t=a\right]+\gamma E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_t=s,a_t=a\right]& i<T-1\end{array}$$
当$i<T-1$时
$$\begin{array}{rl}& E^{\pi }\left[R(s,a,s_{i+1})|s_t=s,a_t=a\right]=\sum _{s_{i+1}\in \mathcal{S}}\left[p(s_{i+1}|s,a)R(s,a,s_{i+1})\right]\\ & E^{\pi }\left[\sum _{t=i+1}^{T-1}{\gamma }^{t-i-1}R(s_t,a_t,s_{t+1})|s_t=s,a_t=a\right]=\sum _{s_{i+1}\in \mathcal{S}}p(s_{i+1}|s,a)V_{i+1}^{\pi }(s_{i+1})=E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=a,a_i=a\right]\end{array}$$
这样就得到
$$Q_i^{\pi }(s,a)=E^{\pi }\left[R(s,a,s_{i+1})|s_t=s,a_t=a\right]+\gamma E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=a,a_i=a\right]$$
也就是状态-行为价值函数是可以利用行为价值函数计算得到的。我们假定策略是随机策略，记
$${\pi }_i(a|s)=P(a_i=a|s_i=s)$$
于是
$$\begin{array}{rl}{V}_i(s)=& E^{\pi }\left[\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})|s_i=s\right]\\ =& \sum _{a_i,s_{i+1},\cdots ,s_T}{p}^{\pi }(a_i,s_{i+1},\cdots ,s_T|s_i=s)\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})\\ p^{\pi }(a_i,s_{i+1},\cdots ,s_T|s_i=s)=& p^{\pi }(s_i,a_i,s_{i+1},\cdots ,s_T|a_i,s_i=s){\pi }_i(a_i|s)\end{array}$$
代入可以得到
$$\begin{array}{rl}{V}_i^{\pi }(s)=& \sum _{a\in \mathcal{A}(s)}{\pi }_i(a|s)\sum _{s_i=s,a_i=a,s_{i+1},\cdots ,s_T}p(s_i,a_i,s_{i+1},\cdots ,s_T|a_i,s_i=s)\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})\\ =& \sum _{a\in \mathcal{A}(s)}{\pi }_i(a|s)E^{\pi }\left[\sum _{t=i}^{T-1}{\gamma }^{t-i}R(s_t,a_t,s_{t+1})|s_i=s,a_i=a\right]\\ =& \sum _{a\in \mathcal{A}(s)}{\pi }_i(a|s)Q_i(s,a)\end{array}$$
我们就得到两个$Q$和$V$的相互关系式
$$\begin{array}{rl}& Q_i^{\pi }(s,a)=\sum _{s_{i+1}\in \mathcal{S}}\left[p(s_{i+1}|s,a)R(s,a,s_{i+1})\right]+\gamma E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=a,a_i=a\right]\\ & V_i^{\pi }(s)=\sum _{a\in \mathcal{A}(s)}{\pi }_i(a|s)Q_i^{\pi }(s,a)\end{array}$$

最优策略

有了状态价值函数和行为-状态价值函数，我们能否得到一个最优策略呢？要回答这个问题，我们首先要解决一个问题，就是何为"更优"？

如果我们有两个策略${\pi }^{(1)},{\pi }^{(2)}$，那么何谓${\pi }^{(2)}$优于${\pi }^{(1)}$？在一开始代理人就要指定一个策略，确定了一个策略后，进入初始状态$s_0$，初始状态对应不同的策略有不同的状态价值函数$V_0^{{\pi }^{(1)}}(s_0),V_0^{{\pi }^{(2)}}(s_0)$，如果$V_0^{{\pi }^{(1)}}(s_0)\ge V_0^{{\pi }^{(2)}}(s_0)$，说明${\pi }^{(1)}$在初始状态为$s_0$时优于${\pi }^{(2)}$，如果不论什么初始状态，${\pi }^{(1)}$都优于${\pi }^{(2)}$，那么我们就称${\pi }^{(1)}$优于${\pi }^{(2)}$，记为${\pi }^{(1)}\ge {\pi }^{(2)}$，即
$${\pi }^{(1)}\ge {\pi }^{(2)}\iff V_0^{{\pi }^{(1)}}(s)\ge V_0^{{\pi }^{(2)}}(s)\forall s\in \mathcal{S}$$
${\pi }^{\prime }$是最优策略当且仅当${\pi }^{\prime }\ge \pi \forall \pi$

最优策略是否存在呢？答案是肯定的。那如果存在该如何求解呢？显然我们可以求助于状态-行为价值函数$Q$。对于有限期的MDP，我们可以采取分步优化的方式，实际上，我们不难看出$V_i^{\pi }(s)$的计算只决定于子策略$({\pi }_i,{\pi }_{i+1},\cdots ,{\pi }_{T-1})$。我们可以先优化${\pi }_{T-1}$，使得每个$V_{T-1}^{\pi }$都达到最优，然后在优化${\pi }_{T-2}$，使得每个$V_{T-2}^{\pi }$都最优，依次倒推，就可以使得每个$V_0^{\pi }$都达到最优，这样我们就得到一个最优策略。

1. 先优化最后阶段的$V$函数：我们已经有最后一阶段的$Q$函数的计算方法
   $$Q_{T-1}^{\pi }(s,a)=\sum _{s_{t+1}\in \mathcal{S}}\left[p(s_T|s,a)R(s,a,s_T)\right]$$再有$Q$函数和$V$函数的关系
   $$V_{T-1}^{\pi }(s)=\sum _{a_{T-1}\in \mathcal{A}(s)}{\pi }_{T-1}(a|s)Q_{T-1}^{\pi }(s,a)$$可以看出$Q_{T-1}^{\pi }(s,a)$的计算只与模型设定有关，与策略$\pi$无关，那么问题转化成
   $$\begin{array}{rlr}& \underset{{\pi }_{T-1}}{\max }{V}_{T-1}^{\pi }(s)& \\ s.t.& \sum _{a\in \mathcal{A}(s)}{\pi }_{T-1}(a|s)=1& \\ & 0\le {\pi }_{T-1}(a|s)\le 1& \forall a\in \mathcal{A}(s)\end{array}$$那么很显然，使得把所有的权重都分配给$Q$值最大的行为即可，所以
   $${\pi }_{T-1}^{\prime }(a|s)=\{\begin{array}{ll}1& a=argma{x}_{a\in \mathcal{A}(s)[Q_{T-1}^{\pi }(s,a)}]\\ 0& otherwise\end{array}$$
2. 同样地，如果${\pi }_{i+1},\cdots ,{\pi }_{T-1}$都得到了优化，现在要优化${\pi }_i$，我们已经计算了最优的$V_{i+1}^{{\pi }^{\prime }}(s)\forall s\in \mathcal{S}$，我们有$Q$函数和$V$函数的关系，通过$V_{i+1}^{{\pi }^{\prime }}$可以计算得到$Q_i^{\pi }$
   $$Q_i^{\pi }(s,a)=\sum _{s_{i+1}\in \mathcal{S}}\left[p(s_{i+1}|s,a)R(s,a,s_{i+1})\right]+\gamma E^{\pi }\left[V_{i+1}^{\pi }(s_{i+1})|s_i=a,a_i=a\right]$$
   可以看到第一项的计算完全与策略$\pi$无关，第二项只依赖于$V_{i+1}^{{\pi }^{\prime }}$，这样$Q_i^{\pi }$的计算完全与${\pi }_i$无关，再优化
   $$V_i^{\pi }(s)=\sum _{a_{T-1}\in \mathcal{A}(s)}{\pi }_{T-1}(a|s)Q_{i+1}^{\pi }(s,a)$$
   显然也是选择一个使得$Q$函数最大化的行为，然后将所有的权重分配给它即可，也是退化成一个确定性的策略。

总而言之，我们从最后一个阶段向第一个阶段进行优化，就可以得到最优策略${\pi }^{\prime }$，${\pi }^{\prime }$在任何一个阶段的任何状态的价值函数都优于其他策略，所以是最优策略，并且从上面的求解过程也可以看出，最优策略一定是一个确定性的策略，不需要考虑随机策略的情况，这样我们既证明了最优策略的存在性，又证明了最优策略一定是一个确定性的策略，还得到了一个求解最优策略的算法，容易看出这个算法就是动态规划。

有限期MDP最优策略求解算法：

for(int t=T-1;t>=0;t--){
	  for(int s=0;s<states;s++){
        计算所有的行为价值函数Q_t(s,a)
        找一个a使得Q_i(s,a)最大化,policy[t][s]=a
        V_t(s)=Q_t(s,a)
	  }
}

计算单个$Q$函数的复杂度为$O(S)$，计算过程就可以求$Q$的最大值以及使$Q$最大化的行为$a$，于是整个求最优策略的时间复杂度为$O(TAS)$，其中$A$是行为总数（不同状态的同一行为算不同的行为），C++代码如下：

void best_policy(){
            //最后一个状态的策略改进
            for(int s=0;s<states;s++){
                int actions=P[s].size();
                bool f=false;
                for(int a=0;a<actions;a++){
                    double v = 0;
                    for(int j=0;j<states;j++) v+=P[s][a][j]*R[s][a][j];
                    if(!f){ 
	                    value_function[T-1][s]=v;
	                    policy[T-1][s]=a;
	                    f=true;
	                }
                    else if(v>value_function[T-1][s]){
                        value_function[T-1][s]=v;
                        policy[T-1][s]=a;
                    }
                }
            }
            for(int t=T-2;t>=0;t--){
                for(int s=0;s<states;s++){
                    int actions=P[s].size();
                    bool f=false;
                    for(int a=0;a<actions;a++){
                        double v = 0;
                        for(int j=0;j<states;j++) v+=P[s][a][j]*(R[s][a][j]+gamma*value_function[t+1][j]);
                        if(!f){
	                        value_function[t][s]=v;
                            policy[t][s]=a;
                            f=true;
                        }
                        else if(v>value_function[t][s]){
                            value_function[t][s]=v;
                            policy[t][s]=a;
                        }
                    }
                }
            }
        }

对于最优策略${\pi }^{\prime }$，最优策略一定满足Bellman方程，我们记$V_t^{{\pi }^{\prime }}=V_t^{\prime }$，于是就有
$$V_t^{\prime }(s)=\underset{a\in \mathcal{A}(s)}{\max }\left[\sum _{s^{\prime }\in \mathcal{S}}p(s,a,s^{\prime })[R(s,a,s^{\prime })+\gamma V_{t+1}^{\prime }(s^{\prime })]\right]$$
这个方程称为Bellman最优方程，某个策略是最优策略当且仅当其满足Bellman最优方程。

小岛商人问题

现在我们回到第二节给出的小岛商人问题，来看看给定一个策略如何计算所有的状态价值函数以及如何求解最优策略。

我们直接用我们前面已经定义了类MDP

int main(){
	vector<vector<vector<double>>> P={
    {{0.2,0.3,0.5},{0.3,0.3,0.4}},
    {{0.1,0.2,0.7},{0.2,0.1,0.7}},
    {{0.2,0.4,0.4},{0.5,0.3,0.2}}
    };
    vector<vector<vector<double>>> R={
    {{0,2,3},{0,2,3}},
    {{3,0,4},{3,0,4}},
    {{5,3,0},{5,3,0}}
    };
    MDP process(3,0.5,5,0,P,R);
    const vector<vector<double>>& v=process.value();
    const vector<vector<int>>& policy=process.show_policy();
    cout<<"value:"<<endl;
    for(const auto& r:v){
        for(auto v0:r){
            cout<<v0<<' ';
        }
        cout<<endl;
    }
}

运行结果：

按照构造函数，我们初始策略选择的是不论去到哪个岛，就选择第一艘船；现在我们求解最优策略，并给出最优策略的状态价值函数，再加上这段代码，并编译运行

		process.best_policy();//求解最优策略
        cout<<"opt policy:"<<endl;
        for(const auto& r:policy){
            for(auto v0:r){
                cout<<v0<<' ';
            }
            cout<<endl;
        }
        cout<<"value:"<<endl;
        for(const auto& r:v){
            for(auto v0:r){
                cout<<v0<<' ';
            }
            cout<<endl;
        }

运行结果：

最优策略变成如果去了第一个小岛，则坐第一艘船，否则坐第二艘船，我们可以清楚地看到，所有的状态价值函数都得到了提升。