Time-optimal path planning in dynamic flows using level set equations: theory and schemes

Abstract

We develop an accurate partial differential equation-based methodology that predicts the time-optimal paths of autonomous vehicles navigating in any continuous, strong, and dynamic ocean currents, obviating the need for heuristics. The goal is to predict a sequence of steering directions so that vehicles can best utilize or avoid currents to minimize their travel time. Inspired by the level set method, we derive and demonstrate that a modified level set equation governs the time-optimal path in any continuous flow. We show that our algorithm is computationally efficient and apply it to a number of experiments. First, we validate our approach through a simple benchmark application in a Rankine vortex flow for which an analytical solution is available. Next, we apply our methodology to more complex, simulated flow fields such as unsteady double-gyre flows driven by wind stress and flows behind a circular island. These examples show that time-optimal paths for multiple vehicles can be planned even in the presence of complex flows in domains with obstacles. Finally, we present and support through illustrations several remarks that describe specific features of our methodology.

Appendices

Appendix A: Preliminaries

In this Appendix, we describe some of the relevant definitions and terminology needed for the theoretical results. Most of the material presented in this Appendix may be found in the study of Bardi and Capuzzo-Dolcetta (2008), Clarke et al. (1998), Cannarsa and Sinestrari (2004), Frankowska (1989), and Bressan (2011). In what follows, we let $n \in \mathbb {N}$, ${\Omega } \subseteq \mathbb {R}^{n}$ be an open set and $\xi : {\Omega } \rightarrow \mathbb {R}$.

Remark 1

Let $\xi \in \mathcal {C}({\Omega })$. Let ∂ ₊ ξ(x ₀) and ∂ ₋ ξ(x ₀) denote the sets of super- and sub-differentials (Bardi and Capuzzo-Dolcetta 2008; Clarke et al. 1998) of ξ at x ₀. Then, q∈∂ ₊ ξ(x ₀) (resp. ∂ ₋ ξ(x ₀)) if and only if there exists a function $\gamma \in \mathcal {C}^{1}({\Omega })$ such that γ(x ₀)=ξ(x ₀), ∇γ(x ₀)=q and the function γ−ξ has a strict local minima (resp. maxima) at x ₀.

Definition 1 (Generalized gradient)

Let ξ be locally Lipschitz at x ₀. For any $\mathbf {u} \in \mathbb {R}^{n}$, let $\xi ^{g}(\mathbf {x}_{0};\mathbf {u} )$ denote the generalized directional derivative of ξ at x ₀ (Clarke et al. 1998). The set of generalized gradients of ξ at x ₀ is the non-empty set

$$\begin{array}{@{}rcl@{}} \partial \xi(\mathbf{x}_{0}) = \left \{ \mathbf{q} \in \mathbb{R}^{n}: \forall \, \mathbf{u} \in \mathbb{R}^{n}, \mathbf{q} \cdot \mathbf{u} \le \xi^{g}(\mathbf{x}_{0};\mathbf{u}) \right \}\, . \end{array} $$

(20)

Definition 2 (Regular function)

ξ is said to be regular at x ₀∈Ω if it is Lipschitz near x ₀ and admits directional derivatives $\xi ^{d}(\mathbf {x}_{0};\mathbf {u})$ for all $\mathbf {u} \in \mathbb {R}^{n}$, with ξ ^g(x ₀;u)=ξ ^d(x ₀;u).

Properties of regular functions

• 1. If ξ is continuously differentiable at x ₀, then it is regular at x ₀. Furthermore, ξ ^d(x ₀;u)=∇ξ(x ₀)⋅u=ξ ^g(x ₀;u) for all $\mathbf {u} \in \mathbb {R}^{n}$.

• 2. If ξ is convex and Lipschitz near x ₀, then it is regular at x_0.

• 3. Let ξ be regular at x ₀∈Ω. Then,

  $$ \partial_{-}\xi(\mathbf{x}_{0}) = \partial \xi(\mathbf{x}_{0}) \, . $$

  (21)

Definition 3 (Viscosity solution)

Let F≥0 and let V(x,t) satisfy assumptions 3–4. Consider the Hamilton-Jacobi equation

$$ \frac{\partial \phi}{\partial t} + F |\nabla \phi| + \mathbf{V}(\mathbf{x,t}) \cdot \nabla \phi = 0\quad\text{in}\quad{\Omega} \times (0,\infty) \, , $$

(22)

with initial conditions

$$ \phi(\mathbf{x},0) = \nu(\mathbf{x}) \, , $$

(23)

where $\nu :{\Omega }\rightarrow \mathbb {R}$ is Lipschitz continuous. A function $\phi \in \mathcal {C}({\Omega } \times [0,\infty ))$ is a viscosity subsolution of Eq. 22 if for every $(\mathbf {x,t})\in {\Omega } \times (0,\infty )$ and (q , p)∈∂ ₊ ϕ(x , t),

$$ p + F |\mathbf{q}| + \mathbf{V}(\mathbf{x},t) \cdot \mathbf{q} \le 0 \, . $$

(24)

A function $\phi \in \mathcal {C}({\Omega } \times [0,\infty ))$ is a viscosity supersolution of Eq. 22 if for every $(\mathbf {x,t})\in {\Omega }\times (0,\infty )$ and (q , p)∈∂ ₋ ϕ(x , t),

