Safe Exploration and Optimization of Constrained MDPs
using Gaussian Processes
Akifumi Wachi

Yanan Sui and Yisong Yue

Masahiro Ono

University of Tokyowachi@space.t.u-tokyo.ac.jp

Calfornia Institute of Technology{ysui, yyue}@caltech.edu

Jet Propulsion Laboratory,
California Institute of Technology
ono@jpl.nasa.gov

Abstract
We present a reinforcement learning approach to explore
and optimize a safety-constrained Markov Decision Process(MDP). In this setting, the agent must maximize discounted
cumulative reward while constraining the probability of entering unsafe states, defined using a safety function being
within some tolerance. The safety values of all states are notknown a priori, and we probabilistically model them via a
Gaussian Process (GP) prior. As such, properly behaving insuch an environment requires balancing a three-way trade-offof exploring the safety function, exploring the reward function, and exploiting acquired knowledge to maximize reward.
We propose a novel approach to balance this trade-off. Specifically, our approach explores unvisited states selectively; thatis, it prioritizes the exploration of a state if visiting that statesignificantly improves the knowledge on the achievable cumulative reward. Our approach relies on a novel information
gain criterion based on Gaussian Process representations of
the reward and safety functions. We demonstrate the effectiveness of our approach on a range of experiments, including
a simulation using the real Martian terrain data.

Introduction
In many environments, an autonomous agent must discoverboth potential hazards as well as rewards on-the-fly. For
example, in case of NASA’s rovers exploring the surfaceof Mars, the driving safety and scientific gain of unvisitedplaces cannot be perfectly known a priori. Hence, at the endof each Martian day, called Sol (approximately 24 hours and
40 minutes), the rover takes a set of panorama images andsend it to the Earth. The ground operators identify sciencetargets as well as hazards in the images, and plan a path forthe next Sol, which is typically less than 100 meters. This
process is exploratory by nature; the rover needs to explore
unvisited regions in order to make discoveries, but in a way
to assure safety. However, the agent is not fully autonomousin this case because there does not exist a good algorithm tofully drive the decision-making process in situ.
Traditional exploration approaches in reinforcement
learning (Sutton and Barto 1998; Bertsekas and Tsitsiklis
1995) are agnostic to safety, and instead focus on efficient
exploration under the assumption of ergodicity of the state
Copyright c 2018, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

transitions, i.e., that all states are assumed to be reachable.
In practical applications, an agent is not necessarily able
to reach all (safe) states because of (blocking) hazardous
states. On the other hand, especially in the field of robotics,
there is a significant body of work on safety-constrained
decision-making (Fleming and McEneaney 1995; Schwarm
and Nikolaou 1999; Blackmore et al. 2010; Ono et al. 2015).
While these methods consider uncertainty what state one
is in, they assume that unsafe states are known a priori.
Removing this assumption requires exploration of safety,
which is a central focus of our work.
To develop a tractable approach for exploring the safe
region, we make regularity assumptions that similar states
have similar safety levels. In particular, we model the safetyfunction as a Gaussian process (GP) (Rasmussen 2006), with
an appropriate kernel to capture similarity between states.
Previous work on safe optimization using GP have largely
focused on the state-less setting, such as in Bayesian optimization (Mockus 2012). Another typical limitation is when
the reward function and the safety function are identical (Suiet al. 2015), e.g., an unsafe state is one where the reward
function is too low.
In stateful settings, such as Markov Decision Processes(MDPs), previous work on safe optimization either ignore additional structure (Moldovan and Abbeel 2012), or
can only exploit structure for exploring the safe region(Turchetta, Berkenkamp, and Krause 2016) (i.e., pure exploration). The former suffers from longer convergence times
since it cannot exploit the GP structure, and the latter suffersfrom not actually optimizing the reward function.
Our Contributions. In this paper, we propose a novelsafe-constrained exploration and optimization approach thatmaximizes discounted cumulative reward while guaranteeing safety. As demonstrated in Figure 1, we optimize overconstrained MDPs with a priori unknown two functions, one
for reward and the other for safety. A state is considered safe
if the safety function value is above a threshold.
As alluded to above, compared to previous work on
jointly exploring and optimizing safety-constrained MDPs
(Moldovan and Abbeel 2012), our approach is unique in two
aspects. First, (Moldovan and Abbeel 2012) puts no assumption on the structure of reward and safety functions. Whilethat approach is more general, in practice, nearby states often have the similar levels of safety. For example, if a Mars

