text:
# List of Figures

## Chapter 1 Introduction

We solve two equations: the Hamilton-Jacobi-Bellman (HJB) PDE and the Dynamic Programming Equations (DPE). Both equations are associated with the infinite horizon optimal control problem of minimizing a running cost subject to a nonlinear control dynamics. Since the optimal control problems arise ubiquitously in engineering, economics and biological science where models of these types surface naturally, this is the motivation to conjure up methods in solving them. It is actually an economist who first studied optimal control problems in modelling capital accumulation . However, it is the advent of modern control theory that had considerable impact on the treatment of these problems. The classic optimal control problem, the nonlinear regulator problem, is what we consider in this dissertation. We seek a feedback or a control which minimizes the running cost under the nonlinear control system for the regulator problem. The optimal control maintains the dynamics close to $`0`$ while keeping the expenditure of the cost at a minimum. The Dynamic Programming technique is used to derive the HJB PDE and the DPE from the infinite horizon optimal control problem. Thus, our main reason for solving the HJB PDE and DPE is to find a stabilizing feedback and a Lyapunov function to confirm local asymptotic stability of the given nonlinear control system. The HJB PDE corresponds to the optimal control problem with continuous dynamics and cost while the DPE is the discrete-time analogue.

Around 1961 one of the earliest results was made by Albrecht  on the nonlinear regulator problem. He discovered the formal power series solutions to the HJB equations while studying analytical nonlinear control systems. In 1969, Lukes proves the existence and uniqueness of these power series representations of the optimal control feedback and the optimal cost. To our knowledge there has been no extension of these results to the discrete-time case. So we look at the discrete-time version of the optimal control problem with the discrete-time nonlinear dynamics and cost. Using the method of Albrecht, we find the formal power series solutions to the DPE. Then we prove the existence and uniqueness of these power series solutions in Theorem 3.1.1. The Pontryagin Maximum Principle (PMP) ,  gives the necessary condition for a control to be optimal. The PMP also provides the Hamiltonian dynamics satisfied by the optimal state and the costate trajectories. The existence of the optimal cost is proved in two cases. One case is when we assume that the Hamiltonian map is invertible. The invertibility of the Hamiltonian matrix allows us to rewrite the dynamics where both the state and the costate dynamics propagate in the direction where time approaches infinity. The other case is when the Hamiltonian map is not a diffeomorphism; i.e., we have a nonlinear bidirectional Hamiltonian discrete-time dynamics. Finding the stable manifold for the Hamiltonian dynamics in both cases is key to finding the existence of the optimal cost. For invertible maps, Hartmann  shows the existence of a stable manifold by the method of successive approximations on the implicit functional equation. There is the technique of using the Contraction Mapping Theorem on a complete space developed by Kelley  in 1966. Applications of this method can be found in the book by Carr  and the paper of Krener . There is another technique by Irwin  based on an application of the inverse function theorem on a Banach space of sequences. One of our main results is Theorem 4.2.1; we show that there exists a local stable manifold for the bidirectional discrete map with a hyperbolic fixed point. After a two-step process of diagonalizing the dynamics we apply the technique of Kelly  on a complete space of Lipschitz functions with the supremum norm.

Computing the solutions of the HJB equations in higher dimension is major challenge. A standard approach requires the temporal and spacial discretizations of the optimal control problem and then solves the corresponding nonlinear program, see  and the appendix by Falcone in . Other methods for solving the HJB PDE are similar to those for conservation laws  and marching methods . Although these standard approaches work for such equations for low dimensional systems, these numerical methods become infeasible for real world problems which are higher dimensional systems. Our method is a higher order approach which requires few discretizations at each dimension. We start with Albrecht’s solutions as the initial approximations for the HJB PDE. We then truncate the computed solutions at the point where the optimal cost satisfies the optimality and stability constraints and begin new approximations at the same point. The idea is to patch together successive approximations to obtain a larger domain for which the numerically computed optimal cost is still a Lyapunov function. The numerical procedure is described in .

The outline of the dissertation is as follows: In Chapter 2, we derive the DPE. We first discuss the linear-quadratic regulator problem followed by the description of the method of finding formal power series solutions of the DPE. The coefficients of the polynomials are expressed in linear equations and the discrete-time algebraic Riccati equation (DTARE). We discuss the solvability of these equations as well. In Chapter 3, we prove the existence of the optimal solutions of the DPE. We show some properties of the Hamiltonian dynamics that are of importance in the proof of the existence of the optimal cost. Chapter 4 is where we prove the existence of a stable manifold for the noninvertible Hamiltonian map. It is actually the generalization of the results in Chapter 3. In Chapter 5, we summarize our numerical approach for solving the HJB PDE. We include numerical results from a 1-d example and describes the actual algorithm in more detail.

## Chapter 2 Formal Solution of the Dynamic Programming Equations in Discrete-Time

We present a method for the solution of the Dynamic Programming Equations, the discrete-time Hamilton Jacobi Bellman PDE that arises in an infinite time optimal control problem. The method extends to the discrete time Al’brecht’s procedure for locally approximating the solution. Assuming that the dynamics and cost are $`𝒞^{k1}\left(^{n+m}\right)`$ and $`𝒞^k\left(^{n+m}\right)`$ smooth, we explicitly find expansions of the optimal control and cost in a term by term fashion. These finite expansions are the first few terms the power series expansions of the optimal control and cost.

### 2.1 Preliminaries

First we discuss some definitions and notions about control systems. Consider the $`n`$-dimensional input and output state equations

$`x^+`$ $`=`$ $`Ax+Bu`$
$`y`$ $`=`$ $`Cx+Du.`$
where $`A,B,C,`$ and $`D`$ are $`n\times n,n\times m,l\times n,l\times m`$ matrices with $`x^n`$, $`y^l`$, and $`u^m`$. We use the notation $`x^+=x_{k+1}`$ and $`x=x_k`$.

###### Definition 2.1.1

A system is *controllable* when any inital state $`x_0`$ can be driven to the final state $`x_F`$ in finite number of steps. Equivalently, if the $`n\times nm`$ controllability matrix

$$𝒞=[B,AB,\mathrm{},A^{n1}B]$$

has rank $`n`$, then $`(A,B)`$ is a controllable pair. Otherwise (A,B) is said to be uncontrollable. The pair (A,B) is stabilizable if all uncontrollable modes are asymptotically stable.

If a state equation is controllable, then all eigenvalues can be assigned arbitrarily by introducing a state feedback. Moreover, every uncontrollable system can be diagonalized into

$`\left[\begin{array}{c}x_c^+\\ x_u^+\end{array}\right]=\left[\begin{array}{cc}A_c& \overline{A}\\ 0& A_u\end{array}\right]\left[\begin{array}{c}x_c\\ x_u\end{array}\right]+\left[\begin{array}{c}b_c\\ 0\end{array}\right]`$ (2.1.9)
where $`(A_c,b_c)`$ is controllable. If $`A_u`$ is stable and $`(A_c,b_c)`$ is controllable, then the system (2.1.9) is stabilizable. As in the definition above, we refer to the eigenvalues of $`A_u`$ as the uncontrollable modes.

Analogously, we define observability and dectectibility, the duals of controllabilily and stabilizability, respectively.

###### Definition 2.1.2

A sytem is observable if for any unknown $`x_0`$, $`T`$ s.t. the input $`u_k`$, and output $`y_k`$ uniquely determine $`x_0`$. Equivalently, if the $`nl\times n`$ observability matrix

$$𝒪=[C,CA,\mathrm{},CA^{n1}]^T$$

has rank $`n`$. Otherwise, $`(A,C)`$ is said to be unobservable. The pair $`(A,C)`$ is detectable if the unobservable eigenvalues are stable.

Given a control dynamics

$`x^+`$ $`=`$ $`Ax+Bu`$
$`x\left(0\right)`$ $`=`$ $`x_0`$
where the state $`x^n`$, the control $`u^m`$, $`A`$, $`B`$ are $`n\times n`$ and $`n\times m`$ matrices, respectively. We want to find $`u=Kx`$ such that the system is driven to 0; i.e., the system

$$x^+=\left(A+BK\right)x$$

is asymptotically stable around 0; i.e., the spectrum of $`A+BK`$ lies inside the unit circle. We call such control a stabilizing feedback. One way to solve this stabilization problem is to set up an optimal control problem. This will be discuss in the next sections in details.

### 2.2 Discrete-Time Optimal Control Problem

We formulate a discrete in time infinite horizon optimal control problem of minimizing the cost functional,

$$\underset{u}{\mathrm{min}}\underset{k=0}{\overset{\mathrm{}}{}}l(x_k,u_k)$$

subject to the dynamics

$`x^+`$ $`=`$ $`f(x,u)`$
$`x\left(0\right)`$ $`=`$ $`x_0`$
where the state vector $`x^n`$, the control $`u^m`$, and

$`f(x,u)`$ $`=`$ $`Ax+Bu+f^{\left[2\right]}(x,u)+f^{\left[3\right]}(x,u)+\mathrm{}`$
$`l(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Qx+x^{}Su+{\displaystyle \frac{1}{2}}u^{}Ru+l^{\left[3\right]}(x,u)+\mathrm{}`$
Here we denote $`f^{\left[m\right]}(x,u)`$ and $`l^{\left[m\right]}(x,u)`$ as homogeneous polynomials in $`x`$ and $`u`$ of degree $`m`$.

To solve the optimal control problem is to look for an optimal feedback $`u^{}=\kappa \left(x\right)`$ such that the cost functional is kept at its minimum, namely the optimal cost $`\pi \left(x_0\right)`$ of starting the system at $`x_0`$, while driving the dynamics to $`0`$.

Here are the assumptions: We assume $`l(x,u)`$ is convex in $`x`$ and $`u`$ so that

$`\left[\begin{array}{cc}Q& S\\ S^{}& R\end{array}\right]0`$
and let $`R>0.`$ In addition, it is assumed that the pair $`(A,B)`$ is stabilizable and the pair $`(A,Q^{1/2})`$ is detectable.

### 2.3 Derivation of the Dynamic Programming Equations

Given that $`x\left(0\right)=x_0`$ the optimal value function $`\pi \left(x\right)`$ is defined by

$$\pi \left(x_0\right)=\underset{u}{\mathrm{min}}\underset{k=0}{\overset{\mathrm{}}{}}l(x_k,u_k).$$

This value function satisfies a functional equation, called the dynamic programming equation. The optimal feedback $`\kappa \left(x\right)`$ is constructed from the dynamic programming equation . First, we state the optimality principle:

###### Theorem 2.3.1

Discrete-Time Optimality Principle:

$`\pi \left(x\right)=\underset{u}{\mathrm{min}}\left\{\pi \left(f(x,u)\right)+l(x,u)\right\}`$ (2.3.11)
Proof:  We have that

$`\pi \left(x_0\right)`$ $`=`$ $`\underset{u_0}{\mathrm{min}}\left\{{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}l(x_k,u_k)\right\}`$
$`=`$ $`\underset{u_0}{\mathrm{min}}\left\{l(x_0,u_0)+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}l(x_k,u_k)\right\}`$
$`=`$ $`\underset{u_0}{\mathrm{min}}\left\{l(x_0,u_0)+\pi \left(x_1\right)\right\}`$
Generalizing the optimality principle at the $`k^{th}`$-step, we have

$`\pi \left(x\right)=\underset{u}{\mathrm{min}}\left\{\pi \left(x^+\right)+l(x,u)\right\}.`$ (2.3.12)
The optimality equation (2.3.11) is the first equation of the dynamic programming equations. An optimal policy $`u^{}=\kappa \left(x\right)`$ must satisfy

$`\pi \left(x\right)\pi \left(f(x,u^{})\right)l(x,u^{})=0`$
if we assume convexity of the LHS of (2.3.12). We can find $`u^{}`$ through

