Title: CONSTRAINED MARKOV DECISION PROCESSES VIA BACKWARD VALUE FUNCTIONS

Abstract: Although Reinforcement Learning (RL) algorithms have found tremendous success in simulated domains, they often cannot directly be applied to physical systems, especially in cases where there are hard constraints to satisfy (e.g. on safety or resources). In standard RL, the agent is incentivized to explore any behavior as long as it maximizes rewards, but in the real world undesired behavior can damage either the system or the agent in a way that breaks the learning process itself. In this work, we model the problem of learning with constraints as a Constrained Markov Decision Process, and provide a new on-policy formulation for solving it. A key contribution of our approach is to translate cumulative cost constraints into state-based constraints. Through this, we define a safe policy improvement method which maximizes returns while ensuring that the constraints are satisfied at every step. We provide theoretical guarantees under which the agent converges while ensuring safety over the course of training. We also highlight computational advantages of this approach. The effectiveness of our approach is demonstrated on safe navigation tasks and in safety-constrained versions of MuJoCo environments, with deep neural networks.

1. INTRODUCTION:

Reinforcement Learning (RL) provides a sound decision-theoretic framework to optimize the behavior of learning agents in an interactive setting (Sutton & Barto, 2018) . Recently, the field of RL has found success in many high-dimensional domains, like video games, Go, robot locomotion and navigation. However, most of the success of RL algorithms has been limited to simulators, where the learning algorithm has the ability to reset the simulator. In the physical world, an agent will need to avoid harmful behavior (e.g. damaging the environment or the agent's hardware) while learning to explore behaviors that maximize the reward.

A few popular approaches for avoiding undesired behaviors for high-dimensional systems include reward-shaping (Moldovan & Abbeel, 2012) , reachability-preserving algorithms (Mitchell, 2003; Eysenbach et al., 2017) , state-level surrogate constraint satisfaction algorithms (Dalal et al., 2018) , risk-sensitive algorithms (Tamar et al., 2013; Chow et al., 2015) and apprenticeship learning (Abbeel & Ng, 2004) . There also exists model-based Bayesian approaches that are focused on imposing the constraints via the dynamics (such as classifying parts of state space as unsafe) and then using model predictive control to incorporate the constraints in the policy optimization and planning (Turchetta et al., 2016; Berkenkamp et al., 2017; Wachi et al., 2018; Koller et al., 2018) . A natural way to model safety is via constraint satisfaction. A standard formulation for adding constraints to RL problems is the Constrained Markov Decision Process (CMDP) framework (Altman, 1999) , wherein the environment is extended to also provide feedback on constraint costs. The agent must then attempt to maximize its expected return while also satisfying cumulative constraints.

A few algorithms have been proposed to solve CMDPs for high-dimensional domains with continuous action spaces -however they come with their own caveats. Reward Constrained Policy Optimization (Tessler et al., 2018) and Primal Dual Policy Optimization (Chow et al., 2015) do not guarantee constraint satisfaction during the learning procedure, only on the final policy. Constrained Policy Optimization (Achiam et al., 2017) provides monotonic policy improvement but is computationally expensive due to requiring a backtracking line-search procedure and conjugate gradient algorithm for approximating the Fisher Information Matrix. Lyapunov-based Safe Policy Optimization (Chow et al., 2019) requires solving a Linear Program (LP) at every step of policy evaluation, although they show that there exists heuristics which can be substituted for the LP at the expense of theoretical guarantees.

In this work, we propose an alternate formulation for solving CMDPs that transforms trajectory-level constraints into localized state-dependent constraints, through which a safe policy improvement step can be defined. In our approach, we define a notion of Backward Value Functions, which act as an estimator of the expected cost collected by the agent so far and can be learned via standard RL bootstrap techniques. We provide conditions under which this new formulation is able to solve CMDPs without violating the constraints during the learning process. Our formulation allows us to define state-level constraints without explicitly solving a LP or the Dual problem at every iteration. Our method is implemented as a reduction to any model-free on-policy bootstrap based RL algorithm, both for deterministic and stochastic policies, and discrete and continuous action spaces. We provide the empirical evidence of our approach with Deep RL methods on various safety benchmarks, including 2D navigation grid worlds (Leike et al., 2017; Chow et al., 2018) , and MuJoCo tasks (Achiam et al., 2017; Chow et al., 2019) .

2. CONSTRAINED MARKOV DECISION PROCESSES:

We write P(Y ) for the set of probability distributions on a space Y . A Markov Decision Process (MDP) (Puterman, 2014) is a tuple (X , A, P, r, x 0 ), where X is a set of states, A is a set of actions, r : X × A → [0, R M AX ] is a reward function, P : X × A → P(X ) is a transition probability function, and x 0 is a fixed starting state. For simplicity we assume a deterministic reward function and starting state, but our results generalize.

A Constrained Markov Decision Process (CMDP) (Altman, 1999 ) is a MDP with additional constraints that restrict the set of permissible policies for the MDP. Formally, a CMDP is a tuple (X , A, P, r, x 0 , d, d 0 ), where d : X → [0, D M AX ] is the cost function1 and d 0 ∈ R ≥0 is the maximum allowed cumulative cost. The set of feasible policies that satisfy the CMDP is the subset of stationary policies Π D := {π :

X → P(A) E[ T t=0 d(x t ) | x 0 , π] ≤ d 0 }.

We consider a finite time horizon T after which the episode terminates. The expected sum of rewards following a policy π from an initial state x is given by the value function

V π (x) = E[ T t=0 r(x t , a t ) | π, x].

Analogously, the expected sum of costs is given by the cost value function

V π D (x) = E[ T t=0 d(x t ) | π, x].

The RL problem in the CMDP is to find the feasible policy which maximizes expected returns from the initial state x 0 , i.e.

π * = arg max

π∈Π D V π (x 0 )

An important point to note about CMDPs is that, in the original formulation, the cost function depends on immediate states but the constraint is cumulative and thus depends on the entire trajectory.

In the case of MDPs, where a model of the environment is not known or is not easily obtained, it is still possible for the agent to find the optimal policy using Temporal Difference (TD) methods (Sutton, 1988) . Broadly, these methods update the estimates of the value functions via bootstraps of previous estimates on sampled transitions (we refer the reader to Sutton & Barto (2018) for more information). In the on-policy setting, we alternate between estimating the state-action value function Q π for a given π and updating the policy to be greedy with respect to the value function.

3. SAFE POLICY ITERATION VIA BACKWARD VALUE FUNCTIONS:

Our approach proposes to convert the trajectory-level constraints of the CMDP into single-step state-wise constraints in such a way that satisfying the state-wise formulation will entail satisfying the original trajectory-level problem. The advantages of this approach are twofold: i) working with single-step state-wise constraints allows us to obtain analytical solutions to the optimization problem, and ii) the state-wise constraints can be defined via value-function-like quantities and can thus be estimated with well-studied value-based methods. The state-wise constraints are defined via Backward Value Functions, in Section 3.2, and in Section 3.3 we provide a safe policy iteration procedure which satisfies said constraints (and thus the original problem).

