A Risk-Sensitive Finite-Time Reachability Approach for Safety of Stochastic Dynamic Systems

TL;DR

This paper introduces a risk-sensitive reachability framework for stochastic systems that quantifies and reduces the probability of constraint violations over finite time horizons using CVaR, with practical algorithms and an illustrative example.

Contribution

It develops a novel risk-sensitive safe set concept, reduces its computation to a CVaR-based MDP, and provides a tractable algorithm with theoretical guarantees for stochastic system safety.

Findings

The approach allows tuning risk sensitivity from worst-case to risk-neutral.

A tractable algorithm for approximating risk-sensitive safe sets is proposed.

Demonstration on stormwater catchment design shows practical utility.

Abstract

A classic reachability problem for safety of dynamic systems is to compute the set of initial states from which the state trajectory is guaranteed to stay inside a given constraint set over a given time horizon. In this paper, we leverage existing theory of reachability analysis and risk measures to devise a risk-sensitive reachability approach for safety of stochastic dynamic systems under non-adversarial disturbances over a finite time horizon. Specifically, we first introduce the notion of a risk-sensitive safe set as a set of initial states from which the risk of large constraint violations can be reduced to a required level via a control policy, where risk is quantified using the Conditional Value-at-Risk (CVaR) measure. Second, we show how the computation of a risk-sensitive safe set can be reduced to the solution to a Markov Decision Process (MDP), where cost is assessed according to CVaR. Third, leveraging this reduction, we devise a tractable algorithm to approximate a risk-sensitive safe set, and provide theoretical arguments about its correctness. Finally, we present a realistic example inspired from stormwater catchment design to demonstrate the utility of risk-sensitive reachability analysis. In particular, our approach allows a practitioner to tune the level of risk sensitivity from worst-case (which is typical for Hamilton-Jacobi reachability analysis) to risk-neutral (which is the case for stochastic reachability analysis).

Tables

Table 1. TABLE I:

Sample moment	Value
Mean	12.16 ft3/s
Variance	3.22 ft6/s2
Skewness	1.68 ft9/s3
Disturbance sample, $d_j$ ft3/s	Probability, $\mathbb{P}\left[w_k=d_j\right]$
$8.57$	$0.0236$
$9.47$	${10}^{-4}$
$10.37$	${10}^{-4}$
$11.26$	$0.5249$
$12.16$	$0.3272$
$13.06$	${10}^{-4}$
$13.95$	${10}^{-4}$
$14.85$	${10}^{-4}$
$15.75$	${10}^{-4}$
$16.65$	$0.1237$

Equations

$$\text{CVaR}_{\alpha}[Z]:={\underset{t\in\mathbb{R}}{\min}}\text{ }\Big{\{}t+\frac{1}{\alpha}\mathbb{E}\big{[}\max\{Z-t,0\}\big{]}\Big{\}}.$$

$$x_{k+1}=f(x_k,u_k,w_k),k=0,1,\dots ,N-1,$$

$$\Pi:=\big{\{}(\mu_{0},\mu_{1},\dots,\mu_{N-1})\mid\mu_{k}:H_{k}\rightarrow U\big{\}},$$

$$x\in \mathcal{K}⟺g(x)<0.$$

$${\mathcal{S}}_{\alpha }^r:=\{x\in \mathcal{X}\mid W_0^{\ast }(x,\alpha )\le r\},$$

$$\displaystyle W_{0}^{*}(x,\alpha):={\underset{\pi\in\Pi}{\min}}\text{ CVaR}_{\alpha}\big{[}Z^{\pi}_{x}\big{]},$$

$$\displaystyle Z^{\pi}_{x}:={\underset{k=0,\dots,N}{\max}}\big{\{}g(x_{k})\big{\}},$$

$$\begin{split}\text{ CVaR}_{\alpha}\big{[}Z^{\pi^{*}}_{x}\big{]}={\underset{\pi\in\Pi}{\inf}}\text{ CVaR}_{\alpha}\big{[}Z^{\pi}_{x}\big{]}={\underset{\pi\in\Pi}{\min}}\text{ CVaR}_{\alpha}\big{[}Z^{\pi}_{x}\big{]}.\end{split}$$