$$\frac{(\pi \left(x\right)\pi \left(f(x,u)\right)l(x,u)}{u}=0$$

which by the chain rule becomes

$$\frac{\pi }{x}\left(f(x,u)\right)\frac{f}{u}(x,u)+\frac{l}{u}(x,u)=0$$

Thus, $`\pi \left(x\right)`$ and $`\kappa \left(x\right)`$ satisfy these equations, the Dynamic Programming Equations (DPE):

$`\pi \left(x\right)\pi \left(f(x,u)\right)l(x,u))`$ $`=`$ $`0`$ (2.3.13)
$`{\displaystyle \frac{\pi }{x}}\left(f(x,u)\right){\displaystyle \frac{f}{u}}(x,u)+{\displaystyle \frac{l}{u}}(x,u)`$ $`=`$ $`0`$ (2.3.14)
We now introduce a method for solving the dynamic programming equations for $`\kappa \left(x\right)`$ and $`\pi \left(x\right)`$.

### 2.4 Power Series Expansion

Our method of solving the DPE is an extension of Al’brecht idea for continuous-time systems . We require that $`f(x,u)𝒞^{k1}\left(^{n+m}\right)`$ and $`l(x,u)𝒞^k\left(^{n+m}\right)`$ are expressible in Taylor’s form around 0 in $`N_\epsilon \left(0\right)`$. Then, all the series representation of the given $`f(x,u),l(x,u)`$ and the unknown $`\kappa \left(x\right),\pi \left(x\right)`$ are substituted into the DPE. With the exception of the first level, gathering the terms of the same degree in each equation will yield linear equations. At the first level, we obtain the discrete-time algebraic Riccati equation. The advantage of our process is the reduction of DPE which is a nonlinear system of equation into a system of one Riccati equation and many linear equations. We will see that the set of linear equations has a triangular structure within the levels of degree, a structure that allows the system to be easily solvable.

On some neighborhood around 0, $`N_\epsilon \left(0\right)^{n+m}`$, the dynamics and cost assume the following form:

$`f(x,u)`$ $`=`$ $`Ax+Bu+f^{\left[2\right]}(x,u)`$
$`+f^{\left[3\right]}(x,u)+\mathrm{}+f^{\left[d\right]}(x,u)`$
$`l(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Qx+x^{}Su+{\displaystyle \frac{1}{2}}u^{}Ru`$
$`+l^{\left[3\right]}(x,u)+\mathrm{}+l^{\left[d+1\right]}(x,u)`$
where $`f𝒞^{k1}`$ and $`l𝒞^k`$.

In consequence, we expect the unknowns to have power series expansions,

$`\pi \left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Px+\pi ^{\left[3\right]}\left(x\right)+\mathrm{}`$ (2.4.15)
$`\kappa \left(x\right)`$ $`=`$ $`Kx+\kappa ^{\left[2\right]}\left(x\right)+\mathrm{}`$ (2.4.16)
as well. The known matrices are $`A,B,f^{\left[2\right]},\mathrm{}`$ and $`Q,S,R,l^{\left[3\right]},\mathrm{}`$ while $`P,\pi ^{\left[3\right]},\mathrm{}`$ and $`K,\kappa ^{\left[2\right]},\mathrm{}`$ are the unknowns.

#### 2.4.1 Special Case: Linear-Quadratic Regulator (LQR)

The linear-quadratic regulator is an infinite-time horizon optimal control problem with

$`f(x,u)`$ $`=`$ $`Ax+Bu`$ (2.4.17)
$`l(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Qx+x^{}Su+{\displaystyle \frac{1}{2}}u^{}Ru;`$ (2.4.18)
i.e., the higher degree homogeneous polynomials are zero. Also, in a linear-quadratic regulator, we expect that the optimal cost and the optimal control will be quadratic and linear, respectively; i.e.,

$`\pi \left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Px`$
$`\kappa \left(x\right)`$ $`=`$ $`Kx.`$
Thus, we look for P and K. After substituting all the expansions in (2.3.13) and collecting quadratic terms, we obtain

$`{\displaystyle \frac{1}{2}}x^{}\left[PA^{}PAK^{}B^{}PAA^{}PBKK^{}B^{}PBKQ2SKK^{}RK\right]x=0`$
Meanwhile gathering linear terms of (2.3.14), we get

$`\left(Ax+Bu\right)^{}PB+x^{}S+u^{}R=0`$ (2.4.20)
that reduces to

$`K=\left(B^{}PB+R\right)^1\left(A^{}PB+S\right)^{}.`$ (2.4.21)
It follows that (2.4.1) can be simplified to

$`{\displaystyle \frac{1}{2}}x^{}\left[PA^{}PA+\left(A^{}PB+S\right)\left(B^{}PB+R\right)^1\left(A^{}PB+S\right)^{}Q\right]x=0`$ (2.4.22)
by (2.4.21). Thus, (2.4.22) and (2.4.20) are the pair of equations obtained by collecting the quadratic terms of (2.3.13) and the linear terms of (2.3.14):

$`0`$ $`=`$ $`PA^{}PA+\left(A^{}PB+S\right)\left(B^{}PB+R\right)^1\left(A^{}PB+S\right)^{}Q`$ (2.4.23)
$`K`$ $`=`$ $`\left(B^{}PB+R\right)^1\left(A^{}PB+S\right)^{}`$ (2.4.24)
Equation (2.4.23) is known as the Discrete Algebraic Riccati Equation (DTARE). Since $`\left(B^{}PB+R\right)`$ is positive definite, the matrix $`K`$ is well-defined once P is known. The next theorem gives the necessary condition for the DTARE to have a unique positive definite solution P.

###### Theorem 2.4.1

If the pair $`(A,B)`$ is stabilizable and the pair $`(A,Q^{1/2})`$ is detectable, then there exists a unique positive definite matrix P which satisfies the Riccati equation and

$$\left|\sigma \left(A+BK\right)\right|<1.$$

where $`\sigma (A+BK)`$ is spectrum of (A+BK).

Moreover, the resulting feedback $`u=Kx`$ is asymptotically stabilizing for the system

$$x^+=\left(A+BK\right)x.$$

The above theorem can be found in .

#### 2.4.2 First Level

When the problem has nonlinear dynamics and cost with higher order terms, the first step is to look for the first terms of (2.4.15) and (2.4.16),i.e.,

$`\pi \left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Px+\mathrm{}`$
$`u\left(x\right)`$ $`=`$ $`Kx+\mathrm{}.`$
The homogeneous polynomial of higher degrees of (2.4.15) and (2.4.16) are eliminated as the only terms gathered are the quadratic terms in (2.3.13) and the linear terms in (2.3.14) when (2.4.15) and (2.4.16) are substituted in the DPE. Then, the first level is exactly the linear-quadratic regulator case. Hence, we get the solutions (2.4.23) and (2.4.24).

#### 2.4.3 Higher Levels: Higher Degree Terms

In the second level, we look for the second terms of the expansions,

$`\pi \left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Px+\pi ^{\left[3\right]}\left(x\right)+\mathrm{}`$
$`u\left(x\right)`$ $`=`$ $`Kx+\kappa ^{\left[2\right]}\left(x\right)+\mathrm{},`$
namely $`\pi ^{\left[3\right]}\left(x\right)`$ and $`\kappa ^{\left[2\right]}\left(x\right)`$. Before the method is applied, we reexpressed the $`f(x,u)`$ and $`l(x,u)`$ with the first level feedback solution; i.e.,

$`\overline{f}(x,u)`$ $`=`$ $`f(x,Kx+u)`$
$`\overline{l}(x,u)`$ $`=`$ $`l(x,Kx+u)`$
These functions have power series representation through terms of $`3^{rd}`$ and $`4^{th}`$ of the form

$`\overline{f}(x,u)`$ $`=`$ $`\left(A+BK\right)x+Bu+\overline{f}^{\left[2\right]}(x,u)+\mathrm{}`$
$`\overline{l}(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(x^{}Qx+2x^{}SKx+x^{}K^{}RKx\right)+x^{}Su+u^{}Ru+\overline{l}^{\left[3\right]}(x,u)+\mathrm{}`$
Again repeating the same process but this time grouping all the cubic terms of (2.3.13) gives

$`\pi ^{\left[3\right]}\left(x\right)\pi ^{\left[3\right]}\left(\left(A+BK\right)x\right)={\displaystyle \frac{1}{2}}x^{}\left(A+BK\right)^{}\overline{f}^{\left[2\right]}(x,0)`$
$`+{\displaystyle \frac{1}{2}}\overline{f}^{\left[2\right]^{}}(x,0)P\left(A+BK\right)x+\overline{l}^{\left[3\right]}(x,0)`$ (2.4.25)
where the terms

$`x^{}Su+{\displaystyle \frac{1}{2}}x^{}K^{}Ru+{\displaystyle \frac{1}{2}}u^{}RKx+{\displaystyle \frac{1}{2}}x^{}(A+BK)^{}PBu)+{\displaystyle \frac{1}{2}}u^{}B^{}P(A+BK)x`$ (2.4.26)
drops out of (2.4.3) as it is zero by (2.4.21) and quadratic terms of (2.3.14) yields

$`\kappa ^{\left[2\right]}\left(x\right)`$ $`=`$ $`(B^{}PB+R)^1`$
$`[{\displaystyle \frac{\overline{f}^{\left[2\right]^{}}}{u}}(x,0)P(A+BK)x+B^{}P\overline{f}^{\left[2\right]}(x,0)+B^{}{\displaystyle \frac{\pi ^{\left[3\right]}}{f}}\left((A+BK)x\right)`$
$`+{\displaystyle \frac{l^{\left[3\right]}}{x}}(x,0)]`$
Cancelling these terms (2.4.26) is crucial because it introduces a block triangular structure within the levels. If $`\kappa ^{\left[2\right]}`$ appears in (2.4.3), then the system is coupled at the second level. Similar condition applies for any $`d`$ level. The block triangular form facilitates solving the system.

The $`d^{th}`$ Level:

We discuss the system of equations at the $`d^{th}`$ level. In the $`d^{th}`$, the two linear systems are obtained by collecting the $`d+1`$ degree terms of the 2.3.13 and $`d`$ degree terms of the 2.3.14. Suppose we have solved through the $`d1^{th}`$ level, we incorporate the feedback up to the degree $`d1`$ into the dynamics and cost to obtain the updated $`\overline{f}(x,u)`$ and $`\overline{l}(x,u)`$; i.e.,

$`\overline{f}(x,u)`$ $`=`$ $`f(x,Kx+\kappa ^{\left[2\right]}\left(x\right)+\kappa ^{\left[3\right]}\left(x\right)+\mathrm{}+\kappa ^{\left[d1\right]}\left(x\right)+u)`$
$`\overline{l}(x,u)`$ $`=`$ $`l(x,Kx+\kappa ^{\left[2\right]}\left(x\right)+\kappa ^{\left[3\right]}\left(x\right)+\mathrm{}+\kappa ^{\left[d1\right]}\left(x\right)+u).`$
In the $`d^{th}`$ level, the corresponding cancellation is

$`x^{}Su+{\displaystyle \frac{1}{2}}x^{}K^{}Ru+{\displaystyle \frac{1}{2}}u^{}RKx+{\displaystyle \frac{1}{2}}x^{}\left(A+BK\right)^{}PBu+{\displaystyle \frac{1}{2}}u^{}B^{}P\left(A+BK\right)x=0`$ (2.4.28)
After the substitution and collection of all the $`d+1`$ degree terms of (2.3.13) and $`d`$ degree terms of (2.3.14), we have the equations

$`\pi ^{\left[d+1\right]}\left(x\right)\pi ^{\left[d+1\right]}\left(\left(A+BK\right)x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[x^{}(A+BK)^{}\overline{f}^{\left[d\right]}(x,0)`$
$`+{\displaystyle \underset{j=2}{\overset{d1}{}}}\left[\kappa ^{\left[j\right]^{}}B^{}+\overline{f}^{\left[j\right]^{}}(x,0)\right]P\left[B\kappa ^{\left[d1j\right]}+\overline{f}^{\left[d1j\right]}(x,0)\right]`$
$`+\overline{f}^{\left[d\right]^{}}(x,0)P(A+BK)x]+\stackrel{~}{\pi }^{\left[d+1\right]}(\overline{f})+\overline{l}^{\left[d+1\right]}(x,0)`$
and

$`\kappa ^{\left[d\right]}\left(x\right)`$ $`=`$ $`(B^{}PB+R)^1`$
$`[{\displaystyle \frac{f^{\left[d\right]^{}}}{u}}(x,0)P(A+BK)x`$
$`+{\displaystyle \underset{j=2}{\overset{d}{}}}{\displaystyle \frac{\overline{f}^{\left[j\right]}}{u}}(PB\kappa ^{\left[dj\right]}+{\displaystyle \frac{\pi ^{\left[d+1j\right]}}{f}}+{\displaystyle \frac{\overline{l}^{\left[d+1\right]}}{u}}(x,0)]`$
where $`\stackrel{~}{\pi }^{\left[d+1\right]}\left(\overline{f}\right)`$ are terms of degree $`d+1`$ of the homogeneous polynomials $`\pi ^{\left[m\right]}`$ where $`2md+1`$. These terms are generated by the input of the feedback control into the DPE.

#### 2.4.4 Solvability of the System of Equations

We know discuss the solutions of the system of equations obtained from level $`2`$. For the case $`d=2`$, we show that $`\pi ^{\left[3\right]}`$ exists as a solution for (2.4.3). First, the equation (2.4.3) is linear since

$`\pi ^{\left[3\right]}\left(x\right)\pi ^{\left[3\right]}\left(x\right)+\pi ^{\left[3\right]}\left(\left(A+BK\right)x\right)`$ (2.4.31)
is linear

###### Lemma 2.4.2

Given that $`|\sigma (A+BK)|<1`$. There exists a unique solution of the homogeneous polynomial $`\pi ^{[3]}`$.

Proof:  For simplicity, we assume that $`\left(A+BK\right)`$ has simple eigenvalues $`\mu _i`$, i.e.

$$w_i\left(A+BK\right)=\mu _iw_i.$$

where $`w_i^{1\times n}`$. Since the system is stable,

$`\left|\mu _i\right|<1.`$ (2.4.32)
The polynomial $`\pi ^{\left[3\right]}\left(x\right)`$ has the following form:

$$\underset{i_1}{\overset{n}{}}\underset{i_2}{\overset{n}{}}\underset{i_3}{\overset{n}{}}c_{i_1,i_2,i_3}\left(w_{i_1}x_{i_1}\right)\left(w_{i_2}x_{i_2}\right)\left(w_{i_3}x_{i_3}\right).$$

where $`\left(w_{i_1}x_{i_1}\right)\left(w_{i_2}x_{i_2}\right)\left(w_{i_3}x_{i_3}\right)`$ are the basis for all cubic polynomials. It follows that

$`\pi ^{\left[3\right]}\left(x\right)\pi ^{\left[3\right]}\left(\left(A+BK\right)x\right)=`$
$`{\displaystyle \underset{i_1}{\overset{n}{}}}{\displaystyle \underset{i_2}{\overset{n}{}}}{\displaystyle \underset{i_3}{\overset{n}{}}}c_{i_1,i_2,i_3}\left(1\mu _{i_1}\mu _{i_2}\mu _{i_3}\right)\left(w_{i_1}x_{i_1}\right)\left(w_{i_2}x_{i_2}\right)\left(w_{i_3}x_{i_3}\right)`$
The RHS of (2.4.3) is

$`{\displaystyle \underset{i_1}{\overset{n}{}}}{\displaystyle \underset{i_2}{\overset{n}{}}}{\displaystyle \underset{i_3}{\overset{n}{}}}d_{i_1,i_2,i_3}\left(w_{i_1}x_{i_1}\right)\left(w_{i_2}x_{i_2}\right)\left(w_{i_3}x_{i_3}\right).`$ (2.4.34)
Together with (2.4.4) and (2.4.34), we express the coefficients as

$$c_{i_1,i_2,i_3}=\frac{d_{i_1,i_2,i_3}}{\left(1\mu _{i_1}\mu _{i_2}\mu _{i_3}\right)}$$

where

$$\left(1\mu _{i_1}\mu _{i_2}\mu _{i_3}\right)0$$

by (2.4.32). Hence, there exists a unique homogeneous polynomial $`\pi ^{\left[3\right]}\left(x\right)`$.   It follows that the linear map (2.4.31) is invertible. For any level d, the homogeneous polynomial $`\pi ^{\left[d\right]}\left(x\right)`$ exists. The proof on the existence for higher order levels follows in the same fashion as above.

For any $`d`$ level, the RHS of equation (2.4.3) contains a multiplicative term

$$\left(B^{}PB+R\right)^1.$$

The inverse matrix exists because

$$\left(B^{}PB+R\right)$$

is positive definite since $`R`$ is positive definite. The rest of the terms on the RHS of (2.4.3) are known from the previous levels. Thus, $`\kappa ^{\left[d\right]}`$ exists for any level $`d`$. Therefore, all of the linear equations for degree 3 or higher are solvable.

## Chapter 3 Existence of the Local Solution of the DPE

In the previous chapter, we found polynomials of degree $`\left(r1\right)`$ and $`r`$, namely the optimal control and optimal cost, respectively, that satisfy the DPE. We will prove the existence of smooth solutions to the DPE which have the same Taylor series expansions as of the formal solutions previously discussed in Chapter $`3`$.

### 3.1 The Main Result

Here we state our main result.

###### Theorem 3.1.1

Suppose that the dynamics and cost are $`C^{r1}`$ and $`C^r`$, respectively and the linear part of the nonlinear Hamiltonian system (3.4.1) is stabilizable and detectable around zero. Then there exists $`C^{r1}(^n)`$ optimal control and $`C^r(^n)`$ optimal cost solutions to the DPE locally around zero.

If we let $`k\mathrm{}`$ in Theorem (3.1.1), we then have the following corollary.

###### Corollary 3.1.2

Suppose that the dynamics and cost are real analytic functions on $`^{2n}`$ and the linear part of the nonlinear Hamiltonian system (3.4.1) is stabilizable and detectable around zero. Then the power series solutions converges to solution of the DPE locally around zero.

### 3.2 PMP and Hamiltonian Dynamics

We look at the Pontryagin Maximum Principle (PMP) as it presents the Hamilton difference equations satisfied by the optimal trajectories. The PMP gives a necessary condition for a control to be optimal. Along with the optimal control problem formulation, there is an associated nonlinear Hamiltonian,

$`H(x,u,\lambda ^+)=\lambda ^+^{}f(x,u)+l(x,u)`$ (3.2.1)
where $`\lambda ^+=\lambda _{k+1}`$. The PMP states the following:

###### Theorem 3.2.1

If $`x_k`$ and $`u_k`$ are optimal for $`k0,1,2,\mathrm{}`$, then there exists $`\lambda _k0`$ for $`k0,1,2,\mathrm{}`$ such that

$`x^+`$ $`=`$ $`{\displaystyle \frac{H}{\lambda ^+}}(x,u,\lambda ^+)`$ (3.2.2)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{H}{x}}(x,u,\lambda ^+)`$ (3.2.3)
$`u^{}`$ $`=`$ $`\text{arg}\underset{u}{\mathrm{min}}H(x,u,\lambda ^+).`$ (3.2.4)
Thus, the minimizer of the nonlinear Hamiltonian evaluated at the optimal $`x`$ and $`\lambda ^+`$ is the optimal control $`u`$ amongst all admissible controls $`v`$. Note that $`u^{}(x,\lambda ^+)𝒞^{r1}`$ since $`f𝒞^{r1}`$ and $`l𝒞^r`$ in (3.2.1). We assume that $`H`$ is convex in u to guarantee a unique optimal control. So the PMP gives the existence of the optimal control. However, the equation (3.2.4) is an optimal control that has $`x`$ and $`\lambda `$ as its independent variables. We must show that $`\lambda `$ is function of $`x`$ to complete this part of the proof.

### 3.3 Forward Hamiltonian Dynamics

If we take the (3.2.2), (3.2.3), and (3.2.4) of the PMP and linearize the equations around zero. The linearized sytem that we obtain is exactly the linear Hamiltonian system. The corresponding Hamiltonian is

$$H(x,\lambda ^+,u)=\lambda ^+^{}\left(Ax+Bu\right)+\frac{1}{2}x^{^{}}Qx+x^{^{}}Su+\frac{1}{2}u^{^{}}Ru$$

and the system is

$`\left[\begin{array}{c}x^+\\ \lambda \end{array}\right]=\left[\begin{array}{c}x\\ \lambda ^+\end{array}\right]`$ (3.3.9)
where

$`=\left[\begin{array}{cc}ABR^1S^{}& BR^1B^{}\\ QSR^1S^{}& A^{}SR^1B^{}\end{array}\right]`$
is the associated Hamiltonian matrix.

Notice the opposing directions of the propagation of the state and costate dynamics in (3.3.9). This presents difficulty in studying some Hamiltonian properties, but it is addressed in Chapter 4. By assuming invertibility of the matrix $`ASR^1B^{}`$, we can rewrite the Hamilton equations (3.3.9) so that both equations propagate in the same direction, in particular the direction of increasing time. The forward recursion,

$`\left[\begin{array}{c}x^+\\ \lambda ^+\end{array}\right]=^F\left[\begin{array}{c}x\\ \lambda \end{array}\right]`$ (3.3.15)
where

$`^F=\left[\begin{array}{cc}\left(ABR^1S^{}\right)BR^1B^{}\left(A^{}SR^1B^{}\right)^1\left(QSR^1S^{}\right)& BR^1B^{}\left(A^{}SR^1B^{}\right)^1\\ \left(A^{}SR^1B^{}\right)^1\left(QSR^1S^{}\right)& \left(A^{}SR^1B^{}\right)^1\end{array}\right]`$ (3.3.18)
is directly derived from the system (3.3.9). The forward recursion describes the flow of the state and costate as $`t`$ approaches to infinity. We refer to $`^F`$ as the forward Hamiltonian matrix. Thus, the dynamics (3.4.41) rewritten in the forward direction relies on the invertibility of $`A^{}SR^1B^{}`$.

The existence of $`\left(A^{}SR^1B^{}\right)^1`$ has the following consequences:

###### Lemma 3.3.1

If $`A^{}SR^1B^{}`$ is invertible, then $`^F`$ exists and is invertible.

Proof:  First, we observe in (3.3.18) that $`\left(A^{}SR^1B^{}\right)^1`$ appears. Then, if $`\left(A^{}SR^1B^{}\right)^1`$ exists so does $`^F`$. Explicit calculation gives

$`\left(^F\right)^1=\left[\begin{array}{cc}\left(ABR^1S^{}\right)^1& \left(ABR^1S^{}\right)BR^1B^{}\\ \left(QSR^1S^{}\right)\left(ABR^1S^{}\right)^1& \left(QSR^1S^{}\right)\left(ABR^1S^{}\right)BR^1B^{}+\left(ABR^1S^{}\right)^{}\end{array}\right]`$ (3.3.21)
The matrix (3.3.21) satisfies the backward hamiltonian dynamics. Thus, $`^B=\left(^F\right)^1`$ where

$`\left[\begin{array}{c}x\\ \lambda \end{array}\right]=^B\left[\begin{array}{c}x^+\\ \lambda ^+\end{array}\right].`$ (3.3.26)
We see that $`^B`$ exists since $`\left(ABR^1S^{}\right)`$ is invertible. Hence, $`^F`$ is invertible.   It follows that $`^F`$ does not have zero as an eigenvalue if $`A^{}SR^1B^{}`$ is invertible. Also, the Lemma (3.3.1) implies that if $`^F`$ has a zero eigenvalue, then so does the spectrum of $`ASR^1B^{}`$. In , if $`(A,B)`$ is stabilizable and $`(Q^{1/2},A)`$ is detectable, then the system is hyperbolic; i.e., none of the eigenvalues lie on the unit circle. If zero is an eigenvalue, then infinity is also an eigenvalue because of the hyperbolicity of the system. We refer to these eigenvalues as infinite eigenvalues. This issue of infinite eigenvalues is discussed in Chapter 4. For the time being we will require that the forward Hamiltonian matrix exists and be invertible. For the rest of the chapter, we assume the invertibility of $`A^{}SR^1B^{}`$ since the existence of $`^F`$ depends on the inverse of $`A^{}SR^1B^{}`$.

The forward Hamiltonian matrix has $`2n`$ eigenvalues; n are stable, i.e., these eigenvalues lie inside the unit circle. The other $`n`$ eigenvalues are unstable and are positioned outside the unit circle. We will see in Section 3.4.3 that the eigenstructure is hyperbolic.

Recall the LQR problem in Chapter $`2`$ of minimizing

$`\underset{u}{\mathrm{min}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{2}}x_j^{}Qx_j+x_j^{}Su_j+{\displaystyle \frac{1}{2}}u_j^{}Ru_j`$
subject to the dynamics (3.3.27)
$`x^+=Ax+Bu`$
$`x_0=x\left(0\right).`$
###### Lemma 3.3.2

Suppose $`u=Kx`$ is the optimal control. If zero is an eigenvalue of $`(ABR^1S^{})`$, then zero is in $`\sigma (A+BK)`$.

Proof:  With the transformation $`u=Lx+v`$ the optimal control problem is changed to the following problem:

$`\underset{u}{\mathrm{min}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{2}}x_j^{}\left(Q+L^{}RL+S\right)x_j+{\displaystyle \frac{1}{2}}v_j^{}Rv_j`$
subject to the dynamics (3.3.28)
$`x^+=\left(ABR^1S^{}\right)x+Bv`$
$`x_0=x\left(0\right)`$
where $`L=R^1S^{}`$. The transformation removes the cross term in the cost. For the modified problem (3.3), we have $`v=\overline{K}x`$ is the optimal control, then $`K=L+\overline{K}`$. The transformation removes the cross term in the cost. Then, we choose the eigenvector $`x_00`$ of $`ABR^1S^{}`$ such that its corresponding eigenvalue is zero; i.e.,

$`\left(ABR^1S^{}\right)x_0=0`$ (3.3.29)
The cost starting from $`x_0`$ is

$`\pi \left(x_0\right)=x_0^{}\left(Q+L^{}RL+S\right)x_0.`$ (3.3.30)
No additional control $`v`$ is exercised since the dynamics reaches zero at one time step. Therefore, $`v=\overline{K}x_0=0`$ and the optimal cost is the cost (3.3.30) at the zero time. We have that

$`Lx_0=Kx_0`$ (3.3.31)
for $`x_00`$ that satisfies (3.3.29) since $`Kx_0=\left(L+\overline{K}\right)x_0`$.

Returning to (3.3), the closed-looped spectrum is $`\sigma \left(A+BK\right)`$. Then, from (3.3.31)

$`\left(A+BL\right)x_0=\left(A+BK\right)x_0.`$ (3.3.32)
Since the LHS of (3.3.32) is equal to 0, we have

$$\left(A+BK\right)x_0=0$$

for $`x_00`$. Thus, $`0\sigma \left(A+BK\right)`$. Hence, zero is a closed-loop eigenvalue.   By Lemma (3.3.2), zero is not a closed-loop eigenvalue. This implies that zero is not an eigenvalue of $`ABR^1S^{}`$ if $`u=Kx`$ is assumed to be the optimal control.

### 3.4 Properties of the Nonlinear Hamiltonian Dynamics

It is imperative to look at some of the properties that the nonlinear Hamiltonian dynamics possesses as it is needed in the proof of the existence of the optimal cost. We look at its tangent dynamics, sympletic form, and eigenstructure.

Recall that the PMP gives the nonlinear Hamilton difference equations (3.2.2), (3.2.3) satisfied by the optimal trajectories. These equations are the nonlinear Hamiltonian dynamics corresponding to the optimal control problem. We calculate the forward nonlinear Hamiltonian dynamics from (3.2.2), (3.2.3). Since the forward linear Hamiltonian dynamics exists at zero and by the Implict Function Theorem, the forward nonlinear dynamics also exists in the neigborhood around 0. The forward nonlinear Hamiltonian dynamics,

$`x^+`$ $`=`$ $`G_1(x,\lambda )`$ (3.4.33)
$`\lambda ^+`$ $`=`$ $`G_2(x,\lambda )`$ (3.4.34)
evolves in the direction of increasing time.

#### 3.4.1 Tangent Dynamics

We now introduce the idea of tangent vectors to $`=\left\{(x,\lambda )\right|x^n,\lambda +^n)\}`$; i.e., $`=^{2n}`$. First, we denote the tangent vectors

$`v=\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right].`$
The set of tangent vector to $``$ at $`(x,\lambda )`$ forms a vector space $`T_{x,\lambda }`$, the tangent space to $``$ at $`(x,\lambda )`$. In addition, the tangent bundle of $``$, denoted by $`T`$, is a differential manifold where

$`T={\displaystyle \underset{x,\lambda ^+}{}}T_{x,\lambda ^+}.`$ (3.4.36)
The local coordinate system on $`T`$ is $`4n`$ numbers $`x_1,\mathrm{},x_n,\lambda _1,\mathrm{},\lambda _n`$ and
$`\delta x_1,\mathrm{},\delta x_n,\delta \lambda _1,\mathrm{},\delta \lambda _n`$ where the former is the local coordinate on $``$ and the latter is the components of the tangent vector.

The tangent dynamics describes the flow of the tangent vector at each sequential point of a map.

###### Definition 3.4.1

Suppose $`x_{k+1}=f(x_k)`$ with known $`x_0`$. Let $`\delta x_0`$ be a tangent vector at $`x_0`$. Then we define the tangent dynamics around the trajectory $`x_k`$ as

$$\delta x_{k+1}=\frac{f}{x}\left(x_k\right)\delta x_k.$$

When the dynamics is linear; i.e., $`x_{k+1}=Ax_k`$, then the tangent dynamics is $`v_{k+1}=Av_k`$ because $`\frac{f}{x}\left(x_k\right)=A`$. Thus, the linear Hamiltonian dynamics (3.4.41) has

$`\left[\begin{array}{c}\delta x^+\\ \delta \lambda ^+\end{array}\right]=^F\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]`$ (3.4.41)
as it tangent dynamics. We now find the nonlinear tangent dynamics as we have a nonlinear system. Note the increase in difficulty in finding the tangent dynamics of the forward nonlinear Hamilton dynamics because of the mixed directions of the propagation. Linearizing the the nonlinear system at the trajectories $`(x,\lambda ^+)`$ invokes a linear perturbation. We perturb $`x`$ and $`\lambda `$ and perform the first variation on the

$$\frac{H}{(x,\lambda ^+)}(x,\lambda ^+).$$

The nonlinear Hamiltonian dynamics from PMP is

$`x^+`$ $`=`$ $`{\displaystyle \frac{H}{\lambda ^+}}(x,\lambda ^+)`$
$`\lambda `$ $`=`$ $`{\displaystyle \frac{H}{x}}(x,\lambda ^+).`$ (3.4.42)
Suppose we replace the parameters

$`x`$ $``$ $`x+\delta x`$
$`\lambda ^+`$ $``$ $`\lambda ^++\delta \lambda ^+.`$
Then,

$`x^++\delta x^+`$ $``$ $`{\displaystyle \frac{H}{\lambda ^+}}(x+\delta x,\lambda ^++\delta \lambda ^+)`$
$`\lambda +\delta \lambda `$ $``$ $`{\displaystyle \frac{H}{x}}(x+\delta x,\lambda ^++\delta \lambda ^+).`$ (3.4.43)
Subtracting (3.4.1) from (3.4.1), we have

$`\delta x^+`$ $``$ $`{\displaystyle \frac{H}{\lambda ^+}}(x+\delta x,\lambda ^++\delta \lambda ^+){\displaystyle \frac{H}{\lambda ^+}}(x,\lambda ^+)`$ (3.4.44)
$`\delta \lambda `$ $``$ $`{\displaystyle \frac{H}{x}}(x+\delta x,\lambda ^++\delta \lambda ^+){\displaystyle \frac{H}{x}}(x,\lambda ^+).`$ (3.4.45)
Expanding the RHS of (3.4.44) and (3.4.45) in Taylor series, we get

$`\delta x^+`$ $`=`$ $`{\displaystyle \frac{^2H}{\lambda ^+x}}(x,\lambda ^+)\delta x+{\displaystyle \frac{^2H}{^2\lambda ^+}}(x,\lambda ^+)\delta \lambda ^+`$
$`\delta \lambda `$ $`=`$ $`{\displaystyle \frac{^2H}{^2x}}(x,\lambda ^+)\delta x+{\displaystyle \frac{^2H}{\lambda ^+x}}(x,\lambda ^+)\delta \lambda ^+.`$ (3.4.46)
Equivalently,

$`\left[\begin{array}{c}\delta x^+\\ \delta \lambda \end{array}\right]=_{\delta ,k}(x,\lambda ^+)\left[\begin{array}{c}\delta x\\ \delta \lambda ^+\end{array}\right].`$ (3.4.51)
where

$`_{\delta ,k}(x,\lambda ^+)=\left[\begin{array}{cc}H_{\lambda ^+x}& H_{\lambda ^+\lambda ^+}\\ H_{xx}& H_{x\lambda ^+}\end{array}\right](x,\lambda ^+).`$
The subscripts of $`H_{x\lambda ^+}`$ denote partial derivatives. We look at the forward perturbed system. The tangent dynamics in the forward time is

$`\left[\begin{array}{c}\delta x^+\\ \delta \lambda ^+\end{array}\right]=_{\delta ,k}^F(x,\lambda ^+)\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]`$ (3.4.57)
where

$`_{\delta ,k}^F(x,\lambda ^+)=\left[\begin{array}{cc}H_{\lambda ^+x}H_{\lambda ^+\lambda ^+}H_{\lambda ^+\lambda ^+}^1H_{xx}& H_{\lambda ^+\lambda ^+}H_{x\lambda }^1\\ H_{x\lambda ^+}^1H_{xx}& H_{x\lambda ^+}^1\end{array}\right](x,\lambda ^+).`$
To have the tangent dynamics written in forward time, we assume the invertibility of $`H_{\lambda ^+x}`$. When $`H_{\lambda ^+x}(x,\lambda ^+)`$ is evaluated at $`\stackrel{}{0}`$, we get the matrix $`A^{}SR^1B^{}`$ which we have assumed to be invertible. Hence, $`H_{\lambda ^+x}`$ is invertible for small $`x`$ and $`\lambda ^+`$.

#### 3.4.2 The Standard Symplectic Form

We define a nondegenerate and a bilinear symplectic two-form $`\mathrm{\Omega }:T_{(x,\lambda )}\times T_{(x,\lambda )}`$,

$`\mathrm{\Omega }(v,w)=v^{}Jw\text{ and }J=\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right]`$
where

$$\mathrm{\Omega }(v,w)=\mathrm{\Omega }(w,v)$$

and

$$v=\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right],w=\left[\begin{array}{c}\stackrel{~}{\delta x}\\ \stackrel{~}{\delta \lambda }\end{array}\right]$$

$`(x,\lambda ^+)`$ and $`(v,w)T_{(x,\lambda )}`$. The matrix $`J`$ is called the symplectic matrix.

With the system (3.4.57), it can be shown through calculations that

$$_{\delta ,k}^F^{}J_{\delta ,k}^F=J.$$

Then, we have that

$`\left[\begin{array}{cc}H_{\lambda ^+x}^TH_{xx}^TH_{x\lambda ^+}^1H_{\lambda ^+\lambda ^+}^T& H_{xx}^TH_{x\lambda ^+}^T\\ H_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}^T& H_{x\lambda ^+}^T\end{array}\right]J\left[\begin{array}{cc}H_{\lambda ^+x}H_{\lambda ^+\lambda ^+}H_{\lambda ^+\lambda ^+}^1H_{xx}& H_{\lambda ^+\lambda ^+}H_{x\lambda }^1\\ H_{x\lambda ^+}^1H_{xx}& H_{x\lambda ^+}^1\end{array}\right](x,\lambda ^+)`$
is equal to

$`\left[\begin{array}{cc}𝔸_{11}& 𝔸_{12}\\ 𝔸_{21}& 𝔸_{22}\end{array}\right]`$ (3.4.63)
where

$`𝔸_{11}`$ $`=`$ $`H_{xx}^TH_{x\lambda ^+}^TH_{\lambda ^+x}H_{xx}^TH_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}H_{x\lambda ^+}^1H_{xx}H_{\lambda ^+x}^TH_{x\lambda ^+}^1H_{xx}+H_{xx}^TH_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}^TH_{x\lambda ^+}^1H_{xx}`$
$`𝔸_{12}`$ $`=`$ $`H_{xx}^TH_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}H_{x\lambda ^+}^1+H_{\lambda ^+x}^TH_{x\lambda ^+}^1H_{xx}^TH_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}^TH_{x\lambda ^+}^1`$
$`𝔸_{21}`$ $`=`$ $`H_{x\lambda ^+}^TH_{\lambda ^+x}+H_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}H_{x\lambda ^+}^1H_{xx}H_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}H_{x\lambda ^+}^1H_{xx}`$
$`𝔸_{22}`$ $`=`$ $`H_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}H_{x\lambda ^+}^1+H_{x\lambda ^+}^TH_{\lambda ^+\lambda ^+}^TH_{x\lambda ^+}^1.`$
Since the matrices $`H_{\lambda ^+x}=H_{x\lambda ^+},H_{xx},\text{and}H_{\lambda ^+\lambda ^+}`$ are symmetric, then (3.4.63) reduces to $`J`$. Thus,

$`_{\delta ,k}^F^{}J_{\delta ,k}^F=J`$ (3.4.64)
Hence, the nonlinear Hamilton dynamics preserves the standard symplectic form. Suppose that $`(v,w)`$ are two tangent vectors that satisfy (3.4.57); i.e., $`v^+=_{\delta ,k}^Fv`$ and $`w^+=_{\delta ,k}^Fw`$. Let $`v`$, $`w`$ be two sequences propagating according to(3.4.57). Then,

$`\mathrm{\Omega }(v^+,w^+)`$ $`=`$ $`v^{}_{}{}^{}+Jw^+`$
$`=`$ $`v^{}_{\delta ,k}^F^{}J_{\delta ,k}^Fw`$
$`=`$ $`v^{}Jw`$
$`=`$ $`\mathrm{\Omega }(v,w).`$
In other words, for any two sequences of tangent vector under the tangent dynamics (3.4.51), the value of the $`\mathrm{\Omega }`$ is unchanged at every point of the map. See fig.(3.1).

As a consequence, if we evaluate (3.4.64) at 0, then at the linear level

$`^F^{}J^F=J.`$ (3.4.65)
where $`^F^{}`$ is the matrix in (3.4.41). The calculations are as follows: For simplicity we let

$`\alpha `$ $`=`$ $`ABR^1S^{}`$
$`\beta `$ $`=`$ $`BR^1B^{}`$
$`\gamma `$ $`=`$ $`QSR^1S^{}.`$
Then,

$`^F=\left[\begin{array}{cc}\alpha \beta \alpha ^T\gamma & \beta \alpha ^T\\ \alpha ^T\gamma & \alpha ^T\end{array}\right]`$
and

$`^F^{}=\left[\begin{array}{cc}\alpha ^T\gamma \alpha ^1\beta & \gamma \alpha ^T\\ \alpha ^1\beta & \alpha ^1\end{array}\right].`$
We denote $`\alpha ^T`$ as the inverse of the transpose of $`\alpha `$. Also, the matrices $`\beta `$ and $`\gamma `$ are symmetric. The equation (3.4.65) is

$`\left[\begin{array}{cc}\alpha ^T\gamma \alpha ^T\beta & \gamma \alpha ^T\\ \alpha ^1\beta & \alpha ^1\end{array}\right]\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right]\left[\begin{array}{cc}\alpha \beta \alpha ^T\gamma & \beta \alpha ^T\\ \alpha ^T\gamma & \alpha ^T\end{array}\right]`$
which is simplified to

$`\left[\begin{array}{cc}\gamma \gamma \alpha ^1\beta \alpha ^T\gamma \gamma +\gamma \alpha ^1\beta \alpha ^T\gamma & I+\gamma \alpha ^1\beta \alpha ^T\gamma \alpha ^1\beta \alpha ^T\\ I+\alpha ^1\beta \alpha ^T\gamma \alpha ^1\beta \alpha ^T\gamma & \alpha ^1\beta \alpha ^T+\alpha ^1\beta \alpha ^T\end{array}\right]=J.`$
Therefore,

$$^F^{}J^F=J$$

holds.

Similarly, if we take tangent vectors satisfying the linear tangent dynamics

$$v^+=^Fv\text{ and }w^+=^Fw,$$

then

$`\mathrm{\Omega }(v^+,w^+)=\mathrm{\Omega }(v,w).`$
Thus, the two-form $`\mathrm{\Omega }(v,w)`$ is unchanged when evaluated at the tangent vectors $`v`$ and $`w`$ of every point of the map under the linear Hamiltonian dynamics.

#### 3.4.3 Eigenstructure

Now with (3.4.65), the eigenstructure of $`^F`$ is as follow:

###### Theorem 3.4.2

Suppose $`𝒜^{^{}}J𝒜=J`$. If $`\mu \sigma (𝒜)`$, then so do $`\frac{1}{\mu },\overline{\mu },`$ and $`\frac{1}{\overline{\mu }}\sigma (𝒜)`$.

Proof:  Suppose $`(\delta x,\delta \lambda )^{}`$ is an eigenvector of $`𝒜`$ for some eigenvalue $`\mu `$.

$`𝒜\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]=\mu \left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]`$
Also, note that

$`J\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]=\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right].`$
Then,

$`𝒜{}_{}{}^{}J𝒜\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]=J\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]`$ $``$ $`\mu 𝒜^{^{}}J\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right]=\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right]`$
$``$ $`\mu 𝒜^{^{}}\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right]=\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right]`$
$``$ $`𝒜^{^{}}\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right]={\displaystyle \frac{1}{\mu }}\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right].`$
Since $`\mu 0`$, we have

$$(\frac{1}{\mu },\left[\begin{array}{c}\delta \lambda \\ \delta x\end{array}\right])$$

as the eigenpair of $`𝒜`$. Since $`\sigma \left(𝒜^{^{}}\right)=\sigma \left(𝒜\right)`$, the eigenvalue $`\frac{1}{\mu }`$ is also $`\sigma \left(𝒜\right)`$. It follows that the complex conjugates $`\overline{\mu }`$ and $`\frac{1}{\overline{\mu }}\sigma \left(𝒜\right)`$ as the characteristic polynomial $`p\left(\mu \right)=p\left(\overline{\mu }\right)`$ implies that $`p\left(\overline{\mu }\right)=0`$. Thus, the eigenvalues of $`𝒜`$ come in reciprocal pairs and complex conjugate pairs.

Since $`^F^{}J^F`$, then the theorem above verifies the claim of the partitioning of $`2n`$ eigenvalues in the complex plane. None of the eigenvalues of $`^F`$ are on the unit circle since the system is linearly stabilizable and detectable.

Moreover, the eigenvectors corresponding to the stable eigenvalues sitting inside the unit circle spans subspace $`E_s`$. Similarly, the subspace $`E_u`$ is spanned by the unstable eigenvectors whose corresponding eigenvectors lie outside the unit circle. In Section 3.6, we show that the subspace $`E_s`$ is invariant.

### 3.5 Local Stable Manifold

#### 3.5.1 Local Stable Manifold Theorem

###### Theorem 3.5.1

Given the dynamics

$$x_{k+1}=G\left(x_k\right)$$

$$G\left(0\right)=0$$

Let $`G:U^n`$ be a $`𝒞^{r1}(^{2n})`$ map with a hyperbolic fixed point $`0`$. Then there is a local stable manifold, $`W^s(0)𝒞^{r1}(^{2n})`$, that is tangent to the eigenspace $`E_0^s`$ of the Jacobian of $`G`$ at $`0`$. Define

$$W^s\left(0\right)=\left\{xU\right|\underset{t\mathrm{}}{lim}G^t\left(x\right)=0\}$$

and
$`E_0^s=\{`$span of eigenvectors whose corresponding eigenvalues are such that $`|\lambda |<1\}`$.
$`W^s(0)`$ is a smooth manifold.

The statement of this theorem is found in . The following theorem was proven by Hartman , but a more modern technique can be found in .

Along with our dynamics (3.4.33) and a critical point $`0`$, the local stable manifold theorem gives the existence of a local stable manifold $`W^s\left(0\right)`$. The tangency of the linear subspace $`E_s`$ to $`W_s`$ implies that the lowest degree term of the local stable manifold is a quadratic homegeneous polynomial.

#### 3.5.2 Construction of the Local Stable Manifold

Through the Taylor approximation technique, we explicitly construct the local stable manifold term by term. These calculations coincide with the power series solutions found in Chapter 2. Having a system that has $`n`$ finite eigenvalues outside and $`n`$ finite eigenvalues inside the unit circle, we can find a linear transformation that block diagonalizes our system (3.4.41). We make a linear change of coordinate such that the dynamics (3.4.41) is block diagonalized,

$`\left[\begin{array}{c}z_s\\ z_u\end{array}\right]=T\left[\begin{array}{c}x\\ \lambda \end{array}\right].`$
Matrix $`T`$ transforms the dynamics (3.4.41) into

$`\left[\begin{array}{c}z_s^+\\ z_u^+\end{array}\right]=\left[\begin{array}{cc}A_s& 0\\ 0& A_u\end{array}\right]\left[\begin{array}{c}z_s\\ z_u\end{array}\right]+\left[\begin{array}{c}f_s(z_s,z_u)\\ f_u(z_s,z_u)\end{array}\right]`$ (3.5.84)
where $`A_s`$ is the stable matrix, $`A_u`$ is the unstable matrix, and

$`f_s(z_s,z_u)`$ $`=`$ $`f_s^{\left[2\right]}(z_s,z_u)+f_s^{\left[3\right]}(z_s,z_u)+\mathrm{}`$
$`f_u(z_s,z_u)`$ $`=`$ $`f_u^{\left[2\right]}(z_s,z_u)+f_u^{\left[3\right]}(z_s,z_u)+\mathrm{}`$
where $`f_s^{\left[d\right]}(z_s,z_u)`$ and $`f_u^{\left[d\right]}(z_s,z_u)`$ are nonlinear stable and unstable homogeneous polynomials of degree $`d`$. We look for the Taylor expansion of local stable manifold $`z_u=\varphi \left(z_s\right)𝒞^{r1}\left(^n\right)`$; thus, we seek the form

$`\varphi \left(z_s\right)=\varphi ^{\left[2\right]}\left(z_s\right)+\varphi ^{\left[3\right]}\left(z_s\right)+\mathrm{}+\varphi ^{\left[r1\right]}\left(z_s\right).`$ (3.5.85)
The linear term $`\varphi ^{\left[1\right]}\left(z_s\right)`$ is zero as stated in the local stable manifold theorem. The graph $`z_u=\varphi \left(z_s\right)`$ is said to be an *invariant manifold* if

$`z_u^+=\varphi \left(z_s^+\right)\text{ whenever }z_u=\varphi \left(z_s\right);`$ (3.5.86)
i.e., the solution of (3.5.84) lies in $`z_u=\varphi \left(z_s\right)`$ for all time. By invariance (3.5.86), we have

$`z_u^+=\varphi \left(z_s^+\right)`$ $``$ $`A_uz_u+f_u(z_s,z_u)=\varphi \left(A_sz_s+f_s(z_s,z_u)\right)`$
$``$ $`A_u\varphi \left(z_u\right)+f_u(z_s,\varphi \left(z_s\right))=\varphi \left(A_sz_s+f_s(z_s,\varphi \left(z_s\right))\right).`$
Thus,

$`\varphi \left(A_sz_s+f_s(z_s,\varphi \left(z_s\right))\right)=A_u\varphi \left(z_s\right)+f_u(z_s,\varphi \left(z_s\right)).`$ (3.5.87)
Taking all the second degree terms of (3.5.87), we see

$`A_u\varphi ^{\left[2\right]}\left(z_s\right)\varphi ^{\left[2\right]}\left(A_sz_s\right)=f_u^{\left[2\right]}\left(z_s\right).`$ (3.5.88)
The map

$`\varphi ^{\left[2\right]}\left(z_s\right)A_u\varphi ^{\left[2\right]}\left(z_s\right)\varphi ^{\left[2\right]}\left(A_sz_s\right)`$ (3.5.89)
is linear. Now, we question the invertibility of the map. The quadratic term $`\varphi ^{\left[2\right]}\left(z_s\right)`$ of (3.5.88) as the spectrum of the map is

$$\xi \mu _i\mu _j0$$

where

$$A_uv_k=\xi v_k\text{ and }w_iA_s=\mu _iw_i$$

and

$`\varphi ^{\left[2\right]}\left(z_s\right)={\displaystyle \underset{i}{\overset{n}{}}}{\displaystyle \underset{j}{\overset{n}{}}}c_{s_i,s_j}v_k\left(w_iz_s\right)\left(w_jz_s\right).`$ (3.5.90)
The details of the proof mirror the steps in Lemma (2.4.2). Moreover, the polynomial $`\varphi ^{\left[d\right]}\left(x\right)`$ for all $`2dk`$ also solve the following form:

$`A_u\varphi ^{\left[d\right]}\left(z_s\right)\varphi ^{\left[d\right]}\left(A_sz_s\right)=f_u^{\left[d\right]}\left(z_s\right).`$ (3.5.91)
Thus, $`\varphi ^{\left[d\right]}\left(z_s\right)`$ exists for $`2dn`$ since it follows the same arguments as above. Hence, $`\varphi \left(z_s\right)`$ that is $`𝒞^{r1}\left(^n\right)`$ smooth has been constructed.

#### 3.5.3 Lagrangian Submanifold

In consequence, if $`v`$ is in the stable subspace of $`_{\delta ,k}^F`$, then so does $`_{\delta ,k}^Fv`$. Since $`\left(_{\delta ,k}^F\right)^kv0`$ as $`k\mathrm{}`$, in the limit $`\mathrm{\Omega }(v,w)=0`$ under the dynamics (3.4.57). Thus, $`\mathrm{\Omega }(v,w)=0`$ if the tangent vectors are restricted in the stable subspace of $`_{\delta ,k}^F`$.

Recall the two-form $`\mathrm{\Omega }:T_{(x,\lambda )}M\times T_{(x,\lambda )}`$ in $`\mathrm{\S }4.4`$. If $`v`$ is $`W^s_{\delta ,k}^F`$, then so does $`_{\delta ,k}^Fv`$. Since $`\left(_{\delta ,k}^F\right)^kv0`$ as $`k\mathrm{}`$, in the limit $`\mathrm{\Omega }(v,w)=0`$ under the map (3.4.57). Thus, the two-form restricted to $`W^s\left(0\right)`$ is zero; i.e.

$`\mathrm{\Omega }(v,w)=0\text{for every v and w}TW^s.`$ (3.5.92)
In addition, $`W^s`$ has the maximal dimension of $`n`$. Hence, $`W^s`$ is a Lagrangian submanifold.

Recall that we represent $`W^s`$ as the graph $`z_u=\varphi \left(z_s\right)`$. The basis of $`T_{z_s,\varphi \left(z_s\right)}W^s`$ are of the form

$`{\displaystyle \frac{}{z_{si}}}\left[\begin{array}{c}z_{s1}\\ \mathrm{}\\ z_{sn}\\ \varphi _1\\ \mathrm{}\\ \varphi _n\end{array}\right]\text{where}i=1,\mathrm{},n.`$
Illuminating (3.5.92), we get

$`\left[\begin{array}{c}0\\ \mathrm{}\\ 1\\ \mathrm{}\\ 0\\ \frac{\varphi _1}{z_{si}}\\ \mathrm{}\\ \frac{\varphi _n}{z_{si}}\end{array}\right]^{}J\left[\begin{array}{c}0\\ \mathrm{}\\ \mathrm{}\\ 1\\ 0\\ \frac{\varphi _1}{z_{sj}}\\ \mathrm{}\\ \frac{\varphi _n}{z_{sj}}\end{array}\right]`$ (3.5.110)
where $`1`$ on the first vector is placed on the *ith* row and $`1`$ on second vector is on the *jth* row. Multiplying (3.5.110) out, we have the condition

$`{\displaystyle \frac{\varphi _i}{z_{sj}}}{\displaystyle \frac{\varphi _j}{z_{si}}}=0\text{ for }i,j=1,\mathrm{},n.`$ (3.5.111)
The equation (3.5.111) implies that $`\varphi \left(z_s\right)`$ is closed. Then, by the *Stokes’ Theorem* there exists $`\psi 𝒞^r\left(^n\right)`$ such that

$`\varphi \left(z_s\right)={\displaystyle \frac{\psi }{z_s}}\left(z_s\right)\text{ where }\psi \left(0\right)=0`$ (3.5.112)
locally on $`𝒩_ϵ\left(0\right)`$. Thus, we have shown that the Lagrangian submanifold is the gradient of $`\psi \left(z_s\right)`$.

### 3.6 The Optimal Cost

From the previous section, we have found that the $`z_u=\varphi \left(z_s\right)`$ is the gradient of the $`\psi \left(z_s\right)`$. In this section we find that the local stable manifold is also described by $`\lambda =\overline{\varphi }\left(x\right)`$, its original coordinates. We start by finding the linear term.

Consider the optimal control problem of minimizing

$$\underset{u}{\mathrm{min}}\frac{1}{2}\underset{j=0}{\overset{M}{}}\left(x_j^{}Q_kx_j+2x_j^{}S_ku_k+u_k^{}Ru_k\right)+x_T^{}Px_T$$

subject to the dynamics

$`x_{k+1}`$ $`=`$ $`A_kx_k+B_ku_k`$ (3.6.113)
$`x\left(0\right)`$ $`=`$ $`x_0`$
where the state vector $`x^n`$, the control $`u^m`$, and $`x_T^{}Px_T`$ is the terminal cost. We denote the fundamental solution matrix,

$`\left[\begin{array}{c}X_k\\ \mathrm{\Lambda }_k\end{array}\right]=\left[\begin{array}{cccc}x^1& x^2& \mathrm{}& x^n\\ \lambda ^1& \lambda ^2& \mathrm{}& \lambda ^n\end{array}\right]`$ (3.6.118)
of the dynamics

$`\left[\begin{array}{c}X_{k+1}\\ \mathrm{\Lambda }_{k+1}\end{array}\right]`$ $`=`$ $`_F\left[\begin{array}{c}X_k\\ \mathrm{\Lambda }_k\end{array}\right]`$ (3.6.123)
$`X\left(0\right)`$ $`=`$ $`Ix_0.`$
The columns

$`\left[\begin{array}{c}x^j\\ \lambda ^j\end{array}\right],j=1,\mathrm{},n`$
are $`n`$ linearly independent solutions and form a basis for the space of solutions.

Through the transformation

$$\mathrm{\Lambda }_k^{}=P_kX_k,$$

the time-varying discrete Riccati equation,

$$P_{k+1}=A_k^{}P_{k+1}A_k+\left(A_k^{}P_{k+1}B_k+S_k^{}\right)K_k+Q_k,$$

and the time-varying discrete control feedback,

$$K_k=\left(B_k^{}P_{k+1}B_k+R_k\right)^1\left(A_k^{}P_{k+1}B_k+S_k\right)^{},$$

are derived as shown in the proof below.

###### Theorem 3.6.1

Suppose ($`X_k,\mathrm{\Lambda }_k`$) is solution to (3.6.123) and $`X_k`$ is invertible, then $`P_k=\mathrm{\Lambda }_k^{}X_k^1`$ satisfies the time-varying Riccati equation.

Proof:  From the Pontryagin Maximum Principle, we have

$`U_k=R_k^1\left(B_k^{}\mathrm{\Lambda }_{k+1}^{}+S_k^{}X_k\right).`$ (3.6.125)
By letting $`\mathrm{\Lambda }_k^{}=P_kX_k`$ and substituting the dynamics (3.6.113),

$`U_k=R_k^1B_k^{}P_{k+1}\left(A_kX_k+B_kU_k\right)R_k^1S_k^{}X_k`$ (3.6.126)
which reduces to

$`\left(I+R_k^1B_k^{}P_{k+1}B\right)U_k=R_k^1B_k^{}P_{k+1}A_kX_kR_k^1S_k^{}X_k.`$ (3.6.127)
Multiplying (3.6.127) with $`R_k`$ becomes

$`U_k=K_kX_k.`$ (3.6.128)
where

$$K_k=\left(B_k^{}P_{k+1}B_k+R_k\right)^1\left(B_k^{}P_{k+1}A_k+S_k^{}\right),$$

the time-varying control feedback. We recall the costate equation for $`\mathrm{\Lambda }`$ in the Pontryagin Maximum Principle,

$`\mathrm{\Lambda }_k^{}=A_k^{}\mathrm{\Lambda }_{k+1}^{}+Q_kX_k+S_k^{}U_k`$ (3.6.129)
Substituting the dynamics (3.6.113) and $`U_k=K_kX_k`$, the costate equation becomes

$`\mathrm{\Lambda }_k^{}=A_k^{}P_{k+1}A_kX_k+\left(A_k^{}P_{k+1}B_k+S_k^{}\right)K_kX_k+Q_kX_k`$ (3.6.130)
Since the LHS of (3.6.130) is $`\mathrm{\Lambda }_k=P_kX_k`$, we obtain the time-varying Riccati equation as $`X_k`$ is cancelled out on both sides,

$`P_k=A_k^{}P_{k+1}A_k+\left(A_k^{}P_{k+1}B_k+S_k^{}\right)K_k+Q_k.`$ (3.6.131)
Thus, $`\mathrm{\Lambda }_k^{}=P_kX_k`$ solves the time-varying Riccati equation once $`(X_k,\mathrm{\Lambda }_k)`$ is found.

It follows that the DTARE (2.4.23) can be solved. By letting $`k\mathrm{}`$ in the dynamics (3.6.123), each linearly independent column vectors of $`(X_k,\mathrm{\Lambda }_k)^{}`$ converges to the stable direction and thus, forming the basis for the stable subspace. If $`X^1`$ exists, then fundamental matrix solution has the same span as

$`\left[\begin{array}{c}I\\ P_k\end{array}\right].`$
Hence, $`P_kP`$ as $`k\mathrm{}`$. Therefore, we solve the DTARE,

$$P=A^{}PA+\left(A^{}PB+S^{}\right)K+Q$$

as $`k\mathrm{}`$. Moreover, the stable linear subspace is

$`E_s=\left[\begin{array}{c}I\\ P\end{array}\right].`$ (3.6.135)
Thus, $`\lambda =Px`$ is the solution of the linear Hamiltonian equations. The linear term solution of the nonlinear Hamiltonian system is also $`\lambda =Px`$. We now show that there exists $`\lambda =\overline{\varphi }\left(x\right)`$ that describes the local stable manifold of the dynamics which has $`\lambda =Px`$ is its linear part.

Define the mapping $`F:𝒩_ϵ(0,0)^{2n}^n`$ by

$$F(x,\lambda )=\lambda Px+g(x,\lambda )$$

where $`g(x,\lambda )`$ contains only the higher order terms. Then, for every $`(x,\lambda )𝒩_ϵ(0,0)`$ $`F(x,\lambda )=0`$. Also the differential operator $`\frac{F}{\lambda }\left(0\right)=I`$ is invertible since $`\frac{g}{\lambda }\left(0\right)=0`$ because $`g`$ only contains the nonlinearities. We then invoke the IFT. It follows that there exists $`\lambda =\overline{\varphi }\left(x\right)`$ where $`\overline{\varphi }\left(x\right)U`$ where $`U`$ is some subset of $`^n`$ such that $`F(x,\overline{\varphi }\left(x\right))=0`$ for $`(x,\overline{\varphi }\left(x\right))𝒩_ϵ(0,0)`$. Thus, $`\lambda =\overline{\varphi }\left(x\right)`$ is the graph of the local stable manifold in the $`x,\lambda `$ coordinates.

### 3.7 Proof of the Main Result (Theorem 3.1.1)

Proof:  Recall the nonlinear Hamiltonian system (3.4.1). Assuming invertibility of $`H_{\lambda ,x}`$ at around 0 gives the nonlinear Hamiltonian system (3.4.33) that propagates in the direction of increasing time. When the nonlinear Hamiltonian sytem is linearized around the critical point zero, the system resembles that of the linear Hamiltonian system (3.4.41). We know that the linear Hamiltonian system is hyperbolic from (3.4.65) and (3.4.2). We then invoke the Local Stable Manifold Theorem as we have the dynamics (3.4.33) with a hyperbolic point $`0`$. Thus, there exists a local stable manifold $`W_s`$ such that the linear invariant subspace (3.6.135) is tangent at $`0`$. Since the linear and nonlinear Hamiltonian structure preserve symplectic form (3.4.65) and (3.4.64), it follows that if $`\mathrm{\Omega }`$ is restricted on $`W^s`$ then $`\mathrm{\Omega }`$ is zero (3.5.92); i.e., for $`u,v`$ are the basis in $`T_pW_{\left(0\right)}^s`$ then $`\mathrm{\Omega }(u,v)=0`$ as $`k\mathrm{}`$. Then, the manifold is a Lagrangian submanifold. The Lagrangian submanifold described by

$$z_u=\varphi \left(z_s\right)$$

satisfies (3.5.111). By Stokes’ Theorem there exists $`\psi 𝒞^r\left(^n\right)`$ such that

$`\varphi \left(z_s\right)={\displaystyle \frac{\psi }{z_s}}\left(z_s\right)\text{ where }\psi \left(0\right)=0.`$
We have also shown that $`\lambda =\overline{\varphi }\left(x\right)`$ describes the local stable manifold by invoking IFT. Moreover, the linear term of $`\lambda =\overline{\varphi }\left(x\right)`$ is $`\lambda =Px`$ which solves DTARE (2.4.23). It follows that

$`\lambda =\overline{\varphi }\left(x\right)={\displaystyle \frac{\pi }{x}}\left(x\right).`$ (3.7.136)
Thus, the graph of the gradient of the optimal cost is the local stable manifold of the associated Hamiltonian dynamics (3.4.33).

Lastly, we show that $`\pi \left(x\right)`$ and $`\kappa \left(x\right)`$ solves the DPE (2.3.13), (2.3.14). From the Pontryagin Maximum Principle and (3.7.136), we have

$$u^{}=\kappa \left(x\right)=\text{arg}\underset{v}{\mathrm{min}}H(x,\frac{\pi }{x}\left(x\right),v).$$

Equivalently,

$$\frac{H}{u}(x,\frac{\pi }{x}\left(x\right),\kappa \left(x\right))=0.$$

Then, $`\frac{\pi }{x}\left(x\right)`$ and $`\kappa \left(x\right)`$ solve DPE2. Since

$`\lambda ={\displaystyle \frac{H}{x}},`$ (3.7.137)
and (3.7.136), we have

$`{\displaystyle \frac{\pi }{x}}={\displaystyle \frac{\pi }{x}}\left(x\right){\displaystyle \frac{f}{x}}(x,u^{})+{\displaystyle \frac{l}{x}}(x,u^{})`$ (3.7.138)
Integrating (3.7.138) w.r.t. $`x`$, we get

$$\pi \left(x\right)\pi \left(f(x,\kappa \left(x\right))\right)l(x,\kappa \left(x\right))=0$$

which DPE1 (2.3.13). Thus, $`\pi `$ and $`\kappa `$ solve the DPE (2.3.13), (2.3.14). Furthermore $`\pi 𝒞^r`$ and $`\kappa 𝒞^{r1}`$ since $`l𝒞^r`$ and $`f𝒞^{r1}`$ in the Hamiltonian.

## Chapter 4 Local Stable Manifold Theorem for the Bidirectional Discrete-Time Dynamics

In this chapter, we look at the case where the linear map (3.3.9) is not a diffeomorphism. This is the case where zero is a closed loop eigenvalue and therefore the Hamiltonian matrix is not invertible. As in the previous case, we study the eigenstructure and symplectic properties of the mixed direction nonlinear Hamiltonian dynamics. Ultimately, we generalize the Local Stable Manifold Theorem for a bidirectional discrete map with a hyperbolic fixed point.

### 4.1 Discrete-Time Version of Gronwall’s Inequalities

We begin with a discrete-time version of Gronwall’s inequality. These lemmas will be useful in the proof of the local existence of a stable manifold.

###### Lemma 4.1.1

(Finite Difference Form) Suppose the sequence of scalars $`\{u_j\}_{j=0}^{\mathrm{}}`$ satisfies the difference inequality

$`u_{k+1}\delta u_k+L,`$ (4.1.1)
then

$$u_k\delta ^ku_0+L\underset{j=0}{\overset{k1}{}}\delta ^{k1j}.$$

The proof of the lemma above clearly follows from summing the equation (4.1.1) $`k`$-times.

###### Lemma 4.1.2

(Summation Form) Suppose $`\{\xi \}_{j=0}^{\mathrm{}}`$ is a sequence that satisfies

$$\left|\xi _k\right|C_1\underset{j=0}{\overset{k1}{}}\left|\xi _j\right|+C_2$$

and constants $`C_1,C_20`$, then

$$\left|\xi _k\right|C_2\underset{j=1}{\overset{k}{}}\left(1+C_1\right)^j$$

Proof:  Let $`s_k=_{j=0}^{k1}\left|\xi _j\right|`$. Then, the sequence $`\left\{s\right\}_{j=0}^{\mathrm{}}`$ satisfies

$$s_{k+1}\left(1+C_1\right)s_k+C_2$$

where $`C_1,C_20`$. By the discrete time form of Gronwall’s inequality,

$$\left|s_k\right|\left(1+C_1\right)^k\left|s_0\right|+C_2\underset{j=0}{\overset{k1}{}}\left(1+C_1\right)^{k1j}.$$

It follows that

$`\left|\xi _k\right|`$ $``$ $`\left(1+C_1\right)\left|s_k\right|+C_2`$
$``$ $`\left(1+C_1\right)\left[\left(1+C_1\right)^k\left|s_0\right|+C_2{\displaystyle \underset{j=0}{\overset{k1}{}}}\left(1+C_1\right)^j\right]`$
$``$ $`C_2{\displaystyle \underset{j=0}{\overset{k1}{}}}\left(1+C_1\right)^{kj}`$
$``$ $`C_2{\displaystyle \underset{j=1}{\overset{k}{}}}\left(1+C_1\right)^j`$
since $`\left|s_0\right|=0`$.

### 4.2 Local Stable Manifold for the Bidirectional Discrete-Time Dynamics

In this section we prove the existence of a local stable manifold $`\lambda =\varphi \left(x\right)`$ for the Hamiltonian dynamics,

$`\left[\begin{array}{c}x^+\\ \lambda \end{array}\right]`$ $`=`$ $`\left[\begin{array}{cc}A& BR^1B^{}\\ Q& A^{}\end{array}\right]\left[\begin{array}{c}x\\ \lambda ^+\end{array}\right]+\left[\begin{array}{c}F(x,\lambda ^+)\\ G(x,\lambda ^+)\end{array}\right]`$ (4.2.10)
where $`x,\lambda ^n`$ and zero is an eigenvalue of A. The nonlinear terms, $`F`$ and $`G`$, are $`C^k`$ functions for $`k1`$ such that

$`F(0,0)=0,`$ $`G(0,0)=0`$ (4.2.11)
$`{\displaystyle \frac{F}{(x,\lambda )}}(0,0)=0,`$ $`{\displaystyle \frac{G}{(x,\lambda )}}(0,0)=0.`$
The proof of the existence of a local stable manifold requires the discussion on the stability of the nonlinear state dynamics. The next subsection deals with the local asymptotic stability of the state dynamics. Consequently, we describe the diagonalization of the bidirectional Hamilton system. Finally, we show the existence of $`\lambda =\varphi \left(x\right)`$.

#### 4.2.1 Preliminaries

First, we introduce a $`C^{\mathrm{}}`$ cut-off function $`\rho \left(y\right):^n[0,1]`$ such that

$`\rho \left(y\right)=\{\begin{array}{cc}1,\hfill & \text{if }0\left|y\right|1\hfill \\ 0,\hfill & \text{if }\left|y\right|>2\hfill \end{array}`$
and $`0\rho \left(y\right)1`$ otherwise. Then we define the functions

$`F(x,\lambda ^+;ϵ)`$ $`:=`$ $`F(x\rho \left({\displaystyle \frac{x}{ϵ}}\right),\lambda ^+\rho \left({\displaystyle \frac{\lambda ^+}{ϵ}}\right))`$ (4.2.12)
$`G(x,\lambda ^+;ϵ)`$ $`:=`$ $`G(x\rho \left({\displaystyle \frac{x}{ϵ}}\right),\lambda ^+\rho \left({\displaystyle \frac{\lambda ^+}{ϵ}}\right))`$
for $`x,\lambda ^+^n`$. Since the $`F(x,\lambda ^+)`$ and $`G(x,\lambda ^+)`$ agree with $`F(x,\lambda ^+;ϵ)`$ and $`G(x,\lambda ^+;ϵ)`$, respectively, for $`\left|x\right|,\left|\lambda ^+\right|ϵ`$, it suffices to prove the existence of a stable manifold for some $`ϵ>0`$.

Now, we show that there exists $`N_1>0`$ and $`N_2>0`$ such that

$`\left|F(x,\lambda ;ϵ)F(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.13)
$`\left|G(x,\lambda ;ϵ)G(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.14)
and

$`\left|{\displaystyle \frac{F}{(x,\lambda )}}(x,\lambda ;ϵ){\displaystyle \frac{F}{(x,\lambda )}}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_2\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.15)
$`\left|{\displaystyle \frac{G}{(x,\lambda )}}(x,\lambda ;ϵ){\displaystyle \frac{G}{(x,\lambda )}}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_2\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right].`$ (4.2.16)
Since $`\rho \left(y\right)`$ and its partial derivatives are continuous functions with compact support there exists $`M>0`$ such that

