Scalable Underapproximation for the Stochastic Reach-Avoid Problem for
High-Dimensional LTI Systems using Fourier Transforms

arXiv:1703.02135v2 [] 16 May 2017

Abraham P. Vinod and Meeko M. K. Oishi
Abstract— We present a scalable underapproximation of the
terminal hitting time stochastic reach-avoid probability at a
given initial condition, for verification of high-dimensional
stochastic LTI systems. While several approximation techniques
have been proposed to alleviate the curse of dimensionality
associated with dynamic programming, these techniques are
limited and cannot handle larger, more realistic systems. We
present a scalable method that uses Fourier transforms to compute an underapproximation of the reach-avoid probability for
systems with disturbances with arbitrary probability densities.
We characterize sufficient conditions for Borel-measurability
of the value functions. We exploit fixed control sequences
parameterized by the initial condition (an open-loop control
policy) to generate the underapproximation. For Gaussian
disturbances, the underapproximation can be obtained using
existing efficient algorithms by solving a convex optimization
problem. Our approach produces non-trivial lower bounds and
is demonstrated on a chain of integrators with 40 states.
Index Terms— Stochastic reachability; Stochastic optimal
control; Open-loop control; Convex optimization.

I. I NTRODUCTION
Reachability analysis of discrete-time stochastic dynamical systems is an established verification tool that provides
probabilistic guarantees of safety or performance, and has
been applied to problems in fishery management and mathematical finance [1], motion planning in robotics [2]–[4],
spacecraft docking [5], and autonomous survelliance [6].
In [1], two classes of problems characterize verification over
a finite horizon — first hitting time and terminal hitting time
– and dynamic programming approaches are formulated to
solve both (similarly to [7], [8]). We focus on the finite
time horizon terminal hitting time stochastic reach-avoid
problem (referred to here as the terminal time problem for
convenience), that is, computing the probability of hitting a
target set at the terminal time, while avoiding an unsafe set
during all the preceding time steps. Specifically, we construct
an underapproximation to the terminal time problem from
a known initial point, as opposed to the typical stochastic
reach-avoid problem. This could be used as a query, for
example, in evaluating feasibility of an initial trajectory an
optimization problem.
The dynamic programming-based discretization approach
(DPBDA), proposed in [8], approximately computes value
This material is based upon work supported by the National Science
Foundation under Grant Number IIS-1528047, CMMI-1254990 (Oishi,
CAREER), and CNS-1329878. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the authors and do
not necessarily reflect the views of the National Science Foundation.
Abraham Vinod and Meeko Oishi are with Electrical and Computer
Engineering, University of New Mexico, Albuquerque, NM; e-mail:
aby.vinod@gmail.com, oishi@unm.edu (corresponding author)

functions for the terminal time problem, but relies on gridding, and hence suffers from the well-known curse of dimensionality. Attempts to circumvent this problem, via approximate dynamic programming [9]–[11], Gaussian mixtures
[10], particle filters [5], [11], and convex chance-constrained
optimization [5], [6], have been applied to systems that are
at most 10-dimensional – far beyond the scope of what is
possible with DPBDA, but not scalable to larger problems.
In this paper, we first characterize sufficient conditions
for Borel-measurability of the value functions for the terminal time problem (characterized so far only for the first
hitting time problem [12]). Using conditional expectations,
we then establish that an open-loop formulation provides an
underapproximation of the stochastic reach-avoid probability for linear systems [5]. We propose a scalable Fourier
transform-based underapproximation (FTBU), for the terminal time problem, exploiting our prior work on uncontrolled
stochastic reachable sets [2]. For an arbitrary probability
density, the FTBU solves an optimization problem with
a multi-dimensional integration as the objective function.
For Gaussian disturbances, the objective function can be
computed efficiently via existing algorithms [13], and the
optimization problem is log-concave. Our approach does not
require gridding of state, input, or disturbance spaces, and
has low memory requirements in contrast to DPBDA.
Our main contribution is twofold: 1) a Fourier transformbased underapproximation of the terminal hitting time
stochastic reach-avoid probability from a known initial condition, based on open-loop control sequences, and 2) the
underlying theory that enables us to exploit measurability and
convexity properties to assure a computationally feasible approach. We extend our previous work on Fourier transformbased stochastic reachable sets for uncontrolled systems [2]
to systems with control inputs, although here we do not seek
to compute the stochastic reach-avoid set [14].
In Section II, we describe the terminal time problem,
its open-loop approximation, and relevant properties from
probability theory and Fourier analysis. Section III presents
sufficient conditions for Borel-measurability, and establishes
the underapproximation result linking the problems in [1]
and [5]. Section IV presents the FTBU and specialized results
for Gaussian disturbances. We demonstrate scalability in
Section V, through application to a 40D chain of integrators.
Section VI concludes the work.
II. P RELIMINARIES AND P ROBLEM F ORMULATION
We denote the Borel σ-algebra by B(·), a discrete-time
time interval by N[a,b] for a, b ∈ N and a ≤ b, which