$${\text{1}}_{\bar{\mathcal{K}}}(x):=\{\begin{array}{l}1\text{ if }x\notin \mathcal{K}\\ 0\text{ if }x\in \mathcal{K}\end{array}.$$

$$\Big{\{}x\in\mathcal{X}\text{ }\Big{|}\text{ }{\underset{\pi\in\Pi}{\min}}\text{ }\mathbb{E}\big{[}\max_{k=0,\dots,N}\textbf{1}_{\mathcal{\bar{K}}}(x_{k})\big{]}\leq\epsilon\Big{\}},$$

$$y_{k+1}={\bar{R}}_{x_k,y_k}(w_k)\cdot y_k,$$

$${\bar{\Pi }}_t:=\{({\bar{\mu }}_t,{\bar{\mu }}_{t+1},\dots ,{\bar{\mu }}_{N-1})\mid {\bar{\mu }}_k:\mathcal{X}\times \mathcal{Y}\to U\},$$

$$t=0,\dots ,N-1.$$

$${\mathcal{U}}_{\alpha }^r:=\{x\in \mathcal{X}\mid J_0^{\ast }(x,\alpha )\le \beta e^{m\cdot r}\},$$

$$\displaystyle J_{0}^{*}(x,\alpha):={\underset{\pi\in\bar{\Pi}_{0}}{\min}}\text{ CVaR}_{\alpha}\big{[}Y_{x}^{\pi}\big{]},$$

$$Y_x^{\pi }:={\sum }_{k=0}^{N}c(x_k),$$

$${\text{CVaR}}_{y_k}[Z|x_k,{\pi }_k]=\underset{R\in \mathcal{R}(y_k,\mathbb{P})}{\max }C(R,Z;x_k,y_k,{\pi }_k),$$

$$C(R,Z;x_k,y_k,{\pi }_k):=$$

$$\displaystyle\mathbb{E}_{w_{k}\sim\mathbb{P}}\big{[}R(w_{k})\cdot\text{CVaR}_{y_{k}R(w_{k})}[Z|x_{k+1},\pi_{k+1}]\big{|}x_{k},\mu_{k}\big{]},$$

$$\displaystyle\mathcal{R}(y_{k},\mathbb{P}):=\big{\{}R:D\to\big{(}0,\textstyle\frac{1}{y_{k}}\big{]}\mathrel{\big{|}}\mathbb{E}_{w_{k}\sim\mathbb{P}}\big{[}R(w_{k})\big{]}=1\big{\}},$$

$$Z:={\sum }_{i=k+1}^{N}c(x_i),$$

$${\bar{R}}_{x_k,y_k}(w_k)=\underset{R\in \mathcal{R}(y_k,\mathbb{P})}{\mathrm{argmax}}C(R,Z;x_k,y_k,{\pi }_k).$$

$$J_k(z_k)$$

$$\displaystyle:={\underset{u\in U}{\min}}\Big{\{}c(x_{k})+{\underset{R\in\mathcal{R}(y_{k},\mathbb{P})}{\max}}\mathbb{E}_{w_{k}\sim\mathbb{P}}\big{[}RJ_{k+1}(x^{\prime},y_{k}R)\big{|}z_{k},u\big{]}\Big{\}},$$

$$\text{ for }k=N-1,\dots ,0,$$

$$J_0(x,\alpha )=J_0^{\ast }(x,\alpha ),$$

$$x_{k+1}=x_k+\frac{△t}{A}(w_k-q_p(x_k,u_k)),k=0,\dots ,N-1,$$

$$q_p(x_k,u_k):=\{\begin{array}{ll}{C}_d\pi r^2{u}_k\sqrt{2\eta (x-E)}& \text{ if }{x}_k\ge E\\ 0& \text{ if }{x}_k<E,\end{array}$$

$$\max \{y_1,\dots ,y_p\}$$

$$\le \max \{y_1,\dots ,y_p\}+\frac{\log p}{m},$$

$$\frac{1}{m}\log (e^{m{y}_1}+\cdots +e^{m{y}_p})\to \max \{y_1,\dots ,y_p\}.$$

$${\text{CVaR}}_{\alpha }[\log (Z)]\le \log \left({\text{CVaR}}_{\alpha }[Z]\right),$$

$$\displaystyle\text{CVaR}_{\alpha}\big{[}Z_{x}^{\pi}\big{]}$$

$$\displaystyle\leq\textstyle\frac{1}{m}\log\left(\text{CVaR}_{\alpha}\big{[}\bar{Y}_{x}^{\pi}\big{]}\right),$$

$$e^{m\cdot r}\geq{\underset{\pi\in\bar{\Pi}_{0}}{\min}}\text{ CVaR}_{\alpha}\big{[}Y_{x}^{\pi}/\beta\big{]}\geq{\underset{\pi\in\Pi}{\min}}\text{ CVaR}_{\alpha}\big{[}Y_{x}^{\pi}/\beta\big{]},$$

$$r\geq\textstyle\frac{1}{m}\log\left(\text{CVaR}_{\alpha}\big{[}Y_{x}^{\pi}/\beta\big{]}\right)\geq\text{CVaR}_{\alpha}\big{[}Z_{x}^{\pi}\big{]},\\$$

$$c(x_k)+\underset{R\in \mathcal{R}(y_k,\mathbb{P})}{\max }\mathbb{E}[R{J}_{k+1}(x_{k+1},y_{k}R)|z_k,{\mu }_k^{ϵ}]\le J_k(z_k)+ϵ.$$

$$J_{k}^{\epsilon}(z_{k}):=\text{CVaR}_{y_{k}}\big{[}\textstyle\sum_{i=k}^{N}c(x_{i})\big{|}z_{k},\pi_{k}^{\epsilon}\big{]},$$

Full text

I Introduction

Reachability analysis is a formal verification method based on optimal control theory that is used to prove safety or performance properties of dynamic systems [1]. A classic reachability problem for safety is to compute the set of initial states from which the state trajectory is guaranteed to stay inside a given constraint set over a given time horizon. This problem was first considered for discrete-time dynamic systems by Bertsekas and Rhodes under the assumption that disturbances are uncertain but belong to known sets [2], [3], [4]. In this context, the problem is solved using a minimax formulation, in which disturbances behave adversarially and safety is described as a binary notion based on set membership [2], [3], [4, Sec. 3.6.2].

In practice, minimax formulations can yield overly conservative solutions, particularly because disturbances are usually non-adversarial. Most storms do not cause major floods, and most vehicles are not involved in pursuit-evasion games. If there are enough observations of the system, one can estimate a probability distribution for the disturbance (e.g., see [5]), and then assess safety properties of the system in a more realistic context. For stochastic discrete-time dynamic systems, Abate et al. [6] developed an algorithm to compute a set of initial states from which the probability of safety of the state trajectory can be increased to a required level by a control policy.111Safety of the state trajectory is the event that the state trajectory stays in the constraint set over a finite time horizon. Summers and Lygeros [7] extended the algorithm of Abate et al. to quantify the probability of safety and performance of the state trajectory, by specifying that the state trajectory should also reach a target set.

Both the stochastic reachability methods [6], [7] and the minimax reachability methods [2], [3], [4] for discrete-time dynamic systems describe safety as a binary notion based on set membership. In Abate et al., for example, the probability of safety to be optimized is formulated as an expectation of a product (or maximum) of indicator functions, where each indicator encodes the event that the state at a particular time point is inside a given set [6]. The stochastic reachability methods [6], [7] do not generalize to quantify the distance between the state trajectory and the boundary of the constraint set, since they use indicator functions to convert probabilities to expectations to be optimized.

In contrast, Hamilton-Jacobi (HJ) reachability methods quantify the deterministic analogue of this distance for continuous-time systems subject to adversarial disturbances (e.g., see [1], [8], [9], [10]). Quantifying the distance between the state trajectory and the boundary of the constraint set in a non-binary fashion may be important in applications where the boundary is not known exactly, or where mild constraint violations are inevitable, but extreme constraint violations must be avoided.

It is imperative that reachability methods for safety take into account the possibility that rare events can occur with potentially damaging consequences. Reachability methods that assume adversarial disturbances (e.g., [1], [3]) suppose that harmful events will always occur, which may yield solutions with limited practical utility, especially in applications with large uncertainty sets. Stochastic reachability methods [6], [7] do not explicitly account for rare high-consequence events, as costs are evaluated in terms of an expectation.

In contrast, in this paper, we harness risk measure theory to formulate a reachability analysis approach that explicitly accounts for the possibility of rare events with negative consequences: harmful events are likely to occur at some time, but they are unlikely to occur all the time. Specifically, a risk measure is a function that maps a random variable $Z$ representing a loss, or a cost, into the real line, according to the possibility of danger associated with $Z$ [11, Sec. 6.3], [12, Sec. 2.2]. Risk-sensitive optimization has been studied in applied mathematics [13], reinforcement learning [14], [15], [16], and optimal control [17], [18]. A risk-sensitive method may provide more practical and protective decision-making machinery (versus stochastic or minimax methods) by encoding a flexible degree of conservativeness.

In this paper, we use a particular risk measure, called Conditional Value-at-Risk (CVaR). If $Z$ is a random variable representing cost with finite expectation, then the Conditional Value-at-Risk of $Z$ at the confidence level $\alpha\in(0,1]$ is defined as [11, Equation 6.22],222Conditional Value-at-Risk is also called Average Value-at-Risk, which is abbreviated as AV@R in [11].

[TABLE]

CVaR captures a full spectrum of risk assessments from risk-neutral to worst-case, since $\text{CVaR}_{\alpha}[Z]$ increases from $\mathbb{E}[Z]$ to $\text{ess}\sup Z$, as $\alpha$ decreases from 1 to 0. CVaR has desirable mathematical properties for optimization [19] and chance-constrained stochastic control [20]. There is a well-established relationship between CVaR and chance constraints that we will use to obtain probabilistic safety guarantees in this paper. Please see [12] and [21] for additional background on CVaR.

Statement of Contributions. This paper introduces a risk-sensitive reachability approach for safety of stochastic dynamic systems under non-adversarial disturbances over a finite time horizon. Specifically, the contributions are four-fold. First, we introduce the notion of a risk-sensitive safe set as a set of initial states from which the risk of large constraint violations can be reduced to a required level via a control policy, where risk is quantified using the Conditional Value-at-Risk (CVaR) measure. Our formulation explicitly assesses the distance between the boundary of the constraint set and the state trajectory of a stochastic dynamic system. This is an extension of stochastic reachability methods (e.g., [6], [7]), which replace this distance with a binary random variable. Further, in contrast to stochastic reachability methods, our formulation explicitly accounts for rare high-consequence events, by posing the optimal control problem in terms of CVaR, instead of a risk-neutral expectation. Second, we show how the computation of a risk-sensitive safe set can be reduced to the solution to a Markov Decision Process (MDP), where cost is assessed according to CVaR. Third, leveraging this reduction, we devise a tractable algorithm to approximate a risk-sensitive safe set and provide theoretical arguments to justify its correctness. Finally, we present a realistic example inspired from stormwater catchment design to demonstrate the utility of risk-sensitive reachability analysis.

Organization. The rest of this paper is organized as follows. We present the problem formulation and define risk-sensitive safe sets in Sec. II. In Sec. III, we show how the computation of a risk-sensitive safe set can be reduced to the solution to a CVaR-MDP problem, i.e., an MDP where cost is assessed according to CVaR. In Sec. IV, we present a value-iteration algorithm to approximate risk-sensitive safe sets, along with theoretical arguments that support its correctness. In Sec. V, we provide numerical experiments on a realistic example inspired from stormwater catchment design. Finally, in Sec. VI, we draw conclusions and discuss directions for future work.

II Problem Formulation

II-A System Model

We consider a fully observable stochastic discrete-time dynamic system over a finite time horizon [4, Sec. 1.2],

[TABLE]

such that $x_{k}\in\mathcal{X}\subseteq\mathbb{R}^{n}$ is the state of the system at time $k$, $u_{k}\in U$ is the control at time $k$, and $w_{k}\in D$ is the random disturbance at time $k$. The control space $U$ and disturbance space $D$ are finite sets of real-valued vectors. The function $f:\mathcal{X}\times U\times D\rightarrow\mathcal{X}$ is bounded and Lipschitz continuous. The probability that the disturbance equals $d_{j}\in D$ at time $k$ is $\mathbb{P}[w_{k}=d_{j}]=p_{j}$, where $0\leq p_{j}\leq 1$ and $\sum_{j=1}^{W}p_{j}=1$. We assume that $w_{k}$ is independent of $x_{k}$, $u_{k}$, and disturbances at any other times. The only source of randomness in the system is the disturbance. In particular, the initial state $x_{0}$ is not random. The set of admissible, deterministic, history-dependent control policies is,

[TABLE]

where $H_{k}:=\underbrace{\mathcal{X}\times\ldots\times\mathcal{X}}_{\text{(k+1) times}}$ is the set of state histories up to time $k$. We are given a constraint set $\mathcal{K}\subseteq\mathcal{X}$, and the safety criterion that the state of the system should stay inside $\mathcal{K}$ over time. For example, if the system is a pond in a stormwater catchment, then $x_{k}$ may be the water level of the pond in feet at time $k$, and $\mathcal{K}=[0,5)$ indicates that the pond overflows if the water level exceeds 5 feet. We quantify the extent of constraint violation/satisfaction using a surface function that characterizes the constraint set. Specifically, similar to [9, Eq. 2.3], let $g:\mathcal{X}\rightarrow\mathbb{R}$ satisfy,

[TABLE]

For example, we may choose $g(x)=x-5$ to characterize $\mathcal{K}=[0,5)$ on the state space $\mathcal{X}=[0,\infty)$.

II-B Risk-Sensitive Safe Sets

A risk-sensitive safe set is a set of initial states from which the risk of large constraint violations can be reduced to a required level via a control policy, where risk is quantified using the CVaR measure. We use the term risk level to mean the allowable level of risk of constraint violations. Formally, the risk-sensitive safe set at the confidence level $\alpha\in(0,1]$ and the risk level $r\in\mathbb{R}$ is defined as,

[TABLE]

such that the state trajectory, $(x_{0},x_{1},...,x_{N})$, evolves according to the dynamics model (2) with the initial state $x_{0}=x$ under the policy $\pi\in\Pi$. The surface function $g$ characterizes distance to the constraint set $\mathcal{K}$ according to (4). Note that the minimum in the definition of $W_{0}^{*}(x,\alpha)$ is attained, as the next lemma states.

Lemma 1** (Existence of a minimizer)**

For any initial state $x_{0}\in\mathcal{X}$ and any confidence level $\alpha\in(0,1]$, there exists a policy $\pi^{*}\in\Pi$ such that

[TABLE]

Proof:

Fix the initial state $x_{0}$. Since the control and disturbance spaces are finite, the set of states that could be visited (starting from $x_{0}$) is finite. Therefore, the set of policies restricted to realizable histories from $x_{0}$ is finite. Hence, the infimum must be attained by some policy $\pi^{*}$. ∎

In the next sections, we will present a tractable algorithm to approximately compute risk-sensitive safe sets at different levels of confidence and risk.

II-C Discussion

Computing risk-sensitive safe sets, as defined by (5), is well-motivated for several reasons. Our formulation incorporates different confidence levels and non-binary distance to the constraint set. In contrast, the stochastic reachability problem addressed by Abate et al. [6] uses a single confidence level and an indicator function to measure distance to the constraint set, in order to quantify the probability of constraint violation. Specifically, let $\epsilon\in[0,1]$ be the maximum tolerable probability of constraint violation (called safety level in [6]), and choose $\alpha:=1$, $r:=\epsilon-\frac{1}{2}$, and $g(x):=\textbf{1}_{\bar{\mathcal{K}}}(x)-\frac{1}{2}$, where

[TABLE]

which is the maximal probabilistic safe set at the $\epsilon$-safety level [6, Eqs. 11 and 13], if we consider non-hybrid dynamic systems that evolve under history-dependent policies.333Abate et al. [6] considers hybrid dynamic systems that evolve under Markov policies.

Risk-sensitive safe sets have two desirable mathematical properties. The first property is that $\mathcal{S}_{\alpha}^{r}$ shrinks as the risk level $r$ or the confidence level $\alpha$ decreases. Since $\mathcal{S}_{\alpha}^{r}$ is an $r$-sublevel set and $\text{CVaR}_{\alpha}$ increases as $\alpha$ decreases, one can show that $\mathcal{S}_{\alpha_{2}}^{r_{2}}\subseteq\mathcal{S}_{\alpha_{1}}^{r_{2}}\subseteq\mathcal{S}_{\alpha_{1}}^{r_{1}}$ and $\mathcal{S}_{\alpha_{2}}^{r_{2}}\subseteq\mathcal{S}_{\alpha_{2}}^{r_{1}}\subseteq\mathcal{S}_{\alpha_{1}}^{r_{1}}$ hold for any $r_{1}\geq r_{2}$ and $1\geq\alpha_{1}\geq\alpha_{2}>0$. In other words, as the allowable level of risk of constraint violation ($r$) decreases, or as the fraction of damaging outcomes that are not fully addressed ($\alpha$) decreases, $\mathcal{S}_{\alpha}^{r}$ encodes a higher degree of safety.

The second property is that risk-sensitive safe sets at the risk level, $r:=0$, have probabilistic safety guarantees.

Lemma 2** (Probabilistic safety guarantee)**

If $x\in\mathcal{S}_{\alpha}^{0}$, then the probability that the state trajectory initialized at $x$ exits $\mathcal{K}$ can be reduced to $\alpha$ by a control policy.

Proof:

The proof follows from the fact that $\text{CVaR}_{\alpha}[Z_{x}^{\pi}]\leq 0\implies\mathbb{P}[Z_{x}^{\pi}\geq 0]\leq\alpha$  [11, Sec. 6.2.4, pp. 257-258]. The event $Z_{x}^{\pi}\geq 0$ is equivalent to the event that there is a state $x_{k}$ of the associated trajectory that exits the constraint set, since $g(x_{k})\geq 0\iff x_{k}\notin\mathcal{K}$. ∎

Lemma 2 indicates that $\mathcal{S}_{\alpha}^{0}$ is a subset of the maximal probabilistic safe set at the safety level $\alpha\in(0,1]$, if we consider non-hybrid dynamic systems that evolve under history-dependent policies [6, Eqs. 9 and 11].

A key difference between our risk-sensitive safe set (5) and the risk-constrained safe set in [18] is that we specify the CVaR of the worst constraint violation of the state trajectory $(x_{0},\dots,x_{N})$ to be below a required threshold, while ref. [18] specifies the CVaR of the constraint violation of $x_{k}$ to be below a required threshold for each $k$.

III Reduction of Risk-Sensitive Safe Set Computation to CVaR-MDP

Computing risk-sensitive safe sets is challenging since the computation involves a maximum of costs (as opposed to a summation of costs) and the Conditional Value-at-Risk measure (as opposed to an expectation). In this section, we show how computing an under-approximation of a risk-sensitive safe set can be reduced to solving a CVaR-MDP, which has been studied, for example, by [15] and [22]. Such a reduction will be leveraged in Section IV to devise a value-iteration algorithm to compute tractable approximations of risk-sensitive safe sets.

III-A Preliminaries

The reduction procedure is inspired by Chow et al. [15]. Specifically, we consider an augmented state space, $\mathcal{X}\times\mathcal{Y}$, that consists of the original state space, $\mathcal{X}$, and the space of confidence levels, $\mathcal{Y}:=(0,1]$. The under-approximations of risk-sensitive safe sets will be defined in terms of the dynamics of the augmented state, $(x,y)\in\mathcal{X}\times\mathcal{Y}$.

Let $(x_{0},y_{0})=(x,\alpha)$ be a given initial condition. The augmented state at time $k+1$, $(x_{k+1},y_{k+1})$, depends on the augmented state at time $k$, $(x_{k},y_{k})$, as follows. Given a control $u_{k}\in U$ and a sampled disturbance $w_{k}\in D$, the next state $x_{k+1}\in\mathcal{X}$ satisfies the dynamics model (2). The next confidence level $y_{k+1}\in\mathcal{Y}$ is given by,

[TABLE]

where $\bar{R}_{x_{k},y_{k}}:D\to(0,\frac{1}{y_{k}}]$ is a known deterministic function, which we will specify in Lemma 3. The augmented state space $\mathcal{X}\times\mathcal{Y}$ is fully observable. Indeed, the history of states and actions, $(x_{0},u_{0},\ldots,x_{k-1},u_{k-1},x_{k})$, is available at time $k$ by (2). Also, the history of confidence levels, $(y_{0},\ldots,y_{k})$, is available at time $k$, since the functions $\bar{R}_{x_{k},y_{k}}$ and the initial confidence level, $y_{0}=\alpha$, are known.

We define the set of deterministic, Markov control policies in terms of the augmented state space as follows,

[TABLE]

There is an important distinction between the set of policies $\bar{\Pi}_{0}$ as defined above, and the set of policies $\Pi$ as defined in (3). Given $\bar{\pi}_{0}\in\bar{\Pi}_{0}$, the control law at time $k$, $\bar{\mu}_{k}\in\bar{\pi}_{0}$, only depends on the current state $x_{k}\in\mathcal{X}$ and the current confidence level $y_{k}\in\mathcal{Y}$. However, given $\pi\in\Pi$, the control law at time $k$, $\mu_{k}\in\pi$, depends on the state history up to time $k$, $(x_{0},\dots,x_{k})\in H_{k}$. In particular, the set of policies $\bar{\Pi}_{0}$ is included in the set of policies $\Pi$. This is because the augmented state at time $k$ is uniquely determined by the initial confidence level and the state history up to time $k$.

The benefits of considering $\bar{\Pi}_{0}$ instead of $\Pi$ are two-fold. First, the computational requirements are reduced when the augmented state at time $k$, $(x_{k},y_{k})$, is processed instead of the initial confidence level and the state history up to time $k$, $(y_{0},x_{0},x_{1},\ldots,x_{k})$. Second, we are able to define an under-approximation of the risk-sensitive safe set given by (5) using $\bar{\Pi}_{0}$, which we explain below.

III-B Under-Approximation of Risk-Sensitive Safe Set

Define the set $\mathcal{U}_{\alpha}^{r}\subseteq\mathcal{X}$, at the confidence level $\alpha\in(0,1]$ and the risk level $r\in\mathbb{R}$,

[TABLE]

where

[TABLE]

such that $c:\mathcal{X}\to\mathbb{R}$ is a stage cost, and the augmented state trajectory, $(x_{0},y_{0},\dots,x_{N-1},y_{N-1},x_{N})$, satisfies (2) and (7) with the initial condition $(x_{0},y_{0})=(x,\alpha)$ under the policy $\pi\in\bar{\Pi}_{0}$. The next theorem, whose proof is provided in the Appendix, states that if the stage cost takes a particular form, then $\mathcal{U}_{\alpha}^{r}$ is an under-approximation of the risk-sensitive safe set $\mathcal{S}_{\alpha}^{r}$.

Theorem 1** (Reduction to CVaR-MDP)**

Choose the stage cost $c(x):=\beta e^{m\cdot g(x)}$, where $\beta>0$ and $m>0$ are constants, and $g$ satisfies (4). Then, $\mathcal{U}_{\alpha}^{r}$ as defined in (9) is a subset of $\mathcal{S}_{\alpha}^{r}$ as defined in (5). Further, the gap between $\mathcal{U}_{\alpha}^{r}$ and $\mathcal{S}_{\alpha}^{r}$ can be reduced by increasing $m$.                  $\blacksquare$

In the definition of the stage costs, the parameter $\beta$ is included to address numerical issues that may arise, if $m$ is set to a very large number.

IV A Value-Iteration Algorithm

to Approximate Risk-Sensitive Safe Sets

By leveraging Theorem 1, one can use existing CVaR-MDP algorithms to compute under-approximations of risk-sensitive safe sets. In this paper, we adapt a value-iteration algorithm from Chow et al. [15] to compute tractable approximations of the risk-sensitive safe set under-approximations $\{\mathcal{U}_{\alpha}^{r}\}$. We start by stating an existing temporal decomposition result for CVaR that will be instrumental to devising the value-iteration algorithm.

IV-A Temporal Decomposition of Conditional Value-at-Risk

In this section, we present an existing result (namely, Lemma 22 in [23]) that specifies how the Conditional Value-at-Risk of a sum of costs can be partitioned over time, and how the confidence level evolves over time, which motivates the choice of the update function (7).

Lemma 3** (Temporal decomposition of CVaR)**

At time $k$, suppose that the system (2) is at the state $x_{k}\in\mathcal{X}$ with the confidence level $y_{k}\in\mathcal{Y}$ and is subject to a policy $\pi_{k}:=(\mu_{k},\pi_{k+1})\in\bar{\Pi}_{k}$. Then,

[TABLE]

such that $c:\mathcal{X}\to\mathbb{R}$ is a stage cost. Further, given the current state $(x_{k},y_{k})$, the current control $u_{k}:=\mu_{k}(x_{k},y_{k})$, and the next state $x_{k+1}$, the function that was introduced in (7) $\bar{R}_{x_{k},y_{k}}:D\to(0,\frac{1}{y_{k}}]$ is defined as,

[TABLE]

Remark 1

The proof of Lemma 3 is a consequence of Lemma 22 in [23], and its proof is omitted for brevity.

Remark 2

If we do not have access to $w_{k}$, but only to $(x_{k},y_{k},u_{k},x_{k+1})$, then the next confidence level is defined as $y_{k+1}:=\bar{R}_{x_{k},y_{k}}(w)$, where $w\in D$ is any disturbance that satisfies $x_{k+1}=f(x_{k},u_{k},w)$.

Remark 3

$\text{CVaR}_{y_{k}}[Z|x_{k},\pi_{k}]$* is the risk of the cumulative cost of the trajectory, $(x_{k+1},\dots,x_{N})$, that is initialized at the state $x_{k}$ with the confidence level $y_{k}$ and is subject to the policy $\pi_{k}\in\bar{\Pi}_{k}$.*

IV-B Value-Iteration Algorithm

Using Lemma 3, we will devise a dynamic programming value-iteration algorithm to compute an approximation $J_{0}$ of $J_{0}^{*}$, and thus, an approximation of $\mathcal{U}_{\alpha}^{r}$ at different levels of confidence $\alpha$ and risk $r$.

Specifically, compute the functions $J_{N-1}$, …, $J_{0}$ recursively as follows: for all $z_{k}:=(x_{k},y_{k})\in\mathcal{X}\times\mathcal{Y}$,

[TABLE]

where $J_{N}(x_{N},y_{N}):=c(x_{N})$, $c(x):=\beta e^{m\cdot g(x)}$, $x^{\prime}:=x_{k+1}$ satisfies (2), and $\mathcal{R}(y_{k},\mathbb{P})$ is defined in (11).

Then, we approximate the set $\mathcal{U}_{\alpha}^{r}$ as $\widehat{\mathcal{U}}_{\alpha}^{r}:=\left\{x\in\mathcal{X}\mid J_{0}(x,\alpha)\leq\beta e^{m\cdot r}\right\}$, where we have replaced $J_{0}^{*}$ in (9) with $J_{0}$. The function, $J_{0}$, is obtained from the last step of the value iteration (13).

In the Appendix, we present theoretical arguments inspired by [15] and [4, Sec. 1.5] that justify such an approximation. In particular, we provide theoretical evidence for the following conjecture.

Conjecture (C): Assume that the functions $J_{N-1}$, …, $J_{0}$ are computed recursively as per the value-iteration algorithm (13). Then, for any $(x,\alpha)\in\mathcal{X}\times\mathcal{Y}$,

[TABLE]

where $J_{0}^{*}$ is given by (10).

This conjecture is further supported by numerical experiments presented next.

V Numerical Experiments

In this section, we provide empirical results that demonstrate the following: 1) our value-iteration estimate of $J_{0}$ is close to a Monte Carlo estimate of $J_{0}^{*}$, 2) our value-iteration estimate of $\widehat{\mathcal{U}}_{y}^{r}$ is an under-approximation of a Monte Carlo estimate of $\mathcal{S}_{y}^{r}$, and 3) estimating $J_{0}$ (and $\widehat{\mathcal{U}}_{y}^{r}$) via the value-iteration algorithm is tractable on a realistic example inspired from stormwater catchment design. Item 1 provides empirical support for the Conjecture. Items 2 and 3 provide empirical support for reducing the computation of risk-sensitive safe sets to a CVaR-MDP. In our experiments, we used MATLAB R2016b (The MathWorks, Inc., Natick, MA) and MOSEK (Copenhagen, Denmark) with CVX [24] on a standard laptop (64-bit OS, 16GB RAM, Intel® CoreTM i7-4700MQ CPU @ 2.40GHz). Our code is available at https://github.com/chapmanmp/ACC_2019_Github.

