Equivalent extensions of Hamilton-Jacobi-Bellman equations on hypersurfaces

TL;DR

This paper introduces a new method to solve Hamilton-Jacobi-Bellman equations on surfaces by extending the problem to a narrow band around the surface, enabling the use of Cartesian grid methods for efficient computation.

Contribution

The authors develop an equivalent formulation of HJB equations on surfaces by extending the control problem to a narrow band, facilitating the use of existing numerical methods.

Findings

Unique viscosity solutions are obtained in the narrow band.

The method allows efficient computation using Cartesian grid techniques.

Applicable to unstructured point clouds sampled from surfaces.

Abstract

We present a new formulation for the computation of solutions of a class of Hamilton Jacobi Bellman (HJB) equations on closed smooth surfaces of co-dimension one. For the class of equations considered in this paper, the viscosity solution of the HJB equation is equivalent to the value function of a corresponding optimal control problem. In this work, we extend the optimal control problem given on the surface to an equivalent one defined in a sufficiently thin narrow band of the co-dimensional one surface. The extension is done appropriately so that the corresponding HJB equation, in the narrow band, has a unique viscosity solution which is identical to the constant normal extension of the value function of the original optimal control problem. With this framework, one can easily use existing (high order) numerical methods developed on Cartesian grids to solve HJB equations on surfaces, with a computational cost that scales with the dimension of the surfaces. This framework also provides a systematic way for solving HJB equations on the unstructured point clouds that are sampled from the surface.

Tables

Table 1. Table 1. Errors and orders of accuracy for the distance function on a sphere, computed by the first order Lax-Friedrichs sweeping algorithm. We use h = 2 / ( N − 1 ) ℎ 2 𝑁 1 h=2/(N-1) and ϵ = 4 ​ h italic-ϵ 4 ℎ \epsilon=4h .

$N$	number of points in $T_{ϵ}$	$L_1$ error	order	$L_{\infty }$ error	order	iterations
101	63302	0.02319		0.03741		32
201	252310	0.01138	1.0345	0.02213	0.76291	49
301	566458	7.5443e-3	1.0158	0.01611	0.78624	65
401	1571974	5.6562e-3	1.0041	0.01280	0.79978	80

Table 2. Table 2. Errors and orders of accuracy for the distance function on a sphere, computed by the third order Lax-Friedrichs sweeping algorithm. We use h = 2 / ( N − 1 ) ℎ 2 𝑁 1 h=2/(N-1) and ϵ = 11 ​ h italic-ϵ 11 ℎ \epsilon=11h .

$N$	number of points in $T_{ϵ}$	$L_1$ error	order	$L_{\infty }$ error	order	iterations
101	183810	5.1129e-4		0.02779		65
201	702626	6.6194e-5	2.9706	0.01137	1.2987	92
301	1566014	2.0373-5	2.9182	7.1199e-3	1.1587	100
401	2776370	8.5992e-6	3.0069	5.1806e-3	1.1085	119

Table 3. Table 3. L 1 subscript 𝐿 1 L_{1} errors for the first order Lax-Friedrichs scheme using different densities of uniformly sampled point cloud representations of the sphere and different noises applied to the point clouds. The noise is applied by adding δ ∗ ( − 1 + 2 ∗ rand ​ ( ) ) / 2 𝛿 1 2 rand 2 \delta*(-1+2*\text{rand}())/2 to each point in the point clouds. We use h = 2 / ( N − 1 ) ℎ 2 𝑁 1 h=2/(N-1) , ϵ = 4 ​ h italic-ϵ 4 ℎ \epsilon=4h .

		$L_1$ errors using ${\Gamma }_K$
$\delta$	N	$K=400$	$K=800$	$K=1600$
0	101	0.02332	0.02317	0.02318
	201	0.01172	0.01144	0.01139
	301	8.0413e-3	7.6406e-3	7.5649e-3
0.001	101	0.02349	0.02321	0.02313
	201	0.01192	0.01153	0.01148
	301	8.2866e-3	7.7936e-3	7.75321e-3
0.005	101	0.02464	0.02364	0.02404
	201	0.01389	0.01339	0.01531
	301	0.01082	0.01265	0.01279

Table 4. Table 4. Comparison of different choices of α 𝛼 \alpha for the sphere on a 101 3 superscript 101 3 101^{3} grid. Here, ϵ = 11 ​ h italic-ϵ 11 ℎ \epsilon=11h and h = 2 / 100 . ℎ 2 100 h=2/100.

$\alpha$	$d\left(𝐳\right)-\left(11-2\sqrt{3}\right)h$	$d\left(𝐳\right)-3h$	$d\left(𝐳\right)$	$d\left(𝐳\right)+3h$
$L_1$ error	0.02298	0.02303	0.02320	0.02349

Table 5. Table 5. Comparison of methods on a 101 3 superscript 101 3 101^{3} grid: (1) Our method using exact closest point mapping and exact singular values and vectors; (2) Our method using closest point mapping sampled from a an equally distributed point cloud of 400 points and finite differences to compute singular values and vectors; (3) method in [ 27 ] . Here ϵ = 2 ​ h 0.7 italic-ϵ 2 superscript ℎ 0.7 \epsilon=2h^{0.7} and h = 2 / 100 ℎ 2 100 h=2/100 .

Method	(1)	(2)	(3)
$L_{\infty }$ error	0.03226	0.03334	0.17072
$L_1$ error	0.07787	0.07667

Equations

$$\displaystyle\min_{\mathbf{a}\in A_{\mathbf{x}}}\bigg{\{}r(\mathbf{x},\mathbf{a})+\nabla_{\Gamma}u(\mathbf{x})\cdot{f}(\mathbf{x},\mathbf{a})\mathbf{a}\bigg{\}}$$

$$u(\mathbf{x})$$

$$T_{ϵ}:=\{\mathbf{z}\in {\mathbb{R}}^n:\underset{\mathbf{x}\in \Gamma }{\min }||\mathbf{x}-\mathbf{z}||)<ϵ\}.$$

$$\overline{H}(\mathbf{z},\nabla v(\mathbf{z}))$$

$$v(\mathbf{z})$$

$$P_{\Gamma }(\mathbf{z}):=\arg \underset{\mathbf{y}\in \Gamma }{\min }||\mathbf{z}-\mathbf{y}||.$$