$`\left|{\displaystyle \frac{\rho }{y}}\left(y\right)\right|`$ $``$ $`M`$
$`\left|{\displaystyle \frac{^2\rho }{y^2}}\left(y\right)\right|`$ $``$ $`M`$
for all $`\lambda ^n`$. We also choose $`M>0`$ large enough that

$`\left|{\displaystyle \frac{F}{x}}(x,\lambda )\right|`$ $``$ $`M\left|x\right|`$ (4.2.17)
$`\left|{\displaystyle \frac{F}{\lambda }}(x,\lambda )\right|`$ $``$ $`M\left|\lambda \right|`$ (4.2.18)
$`\left|{\displaystyle \frac{^2F}{x^i\lambda ^j}}(x,\lambda )\right|`$ $``$ $`M,i,j=1,2`$ (4.2.19)
because of the condition (4.2.11) for $`\left|x\right|,\left|\lambda \right|<1`$. By the Mean Value Theorem,

$`\left|F(x,\lambda ;ϵ)F(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`\left|F(x,\lambda ;ϵ)F(\stackrel{~}{x},\lambda ;ϵ)+F(\stackrel{~}{x},\lambda ;ϵ)F(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$
$``$ $`\left|F(x,\lambda ;ϵ)F(\stackrel{~}{x},\lambda ;ϵ)\right|+\left|F(\stackrel{~}{x},\lambda ;ϵ)F(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$
$``$ $`\left|{\displaystyle \frac{F}{x}}(\xi _1,\lambda ;ϵ)\right|\left|x\stackrel{~}{x}\right|+\left|{\displaystyle \frac{F}{\lambda }}(x,\xi _2;ϵ)\right|\left|\lambda \stackrel{~}{\lambda }\right|`$
where $`\xi _1`$ is between $`x`$ and $`\stackrel{~}{x}`$ and $`\xi _2`$ is between $`\lambda `$ and $`\stackrel{~}{\lambda }`$. Similarly,

$`\left|{\displaystyle \frac{F}{x}}(x,\lambda ;ϵ){\displaystyle \frac{F}{x}}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`\left|{\displaystyle \frac{^2F}{x^2}}(\xi _1,\lambda )\right|\left|x\stackrel{~}{x}\right|+\left|{\displaystyle \frac{^2F}{x\lambda }}(x,\xi _2)\right|\left|\lambda \stackrel{~}{\lambda }\right|`$
$`\left|{\displaystyle \frac{F}{\lambda }}(x,\lambda ;ϵ){\displaystyle \frac{F}{\lambda }}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`\left|{\displaystyle \frac{^2F}{\lambda x}}(\xi _1,\lambda )\right|\left|x\stackrel{~}{x}\right|+\left|{\displaystyle \frac{^2F}{\lambda ^2}}(x,\xi _2)\right|\left|\lambda \stackrel{~}{\lambda }\right|.`$
Next we estimate for $`0ϵ<\frac{1}{2}`$

$`\left|{\displaystyle \frac{F}{x}}(\xi _1,\xi _2;ϵ)\right|`$ $`=`$ $`\left|{\displaystyle \frac{F}{x}}(\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)\xi _1,{\displaystyle \frac{\xi _2}{ϵ}})\xi _2\right)\left|\right|{\displaystyle \frac{\rho }{y}}\left({\displaystyle \frac{\xi _1}{ϵ}}\right){\displaystyle \frac{\xi _1}{ϵ}}+\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)|`$
$``$ $`M\left|\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)\right|\left|\xi _1\right|\left(\left|{\displaystyle \frac{\rho }{y}}\left({\displaystyle \frac{\xi _1}{ϵ}}\right){\displaystyle \frac{\xi _1}{ϵ}}\right|+\left|\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)\right|\right)`$
$``$ $`M\left[M+1\right]ϵ`$
and

$`\left|{\displaystyle \frac{^2F}{x^2}}(\xi _1,\xi _2;ϵ)\right|`$ $`=`$ $`\left|{\displaystyle \frac{^2F}{x^2}}(\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)\xi _1,{\displaystyle \frac{\xi _2}{ϵ}})\xi _2\right)\left|\right|{\displaystyle \frac{\rho }{y}}\left({\displaystyle \frac{\xi _1}{ϵ}}\right){\displaystyle \frac{\xi _1}{ϵ}}+\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)|`$
$`+\left|{\displaystyle \frac{F}{x}}(\rho \left({\displaystyle \frac{\xi _1}{ϵ}}\right)\xi _1,{\displaystyle \frac{\xi _2}{ϵ}})\xi _2\right)\left|\right|{\displaystyle \frac{\rho }{\xi _1}}\left({\displaystyle \frac{\xi _1}{ϵ}}\right){\displaystyle \frac{2}{ϵ}}+{\displaystyle \frac{^2\rho }{y^2}}\left({\displaystyle \frac{\xi _1}{ϵ}}\right){\displaystyle \frac{\xi _1}{ϵ^2}}|`$
$``$ $`M\left(1+2M\right)+8M^2`$
for $`\left|\xi _1\right|,\left|\xi _2\right|ϵ`$.

Let

$`N_1`$ $`=`$ $`M^2\left[M+1\right]`$
$`N_2`$ $`=`$ $`M\left(1+2M\right)+8M^2.`$
It follows that

$`\left|{\displaystyle \frac{F}{(x,\lambda )}}(\xi _1,\xi _2;ϵ)\right|`$ $``$ $`N_1ϵ`$ (4.2.20)
$`\left|{\displaystyle \frac{^2F}{x^i\lambda ^j}}(\xi _1,\xi _2;ϵ)\right|`$ $``$ $`N_2,i,j=1,2`$ (4.2.21)
for $`\left|\xi _1\right|,\left|\xi _2\right|<ϵ`$. The inequalities above also hold $`G`$ as well.

Thus,

$`\left|F(x,\lambda ;ϵ)F(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.22)
$`\left|G(x,\lambda ;ϵ)G(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.23)
and

$`\left|{\displaystyle \frac{F}{(x,\lambda )}}(x,\lambda ;ϵ){\displaystyle \frac{F}{(x,\lambda )}}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_2\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right]`$ (4.2.24)
$`\left|{\displaystyle \frac{G}{(x,\lambda )}}(x,\lambda ;ϵ){\displaystyle \frac{G}{(x,\lambda )}}(\stackrel{~}{x},\stackrel{~}{\lambda };ϵ)\right|`$ $``$ $`N_2\left[\left|x\stackrel{~}{x}\right|+\left|\lambda \stackrel{~}{\lambda }\right|\right].`$ (4.2.25)
Henceforth we suppress $`ϵ`$ and write $`F(x,\lambda ),G(x,\lambda )`$ for $`F(x,\lambda ;ϵ),G(x,\lambda ;ϵ)`$.

#### 4.2.2 Stability of the Nonlinear Dynamics

The stable invariant manifold is described by $`\lambda =Px`$ for the linear bidirectional Hamiltonian dynamics where $`P`$ is the solution to the discrete algebraic Riccati equation. Then, we know from Chapter 3 that the linear term of the stable manifold for the nonlinear bidirectional Hamiltonian dynamics, $`\lambda =\varphi \left(x\right)`$, is $`Px`$. Thus, the local stable manifold is of the form

$`\lambda =\varphi \left(x\right)=Px+\psi \left(x\right)`$ (4.2.26)
where $`\psi \left(x\right)`$ contains all the nonlinear terms.

Suppose we substitute (4.2.26) into the state dynamics in (4.2.10), then the nonlinear state dynamics becomes

$`\left(I+BR^1B^{}P\right)x^+=AxBR^1B^{}\psi \left(x^+\right)+F(x,Px^++\psi \left(x^+\right)).`$
By the Matrix Inversion Lemma (), we have that

$$\left(I+BR^1B^{}P\right)^1=\left(IB\left(B^{}PB+R\right)^1B^{}P\right).$$

Then, it follows that

$`x^+`$ $`=`$ $`\left(A+BK\right)x+f_\psi (x,x^+)`$ (4.2.27)
$`x\left(0\right)`$ $`=`$ $`x_0`$
where $`K=\left(B^{}PB+R\right)^1B^{}P`$ and $`f_\psi (x,x^+)=\left(I+BR^1B^{}P\right)^1\left(F(x,Px^++\psi \left(x^+\right))BR^1B^{}\psi \left(x^+\right)\right)`$.

The implicit equation above can be solved. Let $`:𝒩_ϵ\left(0\right)^{2n}^n`$ such that

$$(x,x^+)=x^+\left(A+BK\right)xf_\psi (x,x^+)=0$$

for $`x,x^+𝒩_ϵ\left(0\right)`$. Then, for $`0𝒩_ϵ\left(0\right)`$ the Jacobian

$`{\displaystyle \frac{}{x^+}}\left(0\right)`$ $`=`$ $`I{\displaystyle \frac{f_\psi }{x^+}}\left(0\right)`$
$`=`$ $`I\left(I+BR^1B^{}P\right)^1\left({\displaystyle \frac{F}{\lambda ^+}}(0,\varphi \left(0\right)){\displaystyle \frac{\psi }{x^+}}\left(0\right)BR^1B^{}{\displaystyle \frac{\psi }{x^+}}\left(0\right)\right)`$
$`=`$ $`I,`$
because of the condition (4.2.11) and $`\psi \left(x^+\right)`$ only contains nonlinear terms. Then, by the Implicit Function Theorem there exists $`𝔽\left(x\right)`$ such that

$`x^+`$ $`=`$ $`𝔽\left(x\right)`$ (4.2.28)
is equivalent to the earlier state dynamics (4.2.27). Moreover, the linear term of $`𝔽\left(x\right)`$ is $`\left(A+BK\right)x`$, i.e.;

$`𝔽\left(x\right)=\left(A+BK\right)x+F_\psi \left(x\right)`$ (4.2.29)
because

$`{\displaystyle \frac{𝔽}{x}}\left(0\right)`$ $`=`$ $`\left({\displaystyle \frac{}{x^+}}\left(0\right)\right)^1{\displaystyle \frac{}{x}}\left(0\right)`$
$`=`$ $`I\left[\left(A+BK\right)\right]`$
$`=`$ $`A+BK.`$
It follows that $`F_\psi \left(x\right)`$ contains only the nonlinear terms and thus,

$`F_\psi \left(0\right)=0`$ (4.2.30)
and

$$\frac{F_\psi }{x_i}\left(0\right)=0i=1,\mathrm{},n.$$

The linear part of (4.2.27) is

$`x^+=\left(A+BK\right)x.`$ (4.2.31)
Since the eigenvalues of $`\left(A+BK\right)`$ lie strictly inside the unit circle, the term $`\left(A+BK\right)^kx_00`$ as $`k\mathrm{}`$. Thus, the system (4.2.31) is asymptotically stable. Also, it implies that there exists a unique positive definite $`P`$ that satisfies the Lyapunov equation

$$\left(A+BK\right)^{}P\left(A+BK\right)P=I.$$

Now we show the stability of the nonlinear dynamics

$`x_k`$ $`=`$ $`\left(A+BK\right)x+F_\psi \left(x\right)`$
$`x\left(0\right)`$ $`=`$ $`0.`$
We must prove that

$`\underset{x0}{lim}{\displaystyle \frac{\left|F_\psi \left(x\right)\right|}{\left|x\right|}}=0;`$
i.e., given any $`\epsilon >0`$ and any $`\psi \left(x\right)`$ satisfying the conditions

$`\psi \left(0\right)`$ $`=`$ $`0`$ (4.2.32)
$`\left|\psi \left(x\right)\psi \left(\overline{x}\right)\right|`$ $``$ $`l\left(ϵ\right)\left|x\overline{x}\right|`$ (4.2.33)
where $`l\left(ϵ\right)0`$ as $`ϵ0`$, there exists $`\delta >0`$ such that

$$\frac{\left|F_\psi \left(x\right)\right|}{\left|x\right|}<\epsilon \text{whenever}\left|x\right|<\delta .$$

We define $`\psi \left(x\right)`$ to be the nonlinear term of the stable manifold in (3.5.90). The conditions (4.2.32) and (4.2.33) will be necessary for the proof of the local stable manifold theorem.

Recall that

$`x^+=𝔽\left(x\right)=\left(A+BK\right)x+F_\psi \left(x\right).`$
Then,

$`0=(x,x^+)`$ $`=`$ $`x^+\left(A+BK\right)xf_\psi (x,x^+)`$
$`=`$ $`\left(A+BK\right)x+F_\psi \left(x\right)\left(A+BK\right)xf_\psi (x,\left(A+BK\right)x+F_\psi \left(x\right))`$
$`=`$ $`F_\psi \left(x\right)f_\psi (x,\left(A+BK\right)x+F_\psi \left(x\right)).`$
It follows that

$`F_\psi \left(x\right)`$ $`=`$ $`(I+BR^1B^{}P)^1[F(x,P((A+BK)x+F_\psi \left(x\right))+\psi ((A+BK)x+F_\psi \left(x\right)))`$ (4.2.34)
$`BR^1B^{}\psi ((A+BK)x+F_\psi \left(x\right))].`$
Let $`𝔹_1=\left(I+BR^1B^{}P\right)^1`$, $`𝔹_2=BR^1B^{}`$, $`=P`$ and $`\alpha =\mathrm{max}_i\left|\lambda _i\right|`$ where $`\lambda _i\sigma \left(A+BK\right)`$ and $`\left|\lambda _i\right|<1`$. We have the following from (4.2.34),

$`\left|F_\psi \left(x\right)\right|`$ $``$ $`𝔹_1N_1ϵ\left[\right|x|+|P(A+BK)x+PF_\psi \left(x\right)+\psi ((A+BK)x+F_\psi \left(x\right)))\left|\right]`$
$`+𝔹_1𝔹_2\left|\psi ((A+BK)x+F_\psi \left(x\right))\right)|.`$
because of (4.2.22). Using the Lipschitz condition (4.2.33) for $`\psi \left(x\right)`$,

$`\left|F_\psi \left(x\right)\right|`$ $``$ $`\left[𝔹_1N_1ϵ+\alpha \left(𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_1𝔹_2l\left(ϵ\right)\right)\right]\left|x\right|`$
$`+\left[𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_1𝔹_2l\left(ϵ\right)\right]\left|F_\psi \left(x\right)\right|.`$
Solving for $`\left|F_\psi \left(x\right)\right|`$,

$`\left|F_\psi \left(x\right)\right|`$ $``$ $`{\displaystyle \frac{𝔹_1N_1ϵ+\alpha \left(𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}{1\left(𝔹_1N_1ϵ+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}}.`$
Let $`\delta =\frac{1\left(𝔹_1N_1ϵ+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}{𝔹_1N_1ϵ+\alpha \left(𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}\epsilon `$. For some $`ϵ>0`$ and $`\epsilon >0`$, we have that $`\delta >0`$. Then,

$`\left|F_\varphi \left(x\right)\right|`$ $``$ $`{\displaystyle \frac{𝔹_1N_1ϵ+\alpha \left(𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}{1\left(𝔹_1N_1ϵ+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}}\left|x\right|`$
$``$ $`{\displaystyle \frac{𝔹_1N_1ϵ+\alpha \left(𝔹_1N_1ϵP+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}{1\left(𝔹_1N_1ϵ+𝔹_1N_1ϵl\left(ϵ\right)+𝔹_2l\left(ϵ\right)\right)}}\delta `$
$``$ $`\epsilon .`$
Thus,

$$F_\varphi \left(x\right)=o\left(\left|x\right|\right).$$

Now, we use the Lyapunov argument. Let $`v\left(x\right)=x^{}Px`$. Then,

$`\mathrm{\Delta }v\left(x\right)`$ $`=`$ $`v\left(x^+\right)v\left(x\right)`$
$`=`$ $`x_{}^{+}{}_{}{}^{}Px^+x^{}Px`$
$`=`$ $`\left[\left(A+BK\right)xF_\psi \left(x\right)\right]^{}P\left[\left(A+BK\right)xF_\psi \left(x\right)\right]x^{}Px`$
$`=`$ $`x^{}\left(\left(A+BK\right)^{}P\left(A+BK\right)P\right)x+2x^{}\left(A+BK\right)^{}PF_\psi \left(x\right)`$
$`=`$ $`\left|x\right|^2+2x^{}\left(A+BK\right)^{}PF_\psi \left(x\right).`$
since

$`\left|F_\psi \left(x\right)\right|{\displaystyle \frac{1}{3p}}\left|x\right|`$
and

$$\left|2x^{}\left(A+BK\right)^{}PF_\psi \left(x\right)\right|\frac{2}{3}\left|x\right|^2$$

for some $`p>0`$. Thus,

$`\mathrm{\Delta }v\left(x\right)={\displaystyle \frac{\left|x\right|^2}{3}}<0.`$
Therefore, the nonlinear dynamics is locally asymptotically stable uniform for all $`\psi 𝕏`$.

#### 4.2.3 Diagonalization of the Hamiltonian Matrix

Recall the bidirectional nonlinear dynamics in (4.2.10)

$`\left[\begin{array}{c}x^+\\ \lambda \end{array}\right]`$ $`=`$ $`\left[\begin{array}{cc}A& BR^1B^{}\\ Q& A^{}\end{array}\right]\left[\begin{array}{c}x\\ \lambda ^+\end{array}\right]+\left[\begin{array}{c}F(x,\lambda ^+)\\ G(x,\lambda ^+)\end{array}\right]`$
where $`x,\lambda ^n`$ and zero is an eigenvalue of A. The nonlinear terms, $`F`$ and $`G`$, are $`C^k`$ functions for $`k1`$ such that

$`F(0,0)=0,`$ $`G(0,0)=0`$
$`{\displaystyle \frac{F}{(x,\lambda )}}(0,0)=0,`$ $`{\displaystyle \frac{G}{(x,\lambda )}}(0,0)=0.`$
By substituting

$`\lambda =Px+\psi \left(x\right)`$ (4.2.36)
into the state dynamics above (4.2.10), we get

$`x^+`$ $`=`$ $`\left(A+BK\right)x+f_\psi (x,x^+)`$ (4.2.37)
where $`f_\psi (x,x^+)=\left(I+BR^1B^{}P\right)^1\left(F(x,Px^++\psi \left(x^+\right))BR^1B^{}\psi \left(x^+\right)\right)`$.

As we substitute (4.2.36) and (4.2.37) into the costate dynamics in (4.2.10), we also add $`0=\left(BK\right)^{}\lambda ^+\left(BK\right)^{}\lambda ^+`$. Then, the costate dynamics becomes

$`\lambda =\left(A+BK\right)^{}\lambda ^++\overline{Q}x+g_\psi (x,x^+).`$
where $`\overline{Q}=QK^{}B^{}P\left(A+BK\right)`$ and $`g_\psi (x,x^+)=G(x,Px^++\psi \left(x^+\right))+K^{}B^{}(\psi \left(x^+\right)Pf_\psi (x,x^+)`$.

Thus, the substitution of

$$\lambda =Px+\psi \left(x\right)$$

into the dynamics (4.2.10) results in a new nonlinear dynamics

$`\left[\begin{array}{c}x^+\\ \lambda \end{array}\right]`$ $`=`$ $`\left[\begin{array}{cc}A+BK& 0\\ \overline{Q}& \left(A+BK\right)^{}\end{array}\right]\left[\begin{array}{c}x\\ \lambda ^+\end{array}\right]+\left[\begin{array}{c}f_\psi (x,x^+)\\ g_\psi (x,x^+)\end{array}\right]`$ (4.2.46)
where

$`f_\psi (x,x^+)`$ $`=`$ $`\left(I+BR^1B^{}P\right)^1\left(F(x,\psi \left(x^+\right))BR^1B^{}\psi \left(x^+\right)\right)`$
and

$`g_\psi (x,x^+)`$ $`=`$ $`G(x,\psi \left(x^+\right))+K^{}B^{}(\psi \left(x^+\right)Pf_\psi (x,x^+).`$
The nonlinear terms $`f_\psi `$ and $`g_\psi `$ are $`C^k`$ functions for $`k1`$ such that

$`f_\psi (0,0)=0,`$ $`g_\psi (0,0)=0`$ (4.2.47)
$`{\displaystyle \frac{f_\psi }{(x,x^+)}}(0,0)=0,`$ $`{\displaystyle \frac{g_\psi }{(x,x^+)}}(0,0)=0.`$ (4.2.48)
because of (4.2.11), $`\psi \left(x\right)`$ only contains nonlinear terms and

$$\frac{\psi }{x}\left(0\right)=0.$$

Now we introduce the $`z`$ coordinate by the transformation

$`\lambda =z+Sx`$ (4.2.49)
for some matrix $`S`$ to block diagonalize the block lower triangular Hamiltonian matrix in (4.2.46). By substitution, the system (4.2.46) becomes

$`x^+`$ $`=`$ $`\left(A+BK\right)x+f_\psi (x,x^+)`$
$`z`$ $`=`$ $`\left(A+BK\right)^{}z^++\left(A+BK\right)^{}S\left(A+BK\right)xSx+\overline{Q}x+h_\psi (x,x^+)`$
where

$$h_\psi (x,x^+)=\left(A+BK\right)^{}Sf_\psi (x,x^+)+g_\psi (x,x^+).$$

Observe from the $`z`$ dynamics above that the terms

$$\left(A+BK\right)^{}S\left(A+BK\right)xSx+\overline{Q}x=0$$

and recall that $`\overline{Q}=QK^{}B^{}P\left(A+BK\right)`$. Indeed,

$`S+A^{}S\left(A+BK\right)+K^{}B^{}S\left(A+BK\right)=Q+K^{}B^{}P\left(A+BK\right).`$ (4.2.50)
We know that

$`S+A^{}S\left(A+BK\right)=Q,`$ (4.2.51)
is the discrete-time algebraic Riccati equation (DTARE). Subtracting (4.2.51) from (4.2.50), we have

$`K^{}B^{}S\left(A+BK\right)`$ $`=`$ $`K^{}B^{}P\left(A+BK\right).`$
Thus,

$`S=P`$ (4.2.52)
and $`S`$ satisfies the DTARE. Therefore, we have a diagonalized system

$`\left[\begin{array}{c}x^+\\ z\end{array}\right]`$ $`=`$ $`\left[\begin{array}{cc}A+BK& 0\\ 0& \left(A+BK\right)^{}\end{array}\right]\left[\begin{array}{c}x\\ z^+\end{array}\right]+\left[\begin{array}{c}f_\psi (x,x^+)\\ g_\psi (x,x^+)\end{array}\right]`$ (4.2.61)
where

$`f_\psi (x,x^+)`$ $`=`$ $`\left(I+BR^1B^{}P\right)^1\left(F(x,\psi \left(x^+\right))BR^1B^{}\psi \left(x^+\right)\right)`$
and

$`h_\psi (x,x^+)=\left(A+BK\right)^{}Sf_\psi (x,x^+)+g_\psi (x,x^+).`$
The nonlinear terms $`f_\psi `$ and $`h_\psi `$ are $`C^k`$ functions for $`k1`$ such that

$`f_\psi (0,0)=0,`$ $`h_\psi (0,0)=0`$ (4.2.62)
$`{\displaystyle \frac{f_\psi }{(x,x^+)}}(0,0)=0,`$ $`{\displaystyle \frac{h_\psi }{(x,x^+)}}(0,0)=0.`$ (4.2.63)
because of (4.2.47) and (4.2.48).

#### 4.2.4 The Local Stable Manifold Theorem

Given the original dynamics (4.2.10) we look for the local stable manifold described by $`\lambda =\varphi \left(x\right)`$. From the Theorem (3.5.90) we have already proven that the linear term of the stable manifold for the system (4.2.10) is $`Px.`$ where $`P`$ is the solution to DTARE. Therefore, the local stable manifold is

$`\lambda =Px+\psi \left(x\right)`$ (4.2.64)
where $`\psi \left(x\right)`$ only contains the nonlinear term of $`\varphi \left(x\right)`$. In the two-step process of diagonalization of the system (4.2.10), we introduce the $`z`$ coordinate through the transformation

$$\lambda =z+Sx.$$

Since $`S=P`$, it must be that

$`z=\lambda Px=Px+\psi \left(x\right)Px=\psi \left(x\right).`$
Then, it suffices to prove existence of the local stable manifold $`z=\psi \left(x\right)`$ for the diagonalized system (4.2.61). In order to show the existence of $`z=\psi \left(x\right)`$, we use the Contraction Mapping Principle (CMP). To invoke the CMP, we will need a map $`T:𝕏𝕏`$ that is a contraction on a complete metric space $`𝕏`$.

###### Theorem 4.2.1

Given the dynamics in (4.2.61) with the nonlinear terms $`f_\psi `$ and $`h_\psi `$ are $`𝒞^k`$ functions satisfying the conditions (4.2.62) and (4.2.63) and a hyperbolic fixed point $`0^{2n}`$, there exists a local stable manifold $`z=\psi (x)`$ around the fixed point $`0`$ where $`\psi `$ is a $`C^k`$ function.

Proof:

First notice that $`f_\psi `$ and $`g_\psi `$ are cut-off functions. It follows that $`h_\psi `$ is also a cut-off function. It suffices to prove the theorem for some $`ϵ>0`$ since the cut-off functions $`f_\psi (x,x^+;ϵ)`$ and $`h_\psi (x,x^+;ϵ)`$ agree with $`f_\psi (x,x^+)`$ and $`h_\psi (x,x^+)`$ for $`\left|x\right|,\left|x^+\right|ϵ`$.

By (4.2.22)-(4.2.25), (4.2.32)-(4.2.33), and (4.2.62)-(4.2.63), there exists $`\overline{N}_1,\overline{N}_2>0`$ such that

$`\left|f_\psi (x,y;ϵ)f_\psi (\stackrel{~}{x},\stackrel{~}{y};ϵ)\right|`$ $``$ $`\overline{N}_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|y\stackrel{~}{y}\right|\right]`$ (4.2.65)
$`\left|h_\psi (x,y;ϵ)h_\psi (\stackrel{~}{x},\stackrel{~}{y};ϵ)\right|`$ $``$ $`\overline{N}_1ϵ\left[\left|x\stackrel{~}{x}\right|+\left|y\stackrel{~}{y}\right|\right]`$ (4.2.66)
and

$`\left|{\displaystyle \frac{f_\psi }{(x,y)}}(x,y;ϵ){\displaystyle \frac{f_\psi }{(x,y)}}(\stackrel{~}{x},\stackrel{~}{y};ϵ)\right|`$ $``$ $`\overline{N}_2\left[\left|x\stackrel{~}{x}\right|+\left|y\stackrel{~}{y}\right|\right]`$ (4.2.67)
$`\left|{\displaystyle \frac{h_\psi }{(x,y)}}(x,y;ϵ){\displaystyle \frac{h_\psi }{(x,y)}}(\stackrel{~}{x},\stackrel{~}{y};ϵ)\right|`$ $``$ $`\overline{N}_2\left[\left|x\stackrel{~}{x}\right|+\left|y\stackrel{~}{y}\right|\right].`$ (4.2.68)
Moreover, from the bound (4.2.20) we know the following

$`\left|{\displaystyle \frac{f_\psi }{(x,y)}}(\xi _1,\xi _2;ϵ)\right|\overline{N}ϵ`$ (4.2.69)
and

$`\left|{\displaystyle \frac{h_\psi }{(x,y)}}(\xi _1,\xi _2;ϵ)\right|\overline{N}ϵ.`$ (4.2.70)
for $`\overline{N}>0`$ and $`\left|\xi _1\right|,\left|\xi _2\right|<ϵ`$.

Henceforth we suppress $`ϵ`$ and write $`f_\psi (x,y),h_\psi (x,y)`$ for $`f_\psi (x,y;ϵ),h_\psi (x,y;ϵ)`$.

Let $`l\left(ϵ\right)`$ with $`l\left(0\right)=0`$ and $`\psi 𝕏^0\left(\left|x\right|ϵ\right)`$ where $`𝕏`$ is space of $`\psi :^n^n`$ such that

$`\psi \left(0\right)`$ $`=`$ $`0`$ (4.2.71)
$`\left|\psi \left(x\right)\psi \left(\overline{x}\right)\right|`$ $``$ $`l\left(ϵ\right)\left|x\overline{x}\right|`$ (4.2.72)
for $`x,\overline{x}\overline{}_ϵ\left(0\right)^n`$. We define

$`\psi =\underset{\left|x\right|ϵ}{sup}\left|{\displaystyle \frac{\psi \left(x\right)}{x}}\right|.`$ (4.2.73)
Since $`𝕏^0\left(\left\{\left|x\right|ϵ\right\}\right)`$, to show $`𝕏`$ is a complete metric space it suffices to show that $`𝕏`$ is closed. We take a sequence $`\left\{\psi _n\right\}𝕏`$ such that $`\psi _n\psi `$ in $`^0`$ norm. For large $`N>n`$, $`\left|\psi _n\left(x\right)\psi \left(x\right)\right|\frac{ϵ}{2}`$ for all $`x\overline{}_ϵ\left(0\right)`$. Then,

$`\left|\psi \left(x\right)\psi \left(\overline{x}\right)\right|`$ $``$ $`\left|\psi \left(x\right)\psi _n\left(x\right)\right|+\left|\psi _n\left(x\right)\psi _n\left(\overline{x}\right)\right|+\left|\psi _n\left(\overline{x}\right)\psi \left(\overline{x}\right)\right|`$
$``$ $`{\displaystyle \frac{ϵ}{2}}+l\left(ϵ\right)\left|x\overline{x}\right|+{\displaystyle \frac{ϵ}{2}}.`$
By letting $`ϵ0`$, we have that $`\left|\psi \left(x\right)\psi \left(\overline{x}\right)\right|l\left(ϵ\right)\left|x\overline{x}\right|`$. Thus, $`\psi `$ is a Lipschitz function. Similarly, the condition (4.2.71) is easily satisfied. It follows that

$`\left|\psi \left(0\right)0\right|`$ $``$ $`\left|\psi \left(0\right)\psi _n\left(0\right)\right|+\left|\psi _n\left(0\right)0\right|ϵ.`$
Thus, $`\psi 𝕏`$. Hence $`𝕏`$ is closed. Moreover, $`𝕏`$ is a complete metric space with the norm defined on (4.2.73).

Solving the $`z`$ dynamics in (4.2.61) via the variation of constants formula, we have

$`z_j`$ $`=`$ $`\left(A^{}+K^{}B^{}\right)^{kj}z_k+{\displaystyle \underset{l=j}{\overset{k1}{}}}\left(A^{}+K^{}B^{}\right)^{lj}h_\psi (x_l,x_{l+1})`$
for $`j<k`$. Let $`j=0`$ and $`k=\mathrm{}`$, then (4.2.4) changes to

$`z_0`$ $`=`$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left(A^{}+K^{}B^{}\right)^lh_\psi (x_l,x_{l+1})`$
We define a mapping $`T:𝕏𝕏`$ by

$`\left(T\psi \right)\left(x_0\right)={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left(A^{}+K^{}B^{}\right)^lh_\psi (x_l,x_{l+1}).`$ (4.2.76)
From this fixed point equation, we look for the solution

$`\left(T\psi \right)\left(x_0\right)=\psi \left(x_0\right).`$
We must show $`T\psi 𝕏`$ and prove $`T`$ is a contraction on $`𝕏`$.

Suppose $`x_0=0`$ is the initial condition. Clearly from the equation (4.2.28)-(4.2.30), $`x_k=0`$ for all $`k`$. Together with $`x_k=0`$ for all $`k`$ and the condition (4.2.63) $`h_\psi (0,0)=0`$, we have

$`T\psi \left(0\right)=0.`$ (4.2.77)
Hence $`T\psi `$ satisfies the condition (4.2.71).

We now prove the Lipschitz condition (4.2.72) for $`T\psi `$.

For $`\psi 𝕏`$ and the initial conditions $`x_0,\overline{x}_0^n`$, we denote $`x_k=x(k,x_0,\psi )`$ to be the solution of the state dynamics.

$`x^+`$ $`=`$ $`\left(A+BK\right)x+f_\psi (x,x^+)`$
$`x\left(0\right)`$ $`=`$ $`x_0.`$
Similarly, for $`\psi 𝕏`$ and the initial conditions $`\overline{x}_0^n`$, let $`x_k=x(k,\overline{x}_0,\psi )`$ tbe the solution of

$`x^+`$ $`=`$ $`\left(A+BK\right)x+f_\psi (x,x^+)`$
$`x\left(0\right)`$ $`=`$ $`\overline{x}_0.`$
Recall $`\alpha =\mathrm{max}_j\left|\lambda _j\right|`$ where $`\lambda _j\sigma \left(A+BK\right)`$ and $`\left|\lambda _j\right|<1`$. Using the estimate (4.2.65), at one-time step

$`\left|x_{k+1}\overline{x}_{k+1}\right|\alpha \left|x_k\overline{x}_k\right|+\overline{N}_1\left[\left|x_k\overline{x}_k\right|+\left|x_{k+1}\overline{x}_{k+1}\right|\right].`$
For $`1\overline{N}_1ϵ>0`$,

$`\left|x_{k+1}\overline{x}_{k+1}\right|{\displaystyle \frac{\alpha +\overline{N}_1ϵ}{1\overline{N}_1ϵ}}\left|x_k\overline{x}_k\right|`$
and recursively,

$`\left|x_k\overline{x}_k\right|\left({\displaystyle \frac{\alpha +\overline{N}_1ϵ}{1\overline{N}_1ϵ}}\right)^k\left|x_0\overline{x}_0\right|`$
As long as

$$ϵ<\frac{1\alpha }{2\overline{N}_1},$$

then

$$\frac{\alpha +\overline{N}_1ϵ}{1\overline{N}_1ϵ}<1.$$

Thus,

$`\left|x_k\overline{x}_k\right|\left|x_0\overline{x}_0\right|.`$ (4.2.78)
Using the bounds (4.2.66) and (4.2.78),

$`\left|T\psi \left(x_0\right)T\psi \left(\overline{x}_0\right)\right|`$ $``$ $`\left|{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left(A^{}+K^{}B^{}\right)^l\left(h_\psi (x_l,x_{l+1})h_\psi (\overline{x}_l,\overline{x}_{l+1})\right)\right|`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l\left|h_\psi (x_l,x_{l+1})h_\psi (\overline{x}_l,\overline{x}_{l+1})\right|`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l\overline{N}_1ϵ\left[\left|x_l\overline{x}_l\right|+\left|x_{l+1}\overline{x}_{l+1}\right|\right]`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l2\overline{N}_1ϵ\left|x_0\overline{x}_0\right|`$
$``$ $`{\displaystyle \frac{2\overline{N}_1ϵ}{1\alpha }}\left|x_0\overline{x}_0\right|.`$
Let

$$l\left(ϵ\right)=\frac{2\overline{N}_1ϵ}{1\alpha }.$$

Notice that $`l\left(ϵ\right)0`$ as $`ϵ0`$. Thus,

$`\left|T\psi \left(x_0\right)T\psi \left(\overline{x}_0\right)\right|`$ $``$ $`l\left(ϵ\right)\left|x_0\overline{x}_0\right|`$
and so $`\left(T\psi \right)`$ satisfies the condition (4.2.72) for $`ϵ>0`$ sufficiently small.

Hence $`T`$ maps from $`𝕏𝕏.`$

Next we show T is a contraction on $`𝕏`$.

We express the solutions to the state dynamics

$$x_k=x(k,x_0,\psi )$$

and

$$\overline{x}_k=\overline{x}(k,\overline{x}_0,\psi )$$

in the implicit form,

$`x_k=\left(A+BK\right)^kx_0+{\displaystyle \underset{j=0}{\overset{k1}{}}}\left(A+BK\right)^{k1j}f_\psi (x_j,x_{j+1})`$
and

$`\overline{x}_k=\left(A+BK\right)^k\overline{x}_0+{\displaystyle \underset{j=0}{\overset{k1}{}}}\left(A+BK\right)^{k1j}f_\psi (\overline{x}_j,\overline{x}_{j+1}),`$
respectively.

We now denote $`x_j=x(j,x_0,\psi )`$ and $`\overline{x}_j=\overline{x}(j,x_0,\overline{\psi })`$ be the solutions to the state dynamics and satisfy the implicit form equations

$$x_k=\left(A+BK\right)^kx_0+\underset{j=0}{\overset{k1}{}}\left(A+BK\right)^{k1j}f_\psi (x_j,x_{j+1})$$

and

$$\overline{x}_k=\left(A+BK\right)^kx_0+\underset{j=0}{\overset{k1}{}}\left(A+BK\right)^{k1j}f_{\overline{\psi }}(\overline{x}_j,\overline{x}_{j+1}),$$

respectively.

The estimates (4.2.65)-(4.2.66) with the trajectories $`x(j,x_0,\psi )`$ and $`x(j,x_0,\overline{\psi })`$ becomes

$`\left|f_\psi (x,y)f_\psi (\stackrel{~}{x},\stackrel{~}{y})\right|`$ $``$ $`r_1\left(ϵ\right)\left|y\overline{y}\right|+r_2\left(ϵ\right)\psi \overline{\psi }+r_3\left(ϵ\right)\left|x\overline{x}\right|`$ (4.2.79)
$`\left|h_\psi (x,y)h_\psi (\stackrel{~}{x},\stackrel{~}{y})\right|`$ $``$ $`r_1\left(ϵ\right)\left|y\overline{y}\right|+r_2\left(ϵ\right)\psi \overline{\psi }+r_3\left(ϵ\right)\left|x\overline{x}\right|`$ (4.2.80)
where

$`r_1\left(ϵ\right)`$ $`=`$ $`n_{1,1}l\left(ϵ\right)+n_{1,2}ϵ+n_{1,3}l\left(ϵ\right)ϵ`$
$`r_2\left(ϵ\right)`$ $`=`$ $`n_{2,1}ϵ+n_{2,2}ϵ^2`$
$`r_3\left(ϵ\right)`$ $`=`$ $`n_3ϵ`$
and $`n_{i,j}`$ are positive constants. Observe that $`r_i\left(ϵ\right)0`$ as $`ϵ0`$.

At one-time step,

$`\left|x_{k+1}\overline{x}_{k+1}\right|`$ $``$ $`m_2\left(ϵ\right)\left|x_k\overline{x}_k\right|+m_3\left(ϵ\right)\psi \overline{\psi }`$
where

$$m_2\left(ϵ\right)=\frac{\alpha +r_3\left(ϵ\right)}{1r_1\left(ϵ\right)}$$

and

$$m_3\left(ϵ\right)=\frac{r_2\left(ϵ\right)}{1r_1\left(ϵ\right)}.$$

By invoking Gronwall’s inequality in finite difference form (4.1.1) and assuming that for some small $`ϵ>0`$

$$m_2\left(ϵ\right)<1,$$

then

$`\left|x_k\overline{x}_k\right|`$ $``$ $`m_3\left(ϵ\right)\left[{\displaystyle \underset{j=0}{\overset{k1}{}}}m_2\left(ϵ\right)^{k1j}\right]\psi \overline{\psi }`$ (4.2.81)
$``$ $`m_3\left(ϵ\right){\displaystyle \frac{1m_2\left(ϵ\right)^k}{1m_2\left(ϵ\right)}}\psi \overline{\psi }`$
$``$ $`{\displaystyle \frac{m_3\left(ϵ\right)}{1m_2\left(ϵ\right)}}\psi \overline{\psi }`$
With the bounds (4.2.80) and (4.2.81), we get

$`\left|T\psi \left(x_0\right)T\overline{\psi }\left(x_0\right)\right|`$ $``$ $`\left|{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\left(A^{}+K^{}B^{}\right)^l\left(h_\psi (x_l,x_{l+1})h_{\overline{\psi }}(\overline{x}_l,\overline{x}_{l+1})\right)\right|`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l\left|h_\psi (x_l,x_{l+1})h_{\overline{\psi }}(\overline{x}_l,\overline{x}_{l+1})\right|`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l\left[r_1\left(ϵ\right)\left|x_{l+1}\overline{x}_{l+1}\right|+r_2\left(ϵ\right)\psi \overline{\psi }+r_3\left(ϵ\right)\left|x_l\overline{x}_l\right|\right]`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\alpha ^l\left[2\left(r_1\left(ϵ\right)+r_3\left(ϵ\right)\right){\displaystyle \frac{m_3\left(ϵ\right)}{1m_2\left(ϵ\right)}}\psi \overline{\psi }+r_2\left(ϵ\right)\psi \overline{\psi }\right].`$
It follows that

$`\left|T\psi \left(x_0\right)T\overline{\psi }\left(\overline{x}_0\right)\right|`$ $``$ $`c\left(ϵ\right)\psi \overline{\psi }`$
where

$$c\left(ϵ\right)=\frac{2m_3\left(ϵ\right)\left(r_1\left(ϵ\right)+r_3\left(ϵ\right)\right)}{\left(1m_2\left(ϵ\right)\right)\left(1\alpha \right)}+\frac{r_2\left(ϵ\right)}{1\alpha }.$$

Notice that $`c\left(ϵ\right)0`$ as $`ϵ0`$. Thus, for sufficiently small $`ϵ>0`$

$$c\left(ϵ\right)<1.$$

Hence, $`T`$ is a contraction on $`𝕏`$ for $`ϵ`$ sufficiently small. Hence there exists a unique $`\psi 𝕏`$ such that

$$\psi =T\psi .$$

### 4.3 Some Properties

#### 4.3.1 Eigenstructure

Recall from the previous chapter the bidirectional linear Hamiltonian dynamics

$`\left[\begin{array}{c}x^+\\ \lambda \end{array}\right]`$ $`=`$ $`\left[\begin{array}{c}x\\ \lambda ^+\end{array}\right]`$ (4.3.86)
where

$`=\left[\begin{array}{cc}A& BR^1B^{}\\ Q& A^{}\end{array}\right].`$
###### Definition 4.3.1

Suppose

$`\left(\begin{array}{c}\delta x\\ \mu \delta \lambda \end{array}\right)=\left(\begin{array}{c}\mu \delta x\\ \delta \lambda \end{array}\right).`$ (4.3.92)
Then we call $`\mu `$ the eigenvalue of the dynamics (4.3.86).

In the case where we do not assume $`0`$ as the eigenvalue of $``$, we found that the forward Hamiltonian $`^F`$ is a hyperbolic system in Theorem 4.4.2. Similarly, we would like to show that the linear bidirectional Hamiltonian matrix $``$ is hyperbolic; i.e., the eigenvalue of $``$ lies strictly inside and outside the unit circle.

###### Theorem 4.3.2

If $`\mu `$ is an eigenvalue of the dynamics (4.3.86) satisfying the relation (4.3.92), then $`\frac{1}{\mu }`$ is also an eigenvalue of $``$.

Proof:  First, we decompose $``$ into

$`\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right]\left[\begin{array}{cc}Q& A^{}\\ A& BR^1B^{}\end{array}\right]=.`$
Let’s call

$`𝕊=\left[\begin{array}{cc}Q& A^{}\\ A& BR^1B^{}\end{array}\right]`$
and notice that $`𝕊`$ is symmetric.

From the equation (4.3.92), we have the following

$`\left(\begin{array}{cc}I& 0\\ 0& \mu I\end{array}\right)\left(\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right)=\left(\begin{array}{cc}\mu I& 0\\ 0& I\end{array}\right)\left(\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right)`$
which is equivalent to

$`\left(\begin{array}{cc}\frac{1}{\mu }I& 0\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ 0& \mu I\end{array}\right)\left(\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right)=\left(\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right).`$
We denote

$`_\mu =\left(\begin{array}{cc}\frac{1}{\mu }I& 0\\ 0& I\end{array}\right)\left(\begin{array}{cc}I& 0\\ 0& \mu I\end{array}\right).`$
Observe that $`1`$ is an eigenvalue of $`_\mu `$. It follows that

$$det\left[I_\mu \right]=0det\left[I_\mu ^{}\right]=0$$

where

$`_\mu ^{}=\left(\begin{array}{cc}I& 0\\ 0& \mu I\end{array}\right)𝕊\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\left(\begin{array}{cc}\frac{1}{\mu }I& 0\\ 0& I\end{array}\right).`$
Again, $`1`$ is an eigenvalue of $`_\mu `$; i.e.,

$`_\mu ^{}\left(\begin{array}{c}\stackrel{~}{\delta x}\\ \stackrel{~}{\delta \lambda }\end{array}\right)=\left(\begin{array}{c}\stackrel{~}{\delta x}\\ \stackrel{~}{\delta \lambda }\end{array}\right).`$
From above, we have

$`𝕊\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\left(\begin{array}{c}\frac{1}{\mu }\stackrel{~}{\delta x}\\ \stackrel{~}{\delta \lambda }\end{array}\right)=\left(\begin{array}{c}\stackrel{~}{\delta x}\\ \frac{1}{\mu }\stackrel{~}{\delta \lambda }\end{array}\right)`$
which is equal to

$`𝕊\left(\begin{array}{c}\stackrel{~}{\delta \lambda }\\ \frac{1}{\mu }\stackrel{~}{\delta x}\end{array}\right)=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\left(\begin{array}{c}\frac{1}{\mu }\stackrel{~}{\delta \lambda }\\ \stackrel{~}{\delta x}\end{array}\right).`$
Thus,

$`\left(\begin{array}{c}\stackrel{~}{\delta \lambda }\\ \frac{1}{\mu }\stackrel{~}{\delta x}\end{array}\right)=\left(\begin{array}{c}\frac{1}{\mu }\stackrel{~}{\delta \lambda }\\ \stackrel{~}{\delta x}\end{array}\right).`$
Hence, $`\frac{1}{\mu }`$ is an eigenvalue of $``$.   Note that $`\mu 1`$. The eigenvalue $`\mu =1`$ corresponds to the trivial eigenvector 0. Also, the infinite eigenvalues $`0`$ and $`\mathrm{}`$ satisfy (4.3.92) due to the singularity of A.

#### 4.3.2 Symplectic Form

The nonlinear tangent dynamics of the bidirectional Hamiltonian system (3.4.1) that we derive in Chapter 3 is

$`\left[\begin{array}{c}\delta x^+\\ \delta \lambda \end{array}\right]=_{\delta ,k}(x,\lambda ^+)\left[\begin{array}{c}\delta x\\ \delta \lambda ^+\end{array}\right]`$ (4.3.107)
where

$`_{\delta ,k}(x,\lambda ^+)=\left[\begin{array}{cc}H_{\lambda ^+x}& H_{\lambda ^+\lambda ^+}\\ H_{xx}& H_{x\lambda ^+}\end{array}\right](x,\lambda ^+).`$
We denote $`H_{x\lambda ^+},H_{xx},H_{\lambda ^+\lambda ^+}`$ as partial derivatives. Recall the nondegenerate and bilinear symplectic two-from $`\mathrm{\Omega }:T_{(x,\lambda )}\times T_{(x,\lambda )}`$,

$`\mathrm{\Omega }(v,w)=v^{}Jw\text{ and }J=\left[\begin{array}{cc}0& I\\ I& 0\end{array}\right]`$
where

$$\mathrm{\Omega }(v,w)=\mathrm{\Omega }(w,v)$$

and

$$v=\left[\begin{array}{c}\delta x\\ \delta \lambda \end{array}\right],w=\left[\begin{array}{c}\stackrel{~}{\delta x}\\ \stackrel{~}{\delta \lambda }\end{array}\right]$$

$`(x,\lambda ^+)`$ and $`(v,w)T_{(x,\lambda )}`$. Also, $`=T^{}𝒩`$ where $`x𝒩`$. We would like to show that under the tangent dynamics (4.3.107), the two-form $`\mathrm{\Omega }`$ is invariant; i.e.,

$`\mathrm{\Omega }(v,w)=\mathrm{\Omega }(v^+,w^+).`$ (4.3.110)
Then,

$`\mathrm{\Omega }(v,w)`$ $`=`$ $`v^{}Jw`$ (4.3.114)
$`=`$ $`\left(\begin{array}{cc}\delta x^{}& \delta x^{}H_{xx}+\delta \lambda ^+^{}H_{x\lambda ^+}\end{array}\right)J\left(\begin{array}{c}\stackrel{~}{\delta x}^{}\\ H_{xx}\stackrel{~}{\delta x}+H_{x\lambda ^+}\stackrel{~}{\delta \lambda }^+\end{array}\right)`$
$`=`$ $`\delta \lambda ^+^{}H_{x\lambda ^+}\stackrel{~}{\delta x}+\delta x^{}H_{x\lambda ^+}\stackrel{~}{\delta \lambda }^+`$ (4.3.115)
and

$`\mathrm{\Omega }(v^+,w^+)`$ $`=`$ $`v^+^{}Jw^+`$ (4.3.119)
$`=`$ $`\left(\begin{array}{cc}\delta x^{}H_{\lambda ^+x}+\delta \lambda ^+^{}H_{\lambda ^+\lambda ^+}^{}& \delta \lambda ^+\end{array}\right)J\left(\begin{array}{c}H_{\delta ^+x}\stackrel{~}{\delta x}+H_{\lambda ^+\lambda ^+}\stackrel{~}{\delta \lambda }^+\\ \stackrel{~}{\delta \lambda }^+\end{array}\right)`$
$`=`$ $`\delta \lambda ^+^{}H_{\lambda ^+x}\stackrel{~}{\delta x}+\delta x^{}H_{\lambda ^+x}\stackrel{~}{\delta \lambda }^+`$ (4.3.120)
Since (4.3.119) and (4.3.114) are equal, then

$$\mathrm{\Omega }(v,w)=\mathrm{\Omega }(v^+,w^+).$$

Thus, for any two tangent vectors satisfying the dynamics (4.3.107), the value of $`\mathrm{\Omega }`$ does not change.

### 4.4 Lagrangian Submanifold

We have two-form

$`\mathrm{\Omega }(v,w)=\delta \lambda ^+^{}H_{\lambda ^+x}\stackrel{~}{\delta x}+\delta x^{}H_{\lambda ^+x}\stackrel{~}{\delta \lambda }^+.`$ (4.4.121)
Using the technique in Chapter 3, we have the following tangent dynamics around the trajectories $`(x_j,x_{j+1})`$

$`\delta x_{j+1}`$ $`=`$ $`\left(\left(A+BK\right)+{\displaystyle \frac{f_\psi }{x_j}}(x_j,x_{j+1})\right)\delta x_j`$
$`+{\displaystyle \frac{f_\psi }{x_{j+1}}}(x_j,x_{j+1})\delta x_{j+1}`$
$`\delta \lambda _j`$ $`=`$ $`\left(A^{}+K^{}B^{}\right)\delta \lambda _{j+1}+{\displaystyle \frac{h_\psi }{x_{j+1}}}(x_j,x_{j+1})\delta x_{j+1}`$
By the Inverse Function Theorem, we can choose $`ϵ>0`$ small enough so that $`I\frac{f_\psi }{x_{k+1}}(x_k,x_{k+1})`$ is invertible for $`\left|x_k\right|,\left|x_{k+1}\right|<ϵ`$ because

$$I\frac{f_\psi }{x_{k+1}}(0,0)=I$$

since (4.2.47).

It follows that the tangent state dynamics is equivalent to

$`\delta x_{k+1}=\left(I{\displaystyle \frac{f_\psi }{x_{k+1}}}(x_k,x_{k+1})\right)^1\left(\left(A+BK\right)+{\displaystyle \frac{f_\psi }{x_k}}(x_k,x_{k+1})\right)\delta x_k`$
and

$`\delta x_{k+1}`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{k}{}}}(I{\displaystyle \frac{f_\psi }{x_{i+1}}}(x_i,x_{i+1}))^1`$
$`\left(\left(A+BK\right)+{\displaystyle \frac{f_\psi }{x_i}}(x_i,x_{i+1})\right)\delta x_0.`$
for $`\left|x_k\right|,\left|x_{k+1}\right|<ϵ`$.

As we let $`k\mathrm{}`$, we have $`x_k0`$, $`\frac{f_\psi }{(x,x^+)}\left(0\right)0`$. It follows that from (4.4)

$$\delta x_{k+1}\left(A+BK\right)^k\delta x_00.$$

Since $`\delta x_k0`$ as $`k\mathrm{}`$, we have that

$$\mathrm{\Omega }(v,w)0$$

for $`v,w`$ restricted to the tangent dynamics. Thus, the local stable manifold is a Lagrangian submanifold. Similarly as in Chapter 3, we have

$`{\displaystyle \frac{\psi _i}{x_j}}{\displaystyle \frac{\psi _j}{x_i}}=0\text{ for }i,j=1,\mathrm{},n.`$ (4.4.124)
The equation (4.4.124) implies that $`\psi \left(x\right)`$ is closed. Then, by the *Stokes’ Theorem* there exists $`\overline{\pi }𝒞^r\left(^n\right)`$ such that

$`\psi \left(x\right)={\displaystyle \frac{\overline{\pi }}{x}}\left(x\right)\text{ where }\psi \left(0\right)=0`$ (4.4.125)
locally on some neighborhood of 0. Since $`z=\psi \left(x\right)=\lambda Px`$, we have that

$`\lambda `$ $`=`$ $`{\displaystyle \frac{\overline{\pi }}{x}}\left(x\right)+Px`$
$`=`$ $`{\displaystyle \frac{\overline{\pi }}{x}}\left(x\right)+{\displaystyle \frac{}{x}}\left({\displaystyle \frac{1}{2}}x^{}Px\right)`$
because $`P`$ is symmetric. Thus, there exists $`\pi 𝒞^r`$ such that

$$\lambda =\frac{\pi }{x}\left(x\right)$$

where

$$\frac{\pi }{x}\left(x\right)=\frac{\overline{\pi }}{x}\left(x\right)+\frac{}{x}\left(\frac{1}{2}x^{}Px\right).$$

Hence, the local stable manifold $`\lambda `$ is the gradient of the optimal cost for the bidirectional Hamiltonian dynamics.

## Chapter 5 Solution of the Hamilton Jacobi Bellman Equations

Our procedure for solving HJB PDE numerically is described in , but we give a summary of the method and some numerical results in this chapter.

### 5.1 A Method for HJB

We solve the Hamilton Jacobi Bellman (HJB) Partial Differential Equation that arises in many control problems. We consider the infinite horizon optimal control problem of minimizing the cost

$`{\displaystyle _t^{\mathrm{}}}l(x,u)𝑑t`$ (5.1.1)
subject to the dynamics

$`\dot{x}=f(x,u)`$ (5.1.2)
and initial condition

$`x\left(t\right)=x^0.`$ (5.1.3)
The state vector $`x`$ is an $`n`$ dimensional column vector, the control $`u`$ is an $`m`$ dimensional column vector and the dynamics $`f(x,u)`$ and Lagrangian $`l(x,u)`$ are assumed to be sufficiently smooth.

If the minimum exists and is a smooth function $`\pi \left(x^0\right)`$ of the initial condition then it satisfies the HJB PDE

$`\underset{u}{\mathrm{min}}\left\{{\displaystyle \frac{\pi }{x}}\left(x\right)f(x,u)+l(x,u)\right\}=0`$ (5.1.4)
and the optimal control $`\kappa \left(x\right)`$ satisfies

$`\kappa \left(x\right)=\text{arg}\underset{u}{\mathrm{min}}\left\{{\displaystyle \frac{\pi }{x}}\left(x\right)f(x,u)+l(x,u)\right\}=0`$ (5.1.5)
These are expressed in terms of the Hamiltonian

$`H(p,x,u)=pf(x,u)+l(x,u)`$ (5.1.6)
where the argument $`p`$ is an $`n`$ dimensional row vector. The HJB PDE becomes

$`0`$ $`=`$ $`\underset{u}{\mathrm{min}}H({\displaystyle \frac{\pi }{x}}\left(x\right),x,u)`$ (5.1.7)
$`\kappa \left(x\right)`$ $`=`$ $`\text{arg}\underset{u}{\mathrm{min}}H({\displaystyle \frac{\pi }{x}}\left(x\right),x,u)`$ (5.1.8)
We assume that the Hamiltonian $`H(p,x,u)`$ is strictly convex in $`u`$ for all $`p,x`$. Then (5.1.4, 5.1.5) become

$`{\displaystyle \frac{\pi }{x}}\left(x\right)f(x,\kappa \left(x\right))+l(x,\kappa \left(x\right))=0`$ (5.1.9)
and

$`{\displaystyle \frac{\pi }{x}}\left(x\right){\displaystyle \frac{f}{u}}(x,\kappa \left(x\right))+{\displaystyle \frac{l}{u}}(x,\kappa \left(x\right))=0`$ (5.1.10)
These are the equations (5.1.9),(5.1.10) that we solve for the optimal cost $`\pi `$ and optimal control $`\kappa `$.

#### 5.1.1 Al’brecht’s Method Revisited

Al’brecht  solved the HJB PDE locally around zero by expanding the problem in a power series,

$`f(x,u)`$ $`=`$ $`Ax+Bu+f^{\left[2\right]}(x,u)`$ (5.1.11)
$`+f^{\left[3\right]}(x,u)+\mathrm{}`$
$`l(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(x^{}Qx+2x^{}Su+u^{}Ru\right)`$ (5.1.12)
$`+l^{\left[3\right]}(x,u)+l^{\left[4\right]}(x,u)+\mathrm{}`$
$`\pi \left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}x^{}Px+\pi ^{\left[3\right]}\left(x\right)`$ (5.1.13)
$`+\pi ^{\left[4\right]}\left(x\right)+\mathrm{}`$
$`\kappa \left(x\right)`$ $`=`$ $`Kx+\kappa ^{\left[2\right]}\left(x\right)+\kappa ^{\left[3\right]}\left(x\right)+\mathrm{}`$ (5.1.14)
where $`^{\left[d\right]}`$ denotes a homogeneous polynomial of degree $`d`$. As in Chapter 3, we subsitute these equations into the HJB PDE (5.1.9, 5.1.10) and equate the terms of like degree to obtain a sequence of algebraic equations for the unknowns.

The first level is the pair of equations obtained by collecting the quadratic terms of (5.1.9) and the linear terms of (5.1.10). We denote the $`d^{th}`$ level as the pair of equations garnered from the $`\left[d+1\right]^{th}`$ degree terms of (5.1.9) and the $`d^{th}`$ degree terms of (5.1.10).

Here are the first set of equations:

$`0`$ $`=`$ $`A^{}P+PA+Q`$ (5.1.15)
$`\left(PB+S\right)R^1\left(PB+S\right)^{}`$
$`K`$ $`=`$ $`R^1\left(PB+S\right)^{}`$ (5.1.16)
The quadratic terms of the HJB PDE reduce to the familiar Riccati equation (5.1.15) and the linear optimal feedback (5.1.16).

We assume $`A,B`$ is stabilizable and $`Q,A`$ is detectable then the Riccati equation has a unique positive definite solution $`P`$ and the linear feedback locally exponentially stabilizes the closed loop system. Moreover the optimal quadratic cost is a local Lyapunov function for the closed loop system.

The $`d^{th}`$ Level

Suppose we have solved through the $`d1^{th}`$ level after repeating the process $`d1`$ times. It is convenient to incorporate this solution into the dynamics and the cost. Let $`\kappa ^{k]}\left(x\right)=Kx+\kappa ^{\left[2\right]}\left(x\right)+\kappa ^{\left[3\right]}\left(x\right)+\mathrm{}+\kappa ^{\left[k\right]}\left(x\right)`$ and define

$`\overline{f}(x,u)`$ $`=`$ $`f(x,\kappa ^{d1]}\left(x\right)+u)`$ (5.1.17)
$`\overline{l}(x,u)`$ $`=`$ $`l(x,\kappa ^{d1]}\left(x\right)+u)`$ (5.1.18)
These have power series expansions through terms of degree $`d`$ and $`d+1`$ of the form

$`\overline{f}(x,u)`$ $`=`$ $`\left(A+BK\right)x+Bu+\overline{f}^{\left[2\right]}(x,u)+\mathrm{}\overline{f}^{\left[d\right]}(x,u)+\mathrm{}`$
$`\overline{l}(x,u)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(x^{}Qx+2x^{}SKx+x^{}K^{}RKx\right)`$
$`+x^{}Su+u^{}Ru+`$
$`\overline{l}^{\left[3\right]}(x,u)+\mathrm{}+\overline{l}^{\left[d+1\right]}(x,u)+\mathrm{}`$
We plug these into the HJB PDE (5.1.9, 5.1.10) and find that they are satisfied through the $`d1`$ level and don’t involve $`u`$. Since $`u=\kappa ^{\left[d\right]}\left(x\right)+\mathrm{}`$, the $`d`$ level equations are

$`0`$ $`=`$ $`{\displaystyle \frac{\pi ^{\left[d+1\right]}}{x}}\left(x\right)\left(A+BK\right)x`$
$`+{\displaystyle \underset{i=2}{\overset{d1}{}}}{\displaystyle \frac{\pi ^{\left[d+2i\right]}}{x}}\left(x\right)\overline{f}^{\left[i\right]}(x,0)+x^{}PBu`$
$`+\overline{l}^{\left[d+1\right]}(x,0)+x^{}Su+{\displaystyle \frac{1}{2}}x^{}K^{}Ru`$
and

$`0`$ $`=`$ $`{\displaystyle \frac{\pi ^{\left[d+1\right]}}{x}}\left(x\right)B+{\displaystyle \underset{i=2}{\overset{d}{}}}{\displaystyle \frac{\pi ^{\left[d+2i\right]}}{x}}{\displaystyle \frac{\overline{f}^{\left[i\right]}}{u}}(x,0)`$ (5.1.20)
$`+{\displaystyle \frac{\overline{l}^{\left[d+1\right]}}{u}}(x,0)+u^{}R`$
Because of (5.1.16), $`u`$ drops out of the first equation which becomes

$`0`$ $`=`$ $`{\displaystyle \frac{\pi ^{\left[d+1\right]}}{x}}\left(x\right)\left(A+BK\right)x`$
$`+{\displaystyle \underset{i=2}{\overset{d1}{}}}{\displaystyle \frac{\pi ^{\left[d+2i\right]}}{x}}\left(x\right)\overline{f}^{\left[i\right]}(x,0)`$
$`+\overline{l}^{\left[d+1\right]}(x,0).`$
After this has been solved for $`\pi ^{\left[d+1\right]}\left(x\right)`$, we can solve the second for $`\kappa ^{\left[d\right]}\left(x\right)`$,

$`\kappa ^{\left[d\right]}\left(x\right)`$ $`=`$ $`R^1({\displaystyle \frac{\pi ^{\left[d+1\right]}}{x}}\left(x\right)B`$
$`+{\displaystyle \underset{i=2}{\overset{d}{}}}{\displaystyle \frac{\pi ^{\left[d+2i\right]}}{x}}{\displaystyle \frac{\overline{f}^{\left[i\right]}}{u}}(x,0)+{\displaystyle \frac{\overline{l}^{\left[d+1\right]}}{u}}(x,0))`$
These equations admit a unique solution up to the smoothness of $`f`$ and $`l`$ since $`A+BK`$ has all its eigenvalues in the left half plane. If $`f`$ and $`l`$ are real analytic, the power series converges to the solution of the HJB PDE locally around $`x=0`$ , . The higher degree equations are linear; hence they are easily solvable.

Albrecht’s approach develops the Taylor series expansion along a point. The Taylor series only converges on a small neigbhorhood and may quickly diverges outside that neighborhood. Also increasing the degree of the polynomial solutions does not guarantee a larger region where the approximated solutions are true. In our approach, we use Albrecht’s method to generate the initial approximations around $`0`$. Then the initial approximations are then improve by piecing successive approximations of smooth solutions around the points on the level sets, see Fig.(5.1). Instead we truncate the computed solutions at the point where the optimal cost satisfies the optimality and stability constraints and begin new approximations at the same point.

#### 5.1.2 The Lyapunov Criterion

Let $`\pi ^{d+1]}`$ and $`\kappa ^{d]}`$ denote the approximate cost and feedback to degrees $`d+1`$ and $`d`$ respectively. We assume $`l(x,u)0`$. We find the largest sublevel set

$`\pi ^{d+1]}\left(x\right)c`$ (5.1.23)
such that

$`{\displaystyle \frac{\pi ^{d+1]}}{x}}\left(x\right)f(x,\kappa ^{d]}\left(x\right))\left(1ϵ_1\right)l(x,\kappa ^{\left[d\right]}\left(x\right))`$ (5.1.24)
$`\left(1ϵ_2\right)l(x,\kappa ^{\left[d\right]}\left(x\right)){\displaystyle \frac{\pi ^{d+1]}}{x}}\left(x\right)f(x,\kappa ^{d]}\left(x\right))`$ (5.1.25)
The parameter $`ϵ_1`$ controls the rate of exponential stability of the closed loop system while $`ϵ_2`$ dictates the rate of the optimality of the feedback. We determine a sublevel set of the approximate cost on which it is an acceptable Lyapunov function and the approximate feedback is stabilizing and satisfying the optimality conditions. We emphasize the stabilizing property of the control law rather than its optimality because usually optimality is only a tool to find a stabilizing feedback. Typically the goal is to find a control law to stabilize the system and the optimal control problem is formulated as a way of finding one.

#### 5.1.3 Power Series Expansions

We specialize to problems with dynamics that is affine in $`u`$ and with a Lagrangian that is quadratic in $`u`$,

$`f(x,u)`$ $`=`$ $`g_0\left(x\right)+g_1\left(x\right)u`$ (5.1.26)
$`l(x,u)`$ $`=`$ $`l_0\left(x\right)+l_1\left(x\right)u+u^{}l_2\left(x\right)u`$ (5.1.27)
where $`l_2\left(x\right)`$ is an invertible $`m\times m`$ matrix for all $`x`$.

Assume that we have solved the HJB PDE locally around 0 and have choosen a sublevel set of value $`c`$ subject to the condition (5.1.23, 5.1.24). We generate a power series solution around a point $`\overline{x}`$ on the level set $`\pi ^{d+1]}\left(x\right)=c`$.

We introduce some notation. Let $`\alpha =(\alpha _1,\mathrm{},\alpha _n)`$ be a multi-index of nonnegative integers and $`\left|\alpha \right|=_i\alpha _i`$. Let $`\beta =(\beta _1,\mathrm{},\beta _n)`$. We say $`\beta \alpha `$ if $`\beta _i\alpha _i,i=1,\mathrm{},n`$ and $`\beta <\alpha `$ if $`\beta \alpha `$ and for at least one $`i`$, $`\beta _i<\alpha _i`$. Let $`\mathrm{𝟎}=(0,\mathrm{},0)`$.

Define the differential operator

$$D^\alpha =\left(\frac{}{x_1}\right)^{\alpha _1}\mathrm{}\left(\frac{}{x_n}\right)^{\alpha _n},$$

the multifactorial

$$\alpha !=\alpha _1!\mathrm{}\alpha _n!,$$

the monomial

$$x^\alpha =x_1^{\alpha _1}\mathrm{}x_n^{\alpha _n},$$

and the coefficient

$$C(\alpha ,\beta )=\left(\begin{array}{c}\alpha _1\\ \beta _1\end{array}\right)\mathrm{}\left(\begin{array}{c}\alpha _n\\ \beta _n\end{array}\right)$$

where $`\beta \alpha `$.

We derive a system of equations for

$$D^\alpha \pi \left(\overline{x}\right),D^\alpha \kappa \left(\overline{x}\right).$$

We already know that if $`\alpha `$ is the $`i^{th}`$ unit vector

$`D^\mathrm{𝟎}\pi \left(\overline{x}\right)`$ $`=`$ $`\pi \left(\overline{x}\right)`$
$`D^\mathrm{𝟎}\kappa \left(\overline{x}\right)`$ $`=`$ $`\kappa \left(\overline{x}\right)`$
$`D^\alpha \pi \left(\overline{x}\right)`$ $`=`$ $`{\displaystyle \frac{\pi }{x_i}}\left(\overline{x}\right).`$
Let the Cauchy data $`D^\alpha \pi \left(\overline{x}\right)`$ and $`D^\alpha \kappa \left(\overline{x}\right)`$ be known for $`\alpha =(0,\alpha _2,\alpha _3,\mathrm{},\alpha _n)`$ and $`f_1(\overline{x},\kappa \left(\overline{x}\right))0`$.

Assume that we have derived algebraic equations $`D^\beta \pi \left(\overline{x}\right)`$ for $`D^\beta \pi \left(\overline{x}\right)`$ and $`D^\beta \kappa \left(\overline{x}\right)`$ for $`\beta <\alpha `$. We apply $`D^\alpha `$ to (5.1.9) to obtain

$`0`$ $`=`$ $`{\displaystyle \frac{}{x}}\left(D^\alpha \pi \left(\overline{x}\right)\right)f(\overline{x},\kappa \left(\overline{x}\right))`$
$`+{\displaystyle \underset{\mathrm{𝟎}<\beta \alpha }{}}C(\alpha ,\beta )\left(D^{\alpha \beta }{\displaystyle \frac{\pi }{x}}\left(\overline{x}\right)\right)D^\beta f(\overline{x},\kappa \left(\overline{x}\right))`$
$`+D^\alpha l(\overline{x},\kappa \left(\overline{x}\right)).`$
This yields an equation for $`\frac{D^\alpha \pi }{x_1}\left(\overline{x}\right)`$ because $`f_1(\overline{x},\kappa \left(\overline{x}\right))0`$ and all the other terms $`\frac{D^\alpha \pi }{x_i}\left(\overline{x}\right)`$ for $`i1`$ are known from the Cauchy data.

We apply $`D^\alpha `$ to (5.1.10) to obtain

$`0`$ $`=`$ $`{\displaystyle \underset{\mathrm{𝟎}\beta \alpha }{}}C(\alpha ,\beta )\left(D^{\alpha \beta }{\displaystyle \frac{\pi }{x}}\left(\overline{x}\right)\right)D^\beta g_1\left(\overline{x}\right)`$ (5.1.29)
$`+D^\alpha l_1\left(\overline{x}\right)+\left(D^\alpha \kappa \left(\overline{x}\right)\right)^{}l_2\left(\overline{x}\right)`$
$`{\displaystyle \underset{\mathrm{𝟎}<\beta \alpha }{}}C(\alpha ,\beta )\left(D^{\alpha \beta }\kappa \left(\overline{x}\right)\right)^{}D^\beta l_2\left(\overline{x}\right)`$
Notice that this equation (5.1.29) only contains $`D^\alpha \kappa \left(x\right)`$ in one term multiplied by an invertible matrix so we can express $`D^\alpha \kappa \left(x\right)`$ as a function of $`D^\beta \kappa \left(x\right)`$ for $`\beta <\alpha `$ and $`D^\gamma \pi \left(x\right)`$.

In this way we obtain the approximations

$`\pi \left(x\right)`$ $``$ $`{\displaystyle \underset{d}{}}{\displaystyle \underset{\left|\alpha \right|d}{}}{\displaystyle \frac{1}{\alpha !}}D^\alpha \pi \left(\overline{x}\right)(x\overline{x}))^\alpha `$ (5.1.30)
$`\kappa \left(x\right)`$ $``$ $`{\displaystyle \underset{d}{}}{\displaystyle \underset{\left|\alpha \right|d}{}}{\displaystyle \frac{1}{\alpha !}}D^\alpha \kappa \left(\overline{x}\right)\left(x\overline{x}\right)^\alpha `$ (5.1.31)
where $`x`$ is close to $`\overline{x}`$.

### 5.2 An Example

Consider the optimal control problem:

$$\underset{u}{\mathrm{min}}_0^{\mathrm{}}\mathrm{ln}^2\left(x+1\right)+u^2dt$$

subject to (5.2.32)
$`\dot{x}`$ $`=`$ $`xu+u`$
$`x\left(0\right)`$ $`=`$ $`x_0`$
where the domain $`𝒟=\{x|x>1)\}`$. Given the problem (5.2.32) and the schemes described above, we solve the optimal control $`\kappa \left(x\right)`$ and optimal cost $`\pi \left(x\right)`$ up to degree $`d`$ and $`d+1`$, respectively; i.e.

$`\pi \left(x\right)`$ $`=`$ $`\pi ^{\left[2\right]}\left(x\right)+\pi ^{\left[3\right]}\left(x\right)+\mathrm{}+\pi ^{\left[d+1\right]}\left(x\right)`$
$`\kappa \left(x\right)`$ $`=`$ $`\kappa ^{\left[1\right]}\left(x\right)+\kappa ^{\left[2\right]}\left(x\right)+\mathrm{}+\kappa ^{\left[d\right]}\left(x\right)`$
We fix the degree at $`d=3`$ for this example. The interval is $`𝒟=(1,4]`$. Note that the analytic solutions are $`\kappa _{}\left(x\right)=\mathrm{ln}\left(x+1\right)`$ and $`\pi _{}\left(x\right)=\mathrm{ln}^2\left(x+1\right)`$.

#### 5.2.1 Approximation I: Albrecht’s Method

We implement Albrecht’s method by using the MATLAB code hjb.m in  to find the coefficients of $`\pi `$ and $`\kappa `$. We then set-up the polynomials. At each point $`x_j𝒟`$, we assign the polynomial approximations $`\pi _{0_j}=\pi _0\left(x_j\right)`$ and $`\kappa _{0_j}=\kappa _0\left(x_j\right)`$. See Fig.(5.2). We solve the optimization problem (5.1.23). We use the MATLAB code, fmincon.m, iteratively to find a point on the level set $`\pi ^{d+1]}\left(x\right)<c`$; i.e., $`\overline{x}𝒟_l`$. For this example, we march along the $`x`$-axis in the both directions. For the interval $`[0,4]`$, we march on the x-axis towards $`\mathrm{}`$, while on interval $`(1,0]`$ we move towards -$`1`$. See Fig.(5.3) for the psuedocode of the algorithm DRIVER1.m and the actual codes in the appendix.

#### 5.2.2 Approximation II

The scheme in the second part is to improve the smooth solutions of Albrecht. See Fig.(5.4). The analogue of hjb.m is kovalesky.m in DRIVER2.m, the codes generating the coefficients of the polynomials. The coefficients of the polynomials are the derivatives and higher order derivatives of $`\pi `$ and $`\kappa `$ evaluated at $`\overline{x}`$. First, we solve the equations derived from the problem (5.2.32) that is described in Section $`\mathrm{5.1.3}`$. We obtain the equations for the derivatives and its higher order derivatives by using MAPLE. Consequently, we write these equations in MATLAB to get the coefficients evaluated at $`\overline{x}`$. Then, we piece together the old and the new solutions $`\pi _{l1}`$ and $`\pi _l`$ and $`\kappa _{l1}`$ and $`\kappa _l`$ on the interval $`𝒟`$. We denote $`\pi _l`$ as the $`l`$th approximated solution and $`\pi _{new}`$ as the union of the truncated $`\pi _l`$ for $`l=0,1,\mathrm{}`$ on the interval $`[\overline{x}_l,\overline{x}_{l+1}]`$. We see that solutions overlapped on some interval around $`\overline{x}`$; i.e. $`\pi _{l1}\left(x\right)=\pi _l\left(x\right)`$ for $`x[\overline{x}_l\epsilon ,\overline{x}_l]`$ in the case $`x\mathrm{}`$. In this example, the solutions $`\pi _{l1}`$ and $`\pi _l`$ and $`\kappa _{l1}`$ and $`\kappa _l`$ coincide on some small interval around $`\overline{x}`$. However, we do not expect the same result for other examples or in general. In these cases, we take

$$\pi _l\left(x\right)=\underset{x}{\mathrm{min}}\{\pi _{l1}\left(x\right),\pi _l\left(x\right)\}$$

for some point $`x`$ on some interval $`𝒟`$ since the lower approximation is the optimal one. The process is then repeated at the next $`\overline{x}_{l+1}`$. In Fig.(5.5), we compare the new approximations, $`\pi _{new}`$ and $`\kappa _{new}`$, with the real solutions, $`\pi _{}`$ and $`\kappa _{}`$ and Albrecht’s solutions, $`\pi _0`$ and $`\kappa _0`$. The new polynomial, $`\pi _{new}`$, is the union of $`\pi _0`$ on $`[0.6094,0]`$ and $`\pi _1`$ on $`[1,0.6094]`$. Similarly, $`\kappa _{new}`$ is the outcome of glueing $`\kappa _0`$ and $`\kappa _1`$. See Figs.(5.6) and (5.7)

## Appendix A Codes

### A.1 Driver 1

driver1.m

```
% generates polynomials around 0 via albrecht
%
% Prager Example
%    Dynamics:        xdot=xu + u; x(0)=x0
%    Cost:            (ln x+1)^2 + u^2
%    Taylorized Cost: x^2 -x^3+(11/12)x^4-...+u^2
%-------------------------------------------------------------------
% Dimension Guide
% d=1,n=1,m=1: input f(1,2), l(1,3) (linear, starts with quadratic)
%              output ka=(1,1), py=(1,1) (linear,st w quadratic)
% d=2,n=1,m=1: input f(1,5), l(1,7) (up to 2nd, up to 3rd)
%              output ka=(1,2), py=(1,2) (up to 2nd, up to 3rd)
% d=3,n=1,m=1: input f(1,9), l(1,12) (up to 3rd, up to 4th)
%              output ka=(1,3), py=(1,3) (up to 3rd, up to 4th)
% d=4,n=1,m=1: input f(1,14), l(1,18) (up to 4th, up to 5th)
%              output ka=(1,4), py=(1,4) (up to 4th, up to 5th)
%--------------------------------------------------------------------
%Subroutines:
%           poly1.m
%           endpoint.m
%           theG.m
%           theF.m
%Find coeffients of Taylor’’s Expansion