inclusively enumerates all integers in between a and b,
random vectors with bold case, and non-random vectors with
an overline. The indicator function of a non-empty set S is
denoted by 1S (ȳ), such that 1S (ȳ) = 1 if ȳ ∈ S and is zero
otherwise. We denote the p-dimensional identity matrix by
Ip , and the matrix with all entries as ones by 1̄p×q ∈ Rp×q .
A. Probability theory
A random vector y is a measurable transformation defined
in the probability space (Ω, Y , P) with sample space Ω, σalgebra Y , and probability measure over Y , P. A sub-σalgebra of Y is a σ-algebra whose members also belong
to Y . The minimal σ-algebra of y, the smallest sub-σalgebra of Y over which y is measurable, is denoted by
σ(y) ⊂ Y . We typically consider Borel-measurable random
vectors, y : Rp → Rp with Ω = Rp and Y = σ(y) =
B(Rp ). For N ∈ N, a random process is a sequence
N
of random vectors {y k }k=0 where the random vectors y k
are defined in the probability space (Ω, Y , P). The random
>
vector Y = [y 0 y 1 . . . y N ] is defined in the probability
N
space (ΩN +1 , σ( k=0 Yk ), PY ), with PY induced from P.
See [15], [16] for details.
Conditional expectations transforms random variable p
in (Ω, P, P) whose mean exists (|E [p] | ≤ ∞) to a subσ-algebra Q, i.e, Q is a σ-algebra with all its members
containingRin P. ForR a Q-measurable random variable q
such that S pdP = S qdP for all S ∈ Q, q = E [p|Q]
almost surely (a.s.) [16, Sec. 7.1, Thm. 1].

×

P1) For Y -measurable random variables y 1 , y 2 with finite
means, a sub-σ-algebra G ⊂ Y , if y 1 ≤ y 2 a.s., then
E [y 1 |G ] ≤ E [y 2 |G ] a.s. [16, Sec. 7.1, eq. (14, i and
iii)].
P2) For bounded random variables y, z that are Y , Z measurable respectively, if |E [yz] | ≤ ∞, then
yE [z|Y ] = E [yz|Y ] a.s. [16, Sec. 7.1, Thm. 3].
P3) E [E [y|G2 ]|G1 ] = E [y|G1 ] a.s. if G1 ⊂ G2 ⊂ σ(y) and
E [y] ≤ ∞ [16, Sec. 7.1, eq. (14, v)].
The characteristic function (CF) of a random vector y ∈
R with probability density function (PDF) ψy (z̄) is


Ψy (ᾱ) , Ey exp j ᾱ> y
Z
>
=
ej ᾱ z̄ ψy (z̄)dz̄ = F {ψy (·)} (−ᾱ) (1)
p

Rp

where F {·} denotes the Fourier transformation operator and
ᾱ ∈ Rp . Given a CF Ψy (ᾱ), the PDF can be computed as
ψy (z̄) = F −1 {Ψy (·)} (−z̄)
 p Z
>
1
e−j ᾱ z̄ Ψy (ᾱ)dᾱ
=
2π
p
R

B. Terminal stochastic reach-avoid analysis
Consider the discrete-time stochastic LTI system,
xk+1 = Axk + B ūk + wk

(3)

with state xk ∈ X = Rn , input ūk ∈ U ⊆ Rm , disturbance
wk ∈ W ⊆ Rn , and matrices A, B assumed to be of
appropriate dimensions. We assume that U is compact, wk
is absolutely continuous with a known PDF ψw , and the
random process w[·] is independent and identical distributed
(IID). Let N be a finite time horizon. For any given sequence
of (non-random) inputs ū[·] and an initial condition x̄0 ∈ X ,
the state xk is a random vector for all k ∈ N[1,N ] via (3).
The system (3) can be equivalently described by a Markov
control process with stochastic kernel that is a Borelmeasurable function Q : B(X ) × X × U → [0, 1], which
assigns to each x̄ ∈ X and ū ∈ U a probability measure on
the Borel space (X , B(X )). For (3),
Q(dȳ|x̄, ū) = ψw (ȳ − Ax̄ − B ū)dȳ.

(4)

We define a Markov policy π = (µ0 , µ1 , . . . , µN −1 ) ∈ M
as a sequence of universally measurable maps µ[·] : X →
>
> >
U. The random vector X = [x>
1 x2 . . . xN ] , defined
x̄0 ,π
N
N
in (X , B(X ), PX ) [1], has probability measure Px̄X0 ,π
defined using Q [19, Prop. 7.45].
Let S, T ∈ B(X ). Define the terminal time probability,
r̂x̄π0 (S, T ), for known x̄0 and π, as the probability that the
execution with policy π is inside the target set T at time N
and stays within the safe set S for all time up to N . From [1],
x̄0 ,π 
r̂x̄π0 (S, T ) = PX
xN ∈ T ∧ xk ∈ S ∀k ∈ N[0,N −1] .
From [1, Def. 10], a Markov policy π ∗ is a maximal reachavoid policy in the terminal sense if and only if it is the
optimal solution of Problem A, defined as
∗

A : r̂x̄π0 (S, T ) = supπ∈M r̂x̄π0 (S, T )

(5)

The solution of Problem A is characterized via dynamic
programming [1, Thm. 11]. Define V̂k∗ : X → [0, 1], k ∈
N[0,N ] , by the backward recursion for x̄ ∈ X ,
V̂N∗ (x̄) = 1T (x̄)

(6)

V̂k∗ (x̄) = sup 1S (x̄)
ū∈U

Z

∗
V̂k+1
(ȳ)Q(dȳ|x̄, ū).

(7)

X
∗

Then, the optimal value to Problem A is r̂x̄π0 (S, T ) =
V̂0∗ (x̄0 ) for every x̄0 ∈ X .
Lemma 1. [1, Thm. 11] A sufficient condition for existence
of a maximal Markov policy for Problem A is
Z
∗
Uk (x̄, λ) = {ū ∈ U :
V̂k+1
(ȳ)Q(ȳ|x̄, ū)dȳ ≥ λ} (8)
X

