# A Second-Order Reachable Sets Computation Scheme via a Cauchy-Type
Variational Hamilton-Jacobi-Isaacs Equation.

Lekan Molu⋆, Ian Abraham‡, and Sylvia Herbert${}^{\dagger}.$ ⋆Microsoft
Research, 300 Lafayette Street, New York, NY10012, USA. ‡Department of
Mechanical Engineering, Yale University, New Haven, Conn, USA, †Department of
Mechanical and Aerospace Engineering, UC San Diego, La Jolla, CA, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>}.

###### Abstract

Motivated by the scalability limitations of a popular Eulerian variational
Hamilton-Jacobi-Isaacs (HJI) computational method for backward reachability
analysis – essentially providing a least restrictive controller in problems
that involve state or input constraints under a worst-possible disturbance, we
introduce a second-order approximation algorithm for computing the zero
sublevel sets of the variational HJI equation. By continuously approximating
the value functions of all possible trajectories within a dynamical system’s
state space, under sufficient HJI partial differential equation regularity and
continuity conditions throughout the state space, we show that with proper
line search for step size adjustments and state feedback control under the
worst-possible disturbance, we can compute the state set that are reachable
within a prescribed verification time bound.

##  I Introduction

At issue is returning point sets for dynamical systems’ behavioral evolution
(or trajectory) throughout a state space. Such sets must satisfy input or
state constraints in a least restrictive sense under a worst-possible
disturbance [1], and over a time period. This is a classic reachability theory
problem. In this letter, we shall provide an approximation scheme for
constructing the discriminating kernel for a controller’s “constraints-
satisfaction” design specifications problem. Our study is motivated by the
exponential computational requirement of Eulerian methods [2, 3]. Ours is an
extension of differential dynamic programming methods in Bolza-like objective
functions inspired by [4, 5], and [6]’s second-order variational methods for
solving general optimal control problems. In emphasis, we focus on pursuit-
evasion games [7], where in an iterative dynamic game fashion [8], each
agent’s strategy depends on its opponent’s control law as agents locally
approximate successive trajectories and their values that emanate from a state
space. We adopt a formal treatment, indicating how to resolve these problems
computationally.

Reachable sets can be analyzed in a  (i) forwardsense, where system
trajectories are examined to determine if they enter certain states from an
initial set;  (ii) backwardsense, where system trajectories are examined to
determine if they enter certain target sets;  (iii) reach setsense, in which
they are examined to see if states reach a set at a particular time; or  (iv)
reach tubesense, in which they are examined to determine that they reach a set
during a time interval.

We focus on the backward construction. This consists in avoiding an unsafe
state set under the worst-possible disturbance within a given time-bound. The
state sets of a backward reachable problem constitute the discriminating
kernel that is “safety-preserving” for a finite-time horizon control problem.
This discriminating kernel’s boundary satisfies a generalized Isaac’s equation
i.e. if the state $\bm{x}$ is within the set’s boundary, but not the target’s
boundary, then for a co-state $p$ that is an exterior proximal normal at
$\bm{x}$, $\bm{H}(\bm{x},p)\leq 0$, where $\bm{H}(\bm{x},p)$ is the problem’s
Hamiltonian. The set associated with the kernel computed in a backward
reachability framework is termed the backward reachable set or BRS.

Eulerian methods [9, 3] resolve the BRS as the zero-level set of an
implicitly-defined value function on a state space. Constructed as an initial
value problem for a Cauchy-type Hamilton-Jacobi-Isaacs (HJI) partial
differential equation (P.D.E.) [2, 10], the BRS in one dimension is equivalent
to the conservation of momentum equation [3]. For a multidimensional problem,
by resolving the P.D.E. on a dimension-by-dimension basis in a consistent and
monotone fashion [11], a numerically precise and accurate solution to the HJ
P.D.E can be determined [9, 3]. However, as state dimensions increase,
spacetime discretization methods become impractical owing to their exponential
complexity. Contrary to Eulerian methods, our approximation method does not
require state space discretization and hence reduces the exponential
complexity to polynomial time.

Contributions: We describe a computational scheme, provide a summary of the
complexity and convergence behavior for the overapproximated BRS using [2]’s
standard variational Hamilton-Jacobi equation. Iteratively approximating
nonlinear trajectories about “near” local paths and under positive-
definiteness requirement of the stagewise cost function’s second-derivatives
(in control terms), we compute the extrema of the nominal state and control-
disturbance policy pair’s cost to find a cost improvement per iteration.
Trajectories are updated until we sweep all possible initial conditions into
the space of all trajectory costs. Within the bounds here set, the zero
isocontour of all “unionized costs” of all trajectories swept along paths that
assure cost improvement per iteration of the algorithm constitute the zero
level set. This work is the first to systematically provide a polynomial time
complexity computational scheme for backward reachable sets to our knowledge.

The body of this letter is structured as follows: $\S$ II describes the
notations used . In $\S$ III, we introduce methods that will enable us
formulate our proposal in $\S$ IV. We describe the computational and
complexity analysis in $\S$ V.

##  II Notations, Terminology, and Assumptions

We employ standard vector-matrix notations throughout. Conventions: small
Latin and Greek letters are scalars; in bold-font they are vectors.
Exceptions: Players in a differential game e.g. $\bm{P},\,\bm{E}$, are
individual entities. All vectors are column-stacked. Capital Greek and
Calligraphic letters are sets. The scalar product of vectors $\bm{x}\text{ and
}\bm{y}$ is written as $\langle\bm{x},\bm{y}\rangle$. Arbitrary real variables
e.g., $t,t_{0},t_{f},\tau,T$ denote time. For a scalar function $\bm{V}$,
$\bm{V}_{\bm{x}}$ is the gradient vector, and $\bm{V}_{\bm{x}\bm{x}}$ is the
jacobian matrix.

Differential Games: At issue are conflicting objectives between two players in
pursuit-evasion games with each player being a pursuer ($\bm{P}$) or an evader
($\bm{E}$). The set of all strategies executed by $\bm{P}$ (resp. $\bm{E}$)
during a game (beginning at a time $t$) is denoted as $\mathcal{B}(t)$ (resp.
$\mathcal{A}(t)$).

System Description: For the dynamical system

$\displaystyle\dot{\bm{x}}(t)$
$\displaystyle=f(t,\bm{x}(t),\bm{u}(t),\bm{v}(t)),\,\,\bm{x}(0)=\bm{x}_{0},\,\,-T\leq
t\leq 0,$ (1)

where the state $\bm{x}(t)$ evolves from some initial negative time $-T$ to a
final time $0$. Assumptions: The flow field $f(t,\cdot,\cdot,\cdot)$ and
$\bm{x}(t)$ are bounded and Lipschitz continuous for fixed controls
$\bm{u}(t)$ and $\bm{v}(t)$; $f(t,\cdot,\cdot,\cdot)$ possesses
$\mathcal{C}^{2}$ continuity in $\bm{x}$111This bounded Lipschitz continuity
property assures uniqueness of the system response $\bm{x}(t)$ to controls
$\bm{u}(t)$ and $\bm{v}(t)$ [12].. We take $\bm{x}(t)$ to belong in the open
set $\Omega\subset\mathbb{R}^{n}$ and the pair $(\bm{x},t)$ as the system’s
phase. The Cartesian product of $\Omega$ and the space $T=\mathbb{R}^{1}$ of
all time values is termed the phase space, $\Omega\times T$. The closure of
$\Omega$ is denoted $\bar{\Omega}$ and we let $\delta\Omega$ denote the
boundary of $\Omega$. Unless otherwise stated, vectors
$\bm{u}(t)\in\mathbb{R}^{n_{u}}$ and $\bm{v}(t)\in\mathbb{R}^{n_{v}}$ are
reserved for admissible control (resp. disturbance) at time $t$ for some
$n_{u},\,n_{v}>0$. Controls $\bm{u}(t)$ (resp. $\bm{v}(t)$) are piecewise
continuous in $t$, if for each $t$, $\bm{u}\in\mathcal{U}$ (resp.
$\bm{v}\in\mathcal{V}$) and $\mathcal{U}(\text{resp. }\mathcal{V})$ are
Lebesgue measurable and compact sets. At all times, any of $\bm{u}(t)$ or
$\bm{v}(t)$ will be under the influence of a player such that the motion of
$\bm{x}(t)$ will be influenced by the will of that player. When a control law
or value function is optimal, it shall be signified by an asterisk superscript
e.g. $\bm{u}^{\star}$. For the phase space $\left(\Omega\times T\right)$, the
set of all controls for players $\bm{P}$ and $\bm{E}$ are respectively drawn
from

