强化学习丨动态规划（DP）学习总结及相关问题的仿真

最新推荐文章于 2023-06-22 10:45:05 发布

XD_MaoHai

最新推荐文章于 2023-06-22 10:45:05 发布

阅读量1.7k

点赞数 10

分类专栏：强化学习文章标签：算法 matlab python

本文链接：https://blog.csdn.net/qq_56937808/article/details/120812857

版权

强化学习专栏收录该内容

6 篇文章 13 订阅

订阅专栏

一、前言

二、动态规划算法（Dynamic Programming）介绍

三、策略迭代（Policy Iteration）

2.1 策略评估（Policy Evaluation）

网格问题

2.2 策略改进（Policy Improvement）

杰克租车问题

四、价值迭代（Value Iteration）

赌徒问题

五、总结

一、前言

本文主要对动态规划算法做以介绍与总结，并利用该算法实现对网格问题、杰克租车问题以及赌徒问题的仿真以及实例分析，其中尤其对于赌徒问题笔者进行了几个拓展思考。本文篇幅较长，含有笔者的自我总结与思考，还望各位批评斧正。文章的开头推荐几个笔者认为不错的博客供各位参考学习：

Jabes：强化学习基础篇（三）动态规划之基础介绍https://www.jianshu.com/p/ac52910f19a1https://www.jianshu.com/p/ac52910f19a1 Jabes：强化学习基础篇（二十三）策略迭代之租车问题https://www.jianshu.com/p/338b173a1dd2https://www.jianshu.com/p/338b173a1dd2 Jabes：强化学习基础篇（二十四）价值迭代之gamblers问题https://www.jianshu.com/p/bfc33ca7b1e5

二、动态规划算法（Dynamic Programming）介绍

在笔者的上篇文章《强化学习丨有限马尔可夫决策过程学习总结》中已经对有限马尔可夫决策过程（有限MDP）及贝尔曼方程的相关概念做了简介，可知对有限MDP问题的求解就是对关于价值函数的贝尔曼方程的求解，而本文要介绍的动态规划算法就是该方程的一个基础解法。

在讲解动态规划具体算法流程前，先介绍《计算方法》中一个常用来求解线性方程组的算法——雅可比迭代法：

设一个线性方程组用矩阵表示为 $Ax = b$ , $A$ 非奇异，且对角元素不为零，则通过移项以及化系数为1可将该方程组改为如下形式：

$\left\{\begin{matrix} x_{1}=\frac{1}{a_{11}}(b_{1}-a_{12}x_{2}-\cdots -a_{1n}x_{n})\\ x_{2}=\frac{1}{a_{22}}(b_{2}-a_{21}x_{1}-\cdots -a_{2n}x_{n})\\ \cdots\cdots\\ x_{n}=\frac{1}{a_{nn}}(b_{n}-a_{n1}x_{1}-\cdots -a_{n,n-1}x_{n-1}) \end{matrix}\right.$

选取初始向量:

$x^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)})$

代入上面方程组右端得：

$x^{(1)}=(x_1^{(1)},x_2^{(1)},\cdots ,x_n^{(1)})$

又将上式代入到方程组右端，得到新的解序列，如此反复继续，得到迭代格式：

$\left\{\begin{matrix} x_{1}^{(k+1)}=\frac{1}{a_{11}}(b_{1}-a_{12}x_{2}^{(k)}-\cdots -a_{1n}x_{n}^{(k)})\\ x_{2}^{(k+1)}=\frac{1}{a_{22}}(b_{2}-a_{21}x_{1}^{(k)}-\cdots -a_{2n}x_{n}^{(k)})\\ \cdots\cdots\\ x_{n}^{(k+1)}=\frac{1}{a_{nn}}(b_{n}-a_{n1}x_{1}^{(k)}-\cdots -a_{n,n-1}x_{n-1}^{(k)}) \end{matrix}\right.$

当迭代前后的解序列变化的最大差值小于要求精度θ时迭代即可结束，这就是雅可比迭代法的过程。

雅可比迭代算法又叫做同步算法，因为该过程需要两个数组，一个存储旧值，一个存储新值，只有进行一次完整的迭代后，新的解序列才会同步更新。但注意到，在一次迭代中，前面几个方程的解已经进行了更新，后面的方程还在用着旧值更新，而为了加快收敛速度，不妨将前面已经更新的解应用到后面方程的更新中，这也是求解线性方程组的另个一种算法——高斯-赛德尔迭代法，也叫做就地更新算法，其迭代公式如下：

$\left\{\begin{matrix} x_{1}^{(k+1)}=\frac{1}{a_{11}}(b_{1}-a_{12}x_{2}^{(k)}-\cdots -a_{1n}x_{n}^{(k)})\\ x_{2}^{(k+1)}=\frac{1}{a_{22}}(b_{2}-a_{21}x_{1}^{(k+1)}-\cdots -a_{2n}x_{n}^{(k)})\\ \cdots\cdots\\ x_{n}^{(k+1)}=\frac{1}{a_{nn}}(b_{n}-a_{n1}x_{1}^{(k+1)}-\cdots -a_{n,n-1}x_{n-1}^{(k+1)}) \end{matrix}\right.$

回到动态规划算法，其核心思想是使用价值函数来结构化的组织对最优策略的搜索，通过将贝尔曼方程转化成为近似逼近理想价值函数的递归更新公式，通过高斯-赛德尔的就地迭代更新算法对其进行精度范围内的求解，这就是动态规划算法。

三、策略迭代（Policy Iteration）

求解有限MDP问题的目的就是寻找最优策略，以达到最高的收益期望，那么这个寻找过程就必然包括了对现阶段下策略的评估（策略评估），以及对该策略的改进（策略改进），在《强化学习》一书中，两者相互作用过程的频率成为“粒度”，而两者不断相互作用与迭代更新的过程就称为策略迭代，不妨将该过程用下面的链式结构进行表示：

$\pi_0\overset{E}{\rightarrow }v_{\pi_0}\overset{I}{\rightarrow }\pi_{1}\overset{E}{\rightarrow }v_{\pi_1}\overset{I}{\rightarrow }\pi_{2}\overset{E}{\rightarrow }\cdots\overset{I}{\rightarrow }\pi_{*}\overset{E}{\rightarrow }v_{*}$

由于一个有限MDP必然只有有限个策略，所以在有限次的迭代后，这种方法一定可以收敛到一个最优的策略与最优价值函数中。

2.1 策略评估（Policy Evaluation）

当确定了一个策略后，我们希望计算该策略下状态的价值函数，这一过程就叫做策略评估。回顾上篇文章介绍的关于状态价值函数的状态方程，对于任意的 $s\in S$ :

$v_{\pi}(s)\doteq \sum_{a}\pi(a|s)\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{\pi}(s^{'})]$

假设一个有限MDP模型中有n个状态，则上式则是一个有着n个未知数以及n个等式联系的线性方程组，利用动态规划算法对该方程组进行迭代求解即可得到迭代公式为：

$v_{k+1}(s)\doteq \sum_{a}\pi(a|s)\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{k}(s^{'})]$

只要折扣因子 $\gamma < 1$ 或者任何状态在 $\pi$ 下都能保证最后终止，那么 $v_{k}$ 在 $k{\rightarrow}\infty$ 时将会收敛到该策略下的状态价值函数 $v_{\pi}$ ，为了减少不必要的迭代次数，我们设置一个收敛精度θ，当更新的解序列与旧序列值之差的最大值Δ小于θ时，迭代结束。以下给出策略评估算法过程的伪代码：

Step1:输入待评估的策略 $\pi$

Step2:输入要求的迭代收敛精度θ，初始化Δ为大于θ的任意实数，初始化状态价值函数 $v(s)$ 使其均为0

Step3:当 $\Delta >\theta$ 时，进行Step4，否则结束程序输出 $v$ 即为 $v_{\pi}$

Step4:记录旧值 $v_{old}=v$ ，对每一个状态 $s\in S$ 进行循环价值更新：

$v(s)= \sum_{a}\pi(a|s)\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v(s^{'})]$

计算误差：

$\Delta=max(|v-v_{old}|)$

返回Step3

现举例网格问题对该算法进行实例应用：

网格问题：

现有一个4×4的网格图（图源自《强化学习》一书，非读者原创）：

非终止状态集合为 $S=\{{}1,2,\cdots,14\}$ ，每个状态有四种可能的动作， $A=\{up,down,right,left\}$ ,每个动作都会导致状态的转移，但当动作会导致智能体移出网格时，状态保持不变，即继续呆在原处。由于存在终止状态（图中阴影），所以这是一个无折扣的分幕式任务，即 $\gamma =1$ 。假设智能体采用等概率随机策略，即每个动作的选择概率均为1/4，在到达终止状态之前，所有的动作收益均为-1。现假设终止状态价值为0，求解其余状态价值函数。

由题意不难得到，对任意的 $s\in S,a\in A$ ，都有：

$\pi(a|s)=1/4,p(s^{'},-1|s,a)=1$

由此不难得出策略评估的迭代更新公式：

$v(s) = -1+\frac{1}{4}(v(s_{1})+v(s_{2})+v(s_{3})+v(s_{4}))$

继而不难写出MATLAB代码如下：
clc,clear;          %清屏
delta = 1e-03;      %估计精度
values = zeros(4,4);%初始化价值矩阵
error = 1;          %初始误差
iteN = 0;           %迭代次数
while(error>delta)
    values_st = values;  %保留原价值矩阵
    for m = 1:4          %行遍历
        for n = 1:4      %列遍历
            flag = m + n;%判断是否为终止状态
            if(flag<=2 || flag>=8)
                continue;
            else
                values(m,n) = -1+1/4*(values(m-sign(m-4),n)+values(m,n-sign(n-4))+values(m-sign(m-1),n)+values(m,n-sign(n-1)));%价值矩阵更新
            end
        end
    end
    iteN = iteN + 1; %迭代次数递增
    error = max(max(abs(values-values_st)));%更新误差
end
最后在88次迭代后，可得到收敛到精度为 $10^{-3}$ 以内的状态价值函数矩阵（保留一位小数结果）：

0.0 -14.0 -20.0 -22.0
-14.0 -18.0 -20.0 -20.0
-20.0 -20.0 -18.0 -14.0
-22.0 -20.0 -14.0 0.0

2.2 策略改进（Policy Improvement）

2.1节已经阐述了策略评估的概念，而策略评估的目的就是以策略的状态价值函数为标尺在每个状态下寻找相应最优的动作，整个状态集合对应的最优动作的集合就是要寻找的最优策略，而这一寻找的过程就叫做策略改进。

如果一个有限MDP的动态函数以及收益分布已经给出，我们利用策略评估很容易能够评估一个状态中某个特定动作的改变会产生什么后果，我们将之延伸到所有的状态和所有可能的动作，在每个状态下根据动作价值函数 $q_{\pi}(s,a)$ 选择一个最优的动作用来更新原策略（若存在多个最优动作可随机选取其一，但随机选择下可能会出现多个动作之间不断切换而算法无法收敛终止的情况，因此最好选择动作标号较小的最优动作），即：

$\pi^{'}(s)\doteq \underset{a}{argmax}q_{\pi}(s,a)=\underset{a}{argmax}\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{\pi}(s^{'})]$

这样构造出的贪心策略满足策略改进定理：

$q_{\pi}(s,\pi^{'}(s))\geq v_{\pi}(s)$

完成对策略的一次更新后，即可对新构造的策略 $\pi^{'}$ 进行新的一轮策略评估，而评估后，对任意状态 $s \in S$ ，改进后的策略能得到一样或者更好的期望回报：

$v_{\pi^{'}}(s)\geq v_{\pi}(s)$

此时我们就称策略 $\pi^{'}$ 相比于原策略 $\pi$ 一样好或者更好，也即完成了一次策略改进。

如前所述，一次策略迭代包含了一次策略评估和策略改进，同时前文也给出了两者的链式关系，先对策略迭代的完整具体算法给以下面的伪代码进行说明：

Step1:初始化

设定初始状态价值矩阵 $v_{\pi}$ 为零矩阵，初始任意策略 $\pi(s)\in A(s)$

Step2:策略评估

①输入要求的迭代收敛精度θ，初始化Δ为大于θ的任意实数

② 当 $\Delta >\theta$ 时，进行步骤③，否则结束程序输出更新的价值矩阵 $v_{\pi}$

        ③ 记录旧值 $v_{old}=v_{\pi}$ ，对每一个状态 $s\in S$ 进行循环价值更新：

$v_{\pi}(s)= \sum_{s^{'},r}p(s^{'},r|s,\pi(s))[r+\gamma v_{\pi}(s^{'})]$

             计算误差：

$\Delta=max(|v_{\pi}-v_{old}|)$

               返回②

Step3：策略改进

① $policy_{stable}\leftarrow True$

② 对每一个状态 $s\in S$ 进行循环：

记录原动作（策略）：

$action_{old}\leftarrow\pi(s)$

策略改进：

$\pi(s)\leftarrow \underset{a}{argmax}\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{\pi}(s^{'})]$

          如果改进后的动作不与原动作相同：

$action_{old}\neq\pi(s)$

则：

$policy_{stable}\leftarrow False$

③ 如果 $policy_{stable}$ 为 $True$ ，则说明最优策略已收敛，结束程序返回最后得到的 $\pi$ 与 $v_{\pi}$ ，否则跳转到Step2

现就杰克租车问题对策略迭代算法做以实例解析：

杰克租车问题：

①问题描述

杰克管理一家有两个地点的租车公司。每一天，一些用户会到某个地点租车，如果该地点有车则杰克会获得10美元/辆的收益，如没有车辆可租则失去这一业务。为了保证每辆车在需要的地方可用，杰克会在夜间在两个地点之间移动车辆，而耳洞每辆车子的代价是2美元/辆。假设每个地点租车与还车的数量服从泊松分布，即数量为n的概率为

$Pr(n)=\frac{\lambda^{n}}{n!}e^{-\lambda}$

其中 $\lambda$ 是期望值，假设租车的 $\lambda$ 在两地的值分别为3和4，还车的 $\lambda$ 在两地的值为3和2。为了简化问题，假设任何一个地点不超过20辆车，多余的车会移送到总公司，即会自动消失，并且每天至多移动5辆车。将该问题视为折扣因子 $\gamma$ 为0.9的持续性任务，则杰克应该怎么调动车辆才能获得最大期望收益？

②问题分析

对有限MDP问题的分析首先要确定其状态和动作集合：

我们将一天结束时作为状态节点，两个地点车辆数作为状态，由于每个地点的车辆有0~20辆，因此问题状态数为21×21=441个。

由于最多可调配5辆车，若用正数表示从地点1移入地点2，负数表示从地点2移入地点1，则动作集合为：

$A=\{-5,-4,-3,-2,-1,0,1,2,3,4,5\}$

接下来就是确定动态函数 $p(s^{'},r|s,a)$ :

设地点1租车与还车数量分别为 $n_{1ren},n_{1ret}$ ，地点2租车与还车数量分别为 $n_{2ren},n_{2ret}$ ,则利用全概率公式有：

$p(s^{'},r|s,a)=\sum_{\begin{matrix} n_{1ren}, & n_{1ret}\\ n_{2ren},& n_{2ret} \end{matrix}}p(n_{1ren}, n_{1ret},n_{2ren} ,n_{2ret} |s,a)\cdot p(s^{'},r|s,a,n_{1ren}, n_{1ret},n_{2ren} ,n_{2ret} )$

又因为 $n_{1ren},n_{1ret},n_{2ren},n_{2ret}$ 互相独立且与 $s,a$ 独立，且当 $n_{1ren},n_{1ret},n_{2ren},n_{2ret}$ 确定时，下一状态和即时收益即可确定，因此上式等式右边乘号右边项为1，进而上式可化为：

$\begin{aligned} p(s^{'},r|s,a)&=\sum_{\begin{matrix} n_{1ren}, & n_{1ret}\\ n_{2ren},& n_{2ret} \end{matrix}}p(n_{1ren}, n_{1ret},n_{2ren} ,n_{2ret} )\cdot 1_{s-n_{1ren}+n_{1ret}-n_{2ren} +n_{2ret}=s^{'}}\\&=\sum_{\begin{matrix} n_{1ren}, & n_{1ret}\\ n_{2ren},& n_{2ret} \end{matrix}}p(n_{1ren})\cdot p(n_{1ret})\cdot p(n_{2ren}) \cdot p(n_{2ret}) \cdot 1_{s-n_{1ren}+n_{1ret}-n_{2ren}+n_{2ret}=s^{'}} \end{aligned}$

其中 $1_{predicate}$ 表示指示函数，当predicate为真时指示函数为1，否则为0。

得到动态函数后，也就不难得到策略评估的迭代公式了，整个问题的重点也就迎刃而解了，之后的工作就按照上述策略迭代的步骤来编写代码即可。

③程序实现

（声明：本问题的代码均参照了本文前言中提到的第三个博客，在此非常感谢博主Jabes的帮助！）

A：两个地点每日还车数为定值（MATLAB）

为了减少运算量先得到初步的结果，不妨先将两个地点每日还车的量设定为定值，且为其原本泊松分布的均值（若没有此操作而采用原泊松分布，则MATLAB平台上这个代码要跑好久）。

先编写获得动作价值函数的函数：

function value = get_value(state,action,values)
%     state :现状态下两个地点的车辆数
%     action:从地点1移至地点2的车辆数
%     values:原状态价值矩阵
%% 初始化
    gamma = 0.9;                  %折扣因子
    MAX_CARS = 20;                %最大存有车辆数
    RENTAL_REQUEST_FIRST_LOC  = 3;%地点1每日租车期望 
    RENTAL_REQUEST_SECOND_LOC = 4;%地点2每日租车期望
    RETURNS_FIRST_LOC = 3;        %地点1每日还车期望
    RETURNS_SECOND_LOC = 2;       %地点2每日还车期望
    RENTAL_CREDIT = 10;           %租车单位收益
    MOVE_CAR_COST = 2;            %移动车所需支出
    POISSON_UPPER_BOUND = 10;     %最高租车容限
    value = 0;                    %状态价值初始化
    value = value - MOVE_CAR_COST * abs(action);           %车辆移动产生的负收益
    NUM_OF_CARS_FIRST_LOC  = min(state(1)-action,MAX_CARS);%移车后地点1的车辆数
    NUM_OF_CARS_SECOND_LOC = min(state(2)+action,MAX_CARS);%移车后地点2的车辆数
%% 遍历两地租车与还车数目求解状态价值函数
    %计算租车概率
    for rental_request_first_loc = 0:POISSON_UPPER_BOUND
        for rental_request_second_loc = 0:POISSON_UPPER_BOUND
            num_of_cars_first_loc  = NUM_OF_CARS_FIRST_LOC ;%地点1原本车辆数
            num_of_cars_second_loc = NUM_OF_CARS_SECOND_LOC;%地点2原本车辆数
            %根据泊松分布计算两地租车的联合概率
            prob_rental = poisspdf(rental_request_first_loc,RENTAL_REQUEST_FIRST_LOC)*poisspdf(rental_request_second_loc,RENTAL_REQUEST_SECOND_LOC);
            valid_rental_first_loc  = min(num_of_cars_first_loc , rental_request_first_loc) ;%一天内地点1实际租车数
            valid_rental_second_loc = min(num_of_cars_second_loc, rental_request_second_loc);%一天内地点2实际租车数
            reward = (valid_rental_first_loc + valid_rental_second_loc) * RENTAL_CREDIT;%租车所得收益
            %根据租车量更新两地车辆数
            num_of_cars_first_loc  = num_of_cars_first_loc  - valid_rental_first_loc ;
            num_of_cars_second_loc = num_of_cars_second_loc - valid_rental_second_loc;
            %当每日还车量一定时
            num_of_cars_first_loc  = min(num_of_cars_first_loc  + RETURNS_FIRST_LOC, MAX_CARS) ;
            num_of_cars_second_loc = min(num_of_cars_second_loc + RETURNS_SECOND_LOC, MAX_CARS);
            value = value + prob_rental * (reward + gamma * values(num_of_cars_first_loc+1,num_of_cars_second_loc+1));%更新状态价值函数
            %计算还车概率
%             for returned_cars_first_loc = 0:POISSON_UPPER_BOUND
%                 for returned_cars_second_loc = 0:POISSON_UPPER_BOUND
%                         %根据泊松分布计算两地还车的联合概率
%                         prob_returned = poisspdf(returned_cars_first_loc,RETURNS_FIRST_LOC)*poisspdf(returned_cars_second_loc,RETURNS_SECOND_LOC);
%                         %根据租车量更新两地车辆数
%                         num_of_cars_first_loc  = min(num_of_cars_first_loc  + returned_cars_first_loc, MAX_CARS) ;
%                         num_of_cars_second_loc = min(num_of_cars_second_loc + returned_cars_second_loc, MAX_CARS);
%                         prob_joint = prob_rental * prob_returned;%计算联合概率，即动态函数
%                         value = value + prob_joint * (reward + gamma * values(num_of_cars_first_loc+1,num_of_cars_second_loc+1));%更新状态价值函数
%                 end
%             end
        end
    end
end

之后编写策略迭代算法的程序：

%% 初始化 
clc,clear;                         %清屏
MAX_CARS = 20;                     %最大存有车辆数
MAX_MOVE_OF_CARS = 5;              %最大可移动车辆数
values  = zeros(MAX_CARS+1,MAX_CARS+1);%状态价值矩阵初始化
policy  = zeros(MAX_CARS+1,MAX_CARS+1);%策略初始化：即不移动车辆
iteN    = 0;                       %迭代次数
delta   = 1e-3;                    %估计精度
error   = 1;                       %初始误差
actions = -MAX_MOVE_OF_CARS:MAX_MOVE_OF_CARS;%动作空间
%% 策略迭代
while(1)
    %策略评估
    while(error>delta)
        values_st = values;  %保留原价值矩阵
        for i = 1:(MAX_CARS+1)
            for j = 1:(MAX_CARS+1)
                %状态价值更新
                value_new   = get_value([i-1,j-1],policy(i,j),values);
                values(i,j) = value_new; 
            end
        end
        error = max(max(abs(values-values_st)));%更新误差
    end
    %策略改进
    policy_stable_flag = 1;%策略稳定标志 
    for i = 1:(MAX_CARS+1)
        for j = 1:(MAX_CARS+1)
            action_old = policy(i,j);%保存原策略下动作
            action_returns = [];     %用以存储各动作的价值函数
            for action = actions     %寻找最优动作
                if(((action<=i-1) && (action>=0)) || ((action<=0) && (action>=-j+1)))
                    %得到各动作的价值函数
                    action_returns = [action_returns get_value([i-1,j-1],action,values)];
                else
                    action_returns = [action_returns -inf];
                end
            end
            %通过最大价值函数确定最优动作
            [returns_max,action_new] = max(action_returns);
            action_new  = action_new-MAX_MOVE_OF_CARS-1;
            policy(i,j) = action_new;
            %根据新旧动作是否一致来判断策略是否稳定
            if(policy_stable_flag && (action_old ~= action_new))
                policy_stable_flag = 0; 
            end
        end
    end
    iteN = iteN + 1;%迭代次数递增
    if(policy_stable_flag)
        break;
    end
end
%% 绘图
figure(1);
imagesc((0:20),(0:20),policy);%绘制最优策略热力图
xlabel('地点2车辆数');ylabel('地点1车辆数');title('最优策略');
figure(2);
imagesc((0:20),(0:20),values);%绘制最优策略下价值函数热力图
xlabel('地点2车辆数');ylabel('地点1车辆数');title('最优策略下状态价值函数');

运行程序可得到最优策略如下：

以及最优策略下状态价值函数如下

B：两个地点每日还车数服从泊松分布（Python）

将两个地点每日还车数从定值改为服从泊松分布，仅需要在获得动作价值函数的函数中添加对两个地点还车数的两重循环即可，但这已经使程序的时间复杂度提升了好多，在MATLAB平台上跑这样的程序会耗费大量时间（笔者亲测在评估精度为 $10^{-3}$ 时1个小时跑不下来），因此改用Python编写程序。

在此再啰嗦一句，看过简书博主Jabes编写的此问题的代码后，深感相形见绌，无出其右，因此就不打算自己编写了，这里直接给出其编写的代码（代码建议看原文，问题分析建议看本文，另原程序选择的是两个地点每日还车数为定值，下面程序选择将其改为泊松分布，即constant_returned_cars=False ）：

导入工具库：

# 载入需要的工具库
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy.stats import poisson

计算泊松分布概率的函数：

# 泊松分布函数
def poisson_probability(n, lam):
    global poisson_cache
    key = n * 10 + lam # 定义唯一key值，除了索引没有实际价值
    if key not in poisson_cache:
        # 计算泊松概率，这里输入为n与lambda，输出泊松分布的概率质量函数，并保存到poisson_cache中
        poisson_cache[key] = poisson.pmf(n, lam)
    return poisson_cache[key]

获得动作价值函数的函数：

# 计算状态价值函数
def expected_return(state, action, state_value, constant_returned_cars):
    """
    @state: 状态定义为每个地点的车辆数

    @action: 车辆的移动数量【-5，5】，负：2->1，正：1->2

    @stateValue: 状态价值矩阵

    @constant_returned_cars:  将还车的数目设定为泊松均值，替换泊松概率分布
    """

    # initialize total return
    returns = 0.0

    # 移动车辆产生负收益
    returns -= MOVE_CAR_COST * abs(action)

    # 移动后的车辆总数不能超过20
    NUM_OF_CARS_FIRST_LOC = min(state[0] - action, MAX_CARS)
    NUM_OF_CARS_SECOND_LOC = min(state[1] + action, MAX_CARS)

    # 遍历两地全部的可能概率下（<11）租车请求数目
    for rental_request_first_loc in range(POISSON_UPPER_BOUND):
        for rental_request_second_loc in range(POISSON_UPPER_BOUND):
            # prob为两地租车请求的联合概率,概率为泊松分布
            # 即：1地请求租车rental_request_first_loc量且2地请求租车rental_request_second_loc量
            prob = poisson_probability(rental_request_first_loc, RENTAL_REQUEST_FIRST_LOC) * \
                   poisson_probability(rental_request_second_loc, RENTAL_REQUEST_SECOND_LOC)

            # 两地原本的车的数量
            num_of_cars_first_loc = NUM_OF_CARS_FIRST_LOC
            num_of_cars_second_loc = NUM_OF_CARS_SECOND_LOC

            # 有效的租车数目必须小于等于该地原有的车辆数目
            valid_rental_first_loc = min(num_of_cars_first_loc, rental_request_first_loc)
            valid_rental_second_loc = min(num_of_cars_second_loc, rental_request_second_loc)

            # 计算回报，更新两地车辆数目变动
            reward = (valid_rental_first_loc + valid_rental_second_loc) * RENTAL_CREDIT
            num_of_cars_first_loc -= valid_rental_first_loc
            num_of_cars_second_loc -= valid_rental_second_loc

            # 如果还车数目为泊松分布的均值
            if constant_returned_cars:
                # 两地的还车数目均为泊松分布均值
                returned_cars_first_loc = RETURNS_FIRST_LOC
                returned_cars_second_loc = RETURNS_SECOND_LOC
                # 还车后总数不能超过车场容量
                num_of_cars_first_loc = min(num_of_cars_first_loc + returned_cars_first_loc, MAX_CARS)
                num_of_cars_second_loc = min(num_of_cars_second_loc + returned_cars_second_loc, MAX_CARS)
                # 核心：
                # 策略评估：V(s) = p(s',r|s,π(s))[r + γV(s')]
                returns += prob * (reward + DISCOUNT * state_value[num_of_cars_first_loc, num_of_cars_second_loc])

            # 否则计算所有泊松概率分布下的还车空间
            else:
                for returned_cars_first_loc in range(POISSON_UPPER_BOUND):
                    for returned_cars_second_loc in range(POISSON_UPPER_BOUND):
                        prob_return = poisson_probability(
                            returned_cars_first_loc, RETURNS_FIRST_LOC) * poisson_probability(returned_cars_second_loc,
                                                                                              RETURNS_SECOND_LOC)
                        num_of_cars_first_loc_ = min(num_of_cars_first_loc + returned_cars_first_loc, MAX_CARS)
                        num_of_cars_second_loc_ = min(num_of_cars_second_loc + returned_cars_second_loc, MAX_CARS)
                        # 联合概率为【还车概率】*【租车概率】
                        prob_ = prob_return * prob
                        returns += prob_ * (reward + DISCOUNT *
                                            state_value[num_of_cars_first_loc_, num_of_cars_second_loc_])
    return returns

策略迭代主函数：

# 策略迭代
def main(constant_returned_cars=False):
    # 初始化价值函数为0
    value = np.zeros((MAX_CARS + 1, MAX_CARS + 1))
    # 初始化策略为0
    policy = np.zeros(value.shape, dtype=np.int)
    # 设置迭代参数
    iterations = 0

    # 准备画布大小，并准备多个子图
    _, axes = plt.subplots(2, 3, figsize=(40, 20))
    # 调整子图的间距，wspace=0.1为水平间距，hspace=0.2为垂直间距
    plt.subplots_adjust(wspace=0.1, hspace=0.2)
    # # 这里将子图形成一个1*6的列表
    axes = axes.flatten()
    while True:
        # 使用seaborn的heatmap作图
        # flipud为将矩阵进行垂直角度的上下翻转，第n行变为第一行，第一行变为第n行，如此。
        # cmap:matplotlib的colormap名称或颜色对象；
        # 如果没有提供，默认为cubehelix map (数据集为连续数据集时) 或 RdBu_r (数据集为离散数据集时)
        fig = sns.heatmap(np.flipud(policy), cmap="Wistia", ax=axes[iterations])
        #
        # 定义标签与标题
        fig.set_ylabel('The num of cars at first location', fontsize=25)
        fig.set_yticks(list(reversed(range(MAX_CARS + 1))))
        fig.set_xlabel('The num of cars at second location', fontsize=25)
        fig.set_title('policy {}'.format(iterations), fontsize=30)

        # policy evaluation (in-place) 策略评估（in-place）
        # 未改进前，第一轮policy全为0，即[0，0，0...]
        while True:
            old_value = value.copy()
            for i in range(MAX_CARS + 1):
                for j in range(MAX_CARS + 1):
                    # 更新V（s）
                    new_state_value = expected_return([i, j], policy[i, j], value, constant_returned_cars)
                    # in-place操作
                    value[i, j] = new_state_value
            # 比较V_old(s)、V(s)，收敛后退出循环
            max_value_change = abs(old_value - value).max()
            print('max value change {}'.format(max_value_change))
            if max_value_change < 1e-2:
                break

        # policy improvement
        # 在上一部分可以看到，策略policy全都是0，如不进行策略改进，其必然不会收敛到实际最优策略。
        # 所以需要如下策略改进
        policy_stable = True
        # i、j分别为两地现有车辆总数
        for i in range(MAX_CARS + 1):
            for j in range(MAX_CARS + 1):
                old_action = policy[i, j]
                action_returns = []
                # actions为全部的动作空间，即[-5、-4...4、5]
                for action in actions:
                    if (0 <= action <= i) or (-j <= action <= 0):
                        action_returns.append(expected_return([i, j], action, value, constant_returned_cars))
                    else:
                        action_returns.append(-np.inf)
                # 找出产生最大动作价值的动作
                new_action = actions[np.argmax(action_returns)]
                # 更新策略
                policy[i, j] = new_action
                if policy_stable and old_action != new_action:
                    policy_stable = False
        print('policy stable {}'.format(policy_stable))

        if policy_stable:
            fig = sns.heatmap(np.flipud(value), cmap="Wistia", ax=axes[-1])
            fig.set_ylabel('The num of cars at first location', fontsize=25)
            fig.set_yticks(list(reversed(range(MAX_CARS + 1))))
            fig.set_xlabel('The num of cars at second location', fontsize=25)
            fig.set_title('V(s) of Optimal Policy', fontsize=30)

            break
        iterations += 1
        print('Now the iteration is  {}'.format(iterations))

        plt.savefig('ResultFigure')

if __name__ == '__main__':
    main()

运行程序后得到迭代过程的策略变化以及最优策略下的状态价值函数热力图如下：

（注意上图中纵轴值从上至下依次减小）

四、价值迭代（Value Iteration）

注意到在策略迭代的算法中，每一次迭代都有一次策略评估，这是一个需要多次遍历状态集合的迭代过程，并且策略评估和策略改进是分开进行的，这样策略评估时选择的动作都是上个迭代周期中策略改进过程中算得到的最优动作，而在本迭代周期中前面的状态价值更新后，尚未更新的状态对应的最优动作可能会改变，对应的价值函数也会改变。

因此现在我们考虑减小策略评估和策略改进的粒度，即将两者融合，在一次迭代中，在根据现有状态价值函数选择最优动作的基础上同时利用该最优动作的价值对状态的价值进行更新，这种算法就称为价值迭代，其算法伪代码如下：

Step1： 设定初始状态价值矩阵 $v$ 为零矩阵，输入要求的迭代收敛精度θ，初始化Δ为大于θ的任意实数，初始任意策略 $\pi(s)\in A(s)$

Step2：  当 $\Delta >\theta$ 时，进行Step3，否则结束程序返回最优价值函数 $v_{\pi}$ 并转入Step4

Step3：记录旧值 $v_{old}=v_{\pi}$ ，对每一个状态 $s\in S$ 进行循环：

                  $v_{\pi}(s)= \underset{a}{max} \sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{\pi}(s^{'})]$



              计算误差：

$\Delta=max(|v_{\pi}-v_{old}|)$

              返回Step2

Step4：输出最优策略：

$\pi(s)= \underset{a}{argmax} \sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v_{\pi}(s^{'})]$

现就赌徒问题对价值迭代算法进行实例分析：

赌徒问题：

①问题描述

一个赌徒利用赌资进行赌注，每次所下赌注是少于自己赌资的一个整数，赢则收获和赌注一样的钱数，败则失去所压之资。当其钱数积累到100时即退出赌局，完成任务，赌徒收益只有完成任务时为1，其余情况均为0，因此此时状态价值函数给出赌徒完成任务的概率。分析每盘胜率为不同的 $P_h(0<P_h<1)$ 情况下各个状态，也即赌徒赌资对应的价值函数以及完成任务的最优策略。

②问题分析

根据题意不难确定这是一个非折扣的可分幕式任务，其状态集合为：

$S=\{0,1,2,3,\cdots ,100\}$

其中 $s=0$ 与 $s=100$ 都是终止状态，其价值分别为0，1。每个状态对应的动作集合为：

$a\in \{0,1,\cdots,min(s,100-s)\}$

由于任务无动作收益，每个状态采取一定的动作后其邻接状态只有两个（输钱和赢钱），因此不难确定动态函数，继而写出其价值迭代公式如下：

$\begin{aligned} v(s)&= \underset{a}{max}\sum_{s^{'},r}p(s^{'},r|s,a)[r+\gamma v(s^{'})]\\ &=\underset{a}{max}[P_{h}\cdot v(s+a)+(1-P_{h})\cdot v(s-a)] \end{aligned}$

接下来就将上式代入到价值迭代公式之中完成对状态价值的更新即可。

③程序实现（MATLAB）

先编写得到动作价值函数的函数：
function value = get_value(state,action,values)
%   state :现状态赌徒赌资
%   action:赌徒投注数
%   values:原状态价值矩阵
%% 初始化
    P_h    = 0.4 ;%正面朝上（获胜）概率
%% 价值评估
    value = P_h*values(state+action+1) + (1-P_h)*values(state-action+1);
end
之后编写价值迭代部分：
%% 初始化
clc,clear    ;%清屏
delta  = 1e-9;%迭代精度
error  = 1   ;%初始化误差
values = zeros(1,101);%价值函数初始化
values(101) = 100      ;%将状态100价值设置为1
policy = zeros(1,99) ;%各个状态对应的策略空间初始化
iteN = 0;
%% 价值迭代
figure(1);
while(error>delta)
    values_old = values();
    for state = 1:99
        action_return = [];
        for action = 1:(min([state 100-state]))
            action_return = [action_return get_value(state,action,values)];
        end
        values(state+1) = max(action_return);%更新价值矩阵
    end
    plot((0:100),values);
    grid on;hold on;
    error = max(abs(values_old-values));%更新误差
    iteN = iteN + 1; %更新迭代次数
end
xlabel('赌资');ylabel('价值估计');title('各次迭代的价值估计曲线');
hold off;
最后迭代达到要求精度后输出稳定最优策略：
%% 最优策略 
figure(2);
optimal_set = zeros(99,50)-1;%各状态对应的最优动作矩阵
for state = 1:99
        action_return = [];  %动作价值向量
        for action = 1:(min([state 100-state]))
            action_return = [action_return get_value(state,action,values)];
        end
        [values_max,policy(state)] = max(round(action_return,5));
        %当某个状态有多个最优动作价值函数时，选择投注数比较小的动作。 
end 
%% 绘制最优策略图
scatter((1:99),policy,'filled');
xlabel('赌资');ylabel('赌注');title('最优策略');
grid on;
需要注意的是，因为在小数位数比较多时计算机计算过程中可能会出现不确定因素，进行比较时可能会出现选择随机的现象，绘制出的policy的点迹也会比较杂乱，因此这里选择action_return的5位小数值进行比较。另外，这里假设赌徒必须做出一定数目的投资，因此选择其动作下限为1。

④结果分析

运行程序得到几个不同 $P_h$ 值下最优状态价值函数曲线以及最优策略：

(a) $P_{h}=0.25$ :





(b) $P_h = 0.4$ :

(c) $P_h = 0.5$ :

(d) $P_h = 0.7$ :

通过观察上图可得到如下结论：

        随着获胜概率的增高，状态价值曲线的会越往上突，即趋近最高值（=1）的速度会越快。更确切的说，当 $P_h$ 小于0.5时，价值曲线大致为下凹函数,说是大致，是因为可以观察到价值曲线中部分状态(如25，50，75)的价值不满足下凹或者上凸特性，但大体上是却是如此；当 $P_h$ 大于0.5时，价值曲线呈上凸函数；等于0.5时为线性函数。现对其整体凹凸性证明如下：

        最优策略下各状态价值函数关系为：

$v_{*}(s)=P_{h}\cdot v_{*}(s+a)+(1-P_{h})\cdot v_{*}(s-a)$

又因为状态价值函数是单调递增的，将上式带入下式有：

$v_{*}(s)-\frac{v_{*}(s+a)+v_{*}(s-a)}{2}\left\{\begin{matrix} <0 & v_{*}(s)\ \text{is generally concave} & P_h<0.5\\ =0&v_{*}(s)\ \text{is linear} & P_h=0.5\\ >0 & v_{*}(s)\ \text{is convex} & P_h>0.5 \end{matrix}\right.$

注意到上面表达式中，当 $P_h$ 小于0.5时，价值曲线大致为下凹函数，而当 $P_h$ 大于0.5时，价值曲线却为严格的凸函数，等于0.5时为严格的线性函数，以下求解出其显式表达式：

（a） $P_h$ 大于0.5

上面已经给出了 $P_h$ 大于0.5情况下的最优策略，即每次都选择投注数为1。我们先不论有没有其他最优策略，前面已经讲述，最优策略共享同一最优价值曲线，因此我们不妨就用已经得到的该最优策略求解价值函数曲线，此时各状态价值函数关系为：

$v_{*}(s)=P_{h}\cdot v_{*}(s+1)+(1-P_{h})\cdot v_{*}(s-1)$

                不难看出该式满足二阶线性递推关系，特征方程的两特征根为1和 $\frac{1-P_h}{P_h}$ ，因为两特征根大于0.5，因此两特征根不相等，其通解为：

$v_{*}(s)=A\cdot1^{n}+B\cdot(\frac{1-P_h}{P_h})^n$

将 $v_{*}(0)=0,v_{*}(100)=1$ 代入可得：

$v_{*}=\frac{1-(\frac{1-P_h}{P_h})^n}{1-(\frac{1-P_h}{P_h})^{100}}$

有了显式公式，就不难解释其单调性和上凸性了，由于篇幅有限，该证明就留给读者吧。

（b） $P_h$ 等于0.5

同上可得，该情况下最优策略依然包括每次都选择投注数为1，因此其特征根依旧为1和 $\frac{1-P_h}{P_h}$ ，但因为 $P_h$ 等于0.5，两个特征根相等，因此其通解为：

$v_{*}(s)=(An+B)\lambda^n$

                 将 $v_{*}(0)=0,v_{*}(100)=1$ 代入可得：

$v_{*}(s)=\frac{n}{100}$

由此可证明其线性性质。

        当 $P_h$ 小于0.5时，最优策略大致呈锯齿状，在赌资为25，50，75时有凸起。个人猜想这是因为在P $P_h$ 小于0.5的情况下，小投多赌的形式会使获胜概率收敛到 $P_h$ ，在此情况下赌徒完成目标任务的概率会越发的小，因此为了减少赌的次数，他会每次都会赌上能投出的赌资的上限，以期最快完成任务。

  但注意到在部分状态下，最优策略并不是投上最高赌资，这是因为出现了多个最优策略时，MATLAB 的max函数只取序数最小的一个，我们不妨求解出每个状态对应的最优策略空间，实现方法仅将输出稳定最优策略部分的代码稍作修改即可：
%% 最优策略 
figure(2);
optimal_set = zeros(99,50)-1;%各状态对应的最优动作矩阵
for state = 1:99
        action_return = [];  %动作价值向量
        for action = 1:(min([state 100-state]))
            action_return = [action_return get_value(state,action,values)];
        end
        %存储各状态对应的最优动作
        [return_sort,index] = sort(round(action_return,5));
        [return_sort_max,index_max] = max(return_sort);
        for i = index(index_max:end)
           optimal_set(state,i) = 1;%最优动作赋值为1 
        end

end 
%% 绘图
imagesc(flipud(optimal_set'));%绘制状态-最优动作矩阵热力图
运行得到 $P_h$ 小于0.5时对应最优策略的热力图如下：

                 可见， $P_h$ 小于0.5时，最优策略中包含了该状态下可投的最大赌注。注意到除此之外，部分状态还对应其他最优策略，不难想到在最优状态价值曲线收敛情况下，这些策略均是满足各状态价值函数关系方程的解。此外，这些解还关于12.5或者25对称，这个性质背后的数学机理笔者还没弄明白，但猜想是因为状态为25，50，75时价值函数斜率出现较大突变的原因。

         当 $P_h$ 大于0.5时，最优策略均为投注1，笔者猜想这是因为在只有任务完成与否为导向的情况下，不考虑完成速度，则在 $P_h$ 大于0.5时，任务完成的概率会比较大，因此稳扎稳打期望概率一定会收敛到 $P_h$ ，任务也会更大概率完成。

为了进一步研究，我们不妨也画出 $P_h$ 大于0.5时最优策略矩阵热力图：

（a） $P_h=0.6$

（b） $P_h=0.7$

（c） $P_h=0.8$

（d） $P_h=0.9$

通过观察上图可知， $P_h$ 值不断上升时，每个状态对应的最优策略数也会不断上升，并且最大投注数会升高。笔者猜想这是因为随着获胜概率的增高，赌徒的容错率会不断升高，因此相应的投注数也会不断升高。

         当 $P_h$ 等于0.5时，同样画出最优策略矩阵热力图如下：

可以看到，  当 $P_h$ 等于0.5时，状态的所有动作选择空间都是该状态下的最优策略，换言之，该情况下不存在最优策略，不难想到这是状态价值函数曲线是线性的原理。

⑤拓展分析——考虑速度因素

         现假设我们想让赌徒以更快的速度完成任务，则可以为每次状态的转移设置一个恒定损耗，该损耗可理解为每次赌局的参赌费（服务费）。为了贴合实际，将赌资为100时状态的价值设定为100，每次赌局的参赌费为1，则改变迭代公式如下：
value = -1+P_h*values(state+action+1) + (1-P_h)*values(state-action+1);
重新运行程序后有如下不同 $P_h$ 对应的价值迭代曲线和最优策略：

（a） $P_h=0.4$

（b） $P_h=0.5$

（c） $P_h=0.7$

（d） $P_h=0.9$

可见，在加入速度限制因素后， $P_h$ 小于等于0.5的情况下最优策略仅为可投最大赌注数。而 $P_h$ 大于0.5时，后半部分状态选择可投最大赌注数，其余状态成间断离散的形式，但每个状态对应的最优策略会大大减小，这也是在加入每局损耗因子后需要尽快完成任务的原因，而至于其前半部分呈间断离散的状态的数学原因笔者还有待考究。

五、总结

总的来说，动态规划是解决贝尔曼方程的一类迭代算法，其基本过程就是策略评估和策略改进的相互作用的过程，根据两者作用的频率（粒度）可将该算法分为策略迭代算法以及价值迭代算法，而由于两者都等次数遍历了整个状态集，同时状态价值在一次迭代后同步更新，两者也叫做同步动态规划，但是这会导致当问题状态数较多时单次遍历成本较高，。

这时可以选择采用异步动态规划算法，即可以以任意顺序和概率更新状态集，而为了保证算法收敛，需要保证至少把每个状态遍历一次，这样虽不一定意味着减少计算量，但保证了一个算法在改进策略前，不需要陷入漫长而无望的遍历，并且可以通过调整状态遍历顺序使得价值信息有效的在状态间进行传播，也就是异步动态规划算法灵活性优点。

XD_MaoHai

关注

10
点赞
踩
30

收藏

觉得还不错? 一键收藏
打赏
0
评论
强化学习丨动态规划（DP）学习总结及相关问题的仿真

一、前言本文主要对动态规划算法做以介绍与总结，并利用该算法实现对网格问题、杰克租车问题以及赌徒问题的仿真，其中尤其对于赌徒问题笔者进行了几个拓展思考，但毕竟水平有限，不尽人意之处还望各位海涵斧正。文章的开头推荐两个笔者认为不错的博客供各位参考学习：Jabes：强化学习基础篇（三）动态规划之基础介绍https://www.jianshu.com/p/ac52910f19a1Jabes：强化学习基础篇（二十三）策略迭代之租车问题https://www.jianshu.com/p/338b...
复制链接

扫一扫