(2)

and Uk is compact for all λ ∈ R, x̄ ∈ X and k ∈ N[0,N −1] .

where F −1 {·} denotes the inverse Fourier transformation
operator and dᾱ is short for dα1 dα2 . . . dαp . Since PDFs are
absolutely integrable, every PDF has a unique CF. See [17,
Sec. 1], [18, Sec. 22.6], [15, Sec. 7.2, 8.2], [2, Sec. 2.1] for
more details about CFs.

Lemma 1 assures universal measurability of V̂k∗ (·), and
that the Markov policy π ∗ consists of universally measurable
maps µ∗k [7, Thm. 1 proof]. However, evaluating (8) is
difficult. We propose alternative sufficient conditions, which
are easier to evaluate, and guarantee Borel-measurability
(stronger than universal measurability [19, Defn. 7.20]).

C. Open-loop stochastic reach-avoid analysis
>
>
>
With Ū = [ū>
∈
0 ū1 . . . ūN −1 ]
>
>
>
>
N
[w0 w1 . . . wN −1 ] ∈ W , we obtain

U

N

and W =

X = Āx̄0 + H̄ Ū + ḠW .

(9)

The matrices Ā, H̄, Ḡ are given by specific combinations of
the matrices A and B (see [20, Sec. 2]).
Consider an open-loop policy ρ : X → U N which provides an open-loop sequence of inputs ρ(x̄0 ) for every initial
condition x̄0 . Then X, defined in (9) under the action of
x̄ ,ρ(x̄ )
ρ(x̄0 ), lies in the probability space (X N , B(X N ), PX0 0 ),
x̄0 ,ρ(x̄0 )
with PX
defined using Q [19, Prop. 7.45]. Note that
ρ(x̄0 ) 6∈ M, since universally measurable maps µk (·) are
functions of xk , not x̄0 . Consequently, a Markov policy with
µk (·) as constants is a special case of ρ(·).
In [5], the authors approximate Problem A, without establishing the direction of approximation, with Problem B,
maximize
B:

subject to

ρ(x̄0 )

r̂
x̄0

(S, T )
x̄ ,ρ(x̄ )
X ∼ PX0 0
ρ(x̄0 ) ∈ U N

with decision variable ρ(x̄0 ), and
x̄ ,ρ(x̄ ) 
ρ(x̄ )
r̂x̄0 0 (S, T ) = PX0 0 xN ∈ T ∧ xk ∈ S ∀k ∈ N[0,N −1]
The optimal solution to Problem B is ρ∗ (x̄0 ). Since ρ(x̄0 ) 6∈
M, the relation between Problems A and B, apart from structural similarity, is not evident. Problem B was solved in [5]
approximately via particle filter and chance-constrained optimization methods.
We first demonstrate that Problem B underapproximates
Problem A, then use a Fourier transform-based approach that
enables an exact solution to Problem B.
Problem 1. Characterize the sufficient conditions under
which V̂k∗ (·) and µ∗k (·) are Borel-measurable for the terminal
time problem.
Problem 2. Show that the terminal time problem (Problem A) is underapproximated by the open-loop formulation
(Problem B).
Problem 3. a) Construct a scalable method for solving
Problem B by characterizing the forward stochastic reach
probability density for stochastic linear systems controlled
by ρ(·) when w has an arbitrary PDF. Additionally, b)
formulate Problem B as a convex optimization problem when
w is Gaussian.
III. T HEORETICAL RESULTS
A. Sufficient conditions for Borel-measurability of V̂k∗ (·)
Definition 1. [19, Defn. 7.12] A stochastic kernel Q(·|x̄, ū)
is continuous if for every (x̄, ū) ∈ X ×U and every sequence
i→∞
(x̄i , ūi ) −−−→ (x̄, ū),
lim Q(dȳ|(x̄i , ūi )) = Q(dȳ|(x̄, ū)).

i→∞

Lemma 2 follows from the fact that continuity is preserved
by composition [21, Cor. 13.1.7]. We have the following
theorem, similar to [12, Prop. 3].
Theorem 1. If U is compact and Q(·|x̄, ū) is continuous,
then V̂k∗ (·) are Borel-measurable functions for k ∈ N[0,N ]
and π ∗ , comprised of Borel-measurable maps µ∗k (·), exists.
Proof: (By induction) Since S, T are Borel sets, 1S (·)
and 1T (·) are Borel-measurable functions, and the result for
k = N follows trivially. Consider the base case k = N − 1.
Since V̂N∗ (·) is a bounded Borel-measurable function and
Q(·|x̄, ū)R is a Borel-measurable function, continuous over
X × U, X V̂N∗ (ȳ)Q(dȳ|x̄, ū) is continuous over X × U [22,
Fact 3.9]. Since continuity implies upper semi-continuity [19,
Lem. 7.13 (b)] and Borel-measurablity [16, Sec. 1.4], and U
is compact, an optimal Borel-measurable input map µ∗N −1 (·)
R
exists and X V̂N∗ (ȳ)Q(dȳ|x̄, µ∗N −1 (x̄)) is Borel-measurable
over X [23, Thm. 2]. Finally, V̂N∗ −1 (·) is Borel-measurable
since the product operator preserves Borel-measurability [21,
Cor. 18.5.6]. For the case k = t, assume for induction
∗
that V̂t+1
(·) is Borel-measurable. By the same arguments as
above, a Borel-measurable µ∗t (·) exists and V̂t∗ (·) is Borelmeasurable, completing the proof.
. Theorem 1 addresses Problem 1. Since Borelmeasurability implies universal measurability [19, Defn.
7.20], the hypotheses of Theorem 1 is stricter than Lemma 1,
but can be easily checked, and implies that V̂k∗ (xk ) is a
B([0, 1])-measurable random variable ∀k ∈ N[0,N ] . The
continuity requirements in Theorem 1 and Lemma 2 may
be weakened to include exponential densities [19, Sec. 8.3].
B. Problem B underapproximates Problem A
Next, we address Problem 2. For xk = x̄, denote the
expectation defined by Q(·|x̄, ū) as Eūx . Under the conditions
∗
proposed by Theorem 1, we know that V̂k+1
(xk+1 ) is a
Borel-measurable random variable for all k ∈ N[0,N ] . From
(7), for any k ∈ N[0,N −1] , we have almost surely (a.s.)1
h
i
∗
V̂k∗ (x̄) = sup 1S (x̄)Eūx V̂k+1
(xk+1 ) xk = x̄ .
(11)
ū∈U