%d -- degree up to
%xbar -- Taylor expand around xbar

function driver1(d)

mesh=256;

[x1,pie,U,dom1,dom2,rrealu,rrealpy]=poly1(d);
[root]=endpoint(pie);

subplot(2,1,1),plot(x1,pie,’m’);
hold;
plot(x1,rrealpy,’c--’);
legend(’\pi^0(x)’,’\pi^{*}(x)’);
hold;
title(’\pi(x) of degree d+1=4’);
xlabel(’x’);
ylabel(’\pi’);
axis([0 4 -2 20]);

subplot(2,1,2),plot(x1,U);
hold;
plot(x1,rrealu,’c--’);
legend(’\kappa^0(x)’,’\kappa^{*}(x)’);
hold;
title(’\kappa(x) of degree d=3’);
xlabel(’x’);
ylabel(’\kappa’);
axis([0 4 -10 10]);

```

Subroutines

poly1.m

```
% make polynomial functions
function [x1,pie,U,dom1,dom2,rrealu,rrealpy]=poly1(d)

n=1;
if (d==1)
  f=[0 1];
  l=[1 0 1];
elseif (d==2)
  f=[0 1 0 1 0];
  l=[1 0 1 -1 0 0 0];
elseif (d==3)
  f=[0 1 0 1 0 zeros(1,4)];
  l=[1 0 1 -1 0 0 0 11/12  zeros(1,4) ];
