Symplectic algorithms for stable manifolds in control theory

It is well known that an optimal feedback control can be given by solving an associated Hamilton-Jacobi (HJ) equation (see e.g. [23]) and $H^{\infty}$ feedback control can be obtained from solutions of one or two Hamilton-Jacobi equations (see e.g. [37, 20, 38, 6]). Unfortunately, the Hamilton-Jacobi equation in general can not be solved analytically. Hence numerical method becomes important. Seeking approximate solutions of Hamilton-Jacobi equations from control theory has been studied extensively. There are several approaches: Taylor series method, Galerkin method, state-dependent Riccati equation method, algebraic method, etc. See e.g. [3, 24, 21, 8, 7, 27, 25, 9, 4, 5, 26, 28, 19, 29, 2] and the references therein. These methods may have good performance for concrete control systems. However, in general, they may have various disadvantages such as heavy computation cost for higher dimensional state spaces, restriction on simple nonlinearity of the systems, etc.

For the stationary Hamilton-Jacobi equations which are related to infinite horizon optimal control and $H^{\infty}$ control problems, [32] developed an iterative procedure to construct an approximate sequence that converges to the exact solution of the associated Hamiltonian system on the stable manifold. It is based on the fact that the stabilizing solutions of stationary Hamilton-Jacobi equations correspond to the generating functions of the stable manifolds (Lagrangian) of the associated Hamiltonian systems at certain equilibriums (cf. e.g. [25, 30, 32]). This approach has better performances for various nonlinear feedback control systems, especially for the ones with more complicated nonlinearities, see e.g. [31, 18, 17].

We should note that the computation approach in [32] (as well as [31, 18, 17]) depends essentially on the radius of convergence of the iterative procedure which is not estimated analytically. Moreover, since the errors of less iterative steps are tremendous when we enlarge the local approximate stable manifolds in unstable direction, to obtain a stable manifold with proper size for applications, the number of iterative steps should be large. This may make the computation time-consuming.

In this note, inspired by [32], we construct a sequence of local approximate stable manifolds of stationary Hamilton-Jacobi equations by iterative procedure with a precise estimate of radius of convergence, and enlarge the local manifolds to larger ones by symplectic algorithms. To be more precise, as in [32] we firstly generate a sequence of local approximate stable manifolds near the equilibrium by an iterative procedure to solve the associated Hamiltonian system of the Hamilton-Jacobi equation. Then we extend the local approximate stable manifolds to large ones by symplectic algorithms for solving initial value problem of the associated Hamiltonian system (Section 4 below). We emphasize that in our approach, the radius of convergence of the iterative sequences are estimated precisely and the error of the local approximate stable manifold can be controlled as small as possible (Theorem 3.1 below). Therefore the significant point is how to extend the local stable manifolds. There are various numerical methods for general ODEs, e.g. Runge-Kutta methods of various orders. For our applications in Hamiltonian systems, we will use symplectic algorithms. The symplectic structure plays a significant role in design of numerical methods for Hamiltonian systems. Generally speaking, the symplectic algorithms are designed to preserve the symplectic structure of the Hamiltonian systems. Hence, for Hamiltonian systems, symplectic algorithm improves qualitative behaviours, and gives a more accurate long-time integration comparing with general-purpose methods such as Runge-Kutta schemes. Various kinds of symplectic algorithms, e.g., symplectic Euler, Störmer-Verlet, symplectic Runge-Kutta, were constructed since 1950s. A detailed history of symplectic methods and related topics can be found in [14]. We refer the readers to the books [16, 13, 10, 34] for a complete presentation of various symplectic algorithms for Hamiltonian systems.

The note is organized as follows. In Section 2, basic notations for the stable manifolds of Hamilton-Jacobi equations in control theory are given. Section 3 is devoted to constructing the iterative procedure, and proves the precise estimate of radius of convergence as well as the error of the approximate solutions. The symplectic scheme which extends the local approximate stable manifold is described in Section 4. In Section 5, the effectiveness of our algorithm is illustrated by application to an optimal control problem. The Appendix gives the details of the proof of Theorem 3.1.

In this section, for the convenience of the reader we recall the relevant material without proofs, thus making our exposition self-contained. For more details, see [32] and [31].

Let $\Omega\subset\mathbb{R}^{d}$ be a domain containing $0$. We consider the following Hamilton-Jacobi equation

