# A semi-Lagrangian scheme for Hamilton-Jacobi-Bellman equations with oblique boundary conditions

We investigate in this work a fully-discrete semi-Lagrangian approximation of second order possibly degenerate Hamilton-Jacobi-Bellman (HJB) equations on a bounded domain with oblique boundary conditions. These equations appear naturally in the study of optimal control of diffusion processes with oblique reflection at the boundary of the domain. The proposed scheme is shown to satisfy a consistency type property, it is monotone and stable. Our main result is the convergence of the numerical solution towards the unique viscosity solution of the HJB equation. The convergence result holds under the same asymptotic relation between the time and space discretization steps as in the classical setting for semi-Lagrangian schemes. We present some numerical results that confirm the numerical convergence of the scheme.

## Authors

• 2 publications
• 2 publications
• 2 publications
• 1 publication
06/10/2017

### A fully semi-Lagrangian discretization for the 2D Navier--Stokes equations in the vorticity--streamfunction formulation

A numerical method for the two-dimensional, incompressible Navier--Stoke...
02/16/2021

### A Hybrid Semi-Lagrangian Cut Cell Method for Advection-Diffusion Problems with Robin Boundary Conditions in Moving Domains

We present a new discretization for advection-diffusion problems with Ro...
06/15/2021

### Convergence of a Lagrangian-Eulerian scheme via the weak asymptotic method

This work presents a suitable mathematical analysis to understand the pr...
11/11/2019

### The waiting time phenomenon in spatially discretized porous medium and thin film equations

Various degenerate diffusion equations exhibit a waiting time phenomenon...
03/23/2022

### A convergent SAV scheme for Cahn–Hilliard equations with dynamic boundary conditions

The Cahn-Hilliard equation is one of the most common models to describe ...
04/06/2020

### Stable Boundary Conditions and Discretization for PN Equations

A solution to the linear Boltzmann equation satisfies an energy bound, w...
09/20/2021

### How to train your solver: A method of manufactured solutions for weakly-compressible SPH

The Weakly-Compressible Smoothed Particle Hydrodynamics (WCSPH) method i...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

In this work we deal with the numerical approximation of the following parabolic Hamilton-Jacobi-Bellman (HJB) equation

 (1.1)

In the system above, , is a nonempty smooth bounded open set and and are nonlinear functions having the form

 (1.2) H(t,x,p,M) = supa∈A{−12{\rm Tr}(σ(t,x,a)σ(t,x,a)⊤M)−⟨μ(t,x,a),p⟩−f(t,x,a)}, (1.3) L(t,x,p) = supb∈B{⟨γ(x,b),p⟩−g(t,x,b)},

where and are nonempty compact sets, , with , , , , with being an open set containing , , and .

If and , for some and , and , with

being the unit outward normal vector to

at , then (1.1) reduces to a standard linear parabolic equation with Neumann boundary conditions. In the general case, and after a simple change of the time variable in order to write (1.1) in backward form, the HJB equation (1.1) appears in the study of optimal control of diffusion processes with controlled reflection on the boundary (see e.g. [27] for the first order case, i.e. , and [26, 11] for the general case). Since the HJB equation (1.1) is possibly degenerate parabolic, one cannot expect the existence of classical solutions and we have to rely on the notion of viscosity solution (see e.g. [16]). Moreover, as it has been noticed in [25, 27], in general the boundary condition in (1.1) does not hold in the pointwise sense and we have to consider a suitable weak formulation of it. We refer the reader to [27, 4] and [16, 6, 7, 24, 12], respectively, for well-posedness results for HJB equations with oblique boundary condition in the first and second order cases.

The study of the numerical approximation of solutions to HJB and, more generally, fully nonlinear second order Partial Differential Equations (PDEs), has made important progress over the last few decades. Most of the related literature consider the case where

, or where a Dirichlet boundary condition is imposed on the boundary . We refer the reader to [19, 20, 30] and the references therein for the state of the art on this topic. By contrast, the numerical approximation of solutions to (1.1) has been much less explored. Indeed, to the best of our knowledge only the methods in [31, 1] can be applied to approximate (1.1) in the particular first order case (). Moreover, in [31], where a finite difference scheme is proposed, the function defining the boundary condition has the particular form . On the other hand, both references consider Hamiltonians which are not necessarily convex with respect to . Let us also mention the reference [2], where, in the context of mean curvature motion with nonlinear Neumann boundary conditions, the authors propose a discretization that combines a Semi-Lagrangian (SL) scheme in the main part of the domain with a finite difference scheme near the boundary.

