$\pi$. We reached a contradiction with the choice of $\sigma$.

For case (b), the only difference is that now the resource level after
$\alpha_{..j}\odot\beta$ can be higher then the one of
$\alpha_{..l}\odot\beta$ due to the removal of the intermediate non-reloading
cycle. Since we need to show that the energy level never drops below 0, the
same argument works. ∎


# Probabilistic Reachability and Invariance Computation of Stochastic Systems
using Linear Programming

Niklas Schmid  John Lygeros  Automatic Control Laboratory, ETH Zürich, 8092,
Zürich, Switzerland (e-mail: {nikschmid, jlygeros}@ethz.ch).

###### Abstract

We consider the safety evaluation of discrete time, stochastic systems over a
finite horizon. Therefore, we discuss and link probabilistic invariance with
reachability as well as reach-avoid problems. We show how to efficiently
compute these quantities using dynamic and linear programming.

###### keywords:

Stochastic Systems, Reachability Analysis, Linear Programming, Dynamic
Programming, Invariance, Viability, Safety, Optimal Control

††thanks: Work supported by the European Research Council under the Horizon
2020 Advanced under Grant 787845 (OCAL).

##  1 Introduction

The need for safety guarantees when controlling systems arises in many
important fields of application, e.g., in air-traffic management (Prandini and
Hu, 2008) and flight control (Lesser et al., 2013). In these applications,
safety is commonly defined by a certain set of safe states, in the sense that
any state trajectory leaving this set during the mission is considered unsafe.
Obviously, whether safety can be guaranteed depends on the initial state that
the system is in during the beginning of the mission. A natural question
arising is how to evaluate which set of initial states allows for a safe
evolution of the state trajectory of a controlled system. For stochastic
systems, safety might only be guaranteed up to some probability. In such cases
the following two questions might arise: What is the maximum / minimum
achievable probability that the state trajectory remains within a certain safe
set (Invariance)? What is the maximum / minimum achievable probability that
the state trajectory reaches a certain target set (Reachability)? Both
questions can be connected by interpreting a trajectory to be safe whenever it
passes through a target set at a desired point in time. Furthermore, one can
combine the two notions, leading to reach-avoid problems, where the goal is to
reach a target set while staying safe.

The computation of safety for stochastic and deterministic, continuous and
discrete time systems has been extensively studied in the literature (e.g., by
Abate et al. (2008), Mitchell et al. (2005), Liao et al. (2022), Esfahani et
al. (2016)). This paper specifically considers nonlinear, stochastic systems,
where closed form solutions for safety guarantees typically do not exist. In
such cases, despite their computational complexity, dynamic programming (DP)
methods have proven useful, see Abate et al. (2008).

Since these dynamic programming solutions are infinite dimensional,
approximation techniques using gridding, semidefinite programming (Drzajic et
al., 2017) and Lagrangian methods (Gleason et al., 2017) are applied. In the
case of reach-avoid problems, Kariotoglou et al. (2017) propose linear
programming formulations, which allow for the utilization of basis function
approximations. Gao et al. (2020) extend these results to maximum invariance
problems, but assume discrete states and actions. In this paper, we generalize
these linear programming formulations to minimum and maximum reachability and
invariance problems with continuous states and actions. Our contributions are
as follows.

  * •

First, we establish a link between probabilistic reachability and invariance.
This can be considered an extension to the observations by Lygeros (2004) and
Liao et al. (2022) for deterministic systems.

  * •

We then utilize simplified variants of established DP recursions (Abate et
al., 2007; Ding et al., 2013) to formulate linear programs (LP) for the
computation of reachable and invariant sets. Our results allow for a clear
perspective and intuition on probabilistic reachability and invariance and
their connection.

Other notable discussions on infinite-dimensional linear programming for the
invariant set computation are given by Korda et al. (2021) and Miller et al.
(2021), although they follow a slightly different approach.

###  1.1 Notation

We denote by $\mathbb{1}_{A}(x)$ the indicator function in the set $A$, where,
if $x\in A$, then $\mathbb{1}_{A}(x)=1$ and if $x\notin A$, then
$\mathbb{1}_{A}(x)=0$. For two sets $X,Y$ we denote by $X\setminus Y=\\{x\in
X:x\notin Y\\}$. $P$ denotes probability, $\mathbf{E}$ expectation and
$\mathcal{B}(X)$ the Borel $\sigma$-Algebra on a topological space $X$.

###  1.2 Structure

Section 2 acts as an introduction. We define what we mean by stochastic
systems in 2.1, safety in 2.2, invariance and reachability in 2.3. Section 3
discusses how to compute these quantities. We first state DP recursions in
3.1, then reformulate these recursions as an infinite dimensional linear
program in 3.2. In order to solve the infinite dimensional linear program, we
discuss approximation methods in 3.3 and apply them in numerical examples in
4. We summarize our results in 5.

##  2 Safety of stochastic Systems

###  2.1 Stochastic Systems and Policies

A discrete time stochastic system is described by the state space
$\mathcal{X}\subseteq\mathbb{R}^{n}$, a compact Borel set
$\mathcal{U}\subseteq\mathbb{R}^{m}$ denoting the action space, and a Borel-
measurable stochastic kernel
$T:\mathcal{X}\times\mathcal{X}\times\mathcal{U}\rightarrow[0,1]$, which,
given $x\in\mathcal{X},u\in\mathcal{U}$, assigns a probability measure
$T(\cdot|x,u)$ on the set $\mathcal{B}(\mathcal{X})$.

We will denote a state at time-step $k\in\mathbb{N}$ as $x_{k}$ and a sequence
of states $x_{k},\dots,x_{N}$ as $x_{k:N}$. The state evolves
probabilistically according to the transition kernel: For
$x_{k}\in\mathcal{X},u_{k}\in\mathcal{U}$ the state transitions to a state
$x_{k+1}$ in the set $B\subseteq\mathcal{B}(\mathcal{X})$ with probability
$T(B|x_{k},u_{k})$.

