Finite Element Methods with Artificial Diffusion for Hamilton-Jacobi-Bellman Equations

## Abstract

In this short note we investigate the numerical performance of the method of artificial diffusion for second-order fully nonlinear Hamilton-Jacobi-Bellman equations. The method was proposed in (M. Jensen and I. Smears, arxiv:1111.5423); where a framework of finite element methods for Hamilton-Jacobi-Bellman equations was studied theoretically. The numerical examples in this note study how the artificial diffusion is activated in regions of degeneracy, the effect of a locally selected diffusion parameter on the observed numerical dissipation and the solution of second-order fully nonlinear equations on irregular geometries.

## 1 Introduction

Hamilton-Jacobi-Bellman (HJB) equations, which describe the value function in the theory of optimal control, are fully nonlinear partial differential equations, which are of second-order if the underlying control problem is stochastic. One challenge arising in the numerical solution of these equations is the presence of spurious generalised solutions of the PDE which do not coincide with the value function. While these spurious solutions often possess the same regularity as the value function, they violate monotonicity properties exhibited by the value function. These properties lead to the concept of viscosity solutions.

Regarding the numerical solution of HJB equations, we would like to highlight three approaches within the finite element methodology, which have been employed to ensure convergence to the value function. For a review of discrete Markov chain approximations before application of the dynamic programming principle we refer to [5]. For the method of vanishing moments we point to [3]. For the approach by Barles and Souganidis we cite the original source [1] and the recent adaption [4] by the authors to a finite element framework. A more comprehensive outline of the literature is in [2].

## 2 Numerical Method

Let be a bounded Lipschitz domain in , . Let be a compact metric space and

 A→C(¯¯¯¯Ω)×C(¯¯¯¯Ω,Rd)×C(¯¯¯¯Ω)×C(¯¯¯¯Ω),α↦(aα,bα,cα,dα)

be continuous. Consider the bounded linear operators of non-negative characteristic form

 Lα:H2(Ω)∩H10(Ω)→L2(Ω),w↦−aαΔw+bα⋅∇w+cαw

where belongs to . Then

 supα∈A∥(aα,bα,cα,dα)∥C(¯¯¯Ω)×C(¯¯¯Ω,Rd)×C(¯¯¯Ω)×C(¯¯¯Ω)+supα∈A∥Lα∥C2(¯¯¯Ω)→C(¯¯¯Ω)<∞. (1)

We assume that the final-time boundary data is non-negative: on . For smooth let , where the supremum is applied pointwise. The HJB equation considered is

 −vt+Hv =0 in (0,T)×Ω, (2a) v =g on (0,T)×∂Ω, (2b) v =vT on {T}×¯¯¯¯Ω. (2c)

Let , , be a sequence of piecewise linear shape-regular finite element spaces, whose underlying meshes are strictly acute. Let , , denote the nodes of the mesh with associated hat functions . Set . The mesh size is denoted . The set of time steps is . Let the th entry of be .

For each and find an approximate splitting into linear operators

 Eαi :w↦−¯aαiΔw+¯bαi⋅∇w+¯cαiw, Iαi :w↦−¯¯aαiΔw+¯¯bαi⋅∇w+¯¯cαiw,

of the form , , and . To impose monotonicity, select the artificial diffusion parameters and as in [4] such that and .

Define, for , ,

 (Eαiw)ℓ :=¯aαi(yℓi)⟨∇w,∇^ϕℓi⟩+⟨¯bαi⋅∇w+¯cαiw,^ϕℓi⟩, (3a) (Iαiw)ℓ :=¯¯aαi(yℓi)⟨∇w,∇^ϕℓi⟩+⟨¯¯bαi⋅∇w+¯¯cαiw,^ϕℓi⟩, (3b) (Cαi)ℓ :=⟨dαi,^ϕℓi⟩. (3c)

Obtain the numerical solution by interpolation of . Then at time is defined, inductively, by interpolating the boundary data and by

 −divi(ski,⋅)+supα(Eαivi(sk+1i,⋅)+Iαivi(ski,⋅)−Cαi)=0, (4)

where the supremum is understood to be applied component-wise to the collection of vectors

 {Eαivi(sk+1i,⋅)+Iαivi(ski,⋅)−Cαi:α∈A}.

## 3 Selection of the Artificial Diffusion Parameter

That the diffusion coefficients and in (3) are placed outside of the scalar product originates from the non-divergence form of the linear operators of HJB operators [4]. This structure makes it straightforward to implement local artificial diffusion parameters and , that is to implement a dependence on the nodal position . Indeed a change of corresponds in the assembly of to a scalar multiplication of the th row of the stiffness matrix, see (3).

Changes in the optimal artificial diffusion coefficients arise from variations of the local mesh quality and the local mesh Péclet number. The parameters and may be chosen by studying (35) in [4] or by examining the assembled matrices of the unstabilised operators. This is in particular simple for the , since only the signs of off-diagonal entries need be corrected to impose a local monotonicity property, thus leading to an algorithm which can be executed row-by-row.

Notice that in contrast, for problems in divergence form, the stabilised diffusion coefficients and need to be determined in the whole domain and not just at nodal positions.

## 4 An Exact Solution and Convergence Rates

We consider a triangular spatial domain with the vertices , and and the HJB equation

 −vt−12√x2+y2T−t+1Δv+121√T−t+1|∇v|=−12√x2+y2(T−t+1)3/2. (5)