The main purpose of this article is to provide a consistent, stable, monotone and convergent SL scheme to approximate the unique viscosity solution to (1.1). By the results in [6], the latter is well-posed in under the assumptions in Sect. 2 below. Semi-Lagrangian schemes to approximate the solution to (1.1) when (see e.g. [13, 17]) can be derived from the optimal control interpretation of (1.1) and a suitable discretization of the underlying controlled trajectories. These schemes enjoy the feature that they are explicit and stable under an inverse Courant-Friedrichs-Lewy (CFL) condition and, consequentely, they allow large time steps. A second important feature is that they permit a simple treatement of the possibly degenerate second order term in . The scheme that we propose for preserves these two properties and seems to be the first convergent scheme to approximate (1.1) with the rather general asumptions in Sect. 2. In particular, our results cover the stochastic and degenerate case. Consequently, from the stochastic control point of view, our scheme allows to approximate the so-called value function of the optimal control of a controlled diffusion process with possibly oblique reflection on the boundary (see [11]). The main difficulty in devising such a scheme is to be able to obtain a consistency type property at points in the space grid which are near the boundary while maintaining the stability. This is achieved by considering a discretization of the underlying controlled diffusion which suitably emulates its reflection at the boundary in the continuous case. We refer the reader to [28] for a related construction of a semi-discrete in time approximation of a second order non-degenerate linear parabolic equation.

The remainder of this paper is structured as follows. In Sect. 2 we state our assumptions, we recall the notion of viscosity solution to (1.1) and the well-posedness result. In Sect. 3 we provide the SL scheme as well as its probabilistic interpretation (in the spirit of [28]). The latter will play an important role in Sect. 4, which is devoted to show a consistency type property and the stability of the SL scheme. By using the half-relaxed limits technique introduced in [5], we show in Sect. 5 our main result, which is the convergence of solutions to the SL scheme towards the unique viscosity solution to (1.1). The convergence is uniform in and holds under the same asymptotic condition between the space and time steps than in the case . Next, in Sect. 6 we first illustrate the numerical convergence of the SL scheme in the case of a one-dimensional linear equation with homogeneous Neumann boundary conditions. In this case the numerical results confirm that the boundary condition in (1.1) is not satisfied at every , but it is satisfied in the viscosity sense recalled in Sect. 2 below. In a second example, we consider a two dimensional degenerate second order nonlinear equation on a circular domain with non-homogeneous Neumann and oblique boundary conditions. In the last example, we consider a two-dimensional non-degenerate nonlinear equation on a non-smooth domain. Due to the lack of regularity of , our convergence result does not apply. However, the SL scheme can be successfully applied, which suggests that our theoretical findings could hold for more general domains. This extension as well as the corresponding study in the stationary framework remain as interesting subjects of future research. Finally, we provide in the Appendix of this work some theoretical results concerning oblique projections and the regularity of the distance to , which play a key role in the definition of the scheme and in the proof of its main properties.

## Acknowledgements

The first two authors would like to thank the Italian Ministry of Instruction, University and Research (MIUR) for supporting this research with funds coming from the PRIN Project (KKJPX entitled “Innovative numerical methods for evolutionary partial differential equations and applications”). Xavier Dupuis thanks the support by the EIPHI Graduate School (contract ANR-17-EURE-0002). Elisa Calzola, Elisabetta Carlini and Francisco J. Silva were partially supported by KAUST through the subaward agreement OSR--CRG-..

## 2. Preliminaries

As mentioned in the introduction, it will be simpler to describe our approximation scheme when (1.1) is written in backward form. This can be done by a simple change of the time variable and a possible modification of the time dependency of . Let us set and . We consider the HJB equation

 (HJB)

where and are respectively given by (1.2) and (1.3).

For notational convenience, throughout this article, we will write for all and . Our standing assumptions for the data in (HJB) are the following.

• () is a nonemtpy, bounded domain with boundary of class .

• The functions , , , and are continuous. Moreover, for every , the functions and are Lipschitz continuous, with Lipschitz constants independent of .

• The function is of class . We also assume that

 (∀(x,b)∈∂O×B)|γb(x)|=1%and⟨n(x),γb(x)⟩>0,

where, for every , we recall that denotes the unit outward normal vector to at .