This section demonstrates the utility of computing approximate risk-sensitive safe sets in a practical setting: to evaluate the design of a retention pond in a stormwater catchment system. We consider a retention pond from our prior work [25] as a stochastic discrete-time dynamic system,

[TABLE]

where $x_{k}\geq 0$ is the water level of the pond in feet at time $k$, $u_{k}\in U:=\{0,1\}$ is the valve setting at time $k$, and $w_{k}\in D:=\{d_{1},\dots,d_{10}\}$ is the random surface runoff in feet-cubed-per-second at time $k$. ($\eta=32.2\frac{\text{ft}}{\text{s}^{2}}$ is the acceleration due to gravity, $\pi\approx 3.14$, $r=\frac{1}{3}$ft is the outlet radius, $A=28,292$ft2 is the pond surface area, $C_{d}=0.61$ is the discharge coefficient, and $E=1$ft is the elevation of the outlet.) We estimated a finite probability distribution for $w_{k}$ using the surface runoff samples that we previously generated from a time-varying design storm (a synthetic storm based on historical rainfall) [25]. We averaged each sample over time and solved for a distribution that satisfied the empirical statistics of the time-averaged samples (Table I). We set $\triangle t:=300$ seconds, and $N:=48$ to yield a 4-hour time horizon. We chose the constraint set $\mathcal{K}:=[0,5)$, and $g(x):=x-5$.