$$ p + F |\mathbf{q}| + \mathbf{V}(\mathbf{x},t) \cdot \mathbf{q} \ge 0 \, . $$

(25)

ϕ is said to be a viscosity solution of Eq. 22 if it is both a viscosity subsolution and a viscosity supersolution.

Theorem 1

(Frankowska 1989 ) A locally Lipschitz function $\phi :{\Omega } \times (0,\infty )\rightarrow \mathbb {R}$ is a viscosity solution to Eq. 22 if and only if for every $(\mathbf {x,t}) \in {\Omega }\times (0,\infty )$ ,

$$ \max_{(\mathbf{q,p}) \in \partial \phi(\mathbf{x,t})} \{ p + F|\mathbf{q}| + \mathbf{V}(\mathbf{x},t) \cdot \mathbf{q} \}= 0\, $$

(26)

and for all (q,p)∈∂₋ ϕ(x,t),

$$ p + F|\mathbf{q}| + \mathbf{V}(\mathbf{x},t) \cdot \mathbf{q = 0} \, . $$

(27)

Theorem 2

(Clarke et al. 1998 ) [Lebourg’s Mean Value Theorem] Let $\mathbb {S} \subseteq \mathbb {R}$ be an open set. Let $x, y \in \mathbb {S}$ and suppose that $f:\mathbb {S}\rightarrow \mathbb {R}$ is Lipschitz on an open set containing the segment [x,y]. Then, there exists 0<λ<1 such that

$$ f(y) - f(x) = g \times (y-x) \, , $$

(28)

for some g∈∂f(z), where z=λx+(1−λ)y.

Theorem 3

(Clarke et al. 1998 )[Chain Rule] Let ${\Omega }_{1} \subseteq \mathbb {R}^{n}$ and ${\Omega }_{2} \subseteq \mathbb {R}^{m}$ be two open sets with $m, n \in \mathbb {N}$. Let $\mathbf {g}:{\Omega }_{1}\rightarrow {\Omega }_{2}$ be continuously differentiable near x∈Ω₁, and let $F:{\Omega }_{2}\rightarrow \mathbb {R}$ be Lipschitz near g(x). Then, f :=F ∘g is Lipschitz near x and

$$ \partial f(\mathbf{x}) \subseteq \left (\mathbf{g}^{\prime}(\mathbf{x}) \right )^{*} \partial F \left (\mathbf{g}(\mathbf{x}) \right ) \, , $$

(29)

where ^∗ denotes the adjoint.

Appendix B: Theoretical results

We now state a lemma that provides a monotonicity result related to ϕ, the viscosity solution of the Hamilton-Jacobi Equation (22). According to this result, the generalized gradient of ϕ is non-positive on trajectories $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)$, along the direction $\left (\frac {\text {d} \widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)}{\text {d} t},1\right )$ for t>0. This lemma is then used to prove Theorem 4, which establishes the relationship between reachable sets and the viscosity solution of a modified level set equation.

Lemma 1

Let ${\Omega } \subseteq \mathbb {R}^{n}$ be open, F>0 and let V(x,t) satisfy assumptions (3 )–( 4 ). Let ϕ be the viscosity solution to Eq. 22 . Let the trajectory $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)$ satisfy ( 1 ) with initial conditions $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},0) = \mathbf {y}_{\mathbf {s}}$. Then,

• 1.
  $$ p + \frac{\text{d} \widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t} \cdot \mathbf{q} \le 0 \, \quad \forall \, (\mathbf{q,p}) \in \partial \phi(\widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t),t) $$

  (30)
• 2.
  $$ {\phi}^{g} \left (\widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t),t;\left (\frac{\text{d} \widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t},1 \right ) \right ) \le 0 \quad \forall \, t>0 \, . $$

  (31)

The proof of this Lemma may be found in Lolla et al. (2014).

Theorem 4

Let ${\Omega } \subseteq \mathbb {R}^{n}$ be an open set, $\mathbf {V}(\mathbf {x},t):{\Omega } \times [0,\infty ) \rightarrow \mathbb {R}^{n}$ satisfy ( 3 )–( 4 ), and F≥0. Let ${T}^{o}(\mathbf {y}):{\Omega } \rightarrow \mathbb {R}$ denote the optimal first arrival time at y. Let the trajectory $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)$ satisfy (1) with initial conditions $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},0) = \mathbf {y}_{\mathbf {s}}$ . Let ϕ ^o(x , t) be the viscosity solution to the Hamilton-Jacobi equation (9) with initial condition (10). Then,

• 1. $\phi ^{o}(\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t),t) \le 0$ for all t≥0.