We now recall the notion of viscosity solution to (see [6]). We need first to introduce some notation. Given a bounded function , its upper semicontinuous (resp. lower semicontinuous) envelope is defined by

 (2.1) (∀(t,x)∈¯¯¯¯OT)z∗(t,x):=limsup(s,y)∈¯¯¯¯OT,(s,y)→(t,x)z(s,y)⎛⎜ ⎜⎝resp. % z∗(t,x):=liminf(s,y)∈¯¯¯¯OT,(s,y)→(t,x)z(s,y)⎞⎟ ⎟⎠.
###### Definition 2.1.

[Viscosity solution]

(i) An upper semicontinuous function is a viscosity subsolution to (HJB) if for any and such that has a local maximum at , we have

 (2.2) −∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x))≤0,

if ,

 (2.3) min{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),L(t,x,Dϕ(t,x))}≤0,

if and,

 (2.4) u1(t,x)≤Ψ(x),

if .

(ii) A lower semicontinuous function is a viscosity supersolution to (HJB) if for any and such that has a local minimum at , we have

 (2.5) −∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x))≥0,

if ,

 (2.6) max{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),L(t,x,Dϕ(t,x))}≥0,

if and,

 (2.7) u2(t,x)≥Ψ(x),

if .

(iii) A bounded function is a viscosity solution to (HJB) if and , defined in (2.1), are, respectively, sub- and supersolutions to (HJB).

###### Remark 2.1.

As shown in [12, Proposition 6], relation (2.4) can be replaced by

 (2.8) min{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),u1(t,x)−Ψ(x)}≤0,

if , and

 (2.9) min{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),L(t,x,Dϕ(t,x)),u1(t,x)−Ψ(x)}≤0,

if . Similarly, condition (2.7) can be replaced by

 (2.10) max{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),u2(t,x)−Ψ(x)}≥0,

if , and

 (2.11) max{−∂tϕ(t,x)+H(t,x,Dϕ(t,x),D2ϕ(t,x)),L(t,x,Dϕ(t,x)),u2(t,x)−Ψ(x)}≥0,

if .

The following well-posedness result for has been shown in [6, Theorem II.1] (see also [11]).

###### Theorem 2.1.

Assume (H1)-(H3). Then there exists a unique viscosity solution to (HJB).

###### Remark 2.2.

(i) [Comparison principle and uniqueness] The existence of at most one solution to (HJB) follows from the following comparison principle (see [6, Theorem II.1] and also [11, Proposition 3.4]). If is a bounded viscosity subsolution to (HJB) and is a bounded viscosity supersolution to (HJB), then

 u1≤u2in¯¯¯¯OT.

(ii) [Existence] Once a comparison principle has been shown, the existence of a solution to (HJB) follows usually from the existence of sub- and supersolutions to (HJB) and Perron’s method. In Sect. 5, we construct sub- and supersolutions to (HJB) as suitable limits of solutions to the approximation scheme that we present in the next section. Together with the comparison principle, this yields an alternative existence proof of solutions to (HJB).

An different and interesting technique to show the existence of a solution to (HJB) is to consider a suitable stochastic optimal control problem, with controlled reflection of the state trajectory at the boundary , and to show that the associated value function is a viscosity solution to (HJB). This strategy has been followed in [11].
(iii) [Continuity] The continuity of the unique viscosity solution to (HJB) follows directly from the comparison principle and the continuity properties required in the definition of sub- and supersolutions to (HJB). Notice that, as usual for parabolic problems with Neumann type boundary conditions, we do not require any compatibility condition between and the operator at the boundary .

## 3. The fully discrete scheme

We introduce in this section a fully discrete SL scheme that approximates the unique viscosity solution to . Throughout this section, we assume that (H1)-(H3) are fulfilled.

### 3.1. Discretization of the space domain O

Let us fix and consider a polyhedral domain such that

 (3.1) d(O,OΔx)=inf{|x−y||x∈O,y∈OΔx}≤C(Δx)2,

for some . A construction of such a domain can be found in [8, Section 3] for or , which explain the dimension constraint in (H1). However, the results in the remainder of this article can be extended to , provided that a numerical domain satisfying (3.1) exists. Let be a triangulation of consisting of simplicial finite elements with vertices in (for some ). We assume that is the mesh size, i.e. the maximum of the diameters of , all the vertices on belong to , at most one face of each element , with vertices on , intersects , and satisfies the following regularity condition: there exists , independent of , such that each is contained in a ball of radius and contains a ball of radius . As in [18], we introduce an auxiliary exact triangulation of with vertices in . The boundary elements of are allowed to be curved and we have

 ¯¯¯¯O=⋃ˆT∈TΔxˆT.