Figure 1: Conceptual images of the planning problem under uncertainty on reward and safety.
rover experiences high slippage at a certain location, it islikely to experience similar slip one meter away. Properly
modeling this dependency can significantly reduce the costof exploration. Second, our approach is able to estimate thelong-term value of exploring new safe regions. In contrast,
the approach by (Moldovan and Abbeel 2012) assigns uniform additional reward for visiting unvisited states, which
can be substantially less efficient.
At a high level, our approach employs two MDPs, which
we call Optimistic and Pessimistic MDPs, and uses the difference in the value functions as the information gain criterion. Our GP safety function yields three classes of states:
safe, unsafe, and uncertain. The only difference between the
Optimistic and Pessimistic MDPs is that uncertain states are
considered safe former and unsafe in the latter. Using thiscriterion, the agent is motivated to explore uncertain statesthat could result in high cumulative reward if they are determined safe.
We demonstrate empirically that our approach yields
more efficient safe exploration and optimization compared to (Moldovan and Abbeel 2012) and (Turchetta,
Berkenkamp, and Krause 2016). The simulations are performed on both synthetic environments, as well as real environments of two locations on Mars.

Problem Statement
We formulate a constrained MDP, defined as a tuple hS, A,
f (s, a), r(s, a), g(s), γi. S is a set of states {s}, A is a setof actions {a}, f : S × A → S is a deterministic transition
model, r : S × A → R is a (bounded) reward function, g :S → R is a safety function, and γ ∈ R is a discount factor.
We assume that reward and safety functions are not known apriori. Hence, an agent must explore the state space to learn
these functions. At every discrete time step t, the agent isrequired to be in a safe state, that is, the safety function valueg(st ) of state st must be above some threshold h ∈ R. A
policy π : S → A maps a state to actions. The value of apolicy is defined as the discounted cumulative reward overthe next N steps while guaranteeing the satisfaction of thesafety constraint. The optimal planning problem in this MDP
is represented as follows:
" N
#
X
π
t−1
maximize : VN (s) = E
γ r(st , at )t=1

subject to : g(st ) ≥ h, ∀t = [1, N ].

In conventional reinforcement learning, the exploration/exploitation trade-off is a central issue. Compared tothe conventional setting, the key difference of our setting is
the need to also explore the safety function g. In the rest of
the paper, we first describe how we use Gaussian processes
to model the safety function and classify the states as safe,
unsafe, or uncertain. Then, we describe our algorithm forleveraging this information to efficiently explore the safe
region to the extent necessary for optimizing the reward
function.

Characterization of Safety
We employ Gaussian processes (GPs) to represent and estimate the safety and reward functions, g and r. For clarityof exposition, we focus our description of GPs with respect
to modeling g.1 In many practical applications, safety exhibits regularity; that is, similar states possess the similarlevel of safety.2 Such regularities can be naturally modeledusing GPs.
Using the GP model, the value of g at unvisited states
are predicted based on the observations at visited states. The
GP-based predictions have uncertainty, which is represented
by Gaussian distributions. As a result, unvisited states arenaturally categorized into three classes: safe, unsafe, and uncertain. Roughly speaking, an unvisited state is consideredsafe if the safety function value is above the threshold almost
certainly (i.e., with a greater probability than a pre-defined
confidence level), is considered unsafe if the safety function
value is below a threshold almost certainly, and is considered
uncertain otherwise. These concepts are shown in Figure 2.
We assume that state space S is endowed with a positive
semidefinite kernel function k(·, ·), and that the safety function g(s) has bounded norm in the associated Reproducing
Kernel Hilbert Space (RKHS). The kernel function k is used
to capture similarity between states, and formalizes the assumption that similar states have similar levels of safety. The
safety function can thus be modeled as:g(s) = GP(µg (s), k g (s, s0 )).
A GP is fully specified by its mean, µg (s), and covariance, k g (s, s0 ). Without loss of generality, let µg (s) = 0
1

Our approach does not prohibit other regression methods frombeing used for the reward function.
2
We do not focus on the situations where environmental hazards
are expressed as binary function (e.g., rocks and cliffs).

(a) t = 1(b) t = 2(c) t = 8

Figure 2: Conceptual image of the characterization of safety. Green curve and blue band represent mean and confidence interval,
respectively. The magenta line means the safety threshold, h. Blue, yellow, and red bars represent the safe region, uncertain
region, and unsafe region, respectively.
for all s ∈ S. We model observation noise as y g =g(s) + ng , where ng ∼ N (0, ω 2 ). The posterior over g(s)is computed analytically, based on t measurements at states
At = {s1 , · · · , st } with measurements, ytg = [g(s1 ) +ng1 , · · · , g(st )+ngt ]> . The posterior mean, variance, and covariance are given as:µgt (s) = ktg (s)> (Ktg + ω 2 I)−1 ytg
σtg (s) = ktg (s, s)(1)