We computed over a grid of states and confidence levels $G:=G_{s}\times G_{c}$, where $G_{s}:=\{0,0.1,\dots,6.4,6.5\}$, and $G_{c}:=\{0.999,0.95,0.80,0.65,0.5,0.35,0.20,0.05,0.001\}$. Since the initial state $x_{0}$ is non-negative and the smallest realization of $w_{k}$ is about 8.5$\frac{\text{ft}^{3}}{\text{s}}$, $x_{k+1}\geq x_{k}$ for all $k$. If $x_{k+1}>6.5$ft, we set $x_{k+1}:=6.5$ft to stay within the grid.

We were able to empirically assess the accuracy of our proposed approach because an optimal control policy is known a priori for the one-pond system. Since $x_{k+1}\geq x_{k}$ for all $k$, and the only way to exit the constraint set is if $x_{k}\geq 5$ft, an optimal policy is to keep the valve open over all time, regardless of the current state, the current confidence level, or the state history up to the current time.

Our value-iteration estimate of the function $J_{0}$ is shown in Fig. 1, and a Monte Carlo estimate of the function $J_{0}^{*}$ is shown in Fig. 2. The estimates of $J_{0}$ and $J_{0}^{*}$ are similar throughout the grid except near the smaller confidence levels. The average (largest) difference normalized by the estimate of $J_{0}^{*}$ is approximately $1.4$ ($18.7$). The average (largest) difference normalized by the estimate of $J_{0}$ is approximately $0.23$ ($0.95$). These results provide empirical support for the Conjecture.

