Abstract
We study the safe reinforcement learning problem with nonlinear functionapproximation, where policy optimization is formulated as a constrainedoptimization problem with both the objective and the constraint being nonconvexfunctions. For such a problem, we construct a sequence of surrogate convexconstrained optimization problems by replacing the nonconvex functions locallywith convex quadratic functions obtained from policy gradient estimators. Weprove that the solutions to these surrogate problems converge to a stationarypoint of the original nonconvex problem. Furthermore, to extend our theoreticalresults, we apply our algorithm to examples of optimal control and multiagentreinforcement learning with safety constraints.
Quick Read (beta)
Convergent Policy Optimization for Safe Reinforcement Learning
Abstract
We study the safe reinforcement learning problem with nonlinear function approximation, where policy optimization is formulated as a constrained optimization problem with both the objective and the constraint being nonconvex functions. For such a problem, we construct a sequence of surrogate convex constrained optimization problems by replacing the nonconvex functions locally with convex quadratic functions obtained from policy gradient estimators. We prove that the solutions to these surrogate problems converge to a stationary point of the original nonconvex problem. Furthermore, to extend our theoretical results, we apply our algorithm to examples of optimal control and multiagent reinforcement learning with safety constraints.
Convergent Policy Optimization for Safe Reinforcement Learning
Ming Yu ^{†}^{†}thanks: The University of Chicago Booth School of Business, Chicago, IL. Email: [email protected]. Zhuoran Yang ^{†}^{†}thanks: Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ. Mladen Kolar ^{†}^{†}thanks: The University of Chicago Booth School of Business, Chicago, IL. Zhaoran Wang ^{†}^{†}thanks: Department of Industrial Engineering and Management Sciences, Northwestern University, Evanston, IL.
noticebox[b]33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada.\[email protected]
1 Introduction
Reinforcement learning [58] has achieved tremendous success in video games [41, 44, 56, 36, 66] and board games, such as chess and Go [53, 55, 54], in part due to powerful simulators [8, 62]. In contrast, due to physical limitations, realworld applications of reinforcement learning methods often need to take into consideration the safety of the agent [5, 26]. For instance, in expensive robotic and autonomous driving platforms, it is pivotal to avoid damages and collisions [25, 9]. In medical applications, we need to consider the switching cost [7].
A popular model of safe reinforcement learning is the constrained Markov decision process (CMDP), which generalizes the Markov decision process by allowing for inclusion of constraints that model the concept of safety [3]. In a CMDP, the cost is associated with each state and action experienced by the agent, and safety is ensured only if the expected cumulative cost is below a certain threshold. Intuitively, if the agent takes an unsafe action at some state, it will receive a huge cost that punishes risky attempts. Moreover, by considering the cumulative cost, the notion of safety is defined for the whole trajectory enabling us to examine the longterm safety of the agent, instead of focusing on individual stateaction pairs. For a CMDP, the goal is to take sequential decisions to achieve the expected cumulative reward under the safety constraint.
Solving a CMDP can be written as a linear program [3], with the number of variables being the same as the size of the state and action spaces. Therefore, such an approach is only feasible for the tabular setting, where we can enumerate all the stateaction pairs. For largescale reinforcement learning problems, where function approximation is applied, both the objective and constraint of the CMDP are nonconvex functions of the policy parameter. One common method for solving CMDP is to formulate an unconstrained saddlepoint optimization problem via Lagrangian multipliers and solve it using policy optimization algorithms [18, 60]. Such an approach suffers the following two drawbacks:
First, for each fixed Lagrangian multiplier, the inner minimization problem itself can be viewed as solving a new reinforcement learning problem. From the computational point of view, solving the saddlepoint optimization problem requires solving a sequence of MDPs with different reward functions. For a large scale problem, even solving a single MDP requires huge computational resources, making such an approach computationally infeasible.
Second, from a theoretical perspective, the performance of the saddlepoint approach hinges on solving the inner problem optimally. Existing theory only provides convergence to a stationary point where the gradient with respect to the policy parameter is zero [28, 37]. Moreover, the objective, as a bivariate function of the Lagrangian multiplier and the policy parameter, is not convexconcave and, therefore, firstorder iterative algorithms can be unstable [27].
In contrast, we tackle the nonconvex constrained optimization problem of the CMDP directly. We propose a novel policy optimization algorithm, inspired by [38]. Specifically, in each iteration, we replace both the objective and constraint by quadratic surrogate functions and update the policy parameter by solving the new constrained optimization problem. The two surrogate functions can be viewed as firstorder Taylorexpansions of the expected reward and cost functions where the gradients are estimated using policy gradient methods [59]. Additionally, they can be viewed as convex relaxations of the original nonconvex reward and cost functions. In §4 we show that, as the algorithm proceeds, we obtain a sequence of convex relaxations that gradually converge to a smooth function. More importantly, the sequence of policy parameters converges almost surely to a stationary point of the nonconvex constrained optimization problem.
Related work.
Our work is pertinent to the line of research on CMDP [3]. For CMDPs with large state and action spaces, [19] proposed an iterative algorithm based on a novel construction of Lyapunov functions. However, their theory only holds for the tabular setting. Using Lagrangian multipliers, [46, 18, 1, 60] proposed policy gradient [59], actorcritic [33], or trust region policy optimization [51] methods for CMDP or constrained risksensitive reinforcement learning [26]. These algorithms either do not have convergence guarantees or are shown to converge to saddlepoints of the Lagrangian using twotimescale stochastic approximations [10]. However, due to the projection on the Lagrangian multiplier, the saddlepoint achieved by these approaches might not be the stationary point of the original CMDP problem. In addition, [65] proposed a crossentropybased stochastic optimization algorithm, and proved the asymptotic behavior using ordinary differential equations. In contrast, our algorithm and the theoretical analysis focus on the discrete time CMDP. Outside of the CMDP setting, [31, 35] studied safe reinforcement learning with demonstration data, [61] studied the safe exploration problem with different safety constraints, and [4] studied multitask safe reinforcement learning.
Our contribution.
Our contribution is threefold. First, for the CMDP policy optimization problem where both the objective and constraint function are nonconvex, we propose to optimize a sequence of convex relaxation problems using convex quadratic functions. Solving these surrogate problems yields a sequence of policy parameters that converge almost surely to a stationary point of the original policy optimization problem. Second, to reduce the variance in the gradient estimator that is used to construct the surrogate functions, we propose an online actorcritic algorithm. Finally, as concrete applications, our algorithms are also applied to optimal control (in §5) and parallel and multiagent reinforcement learning problems with safety constraints (in supplementary material).
2 Background
A Markov decision process is denoted by $(\mathcal{S},\mathcal{A},P,\gamma ,r,\mu )$, where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space, $P$ is the transition probability distribution, $\gamma \in (0,1)$ is the discount factor, $r:\mathcal{S}\times \mathcal{A}\to \mathbb{R}$ is the reward function, and $\mu \in \mathcal{P}(\mathcal{S})$ is the distribution of the initial state ${s}_{0}\in \mathcal{S}$, where we denote $\mathcal{P}(\mathcal{X})$ as the set of probability distributions over $\mathcal{X}$ for any $\mathcal{X}$. A policy is a mapping $\pi :\mathcal{S}\to \mathcal{P}(\mathcal{A})$ that specifies the action that an agent will take when it is at state $s$.
Policy gradient method.
Let $\{{\pi}_{\theta}:\mathcal{S}\to \mathcal{P}(\mathcal{A})\}$ be a parametrized policy class, where $\theta \in \mathrm{\Theta}$ is the parameter defined on a compact set $\mathrm{\Theta}$. This parameterization transfers the original infinite dimensional policy class to a finite dimensional vector space and enables gradient based methods to be used to maximize (1). For example, the most popular Gaussian policy can be written as $\pi (\cdot s,\theta )=\mathcal{N}(\mu (s,\theta ),\sigma (s,\theta ))$, where the state dependent mean $\mu (s,\theta )$ and standard deviation $\sigma (s,\theta )$ can be further parameterized as $\mu (s,\theta )={\theta}_{\mu}^{\top}\cdot x(s)$ and $\sigma (s,\theta )=\mathrm{exp}\left({\theta}_{\sigma}^{\top}\cdot x(s)\right)$ with $x(s)$ being a state feature vector. The goal of an agent is to maximize the expected cumulative reward
$$R(\theta )={\mathbb{E}}_{\pi}\left[\sum _{t\ge 0}{\gamma}^{t}\cdot r({s}_{t},{a}_{t})\right],$$  (1) 
where ${s}_{0}\sim \mu $, and for all $t\ge 0$, we have ${s}_{t+1}\sim P(\cdot {s}_{t},{a}_{t})$ and ${a}_{t}\sim \pi (\cdot {s}_{t})$. Given a policy $\pi (\theta )$, we define the state and actionvalue functions of ${\pi}_{\theta}$, respectively, as
${V}^{\theta}(s)={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}r({s}_{t},{a}_{t})\right{s}_{0}=s],{Q}^{\theta}(s,a)={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}r({s}_{t},{a}_{t})\right{s}_{0}=s,{a}_{0}=a].$  (2) 
The policy gradient method updates the parameter $\theta $ through gradient ascent
${\theta}_{k+1}={\theta}_{k}+\eta \cdot {\widehat{\nabla}}_{\theta}R({\theta}_{k}),$  (3) 
where ${\widehat{\nabla}}_{\theta}R({\theta}_{k})$ is a stochastic estimate of the gradient ${\nabla}_{\theta}R({\theta}_{k})$ at $k$th iteration. Policy gradient method, as well as its variants (e.g. policy gradient with baseline [58], neural policy gradient [64, 39, 16]) is widely used in reinforcement learning. The gradient ${\nabla}_{\theta}R(\theta )$ can be estimated according to the policy gradient theorem [59],
${\nabla}_{\theta}R(\theta )=\mathbb{E}\left[{\nabla}_{\theta}\mathrm{log}{\pi}_{\theta}(s,a)\cdot {Q}^{\theta}(s,a)\right].$  (4) 
Actorcritic method.
To further reduce the variance of the policy gradient method, we could estimate both the policy parameter and value function simultaneously. This kind of method is called actorcritic algorithm [33], which is widely used in reinforcement learning. Specifically, in the value function evaluation (critic) step we estimate the actionvalue function ${Q}^{\theta}(s,a)$ using, for example, the temporal difference method TD(0) [20]. The policy parameter update (actor) step is implemented as before by the MonteCarlo method according to the policy gradient theorem (4) with the actionvalue ${Q}^{\theta}(s,a)$ replaced by the estimated value in the policy evaluation step.
Constrained MDP.
In this work, we consider an MDP problem with an additional constraint on the model parameter $\theta $. Specifically, when taking action at some state we incur some cost value. The constraint is such that the expected cumulative cost cannot exceed some predefined constant. A constrained Markov decision process (CMDP) is denoted by $(\mathcal{S},\mathcal{A},P,\gamma ,r,d,\mu )$, where $d:\mathcal{S}\times \mathcal{A}\to \mathbb{R}$ is the cost function and the other parameters are as before. The goal of an the agent in CMDP is to solve the following constrained problem
$\underset{\theta \in \mathrm{\Theta}}{\text{minimize}}J(\theta )={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot r({s}_{t},{a}_{t})\right],$  (5)  
$\text{subject to}D(\theta )={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot d({s}_{t},{a}_{t})\right]\le {D}_{0},$ 
where ${D}_{0}$ is a fixed constant. We consider only one constraint $D(\theta )\le {D}_{0}$, noting that it is straightforward to generalize to multiple constraints. Throughout this paper, we assume that both the reward and cost value functions are bounded: $\leftr({s}_{t},{a}_{t})\right\le {r}_{\mathrm{max}}$ and $\leftd({s}_{t},{a}_{t})\right\le {d}_{\mathrm{max}}$. Also, the parameter space $\mathrm{\Theta}$ is assumed to be compact.
3 Algorithm
In this section, we develop an algorithm to solve the optimization problem (5). Note that both the objective function and the constraint in (5) are nonconvex and involve expectation without closedform expression. As a constrained problem, a straightforward approach to solve (5) is to define the following Lagrangian function
$L(\theta ,\lambda )=J(\theta )+\lambda \cdot \left[D(\theta ){D}_{0}\right],$  (6) 
and solve the dual problem
$\underset{\lambda \ge 0}{inf}\underset{\theta}{sup}L(\theta ,\lambda ).$  (7) 
However, this problem is a nonconvex minimax problem and, therefore, is hard to solve and establish theoretical guarantees for solutions [2]. Another approach to solve (5) is to replace $J(\theta )$ and $D(\theta )$ by surrogate functions with nice properties. For example, one can iteratively construct local quadratic approximations that are strongly convex [52], or are an upper bound for the original function [57]. However, an immediate problem of this naive approach is that, even if the original problem (5) is feasible, the convex relaxation problem need not be. Also, these methods only deal with deterministic and/or convex constraints.
In this work, we propose an iterative algorithm that approximately solves (5) by constructing a sequence of convex relaxations, inspired by [38]. Our method is able to handle the possible infeasible situation due to the convex relaxation as mentioned above, and handle stochastic and nonconvex constraint. Since we do not have access to $J(\theta )$ or $D(\theta )$, we first define the sample negative cumulative reward and cost functions as
${J}^{*}(\theta )={\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot r({s}_{t},{a}_{t})\mathit{\hspace{1em}\hspace{1em}}\text{and}\mathit{\hspace{1em}\hspace{1em}}{D}^{*}(\theta )={\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot d({s}_{t},{a}_{t}).$  (8) 
Given $\theta $, ${J}^{*}(\theta )$ and ${D}^{*}(\theta )$ are the sample negative cumulative reward and cost value of a realization (i.e., a trajectory) following policy ${\pi}_{\theta}$. Note that both ${J}^{*}(\theta )$ and ${D}^{*}(\theta )$ are stochastic due to the randomness in the policy, state transition distribution, etc. With some abuse of notation, we use ${J}^{*}(\theta )$ and ${D}^{*}(\theta )$ to denote both a function of $\theta $ and a value obtained by the realization of a trajectory. Clearly we have $J(\theta )=\mathbb{E}\left[{J}^{*}(\theta )\right]$ and $D(\theta )=\mathbb{E}\left[{D}^{*}(\theta )\right]$.
We start from some (possibly infeasible) ${\theta}_{0}$. Let ${\theta}_{k}$ denote the estimate of the policy parameter in the $k$th iteration. As mentioned above, we do not have access to the expected cumulative reward $J(\theta )$. Instead we sample a trajectory following the current policy ${\pi}_{{\theta}_{k}}$ and obtain a realization of the negative cumulative reward value and the gradient of it as ${J}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$, respectively. The cumulative reward value is obtained by MonteCarlo estimation, and the gradient is also obtained by MonteCarlo estimation according to the policy gradient theorem in (4). We provide more details on the realization step later in this section. Similarly, we use the same procedure for the cost function and obtain realizations ${D}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{D}^{*}({\theta}_{k})$.
We approximate $J(\theta )$ and $D(\theta )$ at ${\theta}_{k}$ by the quadratic surrogate functions
$\stackrel{~}{J}(\theta ,{\theta}_{k},\tau )$  $={J}^{*}({\theta}_{k})+\u27e8{\nabla}_{\theta}{J}^{*}({\theta}_{k}),\theta {\theta}_{k}\u27e9+\tau {\parallel \theta {\theta}_{k}\parallel}_{2}^{2},$  (9)  
$\stackrel{~}{D}(\theta ,{\theta}_{k},\tau )$  $={D}^{*}({\theta}_{k})+\u27e8{\nabla}_{\theta}{D}^{*}({\theta}_{k}),\theta {\theta}_{k}\u27e9+\tau {\parallel \theta {\theta}_{k}\parallel}_{2}^{2},$  (10) 
where $\tau >0$ is any fixed constant. In each iteration, we solve the optimization problem
${\overline{\theta}}_{k}=\underset{\theta}{argmin}{\overline{J}}^{(k)}(\theta )\mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}\hspace{1em}}{\overline{D}}^{(k)}(\theta )\le {D}_{0},$  (11) 
where we define
${\overline{J}}^{(k)}(\theta )=(1{\rho}_{k})\cdot {\overline{J}}^{(k1)}(\theta )+{\rho}_{k}\cdot \stackrel{~}{J}(\theta ,{\theta}_{k},\tau ),$  (12)  
${\overline{D}}^{(k)}(\theta )=(1{\rho}_{k})\cdot {\overline{D}}^{(k1)}(\theta )+{\rho}_{k}\cdot \stackrel{~}{D}(\theta ,{\theta}_{k},\tau ),$  (13) 
with the initial value ${\overline{J}}^{(0)}(\theta )={\overline{D}}^{(0)}(\theta )=0$. Here ${\rho}_{k}$ is the weight parameter to be specified later. According to the definition (9) and (10), problem (11) is a convex quadratically constrained quadratic program (QCQP). Therefore, it can be efficiently solved by, for example, the interior point method. However, as mentioned before, even if the original problem (5) is feasible, the convex relaxation problem (11) could be infeasible. In this case, we instead solve the following feasibility problem
${\overline{\theta}}_{k}=\underset{\theta ,\alpha}{argmin}\alpha \mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}\hspace{1em}}{\overline{D}}^{(k)}(\theta )\le {D}_{0}+\alpha .$  (14) 
In particular, we relax the infeasible constraint and find ${\overline{\theta}}_{k}$ as the solution that gives the minimum relaxation. Due to the specific form in (10), ${\overline{D}}^{(k)}(\theta )$ is decomposable into quadratic forms of each component of $\theta $, with no terms involving ${\theta}_{i}\cdot {\theta}_{j}$. Therefore, the solution to problem (14) can be written in a closed form. Given ${\overline{\theta}}_{k}$ from either (11) or (14), we update ${\theta}_{k}$ by
${\theta}_{k+1}=(1{\eta}_{k})\cdot {\theta}_{k}+{\eta}_{k}\cdot {\overline{\theta}}_{k},$  (15) 
where ${\eta}_{k}$ is the learning rate to be specified later. Note that although we consider only one constraint in the algorithm, both the algorithm and the theoretical result in Section 4 can be directly generalized to multiple constraints setting. The whole procedure is summarized in Algorithm 3.
[tb] \[email protected]@algorithmic[1] \STATEInput: Initial value ${\theta}_{0}$, $\tau $, $\{{\rho}_{k}\},\{{\eta}_{k}\}$. \FOR$k=1,2,3,\mathrm{\dots}$ \STATEObtain a sample ${J}^{*}({\theta}_{k})$ and ${D}^{*}({\theta}_{k})$ by MonteCarlo sampling. \STATEObtain a sample ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{D}^{*}({\theta}_{k})$ by policy gradient theorem. \IFproblem (11) is feasible \STATEObtain ${\overline{\theta}}_{k}$ by solving (11). \ELSE\STATEObtain ${\overline{\theta}}_{k}$ by solving (14). \ENDIF\STATEUpdate ${\theta}_{k+1}$ by (15). \ENDFOR
Obtaining realizations ${J}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$.
We detail how to obtain realizations ${J}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$ corresponding to the lines 3 and 4 in Algorithm 3. The realizations of ${D}^{*}({\theta}_{k})$ and ${\nabla}_{\theta}{D}^{*}({\theta}_{k})$ can be obtained similarly.
First, we discuss finite horizon setting, where we can sample the full trajectory according to the policy ${\pi}_{\theta}$. In particular, for any ${\theta}_{k}$, we use the policy ${\pi}_{{\theta}_{k}}$ to sample a trajectory and obtain ${J}^{*}({\theta}_{k})$ by MonteCarlo method. The gradient ${\nabla}_{\theta}J(\theta )$ can be estimated by the policy gradient theorem [59],
${\nabla}_{\theta}J(\theta )={\mathbb{E}}_{{\pi}_{\theta}}\left[{\nabla}_{\theta}\mathrm{log}{\pi}_{\theta}(s,a)\cdot {Q}^{\theta}(s,a)\right].$  (16) 
Again we can sample a trajectory and obtain the policy gradient realization ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$ by MonteCarlo method.
In infinite horizon setting, we cannot sample the infinite length trajectory. In this case, we utilize the truncation method introduced in [48], which truncates the trajectory at some stage $T$ and scales the undiscounted cumulative reward to obtain an unbiased estimation. Intuitively, if the discount factor $\gamma $ is close to $0$, then the future reward would be discounted heavily and, therefore, we can obtain an accurate estimate with a relatively small number of stages. On the other hand, if $\gamma $ is close to $1$, then the future reward is more important compared to the small $\gamma $ case and we have to sample a long trajectory. Taking this intuition into consideration, we define $T$ to be a geometric random variable with parameter $1\gamma $: $\mathrm{Pr}(T=t)=(1\gamma ){\gamma}^{t}$. Then, we simulate the trajectory until stage $T$ and use the estimator ${J}_{\text{truncate}}(\theta )=(1\gamma )\cdot {\sum}_{t=0}^{T}r({s}_{t},{a}_{t})$, which is an unbiased estimator of the expected negative cumulative reward $J(\theta )$, as proved in proposition 5 in [43]. We can apply the same truncation procedure to estimate the policy gradient ${\nabla}_{\theta}J(\theta )$.
Variance reduction.
Using the naive sampling method described above, we may suffer from high variance problem. To reduce the variance, we can modify the above procedure in the following ways. First, instead of sampling only one trajectory in each iteration, a more practical and stable way is to sample several trajectories and take average to obtain the realizations. As another approach, we can subtract a baseline function from the actionvalue function ${Q}^{\theta}(s,a)$ in the policy gradient estimation step (16) to reduce the variance without changing the expectation. A popular choice of the baseline function is the statevalue function ${V}^{\theta}(s)$ as defined in (2). In this way, we can replace ${Q}^{\theta}(s,a)$ in (16) by the advantage function ${A}^{\theta}(s,a)$ defined as
$${A}^{\theta}(s,a)={Q}^{\theta}(s,a){V}^{\theta}(s).$$ 
This modification corresponds to the standard REINFORCE with Baseline algorithm [58] and can significantly reduce the variance of policy gradient.
Actorcritic method.
Finally, we can use an actorcritic update to improve the performance further. In this case, since we need unbiased estimators for both the gradient and the reward value in (9) and (10) in online fashion, we modify our original problem (5) to average reward setting as
$\underset{\theta \in \mathrm{\Theta}}{\text{minimize}}J(\theta )=\underset{T\to \mathrm{\infty}}{lim}{\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \frac{1}{T}}{\displaystyle \sum _{t=0}^{T}}r({s}_{t},{a}_{t})\right],$  (17)  
$\text{subject to}D(\theta )=\underset{T\to \mathrm{\infty}}{lim}{\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \frac{1}{T}}{\displaystyle \sum _{t=0}^{T}}d({s}_{t},{a}_{t})\right]\le {D}_{0}.$ 
Let ${V}_{\theta}^{J}(s)$ and ${V}_{\theta}^{D}(s)$ denote the value and cost functions corresponding to (2). We use possibly nonlinear approximation with parameter $w$ for the value function: ${V}_{w}^{J}(s)$ and $v$ for the cost function: ${V}_{v}^{D}(s)$. In the critic step, we update $w$ and $v$ by TD(0) with step size ${\beta}_{w}$ and ${\beta}_{v}$; in the actor step, we solve our proposed convex relaxation problem to update $\theta $. The actorcritic procedure is summarized in Algorithm 3. Here $J$ and $D$ are estimators of $J({\theta}_{k})$ and $D({\theta}_{k})$. Both of $J$ and $D$, and the TD error ${\delta}^{J}$, ${\delta}^{D}$ can be initialized as 0.
The usage of the actorcritic method helps reduce variance by using a value function instead of MonteCarlo sampling. Specifically, in Algorithm 3 we need to obtain a sample trajectory and calculate ${J}^{*}(\theta )$ and ${\nabla}_{\theta}{J}^{*}(\theta )$ by MonteCarlo sampling. This step has a high variance since we need to sample a potentially long trajectory and sum up a lot of random rewards. In contrast, in Algorithm 3, this step is replaced by a value function ${V}_{w}^{J}(s)$, which reduces the variance.
[tb]
\[email protected]@algorithmic[1]
\FOR$k=1,2,3,\mathrm{\dots}$
\STATETake action $a$, observe reward $r$, cost $d$, and new state ${s}^{\prime}$.
\STATECritic step:
\STATE $w\leftarrow w+{\beta}_{w}\cdot {\delta}^{J}{\nabla}_{w}{V}_{w}^{J}(s),J\leftarrow J+{\beta}_{w}\cdot \left(rJ\right)$.
\STATE $v\leftarrow v+{\beta}_{v}\cdot {\delta}^{D}{\nabla}_{v}{V}_{v}^{J}(s),D\leftarrow D+{\beta}_{v}\cdot \left(dD\right)$.
\STATECalculate TD error:
\STATE ${\delta}^{J}=rJ+{V}_{w}^{J}({s}^{\prime}){V}_{w}^{J}(s)$.
\STATE ${\delta}^{D}=dD+{V}_{v}^{D}({s}^{\prime}){V}_{v}^{D}(s)$.
\STATEActor step:
\STATE Solve ${\overline{\theta}}_{k}$ by (11) or (14) with
${J}^{*}({\theta}_{k})$, ${\nabla}_{\theta}{J}^{*}({\theta}_{k})$ in (9) replaced by $J$ and ${\delta}^{J}\cdot {\nabla}_{\theta}\mathrm{log}{\pi}_{\theta}(s,a)$;
${D}^{*}({\theta}_{k})$, ${\nabla}_{\theta}{D}^{*}({\theta}_{k})$ in (10) replaced by $D$ and ${\delta}^{D}\cdot {\nabla}_{\theta}\mathrm{log}{\pi}_{\theta}(s,a)$.
\STATE$s\leftarrow {s}^{\prime}$.
\ENDFOR
4 Theoretical Result
In this section, we show almost sure convergence of the iterates obtained by our algorithm to a stationary point. We start by stating some mild assumptions on the original problem (5) and the choice of some parameters in Algorithm 3.
Assumption 1
The choice of $\mathrm{\{}{\eta}_{k}\mathrm{\}}$ and $\mathrm{\{}{\rho}_{k}\mathrm{\}}$ satisfy ${\mathrm{lim}}_{k\mathrm{\to}\mathrm{\infty}}\mathit{}{\mathrm{\sum}}_{k}{\eta}_{k}\mathrm{=}\mathrm{\infty}$, ${\mathrm{lim}}_{k\mathrm{\to}\mathrm{\infty}}\mathit{}{\mathrm{\sum}}_{k}{\rho}_{k}\mathrm{=}\mathrm{\infty}$ and $$. Furthermore, we have ${\mathrm{lim}}_{k\mathrm{\to}\mathrm{\infty}}\mathit{}{\eta}_{k}\mathrm{/}{\rho}_{k}\mathrm{=}\mathrm{0}$ and ${\eta}_{k}$ is decreasing.
Assumption 2
For any realization, ${J}^{\mathrm{*}}\mathit{}\mathrm{(}\theta \mathrm{)}$ and ${D}^{\mathrm{*}}\mathit{}\mathrm{(}\theta \mathrm{)}$ are continuously differentiable as functions of $\theta $. Moreover, ${J}^{\mathrm{*}}\mathit{}\mathrm{(}\theta \mathrm{)}$, ${D}^{\mathrm{*}}\mathit{}\mathrm{(}\theta \mathrm{)}$, and their derivatives are uniformly Lipschitz continuous.
Assumption 1 allows us to specify the learning rates. A practical choice would be ${\eta}_{k}={k}^{{c}_{1}}$ and ${\rho}_{k}={k}^{{c}_{2}}$ with $$. This assumption is standard for gradientbased algorithms. Assumption 2 is also standard and is known to hold for a number of models. It ensures that the reward and cost functions are sufficiently regular. In fact, it can be relaxed such that each realization is Lipschitz (not uniformly), and the event that we keep generating realizations with monotonically increasing Lipschitz constant is an event with probability 0. See condition iv) in [67] and the discussion thereafter. Also, see [45] for sufficient conditions such that both the expected cumulative reward function and the gradient of it are Lipschitz.
The following Assumption 3 is useful only when we initialize with an infeasible point. We first state it here and we will discuss this assumption after the statement of the main theorem.
Assumption 3
Suppose $\mathrm{(}{\theta}_{S}\mathrm{,}{\alpha}_{S}\mathrm{)}$ is a stationary point of the optimization problem
$\underset{\theta ,\alpha}{\text{\mathit{m}\mathit{i}\mathit{n}\mathit{i}\mathit{m}\mathit{i}\mathit{z}\mathit{e}}}\alpha \mathit{\hspace{1em}\hspace{1em}}\mathit{\text{subject to}}\mathit{\hspace{1em}\hspace{1em}}D(\theta )\le {D}_{0}+\alpha .$  (18) 
We have that ${\theta}_{S}$ is a feasible point of the original problem (5), i.e. $D\mathit{}\mathrm{(}{\theta}_{S}\mathrm{)}\mathrm{\le}{D}_{\mathrm{0}}$.
We are now ready to state the main theorem.
Theorem 4
Suppose the Assumptions 1 and 2 are satisfied with small enough initial step size ${\eta}_{\mathrm{0}}$. Suppose also that, either ${\theta}_{\mathrm{0}}$ is a feasible point, or Assumption 3 is satisfied. If there is a subsequence $\mathrm{\{}{\theta}_{{k}_{j}}\mathrm{\}}$ of $\mathrm{\{}{\theta}_{k}\mathrm{\}}$ that converges to some $\stackrel{\mathrm{~}}{\theta}$, then there exist uniformly continuous functions $\widehat{J}\mathit{}\mathrm{(}\theta \mathrm{)}$ and $\widehat{D}\mathit{}\mathrm{(}\theta \mathrm{)}$ satisfying
$$\underset{j\to \mathrm{\infty}}{lim}{\overline{J}}^{({k}_{j})}(\theta )=\widehat{J}(\theta )\mathit{\hspace{1em}\hspace{1em}}\text{\mathit{a}\mathit{n}\mathit{d}}\mathit{\hspace{1em}\hspace{1em}}\underset{j\to \mathrm{\infty}}{lim}{\overline{D}}^{({k}_{j})}(\theta )=\widehat{D}(\theta ).$$ 
Furthermore, suppose there exists $\theta $ such that $$ (i.e. the Slater’s condition holds), then $\stackrel{\mathrm{~}}{\theta}$ is a stationary point of the original problem (5) almost surely.
The proof of Theorem 4 is provided in the supplementary material.
Note that Assumption 3 is not necessary if we start from a feasible point, or we reach a feasible point in the iterates, which could be viewed as an initializer. Assumption 3 makes sure that the iterates in Algorithm 3 keep making progress without getting stuck at any infeasible stationary point. A similar condition is assumed in [38] for an infeasible initializer. If it turns out that ${\theta}_{0}$ is infeasible and Assumption 3 is violated, then the convergent point may be an infeasible stationary point of (18). In practice, if we can find a feasible point of the original problem, then we proceed with that point. Alternatively, we could generate multiple initializers and obtain iterates for all of them. As long as there is a feasible point in one of the iterates, we can view this feasible point as the initializer and Theorem 4 follows without Assumption 3. In our later experiments, for every single replicate, we could reach a feasible point, and therefore Assumption 3 is not necessary.
Our algorithm does not guarantee safe exploration during the training phase. Ensuring safety during learning is a more challenging problem. Sometimes even finding a feasible point is not straightforward, otherwise Assumption 3 is not necessary.
Our proposed algorithm is inspired by [38]. Compared to [38] which deals with an optimization problem, solving the safe reinforcement learning problem is more challenging. We need to verify that the Lipschitz condition is satisfied, and also the policy gradient has to be estimated (instead of directly evaluated as in a standard optimization problem). The usage of the ActorCritic algorithm reduces the variance of the sampling, which is unique to Reinforcement learning.
5 Application to Constrained LinearQuadratic Regulator
We apply our algorithm to the linearquadratic regulator (LQR), which is one of the most fundamental problems in control theory. In the LQR setting, the state dynamic equation is linear, the cost function is quadratic, and the optimal control theory tells us that the optimal control for LQR is a linear function of the state [23, 6]. LQR can be viewed as an MDP problem and it has attracted a lot of attention in the reinforcement learning literature [12, 13, 21, 47].
We consider the infinitehorizon, discretetime LQR problem. Denote ${x}_{t}$ as the state variable and ${u}_{t}$ as the control variable. The state transition and the control sequence are given by
${x}_{t+1}$  $=A{x}_{t}+B{u}_{t}+{v}_{t},$  (19)  
${u}_{t}$  $=F{x}_{t}+{w}_{t},$ 
where ${v}_{t}$ and ${w}_{t}$ represent possible Gaussian white noise, and the initial state is given by ${x}_{0}$. The goal is to find the control parameter matrix $F$ such that the expected total cost is minimized. The usual cost function of LQR corresponds to the negative reward in our setting and we impose an additional quadratic constraint on the system. The overall optimization problem is given by
$\text{minimize}J(F)=\mathbb{E}\left[{\displaystyle \sum _{t\ge 0}}{x}_{t}^{\top}{Q}_{1}{x}_{t}+{u}_{t}^{\top}{R}_{1}{u}_{t}\right],$  (20)  
$\text{subject to}D(F)=\mathbb{E}\left[{\displaystyle \sum _{t\ge 0}}{x}_{t}^{\top}{Q}_{2}{x}_{t}+{u}_{t}^{\top}{R}_{2}{u}_{t}\right]\le {D}_{0},$ 
where ${Q}_{1},{Q}_{2},{R}_{1},$ and ${R}_{2}$ are positive definite matrices. Note that even thought the matrices are positive definite, both the objective function $J$ and the constraint $D$ are nonconvex with respect to the parameter $F$. Furthermore, with the additional constraint, the optimal control sequence may no longer be linear in the state ${x}_{t}$. Nevertheless, in this work, we still consider linear control given by (19) and the goal is to find the best linear control for this constrained LQR problem. We assume that the choice of $A,B$ are such that the optimal cost is finite.
Random initial state.
We first consider the setting where the initial state ${x}_{0}\sim \mathcal{D}$ follows a random distribution $\mathcal{D}$, while both the state transition and the control sequence (19) are deterministic (i.e. ${v}_{t}={w}_{t}=0$). In this random initial state setting, [24] showed that without the constraint, the policy gradient method converges efficiently to the global optima in polynomial time. In the constrained case, we can explicitly write down the objective and constraint function, since the only randomness comes from the initial state. Therefore, we have the state dynamic ${x}_{t+1}=(ABF){x}_{t}$ and the objective function has the following expression ([24], Lemma 1)
$$J(F)={\mathbb{E}}_{{x}_{0}\sim \mathcal{D}}\left[{x}_{0}^{\top}{P}_{F}{x}_{0}\right],$$  (21) 
where ${P}_{F}$ is the solution to the following equation
$${P}_{F}={Q}_{1}+{F}^{\top}{R}_{1}F+{(ABF)}^{\top}{P}_{F}(ABF).$$  (22) 
The gradient is given by
${\nabla}_{F}J(F)=2\left(\left({R}_{1}+{B}^{\top}{P}_{F}B\right)F{B}^{\top}{P}_{F}A\right)\cdot \left[{\mathbb{E}}_{{x}_{0}\sim \mathcal{D}}{\displaystyle \sum _{t=0}^{\mathrm{\infty}}}{x}_{t}{x}_{t}^{\top}\right].$  (23) 
Let ${S}_{F}={\sum}_{t=0}^{\mathrm{\infty}}{x}_{t}{x}_{t}^{\top}$ and observe that
$${S}_{F}={x}_{0}{x}_{0}^{\top}+(ABF){S}_{F}{(ABF)}^{\top}.$$  (24) 
We start from some ${F}_{0}$ and apply our Algorithm 3 to solve the constrained LQR problem. In iteration $k$, with the current estimator denoted by ${F}_{k}$, we first obtain an estimator of ${P}_{{F}_{k}}$ by starting from ${Q}_{1}$ and iteratively applying the recursion ${P}_{{F}_{k}}\leftarrow {Q}_{1}+{F}_{k}^{\top}{R}_{1}{F}_{k}+{(AB{F}_{k})}^{\top}{P}_{{F}_{k}}(AB{F}_{k})$ until convergence. Next, we sample an ${x}_{0}^{*}$ from the distribution $\mathcal{D}$ and follow a similar recursion given by (24) to obtain an estimate of ${S}_{{F}_{k}}$. Plugging the sample ${x}_{0}^{*}$ and the estimates of ${P}_{{F}_{k}}$ and ${S}_{{F}_{k}}$ into (21) and (23), we obtain the sample reward value ${J}^{*}({F}_{k})$ and ${\nabla}_{F}{J}^{*}({F}_{k})$, respectively. With these two values, we follow (9) and (12) and obtain ${\overline{J}}^{(k)}(F)$. We apply the same procedure to the cost function $D(F)$ with ${Q}_{1},{R}_{1}$ replaced by ${Q}_{2},{R}_{2}$ to obtain ${\overline{D}}^{(k)}(F)$. Finally we solve the optimization problem (11) (or (14) if (11) is infeasible) and obtain ${F}_{k+1}$ by (15).
Random state transition and control.
We then consider the setting where both ${v}_{t}$ and ${w}_{t}$ are independent standard Gaussian white noise. In this case, the state dynamic can be written as ${x}_{t+1}=(ABF){x}_{t}+{\u03f5}_{t}$ where ${\u03f5}_{t}\sim \mathcal{N}(0,I+B{B}^{\top})$. Let ${P}_{F}$ be defined as in (22) and ${S}_{F}$ be the solution to the following Lyapunov equation
$${S}_{F}=I+B{B}^{\top}+(ABF){S}_{F}{(ABF)}^{\top}.$$  (25) 
The objective function has the following expression ([68], Proposition 3.1)
$J(F)={\mathbb{E}}_{x\sim \mathcal{N}(0,{S}_{F})}\left[{x}^{\top}({Q}_{1}+{F}^{\top}{R}_{1}F)x\right]+tr({R}_{1}),$  (26) 
and the gradient is given by
$${\nabla}_{F}J(F)=2\left(\left({R}_{1}+{B}^{\top}{P}_{F}B\right)F{B}^{\top}{P}_{F}A\right)\cdot {\mathbb{E}}_{x\sim \mathcal{N}(0,{S}_{F})}\left[x{x}^{\top}\right].$$  (27) 
Although in this setting it is straightforward to calculate the expectation in a closed form, we keep the current expectation form to be in line with our algorithm. Moreover, when the error distribution is more complicated or unknown, we can no longer calculate the closed form expression and have to sample in each iteration. With the formulas given by (26) and (27), we again apply our Algorithm 3. We sample $x\sim \mathcal{N}(0,{S}_{F})$ in each iteration and solve the optimization problem (11) or (14). The whole procedure is similar to the random initial state case described above.
Other applications.
Our algorithm can also be applied to constrained parallel MDP and constrained multiagent MDP problem. Due to the space limit, we relegate them to supplementary material.
6 Experiment
We verify the effectiveness of the proposed algorithm through experiments. We focus on the LQR setting with a random initial state as discussed in Section 5. In this experiment we set $x\in {\mathbb{R}}^{15}$ and $u\in {\mathbb{R}}^{8}$. The initial state distribution is uniform on the unit cube: ${x}_{0}\sim \mathcal{D}=\text{Uniform}\left({[1,1]}^{15}\right)$. Each element of $A$ and $B$ is sampled independently from the standard normal distribution and scaled such that the eigenvalues of $A$ are within the range $(1,1)$. We initialize ${F}_{0}$ as an allzero matrix, and the choice of the constraint function and the value ${D}_{0}$ are such that (1) the constrained problem is feasible; (2) the solution of the unconstrained problem does not satisfy the constraint, i.e., the problem is not trivial; (3) the initial value ${F}_{0}$ is not feasible. The learning rates are set as ${\eta}_{k}=\frac{2}{3}{k}^{3/4}$ and ${\rho}_{k}=\frac{2}{3}{k}^{2/3}$. The conservative choice of step size is to avoid the situation where an eigenvalue of $ABF$ runs out of the range $(1,1)$, and so the system is stable. ^{1}^{1} 1 The code is available at https://github.com/ming93/Safe_reinforcement_learning
Figure 1(a) and 1(b) show the constraint and objective value in each iteration, respectively. The red horizontal line in Figure 1(a) is for ${D}_{0}$, while the horizontal line in Figure 1(b) is for the unconstrained minimum objective value. We can see from Figure 1(a) that we start from an infeasible point, and the problem becomes feasible after about 100 iterations. The objective value is in general decreasing after becoming feasible, but never lower than the unconstrained minimum, as shown in Figure 1(b).
Comparison with the Lagrangian method.
We compare our proposed method with the usual Lagrangian method. For the Lagrangian method, we follow the algorithm proposed in [18] for safe reinforcement learning, which iteratively applies gradient descent on the parameter $F$ and gradient ascent on the Lagrangian multiplier $\lambda $ for the Lagrangian function until convergence.
Table 1 reports the comparison results with mean and standard deviation based on 50 replicates. In the second and third columns, we compare the minimum objective value and the number of iterations to achieve it. We also consider an approximate version, where we are satisfied with the result if the objective value exceeds less than 0.02% of the minimum value. The fourth and fifth columns show the comparison results for this approximate version. We can see that both methods achieve similar minimum objective values, but ours requires less number of policy updates, for both minimum and approximate minimum version.
min value 
# iterations 
approx. min value  approx. # iterations 