• 2. If ϕ ^o is regular at $(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ for all t>0 and $\mathbf {X}^{o}_{P}$ satisfies

  $$ \frac{\text{d} \mathbf{X}^{o}_{P}}{\text{d} t} = F \frac{\mathbf{q}^{o}}{|\mathbf{q}^{o}|} + \mathbf{V}(\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t), \quad t > 0, $$

  (32)

  for some $(\mathbf {q}^{o},p^{o}) \in \partial \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ , then

  $$ \phi^{o}(\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t) = 0 \quad \forall \, t\ge 0 \, . $$
• 3.
  $$ {T}^{o}(\mathbf{y}) = \inf_{t\ge 0} \{t:\phi^{o}(\mathbf{y,t}) = 0\} \,, $$

  (33)

  where $\inf$ denotes the infimum.

• 4. The optimal trajectory to y _f ∈Ω satisfies (11) whenever ϕ ^o is differentiable at $(\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ and $|\nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)| \ne 0$.

• 5. If $F > \max \limits _{\mathbf {x} \in {\Omega }, t\ge 0} |\mathbf {V}(\mathbf {x},t)|$, then T ^o(y) is the viscosity solution of the modified Eikonal equation

  $$ F |\nabla {T}^{o}(\mathbf{y})| + \mathbf{V}(\mathbf{y},{T}^{o}(\mathbf{y})) \cdot \nabla {T}^{o}(\mathbf{y}) - 1= 0 \, , \mathbf{y} \in {\Omega} \, . $$

  (34)

Proof

• 1. The viscosity solution to Eq. 9 is locally Lipschitz (see Tonon (2011), Bianchini and Tonon (2012), and Cannarsa and Sinestrari (2004)). We now argue that ${\phi ^{o}_{P}}(t) := \phi ^{o}(\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ is locally Lipschitz for all t≥0. Observe that ${\phi ^{o}_{P}}(t) = \phi ^{o}(\mathbf {g}_{P}(t))$ where $\mathbf {g}_{P}(t) := (\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$. Since g _P (t) is continuously differentiable in $(0,\infty )$ with $\frac {\text {d} \mathbf {g}_{P}(t)}{\text {d} t} = \left (\frac {\text {d} \widetilde {\mathbf {X}}_{P}}{\text {d} t},1 \right )$ and ϕ ^o is locally Lipschitz, ${\phi ^{o}_{P}}(t)$ is also locally Lipschitz in $(0,\infty )$ by the chain rule stated in Theorem 3.

  Let t ₁ > 0 be fixed. Since ${\phi ^{o}_{P}}$ is locally Lipschitz, there exists an open interval around t ₁ in which ${\phi ^{o}_{P}}$ is Lipschitz. Thus, for any t ₂>t ₁ in this interval, Lebourg’s Mean Value Theorem (Theorem 2) implies there exist t ₃∈(t ₁,t ₂) and $s \in \partial {\phi ^{o}_{P}} (t_{3})$ such that

  $$ {\phi^{o}_{P}}(t_{2}) - {\phi^{o}_{P}}(t_{1}) = s \times (t_{2} - t_{1}) \, . $$

  (35)

  Using the chain rule of Theorem 3 again (^∗ denotes the adjoint),

  $$\begin{array}{@{}rcl@{}} \partial {\phi^{o}_{P}}(t_{3}) &\subseteq \, \left (\mathbf{g}_{P}^{\prime}(t_{3}) \right )^{*} \partial \phi^{o} \left (\mathbf{g}_{P}(t_{3}) \right ) \\ & = \left \{ p + \mathbf{q}\cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t_{3})}{\text{d} t} \right . \\ &\qquad :\left . (\mathbf{q,p}) \in \partial \phi^{o}(\widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t_{3}),t_{3}) \right \}\, . \end{array} $$

  (36)

  Hence, any $s \in \partial {\phi ^{o}_{P}}(t_{3})$ can be written as

  $$ s = p + \mathbf{q} \cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t_{3})}{\text{d} t} \, . $$

  (37)

  for some $(\mathbf {q,p}) \in \partial \phi ^{o}(\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t_{3}),t_{3})$. From Eq. 30,

  $$\begin{array}{@{}rcl@{}} p + \mathbf{q} \cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t_{3})}{\text{d} t} &\le 0 \, \end{array} $$

  implying that for any $s \in \partial {\phi ^{o}_{P}}(t_{3})$, s≤0. Using this result in Eq 35 yields ${\phi ^{o}_{P}}(t_{2}) \le {\phi ^{o}_{P}}(t_{1})$ for all t ₁,t ₂. Since ${\phi ^{o}_{P}}$ is locally Lipschitz in $(0,\infty )$, we conclude that ${\phi ^{o}_{P}}(t)$ is non-increasing on $(0,\infty )$. Moreover, since ${\phi ^{o}_{P}}$ is continuous on $[0,\infty )$, with ${\phi ^{o}_{P}}(0) = 0$ (from Eq. 10) and non-increasing in $(0,\infty )$, we have ${\phi ^{o}_{P}}(t) = \phi (\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t),t) \le 0$ for all t≥0.