To see that equation (5) is a HJB equation, note that the Euclidean norm of the gradient satisfies

 |∇v|=sup{β⋅∇v:β∈R2 with |β|=1}. (6)

A calculation verifies that the function

 v(x,y,t)=exp⎛⎝−√x2+y2T−t+1⎞⎠+√x2+y2T−t+1

is an exact solution of (5), where boundary and final-time conditions are determined by interpolation of the exact solution . The equation is to be solved on the time interval .

We consider a splitting which discretises the advection term explicitly with the minimum amount of diffusion needed for monotonicity. The remaining (linear) diffusion is incorporated in the implicit term, observing that this leads to a time-step restriction . Figure 2 shows the supremum norm error with a constant ratio , showing the uniform stability of the method in this setting.

Figure 1 illustrates a choice of artificial diffusion that is locally adapted to the Péclet number of a coarse mesh. It is seen that diffusion is only artificially introduced in a neighbourhood of the origin, which is where the operator becomes degenerate. Figure 3 illustrates the rates of convergence of the numerical solution in the -, - and -norms to the exact solution at the initial time, now using the constant ratio .

## 5 Eikonal Equation

The steady-state limit of the time-dependent eikonal equation,

 −vt+|∇v|=1,

equipped with homogeneous boundary and final-time conditions, measures the distance to the boundary. Due to (6) the eikonal equation belongs to the Hamilton-Jacobi-Bellman family. We consider the equation on a domain (Figure 4, top left) whose irregular shape leads to complicated curves on which is not differentiable. The height of the domain is equal to one.

Figure 4 compares locally-adapted choices of the artificial diffusion parameter with global choices, and illustrates the effect of mesh refinement on the quality of the solution. To compare the quality of the numerical solutions we compare their -norms—recalling that excessive numerical dissipation leads to a smearing out of extrema. The coarse grid solution, with internal nodes and with a global diffusion parameter of , only reaches a height of . In contrast, a computation on the same mesh with a local diffusion parameter, with mean value and standard deviation but maximal value , leads to a height of . With one (two) steps of uniform refinement the number of elements increases by a factor of (of ) and the artificial diffusion decreases to an average value of (of ), giving the improved value of (of ) for the -norm in the lower left (right) part of Figure 4. The -norm of a reference solution on a very fine mesh is .

## 6 A Second-Order Fully Nonlinear Equation

The final example, on the same domain as in the previous section, is concerned with the fully nonlinear case of second-order. We examine the equation

 −vt+supα∈[α0,α1]{−αΔv}+|∇v|=−vt+sup(α,β)∈[α0,α1]×∂B(0,1){−αΔv+β⋅∇v}=f,

where is the unit ball in , , , and

 f(x,y) =529(sin(g(x,y))+12sin(2g(x,y))+410sin(8g(x,y)))2, g(x,y) =π2(x−0.63)(y−0.26)/0.07.

The boundary and final-time conditions are homogeneous. As before, the advection term is discretised explicitly, with the locally minimal diffusion needed for monotonicity. The possibly remaining (fully nonlinear) diffusion is placed in the implicit term.

In this example the control of the second-order term is maximised independently of the first-order control ; in the sense that the optimal may be determined without knowledge of . Furthermore, as takes locally either a positive or negative value only the controls and are ever active in the HJB equation. This is an example of the Bang-bang principle. It is also reflected in the left plot in Figure 5, where the value of is plotted. Black colouring signifies maximises the operator, whereas white colouring corresponds to . Observe that no intermediate grey values can appear. The plot of the control mimics some of the features of the value function, which is plotted in the same figure on the right—in part because the Laplacian contains information about the curvature of the solution.

At each time-step of the method, a semi-smooth Newton method [4] was used to solve the nonlinear discrete equation in (4), where each iteration of the algorithm involves solving a linear system. To study the performance of the algorithm, the HJB equation was solved on a sequence of successively refined meshes, with a constant set of tolerances. The linear systems were solved to a tolerance of by GMRES. The stopping criterion for Newton’s method was a relative residual tolerance of in the maximum norm and a convergence of the iterations requirement of in the maximum norm. The sizes of the linear systems for the sequence of meshes were , , , and . The respective average number of Newton iterations for a single time step were , , , and , with respective standard deviations , , , and . This demonstrates a weak dependence of the number of iterations needed for an individual time step on the system size, thus showing that the total number of linear systems to be solved for a complete computation depends principally on the number of time-steps.

### References

1. G. Barles and P.E. Souganidis, Convergence of approximation schemes for fully nonlinear second order equations, J. Asymptotic Analysis, 4:271–283, 1991.
2. X. Feng and R. Glowinsky and M. Neilan, Recent developments in numerical methods for fully nonlinear second order partial differential equations, submitted to SIAM Review.
3. X. Feng and M. Neilan, The vanishing moment method for fully nonlinear second order partial differential equations: formulation, theory, and numerical analysis, arXiv:1109.1183, 2011.
4. M. Jensen and I. Smears, On the convergence of finite element methods for Hamilton-Jacobi-Bellman equations, arXiv:1111.5423, 2011.
5. H.J. Kushner and P. Dupuis, Numerical methods for stochastic control problems in continuous time, Applications of Mathematics 24, Springer-Verlag, 2001.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters