Equivalent extensions of Hamilton-Jacobi-Bellman equations on hypersurfaces
Lindsay Martin Richard Tsai

Equivalent extensions of HJB equations on hypersurfaces]Equivalent extensions of Hamilton-Jacobi-Bellman equations on hypersurfaces



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




[2010]Primary 70H20,49J20, 35R01 58E10, 65N06

[
    Richard Tsai
    December 30, 2023
=====================




 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.




§ 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 <cit.>, path planning <cit.>, and brain mapping <cit.>. 
The equations are fully nonlinear and classical solutions typically do not exist. The unique viscosity solution <cit.> is often sought after.
Sophisticated algorithms have been developed for computing the viscosity solution of HJB equations defined in Euclidean space.

Let Ω⊂ℝ^n where n=2 or 3 be a bounded and connected open set with smooth closed boundary Γ=∂Ω. 
Our goal is to compute solutions of the following HJB equation defined on smooth surfaces, with given Dirichlet boundary conditions:

    min_𝐚∈ A_𝐱{ r(𝐱,𝐚)+∇_Γ u(𝐱) ·f(𝐱,𝐚)𝐚}   =0,    𝐱∈Γ\𝒯
    
    u(𝐱)   =g(𝐱),   𝐱∈𝒯.

The set A_𝐱:=S^n-1∩ T_𝐱Γ is a compact set where S^n-1 is the (n-1)-dimensional unit sphere in ℝ^n and T_𝐱Γ is the tangent space of Γ at the point 𝐱∈Γ. We define A:=S^n-1∩ TM(Γ) to be the unit tangent bundle of Γ where TM(Γ) is the tangent bundle of Γ. Here r,f:Γ× A →ℝ, 𝒯 is a closed subset of Γ, and ∇_Γ u is the surface gradient on Γ <cit.>. 

For ϵ>0, define the narrow band of Γ by

    T_ϵ : = {𝐳∈ℝ^n : min_𝐱∈Γ ||𝐱-𝐳||)< ϵ}.

In this paper, we shall derive a Hamiltonian, H, and extensions 𝒯 and g of 𝒯 and g, respectively, such that the viscosity solution to 

    H(𝐳,∇ v (𝐳))   =0,    𝐳∈ T_ϵ\𝒯
    
    v(𝐳)   =g(𝐳),   𝐳∈𝒯

is the constant normal extension of the solution to (<ref>)-(<ref>), for any positive and sufficiently small ϵ.

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_Γ: T_ϵ→Γ, by

    P_Γ(𝐳) := min_𝐲∈Γ ||𝐳-𝐲||.
  Here, ϵ must be chosen so that the closest point mapping is unique.
Then we define the signed distance function to Γ, d_Γ:T_ϵ→ℝ,  by 

    d_Γ(𝐳)=   ||𝐳-P_Γ𝐳||      𝐳∈Ω^c.
    
     -||𝐳-P_Γ𝐳||      𝐳∈Ω

For |η|<ϵ, we will define the parallel surface, Γ_η, by 

    Γ_η:={𝐳∈ T_ϵ  |   d_Γ(𝐳)=η}.

Closest point mappings are easily derived in the context of level set methods <cit.>, and there are a variety of fast and high order methods available to compute them from the distance function to Γ <cit.>, i.e., 

    P_Γ(𝐳)=𝐳-d_Γ(𝐳) ∇ d_Γ(𝐳).



Finally, we define the constant normal extension of function on Γ, u:Γ→ℝ, to be u:T_ϵ→ℝ given by 
    u(𝐳):=u(P_Γ𝐳).
 Again, our goal is to define an extended HJB equation on T_ϵ  
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 <ref>, 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_ϵ. In Section <ref>, 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 <ref>.





 §.§ Previous work: solving HJB equations in Euclidean space
 


There is extensive work on developing numerical methods to compute solutions to (<ref>) where the domain is a bounded open set in ℝ^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 <cit.>. 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 <cit.> to solve Eikonal equations which are a special case of (<ref>) when the functions r and f are isotropic, i.e., r(𝐱,𝐚)=r(𝐱) and f(𝐱,𝐚)=f(𝐱). The solution to the Eikonal equation is well known to be the distance function to the set 𝒯. FMMs have complexity 𝒪(N log N), where N is the total number of unknowns <cit.>.  Based on the  Godunov upwind numerical scheme developed to solve Eikonal equations on triangulated domains in <cit.>,  FMMs were extended to acute triangulated meshes in <cit.> . 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 <cit.> on orthogonal grids, and ordered upwind methods (OUMs) were developed to solve a class of general HJB equations on any triangulated surface <cit.>. 
In <cit.>, related PDEs are solved using similar fast methods.

FSMs were developed in <cit.> 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 <cit.>, 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 <cit.> on Cartesian grids. The resulting Lax-Friedrichs sweeping scheme is typically more diffusive and requires more iterations. 
FSMs have complexity 𝒪(N) with the caveat that the constant hidden in the notation may be large depending on the characteristics of the PDE. 




 §.§ 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 <cit.> only works for solving Eikonal equations on acute triangulated surfaces. In <cit.>, 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 <cit.> and the FMM in <cit.>

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 <cit.>. 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 <ref>) can then be used. 