Using Theorem 1 and properties of conditional expectations,
we can show the following theorem. See [24] for the proof.
Note that the state xk is not an independent random vector,
but part of a Markov control process controlled by a sequence
of actions. Therefore, the conditional expectation in (11) is
defined on the σ-algebra σ(xk ) = σ(x0 , x1 , . . . , xk ). For
any k ∈ N[0,N −1] , we have from (11)
i
h
∗
V̂k∗ (xk ) = sup 1S (xk )Eūx V̂k+1
(xk+1 ) σ(xk ) a.s. . (12)
ū∈U

Also, by the definition of a stochastic process [16, Sec. 5.3],
σ(xk ) ⊂ σ(xk+1 ) ∀k ∈ N[0,N −1] .

(13)

(10)

Lemma 2. If the PDF of the disturbance w ψw is continuous, then Q(·|x̄, ū) defined in (4) is continuous.

1 The a.s. equality arises because the conditional expectation of
∗ (x
V̂k+1
k+1 ) is defined only within an equivalence (can differ in sets of
zero probability measure) [16, Ch. 7].

Theorem
2. If U is compact and Q(·|x̄, ū) is continuous,
∗
ρ∗ (x̄ )
then r̂x̄0 0 (S, T ) ≤ r̂x̄π0 (S, T ) a.s. in x̄0 ∈ X .
Proof: For notational brevity, given a ∈ N[0,N ] , b ∈
4
Uab =

N[0,N −1] , we define
(ūa , ūa+1 , . . . , ūb ) for some U ∈
N −1
U N with U0N −1 = U and UN
as empty. We will later use
a similar definition for X.
For k ∈ N[0,N −1] , define Ŵk : X × U N −k → [0, 1] based
on (11), (a.s.)
h
i
N −1
Ŵk (x̄, UkN −1 ) , 1S (x̄)Eūxk Ŵk+1 (xk+1 , Uk+1
) xk = x̄

From Property P3 and (13), we have (a.s.)
" 1
!
#
Y
U01
N −1
Ŵ0 (x0 , U ) = EX 2
1S (xk ) Ŵ2 (x2 , U2 ) σ(x0 )
1

k=0

We repeatedly expand Ŵk (·) for k ∈ N[2,N ] and apply the
arguments presented above to obtain
" N −1
!
#
Y
U0N −1
Ŵ0 (x0 , U ) = EX N
1S (xk ) 1T (xN ) σ(x0 ) a.s.
1

k=0
x̄ ,ρ(x̄0 )

By definition of PX0

N −1
ŴN (x̄, UN
) , 1T (x̄).

N −1
, (X N
= U , ρ(x̄0 ))
1 = X, U0
ρ(x̄0 )

We see that Ŵk (·) are Borel measurable by a straightforward proof by induction. Next, we prove:
S1:
S2:

Ŵk (x̄, UkN −1 )
supρ(x̄0 )∈U N

Ŵ0 (x̄0 , ρ(x̄0 )) = r̂x̄0

(S, T ) a.s.,

and the definition of ρ∗ (x̄0 ) completes the proof of S2.
We denote the optimal value of Problem B as Ŵ0∗ (x̄0 ).

≤ V̂k∗ (x̄) a.s. in x̄ ∈ X , ∀k ∈ N[0,N −1]
ρ∗ (x̄ )
Ŵ0 (x̄0 , ρ(x̄0 )) = r̂x̄0 0 (S, T ) a.s. in x̄0 . IV. U NDER - APPROXIMATION VIA F OURIER TRANSFORMS

S1 implies supρ(x̄0 )∈U N Ŵ0 (x̄0 , ρ(x̄0 )) ≤ V̂0∗ (x̄0 ) a.s. in
x̄0 ∈ X . Thus, the proof is complete via S2 and the fact
∗
that V̂0∗ (x̄0 ) = r̂x̄π0 (S, T ) for every x̄0 ∈ X .
Proof of S1: (By induction) The base case is k = N − 1.
N −1
Since UN
−1 = ūN −1 ,
ūN −1
ŴN −1 (x̄, ūN −1 ) , 1S (x̄)Ex
[1T (xN )|xN −1 = x̄] a.s. .

From (11) and (6), we have (a.s.)
ŴN −1 (x̄, ūN −1 ) ≤

sup ŴN −1 (x̄, ūN −1 ) =
ūN −1 ∈U

ρ(x̄0 )

A. FTBU using an analytical expression for r̂x̄0