Our value-iteration estimate of the set $\widehat{\mathcal{U}}_{\alpha}^{r}$ and a Monte Carlo estimate of the set $\mathcal{S}_{y}^{r}$ are shown in Fig. 3 at different levels of confidence $y$ and risk $r$. The empirical results indicate that $\widehat{\mathcal{U}}_{\alpha}^{r}$ is an under-approximation of $\mathcal{S}_{y}^{r}$. We estimated the sets $\{\mathcal{S}_{y}^{r}\}$ using a Monte Carlo estimate of the function $W_{0}^{*}$, which is shown in Fig. 4.

The computation time for our value-iteration estimate of $J_{0}$ was roughly 3h 6min. We deem this performance to be acceptable because

1. computations to evaluate design choices are performed off-line,
2. the problem entailed a realistically sized state space ($|G_{s}|\cdot|G_{c}|=594$ grid points) and time horizon ($N=48$ time points), and 3) our implementation is not yet optimized. Further, there is recent work in scalable approximations of reachable sets (e.g., [9]) that we will investigate for possible extensions to the risk-sensitive case.

Value-iteration implementation. To implement the value-iteration algorithm, we used the interpolation method over the confidence levels proposed by Chow et al. [15] to approximate the expectation in (13) as a piecewise linear concave function, which we maximized by solving a linear program. Further, at each $\alpha\in G_{c}$, we used multi-linear interpolation to approximate the value of $J_{k+1}(x_{k+1},\alpha)$. We set $J_{k+1}(x_{k+1},\alpha):=\frac{(x_{k+1}-x_{i})\cdot J_{k+1}(x^{i+1},\alpha)+(x^{i+1}-x_{k+1})\cdot J_{k+1}(x^{i},\alpha)}{x_{i+1}-x_{i}}$, where $x^{i}\in G_{s}$ and $x^{i+1}\in G_{s}$ are the two nearest grid points to $x_{k+1}$ that satisfy $x^{i}\leq x_{k+1}\leq x^{i+1}$.