ktg (s, s0 ) = k g (s, s0 ) − ktg (s)> (Ktg + ω 2 I)−1 ktg (s0 ),where ktg (s) = [k g (s1 , s) · · · , k g (st , s)]> and Ktg is the
positive semidefinite kernel matrix, [k g (s, s0 )]s,s0 ∈At . We
further assume that the safety function g is Lipschitz continuous, with Lipschitz constant L, with respect to some metric d(·, ·) on S. This assumption is naturally satisfied using
common kernel functions (Ghosal and Roy 2006).
Safe space We employ the same definition of safe space as(Turchetta, Berkenkamp, and Krause 2016). First, let S̃tsafe
be the set of states that satisfy the safety constraint. The safe
space is a subset of S̃tsafe that is reachable and returnable.
Intuitively, a state is reachable if it can be reached from avisited safe state. For example, an “island” of unvisited safe
states that is isolated by unsafe states is not reachable. A
state is returnable if the agent can return to a visited safe
state from the state. For example, on a directed graph, a basinof unvisited safe states that “trap” the agent because there areno outgoing edges is not returnable.
First, we formally define S̃tsafe , the set of states that satisfy
the safety constraint. A GP model enables an agent to judge
safety by providing it the√ confidence interval that has a form
Qt (s) = [µt−1 (s) ± βt σt−1 (s)], where βt is a positive
scalar specifying the required level of safety. In other words,
β inherently specifies the probability of violating the safety
constraint. Following (Srinivas et al. 2010) and (Sui et al.
2015), for our experiments we use βt = 2, ∀t ≥ 0, andrefer to the previous work for a more detailed discussion. We
then consider the intersection of Qt up to iteration t, whichis recursively defined as Ct (s) = Qt (s)∩Ct−1 (s), C0 (s) =[h, ∞], ∀s ∈ S̃0safe . We denote the lower and upper bounds
on Ct (s) by lt (s) := min Ct (s) and ut (s) := max Ct (s),respectively. With the Lipschitz constant L, S̃tsafe is defined
as follows:safe
S̃tsafe = {s ∈ S | ∃s0 ∈ St−1
: lt (s0 ) − Ld(s, s0 ) ≥ h}.

Since the prediction by the GP is probabilistic, being in S̃tsafe
means the safety constraint is satisfied with high probability.
The probability of entering the safe states is guaranteed to be
Pr[g(s) > h | l(s0 ) − Ld(s, s0 ) > h]
≥ Pr[g(s) > h | l(s) > h]
≥ Pr[g(s) > l(s) | l(s) > h]
= Pr[g(s) > l(s)]The last equality exploits the conditional independence.
Note that this probability can be tuned by β.
Next, we formally define the reachable and returnable
sets. Given S, the set that is reachable from S in one step is
Rreach (S) = S ∪ {s ∈ S | ∃s0 ∈ S, a ∈ A : s = f (s0 , a)}.
Given sets S and S̄, the set that is returnable to S̄ with
one step is given as Rret (S, S̄) = S̄ ∪ {s ∈ S | ∃a ∈A : f (s, a) ∈ S̄}. The n-step returnability operator canbe represented as Rnret (S, S̄) = Rret (S, Rret (S, · · · )). Using this operator, the returnable set is given by R̄ret (S, S̄) =limn→∞ Rnret (S, S̄). Finally, the set of safe states that fulfill
reachability and returnability, Stsafe , is defined as:
safe
safe
Stsafe = {s ∈ S̃tsafe | s ∈ Rreach (St−1) ∩ R̄ret (S̃tsafe , St−1)}.