A Markov policy $\pi$ is a sequence $\pi=(\mu_{0},\mu_{1},\dots)$ of
universally measurable maps
$\mu_{k}:\mathcal{X}\rightarrow\mathcal{U},k=0,1,\dots$. Starting from an
initial state $x_{0}\in\mathcal{X}$ and under the policy $\pi$ the state
evolves as $x_{k+1}\sim T(\cdot|x_{k},\mu_{k}(x_{k}))$. We denote the set of
Markov policies by $\Pi$.

The transition kernel $T$, initial state $x_{0}\in\mathcal{X}$ and policy
$\pi\in\Pi$ define a unique probability measure over
$\mathcal{B}(\mathcal{X}^{N+1})$ for the trajectories (see, for example, Abate
et al. (2008)).

###  2.2 Safety of a Policy

Given an initial state $x_{0}\in\mathcal{X},\pi\in\Pi$ and a Borel set
$A\subseteq\mathcal{B}(\mathcal{X})$, which we will call the safe set, we wish
to compute the probability $p^{\pi}_{x_{0}}=P(x_{0:N}\in A|\pi,x_{0})$ that
the system trajectory remains within $A$ for $k=0,\dots,N$. By abuse of
notation, we interpret $x_{0:N}\in A$ to mean $x_{k}\in A$ for $k=0,\dots N$.
As shown by Abate et al. (2008), the safety of a trajectory can be encoded as

$\prod_{k=0}^{N}\mathbb{1}_{A}(x_{k})=\begin{cases}1&\text{if $x_{0:N}\in
A$,}\\\ 0&\text{otherwise}.\end{cases}$

Then, for a stochastic evolution of the state trajectory the probability of
safety is defined by

$\displaystyle p^{\pi}_{x_{0}}$
$\displaystyle=\int_{\mathcal{X}^{N}}\prod_{k=0}^{N}\mathbb{1}_{A}(x_{k})P(dx_{1},\dots,dx_{N}|x_{0},\pi)$
$\displaystyle=\underset{x_{1},\dots,x_{N}}{\mathbf{E}}\left[\prod_{k=0}^{N}\mathbb{1}_{A}(x_{k})\middle|x_{0},\pi\right],$

which can be computed through a dynamic programming recursion. Therefore, we
define functions $V_{k}^{\pi}:\mathcal{X}\rightarrow[0,1]$ to denote
$V_{k}^{\pi}(x_{k})=P(x_{k:N}\in A|x_{k},\pi)$, then at time step $N$ we
trivially obtain $V_{N}^{\pi}(x_{N})=\mathbb{1}_{A}(x_{N})$. Moreover,

$\displaystyle P(x_{k+1:N}\in A|x_{k},\pi)$
$\displaystyle\qquad=\int_{x_{k+1}\in\mathcal{X}}P(x_{k+1:N}\in
A|x_{k+1},x_{k},\pi)P(dx_{k+1}|x_{k},\pi)$
$\displaystyle\qquad=\int_{x_{k+1}\in\mathcal{X}}P(x_{k+1:N}\in
A|x_{k+1},\pi)P(dx_{k+1}|x_{k},\pi)$
$\displaystyle\qquad=\int_{x_{k+1}\in\mathcal{X}}V_{k+1}^{\pi}(x_{k+1})T(dx_{k+1}|x_{k},\mu_{k}(x_{k})),$

where the second equality follows from the Markov property of the system and
the last equality follows from the definitions of $V_{k}^{\pi}$ and $T$.
Consequently,

$\displaystyle V_{k}^{\pi}(x_{k})$ $\displaystyle=P(x_{k:N}\in A|x_{k},\pi)$
$\displaystyle=P(x_{k}\in A|x_{k},\pi)P(x_{k+1:N}\in A|x_{k},\pi)$
$\displaystyle=\mathbb{1}\\!_{A}(x_{k})\\!\int_{x_{k+1}\in\mathcal{X}}V_{k+1}^{\pi}(x_{k+1})T(dx_{k+1}|x_{k},\mu_{k}(x_{k})).$

Given a policy $\pi$ we can thus compute $V_{0:N}^{\pi}(x)$ recursively
backwards in time via

$\displaystyle V_{N}^{\pi}(x_{N})\\!$
$\displaystyle=\\!\mathbb{1}_{A}(x_{N}),$ $\displaystyle
V_{k}^{\pi}(x_{k})\\!$
$\displaystyle=\\!\mathbb{1}_{A}(x_{k})\\!\int_{x_{k+1}\in\mathcal{X}}\\!V_{k+1}^{\pi}(x_{k+1})T(dx_{k+1}|x_{k},\mu_{k}(x_{k})),$

which finally yields the desired probability $p_{x_{0}}^{\pi}=P(x_{0:N}\in
A|x_{0},\pi)=V_{0}^{\pi}(x_{0})$.

###  2.3 Invariance and Reachability

We now define the notion of invariance (probability that the trajectory
remains in the set $A$) and reachability (probability that the trajectory
enters the target set $A$ at least once).

######  Definition 2.1.

The probabilistic maximum invariant set [$I^{\uparrow}$], minimum invariant
set [$I_{\downarrow}$], maximum reachable set [$R^{\uparrow}$] and minimum
reachable set [$R_{\downarrow}$], are defined as