Monte Carlo implementations. For each $(x,\alpha)\in G$, we sampled 100,000 trajectories starting from $x_{0}=x$, subject to keeping the valve open over time, as this is an optimal policy. For each trajectory sample $i$, we computed the cost sample $z_{i}:=\max_{k}\{g(x_{k}^{i})\}$, and estimated the Conditional Value-at-Risk of the 100,000 cost samples at the confidence level $\alpha$. We used the CVaR estimator $\widehat{\text{CVaR}}_{\alpha}[Z]:=\frac{1}{\alpha M}\sum_{i=1}^{M}z_{i}\textbf{1}_{\{z_{i}\geq\hat{Q}_{\alpha}\}}$, where $\hat{Q}_{\alpha}$ is the $(1-\alpha)$-quantile of the empirical distribution of the samples $\{z_{i}\}_{i=1}^{M}$, and $M:=100,000$ is the number of samples [11, Sec. 6.5.1]. Since this estimator is designed for continuous distributions, we added zero-mean Gaussian noise with a small standard deviation, $\sigma:=10^{-12}$, to each cost sample prior to computing the CVaR. Fig. 4 provides a Monte Carlo estimate of $W_{0}^{*}$. To obtain a Monte Carlo estimate of $J_{0}^{*}$, we used the same procedure but with the cost sample $z_{i}:=\xi_{i}+\sum_{k=0}^{N}\beta e^{m\cdot g(x_{k}^{i})}$, $\xi_{i}\sim\mathcal{N}(0,\sigma:=10^{-7})$, $m:=10$, and $\beta:=10^{-3}$; see Fig. 2.