Unsafe space We extend the setting of (Turchetta,
Berkenkamp, and Krause 2016) to also define the unsafe
space, Stunsafe . As in the safe set, we first define the set ofstates that violate the safety constraint with a high confidence as:S̃tunsafe = {s ∈ S | ∃s0 ∈ S : ut (s0 ) + Ld(s, s0 ) ≤ h}.
Using the returnable set defined above, the unsafe space, denoted by Stunsafe , is defined as:
safe
Stunsafe = {s ∈ S | s ∈ S̃tunsafe ∪ (S̃tsafe ∩ R̄cret (S̃tsafe , St−1)},

where Xc represents the complementary set of a set X. Notethat we do not take into account reachability. In other words,
an “island” of unvisited states should be regarded as unsafe.

Uncertain Space Finally, the uncertain space is simply theset of remaining states that belong to neither the safe nor
unsafe spaces:
Stuncertain = S \ Stsafe ∪ Stunsafe .
Intuitively, uncertain states can be either safe or unsafe, andan agent must explore the uncertain space to find new safe
states with high reward. However, uncertain states cannot be
guaranteed to be safe, so the agent must leverage the structure of the safety function (i.e., the similarity kernel) and
explore the boundary of the safe space to collect measurements that are informative of uncertain states.

MethodOur proposed algorithm is named S AFE E XP O PT-MDP ,which is summarized in Algorithm 1. At each iteration, theagent observes the reward and safety function values of thecurrent state, and update the GP models. Using the updated
GPs, the agent then updates the safe, unsafe, and uncertain
sets (Line 3). Given these sets, the agent solves Optimistic
and Pessimistic MDPs (explained in the following subsections) using a PAC-MDP or Near-Bayesian Optimal algorithm, and obtain their value functions (Line 4). Then, theagent chooses the next action such that the linear combination of the two value functions is maximized (Line 5). Notethat the linearly combined value functions in Line 5 can also
be represented as J¯N + η(JˆN − J¯N ).
Therefore, the resulting policy favors an action to visit
states that have a large difference in value between the Optimistic and Pessimistic MDPs. We interpret this difference
corresponds to an information gain criterion, because it becomes zero when the level of uncertainty on safety is similar along with the paths obtained in two MDPs. Hence, thismechanism motivates the agent to selectively explore thestates which could result in significant increase in the cumulative reward.

Value of Exploration of Safety
Ultimately, the agent’s objective is to maximize the amountof reward collected. Thus, the value of exploration should bedirectly tied to the potential reward that could be collected.
Then, how should we evaluate the value of exploration of agiven unvisited state? As alluded to above, we quantify thevalue of exploration as the difference between the optimistic
case (where uncertain states are all safe) and the pessimistic
case (where uncertain states are all unsafe), as described inLine 5 of Algorithm 1. The Bayes-optimal, optimistic, andpessimistic value functions are obtained by solving the following Bellman equation:hX
VN∗ (s, br , bg ) = max
I(s0 | s, a, bg ) · Ps,a,s0 ·as0i{E[r(s0 )] + γVN∗ −1 (s0 , brs,a,s0 , bg )} ,where I is a safety indicator function (explained shortly), P
is the deterministic transition probability, br and bg are thebelief of reward and safety functions, respectively, at time

Algorithm 1 S AFE E XP O PT-MDP
1: loop2:
Observe reward and safety values of the current state, andupdate GP models
3:
Update Stsafe , Stuncertain , and Stunsafe
4:Compute JˆN for Optimistic MDP and J¯N for Pessimistic
MDP by Eq. (5) and Eq. (6)5:Derive the optimal policy to maximize η JˆN + (1 − η)J¯N
by Eq. (8)6:Execute the next optimal action
7: end loop

step N , and brs,a,s0 is the posterior belief of reward after taking the action a and moving to s0 .
Optimistic MDP assumes that all the uncertain states
are safe. Hence, the safety indicator function for Optimistic
MDP is defined as:
1 if s0 ∈ S safe ∪ S uncertain
ˆ 0 | s, a, bg ) :=I(s(2)
0 otherwise
We denote the value function of the Optimistic MDP by V̂N .
Pessimistic MDP assumes that all the uncertain state areunsafe. Hence, the safety indicator function for Pessimistic
MDP is defined as:
1 if s0 ∈ S safe
¯ 0 | s, a, bg ) :=I(s(3)
0 otherwise
We denote the value function of the Pessimistic MDP by V̄N .
The two value functions always satisfy the relationship,∗V̄N ≤ VN∗ ≤ V̂N∗ ,where VN∗ is the value function of theBayes optimal policy for the original problem. Figure 3
shows an example of the classification of the state space,and Table 1 concludes the definition of value functions.

Incentivizing Exploration of Reward Function
Practically, the above Bellman equations for the Optimistic
and Pessimistic MDPs are very difficult to solve. Instead,
we solve them approximately. Previous work proposed algorithms that have a property called Probably Approximate
Correct in Markov Decision Process (PAC-MDP) (Kearns
and Singh 2002; Brafman and Tennenholtz 2002; Kakade
and others 2003; Strehl et al. 2006) or Near-Bayes Optimal (Strens 2000; Kolter and Ng 2009; Araya, Buffet, andThomas 2012). The general idea of PAC-MDP and NearBayesian Optimal is that, with probability of 1 − δ, an agent
following the policy π can obtain the reward such that
Vπ >V∗−for all but m time steps. For example, the Bayesian Exploration Bonus (BEB) algorithm (Kolter and Ng 2009) gives m
6that is O( |S||A|N
log |S||A|ε2δ ). In PAC-MDP and Near-Bayes
Optimal approaches, the approximately optimal policy π is
obtained by Delayed Q-learning or the exploration bonus algorithm, which adds “bonus” to the original reward function
for exploring unexplored states. As a result, the value function including the bonus, denoted by J π , is always greater

(b) Optimistic: 10%(c) Optimistic: 30%(d) Optimistic: 60%(e) Pessimistic: 10%(f) Pessimistic: 30%(g) Pessimistic: 60%(a) True classification

Figure 3: An example of the classification of state space. Blue region represents S safe ∪ S uncertain for optimistic case and S safefor pessimistic. The percentage for each figure means the time course against the time step large enough to explore S.
Table 1: Definition of value (objective) functions. As for reward, we can consider three cases: actually obtainable, Bayes
optimal, and value by algorithm. In the case of safety, we consider pessimistic, Bayes optimal, and optimistic cases.

RewardSafetyV̄NπActuallyobtainable

V̄N∗
Bayes
optimalπ
J¯N
Value byalgorithm

Pessimistic

Pessimistic

Pessimistic

than the actual value function and satisfies the following
with the probability of 1 − δ:
V

π

∗

≤V ≤J

π

and V

π

π

∗

≥ J −  ≥ V − ,(4)for all but m time steps. With an exploration bonus algorithm, such as BEB, the value functions with bonus for the
Optimistic and Pessimistic MDPs, Jˆπ and J¯π , are given as
follows:hX
Jˆπ (s, br , bg ) = maxIˆ · Ps,a,s0 ·
a
s0
i (5)
π
r g
ˆ
{E[r] + R̃ + γ J (s, b , b )} ,

VN∗
Bayes
optimalBayes
optimalV̂Nπactuallyobtainable

V̂N∗
Bayes
optimalπ
JˆN
Value byalgorithm

Optimistic

Optimistic

Optimistic

Overall Policy
As mentioned above, our algorithm evaluates the value of
exploration by the difference of the approximate value functions of the Optimistic and Pessimistic MDPs. Let ∆J be
the difference, and define∆JN := JˆN − J¯N .
Intuitively, the agent obtains a high information gain on
safety when visiting a state with a large ∆J. After computing JˆN and J¯N , the agent chooses the next action by following the policy given as:
πN = arg max {J¯N + η∆JN }(7)
s0 ∈S safe

= arg max {η JˆN + (1 − η)J¯N },J¯π (s, br , bg ) = max
a

hX

(8)

s0 ∈S safe

I¯ · Ps,a,s0 ·

s0

i (6)
π
r g
¯
{E[r] + R̃ + γ J (s, b , b )} ,where R̃ is the term for the exploration bonus. Note that the
belief of reward br and safety bg are not updated in theseequations, which means that we can solve the equation using the standard value iteration or policy iteration algorithm.
Our proposed method optimizes two MDPs, so has a similar
tractability compared with conventional MDP algorithms.

where η ∈ [0, 1] is a weight coefficient. Note that our approach optimizes the interpolating point between the optimistic and pessimistic policies. One can view η as controlling for the degree of optimism versus pessimism.3 Also,
note that the next state is chosen from S safe in Equation (8).

Proof of Consistency
We prove a consistency result that the objective functionin our approach converges as time horizon goes to infinity.
3

η often has intuitive, physical meaning, which helps in tuning.

Table 2: Simulation result for a simple grid world (100 times Monte-Carlo simulation). Our approach achieves higher expected
reward with less exploration, which suggests that our approach is able to better trade off between exploration and exploitation
and avoid unnecessary exploration.
S AFE E XP O PT-MDP
Moldovan & Abbeel (2012)
Turchetta et al. (2016)
Safe known
Safe/reward known

Normalized cumulative reward
0.774 ± 0.183
0.740 ± 0.213
0.701 ± 0.248
0.841 ± 0.089
1.00 ± 0.000

More specifically, the following theorem holds:
Theorem 1. For given δ, there exists the upper bound oferror, ε = η ε̂ + (1 − η)ε̄ for the value function (8) with aprobability of 1 − δ for all but the first m time steps suchthat:
0 ≤ (η JˆN + (1 − η)J¯N ) − (η V̂N∗ + (1 − η)V̄N∗ ) ≤ ε.
Proof. After the sufficiently large number of time steps, (4)holds with high probability for both the optimistic and pessimistic policies, and the following inequalities are satisfied:η V̂N∗ + (1 − η)V̄N∗ ≤ η JˆN + (1 − η)J¯N≤ η(V̂Nπ + ε̂) + (1 − η)(V̄Nπ + ε̄)≤ η(V̂N∗ + ε̂) + (1 − η)(V̄N∗ + ε̄),where ε̂ and ε̄ are confidence intervals in the optimistic andpessimistic policies, respectively. The inequalities are satisfied with the high probability of 1 − δ = (1 − δ̂)(1 − δ̄),
where δ̂ and δ̄ are the probability that Equation (4) does nothold against each policy. Hence, a novel value function in
our paper approaches to η V̂N∗ + (1 − η)V̄N∗ as closely as
ε = η ε̂ + (1 − η)ε̄ with high probability of 1 − δ for all but
a certain time steps.

Experiments
We evaluate our approach on two problem settings. One iswith randomly generated environments on which we perform Monte-Carlo simulations, and the other is with environments created from real Martian terrain data. For eachsimulation, GP hyper-parameters are obtained through trial
and error. As for the two Mars simulations, we refer tothe previous work by (Turchetta, Berkenkamp, and Krause2016) and use similar parameters.

Monte-Carlo Simulations with Randomly
Generated Problems
We consider a 30 × 30 rectangular grid, where reward andsafety values are randomly generated and assigned to states.
At every state (except for the boundary), the agent can takeone of five actions: stay, up, down, left, and right.
The exploring agent predicts the safety function g via GP
with a Radial Basis Function (RBF) kernel with the lengthscales being 2.0 and the prior variance of safety being 1.5.

Percentage explored
7.4 ± 2.9
8.6 ± 3.1
10.6 ± 2.1-

The agent also models the reward function r as a GP with
RBF kernel with the lengthscales of 2.0 and prior varianceof 1.0. Throughout the simulation, the discount rate is set to0.90 and beta is set as βt = 2, ∀t ≥ 0.
The Optimistic and Pessimistic MDPs are solved by the
BEB algorithm with the weight coefficient of 2.5 for the
exploration bonus. In the MDPs, the estimated mean of reward function is used. The safety constraint is imposed byremoving the unsafe states from S in the Optimistic MDP,
and by removing unsafe and uncertain states from S in thePessimistic MDP. The weight coefficient between optimistic
and pessimistic policies is set to η = 0.50.
Our approach is compared with the ones by (Moldovanand Abbeel 2012) and (Turchetta, Berkenkamp, and Krause2016), as well as with two non-exploratory cases where 1)
safety function is known and 2) safety and reward functions
are known a priori. In order to make a fair comparison, the
MDPs are solved by the BEB algorithm in all the targets ofcomparison ((Moldovan and Abbeel 2012) used R-MAX intheir paper). We performed a Monte-Carlo simulation with
100 samples of randomly generated reward and safety. The
number of iteration is fixed and defined as 100.
Table 2 shows the simulation results. We see that our approach outperforms the competing baselines, although there
is a sizable gap compared to more omniscient agents. Notealso that our approach explored the fewest uncertain statescompared to the baselines, and instead spent more effort
finding a good reward function within a high quality safe region. These results show that our approach is able to betterbalance the three-way trade-off between exploring the safe
region, exploring the reward function, and exploiting known
rewards.

