1 Introduction↩︎

Physics-Informed Neural Network(PINN) is a powerful tool for solving Partial Differential Equation(PDE) and even for identifying the missing parameters. Since the conception of PINN was first proposed by M. Raissi etc. in [1], lots of PINN-based applications have been created and showed us its great power in solving PDEs that come from various practical problems, e.g. fluid mechanics, medical diagnosis, etc.[2], [3]. Among all the problems, the one that draws great intension for us is the solving of the Hamiltonian-Jacobi-Equation (HJ equation). It is worth noticing that HJ equation has a deep connection with the Regular equations in classical mechanics. On the one hand, the existence of the HJ equation comes from the transformation of Hamilton quantity in order to solve the regular equations in Hamiltonian system, on the other hand, the regular equations are the characteristic lines of the HJ equation.

Another important application of the HJ equation originates from optimal control, where the Hamiltonian has the form of minimum of a family of linear functions. In this case the HJ equation has the name of HJB where the B represents the Bellman in honor of his contribution to principle of optimality. With the HJB equation, it is easy to find the optimal cost for a finite or infinite horizon optimal control problem. Actually, the difference between the HJB and the popular reinforcement learning technique lies in the knowledge of system or environment model, and also if we know the initial condition when start learning. With an accurate model of the system, we can theoretically solve the minimum cost and the corresponding optimal control very easily.

But there is a huge gap between reality and theory: the difficulty of solving the HJB by using traditional numerical methods such as Finite Difference (FD) method and Semi-Lagrangian(SL) method will increase exponentially with the increase of the state dimension of the system. This problem is called the "curse of dimensionality"[4]. Since PINN doesn’t need to discretize the differential (FD method) or search step by step through discrete dynamic programming (SL method), PINN has theoretically a better performance especially when come across high dimensional system. Some of the existing work using PINN for solving HJB equations could be seen in [5]–[7]

Among all these research, a problem is that none of those research considered the solution in viscosity sense. Actually there exists more than one locally Lipschitz continuous function that satisfy the HJB almost everywhere. The key point is to find the viscosity solution that corresponds to the optimal control. Otherwise, if the PINN is trained directly with the loss term from the HJB, the neural network will converge to other solutions instead the viscosity one, unless the start point of training happens to be close to the answer. The existing solution is to first train the NN to a "guess" solution, which is usually based on the quadratic form of system states. However, it is obvious that this guess method can not guarantee to be correct every time.

In this paper, two PINN-based viscosity solutions for HJB equation are proposed. One is adding an extra term in the loss function, the other is using directly a convex network. Section 2 and Section 3.2 gives the motivation of this paper and some basic introduction of the viscosity solution. Section 4 proves the equivalence of the unique solution of HJB and the convex property of the solution. Section 5 provides two methods to guarantee the uniqueness of solution from neural network and the constructive proof of the convex PINN. Section 6 gives the simulation results of the both methods.

2 Motivation↩︎

Before formally introducing the scheme of the viscosity solution of the HJB equation, it is of great importance to introduce the motivation of studying the solution in viscosity sense instead of classical sense. On the one hand, it is too strong a condition to require differentiable at every point in the domain of a function. Actually, the boundary value problem for a first order P.D.E. does not admit a global \(C^1\) solution in general. Hence it is interesting for us to study the solution in more generalised case, differentiable except limited points. On the other hand, there could be more than one solution even if we could find a smooth solution for our problem. An example is given below.

\[\begin{align} \dot{x} &= ax+bu \\ V^{*} &= \min_{u \in U} \int_{0}^{\infty} \left( x^2 +u^2\right) dt \end{align}\] where \(x\), \(u\) are the system states and control inputs in \(\mathbb{R}^{1}\). The HJB equation corresponding is given below \[H=\min_{u \in U}\left\{ V_x\dot{x}+x^2 +u^2 \right\}=0\]

with the boundary condition \(V(0)=0\). H could be rewritten as \[H=-\frac{1}{4} {V_x}^2b^2 + axV_x + x^2=0\]

From the above equation, it could be imagined that there are at lease two analytical solutions. Actually these solutions are

\[\begin{align} V_x &=2\frac{a \pm \sqrt{a^2+b^2}}{b^2}x \\ V(x) &= \frac{a \pm \sqrt{a^2+b^2}}{b^2} x^2 \end{align}\] we could see that both solutions satisfy the boundary condition, which leads to a natural question that which is the true solution corresponding to our optimal control? It is noticeable that in these two solutions, one is convex, and the other is concave. In our example, it is easy to prove that, the convex one \[\begin{align} V(x) &= \frac{a + \sqrt{a^2+b^2}}{b^2} x^2 \end{align}\] is our true solution.