VI Conclusion

In this paper, we propose the novel idea of a risk-sensitive safe set to encode safety of a stochastic dynamic system in terms of an allowable level of risk of constraint violations $r$ in the $\alpha$-fraction of the most damaging outcomes. We show how the computation of a risk-sensitive safe set can be reduced to the solution to a Markov Decision Process, where cost is assessed according to the Conditional Value-at-Risk measure. Further, we devise a tractable algorithm to approximate a risk-sensitive safe set, and provide theoretical and empirical arguments about its correctness.

Risk-sensitive safe sets have the potential to inform the design of safety-critical infrastructure systems, by revealing trade-offs between the risk of damaging outcomes and design choices at different levels of confidence. We illustrate our risk-sensitive reachability approach on a stormwater retention pond that must be designed to operate safely in the presence of uncertain rainfall. Our results reveal that the current design of the pond is likely undersized: even if the pond starts empty, there is a risk of at least 0.25ft of overflow at most levels of confidence under the random surface runoff of the design storm (see Fig. 3, $r=0.25$ plot at $x=0$).

On the methodological side, future steps are: 1) to formally prove the correctness of the value-iteration algorithm, 2) to devise approximate value-iteration algorithms to improve scalability, and 3) to consider a broader class of risk measures. On the applications side, future steps are: 1) to adjust the parameters of the dynamics model (e.g., outlet radius) to reduce the risk of extreme overflows, 2) to apply our method to a more realistic stormwater system that consists of two ponds in series on a larger grid, and 3) to develop an optimized toolbox for the computation of risk-sensitive safe sets. We are hopeful that with further development, the concept of risk-sensitive reachability will become a valuable tool for the design of safety-critical systems.

ACKNOWLEDGMENTS

We thank Sumeet Singh, Dr. Mo Chen, Dr. Murat Arcak, Dr. Alessandro Abate, and Dr. David Freyberg for discussions. MC is supported by an NSF Graduate Research Fellowship and was supported by a Berkeley Fellowship for Graduate Studies. This work is supported by NSF CPS 1740079 and NSF PIRE UNIV59732.

APPENDIX

Proof:

The proof relies on two facts. The first fact is,

[TABLE]

The second fact is that CVaR is a coherent risk measure, so it satisfies certain properties. CVaR is positively homogeneous, $\text{CVaR}_{\alpha}[\lambda Z]=\lambda\text{CVaR}_{\alpha}[Z]$ for any $\lambda\geq 0$, and monotonic, $\text{CVaR}_{\alpha}[Y]\leq\text{CVaR}_{\alpha}[Z]$ for any random variables $Y\leq Z$ [12, Sec. 2.2]. Also, CVaR can be expressed as the supremum expectation over a particular set of probability density functions [11, Eqs. 6.40 and 6.70]. Using this property and $\mathbb{E}[\log(Z)]\leq\log\left(\mathbb{E}[Z]\right)$, one can show,

[TABLE]

for any random variable $Z$ with finite expectation.