(S, T )

Let the PDF of the random vector X parameterized by the
initial condition x̄0 and the input vector Ū be ψX (X̄; x̄0 , Ū ).
ρ(x̄ )
The objective of Problem B r̂x̄0 0 (S, T ) is
Z Z
Z
ρ(x̄ )
ψX (X̄; x̄0 , ρ(x̄0 ))dX̄ (14)
r̂x̄0 0 (S, T ) =
...
T S
| {z S}
N − 1 times

V̂N∗ −1 (x̄).

For the case k = t, assume for induction that
N −1
∗
Ŵt+1 (x̄, Ut+1
) ≤ V̂t+1
(x̄) a.s. in x̄ ∈ X . We have to show
N −1
∗
Ŵt (ȳ, Ut
) ≤ V̂t (ȳ) a.s. in ȳ ∈ X . By Property P1 and
(12), we have (a.s.)

>

N
>
>
where X̄ = [x̄>
1 x̄2 . . . x̄N ] ∈ X , x̄k ∈ X ∀k ∈ N[1,N ] ,
and dX̄ is short for dx̄1 dx̄2 . . . dx̄N . Therefore, if ψX is
ρ(x̄ )
known, then r̂x̄0 0 (S, T ) is a nN -dimensional integral of
a PDF ψX over S × S × . . . × T . Determining ψX for a
known Ū can be posed as a forward stochastic reachability
problem using the CF of W [2, Prop. P3] defined as