$$d_{\Gamma }(\mathbf{z})=\{\begin{array}{ll}||\mathbf{z}-P_{\Gamma }\mathbf{z}||& \text{mbox}if\mathbf{z}\in {\overline{\Omega }}^c.\\ -||\mathbf{z}-P_{\Gamma }\mathbf{z}||& \text{mbox}if\mathbf{z}\in \Omega \end{array}$$

$${\Gamma }_{\eta }:=\{\mathbf{z}\in T_{ϵ}|d_{\Gamma }(\mathbf{z})=\eta \}.$$

$$P_{\Gamma }(\mathbf{z})=\mathbf{z}-d_{\Gamma }(\mathbf{z})\nabla d_{\Gamma }(\mathbf{z}).$$

$$\overline{u}(\mathbf{z}):=u(P_{\Gamma }\mathbf{z}).$$

$$A_{\mathbf{x}}=S^{n-1}\cap T_{\mathbf{x}}\Gamma .$$

$$\mathcal{A}=\{\text{mbox}measurablefunctions\mathbf{a}:[0,\infty )\to A\}.$$

$$\frac{d\mathbf{y}}{dt}(t)$$

$$\mathbf{y}(0)$$

$${\mathcal{A}}_{\mathbf{x}}:=\{\mathbf{a}(\cdot )\in \mathcal{A}|{\mathbf{y}}_{\mathbf{x}}(t,\mathbf{a}(t))\in \Gamma \text{mbox}forallt\ge 0\}.$$

$$C(\mathbf{x},\mathbf{a}(\cdot ))={\int }_0^{T}r({\mathbf{y}}_{\mathbf{x}}(t),\mathbf{a}(t))dt+g({\mathbf{y}}_{\mathbf{x}}(T)),$$

$$\{\begin{array}{l}f\text{mbox}andr\text{mbox}areLipschitzcontinuousintheirfirstarguments.\\ g\text{mbox}islowersemicontinuousand{\min }_{\mathcal{T}}g<\infty .\\ \text{mbox}Thereexistsconstants{R}_1,R_2,F_1,\text{mbox}and{F}_2\text{mbox}suchthat:\\ 0<R_1\le r(\mathbf{x},\mathbf{a})\le R_2,\text{mbox}forall\mathbf{x}\in \Gamma \text{mbox}and\mathbf{a}\in A,\\ 0<F_1\le f(\mathbf{x},\mathbf{a})\le F_2,\text{mbox}forall\mathbf{x}\in \Gamma \text{mbox}and\mathbf{a}\in A.\end{array}$$

$$u(\mathbf{x})=\underset{\mathbf{a}(\cdot )\in {\mathcal{A}}_{\mathbf{x}}}{\inf }C(\mathbf{x},\mathbf{a}(\cdot )).$$

$$u(\mathbf{x})=\inf_{\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{x}}}\bigg{\{}\int_{0}^{\tau}r(\mathbf{y}_{\mathbf{x}}(t),\mathbf{a}(t))dt+u(\mathbf{y}_{\mathbf{x}}(\tau))\bigg{\}}.$$

$$\min_{\mathbf{a}\in A_{\mathbf{x}}}\bigg{\{}r(\mathbf{x},\mathbf{a})+\nabla_{\Gamma}u(\mathbf{x})\cdot{f}(\mathbf{x},\mathbf{a})\mathbf{a}\bigg{\}}=0,\;\mathbf{x}\in\Gamma\backslash\mathcal{T}.$$

$$u(\mathbf{x})=g(\mathbf{x}),\mathbf{x}\in \mathcal{T}.$$

$$||{\nabla }_{\Gamma }u(\mathbf{x})||=\frac{r(\mathbf{x})}{f(\mathbf{x})}.$$

$$\overline{H}(\mathbf{z},\nabla v(\mathbf{z}))=0$$

$$v(\mathbf{z})=\overline{u}(\mathbf{z}):=u(P_{\Gamma }\mathbf{z}),$$

$$B(\mathbf{z},\mu )=B_0(\mathbf{z})+\mu B_1(\mathbf{z}),$$

$$B_0(\mathbf{z})$$

$$B_1(\mathbf{z})$$

$${\sigma }_i=1-d_{\Gamma }(\mathbf{z}){\kappa }_i(\mathbf{z})$$

$$\frac{d\mathbf{y}}{dt}(t)$$

$$\mathbf{y}(0)$$

$${\mathcal{A}}_{\mathbf{z}}:=\{\mathbf{a}\in \mathcal{A}|{\mathbf{y}}_{\mathbf{z}}(t,\mathbf{a}(t))\in {\Gamma }_{\eta }\}.$$

$$B({\mathbf{y}}_{\mathbf{z}}(t),\mu )\mathbf{a}(t)=\mathbf{a}(t),$$

$$\overline{C}(\mathbf{z},\mathbf{a}(\cdot ))$$

Full text

Department of Mathematics, The University of Texas at Austin, TX 78712, USA. E-mail: lmartin@math.utexas.edu

Department of Mathematics and Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, TX 78712, USA. E-mail: ytsai@math.utexas.edu

1. Introduction

Hamilton-Jacobi-Bellman equations have many applications in optimal control, seismology, geometrical optics, etc. From the solutions of HJB equations on surfaces, the corresponding characteristics curves can be extracted and used in many applications. Some examples include mesh generation [15], path planning [19, 37], and brain mapping [4, 21]. The equations are fully nonlinear and classical solutions typically do not exist. The unique viscosity solution [10] is often sought after. Sophisticated algorithms have been developed for computing the viscosity solution of HJB equations defined in Euclidean space.

Let $\Omega\subset\mathbb{R}^{n}$ where $n=2$ or $3$ be a bounded and connected open set with smooth closed boundary $\Gamma=\partial\Omega$. Our goal is to compute solutions of the following HJB equation defined on smooth surfaces, with given Dirichlet boundary conditions:

[TABLE]

The set $A_{\mathbf{x}}:=S^{n-1}\cap T_{\mathbf{x}}\Gamma$ is a compact set where $S^{n-1}$ is the $(n-1)$-dimensional unit sphere in $\mathbb{R}^{n}$ and $T_{\mathbf{x}}\Gamma$ is the tangent space of $\Gamma$ at the point $\mathbf{x}\in\Gamma$. We define ${A:=S^{n-1}\cap TM(\Gamma)}$ to be the unit tangent bundle of $\Gamma$ where $TM(\Gamma)$ is the tangent bundle of $\Gamma$. Here $r,f:\Gamma\times A\to\mathbb{R}$, $\mathcal{T}$ is a closed subset of $\Gamma$, and $\nabla_{\Gamma}u$ is the surface gradient on $\Gamma$ [11].

For $\epsilon>0$, define the narrow band of $\Gamma$ by

[TABLE]

In this paper, we shall derive a Hamiltonian, $\overline{H}$, and extensions $\overline{\mathcal{T}}$ and $\overline{g}$ of $\mathcal{T}$ and $g$, respectively, such that the viscosity solution to

[TABLE]

is the constant normal extension of the solution to (1)-(2), for any positive and sufficiently small $\epsilon$.

Our contribution includes a theory for how optional control problems on surfaces can be extended into “equivalent” ones defined in a narrow band of the surface. Depending on whether the problem is isotropic or anisotropic, we must take careful consideration of extending the control space. We show for the anisotropic or most general case, the same control space used in the surface problem must be used in defining the extended problem. However, in the isotropic case the extension is equivalent even when we appropriately extend the control space. After defining the extension of optimal control problems, we then extend the corresponding HJB equation defined on the narrow band of the surface. The main advantage of this approach is that to compute solutions of HJB equations on surfaces, we are able to use Cartesian grids and existing methods (including high order methods) which solve HJB equations in Euclidean space. This allows us to avoid unnecessary patching and triangulation for solving the surface HJB equation. We show that in fact the narrow band can be very thin, i.e., it’s thickness is of order $h$, where $h$ is the grid size.

Before proceeding, we will first define the closest point mapping that is used in the extensions. Define the closest point mapping, $P_{\Gamma}:T_{\epsilon}\to\Gamma$, by

[TABLE]

Here, $\epsilon$ must be chosen so that the closest point mapping is unique. Then we define the signed distance function to $\Gamma$, $d_{\Gamma}:T_{\epsilon}\to\mathbb{R}$, by

[TABLE]

For $|\eta|<\epsilon$, we will define the parallel surface, $\Gamma_{\eta}$, by

[TABLE]

Closest point mappings are easily derived in the context of level set methods [30], and there are a variety of fast and high order methods available to compute them from the distance function to $\Gamma$ [37, 33, 8, 36, 39], i.e.,

[TABLE]

Finally, we define the constant normal extension of function on $\Gamma$, $u:\Gamma\to\mathbb{R}$, to be $\overline{u}:T_{\epsilon}\to\mathbb{R}$ given by

[TABLE]

Again, our goal is to define an extended HJB equation on $T_{\epsilon}$ so that the solution is the constant normal extension to the solution of the HJB equation on the surface.

The rest of the paper is organized as follows: In Section 2, we present the optimal control problem which models controlled motion on a surface and define the related class of HJB equations. We then extend the control problem and HJB equations on the narrow band, $T_{\epsilon}$. In Section 3, we present a special case of the problem when the HJB equation reduces to the Eikonal equation and the distance function on the surface is desired. Finally, we give some numerical simulations in Section 4.

1.1. Previous work: solving HJB equations in Euclidean space

There is extensive work on developing numerical methods to compute solutions to (1) where the domain is a bounded open set in $\mathbb{R}^{n}$. In the view of a optimal control problem, one can obtain a semi-Lagrangian discretization for HJB equations which gives a large system of coupled nonlinear equations. Extensive studies of semi-Lagrangian techniques are found in [12, 3, 13, 14]. Solving the discretized system using fixed-point iterations is expensive, thus fast marching methods (FMMs) and fast sweeping methods (FSMs) were developed.

FMMs are a variant of Djikstra’s method and use a heapsort algorithm to determine the order in which the grid nodes are updated. The algorithm was first developed on Cartesian grids [37, 33] to solve Eikonal equations which are a special case of (1) when the functions $r$ and $f$ are isotropic, i.e., $r(\mathbf{x},\mathbf{a})=r(\mathbf{x})$ and $f(\mathbf{x},\mathbf{a})=f(\mathbf{x})$. The solution to the Eikonal equation is well known to be the distance function to the set $\mathcal{T}$. FMMs have complexity $\mathcal{O}(N\log N)$, where $N$ is the total number of unknowns [31, 37]. Based on the Godunov upwind numerical scheme developed to solve Eikonal equations on triangulated domains in [5], FMMs were extended to acute triangulated meshes in [17] . However, FMMs do not directly apply to the anisotropic case. FMMs were applied to compute solutions of a class of “axis-aligned” anisotropic HJ equations in [1] on orthogonal grids, and ordered upwind methods (OUMs) were developed to solve a class of general HJB equations on any triangulated surface [34]. In [6, 7], related PDEs are solved using similar fast methods.

FSMs were developed in [36, 40] to compute solutions of a class of strictly convex HJ equations, including Eikonal equation, on Cartesian grids. The FSM is an iterative method that relies on an upwind discretization of the PDE and Gauss-Seidel style updates with different orderings of grid nodes. The idea is to avoid using heapsort and, instead, update the grid nodes in different orderings (sweeps). In each sweep, characteristics that go in similar directions are propagated automatically. In [22], a FSM was developed for Eikonal equations using a discontinuous Galerkin (DG) local solver for computing the distance function. For more general Hamilton-Jacobi equations, one may use a Lax-Friedrichs type numerical Hamiltonian [16] on Cartesian grids. The resulting Lax-Friedrichs sweeping scheme is typically more diffusive and requires more iterations. FSMs have complexity $\mathcal{O}(N)$ with the caveat that the constant hidden in the notation may be large depending on the characteristics of the PDE.

1.2. Previous work: solving HJB equations on surfaces

Several versions of FMMs and FSMs can be applied to computing the distance function on surfaces. As already mentioned, the FMM in [17] only works for solving Eikonal equations on acute triangulated surfaces. In [35], FMMs equation were extended to solve Eikonal equations on parametric surfaces where the discretization is done on a Cartesian grid in the parametric plane. The FSM in [36] and the FMM in [29] can compute the distance function on surfaces that are defined as the graph of a smooth function. Another approach to computing the distance function on surfaces is to solve a corresponding evolutionary HJ equation where the distance function is the steady state solution [8]. However, with the addition of the time variable this method is not computationally optimal.

OUMs are applicable to more general anisotropic HJB equations on triangulated surfaces. However, if the surface is given in implicit form, triangulation may be unnecessary. We want to avoid triangulation and have the ability to handle general surface representations. Our motivation is to derive equivalent extended HJB equations on a narrow band of the surface such that a variety of meshes and methods can be used to solve the equivalent equations in the narrow band. In particular, the methods developed for Cartesian grids (surveyed in Section 1.1) can then be used.

In [27], the idea of extending the surface by a small offset, $\epsilon$, to a narrow band is used to compute distance functions on surface. On the narrow band, an extended distance function is computed, but it is not the constant normal extension of the distance function on the surface. Therefore, the formulation introduces an analytical error. It is proven that the error between the two distance functions is controlled by $\epsilon$, the width of the narrow band. In order for the method to converge as the grid size goes to zero, it is required that $\epsilon=Ch^{\gamma}$ where $h$ is the grid size and $\gamma\in(0,1).$ The advantage of our framework is that there no analytical error introduced. Our formulation is convergent with respect to the grid size, $h$, and for narrow bands of order $h$. We are also able to perform high order computations for the distance function on surfaces.

The idea of extending functions or differential operators defined on surfaces, via closest point mapping, to the embedding Euclidean space, is used in [32, 24, 9] to compute solutions to parabolic and elliptic PDEs. Our approach is inspired by the work in [9] where variational integrals defined originally on surfaces are extended to ones defined in the narrow band, and the Euler-Lagrange equations of the extended problem are solved by standard numerical methods. Due to the way the extensions are defined, the resulting solution to each equation is automatically the constant normal extension of the solution to the variational problem on the surface.

1.3. Computing PDEs “on point clouds”

Point clouds arise in many applications including facial recognition, manufacturing, medicine, and geosciences and can be acquired easily through modern sensing devices such as laser scanners and cell phones. A typical approach, particularly predominant in the computer graphics community, is to create triangular meshes from the data points and seek to define partial derivatives in certain suitable ways. Among myriad of papers using triangular surface meshes, we mention [5], an early paper which is the most related to the topic of this paper. Another school of approaches aim at solving differential equations on surfaces, using only a finite set of points sampled from the surfaces, without globally reconstructing the surface.

In [28], the algorithm developed in [27] is extended to compute distance functions on surfaces represented as point clouds. In [23], a framework for computing solutions to PDEs on point clouds based on a local approximation of the manifold was presented. A local mesh algorithm was developed in [20] that solves PDEs on point clouds where any of the existing methods valid on triangular meshes discussed in Section 1.2 can be used to compute the solutions to HJB equations once the local mesh is determined.

Our new proposed non-parametric formulation also provides a convenient way to solve equations “on point clouds.” That is we can solve HJB equations which are defined on a smooth closed surface, but the only available information about the surface is a finite set of points that are reasonably distributed over the surface. The discussion of the implementation of our method applied to point clouds is in Section 4.1.4. There, it is assumed that there is no noise perturbing the point cloud in the surface normal directions, and that the point clouds are evenly distributed over the surfaces. The generalization of the proposed algorithm to more general point clouds is the subject of another paper.

2. Static Hamilton-Jacobi-Bellman equations on surfaces

2.1. Modeling controlled motion on a surface

First, we define the optimal control problem from which the HJB equations can be derived. For each $\mathbf{x}\in\Gamma$, let $A_{\mathbf{x}}$ be set of all unit tangent vectors of $\Gamma$ at $\mathbf{x}$, i.e.,

[TABLE]

Note each $A_{\mathbf{x}}$ is a compact set in $\mathbb{R}^{n}$ for each $\mathbf{x}\in\Gamma$. Define $A:=\bigcup_{\mathbf{x}\in\Gamma}A_{\mathbf{x}}$, then $A=S^{n-1}\cap TM(\Gamma)$ where $TM(\Gamma)$ is the tangent bundle of $\Gamma$. The set of admissible controls is given by:

[TABLE]

We are interested in the trajectories governed by the dynamical system:

[TABLE]

where ${f}:\Gamma\times A\to\mathbb{R}$ represents the velocity and $\mathbf{a}(t)\in A_{\mathbf{y}(t)}$, which ensures that $\mathbf{y}(t)\in\Gamma$ for all $t\geq 0$. Denote the solution to (5)-(6) (which exists under the assumption (A1) below) by $\mathbf{y}_{\mathbf{x}}(t,\mathbf{a}(t))$. For each $\mathbf{x}\in\Gamma$ define the following subset of the admissible controls:

[TABLE]

Then the requirement that $\mathbf{a}(t)\in A_{\mathbf{y}(t)}$ for all $t\geq 0$ in the dynamical system (5)-(6) is equivalent to $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{x}}.$ Refer to $\mathbf{y}(\cdot)$ that satisfy (5)-(6) where $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{x}}$ as an admissible path on $\Gamma$. For the rest of the paper, we will suppress writing the dependence of $\mathbf{y}_{\mathbf{x}}(t,\mathbf{a}(t))$ on $\mathbf{a}(\cdot)$ when there is no confusion in the context.