$\displaystyle\Omega^{A}_{I^{\uparrow},k}(p)\\!$
$\displaystyle=\\!\\{x_{k}\\!\in\\!\mathcal{X}|\exists\pi\\!\in\\!\Pi\\!:\\!P(x_{k:N}\\!\in\\!A|x_{k},\pi)\\!\geq\\!p\\},$
$\displaystyle\Omega^{A}_{I_{\downarrow},k}(p)\\!$
$\displaystyle=\\!\\{x_{k}\\!\in\\!\mathcal{X}|\forall\pi\\!\in\\!\Pi\\!:\\!P(x_{k:N}\\!\in\\!A|x_{k},\pi)\\!\geq\\!p\\},$
$\displaystyle\Omega^{A}_{R^{\uparrow},k}(p)\\!$
$\displaystyle=\\!\\{x_{k}\\!\in\\!\mathcal{X}|\exists\pi\\!\in\\!\Pi\\!:\\!P(\exists
x_{i}\\!\in\\!x_{k:N}\\!:\\!x_{i}\\!\in\\!A|x_{k},\\!\pi)\\!\geq\\!p\\},$
$\displaystyle\Omega^{A}_{R_{\downarrow},k}(p)\\!$
$\displaystyle=\\!\\{x_{k}\\!\in\\!\mathcal{X}|\forall\pi\\!\in\\!\Pi\\!:\\!P(\exists
x_{i}\\!\in\\!x_{k:N}\\!:\\!x_{i}\\!\in\\!A|x_{k},\\!\pi)\\!\geq\\!p\\}.$

This definition is in correspondence to the deterministic quantities by Liao
et al. (2022). Maximum invariant sets are often referred to as viability sets
or controlled invariant sets. Minimum invariant and maximum reachable sets are
often simply referred to as invariant and reachable sets.

Assuming that the maximum and minimum are attained - we will state sufficient
conditions later - it holds that

$\displaystyle\exists\pi\in\Pi:P(\cdot|\cdot,\pi)$ $\displaystyle\geq
p\Leftrightarrow\max_{\pi\in\Pi}P(\cdot|\cdot,\pi)\geq p,$ (1)
$\displaystyle\forall\pi\in\Pi:P(\cdot|\cdot,\pi)$ $\displaystyle\geq
p\Leftrightarrow\min_{\pi\in\Pi}P(\cdot|\cdot,\pi)\geq p,$ (2)

and the respective sets correspond to level sets of the functions
$V^{*}_{I^{\uparrow},k,A}(x_{k})=\max_{\pi\in\Pi}P(x_{k:N}\in A|x_{k},\pi)$,
$V^{*}_{R^{\uparrow},k,A}(x_{k})=\max_{\pi\in\Pi}P(\exists x_{i}\in
x_{k:N}:x_{i}\in A|x_{k},\pi)$,
$V^{*}_{I_{\downarrow},k,A}(x_{k})=\min_{\pi\in\Pi}P(x_{k:N}\in A|x_{k},\pi)$
and $V^{*}_{R_{\downarrow},k,A}(x_{k})\allowbreak=\min_{\pi\in\Pi}P(\exists
x_{i}\in x_{k:N}:x_{i}\in A|x_{k},\pi)$.

In the following sections we will discuss how to compute these functions using
linear and dynamic programming. We only consider invariance problems, since
reachability problems turn out to be duals to the respective invariance
problems, as shown in figure 1, and formalized in the following statement.

######  Theorem 2.2.

Let $A^{c}$ denote the complement of $A$. Then

$\displaystyle V^{*}_{R^{\uparrow},k,A}(x_{k})$
$\displaystyle=1-V^{*}_{I_{\downarrow},k,A^{c}}(x_{k}),$ $\displaystyle
V^{*}_{R_{\downarrow},k,A}(x_{k})$
$\displaystyle=1-V^{*}_{I^{\uparrow},k,A^{c}}(x_{k})$

###### Proof.

We start by showing the first relation. Note that $x_{k:N}\in
A\Leftrightarrow\nexists x_{i}\in x_{k:N}:x_{i}\notin A$. Thus,

$\displaystyle V^{*}_{R^{\uparrow},k,A}(x_{k})$
$\displaystyle=\max_{\pi\in\Pi}P(\exists x_{i}\in x_{k:N}:x_{i}\in
A|x_{k},\pi)$ $\displaystyle=\max_{\pi\in\Pi}1-P(\nexists x_{i}\in
x_{k:N}:x_{i}\in A|x_{k},\pi)$ $\displaystyle=1-\min_{\pi\in\Pi}P(\nexists
x_{i}\in x_{k:N}:x_{i}\in A|x_{k},\pi)$
$\displaystyle=1-\min_{\pi\in\Pi}P(\nexists x_{i}\in x_{k:N}:x_{i}\notin
A^{c}|x_{k},\pi)$ $\displaystyle=1-\min_{\pi\in\Pi}P(x_{k:N}\in
A^{c}|x_{k},\pi)$ $\displaystyle=1-V^{*}_{I_{\downarrow},k,A^{c}}(x_{k})$

which yields the desired equality and allows to compute

$\displaystyle\Omega^{A}_{R^{\uparrow},k}(p)$
$\displaystyle=\\{x_{k}\in\mathcal{X}|V^{*}_{R^{\uparrow},k,A}(x_{k})\geq
p\\},$
$\displaystyle=\\{x_{k}\in\mathcal{X}|1-V^{*}_{I_{\downarrow},k,A^{c}}(x_{k})\geq
p\\}.$

The proof for the second relation is similar and has already been shown in
Abate et al. (2008) ∎

Maximum Invariance$\max_{\pi}P(x_{k:N}\\!\in\\!A)$Minimum
Reachability$\min_{\pi}P(\exists
x_{i}\\!\in\\!x_{k:N}\\!:\\!x_{i}\\!\in\\!A)$Minimum
Invariance$\min_{\pi}P(x_{k:N}\\!\in\\!A)$Maximum
Reachability$\max_{\pi}P(\exists
x_{i}\\!\in\\!x_{k:N}\\!:\\!x_{i}\\!\in\\!A)$$V^{*}\leftrightarrow
1-V^{*}$$A\leftrightarrow
A^{c}$$\max\leftrightarrow$$\min$$V^{*}\leftrightarrow
1-V^{*}$$A\leftrightarrow A^{c}$$\max$$\leftrightarrow\min$

Figure 1: Links between invariance and reachability.

From now the dependency on the set $A$ will be omitted wherever possible to
simplify the notation.

##  3 Safety computation

We now show how to solve for maximum/minimum invariance using DP recursions.
Following Abate et al. (2008) we introduce an assumption to guarantee
existence of the solutions $V_{I^{\uparrow},k}^{*}(x_{k})$ and
$V_{I_{\downarrow},k}^{*}(x_{k})$.

######  Assumption 3.1.

The set

$\displaystyle\left\\{u_{k}\\!\in\\!\mathcal{U}\\!:\\!\int_{x_{k+1}\in\mathcal{X}}V_{I^{\uparrow},k+1}^{*}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})\geq\lambda\right\\}$

and the set

$\displaystyle\\!\left\\{u_{k}\\!\in\\!\mathcal{U}\\!:\\!\int_{x_{k+1}\in\mathcal{X}}V_{I_{\downarrow},k+1}^{*}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})\leq\lambda\right\\}$

are compact for all $x_{k}\in\mathcal{X},\lambda\in\mathbb{R},k\in[0,N-1]$.

If Assumption 3.1 holds, the maximum in equation (1) and minimum in equation
(2) are attained (Lemma 3.1 in Bertsekas and Shreve (1996)). A more intuitive,
but also more restrictive condition is the continuity of
$T(\cdot|\cdot,u_{k})$ with respect to $u_{k}$ (see Kariotoglou et al.
(2017)).

###  3.1 Computation via Dynamic Programming

######  Theorem 3.2.

Under Assumption 3.1, the maximum invariance problem can be solved by the DP
recursion (Abate et al. (2008))

$\displaystyle V_{I^{\uparrow},N}^{*}(x_{N})$
$\displaystyle\\!=\\!\mathbb{1}_{A}(x_{N}),$ $\displaystyle
V_{I^{\uparrow},k}^{*}(x_{k})$
$\displaystyle\\!=\\!\max_{u_{k}}\mathbb{1}_{A}(x_{k})\\!\int_{x_{k+1}\in\mathcal{X}}V_{I^{\uparrow},{k+1}}^{*}(x_{k+1})T(dx_{k+1}|x_{k},\\!u_{k}).$

The minimum invariance problem can be solved by the DP recursion

$\displaystyle V_{I_{\downarrow},N}^{*}(x_{N})\\!$
$\displaystyle=\\!\mathbb{1}_{A}(x_{N}),$ $\displaystyle
V_{I_{\downarrow},k}^{*}(x_{k})\\!$
$\displaystyle=\\!\min_{u_{k}}\mathbb{1}_{A}(x_{k})\\!\int_{x_{k+1}\in\mathcal{X}}V_{I_{\downarrow},{k+1}}^{*}(x_{k+1})T(dx_{k+1}|x_{k},\\!u_{k}).$

###### Proof.

The first statement is shown in Abate et al. (2008). The proof for the second
statement is similar. ∎

For the maximum invariance problem, we can ease the computation by the
following observation from Abate et al. (2007): For all $x_{k}\in
A^{c},\pi\in\Pi,V_{I^{\uparrow},k}^{\pi}(x_{k})=P(x_{k:N}\in A|x_{k},\pi)=0$.
Consequently, we only need to compute $V_{I^{\uparrow},k}^{*}(x_{k})$ for
every $x_{k}\in A$, for which

$\displaystyle V_{I^{\uparrow},k}^{*}(x_{k})\\!$
$\displaystyle=\\!\max_{u_{k}}\mathbb{1}_{A}(x_{k})\\!\int_{x_{k+1}\in\mathcal{X}}V_{I^{\uparrow},{k+1}}^{*}(x_{k+1})T(dx_{k+1}|x_{k},\\!u_{k})$
$\displaystyle=\\!\max_{u_{k}}\int_{x_{k+1}\in
A}V_{I^{\uparrow},k+1}^{*}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}),$

where the last equality follows since for any $x_{k+1}\in\mathcal{X}\setminus
A$ we have that $V_{I^{\uparrow},k+1}^{\pi}(x_{k+1})=0$. Thus restricting
integration to $x_{k+1}\in A$ does not change the value of the integral. It is
easy to see that a similar statement holds for the minimum invariance problem,
which simplifies to

$\displaystyle V_{I_{\downarrow},k}^{*}(x_{k})=\min_{u_{k}}\int_{x_{k+1}\in
A}V_{I_{\downarrow},{k+1}}^{*}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}).$

###  3.2 Computation via Linear Programming

An interesting theoretical property of the simplified recursions is that they
allow us to formulate the problems as infinite-dimensional LPs. Under
assumption 3.1, we relax seeking the maximum/minimum of the objective function
into seeking a mimum/maximum value bounded by the set of possible objective
values. For instance we can solve for the maximum invariance probability at
time $k$ as

$\displaystyle\min_{V_{k}}$ $\displaystyle\int_{x_{k}\in
A}V_{k}(x_{k})c(dx_{k})$ s.t. $\displaystyle V_{k}(x_{k})\geq\int_{x_{k+1}\in
A}V_{I^{\uparrow},k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})$

where $c(\cdot)$ is a nonnegative measure that assigns positive mass to all
open subsets of $\mathcal{X}$ and the constraints must hold for all $x_{k}\in
A,u_{k}\in\mathcal{U}$. A similar relaxation has been proposed in Kariotoglou
et al. (2017) to formulate reach-avoid problems as LPs. Concatenating the
constraints for different times $k$, we can solve invariance problems for all
time steps in a single linear program, leading to