$\displaystyle\mathcal{\bar{U}}$
$\displaystyle\equiv\\{\bm{u}:[-T,0]\rightarrow\mathcal{U}|\bm{u}\text{
measurable},\,\mathcal{U}\subset\mathbb{R}^{n_{u}}\\},$ (2)
$\displaystyle\mathcal{\bar{V}}$
$\displaystyle\equiv\\{\bm{v}:[-T,0]\rightarrow\mathcal{V}|\bm{v}\text{
measurable},\,\mathcal{V}\subset\mathbb{R}^{n_{v}}\\},$ (3)

with $\mathcal{U},\,\mathcal{V}$ being compact.

Existence and Uniqueness of Value and Trajectories: For any admissible
control-disturbance pair $(\bm{u}(t),\bm{v}(t))$ and initial phase
$(\bm{x},-T)$, given a pursuer’s non-anticipative strategy, the game admits a
value [13] and there exists a unique trajectory $\bm{\xi}(t)$ such that the
motion of (1) passing through phase $(\bm{x},-T)$ under the action of control
$\bm{u}(t)$, and disturbance $\bm{v}(t)$, and observed at a time $t$
afterwards, given by,

$\displaystyle\bm{\xi}(t)=\bm{\xi}(t;-T,\bm{x},\bm{u}(\cdot),\bm{v}(\cdot))$
(4)

satisfies (1) almost everywhere (a.e.) [2]. Note that under feedback
strategies, the game may not admit a value necessarily [14]; however, if the
target sets are initialized such that there is no overlap of states between
players, and that the evader(s) are not captured by the pursuer(s) at the
start of the game, then one can guarantee a pursuit winning strategy [15,
Assumption 2.1].

##  III Background.

Reachable sets in the context of dynamic programming and two person games is
here introduced. We restrict attention to computing the backward reachable
sets of a dynamical system. We establish the viscosity solution P.D.E. to the
terminal HJI P.D.E, and then describe the formulation of the BRS and backward
reachable tube (BRT). Let us enquire.

###  III-A The Backward Reachable (Target) Set

For any optimal control problem, a value function is constructed based on a
user-defined optimal cost that is bounded and uniformly continuous for any
input phase $(\bm{x},-T)$ e.g. reach goal at the end of a time horizon i.e.,

$\displaystyle|g(0;\bm{x})|$ $\displaystyle\leq
k,\,\,|g(0;\bm{x})-g(t;\hat{\bm{x}})\mid\leq k|\bm{x}-\hat{\bm{x}}\mid$ (5)

for constant $k$ and all $-T\leq t\leq 0$,
$\hat{\bm{x}},\,\bm{x}\in\mathbb{R}^{n}$. The set

$\displaystyle\mathcal{L}_{0}=\\{\bm{x}\in\bar{\Omega}\,|\,g(0;\bm{x})\leq
0\\},$ (6)

is the target set in the phase space $\Omega\times\mathbb{R}$ (proof in [2]).
This target set can represent the failure set (to avoid) or a goal set (to
reach) in the state space.

###  III-B The Backward Reachable Tube

Backward reachability analysis seeks to capture all conditions under which
trajectories of the system may enter a user-defined target set cf. (6). This
could be desirable in goal-regions of the state (safe sets) or undesirable
state configurations (unsafe sets). For a target set construction problem, a
differential game’s (lower) value is equivalent to a solution of (1) for
$\bm{u}(t)$ and $\bm{v}(t)=\beta[\bm{u}](t)$ i.e.,

$\displaystyle\bm{V}(\bm{x},t)=\inf_{\beta\in\mathcal{B}(t)}\sup_{\bm{u}\in\mathcal{U}(t)}\min_{t\in[-T,0]}g\left(0;t,\bm{x},\bm{u}(\cdot),\bm{v}(\cdot)\right).$
(7)

Optimal trajectories emanating from an initial phase $(\bm{x},-T)$ where the
value function is non-negative will maintain non-negative cost over an entire
time horizon, thereby avoiding the target set and vice versa. For the safety
problem setup in (7), we can define the corresponding robustly controlled
backward reachable tube as the closure of the open set

$\displaystyle\mathcal{L}([\tau,0],\mathcal{L}_{0})$
$\displaystyle=\\{\bm{x}\in\Omega\,|\,\exists\,\beta\in\mathcal{\bar{V}}(t)\,\forall\,\bm{u}\in\mathcal{U}(t),\exists\,\tau\in[-T,0],$
$\displaystyle\qquad\qquad\bm{\xi}(\bar{t})\in\mathcal{L}_{0}\\}.$ (8)

Read: The set of states from which the strategies of $\bm{P}$ and for all
controls of $\bm{E}$ imply that we reach and remain in the target set in the
interval $[-T,0]$. Following Lemma 2 of [2], the states in the reachable set
admit the following properties w.r.t the value function $\bm{V}$:

$\displaystyle\bm{x}(t)\in\mathcal{L}(\cdot)\implies\bm{V}(\bm{x},t)\leq
0,\bm{V}(\bm{x},t)\leq 0\implies\bm{x}(t)\in\mathcal{L}(\cdot).$ (9)

Player $\bm{P}$ is minimizing (the game’s termination time c.f. (6)), seeking
to drive system trajectories into the unsafe set; and $\bm{E}$ is maximizing
(the game’s termination time) i.e. is seeking to avoid the unsafe set-222For
the goal-satisfaction (or liveness) problem setups, the strategies are
reversed and the backward reachable tube are the states from which the evader
$\bm{E}$ can successfully reach the target set despite worst-case efforts of
the pursuer $\bm{P}$..

###  III-C The Terminal HJI Value Function for Reachability

The non-anticipative strategy used in level set methods follows from [12]
viz., suppose that the pursuer’s strategy starting at a time $-T$, is
$\beta:\mathcal{\bar{U}}\rightarrow\mathcal{\bar{V}}$. Suppose further that
$\beta$ is provided for each $-T\leq\tau\leq 0$ and
$\bm{u},\hat{\bm{u}}\in\mathcal{\bar{U}}$. Then,

$\displaystyle\begin{cases}&\bm{u}(\bar{t})=\hat{\bm{u}}(\bar{t})\,\,\text{
a.e. on }-T\leq\bar{t}\leq\tau,\\\ &\text{implies
}\beta[\bm{u}](\bar{t})=\beta[\hat{\bm{u}}](\bar{t})\,\,\text{ a.e. on
}-T\leq\bar{t}\leq\tau.\end{cases}$ (10)

In our work, however, we use feedback strategies in an iterative dynamic game
framework while initializing the target sets to ensure the existence of a
value as discussed in §II. It is well-known that a differential game’s value
function admits a “viscosity” (generalized) solution [16, 17] of the
associated HJ-Isaacs (HJI) PDE given by

$\displaystyle\frac{\partial\bm{V}}{\partial t}{(\bm{x},t)}$
$\displaystyle+\min\\{0,\bm{H}(t;\bm{x},p)\\}=0,\,\,V(\bm{x},0)=g(\bm{x})$
(11)

where the vector field $p\,(\equiv\bm{V}_{\bm{x}})$ is known in terms of the
game’s terminal conditions so that the overall game is akin to a two-point
boundary-value problem; and the Hamiltonian
$\bm{H}(t;\bm{x},\bm{u},\bm{v},\bm{V}_{\bm{x}})$ is defined as

$\displaystyle\bm{H}(t;\bm{x},\bm{u},\bm{v},\bm{V}_{\bm{x}})=\max_{u\in\mathcal{U}}\min_{v\in\mathcal{V}}\,\langle
f(t;\bm{x},\bm{u},\bm{v}),\bm{V}_{\bm{x}}\rangle.$ (12)

Equation (7) only allows the game to determine that a trajectory belongs in
the target set at exactly time zero. That is, the evader can chase a
trajectory into the target set and escape it before reaching the final time,
$0$. The $\min$ operation between the scalar $0$ and the Hamiltonian allows
the pursuer to “freeze” trajectories that may be under the evader’s willpower
when the evader tries to evolve such trajectories outside of the target set.
For more details on the construction of this P.D.E., see [2].

##  IV Successive Approximation Scheme

Throughout this section, all feasible trajectories
$\\{\bm{x}_{i},\bm{x}_{i+1},\cdots,\bm{x}_{n}\\}\in\Omega$ are scheduled for
iterative second-order expansion along the nonlinear state variations
$\delta\bm{x}_{i}$ and about nominal states $\bar{\bm{x}}_{i}$; we derive the
variational linear differential equation form of (11) and describe the
integration scheme of the backward and forward steps of standard DDP on the
variational HJI equation. In what follows, for ease of notation we drop the
subscripts that denote the respective trajectories. For clarity’s sake, we
shall drop the subscripts and derive the optimal value for a generic
trajectory $\bm{x}$.

###  IV-A Local Approximations to Nonlinear Trajectories

Suppose that system (1) is controllable everywhere on an interval $[-T,0]$
along a local trajectory $\bm{x}_{r}$ generated by admissible control
disturbance pair $\bm{u}_{r}$, $\bm{v}_{r}$ starting from an initial phase
$(\bm{x}_{r},-T)$ so that

$\displaystyle\dot{\bm{x}}_{r}(\tau)$
$\displaystyle=f(t;\bm{x}_{r}(\tau),\bm{u}_{r}(\tau),\bm{v}_{r}(\tau)).$ (13)

First, we apply local controls $\bm{u}_{r}(t)\,\text{ and }\bm{v}_{r}(t)$ on
(1) so that the nominal value is $\bm{V}(t;\cdot,\cdot,\cdot)$ for a resulting
nominal state $\bm{x}_{r}(\tau)$; $\tau\in\left[-T,0\right]$. System dynamics
that describe variations from the nonlinear system cf. (1) with state and
control pairs $\delta\bm{x}(t),\,\delta\bm{u}(t),\,\delta\bm{v}(t)$
respectively333Note that $\delta\bm{x}(t),\,\delta\bm{u}(t),\,\text{ and
}\delta\bm{v}(t)$ are respectively measured with respect to
$\bm{x}(t),\bm{u}(t),\bm{v}(t)$ and are not necessarily small. can then be
written as

$\displaystyle\bm{x}(t)$
$\displaystyle=\bm{x}_{r}(t)+\delta\bm{x}(t),\,\,\bm{u}(t)=\bm{u}_{r}(t)+\delta\bm{u}(t),$
(14a) $\displaystyle\bm{v}(t)$
$\displaystyle=\bm{v}_{r}(t)+\delta\bm{v}(t),\,\,t\in\left[-T,0\right].$ (14b)

Abusing notation, we drop the templated time arguments in the variations (14)
so that the canonical problem is now

$\displaystyle\dfrac{d}{dt}\left(\bm{x}_{r}+\delta\bm{x}\right)$
$\displaystyle=f(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v}),\,$
(15) $\displaystyle\bm{x}_{r}(-T)$
$\displaystyle+\delta\bm{x}(-T)=\bm{x}(-T),$ (16)

whose terminal value admits the optimal form (see (11)):

$\displaystyle\begin{split}-\frac{\partial\bm{V}^{\ast}}{\partial
t}(\bm{x}_{r}+\delta\bm{x},t)&=\min\left\\{\bm{0},\max_{\delta\bm{u}\in\mathcal{U}}\,\min_{\delta\bm{v}\in\mathcal{V}}\left\langle
f(t;\bm{x}_{r}+\delta\bm{x},\right.\right.\\\
&\left.\left.\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v}),\dfrac{\partial\bm{V}^{\ast}}{\partial\bm{x}}\left(\bm{x}_{r}+\delta\bm{x},t\right)\right\rangle\right\\},\end{split}$
$\displaystyle\bm{V}^{\ast}(\bm{x}_{r}+\delta\bm{x},0)$
$\displaystyle=g(0;\bm{x}_{r}(0)+\delta\bm{x}(0));$ (17)

and state trajectory

$\displaystyle\bm{\xi}(t)=\bm{\xi}(t;-T,\bm{x}_{r}+\delta\bm{x},\bm{u}+\delta\bm{u},\bm{v}+\delta\bm{v}),\quad
t\in\left[-T,0\right].$ (18)

For $-T<t\leq 0$ and a $\tau\in[t,0]$, let the optimal cost for using the
optimal control
$\bm{u}^{\star}(\tau)=\bm{u}_{r}(\tau)+\delta\bm{u}^{\star}(\tau)$ and optimal
disturbance $\bm{v}^{\star}(\tau)=\bm{v}_{r}(\tau)+\delta\bm{v}^{\star}(\tau)$
be $\bm{V}^{\star}(\bm{x}_{r},\tau)$; and the nominal cost for using
$\bm{u}_{r}(\tau)$ be $\bm{V}_{r}(\bm{x}_{r},t)$, so that on the phase
$\left(\bm{x}_{r},t\right)$, we have

$\displaystyle\mathbf{\tilde{V}}^{\star}(\bm{x}_{r},t)=\bm{V}^{\star}(\bm{x}_{r},t)-\bm{V}_{r}(\bm{x}_{r},t).$
(19)

Next, we expand $\bm{V}^{\ast}$ in (17) under sufficient regularity
assumptions.

###  IV-B Power Series Expansion Scheme

Suppose that the optimal terminal cost, $\bm{V}^{\ast}(\cdot)$, is
sufficiently smooth to allow a power series expansion in the state variation
$\delta\bm{x}$ about the nominal state, $\bm{x}_{r}$, then we must have

$\displaystyle\bm{V}^{\star}(\bm{x}_{r}+\delta\bm{x},t)$
$\displaystyle=\mathbf{\tilde{V}}^{\star}(\bm{x}_{r},t)+\bm{V}_{r}(\bm{x}_{r},t)+\left\langle\bm{V}_{\bm{x}}^{\ast}(\bm{x}_{r},t),\delta\bm{x}\right\rangle$
$\displaystyle\quad+\dfrac{1}{2}\left\langle\delta\bm{x},\bm{V}^{\star}_{\bm{x}\bm{x}}(\bm{x}_{r},t)\delta\bm{x}\right\rangle+\text{h.o.t
in }\delta\bm{x},$ (20)

where h.o.t. signifies higher order terms. The smoothness assumption is
necessary for admitting the linear differential equation expansion of (20).
This expansion is consistent with differential dynamic programming schemes [6,
4, 18].

Observe: If $\delta\bm{x}$ is not constrained to be small, (20) may require
huge memory for storage owing to the large dimensionality of terms beyond
second order. However,  (i) if $\bm{x}_{r}$ is constrained to be sufficiently
close to $\bm{x}$, the state variation $\delta\bm{x}$ will be small, resulting
in $\bm{x}\approx\bm{x}_{r}$ cf. (14)444A scheme to keep the variation
$\delta\bm{x}$ small is discussed in section V. Ultimately, the
$\delta\bm{u},\delta\bm{v}$ terms will be quadratic in $\delta\bm{x}$ if we
neglect h.o.t.;  (ii) for small $\delta\bm{x}$, the expansion of (20) becomes
consistent with second-order methods where quadratic terms in $\delta\bm{x}$
dominate higher-order terms.  Hence, we can avoid the infinite data storage
requirement by truncating the expansion in (20) at second-order terms in
$\delta\bm{x}$. This will incur an $O(\delta\bm{x}^{3})$ approximation error,
affording us realizable control laws that can be executed on the system (1).
Thus, we rewrite (20) as

$\displaystyle\bm{V}^{\star}(\bm{x}_{r}+\delta\bm{x},t)$
$\displaystyle=\mathbf{\tilde{V}}^{\star}(\bm{x}_{r},t)+\bm{V}_{r}(\bm{x}_{r},t)+\left\langle\bm{V}_{\bm{x}}^{\ast}(\bm{x}_{r},t),\delta\bm{x}\right\rangle$
$\displaystyle\quad+\dfrac{1}{2}\left\langle\delta\bm{x},\bm{V}^{\star}_{\bm{x}\bm{x}}(\bm{x}_{r},t)\delta\bm{x}\right\rangle.$
(21)

Denote by $\bm{V}_{\bm{x}}^{\star}\left(\bm{x}_{r}+\delta\bm{x},t\right)$ the
optimal value of the co-state
$\frac{\partial\bm{V}^{\star}}{\partial\bm{x}}\left(\bm{x}_{r}+\delta\bm{x},t\right)$
on the r.h.s of (17). Then, expanding up to second order we find that

$\displaystyle\bm{V}_{\bm{x}}^{\star}\left(\bm{x}_{r}+\delta\bm{x},t\right)=\bm{V}^{\star}_{\bm{x}}\left(\bm{x}_{r},t\right)+\langle\bm{V}^{\star}_{\bm{x}\bm{x}}\left(\bm{x}_{r},t\right),\delta\bm{x}\rangle,$
(22)

where we have again omitted the quadratic term
$\bm{V}^{\star}_{\bm{x}\bm{x}\bm{x}}\delta\bm{x}\delta\bm{x}$ owing to the
foregoing reason. Substituting (21) and (22) into (17), abusing notation by
dropping the templated phase arguments, we find that

$\displaystyle\begin{split}-\frac{\partial\mathbf{\tilde{V}}^{\ast}}{\partial
t}-\frac{\partial\bm{V}_{r}}{\partial
t}-\left\langle\dfrac{\partial\bm{V}_{\bm{x}}^{\ast}}{\partial
t},\delta\bm{x}\right\rangle-\dfrac{1}{2}\left\langle\delta\bm{x},\dfrac{\partial\bm{V}_{\bm{x}\bm{x}}^{\ast}}{\partial
t}\delta\bm{x}\right\rangle&=\\\
\min\left\\{0,\max_{\delta\bm{u}\in\mathcal{U}}\,\min_{\delta\bm{v}\in\mathcal{V}}\left\langle
f^{T}(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v}),\right.\right.\\\
\left.\left.\bm{V}^{\ast}_{\bm{x}}+\bm{V}^{\ast}_{\bm{x}\bm{x}}\,\delta\bm{x}\right\rangle\right\\}.\end{split}$
(23)

From (19), we may write

$\displaystyle\frac{d}{dt}\left(\bm{V}_{r}+\mathbf{\tilde{V}}^{\ast}\right)=\dfrac{\partial}{\partial
t}\left(\bm{V}_{r}+\mathbf{\tilde{V}}^{\ast}\right)+\left\langle
f^{T}(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r}),\bm{V}_{\bm{x}}^{\ast}\right\rangle$
(24a) $\displaystyle\dot{\bm{V}}_{\bm{x}}^{\ast}$
$\displaystyle=\dfrac{\partial\bm{V}_{\bm{x}\bm{x}}^{\ast}}{\partial
t}+\langle
f^{T}(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r}),\bm{V}_{\bm{x}\bm{x}}^{\ast}\rangle,\,\,\dot{\bm{V}}_{\bm{x}\bm{x}}^{\ast}=\dfrac{\partial\bm{V}_{\bm{x}\bm{x}}^{\ast}}{\partial
t}.$ (24b)

Therefore, (23) in light of (24) becomes

$\displaystyle\begin{split}-\frac{\partial\mathbf{\tilde{V}}^{\ast}}{\partial
t}-\frac{\partial\bm{V}_{r}}{\partial
t}-\left\langle\dfrac{\partial\bm{V}^{\ast}_{\bm{x}}}{\partial
t},\delta\bm{x}\right\rangle-\dfrac{1}{2}\left\langle\delta\bm{x},\dfrac{\partial\bm{V}^{\ast}_{\bm{x}\bm{x}}}{\partial
t}\delta\bm{x}\right\rangle&=\\\
\min\left\\{\bm{0},\max_{\delta\bm{u}\in\mathcal{U}}\,\min_{\delta\bm{v}\in\mathcal{V}}\left[\bm{H}(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v},\bm{V}_{\bm{x}}^{\ast})+\right.\right.\\\
\left.\left.\left\langle\bm{V}_{\bm{x}\bm{x}}^{\ast}\,\delta\bm{x},f(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v})\right\rangle\right]\right\\}\end{split}$
(25)

where
$\bm{H}(t;\bm{x},\bm{u},\bm{v},\bm{V}_{\bm{x}}^{\ast})=\langle\bm{V}_{\bm{x}}^{\ast},f(t;\bm{x},\bm{u},\bm{v})\rangle$.

###  IV-C Variational HJI Linear Differential Equation

Let the respective maximizing and minimizing control-disturbance pair when
$\bm{x}=\bm{x}_{r}$ at time $t$ be

$\displaystyle\bm{u}^{\ast}$
$\displaystyle=\bm{u}_{r}+\delta\bm{u}^{\ast},\,\,\bm{v}^{\ast}=\bm{v}_{r}+\delta\bm{v}^{\ast}.$
(26)

Given the variation $\delta\bm{x}$ on the r.h.s. of (25), write the respective
maximizing evader’s and minimizing pursuer’s controls for state
$\bm{x}=\bm{x}_{r}+\delta\bm{x}$

$\displaystyle\bm{u}$
$\displaystyle=\bm{u}_{r}+\delta\bm{u}^{\ast}+\delta\bm{u},\,\,\bm{v}=\bm{v}_{r}+\delta\bm{v}^{\ast}+\delta\bm{v}.$
(27)

Note here that $\delta\bm{u}$ and $\delta\bm{v}$ are as yet to be determined.
Henceforward, we adopt the abbreviations
$\bm{H}=\bm{H}(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v},\bm{V}_{\bm{x}}^{\ast})$
and
$f=f(t;\bm{x}_{r}+\delta\bm{x},\bm{u}_{r}+\delta\bm{u},\bm{v}_{r}+\delta\bm{v})$.
Expanding the r.h.s. of (25) about $\bm{x}_{r},\bm{u}^{\ast},\bm{v}^{\ast}$ at
time $t$ up to second-order as before, we have

