%\setcounter{chapter}{3}
\chapter{Optimal Control}\label{chpt:lqr}
\begin{flushright}
\parbox{5.5in}{Nobel laureate H. A. Simon revolutionized management by identifying a managerial phenomenon that causes decision-making failures. Simon~(1957) claimed that executives, engineers, and decision makers are trying to become optimizers while making decisions.
\begin{center}
------------------------
\end{center}
Optimizer: {\em Executive, engineer, or decision maker who wants to make the best possible decision, without consideration of time constraints.}
\begin{center}
------------------------
\end{center}

%To achieve the best possible decision, we must generate all the alternatives, gather all the information, and build a model that will evaluate the alternatives and choose the best among them. All this requires time, effort, and money. In the real world, there is no limit to the number of alternatives we can evaluate.

Simon proposes an alternative approach and suggests that decision makers behave as ``satisficres,'' that is, aspire to reach a ``satisfactory solution,'' a ``good enough solution,'' and not necessarily a perfect one (the ``optimal'' one).
\begin{center}
------------------------
\end{center}
Satisficer: {\em An executive or decision maker who is satisfied with a reasonable solution that will significantly improve the system and does not strive for the optimal solution.}
\begin{center}
------------------------
\end{center}

The satisficer sets a level of aspiration, a threshold he aspires to achieve. The objective is no longer to maximize or minimize some performance measure, but to achieve a solution that will improve the measure beyond the predefined level of aspiration.---Ronen and Pass~\cite[pp.15--16]{ronen08}}
\end{flushright}

\bigskip

One of the essential elements of the control problem\index{control problem!essential elements}, discussed in Section~\ref{sec:whatisasys}, is a means of testing the performance of any proposed control law\index{control law!performance of}. Whenever we use the term ``best'' or ``optimal''\index{optimal!control} to describe the effectiveness of a given control strategy\index{control strategy!optimal}, we do so with respect to some numerical index of performance called the performance index. We assume that the value of the performance index decreases as the quality of the given admissible control law increases. The admissible controller that ensures the completion of the system objective and at the same time minimizes the performance index is called an {\em optimal controller\/}\index{optimal controller} for the system. In this chapter, we discuss various methods of constructing optimal controllers.
\section{Performance Indices}
Constructing a performance index, that is, choosing a means to measure the system performance, can be considered as a part of the system modeling\index{system!modeling}. Here, we discuss some typical choices of performance indices.

First, suppose that the objective is to control a dynamical system modeled by the equations,
\begin{eqnarray}
\dot{\b x}(t) &=& \b A\b x(t) + \b B\b u(t),\quad \b x(t_0)=\b x_0\label{eq:staeqfoperf}\\
\b y(t)&=& \b C\b x(t),
\label{eq:outforperfwq}
\end{eqnarray}
on a fixed interval $[t_0,t_f]$ so that the components of the state vector are ``small''. A suitable performance index to be minimized would be
\begin{equation}
J_1=\int_{t_0}^{t_f}\b x(t)^{\top}\b x(t) dt.
\label{eq:xnormperin}
\end{equation}
Obviously, if $J_1$ is small, then the state vector norm, $\|\b x(t)\|$, is small in the sense of the above performance index.

If the objective is to control the system so that the components of the output, $\b y(t)$, are to be small, then we could use the performance index\index{performance index}
\begin{eqnarray}
J_2&=&\int_{t_0}^{t_f}\b y(t)^{\top}\b y(t) dt\nonumber\\
&=& \int_{t_0}^{t_f} \b x(t)^{\top}\b C^{\top}\b C\b x(t) dt\nonumber\\
&=& \int_{t_0}^{t_f}\b x(t)^{\top}\b Q\b x(t) dt,
\label{eq:outperin}
\end{eqnarray}
where the weight matrix $\b Q=\b C^{\top}\b C$ is symmetric positive semi-definite.

In the case when we wish to control the system in such a manner that the components of the input, $\b u(t)$, are ``not too large'', a suitable performance index to be minimized is
\begin{equation}
J_3=\int_{t_0}^{t_f}\b u(t)^{\top}\b u(t) dt,
\label{eq:miunoonpi}
\end{equation}
or
\begin{equation}
J_4=\int_{t_0}^{t_f}\b u(t)^{\top}\b R\b u(t) dt,
\label{eq:miucomin}
\end{equation}
where the weight matrix $\b R$ is symmetric positive definite. There is no loss of generality in assuming the weight matrix $\b R$ to be symmetric. If $\b R$ was not symmetric, we could represent the quadratic term $\b u^{\top}\b R\b u$ equivalently as
\[
\b u^{\top}\b R\b u=\b u^{\top}\left(\frac{\b R+\b R^{\top}}{2} \right)\b u,
\]
where the matrix $\frac{1}{2}\left(\b R+\b R^{\top} \right)$ is symmetric.

As pointed out by Owens~\cite[page 186]{owens81}, we cannot simultaneously minimize the performance indices~(\ref{eq:xnormperin}) and~(\ref{eq:miunoonpi}) because minimization of~(\ref{eq:xnormperin}) requires large control signals, while minimization of~(\ref{eq:miunoonpi}) requires small control signals. To solve the dilemma, we could compromise between the two conflicting objectives\index{conflicting objectives!convex combination} by minimizing the performance index that is a convex combination of $J_1$ and $J_3$\index{performance index!convex combination},
\begin{eqnarray}
J&=&\lambda J_1 + (1-\lambda)J_3\nonumber\\
&=&\int_{t_0}^{t_f}\left( \lambda\b x(t)^{\top}\b x(t)+(1-\lambda) \b u(t)^{\top}\b u(t)\right)dt,
\label{eq:cocomofpis}
\end{eqnarray}
where $\lambda$ is a parameter in the range $[0,1]$. If $\lambda=1$, then $J=J_1$ and if $\lambda=0$, then $J=J_3$. By trial and error, we could select $\lambda$ from the interval $[0,1]$ to compromise between the two extremes. A generalization of the performance index~(\ref{eq:cocomofpis}) is
\[
J=\int_{t_0}^{t_f}\left(\b x(t)^{\top}\b Q\b x(t)+\b u(t)^{\top}\b R\b u(t) \right) dt.
\]

In certain applications, we may wish the final state $\b x(t_f)$ to be as close as possible to $\b 0$. Then, a possible performance measure to be minimized is
\begin{equation}
\b x(t_f)^{\top}\b F\b x(t_f),
\label{eq:termstsperi}
\end{equation}
where $\b F$ is a symmetric positive definite matrix\index{symmetric positive definite matrix}\index{matrix!symmetric positive definite}.

We can combine the performance measures~(\ref{eq:outperin}), (\ref{eq:miucomin}), and~(\ref{eq:termstsperi}) when our control aim\index{control aim!control not too large} is to keep the state ``small''\index{control aim!state ``small}, the control ``not too large'', and the final state\index{final state!near to $\b 0$} as near to $\b 0$ as possible. The resulting performance index is
\begin{equation}
J=\frac{1}{2}\b x(t_f)^{\top}\b F\b x(t_f)+\frac{1}{2}\int_{t_0}^{t_f}\left(\b x(t)^{\top}\b Q\b x(t)+\b u(t)^{\top}\b R\b u(t) \right) dt,
\label{eq:combperiwoh}
\end{equation}
where the factor $1/2$ is placed to simplify subsequent algebraic manipulations. Minimizing~(\ref{eq:combperiwoh}) subject to~(\ref{eq:staeqfoperf}) is called the {\em linear quadratic regulator\/}\index{linear quadratic regulator (LQR)} problem, or the LQR problem\index{LQR|see{linear quadratic regulator}}  for short.

In some cases, the controller's goal is to force the system state to track a desired state trajectory\index{desired state trajectory}, $\b x_d(t)$, throughout the interval $[t_0,t_f]$ while maintaining the deviations of the actual state $\b x(t)$ ``small'' from the desired trajectory with the control effort $\b u(t)$ ``not too large'' and the final state $\b x(t_f)$ being as near as possible to some desired state $\b x_d(t_f)$. A suitable performance index to be minimized in this case might be
\begin{eqnarray}
J&=&\frac{1}{2}\left(\b x(t_f)-\b x_d(t_f)\right)^{\top}\b F\left(\b x(t_f)-\b x_d(t_f) \right)\nonumber\\
&{}&\quad+\frac{1}{2}\int_{t_0}^{t_f}\left(\left(\b x(t)-\b x_d(t)\right)^{\top}\b Q \left(\b x(t)-\b x_d(t)\right)+\b u(t)^{\top}\b R\b u(t) \right) dt.
\label{eq:trackperfin}
\end{eqnarray}

\section{A Glimpse at the Calculus of Variations}\label{sec:calcofvar}
The goal of this section is to provide a brief introduction to the subject of the {\em calculus of variations}\index{calculus of variations}. The calculus of variations deals with minimization or maximization of {\em functionals}. A functional\index{functional} is a mapping or a transformation that depends on one or more functions and the values of the functional\index{functional values!numbers} are numbers. Examples of functionals are performance indices that we introduced in the previous section. We now give two more examples of functionals.
\begin{example} \rm
The length $l$ of the arc connecting two given points $A(x_0,y_0)$ and $B(x_1,y_1)$ in the $x$-$y$ plane, as depicted in Figure~\ref{fig:arc_len}, is an example of a functional.
\begin{figure}
\centering
\includegraphics[scale=0.5]{arc_len.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/arc_len.bmp}
%\hskip1in\special{bmp:c:/mybook/arc_len.bmp x=4in y=3in}
\caption{The length of arcs joining two given points in the plane is an example of a functional.}
\label{fig:arc_len}
\end{figure}
We know that
\[
l=l(y(x))=\int_{x_0}^{x_1}\sqrt{1+\left(dy(x)/dx \right)^2}dx,
\]
subject to the constraints $y(x_0)=y_0$ and $y(x_1)=y_1$.
\end{example}
\markex
\begin{example} \rm
Another example of a functional is the area of the surface of revolution\index{surface of revolution!example of a functional}\index{functional example!surface of revolution!example of a functional} formed by a surface connecting two circular disks as illustrated in Figure~\ref{fig:surf_rev}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{surf_rev.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/surf_rev.bmp}
%\hskip1in\special{bmp:c:/mybook/surf_rev.bmp x=4in y=3in}
\caption{The area of the surface of revolution connecting two circular disks is an example of a functional.}
\label{fig:surf_rev}
\end{figure}
In particular, the surface of revolution can be formed by a soap film\index{surface of revolution!formed by a soap film} clinging to two circular disks. When a wire circle is dipped into a soap solution and withdrawn, a circular disk of soap film bounded by the circle is formed. If a second smaller circle is made to touch the disk and then moved away, the two circles will be joined by a surface of film that is a surface of revolution in the particular case when the circles are parallel and have their centers on the same axis perpendicular to their planes. If we denote the arc of intersection with the $x$-$y$ plane by $y=y(x)$, then the area of the surface of revolution is
\[
v=v\left(y(x)\right)=2\pi \int_{x_0}^{x_1} y(x) \sqrt{1+\left(dy(x)/dx \right)^2}dx,
\]
subject to the constraints $y(x_0)=y_0$, $y(x_1)=y_1$. One may then want to minimize $v=v\left( y(x)\right)$ to find $y=y(x)$ resulting in the surface of revolution\index{surface of revolution!area of} of minimum area.
\end{example}
\markex
The problem of optimization of functionals dates back to Johann Bernoulli. In 1696 he published a letter in which he posed the problem, which is now known as the {\em brachistochrone problem}\index{brachistochrone problem}\index{brachistochrone!from Greek; shortest time}. Brachistochrone translates from Greek as the shortest time. The brachistochrone problem can be described as follows. A particle of mass $m$ is to slide down along a frictionless track connecting two given points. The motion of the particle is to be constrained only by gravity and the shape of the curve of the track. Find the curve of the track that minimizes the time travel between the two given points\index{brachistochrone problem!minimize time travel between two points}. We illustrate the problem in Figure~\ref{fig:brachist}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.55]{brachist.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/brachist.bmp}
%\hskip1in\special{bmp:c:/mybook/brachist.bmp x=4in y=3in}
\caption{Illustration of the brachistochrone problem.}
\label{fig:brachist}
\end{figure}
The brachistochrone problem was solved by Johann Bernoulli himself,\index{brachistochrone problem!posed and solved by} his brother Jacob, Isaac Newton, the Marquis de~L'Hospital, Gottfried Leibniz and Tschirnhaus. A solution to the brachistochrone problem is presented in Example~\ref{ex:brachisto}. In 1697 Johann Bernoulli posed and solved the {\em geodesics\/} problem\index{geodesics}\index{geodesics problem!points on a surface} which can be formulated as follows. Find a line of shortest length that connects two given points on a given surface---see Figure~\ref{fig:geodesic} for an illustration of the problem. The curves of the least length connecting pairs of points on a surface are called {\em geodesics}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{geodesic.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/geodesic.bmp}
%\hskip1in\special{bmp:c:/mybook/geodesic.bmp x=4in y=3in}
\caption{The problem of finding geodesics.}
\label{fig:geodesic}
\end{figure}
A general procedure for solving\index{geodesics problem!general procedure for solving} the above type of problems, using the calculus of variations, was given by Leonhard Euler (who flourished between 1707 and 1783).

To acquaint ourselves with the methods of the calculus of variations, we need first to introduce concepts concerning functionals.

\subsection{Variation and Its Properties}
A {\em function\/} can be viewed as a rule,\index{function!a rule} or a mapping, that assigns to each element of some set a unique element of a possibly different set. In particular, a function $x:\mb{\Bb R} \to \mb{\Bb R}$ of a real variable $t$,  assigns to each real number a unique real number. An increment of the argument\index{increment of argument} of a function\index{function!increment of argument} of one variable $t$ is $\Delta t = t-t_1$.


Similarly, a {\em functional\/} is a mapping that assigns to each function, from some class of functions, a unique number.\index{functional!mapping on space of functions} Thus, we can say that a functional is a ``function of a function.''\index{functional!function of a function} Let $x:\mb{\Bb R} \to \mb{\Bb R}$ be an argument of a functional. By a {\em variation \/}\index{variation}\index{functional!variation} $\delta x$ of an argument $x$ of a functional $v$ we mean the difference of two functions
\[
\delta x = x-x_1.
\]
We assume that $x$ can change in an arbitrary way in some class of functions.

Recall, that a function $x:\mb{\Bb R} \to \mb{\Bb R}$ is continuous if a small change of its argument $t$ corresponds to a small change of the value $x(t)$ of the function. Similarly, a functional $v$ is said to be continuous\index{continuous functional}\index{functional!continuous} if a ``small'' change of its argument $x$ corresponds to a small change of the value of the functional.

A functional $v$ is called linear\index{linear functional}\index{functional!linear} if
\[
v\left(a x_1 + x_2 \right)=av\left(x_1 \right) + v \left( x_2 \right),
\]
where $a$ is a constant.

We will now introduce the notion of the {\em variation of a functional\/}.\index{variation of a functional}\index{functional!variation of} This notion is analogous to the notion of a function differential.\index{function differential}\index{variation of a functional!analog of function differential} To connect the two, we first consider the case of a function of one variable. Let $x$ be a differentiable function defined on an open interval $U$, and let $t\in U$. The derivative of $x$ at $t$ is defined as
\begin{equation}
x'(t)=\lim_{\Delta t \to 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}.
\label{eq:foovard}
\end{equation}
Let
\begin{equation}
\varphi(\Delta t)=\frac{x(t+\Delta t)-x(t)}{\Delta t}-x'(t).
\label{eq:defphifov}
\end{equation}
The function $\varphi$ is not defined at $\Delta t=0$, however
\[
\lim_{\Delta t \to 0}\varphi(\Delta t)=0.
\]
We rewrite~(\ref{eq:defphifov}) as
\begin{equation}
x(t+\Delta t)-x(t)=x'(t) \Delta t +\varphi(\Delta t) \Delta t.
\label{eq:rewusphi}
\end{equation}
The above relation has meaning only when $\Delta t \ne 0$. To make it hold at $\Delta t = 0$, we define $\left.\varphi(\Delta t)\right\arrowvert_{\Delta t=0}
 = 0$. To proceed, let
\begin{equation}
\begin{array}{lcc}
\beta(\Delta t)=\varphi(\Delta t) & \mb{if} & \Delta t >0\\
\beta(\Delta t)=-\varphi(\Delta t) & \mb{if} & \Delta t <0.
\end{array}
\end{equation}
Then, if $x$ is differentiable, there exists a function $\beta$ such that
\begin{eqnarray}
x(t+\Delta t)-x(t)=\Delta x&=&x'(t)\Delta t +\beta(\Delta t)|\Delta t|\nonumber\\
&=& L(t,\Delta t) +\beta(\Delta t)|\Delta t|,
\label{eq:fordeofx}
\end{eqnarray}
where $\lim_{\Delta t \to 0}\beta(\Delta t)=0$, and $L(t,\Delta t)=x'(t)\Delta t$ is a linear function in $\Delta t$.

For a real-valued function\index{real-valued function}\index{function!real-valued} $f=f(\b x):\mb{\Bb R}^n \to \mb{\Bb R}$, relation~(\ref{eq:fordeofx}) takes the form
\begin{eqnarray}
f(\b x+\Delta \b x)-f(\b x)=\Delta f&=&D f(\b x)\Delta \b x +\beta(\Delta \b x)\|\Delta \b x\|\nonumber\\
&=& L(\b x,\Delta \b x) +\beta(\Delta \b x)\|\Delta \b x\|,
\label{eq:fordeofmnx}
\end{eqnarray}
where $\lim_{\Delta \b x \to \b 0}\beta(\Delta \b x)=0$, and
\[
L(\b x,\Delta \b x)=D f(\b x)\Delta \b x=\nabla f(\b x)^{\top}\Delta \b x
\]
is a linear function\index{linear function}\index{function!linear} in $\Delta \b x$. Relation~(\ref{eq:fordeofmnx}) holds for operators defined over Banach spaces with scalars being either real or complex numbers---see Definition~\ref{defn:banachspa} of Banach space on page~\pageref{defn:banachspa}. The functionals we consider are operators over a Banach space\index{Banach space!of continuous functions}\index{functional!operator over a Banach space} of continuous functions, ${\cal C}([t_0,t_1])$, with the norm given by~(\ref{eq:rvfofone}) on page~\pageref{eq:rvfofone}.

If an increment $\Delta v = v \left( x+ \delta x \right) - v\left(x \right)$ of a functional $v$ can be represented as
\[
\Delta v = L\left( x, \delta x \right) + \beta \left(\delta x \right) \|\delta x \|,
\]
where $L\left( x, \delta x \right)$ is a linear functional with respect to $\delta x$, the term $\|\delta x\|=\max_{t\in[t_0,t_1]}|\delta x |$ denotes the maximal value of $| \delta x |$, and $\beta \left(\delta x \right) \to 0$ if $\|\delta x \| \to 0$, then the linear part $L$ of $\Delta v$ is called the {\em variation of the functional\/}\index{variation of the functional}\index{functional!variation of} and is denoted $\delta v$, that is,
\[
\delta v = L\left( x, \delta x \right).
\]
\begin{example}
\label{ex:variation}
\rm
We find the variation of the functional
\[
v=\int_{0}^{1}\left( 2 x^2(t) + x(t)\right)dt.
\]
For this we first calculate its increment to get
\begin{eqnarray*}
\Delta v &=& v\left( x + \delta x \right) - v \left( x \right)\\
&=& \int_{0}^{1}\left(2(x+\delta x )^2 + (x + \delta x) \right)dt - \int_{0}^{1}\left( 2x^2 + x \right)dt\\
&=& \int_{0}^{1} \left( 2x^2 + 4x \delta x + 2 (\delta x)^2 +x + \delta x -2x^2 - x \right)dt\\
&=& \int_{0}^{1} ( 4x +1 ) \delta x dt + 2 \int_{0}^{1} (\delta x)^2 dt.
\end{eqnarray*}
The linear part of $\Delta v$ is
\[
\delta v = \int_{0}^{1} ( 4x +1 ) \delta x dt,
\]
which is the variation of the given functional.
\end{example}
\markex

To proceed, note that the linear part on the right-hand side of~(\ref{eq:fordeofmnx}) can be computed as
\begin{equation}
D f(\b x)\Delta \b x=L(\b x,\Delta \b x)=\left.\frac{d}{d\alpha}f(\b x+\alpha \Delta \b x)\right\arrowvert_{\alpha=0}
\label{eq:dirdernon}
\end{equation}
As it was mentioned before, the above also holds for operators over Banach spaces other than Euclidean spaces\index{Euclidean space!Banach space}. Using this observation, we state and prove a lemma that allows us to calculate the variation of a given functional more efficiently than the method we used in Example~\ref{ex:variation}.
\begin{lemma}
\label{lemma:variation}
\[
\delta v = \left.\frac{d}{d \alpha}v\left( x + \alpha \delta x \right) \right\arrowvert _{\alpha = 0}.
\]
\end{lemma}
{\bf Proof}\quad Suppose that for a given functional $v$ there exists its variation. This means that we can represent $\Delta v$ as
\[
\Delta v = v\left(x + \alpha \delta x \right) - v \left(x\right)
= L ( x,\alpha \delta x) + \beta(\alpha \delta x)|\alpha|\,\|\delta x\|.
\]
Then, the derivative of $v\left(x + \alpha \delta x \right)$ with respect to $\alpha$ evaluated at $\alpha = 0$ is equal to
\begin{eqnarray*}
\lim _{\alpha \to 0} \frac{\Delta v}{\alpha} &=& \lim _{\alpha \to 0} \frac{L ( x,\alpha \delta x) + \beta(\alpha \delta x)|\alpha|\,\|\delta x\|}{\alpha}\\
&=& \lim _{\alpha \to  0} \frac{L(x,\alpha \delta x)}{\alpha}+\lim_{\alpha \to 0}\frac{\beta(\alpha \delta x)|\alpha|\,\|\delta x\|}{\alpha}\\
&=& L(x, \delta x),
\end{eqnarray*}
because $L(\cdot,\cdot)$ is linear with respect to the second argument and hence,
\[
L(x,\alpha \delta x)=\alpha L(x, \delta x).
\]
Furthermore,
\[
\lim_{\alpha \to 0}\frac{\beta(\alpha \delta x)|\alpha|\,\|\delta x\|}{\alpha} = \lim_{\alpha \to 0}\beta(\alpha \delta x)\|\delta x\|=0.
\]
This completes the proof.\\
\hspace*{5.5in} $\Box$

\begin{example}\rm
We use Lemma~\ref{lemma:variation} to find the variation of the functional considered in Example~\ref{ex:variation}. We have
\begin{eqnarray*}
\delta v &=& \frac{d}{d \alpha}v\left( x + \alpha \delta x \right) \arrowvert _{\alpha = 0}\\
&=& \left.\frac{d}{d \alpha} \left( \int_{0}^{1}\left(2(x+ \alpha \delta x )^2 + (x + \alpha \delta x) \right)dt\right) \right\arrowvert _{\alpha =0}\\
&=& \left.\int_0^1\left( 4(x+\alpha \delta x ) \delta x + \delta x \right) dt \right\arrowvert_{\alpha =0}\\
&=& \int_0^1 ( 4x + 1) \delta x dt,
\end{eqnarray*}
which is the variation of the given functional in Example~\ref{ex:variation}.
\end{example}\markex

\begin{example}\rm
We find the variation of the functional\index{functional!variation of}
\[
v(x_1,x_2)=v(x_1(t),x_2(t))=\int_0^1 x_1e^{-x_2}dt,
\]
using the formula
\[
\delta v=\left.\frac{d}{d \alpha}v(\b x+\alpha\delta \b x)\right\arrowvert_{\alpha=0}.
\]
We have
\begin{eqnarray*}
\delta v&=&\left.\frac{d}{d \alpha}\left( \int_0^1\left(x_1+\alpha\delta x_1 \right)e^{-\left(x_2+\alpha\delta x_2 \right)}dt\right)\right\arrowvert_{\alpha =0}\\
&=& \left.\int_0^1 \frac{d}{d \alpha}\left(\left(x_1+\alpha\delta x_1 \right)e^{-\left(x_2+\alpha\delta x_2 \right)}dt\right)\right\arrowvert_{\alpha =0}\\
&=& \left.\int_0^1 \left(\delta x_1 e^{-(x_2+\alpha \delta x_2)} - \delta x_2 e^{-(x_2+\alpha \delta x_2)}(x_1+\alpha \delta x_1)  \right)dt\right\arrowvert_{\alpha =0}\\
&=& \int_0^1 \left(e^{-x_2}\delta x_1 - x_1e^{-x_2} \delta x_2\right)dt.
\end{eqnarray*}
Thus,
\[
\delta v= \int_0^1 \left(e^{-x_2}\delta x_1 - x_1e^{-x_2} \delta x_2\right)dt .
\]
\end{example}\markex
\begin{defn}
A functional $v$ is maximal on some function\index{functional!maximal on a function} $x_0$ if the value of the functional on any function $x$ close to $x_0$ is less than or equal to that on $x_0$, that is,
\[
\Delta v = v\left( x \right) - v \left( x_0 \right) \le 0.
\]
\end{defn}\markdefn

We now state a first order necessary condition for a given functional to achieve a maximum or minimum\index{first order necessary condition!for a functional} on a curve.
\begin{thm}\label{thm:fonc}
If for a functional $v$ there exists a variation and the functional reaches a maximum or a minimum on a curve $x_0$, then
\[
\delta v\left( x_0 \right) = 0.
\]
\end{thm}
{\bf Proof}\quad Suppose that $x_0$ and $\delta x$ are fixed. Then
\[
v\left( x_0 + \alpha \delta x \right) = \Phi(\alpha ),
\]
that is, $v\left( x_0 + \alpha \delta x \right)$ is a function of $\alpha$ that for $\alpha = 0$ achieves its maximum or minimum. This implies that $\frac{d}{d \alpha} \Phi = 0$, or equivalently,
\[
\left.\frac{d}{d \alpha}v\left( x_0 + \alpha \delta x \right) \right\arrowvert _{\alpha = 0}=0.
\]
Hence, by Lemma~\ref{lemma:variation}, we have $\delta v = 0$. Therefore, a functional that achieves its maximum or a minimum on a function $x_0$, has its variation evaluated at $x_0$ equal to zero.\\
\hspace*{5.5in} $\Box$

The above results are used in the next section to determine conditions that a curve has to satisfy to be a maximizer or a minimizer, that is, to be an extremal\index{extremal!of a functional} of a given functional.

\subsection{Euler-Lagrange Equation}

In this section, we consider the following optimization problem.

{\bf Problem 1}\label{problem_1}\quad Suppose that we are given two points $(t_0,x_0)$ and $(t_1,x_1)$ in the $(t,x)$-plane. We wish to find a curve, trajectory, joining the given points such that the functional
\begin{equation}
v\left( x \right)=\int_{t_0}^{t_1}F\left(t,x,\dot{x}\right)dt
\label{eq:funcfeleq}
\end{equation}
along this trajectory can achieve its extremal, that is, maximal or minimal, value. The problem is illustrated in Figure~\ref{fig:prob_1}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{prob_1.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/prob_1.bmp}
%\hskip1in\special{bmp:c:/mybook/prob_1.bmp x=4in y=3in}
\caption{Illustration of Problem~1.}
\label{fig:prob_1}
\end{figure}
In solving Problem~1, we use the following lemma.
\begin{lemma}
\label{lemma:flcv}
Suppose that a function $\phi:\mb{\Bb R} \to \mb{\Bb R}$ is continuous on the interval $[t_0,\,t_1]$. Then,
\[
\int_{t_0}^{t_1}\phi(t) \delta x(t) dt =0
\]
if and only if $\phi(t)=0$ at every point of the interval $[t_0,\,t_1]$.
\end{lemma}
{\bf Proof}\quad Sufficiency is obvious because we can satisfy $\int_{t_0}^{t_1}\phi(t) \delta x(t) dt =0$ by choosing $\phi$ that vanishes at every point of the interval $[t_0,\,t_1]$.