By monotonicity, positive homogeneity, (16), and (17),

[TABLE]

where $\bar{Y}_{x}^{\pi}:=Y_{x}^{\pi}/\beta$. Now, if $x\in\mathcal{U}_{\alpha}^{r}$, then

[TABLE]

since $\bar{\Pi}_{0}$ is included in $\Pi$. By Lemma 1, there exists $\pi\in\Pi$ such that

[TABLE]

where the second inequality holds by (18). So, $x\in\mathcal{S}_{\alpha}^{r}$. ∎

Theoretical Justification of Conjecture (C): Let $\epsilon>0$. For all $k=0,\dots,N-1$ and $z_{k}:=(x_{k},y_{k})\in\mathcal{X}\times\mathcal{Y}$, let $\mu_{k}^{\epsilon}:\mathcal{X}\times\mathcal{Y}\to U$ satisfy,

[TABLE]

Let $J_{k}^{\epsilon}$ be a sub-optimal cost-to-go starting at time $k$,

[TABLE]

where $\pi_{k}^{\epsilon}:=(\mu_{k}^{\epsilon},\dots,\mu_{N-1}^{\epsilon})=(\mu_{k}^{\epsilon},\pi_{k+1}^{\epsilon})$. Recall $J_{k}$, as defined in (13). Define J_{k}^{*}(z_{k}):=\min_{\pi\in\bar{\Pi}_{k}}\text{CVaR}_{y_{k}}\big{[}\textstyle\sum_{i=k}^{N}c(x_{i})\big{|}z_{k},\pi\big{]}. To prove the Conjecture, we would like to show by induction that for all $z_{k}:=(x_{k},y_{k})\in\mathcal{X}\times\mathcal{Y}$ and $k=N-1,\dots,0$,

[TABLE]

which is the proof technique in [4, Sec. 1.5]. One can show (21) for the base case, $k:=N-1$, since $J_{N}$ is known. Assuming (21) holds for index $k+1$ (induction hypothesis), we want to show that (21) holds for index $k$ (induction step). The key idea is to use the recursion,

[TABLE]

which we justify next. Let $Z:=\sum_{i=k+1}^{N}c(x_{i})$.

[TABLE]

where we use (11), (20), and translation invariance ($\text{CVaR}_{\alpha}[a+Z]=a+\text{CVaR}_{\alpha}[Z]$ for $a\in\mathbb{R}$, see [12, Sec. 2.2]). The last equality (a) is the crux of the Conjecture, as one needs to justify why the worst-case density $R$ is equal to the a priori chosen density $\bar{R}$ that defines the dynamics of the confidence level. Based on [15], we believe this equality to be correct, but we leave its formal proof for future research. Assuming the aforementioned equality is correct, then we show (21a) for index $k$ using (22) and the induction hypothesis. Let $\bar{\epsilon}_{k}:=(N-k-1)\epsilon$, and $x^{\prime}:=x_{k+1}$.

[TABLE]

since $\mathbb{E}[R]=1$, and by (19). By (13), sub-optimality of $\mu_{k}^{\epsilon}(z_{k})\in U$, $J_{k+1}\leq J^{\epsilon}_{k+1}$, and (22),

[TABLE]

which would complete the induction step for (21a), if (22) holds. Next, we show (21b) for index $k$. By definition, $J_{k}^{*}$ is the optimal risk-sensitive cost-to-go from stage $k$, thus, $J_{k}^{*}\leq J_{k}^{\epsilon}$. Let $\hat{\epsilon}_{k}:=(N-k)\epsilon$, $x^{\prime}:=x_{k+1}$, $y^{\prime}:=y_{k}R$, and $Z:=\sum_{i=k+1}^{N}c(x_{i})$. For any $\pi_{k}:=(\mu_{k},\pi^{\prime})\in\bar{\Pi}_{k}$,

[TABLE]

The above statement implies

[TABLE]

which completes the induction step for (21b).

Thus, if (22) holds, then (21a) and (21b) hold for index $k$ for any $\epsilon>0$. So, (21c) would hold for index $k$. Assuming the conjectured equality (a) is correct, this would complete the proof of the Conjecture.

AUTHOR CONTRIBUTIONS

MC: formulation of risk-sensitive safe set based on CVaR so that it generalized stochastic reachability, formulation and proof of reduction to CVaR-MDP theorem (Theorem 1), identified minimax error in Chow 2015 paper, theoretical justification for value iteration using techniques from [4] (see extended manuscript), implementation of numerical example, connected experts in controls, risk measures, stormwater catchment design together to create the paper; wrote initial draft and revised to incorporate co-authors’ feedback.

JL: careful in-depth study of the value iteration algorithm to identify why it is an approximation, important technical revisions (e.g., the intuition behind the confidence level, the reason for the approximate algorithm, clearly specifying the history-dependent policy space), proof of Lemma 1, statement of Lemma 3, expertise of Chow 2015 paper; reviewed drafts of the paper and proofs.

AT: explained Chow 2015 paper to MC during numerous discussions; reviewed drafts of the paper and proofs.

DL: integration of the parameter $m$ into the log-sum-exp approximation in proof of Theorem 1.

KS: expertise in PCSWMM, used PCSWMM to generate realistic simulations of the stormwater catchment in response to the design storm; interpreted the meaning of risk-sensitive safe sets in the context of the stormwater design application, which is a new application for reachability theory; reviewed drafts of the paper and proofs.

VC: generated the probability distribution for the surface runoff using empirical moments.

JF: provided initial golden nugget of “can we formulate a risk-sensitive version of reachability-based control?”; worked with MC on numerous occasions to study risk measures and measure theory; reviewed drafts of the paper and proofs.

SJ: MC’s research mentor and PI at SRI international over the summer when this paper was developed, through numerous discussions helped M.C. convey mathematical ideas more clearly, improved the proof of Lemma 2; reviewed drafts of the paper and proofs.

MP: JL’s adviser, made numerous important revisions to the framing of the paper to take into account the approximate value-iteration method; reviewed drafts of the paper and proofs.

CT: MC’s adviser, numerous discussions with M.C. during the development of the paper to improve accessibility of the ideas and techniques; reviewed drafts of the paper and proofs.