else
  f=[0 1 0 1 0 zeros(1,9)];
  l=[ 1 0 1 -1 0 0 0 11/12 zeros(1,4) -5/6 zeros(1,5)];
  % l=[ 1 0 1 -1 0 0 0 -1 zeros(1,10) ];
end% ifloop

[ka,fk,py,lk]= hjb(f,l,1,1,d);

%generating vectors X, Y (py & ka)
xslots=nchoosek(n+(1+1)-1,1+1);
yslots=nchoosek(n+(1)-1,1);
X=zeros(xslots,1);
Y=zeros(yslots,1);

if (d>=2)
  for i=2:d
   xslots=nchoosek(n+(d+1)-1,d+1);
   yslots=nchoosek(n+(d)-1,d);
   X=[[X] zeros(xslots,1)];
   Y=[[Y] zeros(yslots,1)];
  end% dloop
end %ifloop
X=X’;
Y=Y’;
[xsize,dummy]=size(X);
[ysize,dummy]=size(Y);

aa=0;
bb=4;
mesh=256;
dx=(bb-aa)/mesh;
dy=dx;

for i=1:mesh
  x1(i)= aa + i*dx;
  a=x1(i);
  X=[a^2; a^3; a^4];
  pie(i)=py*X;
  PY=pie(i);
  realpy(i)=x1(i)^2 -x1(i)^3 +(11/12)*x1(i)^4;
  rrealpy(i)=(log(x1(i)+1))^2;
  P=py*X;
  Y=[a; a^2; a^3];
  U(i)=ka*Y;
  realu(i)=-(x1(i)-(1/2)*x1(i)^2+(1/3)*x1(i)^3);
  rrealu(i)=-log(1+x1(i));
  FK=pragerf(x1(i),U(i));
  LK=pragercost(x1(i),U(i));
  gradP=0;