We prove necessity by contraposition. We thus prove that if $\phi$ is nonzero on some subinterval of $[t_0,\, t_1]$, then $\int_{t_0}^{t_1}\phi(t) \delta x(t) dt \ne 0$. Specifically, assume, without loss of generality, that $\phi$ is positive on the interval $[a,\, b]$, where $t_0\le a \le b \le t_1$. Note that if $\phi$ were nonzero at only one point, then because it is continuous, it would have to be nonzero in some neighborhood\index{neighborhood!of a point} of that point. Hence, we can assume that $a<b$. Consider a variation given by
\begin{equation}
\delta x = \delta x(t) = \left\{\begin{array}{ccc}
k^2(t-a)^2(t-b)^2 & \mb{for}& a \le t \le b\\
0 & & \mb{elsewhere}.
\end{array}
\right.
\label{eq:vardef}
\end{equation}
If we assumed $\phi$ to be negative on the interval $[a,\, b]$, then we would use $\delta x$ to be negative of~(\ref{eq:vardef}).  We illustrate $\delta x$ in Figure~\ref{fig:flcv}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{flcv.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/flcv.bmp}
%\hskip1in\special{bmp:c:/mybook/flcv.bmp x=4in y=3in}
\caption{Illustration of the proof of Lemma~\ref{lemma:flcv}.}
\label{fig:flcv}
\end{figure}
Because $\delta x$ is positive in $[a,\,b]$ and zero elsewhere in the interval $[t_0,\, t_1]$, the integrand of $\int_{t_0}^{t_1}\phi(t) \delta x(t) dt$ will have the same sign as $\phi$ in $[a,\, b]$ and will be zero outside this subinterval. Thus
\[
\int_{t_0}^{t_1}\phi(t) \delta x(t) dt =\int_{a}^{b}\phi(t) \delta x(t) dt \ne 0,
\]
and the proof is complete.\\
\hspace*{5.5in} $\Box$

{\bf Solution to Problem 1}\label{soln_to_prob1}\quad We assume that there is an optimal trajectory $x=x(t)$ joining the points $(t_0,x_0)$ and $(t_1,x_1)$ such that $\delta v(x(t))=0$. Let us then consider an arbitrary acceptable curve $x^{*}=x^{*}(t)$ that is close to $x$ and from a family of curves\index{family of curves}
\[
x(t,\alpha)=x(t)+\alpha\left(x^{*}(t)-x(t) \right)=x(t)+\alpha \delta x(t).
\]
Note that for $\alpha = 0$ we get $x(t)$, and for $\alpha =1$ we obtain $x^{*}(t)$. Furthermore,
\[
\frac{d}{dt}(\delta x(t))=\frac{d}{dt}\left(x^{*}(t)-x(t) \right)= \delta x',
\]
and analogously
\[
\frac{d^2}{dt^2} (\delta x(t)) =\frac{d^2}{dt^2}\left(x^{*}(t)-x(t) \right)= \delta x''.
\]
We are now ready to investigate the behavior of the functional $v$ given by~(\ref{eq:funcfeleq}) on the curves from the family $x(t,\alpha)$. Note that the functional $v$ can be considered as a function of $\alpha$, that is, $v\left(x(t,\alpha)\right)=\Phi(\alpha)$. It follows from the first order necessary condition for a point to be an extremizer\index{extremizer!of a function} of a function of one variable that $\alpha =0$ is a candidate to be an extremizer of $\Phi$ if
\[
\left.\frac{d \Phi (\alpha)}{d \alpha}\right\arrowvert _{\alpha = 0}=0,
\]
where
\[
\Phi (\alpha ) = \int_{t_0}^{t_1}F\left(t, x(t, \alpha),x'(t, \alpha) \right)dt,
\]
and $x'(t,\alpha)=\frac{d}{dt}x(t,\alpha)$. We evaluate $\frac{d \Phi (\alpha)}{d \alpha}$ to obtain
\[
\frac{d \Phi(\alpha)}{d \alpha}=\int_{t_0}^{t_1} \left( F_{x}\frac{d}{d \alpha} x(t, \alpha) + F_{x'}\frac{d}{d \alpha}x'(t, \alpha) \right) dt,
\]
where $F_{x}=\frac{\partial}{\partial x} F \left(t, x(t,\alpha ),x'(t,\alpha )\right)$, and $F_{x'}=\frac{\partial}{\partial x'} F \left(t, x(t,\alpha ),x'(t,\alpha) \right)$. Because
\[
\frac{d}{d \alpha}x(t, \alpha) = \frac{d}{d \alpha}\left( x(t) + \alpha \delta x(t) \right)= \delta x(t),
\]
and
\[
\frac{d}{d \alpha}x'(t, \alpha) = \frac{d}{d \alpha}\left( x'(t) + \alpha \delta x'(t) \right)= \delta x'(t),
\]
we can write
\[
\frac{d}{d \alpha}\Phi (\alpha )= \int_{t_0}^{t_1} \left( F_{x}(t, x(t,\alpha), x'(t, \alpha)) \delta x(t) + F_{x'}(t, x(t,\alpha), x'(t, \alpha)) \delta x'(t) \right) dt.
\]
By Lemma~\ref{lemma:variation},
\[
\frac{d}{d \alpha}\Phi(0)=\delta v.
\]
Therefore,
\[
\delta v = \int_{t_0}^{t_1} \left( F_x \delta x + F_{x'} \delta x' \right) dt,
\]
where for convenience we dropped the argument $t$. Integrating by parts the last term, and taking into account that $\delta x' = (\delta x)'$ yields
\[
\delta v = \left.\left( F_{x'} \delta x \right) \right\arrowvert ^{t_1}_{t_0} + \int_{t_0}^{t_1} \left( F_x - \frac{d}{dt} F_{x'} \right) \delta x dt.
\]
The end points $(t_0, x_0)$ and $(t_1, x_1)$ are fixed. This implies that
\[
\delta x \arrowvert_{t=t_0}=x^{*}(t_0)-x(t_0)=0,
\]
 and
\[
\delta x \arrowvert_{t=t_1}=x^{*}(t_1)-x(t_1)=0.
\]
Hence
\[
\delta v =  \int_{t_0}^{t_1} \left( F_x - \frac{d}{dt} F_{x'} \right) \delta x dt.
\]
By Lemma~\ref{lemma:flcv},
\[
\delta v = 0 \quad\mb{if and only if}\quad F_x-\frac{d}{dt}F_{x'}=0.
\]

The equation,
\begin{equation}
\fbox{$\displaystyle F_x-\frac{d}{dt}F_{x'}=0 $}
\label{eq:eleqincalv}
\end{equation}
is known as the Euler-Lagrange equation\index{Euler-Lagrange equation}. The curves $x=x(t,C_1,C_2)$ that are solutions to the Euler-Lagrange equations are called {\em extremals\/}.\index{extremals!Euler-Lagrange equation solutions}\index{Euler-Lagrange equation solutions!extremals} The constants $C_1$ and $C_2$ are the integration constants. Thus, a functional can only achieve its extremum on the extremals.\index{functional!extremum only on extremals}

The Euler-Lagrange equation $F_x-\frac{d}{dt}F_{x'}=0$ can also be represented as
\[
F_x-F_{x' t}-F_{x' x}x'-F_{x'x'}x''=0.
\]
\index{Euler-Lagrange equation!alternative representation}

\begin{example}\rm
We find the extremal of the functional\index{functional!extremal}\index{extremal of!functional}
\begin{equation}
v=\int_{0}^{\frac{\pi}{2}}\left( (x')^2-x^2 \right)dt,
\label{eq:extreforv}
\end{equation}
where $x(0)=0$, and $x\left(\frac{\pi}{2} \right)=1$. The Euler-Lagrange equation for the problem takes the form
\begin{eqnarray*}
F_x-\frac{d}{dt}F_{x'} &=& \frac{\partial}{\partial x}\left((x')^2-x^2 \right) - \frac{d}{dt}\left( \frac{\partial}{\partial x'}\left((x')^2-x^2 \right)\right)\\
&=& -2x -\frac{d}{dt} \left(2 x'\right)\\
&=& -2x -2 x''\\
&=& 0.
\end{eqnarray*}
Equivalently, the Euler-Lagrange equation for the functional~(\ref{eq:extreforv}) can be represented as
\[
x''+x=0.
\]
The characteristic equation of the above differential equation is $s^2+1=0$. Therefore, the general solution of the Euler-Lagrange equation has the form
\[
x(t)=C_1 \cos t + C_2 \sin t.
\]
Taking into account the end-point conditions gives
\begin{eqnarray*}
x(0) &=& C_1 =0,\\
x\left( \frac{\pi}{2}\right) &=& C_2=1.
\end{eqnarray*}
Thus, we have one extremal
\[
x(t)=\sin t.
\]
The functional~(\ref{eq:extreforv}) can be minimized or maximized only on the extremal $x(t)= \sin t$.
\end{example}
\markex

\begin{example}\label{ex:fixedpt_fun}\rm
We find the extremal of the functional\index{functional!extremal}\index{extremal of!functional}
\begin{equation}
v=\int_0^1\left((x')^2 + 12 xt \right)dt,
\label{eq:nextvfore}
\end{equation}
where $x(0)=0$, and $x(1)=1$. The Euler-Lagrange equation for the this functional is
\begin{eqnarray*}
F_x-\frac{d}{dt}F_{x'} &=& \frac{\partial}{\partial x}\left((x')^2 + 12 xt \right)-\frac{d}{dt}\left(\frac{\partial}{\partial x'}\left((x')^2 +12 xt \right)\right)\\
&=& 12t -\frac{d}{dt}\left(2x'\right)\\
&=& 12t - 2x''\\
&=& 0.
\end{eqnarray*}
Equivalently, the Euler-Lagrange equation for the functional~(\ref{eq:nextvfore}) can be represented as
\[
x''=6t.
\]
The general solution to the above differential equation is
\[
x(t)=C_1 t +C_2 +t^3.
\]
From the end-point conditions, we get
\begin{eqnarray*}
x(0) &=& C_2=0,\\
x(1) &=& C_1 +1=1.
\end{eqnarray*}
Hence, $C_1=0$, and the extremal has the form
\[
x(t)=t^3.
\]
\end{example}
\markex

Suppose now that $A=A(t_0,x_{10},\ldots,x_{n0})$ and $B=B(t_1,x_{11},\ldots,x_{n1})$ are two given points in the $(n+1)$-dimensional space of $(n+1)$-tuples\index{(n+1)-tuples}, $(t,x_1,\ldots,x_n)$, along with the functional
\[
v=v(x_1,\ldots,x_n)=\int_{t_0}^{t_1}F(t,x_1,\ldots,x_n,\dot{x}_1,\ldots,\dot{x}_n)dt.
\]
To find a trajectory along which $\delta v=0$, we consider only the variation in one coordinate\index{variation in one coordinate}, say $x_i$, while other coordinates are fixed. We perform this process for each coordinate. The functional $v$ will depend only upon $x_i$, and we write
\[
v(x_1,\ldots,x_i,\ldots,x_n)=\tilde{v}(x_i).
\]
Thus, we reduced the given $n$-dimensional problem to $n$ one-dimensional problems. Proceeding as before, we conclude that the functional $v=v(x_1,\ldots,x_n)$ can only be extremized on the trajectories that satisfy $n$ Euler-Lagrange equations,
\begin{equation}
F_{x_i}-\frac{d}{dt}F_{\dot{x}_i}=0 ,\quad i=1,\ldots,n.
\label{eq:manyeleq}
\end{equation}

Observe that if in~(\ref{eq:manyeleq}) we set $F=L=K-U$, where $L$ is the Lagrangian,\index{Lagrangian} $K=K(\dot{q}_1,\ldots,\dot{q}_n)$ is the kinetic energy of the particle,\index{kinetic energy!of particle} $U=U(q_1,\ldots,q_n)$ is its potential energy,\index{potential energy!of particle} and $x_i=q_i$, then~(\ref{eq:manyeleq})  becomes
\[
\frac{\partial L}{\partial q_i}-\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i} \right)=0,\quad i=1,2,\ldots,n,
\]
which are the Lagrange equations of motion for a particle. Referring to the solution to Problem~1 on page~\pageref{soln_to_prob1}, we conclude that a particle moves along a path on which the functional
\[
J=\int_{t_1}^{t_2}L dt,
\]
is extremized, which agrees with the Hamilton principle\index{Hamilton principle} saying that a particle moves along the path that minimizes the integral $J$.

In our discussion up to now, both end-points were fixed. We would like to solve a more general optimization problem, where the end-points need not be fixed. For simplicity, we assume first that one of the end-points is fixed, while the other is free. After we solve this problem, we can easily generalize to the case when both end-points are free.

{\bf Problem 2}\label{problem_2}\quad In the $(t,x)$-plane the left-end point $A$ is fixed,\index{end point fixed} while the right end-point $B$ is free---see Figure~\ref{fig:fix_end}.\index{end point free}  It is required to find a trajectory $x=x(t)$ starting at $A$ on which the functional $v=\int_{t_0}^{t_1}F(t,x,x')dt$ can achieve its extremal value.
\begin{figure}
\centering
\includegraphics[scale=0.6]{fix_end.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/fix_end.bmp}
%\hskip1.2in\special{bmp:c:/mybook/fix_end.bmp x=4in y=3in}
\caption{Illustration of Problem 2.}
\label{fig:fix_end}
\end{figure}

{\bf Solution to Problem 2}\label{soln_to_prob2}\quad Following Elsgolc~\cite[pages~64--72]{elsgolc61}, we begin by assuming that a solution to Problem 2 exists, that is, there is a trajectory $x=x(t)$ joining $(t_0,x_0)$ and $(t_1,x_1)$ such that $\delta v=0$. This trajectory must satisfy the Euler-Lagrange equation.\index{Euler-Lagrange equation} Consider an adjoining trajectory specified by $x+\delta x$ joining $(t_0,x_0)$ and $(t_1+\delta t_1, x_1+\delta x_1)$. Then,
\begin{eqnarray*}
\Delta v &=& v(x+\delta x)-v(x)\\
&=& \int_{t_0}^{t_1 + \delta t_1} F(t,x+\delta x,x'+\delta x') dt - \int_{t_0}^{t_1}F(t,x,x')dt\\
&=& \int_{t_1}^{t_1+\delta t_1}F(t,x+\delta x,x'+\delta x') dt + \int_{t_0}^{t_1}\left( F(t,x+\delta x,x'+\delta x') - F(t,x,x') \right) dt.
\end{eqnarray*}
Our goal is to extract from $\Delta v$ its linear portion in $\delta x$, which we know is the variation of the given functional. We begin by approximating the first integral as
\[
\int_{t_1}^{t_1+\delta t_1}F(t,x+\delta x,x'+\delta x') dt \approx \left.F(t,x,x')\right\arrowvert_{t=t_1}\delta t_1.
\]
Expanding the integrand of the second integral into a Taylor series\index{Taylor series} and retaining only first-order terms yields
\[
\int_{t_0}^{t_1}\left( F(t,x+\delta x,x'+\delta x') - F(t,x,x') \right) dt=\int_{t_0}^{t_1}\left(F_x(t,x,x')\delta x + F_{x'}(t,x,x')\delta x' \right) dt.
\]
Integrating the second term by parts and combining with the first term gives
\[
\int_{t_0}^{t_1}\left(F_x\delta x + F_{x'}\delta x' \right) dt=\left(F_{x'}\delta x \right)\arrowvert_{t_0}^{t_1}+\int_{t_0}^{t_1}\left(F_x -\frac{d}{dt} F_{x'}\right)\delta x dt,
\]
where for simplicity of notation we suppressed the arguments of $F_x$ and $F_{x'}$. We assumed that a solution to the problem exists. Hence, an extremal must satisfy the Euler-Lagrange equation $F_x -\frac{d}{dt} F_{x'}=0$. Because the left end-point is fixed, we obtain
\[
\delta x \arrowvert_{t=t_0}=0.
\]
Taking into account the results of the above manipulations yields
\begin{equation}
\Delta v \approx \left.F(t,x,x')\right\arrowvert_{t=t_1}\delta t_1 + \left.\left(F_{x'}\delta x \right)\right\arrowvert_{t=t_1}.
\label{eq:deltavgetlin}
\end{equation}
The final step in obtaining $\delta v$ is to evaluate $\left(F_{x'}\delta x \right)\arrowvert_{t_1}$. The problem is two dimensional so we can illustrate it graphically as shown in Figure~\ref{fig:deltas}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{deltas.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/deltas.bmp}
%\hskip1.2in\special{bmp:c:/mybook/deltas.bmp x=4in y=3in}
\caption{A geometric interpretation of calculating the variation of the functional for Problem~2.}
\label{fig:deltas}
\end{figure}
Referring to Figure~\ref{fig:deltas}, we obtain
\begin{eqnarray*}
DB &=&\delta x \arrowvert_{t=t_1},\\
FC &=& \delta x_1,\\
EC &\approx & x'(t_1)\delta t_1,\\
DB &=& FC-EC.
\end{eqnarray*}
Hence,
\[
\delta x \arrowvert_{t=t_1} \approx \delta x_1 - x'(t_1) \delta t_1.
\]
Substituting the above into~(\ref{eq:deltavgetlin}) gives
\begin{eqnarray*}
\delta v &=& F(t,x,x')\arrowvert_{t=t_1}\delta t_1 + F_{x'} \arrowvert_{t=t_1}\left(\delta x_1 - x'(t_1) \delta t_1 \right)\\
&=& \left(F-x'F_{x'}\right) \arrowvert _{t=t_1}\delta t_1 + F_{x'} \arrowvert_{t=t_1} \delta x_1.
\end{eqnarray*}
If the variations $\delta t_1$ and $\delta x_1$ are arbitrary, we obtain $\delta v = 0$ if in addition to the Euler-Lagrange equation the following end-point conditions\index{end-point conditions} are also satisfied,
\[
\fbox{$\displaystyle \left(F-x'F_{x'}\right) \arrowvert _{t=t_1} = 0\quad\mb{and}\quad F_{x'} \arrowvert_{t=t_1} = 0 $}
\]
The above end-point relations are called the {\em transversality conditions}.\index{end-point conditions!transversality conditions}\index{transversality conditions!end-point conditions} In conclusion, the optimal trajectory must satisfy both the Euler-Lagrange equation and the transversality conditions.

\begin{example}\label{ex:freeedpt_xarb}\rm Consider the functional of Example~\ref{ex:fixedpt_fun}, where the left end-point is fixed as before, that is, $x(0)=0$, while the right end-point $x(1)$ is free. Here, we have $t_1=1$. Because $\delta t_1=0$ and $\delta x_1$ is arbitrary, the transversality conditions for the right end-point are reduced to one condition,
\[
F_{x'} \arrowvert_{t=t_1} = 0.
\]
Taking into account the solution to the Euler-Lagrange equation obtained in Example~\ref{ex:fixedpt_fun}, and then using the fact that the left end-point is fixed, we obtain
\begin{eqnarray*}
0&=&F_{x'} \arrowvert_{t=t_1}\\
 &=& 2 x'(t_1)\\
&=& 2(C_1+3t^2)\arrowvert_{t=1}\\
&=& 2(C_1+3).
\end{eqnarray*}
Hence, $C_1=-3$ and the extremal for this case is
\[
x(t)=-3t +t^3.
\]
In Figure~\ref{fig:x_arbit}, we show a plot of the above extremal.
\begin{figure}
\centering
\includegraphics[scale=0.45]{x_arbit.eps}
%\centering \epsfig{file=x_arbit.eps,height=2.75in}
\caption{A plot of the extremal, $x(t)=-3t+t^3$, of Example~\ref{ex:freeedpt_xarb}.}
\label{fig:x_arbit}
\end{figure}
\end{example}
\markex

\begin{example}\label{ex:freeedpt_trb}\rm
We consider again the same functional as in Example~\ref{ex:fixedpt_fun}, where, this time, the left end-point is fixed, as before, while at the right-end point $t_1$ is free and $x(t_1)=1$. Hence, in this case, $\delta x_1=0$ while $\delta t_1$ is arbitrary. The transversality conditions\index{transversality conditions} reduce to one condition,
\[
\left(F-x'F_{x'}\right) \arrowvert _{t=t_1} = 0.
\]
We perform simple calculations to obtain,
\[
F-x'F_{x'} = 12 x t -(x')^2,
\]
where
\[
x'=\left(C_1 t +t^3 \right)'=C_1+3t^2.
\]
Note that because $x(t_1)=C_1t_1+t_1^3=1$, we have
\[
C_1=\frac{1-t_1^3}{t_1}.
\]
Taking into account the above and the fact that $x(t_1)=1$, we perform simple algebraic manipulations to arrive at the following equation for $t_1$,
\[
\left(F-x'F_{x'}\right) \arrowvert _{t=t_1}=\left( 12 x(t_1) t_1-(x'(t_1))^2\right)=-4t_1^6+8t_1^3-1=0.
\]
Solving the above equation, we obtain the following two solutions,
\[
t_1=\left( 1\pm \sqrt{\frac{3}{4}}\right)^{1/3}.
\]
Plots of the two extremals for this example are shown in Figure~\ref{fig:t_arbit}.
\begin{figure}
\centering
\includegraphics[scale=0.45]{t_arbit.eps}
%\centering \epsfig{file=t_arbit.eps,height=2.75in}
\caption{Plots of the extremals of Example~\ref{ex:freeedpt_trb}.}
\label{fig:t_arbit}
\end{figure}
\end{example}
\markex


If the position of the end point $B$ is restricted to a curve $x=g_1(t)$, then
\[
\delta x_1 \approx g_1'(t_1) \delta t_1,
\]
and hence
\[
\delta v = \left( F +\left( g_1'-x'\right) F_{x'} \right)\delta t_1 =0
\]
if
\begin{equation}
\left( F +\left( g_1'-x'\right) F_{x'} \right) \arrowvert _{t=t_1}=0.
\label{eq:transcooncur}
\end{equation}
The above is the transversality condition\index{transversality condition} for the case when the end-point $B$ is restricted to the curve $x=g_1(t)$.

In the next section, we discuss a specific problem of minimizing a functional subject to constraints in the form of differential equations.\index{functional minimization!subject to constraints}

\section{Linear Quadratic Regulator}
\subsection{Algebraic Riccati Equation (ARE)}
Consider a linear system modeled by
\[
\dot{\b x}=\b A\b x+\b B\b u,\quad \b x(0)=\b x_0,
\]
and the associated performance index
\[
J=\int_0^{\infty}\left(\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)dt,
\]
where $\b Q=\b Q^{\top}\ge0$ and $\b R=\b R^{\top}>0$. Our goal is to construct a stabilizing linear state-feedback\index{stabilizing linear state-feedback}\index{linear state-feedback!stabilizing} controller of the form $\b u =-\b K\b x$ that minimizes the performance index $J$. We denote such a linear control law by $\b u^{*}$.

First, we assume that a linear state-feedback optimal controller exists such that the optimal closed-loop system
\[
\dot{\b x}=(\b A-\b B\b K)\b x
\]
is asymptotically stable. The above assumption implies that there is a Lyapunov function $V=\b x^{\top}\b P\b x$ for the closed-loop system, that is, for some positive definite matrix $\b P$ the time derivative $dV/dt$ evaluated on the trajectories of the closed-loop system is negative definite.
\begin{thm}
If the state-feedback controller $\b u^*=-\b K\b x$ is such that
\begin{equation}
\min_{\b u}\left(\frac{dV}{dt}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right) = 0,
\label{eq:mintoare}
\end{equation}
for some $V=\b x^{\top}\b P\b x$, then the controller is optimal.\index{controller!optimal}\index{optimal controller}
\end{thm}
{\bf Proof}\quad We can rewrite the condition of the theorem as
\[
\left.\frac{dV}{dt}\right\arrowvert_{\b u=\b u^*}+ \b x^{\top}\b Q\b x+\b u^{*\top}\b R\b u^*=0.
\]
Hence,
\[
\left.\frac{dV}{dt}\right\arrowvert_{\b u=\b u^*}=-\b x^{\top}\b Q\b x-\b u^{*\top}\b R\b u^*.
\]
Integrating both sides of the resulting equation with respect to time from $0$ to $\infty$, we obtain
\[
V(\b x(\infty))-V(\b x(0))=-\int_0^{\infty}\left(\b x^{\top}\b Q\b x+\b u^{*\top}\b R\b u^* \right)dt.
\]
Because by assumption the closed-loop system is asymptotically stable, $\b x(\infty)=\b 0$, and therefore
\[
V(\b x(0)) = \b x^{\top}_0\b P\b x_0= \int_0^{\infty}\left(\b x^{\top}\b Q\b x+\b u^{*\top}\b R\b u^* \right)dt.
\]
Thus, we showed that if a linear state-feedback controller satisfies the assumption of the theorem, then the value of the performance index for such a controller is
\[
J\left( \b u^* \right)=\b x^{\top}_0\b P\b x_0.
\]
To show that such a controller is indeed optimal, we use a proof by contradiction. We assume that~(\ref{eq:mintoare}) holds and that $\b u^*$ is not optimal. Suppose that $\tilde{\b u}$ yields a smaller value of $J$, that is,
\[
J\left(\tilde{\b u}\right) < J\left(\b u^*\right).
\]
It follows from~(\ref{eq:mintoare}) that
\[
\left.\frac{dV}{dt}\right\arrowvert_{\b u=\tilde{\b u}}+\b x^{\top}\b Q\b x+\tilde{\b u}^{\top}\b R\tilde{\b u}\ge 0,
\]
that is,
\[
\left.\frac{dV}{dt}\right\arrowvert_{\b u=\tilde{\b u}}\ge -\b x^{\top}\b Q\b x-\tilde{\b u}^{\top}\b R\tilde{\b u}.
\]
Integrating the above with respect to time from $0$ to $\infty$ yields
\[
V(\b x(0))\le \int_0^{\infty}\left( \b x^{\top}\b Q\b x+\tilde{\b u}^{\top}\b R\tilde{\b u}\right)dt
\]
implying that
\[
J(\b u^*)\le J(\tilde{\b u}),
\]
which is a contradiction, and the proof is complete.\\
\hspace*{5.5in} $\Box$

It follows from the above theorem that the synthesis of the optimal control law involves finding an appropriate Lyapunov function,\index{optimal control!Lyapunov function} or equivalently, the matrix $\b P$. The appropriate $\b P$ is found by minimizing
\begin{equation}
f(\b u)=\frac{dV}{dt}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u,
\label{eq:miniobfutfp}
\end{equation}
We first apply to~(\ref{eq:miniobfutfp}) the necessary condition for unconstrained minimization,
\[
\frac{\partial}{\partial \b u}\left.\left( \frac{dV}{dt}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)\right\arrowvert_{\b u=\b u^*}=\b 0.
\]
Differentiating the above yields
\begin{eqnarray*}
\frac{\partial}{\partial \b u}\left( \frac{dV}{dt}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)&=&\frac{\partial}{\partial \b u}\left( 2\b x^{\top}\b P\dot{\b x}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)\\
&=&\frac{\partial}{\partial \b u}\left( 2\b x^{\top}\b P\b A\b x +2\b x^{\top}\b P\b B\b u+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)\\
&=& 2\b x^{\top}\b P\b B+2\b u^{\top}\b R\\
&=&\b 0.
\end{eqnarray*}
Hence, the candidate for an optimal control\index{optimal control!candidate for} law has the form
\begin{equation}
\b u^*=-\b R^{-1}\b B^{\top}\b P\b x= -\b K\b x,
\label{eq:expforthegain}
\end{equation}
where $\b K=\b R^{-1}\b B^{\top}\b P$. Note that
\begin{eqnarray*}
\frac{\partial ^2}{\partial \b u^2}\left( \frac{dV}{dt}+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right) &=& \frac{\partial ^2}{\partial \b u^2}\left( 2\b x^{\top}\b P\b A\b x+2\b x^{\top}\b P\b B\b u+\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)\\
&=& \frac{\partial}{\partial \b u}\left( 2\b x^{\top}\b P\b B + 2\b u^{\top}\b R\right)\\
&=& 2\b R\\
&>&0.
\end{eqnarray*}
Thus the second order sufficiency condition for $\b u^*$ to minimize~(\ref{eq:miniobfutfp}) is satisfied.

We next turn our attention to finding an appropriate matrix $\b P$. The optimal closed-loop system\index{optimal!closed-loop system} has the form
\[
\dot{\b x}=\left(\b A-\b B\b R^{-1}\b B^{\top}\b P \right)\b x,\quad \b x(0)=\b x_0.
\]
Our optimal controller satisfies the equation
\[
\left.\frac{dV}{dt}\right\arrowvert_{\b u=\b u^*}+\b x^{\top}\b Q\b x+\b u^{*\top}\b R\b u^*=0,
\]
that is,
\[
2\b x^{\top}\b P\b A\b x+2\b x^{\top}\b P\b B\b u^*+\b x^{\top}\b Q\b x+\b u^{*\top}\b R\b u^*=0.
\]
We substitute the expression for $\b u^*$, given by~(\ref{eq:expforthegain}), into the above equation and represent it as
\[
\b x^{\top}\left( \b A^{\top}\b P+\b P\b A\right)\b x-2\b x^{\top}\b P\b B\b R^{-1}\b B^{\top}\b P\b x + \b x^{\top}\b Q\b x + \b x^{\top}\b P\b B\b R^{-1}\b B^{\top}\b P\b x=0.
\]
Factoring out $\b x$ yields
\[
\b x^{\top}\left( \b A^{\top}\b P+\b P\b A+\b Q-\b P\b B\b R^{-1}\b B^{\top}\b P\right)\b x=0.
\]
The above equation should hold for any $\b x$. For this to be true we have to have
\begin{equation}
\fbox{$\displaystyle \b A^{\top}\b P+\b P\b A+\b Q-\b P\b B\b R^{-1}\b B^{\top}\b P=\b O $}
\label{eq:riccatiale}
\end{equation}
The above equation is referred to as the {\em algebraic Riccati equation\/} (ARE).\index{algebraic Riccati equation (ARE)}\index{ARE|see{algebraic Riccati equation (ARE)}} In conclusion, the synthesis of the optimal linear state feedback controller minimizing the performance index
\[
J=\int_0^{\infty}\left(\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)dt
\]
subject to
\[
\dot{\b x}=\b A\b x+\b B\b u,\quad \b x(0)=\b x_0
\]
requires solving the ARE given by~(\ref{eq:riccatiale}).
\begin{example}\rm
Consider the following model of a dynamical system:
\[
\dot{x}=2u_1+2u_2,\qquad x(0)=3,
\]
and the associated performance index
\[
J=\int_{0}^{\infty}\left(x^2+ru_1^2+ru_2^2\right)dt,
\]
where $r>0$ is a parameter.
\begin{itemize}
\item[(i)] We first find the solution to the ARE corresponding to the linear state-feedback optimal controller. We have
\[
\b A=0, \quad \b B=\left[\begin{array}{cc}
2 & 2
\end{array}\right], \quad \b Q=1, \quad \b R= r \b I_2.
\]
The ARE for this problem is
\[
\b O = \b A^{\top} \b P + \b P \b A +\b Q -\b P \b B \b R^{-1} \b B^{\top} \b P= 1-\frac{8}{r}p^2,
\]
whose solution is
\[
p=\sqrt{\frac{r}{8}}.
\]
\item[(ii)] We write the equation of the closed-loop system driven by the optimal controller,
\[
\b u =-\b R^{-1} \b B^{\top} \b P \b x= -\frac{1}{\sqrt{2r}}\left[\begin{array}{c}
1\\
1
\end{array}\right]x.
\]
The closed-loop optimal system is described by
\[
\dot{x}= \left[\begin{array}{cc}
2 & 2
\end{array}\right]\b u= -\frac{4}{\sqrt{2r}}x.
\]
\item[(iii)] Finally, we find the value of $J$ for the optimal closed-loop system.\index{optimal!closed-loop system} We have
\[
J= \b x(0)^{\top}\b P\b x(0)= \frac{9}{2}\sqrt{\frac{r}{2}}.
\]
We can directly verify that $J=\frac{9}{2}\sqrt{\frac{r}{2}}$ as follows:
\begin{eqnarray*}
\min J &=& \min \int_{0}^{\infty}\left(x^2+ru_1^2+ru_2^2\right)dt\\
&=& \int_{0}^{\infty}9 \exp \left(-\frac{8t}{\sqrt{2r}}\right)\left(1+1 \right)dt= 9\frac{\sqrt{2r}}{4}= \frac{9}{2}\sqrt{\frac{r}{2}}.
\end{eqnarray*}
\end{itemize}
\end{example}
\markex
A major difficulty when solving the ARE is that it is a nonlinear equation. In the next section, we present an efficient method for solving the ARE\index{ARE!solving}\index{solving ARE}.

\subsection{Solving the ARE Using the Eigenvector Method}
We now present a method for solving the ARE referred to as the MacFarlane and Potter method, or the eigenvector method.\index{solving ARE!eigenvector method}\index{solving ARE!MacFarlane and Potter method}\index{eigenvector method for solving ARE}\index{MacFarlane and Potter method|see{eigenvector method for solving ARE}} In our discussion of the method, we follow the development of Power and Simpson~\cite{power78}. We begin by representing the ARE in the form
\begin{equation}
\left[\begin{array}{cc}
\b P & -\b I_n
\end{array}\right]\left[\begin{array}{cc}
\b A & \b -\b B\b R^{-1}\b B^{\top}\\
-\b Q & -\b A^{\top}
\end{array}\right]\left[\begin{array}{c}
\b I_n\\
\b P
\end{array}\right]=\b O.
\label{eq:areinmatform}
\end{equation}
The $2n\times 2n$ matrix in the middle is called the {\em Hamiltonian matrix\/}.\index{Hamiltonian matrix}\index{solving ARE!Hamiltonian matrix} We use the symbol $\b H$ to denote the Hamiltonian matrix, that is,
\[
\b H=\left[\begin{array}{cc}
\b A & \b -\b B\b R^{-1}\b B^{\top}\\
-\b Q & -\b A^{\top}
\end{array}\right].
\]
Thus, the ARE can be represented as
\[
\left[\begin{array}{cc}
\b P & -\b I_n
\end{array}\right]\b H\left[\begin{array}{c}
\b I_n\\
\b P
\end{array}\right]=\b O.
\]
We premultiply  the above equation by $\b X^{-1}$, and then postmultiply it by $\b X$, where $\b X$ is a nonsingular $n \times n$ matrix,
\begin{equation}
\left[\begin{array}{cc}
\b X^{-1}\b P & -\b X^{-1}
\end{array}\right]\b H\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]=\b O.
\label{eq:areafpreandp}
\end{equation}
Observe that if we can find matrices $\b X$ and $\b P\b X$ such that
\[
\b H\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]=\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]\b \Lambda,
\]
then equation~(\ref{eq:areafpreandp}) becomes
\[
\left[\begin{array}{cc}
\b X^{-1}\b P & -\b X^{-1}
\end{array}\right]\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]\b \Lambda=\b O.
\]
We have thus reduced the problem of solving the ARE to that of constructing appropriate matrices $\b X$ and $\b P\b X$. To proceed further, let $\b v_i$ be an eigenvector of $\b H$ and $s_i$ the corresponding eigenvalue, then
\[
\b H\b v_i=s_i\b v_i.
\]
In our further discussion, we assume that $\b H$ has at least $n$ distinct real eigenvalues among its $2n$ eigenvalues. (The results obtained can be generalized for the case when the eigenvalues of $\b H$ are complex or non-distinct.) Thus, we can write
\[
\b H\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]=\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]\left[\begin{array}{cccc}
s_1 & 0 & \cdots & 0\\
0 & s_2 & \cdots & 0\\
\vdots & \vdots & {} & \vdots\\
0 & 0 & \cdots & s_n
\end{array}\right].
\]
Let
\[
\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]=\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right],
\]
and
\[
\b \Lambda=\left[\begin{array}{cccc}
s_1 & 0 & \cdots & 0\\
0 & s_2 & \cdots & 0\\
\vdots & \vdots & {} & \vdots\\
0 & 0 & \cdots & s_n
\end{array}\right].
\]
The above choice of $\b X$ and $\b P\b X$ constitutes a possible solution to the equation
\[
\left[\begin{array}{cc}
\b X^{-1}\b P & -\b X^{-1}
\end{array}\right]\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]\b \Lambda=\b O.
\]
To construct $\b P$, we partition the $2n \times n$ eigenvector matrix $\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]$ into two $n \times n$ submatrices as follows
\[
\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]=\left[\begin{array}{c}
\b W\\
\b Z
\end{array}\right].
\]
Thus,
\[
\left[\begin{array}{c}
\b X\\
\b P\b X
\end{array}\right]=\left[\begin{array}{c}
\b W\\
\b Z
\end{array}\right].
\]
Taking $\b X=\b W$, $\b P\b X=\b Z$, and assuming that $\b W$ is invertible we obtain
\begin{equation}
\fbox{$\displaystyle \b P=\b Z\b W^{-1} $}
\label{eq:pzwinv}
\end{equation}
We now have to decide what set of $n$ eigenvalues should be chosen from those of $\b H$ to construct the particular $\b P$. In the case when all $2n$ eigenvalues of $\b H$ are distinct, the number of different matrices $\b P$ generated by the above described method is
\[
\frac{(2n)!}{(n!)^2}.
\]
Let $\b Q=\b C^{\top}\b C$ be a full rank factorization of $\b Q$, as defined on page~\pageref{eq:furafaofq}. It follows from Theorem~3 of Ku\v{c}era~\cite[page~346]{kucera72} that the Hamiltonian matrix $\b H$ has $n$ eigenvalues in the open left-half complex plane and $n$ in the open right-half plane if and only if the pair $(\b A,\b B)$ is stabilizable\index{stabilizable pair} and the pair $(\b A,\b C)$ is detectable\index{detectable pair}.\index{solving ARE!Ku\v{c}era conditions}\index{Ku\v{c}era conditions!for solving ARE} The matrix $\b P$ that we seek corresponds to the asymptotically stable eigenvalues of $\b H$. The eigenvalues of the Hamiltonian matrix come in pairs $\pm s_i$ (see Exercise~\ref{ex:mupomuneinhm}). The characteristic polynomial of $\b H$ contains only even powers of $s$. With such constructed matrix $\b P$, we have the following result:
\begin{thm}
The poles of the closed-loop system
\[
\dot{\b x}(t)=\left( \b A-\b B\b R^{-1}\b B^{\top}\b P\right)\b x(t)
\]
are those eigenvalues of $\b H$ having negative real parts.
\end{thm}
{\bf Proof}\quad Because $\left[\begin{array}{cccc}
\b v_1 & \b v_2 & \cdots & \b v_n
\end{array}\right]=\left[\begin{array}{c}
\b W\\
\b Z
\end{array}\right]$, we can write
\[
\left[\begin{array}{cc}
\b A & \b -\b B\b R^{-1}\b B^{\top}\\
-\b Q & -\b A^{\top}
\end{array}\right]\left[\begin{array}{c}
\b W\\
\b Z
\end{array}\right]=\left[\begin{array}{c}
\b W\\
\b Z
\end{array}\right]\b \Lambda.
\]
Performing appropriate multiplication of $n\times n$ block matrices yields
\[
\b A\b W-\b B\b R^{-1}\b B^{\top}\b Z=\b W \b \Lambda,
\]
or
\[
\b A-\b B\b R^{-1}\b B^{\top}\b Z\b W^{-1} =\b A-\b B\b R^{-1}\b B^{\top}\b P= \b W\b \Lambda \b W^{-1},
\]
because $\b P=\b Z\b W^{-1}$. Thus, the matrix $\b A-\b B\b R^{-1}\b B^{\top}\b P$ is similar to the matrix $\b \Lambda$ whose eigenvalues are the asymptotically stable eigenvalues of $\b H$. The proof is complete.\\
\hspace*{5.5in}$\Box$

\begin{example}\rm
Consider the following model of a dynamical system
\[
\dot{x} =2x + u,
\]
and the associated performance index
\[
J=\int_{0}^{\infty }(x^2 + ru^2) \, dt.
\]
Find the value of $r$ such that the optimal closed-loop system has its pole at $-3$.

We form the associated  Hamiltonian matrix
\[
{\b H}=\left[ \begin{array}{cc}
{\b A} & -\b B\b R^{-1}\b B^{\top} \\
-{\b Q} & -\b A^{\top}
\end{array}
\right] = \left[ \begin{array}{cc}
2 & -\frac{1}{r} \\
-1 & -2
\end{array}
\right].
\]
The characteristic equation of ${\b H}$ is
\[
\det(s \b I_2 - \b H)=s^2 - 4-\frac{1}{r} = 0.
\]
Hence,
\[
r=\frac{1}{5}
\]
results in the optimal closed-loop system having its pole located at $-3$.
\end{example}
\markex

We mention again that the eigenvector method for solving the ARE can be generalized to the case when the eigenvalues of the Hamiltonian matrix $\b H$ are not necessarily distinct or real.\index{eigenvector method!can be generalized}

\subsection{Optimal Systems With Prescribed Poles}
\label{subsec:optprep}
We can always find a state-feedback $u=-\b K\b x$ such that the closed-loop system $\dot{\b x}=(\b A-\b B\b K)\b x$ will have its poles in prescribed locations (symmetric with respect to the real axis) if and only if the pair $(\b A,\b B)$ is reachable. In the multi-input systems the gain matrix $\b K$ that shifts the poles to the prescribed locations is not unique. We would like to use the remaining degrees of freedom to construct a gain matrix $\b K$ that shifts the system poles to the desired locations and at the same time minimizes a performance index. Specifically, we are given a system model,
\[
\dot{\b x}(t)=\b A\b x(t)+\b B\b u(t),
\]
the desired set of closed-loop poles
\[
\{s_1,s_2,\ldots,s_n\},
\]
and a $m\times m$ positive definite matrix $\b R=\b R^{\top}$. We wish to construct a state-feedback control law $\b u(t)=-\b K\b x(t)$ and an $n \times n$ positive semi-definite matrix $\b Q=\b Q^{\top}$ such that the closed-loop system has its poles in the desired locations, and at the same time the performance index\index{poles!in desired locations}
\[
J=\int_0^{\infty}\left(\b x(t)^{\top}\b Q\b x(t)+\b u(t)^{\top}\b R\b u(t)\right) dt
\]
is minimized. In other words, we are interested in designing an optimal control system with prescribed poles.\index{optimal control system!with prescribed poles} We illustrate the problem with a very simple example.
\begin{example}\rm
For a dynamical system modeled by
\begin{eqnarray*}
\dot{x}&=&2x+u\\
&=&\b A\b x+\b B\b u,
\end{eqnarray*}
and the performance index
\begin{eqnarray*}
J&=&\int_0^{\infty}(qx^2+u^2)dt\\
&=& \int_0^{\infty}\left(\b x^{\top}\b Q\b x+\b u^{\top}\b R\b u\right)dt,
\end{eqnarray*}
find $q$ so that the eigenvalue of
\[
\b A-\b B\b R^{-1}\b B^{\top}\b P
\]
is equal to $-3$.

The corresponding Hamiltonian matrix is
\[
\b H=\left[\begin{array}{cc}
2 & -1\\
-q & -2
\end{array}\right].
\]
The characteristic polynomial of $\b H$ is
\[
\det [s\b I_2-\b H] = \det\left[\begin{array}{cc}
s-2 & 1\\
q & s+2
\end{array}\right]= s^2-4-q.
\]
We want
\[
\det [s\b I_2-\b H]=(s+3)(s-3)= s^2-9.
\]
Hence, the desired value of $q$ is $q=5$.
\end{example}
\markex

We now present a method for designing optimal control systems with prescribed poles\index{optimal control systems!with prescribed poles}. In constructing the gain matrix $\b K$ and the weight matrix $\b Q$, we use repeatedly a state variable transformation. To proceed, we analyze the system model, the performance index, and the Hamiltonian matrix\index{Hamiltonian matrix!in a new basis} in a new basis. Let
\[
\b x=\b T\tilde{\b x},
\]
where $\det \b T \ne 0$. Then, the system model in the new coordinates has the form
\begin{eqnarray*}
\dot{\tilde{\b x}}(t) &=& \b T^{-1}\b A\b T \tilde{\b x}(t)+\b T^{-1}\b B \b u(t)\\
&=& \tilde{\b A}\tilde{\b x}(t)+\tilde{\b B}\b u(t).
\end{eqnarray*}
The performance index $J$, after state variable transformation, takes the form
\[
\tilde{J}=\int_0^{\infty}\left(\tilde{\b x}(t)^{\top}\b T^{\top}\b Q\b T\tilde{\b x}(t)+\b u(t)^{\top}\b R\b u(t) \right) dt.
\]
The Hamiltonian matrix is transformed to
\begin{eqnarray*}
\tilde{\b H} &=& \left[\begin{array}{cc}
\b T^{-1}\b A\b T & -\b T^{-1}\b B\b R^{-1}\b B^{\top}\b T^{-\top}\\
-\b T^{T}\b Q\b T & -\b T^{\top}\b A^{\top}\b T^{-\top}
\end{array}\right]\\
&=& \left[\begin{array}{cc}
\tilde{\b A} & -\tilde{\b B}\b R^{-1}\tilde{\b B}^{\top}\\
-\tilde{\b Q} & -\tilde{\b A}^{\top}
\end{array}\right],
\end{eqnarray*}
where $\b T^{-\top}=\left(\b T^{-1}\right)^{\top}=\left(\b T^{\top}\right)^{-1}$. Observe that
\[
\tilde{\b H}=\left[\begin{array}{cc}
\b T^{-1} & \b 0\\
\b 0 & \b T^{\top}
\end{array}\right] \b H \left[\begin{array}{cc}
\b T & \b 0\\
\b 0 & \b T^{-\top}
\end{array}\right],
\]
that is, $\tilde{\b H}$ and $\b H$ are similarity matrices because
\[
\left[\begin{array}{cc}
\b T^{-1} & \b 0\\
\b 0 & \b T^{\top}
\end{array}\right]\left[\begin{array}{cc}
\b T & \b 0\\
\b 0 & \b T^{-\top}
\end{array}\right]=\b I_{2n}.
\]
Hence, the eigenvalues of $\b H$ and $\tilde{\b H}$ are identical. Our next step is to compute the characteristic equation of $\tilde{\b H}$ as a function of $\tilde{\b Q}$. Computing $\det\left(s\b I_{2n}-\tilde{\b H} \right)$ is facilitated after we perform some manipulations on the characteristic matrix $\left[s\b I_{2n}-\tilde{\b H}\right]$ that do not change its determinant. Specifically, we premultiply the characteristic matrix $\left[s\b I_{2n}-\tilde{\b H}\right]$ by a $2n\times 2n$ matrix whose determinant equals $1$. Then, the determinant of the product is the same as $\det\left(s\b I_{2n}-\tilde{\b H} \right)$, but it is much easier to evaluate. We have
\begin{eqnarray*}
&{}&\left[\begin{array}{cc}
\b I_n & \b 0\\
-\tilde{\b Q}\left(s\b I_n-\tilde{\b A}\right)^{-1} & \b I_n
\end{array}\right]\left[\begin{array}{cc}
s\b I_n-\tilde{\b A} & \tilde{\b B}\b R^{-1}\tilde{\b B}^{\top}\\
\tilde{\b Q} & s\b I_n+\tilde{\b A}^{\top}
\end{array}\right]\\
&{}&\qquad= \left[\begin{array}{cc}
s\b I_n-\tilde{\b A} & \tilde{\b B}\b R^{-1}\tilde{\b B}^{\top}\\
\b 0& s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}\left(s\b I_n-\tilde{\b A}\right)^{-1}\tilde{\b B}\b R^{-1}\tilde{\b B}^{\top}
\end{array}\right].
\end{eqnarray*}
Let
\[
\b F=\tilde{\b B}\b R^{-1}\tilde{\b B}^{\top},
\]
then
\begin{eqnarray*}
\det \left(s\b I_{2n}-\tilde{\b H}\right)&=&\det\left[\begin{array}{cc}
s\b I_n-\tilde{\b A} & \b F\\
\b 0& s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}\left(s\b I_n-\tilde{\b A}\right)^{-1}\b F
\end{array}\right]\\
&=& \det \left(s\b I_n-\tilde{\b A} \right) \det\left( s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}\left(s\b I_n-\tilde{\b A}\right)^{-1}\b F\right).
\end{eqnarray*}
Suppose that the matrix $\b A$ has $n$ linearly independent eigenvectors\index{eigenvectors!linearly independent}, hence it is diagonalizable. Let
\[
\b T=\left[\begin{array}{cccc}
\b t_1 & \b t_2 & \cdots & \b t_n
\end{array}\right]
\]
be a matrix composed of the eigenvectors of $\b A$. Therefore,
\[
\tilde{\b A} = \b T^{-1}\b A\b T= \left[\begin{array}{cccc}
\lambda_1 & 0 & \cdots & 0\\
0 & \lambda_2 & \cdots & 0\\
\vdots & \vdots & {} & \vdots\\
0 & 0 & \cdots & \lambda_n
\end{array}\right].
\]
Note that $\tilde{\b A}^{\top}=\tilde{\b A}$. Assume that
\[
\tilde{\b Q}=\tilde{\b Q}_j=\left[\begin{array}{ccccc}
0 & 0 & 0 &\cdots & 0\\
\vdots & \vdots&\vdots & {} & \vdots\\
0 & 0 & 0 &\cdots & 0\\
0 & 0 & \tilde{q}_{jj} & \cdots & 0\\
\vdots & \vdots &\vdots & {} & \vdots\\
0 & 0 & 0 & \cdots & 0
\end{array}\right],
\]
that is, $\tilde{\b Q}_j$ has only one non-zero element $\tilde{q}_{jj}$. Under the above assumptions the matrix $\left[ s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}_j\left(s\b I_n-\tilde{\b A}\right)^{-1}\b F\right]$ takes the form
\begin{eqnarray*}
&{}& s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}_j\left(s\b I_n-\tilde{\b A}\right)^{-1}\b F= \left[\begin{array}{cccc}
s+\lambda_1 & 0 & \cdots & 0\\
0 & s+\lambda_2 & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & s+\lambda_n
\end{array}\right]\\
&{}&\qquad -\left[\begin{array}{ccccc}
0 & 0 & 0&\cdots & 0\\
\vdots & \vdots & \vdots&{} & \vdots\\
0 & 0 & \tilde{q}_{jj} &0& 0\\
0 & 0 & 0&\cdots & 0
\end{array}\right]\left[\begin{array}{cccc}
\frac{1}{s-\lambda_1} & 0 & \cdots & 0\\
0 & \frac{1}{s-\lambda_2} & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & \frac{1}{s-\lambda_n}
\end{array}\right]\left[\begin{array}{cccc}
f_{11} & f_{12} & \cdots & f_{1n}\\
f_{21} & f_{22} & \cdots & f_{2n}\\
\vdots & \vdots & \ddots & \vdots\\
f_{n1} & f_{n2} & \cdots & f_{nn}
\end{array}\right]\\
&=&\left[\begin{array}{ccccc}
s+\lambda_1 & 0 &0& \cdots & 0\\
0 & s+\lambda_2 &0& \cdots & 0\\
\vdots & \vdots & \vdots &{}& \vdots\\
0 & 0 & s+\lambda_j&{} &0\\
\vdots & \vdots & \vdots &\ddots& \vdots\\
0 & 0 & 0&\cdots & s+\lambda_n
\end{array}\right]-\left[\begin{array}{ccccc}
0 & 0 &0& \cdots & 0\\
0 & 0 &0& \cdots & 0\\
\vdots & \vdots&\vdots & {} & \vdots\\
\frac{\tilde{q}_{jj}f_{j1}}{s-\lambda_j}&\frac{\tilde{q}_{jj}f_{j2}}{s-\lambda_j}&\frac{\tilde{q}_{jj}f_{jj}}{s-\lambda_j}& \cdots & \frac{\tilde{q}_{jj}f_{jn}}{s-\lambda_j}\\
\vdots & \vdots&\vdots & {} & \vdots\\
0 & 0 & 0&\cdots & 0
\end{array}\right]\\
&=& \left[\begin{array}{cccccc}
s+\lambda_1 & 0 &0 & \cdots & 0& 0\\
0 & s+\lambda_2 &0 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots& {}&\vdots&\vdots\\
-\frac{\tilde{q}_{jj}f_{j1}}{s-\lambda_j}&-\frac{\tilde{q}_{jj}f_{j2}}{s-\lambda_j}& \cdots & \cdots&s+\lambda_j-\frac{\tilde{q}_{jj}f_{jj}}{s-\lambda_j} & -\frac{\tilde{q}_{jj}f_{jn}}{s-\lambda_j}\\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\
0 & 0 & 0 & \cdots & 0 & s+\lambda_n
\end{array}\right].
\end{eqnarray*}
Hence,
\[
\det \left( s\b I_n+\tilde{\b A}^{\top}-\tilde{\b Q}_j\left(s\b I_n-\tilde{\b A}\right)^{-1}\b F\right)=\left(\prod^n_{\stackrel{i=1}{i\ne j}}(s+\lambda_i)\right)\left(s+\lambda_j-\frac{\tilde{q}_{jj}f_{jj}}{s-\lambda_j} \right).
\]
We can now evaluate the determinant of the Hamiltonian matrix\index{Hamiltonian matrix!determinant of},
\begin{eqnarray*}
\det(s\b I_{2n}-\b H) &=& \det(s\b I_{2n}-\tilde{\b H})\\
&=&\left(\prod^n_{i=1}(s-\lambda_i) \right)\left(\prod^n_{\stackrel{i=1}{i\ne j}}(s+\lambda_i)\right)\left(s+\lambda_j-\frac{\tilde{q}_{jj}f_{jj}}{s-\lambda_j} \right)\\
&=&\prod^n_{\stackrel{i=1}{i\ne j}}(s+\lambda_i)(s-\lambda_i)\left((s+\lambda_j)(s-\lambda_j)-\tilde{q}_{jj}f_{jj} \right).
\end{eqnarray*}
From the form of the determinant of the Hamiltonian matrix, we conclude that selecting an appropriate value of $\tilde{q}_{jj}$ allows us to shift the $j$-th pole of the closed-loop system into a prespecified location while leaving the other poles in their original positions. To see this we represent the term that involves $\tilde{q}_{jj}$ as
\begin{eqnarray*}
(s+\lambda_j)(s-\lambda_j)-\tilde{q}_{jj}f_{jj} &=& s^2-\lambda^2_j-\tilde{q}_{jj}f_{jj}\\
&=& s^2-\left( \lambda^2_j+\tilde{q}_{jj}f_{jj}\right)\\
&=& (s-s_j)(s+s_j),
\end{eqnarray*}
where
\[
s_j=\pm \sqrt{\lambda_j^2+\tilde{q}_{jj}f_{jj}}.
\]
From the above expression for $s_j$, we calculate $\tilde{q}_{jj}$ that results in the $j$-th pole being shifted to $-\sqrt{\lambda_j^2+\tilde{q}_{jj}f_{jj}}$. We obtain
\begin{equation}
\fbox{$\displaystyle \tilde{q}_{jj}=\frac{s^2_j-\lambda^2_j}{f_{jj}} $}
\label{eq:exprefoqjjt}
\end{equation}
Note that while selecting the desired location for $s_j$, we have to make sure that $\tilde{q}_{jj}>0$ to ensure positive semidefiniteness of the weight matrix\index{weight matrix} $\tilde{\b Q}_j$. Having $\tilde{\b Q}_j$, we compute $\b Q_j=\b T^{-\top}\tilde{\b Q}_j\b T^{-1}$ and solve the Riccati equation\index{Riccati equation}
\begin{equation}
\b A^{\top}\b P_j+\b P_j\b A+\b Q_j-\b P_j\b B\b R^{-1}\b B^{\top}\b P_j=\b O
\label{eq:onericcate}
\end{equation}
for $\b P_j$. We compute $\b u=-\b R^{-1}\b B^{\top}\b P_j\b x+\b u_j=-\b K_j\b x+\b u_j$ that shifts only one pole of the closed-loop system
\begin{eqnarray*}
\dot{\b x}&=&(\b A-\b B\b K_j)\b x+\b B\b u_j\\
&=& \b A_j\b x+\b B\b u_j
\end{eqnarray*}
to its specified position. We then start with the new system
\[
\dot{\b x}=\b A_j\b x+\b B\b u_j
\]
and shift the next pole to a specified location, and then combine the obtained feedback gains. The above method is based on the result, established by Solheim~\cite{solheim72}, that different feedback matrices $\b K_j$ and the corresponding weight matrices $\b Q_j$ can be combined to give the final optimal gain\index{optimal gain} $\b K$ and the corresponding weight matrix\index{weight matrix} $\b Q$. Indeed, suppose we apply the above method to the system $\dot{\b x}=(\b A-\b B\b K_j)\b x+\b B\b u_j$. After finding an appropriate weight matrix $\b Q_{j+1}$, we solve the Riccati equation
\begin{equation}
(\b A-\b B\b K_j)^{\top}\b P_{j+1}+\b P_{j+1}(\b A-\b B\b K_j)+\b Q_{j+1}-\b P_{j+1}\b B\b R^{-1}\b B^{\top}\b P_{j+1}=\b O
\label{eq:secricceq}
\end{equation}
for $\b P_{j+1}$ to compute the gain $\b K_{j+1}$. Adding the two Riccati equations~(\ref{eq:onericcate}) and~(\ref{eq:secricceq}), and performing simple manipulations yields
\[
\b A^{\top}(\b P_j+\b P_{j+1})+(\b P_j+\b P_{j+1})\b A+\b Q_j+\b Q_{j+1}-(\b P_j+\b P_{j+1})\b B\b R^{-1}\b B^{\top}(\b P_j+\b P_{j+1})=\b O.
\]
Let $\b P=\b P_j+\b P_{j+1}$ and $\b Q=\b Q_j+\b Q_{j+1}$, then
\begin{eqnarray*}
\b K&=&\b K_j+\b K_{j+1}\\
&=& \b R^{-1}\b B^{\top}\b P_j+\b R^{-1}\b B^{\top}\b P_{j+1}.
\end{eqnarray*}
In a similar manner we can show that the method works for more than just two iterations.
\begin{example}\label{example:onelink}\rm
We consider a model of one-link robot manipulator\index{one-link robot} presented in Subsection~\ref{subsection:one_link_robot}. Applying Kirchhoff's voltage law\index{Kirchhoff's voltage law} to the armature circuit\index{armature circuit} yields
\[
L_a\frac{di_a}{dt}+R_a i_a+K_bN\frac{d\theta_p}{dt}=u,
\]
where $K_b$ is the back emf constant. We assume now that $L_a\approx 0$. Hence,
\begin{equation}
u=R_ai_a+K_bN\frac{d\theta_p}{dt},
\label{eq:kirchhoff_v}
\end{equation}
We next compute $i_a$ from~(\ref{eq:kirchhoff_v}) and substitute the result into~(\ref{eq:newton_t}) to get
\begin{equation}
ml^2\frac{d^2\theta_p}{dt^2}=mgl\sin \theta_p + NK_m\left(\frac{u}{R_a}-\frac{K_bN\frac{d\theta_p}{dt}}{R_a} \right).
\label{eq:second_o}
\end{equation}
We construct a reduced-order state-space model\index{reduced-order model!of one-link robot} of the one-link robot manipulator. For this we choose the following state and output variables:
\[
x_1=\theta_p,\quad x_2=\frac{d\theta_p}{dt}=\omega_p,\quad\mb{and} \quad y=x_1.
\]
Then, using~(\ref{eq:second_o}), we obtain the following simple state-space model of the robot manipulator:
\begin{eqnarray*}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array} \right] &=& \left[\begin{array}{c}
x_2\\
\frac{g}{l}\sin x_1-\frac{K_bK_mN^2}{ml^2R_a}x_2+\frac{NK_m}{ml^2R_a}u
\end{array}\right]\\
y&=& x_1.
\end{eqnarray*}
Reasonable parameters for the robot\index{robot parameters} are: $l=1\,\mb{m}$, $m=1\,\mb{kg}$, $N=10$, $K_m=0.1\,\mb{Nm}/\mb{A}$, $K_b=0.1\,\mb{Vs}/\mb{rad}$, $R_a=1\,\Omega$. With the above parameters the robot model takes the form:
\begin{eqnarray*}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array} \right] &=& \left[\begin{array}{c}
x_2\\
9.8\,\sin x_1-x_2+u
\end{array}\right]\\
y&=& x_1.
\end{eqnarray*}
Time histories of state trajectories of the uncontrolled nonlinear system, $u=0$, are shown in Figure~\ref{fig:hw4p1st} for the initial conditions are $x_1(0)=1$ and $x_2(0)=0$.
\begin{figure}[h]
\centering
\includegraphics[scale=0.6]{hw4p1st.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p1st.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p1st.bmp x=4.67in y=3.5in}
\caption{Plots of $y=x_1$ and $x_2$ versus time of the uncontrolled nonlinear system of Example~\ref{example:onelink}.}
\label{fig:hw4p1st}
\end{figure}
A phase portrait\index{phase portrait} of the uncontrolled nonlinear system is shown in Figure~\ref{fig:hw4p1ph}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p1ph.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p1ph.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p1ph.bmp x=4.67in y=3.5in}
\caption{A phase portrait of the uncontrolled nonlinear system of Example~\ref{example:onelink}.}
\label{fig:hw4p1ph}
\end{figure}
The linearized model about $\b x=\b 0, u=0$ has the form
\begin{equation}
\left.\begin{array}{ccl}
\frac{d}{dt}\Delta \b x &=& \left[\begin{array}{cc}
0 & 1\\
9.8 & -1
\end{array}\right] \Delta \b x + \left[\begin{array}{c}
0\\
1
\end{array}\right]\Delta u\\
\Delta y &=& \left[\begin{array}{cc}
1 & 0
\end{array}\right]\Delta \b x \end{array}\right\}
\label{eq:limoforop}
\end{equation}
Plots of $y=x_1$ and $x_2$ versus time are shown in Figure~\ref{fig:hw4p2st}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p2st.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p2st.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p2st.bmp x=4.67in y=3.5in}
\caption{Plots of $x_1$ and $x_2$ versus time of the linearized system of Example~\ref{example:onelink}.}
\label{fig:hw4p2st}
\end{figure}
A phase portrait of the linearized system is shown in Figure~\ref{fig:hw4p2pp}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p2pp.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p2pp.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p2pp.bmp x=4.67in y=3.5in}
\caption{A phase portrait of the linearized system of Example~\ref{example:onelink}.}
\label{fig:hw4p2pp}
\end{figure}
Let
\[
J=\int_0^{\infty}\left( y^2+u^2\right)dt.
\]
We find the linear state-feedback control law $u=-\b k\b x$ that minimizes $J$ subject to the equations given by~(\ref{eq:limoforop}). We have
\[
\b Q=\b c^{\top}\b c=\left[\begin{array}{cc}
1 & 0\\
0 & 0
\end{array}\right]\quad\mb{and}\quad \b R=[1].
\]
Solving the associated Riccati equation\index{Riccati equation!solving} yields
\begin{equation}
\b P=\left[\begin{array}{cc}
72.3371 & 19.6509\\
19.6509 & 5.3484
\end{array}\right].
\label{eq:solofareinolm}
\end{equation}
Hence,
\begin{equation}
\b k=\left[\begin{array}{cc}
19.6509 & 5.3484
\end{array}\right],
\label{eq:optgafolm}
\end{equation}
and
\[
\b A_c=\left[\begin{array}{cc}
0 & 1\\
-9.8509 & -6.3484
\end{array}\right].
\]
Plots of $x_1$ and $x_2$ versus time of the closed-loop linear system
\begin{eqnarray*}
\dot{\b x} &=& \left(\b A-\b b\b k \right)\b x=\b A_c\b x\\
y&=&x_1,
\end{eqnarray*}
when the initial conditions are $x_1(0)=1$ and $x_2(0)=0$ are shown in Figure~\ref{fig:hw4p3st}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.6]{hw4p3st.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p3st.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p3st.bmp x=4.67in y=3.5in}
\caption{Plots of $x_1$ and $x_2$ versus time of the linear closed-loop system of Example~\ref{example:onelink}.}
\label{fig:hw4p3st}
\end{figure}
A phase portrait of the closed-loop linear system\index{phase portrait!of closed-loop linear system} is shown in Figure~\ref{fig:hw4p3pp}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p3pp.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p3pp.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p3pp.bmp x=4.67in y=3.5in}
\caption{A phase portrait of the linear closed-loop system of Example~\ref{example:onelink}.}
\label{fig:hw4p3pp}
\end{figure}
We now apply the above optimal controller to the nonlinear model. Plots of $x_1$ and $x_2$ versus time of the closed-loop nonlinear systems are shown in Figure~\ref{fig:hw4p4st}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p4st.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p4st.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p4st.bmp x=4.67in y=3.5in}
\caption{Plots of $x_1$ and $x_2$ of the nonlinear closed-loop system of Example~\ref{example:onelink}.}
\label{fig:hw4p4st}
\end{figure}
A phase portrait of the nonlinear closed-loop system\index{phase portrait!of closed-loop nonlinear system} is shown in Figure~\ref{fig:hw4p4pp}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hw4p4pp.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/hw4p4pp.bmp}
%\hskip.5in\special{bmp:c:/ee680/hw4p4pp.bmp x=4.67in y=3.5in}
\caption{A phase portrait of the nonlinear closed-loop system of Example~\ref{example:onelink}.}
\label{fig:hw4p4pp}
\end{figure}
The closed-loop poles of the linearized system, that is, the eigenvalues of $\b A_c$ are $\lambda_1=-2.7003$ and $\lambda_2=-3.6481$.
\end{example}
\markex

