Deep reinforcement learning for approximate policy iteration: convergence analysis and a post-earthquake disaster response case study

Abstract

Approximate policy iteration (API) is a class of reinforcement learning (RL) algorithms that seek to solve the long-run discounted reward Markov decision process (MDP), via the policy iteration paradigm, without learning the transition model in the underlying Bellman equation. Unfortunately, these algorithms suffer from a defect known as chattering in which the solution (policy) delivered in each iteration of the algorithm oscillates between improved and worsened policies, leading to sub-optimal behavior. Two causes for this that have been traced to the crucial policy improvement step are: (i) the inaccuracies in the policy improvement function and (ii) the exploration/exploitation tradeoff integral to this step, which generates variability in performance. Both of these defects are amplified by simulation noise. Deep RL belongs to a newer class of algorithms in which the resolution of the learning process is refined via mechanisms such as experience replay and/or deep neural networks for improved performance. In this paper, a new deep learning approach is developed for API which employs a more accurate policy improvement function, via an enhanced resolution Bellman equation, thereby reducing chattering and eliminating the need for exploration in the policy improvement step. Versions of the new algorithm for both the long-run discounted MDP and semi-MDP are presented. Convergence properties of the new algorithm are studied mathematically, and a post-earthquake disaster response case study is employed to demonstrate numerically the algorithm’s efficacy.

Appendix

Proof

(Lemma 1) Note that since $\gamma >0$ and $t(.,.,.)>0$, there exists a $\lambda \in (0,1)$ such that:

$$\begin{aligned} \max _{i,j\in \mathcal {S}; a\in \mathcal {A}(i)}\exp (-\gamma t(i,a,j))\le \lambda . \end{aligned}$$

(11)

From G(.)’s definition in Eq. (6), for any two vectors, $\textbf{J}^k$ and $\mathbf {\overline{J}}^k$:

$$\begin{aligned} G(J^k)(i)-G(\overline{J}^k)(i)=\sum _{j \in S}p(i,\mu (i),j) \exp (-\bar{\gamma t(i,\mu (i),j)}\left[ J^k(j)-\overline{J}^k(j)\right] ~~~\forall i. \end{aligned}$$

$$\begin{aligned} \text{ Then, }\ \forall i:\left| G(J^k)(i)-G(\overline{J}^k)(i)\right|&\le \sum _{j \in S}p(i,\mu (i),j) \exp (-\bar{\gamma t(i,\mu (i),j)}\max _{j\in \mathcal {S}}\left| J^k(j)-\overline{J}^k(j)\right| \\&\le \sum _{j \in \mathcal {S}}p(i,\mu (i),j)\lambda \max _{j \in \mathcal {S}}\left| J^k(j) -\overline{J}^k(j)\right| \text{ from } Eq. (11)\\&= \sum _{j \in \mathcal {S}}p(i,\mu (i),j)\lambda ||\textbf{J}^k-\mathbf {\overline{J}}^k||_{\infty }\\&=\lambda ||\textbf{J}^k-\mathbf {\overline{J}}^k||_{\infty } \sum _{j \in S}p(i,\mu (i),j) =\lambda ||\textbf{J}^k-\mathbf {\overline{J}}^k||_{\infty }. \end{aligned}$$

Then, $\max _{i \in \mathcal {S}}\left| G(J^k)(i)-G(\overline{J}^k)(i)\right| =||G(\textbf{J}^k)-G(\mathbf {\overline{J}}^k)||_{\infty } \le \lambda ||\textbf{J}^k-\mathbf {\overline{J}}^k||_{\infty }$. $\square$

Proof

(Lemma 2) It is claimed that for every $i\in \mathcal {S}$:

$$\begin{aligned} |J^k(i)|\le M(1+\lambda +\lambda ^2+\cdots +\lambda ^k), \end{aligned}$$

(12)

$$\begin{aligned} \text{ where } M=\max \left\{ \max _{i,j\in \mathcal{S}, a\in \mathcal{A}(i)}|r(i,a,j)|,\max _{i\in \mathcal{S}} |J^1(i)|\right\} \end{aligned}$$

and $\lambda$ is defined via Eq. (11). Then, the result follows from the fact that:

$$\begin{aligned} \limsup _{k\rightarrow \infty }|J^k(i)|\le M \frac{1}{1-\lambda } \text{ for } \text{ all } i\in \mathcal {S}. \end{aligned}$$

The claim in (12) is proved via induction. In asynchronous updating, two cases can occur:

Case 1 The value for a state visited in the kth iteration is updated: $J^{k+1}(i)=(1-\eta ^k)J^k(i)+\eta ^k\left[ r(i,a,j)+\exp (-\gamma t(i,a,j))J^k(j)\right] .$

Case 2 The value for a state not visited in the kth iteration is not updated: $J^{k+1}(i)=J^k(i)$.

When the update is carried out as in Case 1:

$$ \begin{aligned} |J^{2} (i)|\, & \le \,(1 - \eta ^{1} )|J^{1} (i)| + \eta ^{1} |r(i,a,j) + \exp ( - \gamma t(i,a,j)J^{1} (j)| \\ & \le (1 - \eta ^{1} )|J^{1} (i)| + \eta ^{1} |r(i,a,j) + \lambda J^{1} (j)|{\text{ from }}(11) \\ & \le (1 - \eta ^{1} )M + \eta ^{1} M + \eta ^{1} \lambda M{\text{ from }}M^{\prime}s{\text{ definition}} \\ & < (1 - \eta ^{1} )M + \eta ^{1} M + \lambda M = M(1 + \lambda ){\text{ since }}\eta ^{1} < 1 \\ \end{aligned} $$

When the update is carried out as in Case 2: $|J^2(i)| = |J^1(i)|\le M \le M(1+\lambda ).$ The claim thus holds for $k=1$. When the claim holds for $k=m$, one has that:

$$\begin{aligned}|J^{m}(i)|\le M(1+\lambda +\lambda ^2+\cdots +\lambda ^{m})~~~~\forall i. \end{aligned}$$

Under Case 1: $\forall i$:

$$ \begin{aligned} |J^{{m + 1}} (i)|\, & \le (1 - \eta ^{m} )|J^{m} (i)| + \eta ^{m} |r(i,a,j) + \exp ( - \gamma t(i,a,j))J^{m} (j)| \\ & \le (1 - \eta ^{m} )|J^{m} (i)| + \eta ^{m} |r(i,a,j) + \lambda J^{m} (j)| \\ & \le (1 - \eta ^{m} )M(1 + \lambda + \lambda ^{2} + \cdots + \lambda ^{m} ) + \eta ^{m} M + \eta ^{m} \lambda M(1 + \lambda + \lambda ^{2} + \cdots + \lambda ^{m} ) \\ & = M(1 + \lambda + \lambda ^{2} + \cdots + \lambda ^{m} ) + \eta ^{m} M\lambda ^{{m + 1}} \\ & \le M(1 + \lambda + \lambda ^{2} + \cdots + \lambda ^{m} ) + M\lambda ^{{m + 1}} = M(1 + \lambda + \lambda ^{2} + \cdots + \lambda ^{m} + \lambda ^{{m + 1}} ). \\ \end{aligned} $$

Under Case 2: $\forall i$:

$$\begin{aligned} |J^{m+1}(i)|= & {} |J^m(i)| \le M(1+\lambda +\lambda ^2+\cdots +\lambda ^m) \le M(1+\lambda +\lambda ^2+\cdots +\lambda ^m+\lambda ^{m+1}). \end{aligned}$$

$\square$

Data for discounted SMDP with known transition model $\textbf{P}_a$, $\textbf{R}_a$, and $\textbf{T}_a$, the transition probability, reward, and time matrices for action a, respectively, used are:

$$ {\mathbf{P}}_{1} = \left[ {\begin{array}{*{20}l} {0.6,0.2,0.1,0.1} \hfill \\ {0.2,0.3,0.1,0.4} \hfill \\ {0.1,0.6,0.2,0.1} \hfill \\ {0.2,0.4,0.2,0.2} \hfill \\ \end{array} } \right];~{\mathbf{P}}_{2} = \left[ {\begin{array}{*{20}l} {0.5,0.1,0.2,0.2} \hfill \\ {0.2,0.4,0,0.4} \hfill \\ {0.2,0.5,0.1,0.2} \hfill \\ {0.6,0.2,0.1,0.1} \hfill \\ \end{array} } \right]; {\text{ }}{\mathbf{R}}_{1} = \left[ {\begin{array}{*{20}l} {6, - 5,0,12} \hfill \\ {7,120,3,1} \hfill \\ {5,4,12,\,3} \hfill \\ {7,48,10,10} \hfill \\ \end{array} } \right];~{\mathbf{R}}_{2} = \left[ {\begin{array}{*{20}l} {100,17,0,10} \hfill \\ { - 14,13,0,1} \hfill \\ {9,40,7,12} \hfill \\ {12,12,10,14} \hfill \\ \end{array} } \right]; {\text{ }}{\mathbf{T}}_{1} = \left[ {\begin{array}{*{20}l} {1,5,2,5} \hfill \\ {120,60,30,30} \hfill \\ {2,7,12,14} \hfill \\ {135,55,45,30} \hfill \\ \end{array} } \right];~{\mathbf{T}}_{2} = \left[ {\begin{array}{*{20}l} {50,75,20,12} \hfill \\ {7,2,7,2} \hfill \\ {35,65,30,20} \hfill \\ {50,30,50,70} \hfill \\ \end{array} } \right] $$

Other constants in MAP were set to the following values: $\gamma =0.1$, $m_{\max }=k_{\max }=n_{\max }=100$, and $K=5$.

Data for post-earthquake disaster response SMDP $\textbf{PP}=[0.375,0.2,0.35,0.075]$, where PP(s) denotes the probability of going from the stable state to a primary state s. The following values were used for SP matrix:

$$ {\mathbf{SP}} = \left[ {\begin{array}{*{20}l} {0.1,0.3,0,0.5,0,0,0.1} \hfill \\ {0,0.1,0,0.35,0,0.35,0.2} \hfill \\ {0,0,0.1,0,0.5,0.2,0.2} \hfill \\ {0,0,\,0,0.1,0,0,0.9} \hfill \\ \end{array} } \right], $$

where $SP(s_1,s_2)$ denotes the probability of transition from a state, $s_1$, in the set of primary states to a state,$s_2$, in the union of the sets of the primary and secondary states. The following values were used for the hazard scores: $HS(1)=log(2)$, $HS(2)=log(4)$, $HS(3)=log(8)$, $HS(4)=log(6)$, $HS(5)=log(10)$, $HS(6)=log(12)$, and $HS(7)=log(14)$. Natural logarithms were used here for costs, as RL algorithms can suffer from computer overflow with large absolute values for costs, especially in discounted problems [24]. The response time for a state s is defined following [21]: $RT_c(s)=\sum _{d\in \mathcal {I}(s)}RT(d)\phi$, where $\mathcal {I}(s)$ denotes the incidents present in state s, RT(d) denotes the response time for an incident, d, and $\phi$ is a correction factor that equals 1 for a state that contains one incident, 1.2 for a state that contains two incidents, and 1.3 for a state that contains three incidents. The following values were used (all in hours): $RT(G)=7+\frac{5}{X}$, $RT(F)=21+\frac{15}{X}$, and $RT(BC)=35 +\frac{25}{X}$, where $X=1$ for the local agency and $X=2$ for the federal agency. Finally, for the local and federal agencies, the travel time were TRIA(0.5, 1, 1.5) and TRIA(4, 5, 6), respectively, in hrs, where TRIA denotes the triangular distribution. Other constants in MAP were set as follows: $\gamma =0.1$, $m_{\max }=k_{\max }=n_{\max }=10,000$, and $K=10$.

Steps in conservative approximate policy iteration (CAP) CAP uses Q-values instead of the future state-action values and bypasses model building. It has the same structure as MAP (see Algorithm 1 in main text) with the following exceptions: Step A (model-building) is skipped; Step B2 is replaced by evaluation of the Q-values, i.e., Eqn. (3) is replaced by: Update the Q-value for (i, a) as follows:

$$\begin{aligned} Q^{n+1}(i,a)\leftarrow (1-\beta ^n)Q^n(i,a)+\beta ^n\left[ r(i,a,j)+\exp \left( -\gamma t(i,a,j)\right) J^{\infty }(j)\right] ; \end{aligned}$$

Step C, policy is improved, i.e., for each $i\in \mathcal{S}$, select $\mu '(i)\in {argma }x_{a\in \mathcal{A}(i)}Q^{\infty }(i,a),$ where $Q^{\infty }(.,.)$ denotes the final Q-value obtained in Step B2.