Let $r:\Gamma\times A\to\mathbb{R}$ denote a running cost per unit time and $g:\mathcal{T}\to\mathbb{R}$ be an exit time penalty when the state reaches a closed target set $\mathcal{T}\subset\Gamma$. The total cost associated with initial state $\mathbf{x}$ and control $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{x}}$ to reach $\mathcal{T}$ is given by

[TABLE]

where $T=\min\{t\;|\;\mathbf{y}_{\mathbf{x}}(t)\in\mathcal{T}\}.$

We have the following assumptions on $f,r,$ and $g$:

[TABLE]

The value function $u:\Gamma\to\mathbb{R}$ is defined to be the minimal total cost to $\mathcal{T}$ starting at $\mathbf{x}$ and is given by

[TABLE]

Under the assumption (A1), an optimal control does not necessarily exist. If for the given $f$ and $r$, we have that the set $V(\mathbf{x})=\{f(\mathbf{x},\mathbf{a})\mathbf{a}/r(\mathbf{x},\mathbf{a})\;|\;\mathbf{a}\in A_{\mathbf{x}}\}$ is strictly convex for each $\mathbf{x}\in\Gamma$, then an optimal control is guaranteed [2]. This property is trivially satisfied in the case of isotropic $f$ and $r$. Regardless of whether an optimal control exists, we can still derive the corresponding HJB equation.

The dynamic programming principle [2] for the value function states that for sufficiently small $\tau>0$,

[TABLE]

The corresponding Hamilton-Jacobi-Bellman equation on $\Gamma$ is

[TABLE]

The boundary condition is

[TABLE]

Note that in (10) we can take the minimum since $A_{\mathbf{x}}$ is compact.

In general, classical solutions of (10)-(11) do not exist, and the unique viscosity solution is sought after [10]. A detailed discussion of viscosity solutions of Hamilton-Jacobi equations on manifolds can be found in [25]. Under the assumption (A1), it is a classic result that the value function $u$ coincides with the unique viscosity solution of (10)-(11) [2]. In the case of isotropic running cost and speed, then (10) reduces to the Eikonal equation on the surface:

[TABLE]

We will expand more on the Eikonal equation in the next section.

2.2. Extension to $T_{\epsilon}$

We derive a Hamiltonian, $\overline{H}$, on $T_{\epsilon}$ such that if $v:T_{\epsilon}\to\mathbb{R}$ is the unique viscosity solution of

[TABLE]

with boundary conditions defined appropriately, then

[TABLE]

where $u$ is the viscosity solution to (10)-(11).

First, we extend the optimal control problem on the surface to one in the narrow band, $T_{\epsilon}$. We start by defining a tensor that will be essential in formulating an equivalent optimal control problem on $T_{\epsilon}$.

Definition 2.1**.**

Let $\mathbf{t}_{1}(\mathbf{z})$ and $\mathbf{t}_{2}(\mathbf{z})$ be two orthonormal tangent vectors corresponding to the directions that yield the principle curvatures of $\Gamma_{\eta}$ at a point $\mathbf{z}$. Let $\mathbf{n(\mathbf{z})}$ be the unit normal vector to the tangent plane at $\mathbf{z}$. Let $\sigma_{1}(\mathbf{z})$ and $\sigma_{2}(\mathbf{z})$ be the singular values of $P_{\Gamma}^{\prime}(\mathbf{z})$, the Jacobian matrix of $P_{\Gamma}$ at $\mathbf{z}$. Then for any real number $\mu$, define the tensor $B$ by

[TABLE]

where

[TABLE]

In the above definition, $B_{0}$ corresponds to a weighted orthogonal projection onto the tangent plane of $\Gamma_{\eta}$ at $\mathbf{z}$, and $B_{1}$ is the projection along the normal of $\Gamma_{\eta}$ passing through $\mathbf{z}$. In [18], it is proven that

[TABLE]

where $\kappa_{1}(\mathbf{z})$ and $\kappa_{2}(\mathbf{z})$ are the principal curvatures of the parallel surface $\Gamma_{\eta}$. Note that when $\mathbf{z}\in\Gamma$, $\sigma_{1}=\sigma_{2}=1.$ For the sake of notation, we will suppress the dependence of the tangent vectors, normal vectors, and singular values on the point $\mathbf{z}\in T_{\epsilon}$.

Assuming $\epsilon$ small enough, we have that the tangent planes coincide for $\mathbf{z}\in T_{\epsilon}$ and $\mathbf{x}\in\Gamma$ if $P_{\Gamma}\mathbf{z}=\mathbf{x}$. Thus, for the extended optimal control problem, the set of compact control values at each $\mathbf{z}\in T_{\epsilon}$ is $A_{P_{\Gamma}\mathbf{z}}.$ Then $A$ and $\mathcal{A}$ are as in Section 2.1. Now, consider the following extended dynamical system in $T_{\epsilon}$:

[TABLE]

Here $\overline{f}:T_{\epsilon}\times A\to\mathbb{R}$ is given by $\overline{f}(\mathbf{z},\mathbf{a})=f(P_{\Gamma}\mathbf{z},\mathbf{a}),$ and we have that $\mathbf{a}(t)\in A_{\mathbf{y}(t)}$ for all $t\geq 0.$ This last requirement ensures that the trajectory remains in $\Gamma_{\eta}$ for all $t\geq 0$ where $\eta=d_{\Gamma}(\mathbf{z})$. The inclusion of the tensor $B(\mathbf{y}(t),\mu)$ in the dynamical system above adjusts the velocity according to the curvature of the surface. This has the implication that equivalent paths will have the same total cost.

Denote the solution to (14)-(15) (which exists under the assumption (B1) below) by $\mathbf{y}_{\mathbf{z}}(t,\mathbf{a}(t))$. For each $\mathbf{z}\in T_{\epsilon}$ with $d_{\Gamma}(\mathbf{z})=\eta$, we have the following subset of admissible controls:

[TABLE]