where $p:=\nabla V$ for some unknown function $V$, $f:\Omega\to\mathbb{R}^{d}$, $R:\Omega\to\mathbb{R}^{d\times d}$, $q:\Omega\to\mathbb{R}$ are $C^{\infty}$ and $R(x)$ is symmetric matrix for all $x\in\Omega$. Furthermore, we assume that $f(0)=0$ and $q(0)=0$, $\frac{\partial q}{\partial x}(0)=0$. Hence for $x$ near $0$, $f(x)=Ax+O(|x|^{2}),\,q(x)=\frac{1}{2}x^{T}Qx+O(|x|^{3}),$ where $A=\frac{\partial f}{\partial x}(0)$, $Q=\frac{\partial^{2}q}{\partial x^{2}}(0)$ is the Hessian of $q$ at $0$. We say that a solution $V$ of (2.1) is stabilizing if $p(0)=0$ and $0$ is an asymptotically stable equilibrium of the vector field $f(x)-R(x)p(x)$ where $p(x)=\nabla V(x)$.

This kind of problem arises from infinite horizon optimal control. For example, consider the optimal control problem

where $g$ is a smooth $d\times m$ matrix-valued function, $u$ is an $m$-dimensional feedback control function of column form. Assume the instantaneous cost function $L(x,u)=\frac{1}{2}(q(x)+ru^{T}u),$ where $r>0$, $q(x)\geq 0$ is a smooth function. Define the cost functional by $J(x,u)=\int_{0}^{+\infty}L(x(t),u(t))dt.$ Then the corresponding Hamilton-Jacobi equation is of form (2.1) with $R(x):=r^{-1}g(x)g(x)^{T}$. The optimal controller is

where $p(x)=\nabla V(x)$.

From the symplectic geometry point of view, a stabilizing solution $V$ of (2.1) corresponds to a stable (Lagrangian) submanifold. That is, $\Lambda_{V}:=\{(x,p)\,|\,p=\nabla V\}$ is a stable (Lagrangian) manifold which is invariant under the flow of the associated Hamiltonian system of (2.1):

Conversely, if a $d$-dimensional manifold $\Lambda$ in $(x,p)$-space is invariant with respect to the flow (2.6), and at some point $(x_{0},p_{0})$, the projection of $\Lambda$ to the $x$-space is surjective, then $\Lambda$ is a Lagrangian submanifold in a neighborhood of $(x_{0},p_{0})$ and there is a solution $V$ of (2.1) in a neighborhood of $x_{0}$ such that $\Lambda_{V}=\Lambda$. Therefore, if we get the stable manifold, then the optimal controller can be obtained from (2.3). See e.g. [32].

Denote $z=(x,p)$. Let $J$ be the standard symplectic matrix $\left[\begin{array}[]{cc}0&I_{d}\\ -I_{d}&0\end{array}\right]$, where $I_{d}$ denotes the identity matrix of $d$-dimensional. Then the vector field on left side of (2.6) is $J^{-1}\nabla H(z)$. Let $X_{H}(z)=J^{-1}\nabla H(z)$ be the Hamiltonian vector field of $H$. Note that $z_{0}=(0,0)$ is an equilibrium of $X_{H}$. Then the derivative of the Hamiltonian vector field at $z_{0}$ is a Hamiltonian matrix, i.e., $(JDX_{H}(z_{0}))^{T}=JDX_{H}(z_{0})$. We say that the equilibrium $z_{0}$ is hyperbolic if $DX_{H}(z_{0})$ has no imaginary eigenvalues. It is well known that if $DX_{H}(z_{0})$ is hyperbolic, then its eigenvalues are symmetric with respect to the imaginary axis. From the Stable Manifold Theorem, there exists a global stable manifold $\mathcal{S}_{z_{0}}$ of $z_{0}$. Moreover, $\mathcal{S}_{z_{0}}$ is a Lagrangian submanifold of $(\mathbb{R}^{2d},\omega)$ where $\omega$ is the standard symplectic structure (see e.g. [1]). Near $z_{0}$, the Hamiltonian system (2.6) can be rewritten as

where $N(z)$ is the nonlinear term.

A sufficient condition for the existence local stabilizing solution for (2.1) is obtained by van der Schaft [37] based on an observation on the Riccati equation. Assume $(x_{0},p_{0})=(0,0)$ without loss of generality. Let $R:=R(0)$. The Riccati equation