\begin{example}\rm
We use the data from the previous example to find the weight matrix $\b Q_2$ and the corresponding $\b k_2$ so that $\lambda_2$ is shifted to $s_2=-5$, that is, the new poles of the closed-loop system
\[
\dot{\b x} = \left(\b A_c-\b b\b k_2 \right)\b x
\]
are $s_1=\lambda_1$ and $s_2=-5$. For the control law $u=-\b k\b x+r$, the linear closed-loop system is
\begin{eqnarray*}
\dot{\b x}&=& \left[\begin{array}{cc}
0 & 1\\
-9.8509 & -6.3484
\end{array}\right]\b x+\left[\begin{array}{c}
0 \\
1\end{array}\right]u\\
&=&\b A_c\b x+\b b r,
\end{eqnarray*}
where for simplicity we used $\b x$ rather than $\Delta \b x$, and $u$ instead of $\Delta u$. The eigenvalues of $\b A_c$ are $\lambda_1=-2.7003$ and $\lambda_2=-3.6481$. We have $\lambda_1 \ne \lambda_2$ and thus the matrix $\b A_c$ is diagonalizable. Note that different eigenvalues are only a sufficient condition for diagonalizability\index{different eigenvalues!only sufficient for diagonalizability}. We construct the similarity transformation using the eigenvectors of $\b A_c$ The eigenvectors of $\b A_c$ are
\[
\b v_1=\left[\begin{array}{c}
0.3473\\
-0.9378
\end{array}\right]\quad\mb{and}\quad \b v_2=\left[\begin{array}{c}
-0.2644\\
-0.9644
\end{array}\right].
\]
Let $\b x=\b T\b z$, where $\b T=\left[\begin{array}{cc}
\b v_1&\b v_2
\end{array}\right]$. Then, $\b A_c$ in the new coordinates has the form
\[
\tilde{\b A}_c=\b T^{-1}\b A_c \b T=\left[\begin{array}{cc}
-2.7003 & 0\\
0 & -3.6481
\end{array}\right],
\]
and
\[
\tilde{\b b}=\b T^{-1}\b b=\left[\begin{array}{c}
3.0383\\
3.9912
\end{array}\right].
\]
Next, we compute
\[
\b T^{-1}\b b\b R^{-1}\b b^{\top}\b T^{-\top}=\left[\begin{array}{cc}
9.2312 & 12.1264\\
12.1264 & 15.9296
\end{array}\right].
\]
Hence,
\[
\tilde{q}_{22}=\frac{s_2^2-\lambda_2^2}{15.9296}=0.7339.
\]
Therefore,
\[
\tilde{\b Q}_2=\left[\begin{array}{cc}
0 & 0\\
0 & 0.7339
\end{array}\right]=\b T^{\top}\b Q_2\b T,
\]
and
\[
\b Q_2=\b T^{-\top}\tilde{\b Q}_2\b T^{-1}=\left[\begin{array}{cc}
85.2450 & 31.5686\\
31.5686 & 11.6907
\end{array}\right].
\]
Solving the ARE with $\b A_c$, $\b b$, $\b Q_2$, and $\b R$ yields
\[
\b P_2=\left[\begin{array}{cc}
9.8572 & 3.6504\\
3.6504 & 1.3518
\end{array}\right],
\]
and
\[
\b k_2=\left[\begin{array}{cc}
3.6504 & 1.3518
\end{array}\right].
\]
As a check, if we take
\[
\b Q_{new}=\b Q+\b Q_2=\left[\begin{array}{cc}
86.2450 & 31.5686\\
31.5686 & 11.6907
\end{array}\right],
\]
where $\b Q=\b c^{\top}\b c$ as in Example~\ref{example:onelink}, and $\b R_{new}=\b R$, then solving the associated ARE using $\b A$, rather then $\b A_c$, we obtain
\[
\b P_{new}=\left[\begin{array}{cc}
82.1943 & 23.3013\\
23.3013 & 6.7002
\end{array}\right]=\b P+\b P_2,
\]
where $\b P$ is given by~(\ref{eq:solofareinolm}), and
\[
\b k_{new}=\left[\begin{array}{cc}
23.3013 & 6.7002
\end{array}\right]=\b k+\b k_2,
\]
where $\b k$ is the optimal gain given by~(\ref{eq:optgafolm}).
\end{example}
\markex

\begin{example}\rm
For the linearized model given by~(\ref{eq:limoforop}) and $\b R=[2]$, we determine $\b Q$ and the corresponding $\b k$ that result in the closed-loop system $\dot{\b x}=\left(\b A-\b b\b k \right)\b x$ with poles located at $s_1=-4$ and $s_2=-5$, where
\[
\b A=\left[\begin{array}{cc}
0 & 1\\
9.8 & -1
\end{array}\right],\quad \b b=\left[\begin{array}{c}
0\\
1
\end{array}\right].
\]
We shift poles into the desired locations one by one. We first diagonalize $\b A$ using similarity transformation constructed from eigenvectors of $\b A$. The eigenvalues of $\b A$ are $\lambda_1=2.6702$ and $\lambda_2=-3.6702$. We denote the corresponding eigenvectors by $\b v_1$ and $\b v_2$, respectively. We obtain
\[
\b T^{-1}\b A\b T=\left[\begin{array}{cc}
2.6702 & 0\\
0 & -3.6702
\end{array}\right],
\]
where
\[
\b T=\left[\begin{array}{cc}
\b v_1 & \b v_2
\end{array}\right]=\left[\begin{array}{cc}
0.3507 & -0.2629\\
0.9365 & 0.9648
\end{array}\right].
\]
We construct $\b k_1$ so that $\lambda_1$ is shifted to $s_1=-4$, while $\lambda_2$ remains unchanged, that is, the eigenvalues of $\b A-\b b\b k_1$ become $s_1$ and $\lambda_2$. For this, we first compute
\[
\b F=\b T^{-1}\b b\b R^{-1}\b b^{\top}\b T^{-\top}=\left[\begin{array}{cc}
0.1011 &0.1349\\
0.1349 & 0.1800
\end{array}\right].
\]
Hence, $f_{11}=0.1011$,
\[
\tilde{q}_{11}=\frac{4^2-(2.6702)^2}{0.1011}=87.7352,
\]
and
\[
\tilde{\b Q}_1=\left[\begin{array}{cc}
87.7352 & 0\\
0 & 0
\end{array}\right]=\b T^{\top}\b Q_1\b T.
\]
Hence,
\[
\b Q_1=\b T^{-\top}\tilde{\b Q}_1\b T^{-1}=\left[\begin{array}{cc}
239.0024 & 65.1202\\
65.1202 & 17.7431
\end{array}\right].
\]
Solving the associated ARE, using $\b R=[2]$, yields
\[
\b P_1=\left[\begin{array}{cc}
179.7014 & 48.9626\\
48.9626 & 13.3407
\end{array}\right].
\]
Therefore,
\[
\b k_1=\left[\begin{array}{cc}
24.4813 & 6.6703
\end{array}\right],
\]
and
\[
\b A_1=\b A-\b b\b k_1=\left[\begin{array}{cc}
0 & 1\\
-14.6813 & -7.6703
\end{array}\right].
\]
We construct $\b k_2$ that shifts $\lambda_2$ to $s_2$ leaving $s_1$ unchanged. We first diagonalize $\b A_1$ using its eigenvectors. Note that $\b A_1$ and $\b A$ share the common eigenvector $\b v_2$. The other eigenvector of $\b A_1$ corresponding to its eigenvalue $s_1$ is $\b u_1=\left[\begin{array}{cc}
-0.2425 & 0.9701
\end{array}\right]^{\top}$. We have
\[
\b T^{-1}\b A_1\b T=\left[\begin{array}{cc}
s_1 & 0\\
0 & \lambda_2
\end{array}\right],
\]
where
\[
\b T=\left[\begin{array}{cc}
\b u_1 & \b v_2
\end{array}\right]=\left[\begin{array}{cc}
-0.2425 & 0.2629\\
0.9701 & -0.9648
\end{array}\right].
\]
We next compute
\[
\b T^{-1}\b b\b R^{-1}\b b^{\top}\b T^{-\top}=\left[\begin{array}{cc}
78.0607 & 72.0157\\
72.0157 & 66.4389
\end{array}\right].
\]
Thus, $\tilde{f}_{22}=66.4389$,
\[
\tilde{q}_{22}=\frac{5^2-(-3.6702)^2}{66.4389}=0.1735,
\]
and
\[
\tilde{\b Q}_2=\left[\begin{array}{cc}
0 & 0 \\
0 & 0.1735
\end{array}\right]=\b T^{\top}\b Q_2\b T.
\]
From the above, we obtain
\[
\b Q_2=\b T^{-\top}\tilde{\b Q}_2\b T^{-1}=\left[\begin{array}{cc}
368.9003 & 92.2211\\
92.2211 & 23.0543
\end{array}\right].
\]
Solving the associated ARE, we get
\[
\b P_2=\left[\begin{array}{cc}
42.5495 & 10.6369\\
10.6369 & 2.6591
\end{array}\right],
\]
and
\[
\b k_2=\left[\begin{array}{cc}
5.3185 & 1.3296
\end{array}\right].
\]
The final weight matrix\index{final weight matrix} $\b Q$ that results in shifting poles into the desired locations is
\[
\b Q=\b Q_1+\b Q_2=\left[\begin{array}{cc}
607.9027 & 157.3413\\
157.3413 & 40.7974
\end{array}\right].
\]
The corresponding $\b P$ is
\[
\b P=\b P_1+\b P_2=\left[\begin{array}{cc}
222.2509 & 59.5996\\
59.5996 & 15.9998
\end{array}\right].
\]
The corresponding optimal gain vector\index{optimal gain vector} $\b k$ is
\[
\b k=\b k_1+\b k_2=\left[\begin{array}{cc}
29.7998 & 7.9999
\end{array}\right].
\]
\end{example}
\markex

\subsection{Optimal Saturating Controllers}
Asymptotic stability of an equilibrium state of a nonlinear system is a local property.\index{asymptotic stability!of a nonlinear system a local property} This means that the system's trajectories starting in a neighborhood of the equilibrium state will approach this equilibrium asymptotically. Often, we are interested in determining how far from the equilibrium state the system's trajectories can start and still converge to the equilibrium asymptotically. The region around an equilibrium state within which all trajectories approach the equilibrium is called a {\em region of attraction\/}\index{region of attraction}, or {\em region of asymptotic stability\/}\index{region of asymptotic stability (RAS)}\index{RAS|see{region of asymptotic stability}} abbreviated as RAS. Finding the exact RAS analytically can be a challenging task. A tool often used to estimate the exact RAS is a Lyapunov function. Specifically, if $V=V(\b x)$ is a Lyapunov function of the equilibrium state of interest, that is, $V$ is positive definite with respect to the equilibrium state in a domain ${\cal S}$ and the time derivative of $V$ evaluated on the trajectories of the system is negative definite in ${\cal S}$, then the set
\[
\Omega_c=\left\{ \b x:\b x\in {\cal S}, V(\b x) \le c, \dot{V}(\b x)<0 \right\}
\]
is an estimate of the exact RAS\index{estimate!of the exact RAS}. In other words, a region of attraction is a set of states $\b x$ for which $\dot{V}(\b x)<0$ and $V(\b x)\le c$ for some positive constant $c$. The set $\Omega_c$ may be a conservative estimate, often a crude underestimate, of the true RAS. In order to obtain the least conservative estimate of the domain of attraction, corresponding to the given Lyapunov function, we could try to maximize the size of $\Omega_c$ by finding the largest $V$ contour that can just fit on the equilibrium side of the surface described by $\dot{V}(\b x)=0$. This can be accomplished by solving the optimization problem for $c$ such that
\begin{equation}
c=\min_{\dot{V}(\footnotesize{\b x})=0} V(\b x).
\label{eq:maxiconofv}
\end{equation}
Solving the above optimization problem corresponds to maximizing the size of the region $\{\b x:V(\b x)\le c\}$ in which $\dot{V}(\b x)<0$ by making the boundary of $\{\b x:V(\b x)\le c\}$ to be tangent to the boundary of $\{\b x:\dot{V}(\b x) < 0\}$. This is illustrated in Figure~\ref{fig:rasmaxesti}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{rasesti.eps}
%\vspace{3in}
%\hskip1.2in\special{bmp:c:/mybook/rasesti.bmp x=4in y=3in}
\caption{Estimating the region of asymptotic stability.}
\label{fig:rasmaxesti}
\end{figure}

We next estimate the RAS for single-input linear dynamical systems\index{RAS estimation!for single-input systems} with saturating linear optimal control. Our development follows that of Friedland~\cite[Section~5.1]{friedland96}. We will use such controllers in Chapter~\ref{chp:chaos} to stabilize chaotic systems. We consider a single-input linear system model,
\begin{equation}
\dot{\b x}=\b A\b x+\b b \;U_{max}\mb{sat}\left(u/U_{max}\right),
\label{eq:sisywlinsat}
\end{equation}
where $\b x\in \mb{\Bb R}^n$ and
\[
\mb{sat}(z)=\left\{\begin{array}{ccl}
1&\mb{if}&z\ge 1\\
z&\mb{if}&|z|<1\\
-1&\mb{if}&z\le -1.
\end{array}\right.
\]
In the linear region, $|u|\le U_{max}$, (\ref{eq:sisywlinsat}) becomes $\dot{\b x}=\b A\b x+\b bu$. We consider using the control law of the form
\begin{equation}
u=U_{max}\mb{sat}\left( -\b k \b x/U_{max}\right).
\label{eq:satstafbcl}
\end{equation}
When $|\b k\b x|\le U_{max}$, the control law given by~(\ref{eq:satstafbcl}) becomes a linear state-feedback controller, $u=-\b k\b x$. For simplicity of notation, we can assume that $U_{max}=1$. Indeed, we can always scale the input vector $\b b$ so that $U_{max}=1$. The gain vector $\b k$ is obtained by solving the associated LQR problem, with the performance index,
\[
J=\int_0^{\infty}\left( \b x^{\top}\b Q\b x+ru^2\right)dt,
\]
where we assume that the pair $(\b A,\b b)$ is reachable, $\b Q=\b Q^{\top}> 0$ and $r>0$. Solving the algebraic Riccati equation,
\begin{equation}
\b A^{\top}\b P+\b P\b A+\b Q-\b P\b b\b b^{\top}\b P/r=\b O,
\label{eq:ricceforstacon}
\end{equation}
we use its solution, $\b P=\b P^{\top}>0$, to construct the gain vector,
\begin{equation}
\b k=\b b^{\top}\b P/r.
\label{eq:opgafostacon}
\end{equation}
We then proceed with evaluating of $\frac{d}{dt}V=\frac{d}{dt}\b x^{\top}\b P\b x$ on the trajectories of the closed-loop system,
\begin{eqnarray}
\dot{V} &=& 2\b x^{\top}\b P\dot{\b x}\nonumber\\
&=& \b x^{\top}(\b A^{\top}\b P+\b P\b A)\b x- 2\b x^{\top}\b P\b b\:\mb{sat}(\b k\b x).
\label{eq:lyadrotrclsonsts}
\end{eqnarray}
By~(\ref{eq:opgafostacon}),
\begin{equation}
\b b^{\top}\b P=r\b k.
\label{eq:modifromopgasa}
\end{equation}
Using the above and~(\ref{eq:ricceforstacon}) gives
\begin{equation}
\b A^{\top}\b P+\b P\b A=r\b k^{\top}\b k-\b Q.
\label{eq:modriccfosub}
\end{equation}
Substituting~(\ref{eq:modriccfosub}) and~(\ref{eq:modifromopgasa}) into~(\ref{eq:lyadrotrclsonsts}) yields
\begin{equation}
\fbox{$ \displaystyle \dot{V}=-\b x^{\top}\b Q\b x-2r\b k\b x\:\mb{sat}(\b k \b x) +r(\b k\b x)^2 $}
\label{eq:exfovdosa}
\end{equation}
For $\b x$ such that $|\b k\b x|\le 1$, the above becomes
\[
\dot{V}=-\b x^{\top}\b Q\b x-r(\b k\b x)^2,
\]
which is negative definite. Thus, the origin is (locally) asymptotically stable for the closed-loop system.

Recall that an estimate of the RAS is a set of states bounded by the surface $V=\mb{constant}$ for which $\dot{V}<0$. To find as large an estimate as possible, we attempt to find the surface $\dot{V}=0$ and find $c=\min V(\b x)$ on that surface. Then, the region bounded by $V(\b x)=c$ is an estimate of the RAS\index{RAS!estimation of}. Solving the above problem amounts to solving the optimization problem,
\[
\min_{\footnotesize{\b x}}\b x^{\top}\b P\b x
\]
\begin{equation}
\mb{subject to}\quad
-\b x^{\top}\b Q\b x-2r\b k \b x\:\mb{sat}(\b k\b x)+r(\b k\b x)^2=0.
\label{eq:constrforras}
\end{equation}
Solving the above problem is not a simple task. Instead, we solve an easier problem resulting in a more conservative estimate of the RAS. We first rewrite~(\ref{eq:constrforras}) as
\begin{equation}
y^2-2y\:\mb{sat}(y)- \b x^{\top}\b Q\b x/r=0,
\label{eq:rewrittconstras}
\end{equation}
where $y=\b k\b x$. We investigate the case when $|y|\ge 1$ because for $|y|<1$, we have $\dot{V}<0$ and we are interested in the solution to the above when $\dot{V}=0$. When $y\ge 1$, (\ref{eq:rewrittconstras}) becomes
\begin{equation}
y^2-2y- \b x^{\top}\b Q\b x/r=0.
\label{eq:abfoymaggo}
\end{equation}
because for $y\ge 1$, $\mb{sat}(y)=1$. The solution to~(\ref{eq:abfoymaggo}) for $y\ge 1$ is
\begin{equation}
y=1+\sqrt{1+\b x^{\top}\b Q\b x/r}>2,
\label{eq:firsotoseo}
\end{equation}
while for $y\le -1$, the solution to~(\ref{eq:rewrittconstras}) has the form
\begin{equation}
 y=-1-\sqrt{1+\b x^{\top}\b Q\b x/r}<-2.
\label{eq:secsoltosat}
\end{equation}
The above solutions satisfy the inequality,
\[
|y|\ge 2.
\]
Thus, we have $\dot{V}<0$ for $|y|=|\b k \b x|\le 2$ and hence the RAS includes the region
\[
\{\b x:V(\b x)\le \tilde{c}\},
\]
where $\tilde{c}$ is obtained by solving the optimization problem
\[
\fbox{$\displaystyle \tilde{c}=\max_{|\footnotesize{\b k\b x}|\le 2}\b x^{\top}\b P\b x $}
\]
We use the Lagrange multipliers method\index{Lagrange multipliers method!RAS estimate} to solve the above constrained optimization problem\index{constrained optimization problem!RAS estimate}. We solve the equivalent problem,
\[
\max \b x^{\top}\b P\b x
\]
\[
\mb{subject to}\quad(\b k\b x)^2=4.
\]
The Lagrangian\index{Lagrangian} for the above problem is
\[
l(\b x,\lambda)=\b x^{\top}\b P\b x+\lambda\left(4-(\b k\b x)^2\right).
\]
Applying to the Lagrangian $l$ the first-order necessary conditions for unconstrained minimization\index{first-order necessary conditions!for unconstrained minimization} gives
\begin{eqnarray}
D_{\footnotesize{\b x}}l(\b x,\lambda) &=& 2\b x^{\top}\b P-2\lambda(\b k\b x)\b k=\b 0, \label{eq:fineccofol}\\
D_{\lambda}l(\b x,\lambda) &=& 4-(\b k\b x)^2=0.
\label{eq:secneccofl}
\end{eqnarray}
From~(\ref{eq:fineccofol}),
\[
\b x^{\top}\b P\b x=\lambda(\b k\b x)^2.
\]
Because $(\b k\b x)^2=4$,
\[
\lambda=\frac{1}{4}\b x^{\top}\b P\b x,
\]
and therefore,
\[
\tilde{c}=\b x^{\top}\b P\b x=4\lambda.
\]
Again, using~(\ref{eq:fineccofol}), we obtain
\[
\b x^{\top}\b P=\lambda (\b k\b x)\b k.
\]
Postmultiplying both sides of the above by $\b P^{-1}\b k^{\top}$ gives
\begin{equation}
\b x^{\top}\b k^{\top}=\b k\b x=\lambda (\b k\b x)\b k\b P^{-1}\b k^{\top}.
\label{eq:comlamuspk}
\end{equation}
Note that $\b k\b x\ne 0$ because $(\b k\b x)^2=4$. Hence, from~(\ref{eq:comlamuspk}),
\[
\lambda=\frac{1}{\b k\b P^{-1}\b k^{\top}}.
\]
Substituting~(\ref{eq:opgafostacon}) into the above gives
\[
\lambda=\frac{1}{\b k\b P^{-1}\b k^{\top}}=\frac{r^2}{\b b^{\top}\b P\b b}.
\]
Hence,
\[
\fbox{$\displaystyle \tilde{c}=4\lambda=\frac{4r^2}{\b b^{\top}\b P\b b} $}
\]


\subsection{Linear Quadratic Regulator for Discrete Systems on an Infinite Time Interval}
Consider a discrete linear plant\index{discrete linear plant} model,
\begin{equation}
\b x[k+1]=\b A\b x[k]+\b B\b u[k],\quad k=0,1,2,\ldots,
\label{eq:ifhodlqr}
\end{equation}
with a specified initial condition $\b x(0)=\b x_0$, where $\b x[k]\in \mb{\Bb R}^n$ and $\b u[k]\in \mb{\Bb R}^m$. We assume that the pair $(\b A,\b B)$ is reachable. Our objective is to construct a stabilizing linear state-feedback controller,
\[
\b u[k]=-\b K\b x[k]
\]
that minimizes the quadratic performance index
\begin{equation}
J(\b u)=\sum_{k=0}^{\infty}\left\{\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\right\},
\label{eq:miperdisind}
\end{equation}
where $\b Q=\b Q^{\top}\ge 0$, and $\b R=\b R^{\top}>0$. We assume that the components of the control vector are unconstrained\index{control!unconstrained}. An optimal control\index{optimal control} law will be denoted $\b u^*$.

Because we seek the controller among the ones that stabilize the system~(\ref{eq:ifhodlqr}), the closed-loop system,
\begin{equation}
\b x[k+1]=\left(\b A-\b B\b K\right)\b x[k],
\label{eq:clsydlqr}
\end{equation}
is asymptotically stable, that is, the eigenvalues of the matrix $\b A-\b B\b K$ are in the open unit disk. This implies that there exists a Lyapunov function, $V(\b x[k])= \b x[k]^{\top}\b P\b x[k]$, for the system~(\ref{eq:clsydlqr}). Therefore the first forward difference\index{first forward difference}, $\Delta V(\b x[k])=V(\b x[k+1])-V(\b x[k])$ evaluated on the trajectories of~(\ref{eq:clsydlqr}) is negative definite. We now state and prove a sufficient condition for $\b u[k]=-\b K\b x[k]$ to be optimal.
\begin{thm}
\label{thm:dlqrinfhor}
If the state feedback controller, $\b u[k]=-\b K\b x[k]$, is such that
\begin{equation}
\min_{\b u}\left( \Delta V(\b x[k]) +\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\right)=0,
\label{eq:mintodare}
\end{equation}
then it is optimal.
\end{thm}
{\bf Proof}\quad We can represent~(\ref{eq:mintodare}) as
\[
\left.\Delta V(\b x[k])\right\arrowvert_{\b u=\b u^{*}}+\b x[k]^{\top}\b Q\b x[k]+\b u^*[k]^{\top}\b R\b u^*[k]=0.
\]
Hence,
\[
\left.\Delta V(\b x[k])\right\arrowvert_{\b u=\b u^{*}}=\left.\left(V(\b x[k+1])-V(\b x[k]\right)\right\arrowvert_{\b u=\b u^*}=-\b x[k]^{\top}\b Q\b x[k]-\b u^*[k]^{\top}\b R\b u^*[k].
\]
We sum both sides from $k=0$ to $k=\infty$ to get
\[
V(\b x[\infty])-V(\b x[0])=-\sum_{k=0}^{\infty}\left( \b x[k]^{\top}\b Q\b x[k]+\b u^*[k]^{\top}\b R\b u^*[k]\right).
\]
Because by assumption the closed-loop system is asymptotically stable, we have $\b x[\infty]=\b 0$. Hence,
\[
V(\b x[0])=\b x[0]^{\top}\b P\b x[0]=\sum_{k=0}^{\infty}\left( \b x[k]^{\top}\b Q\b x[k]+\b u^*[k]^{\top}\b R\b u^*[k]\right).
\]
Thus, if a linear stabilizing controller satisfies the hypothesis of the theorem, then the value of the performance index~(\ref{eq:miperdisind}) resulting from applying this controller is
\[
J=\b x[0]^{\top}\b P\b x[0].
\]
To show that such a controller is an optimal one, we use a proof by contradiction. We assume that the hypothesis of the theorem holds but the controller $\b u^*$ that satisfies this hypothesis is not optimal, that is, there is a controller $\tilde{\b u}$ such that
\begin{equation}
J\left(\tilde{\b u}\right)<J\left(\b u^*\right).
\label{eq:jtiulejus}
\end{equation}
The hypothesis of the theorem implies that
\begin{equation}
\left.\Delta V(\b x[k])\right\arrowvert_{\b u=\tilde{\b u}}+\b x[k]^{\top}\b Q\b x[k]+\tilde{\b u}[k]^{\top}\b R\tilde{\b u}[k]\ge 0.
\label{eq:getoeqfuti}
\end{equation}
We represent~(\ref{eq:getoeqfuti}) as
\[
\left.\Delta V(\b x[k])\right\arrowvert_{\b u=\tilde{\b u}}\ge -\b x[k]^{\top}\b Q\b x[k]-\tilde{\b u}[k]^{\top}\b R\tilde{\b u}[k].
\]
Summing the above from $k=0$ to $k=\infty$ yields
\[
V(\b x[0])\le \sum_{k=0}^{\infty}\left( \b x[k]^{\top}\b Q\b x[k]+\tilde{\b u}[k]^{\top}\b R\tilde{\b u}[k]\right),
\]
that is,
\[
J\left(\tilde{\b u}\right)\ge J\left(\b u^*\right),
\]
which contradicts~(\ref{eq:jtiulejus}). The proof is complete.\\
\hspace*{5.5in} $\Box$

Finding an optimal controller involves finding an appropriate quadratic Lyapunov function $V(\b x)=\b x^{\top}\b P\b x$, which then is used to construct the optimal controller. We first find $\b u^*$ that minimizes the function
\begin{equation}
f=f(\b u[k])=\Delta V(\b x[k]) +\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k].
\label{eq:futbemin}
\end{equation}
To proceed, we perform preliminary manipulations on the function $f$. In particular, we use the definition of $\Delta V$ and substitute into~(\ref{eq:futbemin}) the system equation~(\ref{eq:ifhodlqr}) to obtain
\begin{eqnarray}
f(\b u[k])&=&\b x[k+1]^{\top}\b P\b x[k+1]-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\nonumber\\
&=&\left( \b A\b x[k]+\b B\b u[k]\right)^{\top}\b P\left( \b A\b x[k]+\b B\b u[k]\right)-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k]\nonumber\\
&{}&\quad+\b u[k]^{\top}\b R\b u[k].
\label{eq:mifuafsubs}
\end{eqnarray}
Applying to the above function the first-order necessary condition for a relative minimum gives
\begin{eqnarray*}
\frac{\partial f(\b u[k])}{\partial \b u[k]}&=&2\left(\b A \b x[k] + \b B \b u[k]\right)^{\top} \b P \b B+2\b u[k]^{\top}\b R\\
&=&2\b x[k]^{\top}\b A^{\top}\b P\b B+2\b u[k]^{\top}\left( \b R+\b B^{\top}\b P\b B\right)\\
&=&\b 0.
\label{eq:fonctofbu}
\end{eqnarray*}
The matrix $\b R+\b B^{\top}\b P\b B$ is positive definite because $\b R$ is, and therefore the matrix $\b R+\b B^{\top}\b P\b B$ is invertible. Hence,
\begin{equation}
\b u^*[k]=-\left( \b R+\b B^{\top}\b P\b B\right)^{-1}\b B^{\top}\b P\b A\b x[k]=-\b K\b x[k],
\label{eq:leforoptdlqrco}
\end{equation}
where
\[
\b K=\left( \b R+\b B^{\top}\b P\b B\right)^{-1}\b B^{\top}\b P\b A.
\]
Let
\[
\b S=\b R+\b B^{\top}\b P\b B.
\]
Then, we represent~(\ref{eq:leforoptdlqrco}) as
\begin{equation}
\b u^*[k]=-\b S^{-1}\b B^{\top}\b P\b A\b x[k].
\label{eq:optdlqrco}
\end{equation}
We next check if $\b u^*$ satisfies the second-order sufficient condition for a relative minimizer\index{second-order sufficient condition!for a relative minimizer} of the function~(\ref{eq:futbemin}). We have,
\begin{eqnarray}
&{}&\frac{\partial ^2}{\partial \b u^2(k)}\left( \b x[k+1]^{\top}\b P\b x[k+1]-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\right)\nonumber\\
&{}&\quad =\frac{\partial}{\partial \b u[k]}\left( 2\b x[k]^{\top}\b A^{\top}\b P\b B+2\b u[k]^{\top}\left( \b R+\b B^{\top}\b P\b B\right)\right)\nonumber\\
&{}&\quad = 2\left( \b R+\b B^{\top}\b P\b B\right)\nonumber\\
&{}&\quad >0,
\label{eq:soscfust}
\end{eqnarray}
that is, $\b u^*$ satisfies the second-order sufficient condition for a relative minimizer of the function~(\ref{eq:futbemin}).

The optimal controller~(\ref{eq:optdlqrco}) can be constructed if we have found an appropriate positive definite matrix $\b P$. Our next task is to devise a method that would allow us to compute the desired matrix $\b P$. For this, we first find the equation describing the closed-loop system driven by the optimal controller given by~(\ref{eq:optdlqrco}),
\begin{equation}
\b x[k+1]=\left( \b A-\b B\b S^{-1}\b B^{\top}\b P\b A \right)\b x[k].
\label{eq:cloopsysdt}
\end{equation}
Our controller satisfies the hypothesis of Theorem~\ref{thm:dlqrinfhor}, that is,
\[
\b x[k+1]^{\top}\b P\b x[k+1]-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k] +\b u^*[k]^{\top}\b R\b u^*[k]=0.
\]
Substituting into the above~(\ref{eq:cloopsysdt}), the expression for $\b u^*$ given by~(\ref{eq:optdlqrco}), and performing simple manipulations gives
\begin{eqnarray}
&{}&\b x[k]^{\top}\left( \b A-\b B\b S^{-1}\b B^{\top}\b P\b A\right)^{\top}\b P\left( \b A-\b B\b S^{-1}\b B^{\top}\b P\b A\right)\b x[k]-\b x[k]^{\top}\b P\b x[k]\nonumber\\
&{}&\quad +\b x[k]^{\top}\b Q\b x[k] +\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b R\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\quad =\b x[k]^{\top}\b A^{\top}\b P\b A\b x[k]-\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\qquad -\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]+\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\qquad -\b x[k]^{\top}\b P\b x[k]+ \b x[k]^{\top}\b Q\b x[k] +\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b R\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\quad =\b x[k]^{\top}\b A^{\top}\b P\b A\b x[k]-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k]\nonumber\\
&{}&\qquad -2\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\qquad +\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\left(\b R+\b B^{\top}\b P\b B \right)\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\quad =\b x[k]^{\top}\b A^{\top}\b P\b A\b x[k]-\b x[k]^{\top}\b P\b x[k]+\b x[k]^{\top}\b Q\b x[k]\nonumber\\
&{}&\qquad -2\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\qquad +\b x[k]^{\top}\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A\b x[k]\nonumber\\
&{}&\quad = \b x[k]^{\top}\left(\b A^{\top}\b P\b A-\b P+\b Q-\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A \right)\b x[k]\nonumber\\
&{}&\quad = 0.
\label{eq:manidare}
\end{eqnarray}
The above equation should hold for any $\b x$. For this to be true we have to have
\begin{equation}
\fbox{$\displaystyle \b A^{\top}\b P\b A-\b P+\b Q-\b A^{\top}\b P\b B\b S^{-1}\b B^{\top}\b P\b A=\b O $}
\label{eq:dtlqrrieq}
\end{equation}
The above equation is called the {\em discrete algebraic Riccati equation}\index{discrete algebraic Riccati equation}\index{Riccati equation!discrete}.


\section{Dynamic Programming}
The objective of this section is to introduce the reader to the method of dynamic programming\index{dynamic programming} invented by Richard E. Bellman (1920--1984) in 1957. Dynamic programming can be used to synthesize optimal controllers for nonlinear, time-varying dynamical systems. We first present dynamic programming for discrete-time systems and then discuss the continuous version of dynamic programming.

\subsection{Discrete-Time Systems}
Dynamic programming is based on the principle of optimality~(PO)\index{principle of optimality (PO)}\index{PO|see{principle of optimality}}\index{Bellman's principle of optimality (PO)}, that was formulated by Bellman as follows:

\framebox{\parbox[t]{6in}{\em An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision.}}

\noindent In other words, an optimal trajectory has the property that at an intermediate point, no matter how it was reached, the rest of the trajectory must coincide with an optimal trajectory as computed from this intermediate point as the initial point. We justify the PO by contradiction. Suppose that we are given an optimal trajectory starting from initial state $\b x_0=\b x(t_0)$ at time $t_0$ and terminating at the final state $\b x(t_f)$ at time $t_f$. Let $\b x(t)$ be an intermediate state on the optimal trajectory and let $J^{*}(t,\b x(t))$ be the optimal cost of traveling from $\b x(t)$ to the final state $\b x(t_f)$. Thus, $J^{*}(t_0,\b x(t_0))$ is the optimal cost of traveling from $\b x_0=\b x(t_0)$ to $\b x(t_f)$. It then follows that the portion of the optimal trajectory from $\b x(t)$ to $\b x(t_f)$ is also optimal. For if some other path, indicated by dashed curve in Figure~\ref{fig_po}, connecting $\b x(t)$  and $\b x(t_f)$ resulted in a smaller cost than the solid path, then this would provide a less expensive route from $\b x_0$ to $\b x(t_f)$, thus contradicting the optimality of the original trajectory.
\begin{figure}
\centering
\includegraphics[scale=0.6]{fig_po.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/fig_po.bmp}
%\hskip.7in\special{bmp:c:/mybook/fig_po.bmp x=4.67in y=3.5in}
\caption{Illustration of the Principle of Optimality.}
\label{fig_po}
\end{figure}

We note that the PO seems to be already known to Johann Bernoulli who stated in 1706 that ``... all curves, which should give a maximum, preserve also in all their sections the laws of the same maximum''. Here, the reader should interpret maximum as optimum.

Bellman's principle of optimality\index{Bellman's principle of optimality (PO)} serves to ``limit the number of potentially optimal control strategies that must be investigated.''~\cite[page~315]{lewis95}. The optimal control strategy is determined, using the PO, by working backward from the final stage.

We illustrate this ``backwardness'' nature of the PO on an example of a sequential decision making\index{sequential decision making process} process. We define a decision as the choice of alternative paths leaving a given node.
\begin{example} \rm
\label{ex_route}
We use the PO to find the minimum cost path that starts at the node {\bf a} and ends at any one of the nodes {\bf k}, {\bf l}, {\bf m} in the routing network depicted in Figure~\ref{fig_route}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{route.eps}
%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/mybook/route.bmp}
%\hskip.7in\special{bmp:c:/mybook/route.bmp x=4.67in y=3.5in}
\caption{A routing network.}
\label{fig_route}
\end{figure}
The travel costs are shown beside each path segment. The line segments can only be traversed from left to right. Applying the PO to this routing problem\index{routing problem}, we find the optimal path by performing a backward pass\index{backward pass} through the network. Let $J^{*}_q$ be the minimum cost from a generic node {\bf q} to any one of the three possible final nodes. Then, we have
\[
J^{*}_g=16,\quad J^{*}_h=8,\quad J^{*}_i=7,\quad J^{*}_j=6.
\]
The next stage, yields
\[
J^{*}_d = \min\{3+J^{*}_g,11+J^{*}_h\}= \min\{19,19\}=19.
\]
Thus, in this case, we have two optimal paths from node {\bf d}. We can take either {\bf dgk} or {\bf dhk} paths. We then compute optimal paths from nodes {\bf e} and {\bf f}. We have
\[
J^{*}_e = \min\{6+J^{*}_h,10+J^{*}_i\}= \min\{6+8,10+7\}=14,
\]
and
\[
J^{*}_f = \min\{4+J^{*}_i,12+J^{*}_j\}= \min\{4+7,12+6\}=11.
\]
Thus, the optimal path emanating from node {\bf e} is {\bf ehk}, and from node {\bf f} is {\bf fil}. The next stage yields
\[
J^{*}_b = \min\{9+J^{*}_d,2+J^{*}_e\}= \min\{9+19,2+14\}=16,
\]
and
\[
J^{*}_c = \min\{5+J^{*}_e,14+J^{*}_f\}= \min\{5+14,14+11\}=19.
\]
Therefore, the optimal path starting from {\bf b} is {\bf behk}, and starting from {\bf c} the optimal path is
{\bf cehk}. Finally, the initial, and hence final, stage yields
\[
J^{*}_a = \min\{4+J^{*}_b,7+J^{*}_c\}= \min\{4+16,7+19\}=20.
\]
Hence, the optimal, minimum total cost, route from {\bf a} to any one of the three possible termination nodes {\bf k}, {\bf l}, or {\bf m} is {\bf abehk}. The total cost of this path is 20.
\end{example}\markex

We next consider a class of discrete-time dynamical systems modeled by the difference equation\index{difference equation}
\begin{equation}
\b x[k+1]=\b f\left(k,\b x[k],\b u[k]\right),
\label{eq:dtdysylqr}
\end{equation}
with a given initial condition $\b x[k_0]=\b x_0$, and the associated performance index
\[
J=J_0=\Phi(N,\b x[N])+\sum_{k=0}^{N-1}F(k,\b x[k],\b u[k]).
\]
Let $J^{*}_k(\b x[k])$ denote the minimum cost of transferring\index{minimum cost!of transferring} the system from $\b x[k]$ to the terminal point $\b x[N]$. Because $\b x[N]$ is the terminal point, therefore $J^{*}_N(\b x[N])=\Phi(N,\b x[N])$. Using the PO, we obtain the cost of transfer from $\b x[N-1]$ to the terminal point,
\begin{eqnarray*}
J^{*}_{N-1}(\b x[N-1]) &=& \min_{\b u[N-1]}\left\{F(N-1,\b x[N-1],\b u[N-1])+J^{*}_N(\b x[N])\right\}\\
&=& \min_{\b u[N-1]}\left\{F(N-1,\b x[N-1],\b u[N-1])+\Phi(N,\b x[N])\right\}.
\end{eqnarray*}
When applying the PO to solving optimal control problems for discrete-time dynamical systems, we usually perform two stage-by-stage passes\index{two stage-by-stage passes} through the time stages. We begin with a backward pass\index{backward pass}. For this, we use~(\ref{eq:dtdysylqr}) to eliminate $\b x[N]$. Next, we carry out the minimization to find $\b u^{*}[N-1]$. We then repeat the process. A typical step of the backward pass is described as,
\begin{eqnarray*}
J^{*}_k(\b x[k]) &=& \min_{\b u[k]}\left\{F(k,\b x[k],\b u[k])+J^{*}_{k+1}(\b x[k])\right\}\\
&=&\min_{\b u[k]}\left\{F(k,\b x[k],\b u[k])+J^{*}_{k+1}\left(\b f(k,\b x[k],\b u[k])\right)\right\}.
\end{eqnarray*}
The backward pass is completed when the initial time $k_0$ is reached. Now, because $\b x_0=\b x[k_0]$ is known, we find $\b u_0=\b u[k_0]$ in terms of this state. We can then proceed with the {\em forward pass\/}\index{forward pass}. We use $\b x_0,\b u_0$, and~(\ref{eq:dtdysylqr}) to compute $\b x[k_0+1]$. The state $\b x[k_0+1]$ is then used to compute $\b u[k_0+1]$  from $\b x[k_0+2]=\b f(k_0+1,\b x[k_0+1],\b u[k_0+1])$, etc.
\begin{example}\label{examp:dlqr} \rm
Consider a scalar discrete-time dynamical system modeled by
\[
x[k+1]=2x[k]-3u[k],\quad x[0]=4,
\]
and the performance index
\[
J=J_0=(x[2]-10)^2+\frac{1}{2}\sum_{k=0}^{1}\left(x[k]^2+u[k]^2\right).
\]
Note that the final state, in this example, is free. We use the PO to find optimal $u[0]$ and $u[1]$. There are no constraints on $u[k]$. We begin with the backward pass. We have
\[
J^{*}(x[2])=(x[2]-10)^2.
\]
Then,
\begin{eqnarray*}
J^{*}(x[1]) &=& \min_{u[1]}\left\{\frac{1}{2}\left(x[1]^2+u[1]^2\right)+J^{*}(x[2])\right\}\\
&=& \min_{u[1]}\left\{\frac{1}{2}\left(x[1]^2+u[1]^2\right)+\left(2x[1]-3u[1]-10\right)^2\right\}.
\end{eqnarray*}
Because there are no constraints on $u[1]$, we find optimal $u[1]$ as a function of $x[1]$ by applying the first-order necessary condition for unconstrained optimization to get
\[
\frac{\partial}{\partial u[1]}\left\{\frac{1}{2}\left(x[1]^2+u[1]^2\right)+\left(2x[1]-3u[1]-10\right)^2\right\}=0.
\]
Hence,
\[
u[1]=\frac{1}{19}\left(12x[1]-60\right),
\]
and therefore
\[
J^{*}(x[1])=\frac{1}{2}\left(x[1]^2+\left(\frac{12x[1]-60}{19}\right)^2\right)+\left(2x[1]-3\frac{12x[1]-60}{19}-10\right)^2.
\]
Next, we compute $u[0]$. For this observe that
\[
J^{*}(x[0])=\min_{u[0]}\left\{\frac{1}{2}\left(x[0]^2+u[0]^2\right)+J^{*}(x[1])\right\}.
\]
Taking into account that
\[
x[1]=2x[0]-3u[0],\quad x[0]=4,
\]
and the fact that there are no constraints on $u[0]$, we again use the first-order necessary condition for unconstrained optimization\index{first-order necessary condition!for unconstrained optimization} to find optimal $u[0]$. We compute
\[
\frac{\partial}{\partial u[0]}\left\{\frac{1}{2}\left(x[0]^2+u[0]^2\right)+J^{*}(x[1])\right\}=0
\]
to get
\begin{eqnarray*}
&{}&\frac{\partial}{\partial u[0]}\left\{8+\frac{1}{2}u[0]^2+ \frac{1}{2}(8-3u[0])^2+\frac{1}{2}\left(\frac{12(8-3u[0])-60}{19}\right)^2\right\}\\
&{}&\qquad  +\frac{\partial}{\partial u[0]}\left\{2(8-3u[0])-3\frac{12(8-3u[0])-60}{19}-10\right\}^2=0.
\end{eqnarray*}
Performing some manipulations yields $13.7895u[0]=27.7895$. Hence, $u[0]=2.0153$. This completes the backward pass\index{backward pass}. We are now ready for the forward pass\index{forward pass}. We use the system difference equation, $x[0]$, and $u[0]$ to compute $x[1]=1.9541$. Using the above, we find
\[
u[1]= \frac{12x[1]-60}{19}= -1.9237.
\]
Therefore,
\[
x[2] = 2x[1]-3u[1]=9.6794,
\]
which completes the forward pass.
\end{example}\markex


\subsection{Discrete Linear Quadratic Regulator Problem}
In Example~\ref{examp:dlqr} above, we solved a simple discrete linear quadratic regulator\index{discrete linear quadratic regulator} problem where we found the open-loop optimal control strategy. We apply now the PO to find the optimal control of a linear state-feedback form for the discrete linear quadratic regulator problem. This problem is stated as follows. Given a
discrete linear plant model,
\[
\b x[k+1]=\b A\b x[k]+\b B\b u[k],\quad 0\le k\le N-1,
\]
with a specified initial condition $\b x[0]=\b x_0$. We wish to calculate the optimal control sequence\index{optimal control sequence},
\[
\left\{\b u^{*}[0],\b u^{*}[1],\ldots,\b u^{*}[N-1]\right\},
\]
that minimizes the quadratic performance index
\[
J=J_0=\frac{1}{2}\b x[N]^{\top}\b H_N\b x[N]+\frac{1}{2}\sum_{k=0}^{N-1}\left\{\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\right\},
\]
where $\b H_N=\b H^{\top}_N \ge 0$, $\b Q=\b Q^{\top}\ge 0$, and $\b R=\b R^{\top}>0$. We assume that the components of the control vector are unconstrained. To solve this problem we use the PO. To begin, we note that
\[
J^{*}(\b x[N])=\frac{1}{2}\b x[N]^{\top}\b H_N\b x[N].
\]
The above is the penalty for being in a state $\b x[N]$ at a time $N$.

We now decrement $k$ to $N-1$ to get
\[
J_{N-1}(\b x[N-1])=\frac{1}{2}\b x[N]^{\top}\b H_N\b x[N]+\frac{1}{2}\b x[N-1]^{\top}\b Q\b x[N-1]+\frac{1}{2}\b u[N-1]^{\top}\b R\b u[N-1].
\]
We use the state equation to eliminate $\b x[N]$ from $J_{N-1}$ to obtain
\begin{eqnarray*}
J_{N-1} &=& \frac{1}{2}\left(\b A\b x[N-1]+\b B\b u[N-1]\right)^{\top}\b H_N\left(\b A\b x[N-1]+\b B\b u[N-1]\right)\\
&{}&\quad +\frac{1}{2}\b x[N-1]^{\top}\b Q\b x[N-1]+\frac{1}{2}\b u[N-1]^{\top}\b R\b u[N-1].
\end{eqnarray*}
Because there are no constraints on control, we can apply the first-order necessary condition from static optimization to find $\b u^{*}[N-1]$ as a function of $\b x[N-1]$,
\begin{eqnarray*}
\frac{\partial J_{N-1}}{\partial \b u[N-1]} &=& \left(\b A\b x[N-1]+\b B\b u[N-1]\right)^{\top}\b H_N\b B+\b u[N-1]^{\top}\b R\\
&=& \b 0.
\end{eqnarray*}
Solving for the control satisfying the first-order necessary condition, we obtain
\[
\b u^{*}[N-1]=-\left(\b R+\b B^{\top}\b H_N\b B\right)^{-1}\b B^{\top}\b H_N\b A\b x[N-1].
\]
Let
\[
\b K_{N-1}=\left(\b R + \b B^{\top}\b H_N\b B\right)^{-1}\b B^{\top}\b H_N\b A,
\]
then
\[
\b u^{*}[N-1]=-\b K_{N-1}\b x[N-1].
\]
Computing the optimal cost of transferring the system from $N-1$ to $N$, yields
\begin{eqnarray*}
J^{*}_{N-1} &=&\frac{1}{2}\b x[N-1]^{\top}(\b A-\b B\b K_{N-1})^{\top}\b H_N(\b A-\b B\b K_{N-1})\b x[N-1]\\
&{}&\quad +\frac{1}{2}\b x[N-1]^{\top}\left(\b K^{\top}_{N-1}\b R\b K_{N-1}\b + \b Q\right)\b x[N-1].
\end{eqnarray*}
Let
\[
\b H_{N-1}=(\b A-\b B\b K_{N-1})^{\top}\b H_N(\b A-\b B\b K_{N-1})+\b K^{\top}_{N-1}\b R\b K_{N-1}+\b Q,
\]
then
\[
J^{*}_{N-1}=\frac{1}{2}\b x[N-1]^{\top}\b H_{N-1}\b x[N-1].
\]
Decrementing $k$ to $N-2$ yields
\[
J_{N-2}=\frac{1}{2}\b x[N-1]^{\top}\b H_{N-1}\b x[N-1]+\frac{1}{2}\b x[N-2]^{\top}\b Q\b x[N-2]+\frac{1}{2}\b u[N-2]^{\top}\b R\b u[N-2].
\]
Note that $J_{N-2}$ has the same form as $J_{N-1}$. Thus, we obtain optimal feedback gain where $N$ is being replaced with  $N-1$. Continuing in this fashion, we get the following results for each $k=N-1,N-2,\ldots,0$:
\begin{eqnarray*}
\b K_k &=& \left(\b R + \b B^{\top}\b H_{k+1}\b B\right)^{-1}\b B^{\top}\b H_{k+1}\b A,\\
\b u^{*}[k] &=& -\b K_k\b x[k],\\
\b H_{k}&=&(\b A-\b B\b K_{k})^{\top}\b H_{k+1}(\b A-\b B\b K_{k})+\b K^{\top}_{k}\b R\b K_{k}+\b Q,
\end{eqnarray*}
and
\[
J^{*}_k=\frac{1}{2}\b x[k]^{\top}\b H_k\b x[k].
\]
The above control scheme can be implemented by computing the sequence of gain matrices $\{\b K_k\}$ off-line and stored. Then, we can implement the controller $\b u^{*}[k] = -\b K_k\b x[k]$. Note that the time-varying nature of the feedback controller may make an implementation of the control strategy\index{control strategy!implementation of}\index{implementation of control strategy} quite complicated, especially for large $n$ and $N$. However, as we go to the limit as $N\to \infty$, we can simplify the solution to the problem. Specifically, let us
consider a (reachable) system model
\[
\b x[k+1]=\b A\b x[k]+\b B\b u[k],\quad k \ge 0,\quad \b x(0)=\b x_0,
\]
and the performance index
\[
J=\frac{1}{2}\sum_{k=0}^{\infty}\left\{\b x[k]^{\top}\b Q\b x[k]+\b u[k]^{\top}\b R\b u[k]\right\}
\]
obtained from the previously considered performance index by letting $N\to \infty$ and setting $\b H_N=\b 0$. Then, the optimal controller takes the form
\[
\b u^{*}[k] = -\b K_{\infty}\b x[k],\quad k\ge0,
\]
where
\begin{equation}
\b K_{\infty} = \left(\b R + \b B^{\top}\b H_{\infty}\b B\right)^{-1}\b B^{\top}\b H_{\infty}\b A
\label{eq:kinfgafee}
\end{equation}
and $\b H_{\infty}$ is defined as the limit obtained by continuing the recursion
\begin{equation}
\b H_{k}=(\b A-\b B\b K_{k})^{\top}\b H_{k+1}(\b A-\b B\b K_{k})+\b K^{\top}_{k}\b R\b K_{k}+\b Q
\label{eq:recurelfohk}
\end{equation}
back to $k=-\infty$, that is,
\[
\lim_{k\to -\infty}\b H_k=\b H_{\infty}.
\]
Alternatively, letting $k\to -\infty$ in~(\ref{eq:recurelfohk}), we obtain
\begin{equation}
\b H_{\infty}=(\b A-\b B\b K_{\infty})^{\top}\b H_{\infty}(\b A-\b B\b K_{\infty})+\b K^{\top}_{\infty}\b R\b K_{\infty}+\b Q,
\label{eq:hinfredist}
\end{equation}
where $\b K_{\infty}$ is given by~(\ref{eq:kinfgafee}). Substituting~(\ref{eq:kinfgafee}) into~(\ref{eq:hinfredist}) yields
\begin{eqnarray}
{}&{}&\!\!\!\b H_{\infty}\nonumber\\
&=& \!\!\!\left(\b A-\b B\left(\b R + \b B^{\top}\b H_{\infty}\b B\right)^{-1}\b B^{\top}\b H_{\infty}\b A\right)^{\top}\b H_{\infty}\left(\b A-\b B\left(\b R + \b B^{\top}\b H_{\infty}\b B\right)^{-1}\b B^{\top}\b H_{\infty}\b A\right)\nonumber\\
&{}& \!\!\!+\left(\left(\b R + \b B^{\top}\b H_{\infty}\b B\right)^{-1}\b B^{\top}\b H_{\infty}\b A\right)^{\top}\b R\left(\left(\b R + \b B^{\top}\b H_{\infty}\b B\right)^{-1}\b B^{\top}\b H_{\infty}\b A\right)+\b Q.
\label{eq:algeeqfohibt}
\end{eqnarray}
Thus, we can compute $\b H_{\infty}$ by solving the above algebraic equation. We note that the above formidable looking equation is just the discrete algebraic Riccati equation\index{discrete algebraic Riccati equation}. Indeed, setting $\b P=\b H_{\infty}$, using the notation,
\[
\b S=\b R + \b B^{\top}\b H_{\infty}\b B=\b R + \b B^{\top}\b P\b B,
\]
and performing simple manipulations, we obtain~(\ref{eq:dtlqrrieq}).

We next discuss a continuous form of dynamic programming\index{dynamic programming!continuous form of}. We begin our discussion by analyzing the minimum time regulator\index{minimum time regulator} problem. Then, we focus our attention on a more general optimal control problem.

\subsection{Continuous Minimum Time Regulator Problem}\label{subsub:con_min_time_r}
In this subsection, we consider a control system modeled by
\[
\dot{\b x}=\b f(t,\b x,\b u)
\]
with a given initial state $\b x(t_0)=\b x_0$. We restrict the controller $\b u(t)$ to be restrained by the compact set $U \subset \mb{\Bb R}^m$. For example, we may restrict the components $u_i$, $i=1,2,\ldots,m$  of the control vector $\b u$ to satisfy the amplitude constraints\index{amplitude constraints}
\[
\left\arrowvert u_i \right\arrowvert \le k_i,
\]
where $k_i$ are fixed real numbers. A piecewise continuous control\index{piecewise continuous control}\index{control!piecewise continuous} vector $\b u$ satisfying the above constraints is called admissible. The set of admissible controllers is denoted $U$. Our goal is to determine an admissible control\index{admissible control}\index{control!admissible} law $\b u^*$ that steers the system's trajectory from a given initial state $\b x(t_0)=\b x_0$ to a specified final state $\b x_f$ in the shortest possible time. This is equivalent to minimizing the performance index
\[
J=\int_{t_0}^{t_f}dt= t_f-t_0.
\]
We call such a controller {\em time optimal}\index{time optimal!controller}. We assume that for any $\b x\ne \b x_f$ there exists admissible, time optimal control transferring the system from $\b x$ to $\b x_f$ in minimum time. We denote by $J^{*}(t,\b x)=J^{*}(t,\b x(t))$ the time of the optimal transfer from $\b x(t)$ to the specified $\b x_f$. Observe that $J^{*}(t,\b x) = J^{*}(t,x_1,x_2\ldots,x_n)$, that is,
\[
J^{*}:\mb{\Bb R}^{n+1}\to \mb{\Bb R},
\]
which means that $J^{*}$ is a real valued function of $n+1$ variables. We assume that the function $J^{*}$ is continuous and has continuous partial derivatives everywhere except at the point $\b x_f$.

Now, let $\b x_0=\b x(t_0)\ne \b x_f$ be an arbitrary initial state in $\mb{\Bb R}^n$. Assume that we applied some admissible $\b u_0\in U$ on the time interval $[t_0,t]$. As a result of this control action, the system is transferred from $\b x(t_0)$ to $\b x(t)$. The time spent on the motion from $\b x(t_0)$ to $\b x(t)$ is $t-t_0$. We then apply an admissible optimal control that transfers the system in minimum time $J^{*}(t,\b x(t))$ from $\b x(t)$ to $\b x_f$. The above considerations are illustrated in Figure~\ref{fig_transfer}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{transfer.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/transfer.bmp}
%\hskip1.2in\special{bmp:c:/mybook/transfer.bmp x=4in y=3in}
\caption{Illustration of the derivation of the condition that a controller must satisfy when transferring a dynamical system from an initial state to the final state in an optimal manner.}
\label{fig_transfer}
\end{figure}
Denoting the optimal time of transferring\index{optimal time!of transferring} the system from $\b x(t_0)$ to $\b x_f$ by $J^{*}(t_0,\b x(t_0))$, we obtain
\[
(t-t_0)+J^{*}(t,\b x(t))\ge J^{*}(t_0,\b x(t_0)).
\]
Rearranging the terms of the above equation and dividing by  $(t-t_0)$ yields
\[
1+\frac{J^{*}(t,\b x(t))- J^{*}(t_0,\b x(t_0))}{t-t_0}\ge 0.
\]
Letting $t\to t_0$, and taking into account that $\dot{x}_i=f_i(t,\b x,\b u)$, $i=1,2,\ldots,n$, we obtain
\begin{eqnarray*}
1+\frac{\partial J^{*}(t,\b x)}{\partial t}+\sum_{i=1}^{n}\frac{\partial J^{*}(t,\b x)}{\partial x_i}\dot{x}_i &=& 1+\frac{\partial J^{*}(t,\b x)}{\partial t}+\sum_{i=1}^{n}\frac{\partial J^{*}(t,\b x)}{\partial x_i}f_i\\
&=& 1+\frac{\partial J^{*}}{\partial t}+\frac{\partial J^{*}}{\partial \b x}\b f\\
&\ge& 0.
\end{eqnarray*}
We conclude that for any $\b x_0 \ne \b x_f$ and for any $\b u\in U$,
\[
1+\frac{\partial J^{*}}{\partial t}+\frac{\partial J^{*}}{\partial \b x}\b f \ge 0.
\]
We next define the Hamiltonian function\index{Hamiltonian function} for the time optimal control problem to be
\[
H(t,\b x,\b u,\b p)=1+\b p^{\top}\b f(t,\b x,\b u),
\]
where
\[
\b p=\left[\begin{array}{cccc}
\frac{\partial J^*}{\partial x_1}&\frac{\partial J^*}{\partial x_2}&\ldots&\frac{\partial J^*}{\partial x_n}
\end{array}\right]^{\top}.
\]
Thus, we can say that for any $\b x_0\ne \b x_f$ and for any $\b u\in U$,
\[
\frac{\partial J^*(t,\b x)}{\partial t}+H(t,\b x,\b u) \ge 0.
\]

We will use the PO to show that
\[
\frac{\partial J^*(t,\b x)}{\partial t}+H(t,\b x,\b u) = 0
\]
for any optimal transfer, that is,
\[
\frac{\partial J^*(t,\b x)}{\partial t}+\min_{\b u\in U}H(t,\b x,\b u) = 0 \quad \mb{for any} \quad \b x\ne \b x_f.
\]
Indeed, suppose that $\b u(t)\in U$ is an optimal controller that transfers the system from the initial state $\b x(t_0) \ne \b x_f$ to the given final state $\b x_f$ in a time $t_f$. Thus, we have $J^{*}(t_0,\b x(t_0))=t_f-t_0$. The transfer from $\b x(t)$ to $\b x_f$ requires $(t_f-t)$ time units. See Figure~\ref{fig_hjb} for an illustration of the arguments.
\begin{figure}
\centering
\includegraphics[scale=0.6]{hjb.eps}
%\vspace{3in}
%\hskip1.1in\special{bmp:c:/mybook/hjb.bmp x=4in y=3in}
\caption{Illustration of the derivation of the Hamilton-Jacobi-Bellman equation.}
\label{fig_hjb}
\end{figure}
Hence,
\[
t_f-t=J^{*}(t_0,\b x(t_0))-(t-t_0),
\]
and it is not possible to go from $\b x(t)$ to $\b x_f$ in a time less than $J^{*}(t_0,\b x(t_0))-(t-t_0)$. For if such a more rapid motion were to exist, indicated by the dashed curve in Figure~\ref{fig_hjb}, then moving from $\b x_0$ to $\b x(t)$ during the time $t-t_0$, and then from $\b x(t)$ to $\b x_f$ during a time less than $J^{*}(t_0,\b x(t_0))-(t-t_0)$, we would accomplish the transfer from $\b x_0$ to $\b x_f$ during a time less than $J^{*}(t_0,\b x(t_0))$, which is impossible. Therefore, we can write
\begin{eqnarray*}
t_f-t &=& J^{*}(t,\b x(t))\\
&=& J^{*}(t_0,\b x(t_0))-(t-t_0).
\end{eqnarray*}
Rearranging terms, dividing by $t-t_0$, and letting $t\to t_0$, yields
\[
\frac{\partial J^*(t,\b x)}{\partial t}+\min_{\b u\in U}H(t,\b x,\b u) = 0 \quad \mb{for any} \quad \b x\ne \b x_f.
\]
We summarize our considerations in the following theorem.
\begin{thm}
Let $\dot{\b x}=\b f(t,\b x,\b u)$ be a control system model, $U$ a set of admissible controllers, and $\b x_f$ a prescribed final state. Suppose that
\begin{itemize}
\item[(i)] for any $\b x\ne \b x_f$ there exists an admissible optimal control transferring\index{optimal control!transferring} $\b x$ to $\b x_f$,
\item[(ii)]  the function $J^{*}(t,\b x(t))$, which is the time of optimal transfer\index{optimal transfer} of the system from $\b x$, at time $t$,  to $\b x_f$, is continuously differentiable except at the point $\b x_f$,
\end{itemize}
then
\begin{eqnarray*}
\frac{\partial J^*(t,\b x)}{\partial t}+H(t,\b x,\b u) &\ge& 0 \quad \mb{for any}\quad \b u\in U \quad \mb{and}\quad \b x\ne \b x_f,\\
\frac{\partial J^*(t,\b x)}{\partial t}+\min_{\b u\in U}H(t,\b x,\b u) &=& 0 \quad \mb{for any optimal transfer, where}\quad \b x\ne \b x_f.
\end{eqnarray*}
\end{thm}
The above theorem is a special case of a more general result known as the continuous form of the dynamic programming\index{dynamic programming!continuous form of}. We now use the PO to derive a continuous form of dynamic programming for a class of problems with more general performance indices.

\subsection{The Hamilton-Jacobi-Bellman (HJB) Equation}

