强化学习——马尔可夫决策过程（MDP）附python代码

焕森

一、马尔可夫过程过程介绍随机过程在某时刻t的状态StS_tSt​通常取决于t时刻之前的状态。状态St+1S_{t+1}St+1​的概率表示为：P(St+1∣S1,...,St)P(S_{t+1}|S_1,...,S_t)P(St+1​∣S1​,...,St​)马尔可夫过程某时刻t的状态只取决于上一个时刻t-1的状态。状态St+1S_{t+1}St+1​的概率表示为：P(St+1∣St)=P(St+1∣S1,...,St)P(S_{t+1}|S_t)=P(S_{t+1}|S_1,...,S_t)P(St+1​∣St​)=P(St+1​∣S1​,...,St​)  马尔可夫过程也被称为马尔可夫链，通常用元组来描述，其中S是有限数量的状态集合，P是状态转移矩阵。假设有n个状态，则S={s1,s2,...,sn}S=\{s_1,s_2,...,s_n\}S={s1​,s2​,...,sn​},P=[P(s1∣s1)⋯P(sn∣s1)⋮⋱⋮P(s1∣sn)⋯P(sn∣sn)]P=⎡⎣⎢⎢P(s1|s1)⋮P(s1|sn)⋯⋱⋯P(sn|s1)⋮P(sn|sn)⎤⎦⎥⎥[P(s1|s1)⋯P(sn|s1)⋮⋱⋮P(s1|sn)⋯P(sn|sn)]P=​P(s1​∣s1​)⋮P(s1​∣sn​)​⋯⋱⋯​P(sn​∣s1​)⋮P(sn​∣sn​)​​  矩阵P中第i行第j列元素P(sj∣si)=P(St+1=sj∣St=si)P(s_j|s_i)=P(S_{t+1}=s_j|S_t=s_i)P(sj​∣si​)=P(St+1​=sj​∣St​=si​)表示从状态sis_isi​转移到状态sjs_jsj​的概率，称P(s′∣s)P(s'|s)P(s′∣s)为状态转移函数。从某个状态出发，到达其他状态的概率和必须为1。即状态转移矩阵P的每一行和为1。示例：S={Si,1≤i≤6}  状态转移矩阵  P=[0.90.100000.500.50000000.600.400000.30.700.20.30.500000001]S=\{S_i,1\lei\le6\}\\\;\\状态转移矩阵\=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢0.90.500000.10000.2000.5000.30000.600.500000.300000.40.701⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥[0.90.100000.500.50000000.600.400000.30.700.20.30.500000001]S={Si​,1≤i≤6}状态转移矩阵P=​0.90.50000​0.10000.20​00.5000.30​000.600.50​0000.300​000.40.701​​二、马尔可夫奖励过程2.1定义  马尔可夫奖励过程由构成，其中：S是有限状态的集合P是状态转移矩阵r是奖励函数，r(s)r(s)r(s)是指转移到该状态时可以获得的奖励期望γ\gammaγ是折扣因子，取值范围为：[0,1][0,1][0,1]。引入折扣因子是因为远期利益具有一定的不确定性，有时希望能尽快获得有些奖励，所以需要对远期利益打一些折扣。接近1则更关注长期的累积奖励，接近0则更关注短期奖励。2.2回报  回报是指从状态StS_tSt​开始，一直到终止状态，所有奖励的衰减之和。具体公式如下：Gt=∑k=0∞γkRt+kG_t=\sum^{\infty}_{k=0}\gamma^kR_{t+k}Gt​=k=0∑∞​γkRt+k​  其中RtR_tRt​是指在时刻t获得的奖励。示例：【状态旁的数字表示进入该状态获得的奖励】新建项目MDP，项目结构如下：在MDP目录下新建文件MRP.py，文件代码如下：importnumpyasnpnp.random.seed(0)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.5defcompute_return(start,chain,gamma):G=0foriinreversed(range(start,len(chain))):G=gamma*G+rewards[chain[i]-1]returnGchain=[1,2,3,6]start=0G=compute_return(start,chain,gamma)print('计算回报为：%s'%G)12345678910111213141516171819202122232425262728运行MRP.py文件，运行结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\MRP.py计算回报为：-2.5进程已结束，退出代码为012342.3价值函数价值：一个状态的期望回报。【就是给回报取个别名叫价值】价值函数：输入为某个状态，输出为该状态的价值。公式为：V(s)=E[Gt∣St=s]V(s)=E[G_t|S_t=s]V(s)=E[Gt​∣St​=s]，等价于V(s)=E[Rt+γV(St+1)∣St=s]V(s)=E[R_t+\gammaV(S_{t+1})|S_t=s]V(s)=E[Rt​+γV(St+1​)∣St​=s]。【即：当前状态的价值=当前状态获得的奖励+折扣因子∗*∗下一个状态的价值。下一个状态不是确定的状态，而是取决于转移矩阵P】一方面，E[Rt∣St=s]=r(s)E[R_t|S_t=s]=r(s)E[Rt​∣St​=s]=r(s)；另一方面，E[γV(St+1)∣St=s]E[\gammaV(S_{t+1})|S_t=s]E[γV(St+1​)∣St​=s]可以从状态s出发的转移概率得到，即：E[γV(St+1)∣St=s]=γ∑s′∈SP(s′∣s)V(s′)E[\gammaV(S_{t+1})|S_t=s]=\gamma\sum_{s'\inS}P(s'|s)V(s')E[γV(St+1​)∣St​=s]=γs′∈S∑​P(s′∣s)V(s′)所以V(s)=r(s)+γ∑s′∈SP(s′∣s)V(s′)V(s)=r(s)+\gamma\sum_{s'\inS}P(s'|s)V(s')V(s)=r(s)+γs′∈S∑​P(s′∣s)V(s′)这就是有名的贝尔曼方程，对每一个状态都成立。马尔可夫奖励过程有n个状态，即S={s1,s2,⋯ ,sn}S=\{s_1,s_2,\cdots,s_n\}S={s1​,s2​,⋯,sn​}，所有状态的价值表示成一个列向量V=[V(s1),V(s2),⋯ ,V(sn)]TV=[V(s_1),V(s_2),\cdots,V(s_n)]^TV=[V(s1​),V(s2​),⋯,V(sn​)]T，同理，奖励函数列向量R=[r(s1),r(s2),⋯ ,r(sn)]TR=[r(s_1),r(s_2),\cdots,r(s_n)]^TR=[r(s1​),r(s2​),⋯,r(sn​)]T。于是，贝尔曼方程写成矩阵的形式：V=R+γPVV=R+\gammaPVV=R+γPV，解得：V=(E−γP)−1RV=(E-\gammaP)^{-1}RV=(E−γP)−1R，其中E是单位矩阵。【解析解复杂度为O(n3)O(n^3)O(n3)，只适合求解小规模的马尔可夫奖励过程，大规模的可以使用DP算法，蒙特卡洛方法，时序差分算法】示例：修改文件MRP.py，文件代码如下：importnumpyasnpnp.random.seed(0)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.5defcompute_return(start,chain,gamma):G=0foriinreversed(range(start,len(chain))):G=gamma*G+rewards[chain[i]-1]returnGdefcompute(P,rewards,gamma,states_num):rewards=np.array(rewards).reshape((-1,1))#将rewards改写为列向量value=np.dot(np.linalg.inv(np.eye(states_num,states_num)-gamma*P),rewards)returnvaluedefcompute_return_sample():chain=[1,2,3,6]start=0G=compute_return(start,chain,gamma)print('计算回报为：%s'%G)defcompute_value_sample():V=compute(P,rewards,gamma,6)print("马尔可夫奖励过程中每个状态的价值分别为：\n",V)compute_value_sample()1234567891011121314151617181920212223242526272829303132333435363738394041424344运行结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\MRP.py马尔可夫奖励过程中每个状态的价值分别为：[[-2.01950168][-2.21451846][1.16142785][10.53809283][3.58728554][0.]]进程已结束，退出代码为012345678910三、马尔可夫决策过程  在马尔可夫奖励过程加上动作，就得到了马尔可夫决策过程（MDP）。MDP由元组构成，其中：S是状态的集合A是动作的集合γ\gammaγ是折扣因子r(s,a)r(s,a)r(s,a)是奖励函数，此时奖励同时取决于状态s和动作aP(s′∣s,a)P(s'|s,a)P(s′∣s,a)是状态转移函数，表示在状态s执行动作a之后到达状态s′s's′的概率3.1基本概念策略π(a∣s)=P(At=a∣St=s)\pi(a|s)=P(A_t=a|S_t=s)π(a∣s)=P(At​=a∣St​=s)，表示在输入状态为s的情况下采取动作a的概率。确定性策略：在每个状态下只输出一个确定性的动作，即只有该动作的概率为1，其他动作的概率为0。随机性策略：在每个状态下输出的是关于动作的概率分布，根据该分布采样得到一个动作。状态价值函数Vπ(s)=Eπ[Gt∣St=s]V^{\pi}(s)=E_{\pi}[G_t|S_t=s]Vπ(s)=Eπ​[Gt​∣St​=s]，表示从状态s出发遵循策略π\piπ能获得的期望回报动作价值函数Qπ(s,a)=Eπ[Gt∣St=s,At=a]Q^{\pi}(s,a)=E_{\pi}[G_t|S_t=s,A_t=a]Qπ(s,a)=Eπ​[Gt​∣St​=s,At​=a]，表示遵循策略π\piπ对当前状态s执行动作a得到的期望回报。关系：状态s的价值等于在该状态下基于策略π\piπ采取所有动作的概率与相应价值的乘积和，即：Vπ(s)=∑a∈Aπ(a∣s)  Qπ(s,a)V^{\pi}(s)=\sum_{a\inA}\pi(a|s)\;Q^{\pi}(s,a)Vπ(s)=a∈A∑​π(a∣s)Qπ(s,a)在状态s下采取动作a的价值等于即时奖励加上经过衰减的所有可能的下一个状态的转移概率与相应价值的乘积，即：Qπ(s,a)=r(s,a)+γ∑s′∈SP(s′∣s,a)  Vπ(s′)Q^{\pi}(s,a)=r(s,a)+\gamma\sum_{s'\inS}P(s'|s,a)\;V^{\pi}(s')Qπ(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)3.2贝尔曼期望方程状态价值函数：Vπ(s)=∑a∈Aπ(a∣s)  Qπ(s,a)=∑a∈Aπ(a∣s)[r(s,a)+γ∑s′∈SP(s′∣s,a)  Vπ(s′)]Vπ(s)=∑a∈Aπ(a|s)Qπ(s,a)=∑a∈Aπ(a|s)[r(s,a)+γ∑s′∈SP(s′|s,a)Vπ(s′)]Vπ(s)=∑a∈Aπ(a|s)Qπ(s,a)=∑a∈Aπ(a|s)[r(s,a)+γ∑s′∈SP(s′|s,a)Vπ(s′)]Vπ(s)​=a∈A∑​π(a∣s)Qπ(s,a)=a∈A∑​π(a∣s)[r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)]​动作价值函数：Qπ(s,a)=r(s,a)+γ∑s′∈SP(s′∣s,a)  Vπ(s′)=r(s,a)+γ∑s′∈SP(s′∣s,a)  ∑a∈Aπ(a′∣s′)  Qπ(s′,a′)Qπ(s,a)=r(s,a)+γ∑s′∈SP(s′|s,a)Vπ(s′)=r(s,a)+γ∑s′∈SP(s′|s,a)∑a∈Aπ(a′|s′)Qπ(s′,a′)Qπ(s,a)=r(s,a)+γ∑s′∈SP(s′|s,a)Vπ(s′)=r(s,a)+γ∑s′∈SP(s′|s,a)∑a∈Aπ(a′|s′)Qπ(s′,a′)Qπ(s,a)​=r(s,a)+γs′∈S∑​P(s′∣s,a)Vπ(s′)=r(s,a)+γs′∈S∑​P(s′∣s,a)a∈A∑​π(a′∣s′)Qπ(s′,a′)​【Vπ(s)V^{\pi}(s)Vπ(s)和Qπ(s,a)Q^{\pi}(s,a)Qπ(s,a)互相代入就可以得到】示例：【白色圆圈表示状态，蓝色圆圈表示动作，动作旁边的数字表示奖励，虚线旁边没有数字表示概率为1】在MDP目录下新建MDP.py文件，文件代码如下：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.5MDP=(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}defjoin(s1,s2):returns1+'-'+s212345678910111213141516171819202122232425262728293031323334将MDP转化为MRP，公式如下：r′(s)=∑a∈Aπ(a∣s)r(s,a)r'(s)=\sum_{a\inA}\pi(a|s)r(s,a)r′(s)=∑a∈A​π(a∣s)r(s,a)p′(s′∣s)=∑a∈Aπ(a∣s)P(s′∣s,a)p'(s'|s)=\sum_{a\inA}\pi(a|s)P(s'|s,a)p′(s′∣s)=∑a∈A​π(a∣s)P(s′∣s,a)示例：修改MDP.py文件，文件代码如下：importnumpyimportnumpyasnpS=['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.5MDP=(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}defjoin(s1,s2):returns1+'-'+s2#从字典中取出对应状态的参数数组#例如：从策略1中取出状态S1的动作概率分布，即[0.5,0.5]defget_state_parameter(x,s):p=[]foriinx.keys():ifi.split('-')[0]==s:p.append(x[i])returnp#计算a,b两个向量的内积defcompute_sum(a,b):a=np.array(a)b=np.array(b)returnnp.inner(a,b)#将MDP的奖励函数转为MRP的奖励函数#思想：对于每个状态，根据策略将所有动作的概率进行加权，得到的奖励和就可以被认为是在MRP中该状态的奖励defR_MDP_to_MRP(R,Pi,S):MRP_R=[]foriinrange(len(S)):MRP_R.append(compute_sum(get_state_parameter(Pi,S[i]),get_state_parameter(R,S[i])))returnMRP_R#计算MDP的转移概率defcompute_P(P,Pi):foriinP.keys()[i]*=Pi[join(i.split('-')[0],i.split('-')[1])]returnP#根据转移概率创建MRP的转移矩阵defset_MRP_P(P,S):MRP_P=np.zeros((len(S),len(S)))foriinP.keys():start_index=S.index(i.split('-')[0])end_index=S.index(i.split('-')[2])MRP_P[start_index][end_index]=P[i]MRP_P[len(S)-1][len(S)-1]=1.0#终止状态设置转移概率returnMRP_P#将MDP的转移函数转为MRP的转移矩阵#思想：对于每个状态转移到其他状态，计算策略的转移概率与状态转移函数的转移概率乘积和作为MRP的转移概率defP_MDP_to_MRP(P,Pi,S):returnset_MRP_P(compute_P(P,Pi),S)1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586使用MRP方法计算每个状态的价值：示例：新建Main.py文件，文件代码如下：fromMDPimport*fromMRPimportcomputeR=R_MDP_to_MRP(R,Pi_1,S)P=P_MDP_to_MRP(P,Pi_1,S)print('使用策略1，将MDP转化为MRP')print('转化后的MRP奖励函数：\n',R)print('\n转化后的MRP状态转移矩阵：\n',P)V=compute(P,R,gamma,len(S))print("\nMDP中每个状态价值分别为\n",V)123456789101112运行Main.py文件，运行结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\Main.py使用策略1，将MDP转化为MRP转化后的MRP奖励函数：[np.float64(-0.5),np.float64(-1.5),np.float64(-1.0),np.float64(5.5),np.float64(0.0)]转化后的MRP状态转移矩阵：[[0.50.50.0.0.][0.50.0.50.0.][0.0.0.0.50.5][0.0.10.20.20.5][0.0.0.0.1.]]MDP中每个状态价值分别为[[-1.22555411][-1.67666232][0.51890482][6.0756193][0.]]进程已结束，退出代码为01234567891011121314151617181920计算在状态S下采取动作A的价值：修改MDP.py文件，文件代码如下：importnumpyimportnumpyasnpS=['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.5MDP=(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}defjoin(s1,s2):returns1+'-'+s2#从字典中取出对应状态的参数数组#例如：从策略1中取出状态S1的动作概率分布，即[0.5,0.5]defget_state_parameter(x,s):p=[]foriinx.keys():ifi.split('-')[0]==s:p.append(x[i])returnp#计算a,b两个向量的内积defcompute_sum(a,b):a=np.array(a)b=np.array(b)returnnp.inner(a,b)#将MDP的奖励函数转为MRP的奖励函数#思想：对于每个状态，根据策略将所有动作的概率进行加权，得到的奖励和就可以被认为是在MRP中该状态的奖励defR_MDP_to_MRP(R,Pi,S):MRP_R=[]foriinrange(len(S)):MRP_R.append(compute_sum(get_state_parameter(Pi,S[i]),get_state_parameter(R,S[i])))returnMRP_R#计算MDP的转移概率defcompute_P(P,Pi)1=P.copy()foriinP1.keys()1[i]*=Pi[join(i.split('-')[0],i.split('-')[1])]returnP1#根据转移概率创建MRP的转移矩阵defset_MRP_P(P,S):MRP_P=np.zeros((len(S),len(S)))foriinP.keys():start_index=S.index(i.split('-')[0])end_index=S.index(i.split('-')[2])MRP_P[start_index][end_index]=P[i]MRP_P[len(S)-1][len(S)-1]=1.0#终止状态设置转移概率returnMRP_P#将MDP的转移函数转为MRP的转移矩阵#思想：对于每个状态转移到其他状态，计算策略的转移概率与状态转移函数的转移概率乘积和作为MRP的转移概率defP_MDP_to_MRP(P,Pi,S):returnset_MRP_P(compute_P(P,Pi),S)#MDP=(0--S,1--A,2--P,3--R,4--gamma)#计算在状态S下采取动作A的价值Qdefcompute_Q(s,a,MDP,V):r=MDP[3][join(s,a)]sum_PV=0foriinrange(len(MDP[0])):p=join(join(s,a),MDP[0][i])ifpinMDP[2].keys():sum_PV+=MDP[2][p]*V[i]returnr+MDP[4]*sum_PV123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100修改Main.py文件，文件代码如下：fromMDPimport*fromMRPimportcomputeR=R_MDP_to_MRP(R,Pi_1,S)P=P_MDP_to_MRP(P,Pi_1,S)print('使用策略1，将MDP转化为MRP')print('转化后的MRP奖励函数：\n',R)print('\n转化后的MRP状态转移矩阵：\n',P)V=compute(P,R,gamma,len(S))print("\nMDP中每个状态价值分别为\n",V)print("\n在状态为S4时采取动作概率前往的价值为：",compute_Q(S[3],A[6],MDP,V))1234567891011121314运行Main.py文件，运行结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\Main.py使用策略1，将MDP转化为MRP转化后的MRP奖励函数：[np.float64(-0.5),np.float64(-1.5),np.float64(-1.0),np.float64(5.5),np.float64(0.0)]转化后的MRP状态转移矩阵：[[0.50.50.0.0.][0.50.0.50.0.][0.0.0.0.50.5][0.0.10.20.20.5][0.0.0.0.1.]]MDP中每个状态价值分别为[[-1.22555411][-1.67666232][0.51890482][6.0756193][0.]]在状态为S4时采取动作概率前往的价值为：[2.15123859]进程已结束，退出代码为012345678910111213141516171819202122  MRP解析解的方法在状态动作集合比较大的时候比较不适用。四、蒙特卡洛方法用蒙特卡洛方法估计MDP中的状态价值  用策略在MDP上采样很多序列，计算从这个状态出发的回报再求其期望即可，公式如下：Vπ(s)=Eπ[Gt∣St=s]≈1N∑i=1NGt(i)V^{\pi}(s)=E_{\pi}[G_t|S_t=s]\approx\frac1N\sum^N_{i=1}G^{(i)}_tVπ(s)=Eπ​[Gt​∣St​=s]≈N1​i=1∑N​Gt(i)​具体过程如下：使用策略π\piπ采样若干条序列：s0(i)⟶a0(i)r0(i),  s1(i)⟶a1(i)r1(i),  ⋯ ,  sT−1(i)⟶aT−1(i)rT−1(i),  sT(i)s^{(i)}_0\stackrel{a^{(i)}_0}{\longrightarrow}r^{(i)}_0,\;s^{(i)}_1\stackrel{a^{(i)}_1}{\longrightarrow}r^{(i)}_1,\;\cdots,\;s^{(i)}_{T-1}\stackrel{a^{(i)}_{T-1}}{\longrightarrow}r^{(i)}_{T-1},\;s^{(i)}_{T}s0(i)​⟶a0(i)​​r0(i)​,s1(i)​⟶a1(i)​​r1(i)​,⋯,sT−1(i)​⟶aT−1(i)​​rT−1(i)​,sT(i)​对每一条序列中的每一时间步t的状态s进行以下操作：更新状态s的计数器N(s)+=1N(s)+=1N(s)+=1更新状态s的总回报M(s)+=GM(s)+=GM(s)+=G每一个状态的价值被估计为回报的期望V(s)=M(s)/N(s)V(s)=M(s)/N(s)V(s)=M(s)/N(s)实际中，采用增量式更新会更好，即：N(s)+=1N(s)+=1N(s)+=1V(s)+=1N(s)(G−V(s))V(s)+=\frac{1}{N(s)}(G-V(s))V(s)+=N(s)1​(G−V(s))示例：新建MonteCarlo.py文件，文件代码如下：importnumpyasnpfromMDPimportjoinclassMonteCarlo:#采样序列@staticmethoddefsample(MDP,Pi,timestep_max,number):S,A,P,R,gamma=MDPepisodes=[]for_inrange(number):episode=[]timestep=0s=S[np.random.randint(len(S)-1)]#随机选择除一个终止状态外的状态作为起点#一次采样【到达终止状态或者到达最大时间步】whiles!=S[len(S)-1]andtimesteprand:a=a_optr=R.get(join(s,a),0)breakrand,temp=np.random.rand(),0#根据状态转移函数得到下一个状态fors_optinS:temp+=P.get(join(join(s,a),s_opt),0)iftemp>rand:s_next=s_optbreakepisode.append((s,a,r,s_next))s=s_nextepisodes.append(episode)returnepisodes#计算价值@staticmethoddefcompute(episodes,V,N,gamma):forepisodeinepisodes:G=0#从后往前计算foriinrange(len(episode)-1,-1,-1)s,a,r,s_next)=episode[i]G=r+gamma*GN[s]+=1V[s]+=(G-V[s])/N[s]12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849修改Main.py文件，文件代码如下：fromMDPimport*fromMRPimportcomputefromMonteCarloimportMonteCarloS=['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.5MDP=(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}defMDP_to_MRP()S,A,P,R,gamma)=MDPR=R_MDP_to_MRP(R,Pi_1,S)P=P_MDP_to_MRP(P,Pi_1,S)print('使用策略1，将MDP转化为MRP')print('转化后的MRP奖励函数：\n',R)print('\n转化后的MRP状态转移矩阵：\n',P)V=compute(P,R,gamma,len(S))print("\nMDP中每个状态价值分别为\n",V)print("\n在状态为S4时采取动作概率前往的价值为：",compute_Q(S[3],A[6],MDP,V))#MDP_to_MRP(MDP,Pi_1)defMC():timestep_max=20num=1000V={'S1':0,'S2':0,'S3':0,'S4':0,'S5':0}N={'S1':0,'S2':0,'S3':0,'S4':0,'S5':0}episodes=MonteCarlo.sample(MDP,Pi_1,timestep_max,num)print("采样前5条序列为：\n")foriinrange(5):ifi>=len(episodes):breakprint("序列%d：%s"%(i+1,episodes[i]))MonteCarlo.compute(episodes,V,N,gamma)print('\n使用蒙特卡洛方法计算MDP的状态价值为\n',V)MC()12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667运行Main.py文件，运行结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\Main.py采样前5条序列为：序列1：[('S1','前往S2',0,'S2'),('S2','前往S3',-2,'S3'),('S3','前往S5',0,'S5')]序列2：[('S4','概率前往',1,'S4'),('S4','前往S5',10,'S5')]序列3：[('S4','前往S5',10,'S5')]序列4：[('S2','前往S1',-1,'S1'),('S1','保持S1',-1,'S1'),('S1','前往S2',0,'S2'),('S2','前往S3',-2,'S3'),('S3','前往S4',-2,'S4'),('S4','前往S5',10,'S5')]序列5：[('S2','前往S3',-2,'S3'),('S3','前往S4',-2,'S4'),('S4','前往S5',10,'S5')]使用蒙特卡洛方法计算MDP的状态价值为{'S1':-1.2250047295780488,'S2':-1.6935084269467038,'S3':0.48519918492538405,'S4':5.97601313133763,'S5':0}进程已结束，退出代码为012345678910111213五、占用度量  定义MDP的初始状态分布为v0(s)v_0(s)v0​(s)，用Ptπ(s)P^{\pi}_t(s)Ptπ​(s)表示采取策略π\piπ使得智能体在t时刻状态为s的概率，所以P0π(s)=V0(s)P^{\pi}_0(s)=V_0(s)P0π​(s)=V0​(s)，然后定义一个策略的状态访问分布：Vπ(s)=(1−γ)∑t=0∞γtPtπ(s)V^{\pi}(s)=(1-\gamma)\sum^{\infty}_{t=0}\gamma^tP^{\pi}_t(s)Vπ(s)=(1−γ)t=0∑∞​γtPtπ​(s)  其中，1−γ1-\gamma1−γ是归一化因子。状态访问概率表示在一个策略下智能体和MDP交互会访问到的状态的分布。其有如下性质：Vπ(s′)=(1−γ)V0(s′)+γ∫P(s′∣s,a)π(a∣s)Vπ(s)dsdaV^{\pi}(s')=(1-\gamma)V_0(s')+\gamma\intP(s'|s,a)\pi(a|s)V^{\pi}(s)dsdaVπ(s′)=(1−γ)V0​(s′)+γ∫P(s′∣s,a)π(a∣s)Vπ(s)dsda【解释说明：∑t=0∞γt    是等比数列求和，即：∑t=0∞γt=lim⁡n→∞1−rn1−r=11−r，所以(1−γ)∑t=0∞γt=1\sum^{\infty}_{t=0}\gamma^t\;\;是等比数列求和，即：\sum^{\infty}_{t=0}\gamma^t=\lim_{n\rightarrow\infty}\frac{1-r^n}{1-r}=\frac{1}{1-r}，所以(1-\gamma)\sum^{\infty}_{t=0}\gamma^t=1t=0∑∞​γt是等比数列求和，即：t=0∑∞​γt=n→∞lim​1−r1−rn​=1−r1​，所以(1−γ)t=0∑∞​γt=1注：这只是解释说明为什么1−γ1-\gamma1−γ是归一化因子，不可认为Vπ(s)=(1−γ)∑t=0∞γtPtπ(s)=1⋅Ptπ(s)=Ptπ(s)V^{\pi}(s)=(1-\gamma)\sum^{\infty}_{t=0}\gamma^tP^{\pi}_t(s)=1\cdotP^{\pi}_t(s)=P^{\pi}_t(s)Vπ(s)=(1−γ)t=0∑∞​γtPtπ​(s)=1⋅Ptπ​(s)=Ptπ​(s)】  此外，定义策略的占用度量为：ρπ(s,a)=Vπ(s)π(a∣s)=(1−γ)∑t=0∞γtPtπ(s)π(a∣s)ρπ(s,a)=Vπ(s)π(a|s)=(1−γ)∑t=0∞γtPπt(s)π(a|s)ρπ(s,a)=Vπ(s)π(a|s)=(1−γ)∑t=0∞γtPtπ(s)π(a|s)ρπ(s,a)​=Vπ(s)π(a∣s)=(1−γ)t=0∑∞​γtPtπ​(s)π(a∣s)​  它表示状态动作对(s,a)(s,a)(s,a)被访问到的概率。它有两个定理：ρπ1=ρπ2⇔π1=π2\rho^{\pi_1}=\rho^{\pi_2}\Leftrightarrow\pi_1=\pi_2ρπ1​=ρπ2​⇔π1​=π2​给定合法的占用度量ρ\rhoρ，可生成该占用度量的唯一策略πρ\pi_{\rho}πρ​：πρ=ρ(s,a)∑a′ρ(s,a′)\pi_{\rho}=\frac{\rho(s,a)}{\sum_{a'}\rho(s,a')}πρ​=∑a′​ρ(s,a′)ρ(s,a)​示例：MDP.py文件新增函数occupancy：#计算状态动作（s,a）出现的频率，以此估计策略的占用度量defoccupancy(episodes,s,a,timestep_max,gamma):rho=0total_times=np.zeros(timestep_max)occur_times=np.zeros(timestep_max)forepisodeinepisodes:foriinrange(len(episode))s_opt,a_opt,r,s_next)=episode[i]total_times[i]+=1ifs==s_optanda==a_optccur_times[i]+=1foriinreversed(range(timestep_max)):iftotal_times[i]:rho+=gamma**i*occur_times[i]/total_times[i]return(1-gamma)*rho123456789101112131415修改Main.py文件，文件代码如下：fromMDPimport*fromMRPimportcomputefromMonteCarloimportMonteCarloS=['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.5MDP=(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}defMDP_to_MRP()S,A,P,R,gamma)=MDPR=R_MDP_to_MRP(R,Pi_1,S)P=P_MDP_to_MRP(P,Pi_1,S)print('使用策略1，将MDP转化为MRP')print('转化后的MRP奖励函数：\n',R)print('\n转化后的MRP状态转移矩阵：\n',P)V=compute(P,R,gamma,len(S))print("\nMDP中每个状态价值分别为\n",V)print("\n在状态为S4时采取动作概率前往的价值为：",compute_Q(S[3],A[6],MDP,V))#MDP_to_MRP(MDP,Pi_1)defMC():timestep_max=20num=1000V={'S1':0,'S2':0,'S3':0,'S4':0,'S5':0}N={'S1':0,'S2':0,'S3':0,'S4':0,'S5':0}episodes=MonteCarlo.sample(MDP,Pi_1,timestep_max,num)print("采样前5条序列为：\n")foriinrange(5):ifi>=len(episodes):breakprint("序列%d：%s"%(i+1,episodes[i]))MonteCarlo.compute(episodes,V,N,gamma)print('\n使用蒙特卡洛方法计算MDP的状态价值为\n',V)#MC()defoccupancy_instance():gamma=0.5timestep_max=1000num=1000s=S[3]a=A[6]episodes_1=MonteCarlo.sample(MDP,Pi_1,timestep_max,num)episodes_2=MonteCarlo.sample(MDP,Pi_2,timestep_max,num)rho_1=occupancy(episodes_1,s,a,timestep_max,gamma)rho_2=occupancy(episodes_2,s,a,timestep_max,gamma)print('策略1对状态动作对(%s,%s)的占用度量为：%f'%(s,a,rho_1))print('策略2对状态动作对(%s,%s)的占用度量为：%f'%(s,a,rho_2))occupancy_instance()123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384运行Main.py文件，结果如下：D:\RL\MDP\.venv\Scripts\python.exeD:\RL\MDP\Main.py策略1对状态动作对(S4,概率前往)的占用度量为：0.109338策略2对状态动作对(S4,概率前往)的占用度量为：0.226629进程已结束，退出代码为012345六、最优策略  强化学习的目标通常是找到一个策略，使得智能体从初始状态出发能获得最多的期望回报。最优策略π∗(s)\pi^*(s)π∗(s)：在MDP中，至少存在一个不差于其他所有策略的策略。当且仅当对任意的状态s都有Vπ(s)≥Vπ′(s)V^{\pi}(s)\geV^{\pi'}(s)Vπ(s)≥Vπ′(s)，记π⪰π′\pi\succeq\pi'π⪰π′，MDP至少存在一个最优策略最优状态价值函数V∗(s)V^*(s)V∗(s)：V∗(s)=max⁡πVπ(s),    ∀s∈SV^*(s)=\max\limits_{\pi}V^{\pi}(s),\;\;\foralls\inSV∗(s)=πmax​Vπ(s),∀s∈S最优动作价值函数Q∗(s,a)Q^*(s,a)Q∗(s,a)：Q∗(s,a)=max⁡πQπ(s,a),    ∀s∈S,a∈AQ^*(s,a)=\max\limits_{\pi}Q^{\pi}(s,a),\;\;\foralls\inS,a\inAQ∗(s,a)=πmax​Qπ(s,a),∀s∈S,a∈A最优状态价值函数与最优动作价值函数的关系：Q∗(s,a)=r(s,a)+γ∑s′∈SP(s′∣s,a)V∗(s′)  V∗(s)=max⁡a∈AQ∗(s,a)Q^*(s,a)=r(s,a)+\gamma\sum_{s'\inS}P(s'|s,a)V^*(s')\\\;\\V^*(s)=\max\limits_{a\inA}Q^*(s,a)Q∗(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)V∗(s′)V∗(s)=a∈Amax​Q∗(s,a)贝尔曼最优方程：Q∗(s,a)=r(s,a)+γ∑s′∈SP(s′∣s,a)max⁡a′∈AQ∗(s′,a′)  V∗(s)=max⁡a∈A{r(s,a)+γ∑s′∈SP(s′∣s,a)V∗(s′)}Q^*(s,a)=r(s,a)+\gamma\sum_{s'\inS}P(s'|s,a)\max\limits_{a'\inA}Q^*(s',a')\\\;\\V^*(s)=\max\limits_{a\inA}\{r(s,a)+\gamma\sum_{s'\inS}P(s'|s,a)V^*(s')\}Q∗(s,a)=r(s,a)+γs′∈S∑​P(s′∣s,a)a′∈Amax​Q∗(s′,a′)V∗(s)=a∈Amax​{r(s,a)+γs′∈S∑​P(s′∣s,a)V∗(s′)}得到最优策略的方法：第4章动态规划算法第5章时序差分算法