is the linearization of (2.1) at the origin. A symmetric matrix $P$ is said to be the stabilizing solution of (2.8) if it is a solution of (2.8) and $A-RP$ is stable. The Riccati equation (2.8) has a stabilizing solution if and only if the following two conditions hold: (1) $DX_{H}(z_{0})$ is hyperbolic; (2) the generalized eigenspace $E_{-}$ for $d$ stable eigenvalues satisfies $E_{-}\oplus{\rm Im}(0,I_{d})^{T}=\mathbb{R}^{2d}.$ See e.g. [22]. If (2.8) has a stabilizing solution $P$, then there exists a local stabilizing solution $V$ of (2.1) around the origin such that $\frac{\partial^{2}V}{\partial x^{2}}(0)=P$. This yields a local solution of Hamilton-Jacobi solution ([32]). Consider the Lyapunov equation

Some direct computations yield the following result ([32]).

From Lemma 2.1, by coordinates transformation $\left[\begin{array}[]{c}\bar{x}\\ \bar{p}\\ \end{array}\right]=T^{-1}\left[\begin{array}[]{c}x\\ p\\ \end{array}\right]$, system (2.7) becomes a separated form

where $N_{s}(\bar{x},\bar{p})$ and $N_{u}(\bar{x},\bar{p})$ are nonlinear terms corresponding to $N(z)$ after the coordinates transformation.

In this section, we shall give an iterative procedure to construct a sequence of local approximate stable manifolds near equilibrium for the Hamiltonian system in the form

Assumption 1: $B$ has eigenvalues with negative real parts. It follows that there exist positive constants $a,b$ such that $\|e^{Bt}\|\leq ae^{-bt}$ for $t\geq 0$.

Assumption 2: $n_{s},n_{u}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ are continuous and satisfy the following conditions: For all $L>0$, $0<l\leq L$, $|x|+|y|\leq l$ and $|x^{\prime}|+|y^{\prime}|\leq l$,

where $M(L)$ is increasing with respect to $L$.

As in [32], to solve equation (3.3), we define the following iterative sequence for $k=1,2,\cdots$,

and $x_{0}=e^{Bt}\xi$, $y_{0}=0$ for an arbitrary $\xi\in\mathbb{R}^{d}$. Equivalent to (3.6), we consider the following ODE:

with boundary conditions $x_{k+1}(0)=\xi,\,y_{k+1}(+\infty)=0$ and $x_{0}=e^{Bt}\xi$, $y_{0}=0$, $t\geq 0$. This form is more convenient to apply numerical methods to ODEs.

Inspired by [32, Theorem 5], we can prove the following convergence result whose proof is included in Appendix.

In this section, the local stable manifolds will be enlarged by symplectic algorithms.

In this section, we apply the symplectic algorithm to a 2-dimensional optimal control problem with exponential nonlinearity. The existed methods may be difficult to applied to the systems with such kind of complicated nonlinearities as showed in [31, 18, 17].

Throughout this section, we use the following notations. $k$: the iterative times for local approximate stable manifold as in (3.9); $\xi$: the initial condition given in Theorem 3.1; $h_{-}$: the step size of the numerical methods for $t<0$ by (4.7). RK45 method is from scipy.integrate.RK45 in python environment. This algorithm is based on error control as 4-order Runge-Kutta method, and its steps are taken using the 5-order accurate formula. For more details this algorithm, see [11, 35]. We should mention that other numerical environment such as Matlab also includes a similar algorithm named by ode45. In what follows, we implement RK45 in iterative procedure (3.9) for the local stable manifold.

Let us consider the following optimal control problem with exponential nonlinearity:

where $u=(u_{1},u_{2})^{T}$ is the feedback control function. Let

Define the cost function by

Then the corresponding Hamilton-Jacobi equation is

where $p=\nabla V$ is the gradient of the value function in column form, $x=(x_{1},x_{2})^{T}$, $R=I_{2}$, $Q=I_{2}$. Then by (2.3) the optimal controller is $u=-p$. The associated Hamiltonian system is

