In this work, we propose a class of numerical schemes for solving semilinearHamilton-Jacobi-Bellman-Isaacs (HJBI) boundary value problems which arisenaturally from exit time problems of diffusion processes with controlled drift.We exploit policy iteration to reduce the semilinear problem into a sequence oflinear Dirichlet problems, which are subsequently approximated by a multilayerfeedforward neural network ansatz. We establish that the numerical solutionsconverge globally in the $H^2$-norm, and further demonstrate that thisconvergence is superlinear, by interpreting the algorithm as an inexact Newtoniteration for the HJBI equation. Moreover, we construct the optimal feedbackcontrols from the numerical value functions and deduce convergence. Thenumerical schemes and convergence results are then extended to HJBI boundaryvalue problems corresponding to controlled diffusion processes with obliqueboundary reflection. Numerical experiments on the stochastic Zermelo navigationproblem are presented to illustrate the theoretical results and to demonstratethe effectiveness of the method.
Quick Read (beta)
A neural network based policy iteration algorithm with global -superlinear convergence for stochastic games on domains\usetikzlibrary
Abstract. In this work, we propose a class of numerical schemes for solving semilinear Hamilton-Jacobi-Bellman-Isaacs (HJBI) boundary value problems which arise naturally from exit time problems of diffusion processes with controlled drift. We exploit policy iteration to reduce the semilinear problem into a sequence of linear Dirichlet problems, which are subsequently approximated by a multilayer feedforward neural network ansatz. We establish that the numerical solutions converge globally in the -norm, and further demonstrate that this convergence is superlinear, by interpreting the algorithm as an inexact Newton iteration for the HJBI equation. Moreover, we construct the optimal feedback controls from the numerical value functions and deduce convergence. The numerical schemes and convergence results are then extended to oblique derivative boundary conditions. Numerical experiments on the stochastic Zermelo navigation problem are presented to illustrate the theoretical results and to demonstrate the effectiveness of the method.
Key words. Hamilton-Jacobi-Bellman-Isaacs equations, neural networks, policy iteration, inexact semismooth Newton method, global convergence, -superlinear convergence.
AMS subject classifications. 82C32, 91A15, 65M12
In this article, we propose a class of numerical schemes for solving Hamilton-Jacobi-Bellman-Isaacs (HJBI) boundary value problems of the following form:
where is an open bounded domain, is the (nonconvex) Hamiltonian defined as
with given nonempty compact sets , and is a boundary operator, i.e., if is the identity operator, (1.1) is an HJBI Dirichlet problem, while if with some functions , (1.1) is an HJBI oblique derivative problem. Above and hereafter, when there is no ambiguity, we shall adopt the summation convention as in , i.e., repeated equal dummy indices indicate summation from to .
It is well-known that the value function of zero-sum stochastic differential games in domains satisfies the HJBI equation (1.1), and the optimal feedback controls can be constructed from the derivatives of the solutions (see e.g.  and references within; see also Section 6 for a concrete example). In particular, the HJBI Dirichlet problem corresponds to exit time problems of diffusion processes with controlled drift (see e.g. [32, 10, 37]), while the HJBI oblique derivative problem corresponds to state constraints (see e.g. [36, 34]). A nonconvex HJBI equation as above also arises from a penalty approximation of hybrid control problems involving continous controls, optimal stopping and impulse controls, where the HJB (quasi-)variational inequality can be reduced to an HJBI equation by penalizing the difference between the value function and the obstacles (see e.g. [27, 47, 39, 40]). As (1.1) in general cannot be solved analytically, it is important to construct effective numerical schemes to find the solution of (1.1) and its derivatives.
The standard approach to solving (1.1) is to first discretize the operators in (1.1) by finite difference or finite element methods, and then solve the resulting nonlinear discretized equations by using policy iteration, also known as Howard’s algorithm, or generally (finite-dimensional) semismooth Newton methods (see e.g. [18, 8, 45, 39]). However, this approach has the following drawbacks, as do most mesh-based methods: (1) it can be difficult to generate meshes and to construct consistent numerical schemes for problems in domains with complicated geometries; (2) the number of unknowns in general grows exponentially with the dimension , i.e., it suffers from Bellman’s curse of dimensionality, and hence this approach is infeasible for solving high-dimensional control problems. Moreover, since policy iteration is applied to a fixed finite-dimensional equation resulting from a particular discretization, it is difficult to infer whether the same convergence rate of policy iteration remains valid as the mesh size tends to zero ([42, 8]). We further remark that, for a given discrete HJBI equation, it can be difficult to determine a good initialization of policy iteration to ensure fast convergence of the algorithm; see  and references therein on possible accelerated methods.
Recently, numerical methods based on deep neural networks have been designed to solve high-dimensional partial differential equations (PDEs) (see e.g. [35, 15, 6, 16, 26, 44]). Most of these methods reformulate (1.1) into a nonlinear least-squares problem:
where is a collection of neural networks with a smooth activation function. Based on collocation points chosen randomly from the domain, (1.3) is then reduced into an empirical risk minimization problem, which is subsequently solved by using stochastic optimization algorithms, in particular the Stochastic Gradient Descent (SGD) algorithm or its variants. Since these methods avoid mesh generation, they can be adapted to solve PDEs in high-dimensional domains with complex geometries. Moreover, the choice of smooth activation functions leads to smooth numerical solutions, whose values can be evaluated everywhere without interpolations. In the following, we shall refer to these methods as the Direct Method, due to the fact that there is no policy iteration involved.
We observe, however, that the Direct Method also has several serious drawbacks, especially for solving nonlinear nonsmooth equations including (1.1). Firstly, the nonconvexity of both the deep neural networks and the Hamiltonian leads to a nonconvex empirical minimization problem, for which there is no theoretical guarantee on the convergence of SGD to a minimizer (see e.g. ). In practice, training a network with a desired accuracy could take hours or days (with hundreds of thousands of iterations) due to the slow convergence of SGD. Secondly, each SGD iteration requires the evaluation of (with respect to and ) on sample points, but is not necessarily defined everywhere due to the nonsmoothness of . Moreover, evaluating the function (again on a large set of sample points) can be expensive, especially when the sets A and B are of high dimensions, as we do not require more regularity than continuity of the coefficients with respect to the controls, so that approximate optimization may only be achieved by exhaustive search over a discrete coverage of the compact control set. Finally, as we shall see in Remark 4.2, merely including an -norm of the boundary data in the loss function (1.3) does not generally lead to convergence of the derivatives of numerical solutions or the corresponding feedback control laws.
In this work, we propose an efficient neural network based policy iteration algorithm for solving (1.1). At the th iteration, , we shall update the control laws by performing pointwise maximization/minimization of the Hamiltonian based on the previous iterate , and obtain the next iterate by solving a linear boundary value problem, whose coefficients involve the control laws . This reduces the (nonconvex) semilinear problem into a sequence of linear boundary value problems, which are subsequently approximated by a multilayer neural network ansatz. Note that compared to Algorithm Ho-3 in  for discrete HJBI equations, which requires to solve a nonlinear HJB subproblem (involving minimization over the set B) for each iteration, our algorithm only requires to solve a linear subproblem for each iteration, hence it is in general more efficient, especially when the dimension of B is high.
Policy iteration (or Successive Galerkin Approximation) was employed in [4, 5, 28, 30] to solve convex HJB equations on the whole space . Specifically, [4, 5, 28] approximate the solution to each linear equation via a separable polynomial ansatz (without concluding any convergence rate), while  assumes each linear equation is solved sufficiently accurately (without specifying a numerical method), and deduces pointwise linear convergence. The continuous policy iteration in  has also been applied to solve HJBI equations on in , which is a direct extension of Algorithm Ho-3 in  and still requires to solve a nonlinear HJB subproblem at each iteration. In this paper, we propose an easily implementable accuracy criterion for the numerical solutions of the linear PDEs which ensures the numerical solutions converge superlinearly in a suitable function space for nonconvex HJBI equations from an arbitrary initial guess.
Our algorithm enjoys the main advantage of the Direct Method, i.e., it is a mesh-free method and can be applied to solve high-dimensional stochastic games. Moreover, by utilizing the superlinear convergence of policy iteration, our algorithm effectively reduces the number of pointwise maximization/minimization over the sets A and B, and significantly reduces the computational cost of the Direct Method, especially for high dimensional control sets. The superlinear convergence of policy iteration also helps eliminate the oscillation caused by SGD, which leads to smoother and more rapidly decaying loss curves in both the training and validation processes (see Figure 7). Our algorithm further allows training of the feedback controls on a separate network architecture from that representing the value function, or adaptively adjusting the architecture of networks for each policy iteration.
A major theoretical contribution of this work is the proof of global superlinear convergence of the policy iteration algorithm for the HJBI equation (1.1) in , which is novel even for HJB equations (i.e., one of the sets A and B is singleton). Although the (local) superlinear convergence of policy iteration for discrete equations has been proved in various works (e.g. [38, 42, 18, 8, 47, 45, 39]), to the best of our knowledge, there is no published work on the superlinear convergence of policy iteration for HJB PDEs in a function space, nor on the global convergence of policy iteration for solving nonconvex HJBI equations.
Moreover, this is the first paper which demonstrates the convergence of neural network based methods for the solutions and their (first and second order) derivatives of nonlinear PDEs with merely measurable coefficients (cf. [22, 23, 26, 44]). We will also prove the pointwise convergence of the numerical solutions and their derivatives, which subsequently enables us to construct the optimal feedback controls from the numerical value functions and deduce convergence.
Let us briefly comment on the main difficulties encountered in studying the convergence of policy iteration for HJBI equations. Recall that at the th iteration, we need to solve a linear boundary value problem, whose coefficients involve the control laws , obtained by performing pointwise maximization/minimization of the Hamiltonian . The uncountability of the state space and the nonconvexity of the Hamiltonian require us to exploit several technical measurable selection arguments to ensure the measurability of the controls , which is essential for the well-definedness of the linear boundary value problems and the algorithm.
Moreover, the nonconvexity of the Hamiltonian prevents us from following the arguments in [42, 18, 8] for discrete HJB equations to establish the global convergence of our inexact policy iteration algorithm for HJBI equations. In fact, a crucial step in the arguments for discrete HJB equations is to use the discrete maximum principle and show the iterates generated by policy iteration converge monotonically with an arbitrary initial guess, which subsequently implies the global convergence of the iterates. However, this monotone convergence is in general false for the iterates generated by the inexact policy iteration algorithm, due to the nonconvexity of the Hamiltonian and the fact that each linear equation is only solved approximately. We shall present a novel analysis technique for establishing the global convergence of our inexact policy iteration algorithm, by interpreting it as a fixed point iteration in .
Finally, we remark that the proof of superlinear convergence of our algorithm is significantly different from the arguments for discrete equations. Instead of working with the sup-norm for (finite-dimensional) discrete equations as in [42, 18, 8, 47, 39], we employ a two-norm framework to establish the generalized differentiability of HJBI operators, where the norm gap is essential as has already been pointed out in [24, 46, 45]. Moreover, by taking advantage of the fact that the Hamiltonian only involves low order terms, we further demonstrate that the inverse of the generalized derivative is uniformly bounded. Furthermore, we include a suitable fractional Sobolev norm of the boundary data in the loss functions used in the training process, which is crucial for the -superlinear convergence of the neural network based policy iteration algorithm.
We organize this paper as follows. Section 2 states the main assumptions and recalls basic results for HJBI Dirichlet problems. In Section 3 we propose a policy iteration scheme for HJBI Dirichlet problems and establish its global superlinear convergence. Then in Section 4, we shall introduce the neural network based policy iteration algorithm, establish its various convergence properties, and construct convergent approximations to optimal feedback controls. We extend the algorithm and convergence results to HJBI oblique derivative problems in Section 5. Numerical examples for two-dimensional stochastic Zermelo navigation problems are presented in Section 6 to confirm the theoretical findings and to illustrate the effectiveness of our algorithms. The Appendix collects some basic results which are used in this article, and gives a proof for the main result on the HJBI oblique derivative problem.
2 HJBI Dirichlet problems
In this section, we introduce the HJBI Dirichlet boundary value problems of our interest, recall the appropriate notion of solutions, and state the main assumptions on its coefficients. We start with several important spaces used frequently throughout this work.
Let and be a bounded domain in , i.e., a bounded open connected subset of with a boundary. For each integer and real with , we denote by the standard Sobolev space of real functions with their weak derivatives of order up to in the Lebesgue space . When , we use to denote . We further denote by and the spaces of traces from and , respectively (see [21, Proposition 1.1.17]), which can be equivalently defined by using the surface measure on the boundaries as follows (see e.g. ):
We shall consider the following HJBI equation with nonhomogeneous Dirichlet boundary data:
where the nonlinear Hamiltonian is given as in (1.1):
Throughout this paper, we shall focus on the strong solution to (2.2), i.e., a twice weakly differentiable function satisfying the HJBI equation (2.2a) almost everywhere in , and the boundary values on will be interpreted as traces of the corresponding Sobolev space. For instance, on in (2.2b) means that the trace of is equal to in , where denotes the trace operator (see [19, Proposition 1.1.17]). See Section 5 for boundary conditions involving the derivatives of solutions.
We now list the main assumptions on the coefficients of (2.2).
Let , be a bounded domain, A be a nonempty finite set, and B be a nonempty compact metric space. Let , satisfy the following ellipticity condition with a constant :
and satisfy that on , and that is continuous, for all and .
As we shall see in Theorem 3.3 and Corollary 3.5, the finiteness of the set A enables us to establish the semismoothness of the HJBI operator (2.2a), whose coefficients involve a general nonlinear dependence on the parameters and . If all coefficients of (2.2a) are in a separable form, i.e., it holds for all that for some functions (e.g. the penalized equation for variational inequalities with bilateral obstacles in ), then we can relax the finiteness of A to the same conditions on B.
Finally, in this work we focus on boundary value problems in a domain to simplify the presentation, but the numerical schemes and their convergence analysis can be extended to problems in nonsmooth convex domains with sufficiently regular coefficients (see e.g. [21, 45]).
We end this section by proving the uniqueness of solutions to the Dirichlet problem (2.2) in . The existence of strong solutions shall be established constructively via policy iteration below (see Theorem 3.7).
Let be two strong solutions to (2.2), we consider the following linear homogeneous Dirichlet problem:
where we define the following measurable functions: for each ,
with the Hamiltonian defined as in (2.3). Note that the boundedness of coefficients implies that , and . Moreover, one can directly verify that the following inequality holds for all parametrized functions : for all ,
which together with (H.1) leads to the estimate that on the set , and hence we have a.e. . Then, we can deduce from Theorem A.1 that the Dirichlet problem (2.4) admits a unique strong solution and . Since satisfies (2.4) a.e. in and , we see that is a strong solution to (2.4) and hence , which subsequently implies the uniqueness of strong solutions to the Dirichlet problem (2.2). ∎
3 Policy iteration for HJBI Dirichlet problems
In this section, we propose a policy iteration algorithm for solving the Dirichlet problem (2.2). We shall also establish the global superlinear convergence of the algorithm, which subsequently gives a constructive proof for the existence of a strong solution to the Dirichlet problem (2.2).
We start by presenting the policy iteration scheme for the HJBI equations in Algorithm 3, which extends the policy iteration algorithm (or Howard’s algorithm) for discrete HJB equations (see e.g. [18, 8, 39]) to the continuous setting.
Choose an initial guess in , and set .
Given the iterate , update the following control laws: for all ,
Solve the linear Dirichlet problem for :
where for .
If , then terminate with outputs and , otherwise increment by one and go to step 2.
The remaining part of this section is devoted to the convergence analysis of Algorithm 3. For notational simplicity, we first introduce two auxiliary functions: for each with , we shall define the following functions
We then recall several important concepts, which play a pivotal role in our subsequent analysis. The first concept ensures the existence of measurable feedback controls and the well-posedness of Algorithm 3.
Let be a measurable space, and let and be topological spaces. A function is a Carathéodory function if:
for each , the function is -measurable, where is the Borel -algebra of the topological space ; and
for each , the function is continuous.
It is well-known that if are two complete separable metric spaces, and is a Carathéodory function, then for any given measurable function , the composition function is measurable (see e.g. [3, Lemma 8.2.3]). Since any compact metric space is complete and separable, it is clear that (H.1) implies that the coefficients are Carathéodory functions (with and ). Moreover, one can easily check that both and are Carathéodory functions, i.e., (resp. ) is continuous in (resp. ) and measurable in (see Theorem A.3 for the measurability of in ).
We now recall a generalized differentiability concept for nonsmooth operators between Banach spaces, which is referred as semismoothness in  and slant differentiability in [12, 24]. It is well-known (see e.g. [8, 45, 39]) that the HJBI operator in (2.2a) is in general non-Fréchet-differentiable, and this generalized differentiability is essential for showing the superlinear convergence of policy iteration applied to HJBI equations.
Let be defined on a open subset of the Banach space with images in the Banach space . In addition, let be a given a set-valued mapping with nonempty images, i.e., for all . We say is -semismooth in if for any given , we have that is continuous near , and
The set-valued mapping is called a generalized differential of in .
As in , we always require that has a nonempty image, and hence the -semismooth of in shall automatically imply that the image of is nonempty on .
Now we are ready to analyze Algorithm 3. We first prove the semismoothness of the Hamiltonian defined as in (2.3), by viewing it as the composition of a pointwise maximum operator and a family of HJB operators parameterized by the control . Moreover, we shall simultaneously establish that, for each iteration, one can select measurable control laws to ensure the measurability of the controlled coefficients in the linear problem (3.3), which is essential for the well-posedness of strong solutions to (3.3), and the well-definedness of Algorithm 3.
The following proposition establishes the semismoothness of a parameterized family of first-order HJB operators, which extends the result for scalar-valued HJB operators in . Moreover, by taking advantage of the fact that the operators involve only first-order derivatives, we are able to establish that they are semismooth from to for some (cf. [45, Theorem 13]), which is essential for the superlinear convergence of Algorithm 3.
Suppose (H.1) holds. Let be a given constant satisfying if and if , and let be the HJB operator defined by
Then is Lipschitz continuous and -semismooth in with a generalized differential
defined as follows: for any , we have
where is any jointly measurable function such that for all and ,
Since A is a finite set, we shall assume without loss of generality that, the Banach space is endowed with the usual product norm , i.e., for all , . Note that the Sobolev embedding theorem shows that the following injections are continuous: , for all , and , for all . Thus for any given satisfying the conditions in Proposition 3.1, we can find such that the injection is continuous. Then, the boundedness of implies that the mappings and are well-defined, and is Lipschitz continuous.
Now we show that the mapping has a nonempty image from to , where we choose such that the injection is continuous, and naturally extend the operators and from to . For each , we consider the Carathéodory function such that for all , where is defined by (3.5). Theorem A.3 shows there exists a function satisfying (3.8), i.e.,
and is jointly measurable with respect to the product -algebra on . Hence is nonempty for all .
We proceed to show that the operator is in fact -semismooth from to , which implies the desired conclusion due to the continuous embedding . For each , we denote by the -th component of , and by the -th component of . Theorem A.4 and the continuity of in show that for each , the set-valued mapping
is upper hemicontinuous, from which, by following precisely the steps in the arguments for [45, Theorem 13], we can prove that is -semismooth. Then, by using the fact that a direct product of semismooth operators is again semismooth with respect to the direct product of the generalized differentials of the components (see [46, Proposition 3.6]), we can deduce that is semismooth with respect to the generalized differential and finishes the proof. ∎
We then establish the semismoothness of a general pointwise maximum operator, by extending the result in  for the max-function .
Let be a given constant, A be a finite set, and be a bounded subset of . Let be the pointwise maximum operator such that for each ,
Then is -semismooth in with a generalized differential
defined as follows: for any , we have
where is any measurable function such that
Moreover, is uniformly bounded (in the operator norm) for all .
Let the Banach space be endowed with the product norm defined as in the proof of Proposition 3.1. We first show the mappings and are well-defined, has nonempty images, and is uniformly bounded for .
The finiteness of A implies that any can also be viewed as a Carathéodory function . Hence for any given , we can deduce from Theorem A.3 the existence of a measurable function satisfying (3.10). Moreover, for any given measurable function and , the function remains Lebesgue measurable (see Remark 3.1). Then, for any given with , one can easily check that , and , which subsequently implies that and are well-defined, and the image of is nonempty on . Moreover, for any , Hölder’s inequality leads to the following estimate:
which shows that for all .
Now we prove by contradiction that the operator is -semismooth. Suppose there exists a constant and functions such that as , and
where for each , is defined with some measurable function . Then, by passing to a subsequence, we may assume that for all , the sequence converges to zero pointwise a.e. in , as .
For notational simplicity, we define for all and . Then for a.e. , we have for all , for all . By using the finiteness of A and the convergence of , it is straightforward to prove by contradiction that for all such , for all large enough .
We now derive an upper bound of the left-hand side of (3.11). For a.e. , we have
from any . Thus, for each , we have for a.e. that,
where, by applying Theorem A.3 twice, we can see that both the set-valued mapping and the function are measurable.
We then introduce the set for each . The measurability of the set-valued mapping implies the associated distance function is a Carathéodory function (see [1, Theorem 18.5]), which subsequently leads to the measurability of for all . Hence we can deduce that
which leads to the following estimate:
where we have used the bounded convergence theorem and the fact that for a.e. , for all large enough . This contradicts to the hypothesis (3.11), and hence finishes our proof. ∎
Now we are ready to conclude the semismoothness of the HJBI operator. Note that the argument in  does not apply directly to the HJBI operator, due to the nonconvexity of the Hamiltonian defined as in (2.3).
Note that we can decompose the HJBI operator into , where is the linear operator , is the HJB operator defined in Proposition 3.1, is the pointwise maximum operator defined in Proposition 3.2, and is a constant satisfying if , and if .
Proposition 3.1 shows that is Lipschitz continuous and semismooth with respect to the generalized differential defined by (3.7), while Proposition 3.2 shows that is semismooth with respect to the uniformly bounded generalized differential defined by (3.9). Hence, we know the composed operator is semismooth with respect to the composition of the generalized differentials (see [46, Proposition 3.8]), i.e., for all . Consequently, by using the fact that is Fréchet differentiable with the derivative , we can conclude from Propositions 3.1 and 3.2 that is semismooth on , and that (3.12) is a desired generalized differential of at . ∎
Note that the above characterization of the generalized differential of the HJBI operator involves a jointly measurable function , satisfying (3.8) for all . We now present a technical lemma, which allows us to view the control law in (3.2) as such a feedback control on and .
Suppose (H.1) holds. Let , be given measurable functions, and be a measurable function such that for all ,
Then there exists a jointly measurable function such that for all , and it holds for all and that
Let be a given measurable function satisfying (3.14) for all (see Remark 3.1 and Theorem A.3 for the existence of such a measurable function). As shown in the proof of Proposition 3.1, there exists a jointly measurable function satisfying the property (3.15) for all . Now suppose that with (see (H.1)), we shall define the function , such that for all ,
The measurability of and the finiteness of A imply that the set is measurable in the product -algebra on , which along with the joint measurability of leads to the joint measurability of the function .
For any given , we have , from which we can deduce from the definition of that for all . Finally, for any given , we shall verify (3.15) for all and . Let be fixed. If , then the fact that and the definition of imply that , which along with (3.14) and shows that (3.15) holds for the point . On the other hand, if , then and satisfies the condition (3.15) due to the selection of . ∎
As a direct consequence of the above extension result, we now present an equivalent characterization of the generalized differential of the HJBI operator.
Let , and let and be given measurable functions satisfying (3.17) (see Remark 3.1 and Theorem A.3 for the existence of such measurable functions). Then by using Lemma 3.4, we know there exists a jointly measurable function such that satisfies (3.8) for all , and for all . Hence we see the linear operator defined in (3.16) is equal to the following operator
which is a generalized differential of the HJBI operator at due to Theorem 3.3. ∎
The above characterization of the generalized differential of the HJBI operator enables us to demonstrate the superlinear convergence of Algorithm 3 by reformulating it as a semismooth Newton method for an operator equation.
Note that the Dirichlet problem (2.2) can be written as an operator equation with the following operator
where is the HJBI operator defined as in (2.2a), and is the trace operator. Moreover, one can directly check that given an iterate , , the next iterate solves the following Dirichlet problem:
with the differential operator defined as in (3.16). Since is -semismooth (see Corollary 3.5) and , we can conclude that Algorithm 3 is in fact a semismooth Newton method for solving the operator equation .
Note that the boundedness of coefficients and the classical theory of elliptic regularity (see Theorem A.1) imply that under condition (H.1), there exists a constant , such that for any and any , the inverse operator is well-defined, and the operator norm is bounded by , uniformly in . Hence one can conclude from [46, Theorem 3.13] (see also Theorem A.5) that the iterates converges superlinearly to in a neighborhood of . ∎
The next theorem strengthens Theorem 3.6, and establishes a novel global convergence result of Algorithm 3 applied to the Dirichlet problem (2.2), which subsequently provides a constructive proof for the existence of solutions to (2.2). The following additional condition is essential for our proof of the global convergence of Algorithm 3:
Let the function in (H.1) be given as: , for all , where is a sufficiently large constant, depending on , , and .
In practice, (H.2) can be satisfied if (2.2) arises from an infinite-horizon stochastic game with a large discount factor (see e.g. ), or if (2.2) stems from an implicit (time-)discretization of parabolic HJBI equations with a small time stepsize.
If Algorithm 3 terminates in iteration , we have and , from which, we obtain from the uniqueness of strong solutions to (2.2) (Proposition 2.1) that is the strong solution to the Dirichlet problem (2.2). Hence we shall assume without loss of generality that Algorithm 3 runs infinitely.
for a.e. , where the function for all , and the modified Hamiltonian is defined as:
Hence, by taking the difference of equations corresponding to the indices and , one can obtain that
for , and on .
It has been proved in Theorem 9.14 of  that, there exist positive constants and , depending only on and , such that it holds for all with , and for all that
which, together with the identity that and the boundedness of , implies that the same estimate also holds for :
Thus, by assuming and using the boundedness of the coefficients, we can deduce from (3.20) that
for some constant independent of and the index .
Now we apply the following interpolation inequality (see [20, Theorem 7.28]): there exists a constant , such that for all and , we have . Hence, for any given , we have
Then, by taking , , and assuming that satisfies , we can obtain for that
which implies that is a Cauchy sequence with the norm .
We end this section with an important remark that, if one of the sets A and B is a singleton, and for all , then Algorithm 3 applied to the Dirichlet problem (2.2) is in fact monotonically convergent with an arbitrary initial guess. Suppose, for instance, that A is a singleton, then for each , we have that
Hence we can deduce that is a weak subsolution to , i.e.,
Thus, the weak maximal principle (see [19, Theorem 1.3.7]) and the fact that a.e. (with respect to the surface measure), leads to the estimate , which consequently implies that for all and a.e. .
4 Inexact policy iteration for HJBI Dirichlet problems
Note that at each policy iteration, Algorithm 3 requires us to obtain an exact solution to a linear Dirichlet boundary value problem, which is generally infeasible. Moreover, an accurate computation of numerical solutions to linear Dirichlet boundary value problems could be expensive, especially in a high-dimensional setting. In this section, we shall propose an inexact policy iteration algorithm for (2.2), where we compute an approximate solution to (3.3) by solving an optimization problem over a family of trial functions, while maintaining the superlinear convergence of policy iteration.
We shall make the following assumption on the trial functions of the optimization problem.
The collections of trial functions satisfies the following properties: for all , and is dense in .
It is clear that (H.3) is satisfied by any reasonable -conforming finite element spaces (see e.g. ) and high-order polynomial spaces or kernel-function spaces used in global spectral methods (see e.g. [4, 5, 13, 28, 29]). We now demonstrate that (H.3) can also be easily satisfied by the sets of multi-layer feedforward neural networks, which provides effective trial functions for high-dimensional problems. Let us first recall the definition of a feedforward neural network.
Definition 4.1 (Artificial neural networks).
Let be given constants, and be a given function. For each , let be an affine function given as for some and . A function defined as
is called a feedforward neural network. Here the activation function is applied componentwise. We shall refer the quantity as the depth of , as the dimensions of the hidden layers, and as the dimensions of the input and output layers, respectively. We also refer to the number of entries of as the complexity of .
Let , be some nondecreasing sequences of natural numbers, we define for each the set of all neural networks with depth , input dimension being equal to , output dimension being equal to 1, and dimensions of hidden layers being equal to . It is clear that if for all , then we have . The following proposition is proved in [25, Corollary 3.8], which shows neural networks with one hidden layer are dense in .
Let be an open bounded starshaped domain, and satisfying for all . Then the family of all neural networks with depth is dense in .
Now we discuss how to approximate the strong solutions of Dirichlet problems by reformulating the equations into optimization problems over trial functions. The idea is similar to least squares finite-element methods (see e.g. ), and has been employed previously to develop numerical methods for PDEs based on neural networks (see e.g. [35, 6, 44]). However, compared to [6, 35], we do not impose additional constraints on the trial functions by requiring that the networks exactly agree with the boundary conditions, due to the lack of theoretical support that the constrained neural networks are still dense in the solution space. Moreover, to ensure the convergence of solutions in the -norm, we include the -norm of the boundary data in the cost function, instead of the -norm used in  (see Remark 4.2 for more details).
For each , let be the unique solution to the Dirichlet problem (3.3):
where and denote the linear elliptic operator and the source term in (3.3), respectively. For each , we shall consider the following optimization problems:
The following result shows that the cost function provides a computable indicator of the error.
Let and . The definition of implies that , and . Then, by using the assumption that solves (3.3), we deduce that the residual term satisfies the following Dirichlet problem:
Hence the boundedness of coefficients and the regularity theory of elliptic operators (see Theorem A.1) lead to the estimate that
where the constants depend only on the -norms of , which are independent of . The above estimate, together with the facts that is dense in and , leads to the desired result that . ∎
We now present the inexact policy iteration algorithm for the HJBI problem (2.2), where at each policy iteration, we solve the linear Dirichlet problem within a given accuracy.
Choose a family of trial functions , an initial guess in , a sequence of positive scalars, and set .
If , then terminate with outputs and , otherwise increment by one and go to step 2.
In practice, the evaluation of the squared residuals in (4.2) depends on the choice of trial functions. For trial functions with linear architecture, e.g. if are finite element spaces, high-order polynomial spaces, and kernel-function spaces (see [9, 13, 28, 29]), one may evaluate the norms by applying high-order quadrature rules to the basis functions involved.
For trial functions with nonlinear architecture, such as feedforward neural networks, we can replace the integrations in by the empirical mean over suitable collocation points in and on , such as pseudorandom points or quasi-Monte Carlo points (see Section 6; see also [35, 6, 44]). In particular, due to the existence of local coordinate charts of the boundaries, we can evaluate the double integral in the definition of the -norm (see (2.1)) by first generating points in and then mapping the samples onto . The resulting empirical least-squares problem for the -th policy iteration step (cf. (4.1)) can then be solved by stochastic gradient descent (SGD) algorithms; see Section 6. We remark that, instead of pre-generating all the collocation points in advance, one can perform gradient descent based on a sequence of mini-batches of points generated at each SGD iteration. This is particularly useful in higher dimensions, where many collocation points may be needed to cover the boundary, and using mini-batches avoids having to evaluate functions at all collocation points in each iteration.
It is well-known (see e.g. [46, 14]) that the residual term is crucial for the superlinear convergence of inexact Newton methods. This next theorem establishes the global superlinear convergence of Algorithm 4.2.
Suppose (H.1), (H.2) and (H.3) hold, and in Algorithm 4.2. Let be the solution to the Dirichlet problem (2.2). Then for any initial guess , Algorithm 4.2 either terminates with for some , or generates a sequence that converges -superlinearly to in , i.e., . Consequently, we have for a.e. , and for all .
Let be an arbitrary initial guess. We first show that Algorithm 4.2 is always well-defined. For each , if is the strong solution to (3.3), then we can choose , which satisfies (4.2) and terminates the algorithm. If does not solve (3.3), the fact that is dense in enables us to find satisfying the criterion (4.2).
Moreover, one can clearly see from (4.2) that if Algorithm 4.2 terminates at iteration , then is the exact solution to the Dirichlet problem (2.2). Hence in the sequel we shall assume without loss of generality that Algorithm 4.2 runs infinitely, i.e., and for all .
as , where the additional high-order terms are due to the residuals and . Then, by using the interpolation inequality and assuming is sufficiently large, we can deduce that converge linearly to in .
We then reformulate Algorithm 4.2 into a quasi-Newton method for the operator equation , with the operator defined in the proof of Theorem 3.6. Let denote the strong dual space of , and denote the dual product on . For each , by using the fact that , we can choose satisfying , and introduce the following linear operators and :
Then, we can apply the identity and rewrite (4.3) as:
Now we are ready to establish the superlinear convergence of . For notational simplicity, in the subsequent analysis we shall denote by the Banach space with the usual product norm for each . By using the semismoothess of (see Theorem 3.6) and the strong convergence of in , we can directly infer from Theorem A.5 that it remains to show that there exists a open neighborhood of , and a constant , such that
as . Moreover, the boundedness of the coefficients shows that is Lipschitz continuous. Finally, the characterization of the generalized differential of in Theorem 3.6 and the regularity theory of elliptic operators (see Theorem A.1) show that for each , we can choose an invertible operator such that , uniformly in . Thus we can conclude from the semismoothness of at that
for all in some neighborhood of , which completes our proof for -superlinear convergence of .
Finally, we establish the pointwise convergence of and their derivatives. For any given , the superlinear convergence of implies that there exists a constant , depending on , such that for all . Taking the summation over the index , we have
where we used the monotone convergence theorem in the first equality. Thus, we have
which leads us to the pointwise convergence of and its partial derivatives with respect to . ∎
We reiterate that merely including the -norm of the boundary data in the cost functional (4.2) in general cannot guarantee the convergence of the derivatives of the numerical solutions , which can be seen from the following simple example. Let be a sequence such that in but not in , and for each , let be the strong solution to in and on .
The fact that in implies that in as . We now show in . Let be the solution to in and on , we can deduce from the integration by parts and the a priori estimate that
which shows that as . Now let be a given family of trial functions, which is dense in . One can find satisfying , and consequently in as . However, we have
Similarly, one can construct functions such that as , but does not converge to in .
We end this section with a convergent approximation of the optimal control strategies based on the iterates generated by Algorithm 4.2. For any given , we denote by and the set of optimal control strategies for all and for a.e. , such that
As an important consequence of the superlinear convergence of Algorithm 4.2, we now conclude that the feedback control strategies and generated by Algorithm 4.2 are convergent to the optimal control strategies.
Suppose the assumptions of Theorem 4.3 hold, and let be the solution to the Dirichlet problem (2.2). Assume further that there exist functions and such that and for a.e. . Then the measurable functions and , , generated by Algorithm 4.2 converge to the optimal feedback control pointwise almost everywhere.
Theorem A.4 implies that the set-valued mappings and are upper hemicontinuous. Then the result follows directly from the pointwise convergence of in Theorem 4.3, and the fact that and are singleton for a.e. . ∎
If we assume in addition that and for some Banach spaces and , then by using the compactness of A and B (see (H.1)), we can conclude from the dominated convergence theorem that in and in , for any .
5 Inexact policy iteration for HJBI oblique derivative problems
In this section, we extend the algorithms introduced in previous sections to more general boundary value problems. In particular, we shall propose a neural network based policy iteration algorithm with global -superlinear convergence for solving HJBI boundary value problems with oblique derivative boundary conditions. Similar arguments can be adapted to design superlinear convergent schemes for mixed boundary value problems with both Dirichlet and oblique derivative boundary conditions.
We consider the following HJBI oblique derivative problem:
where (5.1a) is the HJBI equation given in (2.2a), and (5.1b) is an oblique boundary condition. Note that the boundary condition on involves the traces of and its first partial derivatives, which exist almost everywhere on (with respect to the surface measure).
The following conditions are imposed on the coefficients of (5.1):
Assume , A, B, , satisfy the same conditions as those in (H.1). Let , , on , and assume there exists a constant , such that on , and on , where are the components of the unit outer normal vector field on .
The next proposition establishes the well-posedness of the oblique derivative problem.
We shall establish the uniqueness of strong solutions to (5.1) in this proof, and then explicitly construct the solution in Theorem 5.2 with the help of policy iteration; see also Theorem 3.7. Suppose that are two strong solutions to (5.1), then we can see is a strong solution to the following linear oblique derivative problem:
where is defined as in Proposition 2.1, and
By following the same arguments as the proof of Proposition 2.1, we can show that , and a.e. in , which, along with Theorem A.2, implies that is the unique strong solution to (5.2), and consequently in . ∎
Now we present the neural network based policy iteration algorithm for solving the oblique derivative problem and establish its rate of convergence.
Choose a family of trial functions , an initial guess in , a sequence of positive scalars, and set .
If , then terminate with outputs and , otherwise increment by one and go to step 2.
Note that the -norm of the boundary term is included in the cost function , instead of the -norm as in Algorithm 4.2. It is straightforward to see that Algorithm 5.1 is well-defined under (H.3) and (H.4). In fact, for each , given the iterate , Corollary 3.5 shows that one can select measurable control laws such that the following linear oblique boundary value problem has measurable coefficients:
and hence admits a unique strong solution in (see Theorem A.2). If , then solve the HJBI oblique derivative problem (5.1), and we can select and terminate the algorithm. Otherwise, the facts that and is dense in allows us to choose sufficiently closed to such that the criterion (5.3) is satisfied, and proceed to the next iteration.
The following result is analogue to Theorem 4.3, and shows the global superlinear convergence of Algorithm 5.1 for solving the oblique derivative problem (5.1). The proof follows precisely the lines given in Theorem 4.3, hence we shall only present the main steps in Appendix B for the reader’s convenience. The convergence of feedback control laws can be concluded similarly to Corollary 4.4 and Remark 4.3.
Suppose (H.2), (H.3) and (H.4) hold, and in Algorithm 5.1. Let be the solution to the oblique derivative problem (5.1). Then for any initial guess , Algorithm 5.1 either terminates with for some , or generates a sequence that converges -superlinearly to in , i.e., . Consequently, we have for a.e. , and for all .
6 Numerical experiments: Zermelo’s Navigation Problem
In this section, we illustrate the theoretical findings and demonstrate the effectiveness of the schemes through numerical experiments. We present a two-dimensional convection-dominated HJBI Dirichlet boundary value problem in an annulus, which is related to stochastic minimum time problems.
In particular, we consider the stochastic Zermelo navigation problem (see e.g. ), which is a time-optimal control problem where the objective is to find the optimal trajectories of a ship/aircraft navigating a region of strong winds, modelled by a random vector field. Given a bounded open set and an adaptive control strategy taking values in A, we assume the dynamics of the ship is governed by the following controlled dynamics:
where the drift coefficient is the sum of the velocity of the wind and the relative velocity of the ship, the nondegenerate diffusion coefficient describes a random perturbation of the velocity field of the wind, and is an -dimensional Brownian motion defined on a probability space .
The aim of the controller is to minimize the expected exit time of the region , taking model ambiguity into account in the spirit of . More generally, we consider the following value function:
over all admissible choices of , where denotes the first exit time of the controlled dynamics , the functions and denote the running cost and the exit cost, respectively, which indicate the desired destinations, and is a family of absolutely continuous probability measures with respect to with density , where is a predictable process satisfying for all and a given parameter . In other words, we would like to minimize a functional of the trajectory up to the exit time under the worst-case scenario, with uncertainty arising from the unknown law of the random perturbation.
By using the dual representation of and the dynamic programming principle (see e.g. [41, 10]), we can characterize the value function as the unique viscosity solution to an HJBI Dirichlet boundary value problem of the form (2.2). Moreover, under suitable assumptions, one can further show that is the strong (Sobolev) solution to this Dirichlet problem (see e.g. ).
For our numerical experiments, we assume that the domain is an annulus, i.e., , the wind blows along the positive -axis with a magnitude :
which decreases in terms of the distance from the bank, and the random perturbation of the wind is given by the constant diffusion coefficient . We also assume that the ship moves with a constant velocity , and the captain can control the boat’s direction instantaneously, which leads to the following dynamics of the boat in the region:
where represents the angle (measured counter-clockwise) between the positive -axis and the direction of the boat. Finally, we assume the exit cost on and on , which represents that the controller prefers to exit the domain through the inner boundary instead of the outer one (see Figure 6). Then the corresponding Dirichlet problem for the value function in (6.1) is given by: on , on , and
where and denote the -norm and -norm on , respectively. The optimal feedback control laws can be further computed as
6.1 Implementation details
In this section, we discuss the implementation details of Algorithm 4.2 for solving (6) with multi-layer neural networks (see Definition 4.1) as the trial functions. We shall now introduce the architecture of the neural networks, the involved hyper-parameters, and various computational aspects of the training process.
For simplicity, we shall adopt a fixed set of trial functions for all policy iterations, which contains fully-connected networks with the activation function , the depth , and the dimension of each hidden layer . The hyper-parameters and will be chosen depending on the complexity of the problem, which ensures that admits sufficient flexibility to approximate the solutions within the desired accuracy. More complicated architectures of neural networks with shortcut connections can be adopted to further improve the performance of the algorithm (see e.g. [16, 44]).
We then proceed to discuss the computation of the cost functional in (4.2) for each policy iteration. It is well-known that Sobolev norms of functions on sufficiently smooth boundaries can be explicitly computed via local coordinate charts of the boundaries (see e.g. ). In particular, due to the annulus shaped domain and the constant boundary conditions used in our experiment, we can express the cost functional as follows: for all and ,
where we define the map for . Note that we introduce an extra weighting parameter in (6.4), which helps achieve the optimal balance between the residual of the PDE and the residuals of the boundary data. We set the parameter for all the computations.
The cost functional (6.4) is further approximated by an empirical cost via the collocation method (see [35, 6]), where we discretize and by sets of collocation points and , respectively, and write the discrete form of (6.4) as follows: for all and ,
where , and for are, respectively, the Lebesgue measures of the domain and boundaries. Note that the choice of the smooth activation function implies that every trial function is smooth, hence all its derivatives are well-defined at any given point. For simplicity, we take the same number of collocation points in the domain and on the boundaries, i.e., .
It is clear that the choice of collocation points is crucial for the accuracy and efficiency of the algorithm. Since the total number of points in a regular grid grows exponentially with respect to the dimension, such a construction is infeasible for high-dimensional problems. Moreover, it is well-known that uniformly distributed pseudorandom points in high dimensions tend to cluster on hyperplanes and lead to a suboptimal distribution by relevant measures of uniformity (see e.g. [11, 6]). Therefore, we shall generate collocation points by a quasi-Monte Carlo (QMC) method based on low-discrepancy sequences. In particular, we first define points in from the generalized Halton sequence (see ), and then map those points into the annulus via the polar map , where and for all . The above transformation preserves fractiona