4 The Pontryagin Maximum Principle

This important chapter moves us beyond the linear dynamics assumed in Chapters 2 and 3, to consider much wider classes of optimal control problems, to introduce the fundamental Pontryagin Maximum Principle, and to illustrate its uses in a variety of examples.

4.1 Calculus of Variations, Hamiltonian Dynamics

We begin in this section with a quick introduction to some variational methods. These ideas will later serve as motivation for the Pontryagin Maximum Principle.

Assume we are given a smooth function $L:{\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathbb{R},L=L(x,v);L$ is called the Lagrangian. Let $T>0,x^0,x^1\in {\mathbb{R}}^n$ be given.

Basic Problem of the Calculus of Variations

Find a curve ${\mathit{x}}^{\ast }(\cdot ):[0,T]\to {\mathbb{R}}^n$ that minimizes the functional

$$I[\mathit{x}(\cdot )]:={\int }_0^{T}L(\mathit{x}(t),\dot{\mathit{x}}(t))\mathrm{d}t$$

among all fintions $\mathit{x}(\cdot )$ satisfying $\mathit{x}(0)=x^0$ and $\mathit{x}(T)=x^1$.

Now assume ${\mathit{x}}^{\ast }(\cdot )$ solves our variational problem. The fundamental question is this: how can we characterize ${\mathit{x}}^{\ast }(\cdot )$?

4.1.1 Derivation of Euler-Lagrange Equations

Notation

We write $L=L(x,v)$, and regard the variable $x$ as denoting position, the variable $v$ as denoting velocity. The partial derivatives of $L$ are

$$\frac{\partial L}{\partial x_i}=L_{x_i},\frac{\partial L}{\partial v_i}=L_{v_i}(1\le i\le n),$$

and we write

$${\mathrm{\nabla }}_{x}L:=(L_{x_1},\cdots ,L_{x_n}),{\mathrm{\nabla }}_{v}L:=(L_{v_1},\cdots ,L_{v_n}).$$

Theorem 4.1: Euler-Lagrange Equations

Let ${\mathit{x}}^{\ast }(\cdot )$ solve the calculus of variations problem. Then ${\mathit{x}}^{\ast }(\cdot )$ solves the Euler–Lagrange differential equations:

$$\begin{array}{}\text{(E-L)}& \frac{\mathrm{d}}{\mathrm{d}t}[{\mathrm{\nabla }}_{v}L({\mathit{x}}^{\ast }(t),{\dot{\mathit{x}}}^{\ast }(t))]={\mathrm{\nabla }}_{x}L({\mathit{x}}^{\ast }(t),{\dot{\mathit{x}}}^{\ast }(t)).\end{array}$$

The significance of preceding theorem is that if we can solve the Euler–Lagrange equations (E-L), then the solution of our original calculus of variations problem (assuming it exists) will be among the solutions.

Note that (E-L) is a quasilinear system of $n$ second–order ODE. The $i^{\mathrm{th}}$ component of the system reads

$$\frac{\mathrm{d}}{\mathrm{d}t}[L_{v_i}({\mathit{x}}^{\ast }(t),{\dot{\mathit{x}}}^{\ast }(t))]=L_{x_i}({\mathit{x}}^{\ast }(t),{\dot{\mathit{x}}}^{\ast }(t)).$$

Proof:

1. Select any smooth curve $\mathit{y}[0,T]\to {\mathbb{R}}^n$, satisfying $\mathit{y}(0)=\mathit{y}(T)=0$(Why? You'll know it). Define

$$i(\tau ):=I({\mathit{x}}^{\ast }(\cdot )+\tau {\mathit{y}}^{\ast }(\cdot )).$$

for $\tau \in \mathbb{R}$ and $\mathit{x}(\cdot )={\mathit{x}}^{\ast }(\cdot )$. (To simplify we omit the superscript $\ast$.) Notice that $\mathit{x}(\cdot )+\tau \mathit{y}(\cdot )$ takes on the proper values at the endpoints ($0,T$). Hence, since $\mathit{x}(\cdot )$ is minimizer (for the variational problem), we have

$$i(\tau )\ge I[\mathit{x}(\cdot )]=i(0).$$

Consequently $i(\cdot )$ has a minimum at $\tau =0$, and so

$$i^{\prime }(0)=0.$$

2. We must compute $i^{\prime }(\tau )$. Note first that

$$i(\tau )={\int }_0^{T}L(\mathit{x}(t)+\tau \mathit{y}(t),\dot{\mathit{x}}(t)+\tau \dot{\mathit{y}}(t))\mathrm{d}t.$$

and hence

$$i^{\prime }(\tau )={\int }_0^T(\sum _{i=1}^n{L}_{x_i}(\mathit{x}(t)+\tau \mathit{y}(t),\dot{\mathit{x}}(t)+\tau \dot{\mathit{y}}(t)){\mathit{y}}_i(t)+\sum _{i=1}^n{L}_{v_i}(\cdots ){\dot{\mathit{y}}}_i(t))\mathrm{d}t.$$

(Why sum: Chain rule) Let $\tau =0$. Then

$$0=i^{\prime }(0)=\sum _{i=1}^n{\int }_0^T{L}_{x_i}(\mathit{x}(t),\dot{\mathit{x}}(t)){\mathit{y}}_i(t)+L_{v_i}(\mathit{x}(t),\dot{\mathit{x}}(t)){\dot{\mathit{y}}}_i(t)\mathrm{d}t.$$

This equality holds for all choices of $\mathit{y}:[0,T]\to {\mathbb{R}}^n$, with $\mathit{y}(0)=\mathit{y}(T)=0$.

3. Fix any $1\le j\le n$. Choose $\mathit{y}(\cdot )$ s.t.

$$y_j(t)=\psi (t),y_i(t)\equiv 0,\mathrm{\forall }i\ne j,$$

where $\psi$ is an arbitary function. Use this choice of $\mathit{y}(\cdot )$ above:

$$0={\int }_0^T{L}_{x_j}(\mathit{x}(t),\dot{\mathit{x}}(t))\psi (t)+L_{v_j}(\mathit{x}(t),\dot{\mathit{x}}(t))\dot{\psi }(t)\mathrm{d}t.$$

Integrate by parts, recalling that $\psi (0)=\psi (T)=0\to [L_{v_j}\psi (t)]{|}_0^T=0$:

$$0={\int }_0^T[L_{x_j}(\mathit{x}(t),\dot{\mathit{x}}(t))-\frac{\mathrm{d}}{\mathrm{d}t}(L_{v_j}(\mathit{x}(t),\dot{\mathit{x}}(t)))]\psi (t)\mathrm{d}t.$$

This holds for alll $\psi :[0,T]\to \mathbb{R},\psi (0)=\psi (T)=0$ and therefore

$$L_{x_j}(\mathit{x}(t),\dot{\mathit{x}}(t))-\frac{\mathrm{d}}{\mathrm{d}t}(L_{v_j}(\mathit{x}(t),\dot{\mathit{x}}(t)))=0$$

for all times $0\le t\le T$. To see this, observe that otherwise $L_{x_j}-\frac{\mathrm{d}}{\mathrm{d}t}(L_{v_j})$ would be, say, positive on some subinterval on $I\subseteq [0,T]$. Choose $\psi \equiv 0$ off $I,\psi >0$ on $I$. Then

$${\int }_0^T(L_{x_j}-\frac{\mathrm{d}}{\mathrm{d}t}(L_{v_j}))\psi \mathrm{d}t>0,$$

a contradiction.

$◻$

4.1.2 Conversion to Hamilton's Equations

Definition: generalized momentum

For the given curve $\mathit{x}(\cdot )$, define

$$\mathit{p}(t):={\mathrm{\nabla }}_{v}L(\mathit{x}(t)\dot{\mathit{x}}(t))(0\le t\le T).$$

We call $\mathit{p}(\cdot )$ the generalized momentum

Our intention now is to rewrite the Euler–Lagrange equations as a system of first–order ODE for $\mathit{x}(\cdot ),\mathit{p}(\cdot )$.

Important Hypothesis: Assume that for all $x,p\in \mathbb{R}$, we can solve the equation

$$p={\mathrm{\nabla }}_{v}L(x,v)$$

for $v$ in terms of $x$ and $p$. That is, we suppose we can solve the identity for $v=\mathit{v}(x,p)$.

Definition: dynamical systems Hamiltonian $H$

Define the dynamical systems Hamiltonian $H:{\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathbb{R}$ by the formula

$$H(x,p)=p\cdot \mathit{v}(x,p)-L(x,\mathit{v}(x,p)),$$

where $\mathit{v}$ id defined above.

Notation

The partial derivatives of $H$ are

$$\frac{\partial H}{\partial x_i}=H_{x_i},\frac{\partial H}{\partial p_i}=H_{p_i}(1\le i\le n),$$

and we write

$${\mathrm{\nabla }}_{x}H:=(H_{x_1},\cdots ,H_{x_n}),{\mathrm{\nabla }}_{p}H:=(H_{p_1},\cdots ,H_{p_n}).$$

Theorem 4.2: Hamiltonian Dynamics

Let $\mathit{x}(\cdot )$ solve the Euler-Lagrange equations (E-L) and define $\mathit{p}(\cdot )$ as above. Then the pair $(\mathit{x}(\cdot ),\mathit{p}(\cdot ))$ solves Hamilton's equations:

$$\{\begin{array}{ll}\dot{\mathit{x}}(t)={\mathrm{\nabla }}_{p}H(\mathit{x}(t),\mathit{p}(t))& (ODE)\\ \dot{\mathit{p}}(t)=-{\mathrm{\nabla }}_{x}H(\mathit{x}(t),\mathit{p}(t))& (ADJ)\end{array}$$

Furthermore, the mapping $t\to H(\mathit{x}(t),\mathit{p}(t))$ is constant.

Proof: Recall that $H(x,p)=p\cdot \mathit{v}(x,p)-L(x,\mathit{v}(x,p))$, where $v=\mathit{v}(x,p)$ or, equivalently, $p={\mathrm{\nabla }}_{v}L(x,v)$. then

$$\begin{array}{rl}{\mathrm{\nabla }}_{x}H(x,p)& =p\cdot {\mathrm{\nabla }}_x\mathit{v}-[{\mathrm{\nabla }}_{x}L(x,\mathit{v}(x,p))+{\mathrm{\nabla }}_{v}L(x,\mathit{v}(x,p))\cdot {\mathrm{\nabla }}_x\mathit{v}]\\ & =[p-{\mathrm{\nabla }}_{v}L(x,v)]\cdot {\mathrm{\nabla }}_x\mathit{v}-{\mathrm{\nabla }}_{x}L(x,\mathit{v}(x,p))\\ & =-{\mathrm{\nabla }}_{x}L(x,\mathit{v}(x,p))\end{array}$$

(Why not consider $p$ in interms of $x$? See $H(x,p)$, they are currently divided)

because $p={\mathrm{\nabla }}_{v}L$. Now $\mathit{p}(t)={\mathrm{\nabla }}_{v}L(\mathit{x}(t),\dot{\mathit{x}}(t))$ iff $\dot{\mathit{x}}(t)=\mathit{v}(\mathit{x}(t),\mathit{p}(t))$. Therefore (E-L) implies

$$\begin{array}{rl}\dot{\mathit{p}}(t)& ={\mathrm{\nabla }}_{x}L(\mathit{x}(t),\dot{\mathit{x}}(t))\\ & ={\mathrm{\nabla }}_{x}L(\mathit{x}(t),\mathbf{(}\mathit{x}(t),\mathit{p}(t)))=-{\mathrm{\nabla }}_{x}H(\mathit{x}(t),\mathit{p}(t)).\end{array}$$

Also

$${\mathrm{\nabla }}_{p}H(x,p)=\mathit{v}(x,p)+p\cdot {\mathrm{\nabla }}_p\mathit{v}-{\mathrm{\nabla }}_{v}L\cdot {\mathrm{\nabla }}_p\mathit{v}=\mathit{v}(x,p)$$

since $p={\mathrm{\nabla }}_{v}L(x,\mathit{v}(x,p))$. This implies

$${\mathrm{\nabla }}_{p}H(\mathit{x}(t),\mathit{p}(t))=\mathit{v}(\mathit{x}(t),\mathit{p}(t)).$$

But

$$\mathit{p}(t)={\mathrm{\nabla }}_{v}L(\mathit{x}(t),\dot{\mathit{x}}(t))$$

and so $\dot{\mathit{x}}(t)=\mathit{v}(\mathit{x}(t),\mathit{p}(t))$. Therefore

$$\dot{\mathit{x}}(t)={\mathrm{\nabla }}_{p}H(\mathit{x}(t),\mathit{p}(t)).$$

Finaly note that

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}H(\mathit{x}(t),\mathit{p}(t))& ={\mathrm{\nabla }}_{x}H\cdot \dot{\mathit{x}}(t)+{\mathrm{\nabla }}_{p}H\cdot \dot{\mathit{p}}(t)\\ & ={\mathrm{\nabla }}_{x}H\cdot {\mathrm{\nabla }}_{p}H+{\mathrm{\nabla }}_{p}H\cdot (-{\mathrm{\nabla }}_{x}H)=0.\end{array}$$

$\square$

A Physical Example

We define the Lagrangian

$$L(x,v)=\frac{m|v{|}^2}{2}-V(x),$$

which we interpret as the kinetic energy minus the potential energy $V$. then

$$\{\begin{array}{l}{\mathrm{\nabla }}_{x}L=-\mathrm{\nabla }V(x)\\ {\mathrm{\nabla }}_{v}L=mv.\end{array}$$

Therefore the Euler-Lagrange equation is

$$m\ddot{x}(t)=-\mathrm{\nabla }V(\mathit{x}(t)),$$

which is Newton’s law. Furthermore

$$p={\mathrm{\nabla }}_{v}L(x,v)=mv$$

is the momentum, and the Hamiltonian is

$$\begin{array}{rl}H(x,p)& =p\cdot \frac{p}{m}-L(x,\frac{p}{m})\\ & =\frac{|p{|}^2}{m}-\frac{m}{2}{\left|\frac{p}{m}\right|}^2+V(x)=\frac{|p{|}^2}{2m}+V(x).\end{array}$$

the sum of the kinetic and potential energies. For this example, Hamilton’s equations read

$$\{\begin{array}{l}\dot{\mathit{x}}(t)=\frac{\mathit{p}(t)}{m}\\ \dot{\mathit{p}}(t)=-\mathrm{\nabla }V(\mathit{x}(t)).\end{array}$$

4.2 Review of Lagrange Multipliers

Constraint and Lagrange Multipliers

What first strikes us about general optimal control problems is the occurence of many constraints, most notably that the dynamics be governed by the differential equation

$$\{\begin{array}{l}\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{\alpha }(t))(t>0)\\ \mathit{x}(0)=x^0.\end{array}$$

This is in contrast to standard calculus of variations problems, as discussed in §4.1, where we could take any curve $\mathit{x}(\cdot )$ as a candidate for a minimizer.

Now it is a general principle of variational and optimization theory that “constraints create Lagrange multipliers” and furthermore that these Lagrange multipliers often “contain valuable information”. This section provides a quick review of the standard method of Lagrange multipliers in solving multivariable constrained optimization problems.

Unconstrainted Optimization

Suppose first that we wish to find a maximum point for a given smooth function $f:{\mathbb{R}}^n\to \mathbb{R}$. In this case there is no constraint, and therefore if $f(x\ast )=\underset{x\in {\mathbb{R}}^n}{max}f(x)$, then $x\ast$ is a critical point of $f$:

$$\mathrm{\nabla }f(x\ast )=0.$$

Constrainted Optimization

We modify the problem above by introducing the region

$$R:=\{x\in {\mathbb{R}}^n|g(x)\le 0\},$$

determined by some given function $g:{\mathbb{R}}^n\to \mathbb{R}$. Suppose $x\ast \in R$ and $f(x\ast )=\underset{x\in R}{max}f(x)$. We would like a characterization of $x\ast$ in terms of the gradients of $f$ and $g$.

• Case 1: $x\ast$ lies in the interior of $R$ Then the constraint is inactive, and so

  $$\mathrm{\nabla }f(x\ast )=0.$$

• Case 2: $x\ast$ lies on $\partial R$ We look at the direction of the vector $\mathrm{\nabla }f(x\ast )$. A geometric picture like Figure above is impossible; for if it were so, then $f(y\ast )$ would be greater that $f(x\ast )$ for some other point $y\ast \in \partial R$. So it must be $\mathrm{\nabla }f(x\ast )$ is perpendicular to $\partial R$ at $x\ast$ (Suppose not perpendicular, then there exists a unit vector $\tau$ in the tangential direction s.t. $\tau \cdot \mathrm{\nabla }f(x\ast )>0$, and this is the direction derivative) Since $\mathrm{\nabla }g$ is perpendicular to $\partial R=\{g=0\}$ (gradient is the fastest change direction, so normal to perpendicular), it follows that $\mathrm{\nabla }f(x\ast )$ is parallel to $\mathrm{\nabla }g(x\ast )$. Therefore

  $$\begin{array}{}\text{(4.4)}& \mathrm{\nabla }f(x\ast )=\lambda \mathrm{\nabla }g(x\ast )\end{array}$$

  for some real number $\lambda$, called a Lagrange multiplier.

Critique

The foregoing argument is in fact incomplete, since we implicitly assumed that $\mathrm{\nabla }g(x^{\ast })\ne 0$, in which case the Implicit Function Theorem implies that the set $\{g=0\}$ is an $(n-1)$ -dimension surface near $x^{\ast }$ (as illustrated).

If instead $\mathrm{\nabla }g(x^{\ast })=0$, the set $\{g=0\}$ need not have this simple form near $x^{\ast }$; and the reasoning discussed as Case 2 above is not complete.

The correct statement is this:

There exist real numbers $\lambda ,\mu$ not both equal to $0$, s.t.

$$\begin{array}{}\text{(4.5)}& \mu \mathrm{\nabla }f(x^{\ast })=\lambda \mathrm{\nabla }g(x^{\ast }).\end{array}$$

If $\mu \ne 0$, we can divide by $\mu$ and convert to the formulation (4.4). And if $\mathrm{\nabla }g(x^{\ast })=0$, we can take $\lambda =1,\mu =0$, making assertion (4.5) correct (if not particularly useful).

4.3 Statemet of Pontryagin Maximum Principle

We come now to the key assertion of this chapter, the theoretically interesting and practically useful theorem that if ${\mathit{\alpha }}^{\ast }(\cdot )$ is an optimal control, then there exists a function ${\mathit{p}}^{\ast }(\cdot )$, called the costate, that satisfies a certain maximization principle. We should think of the function ${\mathit{p}}^{\ast }(\cdot )$ as a sort of Lagrange multiplier, which appears owing to the constraint that the optimal curve ${\mathit{x}}^{\ast }(\cdot )$ must satisfy (ODE). And just as conventional Lagrange multipliers are useful for actual calculations, so also will be the costate.

We quote Francis Clarke [C2]: “The maximum principle was, in fact, the culmination of a long search in the calculus of variations for a comprehensive multiplier rule, which is the correct way to view it: $p(t)$ is a “Lagrange multiplier” ... It makes optimal control a design tool, whereas the calculus of variations was a way to study nature.”

4.3.1 Fixed Time, Free Endpoint Problem

Let us review the basic set-up for our control problem.

We are given $A\subseteq {\mathbb{R}}^m$ and also $f:{\mathbb{R}}^n\times A\to {\mathbb{R}}^n,x^0\in {\mathbb{R}}^n$. We as before denote the set of admissible controls by

$$\mathcal{A}=\{\mathit{\alpha }(\cdot ):[0,\mathrm{\infty })\to A\mid \mathit{\alpha }(\cdot )\text{ is measurable}\}.$$

Then given $\mathit{\alpha }(\cdot )\in \mathcal{A}$, we solve for the corresponding evolution of our system:

$$\{\begin{array}{l}\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{\alpha }(t))(t\ge 0)\\ \mathit{x}(0)=x^0\end{array}$$

We also introduce the payoff functional (To show the influence of free endpoint)

$$\begin{array}{}\text{(P)}& P[\mathit{\alpha }(\cdot )]={\int }_0^{T}r(\mathit{x}(t),\mathit{\alpha }(t))\mathrm{d}t+g(\mathit{x}(T))\end{array}$$

where the terminal time $T>0$, running payoff $r:{\mathbb{R}}^n\times A\to \mathbb{R}$ and terminal payoff $g:{\mathbb{R}}^n\to \mathbb{R}$ are given.

Basic Problem

Find a control ${\mathit{\alpha }}^{\ast }(\cdot )$ such that

$$P[{\mathit{\alpha }}^{\ast }(\cdot )]=\underset{\mathit{\alpha }(\cdot )\in \mathcal{A}}{max}P[\mathit{\alpha }(\cdot )]$$

The Pontryagin Maximum Principle, stated below, asserts the existence of a function ${\mathit{p}}^{\ast }(\cdot )$, which together with the optimal trajectory ${\mathit{x}}^{\ast }(\cdot )$ satisfies an analog of Hamilton's ODE from §4.1.2. For this, we will need an appropriate Hamiltonian:

Definition

The control theory Hamiltonian is the function

$$H(x,p,a):=\mathit{f}(x,a)\cdot p+r(x,a)(x,p\in {\mathbb{R}}^n,a\in A)$$

Theorem 4.3 (Pontryagin Maximum Principle)

Assume ${\mathit{\alpha }}^{\ast }(\cdot )$ is optimal for ( ODE ), ( $P$ ), and ${\mathit{x}}^{\ast }(\cdot )$ is the corresponding trajectory.

Then there exists a function ${\mathit{p}}^{\ast }:[0,T]\to {\mathbb{R}}^n$ such that

$$\begin{array}{}\text{(ODE)}& {\dot{\mathit{x}}}^{\ast }(t)={\mathrm{\nabla }}_{p}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)),\text{(ADJ)}& {\dot{\mathit{p}}}^{\ast }(t)=-{\mathrm{\nabla }}_{x}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)),\end{array}$$

and

$$\begin{array}{}\text{(M)}& H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))=\underset{a\in A}{max}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),a)(0\le t\le T)\end{array}$$

In addition,

$$\text{ the mapping }t\mapsto H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))\text{ is constant. }$$

Finally, we have the terminal condition

$${\mathit{p}}^{\ast }(T)=\mathrm{\nabla }g({\mathit{x}}^{\ast }(T))$$

Remarks and Intepretations

1. The identities (ADJ) are the adjoint equations and (M) the maximization principle. Notice that (ODE) and (ADJ) resemble the structure of Hamilton's equations, discussed in §4.1. We also call ( $T$ ) the transversality condition and will discuss its significance later.

2. More precisely, formula (ODE) says that for $1\le i\le n$, we have

$${\dot{x}}^{i\ast }(t)=H_{p_i}({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))=f^i({\mathit{x}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)),$$

which is just the original equation of motion. Likewise, (ADJ) says

$$\begin{array}{rl}{\dot{p}}^{i\ast }(t)& =-H_{x_i}({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))\\ & =-\sum _{j=1}^n{p}^{j\ast }(t)f_{x_i}^j({\mathit{x}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))-r_{x_i}({\mathit{x}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)).\end{array}$$

4.3.2 Free Time, Fixed Endpoint Problem

Let us next record the appropriate form of the Maximum Principle for a fixed endpoint problem.

As before, given a control $\mathit{\alpha }(\cdot )\in \mathcal{A}$, we solve for the corresponding evolution of our system:

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{\alpha }(t))(t\ge 0)\\ \mathit{x}(0)=x^0\end{array}\end{array}$$

Assume now that a target point $x^1\in {\mathbb{R}}^n$ is given. We introduce then the payoff functional

$$P[\mathit{\alpha }(\cdot )]={\int }_0^{\tau }r(\mathit{x}(t),\mathit{\alpha }(t))\mathrm{d}t.$$

Here $r:{\mathbb{R}}^n\times A\to \mathbb{R}$ is the given running payoff, and $\tau =\tau [\mathit{\alpha }(\cdot )]\le \mathrm{\infty }$ denotes the first time the solution of (ODE) hits the target point $x^1$.

As before, the basic problem is to find an optimal control ${\mathit{\alpha }}^{\ast }(\cdot )$ such that

$$\begin{array}{}\text{(P)}& P[{\mathit{\alpha }}^{\ast }(\cdot )]=\underset{\mathit{\alpha }(\cdot )\in \mathcal{A}}{max}P[\mathit{\alpha }(\cdot )].\end{array}$$

Define the Hamiltonian $H$ as in §4.3.1.

Theorem 4.4 (Pontryagin Maximum Principle)

Assume ${\mathit{\alpha }}^{\ast }(\cdot )$ is optimal for (ODE), ($P$) and ${\mathit{x}}^{\ast }(\cdot )$ is the corresponding trajectory.

Then there exists a function ${\mathit{p}}^{\ast }:[0,{\tau }^{\ast }]\to {\mathbb{R}}^n$ such that

$$\begin{array}{}\text{(ODE)}& {\dot{\mathit{x}}}^{\ast }(t)={\mathrm{\nabla }}_{p}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)),\text{(ADJ)}& {\dot{\mathit{p}}}^{\ast }(t)=-{\mathrm{\nabla }}_{x}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t)),\end{array}$$

and

$$\begin{array}{}\text{(M)}& H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))=\underset{a\in A}{max}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),a)(0\le t\le {\tau }^{\ast })\end{array}$$

Also,

$$H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))\equiv 0(0\le t\le {\tau }^{\ast })$$

Here ${\tau }^{\ast }$ denotes the first time the trajectory $x^{\ast }(\cdot )$ hits the target point $x^1$. We call ${\mathit{x}}^{\ast }(\cdot )$ the state of the optimally controlled system and ${\mathit{p}}^{\ast }(\cdot )$ the costate.

Remark and Warning

More precisely, we should define

$$H(x,p,q,a)=\mathit{f}(x,a)\cdot p+r(x,a)q(q\in \mathbb{R})$$

A more careful statement of the Maximum Principle says "there exists a constant $q\ge 0$ and a function ${\mathit{p}}^{\ast }:[0,t^{\ast }]\to {\mathbb{R}}^n$ such that (ODE), (ADJ), and (M) hold".

• If $q>0$, we can renormalize to get $q=1$, as we have done above.
• If $q=0$, then $H$ does not depend on running payoff $r$ and in this case the Pontryagin Maximum Principle is not useful. This is a so-called "abnormal problem".

Compare these comments with the critique of the usual Lagrange multiplier method at the end of §4.2, and see also the proof in §A. 5 of the Appendix.

4.4 Application and Examples

How to Use the Maximum Principle

We mentioned earlier that the costate ${\mathit{p}}^{\ast }(\cdot )$ can be interpreted as a sort of Lagrange multiplier.

Calculations with Lagrange multipliers

Recall our discussion in §4.2 about finding a point $x^{\ast }$ that maximizes a function $f$, subject to the requirement that $g\le 0$. Now $x^{\ast }={(x_1^{\ast },\dots ,x_n^{\ast })}^T$ has $n$ unknown components we must find. Somewhat unexpectedly, it turns out in practice to be easier to solve (4.4) for the $n+1$ unknowns $x_1^{\ast },\dots ,x_n^{\ast }$ and $\lambda$. We repeat this key insight: it is actually easier to solve the problem if we add a new unknown, namely the Lagrange multiplier. Worked examples abound in multivariable calculus books.

Calculations with the costate

This same principle is valid for our much more complicated control theory problems: it is usually best not just to look for an optimal control ${\mathit{\alpha }}^{\ast }(\cdot )$ and an optimal trajectory ${\mathit{x}}^{\ast }(\cdot )$ alone, but also to look as well for the costate ${\mathit{p}}^{\ast }(\cdot )$. In practice, we add the equations (ADJ) and (M) to (ODE) and then try to solve for ${\mathit{\alpha }}^{\ast }(\cdot ),{\mathit{x}}^{\ast }(\cdot )$ and for ${\mathit{p}}^{\ast }(\cdot )$.

The following examples show how this works in practice, in certain cases for which we can actually solve everything explicitly or, failing that, at least deduce some useful information.

4.4.1 Example 1: Linear Time-Optimal Control.

For this example, let $A$ denote the cube $[-1,1{]}^n$ in ${\mathbb{R}}^n$. We consider again the linear dynamics:

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{\mathit{x}}(t)=M\mathit{x}(t)+N\mathit{\alpha }(t)\\ \mathit{x}(0)=x^0\end{array}\end{array}$$

for the payoff functional

$$P[\mathit{\alpha }(\cdot )]=-{\int }_0^{\tau }1dt=-\tau$$

where $\tau$ denotes the first time the trajectory hits the target point $x^1=0$. We have $r\equiv -1$, and so

$$H(x,p,a)=\mathit{f}\cdot p+r=(Mx+Na)\cdot p-1$$

In §3.2 we introduced the Hamiltonian $H=(Mx+Na)\cdot p$, which differs by a constant from the present $H$. We can redefine $H$ in §3.2 to match the present theory: compare then Theorems 3.4 and 4.4.

4.4.2 Example 2: Control of Production and Consumption

We return to Example 1 in Chapter 1, a model for optimal consumption in a simple economy. Recall that

$$\begin{array}{rl}& x(t)=\text{ output of economy at time }t\\ & \alpha (t)=\text{ fraction of output reinvested at time }t\end{array}$$

We have the constraint $0\le \alpha (t)\le 1$; that is, $A=[0,1]\subset \mathbb{R}$. The economy evolves according to the dynamics

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{x}(t)=\alpha (t)x(t)(0\le t\le T)\\ x(0)=x^0\end{array}\end{array}$$

where $x^0>0$ and we have set the growth factor $k=1$. We want to maximize the total consumption

$$P[\alpha (\cdot )]:={\int }_0^T(1-\alpha (t))x(t)\mathrm{d}t$$

How can we characterize an optimal control ${\alpha }^{\ast }(\cdot )$ ?

Introducing the maximum principle

We apply Pontryagin Maximum Principle, and to simplify notation we will not write the superscripts $\ast$ for the optimal control, trajectory, etc. We have $n=m=1$,

$$f(x,a)=xa,g\equiv 0,r(x,a)=(1-a)x$$

and therefore

$$\begin{array}{rl}H(x,p,a)& =f(x,a)p+r(x,a)\\ & =pxa+(1-a)x=x+ax(p-1)\end{array}$$

The dynamical equation is

$$\begin{array}{}\text{(ODE)}& \dot{x}(t)=H_p=\alpha (t)x(t)\end{array}$$

and the adjoint equation is

$$\begin{array}{}\text{(ADJ)}& \dot{p}(t)=-H_x=-1-\alpha (t)(p(t)-1)\end{array}$$

The terminal condition reads

$$\begin{array}{}\text{(T)}& p(T)=g_x(x(T))=0\end{array}$$

Lastly, the maximality principle asserts

$$\begin{array}{}\text{(M)}& H(x(t),p(t),\alpha (t))=\underset{0\le a\le 1}{max}\{x(t)+ax(t)(p(t)-1)\}\end{array}$$

Using the maximum principle

We now deduce useful information from (ODE), (ADJ), (M) and (T).

According to (M), at each time $t$ the control value $\alpha (t)$ must be selected to maximize $a(p(t)-1)$ for $0\le a\le 1$. This is so, since $x(t)>0$ (It's fixed for the control is decided by the current state $x(t)$). Thus ($a$ is a placeholder variable used to denote $\alpha (t)$)

$$\alpha (t)=\{\begin{array}{ll}1& \text{ if }p(t)>1\\ 0& \text{ if }p(t)\le 1\end{array}$$

Hence if we know $p(\cdot )$, we can design the optimal control $\alpha (\cdot )$. So next we must solve for the costate $p(\cdot )$. We know from $(\mathrm{ADJ})$ and $(\mathrm{T})$ that

$$\{\begin{array}{ll}\dot{p}(t)=-1-\alpha (t)[p(t)-1]& (0\le t\le T)\\ p(T)=0& \end{array}$$

Since $p(T)=0$, we deduce by continuity that $p(t)\le 1$ for $t$ close to $T,t<T$. Thus $\alpha (t)=0$ for such values of $t$. Therefore $\dot{p}(t)=-1$, and consequently $p(t)=T-t$ for times $t$ in this interval. So we have that $p(t)=T-t$ so long as $p(t)\le 1$. And this holds for $T-1\le t\le T$

But for times $t\le T-1$, with $t$ near $T-1$, we have $\alpha (t)=1$; and so ( $ADJ$ ) becomes

$$\dot{p}(t)=-1-(p(t)-1)=-p(t)$$

Since $p(T-1)=1$, we see that $p(t)=e^{T-1-t}>1$ for all times $0\le t\le T-1$. In particular there are no switches in the control over this time interval.

Restoring the superscript * to our notation, we consequently deduce that an optimal control is

$${\alpha }^{\ast }(t)=\{\begin{array}{ll}1& \text{ if }0\le t\le t^{\ast }\\ 0& \text{ if }{t}^{\ast }\le t\le T\end{array}$$

for the optimal switching time $t^{\ast }=T-1$.

We leave it as an exercise to compute the switching time if the growth constant $k\ne 1$. $◻$

4.4.3 Example 3: A Simple Linear-Quadratic Regulator

We take $n=m=1$ for this example, and consider the simple linear dynamics

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{x}(t)=x(t)+\alpha (t)\\ x(0)=x^0\end{array}\end{array}$$

with the quadratic cost functional

$${\int }_0^{T}x(t{)}^2+\alpha (t{)}^2\mathrm{d}t,$$

which we want to minimize. So we want to maximize the payoff functional

$$\begin{array}{}\text{(P)}& P[\alpha (\cdot )]=-{\int }_0^{T}x(t{)}^2+\alpha (t{)}^2\mathrm{d}t.\end{array}$$

For this problem, the values of the controls are not constrained; that is, $A=\mathbb{R}$.

Introducing the maximum principle

To simplify notation further we again drop the superscripts *. We have $n=m=1$,

$$f(x,a)=x+a,g\equiv 0,r(x,a)=-x^2-a^2;$$

and hence

$$H(x,p,a)=fp+r=(x+a)p-(x^2+a^2)$$

The maximality condition becomes

$$\begin{array}{}\text{(M)}& H(x(t),p(t),\alpha (t))=\underset{a\in \mathbb{R}}{max}\{-(x(t{)}^2+a^2)+p(t)(x(t)+a)\}\end{array}$$

We calculate the maximum on the right hand side by setting $H_a=-2a+p=0$. Thus $a=\frac{p}{2}$, and so

$$\alpha (t)=\frac{p(t)}{2}.$$

The dynamical equations are therefore

$$\begin{array}{}\text{(ODE)}& \dot{x}(t)=x(t)+\frac{p(t)}{2}\end{array}$$

and

$$\begin{array}{}\text{(ADJ)}& \dot{p}(t)=-H_x=2x(t)-p(t).\end{array}$$

Moreover $x(0)=x^0$, and the terminal condition is

$$\begin{array}{}\text{(T)}& p(T)=0.\end{array}$$

Using the Maximum Principle

So we must look at the system of equations

$$(\genfrac{}{}{0ex}{}{\dot{x}}{\dot{p}})=\underset{M}{\underbrace{\left(\begin{array}{cc}1& 1/2\\ 2& -1\end{array}\right)}}(\genfrac{}{}{0ex}{}{x}{p}),$$

the general solution of which is

$$(\genfrac{}{}{0ex}{}{x(t)}{p(t)})=e^{tM}(\genfrac{}{}{0ex}{}{x^0}{p^0})$$

Since we know $x^0$, the task is to choose $p^0$ so that $p(T)=0$.

Feedback controls

An elegant way to do so is to try to find optimal control in linear feedback form; that is, to look for a function $c(\cdot ):[0,T]\to \mathbb{R}$ for which

$$\alpha (t)=c(t)x(t)$$

We henceforth suppose that an optimal feedback control of this form exists, and attempt to calculate $c(\cdot )$. Now

$$\frac{p(t)}{2}=\alpha (t)=c(t)x(t)$$

whence $c(t)=\frac{p(t)}{2x(t)}$. Define now

$$d(t):=\frac{p(t)}{x(t)}$$

so that $c(t)=\frac{d(t)}{2}$.

We will next discover a differential equation that $d(\cdot )$ satisfies. Compute

$$\dot{d}=\frac{\dot{p}}{x}-\frac{p\dot{x}}{x^2},$$

and recall that

$$\{\begin{array}{l}\dot{x}=x+\frac{p}{2}\\ \dot{p}=2x-p.\end{array}$$

Therefore

$$\dot{d}=\frac{2x-p}{x}-\frac{p}{x^2}(x+\frac{p}{2})=2-d-d(1+\frac{d}{2})=2-2d-\frac{d^2}{2}.$$

Since $p(T)=0$, the terminal condition is $d(T)=0$. So we have obtained a nonlinear first-order ODE for $d(\cdot )$ with a terminal boundary condition:

$$\begin{array}{}\text{(R)}& \{\begin{array}{l}\dot{d}=2-2d-\frac{1}{2}{d}^2(0\le t<T)\\ d(T)=0.\end{array}\end{array}$$

This is called the Riccati equation.

In summary so far, to solve our linear-quadratic regulator problem, we need to first solve the Riccati equation $(\mathrm{R})$ and then set

$$\alpha (t)=\frac{1}{2}d(t)x(t)$$

How to solve the Riccati equation. It turns out that we can convert ( $R$ ) it into a second-order, linear ODE. To accomplish this, write

$$d(t)=\frac{2\dot{b}(t)}{b(t)}$$

for a function $b(\cdot )$ to be found. What equation does $b(\cdot )$ solve? We compute

$$\dot{d}=\frac{2\ddot{b}}{b}-\frac{2(\dot{b}{)}^2}{b^2}=\frac{2\ddot{b}}{b}-\frac{d^2}{2}.$$

Hence $(\mathrm{R})$ gives

$$\frac{2\ddot{b}}{b}=\dot{d}+\frac{d^2}{2}=2-2d=2-2\frac{2\dot{b}}{b}$$

and consequently

$$\{\begin{array}{l}\ddot{b}=b-2\dot{b}(0\le t<T)\\ \dot{b}(T)=0,b(T)=1\end{array}$$

This is a terminal-value problem for second-order linear ODE, which we can solve by standard techniques. We then set $d=\frac{2\dot{b}}{b}$, to derive the solution of the Riccati equation $(\mathrm{R})$.

We will generalize this example later to systems, in §5.2.

4.4.4 Example 4: Moon Lander

This is a much more elaborate and interesting example, already introduced in Chapter 1.

Introduce the notation

$$\begin{array}{rl}h(t)& =\text{ height at time }t\\ v(t)& =\text{ velocity }=\dot{h}(t)\\ m(t)& =\text{ mass of spacecraft (changing as fuel is used up) }\\ \alpha (t)& =\text{ thrust at time }t\end{array}$$

The thrust is constrained so that $0\le \alpha (t)\le 1$; that is, $A=[0,1]$. There are also the constraints that the height and mass be nonnegative: $h(t)\ge 0,m(t)\ge 0$.

The dynamics are

$$\{\begin{array}{rl}\dot{h}(t)& =v(t)\\ \dot{v}(t)& =-g+\frac{\alpha (t)}{m(t)}\\ \dot{m}(t)& =-k\alpha (t)\end{array}$$

with initial conditions

$$\{\begin{array}{l}h(0)=h_0>0\\ v(0)=v_0\\ m(0)=m_0>0\end{array}$$

The goal is to land on the moon safely, maximizing the remaining fuel $m(\tau )$, where $\tau =\tau [\alpha (\cdot )]$ is the first time $h(\tau )=v(\tau )=0$. Since $\alpha =-\frac{\dot{m}}{k}$, our intention is equivalently to minimize the total applied thrust before landing; so that

$$P[\alpha (\cdot )]=-{\int }_0^{\tau }\alpha (t)dt$$

This is so since

$${\int }_0^{\tau }\alpha (t)dt=\frac{m_0-m(\tau )}{k}$$

Introducing the maximum principle

In terms of the general notation, we have

$$\mathit{x}(t)=\left[\begin{array}{c}h(t)\\ v(t)\\ m(t)\end{array}\right],\mathit{f}=\left[\begin{array}{c}v\\ -g+a/m\\ -ka\end{array}\right].$$

Hence the Hamiltonian is

$$\begin{array}{rl}H(x,p,a)& =\mathit{f}\cdot \mathit{p}+r\\ & =(v,-g+a/m,-ka)\cdot (p_1,p_2,p_3)-a\\ & =-a+p_{1}v+p_2(-g+\frac{a}{m})+p_3(-ka)\end{array}$$

We next have to figure out the adjoint dynamics $(\mathrm{ADJ})$. For our particular Hamiltonian,

$$H_{x_1}=H_h=0,H_{x_2}=H_v=p_1,H_{x_3}=H_m=-\frac{p_{2}a}{m^2}$$

Therefore

$$\begin{array}{}\text{(ADJ)}& \{\begin{array}{l}{\dot{p}}^1(t)=0\\ {\dot{p}}^2(t)=-p^1(t)\\ {\dot{p}}^3(t)=\frac{p^2(t)\alpha (t)}{m(t{)}^2}\end{array}\end{array}$$

The maximization condition $(\mathrm{M})$ reads

$$\begin{array}{rl}\text{(M)}& & H(\mathit{x}(t),\mathit{p}(t),\alpha (t))=\underset{0\le a\le 1}{max}H(\mathit{x}(t),\mathit{p}(t),a)=& \underset{0\le a\le 1}{max}\{-a+p^1(t)v(t)+p^2(t)[-g+\frac{a}{m(t)}]+p^3(t)(-ka)\}\\ =& p^1(t)v(t)-p^2(t)g+\underset{0\le a\le 1}{max}\left\{a(-1+\frac{p^2(t)}{m(t)}-k{p}^3(t))\right\}\end{array}$$

Thus the optimal control law is given by the rule:

$$\alpha (t)=\{\begin{array}{ll}1& \text{ if }1-\frac{p^2(t)}{m(t)}+k{p}^3(t)<0\\ 0& \text{ if }1-\frac{p^2(t)}{m(t)}+k{p}^3(t)>0\end{array}$$

Using the maximum principle

Now we will attempt to figure out the form of the solution, and check it accords with the Maximum Principle.

Let us start by guessing that we first leave rocket engine of (i.e., set $\alpha \equiv 0$ ) and turn the engine on only at the end. Denote by $\tau$ the first time that $h(\tau )=v(\tau )=0$, meaning that we have landed. We guess that there exists a switching time $t^{\ast }<\tau$ when we turned engines on at full power (i.e., set $\alpha \equiv 1$ ). Consequently,

$$\alpha (t)=\{\begin{array}{ll}0& \text{ for }0\le t\le t^{\ast }\\ 1& \text{ for }{t}^{\ast }\le t\le \tau \end{array}$$

Therefore, for times $t^{\ast }\le t\le \tau$ our ODE becomes

$$\{\begin{array}{l}\dot{h}(t)=v(t)\\ \dot{v}(t)=-g+\frac{1}{m(t)}\\ \dot{m}(t)=-k\end{array}(t^{\ast }\le t\le \tau )$$

with $h(\tau )=0,v(\tau )=0,m\left(t^{\ast }\right)=m_0$. We solve these dynamics:

$$\{\begin{array}{l}m(t)=m_0+k(t^{\ast }-t)\\ v(t)=g(\tau -t)+\frac{1}{k}\log \left[\frac{m_0+k(t^{\ast }-\tau )}{m_0+k(t^{\ast }-t)}\right]\\ h(t)=\text{ complicated formula }\end{array}$$

Now put $t=t^{\ast }$ :

$$\{\begin{array}{l}m\left(t^{\ast }\right)=m_0\\ v\left(t^{\ast }\right)=g(\tau -t^{\ast })+\frac{1}{k}\log \left[\frac{m_0+k(t^{\ast }-\tau )}{m_0}\right]\\ h\left(t^{\ast }\right)=-\frac{g{(t^{\ast }-\tau )}^2}{2}-\frac{m_0}{k^2}\log \left[\frac{m_0+k(t^{\ast }-\tau )}{m_0}\right]+\frac{t^{\ast }-\tau }{k}\end{array}$$

Suppose the total amount of fuel to start with was $m_1$; so that $m_0-m_1$ is the weight of the empty spacecraft. When $\alpha \equiv 1$, the fuel is used up at rate $k$. Hence

$$k(\tau -t^{\ast })\le m_1$$

and so $0\le \tau -t^{\ast }\le \frac{m_1}{k}$. Before time $t^{\ast }$, we set $\alpha \equiv 0$. Then (ODE) reads

$$\{\begin{array}{l}\dot{h}=v\\ \dot{v}=-g\\ \dot{m}=0\end{array}$$

and thus

$$\{\begin{array}{l}m(t)\equiv m_0\\ v(t)=-gt+v_0\\ h(t)=-\frac{1}{2}g{t}^2+t{v}_0+h_0\end{array}$$

We combine the formulas for $v(t)$ and $h(t)$, to discover

$$h(t)=h_0-\frac{1}{2g}(v^2(t)-v_0^2)(0\le t\le t^{\ast })$$

We deduce that the freefall trajectory $(v(t),h(t))$ therefore lies on a parabola

$$h=h_0-\frac{1}{2g}(v^2-v_0^2).$$

If we then move along this parabola until we hit the soft-landing curve from the previous picture, we can then turn on the rocket engine and land safely.

In the second case illustrated below, we miss switching curve, and hence cannot land safely on the moon switching once.

To justify our guess about the structure of the optimal control, let us now find the costate $\mathit{p}(\cdot )$ so that $\alpha (\cdot )$ and $\mathit{x}(\cdot )$ described above satisfy $(\mathrm{ODE}),(\mathrm{ADJ}),(\mathrm{M})$. To do this, we will have to figure out appropriate initial conditions

$$p^1(0)={\lambda }_1,p^2(0)={\lambda }_2,p^3(0)={\lambda }_3$$

We solve $(\mathrm{ADJ})$ for $\alpha (\cdot )$ as above, and find

$$\{\begin{array}{ll}{p}^1(t)\equiv {\lambda }_1& (0\le t\le \tau )\\ p^2(t)={\lambda }_2-{\lambda }_{1}t& (0\le t\le \tau )\\ p^3(t)=\{\begin{array}{ll}{\lambda }_3& (0\le t\le t^{\ast })\\ {\lambda }_3+{\int }_{t^{\ast }}^t\frac{{\lambda }_2-{\lambda }_{1}s}{(m_0+k(t^{\ast }-s){)}^2}\mathrm{d}s& (t^{\ast }\le t\le \tau ).\end{array}\end{array}$$

Define

$$r(t):=1-\frac{p^2(t)}{m(t)}+p^3(t)k$$

then

$$\dot{r}=-\frac{{\dot{p}}^2}{m}+\frac{p^2\dot{m}}{m^2}+{\dot{p}}^{3}k=\frac{{\lambda }_1}{m}+\frac{p^2}{m^2}(-k\alpha )+\left(\frac{p^2\alpha }{m^2}\right)k=\frac{{\lambda }_1}{m(t)}.$$

Choose ${\lambda }_1<0$, so that $r$ is decreasing. We calculate

$$r\left(t^{\ast }\right)=1-\frac{({\lambda }_2-{\lambda }_1{t}^{\ast })}{m_0}+{\lambda }_{3}k$$

and then adjust ${\lambda }_2,{\lambda }_3$ so that $r\left(t^{\ast }\right)=0$. Then $r$ is nonincreasing, $r\left(t^{\ast }\right)=0$, and consequently $r>0$ on $[0,t^{\ast }),r<0$ on $(t^{\ast },\tau ]$. But $(\mathrm{M})$ says

$$\alpha (t)=\{\begin{array}{ll}1& \text{ if }r(t)<0\\ 0& \text{ if }r(t)>0\end{array}$$

Thus an optimal control changes just once from 0 to 1 ; and so our earlier guess of $\alpha (\cdot )$ does indeed satisfy the Pontryagin Maximum Principle.

$◻$

4.5 Maximum Principle with Transversality Conditions

Consider again the dynamics

$$\begin{array}{}\text{(ODE)}& \dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{\alpha }(t))(t>0)\end{array}$$

In this section we discuss another variant problem, one for which the initial position is constrained to lie in a given set $X_0\subset {\mathbb{R}}^n$ and the final position is also constrained to lie within a given set $X_1\subset {\mathbb{R}}^n$.

So in this model we get to choose the starting point $x^0\in X_0$ in order to maximize

$$\begin{array}{}\text{(P)}& P[\mathit{\alpha }(\cdot )]={\int }_0^{\tau }r(\mathit{x}(t),\mathit{\alpha }(t))\mathrm{d}t,\end{array}$$

where $\tau =\tau [\mathit{\alpha }(\cdot )]$ is the first time we hit $X_1$.

Notation

We will assume that $X_0,X_1$ are in fact smooth surfaces in ${\mathbb{R}}^n$. We let $T_0$ denote the tangent plane to $X_0$ at $x^0$, and $T_1$ the tangent plane to $X_1$ at $x^1$.

Theorem 4.5 (More Transversality Conditions)

Let ${\mathit{\alpha }}^{\ast }(\cdot )$ and ${\mathit{x}}^{\ast }(\cdot )$ solve the problem above, with

$$x^0={\mathit{x}}^{\ast }(0),x^1={\mathit{x}}^{\ast }\left({\tau }^{\ast }\right)$$

Then there exists a function ${\mathit{p}}^{\ast }(\cdot ):[0,{\tau }^{\ast }]\to {\mathbb{R}}^n$, such that $(\mathrm{ODE}),(\mathrm{ADJ})$ and $(\mathrm{M})$ hold for $0\le t\le {\tau }^{\ast }$. In addition,

$$\begin{array}{}\text{(T)}& \{\begin{array}{l}{\mathit{p}}^{\ast }\left({\tau }^{\ast }\right)\text{ is perpendicular to }{T}_1\\ {\mathit{p}}^{\ast }(0)\text{ is perpendicular to }{T}_0\end{array}\end{array}$$

We call $(\mathrm{T})$ the transversality conditions.

Remarks and Intepretations

• If we have $T>0$ fixed and

$$P[\mathit{\alpha }(\cdot )]={\int }_0^{T}r(\mathit{x}(t),\mathit{\alpha }(t))\mathrm{d}t+g(x(T)),$$

then $(\mathrm{T})$ says

$${\mathit{p}}^{\ast }(T)=\mathrm{\nabla }g({\mathit{x}}^{\ast }(T)),$$

in agreement with our earlier form of the terminal/transversality condition.

• Suppose that the surface $X_1$ is the graph $X_1=\{x\mid g_k(x)=0,k=1,\dots ,l\}$. Then $(\mathrm{T})$ says that ${\mathit{p}}^{\ast }\left({\tau }^{\ast }\right)$ belongs to the "orthogonal complement" of the subspace $T_1$. But orthogonal complement of $T_1$ is the span of $\mathrm{\nabla }{g}_k\left(x^1\right)(k=1,\dots ,l)$. Thus

$${\mathit{p}}^{\ast }\left({\tau }^{\ast }\right)=\sum _{k=1}^l{\lambda }_k\mathrm{\nabla }{g}_k\left(x^1\right)$$

for some unknown constants ${\lambda }_1,\dots ,{\lambda }_l$.

4.6 More Applications

4.6.1 Example 1: Distance between Two Sets

As a first and simple example, let

$$\begin{array}{}\text{(ODE)}& \dot{\mathit{x}}(t)=\mathit{\alpha }(t)\end{array}$$

for $A=S^1$, the unit sphere in ${\mathbb{R}}^2:a\in S^1$ if and only if $|a{|}^2=a_1^2+a_2^2=1$. In other words, we are considering only curves that move with unit speed.

We take

$$\begin{array}{rl}P[\mathit{\alpha }(\cdot )]& =-{\int }_0^{\tau }|\dot{\mathit{x}}(t)|\mathrm{d}t=-\text{ the length of the curve }\\ \text{(P)}& & =-{\int }_0^{\tau }\mathrm{d}t=-\text{ time it takes to reach }{X}_1\end{array}$$

We want to minimize the length of the curve and, as a check on our general theory, will prove that the minimum is of course a straight line.

Using the maximum principle

We have

$$\begin{array}{rl}H(x,p,a)& =\mathit{f}(x,a)\cdot \mathit{p}+r(x,a)\\ & =\mathit{a}\cdot \mathit{p}-1=p_1{a}_1+p_2{a}_2-1\end{array}$$

The adjoint dynamics equation $(\mathrm{ADJ})$ says

$$\dot{\mathit{p}}(t)=-{\mathrm{\nabla }}_{x}H(\mathit{x}(t),\mathit{p}(t),\mathit{\alpha }(t))=0$$

and therefore

$$\mathit{p}(t)\equiv \text{ constant }=p^0\ne 0.$$

The maximization principle ( M ) tells us that

$$H(\mathit{x}(t),\mathit{p}(t),\mathit{\alpha }(t))=\underset{a\in S^1}{max}[-1+p_1^0{a}_1+p_2^0{a}_2].$$

The right hand side is maximized by $a^0=\frac{p^0}{\left|p^0\right|}$, a unit vector that points in the same direction of $p^0$. Thus $\mathit{\alpha }(\cdot )\equiv a^0$ is constant in time. According then to $(\mathrm{ODE})$ we have $\dot{\mathit{x}}=a^0$, and consequently $\mathit{x}(\cdot )$ is a straight line.

Finally, the transversality conditions say that

$$\begin{array}{}\text{(T)}& \mathit{p}(0)\perp T_0,\mathit{p}\left(t_1\right)\perp T_1.\end{array}$$

In other words, $p^0\perp T_0$ and $p^0\perp T_1$; and this means that the tangent planes $T_0$ and $T_1$ are parallel.

Now all of this is pretty obvious from the picture, but it is reassuring that the general theory predicts the proper answer.

$◻$

4.6.2 Example 2: Commodity Trading

Next is a simple model for the trading of a commodity, say wheat. We let $T$ be the fixed length of trading period, and introduce the variables

$$\begin{array}{rl}{x}^1(t)& =\text{ money on hand at time }t\\ x^2(t)& =\text{ amount of wheat owned at time }t\\ \alpha (t)& =\text{ rate of buying or selling of wheat }\\ q(t)& =\text{ price of wheat at time }t\text{ (known) }\\ \lambda & =\text{ cost of storing a unit amount of wheat for a unit of time. }\end{array}$$

We suppose that the price of wheat $q(t)$ is known for the entire trading period $0\le t\le T$ (although this is probably unrealistic in practice). We assume also that the rate of selling and buying is constrained:

$$|\alpha (t)|\le M$$

where $\alpha (t)>0$ means buying wheat, and $\alpha (t)<0$ means selling.

Our intention is to maximize our holdings at the end time $T$, namely the sum of the cash on hand and the value of the wheat we then own:

$$\begin{array}{}\text{(P)}& P[\mathit{\alpha }(\cdot )]=x^1(T)+q(T)x^2(T).\end{array}$$

The evolution is

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{rl}{\dot{x}}^1(t)& =-\lambda x^2(t)-q(t)\alpha (t)\\ {\dot{x}}^2(t)& =\alpha (t)\end{array}\end{array}$$

This is a nonautonomous ( $=$ time dependent) case, but it turns out that the Pontryagin Maximum Principle still applies.

Using the maximum principle

What is our optimal buying and selling strategy? First, we compute the Hamiltonian

$$H(x,\mathit{p},t,a)=\mathit{f}\cdot \mathit{p}+r=p_1(-\lambda x_2-q(t)a)+p_{2}a$$

since $r\equiv 0$. The adjoint dynamics read

$$\begin{array}{}\text{(ADJ)}& \{\begin{array}{l}{\dot{p}}^1=0\\ {\dot{p}}^2=\lambda p^1\end{array}\end{array}$$

with the terminal condition

$$\begin{array}{}\text{(T)}& \mathit{p}(T)=\mathrm{\nabla }g(x(T))\end{array}$$

In our case $g(x_1,x_2)=x_1+q(T)x_2$, and hence

$$\begin{array}{}\text{(T)}& \{\begin{array}{l}{p}^1(T)=1\\ p^2(T)=q(T).\end{array}\end{array}$$

We then can solve for the costate:

$$\{\begin{array}{l}{p}^1(t)\equiv 1\\ p^2(t)=\lambda (t-T)+q(T).\end{array}$$

The maximization principle $(\mathrm{M})$ tells us that

$$\begin{array}{}\text{(M)}& \begin{array}{rl}& H(\mathit{x}(t),\mathit{p}(t),t,\alpha (t))\\ =& \underset{|a|\le M}{max}\{p^1(t)(-\lambda x^2(t)-q(t)a)+p^2(t)a\}\\ =& -\lambda p^1(t)x^2(t)+\underset{|a|\le M}{max}\left\{a(-q(t)+p^2(t))\right\}\end{array}\end{array}$$

So

$$\alpha (t)=\{\begin{array}{ll}M& \text{ if }q(t)<p^2(t)\\ -M& \text{ if }q(t)>p^2(t)\end{array}$$

for $p^2(t):=\lambda (t-T)+q(T)$.

Critique

In some situations the amount of money on hand $x^1(\cdot )$ becomes negative for part of the time. The economic problem has a natural constraint $x_2\ge 0$ (unless we can borrow with no interest charges) which we did not take into account in the mathematical model.

4.7 Maximum Principle with State Constraints

We return once again to our usual setting:

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{\alpha }(t))\\ \mathit{x}(0)=x^0\end{array}\text{(P)}& P[\mathit{\alpha }(\cdot )]={\int }_0^{\tau }r(\mathit{x}(t),\mathit{\alpha }(t))\mathrm{d}t\end{array}$$

for $\tau =\tau [\mathit{\alpha }(\cdot )]$, the first time that $x(\tau )=x^1$. This is the fixed endpoint problem.

State Constraints

We introduce a new complication by asking that our dynamics $\mathit{x}(\cdot )$ must always remain within a given region $R\subset {\mathbb{R}}^n$. We will as above suppose that $R$ has the explicit representation

$$R=\{x\in {\mathbb{R}}^n\mid g(x)\le 0\}$$

for a given function $g(\cdot ):{\mathbb{R}}^n\to \mathbb{R}$.

Definition

It will be convenient to introduce the quantity

$$c(x,a):=\mathrm{\nabla }g(x)\cdot \mathit{f}(x,a).$$

Notice that if $\mathit{x}(t)\in \partial R$ for times $s_0\le t\le s_1$, then $c(\mathit{x}(t),\mathit{\alpha }(t))\equiv 0(s_0\le t\le s_1)$. This is so since $\mathit{f}$ is then tangent to $\partial R$, whereas $\mathrm{\nabla }g$ is perpendicular.

Theorem 4.6 (Maximum Principle with State Constraints)

Let ${\mathit{\alpha }}^{\ast }(\cdot ),{\mathit{x}}^{\ast }(\cdot )$ solve the control theory problem above. Suppose also that ${\mathit{x}}^{\ast }(t)\in \partial R$ for $s_0\le t\le s_1$.

Then there exists a costate function ${\mathit{p}}^{\ast }(\cdot ):[s_0,s_1]\to {\mathbb{R}}^n$ such that $(\mathrm{ODE})$ holds. There also exists ${\lambda }^{\ast }(\cdot ):[s_0,s_1]\to \mathbb{R}$ such that for times $s_0\le t\le s_1$ we have

$$\begin{array}{}\text{(ADJ’)}& {\dot{\mathit{p}}}^{\ast }(t)=-{\mathrm{\nabla }}_{x}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))+{\lambda }^{\ast }(t){\mathrm{\nabla }}_{x}c({\mathit{x}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t));\end{array}$$

and

$$\begin{array}{}\text{(M’)}& \left({\mathrm{M}}^{\prime }\right)H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))=\underset{a\in A}{max}\{H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),a)\mid c({\mathit{x}}^{\ast }(t),a)=0\}.\end{array}$$

To keep things simple, we have omitted some technical assumptions really needed for the Theorem to be valid.

Remarks and Intepretations

• Let $A\subset {\mathbb{R}}^m$ be of this form:

$$A=\{a\in {\mathbb{R}}^m\mid g_1(a)\le 0,\dots ,g_s(a)\le 0\}$$

for given functions $g_1,\dots ,g_s:{\mathbb{R}}^m\to \mathbb{R}$. In this case we can use Lagrange multipliers to deduce from $\left({\mathrm{M}}^{\prime }\right)$ that

$$\begin{array}{}\text{(M”)}& {\mathrm{\nabla }}_{a}H({\mathit{x}}^{\ast }(t),{\mathit{p}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))={\lambda }^{\ast }(t){\mathrm{\nabla }}_{a}c({\mathit{x}}^{\ast }(t),{\mathit{\alpha }}^{\ast }(t))+\sum _{i=1}^s{\mu }_i^{\ast }(t){\mathrm{\nabla }}_a{g}_i({\mathit{x}}^{\ast }(t))\end{array}$$

The function ${\lambda }^{\ast }(\cdot )$ here is that appearing in $({\mathrm{ADJ}}^{\prime })$. If ${\mathit{x}}^{\ast }(t)$ lies in the interior of $R$ for say the times $0\le t<s_0$, then the ordinary Maximum Principle holds.

• Jump conditions. In the situation above, we always have

$${\mathit{p}}^{\ast }(s_0-0)={\mathit{p}}^{\ast }(s_0+0)$$

where $s_0$ is a time that ${\mathit{x}}^{\ast }$ hits $\partial R$. In other words, there is no jump in ${\mathit{p}}^{\ast }$ when we hit the boundary of the constraint $\partial R$.

However,

$${\mathit{p}}^{\ast }(s_1+0)={\mathit{p}}^{\ast }(s_1-0)-{\lambda }^{\ast }\left(s_1\right)\mathrm{\nabla }g\left({\mathit{x}}^{\ast }\left(s_1\right)\right)$$

this says there is (possibly) a jump in ${\mathit{p}}^{\ast }(\cdot )$ when we leave $\partial R$.

4.8 More Applications

4.8.1 Example 1: Shortest Distance between Two Points, Avoiding An Obstacle.

What is the shortest path between two points that avoids the disk $B=B(0,r)$, as drawn?

Let us take

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{\mathit{x}}(t)=\mathit{\alpha }(t)\\ \mathit{x}(0)=x^0\end{array}\end{array}$$

for $A=S^1$, with the payoff

$$\begin{array}{}\text{(P)}& P[\mathit{\alpha }(\cdot )]=-{\int }_0^{\tau }|\dot{\mathit{x}}|\mathrm{d}t=-\text{ length of the curve }\mathit{x}(\cdot ).\end{array}$$

We have

$$H(x,\mathit{p},a)=\mathit{f}\cdot \mathit{p}+r=p_1{a}_1+p_2{a}_2-1.$$

Case 1: avoiding the obstacle

Assume $\mathit{x}(t)\notin \partial B$ on some time interval. In this case, the usual Pontryagin Maximum Principle applies, and we deduce as before that

$$\dot{\mathit{p}}=-{\mathrm{\nabla }}_{x}H=0.$$

Hence

$$\begin{array}{}\text{(ADJ)}& \mathit{p}(t)\equiv \text{ constant }=p^0.\end{array}$$

Condition $(\mathrm{M})$ says

$$H(\mathit{x}(t),\mathit{p}(t),\mathit{\alpha }(t))=\underset{a\in S^1}{max}(-1+p_1^0{a}_1+p_2^0{a}_2)$$

The maximum occurs for $\alpha =\frac{p^0}{\left|p^0\right|}$. Furthermore,

$$-1+p_1^0{\alpha }_1+p_2^0{\alpha }_2\equiv 0$$

and therefore $\alpha \cdot p^0=1$. This means that $\left|p^0\right|=1$, and hence in fact $\alpha =p^0$. We have proved that the trajectory $\mathit{x}(\cdot )$ is a straight line away from the obstacle.

Case 2: touching the obstacle

Suppose now $\mathit{x}(t)\in \partial B$ for some time interval $s_0\le t\le s_1$. Now we use the modified version of Maximum Principle, provided by Theorem 4.6.

First we must calculate $c(x,a)=\mathrm{\nabla }g(x)\cdot \mathit{f}(x,a)$. In our case,

$$R={\mathbb{R}}^2-B=\{x\mid x_1^2+x_2^2\ge r^2\}=\{x\mid g:=r^2-x_1^2-x_2^2\le 0\}$$

Then $\mathrm{\nabla }g=(\genfrac{}{}{0ex}{}{-2x_1}{-2x_2})$. Since $\mathit{f}=(\genfrac{}{}{0ex}{}{a_1}{a_2})$, we have

$$c(x,a)=-2a_1{x}_1-2a_2{x}_2$$

Now condition $({\mathrm{ADJ}}^{\prime })$ implies

$$\dot{\mathit{p}}(t)=-{\mathrm{\nabla }}_{x}H+\lambda (t){\mathrm{\nabla }}_{x}c;$$

which is to say,

$$\begin{array}{}\text{(4.6)}& \{\begin{array}{l}{\dot{p}}^1=-2\lambda {\alpha }^1\\ {\dot{p}}^2=-2\lambda {\alpha }^2\end{array}\end{array}$$

Next, we employ the maximization principle $({\mathrm{M}}^{\prime })$. We need to maximize

$$H(\mathit{x}(t),\mathit{p}(t),a)$$

subject to the requirements that $c(\mathit{x}(t),a)=0$ and $g_1(a)=a_1^2+a_2^2-1=0$, since $A=\{a\in {\mathbb{R}}^2\mid a_1^2+a_2^2=1\}$. According to $({\mathrm{M}}^{\prime \prime })$ we must solve

$${\mathrm{\nabla }}_{a}H=\lambda (t){\mathrm{\nabla }}_{a}c+\mu (t){\mathrm{\nabla }}_a{g}_1;$$

that is,

$$\{\begin{array}{l}{p}^1=\lambda (-2x^1)+\mu 2{\alpha }^1\\ p^2=\lambda (-2x^2)+\mu 2{\alpha }^2\end{array}$$

We can combine these identities to eliminate $\mu$. Since we also know that $\mathit{x}(t)\in \partial B$, we have ${\left(x^1\right)}^2+{\left(x^2\right)}^2=r^2$; and also $\mathit{\alpha }={({\alpha }^1,{\alpha }^2)}^{\mathrm{\top }}$ is tangent to $\partial B$. Using these facts, we find after some calculations that

$$\begin{array}{}\text{(4.7)}& \lambda =\frac{p^2{\alpha }^1-p^1{\alpha }^2}{2r}.\end{array}$$

But we also know

$$\begin{array}{}\text{(4.8)}& {\left({\alpha }^1\right)}^2+{\left({\alpha }^2\right)}^2=1\end{array}$$

and

$$H\equiv 0=-1+p^1{\alpha }^1+p^2{\alpha }^2$$

hence

$$\begin{array}{}\text{(4.9)}& p^1{\alpha }^1+p^2{\alpha }^2\equiv 1.\end{array}$$

Solving for the unknowns. We now have the five equations $(4.6)-(4.9)$ for the five unknown functions $p^1,p^2,{\alpha }^1,{\alpha }^2,\lambda$ that depend on $t$. We introduce the angle $\theta$, as illustrated, and note that $\frac{\mathrm{d}}{\mathrm{d}\theta }=r\frac{\mathrm{d}}{\mathrm{d}t}$. A calculation then confirms that the solutions are

$$\begin{array}{c}\{\begin{array}{l}{\alpha }^1(\theta )=-\sin \theta \\ {\alpha }^2(\theta )=\cos \theta \end{array}\\ \lambda =-\frac{k+\theta }{2r}\end{array}$$

and

$$\{\begin{array}{l}{p}^1(\theta )=k\cos \theta -\sin \theta +\theta \cos \theta \\ p^2(\theta )=k\sin \theta +\cos \theta +\theta \sin \theta \end{array}$$

for some constant $k$.

Case 3: approaching and leaving the obstacle

In general, we must piece together the results from Case 1 and Case 2. So suppose now $\mathit{x}(t)\in R={\mathbb{R}}^2-B$ for $0\le t<s_0$ and $\mathit{x}(t)\in \partial B$ for $s_0\le t\le s_1$.

We have shown that for times $0\le t<s_0$, the trajectory $\mathit{x}(\cdot )$ is a straight line. For this case we have shown already that $p=\alpha$ and therefore

$$\{\begin{array}{l}{p}^1\equiv -\cos {\varphi }_0\\ p^2\equiv \sin {\varphi }_0\end{array}$$

for the angle ${\varphi }_0$ as shown in the picture.

By the jump conditions, $\mathit{p}(\cdot )$ is continuous when $\mathit{x}(\cdot )$ hits $\partial B$ at the time $s_0$, meaning in this case that

$$\{\begin{array}{l}k\cos {\theta }_0-\sin {\theta }_0+{\theta }_0\cos {\theta }_0=-\cos {\varphi }_0\\ k\sin {\theta }_0+\cos {\theta }_0+{\theta }_0\sin {\theta }_0=\sin {\varphi }_0\end{array}$$

These identities hold if and only if

$$\{\begin{array}{l}k=-{\theta }_0\\ {\theta }_0+{\varphi }_0=\frac{\pi }{2}\end{array}$$

The second equality says that the optimal trajectory is tangent to the disk $B$ when it hits $\partial B$.

We turn next to the trajectory as it leaves $\partial B$ : see the next picture. We then have

$$\{\begin{array}{l}{p}^1\left({\theta }_1^{-}\right)=-{\theta }_0\cos {\theta }_1-\sin {\theta }_1+{\theta }_1\cos {\theta }_1\\ p^2\left({\theta }_1^{-}\right)=-{\theta }_0\sin {\theta }_1+\cos {\theta }_1+{\theta }_1\sin {\theta }_1.\end{array}$$

Now our formulas above for $\lambda$ and $k$ imply $\lambda \left({\theta }_1\right)=\frac{{\theta }_0-{\theta }_1}{2r}$. The jump conditions give

$$\mathit{p}\left({\theta }_1^{+}\right)=\mathit{p}\left({\theta }_1^{-}\right)-\lambda \left({\theta }_1\right)\mathrm{\nabla }g\left(\mathit{x}\left({\theta }_1\right)\right)$$

for $g(x)=r^2-x_1^2-x_2^2$. Then

$$\lambda \left({\theta }_1\right)\mathrm{\nabla }g\left(\mathit{x}\left({\theta }_1\right)\right)=({\theta }_1-{\theta }_0)(\genfrac{}{}{0ex}{}{\cos {\theta }_1}{\sin {\theta }_1}).$$

Therefore

$$\{\begin{array}{l}{p}^1\left({\theta }_1^{+}\right)=-\sin {\theta }_1\\ p^2\left({\theta }_1^{+}\right)=\cos {\theta }_1,\end{array}$$

and so the trajectory is tangent to $\partial B$. If we apply usual Maximum Principle after $\mathit{x}(\cdot )$ leaves $B$, we find

$$\begin{array}{rl}& p^1\equiv \mathrm{constant}=-\cos {\varphi }_1\\ & p^2\equiv \mathrm{constant}=-\sin {\varphi }_1.\end{array}$$

Thus

$$\{\begin{array}{l}-\cos {\varphi }_1=-\sin {\theta }_1\\ -\sin {\varphi }_1=\cos {\theta }_1\end{array}$$

and so ${\theta }_1-{\varphi }_1=\frac{\pi }{2}$.

Critique

We have carried out elaborate calculations to derive some pretty obvious conclusions in this example. It is best to think of this as a confirmation in a simple case of Theorem 4.6, which applies in far more complicated situations.

$◻$

4.8.2 An Inventory Control Model

Now we turn to a simple model for ordering and storing items in a warehouse. Let the time period $T>0$ be given, and introduce the variables

$$\begin{array}{rl}x(t)& =\text{ amount of inventory at time }t\\ \alpha (t)& =\text{ rate of ordering from manufacturers, }\alpha \ge 0,\\ d(t)& =\text{ customer demand }(\text{ known })\\ \gamma & =\text{ cost of ordering }1\text{ unit }\\ \beta & =\text{ cost of storing }1\text{ unit. }\end{array}$$

Our goal is to fill all customer orders shipped from our warehouse, while keeping our storage and ordering costs at a minimum. Hence the payoff to be maximized is

$$\begin{array}{}\text{(P)}& P[\alpha (\cdot )]=-{\int }_0^T\gamma \alpha (t)+\beta x(t)\mathrm{d}t.\end{array}$$

We have $A=[0,\mathrm{\infty })$ and the constraint that $x(t)\ge 0$. The dynamics are

$$\begin{array}{}\text{(ODE)}& \{\begin{array}{l}\dot{x}(t)=\alpha (t)-d(t)\\ x(0)=x^0>0\end{array}\end{array}$$

Guessing the optimal strategy

Let us just guess the optimal control strategy: we should at first not order anything ( $\alpha =0$ ) and let the inventory in our warehouse fall off to zero as we fill demands; thereafter we should order just enough to meet our demands ( $\alpha =d$ ).

Using the maximum principle

We will prove this guess is right, using the Maximum Principle. Assume first that $x(t)>0$ on some interval $[0,s_0]$. We then have

$$H(x,p,a,t)=(a-d(t))p-\gamma a-\beta x$$

and $(\mathrm{ADJ})$ says $\dot{p}=-{\mathrm{\nabla }}_{x}H=\beta$. Condition $(\mathrm{M})$ implies

$$\begin{array}{rl}H(x(t),p(t),\alpha (t),t)& =\underset{a\ge 0}{max}\{-\gamma a-\beta x(t)+p(t)(a-d(t))\}\\ & =-\beta x(t)-p(t)d(t)+\underset{a\ge 0}{max}\{a(p(t)-\gamma )\}\end{array}$$

Thus

$$\alpha (t)=\{\begin{array}{ll}0& \text{ if }p(t)\le \gamma \\ +\mathrm{\infty }& \text{ if }p(t)>\gamma \end{array}$$

If $\alpha (t)\equiv +\mathrm{\infty }$ on some interval, then $P[\alpha (\cdot )]=-\mathrm{\infty }$, which is impossible, because there exists a control with finite payoff. So it follows that $\alpha (\cdot )\equiv 0$ on $[0,s_0]$ : we place no orders.

According to $(\mathrm{ODE})$, we have

$$\{\begin{array}{l}\dot{x}(t)=-d(t)(0\le t\le s_0)\\ x(0)=x^0\end{array}$$

Thus $s_0$ is first time the inventory hits 0 . Now since $x(t)=x^0-{\int }_0^{t}d(s)\mathrm{d}s$, we have $x\left(s_0\right)=0$. That is, ${\int }_0^{s_0}d(s)\mathrm{d}s=x^0$ and we have hit the constraint. Now use Pontryagin Maximum Principle with state constraint for times $t\ge s_0$

$$R=\{x\ge 0\}=\{g(x):=-x\le 0\}$$

and

$$c(x,a,t)=\mathrm{\nabla }g(x)\cdot f(x,a,t)=(-1)(a-d(t))=d(t)-a.$$

We have

$$\begin{array}{}\text{(M’)}& H(x(t),p(t),\alpha (t),t)=\underset{a\ge 0}{max}\{H(x(t),p(t),a,t)\mid c(x(t),a,t)=0\}\end{array}$$

But $c(x(t),\alpha (t),t)=0$ if and only if $\alpha (t)=d(t)$. Then ( ODE ) reads

$$\dot{x}(t)=\alpha (t)-d(t)=0$$

and so $x(t)=0$ for all times $t\ge s_0$. We have confirmed that our guess for the optimal strategy was right.

4.9 Numerical Solution of the Two Point Boundary Value Problem (TPBVP)

Reference: Chapter 12.6 in Numerical Optimal Control by Moritz Diehl and Sebastien Gros

At first, give the main drawbacks of the indirect approach:

1. it must be possible to eliminate the controls from the problem by algebraic manipulations, which is not always straightforward or might even be impossible;
2. the optimal controls might be a discontinuous function of $x$ and $p$, such that the BVP is possibly given by a non-smooth differential equation;
3. the differential equation might become very nonlinear and unstable and not suitable for a forward simulation.

All these issues of the indirect approach can partially be addressed, and most importantly, it offers an exact and elegant characterization of the solution of optimal control problems in continuous time.

In this section we address the question of how we can compute a solution of the boundary value problem (BVP) in the indirect approach. The remarkable observation is that the only non-trivial unknown is the initial value for the adjoints, $\mathit{p}(0)$. Once this value has been found, the complete optimal trajectory can in principle be recovered by a forward simulation of the combined differential equation. Let us first recall that the BVP that we want to solve is given as

$$\begin{array}{}\text{(4.1a)}& b_0=\mathit{x}(0)-{\mathit{x}}^0& =0,\text{(4.1b)}& b_T=\mathit{p}(T)-\mathrm{\nabla }g(\mathit{x}(T))& =0,\text{(4.1c)}& \dot{\mathit{x}}(t)-{\mathrm{\nabla }}_{\mathit{p}}{H}^{\ast }(\mathit{x}(t),\mathit{p}(t))& =0,t\in [0,T],\text{(4.1d)}& \dot{\mathit{p}}(t)+{\mathrm{\nabla }}_{\mathit{x}}{H}^{\ast }(\mathit{x}(t),\mathit{p}(t))& =0,t\in [0,T].\end{array}$$

Using the shorthands

$$\mathit{y}=\left[\begin{array}{c}\mathit{x}\\ \mathit{p}\end{array}\right],\mathit{\phi }(\mathit{y})=\left[\begin{array}{c}{\mathrm{\nabla }}_{\mathit{p}}{H}^{\ast }(\mathit{x}(t),\mathit{p}(t))\\ -{\mathrm{\nabla }}_{\mathit{x}}{H}^{\ast }(\mathit{x}(t),\mathit{p}(t))\end{array}\right]$$

and

$$\mathit{b}(\mathit{y}(0),\mathit{y}(T),{\mathit{x}}^0)=\left[\begin{array}{c}{b}_0(\mathit{y}(0),{\mathit{x}}^0)\\ b_T(\mathit{y}(T))\end{array}\right]$$

the system of equations can be summarized as:

$$\begin{array}{}\text{(4.2a)}& \mathit{b}(\mathit{y}(0),\mathit{y}(T),{\mathit{x}}^0)& =0,\text{(4.2b)}& \dot{\mathit{y}}(t)-\mathit{\phi }(\mathit{y}(t))& =0,t\in [0,T].\end{array}$$

This BVP has the $2n_{\mathit{x}}$ differential equations $\dot{\mathit{y}}=\phi$, and the $2n_{\mathit{x}}$ boundary conditions $\mathit{b}$ and is therefore usually well-defined. We detail here three approaches to solve this TPBVP numerically, single shooting, collocation, and multiple shooting.

4.9.1 Single shooting

Single shooting starts with the following idea: for any guess of the initial value ${\mathit{p}}_0$, we can use a numerical integration routine in order to obtain the state-costate trajectory as a function of ${\mathit{p}}_0,{\mathit{x}}^0$, i.e. $\mathit{y}(t,{\mathit{p}}_0,{\mathit{x}}^0)$ for all $t\in [0,T]$. This is visualized in Figure 4.9.1. The result is that the differential equation $(4.2b)$ is by construction already satisfied, as well as the initial boundary condition $(4.1a)$. Thus, we only need to enforce the boundary condition $(4.1b)$, which we can do using the terminal trajectory value $\mathit{y}(T,{\mathit{p}}_0,{\mathit{x}}^0)$ :

$$\underset{=:R_T\left({\mathit{p}}_0\right)}{\underbrace{b_T\left(\mathit{y}(T,{\mathit{p}}_0,{\mathit{x}}^0)\right)}}=0.$$

For nonlinear dynamics $\mathit{\phi }$, this equation can generally not be solved explicitly. We then use the Newton's method, starting from an initial guess ${\mathit{p}}_0$, and iterating to the solution, i.e. we iterate

$$\begin{array}{}\text{(4.3)}& {\mathit{p}}_0^{k+1}={\mathit{p}}_0^k-t_k{\left(\frac{\partial R_T}{\partial {\mathit{p}}_0}\left({\mathit{p}}_0^k\right)\right)}^{-1}{R}_T\left({\mathit{p}}_0^k\right).\end{array}$$

for some adequate step-size $t_k\in ]0,1]$. It is important to note that in order to evaluate $\frac{\partial R}{\partial {\mathit{p}}_0}\left({\mathit{p}}_0^k\right)$ we have to compute the ODE sensitivities $\frac{\partial \mathit{y}(T,y_0)}{\partial {\mathit{p}}_0}$.

In some cases, as said above, the forward simulation of the combined ODE might be an ill-conditioned problem so that single shooting cannot be employed. Even if the forward simulation problem is well-defined, the region of attraction of the Newton iteration on $R_T\left({\mathit{p}}_0\right)=0$ can be very small, such that a good guess for ${\mathit{p}}_0$ is often required. However, such a guess is typically unavailable. In the following example, we illustrate these observation on a simple optimal control problem.

Example 4.9.1

We consider the optimal control problem:

$$\begin{array}{}\text{(4.4)}& \begin{array}{rl}\underset{\mathit{x}(.),u(.)}{\min }& {\int }_0^T{x}_1(t{)}^2+10x_2(t{)}^2+u(t{)}^2\mathrm{d}t\\ \text{ subject to }& {\dot{x}}_1(t)=x_1(t)x_2(t)+u(t),& x_1(0)=0\\ & {\dot{x}}_2(t)=x_1(t),& x_2(0)=1\end{array}\end{array}$$

with $T=5$. This example does not hold a terminal cost or constraints, such that the terminal condition reads as $R_T\left({\mathit{p}}_0\right)=\mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)=0$. The state-costate trajectory at the solution is displayed in Figure 4.9.1. It is then interesting to build the function ${\mathit{p}}_0\mapsto \mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)$ for various values of ${\mathit{p}}_0$, see Figure 4.9.2. This function is very nonlinear, making it difficult for the Newton iteration to find the co-states' initial value ${\mathit{p}}_0$ resulting in $\mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)=0$. More specifically, the Newton iteration (full steps or reduced steps) converges only for a specific set of initial guesses ${\mathit{p}}_0^0$ provided to the iteration (4.3), see Figure 4.9.3.

Figure 4.9.1: Illustration of the state and co-state trajectories for example 4.9.1 at the solution delivering $\mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)=0$ for $T=5$, note that ${\lambda }_{\ast }$ in the image is actually $p_{\ast }$.

Figure 4.9.2: Illustration of the map ${\mathit{p}}_0\mapsto \mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)$, in the form of level curves for $T=5$. The black dot represents the solution of the TPBVP problem, where $R_T\left({\mathit{p}}_0\right)=\mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)=0$. One can observe that the map is very nonlinear, such that the Newton method can struggle to converge to the solution ${\mathit{p}}_0$ of the TPBVP, unless a very good initial guess is provided. Note that ${\lambda }_{\ast }$ in the image is actually $p_{\ast }$

Figure 4.9.3: Illustration of the region of convergence of the Newton iteration $(4.3)$ for problem $(4.4)$ (in black, with full Newton steps on the left-hand side graph and with reduced steps on the right-hand side graph). Here we note ${\mathit{p}}_{0,1},{\mathit{p}}_{0,2}$ the initial guess provided to the Newton iteration. The grey dots (at $(3.22,8.48)$) depict the solution to the TPBVP. Only a fairly small, disconnected and highly convoluted set of initial guess for the co-states initial conditions leads to a convergence of the Newton iteration. Note that ${\lambda }_{\ast }$ in the image is actually $p_{\ast }$

A crucial observation that will motivate an alternative to the single-shooting approach is illustrated in Figure 4.9.4, where the map ${\mathit{p}}_0\mapsto \mathit{p}(t,{\mathit{p}}_0,{\mathit{x}}^0)$ is displayed for the integration times $t=3$ and $t=4$. The crucial observation here is that the map is fairly linear up to $t=3$, and becomes increasingly nonlinear for larger integration times. This observation is general and motivates the core idea behind the alternatives to single shooting, namely that integration shall never be performed over long time intervals, so as to avoid creating strongly nonlinear functions in the TPBVP.

Figure 4.9.4: Illustration of the map ${\mathit{p}}_0\mapsto \mathit{p}(t,{\mathit{p}}_0,{\mathit{x}}^0)$, in the form of level curves for different times $t$. The black dot represents the solution of the TPBVP problem, where $\mathit{p}(T,{\mathit{p}}_0,{\mathit{x}}^0)=0$. One can observe that the map is close to linear for "small" integration times $t$ (upper graphs, where $t=3$ ), and becomes increasingly nonlinear as the integration time increases (lower graph, where $t=4$ ), until it reaches the final time $T=5$, see Figure 4.9.2. This observation is general, and holds for most problems.

4.9.2 Multiple shooting

The nonlinearity of the integration map ${\mathit{p}}_0\mapsto \mathit{y}(t,{\mathit{p}}_0,{\mathit{x}}^0)$ for long integration times $t$ motivates the "breaking down" of the full integration in small pieces, so as to avoid creating very nonlinear map in the TPBVP conditions. The idea is originally due to Osborne, and is based on dividing the time interval $[0,T]$ into (typically uniform) shooting intervals $[t_k,t_{k+1}]\subset [0,T]$, where the most common choice is $t_k=k\frac{T}{N}$. Let us then frame the integration over a short time interval $[t_k,t_{k+1}]$ with initial value ${\mathit{\beta }}_k$ as the function ${\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)$, defined as:

$$\begin{array}{}\text{(4.5a)}& {\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)=\mathit{y}\left(t_{k+1}\right)\text{ where }\text{(4.5b)}& \dot{\mathit{y}}(t)-\mathit{\phi }(\mathit{y}(t))=0,t\in [t_k,t_{k+1}]\text{ and }\mathit{y}\left(t_k\right)={\mathit{\beta }}_k\end{array}$$

for $k=0,\dots ,N-1$. We then rewrite the TPBVP conditions (4.2) as:

$$\begin{array}{}\text{(4.6a)}& b({\mathit{\beta }}_0,{\mathit{\beta }}_N,{\mathit{x}}^0)& =0,& & & & \text{ (boundary conditions) }\text{(4.6b)}& {\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)-{\mathit{\beta }}_{k+1}& =0,& & k=0,\dots ,N-1.& & \text{(continuity conditions)}\end{array}$$

One can then rewrite the conditions (4.6) altogether as the function:

$$\begin{array}{}\text{(4.7)}& R_{\mathrm{MS}}(\mathit{\beta },{\mathit{x}}^0)=0\end{array}$$

where we note $\mathit{\beta }=({\mathit{\beta }}_0,\dots ,{\mathit{\beta }}_N)$. A Newton iteration can be then deployed on $(4.7)$ to find the variables $\mathit{\beta }$, it reads as:

$$\begin{array}{}\text{(4.8)}& {\mathit{\beta }}^{k+1}={\mathit{\beta }}^k-t_k{\left(\frac{\partial R_{\mathrm{MS}}}{\partial \mathit{\beta }}({\mathit{\beta }}^k,{\mathit{x}}^0)\right)}^{-1}{R}_{\mathrm{MS}}({\mathit{\beta }}^k,{\mathit{x}}^0).\end{array}$$

for some step-size $t_k\in ]0,1]$. We illustrate the Multiple-Shooting approach in the following example.

Example 4.9.2

We consider the optimal control problem $(4.4)$ from Example 4.9.1 with $T=5$. If we denote ${\mathit{\beta }}_k=({\mathit{x}}_k,{\mathit{p}}_k)$, the boundary conditions for this example then become:

$$x_0=\left[\begin{array}{l}0\\ 1\end{array}\right],{\mathit{p}}_N=0.$$

We illustrate the Multiple-Shooting procedure (12.14) in Figure 4.9.5 for $N=5$.

Figure 4.9.5: Illustration of the state and co-state trajectories for problem $(4.4)$ during the multiple-shooting iterations $(4.8)$, such that the conditions ${\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)-{\mathit{\beta }}_{k+1}=0$ are not yet fulfilled. Here, the discrete times $t_k$ are depicted as grey dashed lines, the discrete state-costates ${\mathit{\beta }}_k=(x_{k,1},x_{k,2},{\mathit{p}}_{k,1},{\mathit{p}}_{k,2})$ are depicted as black dots, and the resulting integrations ${\mathrm{\Phi }}_k=({\mathrm{\Phi }}_{k,1}^{\mathit{x}},{\mathrm{\Phi }}_{k,2}^{\mathit{x}},{\mathrm{\Phi }}_{k,1}^{\mathit{p}},{\mathrm{\Phi }}_{k,2}^{\mathit{p}})$ are depicted as white dots. The black curves represent the state-costate trajectories on the various time intervals $[t_k,t_{k+1}]$. At the solution (12.14), the conditions ${\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)={\mathit{\beta }}_{k+1}$ are enforced for $k=1,\dots ,N-1$, such that the black and white dots coincide on each discrete time $t_k$.

One ought to observe that the time intervals $[t_k,t_{k+1}]$ are of size $\frac{T}{N}$, and hence get shorter as $N$ increases. Because one can "control" the length of the time interval over which the integration is performed via $N$, and because the functions ${\mathrm{\Phi }}_k\left({\mathit{\beta }}_k\right)-{\mathit{\beta }}_{k+1}$ become less nonlinear as the length of the time interval decreases, one can make them "arbitrarily" linear by increasing $N$. It follows that a sufficiently large $N$ typically allows one to solve the Multiple-Shooting conditions $(4.6)$ using a Newton iteration even if no good initial guess is available.

Figure 4.9.6: Illustration of sparsity pattern of the Jacobian matrix $\frac{\partial R_{\mathrm{MS}}}{\partial \mathit{\beta }}$ in the Newton iteration $(4.8)$ for the optimal control problem $(4.4)$ approached via indirect multiple-shooting, for Example 4.9.2. Here we use $N=5$. One can readily observe that the Jacobian matrix is sparse and highly structured. This structure arises via organising the algebraic conditions $(4.7)$ and the variables $\mathit{\beta }$ in time (i.e. in the order $k=0,\dots ,N$ ). Note that here the last variables ${\mathit{\beta }}_N$ where eliminated using the equality ${\mathit{\beta }}_N={\mathrm{\Phi }}_{N-1}\left({\mathit{\beta }}_{k-1}\right)$. In the specific case of Example 4.9.2, the elimination has no impact on the Newton iteration because the boundary conditions $b({\mathit{\beta }}_0,{\mathit{\beta }}_N,{\mathit{x}}^0)$ are linear.

It is important to observe that the set of algebraic conditions $(4.7)$ holds a large number of variables $\mathit{\beta }$, such that the Newton iteration $(4.8)$ is deployed using large Jacobian matrices $\frac{\partial R_{\mathrm{MS}}}{\partial \mathit{\beta }}$. However, these matrices are sparse, and if the algebraic conditions and variables are adequately organised, they are highly structured (see Figure 4.9.6), such that their factorization can be performed efficiently.

The second alternative to single-shooting is the object of the next section, and can be construed as an extreme case of Multiple-Shooting. We detail this next.

4.9.3 Collocation & Pseudo-spectral methods

The second alternative approach to single shooting is to use simultaneous collocation or Pseudo-spectral methods. As we will see next, the two approaches are fairly similar. The key idea behind these methods is to introduce all the variables involved in processing the integration of the dynamics, and the related algebraic conditions into the set of algebraic equations to be processed.

The most common implementation of this idea is based on the Orthogonal Collocation method.

We consider the collocation-based integration of the state-costate dynamics on a time interval $[t_k,t_{k+1}]$ starting from the initial value ${\mathit{\beta }}_k$, as described in equation

$${\mathit{c}}_k({\mathit{v}}_k,t_{k,i},{\mathit{\beta }}_k)=\left[\begin{array}{c}{v}_{k,0}-{\mathit{\beta }}_k\\ {\dot{\mathit{p}}}_k(t_{k,1},{\mathit{v}}_k)-\mathit{f}({\mathit{v}}_{k,1},t_{k,i})\\ ⋮\\ {\dot{\mathit{p}}}_k(t_{k,d},{\mathit{v}}_k)-\mathit{f}({\mathit{v}}_{k,d},t_{k,i})\end{array}\right]=0.$$

where $\mathit{v}$ is the coefficients of polynominal.

The integration is then based on solving a set of collocation equations:

$$\begin{array}{}\text{(4.10a)}& {\mathit{v}}_{k,0}& ={\mathit{\beta }}_k\text{(4.10b)}& \dot{p}(t_{k,i},{\mathit{v}}_k)& =\mathit{\phi }({\mathit{v}}_{k,i},t_{k,i}),i=1,\dots ,d\end{array}$$

for $k=0,\dots ,N-1$, where $t_{k,i}\in [t_k,t_{k+1}]$ for $i=0,\dots ,d$, and the variables ${\mathit{v}}_k\in {\mathbb{R}}^{2n_{\mathit{x}}\cdot (d+1)}$ hold the discretisation of the continuous state-costates dynamics. The TPBVP discretised using orthogonal collocation then holds the variables ${\mathit{v}}_{k,i}$ and ${\mathit{\beta }}_k$ for $k=0,\dots ,N-1$ and $i=1,\dots ,d$, and the following constraints:

$$\begin{array}{}\text{(4.11a)}& \mathit{b}({\mathit{\beta }}_0,{\mathit{\beta }}_N,{\mathit{x}}^0)=0,& & \text{ (boundary condition), }\text{(4.11b)}& \mathit{p}(t_{k+1},{\mathit{v}}_k)-{\mathit{\beta }}_{k+1}=0,& & \text{ (continuity condition), }\text{(4.11c)}& {\mathit{v}}_{k,0}-{\mathit{\beta }}_k=0,& & \text{ (initial values), }\text{(4.11d)}& \dot{\mathit{p}}(t_{k,i},{\mathit{v}}_k)-\mathit{\phi }({\mathit{v}}_{k,i},t_{k,i})=0,& & \text{ (dynamics). }\end{array}$$

One can observe that equations $(4.11b)$ and $(4.11c)$ are linear, while equation $(4.11d)$ is nonlinear when the dynamics are nonlinear. One can also observe that the variables ${\mathit{\beta }}_{0,\dots ,N-1}$ can actually be eliminated from (4.11), to yield a slightly more compact set of equation, with $k=0,\dots ,N-1$ and $i=1,\dots ,d$ :

$$\begin{array}{}\text{(4.12a)}& \mathit{b}({\mathit{v}}_{k,0},{\mathit{v}}_{N,0},{\mathit{x}}_0)& =0,& & \text{ (boundary condition), }\text{(4.12b)}& \mathit{p}(t_{k+1},{\mathit{v}}_k)-{\mathit{v}}_{k+1,0}& =0,& & \text{ (continuity condition), }\text{(4.12c)}& \dot{\mathit{p}}(t_{k,i},{\mathit{v}}_k)-\mathit{\phi }({\mathit{v}}_{k,i},t_{k,i})& =0,& & \text{ (dynamics). }\end{array}$$

This elimination does not modify the behavior of the Newton iteration. We can gather the algebraic conditions $(4.12)$ and the variables ${\mathit{\beta }}_k,{\mathit{v}}_k$ in the compact form

$$\begin{array}{}\text{(4.13)}& R_{\mathrm{IC}}(\mathit{w},{\mathit{x}}^0)=0,\end{array}$$

where $\mathit{w}=\{{\mathit{v}}_{0,0},\dots ,{\mathit{v}}_{0,d},\dots ,{\mathit{v}}_{N-1,0},\dots ,{\mathit{v}}_{N-1,d},{\mathit{v}}_{N,0}\}$. A Newton iteration can be then deployed on $(4.13)$ to find the variables $\mathit{w}$, it reads as

$$\begin{array}{}\text{(4.14)}& {\mathit{w}}^{k+1}={\mathit{w}}^k-t_k{\left(\frac{\partial R_{\mathrm{IC}}}{\partial w}({\mathit{w}}^k,{\mathit{x}}^0)\right)}^{-1}{R}_{\mathrm{IC}}({\mathit{w}}^k,{\mathit{x}}^0),\end{array}$$

for some step-size $t_k\in [0,1]$. We illustrate the indirect collocation approach in the following example.

Figure 4.9.7: Illustration of the state and co-state trajectories for example 4.9.1 using the orthogonal collocation approach with $N=20$. The grey curves display the state-costate trajectories after the first full Newton step of $(4.14)$, while the black curves report the state-costate trajectories at convergence. The discrete times $t_k$ are depicted as grey dashed lines, the discrete state-costates on the time grid $t_{k,i}$ are depicted as dots. Note that the continuity conditions $(4.11b)$ in the collocation method are linear in the variables $\mathit{w}$, such that the trajectories are continuous after the first full Newton step (hence the grey curves are continuous, even though the problem is not solved yet).

Example 4.9.3

We consider the optimal control problem from Example 4.9.1 with $T=5$. We illustrate the Orthogonal Collocation procedure $(4.11)$ in Figure 4.9.7 for $N=10$. The sparsity pattern of the Jacobian matrix $\frac{\partial R_{\text{IC }}}{\partial \mathit{w}}$ from the Newton iteration $(4.14)$ is illustrated in Figure 4.9.8. The variables and constraints were ordered with respect to time. Even though it is large, the complexity of forming factorisations of the Jacobian matrix $\frac{\partial R_{IC}}{\partial \mathit{w}}$ is limited as it is sparse and highly structured.

Figure 4.9.8: Illustration of the sparsity structure for the Jacobian $\frac{\partial R_{\mathrm{IC}}}{\partial w}$ in the Newton iteration $(4.14)$

Pseudo-spectral methods

Pseudo-spectral methods deploy a very similar approach to the one described here, to the exception that they skip the division of the time interval $[0,T]$ into subintervals $[t_k,t_{k+1}]$, and use a single set of basis functions spanning the entire time interval $[0,T]$. Pseudo-spectral methods for the TPBVP problem $(4.2)$ can then be framed as:

$$\begin{array}{}\text{(4.15a)}& \mathit{b}(p(0,\mathit{v}),\mathit{p}(T,\mathit{v}),{\mathit{x}}^0)& =0,\text{(4.15b)}& \dot{\mathit{p}}(t_k,\mathit{v})-\mathit{\phi }(\mathit{p}(t_k,\mathit{v}),t_k)& =0,i=k,\dots ,n\end{array}$$

where $t_k\in [0,T]$, and the variables $v\in {\mathbb{R}}^{n_{\mathit{x}}\cdot n}$ hold the discretization of the continuous dynamics. Because they attempt to capture the state trajectories in a single function $p(t,v)$, with $t\in [0,T]$, the Newton iteration solving constraints $(4.15)$ generally holds a dense Jacobian matrix, for which structure-exploiting linear algebra is generally inefficient.