$\displaystyle\min_{V_{0:N}}$ $\displaystyle\sum_{k=0}^{N-1}\int_{x_{k}\in
A}V_{k}(x_{k})c(dx_{k})$ s.t. $\displaystyle V_{k}(x_{k})\geq\int_{x_{k+1}\in
A}V_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})$ $\displaystyle V_{N}(x_{N})=1$

for the maximum invariance problem and

$\displaystyle\max_{V_{0:N}}$ $\displaystyle\sum_{k=0}^{N-1}\int_{x_{k}\in
A}V_{k}(x_{k})c(dx_{k})$ s.t. $\displaystyle V_{k}(x_{k})\leq\int_{x_{k+1}\in
A}V_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})$ $\displaystyle V_{N}(x_{N})=1$

for the minimum invariance problem, where the constraints must hold for all
$x_{k},x_{N}\in A,u_{k}\in\mathcal{U},k\in[0,N-1]$. We denote by
$V_{I^{\uparrow},0:N}^{LP}(x_{k})$ and $V_{I_{\downarrow},0:N}^{LP}(x_{k})$
the respective minimizers of the LPs.

Before we state our main results we need the following Lemma (see also Theorem
1a in Kariotoglou et al. (2017) for reach-avoid problems).

######  Lemma 3.3.

There is no feasible solution of the LP such that there exists $k\in[0,N]$ and
$x_{k}\in A$ such that
$V_{I^{\uparrow},k}^{LP}(x_{k})<V_{I^{\uparrow},k}^{*}(x_{k})$.

###### Proof.

Assume, for the sake of contradiction, that there exists a feasible solution
to the LP, where, for some $k\in[0,N-1]$,
$V_{I^{\uparrow},k+1}^{LP}(\cdot)=V_{I^{\uparrow},k+1}^{*}(\cdot)$ and there
exists some $x_{k}$ such that
$V_{I^{\uparrow},k}^{LP}(x_{k})<V_{I^{\uparrow},k}^{*}(x_{k})$. Let
$\pi^{*}=\\{\mu_{0}^{*},\dots,\mu_{N-1}^{*}\\}$ be the optimal policy obtained
from the DP recursion. Since the LP solution is feasible we have for all
$k\in[0,N-1],x_{k}\in A,u_{k}^{*}=\mu_{k}^{*}(x_{k})$

$\displaystyle V_{I^{\uparrow},k}^{LP}(x_{k})$
$\displaystyle\geq\int_{x_{k+1}\in
A}V_{I^{\uparrow},{k+1}}^{LP}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}^{*})$
$\displaystyle=\int_{x_{k+1}\in
A}V_{I^{\uparrow},k+1}^{*}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}^{*})$
$\displaystyle=V_{I^{\uparrow},k}^{*}(x_{k})$
$\displaystyle>V_{I^{\uparrow},k}^{LP}(x_{k})$

We obtain $V_{I^{\uparrow},k}^{LP}(x_{k})>V_{I^{\uparrow},k}^{LP}(x_{k})$,
which is a contradiction.

Applying prior result recursively starting from the given terminal value
$V_{I^{\uparrow},N}^{LP}(x_{N})\allowbreak=V^{*}_{I^{\uparrow},N}(x_{N})=1$
yields the claim for all $k\in[0,N]$. ∎

In fact, the opposite relation, that
$V_{I_{\downarrow},k}^{LP}(x_{k})>V_{I_{\downarrow},k}^{*}(x_{k})$ is
impossible, holds for the minimum invariance case. The proof is very similar
so we will skip it for brevity.

We are now ready to state our main result.

######  Theorem 3.4.

Under Assumption 3.1 the DP recursions for the viability and invariance
problem have a solution. Then so do the respective LPs; moreover the solutions
coincide up to a set of c-measure zero.

###### Proof.

For the sake of brevity we only prove the LP for the maximum invariance case.
The proof for the minimum invariance case is similar (see also Kariotoglou et
al. (2017) for reach-avoid problems).

First note that the DP solution is feasible for the LP, since

$\displaystyle V_{I^{\uparrow},k}^{*}(x_{k})=\int_{x_{k+1}\in
A}V_{I^{\uparrow},k+1}^{*}(x_{k+1})T(dx_{k+1}|x_{k},\mu^{*}_{k}(x_{k}))$

for all $k\in[0,N-1]$ and $V^{*}_{I^{\uparrow},N}(x_{N})=1$ for all $x_{N}\in
A$, which meets the constraints of the LP.

Recall that for all $k\in[0,N]$ and $x_{k}\in A$, any feasible solution yields
$V_{k}^{LP}(x_{k})\geq V_{k}^{*}(x_{k})$ by Lemma 3.3. Thus, for all feasible
solutions

$\displaystyle\sum_{k=0}^{N-1}\int_{x_{k}\in
A}V_{k}^{LP}(x_{k})c(dx_{k})\geq\sum_{k=0}^{N-1}\int_{x_{k}\in
A}V^{*}_{k}(x_{k})c(dx_{k}).$

Assume now that for some $k\in[0,N]$ there is a set $\\{x_{k}\in
A:V_{k}^{LP}(x_{k})>V_{k}^{*}(x_{k})\\}$ with non-zero c-measure. Then this
would consequently lead to a suboptimal objective value, as there is no
$V_{k}^{LP}(x_{k})<V_{k}^{*}(x_{k})$ to compensate by Lemma 3.3. Since
$V_{k}^{*}(x_{k})$ is feasible and it yields the lowest possible objective
value it is a globally optimal solution to the LP. ∎