Then the requirement that $\mathbf{a}(t)\in A_{\mathbf{y}_{\mathbf{z}}(t)}$ for all $t\geq 0$ in the dynamical system (14)-(15) is equivalent to $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{z}}.$ Since $A_{\mathbf{x}}=A_{P_{\Gamma}\mathbf{z}}$, we have $\mathcal{A}_{\mathbf{z}}=\mathcal{A}_{P_{\Gamma}\mathbf{z}}$. Refer to $\mathbf{y}(\cdot)$ that satisfy (14)-(15) where $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{z}}$ as an admissible path on $T_{\epsilon}$. Again, for the rest of the paper we will write suppress writing the dependence of $y_{\mathbf{z}}(t,\mathbf{a}(t))$ on $\mathbf{a}(\cdot)$. There is no confusion between the notation for admissible paths on $\Gamma$ and $T_{\epsilon}$. The initial point will indicate which dynamical system the path solves.

Note that if the initial point $\mathbf{z}\in\Gamma$, $\eta=0$ and

[TABLE]

since $\mathbf{a}(t)\in T_{\mathbf{y}_{\mathbf{z}}(t)}$ and $\sigma_{1}=\sigma_{2}=1.$ Thus, (14)-(15) reduces to the dynamical system on $\Gamma$ given in (5)-(6), and an admissible path on $T_{\epsilon}$ with $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{z}}$ is also an admissible path on $\Gamma$ with the same $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{z}}$.

Now, let $\overline{\mathcal{T}}=\{\mathbf{z}\in T_{\epsilon}\;|\;P_{\Gamma}\mathbf{z}\in\mathcal{T}\}$. Let the running cost and exit time penalty on $T_{\epsilon}$ be given by $\overline{r}:T_{\epsilon}\times A\to\mathbb{R}$ where $\overline{r}(\mathbf{z},\mathbf{a})=r(P_{\Gamma}\mathbf{z},\mathbf{a})$ and $\overline{g}:T_{\epsilon}\to\mathbb{R}$ where $\overline{g}(\mathbf{z})=g(P_{\Gamma}\mathbf{z})$, respectively. Then the total cost associated with initial state $\mathbf{z}$ and control ${\mathbf{a}}(\cdot)\in\mathcal{A}_{\mathbf{z}}$ to reach $\overline{\mathcal{T}}$ is given by

[TABLE]

where $\overline{T}=\min\{t\;|\;\mathbf{{y}_{\mathbf{z}}}(t)\in\overline{\mathcal{T}}\}$. Note that if $\mathbf{z}\in\Gamma$, then $\overline{C}(\mathbf{z},\mathbf{a}(\cdot))=C(\mathbf{z},\mathbf{a}(\cdot)).$ The value function $v:T_{\epsilon}\to\mathbb{R}$ is the minimal total cost to $\overline{\mathcal{T}}$ starting at $\mathbf{z}$ and is given by

[TABLE]

We have that the assumption A1 implies:

[TABLE]

Again, while (B1) does not guarantee the existence of the optimal control, the assumption does ensure that the value function and corresponding viscosity solution of the HJB equation coincide. If we also assume the set $V(\mathbf{x})$ is strictly convex, then $\overline{V}(\mathbf{z})=\{\overline{f}(\mathbf{z},\mathbf{a})B(\mathbf{z},\mu)\mathbf{a}/\overline{r}(\mathbf{z},\mathbf{a})\;|\;\mathbf{a}\in A_{\mathbf{z}}\}$ is strictly convex for each $z\in T_{\epsilon}$. Thus, if the optimal control exists in the problem on $\Gamma$, it also exists for extended problem on $T_{\epsilon}$.

Before deriving the corresponding HJB equation, we prove that the value function, $v$, is the constant normal extension of the value function of the optimal control problem defined on $\Gamma$. The following theorem states that two admissible paths, $\mathbf{y}_{\mathbf{z}_{1}}(\cdot)$ and $\mathbf{y}_{\mathbf{z}_{2}}(\cdot)$, on $T_{\epsilon}$ with equivalent initial points in $T_{\epsilon}$ and equivalent controls stay equivalent for all time.

Theorem 2.1**.**

Suppose that $\mathbf{z}_{1},\mathbf{z}_{2}\in T_{\epsilon}$ such that $P_{\Gamma}\mathbf{z}_{1}=P_{\Gamma}\mathbf{z}_{2}$. Let ${\mathbf{a}_{1}}(\cdot)\in{\mathcal{A}_{\mathbf{z}_{1}}}$ and ${\mathbf{a}_{2}}(\cdot)\in{\mathcal{A}_{\mathbf{z}_{2}}}$ such that ${\mathbf{a}_{1}}\equiv{\mathbf{a}_{2}}$. Suppose ${{\mathbf{y}_{\mathbf{z}_{1}}}:[0,\infty)\to T_{\epsilon}}$ solves (14) with ${\mathbf{y}_{\mathbf{z}_{1}}}(0)=\mathbf{z}_{1}$ and ${{\mathbf{y}_{\mathbf{z}_{2}}}:[0,\infty)\to T_{\epsilon}}$ solves (14) with ${\mathbf{y}_{\mathbf{z}_{2}}}(0)=\mathbf{z}_{2}$. Then $P_{\Gamma}{\mathbf{y}_{\mathbf{z}_{1}}}(t)=P_{\Gamma}{\mathbf{y}_{\mathbf{z}_{2}}}(t)$ for all $t\geq 0$.

Proof.

Let $\mathbf{y}_{1}:[0,\infty)\to\Gamma$ and $\mathbf{y}_{2}:[0,\infty)\to\Gamma$ be given by

[TABLE]

and

[TABLE]

for all $t\geq 0$, respectively. We need to show that $\mathbf{y}_{1}(t)=\mathbf{y}_{2}(t)$ for all $t\geq 0$. We compute $\mathbf{y}_{1}^{\prime}$. We have the following singular value decomposition for $P_{\Gamma}^{\prime}$: $P_{\Gamma}^{\prime}=U\Sigma U^{T}$ where

[TABLE]

Since $\mathbf{t}_{1},$ $\mathbf{t}_{2}$, and $\mathbf{n}$ are orthonormal, for $i,j=1,2$,

[TABLE]

and for $i=1,2$

[TABLE]

Using the above, we have

[TABLE]

The last equality is true since $\mathbf{t}_{1}\otimes\mathbf{t}_{1}+\mathbf{t}_{2}\otimes\mathbf{t}_{2}$ is the orthogonal projection of $\mathbf{a}_{1}(t)$ onto the tangent space at $\mathbf{y}_{\mathbf{z}_{1}}(t)$ and $\mathbf{a}_{1}(t)\in A_{\mathbf{y}_{\mathbf{z}_{1}}(t)}\subset T_{\mathbf{y}_{\mathbf{z}_{1}}(t)}\Gamma$. Thus $\mathbf{y}_{1}$ is an admissible path on $\Gamma$ where $\mathbf{y}_{1}(0)=\mathbf{x}=P_{\Gamma}\mathbf{z}_{1}=P_{\Gamma}\mathbf{z}_{2}$ and $\mathbf{a}_{1}(\cdot)\in\mathcal{A}_{\mathbf{x}}=\mathcal{A}_{P_{\Gamma}\mathbf{z}_{1}}.$ By the same reasoning $\mathbf{y}_{2}$ also an admissible path with $\mathbf{y}_{2}(0)=\mathbf{x}$ and $\mathbf{a}_{2}(\cdot)\in\mathcal{A}_{\mathbf{x}}=\mathcal{A}_{P_{\Gamma}\mathbf{z}_{2}}.$ We assume (A1). Thus, the solution to (5)-(6) is unique. Since $\mathbf{a}_{1}\equiv\mathbf{a}_{2}$, $\mathbf{y}_{1}(t)=\mathbf{y}_{2}(t)$ for all $t\geq 0$. ∎

Now, it easily follows that the value function defined on $T_{\epsilon}$ is the constant normal extension of the value function defined on $\Gamma.$

Corollary 2.2**.**

If $u:\Gamma\to\mathbb{R}$ is the value function on $\Gamma$ given by (8), then

[TABLE]

where $v$ is the value function on $T_{\epsilon}$ given by (16).

Proof.

Let $\mathbf{z}\in T_{\epsilon}$ and $\mathbf{a}(\cdot)\in\mathcal{A}_{\mathbf{z}}$. Then $\mathbf{a}(\cdot)\in\mathcal{A}_{P_{\Gamma}\mathbf{z}}$ since $\mathcal{A}_{\mathbf{z}}=\mathcal{A}_{P_{\Gamma}\mathbf{z}}$. Then we have

[TABLE]

and

[TABLE]

where ${\mathbf{y}}_{\mathbf{z}}(\cdot)$ and ${\mathbf{y}}_{P_{\Gamma}\mathbf{z}}(\cdot)$ are admissible paths on $T_{\epsilon}$ with ${\mathbf{y}}_{\mathbf{z}}(0)=\mathbf{z}$, ${\mathbf{y}}_{P_{\Gamma}\mathbf{z}}(0)=P_{\Gamma}\mathbf{z}$, respectively. Recall, that $\mathbf{y}_{P_{\Gamma}\mathbf{z}}$ is also an admissible path on $\Gamma$ so that $\overline{C}(P_{\Gamma}\mathbf{z},{\mathbf{a}}(\cdot))={C}(P_{\Gamma}\mathbf{z},{\mathbf{a}}(\cdot))$.

Theorem 2.1 implies that $P_{\Gamma}\mathbf{y}_{\mathbf{z}}(t)=\mathbf{y}_{P_{\Gamma}\mathbf{z}}(t)$ for all $t\geq 0$. Since

[TABLE]

we have

[TABLE]

Therefore, $\overline{T}=T$, and

[TABLE]

Finally, we have

[TABLE]

∎

Now that we have showed that the value function on the narrow band is equivalent to the value function on the surface, we define the Hamilton-Jacobi-Bellman equation on $T_{\epsilon}$ associated the value function $v$. Again, the dynamic programming principle states that for sufficiently small $\tau>0$, we have

[TABLE]

The corresponding Hamilton-Jacobi-Bellman equation on $T_{\epsilon}$ is

[TABLE]

The boundary condition is

[TABLE]

The assumption (B1) implies that the value function coincides with the unique viscosity solution [2]. Corollary 2.2 implies that the viscosity solution of (18)-(19) is given by

[TABLE]

where $u$ is the viscosity solution of (10)-(11) on $T_{\epsilon}$. An interesting result is that we can prove (20) without Corollary 2.2, i.e., using only the derived HJB equations and properties of viscosity solutions.

We begin with the following lemma that relates the surface gradients for parallel surfaces at equivalent points:

Lemma 2.3**.**

Assume $\epsilon$ is small enough so that $P_{\Gamma}$ is differentiable, and let $|\eta|<\epsilon$. Given $\phi_{0}\in C^{1}(\Gamma)$, define $\phi_{\eta}\in C^{1}(\Gamma_{\eta})$ by $\phi_{\eta}(\mathbf{z}):=\phi_{0}(P_{\Gamma}\mathbf{z})$. Then for $\mathbf{z}\in\Gamma_{\eta}$,

[TABLE]

Proof.

Let $\overline{\phi}_{0}:T_{\epsilon}\to\mathbb{R}$ be the constant normal extension of $\phi_{0}.$ Then $\overline{\phi}_{0}\in C^{1}(T_{\epsilon})$ and {\left.\kern-1.2pt\overline{\phi}_{0}\vphantom{\big{|}}\right|_{\Gamma_{\eta}}}=\phi_{\eta}. From Theorem A.1 in [18], we have that ${\nabla_{\Gamma}\phi_{0}}(P_{\Gamma}\mathbf{z})=B(\mathbf{z},\mu)\nabla\overline{\phi}_{0}(\mathbf{z})$. Since $\overline{\phi}_{0}$ is the constant normal extension of $\phi_{0}$, $\nabla\overline{\phi}_{0}(\mathbf{z})\in T_{\mathbf{z}}\Gamma_{\eta}$. Thus, if $\mathbf{z}\in\Gamma_{\eta}$, $\nabla\overline{\phi}_{0}(\mathbf{z})=\nabla_{\Gamma_{\eta}}\phi_{\eta}(\mathbf{z})$, and (21) holds. ∎

The following theorem proves (20) directly using the definition of viscosity solutions.

Theorem 2.4**.**

If $u:\Gamma\to\mathbb{R}$ is the viscosity solution to (10)-(11), then the constant normal extension of $u$, $\overline{u}:T_{\epsilon}\to\mathbb{R},$ is the viscosity solution to (18)-(19).

Proof.

We have $u(\mathbf{x})=g(\mathbf{x})$ for $\mathbf{x}\in\mathcal{T}$ and $P_{\Gamma}\mathbf{z}\in\mathcal{T}$ if $\mathbf{z}\in\overline{\mathcal{T}}$. Then for $\mathbf{z}\in\overline{\mathcal{T}}$,

[TABLE]

Thus, $\overline{u}$ satisfies the boundary conditions (19). Next, we will show that $\overline{u}$ is a viscosity subsolution of (18).

Assume $\mathbf{z}_{0}\in T_{\epsilon}$ and $\phi\in C^{1}(T_{\epsilon})$ such that $\overline{u}-\phi$ has a strict local maximum at $\mathbf{z}_{0}$ and $(\overline{u}-\phi)(\mathbf{z}_{0})=0.$ We need to show that

[TABLE]

Let $\eta=d_{\Gamma}(\mathbf{z}_{0})$ and $\phi_{0}:\Gamma\to\mathbb{R}$ be the restriction to $\Gamma$ of the normal extension of {\left.\kern-1.2pt\phi\vphantom{\big{|}}\right|_{\Gamma_{\eta}}}. Then we have that $\phi_{0}\in C^{1}(\Gamma)$ and \phi_{\eta}:={\left.\kern-1.2pt\phi\vphantom{\big{|}}\right|_{\Gamma_{\eta}}}\in C^{1}(\Gamma_{\eta}) such that ${\phi_{\eta}(\mathbf{z})=\phi_{0}(P_{\Gamma}\mathbf{z})}$ for $\mathbf{z}\in\Gamma_{\eta}$.

First, we will prove

[TABLE]

Since $\mathbf{a}\in A_{\mathbf{z}_{0}}\subset T_{\mathbf{z}_{0}}\Gamma$, we have $\mathbf{a}\cdot\mathbf{n}=0$ and

[TABLE]

Now, $\nabla_{\Gamma_{\eta}}\phi_{\eta}(\mathbf{z}_{0})=(I-\mathbf{n}\otimes\mathbf{n})\nabla\phi(\mathbf{z}_{0})$. Therefore, (23) holds.

Lemma 2.3 implies

[TABLE]

Let $\mathbf{x}_{0}=P_{\Gamma}\mathbf{z}_{0}.$ Then using the fact that $B(\mathbf{z}_{0},\mu)$ is symmetric, (23), and (24), we have

[TABLE]

Now, $\mathbf{x}_{0}\in\Gamma,\phi_{0}\in C^{1}(\Gamma)$ such that $u-\phi_{0}$ has a strict local maximum at $\mathbf{x}_{0}$ and

[TABLE]

Since $u$ is a viscosity subsolution of (10),

[TABLE]

Thus, (22) holds. The same argument can be applied to show that $\overline{u}$ is a viscosity supersolution of (10). ∎

In the next section, we will present a special case of the above formulations for when the speed and running cost functions are isotropic. The HJB equation on $\Gamma$ reduces to the Eikonal equation. A noteworthy revelation is that if we ensure that the extended speed and cost functions are isotropic as well, the control values in the extended optimal control problem do not have to be restricted to be in the tangent bundle of $\Gamma$, i.e., the extended set of controls is $\widetilde{A}=S^{n-1}$. In this setup, the value function is the constant normal extension of the value function on $\Gamma.$ Before we proceeding, we present a counter example to show that $\overline{u}$ is not necessarily the value function of the extended control problem if the extended speed function $\overline{f}$ is anisotropic and $\widetilde{A}=S^{n-1}$.

Example

Consider the following control problems on $\Gamma$ and $T_{\epsilon}$: Let

[TABLE]

and

[TABLE]

Let $\epsilon>\sqrt{3}/3$, then

[TABLE]

and

[TABLE]

For the control problem on $\Gamma$, the set of controls is $A=\{(-1,0),(1,0)\}$ and $\mathcal{A}$ is the usual set of admissible controls. Define $f:\Gamma\times A\to\mathbb{R}$ by $f(\mathbf{x},\mathbf{a})=\sqrt{3}$ and suppose that we have unit running cost, i.e., $r(\mathbf{x},\mathbf{a})=1$. Then it is clear that $u(0)=\sqrt{3}/3$ which is the minimum travel time from $\mathbf{x}=0$ to $\mathcal{T}$ while traveling at speed $\sqrt{3}$.

Now for the extended control problem, let the control values be the set $\widetilde{A}=S^{1}$. Then the set of admissible controls is given by

[TABLE]

Define $\overline{f}:T_{\epsilon}\times\widetilde{\mathcal{A}}\to\mathbb{R}$ such that $\overline{f}(\mathbf{z},\widetilde{\mathbf{a}})=\overline{f}(\widetilde{\mathbf{a}})$ and the set $\{\overline{f}(\widetilde{\mathbf{a}})\widetilde{\mathbf{a}}\;|\;\widetilde{\mathbf{a}}\in\widetilde{A}\}$ is the ellipse given by $x^{2}+\frac{y^{2}}{9}=1$ rotated at the origin clockwise by the angle $\pi/3$. Then $\overline{f}(\mathbf{a})=\sqrt{3}$ if $\mathbf{a}\in A$ so that $\overline{f}(\mathbf{z},\mathbf{a})=\overline{f}(P_{\Gamma}\mathbf{z},\mathbf{a})$ for $\mathbf{a}\in A$. Again, assume unit running cost.

Let $v$ be the value function for this extended control problem. If $v=\overline{u}$, then we should have $v((0,0))=u(0)$. However, note that $\overline{f}(\mathbf{z},(\sqrt{3}/2,1/2))=3$ and the Euclidean distance between $(0,0)$ and $(1,\sqrt{3}/3)$ is $2\sqrt{3}/3$. Thus,

[TABLE]

when $\widetilde{a}(\cdot)\equiv(\sqrt{3}/2,1/2)$. By the definition of the value function we must have

[TABLE]

Therefore, $\overline{u}$ cannot be the value function for the extended control problem.

3. The Isotropic Case

We derive the optimal control problems and HJB equations corresponding to the case when the speed and cost functions are isotropic. As mentioned at the end of Section 2.2, we will show that value function on $T_{\epsilon}$ is the normal extension of the value function on $\Gamma$ even when we do not restrict the controls to the tangent space in the extended optimal control problem.

3.1. The controlled dynamics on $\Gamma$ for isotropic speed and cost functions

The optimal control problem formulation is exactly the same as in Section 2.1, and (10) reduces to the Eikonal equation on the surface:

[TABLE]

We have the HJB equation

[TABLE]

and the boundary condition is

[TABLE]

Assuming (A1), the value function $u$ coincides with the unique viscosity solution of (26)-(27). We also have that $V(\mathbf{x})$ is strictly convex for all $\mathbf{x}\in\Gamma$. Thus, an optimal control is guaranteed.

3.2. Extension to $T_{\epsilon}$ in the isotropic case

We derive the extended Hamiltonian, $\overline{H}$, on $T_{\epsilon}$ from the extended isotropic optimal control problem such that if $v:T_{\epsilon}\to\mathbb{R}$ is the unique viscosity solution of

[TABLE]

with boundary conditions defined appropriately, then

[TABLE]

where $u$ is the solution to (26)-(27).