for k=1:d
 gradP=gradP + (k+1)*py(k)*x1(i)^(k);
 gP(i)=gradP;
end %gradPloop

end %for loop

```

endpoint.m

```
%finds pt on levelset
function [root]=endpoint(pie)

C=fmincon(’theG’,2,0,0,0,0,0,4,’theF’);

%finding root
root=100;
m=100;
mesh=256;

for i=1:mesh
     if (abs(C-pie(i))==0)
     root=i;
     end
end %for loop

```

theG.m

```
[ka,py]=justgivemepk;
G=-py*[x^2; x^3; x^4];

```

theF.m

```
%contraint func for fmincon
function [F,eqF]=theF(x)

[ka,py]=justgivemepk;
u=ka*[x; x^2; x^3];
FK=pragerf(x,u);
LK=pragercost(x,u);
gradP=0;
d=3;

for k=1:d
 gradP=gradP + (k+1)*py(k)*x^(k);
end %gradPloop

  eps=2^(-6);
  F1=gradP*FK+(1-eps)*LK;
  F2=-(gradP*FK+(1-eps)*LK);
  F=[F1;F2];
  eqF=0;

```

justgivepk.m

```
%output coefficients
function [ka,py]=justgivemepk

d=3;
n=1;
if (d==1)
  f=[0 1];
  l=[1 0 1];