For a more complex multi-dimensional and nonlinear problem, the solution of HJB could be more complicated. Actually for a nonlinear system, the generalized solution, i.e. locally Lipschitz continuous solution that satisfies HJB almost everywhere, doesn’t have uniqueness even the stability[8]. Next two sections would show that the optimal control problem under the context of concave Hamiltonian would lead to the unique semiconvex solution of equation, just like the example shows.

3 Viscosity solution of HJB↩︎

3.1 HJB for optimal control problem↩︎

Consider an affine system \[\begin{align} \label{affine32system} \dot{x} &= a(x)x+b(x)u, \quad u(t)\in U, \quad t\in \left[0,T\right] \end{align}\tag{1}\] with the cost function

\[\begin{align} \label{cost95function} V(t,x) &= \inf_{u \in U} \int_{0}^{T} \left(x^TQx +u^TRu\right) dt +\psi(x(T))=\inf_{u \in U} \int_{0}^{T} L(x,u) dt +\psi(x(T)) \end{align}\tag{2}\] where \(x \in R^n\) and \(u \in R^m\), and \(Q,R\) are positive definite matrix, and the set \(U\) is all of the admissible control values as defined below.

\[\begin{align} U:= \{ u: \mathbb{R} \to \mathbb{R}^m \;\;measureable, \;\;\;u(t)\in U \;for \;a.e. \;t \} \end{align}\]

Given the initial data \(x(s)=y\), consider the control u and the corresponding trajectory x, we have the following theorem.

3.2 Viscosity solution↩︎

Some useful properties of viscosity solution are listed below:

• Existence. For a nonlinear first order PDE in eq.@eq:HJB , the viscosity solution exists.

• Uniqueness. For such a problem, there exists at most one viscosity solution.

• Stability. The viscosity solution is stable in the following sense
  for any disturbance \(d>0\), \(V^d\) is a viscosity solution of \[\begin{align} \label{HJBd} V^d_t+H(V_x^d, x) & =0 & & \text{ in } {R}^n \times(0, \infty) \\ V^d & =g^d & & \text{ on } {R}^n \times\{t=T\} \nonumber \end{align}\tag{3}\] if parameters of the eq.@eq:HJBd converge to eq.@eq:HJB , e.g. \(H^d \to H\), \(g^d \to g\), and \(V^d \to V\) correspondingly, then \(V\) is the viscosity solution of eq.@eq:HJB .

The proof of the above properties could be find in lots of books [8], [9], [11]

It is worth mentioning that, another definition of the viscosity solution which is more commonly used for analysis is given below

4 Uniqueness of semiconvex solution in viscosity sense↩︎

This section proves the uniqueness of a semi convex solution that satisfy eq.@eq:HJB , given the assumption that Hamiltonian is concave.

4.1 Proof of concave Hamiltonian↩︎

4.2 Equivalence of convex solution and viscosity solution↩︎

Most reference gives the proof in the form of opposite assumptions and opposite conclusions[8], [10], e.g. when the Hamiltonian is convex, the unique viscosity solution is concave. This part of proof would try to give the opposite theorem that is more familiar to control community and the entire proof refers to the skill in[10].

The proof is given in [9] and [10]

5 PINN-based numerical viscosity solution of HJB↩︎

This section introduced two methods to ensure the convex structure of solution based on PINN method

5.1 Method 1: Addding convex term in loss function↩︎

The sufficient and necessary condition of the function \(V(t,.)\) being convex is \[\begin{align} Hess(V_x)=\begin{bmatrix} \frac{\partial^2V(t,x)}{\partial x_1^2} & \cdots & \frac{\partial^2V(t,x)}{\partial x_1x_n}\\ \vdots & \ddots & \vdots\\ \frac{\partial^2V(t,x)}{\partial x_nx_1} & \cdots & \frac{\partial^2V(t,x)}{\partial x_n^2} \end{bmatrix}>0 \\ \iff \left|Hess(V_x)_{i,i}\right|>0 \end{align}\] where \(\left|Hess(V_x)_{i,i}\right|\) is the ith order principal minor determinant.