The notations $\mathbf{t}_{1}$, $\mathbf{t}_{2}$, $\mathbf{n}$, $\sigma_{1}$, and $\sigma_{2}$ are as in Section 2.2. The extended set of compact control values is $\widetilde{A}=S^{n-1}$, and the set of admissible controls is given by:

[TABLE]

This is where the formulation differs from the previous section. We do not restrict the controls to belong to the tangent bundle of $\Gamma$. We extend the speed, cost, and exit time penalty functions as above: $\overline{f}(\mathbf{z})=f(P_{\Gamma}\mathbf{z})$, $\overline{r}(\mathbf{z})=r(P_{\Gamma}\mathbf{z})$, and $\overline{g}(\mathbf{z})=g(P_{\Gamma}\mathbf{z})$. The dynamical system is given by

[TABLE]

where $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$. Admissible paths on $T_{\epsilon}$ with extended controls are solutions to (28)-(29) denoted by $y_{\mathbf{z}}(t)$. Note that the admissible paths are not restricted to the parallel surface, $\Gamma_{\eta}$, to which the initial point $\mathbf{z}$ belongs.

The total cost function associated to the initial state $\mathbf{z}\in T_{\epsilon}$ and control $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$ is the same as in Section 2.2. The value function $v:T_{\epsilon}\to\mathbb{R}$ is the minimal total cost to $\overline{\mathcal{T}}$ starting at $\mathbf{z}$ and is given by

[TABLE]

We have that (A1) implies the following is true:

[TABLE]

We have that (C1) implies the value function, $v$, coincides with the unique viscosity solution of (26)-(27) [2]. We also have that $\overline{V}(\mathbf{z})=\{\overline{f}(\mathbf{z})\widetilde{\mathbf{a}}/\overline{r}(\mathbf{z})\;|\;\widetilde{\mathbf{a}}\in\widetilde{A}\}$ is strictly convex for all $\mathbf{z}\in\Gamma$. Thus, an optimal control is guaranteed.

We now prove that the value function, $v$, is a constant extension of a function on $\Gamma$. However, unlike in the previous section, it is not immediate from the following proofs that $v=\overline{u}$, where $u$ is the value function on $\Gamma$.

Theorem 3.1**.**

Given $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$ and $\mathbf{z}_{1},\mathbf{z}_{2}\in T_{\epsilon}$ such that $P_{\Gamma}\mathbf{z}_{1}=P_{\Gamma}\mathbf{z}_{2}$, let

[TABLE]

and

[TABLE]

solve (28) with ${\mathbf{y}_{\mathbf{z}_{1}}}(0)=\mathbf{z}_{1}$ and ${\mathbf{y}_{\mathbf{z}_{2}}}(0)=\mathbf{z}_{2}$, respectively. Then $P_{\Gamma}\mathbf{y}_{\mathbf{z}_{1}}(t)=P_{\Gamma}{\mathbf{y}}_{\mathbf{z}_{2}}(t)$ for all $t\geq 0$.

The proof of Theorem 3.1 is analogous to the proof of Theorem 2.1 except the control $\widetilde{\mathbf{a}}(\cdot)$ belongs to the extended control space, $\widetilde{\mathcal{A}}$. A consequence of the control belonging to the extended space is that now the projected paths $P_{\Gamma}\mathbf{y}_{\mathbf{z}}(\cdot)$ may not solve (5). This is due to the fact that $(\mathbf{t}_{1}\otimes\mathbf{t}_{1}+\mathbf{t}_{2}\otimes\mathbf{t}_{2})\widetilde{\mathbf{a}}(t)\not=\widetilde{\mathbf{a}}(t)$ if $\widetilde{\mathbf{a}}(t)\not\in T_{P_{\Gamma}\mathbf{y}_{\mathbf{z}}(t)}$ where $\mathbf{t}_{1}$ and $\mathbf{t}_{2}$ are the basis tangent vectors of $T_{P_{\Gamma}\mathbf{y}_{\mathbf{z}}(t)}.$ We now prove that the extended value function in the isotropic case is constant along the normals of $\Gamma$.

Corollary 3.2**.**

If $P_{\Gamma}\mathbf{z}_{1}=P_{\Gamma}\mathbf{z}_{2}$ for $\mathbf{z}_{1},\mathbf{z}_{2}\in T_{\epsilon}$, then $v(\mathbf{z}_{1})=v(\mathbf{z}_{2})$.

Proof.

Let $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$. We have

[TABLE]

and

[TABLE]

where ${\mathbf{y}}_{\mathbf{z}_{1}}(\cdot)$ and ${\mathbf{y}}_{\mathbf{z}_{2}}(\cdot)$ are admissible paths on $T_{\epsilon}$ with extended control, $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$, and ${{\mathbf{y}}_{\mathbf{z}_{1}}(0)=\mathbf{z}_{1}}$, ${\mathbf{y}}_{\mathbf{z}_{2}}(0)=\mathbf{z}_{2}$ where ${P_{\Gamma}\mathbf{z}_{1}=P_{\Gamma}\mathbf{z}_{2}}$.

Theorem 3.1 implies that $P_{\Gamma}\mathbf{y}_{\mathbf{z}_{1}}(t)=P_{\Gamma}\mathbf{y}_{\mathbf{z}_{2}}(t)$ for all $t\geq 0$. Since

[TABLE]

we have

[TABLE]

Therefore, $\overline{T}_{1}=\overline{T}_{2}$ and $\overline{C}(\mathbf{z}_{1},\widetilde{\mathbf{a}}(\cdot))=\overline{C}(\mathbf{z}_{2},\widetilde{\mathbf{a}}(\cdot))$ for $\widetilde{\mathbf{a}}(\cdot)\in\widetilde{\mathcal{A}}$. Now,

[TABLE]

∎

In the above proof, it is not immediate that we can take the infimum over $\mathcal{A}\subset\widetilde{\mathcal{A}},$ which would imply that $v$ is the normal extension of $u$, the value function on $\Gamma$. We will use the corresponding HJB equations and Corollary 3.2 to show that $v=\overline{u}$.

We define the HJB equation on $T_{\epsilon}$ associated the value function $v$. Again, the dynamic programming principle states that for sufficiently small $\tau>0$ we have

[TABLE]

We get the Hamilton-Jacobi-Bellman equation on $T_{\epsilon}$:

[TABLE]

Since the speed and cost functions are isotropic, (32) reduces to an anisotropic Eikonal equation:

[TABLE]

The boundary condition is

[TABLE]

Since the value function coincides with the unique viscosity solution, Corollary 3.2 implies that the viscosity solution of (33)-(34) is a constant along the normals of $\Gamma$, i.e., there is a function $w:\Gamma\to\mathbb{R}$ such that for $\mathbf{z}\in T_{\epsilon}$,

[TABLE]

We now prove that $w=u$ where $u$ is the viscosity solution of (26)-(27).

Theorem 3.3**.**

If $u:\Gamma\to\mathbb{R}$ is the viscosity solution to (26)-(27), then the constant normal extension of $u$, $\overline{u}:T_{\epsilon}\to\mathbb{R},$ is the viscosity solution to (33)-(34).

Proof.

Just as in Theorem 2.4, $\overline{u}$ satisfies the boundary conditions (34). Assume for contradiction that $\overline{u}$ is not the viscosity solution to (33). But from Corollary 3.2 we have that viscosity solution of (33)-(34), $v=\overline{w}$ for some $w\in C^{1}(\Gamma).$ Therefore, if $\overline{u}$ is not the unique viscosity solution then we have $u\not=w$. For contradiction, we will show that $w$ is the viscosity solution of (26) and thus, $u=w$.

Let $\mathbf{x}_{0}\in\Gamma$ and $\phi\in C^{1}(\Gamma)$ such that $w-\phi$ has a local maximum at $\mathbf{x}_{0}$. Consider the normal extensions of $w$ and $\phi$, then $\overline{w}-\overline{\phi}$ has a local max at $\mathbf{x}_{0}$. Since $\overline{w}$ is a viscosity subsolution to (33) and $\overline{\phi}\in C^{1}(T_{\epsilon})$, we have that

[TABLE]

Since $\mathbf{x}_{0}\in\Gamma$, $\sigma_{1}=\sigma_{2}=1$ and $\nabla\overline{\phi}(\mathbf{x}_{0})\in T_{\mathbf{x}_{0}}\Gamma$. Therefore,

[TABLE]

and we have that $B(\mathbf{x}_{0},\mu)\nabla\overline{\phi}(\mathbf{x}_{0})=\nabla\overline{\phi}(\mathbf{x}_{0})=\nabla_{\Gamma}\phi(\mathbf{x}_{0})$. Thus,

[TABLE]

and $w$ is a viscosity subsolution. The same argument can be applied to show that $w$ is a viscosity supersolution. Therefore, $w$ is the viscosity solution to (26) and hence $w=u$. ∎

To summarize, we have shown that we may appropriately extend the control space in the case of isotropic speed and running costs. This has the implication that even though the admissible paths on $T_{\epsilon}$ have the “choice” to leave the tangent space of $\Gamma$, the optimal paths in the extended control problem remain on or parallel to $\Gamma.$

4. Numerical implementation and simulations

4.1. Numerical setup

In the implementation of our new framework, we use a uniform Cartesian grids and denote the grid step size by $h$. $T_{\epsilon}^{h}$ denotes the part of the grid in the narrow band of radius $\epsilon$ around the given surface.

We use the Lax-Friedrichs fast sweeping methods in [16] and a high order version in [39] for our simulations. It is also possible to an upwind scheme generalized from [36]. Unless we mention otherwise, we use the standard first order finite differencing in the approximation of the partial derivatives for the Lax-Friedrichs numerical Hamiltonians.111 All code used to produce the numerical simulations can be found at https://github.com/lindsmart/MartinTsaiExtHJB.

We present several examples of our new formulation performed on surfaces of co-dimension one in three dimensions. For the first three examples, we compute the solution to the Eikonal equation on the surface. As shown in Section 3, the equivalent equation on $T_{\epsilon}$ is

[TABLE]