Simulated Mars Surface Exploration
We next demonstrate S AFE E XP O PT-MDP algorithm
in simulated Mars surface exploration scenarios, as in(Moldovan and Abbeel 2012) and (Turchetta, Berkenkamp,
and Krause 2016). Two scenarios are simulated. The firstscenario uses the similar environment as in (Turchetta,
Berkenkamp, and Krause 2016), where the elevation map ofMars over a 40 by 30 meters area at latitude 30◦ 6’ Southand longitude 202◦ 2’ East is used. The elevation map isderived from the digital terrain models (DEMs), createdfrom HiRISE camera on the Mars Reconnaissance Orbiter(McEwen et al. 2007). Figure 4 shows the elevation map
whose step size is 1 meter.

Table 3: Simulation result for Martian terrain. As with Table 2, our approach achieves higher reward with less exploration.
First Mars exploration scenario
Second Mars exploration scenario
Cumulative reward Percentage explored Cumulative reward Percentage explored
S AFE E XP O PT-MDP
36
5.3 %
67
11.8 %
Moldovan & Abbeel (2012)
29
5.8 %
45
12.9 %
Turchetta et al. (2016)
5
6.7 %
0
15.3 %
Safe known
48
78
Safe/reward known
63
458
-

Figure 4: Mars digital elevation model by HiRISE data [m].