$\displaystyle\begin{split}\min\left\\{\bm{0},\max_{\delta\bm{u}\in\mathcal{U}}\,\min_{\delta\bm{v}\in\mathcal{V}}\left[\bm{H}+\left\langle\bm{H}_{\bm{x}}+\bm{V}^{\ast}_{\bm{x}\bm{x}}f,\delta\bm{x}\right\rangle+\langle\bm{H}_{\bm{u}},\delta\bm{u}\rangle+\right.\right.\\\
\left.\left.\langle\bm{H}_{\bm{v}},\delta\bm{v}\rangle+\langle\delta\bm{u},(\bm{H}_{\bm{u}\bm{x}}+f_{\bm{u}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast})\delta\bm{x}\rangle+\left\langle\delta\bm{v},(\bm{H}_{\bm{v}\bm{x}}+f_{\bm{v}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast})\delta\bm{x}\right\rangle\right.\right.\\\
\left.\left.+\dfrac{1}{2}\left\langle\delta\bm{u},\bm{H}_{\bm{u}\bm{v}}\delta\bm{v}+\dfrac{1}{2}\left\langle\delta\bm{v},\bm{H}_{\bm{v}\bm{u}}\delta\bm{u}\right\rangle\right\rangle+\dfrac{1}{2}\left\langle\delta\bm{u},\bm{H}_{\bm{u}\bm{u}}\delta\bm{u}\right\rangle\right.\right.\\\
\left.\left.+\dfrac{1}{2}\left\langle\delta\bm{v},\bm{H}_{\bm{v}\bm{v}}\delta\bm{v}\right\rangle+\dfrac{1}{2}\left\langle\delta\bm{x},\left(\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}+\bm{V}_{\bm{x}\bm{x}}^{\ast}f_{\bm{x}}\right)\delta\bm{x}\right\rangle\right]\right\\}.\end{split}`$
(28)

When capture occurs i.e. when $\bm{E}$’s separation from $\bm{P}$ becomes less
than a pre-specified (capture) radius, we must have

$\displaystyle\bm{H}_{\bm{u}}(t;\bm{x}_{r},\bm{u}^{\star},\bm{v},\bm{V}_{\bm{x}})=0;\,\bm{H}_{\bm{v}}(t;\bm{x}_{r},\bm{u},\bm{v}^{\star},\bm{V}_{\bm{x}})=0$
(29)

for $t\in\left(-T,0\right]$. Seeking variational feedback controllers of the
form:
$\delta\bm{u}=\bm{k}_{\bm{u}}\delta\bm{x},\,\delta\bm{v}=\bm{k}_{\bm{v}}\delta\bm{x}$
and putting (29) into (28) we may write

$\displaystyle\bm{H}_{\bm{u}\bm{u}}\delta\bm{u}+\left(\bm{H}_{\bm{u}\bm{x}}+f_{\bm{u}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}\right)\delta\bm{x}+\dfrac{1}{2}\left(\bm{H}_{\bm{u}\bm{v}}+\bm{H}_{\bm{v}\bm{u}}^{T}\right)\delta\bm{v}=0,$
(30a)
$\displaystyle\bm{H}_{\bm{v}\bm{v}}\delta\bm{v}+\left(\bm{H}_{\bm{v}\bm{x}}+f_{\bm{v}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}\right)\delta\bm{x}+\dfrac{1}{2}\left(\bm{H}_{\bm{v}\bm{u}}+\bm{H}_{\bm{u}\bm{v}}^{T}\right)\delta\bm{u}=0$
(30b)

so that evaluating (30b) at $\bm{x}_{r}$, $\bm{u}^{\ast}$, and $\bm{v}^{\ast}$
we find that

$\displaystyle\bm{k}_{\bm{u}}$
$\displaystyle=-\bm{H}_{\bm{u}\bm{u}}^{-1}\left[\bm{H}_{\bm{u}\bm{v}}\bm{k}_{\bm{v}}+\left(\bm{H}_{\bm{u}\bm{x}}+f_{\bm{u}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}\right)\right],\text{
and that }$ $\displaystyle\bm{k}_{\bm{v}}$
$\displaystyle=-\bm{H}_{\bm{v}\bm{v}}^{-1}\left[\bm{H}_{\bm{v}\bm{u}}\bm{k}_{\bm{u}}+\left(\bm{H}_{\bm{v}\bm{x}}+f_{\bm{v}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}\right)\right].$
(31)

In equation (31), the two players’ control strategies are interdependent: the
optimal control law for the pursuer is to choose a control strategy that is an
optimal response to the evader’s control input. This is akin to a Newton-
Raphson scheme which “correctly” estimates the maximizing
$\delta\bm{u}^{\ast}$ and and minimizing $\delta\bm{v}^{\ast}$ provided that
$|\bm{u}-\bm{u}^{\ast}|,|\bm{v}-\bm{v}^{\ast}|<\epsilon\text{ for }\epsilon>0$
with the added requirement that $\bm{H}_{\bm{u}\bm{u}}$ and
$\bm{H}_{\bm{v}\bm{v}}$ are positive-definite. Upon substitution, (28)
satisfies

$\displaystyle\begin{split}\min\left\\{\bm{0},\bm{H}+\left\langle\bm{H}_{\bm{x}}+\bm{V}_{\bm{x}\bm{x}}^{\ast}f,\delta\bm{x}\right\rangle+\frac{1}{2}\left\langle\delta\bm{x},\right.\right.\\\
\left.\left.\left(\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}+\bm{V}_{\bm{x}\bm{x}}^{\ast}f_{\bm{x}}+\bm{k}_{\bm{u}}^{T}\bm{H}_{\bm{u}\bm{u}}\bm{k}_{\bm{u}}+\bm{k}_{\bm{v}}^{T}\bm{H}_{\bm{v}\bm{v}}\bm{k}_{\bm{v}}\right)\delta\bm{x}\right\rangle\right\\}.\end{split}$
(32)

We discard terms involving the pairs $(\delta\bm{u},\,\delta\bm{x})$ and
$(\delta\bm{v},\,\delta\bm{x})$ beyond second order. Comparing (32) to the
l.h.s. of (25), we find that

$\displaystyle-\dfrac{\partial\tilde{\bm{V}}^{\ast}}{\partial
t}-\dfrac{\partial\bm{V}_{r}}{\partial t}=\min\\{\bm{0},\bm{H}\\},$ (33a)
$\displaystyle-\dfrac{\partial\bm{V}_{\bm{x}}^{\ast}}{\partial t}$
$\displaystyle=\min\left\\{\bm{0},\bm{H}_{\bm{x}}+\bm{V}_{\bm{x}\bm{x}}^{\ast}\,f\right\\},$
(33b) $\displaystyle-\dot{\bm{V}_{\bm{x}\bm{x}}^{\ast}}$
$\displaystyle=\min\left\\{\bm{0},\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}+\bm{V}_{\bm{x}\bm{x}}^{\ast}f_{\bm{x}}+\bm{k}_{\bm{u}}^{T}\bm{H}_{\bm{u}\bm{u}}\bm{k}_{\bm{u}}+\bm{k}_{\bm{v}}^{T}\bm{H}_{\bm{v}\bm{v}}\bm{k}_{\bm{v}}\right\\}.$

Furthermore, comparing (33) with (24), and using the first-order necessary
condition for optimality cf. (29), we have implies

$\displaystyle-\dot{\tilde{\bm{V}}}^{\ast}$
$\displaystyle=\min\\{\bm{0},\bm{H}-\bm{H}(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r},\bm{V}_{\bm{x}}^{\ast})\\}$
(34a) $\displaystyle-\dot{\bm{V}}_{\bm{x}}^{\ast}$
$\displaystyle=\min\left\\{0,\bm{H}_{\bm{x}}+\bm{V}_{\bm{x}\bm{x}}^{\ast}\left(f-f(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r})\right)\right\\}$
(34b) $\displaystyle-\dot{\bm{V}_{\bm{x}\bm{x}}^{\ast}}$
$\displaystyle=\min\left\\{\bm{0},\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\bm{V}_{\bm{x}\bm{x}}^{\ast}+\bm{V}_{\bm{x}\bm{x}}^{\ast}f_{\bm{x}}+\bm{k}_{\bm{u}}^{T}\bm{H}_{\bm{u}\bm{u}}\bm{k}_{\bm{u}}+\bm{k}_{\bm{v}}^{T}\bm{H}_{\bm{v}\bm{v}}\bm{k}_{\bm{v}}\right\\}$

where every quantity in (34) is evaluated at
$t;\bm{x}_{r},\bm{u}^{\ast},\bm{v}^{\ast}$555The “∗” sign implies that the
optimal $\bm{u}^{\ast}=\bm{u}_{r}+\delta\bm{u}^{\ast}$ and
$\bm{v}^{\ast}=\bm{v}_{r}+\delta\bm{v}^{\ast}$ are used – essentially the
optimal control-disturbance pair for $\bm{x}=\bm{x}_{r}+\delta\bm{x}$.. Note
that equation (34) signifies the terms to be computed in the typical backward
pass of DDP-like algorithms namely, start at the final time $t=0$ and proceed
backwards until $-T$; whilst storing the locally linear control gains i.e.
(31) at every step of the backward pass to keep the necessary conditions of
optimality. In a forward pass, the following new controls

$\displaystyle\bm{u}(\tau)$
$\displaystyle=\bm{u}^{\ast}(\tau)+\bm{k}_{\bm{u}}\delta\bm{x}(\tau),$ (35a)
$\displaystyle\bm{v}(\tau)$
$\displaystyle=\bm{v}^{\ast}(\tau)+\bm{k}_{\bm{v}}\delta\bm{x}(\tau),\,\tau\in\left[-t,0\right],\,t>T$
(35b)

are updated proceeding forward in time from $-T$ to $0$ 666Note that gains
$\bm{k}_{\bm{u}},\,\bm{k}_{\bm{v}}$ are only used in computing feedback
controllers for the discriminating kernel in a backward reachability setting..

The linear differential variational HJI equation for (11) thus becomes (for
sufficiently small $\delta\bm{x}$):

$\displaystyle\begin{split}-\frac{\partial\hat{\mathbf{\tilde{V}}}}{\partial
t}-\frac{\partial\bm{V}_{r}}{\partial
t}-\left\langle\dfrac{\partial\hat{\bm{V}}_{\bm{x}}}{\partial
t},\delta\bm{x}\right\rangle-\dfrac{1}{2}\left\langle\delta\bm{x},\dfrac{\partial\hat{\bm{V}}{\bm{x}\bm{x}}}{\partial
t}\delta\bm{x}\right\rangle=\min\left\\{\bm{0},\right.\\\
\left.\left[\bm{H}(t;\bm{x}_{r},\bm{u}^{\ast},\bm{v}^{\ast},\hat{\bm{V}}_{\bm{x}})+\left\langle\bm{H}_{\bm{x}}+\hat{\bm{V}}_{\bm{x}\bm{x}}f,\delta\bm{x}\right\rangle+\langle\bm{H}_{\bm{u}},\bm{k}_{\bm{u}}\delta\bm{x}\rangle\right.\right.\\\
\left.\left.+\langle\bm{H}_{\bm{v}},\bm{k}_{\bm{v}}\delta\bm{x}\rangle\right.\right.\\\
\left.\left.+\langle\bm{k}_{\bm{u}}\delta\bm{x},(\bm{H}_{\bm{u}\bm{x}}+f_{\bm{u}}^{T}\hat{\bm{V}}_{\bm{x}\bm{x}})\delta\bm{x}\rangle+\left\langle\bm{k}_{\bm{v}}\delta\bm{x},(\bm{H}_{\bm{v}\bm{x}}+f_{\bm{v}}^{T}\hat{\bm{V}}_{\bm{x}\bm{x}})\delta\bm{x}\right\rangle\right.\right.\\\
\left.\left.+\dfrac{1}{2}\left\langle\bm{k}_{\bm{u}}\delta\bm{x},\bm{H}_{\bm{u}\bm{v}}\delta\bm{v}+\dfrac{1}{2}\left\langle\delta\bm{v},\bm{H}_{\bm{v}\bm{u}}\bm{k}_{\bm{u}}\delta\bm{x}\right\rangle\right\rangle+\dfrac{1}{2}\left\langle\bm{k}_{\bm{u}}\delta\bm{x},\bm{H}_{\bm{u}\bm{u}}\bm{k}_{\bm{u}}\delta\bm{x}\right\rangle\right.\right.\\\
\left.\left.+\dfrac{1}{2}\left\langle\bm{k}_{\bm{v}}\delta\bm{x},\bm{H}_{\bm{v}\bm{v}}\bm{k}_{\bm{v}}\delta\bm{x}\right\rangle+\dfrac{1}{2}\left\langle\delta\bm{x},\left(\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\hat{\bm{V}}_{\bm{x}\bm{x}}+\hat{\bm{V}}_{\bm{x}\bm{x}}f_{\bm{x}}\right)\delta\bm{x}\right\rangle\right]\right\\}.\end{split}$
(36)

where the hat terms indicate that (35) is used. Note that we have discarded
the $\bm{V}_{\bm{x}\bm{x}\bm{x}}$ terms in our derivations. We also do away
with the application of the penalizing $\epsilon>0$-term on the variational
controls in (35) as proposed in [6] in lieu of standard step-size adjustment
mechanisms to ensure a sufficiently small $\delta\bm{x}$.

In a similar spirit to (34), we have

$\displaystyle-\dot{\hat{\bm{V}}}^{\ast}$
$\displaystyle=\min\\{\bm{0},\bm{H}-\bm{H}(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r},\hat{\bm{V}}_{\bm{x}}^{\ast})\\}$
(37a) $\displaystyle-\dot{\hat{\bm{V}}}^{\ast}_{\bm{x}}$
$\displaystyle=\min\left\\{0,\bm{H}_{\bm{x}}+\hat{\bm{V}}^{\ast}_{\bm{x}\bm{x}}\left(f-f(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r})\right)\right\\}$
(37b) $\displaystyle-\dot{\hat{\bm{V}}}^{\ast}_{\bm{x}\bm{x}}$
$\displaystyle=\min\left\\{\bm{0},\bm{H}_{\bm{x}\bm{x}}+f_{\bm{x}}^{T}\hat{\bm{V}}^{\ast}_{\bm{x}\bm{x}}+\hat{\bm{V}}^{\ast}_{\bm{x}\bm{x}}f_{\bm{x}}+\bm{k}_{\bm{u}}^{T}\bm{H}_{\bm{u}\bm{u}}\bm{k}_{\bm{u}}+\bm{k}_{\bm{v}}^{T}\bm{H}_{\bm{v}\bm{v}}\bm{k}_{\bm{v}}\right\\}$

with the caret symbols signifying predictions upon application of the policies
(35).

##  V Computational Procedure and Discussion.

Suppose that the position of points, $\\{g.\bm{x}_{i}\\}_{i=1}^{N}$, where all
trajectories emerge is known. In what follows, we describe the numerical
scheme for carrying out the integrations (37) and computing the
overapproximated level set of $\bm{V}(\bm{x}_{r}+\delta\bm{x},t)$ i.e. (36).

###  V-A Computational Scheme

First, a schedule of controls $\\{\bm{u}_{r}(t),\,\bm{v}_{r}(t)\\}_{t=-T}^{0}$
needed for computing nominal states $\\{\bm{x}_{r}(t)\\}_{t=-T}^{0}$ is
initialized as a problem-dependent parameter and then used to run
$\\{\bm{x}_{r}(t)\\}_{t=-T}^{0}$; if the nominal control schedules are not
available, they can be set from the system’s passive dynamics. The predicted
cost improvement starting from the final time and going backwards in time is
(cf. (37)):

$\displaystyle\mid\hat{\bm{V}}^{\ast}(\bm{x},\tau)\mid=\int_{0}^{\tau}\min\\{\bm{0},\bm{H}-\bm{H}(t;\bm{x}_{r},\bm{u}_{r},\bm{v}_{r},\hat{\bm{V}}_{\bm{x}}^{\ast})\\},\,\tau\gg-T,$
(38)

while the actual cost improvement is given by (19). Similar to [4], define a
cost improvement criterion $\rho>0$ such that,

$\displaystyle{\mathbf{\tilde{V}}^{\star}(\bm{x}_{r},t)}/{\mid\hat{\bm{V}}^{\ast}(\bm{x}_{r},t)\mid}>\rho$
(39)

determines the “closeness” of the cost improvement to the cost prediction.

Algorithm 1 Successive Approximation Scheme.

1:procedure VarHJIInt($\eta$ $\rho$) $\triangleright$ Given $\eta>0$,
$\rho\in(0,1]$. $\triangleright$ Stop conditions.

2: Initialize Buffer $V_{buf}=\emptyset$.

3: for all initial trajectories
$X=\\{\bm{\bm{x}}_{i}(t_{0})\\}_{i=1}^{N}\in\Omega$ do

4: Generate schedule
$\pi_{r_{i}}=\\{\bm{u}_{r_{i}}(t),\,\bm{v}_{r_{i}}(t)\\}_{t=0}^{t=-T}$;

5: Initialize $\mid\hat{\bm{V}}_{i}(\bm{x}_{r_{i}},t)\mid=\infty$;
$\triangleright$ cf. (39);

6: $\bm{x}_{r_{i}},\pi_{i}^{\ast}$ =
BackwardPass($\bm{x}_{i},\pi_{r_{i}},\mid\hat{\bm{V}}_{i}(\cdot)\mid$);

7:
$\bm{V}_{i}^{\ast}(\bm{x}_{r_{i}}+\delta\bm{x}_{i},t),\,\pi_{r_{i}}^{\ast}$=
ForwardPass($\bm{x}_{r_{i}},\pi_{i}^{\ast})$;

8:
$\bm{V}_{buf}\leftarrow\bm{V}_{buf}\cup\bm{V}_{i}^{\ast}(\bm{x}_{r_{i}}+\delta\bm{x}_{i},t)$;

9: end for

10: Compute zero-levelset of $V_{buf}$ $\triangleright$ Using [19].

11:end procedure

1:function
BackwardPass($\bm{x}_{i},\pi_{r_{i}},\mid\hat{\bm{V}}_{i}^{\ast}\mid$)

2: Initialize $t_{eff}=-T,\,k=1\cdots K$.

3: for $t=0,-k\Delta t,\cdots-T$ do $\triangleright$ $\Delta t=T/(K-1)$;

4: Unpack $\bm{u}_{r_{i}}(t),\,\bm{v}_{r_{i}}(t):=\pi_{r_{i}}(t)$;

5: $\dot{\bm{x}}_{r_{i}}(t)\leftarrow
f(\bm{x}_{r_{i}},\bm{u}_{r_{i}}(t),\,\bm{v}_{r_{i}}(t))$ & compute
$\bm{V}_{r_{i}}(\bm{x}_{r_{i}},t)$;

6: Compute
$\mid\hat{\bm{V}}_{i}^{\ast}(\bm{x},t)\mid\leftarrow\bm{V}_{i}^{\star}(\bm{x}_{r_{i}},t)-\bm{V}_{r_{i}}(\bm{x}_{r_{i}},t)$;

7: if $\mid\hat{\bm{V}}_{i}^{\ast}(\bm{x}_{r_{i}},t)\mid<\eta$ then

8: Update $\bm{x}_{r_{i}}$$\triangleright$ e.g. via Runge-Kutta integration;

9: Accept $\bm{x}_{r_{i}}$ only if condition (39) is satisfied;

10: Set $t_{eff}=t$ and Terminate loop.

11: else Line search for a smaller $\delta\bm{x}_{r_{i}}$; restart line 3.

12: end if

13: $\delta\pi^{\ast}_{i}(t):=(\delta\bm{u}_{i}^{\ast}$,
$\delta\bm{v}_{i}^{\ast})\leftarrow$ extrema of
$\bm{H}_{i}(\bm{x}_{r_{i}},\pi^{\ast}_{r_{i}},p,t)$

14: Update
$\pi_{i}^{\ast}(t):=(\bm{u}_{i}^{\ast}(t),\,\bm{v}_{i}^{\ast}(t))\leftarrow\pi_{r_{i}}(t)+\delta\pi^{\ast}_{i}(t)$;

15: end for

16: return
$\bm{x}_{r_{i}},\pi_{i}^{\ast}:=\\{\bm{u}_{i}^{\ast}(t),\,\bm{v}_{i}^{\ast}(t)\\}_{t=0}^{t=t_{eff}}$.

17:end function

1:function ForwardPass($\bm{x}_{r_{i}},\pi_{i}^{\ast}$)

2: for $t=-T,(K-1)\Delta t,\cdots,-2\Delta t,-\Delta t,0$ do

3: Unpack controllers
$(\bm{u}_{i}^{\ast}(t),\,\bm{v}_{i}^{\ast}(t)):=\pi_{i}^{\ast}(t)$;

4: Compute the extremizing $\delta\bm{u}\,\&\,\delta\bm{v}$ cf. (28), (31);

5: Run $\delta\dot{\bm{x}}_{i}(t)\leftarrow
f(\cdot,\delta\bm{u}_{i}(t),\,\delta\bm{v}_{i}(t))$;

6: $\delta\bm{x}_{i}(t)\leftarrow\delta\dot{\bm{x}}_{r_{i}}(t)$ set
$t_{f}=t_{ext}$ $\triangleright$ e.g. Euler/RK integration;

7: Compute $\bm{V}_{i}^{\ast}(\bm{x}_{r_{i}}+\delta\bm{x}_{i},t)$ cf. (20);

8: Update $\pi_{r_{i}}^{\ast}(t):=\left(\bm{u}_{r}(t),\,\bm{v}_{r}(t)\right)$
$\triangleright$ Using (26);

9: end for

10: return
$\bm{V}_{i}^{\ast}(\bm{x}_{r_{i}}+\delta\bm{x}_{i},0),\,\pi_{r_{i}}^{\ast}$.

11:end function

As seen in Algorithm 1, the procedure proceeds in two passes for all initial
trajectories:  (i) in a backward pass, costs (37) are estimated with the open
loop sequence $\\{\bm{u}_{r}(t),\,\bm{v}_{r}(t)\\}_{0}^{t=-T}$. We generate
$\bm{u}^{\ast}$ and $\bm{v}^{\ast}$ in (26) afterwards;  (ii) in a forward
pass, reversing the order of integration limits, controls
$\bm{u}^{\ast}(t),\,\bm{v}^{\ast}(t),\,\,t\in\left[-T,0\right]$ of (27) are
then applied and the approximation in (20) is computed for every trajectory
that emanates from the state space.  For the line search procedure, an Armijo-
Goldstein condition in a typical backtracking line search can be applied to
iteratively keep $\delta\bm{x}$ small for a valid approximation of
$\bm{V}(\cdot)$ i.e. (33). A regularization scheme similar to our previous
work [8] can also be applied to keep the stagewise Hessians positive definite.

The union operator from line 8 of Algorithm 1 allows the recovery of the
maximal BRS as proposed in [20]. Given the “safe” and “unsafe” sets that
constitute $\bm{V}_{buf}$, the reachable sets can be over-approximated under a
best-response strategy of the two-player game (See [8]). We have released code
that computes the zero isocontour (levelset) of the (union of all
trajectories’) optimal value function stipulated on Line 10 of Algorithm 1. A
method for obtaining this is available in [19], which is an implementation of
[21].

###  V-B Computational Complexity and Convergence Requirement

The problem introduced in (33) can be solved with Newton’s method to find the
extremizing policies $\delta\pi$ for all $t$ in Algorithm 1: with Hessian
$\frac{\partial^{2}V}{\partial\pi_{u}^{2}}$ for the maximizer and
$\frac{\partial^{2}V}{\partial u^{2}},\frac{\partial^{2}V}{\partial v^{2}}$
for the minimizer, the inversion of the Hessian matrix would constitute
$O(T^{3}m^{3})$ CPU flops. Whereas the overall DDP-style computational scheme
we have presented has a CPU cost per iteration of
$O(N)+T\cdot(2n^{3}+\frac{7}{2}n^{2}(n_{u}+n_{v})+2n(n_{u}+n_{v})^{2}+\frac{1}{3}(n_{u}+n_{v})^{3}+O(n^{2})+O((n_{u}+n_{v})^{2}))$
CPU flops (see [22, Appendix II] for details), where $n_{u}$ and $n_{v}$ are
as given in section II. The polynomial time complexity of the presented scheme
makes it more attractive compared to Newton’s method or even level set methods
which are well-known to scale exponentially. We refer readers to Pantoja [23]
for a thorough differentiation between DDP and Newton’s methods. In our
opinion, recent first-order primal-dual algorithms such as Chambolle-Pock [24]
may prove more computationally parsimonious for these problems.

A Note on Convergence: Conditions upon which the algorithm 1 converges is
premised on the standard Hessians’ positive-definiteness (PD) i.e.
$\bm{H}_{\bm{u}\bm{u}}$ and $\bm{H}_{\bm{v}\bm{v}}$ of DDP algorithms. In
addition, for linear dynamical equations (1), the cost function being PD
convex is a sufficient requirement for PD stagewise Hessians. When 1 is
nonlinear, stagewise PD of $\bm{H}_{\bm{u}\bm{u}}$ and $\bm{H}_{\bm{v}\bm{v}}$
is no longer guaranteed [25]. In such situations, one may explore  (i) a
Levenberg-Marquardt scheme, convert the scheme to steepest descent by turning
the stagewise Hessian to an identity or using the active shift method of
[25]’s Theorem IV.

## References

  * [1] J. Lygeros, “On reachability and minimum cost optimal control,” _Automatica_ , vol. 40, no. 6, pp. 917–927, 2004.
  * [2] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A Time-Dependent Hamilton-Jacobi Formulation of Reachable Sets for Continuous Dynamic Games,” _IEEE Transactions on Automatic Control_ , vol. 50, no. 7, pp. 947–957, 2005.
  * [3] S. Osher and R. Fedkiw, “Level Set Methods and Dynamic Implicit Surfaces,” _Applied Mechanics Reviews_ , vol. 57, no. 3, pp. B15–B15, 2004.
  * [4] D. H. Jacobson and D. Q. Mayne, _Differential Dynamic Programming_. American Elsevier Publishing Company, Inc., New York, NY, 1970.
  * [5] A. E. Bryson and W. F. Denham, “A steepest-ascent method for solving optimum programming problems,” 1962.
  * [6] D. H. Jacobson, “Differential Dynamic Programming Methods for Determining Optimal Control of Nonlinear Systems,” _PhD Thesis, University of London_ , vol. 103, no. 10, pp. 626–627, 1967.
  * [7] R. Isaacs, “Differential games: A mathematical theory with applications to warfare and pursuit, control and optimization.” _Kreiger, Huntigton, NY_ , 1965.
  * [8] O. Ogunmolu, N. Gans, and T. Summers, “Minimax iterative dynamic game: Application to nonlinear robot control tasks,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 6919–6925.
  * [9] J. A. Sethian, “Level Set Methods And Fast Marching Methods: Evolving Interfaces In Computational Geometry, Fluid Mechanics, Computer Vision, And Materials Science,” _Robotica_ , vol. 18, no. 1, pp. 89–92, 2000.
  * [10] I. M. Mitchell, “A Toolbox of Level Set Methods.” _UBC Department of Computer Science Technical Report TR-2007-11_ , 2007.
  * [11] M. G. Crandall and P. L. Lions, “Two Approximations of Solutions of Hamilton-Jacobi Equations,” _Mathematics of Computation_ , vol. 43, no. 167, p. 1, 1984.
  * [12] L. Evans and P. E. Souganidis, “Differential Games And Representation Formulas For Solutions Of Hamilton-Jacobi-Isaacs Equations,” _Indiana Univ. Math. J_ , vol. 33, no. 5, pp. 773–797, 1984.
  * [13] R. J. Elliott and N. J. Kalton, “Cauchy Problems for Certain Isaacs-Bellman Equations and Games of Survival,” _Transactions of the American Mathematical Society_ , vol. 198, no. 1970, p. 45, 1974.
  * [14] P. Cardaliaguet, “Nonsmooth Semipermeable Barriers, Isaacs’ Equation, And Application to a Differential Game with One Target and Two Players,” _Applied Mathematics and Optimization_ , vol. 36, no. 2, pp. 125–146, 1997\.
  * [15] R. Yan, X. Duan, Z. Shi, Y. Zhong, and F. Bullo, “Matching-based Capture Strategies for 3D Heterogeneous Multiplayer Reach-Avoid Differential Games,” _Automatica_ , vol. 140, p. 110207, 2019.
  * [16] P.-L. Lions, _Generalized solutions of Hamilton-Jacobi equations_. London Pitman, 1982, vol. 69.
  * [17] M. G. Crandall and P.-L. Lions, “Viscosity solutions of hamilton-jacobi equations,” _Transactions of the American mathematical society_ , vol. 277, no. 1, pp. 1–42, 1983.
  * [18] W. Denham, “Book Reviews: Differential Dynamic Programming,” _Review of Social Economy_ , vol. 36, no. 2, pp. 228–229, 1978.
  * [19] O. Ogunmolu, “A Numerical Toolbox for the Scalable Analysis of Hamilton-Jacobi PDEs.” https://github.com/robotsorcerer/LevelSetPy, 2022, accessed March 11, 2022.
  * [20] M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin, “Decomposition of Reachable Sets and Tubes for a Class of Nonlinear Systems,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 11, pp. 3675–3688, 2018.
  * [21] W. E. Lorensen and H. E. Cline, “Marching cubes: A high resolution 3d surface construction algorithm,” _ACM siggraph computer graphics_ , vol. 21, no. 4, pp. 163–169, 1987.
  * [22] L.-Z. Liao and C. A. Shoemaker, “Convergence in unconstrained discrete-time differential dynamic programming,” _IEEE Transactions on Automatic Control_ , vol. 36, no. 6, pp. 692–706, 1991.
  * [23] J. D. O. Pantoja, “Differential Dynamic Programming and Newton’s Method,” _International Journal of Control_ , vol. 47, no. 5, pp. 1539–1553, 1988\.
  * [24] A. Chambolle and T. Pock, “A First-Order Primal-Dual Algorithm for Convex Problems With Applications to Imaging,” _Journal of Mathematical Imaging and Vision_ , vol. 40, no. 1, pp. 120–145, 2011.
  * [25] L.-Z. Liao and C. A. Shoemaker, “Convergence in unconstrained discrete-time differential dynamic programming,” _IEEE Transactions on Automatic Control_ , vol. 36, no. 6, pp. 692–706, 1991.