where $\left[\frac{\partial f}{\partial x}(x)\right]^{T}=\left[\begin{array}[]{cc}0&-1-x_{1}^{2}\\ e^{x_{2}}&0\\ \end{array}\right]$. Then $(0,0)$ is an equilibrium and the Hamiltonian matrix is given by $\left[\begin{array}[]{cc}A&-R\\ -Q&-A^{T}\\ \end{array}\right],$ where $A=\left[\begin{array}[]{cc}0&1\\ -1&0\\ \end{array}\right].$ The stabilizing solution $P$ of the Riccati equation is $I_{2}$. From Lemma 2.1, we find a matrix $T=\left[\begin{array}[]{cc}I_{2}&-0.5I_{2}\\ I_{2}&0.5I_{2}\\ \end{array}\right]$ such that the Hamiltonian system (5.6) becomes a separated form

where $P$ is the stabilizing solution of the Riccati equation (2.8), $B=A-RP$, $x=x(\bar{x},\bar{p}),p=p(\bar{x},\bar{p})$ are defined by $\left[\begin{array}[]{c}x\\ p\\ \end{array}\right]=T\left[\begin{array}[]{c}\bar{x}\\ \bar{p}\\ \end{array}\right],$ and $\left[\begin{array}[]{c}n_{s}(\bar{x},\bar{p})\\ n_{u}(\bar{x},\bar{p})\\ \end{array}\right]:=T^{-1}\left[\begin{array}[]{c}f(x)-Ax\\ -\left[\frac{\partial f}{\partial x}\right]^{T}p+A^{T}p\\ \end{array}\right].$ Here the eigenvalues of $B$ are $-1\pm i$. Then we consider an iterative procedure as (3.9) with $\bar{x}_{0}=e^{Bt}\xi,\,\bar{p}_{0}=0$, and $\bar{x}_{k}(0)=\xi,\,\bar{p}_{k}(+\infty)=0$ for $k=1,2,3,\cdots$.

For (5.15) we can choose $M(L)=(9/4)L$ for $L>2/3$ and $M(L)=3/2$ for $L\leq 2/3$. Moreover, $a=1,\,b=1$. Hence by Theorem 3.1, if $M(L)L\geq 3/8$, then $|\bar{x}_{k}|+|\bar{p}_{k}|<L$ for all $k=1,2,\cdots$. That is, we can choose $L\geq 1/4$. Hence for $|\xi|\leq\frac{3}{16M(L)}\approx 0.125$, $\{\bar{x}_{k}(t,\xi)\}$ and $\{\bar{p}_{k}(t,\xi)\}$ converge to exact solution $\bar{x}(t,\xi)$ and $\bar{p}(t,\xi)$ respectively. In what following, we shall take $\xi$ in the sphere $\mathbb{S}_{0.12}$. That yields a local approximate stable manifold $\Lambda_{k}=\{(\bar{x}_{k}(t,\xi),\bar{p}_{k}(t,\xi))\,|\,t\geq 0,\,\xi\in\mathbb{S}_{0.12}\}$ and $\partial\Lambda_{k}=\{(\xi,\bar{p}_{k}(0,\xi))\,|\,\xi\in\mathbb{S}_{0.12}\}$.

Next, as in Section 4.3, the local stable manifold $\Lambda_{k}$ will be enlarged to a global one $\Lambda_{k,g}$ from (4.7). We solve the initial value problem (4.7) by the Störmer-Verlet scheme (4.17). In our problem (5.6), algorithm (4.17) becomes

Here the second equation is the only implicit equation which can be solved by Newton’s iteration method.

In this note, using a symplectic algorithm for the stable manifolds of the Hamilton-Jacobi equations from infinite horizon optimal control, we construct a sequence of local approximate stable manifolds for Hamiltonian system at some hyperbolic equilibrium. In our approach, we first construct an iterative sequence of local approximate stable manifolds at equilibrium based on a precise estimates for the radius of convergence. Then, using symplectic algorithms, we enlarge the local approximate stable manifolds by solving initial value problem of the associated Hamiltonian system for $t<0$. This avoids the divergent case of the iterative sequence as in [32] for unstable direction $t<0$, and the computation cost can be reduced. Moreover, the symplectic algorithms have better long-time behaviours than general purpose numerical methods such as Runge-Kutta. The effectiveness of our algorithm is demostrated by an optimal control problem with exponential nonlinearity. Lastly, we should mention that the symplectic algorithm has potentials to be applied to more general Hamilton-Jacobi equations in control theory since such kind of equations have associated Hamiltonian systems as the characteristic systems.