Again,

[TABLE]

Since the desired solution, $v$, has been proven to be constant along the normals of the surface, $\Gamma$, $\mathbf{n}\otimes\mathbf{n}\nabla v(\mathbf{z})=0$. Therefore, $\mu$ can be any real number. We let $\mu=1$ in all of our computations.

Approximating the solution to (37)-(38) requires the computation of the singular values and vectors of the derivative of the closest point mapping, $P_{\Gamma}^{\prime}$. We defer the discussion of these approximations to Section 4.1.4. In the last example, we apply our framework to solve an HJB equation with an anisotropic speed function.

4.1.1. Boundary closure

We note that the analytical formulation of the HJB equations does not require boundary conditions on $\partial T_{\epsilon}$. However, since we are using Cartesian grids we must take careful consideration of the discretization near the boundary, $\partial T_{\epsilon}^{h}$. When approximating the partial derivatives of $v$, a neighboring grid node may lie outside of $T_{\epsilon}^{h}$. We will call these ghost nodes. In our implementation, we provide a boundary closure procedure to enforce the fact that the solution is a constant along normal function. For each ghost node, we perform the following procedure:

1. (1)

Project the node into the narrow band. 2. (2)

The value at the ghost node is then calculated by interpolating grid values surrounding the projected point. Formally, we must use an interpolation scheme of order higher to the discretization of the PDE on $T_{\epsilon}^{h}$.

Next, we describe how the boundary closure procedure affects the numerical accuracy. Let $\mathbf{z}_{i}$ be a ghost node, and $\mathbf{z}^{\alpha}_{i}=\mathbf{z}_{i}-\alpha\mathbf{n}_{i}$ be the projection of the ghost node into $T_{\epsilon}$. Denote the solution at $\mathbf{z}_{i}^{\alpha}$ by $v_{i}^{\alpha}:=v(\mathbf{z}_{i}^{\alpha})$. Suppose that we use a first order scheme to discretize the HJB equation on $T_{\epsilon}^{h}$. Then the interpolation used in the approximation of $v_{i}^{\alpha}$ should yield at least second order in $h$ accuracy, in order to maintain formally the first order accuracy. This is due to the amplification by a factor of $h^{-1}$ in the finite difference scheme errors of the approximation of the values at the neighboring grid nodes of $\mathbf{z}_{i}$.

4.1.2. Depth of projected points and bounding the thickness of the narrow band

When choosing $\alpha$ for the the projected ghost node discussed above, a necessary condition is that the nearby surrounding nodes used to interpolate the value of $v_{i}^{\alpha}$ must be in $T_{\epsilon}^{h}$. In our numerical simulations, we use cubic interpolation to approximate $v({\bf z}^{\alpha})$. Thus, it requires $4^{3}$ nearby grid nodes for three dimensional cases. It can be shown that if $|d_{\Gamma}(\mathbf{z}_{i})-\alpha|\leq\epsilon-2\sqrt{n}h$, where $n$ is the dimension. Then the inner nodes used to interpolate the value at each ghost node are in $T_{\epsilon}^{h}$. We suggest to use a projected node as close to the boundary node as possible. Therefore in our numerical simulations, we choose $\alpha=d_{\Gamma}(\mathbf{z}_{i})-(\epsilon-2\sqrt{n}h)$. We discuss the effect of the choice of $\alpha$ in Section 4.2.

The maximal value of $\epsilon$ is restricted by the curvatures of $\Gamma$. We must have $\epsilon$ smaller than the reach of the surface, which depends on the reciprocal of the maximal curvatures. The minimal value of $\epsilon$ is determined by the finite difference stencils in approximating the partial derivatives of the PDE and the boundary closure procedure. Both stencils require that the narrow band is sufficiently “thick” relative to the mesh size $h$ with $\epsilon>2\sqrt{n}h$.

The resulting numerical discretization of the HBJ equation on $T_{\epsilon}^{h}$ is convergent for any $\epsilon$ in the interval described above; i.e. it is convergent for very thin $\mathcal{O}(h)$ narrow bands as well as for thicker one whose widths are independent of $h$. We remark that this convergent regime is very different that requiring $\epsilon\sim\mathcal{O}(h^{p}),0<p<1$, in the related work of [18], dealing with singular integrals.

4.1.3. Approximation errors

The error computed by the proposed algorithm can be written as the sum of errors corresponding to different approximations:

[TABLE]

where

• •

$E_{model}$ relates to how the surface PDE is approximated in the tubular neighborhood by an extended PDE and the boundary conditions. Our main contribution is in deriving the extension for which $E_{model}\equiv 0$ for any $\epsilon$ that is smaller than the reach of the surface;

• •

$E_{\Delta}$ corresponds to the numerical error for approximating the extended PDE problem. Our extension allows for higher order (in $h$) methods to be applied;

• •

$E_{\Gamma}$ relates to errors in approximating the surface and its geometric properties. This error corresponds to approximation of the closest point mapping to the surface and its curvatures. Since these quantities are computed on based on finite differencing using the grid, the error depends on $h$;

• •

$E_{L_{p}}$ is the numerical error in approximating the $L_{p}$ norm of the computed errors. Our formulation allows for a very accurate approximation to this error term.

4.1.4. Application to solving HJB on point clouds

We now describe how we compute the closest point mapping when the surface is represented as a point cloud. (See Figure 1.) Denote the point cloud by $\Gamma_{K}\subset\Gamma$, with $K$ indicating the total number of points in the set. Here, we shall assume that the $K$ points uniformly distribute on the surface. The purpose of the exposition here is to point out a promising application of the method to point cloud data. A full analysis of the algorithm in this area warrants a separate paper.

We approximate $P_{\Gamma}$ using the following strategy:

1. (1)

For each $\mathbf{z}_{i}\in T_{\epsilon}^{h}$, estimate $\Gamma_{K}$ locally on the grid nodes. Call the local surface associated to the point $\mathbf{z}_{i},$ $\Gamma_{K}(\mathbf{z}_{i})$. 2. (2)

Compute the closest point of $\mathbf{z}_{i}$ on the local surface, $\Gamma_{K}(\mathbf{z}_{i}).$

In our implementation, we use biquadratic interpolation at each $\mathbf{z}_{i}$ to compute $\Gamma_{N}(\mathbf{z}_{i})$. Another option to locally estimate the surface is least squares as in [23]. We then use Newton’s method to obtain the closest point of each $\mathbf{z}_{i}$ on $\Gamma_{N}(\mathbf{z}_{i}).$ Denote the closest point mappings by $P_{\Gamma_{N}}^{\textrm{int}}(\mathbf{z}_{i})$. Then on each grid node $\mathbf{z}_{i}$, the equation is discretized as usual, with $P_{\Gamma}^{\prime}(\mathbf{z}_{i})$ being approximated by finite differences of $P^{\textrm{int}}_{\Gamma_{N}}(\mathbf{z}_{i}).$ We use a fourth-order finite difference scheme to estimate $P_{\Gamma}^{\prime}$ and singular value decomposition to obtain the singular values and singular vectors of $P_{\Gamma}^{\prime}$. We show the effect of this approximation procedure in Section 4.2.

In any case, the term, $E_{\Gamma}$, in (39) now depends on $1/K$, which corresponds to the uniform spacing of the points. Therefore, with a fixed point cloud, refinement of Cartesian grids eventually will not further improve the accuracy in the numerical solutions when compared to the solution on the idealized smooth surface sampled by the point cloud. In Section 4.3, we present some results revealing the effect of different surface sampling density.

If noise is introduced to the point cloud, and if the amount of perturbation in the original ”surface normals” is significant, more sophisticated surface fitting is required; e.g. adopting a similar strategy proposed in [23]. We will investigate this important issue in a future manuscript. Nevertheless, in Section 4.3, we present some results revealing the effect of noise in the point cloud.

4.2. Example 1.1

First, we present a numerical convergence study. We approximate the solution to (37)-(38), where $\Gamma$ is a sphere centered at $(0,0,0),$ with radius, $r_{0}=0.5$, and $\mathcal{T}$ is a point on the sphere. The exact singular values and vectors in $B(\mathbf{z},1)$ are used. In the case of a sphere

[TABLE]

where $r(\mathbf{z})$ is the radius of the sphere going though $\mathbf{z}$ centered at $(0,0,0)$. Figure 2 shows the distance function to a point on the sphere for two different view points.

Since the solution is constant along the normals of $\Gamma$, we can easily estimate the $L_{1}$ error. Using the formulation in [18], the estimation of the $L_{1}$ error is given by the following formula:

[TABLE]

where $u$ is the exact solution and $u^{h}$ is the approximated solution on $\Gamma$ using the approximate solution, $v_{i}^{h}$, of the extended equation with grid size $h$. Here,

[TABLE]

where $A=7.513931532835806$ and

[TABLE]

where $\sigma_{1},\sigma_{2}$ are given in (40).

Because all the characteristics of the solution emanate from the source point on the sphere, the errors at the source point will propagate through the solution in the entire computational domain. Therefore, we initialize the solution near $\overline{\mathcal{T}}$ with the exact solution if $u(\mathbf{z})<0.2$ The $L_{1}$ and $L_{\infty}$ errors are reported in Table 1. We see that first order convergence is achieved in the $L_{1}$ error. However, because of the singularity of the solution at the point opposite the source point on the sphere, the order in the $L_{\infty}$ error is less than one. The results verify that the boundary closure procedure does not influence the overall order of the scheme.

We report the number of iterations it takes for the $L_{\infty}$ norm of the difference in the solution between successive iterations is less than $1.0e-13$. One iteration includes one sweep of the grid in each of the eight sweeping directions. The timing for the numerical experiments scale linearly to the number of grid nodes in the narrow band $T_{\epsilon}$ and the number of iterations used in the Lax-Friedrichs scheme. The computations were performed on a MacBook laptop with a 1.6 GHz Intel CPU on a single core using Julia. One Lax-Friedrichs iteration took on average 0.34260 seconds for the case $N=101$ in Table 1. We note that our code was not optimized for efficiency.