• When the trajectory $\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t)$ is regular, i.e., when ϕ ^o is regular at points $(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ for all t > 0, Eq. 21 implies $\partial _{-} \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t) = \partial \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ for all t > 0. Since ϕ ^o is the viscosity solution to Eq. 9, we obtain from Theorem 1 that for any t > 0, $p + F|\mathbf {q}| + \mathbf {V}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t) \cdot \mathbf {q} = 0$ for all $(\mathbf {q,p}) \in \partial _{-} \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t) = \partial \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$. Specifically, for the member (q ^o,p ^o) of $\partial \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ that satisfies (32), the definition of generalized gradient (20) implies

  $$\begin{array}{@{}rcl@{}} {\phi^{o}}^{g} &\left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t;\left (\frac{\text{d} \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t},1\right ) \right ) \\ &\ge p^{o} + \mathbf{q}^{o} \cdot \frac{\text{d} \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t} \\ &= p^{o} + F |\mathbf{q}^{o}| + \mathbf{V}(\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t) \cdot \mathbf{q}^{o} \\ &= 0 \,. \end{array} $$

  (38)

  Combining this result with Eq. 31, we obtain

  $${\phi^{o}}^{g}\left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t;\left (\frac{\text{d} \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t},1\right ) \right ) = 0 \, . $$

  Since ϕ ^o is regular at $\left (\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t \right )$ by assumption, we then also have

  $$ {\phi^{o}}^{d}\left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t;\left (\frac{\text{d} \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t},1\right ) \right ) = 0 \, . $$

  (39)

  For any h > 0, the definition of ${\phi ^{o}_{P}}$ implies

  $$\begin{array}{@{}rcl@{}} & \left | \frac{{\phi^{o}_{P}}(t+h) - {\phi^{o}_{P}}(t)}{h} \right | \\ &=\left | \frac{\phi^{o} \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t+h),t+h \right ) - \phi^{o} \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t\right )}{h} \right | \, . \end{array} $$

  (40)

  Since ϕ ^o is locally Lipschitz, ∃ C > 0 such that for h > 0 small enough,

  $$\begin{array}{@{}rcl@{}} &\left | \phi \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t+h),t+h \right ) - \phi^{o} \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t) + h\frac{\text{d} \mathbf{X}^{o}_{P}}{\text{d} t},t+h \right ) \right | \\ &\le C \left | \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t+h) - \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t) - h\frac{\text{d} \mathbf{X}^{o}_{P}}{\text{d} t}\right | \\ &= C \left | \mathbf{o}(h) \right | \, , \end{array} $$

  (41)

  where $\mathbf {o}(h) \in \mathbb {R}^{n}$ denotes a vector whose individual terms are o(h). Adding and subtracting $\phi ^{o} \left (\mathbf {X}^{o}_{P} + h\frac {\text {d} \mathbf {X}^{o}_{P}}{\text {d} t},t+h \right )$ from the numerator of Eq. 40 and using the triangle inequality, we obtain

  $$\begin{array}{cl} &\left | \frac{{\phi^{o}_{P}}(t+h) - {\phi^{o}_{P}}(t)}{h} \right | \\ &\le \left | \frac{\phi^{o} \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t) + h\frac{\text{d} \mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t},t+h \right ) - \phi^{o} \left (\mathbf{X}^{o}_{P}(\mathbf{y}_{\mathbf{s}},t),t\right )}{h} \right | \\ &\qquad+ C\left | \frac{\mathbf{o}(h)}{h}\right |. \end{array} $$

  The first term on the right converges to:${\phi ^{o}}^{d}\left (\mathbf {X}^{o}_{P},t;\left (\frac {\text {d} \mathbf {X}^{o}_{P}}{\text {d} t},1\right ) \right )$ as h ↓0 and by Eq. 39, its value is zero. The second term uniformly converges to zero as h ↓0, by definition. This implies

  $$\lim\limits_{h \downarrow 0} \left | \frac{{\phi^{o}_{P}}(t+h) - {\phi^{o}_{P}}(t)}{h} \right | = 0 \, , $$

  and consequently that

  $$ \lim\limits_{h \downarrow 0} \frac{{\phi^{o}_{P}}(t+h) - {\phi^{o}_{P}}(t)}{h} = 0 \, . $$

  (42)

  Since Eq. 42 holds for all t > 0, ${\phi ^{o}_{P}}$ is right differentiable in $(0,\infty )$ and the value of the right-derivative is zero for all t > 0. This implies that ${\phi ^{o}_{P}}$ is constant in $(0,\infty )$. Since ${\phi ^{o}_{P}}(0) = 0$, we obtain ${\phi ^{o}_{P}}(t) = \phi ^{o}(\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t),t)=0$ for all t≥0. Therefore, trajectories $\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t)$ that are regular and satisfy (32) always remain on the zero-level set of ϕ ^o.

• It has been shown in part (1) that $\phi ^{o}(\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t),t) \le 0$ for all t≥0 for any trajectory $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)$ that satisfies (1) and the initial conditions $\widetilde {\mathbf {X}}_{P}(\mathbf {y}_{\mathbf {s}},0) = \mathbf {y}_{\mathbf {s}}$. Therefore, for a trajectory X _P (y _s ,t) that reaches a given end point y∈Ω at time $\widetilde {T}(\mathbf {y})$ (not necessarily optimal),

  $$ \phi^{o}(\mathbf{y}, \widetilde{T}(\mathbf{y})) = \phi^{o}(\widetilde{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},\widetilde{T}(\mathbf{y})),\widetilde{T}(\mathbf{y}))\le 0 \, . $$

  (43)

  Since this inequality holds for any arbitrary arrival time $\widetilde {T}(\mathbf {y})$, it will also hold for the optimal arrival time T ^o(y), implying

  $$ \phi^{o}(\mathbf{y}, {T}^{o}(\mathbf{y})) \le 0 \,~ \text{for all}\quad \mathbf{y} \in {\Omega} \, . $$

  (44)

  For y=y _s , Eq. 33 holds trivially. For any y≠y _s , ϕ ^o(y,0) > 0 by Eq. 10. The continuity of ϕ ^o and Eq. 44 together then yield

  $$ {T}^{o}(\mathbf{y}) \ge \inf_{t\ge 0}\{t: \phi^{o}(\mathbf{y,t}) = 0\} \, . $$

  (45)

  In part (2), we showed the existence of trajectories that always remain on the zero level set of ϕ ^o. Furthermore, any point on the zero level set of ϕ ^o belongs to a characteristics of Eq. 9 emanating from y _s since y _s is the only point in Ω where ϕ ^o is initially zero. Therefore, when the zero level set reaches y for the first time, it implies the existence of a trajectory $\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},t)$ with $\mathbf {X}^{o}_{P}(\mathbf {y}_{\mathbf {s}},0) = \mathbf {y}_{\mathbf {s}}$ that satisfies (1). For this trajectory, (45) holds with an equality, thereby establishing (33). Physically, this means that fastest arrival time at any end point y∈Ω is when the zero level set of ϕ ^o reaches y for the first time, and equivalently that the reachability front $\partial \mathcal {R}(\mathbf {y}_{\mathbf {s}},t)$ coincides with the zero level set of ϕ ^o at time t.