Denoting by the projection on , the projection is defined by

 pΔx(x)=pT(x), if x∈ˆT∈ˆTΔx and the element T∈TΔx has the same vertices than ˆT.

Set and denote by the linear finite element basis function on . More precisely, for each , is a continuous function, affine on each , , , for all , with , and for all . For any

its linear interpolation

on the mesh is defined by

 (3.2) I[ϕ](x):=NΔx∑i=1ψi(pΔx(x))ϕ(xi),for all x∈¯¯¯¯O.
###### Lemma 3.1.

Let and denote by its restriction to . Then there exists a constant , independent of , such that

 (3.3) supx∈¯¯¯¯O∣∣ϕ(x)−I[ϕ|GΔx](x)∣∣≤Cϕ(Δx)2.
###### Proof.

Let and let and be two elements having the same vertices and such that . By the triangular inequality

 |ϕ(x)−I[ϕ|GΔx](x)|≤|ϕ(x)−ϕ(pT(x))|+|ϕ(pT(x))−I[ϕ|GΔx](x)|.

Using that is Lipschitz, we deduce from (3.1) the existence of , independent of and , such that

. In addition, by standard error estimates for

interpolation (see for instance [15]) and (3.2), there exists , independent of and , such that . Relation (3.3) follows from these two estimates. ∎

### 3.2. A semi-Lagrangian scheme

Let , set , and . We define the time grid . Given , , and , we define the discrete characteristics

 (3.4) y±,ℓk,i(a)=xi+Δtμ(tk,xi,a)±√NσΔtσℓ(tk,xi,a).

Let and let be a fixed constant. For any we set

 (∂O)δ:={x∈RN|d(x,∂O)<δ}.

By Proposition 7.1 in the Appendix, there exist and two functions and , uniquely determined, such that

 (3.5) x=pγb(x)+dγb(x)γb(pγb(x)),% for all (x,b)∈(∂O)R×B.