3.1. BACKWARD MARKOV CHAIN:

Unlike in traditional RL, in the CMDP setting the agent needs to take into account the constraints which it has accumulated so far in order to plan accordingly for the future. Intuitively, the accumulated cost so far can be estimated via the cost value function V D running "backward in time". Before giving the details of our approach and formally introducing the Backward Value Functions, we review the main ideas, which are built upon the work of Morimura et al. (2010) , who also considered time-reversed Markov chains but from the standpoint of estimating the gradient of the log stationary distribution; we extend these ideas to TD methods. Assumption 3.1 (Stationarity). The MDP is ergodic for any policy π, i.e., the Markov chain characterized by the transition probability P π (x t+1 |x t ) = at∈A P(x t+1 |x t , a t )π(a t |x t ) is irreducible and aperiodic.

Let M(π) denote the Markov chain characterized by transition probability P π (x t+1 |x t ). The above assumption implies that there exists a unique stationary distribution η π associated with π, such that it satisfies:

η π (x t+1 ) = xt∈X P π (x t+1 |x t )η π (x t ).

We abuse the notation and denote

P π (x t+1 , a t |x t ) = P(x t+1 |x t , a t )π(a t |x t ).

According to Bayes' rule, the probability q(x t-1 , a t-1 |x t ) of a previous state-action pair (x t-1 , a t-1 ) leading to the current state x t is given by:

q(x t-1 , a t-1 |x t ) = P(x t |x t-1 , a t-1 )P r(x t-1 , a t-1 )

xt-1∈X at-1∈A P(x t |x t-1 , a t-1 )P r(x t-1 , a t-1 )

.

From Assumption 3.1, we have that P r(x t-1 , a t-1 ) = η π (x t-1 )π(a t-1 |x t-1 ), and

xt-1∈X at-1∈A P(x t |x t-1 , a t-1 )P r(x t-1 , a t-1 ) = η π (x t ). We denote the posterior q(x t-1 , a t-1 |x t ) as backward (or time-reversed) probability # » P π (x t-1 , a t-1 |x t ), and we have:

EQUATION

The forward Markov chain, characterized by the transition matrix P π (x t+1 |x t ), runs forward in time, i.e., it gives the probability of the next state in which the agent will end up. Analogously, a backward Markov chain is denoted by the transition matrix 

# » P π (x t-1 |x t ) = at-1∈A # » P π (x t

3.2. BACKWARD VALUE FUNCTION:

We define the Backward Value Function (BVF) to be a value function running on the backward Markov chain # » B (π). A BVF is the expected sum of returns or costs collected by the agent so far. We are mainly interested in maintaining estimates of the cumulative cost incurred at a state in order to express the total constraint state-wise.

We note that, since every Markov chain M(π) is ergodic by Assumption 3.1, the corresponding backward Markov chain B(π) is also ergodic (Morimura et al., 2010, Prop. B.1) . In particular, every policy π can reach the initial state via some path in the transition graph of the backward Markov chain. Thus, the backwards Markov chain are also finite-horizon for some T B , with x 0 corresponding to the terminal state. We define a finite-horizon Backward Value Function for cost as:

EQUATION

Proposition 3.1 (Sampling). Samples from the forward Markov chain M(π) can be used directly to estimate the statistics of the backward Markov chain # » B (π) (or the Backward Value Function). We have:

EQUATION

= E M(π) K k=0 d(x t+k )|x t+K , η π (x t )

,

where E M(π) and E # » B (π) are expectations over the forward and backward chains respectively. The Equation (3) holds true even in the limit K → ∞.

The proof is given in Appendix B.1. Using the above proposition, we get an interchangeability property that removes the need to sample from the backward chain. We can use the traditional RL setting and draw samples from the forward chain and still estimate the BVFs. Equation ( 2) can be written recursively as:

# » V π D (x t ) = E # » B (π) d(x t ) + # » V π D (x t-1 ) .

In operator form, the above equation can also be written as:

EQUATION

)

Proposition 3.2 (Fixed point). For a policy π, the associated Backward Value Function vector,

# » V π , satisfies lim k→∞ # » T π k # » V = # » V π for every vector # » V , and # » V π is the unique solution of the equation # » V π = # » T π # » V π .

The proof is given in Appendix B.2. The above proposition allows us to soundly extend the RL methods based on Bellman operators for the estimation of BVFs.

3.3. SAFE POLICY IMPROVEMENT VIA BVF-BASED CONSTRAINTS:

With the Backward Value Function framework, the trajectory-level optimization problem associated with a CMDP can be rewritten in state-wise form. Recall that a feasible policy must satisfy the constraint:

E M(π) T k=0 d(x k ) | x 0 ≤ d 0 .

Alternatively, for each timestep t ∈ [0, T ] of a trajectory:

E t k=0 d(x k ) | x 0 , π + E T k=t d(x k ) | x 0 , π -E d(x t ) | x 0 ≤ d 0 . Via the identities E[ T k=t d(x k ) | x 0 , π] ≤ E xt∼δx 0 (P π ) t [V π D (x t )] and E[ t k=0 d(x k ) | x 0 , π] ≤ E x k ∼δx 0 (P π ) t [ # » V π D (x t )](derived in Appendix C) 2

, we remark that the quantity on the LHS is less than the expectation over k-step trajectories of

# » V π D (x t ) + V π D (x t ) -d(x t ).

In other words, for each t ∈ [0, T ] :

E M(π) T k=0 d(x k ) | x 0 ≤ E xt∼δx 0 (P π ) t # » V π D (x t ) + V π D (x t ) -d(x t ) ≤ d 0 .

These are the state-wise constraints that should hold at each step in a given trajectory -we refer to them as the value-based constraints. Satisfying the value-based constraints will automatically satisfy the given CMDP constraints.

This formulation allows us to introduce a policy improvement step, which maintains a safe feasible policy at every iteration by using the previous estimates of the forward and backward value functions3 . The policy improvement step is defined by a linear program, which performs a greedy update with respect to the current state-action value function subject to the value-based constraints:

π k+1 (•|x) = arg max π∈Π π(•|x), Q π k (x, •) , (SPI) s.t. π(•|x), Q π k D (x, •) + # » V π k D (x) -d(x) ≤ d 0 , ∀x ∈ X .

Our first result is that the policies obtained by the policy improvement step will satisfy the safety constraints. We write TV(•, •) for the total variation metric between distributions.

Theorem 3.1 (Consistent Feasibility). Assume that successive policies are updated sufficiently

slowly, i.e. TV(π k+1 (•|x), π k (•|x)) ≤ d0-V π k D (x0)

2DMAXT 2 . 4 Then the policy iteration step given by (SPI) is consistently feasible, i.e. if π k is feasible at x 0 then so is π k+1 .

It is also possible to consider larger neighbourhoods for updates of successive policies, but at the cost of everywhere-feasibility. For want of space, we present that result in Appendix D.

Next we show that the policy iteration step given by (SPI) leads to monotonic improvement.

Theorem 3.2 (Policy Improvement). Let π n and π n+1 be successive policies generated by the policy iteration step of (SPI). Then V πn+1 (x) ≥ V πn (x) ∀x ∈ X . In particular, the sequence of value functions {V πn } n≥0 given by (SPI) monotonically converges.

Proofs for Theorems 3.1 and 3.2 are given in Appendix D. Finding the sub-optimality gap (if any) remains an interesting question left for future work.

4. PRACTICAL IMPLEMENTATION CONSIDERATIONS:



4.1. DISCRETE ACTION SPACE:

In discrete action spaces, the problem in (SPI) can be solved exactly as a Linear Programming problem. It is possible to approximate its analytical solution by casting it into the corresponding entropy-regularized counterpart (Neu et al., 2017; Chow et al., 2018) . The details of the closed form solution can be found in Appendix E.

Furthermore, if we restrict the set of policies to be deterministic, then it is possible to have an in-graph solution as well. The procedure then closely resembles the Action Elimination Procedure (Puterman, 2014, Chapter 6) , where non-optimal actions are identified as being those which violate the constraints.

4.2. EXTENSION TO CONTINUOUS CONTROL:

For MDPs with only state-dependent costs, Dalal et al. (2018) proposed the use of safety layers, a constraint projection approach, that enables action correction at each step. At any given state, an unconstrained action is selected and is passed to the safety layer, which projects the action to the nearest action (in Euclidean norm) satisfying the necessary constraints. We extend this approach to stochastic policies to handle the corrections for the actions generated by stochastic policies. When the policy is parameterized with a Gaussian distribution, then the safety-layer can still be used by projecting both the mean and standard-deviation vector to the constraint-satisfying hyper-plane5 . In most cases, the standard-deviation vector is kept fixed or independent of the state (Kostrikov, 2018; Dhariwal et al., 2017) , which allows us to formulate the problem as solving the following L2-projection of the mean of the Gaussian in Euclidean space. For µ π (.; θ), at any given state x ∈ X , the safety layer solves the following projection problem:

arg min µ 1 2 (µ -µ π (x) 2 , s.t. Q π D (x, µ) + # » V π D (x) -d(x) ≤ d 0 .

As shown in Dalal et al. (2018) ; Chow et al. (2019) , if the constraints have linear nature then an analytical solution exists. In order to get a linearized version of the constraints (and simplify the projection), we can approximate the constraint with its first-order Taylor series at µ = µ π (x):

EQUATION

s.t. # » V π D (x) -d(x) + Q π D (x, µ π (x)) + (µ -µ π (x)) T (∇Q π D (x, µ)| µ=µπ(x) )

First order Taylor expansion

≤ d 0 .

The above objective function is positive-definite and quadratic, and the constraints are linear. Though this problem can be solved by an in-graph QP solver, there exists an analytical solution (see Appendix G): Proposition 4.1. At a given state x ∈ X , the solution to the Eq. ( 5), µ * is:

µ * = µ π (x) -λ * (x) • g µ,D (x), where, g µ,D (x) = ∇Q π D (x, µ)| µ=µπ(x) , λ * (x) = -(d 0 + d(x) - # » V π D (x) -Q π D (x, µ π (x))) g µ,D (x) T g µ,D (x) + .

5. RELATED WORK:

Lagrangian-based methods: Initially introduced in Altman (1999), more scalable versions of the Lagrangian based methods have been proposed over the years (Moldovan & Abbeel, 2012; Tessler et al., 2018; Chow et al., 2015) . The general form of the Lagrangian methods is to convert the problem to an unconstrained problem via Langrange multipliers. If the policy parameters are denoted by θ, then Lagrangian formulation becomes:

min λ≥0 max θ (L(θ, λ) = min λ≥0 max θ [V π θ (x 0 ) -λ(V π θ D (x 0 ) -d 0 ))] ,

where L is the Lagrangian and λ is the Lagrange multiplier (penalty coefficient). The main problems of the Lagrangian methods are that the Lagrangian multiplier is either a hyper-parameter (without much intuition), or is solved on a lower time-scale. That makes the unconstrained RL problem a three time-scale6 problem, which makes it very difficult to optimize in practice. Another problem is that during the optimization, this procedure can violate the constraints. Ideally, we want a method that can respect the constraint throughout the training and not just at the final optimal policy.

Lyapunov-based methods: In control theory, the stability of the system under a fixed policy is computed using Lyapunov functions (Khalil, 1996) . A Lyapunov function is a type of scalar potential function that keeps track of the energy that a system continually dissipates. Recently, Chow et al. (2018; 2019) provide a method of constructing the Lyapunov functions to guarantee global safety of a behavior policy using a set of local linear constraints. Their method requires the knowledge of T V (π, π * ) to guarantee the theoretical claims. They substitute the ideally required Lyapunov function with an approximate solution that requires solving a LP problem at every iteration. For the practical scalable versions, they use a heuristic, a constant Lyapunov function for all states that only depends on the initial state and the horizon. While our methods also constructs state-wise constraints, there are two notable differences: a) our assumption only rely on the current policy candidate and the baseline policy, instead of the baseline and the optimal policy, b) our method does not require solving an LP at every update step to construct the constraint and as such the only approximation error that is introduced comes from the function approximation.

Conservative Policy Improvement: Constrained Policy Optimization (CPO) (Achiam et al., 2017) extends the trust-region policy optimization (Schulman et al., 2015) algorithm to satisfy constraints during training as well as after convergence. CPO is computationally expensive as it uses an approximation to the Fisher Information Matrix which requires many steps of conjugate gradient descent (n cg steps) followed by a backtracking line-search procedure (n ls steps) for each iteration, so it is more expensive by O(n cg + n ls ) per update. Furthermore, accurately estimating the curvature requires a large number of samples in each batch (Wu et al., 2017) .

6. EXPERIMENTS:

We empirically validate our approach on RL benchmarks to measure the performance of the agent with respect to the accumulated return and cost during training in the presence of neural-networks based function approximators. We compare our approach with the respective Unconstrained versions, and the Lyapunov-based approach (Chow et al., 2018; 2019) in each setting. Even though our formulation is based on the undiscounted case, we use discounting with γ = 0.99 for estimating the value functions in order to be consistent with the baselines.7 6.1 STOCHASTIC GRID WORLD Motivated by the safety in navigation tasks, we first consider a stochastic 2D grid world (Leike et al., 2017; Chow et al., 2018) . The agent (green cell in Fig. 1a ) starts in the bottom-right corner, the safe region, and the objective is to move to the goal on the other side of the grid (blue cell). The agent can only move in the adjoining cells in the cardinal directions. It gets a reward of +1000 on reaching the goal, and a penalty of -1 at every timestep. Thus, the task is to reach the goal in the shortest amount of time. There are a number of pits in the terrain (red cells) that represent the safety constraint and the agent gets a cost of 10 on passing through any pit cell. Occasionally, with probability p = 0.05, a random action will be executed instead of the one selected by the agent. Thus, the task is to reach to the goal in the shortest amount of time, while passing through the red grids at most d 0 /10 times. The size of the grid is 12 × 12 cells, and the pits are randomly generated for each grid with probability ρ = 0.3. The agent starts at (12, 12) and the goal is selected uniformly on (α, 0), where α ∼ U (0, 12). The threshold d 0 = 20 implies the agent can pass at most two pits. The maximum horizon is 200 steps, after which the episode terminates.

We use the action elimination procedure described in Sec 4.1 in combination with n-step SARSA (Rummery & Niranjan, 1994; Peng & Williams, 1994) using neural networks and multiple synchronous agents as in Mnih et al. (2016) . We use -greedy exploration. The results are shown in Fig. 1 (more experimental details can be found in Appendix H). We observe that the agent is able to respect the safety constraints more adequately than the Lyapunov-based method, albeit at the expense of some decrease in return, which is the expected trade-off for satisfying the constraints.

6.2. MUJOCO BENCHMARKS:

Based on the safety experiments in Achiam et al. (2017) ; Chow et al. (2019) , we design three simulated robot locomotion continuous control tasks using the MuJoCo simulator (Todorov et al., 2012) and OpenAI Gym (Brockman et al., 2016) : (1) Point-Gather: A point-mass agent (S ⊆ R 9 , A ⊆ R 2 ) is rewarded for collecting the green apples and constrained to avoid the red bombs; (2) Safe-Cheetah: A bi-pedal agent (S ⊆ R 18 , A ⊆ R 6 ) is rewarded for running at high speed, but at the same time constrained by a speed limit; (3) Point-Circle: The point-mass agent (S ⊆ R 9 , A ⊆ R 2 ) is rewarded for running along the circumference of a circle in counter-clockwise direction, but is constrained to stay within a safe region smaller than the radius of the circle.

We integrate our method on top of the A2C algorithms (Mnih et al., 2016) and PPO (Schulman et al., 2017) , using the procedure described in Section 4.2. More details about the tasks and network architecture can be found in the Appendix I. Algorithmic details can be found in Appendix J. The results with A2C are shown in Fig. 2 and the results with PPO are shown in Fig. 3 . We observe that our Safe method is able to respect the safety constraint throughout most of the learning, and with much greater degree of compliance than the Lyapunov-based method, especially when combined with A2C. The one case where the Safe method fails to respect the constraint is in Point-Circle with PPO (Fig. 3(c) ). Upon further examination, we note that the training in this scenario has one of two outcomes: some runs end with the learner in an infeasible set of states from which it cannot recover; other runs end in a good policy that respects the constraint. We discuss solutions to overcome this in the final section. 

7. DISCUSSION:

We present a method for solving constrained MDPs that respects trajectory-level constraints by converting them into state dependent value-based constraints, and show how the method can be used to handle safety limitations in both discrete and continuous spaces. The main advantage of our approach is that the optimization problem is more easily solved with value-based constraints, while providing similar guarantees and requiring less approximations. The empirical results presented show that our approach is able to solve the tasks with good performance while maintaining safety throughout training. It is important to note that there is a fundamental trade-off between exploration and safety. It is impossible to be 100% safe without some knowledge; in cases where that knowledge is not provided a priori, it must be acquired through exploration. We see this in some of our results (Gridworld, Point-Circle) where our safe policy goes above the constraint in the very early phases of training (all our experiments started from a random policy). We note that the other methods also suffer from this shortcoming. An open question is how to provide initial conditions or a priori knowledge, to avoid this burn-in phase. Another complementary strategy to explore is for cases where an agent is stuck in an unsafe or infeasible policy space, where a recovery method (trained by purely minimizing the constraints) could be useful to help the agent recover (Achiam et al., 2017; Chow et al., 2019) . 

A REPRODUCIBILITY CHECKLIST:

We follow the reproducibility checklist (Pineau, 2018) and point to relevant sections explaining them here. For all algorithms presented, check if you include:

• A clear description of the algorithm. The algorithms are explained in Sec. J. Any additional details for Discrete methods are provided in Sec. 4.1, and for continuous Sec. 4.2.

• An analysis of the complexity (time, space, sample size) of the algorithm. The analytical solution in Eq. ( 5) consists of a few vector arithmetic and relu operator and as such has the same complexity as the baselines. For the discrete case, with deterministic policies the solution again can be implemented as part of the computation graph, consisting of basic vector arithmetic operations, and has very little additional overhead. For discrete actions with stochastic policies, one needs to sovle the LP problem in (SPI). In that case the complexity is same as the baseline safe-methods (Lyapunov), and is higher than the unconstrained versions. In terms of computation time (for Deep-RL experiments) the newly proposed algorithms are almost identical to the baselines due to its parallelizable nature. We do not make any claims about the sample complexity.

• A link to a downloadable source code, including all dependencies. The code will be made available after the acceptance of the paper.

For any theoretical claim, check if you include:

• A statement of the result. See the main paper for all the claims we make. Additional details are provided in the Appendix.

• A clear explanation of any assumptions. See the main paper for all the assumptions.

• A complete proof of the claim. See the main paper. The cross references to the proofs in the Appendix have been included in the main paper.

For all figures and tables that present empirical results, check if you include:

• A complete description of the data collection process, including sample size. For the base agent we standard benchmarks provided in OpenAI Gym (Brockman et al., 2016) , and rllab (Duan et al., 2016) . We use the code from Achiam et al. (2017) for building the Point-Circle and Point-Gather environments.

• A link to downloadable version of the dataset or simulation environment. See: github.com/openai/gym for OpenAI Gym benchmarks, github.com/jachiam/cpo for rllab based Circle and Gather environments.

• An explanation of how samples were allocated for training / validation / testing. We do not use a split as we run multiple runs over random seeds to examine the optimization performance.

• An explanation of any data that were excluded. NA

• The range of hyper-parameters considered, method to select the best hyper-parameter configuration, and specification of all hyper-parameters used to generate results. The default hyper-parameters for the MuJoCo baselines are taken from Kostrikov (2018) . The ranges and parameters for Grid experiments are described in Sec. H, and for MuJoCo are described in Sec. I.

• The exact number of evaluation runs. The number of evaluation runs is mentioned in the caption corresponding to each result.

• A description of how experiments were run. See Experiments Sec. 6 in the main paper and in the Appendix Sec. H and Sec. I.

• A clear definition of the specific measure or statistics used to report results. Undiscounted return and cost using the current policy over the horizon are plotted after every 1000 episodes are plotted. We use a linear-filter with 0.7 weight for smoothing. We use the smooting algorithm provided by TensorBoard (https://github.com/tensorflow/ tensorboard).

• Clearly defined error bars. Standard error used in all cases.

• A description of results with central tendency (e.g. mean) and variation (e.g. stddev).

The bold lines in the figure represent the mean, and the shaded region denotes the 80% confidence interval.

• A description of the computing infrastructure used. We distribute all runs across 10 CPU nodes (Intel(R) Xeon(R) CPU E5-2650 v4) and 1 GPU (GP 100) per run for experiments.

B BACKWARD VALUE FUNCTIONS:

We have the following result from Proposition 1 from Morimura et al. (2010) . We give the proof too for the sake of completeness.

Proposition B.1. Let the forward Markov chain M(π) be irreducible and ergodic, i.e., has a stationary distribution. Then the associated backward Markov chain # » B (π) is also ergodic and has the same unique stationary distribution as M(π):

η π (x) = # » η π (x), (∀x ∈ X )

where η π (x) and # » η π (x) are the stationary distributions of M(π) and # » B (π).

Proof. Multiply both sides of Eq. ( 1) by η π (x t ) and sum over all actions a t-1 ∈ A we obtain detailed balance like equations (with respect to time):

# » P π (x t-1 |x t )η π (x t ) = P π (x t , a t-1 |x t-1 )η π (x t-1 ). (∀x t-1 ∈ X , x t ∈ X )

Sum over all possible x t we have:

xt∈X # » P π (x t-1 |x t )η π (x t ) = η π (x t-1 ).

The above equation indicates that # » B (π) has same stationary distribution as M(π). In the matrix form the above equation can be written as η # » P π = η, that implies that η is stationary distribution with # » P π transition matrix.

B.1 RELATION BETWEEN FORWARD AND BACKWARD MARKOV CHAINS AND BACKWARD VALUE FUNCTIONS:

Proof. We use the technique of Proposition 2 of Morimura et al. (2010) to prove this. Using the Markov property and then substituting Eq. ( 1) for each term we have:

# » P π (x t-1 , a t-1 , . . . , x t-K , a t-K |x t ) = # » P π (x t-1 , a t-1 |x t ) . . . # » P π (x t-K , a t-K |x t-K+1 ), = P π (x t , a t-1 |x t-1 ) . . . P π (x t-K+1 , a t-K |x t-K )η π (x t-K ) η π (x t ) , ∝ P π (x t , a t-1 |x t-1 ) . . . P π (x t-K+1 , a t-K |x t-K )η π (x t-K ).

This proves the proposition for finite K. Using the Prop. B.1, K → ∞ case is proven too:

lim K→∞ E # » B (π) K k=0 d(x t-k )|x t = lim K→∞ E M(π) K k=0 d(x t-k )|x t , η π (x t-K ) , = x∈X a∈A π(a|x)η π (x)d(x).

B.2 TD FOR BVF:

Proof. We use the same technique from Stochastic Shortest Path dynamic programming (Bertsekas et al., 1995, Vol 2, Proposition 1.1) to prove the above proposition. The general outline of the proof is given below, for more details we refer the reader to the textbook.

We have,

# » T π # » V = d + # » P π # » V .

(Eq. ( 4) in matrix notation)

Using induction argument, we have for all # » V ∈ R n and k ≥ 1, we have:

# » T π k # » V = # » P π k # » V + k-1 m=0 # » P π m d,

Taking the limit, and using the result, 1.2), we have:

lim k→∞ # » P π k # » V = 0,

lim k→∞ # » T π k # » V = lim k→∞ k-1 m=0 # » P π m d = # » V π ,

Also we have by definition:

# » T π k+1 # » V = d + # » P π # » T π k # » V ,

and by taking the limit k → ∞, we have:

# » V π = d + # » P π # » V π , which is equivalent to, # » V π = # » T π # » V π .

To show uniqueness, note that if

# » V = # » T π # » V , then # » V = # » T π k # » V for all k and letting k → ∞ we get # » V = # » V π . C VALUE-BASED CONSTRAINT LEMMA Lemma C.1. E T k=t d(x k ) | x 0 , π ≤ E xt∼δx 0 (P π ) t V π D (x t ) and E t k=0 d(x k ) | x 0 , π ≤ E x k ∼δx 0 (P π ) t # » V π D (x k )

Proof. Follows since adding more steps to the trajectory (from T -t steps to T ) can only increase the expected total cost.

E T k=t d(x k ) | x 0 , π = δ x0 (P π ) t T k=t (P π ) k d ≤ δ x0 (P π ) t T +t k=t (P π ) k d = E xt∼δx 0 (P π ) t V π D (x t )

. The backward case is analogous.

D PROPERTIES OF THE POLICY ITERATION (SPI):

Theorem D.1. Let σ(x) := TV(π k+1 (•|x), π k (•|x)) = (1/2) a |π k+1 (a|x) -π k (a|x)

| denote the total variation between policies π k (•|x) and π k+1 (•|x). If the policies are updated sufficiently slowly and π k is feasible, then so is π k+1 . More specifically:

(I) If π k is feasible at x 0 and σ(x) ≤ d0-V π k D (x0)

2T 2 DMAX ∀x then π k+1 is feasible at x 0 .

(II) If π k is feasible everywhere (i.e.

V π k D (x) ≤ d 0 ∀x) and σ(x) ≤ d0-V π k D (x) 2T max x {d0- # » V π k D (x )-d(x )}

∀x then π k+1 is feasible everywhere.

We note that the second case allows the policies to be updated in a larger neighborhood but requires π k to be feasible everywhere. By contrast the first item updates policies in a smaller neighbourhood but only requires feasibility at the starting state.

Proof. Similar to the analysis in Chow et al. (2018) . We aim to show that V π k+1 D (x 0 ) ≤ d 0 . For simplicity we consider k = 0, and by induction the other cases will follow. We write P 0 = P π0 , P 1 = P π1 , ∆(a|x) = π 1 (a|x) -π 0 (a|x), and P ∆ = a∈A ∆(a|x)P (x |x, a) {x ,x} . Note that (I -P 0 ) = (I -P 1 + P ∆ ), and therefore (I -P 1 + P ∆ )(I -P 0 ) -1 = I |X |×|X | . Thus, we find

(I -P 0 ) -1 = (I -P 1 ) -1 (I |X |×|X | + P ∆ (I -P 0 ) -1 ).

Multiplying both sides by the cost vector d one has

V π0 D (x) = E T t=0 d(x t ) + ε(x t ) | π 1 , x , for each x, where ε(x) = a∈A ∆(a|x) x ∈X P (x |x, a)V π0 D (x ).

Splitting the expectation, we have

V π1 D (x) = V π0 D (x) -E T t=0 ε(x t ) | π 1 , x

For case (I) we note that V π0 D (x ) ≤ D MAX T and so -2σ

(x t )D MAX T ≤ ε(x t ) ∀x t . Using σ(x t ) ≤ (d 0 -V π k D )/2D MAX T 2 gives V π1 D (x 0 ) ≤ V π0 D (x 0 )-2D MAX T 2 (d 0 -V π0 D (x 0 ))/(2D MAX T 2 ) = d 0 , i.e. π 0 is feasible at x 0 . For case (II) we note that V π0 D (x) ≤ max x {d 0 - # » V π0 D (x ) -d(x )} =:

Θ since π 0 is feasible at every x. As before, we have -2σ(x t )Θ ≤ ε(x t ) ∀x t and so V π1 D (x) ≤ V π0 D (x) -2ΘT (d 0 -V π0 D (x))/(2ΘT ) = d 0 ∀x, i.e. π 1 is feasible everywhere. Theorem D.2. Let π n and π n+1 be successive policies generated be the policy iteration algorithm of (SPI). Then V πn+1 ≥ V πn .

Proof. Note that π n+1 and π n are both feasible solutions of the LP (SPI). Since π n+1 maximizes V π over all feasible solutions, the result follows.

E ANALYTICAL SOLUTION OF THE UPDATE -DISCRETE CASE:

We follow the same procedure as (Chow et al., 2018, Section E.1) to convert the problem to its Shannon entropy regularized version:

EQUATION

where τ > 0 is a regularization constant. Consider the Lagrangian problem for optimization:

max λ≥0 max π∈∆ Γ x (π, λ) = π(.|x) T (Q(x, .) + λQ D (x, .) + τ log π(.|x)) + λ(d 0 + d(x) - # » V (x))

From entropy-regularized literature (Neu et al., 2017) , the inner λ-solution policy has the form:

π * Γ,λ (.|x) ∝ exp - Q(x, .) + λQ D (x, .) τ

We now need to solve for the optimal lagrange multiplier λ * at x.

max λ≥0 -τ log-sum-exp - Q(x, .) + λQ D (x, .) τ + λ(d 0 + d(x) - # » V D (x)),

where log-sum-exp(y) = log a exp(y a ) is a convex function in y, and objective is a concave function of λ. Using KKT conditions, the ∇ λ gives the solution:

(d 0 + d(x) - # » V D (x)) - a Q D (x, a) exp( -Q(x,a)+λQ D (x,a) τ ) a exp( -Q(x,a)+λQ D (x,a) τ ) = 0

Using parameterization of z = exp(-λ), the above condition can be written as polynomial equation in z:

a d 0 + d(x) - # » V D (x) -Q D (x, a) . exp(- Q(x, a) τ ) z Q D (x,a) τ = 0

The roots to this polynomial will give 0 ≤ z * (x) ≤ 1, using which one can find λ * (x) = -log(z * (x)). The roots can be found using the Newton's method. The final optimal policy of the entropy-regularized process is then:

π * Γ ∝ exp - Q(x, •) + λ * Q D (x, •) τ

F EXTENSION OF SAFETY LAYER TO STOCHASTIC POLICIES WITH GAUSSIAN PARAMTERIZATION:

Consider stochastic gaussian policies parameterized by mean µ(x; θ) and standard-deviation σ(x; φ), and the actions sampled have the form µ(x; θ) + σ(x; φ) , where ∼ N (0, I) is the noise. Here, < µ(x; θ), σ(x; φ) > are both deterministic w.r.t. the parameters θ, φ and x, and as such both of them together can be treated in the same way as deterministic policy (π(x) =< µ(x), σ(x) >). The actual action sampled and executed in the environment is still stochastic, but we have moved the stochasticity fron the policy to the environment. This allows us to define and work with action-value functions Q D (x, µ π (x), σ π (x)). In this case, the corresponding projected actions have the form µ + σ . The main objective of the safety layer (without the constraints) can be further simplified as:

arg min µ ,σ E ∼N (0,I) 1 2 (µ + σ ) -(µ π (x) + σ π (x) ) 2 arg min µ ,σ E ∼N (0,I) 1 2 (µ -µ π (x)) + ((σ -σ π (x)) ) 2 arg min µ ,σ 1 2 E ∼N (0,I)    µ -µ π (x) 2 + (σ -σ π (x)) 2 + 2 < µ -µ π (x), (σ -σ π (x)) > =0,due to linearity of expectation, ∼N (0,I)    arg min µ ,σ 1 2 µ -µ π (x) 2 + E ∼N (0,I) (σ -σ π (x)) 2 arg min µ ,σ 1 2     µ -µ π (x) 2 + (σ -σ π (x)) 2 E ∼N (0,I) 2 =1,second moment of     arg min µ ,σ 1 2 µ -µ π (x) 2 + (σ -σ π (x)) 2

As both µ π (.; θ) and σ π (.; φ) are modelled by independent set of parameters (different neural networks, usually) we can solve each of the safety layer problem independently, w.r.t. only those parameters. • Point-Gather: The environment (Fig. 4c ) is taken from Achiam et al. (2017) , where the point mass agent gets a reward of +10.0 for collecting a green apple, and a cost of 1 for collecting a red bomb. Two apples and eight bombs are spawned randomly at the start of each episode. The constraints are defined over the nmber of bombs collected over the episode. Episode horizon is 15 and threshold d 0 = 4.

• Safe-Cheetah: This environment (Fig. 4b ) is taken from Chow et al. (2019) . A bi-pedal agent (HalfCheetah-v0) is augmented with speed safety constraints. The agent gets the reward based on the speed with which it runs, and the constrain is define on the speed to be less than 1, i.e., it gets a constraint cost based on 1[|v| > 1], where v is the velocity at the state. The maximum length of the episode is 200 and the constraint threshold is d0 = 50.

• Point-Circle: This environment (Fig. 4a ) is taken from Achiam et al. (2017) . The pointmass agent is rewarded for running along the circumference of a circle of radius 15 in counter-clockwise direction, with the reward and cost function:

R(s) = v T [-y, x] 1 + | [x, y] 2 -15| , C(s) = 1[|x| > 2.5],

where x, y are coordinates in the plane and v is the velocity. The length of the episode is 65 and the constraint threshold d 0 = 10.0.

I.2 NETWORK ARCHITECTURE AND TRAINING DETAILS:

The architecture and the training procedure is based on the open-source implementations (Kostrikov, 2018) . All the value based estimators use a network architecture of 2 hidden layers of size 200, 50 hidden units with tanh non-linearity, followed by a linear layer with single output. For the actor, we model mean using a network architecture of 2 hidden layers of size 100, 50 hidden units with tanh non-linearity, followed by a linear layer with dimensions of the action-space and tanh non-linearity. For the Q(x, µ) we also a 2 layer neural network with 200, (50 + action-dimension) hidden units and tanh non-linearity. We concatenate the mean in the second layer, and add a linear layer with single output in the end.

Entropy regularization with β = 0.001 was used for all the experiments and the baselines. The trajectory length for different environments. For PPO GAE with λ = 0.95 was used for every algorithm. 20 parallel actors were used for every algorithm for each experiment. We searched the trajectory length hyper-parameter in the range 5,20,100 for every environment. We used the trajectory length of 1000 over which the samples are collected for PPO, for all environments. For the A2C experiments, for SafeCheetah trajectory length of 5 is used and for the rest 20 is used.

We use Adam Optimizer for training all the estimators. The learning rate of the critic is always 0.5 the learning rate of the actor. For the cost estimators, the same learning rate was used for forward and backward estimators. The same range of learning rate parameters for considered for all estimators i.e. {1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1}.

I.3 OTHER DETAILS:

As we mentioned in Sec. 7, due to exploration the agent can potentially end up being in an infeasible policy space. To prevent that from happening a recovery policy (or safe-guard policy) (Achiam et al., 2017; Chow et al., 2019) is used to recover back to the feasible policy space. We run the experiments with and without the use of recovery policies (in the same procedure as the baselines), and chose the run that performs the best. We noticed that, empirically, for our approach recovery policies are only required for Point-Circle environments, as the agent has much more probability of being stuck in the constraint space.

In order to take error due to function approximation into account, Achiam et al. (2017) use costshaping to smooth out the sparse constraint, and Chow et al. ( 2019) use a relaxed threshold, i.e. d 0 • (1 -δ), instead of d 0 , where δ ∈ (0, 1). We run experiments with δ = {0.0, 0.2} for each algorithms, and use the best among them. We found that empirically, only for Safe-Cheetah δ = 0.2 works better compared to δ = 0.0.

J ALGORITHM DETAILS:

J.1 N-STEP SYNCHRONOUS SARSA

The algorithm for n-step Synchronous SARSA is similar to the n-step Asynchronous Q-learning of Mnih et al. (2016) , except that it uses SARSA instead of Q-learning, is synchronous, and instead of greedy maximization step of -greedy we use (SPI). When working with discrete actions and deterministic policies, this can be solved as part of the computation-graph itself. The algorithm is presented in Alg. 1.

J.2 A2C:

In Actor Critic (Konda & Tsitsiklis, 2000) algorithms, the parameterized policy (actor) is denoted by π(a|x; θ), and is updated to minimizing the following loss:

L(θ) = E[-log π(a t |x t ; θ)(r t + γV π (x t+1 -V xt ))]

The algorithm for A2C with Safety Layer given by Eq. ( 5) is similar to the Synchronous version of Actor-Critic (Mnih et al., 2016) , except that it has estimates for the costs and safety-layer. Note that due to the projection property of the safety layer, it is possible to sample directly from the projected mean. Also, as the projection is a result of vector products and max, it is differentiable and and computed in-graph (via relu). The algorithm is presented in Alg. 2.

J.3 PPO:

The PPO algorithm build on top of the Actor-Critic algorithm and is very similar to Algorithm 2. The main difference is how the PPO loss for the actor is defined as:

L CLIP (θ) = E[min(ρ t (θ)A t , clip(ρ t (θ), 1 -, 1 + )A t )],

where the likelihood ration is ρ t (θ) = π θ (at|xt) π θ old (at|xt) , with π old being the policy parameters before the update, < 1 is a hyper-parameters that controls the clipping and A t is the generalized advantage estimator: 

A GAE(λ,γ) t = T -1 k=0 (λγ) k δ V π t+k ,

dφ D ← dφ D + ∂( # » R D - # » V D (x i ; φ D )) 2 ∂φ D

end for Do synchronous batch update with the accumulated gradients to update θ, θ D , φ D using dθ, dθ D , dφ D .

end while Empty the trajectory buffer, τ end for where T is the maxmimum number of timestamps in an episode trajectory, and δ j denotes the TD error at j. The value function is updated using the γλ-returns from the GAE:

L(φ) = E[(V π (x; φ) -(V π (x; φ old ) + A t )) 2 ].

Similar to the the forward value estimates the backward value estimates are defined in the similar sense. One way to think of it is to assume the trajectories are reversed and we are doing the regular GAE estimation for the value functions.

The GAE updates for the regular value function can be seen in the λ-operator form as:

T π λ v π = (I -γλP π ) -1 (r π + γP π v π -v π ) + v π .

In similar spirit it can be shown that the λ-operator for SARSA has the form:

T π λ q π = (I -λγP π ) -1 (T π q π -q π ) + q π , where (T π q π -q π ) denotes the TD error. Thus, the GAE estimates can be applied for the Q-functions in the similar form, i.e. 

:

Here the cost only depends on states and not state-action pairs.

δx 0 is a Dirac distribution at x0, and δx 0 (P π ) t is the distribution of states at time t.

In general, it is not possible to obtain the expectation E x t ∼δx 0 (P π ) t [•] directly as it may be intractable to compute or we may not have access to the true transition distributions of the model. Thus, we sample a batch of transitions from the current policy and use them for the updates.

This can be enforced, for example, by constraining iterates to a neighborhood D(π, π k ) ≤ δ.

More information about this claim can be found in Appendix F.

Classic Actor Critic is two time-scale(Konda & Tsitsiklis, 2000), and adding a learning schedule over the Lagrangian makes it three time scale.

In practice, the starting states in the episode are unlikely to be distributed according to the stationary distribution η π . We still use the initial trajectories to update the estimates nonetheless in our experiments, but we use n-step updates.

Figures and Tables:

Figure 1: (a) Example of a gridworld environment. (b,c) Performance over the training for Unconstrained (red), Lyapunov-based (green), and our method (blue) all trained with n-step SARSA on 2D GridWorld task over 20 random seeds. The x-axis is the number of episodes in thousands. The dotted black line in (c) denotes the constraint threshold, d0. The bold line represents mean, and the shaded region denotes 80% confidence-intervals.

Figure 2: A2C Performance over the training for Unconstrained (red), Lyapunov-based (green), and our method (blue) all trained with A2C on MuJoCo tasks over 10 random seeds. The x-axis is the number of episodes in thousands. The dotted black line denotes d0. The bold line represents the mean, and the shaded region denotes the 80% confidence-intervals.

regarding proper policies fromBertsekas et al. (1995, Vol 2, Equation

Figure 4: MuJoCo Safety Environments

D ) = E[(Q π D (x, a; θ D ) -(Q θ D D (x, a; θ D old ) + B t )) 2 ].

-1 , a t-1 |x t ), and describes the state and action the agent took to reach the current state.

Table content: <table><tr><td>Definition 3.1 (Backward Markov Chain). A backward Markov chain associated with M(π) is denoted by # » B (π) and is characterized by the transition probability # » P</td></tr></table>

Algorithm 1 Synchronous n-step SARSA Input: θ parameters for Q(x, .; θ), θ D parameters for Q D (x, .; θ D ),φ D parameters for # » V D (x; φ D ); π 0 initial feasible policy. for episode e ∈ 1, ..., M do Add the initial state to the trajectory buffer τ ← {x 0 } t ← 1 while t < T do: t start ← t while t < t + n or t == T do Select a t using (SPI), execute a t , observe x t+1 and reward r t and cost d t . Add experiences to a buffer, i.e., τ ← (a t , r t , d t , x t+1 ). t ← t + 1 end while Calculate the next action for x t+1 using the current policy estimates, a t+1 Bootstrap the targets: t+1 , a t+1; θ) otherwise R D ← 0 if t == T Q D (x t+1 , a t+1 ; θ D ) otherwise tstart-1 ; φ D ) otherwiseCalculate the targets for the transitions in bufferfor i ∈ {t -1, . . . , t start } do R ← r i + γR R D ← d i + γR DAccumulate the gradients wrt θ, θ D :dθ ← dθ + ∂(R -Q(x i , a i ; θ)) 2 ∂θ dθ D ← dθ D + ∂(R D -Q D (x i , a i ; θ D ))2 ∂θ D end for for i ∈ {t start , . . . , t} do # » R D ← d i + γ # » R D Accumulate the gradients wrt φ D :

Table content: <table><tr><td>0 Q(x # R ← » R D ← 0 # » V (x</td><td>if t == T if t == 0</td></tr></table>

Peer review (iclr2020):
This paper presents a new approach for solving Constrained MDPs. Because the cost constraint is cumulative, the best action depends on the cumulative cost so far. They address this issue by learning a backward value function of the estimated cumulative cost so far. Their theoretical results show that the same properties for forward value functions hold here for backwards ones. They are then able to use the forward and backward cost estimates to constraint the actions selection, by adding a safety layer to the algorithm. The results show that the method does a better job of meeting safety constraints than the Lyapunov based method. The backward value function idea is a nice novel way of addressing the cumulative cost constraint problem. The paper is clearly written and the results are nice. The biggest issue with this paper is that too much material has been pushed to the appendix. I think at the least some of the details of the actual implementation with PPO or A2C should be moved into the main text. For the empirical results, it would be great to see something about the computation or wall clock time required. Does this method run more quickly than the Lyapunov based method?