elseif (d==2)
  f=[0 1 0 1 0];
  l=[1 0 1 -1 0 0 0];
elseif (d==3)
  f=[0 1 0 1 0 zeros(1,4)];
  l=[1 0 1 -1 0 0 0 11/12  zeros(1,4) ];
else
%d=4
  f=[0 1 0 1 0 zeros(1,9)];
  l=[ 1 0 1 -1 0 0 0 11/12 zeros(1,4) -5/6 zeros(1,5)];
end% ifloop

%hjb(f,l,n,m,d,f_,n_,m_)
[ka,fk,py,lk]= hjb(f,l,1,1,d);

```

### A.2 driver 2

Driver2.m

```
% generate polynomial around the x_0
% improves albrecht approx
%
% Prager Example
%    Dynamics:        xdot=xu + u; x(0)=x0
%    Cost:            (ln x+1)^2 + u^2
%    Taylorized Cost: x^2 -x^3+(7/12)x^4-...+u^2
%-------------------------------------------------------------------
% Dimension Guide
% d=1,n=1,m=1: input f(1,2), l(1,3) (linear, starts with quadratic)
%              output ka=(1,1), py=(1,1) (linear,st w quadratic)
% d=2,n=1,m=1: input f(1,5), l(1,7) (up to 2nd, up to 3rd)
%              output ka=(1,2), py=(1,2) (up to 2nd, up to 3rd)
% d=3,n=1,m=1: input f(1,9), l(1,12) (up to 3rd, up to 4th)
%              output ka=(1,3), py=(1,3) (up to 3rd, up to 4th)
% d=4,n=1,m=1: input f(1,14), l(1,18) (up to 4th, up to 5th)
%              output ka=(1,4), py=(1,4) (up to 4th, up to 5th)
%--------------------------------------------------------------------
%Subroutines:
%           kovalesky.m
%           poly2.m
%           glue.m
%
%d -- degree up to
%xbar -- Taylor expand around xbar