Hence, the following loss function is thus proposed during training \[\begin{align} Loss = \frac{1}{N_i} \Sigma_{j=1}^{N_i}\left| V_{tj}+H(x_j,V_{xj}) \right|^2 +\frac{1}{N_b} \Sigma_{j=1}^{N_b} \left|V(tj,xj)-V_{NN}(tj,xj)\right|^2 + \frac{1}{N_H} \Sigma_{j=1}^{N_H}\left|\min (0,\left|Hess(V_x)_{j,j} \right|)\right|^2 \end{align}\] where \(N_i\) is the number of inner point, \(N_b\) is the number of boundary point, \(N_H\) is the number of point calculating Hessin matrix. It is worth mentioning that the Tanh activation function would be used in the neural network, hence the smoothness of the solution is guaranteed. Once we can check the positiveness of Hessian matrix of \(V_x\), then the uniqueness of solution is satisfied according to theorem 3.

5.2 Method 2: Design of input convex neural network↩︎

Another method is to construct a partially input convex neural network directly to ensure the convexity of solution in theorem3.

To build a first full input convex neural network, it is necessary to review some basic properties of convex functions

From the proposition above, it is easy to prove that the following neural network is a convex function

\[\begin{align} y = W_1\sigma(W_0x+b_0) + b_1\\ where \left\{\left. w_i>0\right|\forall w_i \in W_1 \nonumber \right\},\sigma(x)=\max(0,x) \nonumber \end{align}\] It is clear that this single neural network represents the sum of a bunch of convex functions. To improve the approximation ability, we add a feed forward channel to the output, so that the single layer NN become \[\begin{align} \label{ConvexNN} y = W_1\sigma(W_0x+b_0) + b_1 + fx\\ where\left\{\left. w_i>0\right| \forall w_i \in W_1 \nonumber \right\} \end{align}\tag{4}\] The structure is shown in Fig.1 and the following theorem holds

Next phase is to build the partially convex network. A partial convex neural network is proposed as Fig.2 shows, where the input output relation is given as

\[\begin{align} y &= W_2\sigma(W_0x + W_1{\mathcal{F}}_{{NN}_2}(\hat{y}) + b_0) + b_1\\ \hat{y} &= {\mathcal{F}}_{{NN}_1}(t) \end{align}\] where \(\left\{\left. w_i>0\right|\forall w_i \in W_2 \nonumber \right\},\sigma(x)=\max(0,x)\), and \({\mathcal{F}}_{{NN}_1},{\mathcal{F}}_{{NN}_2}\) are normal fully connected networks.

On the one hand, for a given \(t\), \(V(t,.)\) could be seen as a convex neural network in Fig.1 with a special bias decided by input t, hence \(V(t,.)\) is a convex neural network. On the other hand, for a given \(x\), \(W_0x\) could be regarded as a special bias for the network of \(V(t,.)\). It is worth mentioning that, although \(W_2\) should be element wise positive, there exists a network \({\mathcal{F}}_{{NN}_2}\) with certain layers, so that the input \(\hat{y}\) approximate output \(y\). Hence given a normal fully connected network \({\mathcal{F}}_{{NN}_1}\) with certain layers, the representability from input \(t\) to output \(y\) is not worth than a normal fully connected neural network.

5.3 Strip training method↩︎

A key issue during the training of PINN, is how to find global minimum. For HJB-based problem, we found a efficient training method that could avoid local minimum theoritically, we call it strip training method, which is similar to the curriculum learning in reinforcement learning.

For a given system \[\begin{align} \dot{x} &= a(x)x+b(x)u, \quad u(t)\in U, \quad t\in \left[0,T\right] \end{align}\] with the cost function \[\begin{align} V(t,x) &= \inf_{u \in U} \int_{0}^{T} \left(x^TQx +u^TRu\right) dt +\psi(x(T))=\inf_{u \in U} \int_{0}^{T} L(x,u) dt +\psi(x(T)) \end{align}\] where the cost function \(V(T,x)\), at end time T, is known. The corresponding HJB is given as

\[\begin{align} V_t + \inf_{u \in U} \{ L(x,u) + V_x\dot{x}\} = V_t + H(x,u,V_x) = 0 \end{align}\]

We have the strip training algorithm which is shown in algorithm 4. It worth mentioning that, the method 2 has a better performance than method 1 with less computation time in practice, so this strip training method is studied based on method 2. Although method 1 could also be realized in similar way. The sketch figure is shown in Fig.3

The main idea of strip training is, the optimal function \(V(t,x)\) is continuous, and if we start from the boundary, and take \(\Delta t\) small enough, the shape of neural network could approximate the true solution as close as we want. Then the training step is first train the boundary point, then under the shape of boundary points (realized by large \(\lambda\)), add the information of HJB equation to the loss (in step 2.2). Then we treat the trained area as a known area, so after adding a new strip to the training area (in step 3), the large \(\lambda\) is used to maintain the shape of known area. Repeat the procedure till the time reaches to 0.

6 Simulation results↩︎

Here we provide three nonlinear examples.

6.1 Example 1↩︎

\[V^{*} = \min_{u \in U} \int_{0}^{\infty} x^Tx +u^2 =0.5x_1^2+x_2^2\] where \[\begin{align} \dot{x_1} &= -x_1+x_2 \\ \dot{x_2} &= -0.5x_1-0.5x_2+0.5x_1^2x_2 +x_1u \\ \end{align}\] The corresponding HJB is \[\begin{align} x_1^2+x_2^2+(-x_1+x_2)V_{x_1}+(-0.5x_1-0.5x_2+0.5x_1^2x_2)V_{x_2} -0.25x_1^2V_{x_2}^2=0\\ V(0)=0 \end{align}\]

The training data is sampled with the range \([-1,-1]\) to \([1,1]\). results with both methods are shown in Fig.5. The orange surface is the Analytical solution, while the blue surface is the prediction of neural network. It could be seen that in this case, the method 2 has a better accuracy.

6.2 Example 2↩︎

\[V^{*} = \min_{u \in U} \int_{0}^{\infty} x^Tx +u^2 =0.5x_1^2+x_2^2\] where \[\begin{align} \dot{x_1} &= -x_1+x_2 \\ \dot{x_2} &= -0.5x_1-0.5x_2(1-(cos(2x_1)+2)^2)+(cos(2x_1)+2)u \\ \end{align}\] The corresponding HJB is \[\begin{align} x_1^2+x_2^2+(-x_1+x_2)V_{x_1}+(-0.5x_1-0.5x_2(1-(cos(2x_1)+2)^2))V_{x_2} -0.25(cos(2x_1)+2)^2V_{x_2}^2=0\\ V(0)=0 \end{align}\]

The training data is sampled with the range \([-1,-1]\) to \([1,1]\). results with both methods are shown in Fig.6 . The orange surface is the Analytical solution, while the blue surface is the prediction of neural network. In this case, because of the nonlinearity of system and the convex optimization term, the method 1 could not find the global optimum. While the method 2 find the solution directly.

6.3 Example 3↩︎

\[V^{*} = \min_{u \in U} \int_{0}^{\infty} x_2^2+u^2 =x_1^2(\frac{\pi}{2}+arctan(5x_1))+x_2^2\] where \[\begin{align} \dot{x_1} &= x_2 \\ \dot{x_2} &= -x_1(\frac{\pi}{2}+arctan(5x_1))-\frac{5x_1^2}{2+50x_1^2}+4x_2 +3u \\ \end{align}\] The corresponding HJB is \[\begin{align} x_2^2+x_2V_{x_1}+(-x_1(\frac{\pi}{2}+arctan(5x_1))-\frac{5x_1^2}{2+50x_1^2}+4x_2)V_{x_2} -\frac{9}{4} V_{x_2}^2=0\\ V(0)=0 \end{align}\] This example is the most hard one to be trained. The results are given in Fig.7. To show that we don’t start at some point that is already close to the true solution, the initial function has a really large distance from the true solution, and by method 2, we could see that the PINN could convergent to the true solution. But the accuracy is not as good as the two example before.

6.4 Example 4↩︎

For a linear system \[\begin{align} \dot{x}(t) &= x(t) + u(t) \end{align}\] the time limited quadratic cost function is given as \[V^{*}(t,x) = \min_{u \in U}V(t,x)= x^2(T) + \int_{t}^{T} u^2 dt=\frac{2x^2}{1+ e^{2t-2T}}\] where \(T= 10\) in this example. and corresponding HJB is \[\begin{align} V_t - 0.25V^2_x + V_xx = 0 \\ V(10,x)=x^2(10) \end{align}\]

The corresponding training result is shown in Fig.8. It could be seen that, for time limited case, the method 2 could well approximate the true solution of the HJB.

References↩︎