Figure 5: Data on the slope angle [deg] and ROI employed
for the landing site selection of Mars 2020 rover.

The reward function is binary: one if a state is within ascientific region of interest (ROI), and zero otherwise. We
assume a hypothetical rectangular ROI at 16 ≤ x ≤ 24 and21 ≤ y ≤ 30. The rover (i.e., the agent) starts from (0, 0),and at each time step it takes one of five actions: stay, up,
down, left, and right. We also assume that any states where
the slope is greater than 25 degrees are unsafe. Hence, thesafety function g is defined as g(s, a) = H(s) − H(s0 ),where H is the elevation of a given state. The safety threshold is h = − tan(25◦ ).
The rover predicts the elevation using a GP with a Matern
kernel with ν = 5/2. The lengthscales are 15.0 m and theprior variance over elevation is 100 m2 . We assume a noise
standard deviation of 0.075 m. To predict the reward function, it uses a GP with Radial Basis Function (RBF) kernel.
The lengthscales are 10.0 m and the prior variance over reward is 50. We set the level of confidence as βt = 2, ∀t ≥ 0.
Optimistic/Pessimistic MDPs are solved by BEB algorithmwith the weight coefficient for exploration bonus at 2.5. The
discount factor is set to 0.90, and the weight coefficient between optimistic and pessimistic policies is η = 0.50.
Our algorithm is compared against the same algorithms
as in the previous simulation. The simulation result is summarized in the left half of Table 3. The result shows that
our algorithm outperforms existing algorithms in terms ofthe cumulative reward even for the binary reward function.
We observe a 24% relative improvement in reward compared to the next best baseline. Note that since the reward
function is more peaked, the pure exploration approach by(Turchetta, Berkenkamp, and Krause 2016) performs substantially worse. We again observe that our approach con-ducts the least amount of exploration of the safety region,
which indicates it is better able to balance the trade-off between exploring the safety region and the reward function.
The second scenario is based on the potential exploration
scenario of Mars 2020 rover. We employ the real data of aregion over 1.75 by 1.75 kilometers around Holden Crater,
one of the previous candidates of the landing site of Mars
2020 rover. The safety function g is the slope angle, and weassume that every state where the slope is greater that 25
degrees is unsafe.
The reward function r is defined as a binary function
meaning scientific gain: one if a state is within a ROI, andzero otherwise. The rover predicts the elevation using a GP
with a Matern kernel with ν = 5/2. The lengthscales are
50.0 m and the prior variance over elevation is 100 deg2 . We
assume a noise standard deviation of 0.075 deg. To predictthe reward function, the agent (rover) uses a GP with RBF
kernel. The lengthscales are 50.0 m and the prior varianceover reward is 50. Other simulation settings are equivalent
with previous Mars simulation in this paper.
Figure 5 represents simulation environment in which thecolor map behind is the slope angle [deg] and the region surrounded by the black close curve represents ROI. The step
size is 25 meters, and we consider 70 × 70 rectangular grid.
The rover starts from (0, 0), and takes a route represented by
red circles. At the initial phase before arriving at the ROI,
all the area around an agent is safe and reward-poor, so anagent considers that all the environment has similar characteristics. Hence an agent explores the state space all evenly
by randomly choosing the next action. When the value functions are exactly equal between more than one state, an agent

randomly chooses the next state and action pair. However,once the agent discovers the ROI, it is likely to prefer to staywithin the ROI. In addition, observe that the rover does not
run through the states in which the slope angle is greaterthan 25 degrees, which shows that the rover succeeded insafely arriving at ROI. The simulation result is summarized
in the right half of Table 3. There is a large gap between
not knowing the safety/reward functions and knowing them(i.e., unrealistic case), which suggests there is still room for
improving the algorithm.

Conclusion
We presented a novel approach for exploring and optimizing
safety constrained MDPs. By modeling a priori unknown reward and safety via GPs, an agent can classify state space
into safe, uncertain, and unsafe regions. The agent then optimizes the linear combination of optimistic policy and pessimistic policy. This algorithm automatically encourages exploration of safety while balancing exploration/exploitationof reward. In addition, the agent can guarantee safety withhigh probability by exploiting GP structure of the safetyfunction. Finally, we demonstrated the effectiveness of our
proposed method by theoretical analysis and numerical simulation. Moving forward, it would be interesting to develop
rigorous finite-time rates of convergence (rather than just an
asymptotic consistency result), as well as extend to multiple
heterogeneous safety and reward functions.

Acknowledgement
The research was carried out in part at the Jet Propulsion
Laboratory, California Institute of Technology. This research
was also supported in part by NSF Awards #1564330 &#1637598, JPL PDF IAMS100224, a Bloomberg Data Science Research Grant, and a gift from Northrop Grumman.

References
Araya, M.; Buffet, O.; and Thomas, V. 2012. Near-optimal
BRL using optimistic local transitions. In Proceedings ofthe 29th International Conference on Machine Learning(ICML), 97–104.
Bertsekas, D. P., and Tsitsiklis, J. N. 1995. Neuro-dynamic
programming: an overview. In Decision and Control, 1995.,
Proceedings of the 34th IEEE Conference on, volume 1,
560–564. IEEE.
Blackmore, L.; Ono, M.; Bektassov, A.; and Williams,
B. C. 2010. A probabilistic particle-control approximationof chance-constrained stochastic predictive control. IEEE
transactions on Robotics 26(3):502–517.
Brafman, R. I., and Tennenholtz, M. 2002. R-max - ageneral polynomial time algorithm for near-optimal reinforcement learning. Journal of Machine Learning Research3(Oct):213–231.
Fleming, W. H., and McEneaney, W. M. 1995. Risksensitive control on an infinite time horizon. SIAM Journal
on Control and Optimization 33(6):1881–1915.

Ghosal, S., and Roy, A. 2006. Posterior consistency of Gaussian process prior for nonparametric binary regression. The
Annals of Statistics 2413–2429.
Kakade, S. M., et al. 2003. On the sample complexity ofreinforcement learning. Ph.D. Dissertation, University of
London.
Kearns, M., and Singh, S. 2002. Near-optimal reinforcement learning in polynomial time. Machine Learning 49(23):209–232.
Kolter, J. Z., and Ng, A. Y. 2009. Near-Bayesian explorationin polynomial time. In Proceedings of the 26th Annual International Conference on Machine Learning, 513–520. ACM.
McEwen, A. S.; Eliason, E. M.; Bergstrom, J. W.; Bridges,
N. T.; Hansen, C. J.; Delamere, W. A.; Grant, J. A.; Gulick,
V. C.; Herkenhoff, K. E.; Keszthelyi, L.; et al. 2007. Mars
reconnaissance orbiter’s high resolution imaging science experiment (HiRISE). Journal of Geophysical Research: Planets 112(E5).
Mockus, J. 2012. Bayesian approach to global optimization:
theory and applications, volume 37. Springer Science &Business Media.
Moldovan, T., and Abbeel, P. 2012. Safe exploration in
Markov decision processes. In International Conference on
Machine Learning (ICML).
Ono, M.; Pavone, M.; Kuwata, Y.; and Balaram, J. 2015.
Chance-constrained dynamic programming with application
to risk-aware robotic space exploration. Auton. Robots
39(4):555–571.
Rasmussen, C. E. 2006. Gaussian processes for machinelearning.
Schwarm, A. T., and Nikolaou, M. 1999. Chanceconstrained model predictive control.
AIChE Journal
45(8):1743–1752.
Srinivas, N.; Krause, A.; Kakade, S. M.; and Seeger, M.
2010. Gaussian process optimization in the bandit setting:
No regret and experimental design. In International Conference on Machine Learning (ICML).
Strehl, A. L.; Li, L.; Wiewiora, E.; Langford, J.; andLittman, M. L. 2006. PAC model-free reinforcement learning. In Proceedings of the 23rd International Conference on
Machine Learning (ICML), 881–888. ACM.
Strens, M. 2000. A bayesian framework for reinforcement
learning. In International Conference on Machine Learning(ICML), 943–950.
Sui, Y.; Gotovos, A.; Burdick, J. W.; and Krause, A. 2015.
Safe exploration for optimization with gaussian processes.
In International Conference on Machine Learning (ICML).
Sutton, R. S., and Barto, A. G. 1998. Reinforcement learning: An introduction, volume 1. MIT press Cambridge.
Turchetta, M.; Berkenkamp, F.; and Krause, A. 2016. Safe
exploration in finite Markov decision processes with gaussian processes. In Advances In Neural Information Processing Systems, 4305–4313.