• Let y _f ∈Ω be fixed. From part (3), the optimal trajectory to y _f satisfies $\phi ^{o}(\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t),t) = 0$ for all t≥0. Hence, ${\phi ^{o}_{P}}(t) := \phi ^{o}(\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t),t)$ equals zero for all $0 \le t \le {T}^{\star }(\mathbf {y}_{\mathbf {f}})$. Let us fix a time $0 < t < {T}^{\star }(\mathbf {y}_{\mathbf {f}})$ such that ϕ ^o is differentiable at $(\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t),t)$. The usual chain rule then yields

  $$ 0 = \frac{\text{d} {\phi^{o}_{P}}(t)}{\text{d} t} = \frac{\partial \phi^{o}}{\partial t} + \nabla \phi^{o} \cdot \frac{\text{d} \mathbf{X}^{\star}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t} \, , $$

  (46)

  where the derivatives of ϕ ^o are evaluated at $(\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t),t)$. Since ϕ ^o is assumed to be differentiable at this point, Eq. 9 holds in the classical sense and $\frac {\partial \phi ^{o}}{\partial t} = - F |\nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)| - \mathbf {V}(\mathbf {X}^{\star }_{P},t)\cdot \nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)$. Substituting this in Eq. 46 gives

  $$ \frac{\text{d} \mathbf{X}^{\star}_{P}}{\text{d} t} \cdot \nabla \phi^{o}(\mathbf{X}^{\star}_{P},t) = F |\nabla \phi^{o} (\mathbf{X}^{\star}_{P},t)| + \mathbf{V}(\mathbf{X}^{\star}_{P},t) \cdot \nabla\phi^{o} (\mathbf{X}^{\star}_{P},t) \, . $$

  (47)

  Using Eq. 1,

  $$\begin{array}{cl} &\frac{\text{d} \mathbf{X}^{\star}_{P}}{\text{d} t} \cdot \nabla \phi^{o}(\mathbf{X}^{\star}_{P},t) \\ &= F_{P}^{\star}(t)\, \hat{h}^{\star}(t) \cdot \nabla \phi^{o}(\mathbf{X}^{\star}_{P},t) + \mathbf{V}(\mathbf{X}^{\star}_{P},t) \cdot \nabla\phi^{o} (\mathbf{X}^{\star}_{P},t) \, \\ &\le F |\nabla \phi^{o} (\mathbf{X}^{\star}_{P},t)| + \mathbf{V}(\mathbf{X}^{\star}_{P},t) \cdot \nabla\phi^{o} (\mathbf{X}^{\star}_{P},t) \, , \end{array} $$

  equality holding iff $F_{P}^{\star }(t) = F$ and $\hat {h}^{\star }(t) = \frac {\nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)}{|\nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)|}$, for $|\nabla \phi ^{o} (\mathbf {X}^{\star }_{P},t)| \ne 0$. Using this result in Eq. 47 yields

  $$ \frac{\text{d} \mathbf{X}^{\star}_{P}}{\text{d} t} = F \, \frac{\nabla \phi^{o} (\mathbf{X}^{\star}_{P},t)}{|\nabla \phi^{o} (\mathbf{X}^{\star}_{P},t)|} + \mathbf{V}(\mathbf{X}^{\star}_{P},t) \, . $$