In <cit.>, the idea of extending the surface by a small offset, ϵ, 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 ϵ, the width of the narrow band. In order for the method to converge as the grid size goes to zero, it is required that ϵ=Ch^γ where h is the grid size and γ∈ (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 <cit.> to compute solutions to parabolic and elliptic PDEs. Our approach is inspired by the work in <cit.> 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.  





 §.§ 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 <cit.>, 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 <cit.>, the algorithm developed in <cit.> is extended to compute distance functions on surfaces represented as point clouds. In <cit.>, 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 <cit.> that solves PDEs on point clouds where any of the existing methods valid on triangular meshes discussed in Section <ref> 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 <ref>. 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.



§ STATIC HAMILTON-JACOBI-BELLMAN EQUATIONS ON SURFACES





 §.§ Modeling controlled motion on a surface

First, we define the optimal control problem from which the HJB equations can be derived. For each 𝐱∈Γ, let A_𝐱 be set of all unit tangent vectors of Γ at 𝐱, i.e., 
    A_𝐱=S^n-1∩ T_𝐱Γ.
 Note each A_𝐱 is a compact set in ℝ^n for each 𝐱∈Γ. Define A:=⋃_𝐱∈Γ A_𝐱, then A=S^n-1∩ TM(Γ) where TM(Γ) is the tangent bundle of Γ.  The set of admissible controls is given by: 

    𝒜={𝐚:[0,∞) → A}.
 We are interested in the trajectories governed by the dynamical system:

    d𝐲/dt(t)   =f(𝐲(t),𝐚(t))𝐚(t),   t>0 
    𝐲(0)   =𝐱,

where f:Γ× A →ℝ represents the velocity and 𝐚(t) ∈ A_𝐲(t), which ensures that 𝐲(t) ∈Γ for all t≥0. Denote the solution to (<ref>)-(<ref>) (which exists under the assumption (<ref>) below) by 𝐲_𝐱(t,𝐚(t)). For each 𝐱∈Γ define the following subset of the admissible controls:


    𝒜_𝐱:={𝐚(·) ∈𝒜  |  𝐲_𝐱(t,𝐚(t)) ∈Γ t≥0}.
 Then the requirement that 𝐚(t) ∈ A_𝐲(t) for all t ≥ 0 in the dynamical system (<ref>)-(<ref>) is equivalent to 𝐚(·)∈𝒜_𝐱. Refer to 𝐲(·) that satisfy (<ref>)-(<ref>) where 𝐚(·) ∈𝒜_𝐱 as an admissible path on Γ. For the rest of the paper, we will suppress writing the dependence of 𝐲_𝐱(t, 𝐚(t)) on 𝐚(·) when there is no confusion in the context. 
 
 Let r : Γ× A  →ℝ denote a running cost per unit time and g:𝒯→ℝ be an exit time penalty when the state reaches a closed target set 𝒯⊂Γ. The total cost associated with initial state 𝐱 and control 𝐚(·)∈𝒜_𝐱 to reach 𝒯 is given by 

    C(𝐱,𝐚(·))=∫_0^T r(𝐲_𝐱(t),𝐚(t))dt +g(𝐲_𝐱(T)),

where T=min{ t   |  𝐲_𝐱(t) ∈𝒯}.

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

    f  r . 
    
    g min_𝒯 g < ∞.
     R_1,R_2,F_1,  F_2 
     0<R_1≤ r(𝐱,𝐚) ≤ R_2,   𝐱∈Γ𝐚∈ A, 
     0<F_1≤ f(𝐱,𝐚) ≤ F_2,  𝐱∈Γ𝐚∈ A. A1
 The value function u: Γ→ℝ is defined to be the minimal total cost to 𝒯 starting at 𝐱 and is given by 

    u(𝐱)=inf_𝐚(·) ∈𝒜_𝐱 C(𝐱, 𝐚(·)).
 






Under the assumption (<ref>), an optimal control does not necessarily exist. If for the given f and r, we have that the set V(𝐱)={f(𝐱,𝐚)𝐚/r(𝐱, 𝐚)   |  𝐚∈ A_𝐱} is strictly convex for each 𝐱∈Γ, then an optimal control is guaranteed <cit.>. 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 <cit.> for the value function states that for sufficiently small τ >0,

    u(𝐱)=inf_𝐚(·) ∈𝒜_𝐱{∫_0^τ r(𝐲_𝐱(t),𝐚(t))dt+u(𝐲_𝐱(τ)) }.


The corresponding Hamilton-Jacobi-Bellman equation on Γ is

    min_𝐚∈ A_𝐱{ r(𝐱,𝐚)+∇_Γ u(𝐱) ·f(𝐱,𝐚)𝐚}=0,  𝐱∈Γ\𝒯.

The boundary condition is 

    u(𝐱)=g(𝐱),   𝐱∈𝒯.

Note that in (<ref>) we can take the minimum since A_𝐱 is compact.


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

    ||∇_Γ u(𝐱)||=r(𝐱)/f(𝐱).

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



 §.§ Extension to T_ϵ

We derive a Hamiltonian, H, on T_ϵ such that if v: T_ϵ→ℝ is the unique viscosity solution of 
    H(𝐳,∇ v(𝐳))=0
 
with boundary conditions defined appropriately, then 
    v(𝐳)=u(𝐳):=u(P_Γ𝐳),
 
where u is the viscosity solution to (<ref>)-(<ref>).

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


Let 𝐭_1(𝐳) and 𝐭_2(𝐳) be two orthonormal tangent vectors corresponding to the directions that yield the principle curvatures of Γ_η at a point 𝐳. Let 𝐧(𝐳) be the unit normal vector to the tangent plane at 𝐳. Let σ_1(𝐳) and σ_2(𝐳) be the singular values of P_Γ'(𝐳), the Jacobian matrix of P_Γ at 𝐳. 
Then for any real number μ, define the tensor B by

    B(𝐳,μ)=B_0(𝐳)+μ B_1(𝐳),

where 

    B_0(𝐳)   =σ_1^-1(𝐳)𝐭_1(𝐳) ⊗𝐭_1(𝐳) +σ_2^-1(𝐳) 𝐭_2(𝐳) ⊗𝐭_2(𝐳),
    
    B_1(𝐳)   =𝐧(𝐳)⊗𝐧(𝐳).
 
 
In the above definition, B_0 corresponds to a weighted orthogonal projection onto the tangent plane of Γ_η at 𝐳, and B_1 is the projection along the normal of Γ_η passing through 𝐳. 
In <cit.>, it is proven that 
    σ_i=1-d_Γ(𝐳)κ_i(𝐳)
 where κ_1(𝐳) and κ_2(𝐳) are the principal curvatures of the parallel surface Γ_η. 
Note that when 𝐳∈Γ, σ_1=σ_2=1. For the sake of notation, we will suppress the dependence of the tangent vectors, normal vectors, and singular values on the point 𝐳∈ T_ϵ. 

Assuming ϵ small enough, we have that the tangent planes coincide for 𝐳∈ T_ϵ and 𝐱∈Γ if P_Γ𝐳=𝐱. Thus, for the extended optimal control problem, the set of compact control values at each 𝐳∈ T_ϵ is A_P_Γ𝐳. Then A and 𝒜 are as in Section <ref>. Now, consider the following extended dynamical system in T_ϵ: 

    d𝐲/dt(t)   =f(𝐲(t),𝐚(t))B(𝐲(t),μ)𝐚(t),   t>0 
    𝐲(0)   =𝐳, 𝐳∈ T_ϵ.
  Here f: T_ϵ× A →ℝ is given by f(𝐳,𝐚)=f(P_Γ𝐳, 𝐚), and we have that 𝐚(t) ∈ A_𝐲(t) for all t≥0. This last requirement ensures that the trajectory remains in Γ_η for all t ≥ 0 where η=d_Γ(𝐳). The inclusion of the tensor B(𝐲(t),μ) 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 (<ref>)-(<ref>) (which exists under the assumption (<ref>) below) by 𝐲_𝐳(t,𝐚(t)). For each 𝐳∈ T_ϵ with d_Γ(𝐳)=η, we have the following subset of admissible controls:

    𝒜_𝐳:={𝐚∈𝒜  | 𝐲_𝐳(t,𝐚(t)) ∈Γ_η}.

Then the requirement that 𝐚(t) ∈ A_𝐲_𝐳(t) for all t ≥ 0 in the dynamical system (<ref>)-(<ref>) is equivalent to 𝐚(·)∈𝒜_𝐳. Since A_𝐱=A_P_Γ𝐳, we have 𝒜_𝐳=𝒜_P_Γ𝐳.   Refer to 𝐲(·) that satisfy (<ref>)-(<ref>) where 𝐚(·) ∈𝒜_𝐳 as an admissible path on T_ϵ. Again, for the rest of the paper we will write suppress writing the dependence of y_𝐳(t,𝐚(t)) on 𝐚(·). There is no confusion between the notation for admissible paths on Γ and T_ϵ. The initial point will indicate which dynamical system the path solves.

Note that if the initial point 𝐳∈Γ, η=0 and 
    B(𝐲_𝐳(t), μ)𝐚(t)=𝐚(t),
  since 𝐚(t) ∈ T_𝐲_𝐳(t) and σ_1=σ_2=1. Thus, (<ref>)-(<ref>) reduces to the dynamical system on Γ given in (<ref>)-(<ref>), and an admissible path on T_ϵ with 𝐚(·) ∈𝒜_𝐳 is also an admissible path on Γ with the same 𝐚(·) ∈𝒜_𝐳.  

Now, let  𝒯={𝐳∈ T_ϵ  |   P_Γ𝐳∈𝒯}. Let the running cost and exit time penalty on T_ϵ be given by r: T_ϵ× A →ℝ where r(𝐳,𝐚)=r(P_Γ𝐳, 𝐚) and g:T_ϵ→ℝ where g(𝐳)=g(P_Γ𝐳), respectively.  Then the total cost associated with initial state 𝐳 and control 𝐚(·)∈𝒜_𝐳 to reach 𝒯 is given by 

    C(𝐳,𝐚(·))   :=∫_0^Tr( 𝐲_𝐳(t),𝐚(t)) dt +g( 𝐲_𝐳(T))

where T=min{ t  |  𝐲_𝐳(t) ∈𝒯}. Note that if 𝐳∈Γ, then C(𝐳,𝐚(·))=C(𝐳,𝐚(·)). 
The value function v:T_ϵ→ℝ is the minimal total cost to 𝒯 starting at 𝐳 and is given by

    v(𝐳)=inf_𝐚(·) ∈𝒜_𝐳C(𝐳,𝐚(·)).



We have that the assumption <ref> implies:



    (𝐳,𝐚) ↦r( 𝐳,𝐚)   (𝐳,𝐚) ↦f( 𝐳,𝐚)B(𝐳,μ)𝐚
    
    gmin_𝒯g < ∞.
    
     0<R_1≤r(𝐳,𝐚) ≤ R_2,   𝐳∈ T_ϵ𝐚∈ A, 
    
      0<G_1≤g(𝐳) ≤ G_2,  𝐳∈𝒯,
    F_1,F_2 
    0<F_1≤ ||f(𝐳,𝐚)B(𝐳,μ)𝐚|| ≤F_2,  𝐳∈ T_ϵ𝐚∈ A.  B1

Again, while (<ref>) 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(𝐱) is strictly convex, then V(𝐳)={f(𝐳,𝐚)B(𝐳,μ)𝐚/r(𝐳, 𝐚)   |  𝐚∈ A_𝐳} is strictly convex for each z ∈ T_ϵ. Thus, if the optimal control exists in the problem on Γ, it also exists for extended problem on T_ϵ.

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 Γ.  The following theorem states that two admissible paths, 𝐲_𝐳_1(·) and 𝐲_𝐳_2(·), on T_ϵ with equivalent initial points in T_ϵ and equivalent controls stay equivalent for all time.  
 Suppose that 𝐳_1,𝐳_2 ∈ T_ϵ such that  P_Γ𝐳_1=P_Γ𝐳_2. Let  𝐚_1(·) ∈𝒜_𝐳_1 and 𝐚_2(·) ∈𝒜_𝐳_2 such that 𝐚_1≡𝐚_2. Suppose 𝐲_𝐳_1:[0,∞) → T_ϵ   solves (<ref>) with  𝐲_𝐳_1(0)=𝐳_1 and 𝐲_𝐳_2:[0,∞) → T_ϵ solves (<ref>) with  𝐲_𝐳_2(0)=𝐳_2.  Then P_Γ𝐲_𝐳_1(t)=P_Γ𝐲_𝐳_2(t) for all t ≥ 0.



Let 𝐲_1:[0,∞) →Γ and 𝐲_2:[0,∞) →Γ be given by 
    𝐲_1(t)=P_Γ𝐲_𝐳_1(t)
 and 
    𝐲_2(t)=P_Γ𝐲_𝐳_2(t)
 for all t ≥ 0, respectively. We need to show that 𝐲_1(t)=𝐲_2(t) for all t≥ 0. We compute 𝐲_1'. We have the following singular value decomposition for P_Γ ': P_Γ'=U Σ U^T where  
    U=
      [ [ 𝐭_1 𝐭_2   𝐧 ]] Σ= [ [ σ_1   0   0;   0 σ_2   0;   0   0   0 ]].


Since 𝐭_1,  𝐭_2, and 𝐧 are orthonormal, for i,j=1,2,


    (𝐭_i ⊗𝐭_i)(𝐭_j ⊗𝐭_j)=
        𝐭_i ⊗𝐭_i,      i=j,
    
        0,      i≠ j,

and for i=1,2

    (𝐭_i ⊗𝐭_i)(𝐧⊗𝐧) =0.

Using the above, we have  

    𝐲_1'(t)     = (P_Γ𝐲_𝐳_1(t))'
       =P_Γ ' d𝐲_𝐳_1/dt(t)
       = P_Γ'f(𝐲_𝐳_1(t),𝐚_1(t)) B(𝐲_𝐳_1(t),μ)𝐚_1(t)
       =P_Γ'B(𝐲_𝐳_1(t),μ)f(𝐲_𝐳_1(t),𝐚_1(t)) 𝐚_1(t)
       = (𝐭_1 ⊗𝐭_1+𝐭_2 ⊗𝐭_2 )f(𝐲_1(t),𝐚_1(t))𝐚_1(t) 
       = f(𝐲_1(t),𝐚_1(t))(𝐭_1 ⊗𝐭_1+𝐭_2 ⊗𝐭_2 )𝐚_1(t)
       =f(𝐲_1(t),𝐚_1(t))𝐚_1(t).
  
The last equality is true since 𝐭_1 ⊗𝐭_1+𝐭_2 ⊗𝐭_2 is the orthogonal projection of 𝐚_1(t) onto the tangent space at 𝐲_𝐳_1(t) and 𝐚_1(t) ∈ A_𝐲_𝐳_1(t)⊂ T_𝐲_𝐳_1(t)Γ. Thus 𝐲_1 is an admissible path on Γ where 𝐲_1(0)=𝐱=P_Γ𝐳_1=P_Γ𝐳_2 and 𝐚_1(·) ∈𝒜_𝐱=𝒜_P_Γ𝐳_1. By the same reasoning 𝐲_2 also an admissible path with 𝐲_2(0)=𝐱 and 𝐚_2(·) ∈𝒜_𝐱=𝒜_P_Γ𝐳_2.  We assume (<ref>). Thus, the solution to (<ref>)-(<ref>) is unique. Since 𝐚_1≡𝐚_2, 𝐲_1(t)=𝐲_2(t) for all t ≥ 0. 


Now, it easily follows that the value function defined on T_ϵ is the constant normal extension of the value function defined on Γ.


If u:Γ→ℝ is the value function on Γ given by (<ref>), then 

    v(𝐳)=u(𝐳):=u(P_Γ𝐳)

where v is the value function on T_ϵ given by (<ref>).




Let 𝐳∈ T_ϵ and 𝐚(·) ∈𝒜_𝐳. Then 𝐚(·) ∈𝒜_P_Γ𝐳 since 𝒜_𝐳=𝒜_P_Γ𝐳.  Then we have


    C(𝐳,𝐚(·))=∫_0^Tr( 𝐲_𝐳(t),𝐚(t)) dt +g( 𝐲_𝐳(T))

and 

    C(P_Γ𝐳,𝐚(·))=∫_0^Tr( 𝐲_P_Γ𝐳(t),𝐚(t)) dt +g( 𝐲_P_Γ𝐳(T)),
  where 𝐲_𝐳(·) and  𝐲_P_Γ𝐳(·) are admissible paths on T_ϵ with 𝐲_𝐳(0)=𝐳, 𝐲_P_Γ𝐳(0)=P_Γ𝐳, respectively. Recall, that 𝐲_P_Γ𝐳 is also an admissible path on Γ so that C(P_Γ𝐳,𝐚(·))=C(P_Γ𝐳,𝐚(·)).


Theorem <ref> implies that P_Γ𝐲_𝐳(t) = 𝐲_P_Γ𝐳(t) for all t ≥ 0.  
Since 
    𝒯 = {𝐳∈ T_ϵ  |  P_Γ𝐳∈𝒯},

we have 

    min{t   |   P_Γ𝐲_𝐳(t)= 𝐲_P_Γ𝐳(t) ∈𝒯}    = min{ t   |  𝐲_𝐳(t) ∈𝒯}.
 

Therefore, T=T, and

    C(𝐳,𝐚(·))
       =C(P_Γ𝐳, 𝐚(·))
       =C(P_Γ𝐳,𝐚(·)).
 
Finally, we have

    v(𝐳)   =inf_𝐚(·) ∈𝒜_𝐳C(𝐳,𝐚(·))
       =inf_𝐚(·) ∈𝒜_P_Γ𝐳C(P_Γ𝐳,𝐚(·))
       =u(P_Γ𝐳).




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_ϵ associated the value function v. Again, the dynamic programming principle states that for sufficiently small τ>0, we have  

    v(𝐳)=inf_𝐚(·) ∈𝒜_𝐳{∫_0^τr( 𝐲_𝐳(t),𝐚(t))dt +v(𝐲_𝐳(τ))}.

The corresponding Hamilton-Jacobi-Bellman equation on T_ϵ is

    min_𝐚∈A_𝐳{r( 𝐳,𝐚)+∇ v(𝐳) ·f( 𝐳,𝐚)B(𝐳,μ)𝐚}=0,   𝐳∈ T_ϵ\𝒯 .

The boundary condition is 

    v(𝐳)=g( 𝐳),   𝐳∈𝒯.
 

The assumption (<ref>) implies that the value function coincides with the unique viscosity solution <cit.>. Corollary <ref> implies that the viscosity solution of (<ref>)-(<ref>) is given by 
    v(𝐳)=u(P_Γ𝐳),
 where u is the viscosity solution of (<ref>)-(<ref>) on T_ϵ. An interesting result is that we can prove (<ref>) without Corollary <ref>, 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:

Assume ϵ is small enough so that P_Γ is differentiable, and let |η| <ϵ. Given ϕ_0 ∈ C^1(Γ), define ϕ_η∈ C^1(Γ_η)  by ϕ_η(𝐳):=ϕ_0(P_Γ𝐳). Then for 𝐳∈Γ_η, 

    ∇_Γϕ_0(P_Γ𝐳)=B(𝐳,μ) ∇_Γ_ηϕ_η (𝐳).



Let ϕ_0: T_ϵ→ℝ be the constant normal extension of ϕ_0. Then ϕ_0 ∈ C^1(T_ϵ) and ϕ_0Γ_η=ϕ_η. From Theorem A.1 in <cit.>,  we have  that ∇_Γϕ_0(P_Γ𝐳)=B(𝐳,μ)∇ϕ_0 (𝐳).  Since ϕ_0 is the constant normal extension of ϕ_0, ∇ϕ_0(𝐳) ∈ T_𝐳Γ_η. Thus, if 𝐳∈Γ_η, ∇ϕ_0(𝐳)=∇_Γ_ηϕ_η(𝐳), and (<ref>) holds.

The following theorem proves (<ref>) directly using the definition of viscosity solutions.
 If u:Γ→ℝ is the viscosity solution to (<ref>)-(<ref>), then the constant normal extension of u, u: T_ϵ→ℝ, is the viscosity solution to (<ref>)-(<ref>).

 
We have  u(𝐱)=g(𝐱) for 𝐱∈𝒯 and P_Γ𝐳∈𝒯 if 𝐳∈𝒯. Then for 𝐳∈𝒯, 

    u(𝐳)   =u(P_Γ𝐳)=g(P_Γ𝐳)=g(𝐳).

Thus, u satisfies the boundary conditions (<ref>). Next, we will show that u is a viscosity subsolution of (<ref>).

Assume 𝐳_0 ∈ T_ϵ and ϕ∈ C^1(T_ϵ) such that u -ϕ has a strict local maximum at 𝐳_0 and  (u-ϕ)(𝐳_0)=0. We need to show that 

    min_𝐚∈A_𝐳{r( 𝐳_0,𝐚)+∇ϕ(𝐳_0) ·f( 𝐳_0,𝐚)B(𝐳_0,μ)𝐚}≤0.


Let η=d_Γ(𝐳_0) and ϕ_0: Γ→ℝ be the restriction to Γ of the normal extension of ϕΓ_η. Then we have that ϕ_0 ∈ C^1(Γ) and ϕ_η:= ϕΓ_η∈ C^1(Γ_η) such that ϕ_η(𝐳)=ϕ_0(P_Γ𝐳) for 𝐳∈Γ_η.

First, we will prove 

    B(𝐳_0,μ)∇ϕ(𝐳_0) ·𝐚 =B(𝐳_0,μ)∇_Γ_ηϕ_η(𝐳_0) ·𝐚.

 Since 𝐚∈ A_𝐳_0⊂ T_𝐳_0Γ, we have 𝐚·𝐧=0 and 
    B(𝐳_0,μ)𝐧⊗𝐧∇ϕ(𝐳_0)·𝐚=0.
 Now, ∇_Γ_ηϕ_η(𝐳_0)=(I-𝐧⊗𝐧)∇ϕ(𝐳_0). Therefore, (<ref>) holds.

Lemma <ref>  implies

    B(𝐳_0,μ)∇_Γ_ηϕ_η(𝐳_0)=∇_Γϕ_0(P_Γ𝐳_0).

Let 𝐱_0=P_Γ𝐳_0. Then using the fact that B(𝐳_0,μ) is symmetric, (<ref>), and (<ref>), we have 

    min_𝐚∈A_𝐳_0{r( 𝐳_0,𝐚)+∇ϕ(𝐳_0) ·   f( 𝐳_0,𝐚)B(𝐳_0,μ)𝐚}
       =min_𝐚∈A_𝐳_0{r( 𝐳_0,𝐚)+B(𝐳_0,μ)∇ϕ(𝐳_0) ·f( 𝐳_0,𝐚)𝐚}
       =min_𝐚∈A_𝐳_0{r( 𝐳_0,𝐚)+B(𝐳_0,μ)∇_Γ_ηϕ_η(𝐳_0) ·f( 𝐳_0,𝐚)𝐚}
       =min_𝐚∈A_𝐳_0{r( 𝐳_0,𝐚)+∇_Γϕ_0(P_Γ𝐳_0) ·f( 𝐳_0,𝐚)𝐚}
       =min_𝐚∈A_𝐱_0{r( 𝐱_0,𝐚)+∇_Γϕ_0(𝐱_0) ·f( 𝐱_0,𝐚)𝐚}.


Now, 𝐱_0 ∈Γ, ϕ_0 ∈ C^1(Γ) such that u-ϕ_0 has a strict local maximum at 𝐱_0 and 
    (u-ϕ_0)(𝐱_0)=(u-ϕ)(𝐳_0)=0.
 Since u is a viscosity subsolution of (<ref>), 
    min_𝐚∈A_𝐱_0{r( 𝐱_0,𝐚)+∇_Γϕ_0(𝐱_0) ·f( 𝐱_0,𝐚)𝐚}≤0.
 Thus, (<ref>) holds. The same argument can be applied to show that u is a viscosity supersolution of (<ref>).


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 Γ 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 Γ, i.e., the extended set of controls is A=S^n-1. In this setup, the value function is the constant normal extension of the value function on Γ. Before we proceeding, we present a counter example to show that u is not necessarily the value function of the extended control problem if the extended speed function f is anisotropic and A=S^n-1.



  §.§.§ Example

Consider the following control problems on Γ and T_ϵ: Let 
    Γ=(-1,1) ⊂ℝ
 and 
    𝒯={1}.
 Let ϵ>√(3)/3, then 
    T_ϵ=(-1,1)×(-ϵ,ϵ)
 and 
    𝒯={1}×(-ϵ,ϵ).
 For the control problem on Γ, the set of controls is A={(-1,0),(1,0)} and 𝒜 is the usual set of admissible controls. Define f:Γ× A →ℝ by f(𝐱,𝐚)=√(3) and suppose that we have unit running cost, i.e., r(𝐱,𝐚)=1. Then it is clear that u(0)=√(3)/3 which is the minimum travel time from 𝐱=0 to 𝒯 while traveling at speed √(3).

Now for the extended control problem, let the control values be the set A=S^1. Then the set of admissible controls is given by 
    𝒜:={a:[0,∞) →A}.
 Define f: T_ϵ×𝒜→ℝ such that f(𝐳,𝐚)=f(𝐚) and the set {f(𝐚) 𝐚  |  𝐚∈A} is the ellipse given by x^2+y^2/9=1 rotated at the origin clockwise by the angle π/3. Then f(𝐚)=√(3) if 𝐚∈ A so that f(𝐳,𝐚)=f(P_Γ𝐳, 𝐚) for 𝐚∈ A.  Again, assume unit running cost.

Let v be the value function for this extended control problem. If v=u, then we should have v((0,0))=u(0). However, note that f(𝐳,(√(3)/2,1/2))=3 and the Euclidean distance between (0,0) and (1,√(3)/3) is 2√(3)/3. Thus, 
    C((0,0), a(·))=2√(3)/9
 when a(·) ≡ (√(3)/2,1/2). By the definition of the value function we must have 
    v((0,0))≤2√(3)/9 < √(3)/3= u(0).
 Therefore, u cannot be the value function for the extended control problem.



§ 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 <ref>, we will show that value function on T_ϵ is the normal extension of the value function on Γ even when we do not restrict the controls to the tangent space in the extended optimal control problem.


 §.§ The controlled dynamics on Γ for isotropic speed and cost functions

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

    ||∇_Γ u(𝐱)||=r(𝐱)/f(𝐱).

We have the HJB equation

    H(𝐱, ∇_Γ u (𝐱))=r(𝐱)-f(𝐱)||∇_Γ u(𝐱)||=0,   𝐱∈Γ,

and the boundary condition is 

    u(𝐱)=g(𝐱),   𝐱∈𝒯.

Assuming (<ref>), the value function u coincides with the unique viscosity solution of (<ref>)-(<ref>). We also have that V(𝐱) is strictly convex for all 𝐱∈Γ. Thus, an optimal control is guaranteed. 



 §.§ Extension to T_ϵ in the isotropic case

We derive the extended Hamiltonian, H, on T_ϵ from the extended isotropic optimal control problem such that if v: T_ϵ→ℝ is the unique viscosity solution of 
    H(𝐱,∇ v(𝐳))=0
 with boundary conditions defined appropriately, then 
    v(𝐳)=u(P_Γ𝐳),
 
where u is the solution to (<ref>)-(<ref>).

The notations 𝐭_1, 𝐭_2, 𝐧, σ_1, and σ_2 are as in Section <ref>. The extended set of compact control values is A=S^n-1, and the set of admissible controls is given by: 
    𝒜={𝐚: [0,∞) →A}.
 This is where the formulation differs from the previous section. We do not restrict the controls to belong to the tangent bundle of Γ. We extend the speed, cost, and exit time penalty functions  as above: f(𝐳)=f(P_Γ𝐳), r(𝐳)=r(P_Γ𝐳), and g(𝐳)=g(P_Γ𝐳). The dynamical system is given by 

    d𝐲/dt(t)   =f(𝐲(t))B(𝐲(t),μ)𝐚(t),   t>0 
    𝐲(0)   =𝐳, 𝐳∈ T_ϵ,
 
where 𝐚(·) ∈𝒜. Admissible paths on T_ϵ with extended controls are solutions to (<ref>)-(<ref>) denoted by y_𝐳(t).  Note that the admissible paths are not restricted to the parallel surface, Γ_η, to which the initial point 𝐳 belongs.

The total cost function associated to the initial state 𝐳∈ T_ϵ and control 𝐚(·) ∈𝒜 is the same as in Section <ref>.
The value function v:T_ϵ→ℝ is the minimal total cost to 𝒯 starting at 𝐳 and is given by 
    v(𝐳)=inf_𝐚(·) ∈𝒜C(𝐳,𝐚(·)).


We have that (<ref>) implies the following is true:


    (𝐳,𝐚) ↦f( 𝐳)B(𝐳,μ)𝐚
    r
    gmin_𝒯g < ∞.
    
     0<R_1≤r(𝐳) ≤ R_2,   𝐳∈ T_ϵ. 
    F_1,F_2 
    0<F_1≤ ||f(𝐳)B(𝐳,μ)𝐚|| ≤F_2,  𝐳∈ T_ϵ𝐚∈A.  C1

We have that (<ref>) implies the value function, v, coincides with the unique viscosity solution of (<ref>)-(<ref>) <cit.>. We also have that V(𝐳)={f(𝐳)𝐚/r(𝐳) | 𝐚∈A} is strictly convex for all 𝐳∈Γ. Thus, an optimal control is guaranteed. 

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

Given 𝐚(·) ∈𝒜 and 𝐳_1,𝐳_2 ∈ T_ϵ such that P_Γ𝐳_1=P_Γ𝐳_2,  let 
    𝐲_𝐳_1:[0,∞) → T_ϵ
 and 
    𝐲_𝐳_2:[0,∞) → T_ϵ
   solve (<ref>)  with  𝐲_𝐳_1(0)=𝐳_1 and  𝐲_𝐳_2(0)=𝐳_2, respectively.  Then P_Γ𝐲_𝐳_1(t)=P_Γ𝐲_𝐳_2(t) for all t ≥ 0.



The proof of Theorem <ref> is analogous to the proof of Theorem <ref> except the control 𝐚(·) belongs to the extended control space, 𝒜. A consequence of the control belonging to the extended space is that now the projected paths P_Γ𝐲_𝐳(·) may not solve (<ref>). This is due to the fact that (𝐭_1 ⊗𝐭_1+𝐭_2 ⊗𝐭_2 )𝐚(t)≠𝐚(t) if 𝐚(t) ∉T_P_Γ𝐲_𝐳(t) where 𝐭_1 and 𝐭_2 are the basis tangent vectors of T_P_Γ𝐲_𝐳(t).  We now prove that the extended value function in the isotropic case is constant along the normals of Γ.


If P_Γ𝐳_1=P_Γ𝐳_2 for 𝐳_1,𝐳_2 ∈ T_ϵ, then v(𝐳_1)=v(𝐳_2).


Let 𝐚(·)∈𝒜. We have  

    C(𝐳_1,𝐚(·))=∫_0^T_1r( 𝐲_𝐳_1(t)) dt +g( 𝐲_𝐳_1(T_1))

and 

    C(𝐳_2,𝐚(·))=∫_0^T_2r( 𝐲_𝐳_2(t)) dt +g( 𝐲_𝐳_2(T_2)),
  where 𝐲_𝐳_1(·) and  𝐲_𝐳_2(·) are admissible paths on T_ϵ with extended control, 𝐚(·)∈𝒜, and 𝐲_𝐳_1(0)=𝐳_1, 𝐲_𝐳_2(0)=𝐳_2 where P_Γ𝐳_1=P_Γ𝐳_2.


Theorem <ref> implies that P_Γ𝐲_𝐳_1(t) = P_Γ𝐲_𝐳_2(t) for all t ≥ 0.  
Since 
    𝒯 = {𝐳∈ T_ϵ  |  P_Γ𝐳∈𝒯},

we have 

    min{t   |   P_Γ𝐲_𝐳_1(t)=P_Γ𝐲_𝐳_2(t) ∈𝒯}    = min{ t   |  𝐲_𝐳_1(t) ∈𝒯} =  min{ t   |  𝐲_𝐳_2(t) ∈𝒯}.
 
Therefore, T_1=T_2 and  C(𝐳_1, 𝐚(·))=C(𝐳_2, 𝐚(·)) for 𝐚(·) ∈𝒜. Now, 

    v(𝐳_1)   =inf_𝐚(·) ∈𝒜C(𝐳_1,𝐚(·)) =inf_𝐚(·) ∈𝒜C(𝐳_2,𝐚(·)) =v(𝐳_2).




In the above proof, it is not immediate that we can take the infimum over 𝒜⊂𝒜, which would imply that v is the normal extension of u, the value function on Γ. We will use the corresponding HJB equations and Corollary <ref> to show that v=u.

We define the HJB equation on T_ϵ associated the value function v. Again, the dynamic programming principle states that for sufficiently small τ>0 we have  

    v(𝐳)=inf_𝐚(·) ∈𝒜{∫_0^τr( 𝐲_𝐳(t))dt +v(𝐲_𝐳(τ))}.

We get the Hamilton-Jacobi-Bellman equation on T_ϵ:

    min_𝐚∈A{r( 𝐳)+∇ v(𝐳) · B(𝐳,μ)f( 𝐳)𝐚}=0.

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

    H(𝐳,∇ v(𝐳))=r( 𝐳)-f( 𝐳)||B(𝐳,μ) ∇ v(𝐳)||=0.

The boundary condition is 

    v(𝐳)=g( 𝐳),   𝐳∈𝒯.



Since the value function coincides with the unique viscosity solution, Corollary <ref> implies that the viscosity solution of (<ref>)-(<ref>) is a constant along the normals of Γ, i.e., there is a function w:Γ→ℝ such that for 𝐳∈ T_ϵ, 
    v(𝐳)=w(𝐳):=w(P_Γ𝐳).
 We now prove that w=u where u is the viscosity solution of (<ref>)-(<ref>). 


 
If u:Γ→ℝ is the viscosity solution to (<ref>)-(<ref>), then the constant normal extension of u, u: T_ϵ→ℝ, is the viscosity solution to (<ref>)-(<ref>).

 
Just as in Theorem <ref>, u satisfies the boundary conditions (<ref>).  Assume for contradiction that u is not the viscosity solution to (<ref>). But from Corollary <ref> we have that viscosity solution of (<ref>)-(<ref>), v=w for some w ∈ C^1(Γ). Therefore, if u is not the unique viscosity solution then we have u ≠ w. For contradiction, we will show that w is the viscosity solution of (<ref>) and thus, u=w.

Let 𝐱_0 ∈Γ and ϕ∈ C^1(Γ) such that w-ϕ has a local maximum at 𝐱_0. Consider the normal extensions of w and ϕ, then w-ϕ has a local max at 𝐱_0. Since w is a viscosity subsolution to (<ref>) and ϕ∈ C^1(T_ϵ), we have that 

    r( 𝐱_0)-f( 𝐱_0)||B(𝐱_0,μ) ∇ϕ(𝐱_0)|| ≤ 0.
 
Since 𝐱_0 ∈Γ, σ_1=σ_2=1 and ∇ϕ(𝐱_0) ∈ T_𝐱_0Γ. Therefore, 
    B(𝐱_0,μ)=𝐭_1⊗𝐭_1+𝐭_2⊗𝐭_2+μ𝐧⊗𝐧,
 and we have that B(𝐱_0,μ)∇ϕ(𝐱_0)=∇ϕ(𝐱_0)=∇_Γϕ(𝐱_0).  Thus, 

    r( 𝐱_0)-f( 𝐱_0)||∇_Γϕ(𝐱_0)|| ≤ 0,
 
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 (<ref>) 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_ϵ have the “choice" to leave the tangent space of Γ, the optimal paths in the extended control problem remain on or parallel to Γ.




§ NUMERICAL IMPLEMENTATION AND SIMULATIONS




 §.§ Numerical setup

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

We use the Lax-Friedrichs fast sweeping methods in <cit.> and a high order version in <cit.> for our simulations. It is also possible to an upwind scheme generalized from <cit.>. 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.[ 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 <ref>, the equivalent equation on T_ϵ is 
    ||B(𝐳,μ) ∇ v(𝐳)||   =1,   𝐳∈ T_ϵ\𝒯
    
    v(𝐳)   =g(𝐳),   𝐳∈𝒯.

Again, 
    B(𝐳,μ)=σ_1^-1𝐭_1 ⊗𝐭_1 +σ_2^-1𝐭_2 ⊗𝐭_2+μ𝐧⊗𝐧.
 Since the desired solution, v, has been proven to be constant along the normals of the surface, Γ, 𝐧⊗𝐧∇ v(𝐳)=0. Therefore, μ can be any real number. We let μ=1 in all of our computations.

Approximating the solution to (<ref>)-(<ref>) requires the computation of the singular values and vectors of the derivative of the closest point mapping, P_Γ'. We defer the discussion of these approximations to Section <ref>. In the last example, we apply our framework to solve an HJB equation with an anisotropic speed function. 




  §.§.§ Boundary closure


We note that the analytical formulation of the HJB equations does not require boundary conditions on ∂ T_ϵ. However, since we are using Cartesian grids we must take careful consideration of the discretization near the boundary, ∂ T_ϵ^h. When approximating the partial derivatives of v, a neighboring grid node may lie outside of T_ϵ^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:

    
  * Project the node into the narrow band.
    
  * 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_ϵ^h. 


Next, we describe how the boundary closure procedure affects the numerical accuracy. Let 𝐳_i be a ghost node, and 𝐳^α_i=𝐳_i-α𝐧_i be the projection of the ghost node into T_ϵ. Denote the solution at 𝐳_i^α by v_i^α:=v(𝐳_i^α). Suppose that we use a first order scheme to discretize the HJB equation on T_ϵ^h. Then the interpolation used in the approximation of v_i^α 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 𝐳_i. 



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

When choosing α 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^α must be in T_ϵ^h. In our numerical simulations, we use cubic interpolation to approximate v( z^α). Thus, it requires 4^3 nearby grid nodes for three dimensional cases.  It can be shown that if |d_Γ(𝐳_i)-α|≤ϵ-2√(n)h, where n is the dimension. Then the inner nodes used to interpolate the value at each ghost node are in T_ϵ^h. We suggest to use a projected node as close to the boundary node as possible. Therefore in our numerical simulations, we choose α=d_Γ(𝐳_i)-(ϵ-2√(n)h). We discuss the effect of the choice of α in Section <ref>.

The maximal value of ϵ is restricted by the curvatures of Γ. We must have ϵ smaller than the reach of the surface, which depends on the reciprocal of the maximal curvatures.

The minimal value of ϵ 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 ϵ>2√(n)h.


 The resulting numerical discretization of the HBJ equation on T_ϵ^h is convergent 
 for any ϵ in the interval described above; i.e. it is convergent for very thin 𝒪(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 ϵ∼𝒪(h^p), 0<p<1, in the related work of <cit.>, dealing with singular integrals.
 




  §.§.§ Approximation errors

The error computed by the proposed algorithm can be written as the sum of errors corresponding to different approximations:
 
    Error = E_model(ϵ)+E_Δ(h)+E_Γ(h)+E_L_p(h,ϵ)   ,

 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≡ 0 for any ϵ that is smaller than the reach of the surface;
     
  * E_Δ corresponds to the numerical error for approximating the extended PDE problem. Our extension allows for higher order (in h) methods to be applied;
     
  * E_Γ 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. 
 
 
 


  §.§.§ 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 <ref>.) Denote the point cloud by Γ_K⊂Γ, 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_Γ using the following strategy:

    
  * For each 𝐳_i ∈ T_ϵ^h, estimate Γ_K locally on the grid nodes. Call the local surface associated to the point 𝐳_i, Γ_K(𝐳_i).
    
  * Compute the closest point of 𝐳_i on the local surface, Γ_K(𝐳_i).


 In our implementation, we use biquadratic interpolation at each 𝐳_i to compute Γ_N(𝐳_i). Another option to locally estimate the surface is least squares as in <cit.>. We then use Newton's method to obtain the closest point of each 𝐳_i on  Γ_N(𝐳_i). Denote the closest point mappings by P_Γ_N^int(𝐳_i). Then on each grid node 𝐳_i, the equation is discretized as usual, with P_Γ^'(𝐳_i) being approximated by finite differences of P^int_Γ_N (𝐳_i). We use a fourth-order finite difference scheme to estimate P_Γ' and singular value decomposition to obtain the singular values and singular vectors of P_Γ'. We show the effect of this approximation procedure in Section <ref>.


In any case, the term, E_Γ, in (<ref>) 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 <ref>, 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 <cit.>.  We will investigate this important issue in a future manuscript. Nevertheless, in Section <ref>, we present some results revealing the effect of noise in the point cloud.




 §.§ Example 1.1

First, we present a numerical convergence study. We approximate the solution to (<ref>)-(<ref>), where Γ is a sphere centered at (0,0,0), with radius, r_0=0.5, and 𝒯 is a point on the sphere. The exact singular values and vectors in B(𝐳,1
) are used. In the case of a sphere 
    σ_1=σ_2=1-r(𝐳)-r_0/r(𝐳)=r_0/r(𝐳),
 where r(𝐳) is the radius of the sphere going though 𝐳 centered at (0,0,0). Figure <ref> shows the distance function to a point on the sphere for two different view points.

Since the solution is constant along the normals of Γ, we can easily estimate the L_1 error. Using the formulation in <cit.>, the estimation of the L_1 error is given by the following formula:

    ||u-u^h||_L_1(Γ)≈∑_𝐳_i ∈ T_ϵ^h (u(P_Γ𝐳_i)-v_i^h)K_ϵ(d_Γ(𝐳_i))J(𝐳_i)h^3,

where u is the exact solution and u^h is the approximated solution on Γ using the approximate solution, v_i^h, of the extended equation with grid size h.
Here, 
    K_ϵ(d_Γ(𝐳_i))=A *exp (2/d_Γ(𝐳_i)^2-1)
 where A =7.513931532835806 
and 

    J(𝐳_i)=σ_1σ_2,
 where σ_1,σ_2 are given in (<ref>).

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 𝒯 with the exact solution if u(𝐳) <0.2 The L_1 and L_∞ errors are reported in Table <ref>. 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_∞ 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_∞ 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_ϵ 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 <ref>, we present the L_1 and L_∞ errors and orders for the third order method given in <cit.>. We initialize the solution near 𝒯 to the exact solution if u(𝐳) < 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_∞ error is first order.









 §.§ 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.[The 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 <ref> 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 𝒯 as in Section <ref>.   By comparing the L_1 error from Table <ref> for N=101 and the L_1 error in Table <ref> 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 <ref> for N=301 and the L_1 error from Table <ref> 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 <ref>, we also consider when noise is introduced to the point clouds of different sampling densities. The noise is applied by adding δ*(-1+2*rand())/2 to each point in the point clouds, where δ = 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.






 §.§ 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 ϵ=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, 𝐳_i, is outside of T_ϵ^h and is projected into T_ϵ along the normal of the surface at that point, i.e., 

    𝐳_i ↦𝐳_i-α𝐧_i.
 Table <ref> shows the L_1 error for the depths: 
    α=d(𝐳)-(10-2√(3))h, d(𝐳)-3h, d(𝐳), d(𝐳)+3h.
 We can see that the error the solution is not very sensitive to the choice of depth. In our simulations, we chose α=ϵ -2√(3)h.




 §.§ 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 <ref>; (2) our method with the point cloud procedure described in <ref> and finite differences to estimate P_Γ'; (3) the method in <cit.>. Since method (3) requires ϵ=Ch^γ where C=√(3) and  γ∈ (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 <cit.>. 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 Γ. Therefore, we can only report the L_∞ error, and we use trilinear interpolation to approximate the solutions on Γ for (3). The errors are reported in Table <ref>. We can also see that estimating P_Γ from a point cloud 
and P_Γ' from finite differences does not greatly influence the L_1 error of the method. 







 §.§ Example 2












In Figures <ref>, <ref>, and <ref>, 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 ϵ=4h and point cloud representations for the surfaces. The number of points in the point clouds of the torus[The point cloud for the torus was generated using the standard parametrization of a torus.]  and bunny[The point cloud for the Stanford bunny is generated from a refinement of the triangulated Stanford bunny from https://casual-effects.com/data/ <cit.>.]  are 178,350 and 228,096, respectively. For the elephant, mask, human skull, and dinosaur skull[The point clouds for the surfaces in Figure <ref> 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 <ref>, 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 <ref> 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 Γ.  Since the solution to (<ref>)-(<ref>) is constant along the normals of Γ, the gradient always belongs to the tangent spaces of Γ or its parallel surfaces, Γ_η. We have that B(𝐱,μ) 𝐚=𝐚 for 𝐱∈Γ and 𝐚∈ A_𝐱=T_𝐱Γ∩ S^n-1 since σ_1=σ_2=1. Thus, if the initial point is 𝐳∈Γ, the geodesic on the surface can be obtained by solving the dynamical system 
    d𝐲/dt(t)   =-∇ v(𝐲(t))/||∇ v(𝐲(t))||,   t>0 
    𝐲(0)   =𝐳, 𝐳∈Γ .



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



 §.§ Example 3




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


    min_𝐚∈ A_𝐳{∇ v(𝐳) ·f( 𝐳,𝐚)B(𝐳,μ)𝐚 +1 }   =0 ,   𝐳∈ T_ϵ\𝒯,
     
     v(𝐳)   =g(𝐳),   𝐳∈𝒯,

where f: Γ× A →ℝ is a curvature based speed function on the surface. We use the speed function proposed in <cit.>. 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 Γ at 𝐱 in the direction 𝐚 is given by 
    κ_𝐚(𝐱)=𝐚^T[ [ κ_1   0;   0 κ_2;     ]]𝐚

  where κ_1 and κ_2 are the principal curvatures of Γ at 𝐱.
Then the curvature-minimizing  speed function is given by
   
    f(𝐱,𝐚)=(-b ||κ_𝐚(𝐱)||),

where b is a positive constant. When b=0, (<ref>) 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 <ref>. 

In this example, the set V(𝐱)={f(𝐱,𝐚) 𝐚  |  𝐚∈ A_𝐳}  is not neccesarily convex for all 𝐱 in Γ. 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(𝐱, 𝐚)𝐚=𝐚 for 𝐱∈Γ. The paths can easily be extracted because we compute the minimizing control at each grid node when solving (<ref>)-(<ref>). Once we have the optimal (or suboptimal) control values, 𝐚^*(·), we then solve dynamical system 
    d𝐲/dt(t)   =f(𝐲(t),𝐚^*(t))𝐚^*(t),   t>0 
    𝐲(0)   =𝐳, 𝐳∈Γ.
 We plot three anisotropic geodesics for each b value in Figure <ref>. 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.





§ 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 <cit.>, 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.




abbrv