Another advantage of our setup is that we can use existing high order methods to compute the solution to HJB equations on surfaces. In Table 2, we present the $L_{1}$ and $L_{\infty}$ errors and orders for the third order method given in [39]. We initialize the solution near $\overline{\mathcal{T}}$ to the exact solution if $u(\mathbf{z})<0.2.$ We see that third order is achieved in the $L_{1}$ error. Again, because of the singularity in the solution at the point opposite the source point on the sphere, the order in the $L_{\infty}$ error is first order.

4.3. Example 1.2

In this example, we compute the distance function on the sphere when the sphere is represented as a point cloud. We show the effect of our method when using different uniform samplings of the sphere.222The uniform samplings of the sphere were computed using the code provided at https://github.com/AntonSemechko/S2-Sampling-Toolbox. The results are tabulated in Table 3 when the sampling of the sphere has $400,800,$ and $1600$ points. We again initialize the solution by exact solutions around the same neighborhood of $\overline{\mathcal{T}}$ as in Section 4.2. By comparing the $L_{1}$ error from Table 1 for $N=101$ and the $L_{1}$ error in Table 3 for $N=101$ and no noise, it appears that the dominating source of error is the numerical discretization of the PDE. As the grid is refined, we see that it appears the error related to approximating the surface begins to dominate. This is apparent by comparing the $L_{1}$ error from Table 1 for $N=301$ and the $L_{1}$ error from Table 3 for $N=301$ with no noise. Notice as the point cloud becomes more dense and the accuracy of the surface approximation increases, it appears the dominating error is again due to the numerical discretization of the PDE.

In Table 3, we also consider when noise is introduced to the point clouds of different sampling densities. The noise is applied by adding $\delta*(-1+2*\text{rand}())/2$ to each point in the point clouds, where $\delta=0.001,0.005$ and rand() is a uniformly distributed random number in $(0,1)$. We can see if the noise is too large the error begins to deteriorate and we need a more sophisticated surface fitting algorithm in our method.

4.4. Example 1.3

Next, we study the affect that the depth of projected ghost nodes has on the overall error for the distance function on the sphere. The set up is the same as in the previous section. Here, we choose a $101^{3}$ sized grid with $\epsilon=11h$ where $h=2/100$. We estimate the $L_{1}$ error for when the boundary closure procedure is carried out at four different depths. Recall that a ghost node, $\mathbf{z}_{i}$, is outside of $T_{\epsilon}^{h}$ and is projected into $T_{\epsilon}$ along the normal of the surface at that point, i.e.,

[TABLE]

Table 4 shows the $L_{1}$ error for the depths:

[TABLE]

We can see that the error the solution is not very sensitive to the choice of depth. In our simulations, we chose $\alpha=\epsilon-2\sqrt{3}h$.

4.5. Example 1.4

Finally, we compare the following: (1) our method on the sphere, using exact singular values and vectors as in the set up for Table 1; (2) our method with the point cloud procedure described in 4.1.4 and finite differences to estimate $P_{\Gamma}^{\prime}$; (3) the method in [27]. Since method (3) requires $\epsilon=Ch^{\gamma}$ where $C=\sqrt{3}$ and $\gamma\in(0,1)$, to fairly compare the error of the methods, we increase the width of the narrow band to $h=2h^{0.7},$ which is used in the convergence analysis of (3) in [27]. We initialize the experiments using methods (1), (2), and (3) with exact values given in a box of length $h$ around the source point/points in the initial sets. We also are not able to compute the $L_{1}$ error for method (3) since the solution is not constant along the normals of $\Gamma$. Therefore, we can only report the $L_{\infty}$ error, and we use trilinear interpolation to approximate the solutions on $\Gamma$ for (3). The errors are reported in Table 5. We can also see that estimating $P_{\Gamma}$ from a point cloud and $P_{\Gamma}^{\prime}$ from finite differences does not greatly influence the $L_{1}$ error of the method.

4.6. Example 2

In Figures 3, 4, and 5, the distance function to a source point is shown on various surfaces. The contours shown are equally spaced and parallel. The solutions for the torus and bunny are computed on a $201^{3}$ grid, and a $101^{3}$ grid is used for the elephant, mask, human skull, and dinosaur skull. All computations use a narrow band width of $\epsilon=4h$ and point cloud representations for the surfaces. The number of points in the point clouds of the torus333The point cloud for the torus was generated using the standard parametrization of a torus. and bunny444The point cloud for the Stanford bunny is generated from a refinement of the triangulated Stanford bunny from https://casual-effects.com/data/ [26]. are 178,350 and 228,096, respectively. For the elephant, mask, human skull, and dinosaur skull555The point clouds for the surfaces in Figure 5 were generated from the triangulations downloaded at https://www.myminifactory.com/scantheworld/., we used point clouds with 65,292, 1,199,988, 234,618, and 139,491 points, respectively. In Figure 3, a cross section of the solution is displayed to show that the solution is indeed constant along the normals of the surface. Our new framework also allows us to sort point clouds. Figure 3 displays level “belts” of the point cloud, i.e., points in the point cloud whose distance from the source point lies in given interval.

Another advantage of our formulation is that when we compute the characteristics, known as geodesics, of the Eikonal equation via the extended Eikonal equation on the narrowband, the paths remain on the surface if the initial point lies on $\Gamma$. Since the solution to (37)-(38) is constant along the normals of $\Gamma$, the gradient always belongs to the tangent spaces of $\Gamma$ or its parallel surfaces, $\Gamma_{\eta}$. We have that $B(\mathbf{x},\mu)\mathbf{a}=\mathbf{a}$ for $\mathbf{x}\in\Gamma$ and $\mathbf{a}\in A_{\mathbf{x}}=T_{\mathbf{x}}\Gamma\cap S^{n-1}$ since $\sigma_{1}=\sigma_{2}=1$. Thus, if the initial point is $\mathbf{z}\in\Gamma$, the geodesic on the surface can be obtained by solving the dynamical system

[TABLE]

In Figure 6, we consider the case when we have two source points and display the solution and some geodesics. Recall that in Theorem 2.1 we showed that paths with equivalent initial points remain equivalent for all time. We show a visualization of this property in Figure 7. In the next section we will compare geodesics of the HJB equation given different anisotropic speed functions.

4.7. Example 3

Finally, we test implement the framework on an anisotropic speed function. We solve

[TABLE]

where $f:\Gamma\times A\to\mathbb{R}$ is a curvature based speed function on the surface. We use the speed function proposed in [38]. The speeds are fast on low curvature areas of the surface and slow on high curvature areas. This means that the corresponding “shortest” paths traverse areas of the lowest curvatures. These paths are typically longer than the geodesics of the surface.

The normal curvature of $\Gamma$ at $\mathbf{x}$ in the direction $\mathbf{a}$ is given by

[TABLE]

where $\kappa_{1}$ and $\kappa_{2}$ are the principal curvatures of $\Gamma$ at $\mathbf{x}$. Then the curvature-minimizing speed function is given by

[TABLE]

where $b$ is a positive constant. When $b=0$, (43) reduces Eikonal equation on the surface. A larger value of $b$ corresponds to a greater difference in speeds between areas of low and high curvature. We display the solutions for varying $b$ values to a source point on the bunny in Figure 8.

In this example, the set $V(\mathbf{x})=\{f(\mathbf{x},\mathbf{a})\mathbf{a}\;|\;\mathbf{a}\in A_{\mathbf{z}}\}$ is not neccesarily convex for all $\mathbf{x}$ in $\Gamma.$ Therefore, an optimal control may not exist. However, we can still extract suboptimal paths, called anisotropic geodesics, whose total cost is arbitrarily close to the value function at the starting point. Just as in the isotropic case, the anisotropic geodesics computed from the extended HJB equation on the narrow band will lie on the surface since $B(\mathbf{x},\mathbf{a})\mathbf{a}=\mathbf{a}$ for $\mathbf{x}\in\Gamma$. The paths can easily be extracted because we compute the minimizing control at each grid node when solving (43)-(44). Once we have the optimal (or suboptimal) control values, $\mathbf{a}^{*}(\cdot)$, we then solve dynamical system

[TABLE]

We plot three anisotropic geodesics for each $b$ value in Figure 8. We can see as $b$ increases the Euclidean distances of the paths are longer, and the paths start to seek out the narrow valleys of the bunny.

5. Summary and conclusion

In this paper, we presented a new formulation to compute solutions of a class of HJB equations on smooth hypersurfaces. We extend the HJB equation’s associated optimal control problem from the surface to an equivalent problem defined in a sufficiently ”thin” narrow band around the surface, in the embedding Euclidean space. The extension was done so that the resulting value function is the constant normal extension of the value function defined in the optimal control problem on the surface. We presented the formulations for the general anisotropic equation and showed that the viscosity solution of the HJB equation on the narrow band is the constant normal extension of the viscosity solution on the surface, independently of the optimal control problems. We also presented the isotropic case and showed there is no need to restrict the control space in order to have an equivalent formulation. The proposed approach is independent of surface representation and can be used to compute and define optimal control problems on uniformly distributed point clouds sampled from some smooth surface. It is also clear that our proposed extension approach can be applied in to time dependent equations arising from the finite horizon control problems. Together with [9], our extension approach provides a good framework for solving to mean field games to high order accuracy on complicated and non-parametrized surfaces.

With this new framework, we are able to use Cartesian grids and the existing methods for computing solutions to HJB equations in Euclidean space on a very thin narrow band to solve HJB equations on surfaces coupled with a simple boundary closure procedure. Our numerical examples verify that the boundary closure procedure does not influence the overall order of the method. We also show that our formulation allows one to easily solve the surface HJB equations to high order accuracy.

Acknowledgment

The authors are supported partially by National Science Foundation Grant DMS-1720171. Tsai also thanks National Center for Theoretical Study, Taipei for hosting his visits, in which some of the ideas presented in this paper originated.