• Under the assumption $F > \sup _{\mathbf {x} \in {\Omega }, t\ge 0} \{ |\mathbf {V}(\mathbf {x},t)| \}$, the start point y _s belongs to the interior of the reachable set $\mathcal {R}(\mathbf {y}_{\mathbf {s}},t)$ for all t>0, i.e., for any t > 0, there exists 𝜖 _t > 0 such that all points $\mathbf {x}^{\prime }$ satisfying $|\mathbf {y}_{\mathbf {s}} - \mathbf {x}^{\prime }| < \epsilon _{t}$ are members of $\mathcal {R}(\mathbf {y}_{\mathbf {s}},t)$. This condition is equivalent to the “small time local controllability” condition discussed in the study of Bardi and Capuzzo-Dolcetta (2008) as a result of which, T ^o is continuous in Ω. See (Bardi and Capuzzo-Dolcetta 2008) for a formal proof of this statement. Let us fix y∈Ω. By definition, T ^o(y) satisfies

  $$ {T}^{o}(\mathbf{y}) = \inf_{h>0} \{ {T}^{o}(\widetilde{\mathbf{y}}) + h\} \, , $$

  (48)

  where $\widetilde {\mathbf {y}} \in {\Omega }$ is a point such that there exists a trajectory $\mathbf {\widetilde {X}}_{P}(\mathbf {y}_{\mathbf {s}},t)$ satisfying (1) and the limiting conditions

  $$ \mathbf{\widetilde{X}}_{P}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\mathbf{\widetilde{y}})) = \mathbf{\widetilde{y}},~~~ \mathbf{\widetilde{X}}_{P}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\mathbf{\widetilde{y}})+h) = \mathbf{y} \, . $$

  (49)

  In order to show that T ^o is a viscosity solution to Eq. 34, we show that it is both a viscosity subsolution and a supersolution to Eq. 34.

  Viscosity subsolution From Definition 3 and Remark 1, ${T}^{o} \in \mathcal {C}({\Omega })$ is a viscosity subsolution to Eq. 34 if at every y∈Ω and for every $\mathcal {C}^{1}$ function $\tau _{s}:{\Omega } \rightarrow \mathbb {R}$ such that $\tau _{s}(\mathbf {y}) = {T}^{o}(\mathbf {y})$ and $\tau _{s} - {T}^{o}$ has a local minima at y,

  $$ F \, |\nabla \tau_{s} (\mathbf{y})| + \mathbf{V}(\mathbf{y}, {T}^{o}(\mathbf{y})) \cdot \nabla \tau_{s}(\mathbf{y}) - 1 \le 0 \, . $$

  (50)

  Since $\tau _{s} \ge {T}^{o}$ in a neighborhood of y, we obtain for h>0 small enough,

  $$\tau_{s}(\mathbf{y}) - \tau_{s}(\widetilde{\mathbf{y}}) \le {T}^{o}(\mathbf{y}) - {T}^{o}(\widetilde{\mathbf{y}}) \, . $$

  Moreover, for this choice of h and the resulting $\widetilde {\mathbf {y}}$, (48) implies

  $${T}^{o}(\mathbf{y}) \le {T}^{o}(\widetilde{\mathbf{y}}) + h \, . $$

  Combining the above two inequalities yields

  $$ \tau_{s}(\mathbf{y}) - \tau_{s}(\widetilde{\mathbf{y}}) \le {T}^{o}(\mathbf{y}) - {T}^{o}(\widetilde{\mathbf{y}}) \le h \, . $$

  (51)

  Since τ _s is differentiable at y, Taylor’s theorem may be used to expand $\tau _{s}(\widetilde {\mathbf {y}})$ near y.

  $$\begin{array}{@{}rcl@{}} \tau_{s}(\widetilde{\mathbf{y}}) &= \tau_{s}(\mathbf{y}) + \nabla \tau_{s}(\mathbf{y}) \cdot (\widetilde{\mathbf{y}} -\mathbf{y}) + \mathbf{o} \left (|\widetilde{\mathbf{y}} - \mathbf{y}| \right ) \, \\ &= \tau_{s}(\mathbf{y}) - {\int}_{{T}^{o}(\widetilde{\mathbf{y}})}^{{T}^{o}(\widetilde{\mathbf{y}})+h} \nabla \tau_{s}(\mathbf{y}) \cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}}{\text{d} t} \, \text{d}t + \mathbf{o} \left (|\widetilde{\mathbf{y}} - \mathbf{y}| \right ).\\ \end{array} $$

  (52)

  Inserting Eq. 52 in Eq. 51 and dividing by h,

  $$\frac{1}{h} {\int}_{{T}^{o}(\widetilde{\mathbf{y}})}^{{T}^{o}(\widetilde{\mathbf{y}})+h} \nabla \tau_{s}(\mathbf{y}) \cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}}{\text{d} t} \, \text{d}t + \frac{\mathbf{o} \left (|\widetilde{\mathbf{y}} - \mathbf{y}| \right )}{h} \le 1 \, . $$

  As h ↓0 and after noting that the second term on the left vanishes under this limit, we obtain

  $$ \nabla \tau_{s}(\mathbf{y}) \cdot \frac{\text{d} \widetilde{\mathbf{X}}_{P}}{\text{d} t}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\mathbf{y})) \le 1 \, . $$

  (53)

  One can see that Eq. 50 is satisfied trivially when |∇τ _s (y)|=0. Thus, we may assume |∇τ _s (y)|≠0. Since (53) holds for any valid choice of $\frac {\text {d} \widetilde {\mathbf {X}}_{P}}{\text {d} t}$, we may choose $\frac {\text {d} \widetilde {\mathbf {X}}_{P}}{\text {d} t} (\mathbf {y}_{\mathbf {s}},{T}^{o}(\mathbf {y}))= F \frac {\nabla \tau _{s}(\mathbf {y})}{|\nabla \tau _{s}(\mathbf {y})|} + \mathbf {V}(\mathbf {y},{T}^{o}(\mathbf {y}))$ to obtain

  $$\begin{array}{cl} &\nabla \tau_{s}(\mathbf{y}) \cdot \left (F \frac{\nabla \tau_{s}(\mathbf{y})}{|\nabla \tau_{s}(\mathbf{y})|} + \mathbf V(\mathbf{y},{T}^{o}(\mathbf{y})) \right ) \\ &= F|\nabla \tau_{s}(\mathbf{y})| + \mathbf{V}(\mathbf{y},{T}^{o}(\mathbf{y})) \cdot \nabla \tau_{s}(\mathbf{y}) \le 1 \, , \end{array} $$

  thereby establishing (50). Therefore, T ^o is a viscosity subsolution to Eq. 34.

  Viscosity supersolution T ^o is a viscosity supersolution to Eq. 34 if at any y∈Ω and for every $\mathcal {C}^{1}$ function $\tau ^{s}:{\Omega } \rightarrow \mathbb {R}$ such that τ ^s(y)=T ^o(y) and τ ^s−T ^o has a local maxima at y,

  $$ F |\nabla \tau^{s} (\mathbf{y})| + \mathbf{V}(\mathbf{y}, {T}^{o}(\mathbf{y})) \cdot \nabla \tau^{s}(\mathbf{y}) - 1 \ge 0 \, . $$

  (54)

  For any 0 < h < T ^o(y), there exists $\widehat {\mathbf y} \in {\Omega }$ satisfying ${T}^{o}(\widehat {\mathbf {y}}) + h = {T}^{o}(\mathbf {y})$ and a trajectory $\widehat {\mathbf X}_{P}(\mathbf {y}_{\mathbf {s}},t)$ satisfying (1) and the limiting conditions

  $$ \widehat{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\widehat{\mathbf y})) = \widehat{\mathbf{y}},~~~ \widehat{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\mathbf{y})) = \mathbf{y} \, . $$

  (55)

  Of course, the optimal trajectory leading to y is a valid choice for $\widehat {\mathbf X}_{P}(\mathbf {y}_{\mathbf {s}},t)$. For h > 0 small enough,

  $$ h = {T}^{o}(\mathbf{y}) - {T}^{o}(\widehat{\mathbf{y}}) \le \tau^{s}(\mathbf{y}) - \tau^{s}(\widehat{\mathbf{y}}) \, . $$

  (56)

  As in the earlier sub-section, we may use Taylor’s theorem to expand $\tau ^{s}(\widehat {\mathbf {y}})$ near y to obtain

  $$\begin{array}{@{}rcl@{}} \tau^{s}(\widehat{\mathbf{y}}) &= \tau^{s}(\mathbf{y}) + \nabla \tau^{s}(\mathbf{y}) \cdot (\widehat{\mathbf{y}} -\mathbf{y}) + \mathbf{o} \left (|\widehat{\mathbf{y}} - \mathbf{y}| \right ) \, \\ &= \tau^{s}(\mathbf{y}) - {\int}_{{T}^{o}(\widehat{\mathbf{y}})}^{{T}^{o}(\widehat{\mathbf{y}})+h} \nabla \tau^{s}(\mathbf{y}) \cdot \frac{\text{d} \widehat{\mathbf{X}}_{P}}{\text{d} t} \, \text{d}t + \mathbf{o} \left (|\widehat{\mathbf{y}} - \mathbf{y}| \right ). \end{array} $$

  (57)

  Inserting (57) in Eq. 56 and dividing by h,

  $$ \frac{1}{h} {\int}_{{T}^{o}(\widehat{\mathbf{y}})}^{{T}^{o}(\widehat{\mathbf{y}})+h} \nabla \tau^{s}(\mathbf{y}) \cdot \frac{\text{d} \widehat{\mathbf{X}}_{P}}{\text{d} t} \, \text{d}t + \frac{\mathbf{o} \left (|\widehat{\mathbf{y}} - \mathbf{y}| \right )}{h} \ge 1 \, . $$

  (58)

  Observe that from Eq. 1,

  $$\nabla \tau^{s}(\mathbf{y}) \cdot \frac{\text{d} \widehat{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\text{d} t} \le F |\nabla \tau^{s}(\mathbf{y})| + \mathbf{V}(\widehat{\mathbf{X}}_{P}(\mathbf{y}_{\mathbf{s}},t),t) \cdot \nabla \tau^{s}(\mathbf{y}) \, . $$

  Taking limits of Eq. 58 as h ↓0 gives

  $$1 \le \nabla \tau^{s}(\mathbf{y}) \cdot \frac{\text{d} \widehat{\mathbf X}_{P}}{\text{d} t}(\mathbf{y}_{\mathbf{s}},{T}^{o}(\mathbf{y})) \le F |\nabla \tau^{s}(\mathbf{y})| + \mathbf{V}(\mathbf{y},{T}^{o}(\mathbf{y})) \cdot \nabla \tau^{s}(\mathbf{y}) \, , $$

  which proves that T ^o is a viscosity supersolution of Eq. 34. Therefore, T ^o is a viscosity solution to Eq. 34. □

Appendix C: Numerical schemes

We now summarize the numerical schemes utilized to discretize and solve (9) and (12). These equations are solved using a finite volume framework implemented in MATLAB. The term |∇ϕ ^o| in (9) is discretized using either a first order (Sethian 1999b; Lolla 2012) or a higher order (Yigit 2011) upwind scheme and V(x,t)⋅∇ϕ ^o is discretized using a second order TVD scheme on a staggered C-grid (Ueckermann and Lermusiaux 2011).

1.1 C.1 Forward level set evolution

We discretize (9) in time using a fractional step method as follows:

$$\begin{array}{@{}rcl@{}} \frac{\bar{\phi} - \phi^{o}(\mathbf{x},t)}{\Delta t/2} &= -F|\nabla \phi^{o}(\mathbf{x},t)| \end{array} $$

(59)

$$\begin{array}{@{}rcl@{}} \frac{\bar{\bar{\phi}} - \bar{\phi}}{\Delta t} = -\mathbf{V}\left (\mathbf{x},t+\frac{\Delta t}{2} \right ) \cdot \nabla \bar{\phi} \end{array} $$

(60)

$$\begin{array}{@{}rcl@{}} \frac{\phi^{o}(\mathbf{x},t+{\Delta} t) - \bar{\bar{\phi}}}{\Delta t/2} &= -F |\nabla \bar{\bar{\phi}}| \end{array} $$

(61)

Equations 59–61 are solved only in the interior nodes of the discretized system. For boundaries that are open inlets/outlets or side walls (i.e., not interior obstacles nor forbidden regions), open boundary conditions are used on ϕ ^o and on the intermediate variables $\bar {\phi }$ and $\bar {\bar {\phi }}$ at each time step. Specifically, a radiation boundary condition with infinite wave speed is assumed, which amounts to an internal zero normal gradient (Neumann) condition, so the boundary values are updated by replacing them with the value of the variable one cell interior to the boundary. Obstacles and forbidden regions in the domain are masked, i.e., (9) is solved only at interior nodes not lying under these regions. For points adjacent to the mask, open boundary conditions are implemented and necessary spatial gradients are evaluated using neighboring nodes that do not lie under the mask. As a result, the value of ϕ ^o under the mask is never used in the computation. We note that in some situations, more complex open boundary conditions could be used as done in regional ocean modeling (Lermusiaux 1997; Haley Jr. and Lermusiaux 2010). We have implemented the narrow-band scheme of Adalsteinsson and Sethian (1995) to solve (59)–(61).

The reachability front $\partial \mathcal {R}_{\phi ^{o}}(t)$ is extracted from the ϕ ^o field at every time step using a contour algorithm. In a 2D problem, the amount of storage required for this is not significant because $\partial \mathcal {R}_{\phi ^{o}}(t)$ is a 1D curve which is numerically represented by a finite number of points. We also note that this contour extraction is not needed: We could simply store the times when the zero contour of ϕ ^o crosses each grid point in order to compute the normals for the backtracking (Yigit 2011).

C.2 Backtracking

(12) is discretized using first order (Lolla 2012) or higher-order (Yigit 2011) time integration schemes. Ideally, it suffices to solve Eq. 9 until the level set front first reaches y _f . However, due to the discrete time steps, a more convenient stopping criterion is the first time, T, when $\phi ^{o}(\mathbf {y}_{\mathbf {f}},T) \le 0$. Due to this, y _f does not lie on the final contour $\partial \mathcal {R}_{\phi ^{o}}(t)$ exactly. Thus, we first project y _f onto $\partial \mathcal {R}_{\phi ^{o}}(t)$. The projected $\hat {\mathbf {n}}_{p}$ is computed as the unit normal to $\partial \mathcal {R}_{\phi ^{o}}(t)$ at the projected point. The discretized form of Eq. 12,

$$ \frac{\mathbf{X}^{\star}_{P}(\mathbf{y}_{\mathbf{s}},t-{\Delta} t) - \mathbf{X}^{\star}_{P}(\mathbf{y}_{\mathbf{s}},t)}{\Delta t} = -\mathbf{V}(\mathbf{X}^{\star}_{P},t) - F\underbrace{\frac{\nabla \phi^{o}(\mathbf{X}^{\star}_{P},t)}{|\nabla \phi^{o}(\mathbf{X}^{\star}_{P},t)|}}_{\hat{\mathbf{n}}_{p}(\mathbf{x},t)} \, , $$

(62)

is marched back in time until we reach a point on the first saved contour and this generates a discrete representation of $\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t)$. Along the way, we project each newly computed trajectory point, $\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t-{\Delta } t)$ onto the corresponding intermediate level set contour (see Lolla (2012)). Instead of performing these projections, one can use the two intermediate discrete level set contours between which an unprojected trajectory point lies to interpolate either the normal $\hat {\mathbf {n}}_{p}$ at the trajectory point or a contour passing through the trajectory point from which $\hat {\mathbf {n}}_{p}$ can be computed. This interpolation should be of sufficiently high order to prevent potential biases that may occur. One can also use a predictor-corrector scheme to compute $\mathbf {X}^{\star }_{P}(\mathbf {y}_{\mathbf {s}},t-{\Delta } t)$ using the normals both at t−Δt and t.

As discussed in Section 4.1 and the Uniqueness remark of Section 3.3, multiple optimal paths exist to end points y _f which lie on shock lines. However, as Eq. 9 is solved numerically, y _f does not lie on shock lines exactly due to discretization errors. In fact, y _f does not even lie exactly on the final level set contour $\partial \mathcal {R}_{\phi ^{o}}(t)$, as mentioned above. Consequently, solving (12) in such cases yields only one of the optimal trajectories, depending on the numerical errors.