N
−1
Ŵt (ȳ, UtN −1 ) ≤ sup Ŵt (ȳ, UtN −1 )
Y
ūt ∈U
ΨW (ᾱ) =
Ψw (ᾱk )
(15)
h
i
N −1
k=0
) xt = ȳ
≤ sup 1S (ȳ)Eūxt Ŵt+1 (xt+1 , Ut+1
ūt ∈U
>
i
h
(nN )
>
and ᾱk ∈ Rn for
where ᾱ = [ᾱ0> ᾱ1> . . . ᾱN
−1 ] ∈ R
∗
≤ sup 1S (ȳ)Eūxt V̂t+1
(xt+1 ) xt = ȳ
all k ∈ N[0,N −1] . We compute ψX via Proposition 1.
ūt ∈U

= V̂t∗ (ȳ).
This completes the proof of S1.
Proof of S2: We use the definition of Ŵk (·) based on
σ-algebra (similar to (12)). By (13), 1S (xk ) is a σ(xt )measurable random variable for t ∈ N[k,N ] . Expanding
Ŵ0 (·) and Ŵ1 (·) and using Property P2, we have (a.s.)
h
i
Ŵ0 (x0 , U ) = 1S (x0 )Eūx0 Ŵ1 (x1 , U1N −1 ) σ(x0 )
h
i
= Eūx0 1S (x0 )Ŵ1 (x1 , U1N −1 ) σ(x0 )
h
= Eūx0 1S (x0 )1S (x1 )
i
i
h
Eūx1 Ŵ2 (x2 , U2N −1 ) σ(x1 ) σ(x0 )
h
h
= Eūx0 Eūx1 1S (x0 )1S (x1 )
i
i
Ŵ2 (x2 , U2N −1 ) σ(x1 ) σ(x0 )

Proposition 1. For initial state x̄0 ∈ X , dynamics as in (9),
and open-loop control vector Ū , the PDF and CF of X are
ΨX (β̄; x̄0 , Ū ) = exp(j β̄ > (Āx̄0 + H̄ Ū ))ΨW (Ḡ> β̄) (16)

ψX (X̄; x̄0 , Ū ) = F −1 ΨX (β̄; x̄0 , Ū ) (−X̄)
(17)
>

>
where β̄ = [β̄1> β̄2> . . . β̄N
]
all k ∈ N[1,N ] .

∈ R(nN ) and β̄k ∈ Rn for

Proof: From (9), (15), and [2, Property P2].
In general, (17) is a nN -dimensional integration (2).
However, when the CF of X is in a standard form, a closedform expression for ψX (·) can be obtained, and we compute
ρ(x̄ )
ρ(x̄ )
r̂x̄0 0 (S, T ) via (14). Else, we can compute r̂x̄0 0 (S, T )
using ΨX if hthe Fourier
transform of 1S (·) and 1T (·) is
i
known and E X > X < ∞ [2, Sec. 4.2]. Thus, we can solve
Problem B, and thereby Problem 3a, using Proposition 1 and
(14) for arbitrary ψw .

Note that while scalability of this approach is contingent
on high-dimensional quadrature, this challenge is far more
tractable than the computational and memory costs associated with DPBDA. In general, we can compute (14) for
arbitrary disturbance densities through Monte-Carlo simulations [25, Sec. 4.8], [26, Ch. 4.2.1] and quasi-Monte Carlo
simulations [26, Ch. 4.2.2].
B. Gaussian disturbance
When w is a Gaussian random vector, the CF of w ∼
N (m̄, Σ) [15, Sec. 9.3] is


ᾱ> Σᾱ
>
.
(18)
Ψw (ᾱ) = exp j ᾱ m̄ −
2
Using (15), (18), and Proposition 1, ψX (·) is described by

X ∼ N m̄X , ΣX
(19a)
m̄X = Ḡ(1̄N ×1 ⊗ m̄) + Āx̄0 + H̄ Ū
>

ΣX = Ḡ(IN ⊗ Σ)Ḡ .

(19b)
(19c)

Proposition 2. For convex U, S, and T , dynamics as in (3)
and a Gaussian disturbance w, Problem B is log-concave.
Proof: From [27, Sec. 2.3], y ∼ N (0, ΣX ) is logconcave with respect to y. By [28, Sec. 3.2.2], ψX (19a) is
log-concave in Ū since it is an affine transformation of ψy by
y−m̄X . From [28, Sec. 2.3.2, Sec. 3.5.2], sets T ×S ×. . .×S
ρ(x̄ )
and U N are convex, and r̂x̄0 0 (S, T ) is log-concave over
N
U . Thus, Problem B is log-concave.
Proposition 2 addresses Problem 3b. For stochastic linear
systems with a Gaussian disturbance and polytopic S and
T , (14) is the integration of a multivariate Gaussian random
variable over a polytope. Efficient computation of (14) and
log-concavity (Proposition 2) enables a scalable solution to
Problem B when w is Gaussian.
C. FTBU implementation for the Gaussian disturbance case
To solve (14) when w is Gaussian, we use Genz’s algorithm [29], which is based on quasi-Monte-Carlo simulations
and Cholesky decomposition [13]. Genz’s algorithm provides
an error estimate that is the result of a trade-off between
accuracy and computation time. We set the number of
particles for the Monte-Carlo simulation so that the error
estimate is less than some  > 0. This results in a runtime
ρ(x̄ )
evaluation of r̂x̄0 0 (S, T ) that is dependent on x̄0 , unlike
typical Monte-Carlo simulations. To take the logarithm of
ρ(x̄ )
ρ(x̄ )
r̂x̄0 0 (S, T ) in Proposition 2, we set r̂x̄0 0 (S, T ) =  if
ρ(x̄0 )
r̂x̄0 (S, T ) < .
While the convexity result in Proposition 2 ensures a
tractable, globally optimal solution to Problem B, the lack
of a closed-form expression for the objective (14) requires
black-box optimization techniques. Further, since Genz’s
algorithm enforces an accuracy of only , the log-concavity
ρ(x̄ )
of r̂x̄0 0 (S, T ) may not be preserved. Hence the ideal solver
for Problem B should handle the “noisy” evaluation of (14)
as an oracle, and solve a constrained optimization problem.

We use MATLAB’s patternsearch to solve Problem B,
because it is based on direct search optimization [30] and
can handle estimation errors in (14) efficiently. The solver
is a derivative-free optimizer and uses evaluations over an
adaptive mesh to obtain feasible descents towards the globally optimal solution. However, it requires a larger number
of function evaluations as compared to fmincon. For linearlyconstrained and bound-constrained optimization problems
(such as Problem B, which is linearly constrained when U
is a polytope), creating the mesh using generating set search
reduces the number of function evaluations [30, Sec. 8].
D. Advantages and limitations of FTBU
The main advantage of FTBU is that it does not require
gridding of the state, input, or disturbance spaces. Unlike
the DPBDA [8], which solves Problem A on a grid over S,
irrespective of the size of the initial set of interest, the FTBU
solves Problem B at a desired x̄0 . By converting the terminal
time problem into an optimization problem involving a multidimensional integral, FTBU achieves higher computational
speed at lower memory cost (Figure 1) for a given initial
condition. Probabilistically verifying a set of initial conditions would require performing FTBU over a grid on the state
space (Figure 2), thereby losing any computational advantage
over DPBDA. An alternative approach for the verification
problem relies on Lagrangian methods [14].
While evaluating (14) can be computationally expensive
for arbitrary disturbances, for Gaussian disturbances we can
compute (14) efficiently (see Section IV-C). Further, since the
dimension of the integral in (14) is nN , large n effectively
limits the time horizon N . Additionally, the lack of feedback
in ρ(·) implies N cannot be large [5], as it may induce
excessive conservatism in the underapproximation.
V. N UMERICAL EXAMPLE
Consider the discrete-time chain of integrators, with state
xk ∈ Rn , input uk ∈ [−1, 1], a Gaussian disturbance wk ∼
N (0, 0.01In ), sampling time Ns = 10, and time horizon
N = 10.


1
Nsn−1
1 Ns 21 Ns2 . . . (n−1)!

 0 1
Ns



 .
.
..
xk+1 =  ..
..
 xk
.



 0 0
0
...
Ns
+

0


0
1
n
n! Ns

0

...

...

1
2 Ns

Ns

1
>

uk + wk

(20)

All computations were performed using MATLAB on an
Intel Core i7 CPU with 3.4GHz clock rate and 16 GB
RAM. MATLAB code for this work is available at http:
//hscl.unm.edu/files/code/LCSS17.zip.
A. Comparison of FTBU and DPBDA runtimes and bounds
We first demonstrate 1) scalability of the underapproximation as compared to the DPBDA, and 2) non-trivial lower
bounds obtained using FTBU. Figure 1 shows how FTBU
and DPBDA scale with state dimension n, for n ≤ 40.

Grid spacing

0.1

0.05

0.01

0.005

V̂0∗ (x̄1 )
0.422 0.506 0.476
0.478
V̂0∗ (x̄2 )
0.527 0.510 0.483
0.479
Computation time (seconds) 11.22 42.68 1206.59 5710.06

TABLE II: Grid spacing in DPBDA (n = 2, x̄1 = −x̄2 =
>
2
2
[0.1 0.9] , T = [−0.5, 0.5] , S = [−1, 1] , Ŵ0∗ (x̄1 ) =
∗
Ŵ0 (x̄2 ) = 0.436,  = 0.001).