We consider a control process model of the form
\begin{equation}
\dot{\b x}=\b f(t,\b x,\b u),\quad \b x(t_0)=\b x_0,
\label{eq:sysmofpoamp}
\end{equation}
and we consider the performance index that we wish to minimize
\begin{equation}
J(t_0,\b x(t_0),\b u)=\Phi(t_f,\b x(t_f))+\int_{t_0}^{t_f}F(\tau,\b x(\tau),\b u(\tau))d\tau.
\label{eq:perfinfgenop}
\end{equation}
Define
\[
J(t,\b x(t),\b u(\tau))=\Phi(t_f,\b x(t_f))+\int_{t}^{t_f}F(\tau,\b x(\tau),\b u(\tau))d\tau
\]
with $t\le \tau\le t_f$, and let
\[
J^{*}(t,\b x(t))=\min_{\b u}J(t,\b x(t),\b u(\tau)).
\]
We subdivide the interval $[t,t_f]$ as
\[
[t,t_f]=[t,t+\Delta t]\cup[t+\Delta t,t_f].
\]
Then, we write
\[
J^{*}(t,\b x(t),\b u(\tau))=\min_{\b u}\left\{\int_{t}^{t+\Delta t}Fd\tau+\int_{t+\Delta t}^{t_f}Fd\tau+\Phi(t_f,\b x(t_f))\right\}.
\]
By the PO the trajectory on the interval $[t+\Delta t, t_f]$ must be optimal. Hence,
\[
J^{*}(t,\b x(t))=\min_{\b u}\left\{\int_{t}^{t+\Delta t}Fd\tau + J^{*}(t+\Delta t,\b x(t+\Delta t))\right\}.
\]
We next expand $J^{*}(t+\Delta t,\b x(t+\Delta t))$ in a Taylor series about the point $(t,\b x(t))$ to get
\[
J^{*}(t,\b x(t))=\min_{\b u}\left\{\int_{t}^{t+\Delta t}Fd\tau + J^*+\frac{\partial J^*}{\partial t}\Delta t+\frac{\partial J^*}{\partial \b x}\left(\b x(t+\Delta t)-\b x(t)\right)+ \mb{H.O.T.}\right\},
\]
where H.O.T. stands for higher order terms. Because $J^*(t,\b x(t))$ is independent of $\b u$, we cancel the terms $J^*(t,\b x(t))$ out and use the fact that $\b x(t+\Delta t)-\b x(t) \approx \dot{\b x}\Delta t$ to obtain
\[
0=\min_{\b u}\left\{\int_{t}^{t+\Delta t}Fd\tau + \frac{\partial J^*}{\partial t}\Delta t+\frac{\partial J^*}{\partial \b x}\dot{\b x}\Delta t+ \mb{H.O.T.}\right\}.
\]
By assumption $\Delta t$ is small, hence we can write
\[
0=\min_{\b u}\left\{F\Delta t+\frac{\partial J^*}{\partial t}\Delta t+\frac{\partial J^*}{\partial \b x}\b f\Delta t + \mb{H.O.T.}\right\}.
\]
Dividing the above equation by $\Delta t$ and letting $\Delta t \to 0$ yields
\begin{equation}
0=\frac{\partial J^*}{\partial t}+\min_{\b u}\left\{F+\frac{\partial J^*}{\partial \b x}\b f\right\}.
\label{eq:hjbbefh}
\end{equation}
Let $H=F+\frac{\partial J^*}{\partial \b x}\b f$ denote the Hamiltonian function\index{Hamiltonian function}. Then, we write the partial differential equation~(\ref{eq:hjbbefh}) as
\begin{equation}
\fbox{$ \displaystyle 0=\frac{\partial J^*}{\partial t}+\min_{\b u}H $}
\label{eq:hjbcoeqder}
\end{equation}
subject to the boundary condition\index{boundary condition}
\[
J^*(t_f,\b x(t_f))=\Phi(t_f,\b x(t_f)).
\]
The partial differential equation~(\ref{eq:hjbcoeqder}) for the optimal cost $J^*(t,\b x(t))$  is called the {\em Hamilton-Jacobi-Bellman\/} (HJB) equation\index{Hamilton-Jacobi-Bellman(HJB) equation}\index{HJB|see{Hamilton-Jacobi-Bellman(HJB)}}. It provides the solution to the optimal control problem for general nonlinear dynamical systems. However, analytical solution to the HJB  equation\index{HJB  equation} is difficult to obtain in most cases. We now illustrate, with three examples, how the HJB equation can be used to construct optimal controllers.

\begin{example}\rm
We consider a general class of dynamical systems modeled by the scalar differential equation
\[
\dot{x}=ax+bu,\quad x(t_0)=x_0,
\]
with the associated cost index
\[
J(t_0)=\frac{1}{2}fx(t_f)^2+\frac{1}{2}\int_{t_0}^{t_f}\left(qx(t)^2+ru(t)^2\right)dt,
\]
where the initial time $t_0=0$, the final time $t_f<\infty$ is fixed, and the final state $x(t_f)$ is free. Our goal is to construct the control $u^{*}$ that minimizes the cost index $J$.

We first form the Hamiltonian function for the problem,
\[
H\left(t,x,u,\frac{\partial J^{*}}{\partial x}\right)=\frac{1}{2}qx^2+\frac{1}{2}ru^2+\frac{\partial J^{*}}{\partial x}(ax+bu).
\]
Because there are no constraints on the control $u$, we can determine the form of the optimal control by applying the first order necessary condition for static optimization. We compute
\[
\frac{\partial H}{\partial u} = ru+b\frac{\partial J^{*}}{\partial x}= 0.
\]
The optimal control, therefore, has the form
\[
u^{*}=-\frac{1}{r}b\frac{\partial J^{*}}{\partial x}.
\]
This control indeed minimizes the Hamiltonian function $H$ because
\[
\frac{\partial ^2 H}{\partial u^2}= r > 0,
\]
that is, the second order sufficient condition for static minimization is satisfied. We thus reduced the problem of constructing the optimal control to that of finding $\frac{\partial J^{*}}{\partial x}$. Substituting $u^{*}$ into the HJB equation yields
\begin{eqnarray*}
&{}&\frac{\partial J^{*}}{\partial t}+\frac{1}{2}qx^2+\frac{b^2}{2r}\left(\frac{\partial J^{*}}{\partial x} \right)^2+ax\frac{\partial J^{*}}{\partial x}-\frac{b^2}{r}\left( \frac{\partial J^{*}}{\partial x}\right)^2\\
&{}&\quad = \frac{\partial J^{*}}{\partial t}+\frac{1}{2}qx^2-\frac{b^2}{2r} \left(\frac{\partial J^{*}}{\partial x} \right)^2+ax \frac{\partial J^{*}}{\partial x}\\
&{}&\quad = 0.
\end{eqnarray*}
Let us assume that $J^{*}$ is a quadratic function of the state for all $t\le t_f$, that is,
\[
J^{*}=\frac{1}{2}p(t)x^2,
\]
where $p(t)$ is a scalar function of time to be determined. Partial derivatives of $J^{*}$ are
\[
\frac{\partial J^{*}}{\partial t}= \frac{1}{2}\dot{p}(t)x^2\quad\mb{and}\quad \frac{\partial J^{*}}{\partial x}=p(t)x.
\]
Substituting the above into the HJB equation and performing some manipulations gives
\[
\left(\frac{1}{2}\dot{p}(t)+\frac{1}{2}q+ap(t)-\frac{b^2p(t)^2}{2r}\right)x^2=0.
\]
The above holds for all $x=x(t)$ if and only if
\[
\frac{1}{2}\dot{p}(t)+\frac{1}{2}q+ap(t)-\frac{b^2p(t)^2}{2r}=0,
\]
or equivalently
\begin{equation}
\dot{p}(t)+q+2ap(t)-\frac{b^2p(t)^2}{r}=0
\label{eq:eqscaldae}
\end{equation}
subject to the boundary condition $p(t_f)=f$ that we obtain by comparing $J^{*}=\frac{1}{2}p(t)x^2$ and $J(t_f)=\frac{1}{2}fx(t_f)^2$. To solve the above equation we assume that
\[
p(t)=\rho\frac{\dot{w}(t)}{w(t)},
\]
where the parameter $\rho$ is to be determined. Note that
\[
\dot{p}(t)=\rho \frac{\ddot{w}(t)w(t)-\dot{w}(t)^2}{w(t)^2}.
\]
Substituting the above into~(\ref{eq:eqscaldae}), we obtain
\begin{eqnarray*}
&{}&\dot{p}+q+2ap-\frac{b^2p^2}{r}\\
&{}& \quad =\rho \frac{\ddot{w}w-\dot{w}^2}{w^2}+q\frac{w^2}{w^2}+2a\rho\frac{\dot{w}}{w}-\frac{b^2\rho ^2\dot{w}^2/r}{w^2}\\
&{}&\quad =\frac{\rho \ddot{w}w-\rho \dot{w}^2+qw^2+2a\rho \dot{w}w-b^2\rho ^2\dot{w}^2/r }{w^2}\\
&{}&\quad =0.
\end{eqnarray*}
We choose $\rho$ so that the terms nonlinear in $\dot{w}$ cancel out. For this it is enough that
\[
-\rho \dot{w}^2-b^2\rho ^2\dot{w}^2/r=0.
\]
Hence, if we set
\[
\rho=-\frac{r}{b^2},
\]
then the equation to be solved takes the form
\[
\frac{\left(-\frac{r}{b^2}\ddot{w}-\frac{2ar}{b^2}\dot{w}+qw \right)w}{w^2}=0.
\]
The above is satisfied if
\[
-\frac{r}{b^2}\ddot{w}-\frac{2ar}{b^2}\dot{w}+qw=0,
\]
or equivalently
\[
\ddot{w}+2a\dot{w}-\frac{qb^2}{r}w=0,
\]
which is just a second order linear differential equation. Its characteristic equation\index{characteristic equation} has two real roots:
\[
s_{1,2}=-a\pm \sqrt{a^2+\frac{qb^2}{r}}.
\]
Hence
\[
w(t)=C_1e^{s_1t}+C_2e^{s_2t},
\]
and
\[
p(t)=\rho\frac{\dot{w}(t)}{w(t)}= \rho \frac{s_1C_1e^{s_1t}+s_2C_2e^{s_2t}}{C_1e^{s_1t}+C_2e^{s_2t}}.
\]
Using the boundary value $p(t_f)=f$, we express $C_1$ in terms of $C_2$ to obtain
\[
C_1=\frac{\rho s_2 C_2 e^{s_2 t_f}-fC_2e^{s_2 t_f}}{fe^{s_1 t_f}-\rho s_1e^{s_1t_f}}.
\]
We substitute the above expression for $C_1$ into $p(t)$ and cancel out $C_2$. Then, we use the obtained result to construct the optimal state-feedback controller\index{optimal state-feedback}\index{state-feedback!optimal}
\[
\fbox{$\displaystyle u^{*} = -\frac{1}{r}b\frac{\partial J^{*}}{\partial x}= -\frac{1}{r}bp(t)x $}
\]
\label{ex:soldifriceq}
\end{example}
\markex
\begin{example}\label{ex:dcmotdynpro}\rm
Consider the schematic of an armature controlled DC motor\index{DC motor!armature controlled}\index{armature controlled!DC motor}, given in Figure~\ref{fig:fins7dc}, where the system parameters are: $R_a=2\,\Omega$, $L_a\approx 0\,\mb{H}$, $K_b=2\,\mb{V/rad/sec}$, $K_i=2\,\mb{Nm/A}$, and the equivalent moment of inertia\index{equivalent moment of inertia}\index{moment of inertia!equivalent} referred to the motor shaft\index{referred to motor shaft!moment of inertia} is $I_{eq}=1$~kg$\cdot$m$^2$. The friction is assumed to be negligible.
\begin{figure}
\centering
\includegraphics[scale=0.6]{fins7dc.eps}
%\vspace{3.5in}
%\hskip1in\special{bmp:c:/mybook/fins7dc.bmp x=4.67in y=3.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/fins7dc.bmp}
\caption{Schematic of an armature-controlled DC motor system of Example~\ref{ex:dcmotdynpro}.}
\label{fig:fins7dc}
\end{figure}
\begin{itemize}
\item[(i)] Construct the first-order differential equation modeling the DC motor;
\item[(ii)] Use the Hamilton-Jacobi-Bellman equation\index{Hamilton-Jacobi-Bellman equation} to find the optimal state feedback controller, $e_a=k(t)\omega$, that minimizes the performance index
\[
J=\frac{1}{2}\omega(10)^2 + \int_0^{10}R_a i_a(t)^2dt,
\]
where $i_a$ is the armature current. The final state is free. There are no constraints on $e_a$. Assume $J^{*}=\frac{1}{2}p(t)\omega^2$.
\end{itemize}

%\markdefn

Applying Kirchhoff's voltage law\index{Kirchhoff's voltage law} to the armature circuit\index{armature circuit} yields
\[
R_a i_a+K_b\omega = e_a.
\]
The torque developed by the motor, $T_m$, is
\[
I_{eq}\dot{\omega}=T_m=K_i i_a.
\]
Substituting $i_a$ from the first equation into the second and dividing both sides by $I_{eq}$ yields
\[
\dot{\omega}=-\frac{K_iK_b}{R_aI_{eq}}\omega +\frac{K_i}{I_{eq}R_a}e_a.
\]
Substituting into the above the parameter values, we obtain
\begin{equation}
\dot{\omega}=-2\omega +e_a.
\label{eq:dcexsimf7hjb}
\end{equation}
We represent the performance index, $J$, in terms of $\omega$ and $e_a$. Applying Ohm's law to the armature circuit, we get
\[
e_a-K_b\omega=R_a i_a.
\]
Hence,
\[
J=\frac{1}{2}\omega(10)^2 + \int_0^{10}\frac{\left(e_a-K_b \omega \right)^2}{R_a}dt.
\]
The Hamiltonian function\index{Hamiltonian function} is
\[
H=\frac{\left(e_a-K_b \omega \right)^2}{R_a}+\frac{\partial J^{*}}{\partial \omega}\left( -\frac{K_iK_b}{R_aI_{eq}}\omega +\frac{K_i}{I_{eq}R_a}e_a\right)=\frac{\left( e_a-2\omega\right)^2}{2}+\frac{\partial J^{*}}{\partial \omega}\left(-2\omega+e_a \right).
\]
Because there are no constraints on $e_a$, we can find the optimal control\index{optimal control} by solving the equation
\[
0=\frac{\partial H}{\partial e_a}=\frac{2\left( e_a-K_b\omega\right)}{R_a}+\frac{K_i}{I_{eq}R_a}\frac{\partial J^{*}}{\partial \omega}.
\]
Hence,
\begin{equation}
e_a^*=2\omega-\frac{\partial J^{*}}{\partial \omega}.
\label{eq:confin7hjb}
\end{equation}
Substituting~(\ref{eq:confin7hjb}) into the HJB equation\index{HJB equation} yields
\[
0=\frac{\partial J^{*}}{\partial t}+\frac{1}{2}\left(-\frac{\partial J^{*}}{\partial \omega} \right)^2+\frac{\partial J^{*}}{\partial \omega}\left( -\frac{\partial J^{*}}{\partial \omega}\right)=\frac{\partial J^{*}}{\partial t}-\frac{1}{2}\left(\frac{\partial J^{*}}{\partial \omega} \right)^2.
\]
By assumption $J^{*}=\frac{1}{2}p(t)\omega^2$. Hence,
\[
\frac{\partial J^{*}}{\partial t}=\frac{1}{2}\dot{p}\omega^2\quad\mb{and}\quad \frac{\partial J^{*}}{\partial \omega}=p\omega.
\]
Substituting the above into the HJB equation, we get
\[
\frac{1}{2}\left(\dot{p}-p^2 \right)\omega^2=0.
\]
We next solve the nonlinear differential equation\index{nonlinear differential equation!HJB equation}
\begin{equation}
\dot{p}-p^2 =0.
\label{eq:fseqforp}
\end{equation}
We assume that $p=\rho\frac{\dot{w}}{w}$. Therefore,
\[
\dot{p}=\rho\frac{\ddot{w}w-\dot{w}^2}{w^2}.
\]
We substitute the expressions for $p$ and $\dot{p}$ into~(\ref{eq:fseqforp}) to get
\[
\frac{\rho \ddot{w}w-\rho \dot{w}^2-\rho^2\dot{w}^2}{w^2}=0.
\]
Selecting $\rho=-1$ eliminates the nonlinear terms in the numerator and we are left with the equation
\[
\ddot{w}=0,
\]
whose solution is
\[
w=w(t)=C_1t+C_2,
\]
where $C_1$ and $C_2$ are integration constants\index{integration constants!eliminate}. Hence,
\[
p=-\frac{\dot{w}}{w}=-\frac{C_1}{C_1t+C_2}.
\]
Because
\[
J^{*}(10)=\frac{1}{2}\omega(10)^2=\frac{1}{2}p(10)\omega(10)^2,
\]
we conclude that $p(10)=1$. We use this information to eliminate one of the integration constants. We obtain
\[
p(10)=1=-\frac{C_1}{10 C_1+C_2}.
\]
Hence,
\[
C_2=-11C_1,
\]
and
\[
p=p(t)=-\frac{C_1}{C_1 t+C_2}=-\frac{C_1}{C_1 t-11C_1}=-\frac{1}{t-11}.
\]
Therefore,
\[
\fbox{$\displaystyle e_a^{*}=2\omega-\frac{\partial J^{*}}{\partial \omega}=2\omega-p\omega=\left(2+\frac{1}{t-11} \right)\omega$}
\]
\end{example}
\markex
In the above examples, the final time $t_f$ was fixed and $t_f<\infty$. In the following example, we consider the linear quadratic regulator design for $t_f=\infty$ using the HJB equation.

\begin{example}\rm
For the dynamical system model,
\[
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array}\right]=\left[\begin{array}{cc}
0 & 1\\
0 & -1
\end{array}\right]\left[\begin{array}{c}
x_1\\
x_2
\end{array}\right]+\left[\begin{array}{c}
0\\
1
\end{array}\right]u,
\]
and the performance index
\[
J=\frac{1}{2}\int_0^{\infty}\left(x^2_1+x^2_2+u^2\right)dt,
\]
we will use the HJB equation to find the optimal state-feedback
controller. We first form the Hamiltonian function
\[
H=\frac{1}{2}\left(x^2_1+x^2_2+u^2\right)+\frac{\partial J^*}{\partial x_1}x_2+\frac{\partial J^*}{\partial x_2}(-x_2+u).
\]
We assume that $J^*$ has the form
\[
J^*=\frac{1}{2}\left(\alpha x^2_1+2\beta x_1 x_2 + \gamma x^2_2\right).
\]
Hence,
\[
\frac{\partial J^*}{\partial t}=0,
\]
and the HJB equation takes the form
\[
\min_uH=0,
\]
with the boundary condition\index{HJB equation!with boundary condition}
\[
J^*(\b x(t_f)) = J^*(\b x(\infty))= 0.
\]
Because there are no constraints on $u$, we apply the first-order necessary condition for unconstrained optimization to find $u$ minimizing $H$,\index{minimizing the Hamiltonian}
\[
\frac{\partial H}{\partial u} = 0.
\]
We denote by $u^*$ the solution to the above equation. Differentiating, we obtain
\[
u^*+\frac{\partial J^*}{\partial x_2}=0.
\]
Hence,
\[
u^*=-\frac{\partial J^*}{\partial x_2}.
\]
We substitute $u^*$ into the HJB equation to get
\[
\frac{1}{2}\left(x^2_1+x^2_2+\left( \frac{\partial J^*}{\partial x_2}\right)^2\right)+\frac{\partial J^*}{\partial x_1}x_2-\frac{\partial J^*}{\partial x_2}x_2-\left( \frac{\partial J^*}{\partial x_2}\right)^2=0.
\]
Note that
\[
\frac{\partial J^*}{\partial x_1} = \alpha x_1 + \beta x_2\quad\mb{and}\quad
\frac{\partial J^*}{\partial x_2} = \beta x_1 + \gamma x_2.
\]
Substituting the above into the HJB equation yields
\[
\frac{1}{2}\left(x^2_1+x^2_2+(\beta x_1+\gamma x_2)^2 \right) + (\alpha x_1 + \beta x_2)x_2 - (\beta x_1 + \gamma x_2)x_2-(\beta x_1 + \gamma x_2)^2=0.
\]
We rearrange the above equation to obtain
\[
\frac{1}{2}x^2_1(1-\beta ^2)+x_1x_2(\alpha -\beta -\beta \gamma)+\frac{1}{2}x^2_2(1-2\gamma -\gamma^2 + 2\beta )=0.
\]
For this equation to hold, we have to have
\[
\left.\begin{array}{lll}
1-\beta ^2 &=& 0\\
\alpha -\beta -\beta \gamma &=& 0\\
1-2\gamma -\gamma^2 + 2\beta &=& 0.
\end{array}\right\}
\]
There are three solutions to the above system of equations:
\begin{enumerate}
\item $\alpha=2$, $\beta = \gamma =1$;
\item $\alpha =-2$, $\beta =1$, $\gamma=-3$;
\item $\alpha =0$, $\beta = \gamma = -1$.
\end{enumerate}
Thus, we have three candidate solutions to the problem:
\begin{eqnarray*}
u^{(1)} &=& -x_1-x_2,\\
u^{(2)} &=& -x_1+3x_2,\\
u^{(3)} &=& x_1+x_2.
\end{eqnarray*}
However, only $u^{(1)}$ stabilizes the plant. Hence, the optimal control law for the problem is
\[
u^*=-x_1-x_2.
\]
Note, that only for this controller the boundary condition, $J^*(\b x(t_f)) = J^*(\b x(\infty))= 0$, is satisfied\index{boundary condition!satisfied}.
\end{example}
\markex

\section{Pontryagin's Minimum Principle}
\label{sec:pontrminpr}
\subsection{Optimal Control With Constraints on Inputs}
In this section, we provide a derivation of the minimum principle of Pontryagin\index{minimum principle of Pontryagin}, which is a generalization of the Euler-Lagrange equations\index{Pontryagin's minimum principle!generalization of Euler-Lagrange equations} that also includes problems with constraints on the control\index{Pontryagin's minimum principle!includes constraints on control} inputs. Only a special case of the minimum principle is stated. However, this special case covers a large class of control problems. We consider the problem of minimizing the performance index given by
\begin{equation}
J=\Phi(\b x(t_f))+\int_{t_0}^{t_f}F(\b x(t),\b u(t))dt
\label{eq:tiinvpeinmp}
\end{equation}
subject to
\begin{equation}
\dot{\b x}=\b f(\b x,\b u),\quad \b x(t_0)=\b x_0,\quad\mb{and}\quad \b x(t_f)=\b x_f.
\label{eq:tiinplfomp}
\end{equation}
We consider two cases: fixed final state and free final state.\index{Pontryagin's minimum principle!fixed final state}\index{Pontryagin's minimum principle!free final state}

%It follows from our discussion of the HJB equation that the optimal controller, $\b u^{*}$, minimizes the Hamiltonian function, that is,
%\[
%\b u^{*}=\arg\min_{\b u\in U}H(\b x,\b u,\b p),
%\]
To proceed we define the Hamiltonian function $H(\b x,\b u,\b p)$ as
\[
H(\b x,\b u,\b p)=H=F+\b p^{\top}\b f,
\]
where the costate vector $\b p$ will be determined in the analysis to follow.
%and
%\[
%\b p=\left[\begin{array}{cccc}
%p_1 & p_2& \cdots & p_n
%\end{array}\right]^{\top}=
%\left[\begin{array}{cccc}
%\frac{\partial J^*}{\partial x_1}&\frac{\partial J^*}{\partial x_2}&\cdots&\frac{\partial J^*}{\partial x_n}\end{array}\right]^{\top}.
%\]
%To proceed, we need one more assumption in addition to the assumptions we made while discussing the HJB equation. We assume that $J^*(\b x)$ has continuous second partial derivatives $\frac{\partial ^2 J^*}{\partial x_i \partial x_j}$, that is, the functions $p_i$ have continuous first derivatives $\frac{\partial p_i}{\partial x_j}$. This implies that
%\begin{equation}
%\frac{\partial \b p}{\partial \b x}=\left( \frac{\partial \b p}{\partial \b x}\right)^{\top}.
%\label{eq:consofeqmder}
%\end{equation}
%We also assume that the components of $\b f$ have continuous first derivatives $\frac{\partial f_i}{\partial x_j}$.
%We then compute
%\begin{eqnarray}
%\frac{\partial H}{\partial \b x} &=& \frac{\partial}{\partial \b x}\left(F+\b p^{\top} \b f\right)\nonumber\\
%&=& \frac{\partial F}{\partial \b x}+\b f^{\top}\frac{\partial \b p}{\partial \b x}+\b p^{\top}\frac{\partial \b f}{\partial \b x},
%\label{eq:derhwxwf}
%\end{eqnarray}
%where $\frac{\partial F}{\partial \b x}=\nabla_{\b x}^{\top} F$. From our previous discussion it follows that the Hamiltonian function attains its minimum on the optimal trajectory. Taking this into account, we conclude that for the optimal process
%\[
%\frac{\partial H}{\partial \b x}=\b 0 \quad \mb{for} \quad (\b x,\b u)=(\b x^*,\b u^*).
%\]
%Hence, for the optimal process $(\b x^*,\b u^*)$, we have
%\begin{equation}
%\frac{\partial F}{\partial \b x}+\b f^{\top}\frac{\partial \b p}{\partial \b x}+\b p^{\top}\frac{\partial \b f}{\partial \b x}=\b 0.
%\label{eq:derhwx}
%\end{equation}
%We now differentiate $\b p$ with respect to time to get
%\begin{equation}
%\dot{\b p} = \frac{\partial \b p}{\partial \b x}\dot{\b x}= \frac{\partial \b p}{\partial \b x}\b f.
%\label{eq:pdotwt}
%\end{equation}
%Combining~(\ref{eq:derhwx}) and~(\ref{eq:pdotwt}) and taking into account~(\ref{eq:consofeqmder}), we obtain
%\begin{equation}
%\dot{\b p}=-\left(\frac{\partial F}{\partial \b x} \right)^{\top}-\left(\frac{\partial \b f}{\partial \b x}\right)^{\top}\b p.
%\label{eq:dieqfopindox}
%\end{equation}
%Thus, $\b p$ is the solution to the differential equation~(\ref{eq:dieqfopindox}), and hence $\b p=\b p(t)$. Therefore $\frac{\partial \b p}{\partial \b x}=\b 0$, and by~(\ref{eq:derhwxwf}), we have
%\begin{equation}
%\frac{\partial H}{\partial \b x}=\frac{\partial F}{\partial \b x}+\b p^{\top}\frac{\partial \b f}{\partial \b x}.
%\label{eq:frahwrtbx}
%\end{equation}
%Comparing~(\ref{eq:dieqfopindox}) and(\ref{eq:frahwrtbx}) gives
%\begin{equation}
%\dot{\b p}=-\left( \frac{\partial H}{\partial \b x}\right)^{\top}
%\label{eq:cooradjeqmp}
%\end{equation}
%The above equation is called in the literature the {\em adjoint\/} or {\em costate\/} equation.

We discuss the case when the final time $t_f$ is fixed. We follow the development of Luenberger~\cite[Section~11.1]{luenberger79}. We ``adjoin'' to $J$ additional terms that sum up to zero. We note that because the state trajectories must satisfy the equation, $\dot{\b x}=\b f(\b x,\b u)$, we have
\begin{equation}
\b f(\b x,\b u)-\dot{\b x}=\b 0.
\label{ex:diff_eqn_zero}
\end{equation}
We introduce the modified objective performance index,
\begin{equation}
\tilde{J}=J+\int_{t_0}^{t_f} \b p(t)^{\top}(\b f(\b x,\b u)-\dot{\b x})\, dt .
\label{eq:modified_J}
\end{equation}
Because for any state trajectory,~\eqref{ex:diff_eqn_zero} is satisfied, for any choice of $\b p(t)$ the value of $\tilde{J}$ is the same as that of $J$. Therefore, we can express $\tilde{J}$ as
\begin{eqnarray}
\tilde{J} &=& \Phi(\b x(t_f))+\int_{t_0}^{t_f}F(\b x(t),\b u(t))\,dt+\int_{t_0}^{t_f} \b p(t)^{\top}(\b f(\b x,\b u)-\dot{\b x})\, dt\nonumber\\
&=& \Phi(\b x(t_f))+\int_{t_0}^{t_f} (H(\b x, \b u, \b p)-\b p^{\top}\dot{\b x})\,dt .
\label{eq:jtideandH}
\end{eqnarray}
Let $\b u(t)$ be a nominal control strategy;\index{nominal control strategy} it determines a corresponding state trajectory $\b x(t)$. If we apply another control strategy, say $\b v(t)$, that is ``close'' to $\b u(t)$, then $\b v(t)$ will produce a state trajectory close to the nominal trajectory. This new state trajectory is just a perturbed version of $\b x(t)$ and it can be represented as
\[
\b x(t) + \delta \b x(t).
\]
The change in the state trajectory yields a corresponding change in the modified performance index\index{modified performance index!change in it}. We represent  this change as $\delta \tilde{J}$; it has the form,
\begin{equation}
\delta \tilde{J}=\Phi(\b x(t_f) + \delta \b x(t_f))-\Phi(\b x(t_f)) + \int_{t_0}^{t_f} (H(\b x+\delta \b x,\b v,\b p) - H(\b x,\b u,\b p) - \b p^{\top}\delta \dot{\b x}  )\,dt .
\label{eq:delta Jtilde}
\end{equation}
We evaluate the above integral. We first perform an integration by parts to obtain
\begin{equation}
\int_{t_0}^{t_f} \b p^{\top}\delta \dot{\b x}\,dt = \b p(t_f)^{\top}\delta \b x(t_f) - \b p(t_0)^{\top}\delta \b x(t_0) - \int_{t_0}^{t_f} \dot{\b p}^{\top}\delta \b x\,dt .
\label{eq:int_by_partlxd}
\end{equation}
Note that $\delta \b x(t_0)=\b 0$ because a change in the control strategy does not change the initial state. Taking into account the above, we represent~\eqref{eq:delta Jtilde} as
\begin{eqnarray*}
\delta \tilde{J} &=& \Phi(\b x(t_f) + \delta \b x(t_f))-\Phi(\b x(t_f)) -\b p(t_f)^{\top}\delta \b x(t_f)\\
&{}&\quad + \int_{t_0}^{t_f} (H(\b x+\delta \b x,\b v,\b p) - H(\b x,\b u,\b p) + \dot{\b p}^{\top}\delta \b x  )\,dt .
\label{eq:delta Jtilde_wpxdo}
\end{eqnarray*}
We next replace $(\Phi(\b x(t_f) + \delta \b x(t_f))-\Phi(\b x(t_f)))$  with its first-order approximation, and add and subtract the term $H(\b x,\b v,\b p)$ under the integral to obtain
\begin{eqnarray*}
\delta \tilde{J} & = & \left(\left.\nabla_{\b x}\Phi\right\arrowvert_{t=t_f} -\b p(t_f)\right)^{\top}\delta \b x(t_f)\\
&{}&\quad + \int_{t_0}^{t_f} (H(\b x+\delta \b x,\b v,\b p) -H(\b x,\b v,\b p) + H(\b x,\b v,\b p)- H(\b x,\b u,\b p) + \dot{\b p}^{\top}\delta \b x  )\,dt .
\label{eq:delta Jtilde_wpaproxphi}
\end{eqnarray*}
Replacing $(H(\b x+\delta \b x,\b v,\b p) - H(\b x,\b v,\b p))$ with its first-order approximation,
\[
H(\b x+\delta \b x,\b v,\b p) - H(\b x,\b v,\b p)=\frac{\partial H}{\partial \b x}\delta \b x,
\]
gives
\begin{eqnarray*}
\delta \tilde{J} &=& \left(\left.\nabla_{\b x}\Phi\right\arrowvert_{t=t_f} -\b p(t_f)\right)^{\top}\delta \b x(t_f)\\
&{}&\quad + \int_{t_0}^{t_f} \left(\left( \frac{\partial H}{\partial \b x} + \dot{\b p}^{\top}\right)\delta \b x + H(\b x,\b v,\b p)- H(\b x,\b u,\b p)  \right)\,dt .
\label{eq:delta Jtilde_wpxdopm}
\end{eqnarray*}
Selecting in the above $\b p$ as the solution to the differential equation
\[
\frac{\partial H}{\partial \b x} + \dot{\b p}^{\top}=\b 0^{\top}
\]
with the final condition
\[
\b p(t_f)=\left.\nabla_{\b x}\Phi\right\arrowvert_{t=t_f},
\]
reduces $\delta \tilde{J}$ to
\[
\delta \tilde{J}=\int_{t_0}^{t_f} \left(H(\b x,\b v,\b p)- H(\b x,\b u,\b p)\right)\, dt.
\]
If the original control $\b u(t)$ is optimal, then we should have $\delta \tilde{J}\ge 0$, that is, $H(\b x,\b v,\b p)\ge H(\b x,\b u,\b p)$. We summarize our development in the following theorem.
\begin{thm}
\label{thm:pontryagins}
Necessary conditions for $\b u\in U$ to minimize~(\ref{eq:tiinvpeinmp}) subject to~(\ref{eq:tiinplfomp}) are:
\[
\dot{\b p}=-\left( \frac{\partial H}{\partial \b x}\right)^{\top},
\]
where $H=H(\b x,\b u,\b p)=F(\b x,\b u)+\b p^{\top}\b f(\b x,\b u)$,
\[
H(\b x^*,\b u^*,\b p^*)=\min_{\b u\in U}H(\b x,\b u,\b p).
\]
If the final state, $\b x(t_f)$, is free, then in addition to the above conditions it is required that the following end-point condition\index{end-point condition} is satisfied,
\[
\b p(t_f)=\left.\nabla_{\b x}\Phi\right\arrowvert_{t=t_f}.
\]
\end{thm}

The equation,
\begin{equation}
\dot{\b p}=-\left( \frac{\partial H}{\partial \b x}\right)^{\top}
\label{eq:cooradjeqmp}
\end{equation}
is called in the literature the {\em adjoint\/} or {\em costate\/} equation.\index{costate equation}\index{adjoint equation|see{costate equation}}

We illustrate the above theorem with the following well-known example.
\begin{example}\label{ex:doubintboundco}\rm
We consider a controlled object modeled by
\begin{eqnarray*}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2
\end{array}\right]&=&\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\left[\begin{array}{c}
{x}_1\\
{x}_2
\end{array}\right]+\left[\begin{array}{c}
0\\
1
\end{array}\right]u\\
&=& \b A\b x + \b bu
\end{eqnarray*}
with the performance index
\[
J=\int_{0}^{t_f}dt.
\]
The control is required to satisfy
\[
|u(t)|\le 1
\]
for all $t\in[0,t_f]$. This constraint means that the control must have magnitude no greater than $1$. Our objective is to find admissible control that minimizes $J$ and transfers the system from a given initial state $\b x_0$ to the origin.