Our method  30.689 $\pm $ 0.114  2001 $\pm $ 1172  30.694 $\pm $ 0.114  604.3 $\pm $ 722.4 
Lagrangian  30.693 $\pm $ 0.113  7492 $\pm $ 1780  30.699 $\pm $ 0.113  5464 $\pm $ 2116 
References
 [1] Joshua Achiam, David Held, Aviv Tamar, and Pieter Abbeel. Constrained policy optimization. In International Conference on Machine Learning, pages 22–31, 2017.
 [2] Leonard Adolphs. Non convexconcave saddle point optimization. Master’s thesis, ETH Zurich, 2018.
 [3] Eitan Altman. Constrained Markov decision processes, volume 7. CRC Press, 1999.
 [4] Haitham Bou Ammar, Rasul Tutunov, and Eric Eaton. Safe policy search for lifelong reinforcement learning with sublinear regret. In International Conference on Machine Learning, pages 2361–2369, 2015.
 [5] Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané. Concrete problems in ai safety. arXiv preprint arXiv:1606.06565, 2016.
 [6] Brian DO Anderson and John B Moore. Optimal control: linear quadratic methods. Courier Corporation, 2007.
 [7] Yu Bai, Tengyang Xie, Nan Jiang, and YuXiang Wang. Provably efficient qlearning with low switching cost. arXiv preprint arXiv:1905.12849, 2019.
 [8] Marc G Bellemare, Yavar Naddaf, Joel Veness, and Michael Bowling. The arcade learning environment: An evaluation platform for general agents. Journal of Artificial Intelligence Research, 47:253–279, 2013.
 [9] Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe modelbased reinforcement learning with stability guarantees. In Advances in Neural Information Processing Systems, pages 908–918, 2017.
 [10] Vivek S Borkar. Stochastic approximation with two time scales. Systems & Control Letters, 29(5):291–294, 1997.
 [11] Craig Boutilier. Planning, learning and coordination in multiagent decision processes. In Proceedings of the 6th conference on Theoretical aspects of rationality and knowledge, pages 195–210. Morgan Kaufmann Publishers Inc., 1996.
 [12] Steven J Bradtke. Reinforcement learning applied to linear quadratic regulation. In Advances in neural information processing systems, pages 295–302, 1993.
 [13] Steven J Bradtke, B Erik Ydstie, and Andrew G Barto. Adaptive linear quadratic control using policy iteration. In Proceedings of the American control conference, volume 3, pages 3475–3475. Citeseer, 1994.
 [14] Leo Breiman. Random forests. Machine learning, 45(1):5–32, 2001.
 [15] Lucian Busoniu, Robert Babuska, and Bart De Schutter. A comprehensive survey of multiagent reinforcement learning. IEEE Transactions on Systems, Man, And CyberneticsPart C: Applications and Reviews, 38 (2), 2008, 2008.
 [16] Qi Cai, Zhuoran Yang, Jason D Lee, and Zhaoran Wang. Neural temporaldifference learning converges to global optima. arXiv preprint arXiv:1905.10027, 2019.
 [17] Tianyi Chen, Kaiqing Zhang, Georgios B Giannakis, and Tamer Başar. Communicationefficient distributed reinforcement learning. arXiv preprint arXiv:1812.03239, 2018.
 [18] Yinlam Chow, Mohammad Ghavamzadeh, Lucas Janson, and Marco Pavone. Riskconstrained reinforcement learning with percentile risk criteria. Journal of Machine Learning Research, 18(167):1–167, 2017.
 [19] Yinlam Chow, Ofir Nachum, Edgar DuenezGuzman, and Mohammad Ghavamzadeh. A lyapunovbased approach to safe reinforcement learning. arXiv preprint arXiv:1805.07708, 2018.
 [20] Christoph Dann, Gerhard Neumann, and Jan Peters. Policy evaluation with temporal differences: A survey and comparison. The Journal of Machine Learning Research, 15(1):809–883, 2014.
 [21] Sarah Dean, Horia Mania, Nikolai Matni, Benjamin Recht, and Stephen Tu. On the sample complexity of the linear quadratic regulator. arXiv preprint arXiv:1710.01688, 2017.
 [22] Nelson Dunford and Jacob T Schwartz. Linear operators part I: general theory, volume 7. Interscience publishers New York, 1958.
 [23] Lawrence C Evans. An introduction to mathematical optimal control theory. Lecture Notes, University of California, Department of Mathematics, Berkeley, 2005.
 [24] Maryam Fazel, Rong Ge, Sham Kakade, and Mehran Mesbahi. Global convergence of policy gradient methods for the linear quadratic regulator. In International Conference on Machine Learning, pages 1466–1475, 2018.
 [25] Jaime F Fisac, Anayo K Akametalu, Melanie N Zeilinger, Shahab Kaynama, Jeremy Gillula, and Claire J Tomlin. A general safety framework for learningbased control in uncertain robotic systems. IEEE Transactions on Automatic Control, 2018.
 [26] Javier Garcıa and Fernando Fernández. A comprehensive survey on safe reinforcement learning. Journal of Machine Learning Research, 16(1):1437–1480, 2015.
 [27] Ian Goodfellow, Jean PougetAbadie, Mehdi Mirza, Bing Xu, David WardeFarley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672–2680, 2014.
 [28] Ivo Grondman, Lucian Busoniu, Gabriel AD Lopes, and Robert Babuska. A survey of actorcritic reinforcement learning: Standard and natural policy gradients. IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews), 42(6):1291–1307, 2012.
 [29] Jingyu He, Saar Yalov, and P Richard Hahn. XBART: Accelerated Bayesian additive regression trees. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1130–1138, 2019.
 [30] Jingyu He, Saar Yalov, Jared Murray, and P Richard Hahn. Stochastic tree ensembles for regularized supervised learning. Technical report, 2019.
 [31] Jessie Huang, Fa Wu, Doina Precup, and Yang Cai. Learning safe policies with expert guidance. arXiv preprint arXiv:1805.08313, 2018.
 [32] John L Kelley. General topology. Courier Dover Publications, 2017.
 [33] Vijay R Konda and John N Tsitsiklis. Actorcritic algorithms. In Advances in neural information processing systems, pages 1008–1014, 2000.
 [34] R Matthew Kretchmar. Parallel reinforcement learning. In The 6th World Conference on Systemics, Cybernetics, and Informatics. Citeseer, 2002.
 [35] Jonathan Lacotte, Yinlam Chow, Mohammad Ghavamzadeh, and Marco Pavone. Risksensitive generative adversarial imitation learning. arXiv preprint arXiv:1808.04468, 2018.
 [36] Dennis Lee, Haoran Tang, Jeffrey O Zhang, Huazhe Xu, Trevor Darrell, and Pieter Abbeel. Modular architecture for starcraft ii with deep reinforcement learning. In Fourteenth Artificial Intelligence and Interactive Digital Entertainment Conference, 2018.
 [37] Yuxi Li. Deep reinforcement learning: An overview. arXiv preprint arXiv:1701.07274, 2017.
 [38] An Liu, Vincent Lau, and Borna Kananian. Stochastic successive convex approximation for nonconvex constrained stochastic optimization. arXiv preprint arXiv:1801.08266, 2018.
 [39] Boyi Liu, Qi Cai, Zhuoran Yang, and Zhaoran Wang. Neural proximal/trust region policy optimization attains globally optimal policy. arXiv preprint arXiv:1906.10306, 2019.
 [40] Volodymyr Mnih, Adria Puigdomenech Badia, Mehdi Mirza, Alex Graves, Timothy Lillicrap, Tim Harley, David Silver, and Koray Kavukcuoglu. Asynchronous methods for deep reinforcement learning. In International conference on machine learning, pages 1928–1937, 2016.
 [41] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. Humanlevel control through deep reinforcement learning. Nature, 518(7540):529, 2015.
 [42] Arun Nair, Praveen Srinivasan, Sam Blackwell, Cagdas Alcicek, Rory Fearon, Alessandro De Maria, Vedavyas Panneershelvam, Mustafa Suleyman, Charles Beattie, Stig Petersen, et al. Massively parallel methods for deep reinforcement learning. arXiv preprint arXiv:1507.04296, 2015.
 [43] Santiago Paternain. Stochastic Control Foundations of Autonomous Behavior. PhD thesis, University of Pennsylvania, 2018.
 [44] Peng Peng, Ying Wen, Yaodong Yang, Quan Yuan, Zhenkun Tang, Haitao Long, and Jun Wang. Multiagent bidirectionallycoordinated nets: Emergence of humanlevel coordination in learning to play starcraft combat games. arXiv preprint arXiv:1703.10069, 2017.
 [45] Matteo Pirotta, Marcello Restelli, and Luca Bascetta. Policy gradient in lipschitz markov decision processes. Machine Learning, 100(23):255–283, 2015.
 [46] LA Prashanth and Mohammad Ghavamzadeh. Varianceconstrained actorcritic algorithms for discounted and average reward mdps. Machine Learning, 105(3):367–417, 2016.
 [47] Benjamin Recht. A tour of reinforcement learning: The view from continuous control. Annual Review of Control, Robotics, and Autonomous Systems, 2018.
 [48] Changhan Rhee and Peter W Glynn. Unbiased estimation with square root convergence for sde models. Operations Research, 63(5):1026–1043, 2015.
 [49] Andrzej Ruszczyński. Feasible direction methods for stochastic programming problems. Mathematical Programming, 19(1):220–229, 1980.
 [50] Joris Scharpff, Diederik M Roijers, Frans A Oliehoek, Matthijs TJ Spaan, and Mathijs Michiel de Weerdt. Solving transitionindependent multiagent mdps with sparse interactions. In AAAI, pages 3174–3180, 2016.
 [51] John Schulman, Sergey Levine, Pieter Abbeel, Michael Jordan, and Philipp Moritz. Trust region policy optimization. In International Conference on Machine Learning, pages 1889–1897, 2015.
 [52] Gesualdo Scutari, Francisco Facchinei, Peiran Song, Daniel P Palomar, and JongShi Pang. Decomposition by partial linearization: Parallel optimization of multiagent systems. IEEE Transactions on Signal Processing, 62(3):641–656, 2013.
 [53] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al. Mastering the game of go with deep neural networks and tree search. nature, 529(7587):484, 2016.
 [54] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, et al. A general reinforcement learning algorithm that masters chess, shogi, and go through selfplay. Science, 362(6419):1140–1144, 2018.
 [55] David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al. Mastering the game of go without human knowledge. Nature, 550(7676):354, 2017.
 [56] Peng Sun, Xinghai Sun, Lei Han, Jiechao Xiong, Qing Wang, Bo Li, Yang Zheng, Ji Liu, Yongsheng Liu, Han Liu, et al. Tstarbots: Defeating the cheating level builtin ai in starcraft ii in the full game. arXiv preprint arXiv:1809.07193, 2018.
 [57] Ying Sun, Prabhu Babu, and Daniel P Palomar. Majorizationminimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3):794–816, 2016.
 [58] Richard S Sutton and Andrew G Barto. Reinforcement learning: An introduction. MIT press, 2018.
 [59] Richard S Sutton, David A McAllester, Satinder P Singh, and Yishay Mansour. Policy gradient methods for reinforcement learning with function approximation. In Advances in neural information processing systems, pages 1057–1063, 2000.
 [60] Chen Tessler, Daniel J Mankowitz, and Shie Mannor. Reward constrained policy optimization. arXiv preprint arXiv:1805.11074, 2018.
 [61] Matteo Turchetta, Felix Berkenkamp, and Andreas Krause. Safe exploration in finite markov decision processes with gaussian processes. In Advances in Neural Information Processing Systems, pages 4312–4320, 2016.
 [62] Oriol Vinyals, Timo Ewalds, Sergey Bartunov, Petko Georgiev, Alexander Sasha Vezhnevets, Michelle Yeo, Alireza Makhzani, Heinrich Küttler, John Agapiou, Julian Schrittwieser, et al. Starcraft ii: A new challenge for reinforcement learning. arXiv preprint arXiv:1708.04782, 2017.
 [63] HoiTo Wai, Zhuoran Yang, Zhaoran Wang, and Mingyi Hong. Multiagent reinforcement learning via double averaging primaldual optimization. arXiv preprint arXiv:1806.00877, 2018.
 [64] Lingxiao Wang, Qi Cai, Zhuoran Yang, and Zhaoran Wang. Neural policy gradient methods: Global optimality and rates of convergence. arXiv preprint arXiv:1909.01150, 2019.
 [65] Min Wen and Ufuk Topcu. Constrained crossentropy method for safe reinforcement learning. In Advances in Neural Information Processing Systems, pages 7461–7471, 2018.
 [66] Sijia Xu, Hongyu Kuang, Zhi Zhuang, Renjie Hu, Yang Liu, and Huyang Sun. Macro action selection with deep reinforcement learning in starcraft. arXiv preprint arXiv:1812.00336, 2018.
 [67] Yang Yang, Gesualdo Scutari, Daniel P Palomar, and Marius Pesavento. A parallel decomposition method for nonconvex stochastic multiagent optimization problems. IEEE Transactions on Signal Processing, 64(11):2949–2964, 2016.
 [68] Zhuoran Yang, Yongxin Chen, Mingyi Hong, and Zhaoran Wang. On the global convergence of actorcritic: A case for linear quadratic regulator with ergodic cost. arXiv preprint arXiv:1907.06246, 2019.
 [69] Kaiqing Zhang, Zhuoran Yang, Han Liu, Tong Zhang, and Tamer Başar. Finitesample analyses for fully decentralized multiagent reinforcement learning. arXiv preprint arXiv:1812.02783, 2018.