Since the LP solution might differ from the DP solution for a set of states of
zero c-measure, the optimal solution to the LP does not have to be unique.
Moreover, if $V^{LP}_{I^{\uparrow},k}(x_{k})$ and
$V^{*}_{I^{\uparrow},k}(x_{k})$ are solutions to the LP and DP recursions,
respectively, then by optimality of the DP solution there does not exist
$\pi\in\Pi,k\in[0,N],x_{k}\in
A:V^{LP}_{I^{\uparrow},k}(x_{k})=V^{\pi}_{I^{\uparrow},k}(x_{k})>V^{*}_{I^{\uparrow},k}(x_{k})$.

Next, we consider the set of optimal policies as

$\Pi^{*}\\!=\\!\left\\{\pi\\!\in\\!\Pi\\!:\\!\sum_{k=0}^{N-1}\int_{x\in
A}V_{I^{\uparrow},k}^{\pi}(x)c(x)\\!=\\!\sum_{k=0}^{N-1}\int_{x\in
A}V_{I^{\uparrow},k}^{*}(x)c(x)\right\\}.$

One might be tempted to think that an optimal policy $\pi^{*}\in\Pi^{*}$ is
given by the active constraints of the LP (if they exist). While this is true
for countable state problems with a c-measure that is strictly positive on
every state, it is not always true for infinite state problems. As an example,
assume that $V_{0:N}^{LP}$ is an optimal solution to the LP. We denote by
$\mathbb{B}_{k}=\\{x_{k}\in A:V_{k}^{LP}(x_{k})>V^{*}_{k}(x_{k})\\}$, which
must have zero c-measure for the objective to be optimal. Interestingly, for
states $x_{k}\in A\setminus\mathbb{B}_{k}$ the optimal inputs
$u_{k}^{*}=\mu_{k}^{*}(x_{k})$ corresponding to the optimal policy obtained
from the DP recursion must be contained within the active constraints, since
for these inputs

$\displaystyle\int_{x_{k+1}\in A}V_{k+1}^{LP}$
$\displaystyle(x_{k+1})T(dx_{k+1}|x_{k},u_{k}^{*})={V}_{k}^{LP}(x_{k})$
$\displaystyle={V}^{*}_{k}(x_{k})=\int_{x_{k+1}\in
A}V^{*}_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}^{*}).$

Note that these inputs must consequently yield zero $T$-measure on
$\mathbb{B}_{k+1}$. However, there may also exist an input
$\tilde{u}_{k}\in\mathcal{U}$ corresponding to a tight constraint with

$\displaystyle\int_{x_{k+1}\in A}V^{*}_{k+1}$
$\displaystyle(x_{k+1})T(dx_{k+1}|x_{k},\tilde{u}_{k})<{V}^{*}_{k}(x_{k})$
$\displaystyle={V}_{k}^{LP}(x_{k})=\int_{x_{k+1}\in
A}V_{k+1}^{LP}(x_{k+1})T(dx_{k+1}|x_{k},\tilde{u}_{k}).$

This is possible since $\mathbb{B}_{k+1}$ can have non-zero $T$-measure.
Consequently, choosing a policy $\pi\in\Pi$ based on active constraints may
result in such suboptimal inputs $\tilde{u}_{k}=\mu_{k}(x_{k})$. If it does so
for a non-zero c-measure set of states $x_{k}\in A$, then the policy is
suboptimal, since the actual incurred invariance is

$\displaystyle{V}^{\pi}_{k}(x_{k})$ $\displaystyle=\int_{x_{k+1}\in
A}V^{\pi}_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},\tilde{u}_{k})$
$\displaystyle\leq\int_{x_{k+1}\in
A}V^{*}_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},\tilde{u}_{k})<{V}^{*}_{k}(x_{k}).$

For the sake of completeness we also provide the following results from
Kariotoglou et al. (2017) dealing with reach-avoid problems, where one tries
to reach a target set $\mathcal{T}\subseteq A$ while avoiding to step out of
the safe set $A$. This means that we want to maximize $P(\exists
x_{i}\\!\in\\!x_{k:N}\\!:\\!x_{i}\\!\in\\!\mathcal{T}\ \text{and}\ x_{k:i}\in
A|x_{k},\pi)$. This probability is one for every state $x_{k}\in\mathcal{T}$
and zero for every state $x_{k}\notin A$. For all $x_{k}\in
A\setminus\mathcal{T}$, the maximum reach-avoid probability is given by

$\displaystyle\min_{V_{0:N}}$ $\displaystyle\sum_{k=0}^{N-1}\int_{x\in
A}V_{k}(x)c(dx)$ s.t. $\displaystyle V_{k}(x_{k})\geq\int_{x_{k+1}\in
A\setminus\mathcal{T}}V_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k})$
$\displaystyle\qquad+\int_{x_{k+1}\in\mathcal{T}}T(dx_{k+1}|x_{k},u_{k})$
$\displaystyle V_{N}(x_{N})=0,$

where the constraints must hold for all $x_{k},x_{N}\in
A\setminus\mathcal{T},u_{k}\in\mathcal{U},k\in[0,N-1]$.

###  3.3 Approximation Techniques

Several techniques have been proposed to approximate the solution to the DP/LP
formulations. The state and input space may be divided into a finite number of
sets, each of which is represented by a discrete point; this effectively
amounts to approximating the value functions in the class of piecewise
constant functions and replacing integration by summation over a finite set
(Abate et al., 2008). The computation then reduces to the case of finite
states and actions. Moreover, in this class Assumption 3.1 is always fulfilled
and the DP and LP solutions always coincide.

