搭建游戏环境
这次我们来用C++实践一下几种算法,我们采用的小游戏是这样的。有一个矿工在一个
n
×
m
n\times m
n×m的网格上进行淘金,每个网格上都有一个数值,是这个格子能掏到的钱的数目,矿工在每一个格子上可以有四种行动,上下左右,决定了行走方向之后,矿工将投一枚骰子来决定走几步,如果走到边界会折返,如下图所示:
红色是矿工所在的初始位置,当矿工决定向上走的时候,投到1、2、3时会走到1、2、3所在的位置,投到4时会折返改变方向向下,4、5、6时分别会走到4、5、6所在的位置。因此可能走到的位置有4个,3、4、5、6号格子,概率分别为 1 6 \frac{1}{6} 61, 1 3 \frac{1}{3} 31, 1 3 \frac{1}{3} 31, 1 6 \frac{1}{6} 61。
游戏的基本设定很简单,下面我们要搭建一个环境,让代理人能够和这个环境进行交互,产生样本。我们先定义一个比较广义的MDP_environment类,这是一个MDP的环境,代理人能和这个environment进行交互,给定状态和行为,环境会返回下一时刻的状态和回报。
其中的interact函数的作用就是,给定state0,action0,返回reward和state1,由于我们要返回的是两个返回值,所以我们在用引用传递的方式进行参数传递。
/*MDP_environment*/
struct MDP_exception{
int exception_id;
string exception_info;
MDP_exception(int id,string info=""):exception_id(id),exception_info(info){}
};
class Environment{
public:
virtual void interact(int state0,int action0,double& reward,int& state1)=0;//交互函数
virtual int actions(int states)=0;//获取states状态的行为数量
virtual int num_states()=0;//获取状态空间的状态数量
};
class MDP_environment:public Environment{
private:
vector<vector<vector<double>>> P,R;
int states;
static bool is_Distribution(const vector<double>& p){
//用来判断p是不是一个分布
double sum=0;
for(auto p_:p){
if(p_<0) return false;
sum+=p_;
}
return abs(sum-1)<=1e-5;
}
static int random(const vector<double>& p){
//给定分布p,产生一个服从分布p的随机数
static default_random_engine e;//Define an random engine
static uniform_real_distribution<double> u(0,1);//Define a distribution
double random_number = u(e);
double p_0=0;
for(int i=0;i<p.size();i++){
if(random_number>=p_0&&random_number<p_0+p[i]) return i;
else p_0+=p[i];
}
return p.size()-1;
}
public:
/*两个构造函数*/
MDP_environment(int states){
if(states<=0) throw MDP_exception(0);
this->states=states;
P=vector<vector<vector<double>>>(states,vector<vector<double>>());
R=vector<vector<vector<double>>>(states,vector<vector<double>>());
}
MDP_environment(const vector<vector<vector<double>>>& P,const vector<vector<vector<double>>>& R){
states=P.size();
if(R.size()!=states) throw MDP_exception(1);
else{
for(int state=0;state<states;state++){
int actions = P[state].size();
if(R[state].size()!=actions) throw MDP_exception(2);
else{
for(int action=0;action<actions;action++){
if(P[state][action].size()!=states||R[state][action].size()!=states) throw MDP_exception(3,"States "+to_string(state)+" action "+to_string(action));
else if(!is_Distribution(P[state][action])) throw MDP_exception(4,"States "+to_string(state)+" action "+to_string(action));
}
}
}
this->P=P;
this->R=R;
}
}
/*交互函数,用来采集样本*/
void interact(int state0,int action0,double& reward,int& state1){
if(state0<0||state0>=states) throw MDP_exception(5);
else if(action0<0||action0>=P[state0].size()) throw MDP_exception(6);
else{
state1=random(P[state0][action0]);
reward=R[state0][action0][state1];
}
}
/*下面是获取环境参数的一些接口*/
int actions(int states){return P[states].size();}
int num_states(){return states;}
void insert_action(int state){
if(state<0||state>=states) throw MDP_exception(7);
else{
vector<double> p(states,0);
p[0]=1;
P[state].push_back(p);
R[state].push_back(vector<double>(states,0));
}
}
void insert_action(int state,const vector<double>& rewards,const vector<double>& probabilities){
if(state<0||state>=states) throw MDP_exception(8);
else if(rewards.size()!=states) throw MDP_exception(9);
else if(probabilities.size()!=states) throw MDP_exception(10);
else if(!is_Distribution(probabilities)) throw MDP_exception(11);
else{
P[state].push_back(probabilities);
R[state].push_back(rewards);
}
}
const vector<vector<vector<double>>>& transition_probability(){return P;}
const vector<vector<vector<double>>>& rewards(){return R;}
void change_transition_probability(const vector<double>& probabilities,int state,int action){
if(state<0||state>=states) throw MDP_exception(12);
else if(action<0||action>=P[state].size()) throw MDP_exception(13);
else if(probabilities.size()!=states) throw MDP_exception(14);
else if(!is_Distribution(probabilities)) throw MDP_exception(15);
else P[state][action]=probabilities;
}
void change_rewards(const vector<double>& rewards,int state,int action){
if(state<0||state>=states) throw MDP_exception(16);
else if(action<0||action>=P[state].size()) throw MDP_exception(17);
else if(rewards.size()!=states) throw MDP_exception(18);
else R[state][action]=rewards;
}
};
我们还定义一个RL_Agent类,表示一个代理人,用来和环境进行交互
class RL_Agent{
private:
Environment* e;
vector<vector<double>> s_policy;
vector<int> d_policy;
bool d;
static bool is_Distribution(const vector<double>& p){
//用来判断p是不是分布
double sum=0;
for(auto p_:p){
if(p_<0) return false;
sum+=p_;
}
return abs(sum-1)<=1e-5;
}
bool legal_stochostic_policy(const vector<vector<double>>& s_policy){
//用来判断给定的随机策略s_policy是不是一个合法的随机策略
if(s_policy.size()!=e->num_states()) return false;
else{
for(int i=0;i<e->num_states();i++){
if(s_policy[i].size()!=e->actions(i)) return false;
if(!is_Distribution(s_policy[i])) return false;
}
return true;
}
}
static int random(const vector<double>& p){
//用来生成服从分布p的一个随机数
static default_random_engine e;//Define an random engine
static uniform_real_distribution<double> u(0,1);//Define a distribution
double random_number = u(e);
double p_0=0;
for(int i=0;i<p.size();i++){
if(random_number>=p_0&&random_number<p_0+p[i]) return i;
else p_0+=p[i];
}
return p.size()-1;
}
public:
//三个构造函数
RL_Agent(Environment& e){
d=true;
d_policy=vector<int>(e.num_states(),0);
this->e=&e;
}
RL_Agent(Environment& e,const vector<int>& d_policy){
d=true;
if(d_policy.size()!=e.num_states()) throw MDP_exception(19);
for(int i=0;i<d_policy.size();i++){
if(d_policy[i]<0||d_policy[i]>=e.actions(i)) throw MDP_exception(20);
}
this->d_policy=d_policy;
this->e=&e;
}
RL_Agent(Environment& e,const vector<vector<double>>& s_policy){
d=false;
if(s_policy.size()!=e.num_states()) throw MDP_exception(21);
for(int i=0;i<s_policy.size();i++){
if(s_policy[i].size()!=e.actions(i)) throw MDP_exception(22);
else if(!is_Distribution(s_policy[i])) throw MDP_exception(23);
}
this->s_policy=s_policy;
this->e=&e;
}
const vector<int>& get_d_policy(){
return d_policy;
}
const vector<vector<double>>& get_s_policy(){
return s_policy;
}
//调整策略的四个函数
void change_policy(const vector<vector<double>>& s_policy){
//将策略调整为随机性策略s_policy
if(s_policy.size()!=e->num_states()) throw MDP_exception(24);
else{
if(!legal_stochostic_policy(s_policy)) throw MDP_exception(25);
this->s_policy=s_policy;
d=false;
}
}
void change_policy(const vector<int>& d_policy){
//将策略调整为确定性策略d_policy
if(d_policy.size()!=e->num_states()) throw MDP_exception(26);
else{
for(int i=0;i<e->num_states();i++){
if(d_policy[i]<0||d_policy[i]>=e->actions(i)) throw MDP_exception(27);
}
s_policy.clear();
this->d_policy=d_policy;
d=true;
}
}
void change_policy(int states,const vector<double>& pi){
//将状态states的随机性策略调整为pi
if(states<0||states>=e->num_states()) throw MDP_exception(28);
else if(pi.size()!=e->actions(states)) throw MDP_exception(29);
else if(!is_Distribution(pi)) throw MDP_exception(30);
else if(d) throw MDP_exception(31);
else s_policy[states]=pi;
}
void change_policy(int states,int actions){
//将状态states的确定性策略调整为actions
if(states<0||states>=e->num_states()) throw MDP_exception(32);
else if(actions<0||actions>=e->actions(states)) throw MDP_exception(33);
else if(!d) throw MDP_exception(34);
else d_policy[states]=actions;
}
int action(int state){
//给定状态state,采取行动
if(state<0||state>=e->num_states()) throw MDP_exception(35);
else if(d) return d_policy[state];
else return random(s_policy[state]);
}
void single_interact(int state0,int& action0,double& reward,int& state1){
//代理人和环境进行单词交互,给定状态state0,采取行为,并同时获得回报和下一时刻的状态
action0=action(state0);
e->interact(state0,action(state0),reward,state1);
}
void trajectory(int state0,int action0,int T,vector<int>& state_series,vector<double>& reward_series,vector<int>& action_series){
//给定初始状态state0,和初始行为action0,产生一条长度为T的轨迹
if(T<=0) throw MDP_exception(36);
state_series.clear();
reward_series.clear();
action_series.clear();
state_series.push_back(state0);
action_series.push_back(action0);
int state1;double reward;
e->interact(state0,action0,reward,state1);
reward_series.push_back(reward);
state0=state1;
for(int t=0;t<T;t++){
state_series.push_back(state0);
action0=action(state0);
action_series.push_back(action0);
e->interact(state0,action0,reward,state1);
reward_series.push_back(reward);
state0=state1;
}
state_series.push_back(state0);
}
void trajectory(int state0,int T,vector<int>& state_series,vector<double>& reward_series,vector<int>& action_series){
//给定初始状态state0,产生一条长度为T的轨迹
state_series.clear();
reward_series.clear();
action_series.clear();
int action0,state1;
double reward;
if(T<=0) throw MDP_exception(36);
for(int t=0;t<T;t++){
state_series.push_back(state0);
action0=action(state0);
action_series.push_back(action0);
e->interact(state0,action0,reward,state1);
reward_series.push_back(reward);
state0=state1;
}
state_series.push_back(state0);
}
};
代理人类最重要的函数就是single_interact和trajectory,分别用来和环境进行交互以及产生一条完整的轨迹,我们后面TD算法主要用single_interact函数来进行采样,而MC算法则用trajectory函数进行采样。
现在,我们来搭建游戏环境,我们搭建一个Find_Money_Game的类用来表示一个矿工游戏,给定一个网格vector<vector<double>>
,自动计算转移概率矩阵和回报矩阵,并搭建一个MDP_environment
。
class Find_Money_Game{
MDP_environment env;
int n,m;
int index(int i,int j){
return i*m+j;
}
void h_generate_probability_and_reward(vector<vector<double>>& R,vector<vector<double>>& P,const vector<vector<double>>& Rewards,int i,int j){
//产生向上走和向下走的转移概率和回报
int i0=i;
bool up=true;
for(int k=1;k<=6;k++){
if(up){
if(i0==0){
up=false;
P[0][index(1,j)]+=1.00/6;
R[0][index(1,j)]=Rewards[1][j];
i0=1;
}else {
P[0][index(--i0,j)]+=1.00/6;
R[0][index(i0,j)]=Rewards[i0][j];
}
}else{
if(i0==n-1){
up=true;
P[0][index(n-2,j)]+=1.00/6;
R[0][index(n-2,j)]=Rewards[n-2][j];
i0=n-2;
}else{
P[0][index(++i0,j)]+=1.00/6;
R[0][index(i0,j)]=Rewards[i0][j];
}
}
}
i0=i,up=false;
for(int k=1;k<=6;k++){
if(up){
if(i0==0){
up=false;
P[1][index(1,j)]+=1.00/6;
R[1][index(1,j)]=Rewards[1][j];
i0=1;
}else {
P[1][index(--i0,j)]+=1.00/6;
R[1][index(i0,j)]=Rewards[i0][j];
}
}else{
if(i0==n-1){
up=true;
P[1][index(n-2,j)]+=1.00/6;
R[1][index(n-2,j)]=Rewards[n-2][j];
i0=n-2;
}else{
P[1][index(++i0,j)]+=1.00/6;
R[1][index(i0,j)]=Rewards[i0][j];
}
}
}
}
void l_generate_probability_and_reward(vector<vector<double>>& R,vector<vector<double>>& P,const vector<vector<double>>& Rewards,int i,int j){
//产生向左走和向右走的转移概率
int j0=j;
bool left=true;
for(int k=1;k<=6;k++){
if(left){
if(j0==0){
left=false;
P[2][index(i,1)]+=1.00/6;
R[2][index(i,1)]=Rewards[i][1];
j0=1;
}else {
P[2][index(i,--j0)]+=1.00/6;
R[2][index(i,j0)]=Rewards[i][j0];
}
}else{
if(j0==m-1){
left=true;
P[2][index(i,m-2)]+=1.00/6;
R[2][index(i,m-2)]=Rewards[i][m-2];
j0=n-2;
}else{
P[2][index(i,++j0)]+=1.00/6;
R[2][index(i,j0)]=Rewards[i][j0];
}
}
}
j0=j,left=false;
for(int k=1;k<=6;k++){
if(left){
if(j0==0){
left=false;
P[3][index(i,1)]+=1.00/6;
R[3][index(i,1)]=Rewards[i][1];
j0=1;
}else {
P[3][index(i,--j0)]+=1.00/6;
R[3][index(i,j0)]=Rewards[i][j0];
}
}else{
if(j0==m-1){
left=true;
P[3][index(i,m-2)]+=1.00/6;
R[3][index(i,m-2)]=Rewards[i][m-2];
j0=n-2;
}else{
P[3][index(i,++j0)]+=1.00/6;
R[3][index(i,j0)]=Rewards[i][j0];
}
}
}
}
public:
Find_Money_Game(const vector<vector<double>>& Rewards):env(5){
n=Rewards.size();
if(n<=1) throw MDP_exception(0);
else{
m=Rewards[0].size();
if(m<=1) throw MDP_exception(1);
else{
for(int i=0;i<n;i++){
if(Rewards[i].size()!=m) throw MDP_exception(2);
}
}
}
vector<vector<vector<double>>> R(n*m,vector<vector<double>>(4,vector<double>(n*m,0)));
vector<vector<vector<double>>> P(n*m,vector<vector<double>>(4,vector<double>(n*m,0)));
for(int i=0;i<n;i++){
for(int j=0;j<m;j++){
//generate R,P
h_generate_probability_and_reward(R[index(i,j)],P[index(i,j)],Rewards,i,j);
l_generate_probability_and_reward(R[index(i,j)],P[index(i,j)],Rewards,i,j);
}
}
env=MDP_environment(P,R);
}
MDP_environment& environment() {return env;}
vector<int> state_to_position(int state){
return {state/m,state%m};
}
};
坐标
(
i
,
j
)
(i,j)
(i,j)表示一个状态,我们用映射
i
∗
m
+
j
i*m+j
i∗m+j映射到一个一维的数组上,
m
m
m是矩阵的列数,state_to_position(state)
函数就可以把状态还原为位置,调用environment()
函数就可以获取生成的这个游戏环境。
现在,我们随机生成一个 n × m n\times m n×m的矩阵,然后再生成一个游戏环境。
int n=20,m=20;
default_random_engine e;
uniform_real_distribution<double> u(-3,3);//回报服从一个-3到3的均匀分布
vector<vector<double>> Rewards(n,vector<double>(m,0));
//随机生成回报
for(auto& row:Rewards){
for(auto& el:row) el=u(e);
}
//搭建游戏环境
Find_Money_Game game(Rewards);
//把搭建的MDP环境提取出来
MDP_environment& env=game.environment();
DP方法
我们先用动态规划方法去求解这个模型,我们把之前的动态规划方法的代码封装成一个DP的类。
class DP{
private:
MDP_environment* environment;//需要求解的MDP环境
double gamma;//折现系数gamma
vector<double> V;//状态价值函数向量V
vector<int> policy;//当前策略policy
void policy_evaluation(const vector<int>& policy,vector<double>& state_value,double e=1e-5){
//策略评估算法的代码
int s=environment->num_states();
const vector<vector<vector<double>>>& P=environment->transition_probability();
const vector<vector<vector<double>>>& R=environment->rewards();
vector<double> R2(s,0);
vector<vector<double>> P2(s,vector<double>(s,0));
//求解线性方程组的系数
for(int i=0;i<s;i++){
for(int k=0;k<s;k++){
P2[i][k]=P[i][policy[i]][k];
R2[i]+=P[i][policy[i]][k]*R[i][policy[i]][k];
}
}
//迭代求解价值函数
double e0=1;
while(e0>=e){
vector<double> v2(s,0);
e0=0;
for(int i=0;i<s;i++){
for(int j=0;j<s;j++) v2[i]+=P2[i][j]*state_value[j];
v2[i]*=gamma;
v2[i]+=R2[i];
e0=max(e0,abs(v2[i]-state_value[i]));
state_value[i]=v2[i];
}
}
}
bool policy_improvement(vector<int>& old_policy,const vector<double>& state_value){
//策略改进算法的代码
const vector<vector<vector<double>>>& P=environment->transition_probability();
const vector<vector<vector<double>>>& R=environment->rewards();
bool i0=false;
int s=environment->num_states();
for(int i=0;i<s;i++){
int new_policy=0;
double maxQ=0;
for(int k=0;k<s;k++) maxQ+=P[i][0][k]*(R[i][0][k]+gamma*state_value[k]);
for(int j=1;j<P[i].size();j++){
double Q=0;
for(int k=0;k<s;k++) Q+=P[i][j][k]*(R[i][j][k]+gamma*state_value[k]);
if(Q>maxQ){
new_policy=j;
maxQ=Q;
}
}
if(new_policy!=old_policy[i]) i0=true;
old_policy[i]=new_policy;
}
return i0;
}
public:
//构造函数
DP(MDP_environment& environment,double gamma=0.5){
this->environment=&environment;
this->gamma=gamma;
V=vector<double>(environment.num_states(),0);
policy=vector<int>(environment.num_states(),0);
}
vector<double> value_function(){
//获得价值函数
return V;
}
void policy_iteration(double e=1e-5){
//策略迭代算法
int s=environment->num_states();
V = vector<double>(s,0);
policy_evaluation(policy,V,e);
while(policy_improvement(policy,V)) {
policy_evaluation(policy,V,e);
}
}
void value_iteration(double e=1e-5){
//价值迭代算法
const vector<vector<vector<double>>>& P=environment->transition_probability();
const vector<vector<vector<double>>>& R=environment->rewards();
int s=environment->num_states();
vector<vector<double>> R2;
//计算系数
for(int i=0;i<s;i++){
R2.push_back(vector<double>(P[i].size(),0));
for(int j=0;j<P[i].size();j++){
for(int k=0;k<s;k++){
R2[i][j]+=R[i][j][k]*P[i][j][k];
}
}
}
//价值迭代
double e0=1;
V = vector<double>(s,0);
vector<double>& value = V;
while(e0>e){
e0=0;
for(int i=0;i<s;i++){
double v=0;
for(int k=0;k<s;k++){
v+=P[i][0][k]*value[k];
}
v*=gamma;
v+=R2[i][0];
policy[i]=0;
for(int j=1;j<P[i].size();j++){
double q=0;
for(int k=0;k<s;k++){
q+=P[i][j][k]*value[k];
}
q*=gamma;
q+=R2[i][j];
if(q>v){
policy[i]=j;
v=q;
}
}
e0=max(e0,abs(v-value[i]));
value[i]=v;
}
}
}
RL_Agent agent(){
//获取一个按当前策略进行行动的代理人
return RL_Agent(*environment,policy);
}
vector<double> policy_evaluation(const vector<int>& policy,double e=1e-5){
//给定确定性策略policy进行策略评估
int s=environment->num_states();
vector<double> state_value(s,0);
RL_Agent agent(*environment,policy);
const vector<vector<vector<double>>>& P=environment->transition_probability();
const vector<vector<vector<double>>>& R=environment->rewards();
vector<double> R2(s,0);
vector<vector<double>> P2(s,vector<double>(s,0));
//求解线性方程组的系数
for(int i=0;i<s;i++){
for(int k=0;k<s;k++){
P2[i][k]=P[i][policy[i]][k];
R2[i]+=P[i][policy[i]][k]*R[i][policy[i]][k];
}
}
//迭代求解价值函数
double e0=1;
while(e0>=e){
vector<double> v2(s,0);
e0=0;
for(int i=0;i<s;i++){
for(int j=0;j<s;j++) v2[i]+=P2[i][j]*state_value[j];
v2[i]*=gamma;
v2[i]+=R2[i];
e0=max(e0,abs(v2[i]-state_value[i]));
state_value[i]=v2[i];
}
}
return state_value;
}
vector<double> policy_evaluation(const vector<vector<double>>& policy,double e=1e-5){
//给定随机策略policy进行策略评估
int s=environment->num_states();
vector<double> state_value(s,0);
RL_Agent agent(*environment,policy);
const vector<vector<vector<double>>>& P=environment->transition_probability();
const vector<vector<vector<double>>>& R=environment->rewards();
vector<double> R2(s,0);
vector<vector<double>> P2(s,vector<double>(s,0));
//求解线性方程组的系数
for(int i=0;i<s;i++){
for(int a=0;a<P[i].size();a++){
for(int k=0;k<s;k++){
P2[i][k]+=P[i][a][k]*policy[i][a];
R2[i]+=P[i][a][k]*R[i][a][k]*policy[i][a];
}
}
}
//迭代求解价值函数
double e0=1;
while(e0>=e){
vector<double> v2(s,0);
e0=0;
for(int i=0;i<s;i++){
for(int j=0;j<s;j++) v2[i]+=P2[i][j]*state_value[j];
v2[i]*=gamma;
v2[i]+=R2[i];
e0=max(e0,abs(v2[i]-state_value[i]));
state_value[i]=v2[i];
}
}
return state_value;
}
};
要求解这个MDP模型的最优策略,用策略迭代算法和价值迭代算法都是可以的,最终得到的最优价值函数和最优策略都一样的。
//策略迭代算法
DP model1(env,0.5);
model1.policy_iteration();
vector<double> value = model1.value_function();
vector<int> policy = model1.agent().get_d_policy();
//价值迭代算法
DP model2(env,0.5);
model2.value_iteration();
vector<double> value2 = model2.value_function();
vector<int> policy2 = model2.agent().get_d_policy();
可以输出两个最优策略和价值函数来看看它们是否有不同。下图是策略迭代的价值函数矩阵(部分):
下图是价值迭代的价值函数矩阵(部分):
我们可以看到两种算法得出的价值函数是一模一样的(在误差不超过0.00001的情况下),下面是两种算法得出的最优策略。下图是策略迭代得到的最优策略:
而下图是价值迭代得到的最优策略
0表示向上走,1表示向下走,2表示向左走,3表示向右走。我们可以看到两种算法得出的最优策略也是一模一样的。
MC方法
预测
我们现在利用MC方法对最优策略进行评估,我们这里采用的方案是First-visit的,也就是对每条轨迹,当前状态第一次出现的时刻才会纳入计算,其他时刻不会纳入计算,这样既能重复利用一条轨迹,也不会破坏样本之间的独立性。
C++代码如下:
vector<double> policy_evaluation(const vector<int>& policy,int T,int N){
RL_Agent agent(e,policy);
int s=e.num_states();
vector<int> states,actions;
vector<double> rewards;
if(T<=0||N<=0) throw MDP_exception(39);
vector<int> N2(s,0);//用来记录每个状态已经被算了多少次
vector<double> V(s,0);//价值函数估计
for(int i=0;i<s;i++){//所有状态都有机会做初始状态
for(int j=0;j<N;j++){
agent.trajectory(i,T,states,rewards,actions);
vector<bool> visit(s,false);//用来记录每个状态是否被访问过
vector<double> V2(T,0);//用来存每个时刻的累积回报函数
V2[T-1]=rewards[T-1];
//倒着计算所有的累积回报函数
for(int t=T-2;t>=0;t--) V2[t]=rewards[t]+gamma*V2[t+1];
for(int t=0;t<T;t++){
if(!visit[states[t]]){
N2[states[t]]++;
V[states[t]]+=(V2[t]-V[states[t]])/N2[states[t]];
visit[states[t]]=true;
}
}
}
/*下面用来输出101号状态,也就是(5,0)的价值函数计算过程
std::cout<<'('<<N2[100]<<',';
std::cout<<V[100]<<"),";
*/
}
return V;
}
我们来看看101号状态经验价值函数在迭代过程中的变化
黄线是状态价值函数的真实值,我们看到,随着经过(5,0)各自次数增多,状态计算越来越准确,最后能够接近真实的状态价值函数。最后所有状态的价值函数估计如下:
估计误差见下图
控制
on policy mc
我们先来实验on policy的MC算法:同样地,我们也采用first visit的方案,我们采用的伪代码如下:
输入:
T:每条轨迹的长度
N:每个状态行为对作为初始状态行为的轨迹数
max_iter:策略改进次数
min_Length:最小轨迹长度
算法:
1. 初始化:epsilon=1,policy[i]=0,N[s][a]=0,Q[s][a]=0
2. 构造epsilon-greedy策略epsilon_greedy
for(int i0=0;i0<max_iter;i0++){
for(int i=0;i<states;i++){
for(int a=0;a<actions(i);a++){
for(int j=0;j<N;j++){
3. visited[i][a]=false;(visited每个值都设为false)
4. 以i为初始状态,a为初始行为抽取一条策略为epsilon_greedy的长度为T的轨迹
5. 从最后一个时刻到第一个时刻,计算每一个时刻的累计回报值
AR[t]=Rewards[t]+gamma*AR[t+1]
6. 从第一个时刻到最后一个时刻计算均值,必须保证到最后一个时刻的时间差至少为min_Length
for(int t=0;t<T-min_Length+1;t++){
if(!visited[states[t]][actions[t]]){
N[states[t]][actions[t]]++;
Q[states[t]][actions[t]]+=(AR[t]-Q[states[t]][actions[t]])/N[states[t]][actions[t]];
}
}
}
}
}
7. 按argmaxQ(a,s)的方式更新策略policy
8. epsilon=1/(i0+2)
9. 按policy构造epsilon-greedy策略
}
输出:最优策略policy
我们来进行简要的讨论:
- 我们这里采用的方案是每个状态行为对都采N条轨迹来计算Q-table,这样做的好处是能够保证Q_table每个格子都有充足的样本量(至少能抽到N个样本来计算Q值),因为我们知道MC方法的基础是大数定理,如果我们随机设定初始状态行为对,那么某些状态行为对的样本量可能会非常少,估计出来的Q函数方差非常大,会使得策略不断的波动,并且也得不到最优策略。
- 我们这里设定了一个最小轨迹长度min_Length,实际上,上节我们的误差分析可以看出,MC的偏差取决于计算Q函数的轨迹长度,轨迹长度越长,偏差越小,因此为了控制偏差,我们尽量不选择轨迹尾部的时刻计算Q函数,这样就能起到误差控制的效果。
- 在这里,我们的epsilon不是一个恒定值,而是一个随着策略改进次数增大不断缩减的值,即 ϵ = 1 , 1 2 , 1 3 , 1 4 , ⋯ \epsilon=1,\frac{1}{2},\frac{1}{3},\frac{1}{4},\cdots ϵ=1,21,31,41,⋯,我们知道,epsilon越大,越趋向于探索,和策略policy偏离的越远,而epsilon越小,越趋向于利用,和policy越靠近,但从采样的角度看,探索偏离策略的状态行为对的机会就变小,可能使得Q函数计算不准确,同样达不到求得最优策略的效果。这是on policy MC面临的探索利用困境,为此,我们第1点已经给出一个在epsilon非常小的情况下也能保证其他状态-行为对有足够样本量的方案,一定程度缓解了探索利用困境。另一方面,在训练一开始,我们当然要倾向于探索,因为在MC的设定下,我们没有关于环境的信息,所以要多探索以全面评估各种状态行为对;但在后期,如果再进行探索,可能会使得我们的epsilon-greedy策略和policy偏离的太远,在后期应当让epsilon-greedy和policy尽可能靠近,使得评估epsilon-greedy一定程度上相当于评估epsilon-greedy策略,这样得出的policy能更靠近最优策略。所以我们采用的epsilon是不断衰减的epsilon。
我们采用的C++代码如下:
class Monte_Carlo{
//其他成员暂略
public:
void on_policy_train(int T,int N,int max_iter=50,int min_Length=20){
if(min_Length<=1) throw MDP_exception(0);
int s=e.num_states();
vector<vector<double>> e_greedy_policy;
epsilon=1;
for(int i=0;i<s;i++){
e_greedy_policy.push_back(vector<double>(e.actions(i)));
for(int a=0;a<e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=epsilon/(e.actions(i))+1-epsilon;
else e_greedy_policy[i][a]=epsilon/(e.actions(i));
}
}
epsilon_greedy.change_policy(e_greedy_policy);
vector<int> states,actions;
vector<double> rewards;
for(int i0=0;i0<max_iter;i0++){
epsilon=1/(i0+2);
vector<vector<double>> Q_table;
vector<vector<int>> N_table;
vector<vector<bool>> visited;
int s=e.num_states();
for(int i=0;i<s;i++) {
Q_table.push_back(vector<double>(e.actions(i),0));
N_table.push_back(vector<int>(e.actions(i),0));
visited.push_back(vector<bool>(e.actions(i),false));
}
for(int i=0;i<s;i++){
for(int a=0;a<e.actions(i);a++){
for(int j=0;j<N;j++){
for(int i=0;i<s;i++) for(int j=0;j<visited[i].size();j++) visited[i][j]=false;
epsilon_greedy.trajectory(i,a,T,states,rewards,actions);
vector<double> q0(T,0);
q0[T-1]=rewards[T-1];
for(int t=T-2;t>=0;t--) q0[t]=rewards[t]+gamma*q0[t+1];
for(int t=0;t<T-min_Length+1;t++){
if(!visited[states[t]][actions[t]]){
N_table[states[t]][actions[t]]++;
Q_table[states[t]][actions[t]]+=(q0[t]-Q_table[states[t]][actions[t]])/N_table[states[t]][actions[t]];
}
}
}
}
}
for(int i=0;i<s;i++){
double maxQ=Q_table[i][0];
policy[i]=0;
for(int a=1;a<e.actions(i);a++){
if(Q_table[i][a]>maxQ){
policy[i]=a;
maxQ=Q_table[i][a];
}
}
}
//construct epsilon greedy policy
vector<vector<double>> e_greedy_policy;
for(int i=0;i<s;i++){
e_greedy_policy.push_back(vector<double>(e.actions(i)));
for(int a=0;a<e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=epsilon/(e.actions(i))+1-epsilon;
else e_greedy_policy[i][a]=epsilon/(e.actions(i));
}
}
epsilon_greedy.change_policy(e_greedy_policy);
}
}
};
我们可以用一张图来表示整个on policy MC的策略改进过程
当然如此改进的前提是在epsilon很小的情况下也能准确评估的Q函数,所以我们引入了上面的第1点的机制,否则在epsilon很小的情况偏离策略状态行为对能被采样的次数非常小,无法准确评估Q函数。这样一来MC算法的计算开销是很大的,因为我们要采样很多的轨迹才能达到准确评估Q函数的目的。
当然,在本例中,通过训练,我们很幸运地策略改进一次就得到比较接近最优解的一个策略,下图是on policy MC改进一次之后得到的策略的价值函数矩阵,这个价值函数矩阵是通过DP方法计算得到的,是这个策略的准确的状态价值函数值。我们可以看到这个价值函数矩阵已经和最优的价值函数矩阵相差无几。
off policy MC
off policy方法直接评估的就是一个确定性策略,但如何解决探索-利用困境呢?off-policy采用另外一个参考策略进行采样,也就是epsilon-greedy策略,然后用参考策略进行采样,通过这采样的轨迹来评估和改进原策略,其原理是在经验累计回报上再乘以一个off policy correction。
在off policy MC下,我们直接评估的是确定性策略的Q函数,此时我们可以不采用epsilon衰减的办法,因为我们直接改进的就是原策略。但是MC方法的基石就是大数定律,因此还是要求每个状态行为对能有大量的样本量以取得Q函数的准确估计;此时策略提升的依据还是MDP的策略提升定理,如果Q函数估计的足够准确,那么能保证每一轮的确定性策略都是提升的。
我们采用的伪代码是
输入:
T:每条轨迹的长度
N:每个状态行为对作为初始状态行为的轨迹数
max_iter:策略改进次数
min_Length:最小轨迹长度
算法:
1. 初始化:policy[i]=0,N[s][a]=0,Q[s][a]=0
2. 构造epsilon-greedy策略epsilon_greedy
for(int i0=0;i0<max_iter;i0++){
for(int i=0;i<states;i++){
for(int a=0;a<actions(i);a++){
for(int j=0;j<N;j++){
3. visited[i][a]=false;(visited每个值都设为false)
4. 以i为初始状态,a为初始行为抽取一条策略为epsilon_greedy的长度为T的轨迹
5. 从最后一个时刻到第一个时刻,计算每一个时刻的累计回报值
rho=rho*(I(policy[states[t]]==actions[t]))/e_greedy_policy[states[t]][actions[t]]
ar=Rewards[t]+gamma*ar;
AR[t]=rho*ar;
6. 从第一个时刻到最后一个时刻计算均值,必须保证到最后一个时刻的时间差至少为min_Length
for(int t=0;t<T-min_Length+1;t++){
if(!visited[states[t]][actions[t]]){
N[states[t]][actions[t]]++;
Q[states[t]][actions[t]]+=(AR[t]-Q[states[t]][actions[t]])/N[states[t]][actions[t]];
}
}
}
}
}
7. 按argmaxQ(a,s)的方式更新策略policy
8. 按policy构造epsilon-greedy策略
}
输出:最优策略policy
我们采用的C++代码如下:
void off_policy_train(int T,int N,int max_iter=50,int min_Length=20){
int s=e.num_states();
vector<vector<double>> e_greedy_policy;
for(int i=0;i<s;i++){
e_greedy_policy.push_back(vector<double>(e.actions(i)));
for(int a=0;a<e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=epsilon/(e.actions(i))+1-epsilon;
else e_greedy_policy[i][a]=epsilon/(e.actions(i));
}
}
epsilon_greedy.change_policy(e_greedy_policy);
vector<int> states,actions;
vector<double> rewards;
for(int i0=0;i0<max_iter;i0++){
vector<vector<double>> Q_table;
vector<vector<int>> N_table;
vector<vector<bool>> visited;
int s=e.num_states();
for(int i=0;i<s;i++) {
Q_table.push_back(vector<double>(e.actions(i),0));
N_table.push_back(vector<int>(e.actions(i),0));
visited.push_back(vector<bool>(e.actions(i),false));
}
for(int i=0;i<s;i++){
for(int a=0;a<e.actions(i);a++){
for(int j=0;j<N;j++){
for(int i=0;i<s;i++) for(int a=0;a<e.actions(i);a++) visited[i][a]=false;
epsilon_greedy.trajectory(i,a,T,states,rewards,actions);
//倒序计算加权回报
vector<double> AR(T,0);
AR[T-1]=rewards[T-1];
double rho=1,ar=AR[T-1];
for(int t=T-2;t>=0;t--){
ar=rewards[t]+gamma*ar;
rho*=((actions[t+1]==policy[states[t+1]])?1:0)/e_greedy_policy[states[t+1]][actions[t+1]];
if(rho==0) break;
else AR[t]=rho*ar;
}
for(int t=0;t<T-min_Length+1;t++){
if(!visited[states[t]][actions[t]]){
N_table[states[t]][actions[t]]++;
Q_table[states[t]][actions[t]]+=(AR[t]-Q_table[states[t]][actions[t]])/N_table[states[t]][actions[t]];
}
}
}
}
}
for(int i=0;i<s;i++){
double maxQ=Q_table[i][0];
policy[i]=0;
for(int a=1;a<e.actions(i);a++){
if(Q_table[i][a]>maxQ){
policy[i]=a;
maxQ=Q_table[i][a];
}
}
}
//construct epsilon greedy policy
for(int i=0;i<s;i++){
for(int a=0;a<e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=epsilon/(e.actions(i))+1-epsilon;
else e_greedy_policy[i][a]=epsilon/(e.actions(i));
}
}
epsilon_greedy.change_policy(e_greedy_policy);
}
}
通过10轮的训练,最终得到的策略的真实价值函数的矩阵如下,也达到了很好的效果。
时序差分
MC算法基于大数定律来进行预测和控制,基本不依赖于对环境的假设,但是基于大数定律进行预测和控制有一个前提:我们要进行大量的采样以保证估计量的精度,否则我们估计出来的Q函数和V函数的准确性无法保证,也不能保证能得到最优的策略。因此,MC算法的计算量是非常大的,是一种离线学习的算法。
时序差分算法TD结合了MC算法和环境的MDP结构,其基本思想是随机逼近,只要满足 ∑ a t = + ∞ , ∑ a t 2 < + ∞ \sum{a_t}=+\infty,\sum{a_t^2}<+\infty ∑at=+∞,∑at2<+∞就能以概率1收敛到真实的V函数和最优的Q函数。
预测
我们的迭代式是:
V
t
+
1
π
(
s
t
)
=
V
t
π
(
s
t
)
+
α
t
(
r
t
+
γ
V
t
π
(
s
t
+
1
)
−
V
t
π
(
s
t
)
)
V_{t+1}^\pi(s_t)=V_{t}^\pi(s_t)+\alpha_t(r_t+\gamma V_t^\pi(s_{t+1})-V_t^\pi(s_t))
Vt+1π(st)=Vtπ(st)+αt(rt+γVtπ(st+1)−Vtπ(st))
现在的问题是
α
t
\alpha_t
αt怎么选取,受到MC算法的启发,MC算法的迭代式是
V
t
+
1
π
(
s
t
)
=
V
t
π
(
s
t
)
+
1
N
t
(
s
t
)
(
G
t
−
V
t
π
(
s
t
)
)
V_{t+1}^\pi(s_t)=V_t^\pi(s_t)+\frac{1}{N_t(s_t)}(G_t-V_t^\pi(s_t))
Vt+1π(st)=Vtπ(st)+Nt(st)1(Gt−Vtπ(st))
其中
G
t
G_t
Gt是累积回报,因此我们可以取
α
t
=
1
N
t
(
s
t
)
\alpha_t=\frac{1}{N_t(s_t)}
αt=Nt(st)1,这样
∑
a
t
=
+
∞
,
∑
a
t
2
<
+
∞
\sum{a_t}=+\infty,\sum{a_t^2}<+\infty
∑at=+∞,∑at2<+∞肯定能得到满足,因为
∑
1
n
=
+
∞
,
∑
1
n
2
=
π
6
\sum{\frac{1}{n}}=+\infty,\sum{\frac{1}{n^2}}=\frac{\pi}{6}
∑n1=+∞,∑n21=6π,于是,其次为了保证每个状态能走足够多次,我们每个状态为初始状态走N条轨迹,走T步进行逐步更新,于是我们的伪代码如下
输入:策略policy
初始化:V[i]=0,N_s[i]=0
for(int i=0;i<s;i++){
1. state0=i;
for(int t=0;t<T;t++){
for(int j=0;j<N;j++){
2.以state0为初始状态,按策略pi和环境交互一次
3.按上面的迭代式进行迭代
N[state0]++
V[state0]+=(Rewards+gamma*V[state1]-V[state0])/N[state0]
}
}
}
返回:V
按照伪代码写出的C++代码如下:
vector<double> policy_evaluation(const vector<int>& policy,int T,int N){
if(T<=0) throw MDP_exception(42);
RL_Agent agent(e,policy);
int s=e.num_states();
vector<double> V(s,0);
vector<int> N_s(s,0);
for(int i=0;i<s;i++){
for(int j=0;j<N;j++){
int state0=i,action0,state1;
double reward;
for(int t=0;t<T;t++){
N_s[state0]++;
agent.single_interact(state0,action0,reward,state1);
V[state0]+=(reward+gamma*V[state1]-V[state0])/N_s[state0];
state0=state1;
}
}
}
return V;
}
下图是V[1]在整个迭代过程中的变化情况:
纵轴是V[1],横轴是经过1号状态的次数,其中的红线是状态价值函数的真实值,随着经过这个状态次数不断增加,最终的价值函数估计V[1]就稳定在红线附近。同样地我们还可以观察399号状态的价值函数变化,也是如此。
控制
随机逼近理论能够证明,在有限状态空间和有限行为空间的情况下,不论是SARSA算法还是Q-learning算法, Q t ( s , a ) Q_t(s,a) Qt(s,a)总会以概率1收敛到 Q ⋆ ( s , a ) Q^\star(s,a) Q⋆(s,a),我们就这两种算法分别进行实验。
SARSA
SARSA的迭代式是
Q
t
+
1
(
s
t
,
a
t
)
←
Q
t
(
s
t
,
a
t
)
+
α
t
(
r
t
+
γ
Q
t
(
s
t
+
1
,
a
t
+
1
)
−
Q
t
(
s
t
,
a
t
)
)
Q_{t+1}(s_t,a_t)\leftarrow Q_t(s_t,a_t)+\alpha_t(r_t+\gamma Q_t(s_{t+1},a_{t+1})-Q_t(s_t,a_t))
Qt+1(st,at)←Qt(st,at)+αt(rt+γQt(st+1,at+1)−Qt(st,at))
其中
a
t
+
1
a_{t+1}
at+1按
ϵ
\epsilon
ϵ-greedy策略采取行动,同时我们每更新一次都要调整一次policy和epsilon-greedy策略,
α
t
=
1
N
t
(
s
t
,
a
t
)
\alpha_t=\frac{1}{N_t(s_t,a_t)}
αt=Nt(st,at)1。因此,伪代码如下:
输入:N,T
初始化:Q[s][a]=0,N[s][a]=0,policy[s]=0,按policy构造epsilon_greedy策略
for(int j=0;j<N;j++){
for(int i=0;i<s;i++){
1.state0=i;
2.按epsilon-greedy策略以state0为初始状态先交互一轮,得到action0,reward0,action1
for(int t=0;t<T;t++){
3.按epsilon-greedy策略以state1为初始状态交互一轮,得到action1,reward1,state2
4.更新Q表:
N[state0][action0]++;
Q[state0][action0]+=(reward0+gamma*Q[state1][action1]-Q[state0][action0])/N[state0][action0];
5.更新policy
6.更新epsilon-greedy策略
7.更新初始状态、行为、回报
state0=state1
action0=action1
reward0=reward1
}
}
}
按照伪代码写的C++代码如下:
void train(int T,int N){
int s=TD_method::e.num_states();
if(T<=0) throw MDP_exception(45);
if(N<=0) throw MDP_exception(46);
//Construct Q-table
vector<vector<double>> e_greedy_policy;
for(int i=0;i<s;i++){
e_greedy_policy.push_back(vector<double>(TD_method::e.actions(i),0));
for(int a=0;a<TD_method::e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=1-epsilon+epsilon/TD_method::e.actions(i);
else e_greedy_policy[i][a]=epsilon/TD_method::e.actions(i);
}
}
RL_Agent ep_greedy_agent(TD_method::e,e_greedy_policy);
vector<vector<double>> Q_table;
vector<vector<int>> N_table;
for(int i=0;i<s;i++){
Q_table.push_back(vector<double>(TD_method::e.actions(i),0));
N_table.push_back(vector<int>(TD_method::e.actions(i),0));
}
for(int j=0;j<N;j++){
for(int s0=0;s0<TD_method::e.num_states();s0++){
int state0=s0,action0,state1,action1,state2;
double reward0,reward1;
ep_greedy_agent.single_interact(state0,action0,reward0,state1);
for(int t=0;t<T;t++){
ep_greedy_agent.single_interact(state1,action1,reward1,state2);
N_table[state0][action0]++;
Q_table[state0][action0]+=(reward0+TD_method::gamma*Q_table[state1][action1]-Q_table[state0][action0])/N_table[state0][action0];
//policy improvement
double maxQ=Q_table[state0][0];
policy[state0]=0;
for(int k=0;k<TD_method::e.actions(state0);k++){
if(Q_table[state0][k]>maxQ){
maxQ=Q_table[state0][k];
policy[state0]=k;
}
}
//update epsilon greedy policy
for(int a=0;a<TD_method::e.actions(state0);a++){
if(a==policy[state0]) e_greedy_policy[state0][a]=1-epsilon+epsilon/TD_method::e.actions(state0);
else e_greedy_policy[state0][a]=epsilon/TD_method::e.actions(state0);
}
ep_greedy_agent.change_policy(state0,e_greedy_policy[state0]);
state0=state1;
action0=action1;
reward0=reward1;
state1=state2;
}
}
}
}
SARSA学到的策略如下:
SARSA学到的策略的真实的状态价值函数如下
Q-learning
Q-learning的迭代式是
Q
t
+
1
(
s
t
,
a
t
)
←
Q
t
(
s
t
,
a
t
)
+
α
t
(
r
t
+
γ
max
a
Q
t
(
s
t
+
1
,
a
)
−
Q
t
(
s
t
,
a
t
)
)
Q_{t+1}(s_t,a_t)\leftarrow Q_t(s_t,a_t)+\alpha_t(r_t+\gamma \max_{a}Q_t(s_{t+1},a)-Q_t(s_t,a_t))
Qt+1(st,at)←Qt(st,at)+αt(rt+γamaxQt(st+1,a)−Qt(st,at))
α
t
=
1
N
t
(
s
t
,
a
t
)
\alpha_t=\frac{1}{N_t(s_t,a_t)}
αt=Nt(st,at)1,这个迭代式相当于每次更新一次Q-table,我们都更新确定性策略policy,然后再迭代式内部按照之前的policy进行行动,是一种off-policy的算法,我们用epsilon-greedy策略进行探索,然后下一步执行的时候用的是确定性策略。
输入:N,T
初始化:Q[s][a]=0,N[s][a]=0,policy[s]=0,按policy构造epsilon_greedy策略
for(int j=0;j<N;j++){
for(int i=0;i<s;i++){
1.state0=i;
for(int t=0;t<T;t++){
2.按epsilon-greedy策略以state0为初始状态交互一轮,得到action0,reward0,state1
3.更新Q表:
N[state0][action0]++;
Q[state0][action0]+=(reward0+gamma*Q[state1][policy[state1]]-Q[state0][action0])/N[state0][action0];
4.更新policy
5.更新epsilon-greedy策略
}
}
}
C++代码如下:
void train(int T,int N){
int s=TD_method::e.num_states();
if(T<=0) throw MDP_exception(45);
//Construct Q-table
vector<vector<double>> e_greedy_policy;
for(int i=0;i<s;i++){
e_greedy_policy.push_back(vector<double>(TD_method::e.actions(i),0));
for(int a=0;a<TD_method::e.actions(i);a++){
if(a==policy[i]) e_greedy_policy[i][a]=1-epsilon+epsilon/TD_method::e.actions(i);
else e_greedy_policy[i][a]=epsilon/TD_method::e.actions(i);
}
}
RL_Agent ep_greedy_agent(TD_method::e,e_greedy_policy);
vector<vector<double>> Q_table;
vector<vector<int>> N_table;
for(int i=0;i<s;i++){
Q_table.push_back(vector<double>(TD_method::e.actions(i),0));
N_table.push_back(vector<int>(TD_method::e.actions(i),0));
}
for(int j=0;j<N;j++){
for(int s0=0;s0<TD_method::e.num_states();s0++){
int state0=s0,action0,state1;
double reward0;
for(int t=0;t<T;t++){
ep_greedy_agent.single_interact(state0,action0,reward0,state1);
N_table[state0][action0]++;
Q_table[state0][action0]+=(reward0+TD_method::gamma*Q_table[state1][policy[state1]]-Q_table[state0][action0])/N_table[state0][action0];
policy[state0]=action0;
double maxQ=Q_table[state0][action0];
for(int a=0;a<TD_method::e.actions(state0);a++){
if(Q_table[state0][a]>maxQ){
policy[state0]=a;
maxQ=Q_table[state0][a];
}
}
//update epsilon greedy policy
for(int a=0;a<TD_method::e.actions(state0);a++){
if(a==policy[state0]) e_greedy_policy[state0][a]=1-epsilon+epsilon/TD_method::e.actions(state0);
else e_greedy_policy[state0][a]=epsilon/TD_method::e.actions(state0);
}
ep_greedy_agent.change_policy(state0,e_greedy_policy[state0]);
state0=state1;
}
}
}
}
Q-learning学到的策略如下:
Q-learning学到的策略的真实价值函数如下:
可以看到Q-learning和SARSA都学到了最优策略,只不过学到的最优策略是不一样的,但其价值函数都等于最优价值函数。