We begin by finding the Hamiltonian function\index{Hamiltonian function} for the problem,
\begin{eqnarray*}
H &=& 1+\b p^{\top}(\b A\b x + \b bu)\\
&=& 1 + \left[\begin{array}{cc}
p_1 & p_2
\end{array}\right]\left(\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\left[\begin{array}{c}
{x}_1\\
{x}_2
\end{array}\right]+\left[\begin{array}{c}
0\\
1
\end{array}\right]u \right)\\
&=& 1 + p_1x_2 + p_2 u.
\end{eqnarray*}
The costate equations are
\[
\left[\begin{array}{c}
\dot{p}_1\\
\dot{p}_2
\end{array}\right] = -\left[\begin{array}{c}
\frac{\partial H}{\partial x_1}\\
\frac{\partial H}{\partial x_2}
\end{array}\right]= \left[\begin{array}{c}
0\\
-p_1
\end{array}\right].
\]
Solving the costate equations\index{costate equations} yields
\[
p_1 = d_1\quad\mb{and}\quad p_2 = -d_1t+d_2,
\]
where $d_1$ and $d_2$ are integration constants. We now find an admissible control\index{admissible control} minimizing the Hamiltonian\index{minimizing the Hamiltonian},
\begin{eqnarray*}
\arg_{u}\min H &=& \arg_u \min (1 + p_1x_2+p_2u)\\
&=& \arg_u \min (p_2u)\\
&=& \left\{\begin{array}{ccc}
u(t)=1 & \mb{if}& p_2<0\\
u(t)=-1 & \mb{if}& p_2>0
\end{array}\right..
\end{eqnarray*}
Hence,
\[
u^*(t) = -\mb{sign}(p_2^*)= -\mb{sign}(-d_1t+d_2),
\]
where
\[
\mb{sign}(z)=\left\{\begin{array}{ccc}
1 & \mb{if}& z>0\\
%0 & \mb{if}& z=0\\
-1 & \mb{if}& z<0.
\end{array}\right.
\]
Thus, the optimal control\index{optimal control} law is piecewise constant taking the values $1$ or $-1$. This control law has at most two intervals of constancy because the argument of the sign function is a linear function, $-d_1t+d_2$, that changes its sign at most once. This type of control is called a {\em bang-bang control\/}\index{bang-bang control} because it switches back and forth between its extreme values. System trajectories for $u=1$ and $u=-1$ are families of parabolas and their equations were derived in Example~\ref{example_sat}. A phase portrait\index{phase portrait} for $u=1$ is shown in Figure~\ref{fig:bbol1}, while for $u=-1$ in Figure~\ref{fig:bbol2}. Note that only one parabola from each family passes through the specified terminal point $\b x=\left[\begin{array}{cc}
0& 0
\end{array}\right]^{\top}$ in the state plane\index{state plane}. Segments of the two parabolas through the origin form the switching curve,\index{switching curve}
\[
x_1=-\frac{1}{2}x^2_2\, \mb{sign}(x_2).
\]
This means that if an initial state is above the switching curve, then $u=-1$ is used until the switching curve is reached. Then, $u=1$ is used to reach the origin. For an initial state below the switching curve, the control $u=1$ is used first to reach the switching curve, and then the control is switched to $u=-1$. We implement the above control action as
\[
u=-\mb{sign}(v)=\left\{\begin{array}{ccc}
1&\mb{if}& v<0\\
%0 &\mb{if}&v=0\\
-1&\mb{if}& v>0,
\end{array}\right.
\]
where $v=v(x_1,x_2)=x_1+\frac{1}{2}x^2_2\mb{sign}(x_2)=0$ is the equation describing the switching curve. We can use the above equation to synthesize a closed-loop system such that starting at an arbitrary initial state in the state plane, the trajectory will always be moving in an optimal fashion towards the origin. Once the origin is reached, the trajectory will stay there.
%A block diagram of the closed-loop time optimal system is shown in %Figure~\ref{fig:bbcloop}.
%\begin{figure}
%\vspace{4.5in}
%\hskip.2in\special{bmp:c:/mybook/bbcloop.bmp x=6in y=4.5in}
%\caption{Time optimal closed-loop system with state feedback.}
%\label{fig:bbcloop}
%\end{figure}
A phase portrait of the closed-loop time optimal system is given in Figure~\ref{fig:bbclp}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{bbclp.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/mybook/bbclp.bmp}
%\hskip1in\special{bmp:c:/mybook/bbclp.bmp x=4in y=3in}
\caption{A phase-plane portrait of the time optimal closed-loop system of Example~\ref{ex:doubintboundco}.}
\label{fig:bbclp}
\end{figure}
\end{example}\markex

\begin{example}\rm
Consider the schematic of an armature-controlled DC motor\index{DC motor!armature-controlled}, given in Figure~\ref{fig:fins7dc}. Use the minimum principle\index{minimum principle} to find the armature voltage $e_a(t)$ such that the motor angular velocity $\omega$ changes from $\omega(0)=0\,\mb{rad/sec}$ at $t=0$ to $\omega(1)=10\,\mb{rad/sec}$ while minimizing the energy dissipated in the armature resistor\index{minimizing energy!dissipated in armature resistor}, $R_a$. There are no constraints on $e_a$. Note that the energy dissipated in the armature resistor is
\[
J=\int_0^{1}R_a i_a(t)^2dt,
\]
where $i_a$ is the armature current.

%\markdefn

The modeling equation of the DC motor was derived in Example~\ref{ex:dcmotdynpro}. It is given by~(\ref{eq:dcexsimf7hjb}) on page~\pageref{eq:dcexsimf7hjb}. We represent the expression for the energy dissipated in the armature resistor in terms of $\omega$ and $e_a$. Applying Ohm's law to the armature circuit, we get
\[
e_a-K_b\omega=R_a i_a.
\]
Hence,
\[
J=\int_0^1R_ai_a(t)^2dt=\int_0^1\frac{\left(e_a-K_b \omega \right)^2}{R_a}dt.
\]
The Hamiltonian function\index{Hamiltonian function!DC motor} is
\[
H=\frac{\left(e_a-K_b \omega \right)^2}{R_a}+p\left( -\frac{K_iK_b}{R_aI_{eq}}\omega +\frac{K_i}{I_{eq}R_a}e_a\right).
\]
Because there are no constraints on $e_a$, we can find the optimal control by solving the equation
\[
0=\frac{\partial H}{\partial e_a}=\frac{2\left( e_a-K_b\omega\right)}{R_a}+\frac{K_i}{I_{eq}R_a}p.
\]
Hence,
\begin{equation}
e_a=2\omega-p.
\label{eq:confin7}
\end{equation}
The costate equation\index{costate equation!DC motor} is
\begin{equation}
\dot{p}=-\frac{\partial H}{\partial \omega}=\frac{2K_b\left( e_a-K_b\omega\right)}{R_a}+\frac{K_iK_b}{R_aI_{eq}}p=2(e_a-2\omega)+2p.
\label{eq:costfin7}
\end{equation}
We substitute~(\ref{eq:confin7}) into~(\ref{eq:costfin7}) to obtain
\[
\dot{p}=0.
\]
Hence,
\[
p(t)=\mb{constant}=c.
\]
Therefore,
\begin{equation}
e_a=2\omega-c.
\label{eq:contfin7}
\end{equation}
Substituting~(\ref{eq:contfin7}) into~(\ref{eq:dcexsimf7hjb}), we get
\[
\dot{\omega}=-c.
\]
Solving the above yields
\[
\omega(t)=-ct+\omega(0)=-ct,
\]
because $\omega(0)=0$. We next use the boundary condition\index{boundary condition!DC motor} $\omega(1)=10$ to find $c=-10$. Hence, the optimal armature voltage\index{optimal armature voltage!DC motor} is
\[
\fbox{$\displaystyle e_a^{*}(t)=20t+10$}
\]
\end{example}
\markex

\begin{example}\label{ex:brachisto}\rm In this example, we will solve the Johann Bernoulli's brachistochrone\index{Johann Bernoulli's brachistochrone problem}\index{brachistochrone problem!of Johann Bernoulli's } problem, introduced in Section~\ref{sec:calcofvar}, using the minimum principle. We first recall the problem as re-stated by Johann Bernoulli in 1697.

``Among the infinitely many curves which join the two given points ... choose the one such that, if the curve is replaced by a thin tube or groove, and a small sphere placed in it and released, then this will pass from one point to the other in the shortest time.'' quoted from W.~Dunham~\cite[page~200]{dunham90}.

It is difficult to resist not to mention about the reward that one to obtain for solving the problem:


``Let who can seize quickly the prize which we have promised to the solver. Admittedly this prize is neither of gold nor silver, for these appeal only to base and venal souls. ... Rather, since virtue itself is its own most desirable reward and a fame is a powerful incentive, we offer the prize fitting for the man of noble blood, compounded of honor, praise, and approbation...'' quoted from W.~Dunham~\cite[page~200]{dunham90}. We add that Johann Bernoulli wrote the above after he already knew the solution to the problem!

\noindent{\bf Model and the Performance Index}

The field is conservative\index{conservative field!brachistochrone problem}, hence the sums of the kinetic and potential energy is constant,\index{kinetic energy!brachistochrone problem}\index{potential energy!brachistochrone problem}
\[
\frac{1}{2}mv(t)^2-mgy(t)= \frac{1}{2}mv^2(t_0)-mgy(t_0)=0,
\]
where we assumed that the particle is at rest at $t_0$. Using the above, we can find the particle's speed at any time $t\ge t_0$,
\[
v(t)=\sqrt{2 g y(t)}.
\]
Then, the brachistochrone problem can be formulated as an optimal control problem of the form,
\begin{eqnarray*}
\mb{minimize}&{}&\; J(t_0)=\int_{t_0}^{T} dt\\
\mb{subject to}&{}& \left\{\begin{array}{lll}
\dot{x}&=& v\cos \theta\\
\dot{y}&=& v\sin \theta,
\end{array}\right.
%{}&{}& \dot{y}=v\sin \theta,
\end{eqnarray*}
where $\theta$ is the angle that the particle velocity vector forms with the positive $x$-axis, while $\dot{x}$ and $\dot{y}$ are the horizontal and vertical components of the particle velocity vector. The path angle $\theta$ is considered as the control input to be determined.

\noindent{\bf Particle Trajectory}

To proceed, we form the Hamiltonian function\index{Hamiltonian function!brachistochrone problem},
\begin{equation}
H=1+p_1 v \cos\theta + p_2 v \sin\theta,
\label{eq:hamilton}
\end{equation}
where $p_1$ and $p_2$ are the costate variables. The costate equations are:
\[
\left[\begin{array}{c}
\dot{p}_1\\
\dot{p}_2
\end{array}\right]=-\left[\begin{array}{c}
\frac{\partial H}{\partial x}\\
\frac{\partial H}{\partial y}
\end{array}\right]=\left[\begin{array}{l}
0\\
-\frac{g}{v}(p_1\cos\theta + p_2 \sin \theta)
\end{array}\right].
\]
It follows from the above that $p_1=K=\mb{constant}$. Because there are no constraints on the control input $\theta$, the necessary condition for the control optimality is\index{condition for the control optimality!brachistochrone problem}
\[
\frac{\partial H}{\partial \theta}=0,
\]
where
\[
\frac{\partial H}{\partial \theta}=-p_1 v \sin \theta + p_2 v \cos \theta.
\]
Hence,
\[
p_2=p_1\tan \theta=K\tan \theta .
\]
We next find an expression for $\dot{\theta}$. First note that
\begin{equation}
\dot{p}_2=\frac{\partial p_2}{\partial \theta} \dot{\theta}=\frac{\partial (K\tan \theta)}{\partial \theta}\dot{\theta}=\frac{K}{\cos ^2 \theta}\dot{\theta}.
\label{eq:theatdottot}
\end{equation}
On the other hand,
\begin{eqnarray}
\dot{p}_2 &=& -\frac{g}{v}(p_1\cos \theta +p_2 \sin \theta)\nonumber\\
&=& -\frac{g}{v}(K\cos \theta +K\tan \theta \sin \theta)\nonumber\\
&=& -\frac{K g}{v}(\cos \theta + \frac{\sin^2 \theta}{\cos \theta} )\nonumber\\
&=& -\frac{K g}{v}\frac{\cos ^2 \theta + \sin^2 \theta}{\cos \theta} )\nonumber\\
&=& -\frac{K g}{v \cos \theta}.
\label{eq:altefordopt}
\end{eqnarray}
Combining~(\ref{eq:theatdottot}) and~(\ref{eq:altefordopt}) yields
\begin{equation}
\dot{\theta}=-\frac{g}{v}\cos\theta.
\label{eq:thetadot}
\end{equation}
We now use the boundary conditions\index{boundary conditions!brachistochrone problem} to express $x$ and $y$ as functions of time and $\theta$. First, using the above we can write
\[
\dot{y}=\frac{d y}{d \theta}\dot{\theta}=\frac{d y}{d \theta}\left(-\frac{g}{v}\cos\theta \right).
\]
On the other hand,
\[
\dot{y}=v\sin \theta.
\]
Equating the right-hand sides of the above two equations yields
\[
\frac{dy}{v^2}=\frac{dy}{2 g y}=-\frac{1}{g}\tan \theta d\theta.
\]
Integrating both sides gives
\[
\frac{1}{2g}\int \frac{dy}{y}=-\frac{1}{g}\int \tan \theta d\theta,
\]
that is,
\[
\frac{1}{2g}\ln |y|=\frac{1}{g}\ln |\cos \theta |.
\]
Performing manipulations, we obtain
\[
y=C \cos ^2 \theta,
\]
where $C$ is an integration constant. Evaluating the above at $t=T$ gives
\[
y(T)=y_1=C \cos ^2 \theta(T).
\]
Hence,
\[
C=\frac{y_1}{\cos ^2 \theta(T)}
\]
and
\begin{equation}
\fbox{$\displaystyle y(t)= \frac{y_1}{\cos ^2 \theta(T)}\cos^2\theta(t) $}
\label{eq:yofcycloid}
\end{equation}
In a similar way, we can obtain an expression for $x$ as a function of time and $\theta$. Indeed, on the one hand,
\[
\dot{x}=\frac{d x}{d \theta}\dot{\theta}=-\frac{d x}{d \theta}\frac{g}{v}\cos\theta,
\]
on the other hand,
\[
\dot{x}=v \cos \theta.
\]
Comparing the right-hand sides of the above two equations gives
\[
\frac{d x}{d \theta}=-\frac{v^2}{g}=-2y=-2\frac{y_1}{\cos ^2 \theta(T)}\cos\theta.
\]
We represent the above as
\[
d x = -2y=-2\frac{y_1}{\cos ^2 \theta(T)}\cos\theta d\theta .
\]
Integrating both sides, using the identity, $\cos ^2 \theta = \frac{1}{2}(1+\cos 2 \theta)$, yields
\[
x=-2\frac{y_1}{\cos ^2 \theta(T)}\left(\frac{\theta}{2} + \frac{\sin 2 \theta}{4}  \right) + \tilde{C},
\]
where $\tilde{C}$ is the integration constant. Evaluating the above at $t=T$ gives
\[
x_1=-2\frac{y_1}{\cos ^2 \theta(T)}\left(\frac{\theta(T)}{2} + \frac{\sin 2 \theta(T)}{4}  \right) + \tilde{C}.
\]
Hence,
\[
\tilde{C}=x_1+\frac{y_1}{\cos ^2 \theta(T)}\left(\theta(T) + \frac{\sin 2 \theta(T)}{2}  \right).
\]
Using the above, we obtain
\begin{equation}
\fbox{$\displaystyle x(t)=x_1+\frac{y_1}{2\cos ^2 \theta(T)}\left(2(\theta(T)-\theta(t)) + \sin 2 \theta(T)-\sin 2\theta(t) \right) $}
\label{eq:xofcycloid}
\end{equation}
We point out that the equation for $x(t)$, $y(t)$ and $\dot{\theta}$ specify an implicit feedback control law\index{implicit feedback control!brachistochrone problem} in the sense that the control $\theta$ can be calculated given the state $(x(t), y(t))$.


We use the boundary conditions\index{boundary conditions!brachistochrone problem} to find the description of the time-optimal closed loop system. We assume, without loss of generality, that $x(t_0)=y(t_0)=x(0)=y(0)=0$. Then, using the expression for $y(t)$ we note that $y(t_0)=0=\frac{y_1}{\cos ^2\theta(T)}\cos ^2 \theta(t_0) $. Thus, the angle of departure  is $\theta(t_0)=90^{\circ}$. With the above initial path angle, $y$, and hence $v$, is increased at the outset as quickly as possible. Using this data, we next evaluate $x$ at $t=t_0$ to form an equation that can be used to find the angle of arrival $\theta(T)$. The equation is
\[
2x_1\cos^2 \theta(T) +y_1\left(2\theta(T)-\pi +\sin 2\theta(T)\right) = 0.
\]
The above equation can be solved easily using MATLAB's Symbolic Toolbox or the MATLAB's \verb!fsolve! function.\index{MATLAB's function!fsolve} Once $\theta(T)$ is found, we use the expression for $y$ to find $\cos\theta$, where
\[
\cos\theta(t)=\sqrt{\frac{y(t)}{y_1}}\cos\theta(T).
\]
Hence,
\[
\sin\theta (t)=\sqrt{1-\cos^2\theta (t)}.
\]
We substitute the above expression into $\dot{x}=v\cos\theta$ and $\dot{y}=v\sin\theta$ to obtain the state-space model of the time optimal motion of the particle\index{optimal motion of the particle!brachistochrone problem} moving under the force of gravity from a given point to another specified point. The equations are,
\[
\left\{\begin{array}{lll}
\dot{x}(t) &=& \sqrt{2 g y(t)}\sqrt{\frac{y(t)}{y_1}}\cos\theta(T)\\
\dot{y}(t)&=& \sqrt{2 g y(t)}\sqrt{1-\frac{y(t)}{y_1}\cos^2\theta(T)}
\end{array}\right.
\]
subject to the boundary conditions $\left[\begin{array}{cc}
x(t_0) & y(t_0)
\end{array}\right]=\left[\begin{array}{cc}
x_0 & y_0
\end{array}\right]$ and $\left[\begin{array}{cc}
x(T) & y(T)
\end{array}\right]=\left[\begin{array}{cc}
x_1 & y_1
\end{array}\right]$.

\noindent{\bf The Time of Quickest Descent}

We derive an expression that can be used to compute the time of the particle travel.\index{travel time!brachistochrone problem} We first establish an auxiliary result that we will use in our derivation. Specifically, we will show that $\dot{\theta}(t)$ is constant on the interval $[t_0, T]$. Differentiating both sides of equation~(\ref{eq:thetadot}) gives
\begin{eqnarray*}
\frac{d \dot{\theta}}{dt} &=& \frac{d}{dt}\left(-\frac{g}{\sqrt{2g y}} \cos \theta \right)\\
&=& -\frac{g}{\sqrt{2g}}\left(\frac{d}{dt}\left(\frac{1}{\sqrt{y}}\cos \theta \right) \right)\\
&=& -\frac{g}{\sqrt{2g}}\left(-\frac{1}{2} \frac{\dot{y} \cos \theta}{y^{3/2}}-\frac{\sin \theta}{\sqrt{y}}\dot{\theta} \right)\\
&=& \frac{g}{\sqrt{2g}}\left(\frac{1}{2} \frac{\dot{y} \cos \theta}{y^{3/2}}+\frac{\sin \theta}{\sqrt{y}}\dot{\theta} \right).
\end{eqnarray*}
Substituting into the above the expressions for $\dot{\theta}$ and $\dot{y}$ gives
\begin{eqnarray*}
\frac{d \dot{\theta}}{dt} &=& \frac{g}{\sqrt{2g}}\left(\frac{1}{2} \frac{\sqrt{2gy}\sin\theta \cos \theta}{y^{3/2}}- g\frac{\sin \theta \cos\theta }{\sqrt{y}\sqrt{2g y}} \right)\\
&=& \frac{g}{\sqrt{2g}}\left( \frac{\sqrt{g}}{y\sqrt{2}}\sin\theta \cos \theta - \frac{\sqrt{g}}{y\sqrt{2}}\sin\theta \cos \theta\right)\\
&=& 0.
\end{eqnarray*}
Because $\dot{\theta}$ is constant on the interval $[t_0, T]$, we use~(\ref{eq:thetadot}) to write
\[
\frac{d \theta}{dt}=-\frac{g}{\sqrt{2 g y(T)}}\cos\theta(T).
\]
Separating the variables yields
\[
d\theta = -\frac{g}{\sqrt{2 g y(T)}}\cos\theta(T) dt.
\]
Integrating both sides of the above, we obtain
\[
\int_{\theta(t)}^{\theta(T)} d\theta = -\frac{g}{\sqrt{2 g y(T)}}\cos\theta(T) \int_{t_0}^{T} dt.
\]
Hence,
\[
\theta(T)-\theta(t)=-\frac{g}{\sqrt{2 g y(T)}}\cos\theta(T)(T-t_0).
\]
To obtain an expression for the particle travel time, we set in the above equation $t=t_0=0$ and perform simple manipulations to obtain
\[
\fbox{$\displaystyle T=\frac{\pi/2-\theta(T)}{\cos\theta(T)}\sqrt{{2 y(T)}{g}} $}
\]

\noindent{\bf It's an Upside-Down Cycloid!}

The path of quickest descent is an upside-down cycloid.\index{path of quickest descent!upside-down cycloid} A cycloid is a curve described by a point on the circumference of a circle that rolls without slipping---see Figure~\ref{fig:brachis}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{brachis.eps}
%\centerline{\psfig{figure=c:/mybook/brachis.eps,height=2in}}
\caption{Constructing an upside-down cycloid---the solution to the brachistochrone problem.}
\label{fig:brachis}
\end{figure}
The equations describing a cycloid are
\[
\left.\begin{array}{lll}
x(\phi) &=& R(\phi-\sin \phi)\\
y(\phi) &=& R(1 - \cos \phi)
\end{array}\right\}
\]
We will now show that equations~(\ref{eq:xofcycloid}) and~(\ref{eq:yofcycloid}) describe a cycloid. To proceed, we define
\begin{equation}
R=\frac{y_1}{2\cos^2 \theta(T)}\quad\mb{and}\quad \phi(t)=\pi-2\theta(t).
\label{eq:rphinotation}
\end{equation}
Using the above notation, and taking into account that $\sin(\pi-\phi)=\sin \phi$, we re-write (\ref{eq:xofcycloid}) as
\[
x(t)-x_1+R(\phi(T)-\sin \phi(T))= R(\phi(t)-\sin\phi(t)).
\]
Using the notation given in~(\ref{eq:rphinotation}) and the fact that $\cos 2\phi = 2\cos^2\phi -1$, we represent~(\ref{eq:yofcycloid}) as
\begin{eqnarray*}
y(t) &=& \frac{y_1}{\cos ^2 \theta(T)}\cos^2\theta(t)\\
&=& \frac{y_1}{\cos ^2 \theta(T)}\frac{\cos 2\theta(t) + 1}{2}\\
&=& R(1 + \cos(\pi - \phi(t))\\
&=& R(1-\cos \phi(t)),
\end{eqnarray*}
because $\cos(\pi - \phi(t))=-\cos \phi$.
\end{example}
\markex

In the following, we discuss a general minimum time regulator\index{minimum time regulator problem} problem for a class of multi-input dynamical systems modeled by
\[
\dot{\b x}=\b A\b x+\b B\b u.
\]
The goal is to transfer the system from some initial state to the origin in minimum time subject to
\[
|u_i|\le k_i,\quad k_i>0,\quad i=1,2,\ldots,m.
\]
Thus, the performance index for the problem is $J=\int_{t_0}^{t_f}dt$. We begin our analysis by forming the associated Hamiltonian function,\index{Hamiltonian function!minimum time regulator problem}
\begin{eqnarray*}
H &=& 1+\b p^{\top}(\b A\b x+\b B\b u)\\
&=& 1+\b p^{\top}\b A\b x+\sum_{i=1}^m\left(\b p^{\top}\b b_i\right)u_i,
\end{eqnarray*}
where $\b b_i$ is the $i$-th column of $\b B$. Let
\[
s_i(t)=\b p^{*\top}(t)\b b_i.
\]
Then, applying Pontryagin's minimum principle\index{Pontryagin's minimum principle!minimum time regulator problem} we obtain the following form of the optimal control,
\[
u^{*}_i(t)=-k_i\mb{sign}(s_i(t)),\quad i=1,2,\ldots,m.
\]
Thus, the optimal control for the $i$-th control channel has the bang-bang form, $u^*_i=\pm k_i$. Note that $s_i(t)=\b p^{*\top}(t)\b b_i=0$ is the switching curve for the $i$-th control channel. We compute the costate vector from the costate equation,\index{costate equation!minimum time regulator problem}
\begin{eqnarray*}
\dot{\b p}^{*\top} &=& -\frac{\partial H}{\partial \b x}\\
&=& -\frac{\partial}{\partial \b x}\left(\b p^{*\top}\b A\b x\right)\\
&=& -\b p^{*\top}\b A,
\end{eqnarray*}
or
\[
\dot{\b p}^* = -\b A^{\top}\b p^*.
\]
The solution to the costate equation is
\[
\b p^*(t)=e^{-\b A^{\top}t}\b p^*(0).
\]
Thus, the $i$-th switching curve\index{switching curve!minimum time regulator problem} can be represented as
\[
s_i(t) = \b p^{*}(0)^{\top}e^{-\b At}\b b_i= 0.
\]
Observe that if $s_i(t)=0$ on some time interval, then $u^*_i(t)$ is indefinite on this interval. We show that under some mild assumptions there is no time interval on which $s_i(t)=0$, that is, there is no time interval on which $u^*_i(t)$ is indefinite. First, we assume that $\b b_i\ne \b 0$. Note that in our problem here, $\min_{\b u\in U}H(\b x,\b u)=0$ because $J^*(t,\b x)=J^*(\b x)$ and therefore $\frac{\partial J^*}{\partial t}=0$. Hence, $1+\b p^{*}(t)^{\top}(\b A\b x^*+\b B\b u^*)=0$ for all $t$ for the optimal process. Thus, the optimal costate vector $\b p^*(t)$ cannot be zero for any value of $t$. Finally, if $s_i=\b p^{*\top}\b b_i$ is zero on some time interval, that is, $s_i$ is constant on a time interval, then this implies that on the above time interval,
\[
\dot{s}_i = \dot{\b p}^{*\top}\b b_i= -\b p^{*\top}\b A\b b_i=0.
\]
Similarly, higher derivatives of $s_i(t)=\b p^{*\top}\b b_i$ will also be zero on the time interval on which $s_i(t)=\b p^{*\top}\b b_i$ is zero, that is, $s^{(j)}_i(t)=(-1)^j\b p^{*\top}\b A^j\b b_i=0$, $j=1,2,\ldots$. We represent the above relations as
\[
\b p^{*\top}\left[\begin{array}{ccccc}
\b b_i& \b A\b b_i&\cdots&\b A^{n-1}\b b_i&\cdots
\end{array}\right]=\b 0.
\]
If the system is controllable\index{controllable system!minimum time regulator problem} by the i-th input acting alone, then
\[
\mb{rank}\left[\begin{array}{ccccc}
\b b_i& \b A\b b_i&\cdots&\b A^{n-1}\b b_i&\cdots
\end{array}\right]=n,
\]
and we would have to have $\b p^*=\b p^*(t)=\b 0$. However, we already ruled out this case and thus $s_i(t)=\b p^{*\top}\b b_i$ cannot be zero on any time interval. In conclusion, if the system is controllable by each input acting alone, then there is no time interval on which optimal control vector is indefinite.

\subsection{A Two-Point Boundary-Value Problem}\label{subsec:TPBVP}
We consider the problem of minimizing the quadratic performance index
\begin{equation}
J=\frac{1}{2}\b x(t_f)^{\top}\b F \b x(t_f)+\frac{1}{2}\int_{t_0}^{t_f}\left(\b x(t)^{\top}\b Q\b x(t)+\b u(t)^{\top}\b R\b u(t) \right)dt
\label{eq:qperintpbv}
\end{equation}
subject to
\begin{equation}
\dot{\b x}(t)=\b A\b x(t)+\b B\b u(t),\quad \b x(t_0)=\b x_0,\quad \b x(t_f)\;\;\;\mb{is free}.
\label{eq:lindiffeq}
\end{equation}
We assume that there are no constraints on the control input $\b u$. The weight matrices\index{weight matrices} $\b F$ and $\b Q$ are real symmetric positive semi-definite, and the matrix $\b R$ is real symmetric positive definite. The Hamiltonian function for the above linear quadratic control problem is
\begin{equation}
H=\frac{1}{2}\b x^{\top}\b Q\b x+\frac{1}{2}\b u^{\top}\b R\b u +\b p^{\top}\left(\b A\b x+\b B\b u \right),
\label{eq:hamilforlqrc}
\end{equation}
where $\b p \in \mb{\Bb R}^n$ is the costate vector. It follows from Theorem~\ref{thm:pontryagins} that the optimal controller must minimize the Hamiltonian function\index{optimal controller!must minimize the Hamiltonian}. Because the control vector is unconstrained, the necessary condition for optimality of the control $\b u$ is
\begin{equation}
\frac{\partial H}{\partial \b u}=\b 0.
\label{eq:necconforopu}
\end{equation}
Evaluating~(\ref{eq:necconforopu}) yields
\begin{equation}
\frac{\partial H}{\partial \b u}=\b u^{\top}\b R+\b p^{\top}\b B=\b 0.
\label{eq:evalderofh}
\end{equation}
Hence, the optimal controller has the form
\begin{equation}
\b u=-\b R^{-1}\b B^{\top}\b p.
\label{eq:optuncco}
\end{equation}
Substituting the above into~(\ref{eq:lindiffeq}), we obtain
\begin{equation}
\dot{\b x}=\b A\b x-\b B\b R^{-1}\b B^{\top}\b p,\quad \b x(t_0)=\b x_0.
\label{eq:aftsubuinxd}
\end{equation}
We next find the costate equation,
\begin{eqnarray}
\dot{\b p} &=& -\left(\frac{\partial H}{\partial \b x} \right)^{\top}\nonumber\\
&=& -\b Q\b x-\b A^{\top}\b p.
\label{eq:costeqforlqc}
\end{eqnarray}
Combining~(\ref{eq:aftsubuinxd}) and the above yields $2n$ linear differential equations
\begin{equation}
\left[\begin{array}{c}
\dot{\b x}\\
\dot{\b p}
\end{array}\right]=\left[\begin{array}{cc}
\b A & -\b B\b R^{-1}\b B^{\top}\\
-\b Q & -\b A^{\top}
\end{array}\right]\left[\begin{array}{c}
\b x\\
\b p
\end{array}\right].
\label{eq:bigeqtn}
\end{equation}
By Theorem~\ref{thm:pontryagins}, the state vector at time $t_f$ must satisfy
\begin{equation}
\b p(t_f)=\frac{1}{2}\left.\left(\nabla\b x^{\top}\b F \b x\right)\right\arrowvert_{t=t_f}=\b F\b x(t_f).
\label{eq:terconfropo}
\end{equation}
In summary, we reduced the linear quadratic control problem to solving $2n$ linear differential equations with mixed boundary conditions, which is an example of a {\em two-point boundary value problem\/}, or TPBVP for short.\index{two-point boundary value problem (TPBVP)}\index{TPBVP|see{two-point boundary value problem}} We solve the above TPBVP. If we had the initial conditions $\b x(t_0)$ and $\b p(t_0)$, then the solution to~(\ref{eq:bigeqtn}) would have the form
\begin{equation}
\left[\begin{array}{c}
{\b x(t_f)}\\
{\b p(t_f)}
\end{array}\right]=e^{\b H(t_f-t_0)}\left[\begin{array}{c}
\b x(t_0)\\
\b p(t_0)
\end{array}\right],
\label{eq:soltotnebp}
\end{equation}
where
\[
\b H=\left[\begin{array}{cc}
\b A & -\b B\b R^{-1}\b B^{\top}\\
-\b Q & -\b A^{\top}
\end{array}\right].
\]
However, in this particular TPBVP, $\b p(t_0)$ is unknown and instead we are given $\b p(t_f)$. To proceed, let
\begin{equation}
\left[\begin{array}{c}
{\b x(t_f)}\\
{\b p(t_f)}
\end{array}\right]=e^{\b H(t_f-t)}\left[\begin{array}{c}
\b x(t)\\
\b p(t)
\end{array}\right],
\label{eq:soltotnwtio}
\end{equation}
where
\[
e^{\b H(t_f-t)}=\left[\begin{array}{cc}
\b{\Phi}_{11}(t_f,t) & \b{\Phi}_{12}(t_f,t)\\
\b{\Phi}_{21}(t_f,t) & \b{\Phi}_{22}(t_f,t)
\end{array}\right].
\]
Each of the blocks $\b{\Phi}_{ij}(t_f,t)$, $i,j=1,2$, is $n\times n$. Thus we have
\begin{equation}
\left[\begin{array}{c}
{\b x(t_f)}\\
{\b p(t_f)}
\end{array}\right]=\left[\begin{array}{cc}
\b{\Phi}_{11}(t_f,t) & \b{\Phi}_{12}(t_f,t)\\
\b{\Phi}_{21}(t_f,t) & \b{\Phi}_{22}(t_f,t)
\end{array}\right]\left[\begin{array}{c}
\b x(t)\\
\b p(t)
\end{array}\right].
\label{eq:solwitphibls}
\end{equation}
The above matrix equation can be represented as
\begin{eqnarray}
\b x(t_f)&=& \b{\Phi}_{11}(t_f,t)\b x(t) + \b{\Phi}_{12}(t_f,t)\b p(t)\label{eq:firsforx},\\
\b p(t_f) &=& \b{\Phi}_{21}(t_f,t)\b x(t)+\b{\Phi}_{22}(t_f,t)\b p(t).
\label{eq:secfoptf}
\end{eqnarray}
Using the boundary condition~(\ref{eq:terconfropo}) and~(\ref{eq:firsforx}) gives
\begin{equation}
\b p(t_f)=\b F\b x(t_f)=\b F\b{\Phi}_{11}(t_f,t)\b x(t) + \b F\b{\Phi}_{12}(t_f,t)\b p(t).
\label{eq:comthbxap}
\end{equation}
Subtracting~(\ref{eq:secfoptf}) from~(\ref{eq:comthbxap}), we obtain
\begin{equation}
\b 0=\left(\b F\b{\Phi}_{11}(t_f,t)-\b{\Phi}_{21}(t_f,t) \right)\b x(t)+\left(\b F\b{\Phi}_{12}(t_f,t)-\b{\Phi}_{22}(t_f,t) \right)\b p(t).
\label{eq:aftsutps}
\end{equation}
Using the above, we compute
\begin{eqnarray}
\b p(t)&=&\left( \b{\Phi}_{22}(t_f,t)-\b F\b{\Phi}_{12}(t_f,t)\right)^{-1}\left( \b F\b{\Phi}_{11}(t_f,t)-\b{\Phi}_{21}(t_f,t)\right)\b x(t)\nonumber\\
&=& \b P(t)\b x(t),
\label{eq:compasfx}
\end{eqnarray}
where,
\begin{equation}
\b P(t)=\left( \b{\Phi}_{22}(t_f,t)-\b F\b{\Phi}_{12}(t_f,t)\right)^{-1}\left( \b F\b{\Phi}_{11}(t_f,t)-\b{\Phi}_{21}(t_f,t)\right).
\label{eq:soldare}
\end{equation}
Substituting~(\ref{eq:compasfx}) into~(\ref{eq:optuncco}) yields the optimal state-feedback controller
\begin{equation}
\b u(t)=-\b R^{-1}\b B^{\top}\b P(t)\b x(t).
\label{eq:optlqrtvar}
\end{equation}

We show that $\b P(t)$ satisfies a matrix differential equation\index{matrix differential equation!of TPBVP}. Indeed, differentiating~(\ref{eq:compasfx}) gives
\begin{equation}
\dot{\b P}\b x+\b P\dot{\b x}-\dot{\b p}=\b 0.
\label{eq:aftdippx}
\end{equation}
Substituting into the above~(\ref{eq:aftsubuinxd}), (\ref{eq:costeqforlqc}), and~(\ref{eq:compasfx}) yields
\begin{equation}
\left(\dot{\b P}+\b P\b A-\b P\b B\b R^{-1}\b B^{\top}\b P+\b Q+\b A^{\top}\b P  \right)\b x=\b 0.
\label{eq:afsuxdpd}
\end{equation}
The preceding must hold throughout $t_0\le t\le t_f$. Hence, $\b P$ must satisfy
\begin{equation}
\fbox{$\displaystyle \dot{\b P}=\b P\b B\b R^{-1}\b B^{\top}\b P-\b A^{\top}\b P-\b P\b A-\b Q $}
\label{eq:mariccdieqtpbvp}
\end{equation}
subject to the boundary condition
\begin{equation}
\b P(t_f)=\b F.
\label{eq:bouncoofdirixe}
\end{equation}
Equation~(\ref{eq:mariccdieqtpbvp}) is called the matrix {\em Riccati differential equation\/}\index{Riccati differential equation}. Note that because $\b F$ is symmetric, so is $\b P=\b P(t)$. Observe that $\b P=\b P(t)$ given by~(\ref{eq:soldare}) is a solution to~(\ref{eq:mariccdieqtpbvp}) subject to the boundary condition~(\ref{eq:bouncoofdirixe}).

The Riccati differential equation~(\ref{eq:mariccdieqtpbvp}) associated with a class of one-dimensional plants was solved in Example~\ref{ex:soldifriceq}.
\section{Notes}
Sussmann and Willems~\cite{sussmann97} made the following statement, in 1997, regarding the origins of optimal control: ``Optimal control was born in 1697---300 years ago---in Groningen, a university town in the north of The Netherlands, when Johann Bernoulli, professor of mathematics at the local university from 1695 to 1705, published his solution of the {\it brachystochrone problem\/}.''\index{brachystochrone problem|see{brachistochrone problem}}

For a serious student of optimal control, the seminal text by Lee and Markus~\cite{eblee86} is recommended. Less advanced than Lee and Markus are Kirk~\cite{kirk70}, Bryson and Ho~\cite{bryson75}, Owens~\cite{owens81}, and Pierre~\cite{pierre86}. The classical text by Kwakernaak and Sivan~\cite{kwakernaak72} is also worth perusing because of its clarity and its treatment of the stochastic aspects of control problems. Informative, and at the same time entertaining, accounts of the origins of the calculus of variations\index{calculus of variations!origins of} are given by Dunham~\cite[Chapter 8]{dunham90}, and Sussmann and Willems~\cite{sussmann97}. Elsgolc~\cite{elsgolc61} provides a lucid introduction, illustrated with many examples, to the calculus of variations. Applications of the calculus of variations in physics and engineering are provided by Weinstock~\cite{weinstock74}. Interrelations between the calculus of variations and optimal control are discussed by Hestenes~\cite{hestenes66}. Frequency-domain aspects of optimal control and robust optimal control\index{robust optimal control} design is discussed by Burl in~\cite{burl99}. Our presentation of the continuous minimum time regulator problem in Subsection~\ref{subsub:con_min_time_r}, was based on Boltyanskii~\cite{boltyanskii71}. For more numerical examples related to optimal control, see Lewis and Syrmos~\cite{lewis95}, and Bryson~\cite{bryson99}. Exciting applications of optimal control to a number of control problems of technology and science can be found in Lee and Markus~\cite{eblee86} .
\section{Exercises}
\begin{exercise}\rm
Find the variations of the functionals\index{variations of the functionals}:
\begin{itemize}
\item[(i)] $v(x_1,x_2)=\int_{0}^{\pi/2}\left(x_1^4 +2x_1^2x_2^2+x_2^2 \right)dt$,
\item[(ii)] $v(x_1,x_2)=\int_{0}^1\left( e^{x_1-x_2}-e^{x_1+x_2}\right) dt$.
\end{itemize}
\end{exercise}
\begin{exercise}\rm
Find the equations of the curves that are extremals of the functional\index{extremals of the functional}
\[
v(x)=\int_0^{t_f}\left( \frac{1}{4}\dot{x}^2-t\dot{x}-x+\frac{1}{2}x\dot{x}\right)dt
\]
for the boundary conditions
\begin{itemize}
\item[(i)] $x(0)=1$, $t_f=2$, and $x(2)=10$;
\item[(ii)] $x(0)=1$, $t_f=10$, and $x(10)$ is free.
\end{itemize}
\end{exercise}
\begin{exercise}\rm
Find the equation of the curve that is an extremal of the functional
\[
v(x)=\int_0^{t_f}\left( \frac{1}{2}\dot{x}^2-t\dot{x}+2x^2\right)dt
\]
subject to the boundary conditions
\[
x(0)=3/4,\; t_f=1,\;\;\mb{and}\;\; x(1)\;\;\mb{is free}.
\]
\end{exercise}
\begin{exercise}\rm
Given a dynamical system model
\begin{eqnarray*}
\dot{\b x}&=&\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\b x+\left[\begin{array}{c}
0\\
1
\end{array}\right]u\\
y &=& \left[\begin{array}{cc}
2 & 0
\end{array}\right]\b x,
\end{eqnarray*}
and the associated performance index
\[
J=\int^{\infty}_0\left( y^2+u^2\right)dt.
\]
Assuming that $u=-\b k\b x$ is such that $J$ is minimized, find the poles of the closed-loop system.
\end{exercise}
\begin{exercise}\rm
Given a dynamical system model
\begin{eqnarray*}
\dot{\b x}&=&\left[\begin{array}{cc}
0 & 1\\
0 & 0
\end{array}\right]\b x+\left[\begin{array}{c}
0\\
1
\end{array}\right]u\\
y &=& \left[\begin{array}{cc}
2 & 0
\end{array}\right]\b x,
\end{eqnarray*}
and the performance index
\[
J=\int^{\infty}_0\left( y^2+u^2\right)dt.
\]
Use the Kronecker product\index{Kronecker product} to represent the associated algebraic Riccati equation\index{algebraic Riccati equation} as the nonlinear vector algebraic equation. Solve it, and construct the optimal state feedback control law $u=-\b k\b x$. Assume that
\[
\b x(0)=\left[\begin{array}{c}
1\\
2
\end{array}\right].
\]
Find the optimal value of $J$.
\end{exercise}
\begin{exercise}\rm
Find a state feedback control law $u$ so as to minimize
\[
J=\int_0^1\left( x^2+u^2\right)dt
\]
subject to
\[
\dot{x}=-x+u,\quad x(0)=2.
\]
Assume $J^{*}=p(t)x^2$.
\end{exercise}
\begin{exercise}\rm
Use the principle of optimality\index{principle of optimality} to find the set of numbers $x_1,x_2,\ldots,x_n$ satisfying $x_1+x_2+\cdots +x_n=c$, where $c$ is a positive constant such that the product $x_1x_2\cdots x_n$ is maximized.
\end{exercise}
\begin{exercise}\label{ex:dyprowithtoll}\rm
Use dynamic programming to find the maximum cost path between node {\bf a} and node {\bf t} in the routing network depicted in Figure~\ref{fig:test2f1}, where the travel costs are shown beside each segment and toll charges are shown by each node. The network can only be traversed from left to right.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{test2f1t.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/test2f1t.bmp}
%\hskip1in\special{bmp:c:/ee680/test2f1t.bmp x=4.67in y=3.5in}
\caption{A routing network of Exercise~\ref{ex:dyprowithtoll}.}
\label{fig:test2f1}
\end{figure}
\end{exercise}
\begin{exercise}\rm
Use dynamic programming\index{dynamic programming} to find $u[0]$ and $u[1]$ that minimize
\[
J=(x[2]-2)^2+\sum_{k=0}^1u[k]^2
\]
subject to
\[
x[k+1]=x[k]+u[k],\quad x[0]=1.
\]
\end{exercise}
\begin{exercise}\rm Consider the following optimal control problem\index{optimal control problem}:
\begin{eqnarray*}
\mb{minimize}&{}&\; J_2(u)=\frac{1}{2}\sum_{k=0}^1 u[k]^2\\
\mb{subject to}&{}& x[k+1]=x[k]+2u[k],\quad x[0]=3,\quad x[2]=0.
\end{eqnarray*}
Convert the above optimal control problem into a quadratic optimization\index{quadratic optimization} problem with an equality constraint using the composite input vector
\[
\b z=\left[\begin{array}{c}
u[0]\\
u[1]
\end{array}\right].
\]
Then, employ Lagrange's multiplier theorem\index{Lagrange's multiplier theorem} to find the optimal sequence $\{u^*[0], u^*[1]\}$.
\end{exercise}
\begin{exercise}\rm Consider the following optimal control problem:
\begin{eqnarray*}
\mb{minimize}&{}&\; J_2(u)=\frac{1}{2}x[2]^2+\frac{1}{2}\sum_{k=0}^1 \left(x[k]^2 + u[k]^2\right)\\
\mb{subject to}&{}& x[k+1]=x[k]+2u[k],\quad x[0]=3.
\end{eqnarray*}
Convert the above optimal control problem into a quadratic optimization problem with an equality constraint using the composite input-state vector
\[
\b z=\left[\begin{array}{c}
u[0]\\
x[1]\\
u[1]\\
x[2]
\end{array}\right].
\]
Then, employ Lagrange's multiplier theorem to find the optimal sequence\index{Lagrange's multiplier theorem!to find optimal sequence} $\{u^*[0], u^*[1]\}$.
\end{exercise}
\begin{exercise}\rm
Find $u=u(t)$ so as to minimize
\[
J=\int^1_0u^2dt
\]
subject to
\[
\dot{x}=-2x+u,\quad x(0)=10,\quad x(1)=0.
\]
\end{exercise}
\begin{exercise}\rm
For the circuit shown in Figure~\ref{fig:lrcirqe}, find the input voltage, $u(t)$, that changes the current in the inductor form $i=10$~A at $t=0$ to $i=0$ at $t=1$~sec while minimizing the performance index,
\[
J=\int_0^1 u(t)^2dt.
\]
\begin{figure}[h]
\centering
\includegraphics[scale=0.6]{lrcirqe.eps}
%\centerbmp{3.33in}{2.5in}{c:/mybook/lrcirqe.bmp}
\caption{A circuit for Exercise~\ref{exercise:rlcirc}.}
\label{fig:lrcirqe}
\end{figure}
\label{exercise:rlcirc}
\end{exercise}
\begin{exercise}\label{pierre_based}\rm
Use the minimum principle\index{minimum principle} to find the input voltage $u(t)$ that charges the capacitor, in Figure~\ref{fig:finalf1}, from $x(0)=2\,V$ at $t=0$ to $x(10)=12\,V$ while minimizing the energy dissipated in the resistor. There are no constraints on $u(t)$. Note that, by Kirchhoff's voltage law\index{Kirchhoff's voltage law},
\[
Ri(t)+x(t)=u(t),
\]
where $i(t)=C\frac{dx(t)}{dt}$, and the energy dissipated is
\[
J=\int_0^{10}Ri(t)^2dt.
\]
\begin{figure}
\centering
\includegraphics[scale=0.6]{finalf1.eps}
%\centerbmp{3in}{2.25in}{c:/mybook/finalf1.bmp}
%\hskip.6in\special{bmp:c:/mybook/finalf1.bmp x=4.67in y=3.5in}
\caption{A capacitor charging circuit of Exercise~\ref{pierre_based}.}
\label{fig:finalf1}
\end{figure}
\end{exercise}
\begin{exercise}\rm
Show that if $\mu$ is an eigenvalue of the Hamiltonian matrix\index{Hamiltonian matrix} $\b H$, then so is $-\mu$.
\label{ex:mupomuneinhm}
\end{exercise}
\begin{exercise}\rm
Consider the Hamiltonian function defined by
\[
H=F+\b p^{\top}\b f.
\]
Show that if the functions $\b f$ and $F$ do not explicitly depend upon $t$, then
\[
\left.H\right\arrowvert_{\b u=\b u^*}=\mb{constant},\quad t_0\le t\le t_f.
\]
\label{ex:conshamfu}
\end{exercise}
\begin{exercise}\rm
Given the Hamiltonian matrix corresponding to a dynamical system and its associated performance index. The Hamiltonian matrix is
\[
\mb{\bm$H$} = \left[ \begin{array}{cc}
2 & -5\\
-1 & -2
\end{array}\right].
\]
The MATLAB command \verb+[V,D] = eig(H)+ \index{MATLAB command!eig(H)}resulted in the following data:
\[
\mb{\bm$V$}=\left[\begin{array}{cc}
0.9806 & 0.7071\\
-0.1961 & 0.7071
\end{array} \right], \qquad \mb{\bm$D$}=\left[\begin{array}{cc}
3 & 0\\
0 & -3
\end{array} \right].
\]
\begin{itemize}
\item[(i)] Write the equation of the closed loop system driven by the optimal controller $u=-kx$.
\item[(ii)] Find the solution to the algebraic Riccati equation\index{algebraic Riccati equation} corresponding to the optimal controller.
\end{itemize}
\end{exercise}
\begin{exercise}\rm\begin{itemize}\item[{\bf (i)}] Suppose that $s_i$ is an eigenvalue of the matrix
\[
\b A-\b B\b R^{-1}\b B^{\top}\b P,
\]
and $\b v_i\ne \b 0$ is the corresponding eigenvector, where $\b P$ is the symmetric positive definite solution of the algebraic Riccati equation. Show that
\[
\left[\begin{array}{c}
\b v_i\\
\b P\b v_i
\end{array}\right]
\]
is an eigenvector of the associated Hamiltonian matrix\index{Hamiltonian matrix} corresponding to the eigenvalue $s_i$.
\item[{\bf (ii)}] Show that if $\b Q=\b Q^{\top}>0$, then the real parts of the eigenvalues of the matrix
\[
\b A-\b B\b R^{-1}\b B^{\top}\b P
\]
are all negative, where $\b P=\b P^{\top}>0$ is the solution of the associated algebraic Riccati equation.\index{algebraic Riccati equation}
\end{itemize}
\end{exercise}
\begin{exercise}\rm
Given the following model of a dynamical system:
\begin{eqnarray*}
\dot{x}_1 &=& x_2+5\\
\dot{x}_2 &=& u,
\end{eqnarray*}
where
\[
|u| \le 1.
\]
The performance index to be minimized is
\[
J=\int_{0}^{t_f}dt.
\]
Find the state-feedback control law $u=u(x_1,x_2)$ that minimizes $J$ and drives the system from a given initial condition $\b x(0)=[x_1(0),x_2(0)]^{\top}$ to the final state $\b x(t_f)= \b 0$. Proceed as indicated below.
\begin{itemize}
\item[{\bf (a)}] Derive the equations of the optimal trajectories.\index{optimal trajectories}
\item[{\bf (b)}] Derive the equation of the switching curve.\index{switching curve}
\item[{\bf (c)}] Write the expression for the optimal state-feedback controller.\index{optimal state-feedback}
\end{itemize}
\end{exercise}
\begin{exercise}[based on~\cite{rugh69}]\rm
Given a dynamical system model
\[
\dot{\b x}=\b A\b x+\b bu,
\]
where
\[
\b A=\left[\begin{array}{cccccc}
0 & 1 & 0 & \cdots & 0 & 0\\
0 & 0 & 1 & \cdots & 0 & 0\\
\vdots & {} & {} &\cdots & {} & \vdots\\
0 & 0 & 0 & \cdots & 0 & 1\\
-a_0 & -a_1 & -a_2& \cdots & a_{n-2}& a_n
\end{array}\right],\quad \b b=\left[\begin{array}{c}
0\\
0\\
\vdots\\
0\\
1
\end{array}\right],
\]
and the associated performance index
\[
J=\frac{1}{2}\int_0^{\infty}\left( \b x^{\top}\b Q\b x +u^2\right)dt,\quad \b Q=\b Q^{\top} \ge 0.
\]
We label the above problem as $(\b A,\b b)$. Consider now the dynamical system model
\[
\dot{\tilde{\b x}}=\tilde{\b A}\tilde{\b x}+\tilde{\b b}\tilde{u},
\]
where $\b b=\tilde{\b b}$,
\[
\tilde{\b A}=\left[\begin{array}{cccccc}
0 & 1 & 0 & \cdots & 0 & 0\\
0 & 0 & 1 & \cdots & 0 & 0\\
\vdots & {} & {} &\cdots & {} & \vdots\\
0 & 0 & 0 & \cdots & 0 & 1\\
0 & 0 & 0& \cdots & 0& 0
\end{array}\right],
\]
and the associated performance index has the form
\[
\tilde{J}=\frac{1}{2}\int_0^{\infty}\left( \tilde{\b x}^{\top}\tilde{\b Q}\tilde{\b x} +\tilde{u}^2\right)dt,
\]
where $\tilde{\b Q}=\b Q+\b D\tilde{\b A}+\tilde{\b A}^{\top}\b D+\b D\b b\b b^{\top}\b D$, and
\[
\b D=\left[\begin{array}{cccccc}
0 & 1 & 0 & \cdots & 0 & -a_0\\
0 & 0 & 1 & \cdots & 0 & -a_1\\
\vdots & {} & {} &\cdots & {} & \vdots\\
0 & 0 & 0 & \cdots & 0 & -a_{n-2}\\
-a_0 & -a_1 & -a_2& \cdots & -a_{n-2}& -a_{n-1}
\end{array}\right].
\]
We refer to the second problem as $\left(\tilde{\b A},\tilde{\b b} \right)$.
\begin{itemize}
\item [{\bf (a)}] Show that the symmetric matrix $\b P$ is a solution to the $\left( \b A,\b b\right)$ problem if, and only if, the symmetric matrix $\tilde{\b P}=\b P-\b D$ is a solution to the $\left( \tilde{\b A},\tilde{\b b}\right)$ problem.
\item [{\bf (b)}] Show that $\b A-\b b\b b^{\top}\b P=\tilde{\b A}-\b b\b b^{\top}\tilde{\b P}$.
\end{itemize}
{\bf Hint}: Notice that $\tilde{\b A}=\b A-\b b\b b^{\top}\b D$, and $\tilde{\b b}=\b b$. Write down the algebraic Riccati equation\index{algebraic Riccati equation} for the $\left(\tilde{\b A},\tilde{\b b}\right)$ problem and perform appropriate substitutions.
\end{exercise}

\begin{exercise}\label{ex:brchistochrone}\rm Consider the following variation of the brachistochrone problem:\index{brachistochrone problem!generalized} A particle of mass $m$ is to slide down a frictionless track connecting two given points that are not located on the vertical line. The particle's initial speed is $v_0$. The motion of the particle is constrained only by gravity and the shape of the track. Find the curve of the track that minimizes the time travel between the two given points. Solve this brachistochrone problem using the minimum principle of Pontryagin\index{brachistochrone problem!minimum principle of Pontryagin}. Illustrate your derivations with a numerical example and simulations using MATLAB.
\end{exercise}
\begin{exercise}\label{ex:anim_vehicle}\rm\begin{itemize}\item[{\bf (a)}] Consider the ground vehicle model described in Subsection~\ref{subs:vehmod}. Let the input signal be $u=\delta_f$---the front wheel steer angle, and $\delta_r=0$. Write a MATLAB script that animates the lateral motion $Y$ versus the longitudinal motion $X$ of the vehicle, the steer angle $\delta_f$, and the heading angle $\psi$. Use the steer angle shown in Figure~\ref{fig:steer}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.6]{ster_a.eps}
%\centering \epsfig{file=c:/mybook/ster_a.eps,height=3in}
\caption{Plot of the steer angle $\delta_f$ versus time for Exercise~\ref{ex:anim_vehicle}.}
\label{fig:steer}
\end{figure}
\item[{\bf (b)}] Use the linear model and the data from Table~\ref{tab:par} to design a linear quadratic regulator~(LQR)\index{LQR}. In this part, the control has the form
\[
\delta_f=-\b k\b x+Y_{ref},
\]
where $Y_{ref}$ is the command input shown in Figure~\ref{fig:lqr_in}. That is, the closed-loop system here has the form,
\[
\dot{\b x}=(\b A-\b b\b k)\b x+\b b Y_{ref}.
\]
\begin{figure}
\centering
\includegraphics[scale=0.6]{lqr_in.eps}
%\centering \epsfig{file=c:/mybook/lqr_in.eps,height=3in}
\caption{Plot of the command input $Y_{ref}$ versus time for Exercise~\ref{ex:anim_vehicle}.}
\label{fig:lqr_in}
\end{figure}
You as a designer, will select the weights $\b Q$ and $r$. See how the linear vehicle tracks the command input. Observe the sign difference in the command input and the heading angle. Your LQR gain should be calculated for each chosen $v_x$. In other words, if in your interactive simulation, a specific value of $v_x$ is chosen, then your program should calculate the corresponding gain for this particular value of $v_x$. Implement the LQR state-feedback employing the vehicle model that uses the nonlinear cornering characteristics\index{cornering characteristics} shown in Figure~\ref{fig:dugoffb} and test the performance of the closed-loop system utilizing the command input shown in Figure~\ref{fig:lqr_in}. You may need to impose constraints on $\delta_f$, for example, $|\delta_f |\le 20^{\circ}$. (Do not forget about converting degrees into radians\index{converting!degrees into radians}.) Write a MATLAB script that animates\index{MATLAB animation} the lateral motion $Y$ versus the longitudinal motion\index{longitudinal motion!of vehicle} $X$ of the vehicle, the command input, and the heading angle.\index{heading angle!of vehicle}
\end{itemize}

\end{exercise}

\begin{exercise}\label{ex:one_link_lqr_design}\rm
For the one-link manipulator from Exercise~\ref{exer:one-link_compensat}, use the linearized model from Part~{\bf (b)} of that exercise to
\begin{itemize}
\item[{\bf (a)}] design the optimal controller that minimizes the performance index
\[
J=\frac{1}{2}\Delta \b x(t_f)^{\top}\b F \Delta \b x(t_f) + \frac{1}{2}\int_0^{t_f}\left(\Delta \b x(t)^{\top}\b Q \Delta \b x(t) + \Delta u(t)^2\right) dt,
\]
where  $t_f=2$~sec, $\b F= 3 \b I_3$, and $\b Q= 2 \b I_3$. (Hint: Use the method
presented in Subsection~\ref{subsec:TPBVP});

\item[{\bf (b)}]  Implement and perform animation of the closed-loop system comprised of the optimal state-feedback controller driving the nonlinear model of the robot;

\item[{\bf (c)}]Implement and perform animation of the closed-loop system comprised of the combined optimal controller-observer compensator driving the nonlinear model of the robot.
\end{itemize}


\end{exercise} 