More generally, the approach proposed in Kariotoglou et al. (2013)
approximates the value functions in a subspace spanned by a finite number of
basis functions, for example Radial Basis Functions (RBF). The optimization
variables of the LP then become parameters of the basis functions. By sampling
constraints one can obtain probabilistic guarantees based on the scenario
approach. To reconstruct the optimal policy at a given state one has to sample
inputs and then compute and compare the corresponding invariance or
reachability probabilities. This can by circumvented by computing the state-
action value functions $Q_{k}(x_{k},u_{k})$ in addition to $V_{k}(x_{k})$. The
linear program for the maximum invariance computation is then solved stagewise
for all $k\in[0,N-1]$ and reads as

$\displaystyle\min_{Q_{k}}$ $\displaystyle\int_{x_{k}\in
A,u_{k}\in\mathcal{U}}Q_{k}(x_{k},u_{k})c(dx_{k},du_{k})$ s.t. $\displaystyle
Q_{k}(x_{k},u_{k})\geq\int_{x_{k+1}\in
A}V_{k+1}(x_{k+1})T(dx_{k+1}|x_{k},u_{k}),$

where the constraints must hold for all $x_{k}\in A$ and
$u_{k}\in\mathcal{U}$, and we denote by
$V_{k}(x_{k})=\min_{u_{k}\in\mathcal{U}}Q_{k}(x_{k},u_{k})$, where
$V_{N}(x_{N})=1$ for all $x_{N}\in A$. The optimal policy can then be
recovered as the minimizing $u_{k}$ for the respective $Q_{k}$ at a given
state $x_{k}$. Again, there may be a zero c-measure set of states and inputs
for which the Q-function does not attain its feasible minimum.

An inherent problem of these approximations is that the computational
complexity of the LP grows exponentially with the state and input space
dimensionality.

##  4 Numerical example

As a numerical example we consider an autonomous robot that cannot rotate
while moving. Thus, the robot first rotates to a desired orientation $u_{k}$
at its current location $x_{k}=(x_{k}^{x},x_{k}^{y})^{\top}$ and afterwards
moves for a certain distance, which we fix at $3\text{\,}\mathrm{m}$ for
simplicity. However, due to measurement noise, the true orientation will be
$\theta_{k}=u_{k}+w_{k}$, with $w_{k}\sim\mathcal{N}(0,\pi/5)$, leading to

$\displaystyle\begin{bmatrix}x_{k+1}^{x}\\\ x_{k+1}^{y}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}x_{k}^{x}+3sin(\theta_{k})\\\
x_{k}^{y}+3cos(\theta_{k})\end{bmatrix}.$

unsafeunsafeunsafetarget Figure 2: The target set is marked yellow, edges of
unsafe sets white.

We define the state space as $\mathcal{X}=\mathbb{R}^{2}$, the action space as
any rotation $\mathcal{U}=(0,2\pi]$ and the safe set as a room of
$A=[0,50]\times[0,50]$\text{\,}{\mathrm{m}}^{2}$$ excluding some additional
interior walls. The robot must reach the top left corner of the room without
hitting any walls (see figure 2).

We compute the reach-avoid probability of the robot for $100$ time steps.
Therefore, we discretize the room into
$$1\text{\,}\mathrm{m}$\times$1\text{\,}\mathrm{m}$$ blocks and the action
space into $18$ actions. To compute the transition kernel, we simulate the
state transition at every state action pair $1000$ times with random samples
of $w_{k}$. In the second example we use gaussian radial basis functions
$\phi_{\epsilon,c}(x)=e^{-\epsilon^{2}||x-c||_{2}^{2}}$. We evenly place
$15\times 15$ basis functions across the state space and choose
$\epsilon=\frac{15}{50}$. To generate constraints for the LP, we randomly
sample $800$ states, $30$ inputs per state and $15$ samples of $w_{k}$ for
every state action pair. For simplicity, we compute the results sequentially
stage-wise instead of solving all time-steps in a single LP. The results are
shown in figure 3. The RBF-based solution shows artifacts visible as darker
gaps and tends to be more optimistic than the gridding approach, i.e., it
generally yields higher probability values. This is especially notable in the
time-step $k=N-100$, where the gridding approach computes a reduced reach-
avoid probability at the edges and in between the two lower walls.

DiscretizationRBFs$k=N-3$$k=N-20$$k=N-30$$k=N-40$$k=N-100$

Figure 3: The plots show reach-avoid probabilities, which are computed
backwards in time for $100$ time-steps using an RBF and discretization based
approach. The robot aims to reach the upper left corner while avoiding to hit
any walls. Initial states with high probabilities of achieving this goal until
the terminal time correspond to bright areas. The further back in time (i.e.,
the further right), the higher these probabilities.

To compare both approximation techniques further, the computation has been
carried out with a different number of basis functions and grid densities.
Grid100 to Grid25 denote results when using gridding with $100\times 100$ and
$25\times 25$ representative points, respectively. RBF20 to RBF5 denote
results when using $20\times 20$ to $5\times 5$ basis functions over the state
space. The value of $\epsilon$ and the number of state-samples has been varied
accordingly. Figure 4 shows the two norm distance of the respective results
evaluated at $100\times 100$ grid points to those of Grid100. With increasing
fineness of the grid and with increasing number of basis functions, the
results are expected to converge to the true reach-avoid probabilities.
Indeed, with increasing grid density and number of basis functions the plot
shows a decreasing difference to the results of Grid100. Interestingly, in our
simulation, the gridding based approach has been less sensitive to the
approximation density than the basis function based approach. In addition, the
computation times when using gridding have been significantly smaller (for
$20$ time steps: $90\text{\,}\mathrm{s}$ for Grid100, $4\text{\,}\mathrm{s}$
for Grid25, $27\text{\,}\mathrm{s}$ for RBF5, $5355\text{\,}\mathrm{s}$ for
RBF20 on a Surface Pro 8 i7). This is mainly due to the fact that the
transition kernel is computed once for the gridding based approach and stored
as a matrix at the cost of high memory usage, while the constraints in the
basis function based approach have been resampled at every time-step. However,
storing the full transition matrix is only feasible for low dimensional
systems since its size scales exponentially with the state dimensionality.