Fig. 1: Scalability of DPBDA and FTBU with state dimension n to compute V̂0∗ (x̄0 ) and Ŵ0∗ (x̄0 ) for some x̄0 ∈
X . Average computation time for the FTBU based on 20
randomly chosen points in T at each n.
Initial state of
Runtime (s)
Ŵ0∗ (x̄0 )
Member
V̂0∗ (x̄0 ) fm
interest x̄>
fm
ps
ps
of
0
[0 0 0 . . . 0]
T
1
1
[0.999, 1] 12 302
[2.5 2.5 2.5 . . .]
T
0.984 0.986 [0.985, 1] 798 1196
[-8.5 8 -8.5 8 . . .] S \ T
0.500 0.999 [0.998, 1] 12 441

TABLE I: Non-trivial bounds for V̂0∗ (x̄0 ) (x̄0 ∈ R40 ); fm and
ps is FTBU with fmincon and patternsearch respectively.

S, T , we obtain non-trivial lower bounds even for highdimensional systems. Figure 2(f) and Table I show that
FTBU with patternsearch clearly outperforms fmincon in
the quality of the underapproximation, at the expense of
computation time.
Lastly, note that the expected symmetry about the origin of the terminal time probability for the system (20)
is not evident, unless a fine grid is used (Figure 2(a)).
Table II shows that FTBU can serve as a “certificate” for
the validity of the grid spacing in DPBDA by relying on
the conservativeness established by Theorem 2. That is,
the FTBU underapproximation provides a grid-independent
lower bound on the value function V̂0∗ (x̄0 ) computed using
DBPDA. For example, for the double integrator, a grid
spacing of 0.1 will not give accurate results with DPBDA,
since V̂0∗ (x̄1 ) < Ŵ0∗ (x̄1 ) contradicts Theorem 2.

n

We solve Problems A and B with S = [−10, 10] , T =
n
[−5, 5] , and  = 0.01. For DPBDA, we restrict the grid
over X to S for n ≤ 3. We approximate the disturbance
space as [−0.5, 0.5]n , based on the covariance matrix of wk .
We discretize X , W, and U with grid spacings of 0.05, 0.05,
and 0.1, respectively. As expected, FTBU implemented using
patternsearch is slower than fmincon, but both implementations scale with dimension n much better than DPBDA.
Table I summarizes the bounds on V̂0∗ (·) for various x̄0 at
n = 40. For high n, (20) becomes severely under-actuated
and the influence of the disturbance becomes very strong.
This leads ∗to the open-loop formulation yielding trivial lower
ρ (x̄ )
bounds, r̂x̄0 0 (S, T ) = , for many x̄0 in the original S, T .
40
40
We therefore set T = [−8, 8] and S = [−10, 10] . While
Theorem 2 assures that Ŵ0∗ (·) is a lower bound on V̂0∗ (·),
this bound is subject to , hence the discrepancies between
the numerical values for Ŵ0∗ (·) and the lower bounds on
V̂0∗ (·). We use  = 0.001.

VI. C ONCLUSION
We show the conservativeness of the open-loop formulation of the finite time horizon terminal hitting time stochastic
reach-avoid problem for stochastic linear systems, using
conditional expectations and sufficient conditions for Borelmeasurability of the value functions. The open-loop formulation converts the verification problem into a simpler
optimization problem. The objective function is a multidimensional integral, and an analytical expression of the integrand can be obtained using Fourier transforms. For Gaussian
disturbances, the objective function can be evaluated efficiently and the optimization problem is log-concave. Because
the underapproximation technique does not rely on a grid,
it mitigates the curse of dimensionality, and provides nontrivial lower bounds on the stochastic reach-avoid probability.
The method is demonstrated on 40D dynamical system.
R EFERENCES

B. Conservativeness of FTBU
Figure 2(d) and (e) shows the relative error of FTBU with
2
2
respect to DPBDA, with T = [−0.5, 0.5] , S = [−1, 1] ,
grid spacing of 0.05, and  = 0.01. FTBU implemented using
patternsearch has 77.57% grid points with the relative error
less than 30% as compared to 6.6% grid points for fminconbased FTBU. This is also reflected in Figure 2(a), (b), (c).
The sharp rise in Figure 2(d) is due to points where
Ŵ0∗ (x̄0 ) =  and V̂0∗ (x̄0 ) >> , resulting in a large relative
error. The conservativeness of FTBU highlights the role of
feedback in increasing the terminal time probability for any
x̄0 ∈ X . However, as seen in Table I, for sufficiently large

[1] S. Summers and J. Lygeros, “Verification of discrete time stochastic
hybrid systems: A stochastic reach-avoid decision problem,” Automatica, vol. 46, no. 12, pp. 1951–1961, 2010.
[2] A. P. Vinod, B. HomChaudhuri, and M. Oishi, “Forward stochastic
reachability analysis for uncontrolled linear systems using Fourier
transforms,” in Hybrid Systems: Comp. & Control, pp. 35–44.
[3] B. HomChaudhuri, A. P. Vinod, and M. Oishi, “Computation of
forward stochastic reach sets: Application to stochastic, dynamic
obstacle avoidance,” in American Control Conf., Seattle, WA, 2017.
[4] N. Malone, K. Lesser, M. Oishi, and L. Tapia, “Stochastic reachability
based motion planning for multiple moving obstacle avoidance,” in
Proc. Hybrid Syst.: Comput. and Ctrl., 2014, pp. 51–60.
[5] K. Lesser, M. Oishi, and R. Erwin, “Stochastic reachability for control
of spacecraft relative motion,” in IEEE Conf. Dec. Ctrl., 2013, pp.
4705–4712.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 2: Conservativeness of the underapproximation for a double integrator (n = 2). (a) Dynamic programming (DPBDA)
approximation of V̂0∗ (·); Fourier transform-based underapproximation (FTBU) Ŵ0∗ (·) computed (b) using patternsearch and
V̂ ∗ (·)−Ŵ ∗ (·)
(c) using fmincon for each grid point; Relative error defined as 0 V̂ ∗ (·)0 × 100 for V̂0∗ (·) >  and Ŵ0∗ (·) computed
0
(d) using patternsearch and (e) using fmincon for each grid point; Note that higher errors occur closer to the boundary, as
expected, due to the lack of feedback; (f) Improvement in Ŵ0∗ (·) using patternsearch instead of fmincon;  = 0.01