Therefore, there exists such that for all , , , , and , the reflected characteristic

 (3.6) ~ysk,i(a,b):={ysk,i(a)if ysk,i(a)∈¯¯¯¯O,pγb(ysk,i(a))−¯c√Δtγb(pγb(ysk,i(a)))otherwise

is well-defined. In Figure 1 we illustrate how the reflected characteristic is computed from the projection of onto parallel to . Let us also set

 (3.7) ~dsk,i(a,b) :={0if ysi,k(a)∈¯¯¯¯O,dγb(ysk,i(a))+¯c√Δtotherwise, (3.8) ~gsk,i(a,b) :=⎧⎨⎩0if ysk,i(a)∈¯¯¯¯O,g(tk,pγb(ysk,i(a)),b)% otherwise.

Notice that if , then (3.5), (3.6), and (3.7) imply that

 (3.9) ~ysk,i(a,b)=ysk,i(a)−~dsk,i(a,b)γb(pγb(ysk,i(a))).

For and , let us define by

 (3.10) Sk,i[Φ](a,b):=12Nσ∑s∈I[I[Φ](~ysk,i(a,b))+~dsk,i(a,b)~gsk,i(a,b)]+Δtf(tk,xi,a)

and set

 (3.11) Sk,i[Φ]:=infa∈A,b∈BSk,i[Φ](a,b).

In the remainder of this work, we will consider the following fully discrete SL scheme to approximate the solution to .

 (HJB{\rm\tiny disc}) Uk,i=Sk,i[Uk+1,(⋅)],% for(k,i)∈I∗Δt×IΔx,UNΔt,i=Ψ(xi),fori∈IΔx.

### 3.3. Probabilistic interpretation of the scheme

The fully-discrete SL to approximate the solution to (HJB) in the unbounded case, i.e. , has a natural interpretation in terms of a discrete time, finite state, Markov control process (see e.g. [13, Section 3]). We show below that a similar interpretation holds for (HJB). The latter will play an important role in the stability analysis of (HJB) presented in the next section. Given and , , let us define the controlled transition law

 (3.12) pk,i,j(a,b):=12Nσ∑s∈Iβj(~ysk,i(a,b)),for all i,j∈IΔx.

We say that is a -policy if for all we have . The set of -policies is denoted by . Let us fix and, for notational convenience, set . Associated to and

, there exists a probability measure

on (the powerset of

) and a Markov chain

, with state space , such that

 (3.13) Pk,xi,π(Xk=xi)=1andPk,xi,π(Xm+1=xj|Xm=xi)=pm,i,j(πm(xi)),

for . Now, consider a family of

-valued independent random variables, which are also independent of

, and with common distribution given by

 P(ξm=±eℓ)=12Nσ,for m=k+1,…,NΔt and ℓ=1,…,Nσ,

where denotes the -th canonical vector of . By a slight abuse of notation (see (3.4)), for , , and , let us set

 (3.14) ym(xi,a)=xi+Δtμ(tm,xi,a)+√NσΔtσ(tm,xi,a)ξm+1.

For , , , and , define the random variable

 (3.15) h(tm,xi,a,b)=⎧⎨⎩0if ym(xi,a)∈¯¯¯¯O,(dγb(ym(xi,a))+¯c√Δt)g(tm,pγb(ym(xi,a)),b)otherwise.

For all and , let us define

 Jk,i(π)=EPk,xi,π(∑NΔt−1m=k[Δtf(tm,Xm,αm)+h(tm,Xm,αm,βm)]+Ψ(XNΔt)),JNΔt,i(π)=Ψ(xi),

where, for notational convenience, we have denoted, respectively, by and the first and the last coordinates of . Notice that, by construction and (3.10), we have that

 Jk,i(π)=\SSk,i[Jk+1,(⋅)(π)](αk,βk).

Moreover, setting

 ^Uk,i=infπ∈ΠNΔtJk,i(π),^UNΔt,i=Ψ(xi),

for all , the dynamic programming principle (see e.g. [23, Theorem 12.1.5]) implies that satisfies (HJB). Since the latter has a unique solution, we deduce that for all and .

###### Remark 3.1.

Scheme (HJB) can thus be interpreted as a Markov chain discretization of an stochastic control problem with oblique reflection in the boundary (see e.g. [11]).

## 4. Properties of the fully discrete scheme

In this section, we establish some basic properties of (HJB).

###### Proposition 4.1.

The following hold:
(i) (Monotonicity) For all with , we have

 Sk,i[U]≤Sk,i[V],for k∈I∗Δtand i∈IΔx.

(ii) (Commutation by constant) For any and ,

 Sk,i[U+c]=Sk,i[U]+c,for k∈I∗Δtand i∈IΔx.
###### Proof.

Both assertions follow directly from (3.10) and (HJB). ∎

We show in Proposition 4.2 below a consistency result for (HJB). For this purpose, let us set

 H(t,x,p,M,a) = −12{\rm Tr}(σ(t,x,a)σ(t,x,a)⊤M)−⟨μ(t,x,a),p⟩−f(t,x,a), for (t,x,p,M,a)∈¯¯¯¯OT×RN×RN×Nσ×A, L(t,x,p,b) = ⟨γ(x,b),p⟩−g(t,x,b), for (t,x,p,b)∈[0,T]×∂O×RN×B,

and for all , , , , , and , define

 (4.3) ~Lsk,i(q,a,b):=⎧⎪⎨⎪⎩0if ysk,i(a)∈¯¯¯¯O,L(tk,pγb(ysk,i(a)),q,b)% otherwise.
###### Proposition 4.2 (Consistency).

Let and denote by its restriction to . Then the following hold:

1. For all , , , and , we have

 Sk,i[ϕ|GΔx](a,b)−ϕ(xi)=−ΔtH(tk,xi,Dϕ(xi),D2ϕ(xi),a)−12Nσ∑s∈I~dsk,i(a,b)(~Lsk,i(Dϕ(xi),a,b)−√ΔtKsk,i(a,b))+O(Δt√Δt+(Δx)2),

where the set of constants is bounded, independently of .

2. For all and , we have

 Sk,i[ϕ|GΔx]−ϕ(xi)=−supa∈A,b∈B{ΔtH(tk,xi,Dϕ(xi),D2ϕ(xi),a)+12Nσ∑s∈I~dsk,i(a,b)(~Lsk,i(Dϕ(xi),a,b)−√ΔtKsk,i(a,b))}+O(Δt√Δt+(Δx)2).
###### Proof.

In what follows, we denote by a generic constant, which is independent of , , , , and . Since assertion (ii) follows directly from , we only show the latter.

For every , (3.4) and (3.7) imply that . Thus, by (3.4), (3.9), and a second order Taylor expansion of around , for every , we have

 ϕ(~y±,ℓk,i(a,b))=ϕ(xi)+Δt⟨Dϕ(xi),μ(tk,xi,a)⟩+NσΔt2⟨D2