$-2$$0$$2$$4$$6$$8$$10$$12$$14$$16$$18$$20$$22$$0$$20$$40$$60$Timestep
$N-k$Two norm of distance to Grid100RBF5RBF10RBF15RBF20Grid25Grid50Grid75
Figure 4: We computed reach-avoid probabilities using gridding and basis
function approximations and varied the gridding density/number of basis
functions. The plot shows an increasing deviation over time between the
computed reach-avoid probabilities to those obtained using gridding with
$100\times 100$ representative points.

##  5 Conclusion

We established a link between probabilistic invariance and reachability for
stochastic systems and proposed infinite dimensional DP and LP formulations to
solve for these quantities. Approximate formulations have been evaluated in a
numerical example. We are aware that the LP formulations are not
computationally more efficient than classical DP formulations. However, they
may allow for different analysis and provide a new perspective on the problem.

## References

  * Abate et al. (2007) Abate, A., Amin, S., Prandini, M., Lygeros, J., and Sastry, S. (2007).  Computational approaches to reachability analysis of stochastic hybrid systems.  In _Hybrid Systems: Computation and Control_ , 4–17.  10.1007/978-3-540-71493-4_4.
  * Abate et al. (2008) Abate, A., Prandini, M., Lygeros, J., and Sastry, S. (2008).  Probabilistic reachability and safety for controlled discrete time stochastic hybrid systems.  _Automatica_ , 44, 2724–2734.  10.1016/j.automatica.2008.03.027.
  * Bertsekas and Shreve (1996) Bertsekas, D. and Shreve, S.E. (1996).  _Stochastic optimal control: the discrete-time case_ , volume 5.  Athena Scientific.
  * Ding et al. (2013) Ding, J., Kamgarpour, M., Summers, S., Abate, A., Lygeros, J., and Tomlin, C. (2013).  A stochastic games framework for verification and control of discrete time stochastic hybrid systems.  _Automatica_ , 49.  10.1016/j.automatica.2013.05.025.
  * Drzajic et al. (2017) Drzajic, D., Kariotoglou, N., Kamgarpour, M., and Lygeros, J. (2017).  A semidefinite programming approach to control synthesis for stochastic reach-avoid problems.  In _ARCH16. 3rd International Workshop on Applied Verification for Continuous and Hybrid Systems_ , volume 43 of _EPiC Series in Computing_ , 134–143.  10.29007/fqg6.
  * Esfahani et al. (2016) Esfahani, P.M., Chatterjee, D., and Lygeros, J. (2016).  The stochastic reach-avoid problem and set characterization for diffusions.  _Automatica_ , 70, 43–56.  10.1016/j.automatica.2016.03.016.
  * Gao et al. (2020) Gao, Y., Johansson, K., and Xie, L. (2020).  Computing probabilistic controlled invariant sets.  _IEEE Transactions on Automatic Control_ , PP, 1–1.  10.1109/TAC.2020.3018438.
  * Gleason et al. (2017) Gleason, J.D., Vinod, A.P., and Oishi, M.M.K. (2017).  Underapproximation of reach-avoid sets for discrete-time stochastic systems via lagrangian methods.  In _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_ , 4283–4290.  10.1109/CDC.2017.8264291.
  * Kariotoglou et al. (2017) Kariotoglou, N., Kamgarpour, M., Summers, T.H., and Lygeros, J. (2017).  The linear programming approach to reach-avoid problems for markov decision processes.  _Journal of Artificial Intelligence Research_ , 60(1), 263–285.  10.1613/jair.5500.
  * Kariotoglou et al. (2013) Kariotoglou, N., Summers, S., Summers, T., Kamgarpour, M., and Lygeros, J. (2013).  Approximate dynamic programming for stochastic reachability.  In _2013 European Control Conference (ECC)_ , 584–589.  10.23919/ECC.2013.6669603.
  * Korda et al. (2021) Korda, M., Henrion, D., and Mezić, I. (2021).  Convex computation of extremal invariant measures of nonlinear dynamical systems and markov processes.  _Journal of Nonlinear Science_ , 31, 1–26.  10.1007/s00332-020-09658-1.
  * Lesser et al. (2013) Lesser, K., Oishi, M., and Erwin, R.S. (2013).  Stochastic reachability for control of spacecraft relative motion.  In _52nd IEEE Conference on Decision and Control_ , 4705–4712.  10.1109/CDC.2013.6760626.
  * Liao et al. (2022) Liao, W., Liang, T., Wei, X., and Lai, J. (2022).  A novel unified framework for solving reachability and invariance problems.  _International Journal of Control_.  10.1080/00207179.2022.2051749.
  * Lygeros (2004) Lygeros, J. (2004).  On reachability and minimum cost optimal control.  _Automatica_ , 40, 917–927.  10.1016/j.automatica.2004.01.012.
  * Miller et al. (2021) Miller, J., Henrion, D., Sznaier, M., and Korda, M. (2021).  Peak estimation for uncertain and switched systems.  In _2021 60th IEEE Conference on Decision and Control (CDC)_ , 3222–3228.  10.1109/CDC45484.2021.9683778.
  * Mitchell et al. (2005) Mitchell, I.M., Bayen, A.M., and Tomlin, C.J. (2005).  A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games.  _IEEE Transactions on automatic control_ , 50(7), 947–957.  10.1109/TAC.2005.851439.
  * Prandini and Hu (2008) Prandini, M. and Hu, J. (2008).  Application of reachability analysis for stochastic hybrid systems to aircraft conflict prediction.  In _2008 47th IEEE Conference on Decision and Control_ , 4036–4041.  10.1109/CDC.2008.4739248.