[6] N. Kariotoglou, D. M. Raimondo, S. J. Summers, and J. Lygeros,
“Multi-agent autonomous surveillance: a framework based on stochastic reachability and hierarchical task allocation,” J. Dyn. Sys., Meas.,
Control, vol. 137, no. 3, pp. 031 008–031 008–14, 2014.
[7] A. Abate, M. Prandini, J. Lygeros, and S. Sastry, “Probabilistic
reachability and safety for controlled discrete time stochastic hybrid
systems,” Automatica, vol. 44, no. 11, pp. 2724–2734, 2008.
[8] A. Abate, S. Amin, M. Prandini, J. Lygeros, and S. Sastry, “Computational approaches to reachability analysis of stochastic hybrid
systems,” in Proc. Hybrid Syst.: Comput. and Ctrl., 2007, pp. 4–17.
[9] N. Kariotoglou, S. Summers, T. Summers, M. Kamgarpour, and
J. Lygeros, “Approximate dynamic programming for stochastic reachability,” in Proc. European Ctrl. Conf., 2013, pp. 584–589.
[10] N. Kariotoglou, K. Margellos, and J. Lygeros, “On the computational
complexity and generalization properties of multi-stage and stage-wise
coupled scenario programs,” Sys. & Ctr. Lett., vol. 94, pp. 63–69, 2016.
[11] G. Manganini, M. Pirotta, M. Restelli, L. Piroddi, and M. Prandini,
“Policy search for the optimal control of Markov Decision Processes:
A novel particle-based iterative scheme,” IEEE Trans. Cybern., pp.
1–13, 2015.
[12] J. Ding, M. Kamgarpour, S. Summers, A. Abate, J. Lygeros, and
C. Tomlin, “A stochastic games framework for verification and control
of discrete time stochastic hybrid systems,” Automatica, vol. 49, no. 9,
pp. 2665–2674, 2013.
[13] A. Genz, “Numerical computation of multivariate normal probabilities,” J. of Comp. and Graph. Stat., vol. 1, no. 2, pp. 141–149, 1992.
[14] J. Gleason, A. Vinod, and M. Oishi, “Underapproximation of reachavoid sets for discrete-time stochastic systems via Lagrangian methods,” in IEEE Conf. Dec. Ctrl., https://arxiv.org/abs/1704.03555.
[15] J. A. Gubner, Probability and random processes for electrical and
computer engineers. Cambridge Univ. Press, 2006.
[16] Y. Chow and H. Teicher, Probability Theory: Independence, Interchangeability, Martingales, 3rd ed. Springer New York, 1997.
[17] E. M. Stein and G. L. Weiss, Introduction to Fourier analysis on
Euclidean spaces. Princeton Univ. Press, 1971, vol. 1.
[18] H. Cramér, Mathematical methods of statistics. Princ. Univ. Pr., 1961.
[19] D. Bertsekas and S. Shreve, Stochastic optimal control: The discrete
time case. Academic Press, 1978.

[20] J. Skaf and S. Boyd, “Design of affine controllers via convex optimization,” IEEE Trans. Auto. Ctr., vol. 55, no. 11, pp. 2476–87, 2010.
[21] T. Tao, Analysis II, 2nd ed. Hindustan Book Agency, 2009.
[22] A. S. Nowak, “Universally measurable strategies in zero-sum stochastic games,” The Annals of Probability, pp. 269–287, 1985.
[23] C. J. Himmelberg, T. Parthasarathy, and F. S. VanVleck, “Optimal
plans for dynamic programming problems,” Mathematics of Operations Research, vol. 1, no. 4, pp. 390–394, 1976.
[24] A. P. Vinod and M. Oishi, “Scalable underapproximation for stochastic
reach-avoid problem for high-dimensional LTI systems using Fourier
transforms,” https://arxiv.org/abs/1703.02135.
[25] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical
recipes: The art of scientific computing. Cambridge Univ. Pr., 2007.
[26] A. Genz and F. Bretz, Computation of multivariate normal and t
probabilities. Springer Science & Business Media, 2009, vol. 195.
[27] S. Dharmadhikari and K. Joag-Dev, Unimodality, convexity, and applications. Elsevier, 1988.
[28] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge Univ.
Press, 2004.
[29] A. Genz, “QSCMVNV.” [Online]. Available: http://www.math.wsu.
edu/faculty/genz/software/matlab/qscmvnv.m
[30] T. G. Kolda, R. M. Lewis, and V. Torczon, “Optimization by direct
search: New perspectives on some classical and modern methods,”
SIAM review, vol. 45, no. 3, pp. 385–482, 2003.

