马尔科夫决策过程(MDP)：汽车租赁问题

汽车租赁问题：

杰克是一家汽车租赁公司的经理，他管理着两个汽车租赁点。每天都有一些人在不同的租赁点租车，每出租一辆杰克会得到200￥的奖励，汽车被归还后第二天才能被再次出租。为了尽可能满足广大客户的需求，同时也为了获得更多的现金奖励，杰克计划在两个租赁点之间进行灵活的资源调配，在夜间进行两个租赁点的车辆转移，以满足第二天的租赁需求，两地之间移动车辆的费用是40￥/辆，每次移动车辆的上限是5(辆)，每个租赁点的容量上限是20，当一个租赁点容量达到上限后，车辆会归还到就近的其他租赁点。杰克是个大聪明，他很快就意识到这是一个连续的马尔可夫决策过程，而且状态和动作空间还是有限的。他当即就决定采用强化学习来帮自己找到一个最优的策略，在控制成本的同时，获得更多的收益。

问题抽象和假设

sss : 每个租赁点结束营业后，可用车辆的数量
A\mathcal AA : [-5, 5],正数:表示租赁点1→\to→租赁点2移动的车辆，负数:租赁点2→\to→租赁点1
ttt : time step 天
γ\gammaγ : discount 0.9
假设：每个租赁点的租车需求量和归还量是服从泊松分布的随机变量。
- λnn!e−λ\frac{\lambda^n}{n!}e^{-\lambda}n!λne−λ
- 租赁点1：λ=(3,4)\lambda = (3, 4)λ=(3,4)(需求，归还）
- 租赁点2：λ=(3,2)\lambda = (3, 2)λ=(3,2)(需求，归还）
- 泊松分布的期望和方差均为λ\lambdaλ

模拟求解

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from matplotlib import cm
from scipy.stats import poisson

# 每个租赁点最大容量
MAX_CARS = 20
# 每次最多移动车辆
MAX_MOVE_CARS = 5
# 租赁点1每日租车需求量的期望和方差
RENTAL1_REQUEST = 3
RENTAL2_REQUEST = 4
# 租赁点1每日还车量的期望和方差
RENTAL1_RETURNS = 3
RENTAL2_RETURNS = 2DISCOUNT = 0.9
# 租车奖励
RENTAL_CREDIT = 200
# 转移车辆价格
MOVE_CAR_COST = 40
# 动作空间
actions = np.arange(-MAX_MOVE_CARS, MAX_MOVE_CARS + 1)proba_cache = {} # global

def poisson_proba(n, lam):key = n * 10 + lamif key not in proba_cache:proba_cache[key] = poisson.pmf(n, lam)return proba_cache[key]

动作值函数

qπ(s,a)=∑s′,rp(s′,r∣s,a)(r+vπ(s′))q_{\pi}(s, a) = \sum_{s',r}p(s',r|s,a)\big(r + v_{\pi}(s')\big)qπ(s,a)=s′,r∑p(s′,r∣s,a)(r+vπ(s′))

def expect_return(state, action, state_value, constant_returned_cars, bounds):'''state = (num_cars_rental1, num_cars_rental2)bounds = (num1, num2),request bound, return bound,租还上限'''# init returnreturns = 0# move costreturns -= MOVE_CAR_COST * abs(action)# move carsNUM_CARS_RENTAL1 = min(state[0] - action, MAX_CARS)NUM_CARS_RENTAL2 = min(state[1] + action, MAX_CARS)# 遍历所有的request,return组合for request_num1 in range(bounds[0]):for request_num2 in range(bounds[1]):# 每一种需求组合的概率req_proba = poisson_proba(request_num1, RENTAL1_REQUEST) * poisson_proba(request_num2, RENTAL2_REQUEST)num_cars_rental1 = NUM_CARS_RENTAL1num_cars_rental2 = NUM_CARS_RENTAL2# 实际可满足的需求valid_rental1 = min(num_cars_rental1, request_num1)valid_rental2 = min(num_cars_rental2, request_num2)# rewardreward = (valid_rental1 + valid_rental2) * RENTAL_CREDITnum_cars_rental1 -= valid_rental1num_cars_rental2 -= valid_rental2# 固定的return数量均值if constant_returned_cars:returned_num1 = RENTAL1_RETURNSreturned_num2 = RENTAL2_RETURNS# new statenum_cars_rental1 = min(num_cars_rental1 + returned_num1, MAX_CARS)num_cars_rental2 = min(num_cars_rental2 + returned_num2, MAX_CARS)returns += req_proba * (reward + DISCOUNT * state_value[num_cars_rental1, num_cars_rental2])else:# 每种return组合的概率for returned_num1 in range(bounds[0]):for returned_num2 in range(bounds[1]):ret_proba = poisson_proba(returned_num1, RENTAL1_RETURNS) * poisson_proba(returned_num1, returned_num2)# new statenum_cars_rental1 = min(num_cars_rental1 + returned_num1, MAX_CARS)num_cars_rental2 = min(num_cars_rental2 + returned_num2, MAX_CARS)# p(s'|s,a)proba = req_proba * ret_probareturns += proba * (reward + DISCOUNT * state_value[num_cars_rental1, num_cars_rental2])return returns

策略迭代

Policy Iteration

π0⟶Evπ0⟶Iπ1⟶Evπ1⟶Iπ2⟶Evπ2⟶I...π∗⟶Evπ∗\pi_0 \stackrel{E}{\longrightarrow}v_{\pi_0}\stackrel{I}{\longrightarrow}{\pi_1}\stackrel{E}{\longrightarrow}v_{\pi_1}\stackrel{I}{\longrightarrow}{\pi_2}\stackrel{E}{\longrightarrow}v_{\pi_2}\stackrel{I}{\longrightarrow}...\pi_*\stackrel{E}{\longrightarrow}v_{\pi_*}π0⟶Evπ0⟶Iπ1⟶Evπ1⟶Iπ2⟶Evπ2⟶I...π∗⟶Evπ∗

E:Policy Evaluation

vk+1(s)=Eπ[Rt+1+γvk(St+1)∣St=s]v_{k+1}(s) = E_{\pi}[R_{t+1} + \gamma v_k(S_{t+1})|S_t = s ]vk+1(s)=Eπ[Rt+1+γvk(St+1)∣St=s]

I:Policy Improvement

πk+1(s)=argmaxaqπk(s,a)\pi_{k+1}(s) = \mathop{argmax}\limits_{a}q_{\pi_k}(s,a)πk+1(s)=aargmaxqπk(s,a)

def policy_evalution(value, policy, constant_returned_cars, bounds=(11, 11)):while True:old_value = value.copy()for i in range(MAX_CARS+1):for j in range(MAX_CARS+1):new_state_value = expect_return((i, j), policy[i, j], value, constant_returned_cars, bounds)# inplace updatevalue[i, j] = new_state_valuedelta = abs(old_value - value).max()if delta <= 1e-4:breakreturn value

def policy_improvement(value, policy, constant_returned_cars, bounds=(11, 11)):policy_stable = Truefor i in range(MAX_CARS+1):for j in range(MAX_CARS+1):old_action = policy[i, j]act_returns = []for action in actions:if (0 <= action <= i) or (-j <= action <= 0):act_returns.append(expect_return((i, j), action, value, constant_returned_cars, bounds))else:act_returns.append(-np.inf)new_action = actions[np.argmax(act_returns)]policy[i, j] = new_actionif policy_stable and old_action != new_action:policy_stable = Falsereturn policy, policy_stable

# init value policy
value = np.zeros((MAX_CARS+1, MAX_CARS+1))
policy = np.zeros_like(value, dtype=int)
iterations = 0
value_history = [value.copy()]
policy_history = [policy.copy()]
while True:# policy evaluation, change valuevalue = policy_evalution(value, policy, True)value_history.append(value.copy())# policy improvementpolicy, policy_stable = policy_improvement(value, policy, True)policy_history.append(policy.copy())iterations += 1if policy_stable:break

动作值可视化

策略的可视化

最优动作值