Appendix A Other applications
A.1 Constrained Parallel Markov Decision Process
We consider the parallel MDP problem [34, 42, 17] where we have a singleagent MDP task and $N$ workers, where each worker acts as an individual agent and aims to solve the same MDP problem. In the parallel MDP setting, each agent is characterized by a tuple $(\mathcal{S},\mathcal{A},P,\gamma ,{r}^{i},{d}^{i},{\mu}^{i})$, where each agent has the same but individual state space, action space, transition probability distribution, and the discount factor. However, the reward function, cost function, and the distribution of the initial state ${s}_{0}\in \mathcal{S}$ could be different for each agent, but satisfy $\mathbb{E}[{r}^{i}(s,a)]=r(s,a)$, $\mathbb{E}[{d}^{i}(s,a)]=d(s,a)$, and $\mathbb{E}[{\mu}^{i}(s,a)]=\mu (s,a)$. Each agent $i$ generates its own trajectory $\{{s}_{0}^{i},{a}_{0}^{i},{s}_{1}^{i},{a}_{1}^{i},\mathrm{\dots}\}$ and collects its own reward/cost value $\{{r}_{0}^{i},{d}_{0}^{i},{r}_{1}^{i},{d}_{1}^{i},\mathrm{\dots}\}$.
The hope is that by solving the singleagent problem using $N$ agents in parallel, the algorithm could be more stable and converge much faster [40]. Intuitively, each agent $i$ may have a different initial state and will explore different parts of the state space due to the randomness in the state transition distribution and the policy. It also helps to reduce the correlation between agents’ behaviors. As a result, by running multiple agents in parallel, we are more likely to visit different parts of the environment and get the experience of the reward/cost function values more efficiently. This mimics the strategy used in treebased supervised learning algorithms [14, 29, 30].
Following the settings in [17], we have $N$ agents (i.e., $N$ workers) and one central controller in the system. The global parameter is denoted by $\theta $, and we consider the constrained parallel MDP problem where the goal is to solve the following optimization problem:
$\underset{\theta}{\text{minimize}}J(\theta )={\displaystyle \sum _{i=1}^{N}}{\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot {r}^{i}({s}_{t}^{i},{a}_{t}^{i})\right],$  (28)  
$\text{subject to}D(\theta )={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot {d}^{i}({s}_{t}^{i},{a}_{t}^{i})\right]\le {D}_{0},i\in \mathcal{N}.$ 
During the estimation step, the controller broadcasts the current parameter ${\theta}_{k}$ to each agent and each agent samples its own trajectory and obtains estimators for function value/gradient of the reward/cost function. Next, each agent uploads its estimators to the central controller and the central controller takes the average over these estimators, and then follow our proposed algorithm to solve for the QCQP problem and update the parameter to ${\theta}_{k+1}$. This process continues until convergence.
A.2 Constrained Multiagent Markov Decision Process
A natural extension of the (singleagent) MDP is to consider a model with $N$ agents termed multiagent Markov decision process (MMDP). Recently this kind of problem has been attracting more and more attention. See [15] for a comprehensive survey. Most of the work on multiagent MDP problems consider the setting where the agents share the same global state space, but each with their own collection of actions and rewards [11, 63, 69]. In each stage of the system, each agent observes the global state and chooses its own action individually. As a result, each agent receives its reward and the state evolves according to the joint transition distribution. An MMDP problem can be fully collaborative where all the agents have the same goal, or fully competitive where the problem consists of two agents with an opposite goal, or the mix of the two.
Here we consider a slightly different setting where each agent has its own state space. The only connection between the agents is that the global reward is a function of the overall states and actions. Furthermore, each agent has its own constraint which depends on its own state and action only. This problem is known as TransitionIndependent Multiagent MDP and is considered in [50]. Specifically, each agent’s task is characterized by a tuple $({\mathcal{S}}^{i},{\mathcal{A}}^{i},{P}^{i},\gamma ,{d}^{i},{\mu}^{i})$ with each component defined as usual. Note that ${P}^{i}:{\mathcal{S}}^{i}\times {\mathcal{A}}^{i}\to \mathcal{P}({\mathcal{S}}^{i})$ and ${d}^{i}:{\mathcal{S}}^{i}\times {\mathcal{A}}^{i}\to \mathbb{R}$ are functions of each agent’s state and action only and do not depend on other agents. Denote $\mathcal{S}={\mathrm{\Pi}}_{i\in \mathcal{N}}{\mathcal{S}}^{i}$ and $\mathcal{A}={\mathrm{\Pi}}_{i\in \mathcal{N}}{\mathcal{A}}^{i}$ as the joint state space and action space. The global reward function is given by $r:\mathcal{S}\times \mathcal{A}\to \mathbb{R}$ that depends on the joint state and action. Here we consider the fully collaborative setting where all the agents have the same goal. Under this setting, the policy set of each agent is parameterized as $\{{\pi}_{{\theta}^{i}}^{i}:{\mathcal{S}}^{i}\to \mathcal{P}({\mathcal{A}}^{i})\}$ and we denote $\theta =[{\theta}^{1},\mathrm{\dots},{\theta}^{N}]$ as the overall parameters and ${\pi}_{\theta}$ as the overall policy. In the following, we use $\mathcal{N}=\{1,2,\mathrm{\dots},N\}$ to denote the $N$ agents. Denote ${a}_{t}^{i}$ as the action chosen by agent $i$ at stage $t$ and ${a}_{t}={\mathrm{\Pi}}_{i\in \mathcal{N}}{a}_{t}^{i}$ as the joint action chosen by all the agents. The goal of this constrained MMDP is to solve the following problem
$\underset{\theta}{\text{minimize}}J(\theta )={\mathbb{E}}_{{\pi}_{\theta}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot r({s}_{t},{a}_{t})\right],$  (29)  
$\text{subject to}{D}^{i}({\theta}^{i})={\mathbb{E}}_{{\pi}_{{\theta}^{i}}}\left[{\displaystyle \sum _{t\ge 0}}{\gamma}^{t}\cdot {d}^{i}({s}_{t}^{i},{a}_{t}^{i})\right]\le {D}_{0}^{i},i\in \mathcal{N}.$ 
Inspired by the parallel implementation ([38], Section V), our algorithm applies naturally to constrained MMDP problem with some modifications. This modified procedure can also be viewed as a distributed version of the original algorithm. The overall problem (29) can be viewed as a large “singleagent" problem where the constraints are decomposable into $N$ parts. In this case, instead of solving a large QCQP problem in each iteration, each agent could solve its own QCQP problem in a distributed manner which is much more efficient. As before, we denote the sample negative reward and cost function as
$${J}^{*}(\theta )=\sum _{t\ge 0}{\gamma}^{t}\cdot r({s}_{t},{a}_{t})\mathit{\hspace{1em}\hspace{1em}}\text{and}\mathit{\hspace{1em}\hspace{1em}}{D}^{i,*}({\theta}^{i})=\sum _{t\ge 0}{\gamma}^{t}\cdot {d}^{i}({s}_{t}^{i},{a}_{t}^{i}).$$ 
In each iteration with ${\theta}_{k}=[{\theta}_{k}^{1},\mathrm{\dots},{\theta}_{k}^{N}]$, we approximate $J(\theta )$ and $D(\theta )$ as
${\stackrel{~}{J}}^{i}({\theta}^{i},{\theta}_{k},\tau )$  $={\displaystyle \frac{1}{N}}{J}^{*}({\theta}_{k})+\u27e8{\nabla}_{{\theta}^{i}}{J}^{*}({\theta}_{k}),{\theta}^{i}{\theta}_{k}^{i}\u27e9+\tau {\parallel {\theta}^{i}{\theta}_{k}^{i}\parallel}_{2}^{2},$  (30)  
${\stackrel{~}{D}}^{i}({\theta}^{i},{\theta}_{k},\tau )$  $={D}^{i,*}({\theta}_{k}^{i})+\u27e8{\nabla}_{{\theta}^{i}}{D}^{i,*}({\theta}_{k}^{i}),{\theta}^{i}{\theta}_{k}^{i}\u27e9+\tau {\parallel {\theta}^{i}{\theta}_{k}^{i}\parallel}_{2}^{2}.$  (31) 
Note that the constraint function is naturally decomposable into $N$ parts. We also “manually" split the objective function into $N$ parts, so that each agent could solve its own QCQP problem in a distributed manner. As before, we define
${\overline{J}}^{i,(k)}({\theta}^{i})=(1{\rho}_{k})\cdot {\overline{J}}^{i,(k1)}({\theta}^{i})+{\rho}_{k}\cdot {\stackrel{~}{J}}^{i}({\theta}^{i},{\theta}_{k},\tau ),$  
${\overline{D}}^{i,(k)}({\theta}^{i})=(1{\rho}_{k})\cdot {\overline{D}}^{i,(k1)}({\theta}^{i})+{\rho}_{k}\cdot {\stackrel{~}{D}}^{i}({\theta}^{i},{\theta}_{k},\tau ).$ 
With this surrogate functions, each agent then solves its own convex relaxation problem
$${\overline{\theta}}_{k}^{i}=\underset{{\theta}^{i}}{argmin}{\overline{J}}^{i,(k)}({\theta}^{i})\mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}\hspace{1em}}{\overline{D}}^{i,(k)}({\theta}^{i})\le {D}_{0}^{i},$$  (32) 
or, alternatively, solves for the feasibility problem if (32) is infeasible
$${\overline{\theta}}_{k}^{i}=\underset{{\theta}^{i},{\alpha}^{i}}{argmin}{\alpha}^{i}\mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}\hspace{1em}}{\overline{D}}^{i,(k)}({\theta}^{i})\le {D}_{0}^{i}+{\alpha}^{i}.$$  (33) 
This step can be implemented in a distributed manner for each agent and is more efficient than solving the overall problem with the overall parameter $\theta $. Finally, the update rule for each agent $i$ is as usual
$${\theta}_{t+1}^{i}=(1{\eta}_{k})\cdot {\theta}_{k}^{i}+{\eta}_{k}\cdot {\overline{\theta}}_{k}^{i}.$$ 
This process continues until convergence.
Appendix B Proof of Theorem 4
According to the choice of the surrogate function (9) and Assumption 2, it is straightforward to verify that the function ${\overline{J}}^{(k)}(\theta )$ defined in (12) is uniformly strongly convex in $\theta $ for each iteration $t$. Moreover, both ${\overline{J}}^{(k)}(\theta )$ and ${\nabla}_{\theta}{\overline{J}}^{(k)}(\theta )$ are Lipschitz continuous functions.
From Lemma 1 in [49] we have
$$\underset{t\to \mathrm{\infty}}{lim}\left{\overline{J}}^{(k)}(\theta )\mathbb{E}\left[\stackrel{~}{J}(\theta ,{\theta}_{k},\tau )\right]\right=0.$$ 
Since the function $\mathbb{E}\left[\stackrel{~}{J}(\theta ,{\theta}_{k},\tau )\right]$ is Lipschitz continuous in ${\theta}_{k}$, we obtain that
$$\left{\overline{J}}^{({k}_{1})}(\theta ){\overline{J}}^{({k}_{2})}(\theta )\right\le {L}_{0}\cdot \parallel {\theta}_{{k}_{1}}{\theta}_{{k}_{2}}\parallel +\u03f5,$$ 
for some constant ${L}_{0}$ and the error term $\u03f5$ that goes to $0$ as ${k}_{1},{k}_{2}$ go to infinity. This shows that the function sequence ${\overline{J}}^{({k}_{j})}(\theta )$ is equicontinuous. Since $\mathrm{\Theta}$ is compact and the discounted cumulative reward function is bounded by ${r}_{\mathrm{max}}/(1\gamma )$, we can apply ArzelaAscoli theorem [22, 32] to prove existence of $\widehat{J}(\theta )$ that converges uniformly. Moreover, since we apply the same operations on the constraint function $D(\theta )$ as to the reward function $J(\theta )$ in Algorithm 3, the above properties also hold for $D(\theta )$.
The rest of the proof follows in a similar way as the proof of Theorem 1 in [38]. Under Assumptions 1  3, the technical conditions in [38] are satisfied by the choice of the surrogate functions (9) and (10). According to Lemma 2 in [38], with probability one we have
$$\underset{k\to \mathrm{\infty}}{limsup}D({\theta}_{k})\le {D}_{0}.$$ 
This shows that, although in some of the iterations the convex relaxation problem (11) is infeasible, and we have to solve the alternative problem (14), the iterates $\{{\theta}_{k}\}$ converge to the feasible region of the original problem (5) with probability one. Furthermore, with probability one, the convergent point $\stackrel{~}{\theta}$ is the optimal solution to the following problem
$\underset{\theta \in \mathrm{\Theta}}{\text{minimize}}\widehat{J}(\theta )\mathit{\hspace{1em}\hspace{1em}}\text{subject to}\mathit{\hspace{1em}\hspace{1em}}\widehat{D}(\theta )\le {D}_{0}.$  (34) 
The KKT conditions for (34) together with the Slater condition show that the KKT conditions of the original problem (5) are also satisfied at $\stackrel{~}{\theta}$. This shows that $\stackrel{~}{\theta}$ is a stationary point of the original problem almost surely.