function driver2(d)

mesh=256;
[x1,pie,U,gP]=poly1(d);

for j=1:4

root=endpoint(pie);
ii=root;
[ka,py]=kovalesky(d,xbar,k1,pi1,pi2);
[x1,pie2,U2,rrealu,rrealpy,gP]=poly2(d,ka,py,xbar,x1);
[newp,newu]=glue(pie2,U2,pie,U,istar,x1,rrealu,rrealpy);

figure;
subplot(2,1,1);
hold;
plot(x1,rrealpy,’c--’);
plot(x1,newp,’k’);
hold;
legend(’\pi_{*}(x)’,’\pi_{new}(x)’);
title(’\pi(x) of degree d+1=4’);
xlabel(’x’);
ylabel(’\pi’);
axis([0 4 -2 20]);

subplot(2,1,2);
hold;
plot(x1,rrealu,’c--’);
plot(x1,newu,’k’);
hold;
legend(’\kappa_{*}(x)’,’\kappa_{new}(x)’);
title(’\kappa(x) of degree d=3’);
xlabel(’x’);
ylabel(’\kappa’);
axis([0 4 -10 10]);

pie=newp;
U=newu;

end %for

```

Subroutines

kovalesky.m

```
%generates coefficients
function [ka,py]=kovalesky(d,xbar,k1,pi1,pi2)

%py -- vector containing Taylor coefficients of py
%ka -- vector containing Taylor coefficients of ka

ka=zeros(1,d+1);
py=zeros(1,d+2);
%degree=d;

     if (d >= 1)
        [ka1,py1]=Coeffd1(xbar,ka,py,d,k1,pi1,pi2);
     end% if loop

     if (d >= 2)
     [ka2,py2]=Coeffd2(xbar,ka1,py1);
     end% if loop

     if (d >= 3)
        [ka,py]=Coeffd3(xbar,ka2,py2);
     end% if loop

```

poly2.m

```
%set up polynomial functions
function [x1,pie,U,rrealu,rrealpy,gP]=poly2(d,ka,py,xbar,x1)

mesh=256;

for i=1:mesh
  a=x1(i);
  X=[1; (a-xbar);(1/factorial(2))*(a-xbar)^2;(1/factorial(3))*(a-xbar)^3;(1/factorial(4))*(a-xbar)^4];
  pie(i)=py*X;
  PY=pie(i);
  rrealpy(i)=(log(x1(i)+1))^2;
  Y=[1; (a-xbar);(1/factorial(2))*(a-xbar)^2;(1/factorial(3))*(a-xbar)^3];
  U(i)=ka*Y;
  rrealu(i)=-log(1+x1(i));
  FK=pragerf(x1(i),U(i));
  LK=pragercost(x1(i),U(i));
  gradP=0;

for k=1:d+1
   gradP=gradP + (k)*(1/factorial(k))*py(k+1)*(x1(i)-xbar)^(k-1);
   gP(i)=gradP;
end %gradPloop

end %for loop i

```

glue.m

```
% attach new polynomial in appropriate interval
function [newp,newu]=glue(p,u,pie,U,istar,x1,rrealu,rrealpy)

mesh=256;
[dum n]=size(p);
ii=istar;

for i=1:ii-1
  newp(i)=pie(i);
  newu(i)=U(i);
end

newp(ii:n)=p(ii:n);
newu(ii:n)=u(ii:n);

```

coeffd1.m

```
%Calculate the deg=1 coefficients of Prager’’s Example
function [k,p]=Coeffd1(a,k,p,deg,k1,pi1,pi2)

[PP,KK]=Tconst(a,deg);
p(1)=pi1;
k(1)=k1;
p(2)=pi2;
l=(log(a+1))^2 + k(1)^2;
f=a*k(1) + k(1);
p(3)=(-1/f)*( p(2)*( k(2) + k(1) +a*k(2) ) +2*log(a+1)/(a+1) +2*k(1)*k(2));
k(2)= (1/2)*(- p(2) - p(3) - a*p(3));

```

coeffd2.m

```
%Calculate the deg=2 coefficients of Prager’’s Example
function [k,p]=Coeffd2(a,k,p)

l=pragercost(a,k(1));
f=pragerf(a,k(1));
p(4)=(-1/f)*( 2*p(3)*(k(2)+ k(1) + a*k(2)) + p(2)*(k(3)+2*k(2)+a*k(3)) +2/(a+1)^2 -2*log(a+1)/(a+1)^2 + 2*k(2)^2 +2*k(1)*k(3));
k(3)=(-1/2)*(p(3) + p(4) +a*p(4));

```

coeff3.m

```
%Calculate the deg=3 coefficients of Prager’’s Example
function [k,p]=Coeffd3(a,k,p)

l=pragercost(a,k(1));
f=pragerf(a,k(1));
p(5)=(-1/f)*( 3*p(4)*(k(2) + k(1) + a*k(2)) + 3*p(3)*(k(3) + 2*k(2) + a*k(3))
+ p(2)*(k(4)+3*k(3)+a*k(4)) -6/(a+1)^3 + 4*log(a+1)/(a+1)^3 + 6*k(2)*k(3) +2*k(1)*k(4));
k(4)=(-1/2)*(3*p(4) + p(5) + a*p(5));

```