\chapter{Dynamical System Modeling}
\label{chpt:Model_Catalog}
\begin{flushright}
\parbox{5.5in}{To catch a mouse alive, try this trick: When you see the beast running around a room, open a black umbrella and hold it on the floor at an angle that allows the mouse to run in (it helps if someone else can chase the mouse). Your prey, drawn to dark places, will probably head straight for the umbrella. Then quickly close it, go outside, and set your prisoner free.~\cite[page~135]{rd91}}
\end{flushright}

\bigskip

In this chapter, we discuss the Lagrange equations of motion that are useful in mathematical modeling of dynamic systems. We first review relevant concepts that are needed to derive the Lagrange equations. We then construct mathematical models of selected dynamical systems.

\section{Mathematical Modeling}\label{sec:mathmo}

Mathematical models of systems that we wish to control are key elements in their analysis and controller design. We devote this section~\ref{sec:mathmo} to mathematical modeling\index{mathematical modeling}. Constructing models using linguistic description of the process to be controlled will be discussed in Subsection~\ref{subsec:fuzzymodeling}.

The goal of {\em mathematical modeling} is to provide a mathematical description of the inter-relations between the system quantities as well as the relations between these quantities and external inputs. A mathematical model can be viewed as ``a mathematical representation of the significant, relevant aspects of an existing or proposed physical system. (Significance and relevance being in relation to an application where the model is to be used).''~\cite[page~xi]{leigh83}. Another definition of mathematical model that we found in~\cite[page~1]{elworthy94} is: ``A mathematical model is a system of equations that purports to represent some phenomenon in a way that gives insight into the origins and consequences of its behavior.''

At this point, we wish to point out that it is important to distinguish between a physical system\index{physical system} and its model. A physical system may be modeled using different models. In particular, it may be advantageous to model the physical system using different models for its different operating regions. For example, it is common to model an electric  circuit using different models for different ranges of frequencies in which the circuit operates. In what follows, we refer to the physical system's model as the system.

In general, the more accurate the mathematical model is, the more complex are the modeling equations. This means that during the mathematical modeling process we may be forced to find a middle ground between accuracy of the model and its complexity. Our goal is to generate the most accurate possible description of the given system while still being able to carry out the model analysis and computations in a reasonable amount of time. At first, we tend to make oversimplifying assumptions in order to obtain a very simple model. However, we should follow the principle of maximum simplicity\index{principle of maximum simplicity}, which states that "everything should be as simple as possible---but no simpler." After extracting some information about the system behavior based on this first, possibly crude model, it may be desirable to construct a more accurate model. We try to create a mathematical description of a given physical process that incorporates realistic assumptions and physical constraints. We then extract information from this model and compare it with the experimental data. We perform this comparative study to determine the value of the obtained model. If the model appears to be inadequate for our needs, then it must be adjusted. The process of mathematical modeling is summarized in the diagram in Figure~\ref{fig:mat_mod}.
\begin{figure}
\centering
\includegraphics[scale=0.57]{mat_mod.eps}
%\vspace{3.5in}
%\centerbmp{4.66in}{3.5in}{c:/mybook/mat_mod.bmp}
%\hskip1.5in\special{bmp:c:/mybook/mat_mod.bmp x=4.6in y=3.5in}
\caption{Diagram of mathematical modeling process and its use in the controller design process.}
\label{fig:mat_mod}
\end{figure}
In the process of mathematical modeling, one can use different types of equations and their combinations like algebraic, ordinary and partial differential equations, difference, integral, or functional equations. In this book, only difference and ordinary differential equations are used to model dynamical systems.

For the controller design purposes, we first construct a truth model\index{truth model}\index{truth model|see{simulation model}} of the given plant and then its design model\index{design model}. The truth model is the simulation model\index{simulation model} that should include all the relevant characteristics of the physical system to be controlled. The control designs are being simulated using the truth model. A design model is developed by eliminating ``all unpleasant nonlinear effects while capturing all the essential features of the process, is much more amenable to use for the design of a control system''~\cite[p.~26]{friedland96}.


\section{Review of Work and Energy Concepts}
The purpose of this section is to review and reacquaint ourselves with some basic notions from physics like work, energy, etc. After a brief review of work and energy concepts, we derive the Lagrange equations of motion and use them to devise mathematical models of dynamic systems.

Suppose we are given a particle of a constant mass\index{particle of a constant mass} $m$ subjected to a force\index{force} $\mb{\bm$F$}$. Then, by Newton's second law\index{Newton's second law} we have
\begin{equation}
\mb{\bm$F$}=m \frac{d\mb{\bm$v$}}{dt}=m\mb{\bm$a$}.
\label{eq:oldnewt}
\end{equation}
Here the force $\mb{\bm$F$}$ and the velocity $\mb{\bm$v$}$ are vectors. Therefore, they can be represented as
\[
\mb{\bm$F$}=\left[\begin{array}{ccc}
F_{x_1} & F_{x_2} & F_{x_3}
\end{array}\right]^{\top}\quad \mb{and} \quad \mb{\bm$v$}=\left[\begin{array}{ccc}
v_{x_1} & v_{x_2} & v_{x_3}
\end{array}\right]^{\top}=\left[\begin{array}{ccc}
\dot{x}_1 & \dot{x}_2 & \dot{x}_3
\end{array}\right]^{\top}.
\]
Using the above notation, we represent equation~(\ref{eq:oldnewt}) as
\begin{equation}
\left[\begin{array}{c}
F_{x_1}\\
F_{x_2}\\
F_{x_3}
\end{array}\right]=m\left[\begin{array}{c}
\frac{d v_{x_1}}{dt}\\
\frac{d v_{x_2}}{dt}\\
\frac{d v_{x_3}}{dt}
\end{array}\right].
\label{eq:newtinvecn}
\end{equation}
Suppose now that a force $\mb{\bm$F$}$ is acting on a particle located at a point $A$ and the particle moves to a point $B$. The work\index{work} $\delta W$ done by $\mb{\bm$F$}$ along infinitesimally small distance\index{infinitesimally small distance} $\delta \mb{\bm$s$}$ is
\begin{eqnarray*}
\delta W &=& \mb{\bm$F$}^{\top} \delta \mb{\bm$s$}\\
&=& \left[\begin{array}{ccc}
F_{x_1} & F_{x_2} & F_{x_3}
\end{array}\right] \left[\begin{array}{c}
\delta x_1\\
\delta x_2\\
\delta x_3
\end{array}\right]\\
&=& F_{x_1} \delta x_1 + F_{x_2} \delta x_2 + F_{x_3} \delta x_3.
\end{eqnarray*}
where $\delta \mb{\bm$s$}=\left[\begin{array}{ccc}
\delta x_1 & \delta x_2 & \delta x_3
\end{array} \right]^{\top}$. The work $W_{AB}$ done on the path from $A$ to $B$ is obtained by integrating the above equation,
\[
W_{AB}=\int ^{B}_{A}\left(F_{x_1}dx_1 + F_{x_2}dx_2 + F_{x_3}dx_3 \right).
\]
We now would like to establish a relation between work and kinetic energy\index{kinetic energy}. For this observe that
\[
\ddot{x}\,dx=\frac{d \dot{x}}{dt}dx=\dot{x} d\dot{x}.
\]
Using the above relation and Newton's second law, we express $W_{AB}$ in a different way as
\begin{eqnarray}
W_{AB} &=& \int ^{B}_{A}\left(m\ddot{x}_1dx_1 + m\ddot{x}_2dx_2 + m\ddot{x}_3 dx_3 \right)\nonumber\\
&=& \int^{B}_{A}m\left(\dot{x}_1d\dot{x}_1 + \dot{x}_2d\dot{x}_2 + \dot{x}_3d\dot{x}_3 \right)\nonumber\\
&=& m\left.\left(\frac{\dot{x}_1^2}{2}+\frac{\dot{x}_2^2}{2}+\frac{\dot{x}_3^2}{2} \right)\right\arrowvert^{B}_{A}\nonumber\\
&=& m\left(\frac{\| \mb{\bm$v$}_{B} \|^2}{2}-\frac{\| \mb{\bm$v$}_{A} \|^2}{2} \right)\nonumber\\
&=& \frac{m\| \mb{\bm$v$}_{B}\|^2}{2}-\frac{m\| \mb{\bm$v$}_{A}\|^2}{2}.
\label{eq:wablen}
\end{eqnarray}
where $\mb{\bm$v$}_{A}=\left[\begin{array}{ccc}
\dot{x}_{1A} & \dot{x}_{2A} & \dot{x}_{3A}
\end{array}\right]^{\top}$, and $\mb{\bm$v$}_{B}$ is defined in an analogous manner, and $\| \b v\|^2=\b v^{\top}\b v$. Observing that $\frac{m\| \mb{\bm$v$}_{B}\|^2}{2}$ is the {\em kinetic energy} of the particle at the point $B$ and $\frac{m\| \mb{\bm$v$}_{A}\|^2}{2}$ is its kinetic energy at the point $A$, we obtain
\[
W_{AB}=K_{B}-K_{A}.
\]
The above relation can be used to define the kinetic energy of a particle as the work required to change its velocity from some value $\mb{\bm$v$}_{A}$  to a final value $\mb{\bm$v$}_{B}$. The relation, $W_{AB}=K_{B}-K_{A}=\Delta K$ is also known as the {\em work-energy theorem for a particle\/}\index{work-energy theorem}. A force is conservative if the work done by the force on a particle that moves through any round trip is zero. In other words, a force is conservative if the work done by it on a particle that moves between two points depends only on these points and not on the path followed.


We next review the notion of {\em potential energy\/}\index{potential energy} or what is the same, the {\em energy of configuration}. Recall that if the kinetic energy $K$ of a particle changes by $\Delta K$, then the potential energy $U$ must change by an equal but opposite amount so that the sum
of the two changes is zero, that is,
\[
\Delta K + \Delta U = 0.
\]
This is equivalent to saying that any change in the kinetic energy of the particle is compensated for by an equal but opposite change in the potential energy $U$ of the particle
so that their sum remains constant, that is,
\[
K + U = \mb{constant}.
\]
See Exercise~\ref{ex:kplusuconst} for a rigorous proof of the above relation. The potential energy of a particle represents a form of stored energy that can be recovered and converted into the kinetic energy. If we now use the work-energy theorem\index{work-energy theorem} and the equation $\Delta K + \Delta U =0$, then we obtain
\[
W=\Delta K = -\Delta U.
\]
The work done by a conservative force depends only on the starting and the end points of motion and not on the path followed between them. Therefore, for motion in one dimension,
\begin{equation}
\Delta U =-W =- \int_{x_0}^{x}F(s)ds.
\label{eq:onedimdelu}
\end{equation}
We thus can write
\begin{equation}
F(x)=-\frac{dU(x)}{dx}.
\label{eq:fisnegu}
\end{equation}
A generalization of equation~(\ref{eq:onedimdelu})  to motion in three dimensions yields
\[
\Delta U = -\int_{x_{10}}^{x_1} F_{x_1}ds -\int_{x_{20}}^{x_2} F_{x_2}ds -\int_{x_{30}}^{x_3} F_{x_3}ds,
\]
and a generalization of equation~(\ref{eq:fisnegu}) to motion in three dimensions is
\begin{equation}
\mb{\bm$F$}\left( \mb{\bm$x$} \right) = -\left[\begin{array}{ccc}
\frac{\partial U}{\partial x_1} & \frac{\partial U}{\partial x_2} & \frac{\partial U}{\partial x_3}
\end{array}\right]^{\top}
= -\nabla U \left( \mb{\bm$x$} \right),
\label{eq:clerdoconf}
\end{equation}
where $\mb{\bm$x$}= \left[\begin{array}{ccc}
x_1 & x_2 & x_3
\end{array}\right]^{\top}$.


For a formal proof of the property of a conservative vector field that the work done by it on a particle that moves between two points depends only on these points and not on the path followed, see Exercise~\ref{ex:potenergy}.

We next represent Newton's equation\index{Newton's equation} in an equivalent format that establishes a connection with the Lagrange equations\index{Lagrange equations} of motion to be discussed in the following section. To proceed, note that
\begin{equation}
K=\frac{m\left\Arrowvert\dot{\b x}\right\Arrowvert^2}{2}=\frac{m\dot{\b x}^{\top}\dot{\b x}}{2}.
\label{eq:kineneqrep}
\end{equation}
Hence,
\begin{equation}
\frac{\partial K}{\partial \dot{x}_i}=m\dot{x}_i,\quad i=1,2,3.
\label{eq:pderoken}
\end{equation}
Combining~(\ref{eq:clerdoconf}) and~(\ref{eq:pderoken}), we represent Newton's equations as
\begin{equation}
\frac{d}{dt}\left( \frac{\partial K}{\partial \dot{x}_i}\right)+(\nabla U)_i=0,\quad i=1,2,3.
\label{eq:fstefen}
\end{equation}
Let
\begin{equation}
L=K-U.
\label{eq:inflagrn}
\end{equation}
The function $L$ defined above is called the {\em Lagrangian function\/}\index{Lagrangian function} or just the Lagrangian\index{Lagrangian}. Note that
\begin{equation}
\frac{\partial L}{\partial \dot{x}_i}=\frac{\partial K}{\partial \dot{x}_i},
\label{eq:doflagwxd}
\end{equation}
and
\begin{equation}
\frac{\partial L}{\partial {x}_i}=-\frac{\partial U}{\partial {x}_i},
\label{eq:doflagwxon}
\end{equation}
Substituting~(\ref{eq:doflagwxd}) and~(\ref{eq:doflagwxon}) into~(\ref{eq:fstefen}), we obtain
\begin{equation}
\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{x}_i}\right)-\frac{\partial L}{\partial x_i}=0,\quad i=1,2,3.
\label{eq:newtlagbr}
\end{equation}
Equations~(\ref{eq:newtlagbr}) are called the Lagrange equations of motion, in Cartesian coordinates\index{Cartesian coordinates}, for a single particle. They are just an equivalent representation of Newton's equations given by~(\ref{eq:newtinvecn}).

\section{The Lagrange Equations of Motion in Generalized Coordinates}\label{lagrange_eqns}
A single particle with no constraints is free to take up any position in the three dimensional space. This is the reason we need three coordinates to describe its position. In the case of a system of $p$ particles\index{system of $p$ particles} there are $3p$ degrees of freedom\index{degree of freedom}. In general, if there are constraints imposed on the particles, then the number of the degrees of freedom is
\[
n=3p - \mb{number of constraints}.
\]
For example, if the system consists of only one particle in space and if the particle is restricted to move along a line, then $n = 3 - 2 = 1$. We can use many different coordinate systems to describe a location of a particle in space. We use the letter $q$ as a symbol for a coordinate regardless of its nature. A generalized coordinate\index{generalized coordinate} may represent, for example, angular displacement\index{angular displacement}\index{generalized coordinate!angular displacement} or electric charge\index{electric charge}\index{generalized coordinate!electric charge}. This indicates that the results we are going to obtain are of general nature and are not be restricted to mechanical systems. We refer to $q$ as a generalized coordinate. Thus, if
\[
\left.\begin{array}{lll}
x &=& x(q_1,q_2,q_3)\\
y &=& y(q_1,q_2,q_3)\\
z &=& z(q_1,q_2,q_3)
\end{array}\right\}
\]
then
\[
\dot{x}=\frac{\partial x}{\partial q_1}\dot{q}_1 + \frac{\partial x}{\partial q_2}\dot{q}_2 + \frac{\partial x}{\partial q_3}\dot{q}_3,
\]
and so on. We assume that the functions $x, y, z$ are twice continuously differentiable, that is, $x, y, z \in {\cal C}^2$, which implies that their corresponding mixed second partial derivatives are equal.

Suppose now that a particle is confined to move on a surface. Therefore, to describe the location of the particle, we need only two generalized coordinates. The position of the
particle on a surface in a three-dimensional space is then described as
\begin{eqnarray*}
x &=& x(q_1,q_2),\\
y &=& y(q_1,q_2),\\
z &=& z(q_1,q_2),
\end{eqnarray*}
and hence
\begin{equation}
\left.\begin{array}{lll}
\delta x &=& \frac{\partial x}{\partial q_1} \delta q_1 + \frac{\partial x}{\partial q_2} \delta q_2\\
\delta y &=& \frac{\partial y}{\partial q_1} \delta q_1 + \frac{\partial y}{\partial q_2} \delta q_2\\
\delta z &=& \frac{\partial z}{\partial q_1} \delta q_1 + \frac{\partial z}{\partial q_2} \delta q_2.
\end{array}\right\}
\label{eq:dxdydzdq}
\end{equation}
We know that the work $\delta W$ done by the force $\mb{\bm$F$}$ is $\delta W=F_{x}\delta x + F_{y} \delta y + F_{z} \delta z$. On the other hand, this work is equal to the change in the kinetic energy. Therefore,
\begin{equation}
\fbox{$\displaystyle m\left( \ddot{x}\delta x +\ddot{y} \delta y + \ddot{z} \delta z\right)= F_{x}\delta x + F_{y} \delta y + F_{z} \delta z $}
\label{eq:dalemertse}
\end{equation}
The above equation is known as {\em D'Alembert's\/} equation\index{D'Alembert's equation}.

If we have a system of particles acted upon by forces $\mb{\bm$F$}_1, \ldots, \mb{\bm$F$}_p$ and if each of the particles undergoes displacements $\delta \mb{\bm$s$}_1, \ldots, \delta \mb{\bm$s$}_p$, then the total work done by the forces is
\begin{equation}
\delta W=\sum_{i=1}^{p}\left( F_{{x}_{i}} \delta x_{i} + F_{{y}_{i}} \delta y_{i} + F_{{z}_{i}} \delta z_{i} \right).
\label{eq:totwdobyf}
\end{equation}
Substituting~(\ref{eq:dxdydzdq}) into D'Alembert's equation and rearranging, we obtain
\begin{eqnarray}
\delta W &=& m \left(\ddot{x}\frac{\partial x}{\partial q_1}+\ddot{y}\frac{\partial y}{\partial q_1}+\ddot{z}\frac{\partial z}{\partial q_1}\right) \delta q_1+m \left(\ddot{x}\frac{\partial x}{\partial q_2}+\ddot{y}\frac{\partial y}{\partial q_2}+\ddot{z}\frac{\partial z}{\partial q_2}\right)\delta q_2\nonumber\\
&=& \left(F_{x}\frac{\partial x}{\partial q_1} + F_{y}\frac{\partial y}{\partial q_1} + F_{z}\frac{\partial z}{\partial q_1}\right)\delta q_1 + \left( F_{x}\frac{\partial x}{\partial q_2} + F_{y}\frac{\partial y}{\partial q_2} + F_{z}\frac{\partial z}{\partial q_2}\right)\delta q_2.
\label{eq:afsuinaaaleme}
\end{eqnarray}
Let
\[
F_{q_{1}}=F_{x}\frac{\partial x}{\partial q_1} + F_{y}\frac{\partial y}{\partial q_1} + F_{z}\frac{\partial z}{\partial q_1},
\]
and
\[
F_{q_{2}}=F_{x}\frac{\partial x}{\partial q_2} + F_{y}\frac{\partial y}{\partial q_2} + F_{z}\frac{\partial z}{\partial q_2}.
\]
Then, (\ref{eq:afsuinaaaleme}) can be represented as,
\[
\delta W=F_{q_1}\delta q_1 + F_{q_2}\delta q_2.
\]
To proceed further, note that
\[
\frac{d}{dt}\left(\dot{x}\frac{\partial x}{\partial q_1}\right)=\ddot{x}\frac{\partial x}{\partial q_1} + \dot{x}\frac{d}{dt}\left(\frac{\partial x}{\partial q_1} \right).
\]
Hence
\begin{equation}
\fbox{$\displaystyle \ddot{x}\frac{\partial x}{\partial q_1}=\frac{d}{dt}\left(\dot{x}\frac{\partial x}{\partial q_1}\right)-\dot{x}\frac{d}{dt}\left(\frac{\partial x}{\partial q_1} \right) $}
\label{eq:oneofthbugb}
\end{equation}
The time derivative of $x=x(q_1,q_2)$ is
\begin{equation}
\dot{x}=\frac{\partial x}{\partial q_1}\dot{q}_1 + \frac{\partial x}{\partial q_2} \dot{q}_2.
\label{eq:timderofxb}
\end{equation}
Differentiating the above partially with respect to $\dot{q}_1$ yields
\begin{equation}
\fbox{$\displaystyle \frac{\partial \dot{x}}{\partial \dot{q}_1} = \frac{\partial x}{\partial q_1} $}
\label{eq:prtdeoxdwq}
\end{equation}
Finally, we show that
\begin{equation}
\fbox{$\displaystyle \frac{d}{dt} \left( \frac{\partial x}{\partial q_1} \right) = \frac{\partial \dot{x}}{\partial q_1} $}
\label{eq:tlontoleqm}
\end{equation}
Indeed, because $x=x(q_1,q_2)$ and the partial derivative ${\partial x}/{\partial q_1}$ is in general a function of $q_1$ and $q_2$, we have
\begin{equation}
\frac{d}{dt}\left(\frac{\partial x}{\partial q_1} \right) = \frac{\partial}{\partial q_1} \left(\frac{\partial x}{\partial q_1} \right) \dot{q}_1 + \frac{\partial}{\partial q_2} \left( \frac{\partial x}{\partial q_1} \right) \dot{q}_2.
\label{eq:ddtofxwrtqo}
\end{equation}
Taking the partial derivative of~(\ref{eq:timderofxb}) with respect to $q_1$, we get
\begin{eqnarray}
\frac{\partial \dot{x}}{\partial q_1} &=& \frac{\partial}{\partial q_1} \left(\frac{\partial x}{\partial q_1} \right) \dot{q}_1 + \frac{\partial}{\partial q_1} \left( \frac{\partial x}{\partial q_2} \right) \dot{q}_2\nonumber\\
&=& \frac{\partial}{\partial q_1} \left(\frac{\partial x}{\partial q_1} \right) \dot{q}_1 + \frac{\partial}{\partial q_2} \left( \frac{\partial x}{\partial q_1} \right) \dot{q}_2,
\label{eq:mixeddrhelp}
\end{eqnarray}
because by assumption $x,y,z \in {\cal C}^2$ and hence
\[
\frac{\partial}{\partial q_1} \left( \frac{\partial x}{\partial q_2} \right)=\frac{\partial}{\partial q_2} \left( \frac{\partial x}{\partial q_1} \right).
\]
With the above in mind, comparing the right hand sides of~(\ref{eq:ddtofxwrtqo}) and~(\ref{eq:mixeddrhelp}) yields~(\ref{eq:tlontoleqm}). We now use~(\ref{eq:oneofthbugb}), (\ref{eq:prtdeoxdwq}), and~(\ref{eq:tlontoleqm}) to obtain the Lagrange equation of motion for a single particle. We first substitute~(\ref{eq:prtdeoxdwq}) and~(\ref{eq:tlontoleqm}) into~(\ref{eq:oneofthbugb}) to get
\begin{equation}
\ddot{x}\frac{\partial x}{\partial q_1}=\frac{d}{dt}\left(\dot{x}\frac{\partial \dot{x}}{\partial \dot{q}_1}\right) -\dot{x}\frac{\partial \dot{x}}{\partial q_1}.
\label{eq:aftsubthe}
\end{equation}
Note that
\begin{equation}
\frac{\partial}{\partial \dot{q}_1} \left(\frac{\dot{x}^2}{2} \right) = \dot{x} \frac{\partial \dot{x}}{\partial \dot{q}_1},
\label{eq:thfiforab}
\end{equation}
and
\begin{equation}
\frac{\partial}{\partial q_1} \left(\frac{\dot{x}^2}{2} \right) = \dot{x} \frac{\partial \dot{x}}{\partial q_1}.
\label{eq:thesecforab}
\end{equation}
Substituting~(\ref{eq:thfiforab}) and~(\ref{eq:thesecforab}) into~(\ref{eq:aftsubthe}) gives
\[
\ddot{x}\frac{\partial x}{\partial q_1} = \frac{d}{dt}\left( \frac{\partial \left( \frac{\dot{x}^2}{2} \right)}{\partial \dot{q}_1} \right) - \frac{\partial \left( \frac{\dot{x}^2}{2} \right)}{\partial q_1}.
\]
We can obtain similar expressions for $y$ and $z$. Taking this into account, the expression~(\ref{eq:afsuinaaaleme}) when $\delta q_2=0$ becomes,
\begin{eqnarray}
\delta W &=&\left( \frac{d}{dt} \left(\frac{\partial}{\partial \dot{q}_1} \frac{m\left(\dot{x}^2 + \dot{y}^2 + \dot{z}^2 \right)}{2}\right)\right. - \left.\frac{\partial}{\partial q_1} \frac{m \left( \dot{x}^2 + \dot{y}^2 + \dot{z}^2 \right)}{2}  \right) \delta q_1\nonumber\\
&=& \left( F_{x} \frac{\partial x}{\partial q_1} + F_{y} \frac{\partial y}{\partial q_1}+F_{z} \frac{\partial z}{\partial q_1}\right) \delta q_1\nonumber\\
&=&F_{q_1} \delta q_1.
\label{eq:lexffqo}
\end{eqnarray}
Let $K=\frac{1}{2}m \left(\dot{x}^2 + \dot{y}^2 + \dot{z}^2 \right)$ denote the kinetic energy of the particle. Then, we represent~(\ref{eq:lexffqo}) as
\begin{equation}
\fbox{$ \displaystyle \frac{d}{dt} \left( \frac{\partial K}{\partial \dot{q}_1} \right) - \frac{\partial K}{\partial q_1} = F_{q_1} $}
\label{eq:lagfqocor}
\end{equation}
This equation is called the {\em Lagrange equation of motion}\index{Lagrange equation of motion} for the $q_1$ coordinate. Using the same arguments as above, we can derive the Lagrange equation of motion for the $q_2$ coordinate. In general there are as many the Lagrange equations of motion as there are degrees of freedom of the particle.

Recall that the expression $\delta W=F_{x} \delta x + F_{y} \delta y + F_{z} \delta z$ is a work done by a force $\mb{\bm$F$}=\left[\begin{array}{ccc}
F_{x} & F_{y} & F_{z}
\end{array} \right]^{\top}$ for a general displacement $\delta \mb{\bm$s$}=\left[\begin{array}{ccc}
\delta x & \delta y & \delta z
\end{array} \right]^{\top}$ . In terms of generalized coordinates and when, say, $\delta q_2 =0$,  we have
\begin{eqnarray*}
\delta W &=& \left( F_{x} \frac{\partial x}{\partial q_1} + F_{y} \frac{\partial y}{\partial q_1}+F_{z} \frac{\partial z}{\partial q_1}\right) \delta q_1\\
&=& F_{q_1} \delta q_1\\
&=& \delta W_{q_1}.
\end{eqnarray*}
A generalized force\index{generalized force} $F_{q_{r}}$ is of such nature that the product $F_{q_{r}} \delta q_{r}$ is the work
done by driving forces when $q_{r}$  alone is changed by $\delta q_{r}$. We mention that a generalized force does not have to be a force in the usual sense. For example, if $q_{r}$ is an angle, then $F_{q_{r}}$ must be a torque in order that $F_{q_{r}} \delta q_{r}$ be work.

We now derive a more general version of the Lagrange equations of motion for a single particle. In our derivation we need the notion of {\em nonconservative\/} forces\index{nonconservative forces} that we discuss next. First recall that we interpret the kinetic energy of a particle as its ability to do work by virtue of its motion. This motion is the result of conservative and nonconservative forces acting upon the particle. A force is nonconservative if the work done by the force on a particle that moves through any round trip is not zero. Thus a force is nonconservative if the work done by that force on a particle that moves between two points depends on the path taken between these points.

Let us now suppose that in addition to the conservative forces, a single nonconservative force due to friction\index{friction} acts on the particle. Denote by $\sum W_{c}$ the sum of the work done by the conservative forces, and by $W_{f}$ the work done by friction. We then have
\begin{equation}
W_{f}+\sum W_{c}= \Delta K.
\label{eq:nonconser}
\end{equation}
Thus, if a nonconservative force, like friction, acts on a particle, then the total mechanical energy\index{total mechanical energy} is not constant, but changes by the amount of work done by a nonconservative force according to the relation~(\ref{eq:nonconser}). Let us denote by $E$ the final mechanical energy\index{mechanical energy} of the particle and by $E_0$ its initial mechanical energy. Then, we can write
\[
E-E_0=\Delta E = W_{f}.
\]
The work done by friction, the nonconservative force, is always negative. Hence, the final mechanical energy $E$ is less than the initial mechanical energy. The ``lost'' mechanical energy\index{lost mechanical energy} is transformed into {\em internal energy\/}\index{internal energy} $U_{int}$. In other words, the loss of mechanical energy equals the gain in internal energy and we write
\[
\Delta E + U_{int} = 0.
\]
In general, we have
\[
\Delta K + \sum \Delta U + U_{int} + \left(\mb{change in other forms of energy}\right)=0.
\]
This means that the total energy, kinetic plus potential plus internal plus all other forms of energy, does not change. This is the {\em principle of the conservation of energy}\index{principle of the conservation of energy}. To proceed further, recall that any generalized force $F_{q_i}$ acting on our particle can be expressed as
\[
F_{q_i}=F_{x}\frac{\partial x}{\partial q_i}+F_{y}\frac{\partial y}{\partial q_i} + F_{z}\frac{\partial z}{\partial q_i}.
\]
Assuming that the forces are conservative, we get
\[
F_{q_i}=-\left( \frac{\partial U}{\partial x}\frac{\partial x}{\partial q_i} + \frac{\partial U}{\partial y}\frac{\partial y}{\partial q_i} + \frac{\partial U}{\partial z}\frac{\partial z}{\partial q_i}\right).
\]
We thus have
\[
F_{q_i}=-\frac{\partial U}{\partial q_i}.
\]
Therefore, the Lagrange equations of motion for a single particle for conservative forces can be represented as
\[
\frac{d}{dt}\left(\frac{\partial K}{\partial \dot{q}_i} \right) - \frac{\partial K}{\partial q_i}=-\frac{\partial U}{\partial q_i}.
\]
We can rearrange the above equations as
\[
\frac{d}{dt}\left(\frac{\partial K}{\partial \dot{q}_i} \right) - \frac{\partial}{\partial q_i}\left( K-U \right)=0.
\]
We now introduce the {\em Lagrangian function\/} $L$, defined as
\[
L=K-U.
\]
Note that because $U$ is not a function of $\dot{q}_i$, we have
\[
\frac{\partial K}{\partial \dot{q}_i}=\frac{\partial L}{\partial \dot{q}_i}.
\]
Taking the above into account, we can represent the Lagrange equations of motion for conservative forces as
\[
\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = 0.
\]
If some of the forces acting on the particle are nonconservative, the Lagrange equations of motion for our particle will take the form
\begin{equation}
\fbox{$\displaystyle \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = \tilde{F}_{q_i} $}
\label{eq:finfoflaeq}
\end{equation}
where $\tilde{F}_{q_i}$ are nonconservative forces\index{nonconservative forces} acting on the particle.
\begin{example}\label{ex:simP_pend_Lagrange}\rm
Consider the simple pendulum shown in Figure~\ref{fig:sim_ped}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{sim_ped.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/sim_ped.bmp}
%\hskip1in\special{bmp:c:/mybook/sim_ped.bmp x=4.6in y=3.5in}
\caption{A simple pendulum.}
\label{fig:sim_ped}
\end{figure}
The simple pendulum\index{simple pendulum} is an idealized body consisting of a point mass $M$, suspended by a weightless inextensible cord of length $l$. The simple pendulum is an example of a one degree of freedom system with a generalized coordinate being the angular displacement $\theta$. The pendulum kinetic energy is
\[
K=\frac{1}{2}Ml^2\dot{\theta}^2.
\]
We assume that the potential energy of the pendulum is zero when the pendulum is at rest. Hence, its potential energy is
\[
U=Mgl ( 1 - \cos \theta ),
\]
where $g=10\,\mb{m/sec}^2$ is the acceleration due to gravity. The Lagrangian function is
\begin{eqnarray*}
L &=& K-U\\
&=& \frac{1}{2}Ml^2 \dot{\theta}^2 - Mgl ( 1 - \cos \theta ).
%&=& 6 \dot{\theta}^2 - 60(1 - \cos \theta ).
\end{eqnarray*}
We compute
\[
\frac{\partial L}{\partial \theta}=-Mgl \sin\theta
\]
and
\[
\frac{\partial L}{\partial \dot{\theta}} = Ml^2 \dot{\theta}.
\]
Hence
\[
\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{\theta}} \right) = Ml^2 \ddot{\theta}.
\]
Combining the above terms together we obtain the Lagrange equation describing the pendulum motion,
\begin{eqnarray*}
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}} \right) - \frac{\partial L}{\partial \theta} &=&  Ml^2 \ddot{\theta} + Mgl \sin\theta\\
&=& 0.
\end{eqnarray*}
Solving for $\ddot{\theta}$, we obtain
\begin{equation}
\ddot{\theta}=-\frac{g}{l}\sin\theta.
\label{eq:gen_fo_sim_pen}
\end{equation}
In our example the above equation evaluates to
\begin{equation}
\ddot{\theta}=-5 \sin \theta.
\label{eq:sec_ord_mode_simp_pen}
\end{equation}
\end{example}
\markex
\begin{example}\rm
\label{ex:thaetarman}
In Figure~\ref{fig_theta_r}, the $\theta$-$r$ robot manipulator is shown. This example was adapted from Snyder~\cite[Chapter~10]{snyder85}.
\begin{figure}
\centering
\includegraphics[scale=0.35]{theta_r.eps}
%\vspace{3in}
%\centerbmp{4in}{3in}{c:/ee675/theta_r.bmp}
%\hskip1in \special{bmp:c:/mybook/theta_r.bmp x=4in y=3in}
\caption{The $\theta$-$r$ robot manipulator.}
\label{fig_theta_r}
\end{figure}
A lumped mass representation\index{lumped mass representation} of the $\theta$-$r$ robot manipulator\index{robot manipulator}\index{robot manipulator!lumped mass representation} is shown in Figure~\ref{fig_t_r_sche}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{t_r_sche.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/ee675/t_r_sche.bmp}
%\hskip.5in \special{bmp:c:/ee675/t_r_sche.bmp x=4.67in y=3.5in}
\caption{A lumped mass representation of the $\theta$-$r$ robot manipulator.}
\label{fig_t_r_sche}
\end{figure}
The mass $m_1=10\,\mb{kg}$ represents the mass of the outer cylinder positioned at its center of mass. The constant distance $r_1=1\mb{m}$ designates the fixed distance between the center of the mass of the outer cylinder and the center of rotation. The mass of the load is represented by $m_2=3\mb{kg}$ and is assumed to be located at the end of a piston of a telescoping arm that is a variable radial distance $r$ from the hub or center of rotation. The angle of rotation of the manipulator arm is denoted by $\theta$. The inputs to the system are assumed to be a torque $T_{\theta}$ applied at the hub in the direction of $\theta$, and a translational force $F_r$ applied in the direction $r$. Finally, $g=10\,\mb{m}/\mb{s}^2$ is the gravitational constant. We construct the Lagrange equations of motion for the system, and then represent the obtained model in state space format. We first find the total kinetic energy of the system. The position of the mass $m_1$ is
\[
\left.\begin{array}{lll}
x_1 &=& r_1\cos\theta\\
y_1 &=& r_1\sin\theta.
\end{array}\right\}
\]
Note that $r_1$ is a constant. Differentiating $x_1$ and $y_1$ with respect to time yields
\[
\left.\begin{array}{lll}
\dot{x}_1 &=& -r_1 \dot{\theta}\sin \theta\\
\dot{y}_1 &=& r_1\dot{\theta}\cos\theta.
\end{array}\right\}
\]
The magnitude squared of the velocity vector of the mass $m_1$ is
\begin{eqnarray*}
v_1^2 &=& \dot{x}_1^2 + \dot{y}_1^2\\
&=& r_1^2\dot{\theta}^2\sin^2\theta+r_1^2\dot{\theta}^2\cos^2\theta\\
&=& r_1^2\dot{\theta}^2.
\end{eqnarray*}
Therefore, the kinetic energy of the mass $m_1$ is
\[
K_1=\frac{1}{2}m_1v_1^2=\frac{1}{2}m_1r_1^2\dot{\theta}^2.
\]
We now derive an expression for the kinetic energy of the second mass. The position of the mass $m_2$ is
\[
\left.\begin{array}{lll}
x_2 &=& r\cos\theta\\
y_2 &=& r \sin \theta.
\end{array}\right\}
\]
Note that the distance $r$ is not a constant. Differentiating $x_2$ and $y_2$ with respect to time yields
\[
\left.\begin{array}{lll}
\dot{x}_2 &=& \dot{r}\cos \theta - r\dot{\theta}\sin \theta\\
\dot{y}_2 &=& \dot{r}\sin \theta + r\dot{\theta}\cos \theta.
\end{array}\right\}
\]
The magnitude squared of the velocity vector of the mass $m_2$ is
\begin{eqnarray*}
v_2^2 &=& \left( \dot{x}_2\right)^2+ \left( \dot{y}_2\right)^2\\
&=& \left( \dot{r}\cos \theta - r\dot{\theta}\sin \theta\right)^2+\left(\dot{r}\sin \theta + r\dot{\theta}\cos \theta \right)^2\\
&=& \dot{r}^2+r^2\dot{\theta}^2.
\end{eqnarray*}
Hence, the kinetic energy of the mass $m_2$ is
\[
K_2=\frac{1}{2}m_2v_2^2=\frac{1}{2}m_2\left(\dot{r}^2+r^2\dot{\theta}^2 \right).
\]
The total kinetic energy of the system is
\[
K=K_1+K_2=\frac{1}{2}m_1r_1^2\dot{\theta}^2+\frac{1}{2}m_2\left(\dot{r}^2+r^2\dot{\theta}^2 \right).
\]
We now find the potential energy of the system. The potential energy of the mass $m_1$ is
\[
U_1=m_1gr_1\sin \theta.
\]
The potential energy of the mass $m_2$ is
\[
U_2=m_2gr\sin \theta.
\]
The total potential energy of the system is
\[
U=U_1+U_2=m_1gr_1\sin \theta+m_2gr\sin \theta.
\]
The Lagrangian function for the system is
\[
L=K-U=\frac{1}{2}m_1r_1^2\dot{\theta}^2+\frac{1}{2}m_2\dot{r}^2+\frac{1}{2}m_2 r^2\dot{\theta}^2-m_1gr_1\sin \theta-m_2gr\sin \theta.
\]
The manipulator has two degrees of freedom. Hence, we have two Lagrange's equations of motion for this system,
\[
\left.\begin{array}{lll}
\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{\theta}}\right)-\frac{\partial L}{\partial \theta} &=& T_{\theta}\\
\frac{d}{dt}\left( \frac{\partial L}{\partial \dot{r}}\right)-\frac{\partial L}{\partial r} &=& F_{r}.
\end{array}\right\}
\]
The first Lagrange equation yields
\[
m_1r_1^2\ddot{\theta}+m_2r^2\ddot{\theta}+2m_2r\dot{r}\dot{\theta}+g\cos\theta(m_1r_1+m_2r)=T_{\theta}.
\]
The second Lagrange equation gives
\[
m_2\ddot{r}-m_2r\dot{\theta}^2+m_2g\sin \theta=F_r.
\]
We next represent the above modeling equations in state space format. (One can use, for example, MATLAB's Symbolic Math Toolbox to perform needed calculations.)
% We use the commands
%\begin{verbatim}
%M=sym('[m1*r1^2+m2*z3^2,0;0,m2]');
%v=sym('[-2*m2*z3*z4*z2-g*cos(z1)*(m1*r1+m2*z3)+u1;
%m2*z3*z2^2-m2*g*sin(z1)+u2]');
%Minv=inverse(M);
%g=symmul(Minv,v);
%simplify(g);
%pretty(ans)
%\end{verbatim}
%Then, we use the obtained results to represent the model in state space %format.
Let the state variables be defined as:
\[
z_1=\theta,\quad z_2=\dot{\theta},\quad z_3=r,\quad z_4=\dot{r}.
\]
Then, we can write
\begin{equation}
\left[\begin{array}{c}
\dot{z}_1\\
\dot{z}_2\\
\dot{z}_3\\
\dot{z}_4
\end{array}\right]=\left[\begin{array}{c}
z_2\\
\frac{-2m_2z_2z_3z_4-g\cos z_1(m_1r_1+m_2z_3)+u_1}{m_1r_1^2+m_2z_3^2}\\
z_4\\
z_2^2z_3-g\sin z_1+u_2/m_2
\end{array}\right].
\label{eq:stamothr}
\end{equation}
\end{example}
\markex
In the following section, we develop models of a few more dynamical systems using basic laws of physics.

\section{Model Catalog}\label{modeling_exampls}
\subsection{Centrifugal Governor}
\label{subsec:watt}
In this subsection, we model the centrifugal governor\index{centrifugal governor}\index{centrifugal governor|see{Watt's governor}}. In our modeling, we follow Pontryagin~\cite[pages 213--220]{pontryagin62}. For an alternative approach to modeling the governor see Maxwell~\cite{maxwell68} and Fuller~\cite{fuller96}. Some, like for example MacFarlane~\cite[page~250]{macfarlane79}, claim that Maxwell founded the theory of automatic control\index{automatic control} systems with his 1868 paper ``On governors''. However, this landmark paper is written in a terse and cryptic style. Fuller~\cite{fuller96} does a wonderful job elucidating, correcting, and expanding Maxwell's paper. Mayr~\cite[Chapter~X]{mayr70} gives a fascinating historical account of the application of the centrifugal governor to the speed regulation of the steam engine and tracks the origins of the flyball governor\index{flyball governor}\index{flyball governor|see{Watt's governor}}.

A schematic of the governor is depicted in Figure~\ref{fig:watts}, while forces acting on the governor are shown in Figure~\ref{fig_watt_s}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.75]{watt_s.eps}
%\vspace{4.5in}
%\centerbmp{6in}{4.5in}{c:/mybook/watt_s.bmp}
%\special{bmp:c:/mybook/watt_s.bmp x=6in y=4.5in}
\caption{Forces acting on the weights in the governor.}
\label{fig_watt_s}
\end{figure}
%For simplicity of notation, we assume that the length of the arms, on which the weights are fastened, are of unity length.
If the weights rotate with angular velocity $\omega_{gov}$, then the arms form an angle $\varphi$ with the vertical position. Let $v_{weight}$ denote the linear velocity of the weights, and $r$ be the radius of the circle along which the weights rotate. Then, $v_{weight}=\omega_{gov}r$. Note that $r=L\sin \varphi$. The resulting centrifugal force\index{centrifugal force} acting on each weight is
\begin{eqnarray*}
m\frac{v_{weight}^2}{r} &=& m\omega_{gov}^2r\\
&=& m\omega_{gov}^2 L\sin \varphi.
\end{eqnarray*}
Simultaneously, a gravitational force $mg$ acts on each weight. In addition, a frictional force in the hinge joints acts on each weight. We assume that the frictional force is proportional to $\dot{\varphi}$ with the friction coefficient $b$. Summing forces along the axis perpendicular to an arm yields an equation of motion of the weight
\begin{eqnarray}
m\ddot{\varphi}L &=& m\omega_{gov}^2 L\sin\varphi\cos\varphi - mg\sin\varphi - b\dot{\varphi}\nonumber\\
&=& \frac{1}{2}m\omega_{gov}^2 L\sin 2\varphi - mg\sin\varphi - b\dot{\varphi}.
\label{eq:eformoofwe}
\end{eqnarray}
%Note that $m\ddot{\varphi}L=m\ddot{\varphi}$ because by assumption $L=1$.
The governor is connected to the engine via a gear with the gear ratio $N$ such that
\begin{equation}
\omega_{gov}=N\omega,
\label{eq:gearratengo}
\end{equation}
where $\omega$ is the angular velocity of the flywheel\index{flywheel}. Substituting~(\ref{eq:gearratengo}) into~(\ref{eq:eformoofwe}) gives
\[
m\ddot{\varphi}L=\frac{1}{2}mN^2\omega^2 L\sin 2\varphi - mg\sin\varphi - b\dot{\varphi}.
\]
Let $I$ be the moment of inertia of the flywheel, $\tau_1$ the torque delivered by the engine, and $\tau_l$ the load torque. Then, the equation modeling the engine is
\begin{equation}
I\dot{\omega}=\tau_1-\tau_l.
\label{eq:emotheeng}
\end{equation}
The collar of the governor\index{collar of the governor} is connected with the valve so that at the nominal value of $\varphi=\varphi_0$ the engine delivers the nominal torque $\tau_d$ that results in the desired angular velocity $\omega=\omega_d$ of the flywheel. We model the interaction of the collar and the valve using the equation
\begin{equation}
\tau_1=\tau_d+\kappa(\cos\varphi-\cos\varphi_0),
\label{eq:collavav}
\end{equation}
where $\kappa>0$ is a proportionality constant that depends on the governor. Substituting~(\ref{eq:collavav}) into~(\ref{eq:emotheeng}) yields
\begin{eqnarray*}
I\dot{\omega} &=& \tau_1-\tau_l\\
&=&\tau_d+\kappa(\cos\varphi-\cos\varphi_0)-\tau_l\\
&=& \kappa\cos\varphi-\tau,
\end{eqnarray*}
where $\tau=\tau_l-\tau_d+\kappa \cos\varphi_0$. We thus have two equations that model the governor,
\begin{eqnarray*}
m\ddot{\varphi}L &=& \frac{1}{2}mN^2\omega^2 L\sin 2\varphi - mg\sin\varphi - b\dot{\varphi},\\
I\dot{\omega} &=& \kappa\cos\varphi-\tau.
\end{eqnarray*}
We use the state variables
\[
x_1=\varphi,\quad x_2=\dot{\varphi},\quad\mb{and}\quad x_3=\omega
\]
to represent the above model in state space format,
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& \frac{1}{2}N^2x_3^2\sin 2x_1-\frac{g}{L}\sin x_1 - \frac{b}{m L}x_2\\
\dot{x}_3 &=& \frac{\kappa}{I}\cos x_1-\frac{\tau}{I}.
\end{array}\right\}
\label{eq:watt_model}
\end{equation}

\subsection{Ground Vehicle}
\label{subs:vehmod}
In this subsection, we model a ground vehicle\index{ground vehicle} in a turning maneuver\index{turning maneuver}\index{turning maneuver|see{ground vehicle}}, such as, a step lane change maneuver\index{step lane change maneuver}. During such a maneuver, the vehicle undergoes translational as well as rotational motion. When modeling a moving vehicle, it is convenient to use reference coordinate frame attached to and moving with the vehicle. This is because with respect to the vehicle coordinate frame moments of inertias of the vehicle are constant. We use the Society of Automotive Engineers~(SAE)\index{Society of Automotive Engineers~(SAE)}\index{SAE|see{Society of Automotive Engineers}} standard coordinate system as shown in Figure~\ref{fig:one}. Our modeling follows that of Wong~\cite{wong93} and Will and \.{Z}ak~\cite{will97}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{sae2.eps}
%\vspace{4in}
%\centerbmp{4.67in}{3.5in}{c:/sliding/sae2.bmp}
%\hskip.2in\special{bmp:c:/sliding/sae2.bmp x=5.33in y=4in}
\caption{The Society of Automotive Engineers (SAE) coordinate system attached to the vehicle.}
\label{fig:one}
\end{figure}
We use lower case letters to label axes of the vehicle coordinate system\index{vehicle coordinate system}. Following the SAE convention, we define the coordinates as follows:
\begin{itemize}
\item the $x$-coordinate is directed forward and on the longitudinal plane of symmetry of the vehicle,
\item the $y$-coordinate has a lateral direction out the right side of the vehicle,
\item the $z$-coordinate is directed downward with respect to the vehicle,
\item the yaw velocity about the $z$-axis is denoted $\dot{\theta}$.
\end{itemize}
In developing our simple model, we neglect the effects of roll\index{roll}index{roll|see{ground vehicle}} and pitch\index{pitch}index{pitch|see{ground vehicle}}. The origin of the vehicle coordinate system is located at the vehicle center of gravity (CG)\index{center of gravity (CG)}\index{CG|see{center of gravity}}. We also use a global, or fixed coordinate system. Its axes are labeled using upper case letters as shown in Figure~\ref{fig:veh_pth}.
\begin{figure}
\centering
\includegraphics[scale=0.55]{VEH_PTH.eps}
%\vspace{3in}
%\hskip1in\special{bmp:c:/sliding/veh_pth.bmp x=4in y=3in}
\caption{Analyzing vehicle motion using coordinate frame attached to the vehicle; snapshots of the vehicle velocity at times $t$ and $t+\Delta t$.}
\label{fig:veh_pth}
\end{figure}
To form the equations of motion of the vehicle, we need an expression for the vehicle CG acceleration. We first compute the change in the longitudinal component of the velocity. Using Figure~\ref{fig:veh_pth} and Figure~\ref{fig:veh_acc}, and assuming small angles, we obtain
\begin{eqnarray*}
\left( v_x+\Delta v_x\right)\cos \Delta \theta -\left( v_y+\Delta v_y\right)\sin \Delta \theta - v_x &\approx& \Delta v_x-v_y \Delta \theta\\
&=& \Delta \dot{x}-\dot{y}\Delta \theta.
\end{eqnarray*}
\begin{figure}
\centering
\includegraphics[scale=0.65]{VEH_ACC.eps}
%\vspace{3in}
%\hskip1in\special{bmp:c:/sliding/veh_acc.bmp x=4in y=3in}
\caption{Computing the change in the longitudinal and lateral components of the vehicle velocity.}
\label{fig:veh_acc}
\end{figure}
We form the Newton's quotient\index{Newton's quotient} using the above. Then, taking the limit, we obtain the total longitudinal acceleration component $a_x$,
\begin{eqnarray}
a_x &=&  \lim_{\Delta t \to 0}\frac{\Delta v_x - v_y \Delta \theta}{\Delta t}\nonumber\\
&=& \dot{v}_x-v_y\dot{\theta}\nonumber\\
&=& \ddot{x}-\dot{y}\omega_z,
\label{eq:axinloncoov}
\end{eqnarray}
where $\omega_z=\dot{\theta}$. Using similar argument to the one above, we obtain the lateral acceleration component $a_y$ to be
\begin{eqnarray}
a_y &=& \dot{v}_y+v_x\dot{\theta}\nonumber\\
&=& \ddot{y}+\dot{x}\omega_z.
\label{eq:axinlaterdir}
\end{eqnarray}
We assume that the vehicle is symmetric with respect to the $x$-$z$ plane. We use the so called bicycle model\index{bicycle model}\index{bicycle model|see{ground vehicle}} that is shown in Figure~\ref{fig:veh_for}, where $\delta_f$ denotes the front wheels' steering angle and $\delta_r$ is the rear wheels' steering angle\index{steering angle}\index{steering angle|see{ground vehicle}}. We denote by $2F_{xf}$ lumped longitudinal force at the front axle and by $2F_{xr}$ lumped longitudinal force\index{longitudinal force} at the rear axle. The lumped lateral force at the front axle is denoted $2F_{yf}$, while lumped lateral force\index{lateral force} at the rear axle is denoted $2F_{yr}$.
\begin{figure}
\centering
\includegraphics[scale=0.6]{VEH_FOR2.eps}
%\vspace{3in}
%\hskip1.2in\special{bmp:c:/sliding/veh_for2.bmp x=4in y=3in}
\caption{A simplified, bicycle model of a vehicle.}
\label{fig:veh_for}
\end{figure}
Summing the forces along the $x$ axis and assuming small angles, we obtain
\begin{eqnarray}
m\left(\ddot{x}-\dot{y}\omega_z \right)&=&2F_{xf}\cos \delta_f-2F_{yf}\sin\delta_f +2F_{xr}\cos \delta_r + 2F_{yr}\sin \delta_r\nonumber\\
&\approx& 2F_{xf}+2F_{xr}.
\label{eq:suoffoinlongd}
\end{eqnarray}
Next, summing the forces along the $y$ axis and assuming small angles yields
\begin{eqnarray}
m\left(\ddot{y}+\dot{x}\omega_z \right)&=&2F_{xf}\sin \delta_f+2F_{yf}\cos\delta_f +2F_{yr}\cos\delta_r-2F_{xr}\sin\delta_r\nonumber\\
&\approx& 2F_{yf}+2F_{yr}.
\label{eq:suoflongdi}
\end{eqnarray}
The third equation of motion is obtained by summing the torques with respect to the $z$-axis. Let $I_z$ denote the moment of inertia of the vehicle about the $z$ axis. Then, referring again to Figure~\ref{fig:veh_for} and assuming small angles, we obtain:
\begin{eqnarray}
I_z\dot{\omega}_z&=&2l_1F_{yf}\cos\delta_f + 2l_1F_{xf}\sin\delta_f-2l_2F_{yr}\cos\delta_r + 2l_2F_{xr}\sin\delta_r\nonumber\\
&\approx& 2l_1F_{yf}-2l_2F_{yr}.
\label{eq:rotnelabom}
\end{eqnarray}
In our further analysis, we assume $v_x=\dot{x}=\mb{constant}$, that is, the driver neither is accelerating nor braking the vehicle. This means that we need to use only the last two equations of motion in our model. To proceed, we need to define tire slip\index{tire slip} angles and cornering stiffness coefficients. The tire slip angle is the difference between the tire's actual travel direction and the tire's heading direction, which is the steer angle. In other words, the tire slip angle is the angle between the desired direction of motion of the tire specified by the steer angle and the actual direction of travel of the center of the tire contact. We denote by $C_{\alpha f}$ the cornering stiffness coefficient\index{cornering stiffness coefficient} of the front tires and by $C_{\alpha r}$ the cornering stiffness coefficient of the rear tires. The coefficients $C_{\alpha f}$ and $C_{\alpha r}$ are measured experimentally. A typical relation between a tire slip angle and a cornering force is shown in Figure~\ref{fig:dugoffb}. It is a nonlinear relation, where the cornering coefficient is
\[
C_{\alpha w}=\frac{d F_{yw}}{d \alpha_w},
\]
where $w=f$ or $w=r$.
\begin{figure}
\centering
\includegraphics[scale=0.7]{DUGOFFB.eps}
%\vspace{3.5in}
%\hskip.9in\special{bmp:c:/mybook/dugoffb.bmp x=4.67in y=3.5in}
\caption{Typical cornering characteristic.}
\label{fig:dugoffb}
\end{figure}
The lateral forces $F_{yf}$ and $F_{yr}$ that appear in~(\ref{eq:suoffoinlongd})--(\ref{eq:rotnelabom}) are functions of the corresponding tire slip angles $\alpha_f$ and $\alpha_r$ and cornering stiffness coefficients $C_{\alpha f}$ and $C_{\alpha r}$. In our simple linearized model, we assume that $C_{\alpha f}=C_{\alpha r} =$~constant, and model dependencies between the lateral forces and tire slip angles as
\begin{equation}
F_{yf}=C_{\alpha f}\alpha_f,
\label{eq:laforacof}
\end{equation}
and
\begin{equation}
F_{yr}=C_{\alpha r}\alpha_r.
\label{eq:relafocor}
\end{equation}
We formulate expressions for the tire slip angles using Figure~\ref{fig:veh_vel}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{VEH_VEL2.eps}
%\vspace{3in}
%\hskip1in\special{bmp:c:/sliding/veh_vel2.bmp x=4in y=3in}
\caption{Simplified vehicle model for computing tire slip angles.}
\label{fig:veh_vel}
\end{figure}
Assuming small angles, we form an expression for the front tire slip angle
\begin{eqnarray}
\alpha_f &=& \delta_f-\tan^{-1}\left( \frac{l_1\omega_z+v_y}{v_x}\right)\nonumber\\
&\approx& \delta_f-\frac{l_1\omega_z+v_y}{v_x}\nonumber\\
&=& \delta_f-\frac{l_1\omega_z+\dot{y}}{\dot{x}}.
\label{eq:alphafexpr}
\end{eqnarray}
The rear tire slip angle is
\begin{eqnarray}
\alpha_r &=& \delta_r-\tan^{-1}\left(-\frac{l_2\omega_z-v_y}{v_x}\right)\nonumber\\
 &=& \delta_r+\tan^{-1}\left(\frac{l_2\omega_z-v_y}{v_x}\right)\nonumber\\
&\approx& \delta_r+\frac{l_2\omega_z-v_y}{v_x}\nonumber\\
&=& \delta_r+\frac{l_2\omega_z-\dot{y}}{\dot{x}}.
\label{eq:reatialphr}
\end{eqnarray}
Hence
\begin{equation}
F_{yf} = C_{\alpha f}\left( \delta_f-\frac{l_1\omega_z+\dot{y}}{\dot{x}}\right),
\label{eq:fyfwafacaf}
\end{equation}
and
\begin{equation}
F_{yr} = C_{\alpha r}\left( \delta_r+\frac{l_2\omega_z-\dot{y}}{\dot{x}}\right).
\label{eq:fyrcar}
\end{equation}
Substituting the above two equations into~(\ref{eq:suoflongdi}), and assuming small angles yields
\[
m\ddot{y}+m\dot{x}\omega_z=2C_{\alpha f}\delta_f-2C_{\alpha f}\frac{l_1\omega_z+\dot{y}}{v_x}+2C_{\alpha r}\delta_r+2C_{\alpha r}\frac{l_2\omega_z-\dot{y}}{v_x}.
\]
Rearranging terms and performing some manipulations, we get
\begin{equation}
\ddot{y}=-\left( \frac{2C_{\alpha f}+2C_{\alpha r}}{mv_x}\right)\dot{y}+\left( -v_x-\frac{2C_{\alpha f}l_1-2C_{\alpha r}l_2}{mv_x}\right)\omega_z +\frac{2 C_{\alpha f}}{m}\delta_f+\frac{2 C_{\alpha r}}{m}\delta_r.
\label{eq:lataxcafr}
\end{equation}
Substituting the expressions for $F_{yf}$ and $F_{yr}$ into~(\ref{eq:rotnelabom}) and assuming small angles, we obtain
\[
I_z\dot{\omega}_z=2l_1C_{\alpha f}\left( \delta_f -\frac{l_1\omega_z + \dot{y}}{v_x}\right) -2l_2C_{\alpha r}\left(\delta_r+\frac{l_2\omega_z-\dot{y}}{v_x}\right).
\]
Rearranging and performing manipulations yields
\begin{equation}
\dot{\omega}_z=-\left(\frac{2l_1C_{\alpha f} - 2l_2C_{\alpha r}}{I_zv_x} \right)\dot{y}-\left( \frac{2l_1^2C_{\alpha f}+2l_2^2C_{\alpha r}}{I_zv_x}\right)\omega_z+\frac{2l_1C_{\alpha f}}{I_z}\delta_f-\frac{2l_2C_{\alpha r}}{I_z}\delta_r.
\label{eq:rotnewaftman}
\end{equation}
Later, we will analyze a vehicle performance in the fixed coordinates. In particular, we will be interested in an emergency lane change maneuver\index{emergency lane change maneuver}. This requires that we transform the vehicle lateral variable into the fixed coordinates. Assuming small angles and referring to Figure~\ref{fig:veh_lat}, we obtain
\begin{eqnarray}
\dot{Y} &=& -\dot{x}\sin \psi -\dot{y}\cos \psi\nonumber\\
&\approx& -v_x\psi - \dot{y}.
\label{eq:lineexfobiy}
\end{eqnarray}
\begin{figure}
\centering
\includegraphics[scale=0.55]{VEH_LAT.eps}
%\vspace{3.5in}
%\hskip.7in\special{bmp:c:/sliding/veh_lat.bmp x=4.67in y=3.5in}
\caption{Computing the velocity of a vehicle with respect to the fixed coordinates.}
\label{fig:veh_lat}
\end{figure}
Note that
\[
\dot{X} = \dot{x}\cos \psi - \dot{y}\sin \psi.
\]
For the step lane change, the yaw rate of the vehicle is equal to the yaw rate error, that is,
\[
\omega_z=\dot{\psi}.
\]
We further assume that
\begin{equation}
\theta=\psi.
\label{eq:theaeqpsi}
\end{equation}
We can now represent the vehicle modeling equations in a state space format. Let
\[
x_1 =\dot{y},\quad x_2=\theta,\quad x_3=\omega_z,\quad x_4=Y.
\]
Then, using~(\ref{eq:lataxcafr}), (\ref{eq:rotnewaftman}), (\ref{eq:theaeqpsi}) and~(\ref{eq:lineexfobiy}), we obtain
\begin{eqnarray}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]&=&\left[\begin{array}{cccc}
-\frac{2C_{\alpha f}+2C_{\alpha r}}{mv_x}& 0 & -v_x-\frac{2C_{\alpha f}l_1-2C_{\alpha r}l_2}{mv_x}& 0\\
0 & 0 & 1 & 0\\
-\frac{2l_1C_{\alpha f} - 2l_2C_{\alpha r}}{I_zv_x}& 0 & -\frac{2l_1^2C_{\alpha f}+2l_2^2C_{\alpha r}}{I_zv_x}& 0\\
-1 & -v_x & 0 & 0
\end{array}\right]\left[\begin{array}{c}
{x_1}\\
{x_2}\\
{x_3}\\
{x_4}
\end{array}\right]\nonumber\\
&{}&\quad+\left[\begin{array}{cc}
\frac{2C_{\alpha f}}{m}&\frac{2C_{\alpha r}}{m} \\
0 & 0\\
\frac{2l_1C_{\alpha f}}{I_z}&-\frac{2l_2C_{\alpha r}}{I_z} \\
0 & 0
\end{array}\right]\left[\begin{array}{c}
\delta_f\\
\delta_r
\end{array}\right].
\label{veh_model}
\end{eqnarray}
A possible vehicle parameters are given in Table~\ref{tab:par}.
\begin{table}
\caption{Typical vehicle data.}
\label{tab:par}
\begin{center}
\begin{tabular}{cc}\hline
$C_{\alpha f}=C_{\alpha r}$ & $30,000.00\,\mb{N/rad}$\\
$l_1$ & $1.20\,\mb{m}$\\
$l_2$ & $1.22\,\mb{m}$\\
$m$ & $1280\,\mb{kg}$\\
$I_z$ & 2500\,kg$\cdot$m$^2$\\
$v_x$ & $18.3\,\mb{m/sec}$\\ \hline
\end{tabular}
\end{center}
\end{table}


\subsection{Permanent Magnet Stepper Motor}
\label{subsec:pmstepmo}
In this subsection, we develop a model of a permanent magnet~(PM) stepper motor\index{permanent magnet~(PM) stepper motor}\index{PM|see{permanent magnet}}. Stepper motors can be viewed as digital positioning devices. Specifically, stepper motors are electromechanical motion devices\index{electromechanical motion device|see{permanent magnet~(PM) stepper motor}} that convert input information in digital form into output that has mechanical form. The stepper motor can be driven by digital pulses, for example of the form: move clockwise or move counterclockwise. We now describe an idea of an operation of a PM stepper motor and then derive its modeling equations.

A simplified organization of a stepper motor is shown in Figure~\ref{fig:step_org}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{STEP_ORG.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/figures/step_org.bmp x=4.67in y=3.5in}
%\hskip0.75in\special{bmp:c:/figures/step_org.bmp x=4.67in y=3.5in}
\caption{A simplified schematic of an operation of a stepper motor.}
\label{fig:step_org}
\end{figure}
The armature\index{armature}, or rotor\index{rotor}\index{rotor|see{armature}}, is a magnet. At any instant of the operation, exactly one pair of windings of the stator is energized. A pair of windings represents a phase. Each pair of the stator windings is wound in such a way that their magnetic fields are collinear, for example, windings A and A'. The rotor then rotates until its magnetic field is aligned with the fields of the energized pair of the stator windings. Figure~\ref{fig:step_org} shows the situation where the pair of windings~A and A' is energized. If the rotor is at rest, as shown in Figure~\ref{fig:step_org}, and the pair~A, A' is de-energized, while the pair~B, B' is energized, the rotor rotates to the new position so that its magnetic field is aligned with the magnetic field of the pair~B, B'. In this example of the three-phase motor\index{three-phase motor}, a single step results in the rotor motion of $360^{\circ}/6=60^{\circ}$. We next give a more detailed description of the electrical organization of the stepper motor and construct its model.

The basic PM stepper motor consists of a slotted stator and a rotor. There can be two sets of teeth on the rotor. These two sets of teeth are out of alignment with each other by a tooth width\index{tooth width}, as shown in Figure~\ref{fig:rot_tee}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{ROT_TEE.eps}
%\vspace{3.5in}
%\hskip0.75in\special{bmp:c:/figures/rot_tee.bmp x=4.67in y=3.5in}
%\centerbmp{4.67in}{3.5in}{c:/figures/rot_tee.bmp x=4.67in y=3.5in}
\caption{A possible rotor structure of a PM stepper motor.}
\label{fig:rot_tee}
\end{figure}
Specifically, the rotor teeth of the left end-cap are displaced one half pitch from the teeth on the right end-cap. A magnet, or an electromagnet\index{electromagnet}, that is mounted on the rotor magnetizes the iron end-caps. The left end-cap is magnetized as a north pole, while the right end-cap is magnetized as a south pole. The slotted stator is equipped with two or more individual coils. Figure~\ref{fig:pm_step} shows a four-pole stator and a five-pole rotor structure.
\begin{figure}
\centering
\includegraphics[scale=0.75]{PM_STEPA.eps}
%\vskip3.5in
%\hskip.5in\special{bmp:c:/figures/pm_stepa.bmp x=4.67in y=3.5in}
%\vskip3.5in
%\hskip.5in\special{bmp:c:/figures/pm_stepb.bmp x=4.67in y=3.5in}
%\centerbmp{4.67in}{3.5in}{c:/figures/pm_step.bmp x=4.67in y=3.5in}
\caption{Construction features of a four-pole stator and a five-pole rotor PM stepper motor.}
\label{fig:pm_step}
\end{figure}

We now discuss the way in which the PM stepper motor can be used to perform precise positioning\index{precise positioning}. Assume that the starting point is as represented in Figure~\ref{fig:pm_start}.
\begin{figure}
\centering
\includegraphics[scale=0.75]{PM_START.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/figures/pm_start.bmp x=4.67in y=3.5in}
%\hskip0.75in\special{bmp:c:/figures/pm_start.bmp x=4.67in y=3.5in}
\caption{Coil A-A' energized for an S-N orientation, while coil B-B' de-energized.}
\label{fig:pm_start}
\end{figure}
Here, coil A-A' is energized so that the upper stator pole is given an S-pole polarity and the corresponding lower stator is given an N-pole polarity. The rotor remains in this position as long as stator coil voltage, and hence stator polarities, remain unchanged. This is because the net torque acting on the rotor is zero. If we now de-energize coil A-A' and energize coil B-B' so that the left pole becomes an S-pole and the right pole assumes an N-polarity, then this results in the new steady-state position shown in Figure~\ref{fig:pm_new}.
\begin{figure}
\centering
\includegraphics[scale=0.75]{PM_NEW.eps}
%\vspace{3.5in}
%\hskip0.75in\special{bmp:c:/figures/pm_new.bmp x=4.67in y=3.5in}
%%\centerbmp{4.67in}{3.5in}{c:/figures/pm_new.bmp x=4.67in y=3.5in}
\caption{Coil A-A' de-energized, while coil B-B' energized for an S-N orientation.}
\label{fig:pm_new}
\end{figure}
We can proceed by alternatively energizing the two stator coils. The PM stepper motor can thus be controlled by digital pulses applied to its two inputs: move clockwise or move counterclockwise.

Our next objective is to devise a mathematical model of the PM stepper motor\index{PM stepper motor!mathematical model}. In our discussion we use the following notation. We denote by $R_t$ the number of rotor teeth. The angular displacement between rotor teeth is called the tooth pitch and denoted $T_p$. The angular displacement of the rotor as we change the input voltage from one pair of windings to the other is called the step length and denoted $S_l$. The number of phases, that is, the number of pairs of stator windings, is denoted $N_p$. In our example, we have $N_p=2$, $R_t=5$, and $T_p=2\pi/R_t=72^{\circ}$. In this example, it takes four input voltage switchings, and hence steps, to advance the rotor one tooth pitch. Indeed, because
\begin{equation}
T_p=2N_pS_l=\frac{2\pi}{R_t},
\label{eq:tooth_pitch}
\end{equation}
hence
\begin{equation}
S_l=\frac{\pi}{R_tN_p}=18^{\circ}.
\label{eq:step_leng}
\end{equation}
We denote the rotor angular displacement by $\theta$.

Following Krause and Wasynczuk~\cite[pages~350--351]{krause89}, we trace the main path of flux linking\index{main path of flux linking} the A-A' windings for the rotor position shown in Figure~\ref{fig:pm_flux}.
\begin{figure}
\centering
\includegraphics[scale=0.75]{PM_FLUXA.eps}
%\vskip3.5in
%\hskip.5in\special{bmp:c:/figures/pm_fluxa.bmp x=4.67in y=3.5in}
%\vskip3.5in
%\hskip.5in\special{bmp:c:/figures/pm_fluxb.bmp x=4.67in y=3.5in}
%%\centerbmp{4.67in}{3.5in}{c:/figures/pm_flux.bmp x=6in y=4.5in}
\caption{The main path of flux linking the A-A' winding.}
\label{fig:pm_flux}
\end{figure}
Flux leaves the left end-cap through the rotor tooth at the top that is aligned with the stator tooth having the A part of the A-A' winding. The flux then travels up through the stator tooth in the stator iron. Next, the flux splits and travels around the circumference of the stator, as shown in Figure~\ref{fig:pm_flux}(a). The flux returns to the south pole of the rotor through the stator tooth positioned at the bottom of the A' part of the A-A' winding, as shown in Figure~\ref{fig:pm_flux}(b). We denote the peak flux linkage produced by the rotor permanent magnet, in the A-A' stator winding by $n\phi_m$, where $n$ is the number of turns of the winding. For $\theta=0^{\circ}$, corresponding to the rotor position in Figure~\ref{fig:pm_flux}, the flux linkage\index{flux linkage} in the A-A' winding of the stator is maximal and is approximately equal to $n\phi_m$. For $\theta=90^{\circ}/R_t$, corresponding to the rotor position in Figure~\ref{fig:pm_new}, the flux linkage in the A-A' winding of the stator is minimal, approximately equal to zero. We model the flux linkage, $\psi_a$, produced by the permanent magnet, in the A-A' stator winding as
\begin{equation}
\psi_a=n\phi_m\cos(R_t\theta).
\label{eq:flux_a}
\end{equation}
The (main) flux, linking the B-B' winding for the rotor position shown in Figure~\ref{fig:pm_new}, would enter the stator tooth on which the B part of the B-B' winding is located from the rotor tooth on the left. The flux then would travel around the circumference of the stator and return to the rotor through the stator pole on which the B' part of the B-B' winding is wound. For $\theta=90^{\circ}/R_t$, corresponding to the rotor position in Figure~\ref{fig:pm_new}, the flux linkage in the B-B' winding of the stator is maximal and is approximately equal to $n\phi_m$. For $\theta=0^{\circ}$, corresponding to the rotor position in Figure~\ref{fig:pm_flux}, the flux linkage in the B-B' winding is minimal, approximately equal to zero.  We model the flux linkage in the B-B' phase of the stator by
\begin{equation}
\psi_b=n\phi_m \sin(R_t\theta).
\label{eq:flux_b}
\end{equation}
\begin{table}
\caption{Symbols used in modeling equations of the PM stepper motor.}
\begin{center}
\begin{tabular}{l}\hline
$T_p$ = tooth pitch\\
$S_l$ = step length\\
$R_t$ = number of teeth on rotor\\
$N_p$ = number of phases, that is, stator winding pairs\\
$\theta$ = angular displacement of the rotor\\
$\omega$ = angular velocity of the rotor\\
$v_a(t)$ = voltage applied to phase $a$, that is, A-A' stator winding\\
$v_b(t)$ = voltage applied to phase $b$, that is, B-B' stator winding\\
$i_a(t)$ = current in phase $a$\\
$i_b(t)$ = current in phase $b$\\
$R$ = stator winding resistance per phase\\
$L$ = stator winding inductance per phase\\
$B$ = viscous friction coefficient\\
$I$ = rotor inertia\\
$K_b$ = back emf constant\\\hline
\end{tabular}
\end{center}
\label{table:pm_step}
\end{table}
In the following development, we use symbols that are described in Table~\ref{table:pm_step}. Applying Kirchhoff's voltage law\index{Kirchhoff's voltage law} to the circuit of phase $a$ yields
\[
v_a(t)= R i_a(t)+L\frac{d i_a(t)}{dt}+e_a(t),
\]
where $e_a$ is the back emf induced in phase $a$, that is,
\begin{eqnarray*}
e_a &=& \frac{d\psi_a}{dt}\\
&=& \frac{d\psi_a}{d\theta}\frac{d\theta}{dt}\\
&=& -nR_t\phi_m \sin(R_t \theta)\frac{d\theta}{dt}\\
&=&-K_b \,\omega \sin(R_t \theta),
\end{eqnarray*}
where $K_b=nR_t\phi_m$. Hence, the voltage equation for phase $a$ can be represented as
\[
v_a=Ri_a+L\frac{di_a}{dt}-K_b \,\omega \sin(R_t \theta).
\]
Applying now Kirchhoff's voltage law to the circuit of phase $b$ yields
\[
v_b(t)= R i_b(t)+L\frac{d i_b(t)}{dt}+e_b(t),
\]
where $e_b$ is the back emf induced in phase $b$, that is,
\begin{eqnarray*}
e_b &=& \frac{d\psi_b}{dt}\\
&=& \frac{d\psi_b}{d\theta}\frac{d\theta}{dt}\\
&=& nR_t\phi_m \cos(R_t \theta)\frac{d\theta}{dt}\\
&=&K_b \,\omega \cos(R_t \theta),
\end{eqnarray*}
where as before $K_b=nR_t\phi_m$. Hence, the voltage equation for phase $b$ can be represented as
\[
v_b=Ri_b+L\frac{di_b}{dt}+K_b \,\omega \cos(R_t \theta).
\]
Let $T_a$ and $T_b$ denote torques produced by the currents $i_a$ and $i_b$, respectively. Then, applying Newton's second law\index{Newton's second law}\index{Newton's second law!applied to rotor} to the rotor yields
\[
I\frac{d^2 \theta}{dt^2}+B\frac{d\theta}{dt}=T_a+T_b.
\]
The PM stepper motor is an electromechanical device in which electric power is transformed into mechanical power\index{electric power!transformed into mechanical power}. The mechanical power, $P_a$, due to torque $T_a$ is
\begin{equation}
P_a=T_a\frac{d\theta}{dt}.
\label{eq:potadth}
\end{equation}
Assuming that there are no losses, this mechanical power is equal to the power developed in phase $a$,
\begin{equation}
P_a=i_a\frac{d\psi_a}{dt}=i_a\frac{d\psi_a}{d\theta}\frac{d\theta}{dt}.
\label{eq:therminpo}
\end{equation}
Equating the right hand sides of~(\ref{eq:potadth}) and~(\ref{eq:therminpo}) and performing simple manipulations, we obtain
\[
T_a = i_a\frac{d\psi_a}{d \theta}= -i_aK_b \sin(R_t\theta).
\]
Applying similar arguments to phase $b$, we get
\[
T_b = i_b\frac{d \psi_b}{d \theta}= i_b K_b \cos(R_t\theta).
\]
Taking the above equations into account, we represent the equation of motion of the rotor as
\[
I\frac{d^2 \theta}{dt^2}+B\frac{d\theta}{dt}=T_a+T_b=-i_aK_b \sin(R_t\theta)+i_b K_b \cos(R_t\theta).
\]
Combining the above equations, we obtain a model of a two-phase, that is, a two input PM stepper motor\index{two input PM stepper motor}. We have
\[
\left.\begin{array}{lll}
\frac{d\theta}{dt}&=&\omega\\
\frac{d\omega}{dt}&=& -\frac{K_b}{I}i_a \sin(R_t\theta) + \frac{K_b}{I}i_b \cos(R_t\theta)-\frac{B}{I}\omega\\
\frac{di_a}{dt} &=& -\frac{R}{L}i_a+\frac{K_b}{L}\omega \sin(R_t\theta)+\frac{1}{L}v_a\\
\frac{di_b}{dt} &=& -\frac{R}{L}i_b-\frac{K_b}{L}\omega \cos(R_t\theta)+\frac{1}{L}v_b.
\end{array}\right\}
\]
We now define the state variables: $x_1=\theta$, $x_2=\omega$, $x_3=i_a$, and $x_4=i_b$. We let $K_1=B/I$, $K_2=K_b/I$, $K_3=R_t$, $K_4=K_b/L$, $K_5=R/L$, $u_1=\frac{1}{L}v_a$, and $u_2=\frac{1}{L}v_b$. Using the above notation, we represent the PM stepper motor model in state space format:
\begin{eqnarray}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]&=&\left[\begin{array}{c}
x_2\\
-K_1x_2-K_2x_3\sin(K_3x_1)+K_2x_4\cos(K_3x_1)\\
K_4x_2\sin(K_3x_1)-K_5x_3\\
-K_4x_2\cos(K_3x_1)-K_5x_4
\end{array}\right]+\left[\begin{array}{cc}
0 & 0\\
0 & 0\\
1 & 0\\
0 & 1
\end{array}\right]\left[\begin{array}{c}
u_1\\
u_2
\end{array}\right].
\label{eq:nomoopmstep}
\end{eqnarray}
\subsection{Stick Balancer}
\label{subsection:truth}
In this subsection, we develop a nonlinear model of the stick balancer\index{stick balancer}\index{stick balancer|see{cart with an inverted pendulum}}, also known as the cart with an inverted pendulum\index{cart with an inverted pendulum}. In our derivation of the model, we use Newton's laws, in a similar fashion as Kwakernaak and Sivan~\cite{kwakernaak72}, and Ogata~\cite{ogata97} did, where linearized models of the stick balancer were developed. An alternative derivation, using the method of D'Alembert, can be found in Cannon~\cite[Section 22.4]{cannon67}.

A free body diagram\index{free body diagram} of a cart on which a stick, that is, an inverted pendulum, is mounted is shown in Figure~\ref{fig:inp_free}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{INP_FREE.eps}
%\vspace{4in}
%\hskip.75in\special{bmp:c:/papers/inp_free.bmp x=5.33in y=4in}
%%\centerbmp{4.67in}{3.5in}{c:/papers/inp_free.bmp x=4.67in y=3.5in}
\caption{Free body diagram of a cart with an inverted pendulum.}
\label{fig:inp_free}
\end{figure}
Let $H=H(t)$ and $V=V(t)$ be horizontal and vertical reaction forces, respectively. The $x$ and $y$ are the coordinates of the fixed, nonrotating coordinate frame\index{nonrotating coordinate frame} $x-y$. The angular displacement of the stick from the vertical position is $\theta=\theta(t)$. The mass of the cart is $M$, while the mass of the stick is $m$. The length of the stick is $2 l$, and its center of gravity is at its geometric center. The control force applied to the cart is labeled as $u$. We assume that the wheels of the cart do not slip. We model the frictional force of the cart wheels on the track by
\begin{equation}
f_c=\mu_c\mb{sign}(\dot{x}),
\label{eq:kinetic_fr}
\end{equation}
where $\mu_c$ is the cart friction coefficient\index{cart friction coefficient}. Friedland~\cite[page 201]{friedland96} refers to this model of friction as the classical Coulomb friction model. Let $(x_G,y_G)$ be the coordinates of the center of gravity of the stick. Then,
\begin{equation}
\left\{\begin{array}{lcl}
x_G & = & x+l\sin(\theta)\\
y_G &=& l\cos(\theta).
\end{array}\right.
\label{eq:cog_stick}
\end{equation}
We are now ready to write down equations modeling the system. The equation that describes the rotational motion of the stick about its center of gravity is obtained by applying the rotational version of Newton's second law. Summing the moments about the center of gravity of the stick, we obtain
\begin{equation}
I\frac{d^2\theta}{dt^2}=Vl\sin(\theta)-Hl\cos(\theta),
\label{eq:newton_rot}
\end{equation}
where
\begin{equation}
I=\int_{-l}^lr^2dm=ml^2/3,
\label{eq:mom_inertia}
\end{equation}
is the moment of inertia of the stick with respect to its center of gravity.

We next write the equation that describes the horizontal motion of the center of gravity of the stick. Applying Newton's second law along the $x$ axis yields
\begin{equation}
m\frac{d^2}{dt^2}\left( x+l\sin(\theta)\right)=H.
\label{eq:newt_hor}
\end{equation}
Performing the required differentiation, we get
\begin{equation}
m\left(\ddot{x}+l\left(-\dot{\theta}^2\sin(\theta)+\ddot{\theta}\cos(\theta)\right)\right)=H.
\label{eq:newt_hort}
\end{equation}
The equation that describes the vertical motion of the center of gravity of the stick is obtained by applying Newton's second law\index{Newton's second law} along the $y$ axis. We have
\begin{equation}
m\frac{d^2}{dt^2}\left( l\cos(\theta)\right)=V-mg.
\label{eq:newt_vert}
\end{equation}
Performing the needed differentiation yields
\begin{equation}
ml\left(-\dot{\theta}^2\cos(\theta)-\ddot{\theta}\sin(\theta)\right)=V-mg.
\label{eq:newt_vertt}
\end{equation}
Finally, we apply Newton's second law to the cart to get
\begin{equation}
M\frac{d^2x}{dt^2}=u-H-f_c.
\label{eq:newt_cart}
\end{equation}
Substituting~(\ref{eq:newt_cart}) into~(\ref{eq:newt_hort}) gives
\begin{equation}
m\ddot{x}+ml\ddot{\theta}\cos(\theta)-ml\dot{\theta}^2\sin(\theta)+f_c=u-M\ddot{x}.
\label{eq:four_two}
\end{equation}
Substituting now equations~(\ref{eq:newt_vertt}) and~(\ref{eq:newt_cart}) into~(\ref{eq:newton_rot}) yields
\begin{equation}
I\ddot{\theta}=\left( mg-ml\dot{\theta}^2\cos(\theta)-ml\ddot{\theta}\sin(\theta)\right)l\sin(\theta)+f_cl\cos(\theta)-\left(u-M\ddot{x}\right)l\cos(\theta).
\label{eq:three_four}
\end{equation}
We next substitute $u-M\ddot{x}$ from~(\ref{eq:four_two}) into~(\ref{eq:three_four}) and perform manipulations to get
\begin{equation}
I\ddot{\theta}=mgl\sin(\theta)-ml^2\ddot{\theta}-m\ddot{x}l\cos(\theta).
\label{eq:almost}
\end{equation}
Let
\[
a=\frac{1}{m+M}.
\]
Then, we can represent~(\ref{eq:four_two}) as
\begin{equation}
\ddot{x}=-mal\ddot{\theta}\cos(\theta)+mal\dot{\theta}^2\sin(\theta)-af_c+au.
\label{eq:xtwodot}
\end{equation}
We substitute~(\ref{eq:xtwodot}) into~(\ref{eq:almost}) to obtain
\begin{equation}
\ddot{\theta}=\frac{mgl\sin(\theta)-m^2l^2a\dot{\theta}^2\sin(2\theta)/2+ mal\cos(\theta)f_c - mal\cos(\theta)u}{I-m^2l^2a\cos^2(\theta)+ml^2}.
\label{eq:theta_fin}
\end{equation}
Let $x_1=\theta$ and $x_2=\dot{\theta}$. Using the expression for $I$, given by~(\ref{eq:mom_inertia}), we represent~(\ref{eq:theta_fin}) in state space format as follows
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1 &=& x_2\\
\dot{x}_2 &=& \frac{g\sin(x_1)-mlax_2^2\sin(2x_1)/2+ a\cos(x_1)f_c}{4l/3-mla\cos^2(x_1)}-\frac{mal\cos(x_1)u}{4ml^2/3-m^2l^2a\cos^2(x_1)}.
\end{array}\right\}
\label{eq:mod_space}
\end{equation}
We next substitute $\ddot{\theta}$, obtained from~(\ref{eq:almost}), into~(\ref{eq:xtwodot}) to get
\begin{equation}
\ddot{x}=\frac{-mag\sin(2x_1)/2+ax_2^2\sin(x_1)4ml/3+(u-f_c)4a/3}{4/3-ma\cos^2(x_1)}.
\label{eq:xddot}
\end{equation}
Let $x_3=x$ and $x_4=\dot{x}$. Then, combining~(\ref{eq:mod_space}) and~(\ref{eq:xddot}) we obtain a state space model of the inverted pendulum on a cart of the form,
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]=\left[\begin{array}{c}
x_2\\
\frac{g\sin(x_1)-mlax_2^2\sin(2x_1)/2}{4l/3-mla\cos^2(x_1)}\\
{x}_4\\
\frac{-mag\sin(2x_1)/2+alx_2^2\sin(x_1)4m/3}{4/3-ma\cos^2(x_1)}\end{array}\right]+\left[\begin{array}{c}
0\\
\frac{-a\cos(x_1)}{4l/3-mla\cos^2(x_1)}\\
0\\
\frac{4a/3}{4/3-ma\cos^2(x_1)}\end{array}\right](u-f_c).
\label{eq:state_stick}
\end{equation}
While designing a controller for the stick balancer system, we will use an equivalent representation of~(\ref{eq:state_stick}),
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]=\left[\begin{array}{c}
x_2\\
\frac{g\sin(x_1)}{4l/3-mla\cos^2(x_1)}\\
{x}_4\\
\frac{-mag\sin(2x_1)/2}{4/3-ma\cos^2(x_1)}\end{array}\right]+\left[\begin{array}{c}
0\\
\frac{-a\cos(x_1)}{4l/3-mla\cos^2(x_1)}\\
0\\
\frac{4a/3}{4/3-ma\cos^2(x_1)}\end{array}\right](u-f_c+mlx_2^2\sin(x_1)).
\label{eq:stick_eqiv}
\end{equation}
\label{ex:fullmo}


\subsection{One-Link Robot Manipulator}
\label{subsection:one_link_robot}


The example comes from~\cite{zak86}. It is a simple model of a one-link robot manipulator\index{robot manipulator!one-link} shown in Figure~\ref{fig:onelinkex}. The motion of the robot arm is controlled by a DC motor\index{DC motor} via a gear\index{gear}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{onelink.eps}
%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/onelink.bmp}
%\hskip1.4in\special{bmp:c:/mybook/onelink.bmp x=4.67in y=3.5in}
\caption{One-link manipulator controlled by a DC motor via a gear.}
\label{fig:onelinkex}
\end{figure}
The DC motor is armature-controlled and its schematic is shown in Figure~\ref{fig:dcmotorex}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{dcmotor.eps}
%\vspace{3in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/dcmotor.bmp}
%\hskip1.5in\special{bmp:c:/mybook/dcmotor.bmp x=4in y=3in}
\caption{Schematic of an armature-controlled DC motor of Example~\ref{example:onelink}.}
\label{fig:dcmotorex}
\end{figure}
We assume that the motor moment of inertia is negligible compared with that of the robot arm. We model the arm as a point mass $m$ attached to the end of a massless rod of length $l$. Hence the arm inertia is $I_a=ml^2$. We assume that the gear train has no backlash\index{backlash}, and all connecting shafts are rigid. As can be seen from Figure~\ref{fig:onelinkex}, counterclockwise rotation of the arm is defined as positive and clockwise rotation\index{clockwise rotation} of the arm as negative, while counter-clockwise rotation\index{counter-clockwise rotation} of the motor shaft is defined to be negative, and clockwise rotation of the shaft as positive. The torque delivered by the motor is
\[
T_m=K_mi_a,
\]
where $K_m$ is the motor-torque constant\index{motor-torque constant}, and $i_a$ is the armature current\index{armature current}. Let $N$ denote the gear ratio\index{gear ratio}. Then we have
\[
\frac{\theta_p}{\theta_m}=\frac{\mb{radius of motor gear}}{\mb{radius of arm gear}}=\frac{\mb{number of teeth of motor gear}}{\mb{number of teeth of arm gear}} = \frac{1}{N}.
\]
This is because the gears are in contact and thus
\[
\theta_p \times\mb{radius of arm gear}=\theta_m \times\mb{radius of motor gear},
\]
and the radii of the gears are proportional to their number of teeth. The work done by the gears must be equal. Let $T_p$ denote the torque applied to the robot arm. Then,
\[
T_p\theta_p = T_m\theta_m.
\]
Thus, the torque applied to the pendulum is
\[
T_p=NT_m=NK_mi_a.
\]
We use Newton's second law to write the equation modeling the arm dynamics,
\begin{equation}
I_a\frac{d^2\theta_p}{dt^2}=mgl\sin \theta_p +T_p.
\label{eq:prenewton}
\end{equation}
Substituting into~(\ref{eq:prenewton}) the expressions for $I_a$ and $T_p$ and then rearranging yields
\begin{equation}
ml^2\frac{d^2\theta_p}{dt^2}=mgl\sin \theta_p + NK_m i_a,
\label{eq:newton_t}
\end{equation}
where $g=9.8\;\mb{m}/\mb{s}^2$ is the gravitational constant. We next apply Kirchhoff's voltage law to the armature circuit to obtain
\[
L_a\frac{d i_a}{dt}+Ri_a+e_b=L_a\frac{d i_a}{dt}+Ri_a+K_b N\frac{d\theta_p}{dt}=u,
\]
where $K_b$ is the back emf constant.  We can now construct a third-order state-space model of the one-link robot. For this we choose the following state variables:
\[
x_1=\theta_p, \quad x_2=\frac{d\theta_p}{dt}=\omega_p,\quad x_3=i_a.
\]
Then, the model in state-space format is:
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array} \right] = \left[\begin{array}{c}
x_2\\
\frac{g}{l}\sin x_1+\frac{NK_m}{ml^2}x_3\\
-\frac{K_bN}{L_a}x_2-\frac{R_a}{L_a}x_3
\end{array}\right]+\left[\begin{array}{c}
0\\
0\\
\frac{1}{L_a}
\end{array}\right]u,
\label{eq:nomooolima}
\end{equation}
Reasonable parameters of the robot 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$, $L_a=100\,\mb{mH}$. With the above parameters the robot model takes the form:
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3
\end{array} \right] = \left[\begin{array}{c}
x_2\\
9.8\,\sin x_1+x_3\\
-10\,x_2-10\,x_3
\end{array}\right]+\left[\begin{array}{c}
0\\
0\\
10
\end{array}\right]u.
\label{eq:nomoolinum}
\end{equation}
We assume that the output, $y$, is
\begin{equation}
y=x_1.
\label{eq:ouoonlima}
\end{equation}

\subsection{Double Pendulum}
\label{subsection:two_link_robot}
We construct a model of the double pendulum\index{double pendulum} shown in Figure~\ref{fig:dbl_pend}.
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{dbl_pend.eps}
\caption{Double pendulum.}
\label{fig:dbl_pend}
\end{figure}
This double pendulum can be viewed as a model of a double-link robot manipulator\index{robot manipulator!double-link}. We assume that the input torques, $T_1$ and $T_2$ are applied at each joint, $\omega_i=\dot{\theta}_i$, $i=1,2$, are the angular velocities of the pendulum arms, $\Delta \theta=\theta_2-\theta_1$, and $m_i$ and $l_i$ are the masses and arm lengths, respectively. The parameter values are given in Table~\ref{tab:dbl_pendpar}.
\begin{table}
\caption{Parameter numerical values of double pendulum system shown in Figure~\ref{fig:dbl_pend}.}
\begin{center}
\begin{tabular}{ll}\hline
$m_1=m_2$ & 1 kg\\
$l_1=l_2$ & 0.5 m\\
%$I_1=I_2$& 1/48 kg$\cdot$m$^2$\\
g& 9.8 m/sec$^2$ \\ \hline
\end{tabular}
\end{center}
\label{tab:dbl_pendpar}
\end{table}
We write the Lagrange equations of motion for the system and represent them in state-space format, where the state variables are defined as, $x_1=\theta_1$, $x_2=\dot{\theta}_1$, $x_3=\theta_2$, and $x_4= \dot{\theta}_2$.

\noindent{\bf Lagrange Equations of Motion of the Double Pendulum System}

Recall that the Lagrangian\index{Lagrangian} $L$ is defined as $L=K-U$, where $K$ is the kinetic energy of the double pendulum\index{double pendulum!kinetic energy} system shown in Figure~\ref{fig:dbl_pend} and $U$ is the potential energy of the system.
%\begin{figure}
%\centerbmp{4in}{3in}{c:/mysolns/dbl_pend.bmp}
%\caption{Double pendulum of Exercise~\ref{exer:dbl_pend}.}
%\label{fig:dbl_pend}
%\end{figure}
The kinetic energy is
\begin{eqnarray*}
K&=&K_{m_1}+K_{m_2}\\
&=& \frac{1}{2}m_1\|\b v_{m_1}\|^2+\frac{1}{2}m_2\|\b v_{m_2}\|^2,
\end{eqnarray*}
where $\b v_{m_i}$, $i=1,2$, are the velocities of the point masses. The velocity vectors can be computed from the position of the point masses. The position of the first mass is described by
\begin{equation}
x_1=l_1 \sin (\theta_1),\quad y_1=l_1 \cos (\theta_1).
\label{eq:massonepos}
\end{equation}
The position of $m_2$ is determined by
\begin{equation}
x_2=l_1 \sin (\theta_1) + l_2 \sin (\theta_2),\quad y_2=l_1 \cos (\theta_1) + l_2 \cos (\theta_2).
\label{eq:masstwopos}
\end{equation}
We next find time derivatives of the expressions in~(\ref{eq:massonepos}) and~(\ref{eq:masstwopos}). The time derivatives of $x_1$ and $y_1$ are,
\[
\dot{x}_1=l_1\dot{\theta}_1 \cos(\theta_1),\quad \dot{y}_1=-l_1\dot{\theta}_1 \sin(\theta_1).
\]
The time derivatives of $x_2$ and $y_2$ are,
\[
\dot{x}_2=l_1\dot{\theta}_1 \cos(\theta_1)+l_2\dot{\theta}_2 \cos(\theta_2),\quad \dot{y}_2=-l_1\dot{\theta}_1 \sin(\theta_1)-l_2\dot{\theta}_2 \sin(\theta_2).
\]
Using the above expressions, we compute
\[
\|\b v_{m_1}\|^2 = \dot{x_1}^2 + \dot{y_1}^2=\left( l_1\dot{\theta}_1\right)^2
\]
and
\begin{eqnarray*}
\|\b v_{m_2}\|^2 &=& \dot{x_2}^2 + \dot{y_2}^2\\
{}&=&\left(l_1 \dot{\theta}_1 \right)^2\cos^2(\theta_1)+2 l_1 l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1)\cos(\theta_1)+ \left(l_2\dot{\theta}_2 \right)^2\cos^2(\theta_2)\\
{}&{}&\quad +\left(l_1 \dot{\theta}_1 \right)^2\sin^2(\theta_1)+2 l_1 l_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1)\sin(\theta_1)+ \left(l_2\dot{\theta}_2 \right)^2\sin^2(\theta_2)\\
&=& \left(l_1 \dot{\theta}_1 \right)^2+\left(l_2 \dot{\theta}_2 \right)^2 + 2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos (\theta_2 - \theta_1),
\end{eqnarray*}
where we used the well known trigonometric identity,
\[
\cos(\theta_2 - \theta_1) = \cos(\theta_2)\cos(\theta_1) + \sin(\theta_2)\sin(\theta_1).
\]
In view of the above, the kinetic energy of the system is
\begin{eqnarray*}
K&=&K_{m_1}+K_{m_2}\\
&=& \frac{1}{2}m_1\|\b v_{m_1}\|^2+\frac{1}{2}m_2\|\b v_{m_2}\|^2\\
&=& \frac{1}{2}m_1 \left(l_1 \dot{\theta}_1 \right)^2 + \frac{1}{2}m_2 \left(\left(l_1 \dot{\theta}_1 \right)^2+\left(l_2 \dot{\theta}_2 \right)^2 + 2l_1 l_2\dot{\theta}_1\dot{\theta}_2\cos (\theta_2 - \theta_1)  \right).
\end{eqnarray*}
To compute the potential energy, we observe the positive direction of the $y$-coordinate and the direction of the gravity force are the same. Therefore, the potential energy is
\begin{eqnarray*}
U &=& -m_1 g y_1 - m_2g y_2\\
&=& -m_1gl_1\cos(\theta_1) - m_2 g \left(l_1\cos(\theta_1) + l_2\cos(\theta) \right)\\
&=&-(m_1+m_2)g l_1 \cos(\theta_1)-m_2gl_2 \cos(\theta_2).
\end{eqnarray*}
Hence, the Lagrangian function of the double pendulum system\index{Lagrangian function!double pendulum} is
\begin{eqnarray*}
L &=& K-U\\
&=& \frac{1}{2}m_1 \left(l_1 \dot{\theta}_1 \right)^2 + \frac{1}{2}m_2 \left(\left(l_1 \dot{\theta}_1 \right)^2+\left(l_2 \dot{\theta}_2 \right)^2 + 2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos (\theta_2 - \theta_1)  \right)\\
{}&{}&\quad + (m_1+m_2)g l_1 \cos(\theta_1) + m_2gl_2 \cos(\theta_2).
\end{eqnarray*}
We next derive the Lagrange equations of motion of the system. For this, we need to determine the generalized forces that constitute the right hand sides of the Lagrange equations of motion. Referring to Figure~\ref{fig:dbl_pwk}, we write,
\begin{eqnarray}
\delta W&=&T_1\delta \theta_1 + T_2\delta (\theta_2 - \theta_1)\nonumber\\
&=& (T_1-T_2)\delta\theta_1 + T_2 \delta\theta_2,
\label{virt_work}
\end{eqnarray}
where $T_1$ and $T_2$ are the torques applied at each joint.
\begin{figure}
\centering
\includegraphics[scale=0.5]{dbl_pwk.eps}
%\centerbmp{4in}{3in}{c:/mysolns/dbl_pwk.bmp}
\caption{Angles in the double pendulum system of Figure~\ref{fig:dbl_pend}.}
\label{fig:dbl_pwk}
\end{figure}
Hence the generalized forces are
\begin{equation}
F_{\theta_1}=T_1-T_2\quad\mb{and}\quad F_{\theta_2}=T_2.
\label{eq:genforces}
\end{equation}
To proceed, we need to find partial derivatives of the Lagrangian function. We first evaluate $\partial L/\partial \dot{\theta}_1$ to get
\begin{equation}
\frac{\partial L}{\partial \dot{\theta}_1}=m_1l_1^2\dot{\theta}_1 + m_2l_1^2 \dot{\theta}_1 + m_2l_1l_2\dot{\theta}_2 \cos(\theta_2 - \theta_1).
\label{eq:pwrtdottheta}
\end{equation}
Next, we find
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}_1}=m_1l_1^2\ddot{\theta_1}+m_2l_1^2\ddot{\theta}_1+m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_2 - \theta_1)-m_2 l_1 l_2 \dot{\theta}_2 (\dot{\theta}_2 - \dot{\theta}_1)\sin(\theta_2 - \theta_1).
\label{eq:tderofl}
\end{equation}
The partial derivative $\partial L/\partial \theta_1$ is
\begin{equation}
\frac{\partial L}{\partial \theta_1}= m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_2 - \theta_1) - (m_1+m_2)gl_1 \sin(\theta_1).
\label{eq:pawrttofl}
\end{equation}
Combining~(\ref{eq:genforces}), (\ref{eq:tderofl}), and~(\ref{eq:pawrttofl}), we obtain the first Lagrange equation of motion,
\begin{equation}
(m_1+m_2)l_1^2 \ddot{\theta}_1 + m_2l_1l_2\ddot{\theta}\cos(\theta_2 - \theta_1) - m_2l_1l_2\dot{\theta}_2^2 \sin(\theta_2 - \theta_1) + (m_1 + m_2) g l_1 \sin(\theta_1) = T_1-T_2.
\label{eq:firslageqmo}
\end{equation}
We next compute partial derivatives of $L$ that we use to write the second Lagrange equation of motion. We first find $\partial L/\partial \dot{\theta}_2$,
\begin{equation}
\frac{\partial L}{\partial \dot{\theta}_2}= m_2l_2^2 \dot{\theta}_2 + m_2 l_2 l_2 \dot{\theta}_1 \cos(\theta_2 - \theta_1).
\label{eq:pawrttht}
\end{equation}
The time derivative of the above expression is
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}_2} = m_2l_2^2\ddot{\theta}_2 + m_2l_1l_2\ddot{\theta}_1 \cos(\theta_2 - \theta_1) - m_2l_1l_2(\dot{\theta}_2 - \dot{\theta}_1)\dot{\theta}_1 \sin(\theta_2 - \theta_1).
\label{eq:tderofttwo}
\end{equation}
The partial derivative $\partial L/\partial \theta_2$ is
\begin{equation}
\frac{\partial L}{\partial \theta_2}= -m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_2 - \theta_1) - m_2gl_2 \sin(\theta_2).
\label{eq:pawrttof2}
\end{equation}
Combining~(\ref{eq:genforces}), (\ref{eq:tderofttwo}), and~(\ref{eq:pawrttof2}), we obtain the second Lagrange equation of motion,
\begin{equation}
m_2l_2^2\ddot{\theta}_2 + m_2l_1l_2\ddot{\theta}_1\cos(\theta_2 - \theta_1) + m_2l_1l_2\dot{\theta}_1^2\sin(\theta_2 - \theta_1) + m_2gl_2\sin(\theta_2)=T_2.
\label{eq:seclageqmo}
\end{equation}

\noindent{\bf State-Space Model of the Double Pendulum System}

Taking into account the definition of the state variables, we obtain
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]=\left[\begin{array}{c}
x_2\\
\ddot{\theta}_1\\
x_4\\
\ddot{\theta}_2
\end{array}\right].
\label{eq:stvecdpsys}
\end{equation}
Thus, to obtain a state-space model of the double pendulum\index{double pendulum!state-space model} system, we need to find the acceleration of each link. For this, we represent Lagrange equations~(\ref{eq:firslageqmo}) and~(\ref{eq:seclageqmo}) in matrix form,
\begin{eqnarray}
&{}&\left[\begin{array}{cc}
(m_1+m_2)l_1^2 & m_2l_1l_2\cos(\theta_2 - \theta_1)\\
m_2l_1l_2\cos(\theta_2 - \theta_1) & m_2l_2^2
\end{array}\right]\left[\begin{array}{c}
\ddot{\theta}_1\\
\ddot{\theta}_2
\end{array}\right]\nonumber\\
&{}&\qquad\qquad = \left[\begin{array}{c}
T_1 - T_2 + m_2l_1l_2\dot{\theta}_2^2 \sin(\theta_2 - \theta_1) - (m_1 + m_2) g l_1 \sin(\theta_1)\\
T_2-m_2l_1l_2\dot{\theta}_1^2\sin(\theta_2 - \theta_1) -m_2gl_2\sin(\theta_2)
\end{array}\right].
\label{eq:matftwlags}
\end{eqnarray}
Note that the square matrix, call it $\b M$, on the left-hand side of the above equation is invertible because its determinant is
\[
\det \b M = (m_1+m_2)m_2l_1^2l_2^2 - m_2^2l_1^2l_2^2 \cos^2(\theta_2 - \theta_1) > 0.
\]
Solving~(\ref{eq:matftwlags}) for $\ddot{\theta}_1$ and $\ddot{\theta}_2$, we obtain
\begin{equation}
\ddot{\theta}_1 = \frac{\begin{array}{l}
m_2l_2^2\left(T_1-T_2+ m_2l_1l_2\dot{\theta}_2^2 \sin(\theta_2 - \theta_1) - (m_1 + m_2) g l_1 \sin(\theta_1)   \right)\\
-m_2l_1l_2 \cos(\theta_2 - \theta_1)\left( T_2-m_2l_1l_2\dot{\theta}_1^2\sin(\theta_2 - \theta_1) -m_2gl_2\sin(\theta_2)
\right)
\end{array}}{(m_1+m_2)m_2l_1^2l_2^2 - m_2^2l_1^2l_2^2 \cos^2(\theta_2 - \theta_1)}
\label{eq:theaonedd}
\end{equation}
and
\begin{equation}
\ddot{\theta}_2 = \frac{\begin{array}{l}
-m_2l_1l_2\cos(\theta_2 - \theta_1)\left(T_1-T_2+ m_2l_1l_2\dot{\theta}_2^2 \sin(\theta_2 - \theta_1) - (m_1 + m_2) g l_1 \sin(\theta_1)   \right)\\
+(m_1+m_2)l_1^2\left( T_2-m_2l_1l_2\dot{\theta}_1^2\sin(\theta_2 - \theta_1) -m_2gl_2\sin(\theta_2)
\right)
\end{array}}{(m_1+m_2)m_2l_1^2l_2^2 - m_2^2l_1^2l_2^2 \cos^2(\theta_2 - \theta_1)}.
\label{eq:theatwodd}
\end{equation}
The input vector is defined as
\[
\b u=\left[\begin{array}{cc}
u_1 & u_2
\end{array}\right]^{\top}=\left[\begin{array}{cc}
T_1 & T_2
\end{array}\right]^{\top}.
\]
Taking into account the above definition of the input, the definition of the state variables given by~(\ref{eq:stvecdpsys}) and relations~(\ref{eq:theaonedd}) and~(\ref{eq:theatwodd}), we obtain the following state-space model of the double pendulum system
\begin{eqnarray*}
&{}&\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]=\b f(\b x,\b u)\\
 &=& \left[\begin{array}{c}
x_2\\
\frac{\begin{array}{l}
m_2l_2^2\left(u_1-u_2+ m_2l_1l_2x_4^2 \sin(x_3 - x_1) - (m_1 + m_2) g l_1 \sin(x_1)   \right)\\
-m_2l_1l_2 \cos(x_3 - x_1)\left( u_2-m_2l_1l_2x_2^2\sin(x_3 - x_1) -m_2gl_2\sin(x_3)
\right)
\end{array}}{(m_1+m_2)m_2l_1^2l_2^2 - m_2^2l_1^2l_2^2 \cos^2(x_3 - x_1)}\\
x_4\\
\frac{\begin{array}{l}
-m_2l_1l_2\cos(x_3 - x_1)\left(u_1-u_2+ m_2l_1l_2x_4^2 \sin(x_3 - x_1) - (m_1 + m_2) g l_1 \sin(x_1)   \right)\\
+(m_1+m_2)l_1^2\left( u_2-m_2l_1l_2x_2^2\sin(x_3 - x_1) -m_2gl_2\sin(x_3)
\right)
\end{array}}{(m_1+m_2)m_2l_1^2l_2^2 - m_2^2l_1^2l_2^2 \cos^2(x_3 - x_1)}
\end{array}\right].
%\label{eq:ssmodpsys}
\end{eqnarray*}

\subsection{Population Dynamics}
\label{subsec:popugromo}
We first model dynamics of a single population\index{single population!dynamics}. Then, we use the obtained model to investigate two interacting populations\index{interacting populations}.

Let $x=\phi(t)$ be the population of the given species at a time $t$. The simplest model of the population dynamics of the given species is the one resulting from the hypothesis that the net rate of change of $x$ is proportional to the current value of $x$. Let $r$ be the rate of growth of the population\index{rate of growth!of the population}. Then, the simple model of the population dynamics is
\begin{equation}
\dot{x}=r x.
\label{eq:simoofpopu}
\end{equation}
Solving~(\ref{eq:simoofpopu}) subject to the initial condition, $x(0)=x_0$,
%\begin{equation}
%x(0)=x_0,
%\label{eq:icforpopugr}
%\end{equation}
yields
\begin{equation}
x=e^{rt}x_0.
\label{eq:solsimopop}
\end{equation}
In the above linear model, the growth rate, $r$, is constant. It does not account for the limitation of space, food supply, etc. Therefore, we modify the unrestricted growth rate\index{unrestricted growth rate} $r$ to take into account the fact that the environment of the population can only support a certain number, say $K$, of the species. The constant $K$ is called the {\em carrying capacity\/}\index{carrying capacity} of the environment or the {\em saturation level}\index{saturation level}. Specifically, if $x>K$, then because there is not enough food or space to support $x$, more species will die than are born, which results in a negative growth rate. On the other hand, if $x<K$, then there is extra food and space and the population growth is positive. One function that satisfies the above requirements has the form
\[
r\left( 1-\frac{x}{K}\right).
\]
Using the above model of the growth rate, we obtain the following model\index{model!population growth} of the population growth\index{population growth},
\begin{equation}
\dot{x}=r\left( 1-\frac{x}{K}\right)x.
\label{eq:prelogeq}
\end{equation}
Let $c=r/K$. Then, we represent~(\ref{eq:prelogeq}) as
\begin{equation}
\fbox{$\displaystyle \dot{x}=\left(r-cx \right)x $}
\label{eq:stafoofloeq}
\end{equation}
Equation~(\ref{eq:stafoofloeq}) is known as the {\em Verhulst equation\/}\index{Verhulst equation}\index{Verhulst equation|see{logistic equation}} or the {\em logistic equation}\index{logistic equation}. According to Boyce and DiPrima~\cite[page~60]{boyce97}, Verhulst, a Belgian mathematician, introduced equation~(\ref{eq:stafoofloeq}) in 1838 to model human population growth\index{human population growth}. For reasons somewhat unclear, Verhulst referred to it as logistic growth, hence~(\ref{eq:stafoofloeq}) is often called the logistic equation.

We use the logistic equation to analyze interacting species. We first model two interacting species that do not prey upon each other, but compete for the same resources, for example for a common food supply. Let $x_1$ and $x_2$ be the populations, at time $t$, of the two competing species. We assume that the population of each of the species in the absence of the other is governed by a logistic equation. Thus, a model of two non-interacting species\index{non-interacting species} has the form
\begin{equation}
\left.\begin{array}{ccl}
\dot{x}_1 & = & \left( r_1-c_1x_1\right)x_1\\
\dot{x}_2 & = & \left( r_2-c_2x_2\right)x_2,
\end{array}\right\}
\label{eq:noninspec}
\end{equation}
where $r_1$ and $r_2$ are the unrestricted growth rates of the two populations, while $r_1/c_1$ and $r_2/c_2$ are their saturation levels. When both species are present, each will affect the growth rate of the other. Since $x_1$ and $x_2$ use the same resources, the growth rate of $x_1$ will be lowered by an amount proportional to the size of $x_2$. We modify the growth equation for species $x_1$ as
\begin{equation}
\dot{x}_1  =\left( r_1-c_1x_1 - a_1x_2\right)x_1,
\label{eq:mogrospxo}
\end{equation}
where $a_1$ is a measure by which species $x_2$ lower the growth rate of $x_1$. Similarly, the growth rate of $x_2$ is reduced by $a_2x_1$, where $a_2$ is a measure of the degree by which species $x_1$ lower the growth rate of $x_2$. Therefore, the growth equation for species $x_2$ becomes
\begin{equation}
\dot{x}_2  = \left( r_2-c_2x_2-a_2x_1\right)x_2.
\label{eq:groeqfsxtm}
\end{equation}
In summary, equations~(\ref{eq:mogrospxo}) and~(\ref{eq:groeqfsxtm}) model two species that interact by competing for a common resource like food or other supply.

We now analyze the situation in which one species, the predator, preys upon the other species, the prey, and there is no competition between the predator and prey for a common resource. Let $x_1$ and $x_2$ denote the populations of the prey and predator, respectively. We make the following modeling assumptions:
\begin{enumerate}
\item When there is no predator, that is, $x_2=0$, the prey enjoys an exponential growth governed by
\begin{equation}
\dot{x}_1=r_1x_1,\quad r_1>0.
\label{eq:nopreexpo}
\end{equation}
\item When there is no prey, that is, $x_1=0$, the predator perishes. We model this scenario using the equation
\begin{equation}
\dot{x}_2=-r_2x_2,\quad r_2>0.
\label{eq:nopreyhun}
\end{equation}
\item The larger the population of the predator, the more prey is eaten. The larger the population of the prey, the easier hunting is for the predator. We model the growth rate decrease of the prey by including in~(\ref{eq:nopreexpo}) the term $-dx_1x_2$, $d>0$.  The growth rate increase of the predator is modeled by adding the term $ex_1x_2$ in~(\ref{eq:nopreyhun}), where $e>0$.
\end{enumerate}
The resulting model of the predator-prey relationship\index{predator-prey relationship} is
\begin{equation}
\fbox{$\displaystyle\begin{array}{ccl}
\dot{x}_1 & =& \left( r_1-dx_2\right)x_1\\
\dot{x}_2 & =& \left( -r_2+ex_1\right)x_2
\end{array} $}
\label{eq:predpreye}
\end{equation}
Equations~(\ref{eq:predpreye}) are known as the {\em Lotka-Volterra equations}\index{Lotka-Volterra equations!predator-prey relationship}. They were published by Lotka in 1925 and by Volterra in 1926.

\subsection{The Hypothalamic-Pituitary-Adrenal Axis Model}
\label{subsec:HPA_axis_model}

The hypothalamic-pituitary-adrenal~(HPA)\index{hypothalamic-pituitary-adrenal (HPA) axis}\index{HPA axis}\index{HPA axis|see{hypothalamic-pituitary-adrenal axis}} axis is the body's primary stress management\index{stress management} system. The role of the HPA axis is to maintain bodily functions in the presence of physical and mental stressors. The (HPA) axis is a part of the endocrine system. The endocrine system as well as its subsystem, the HPA axis, uses hormones to communicate between the regions of the body. The regulation of hormones maintains homeostasis\index{homeostasis}\index{homeostasis|see{HPA axis}}---the process by which bodily functions are maintained at a constant level. This leads to a definition of stress as a state of disharmony in which the homeostasis of the organism is threatened.

Another approach to define stress was proposed by McEwen~\cite{mcewen02} in 2002. To define stress McEwen introduced a notion of allostasis\index{allostasis}---the process by which the body functions change in response to surrounding stimuli. The term allostasis is the opposite to the notion of homeostasis. An example of allostais is the fight-or-flight response in which the sympathetic nervous system as well as the HPA axis are involved. Dysregulation of the HPA axis\index{dysregulation!of the HPA axis} is associated with a number of neuroimmune disorders\index{neuroimmune disorders} such as chronic fatigue syndrome~(CFS)\index{chronic fatigue syndrome (CFS)}\index{CFS|see{HPA axis}}\index{CFS}, depression, Gulf War illness~(GWI)\index{Gulf War illness}\index{GWI|see{Gulf War illness}}\index{GWI|see{neuroimmune disorders}}, or posttraumatic stress disorder~(PTSD)\index{posttraumatic stress disorder}\index{PTSD|see{posttraumatic stress disorder}}\index{PTSD|see{neuroimmune disorders}}, among other stress related diseases~\cite{benzvi09}. At present, it is not clear what causes dysregulation of the HPA axis.


The HPA model presented here was developed by Gupta et al.~\cite{gupta07} and modified by Ben-Zvi et al.~\cite{benzvi09}. A simplified schematic diagram of the HPA is shown in Figure~\ref{fig:HPA_pic3}.
\begin{figure}
  \centering
  \includegraphics[width=2.75in]{HPA_pic3.eps}
  \caption{A simplified schematic diagram of the HPA axis.}
  \label{fig:HPA_pic3}
\end{figure}

The HPA axis is responsible for a rapid response to stress stimuli. An activation of the hypothalamus by a stressor causes the release of the corticotropin releasing hormone~(CRH)\index{corticotropin releasing hormone}\index{CRH|see{corticotropin releasing hormone}}. The hypothalamus is the control center of most of the body's hormonal systems. Upon reaching the pituitary gland, the CRH hormone induces the release of the adrenocorticotropic hormone~(ACTH)\index{adrenocorticotropic hormone}\index{ACTH|see{corticotropin releasing hormone}} by the pituitary into the circulation that reaches the adrenal glands that are located on top of the kidneys. The ACTH stimulates the secretion of cortisol\index{cortisol} by the adrenal glands\index{adrenal glands}. The release of cortisol initiates metabolic effects to fight the harmful effects of stress through negative feedback to the hypothalamus and pituitary---see Figure~\ref{fig:HPA_pic3}. Once the state of stress subsides, the concentration  of ACTH and cortisol decreases. In their model, Gupta et al.~\cite{gupta07} also include the dynamics of the glucocorticoid receptor~(GR)\index{glucocorticoid receptor}\index{GR|see{corticotropin releasing hormone}} to model more accurately the HPA axis dynamics.

We use the variables described in Table~\ref{tab:HPA_variables} in the HPA axis modeling equations.
\begin{table}
%\renewcommand{\arraystretch}{1}
\caption{Description of variables in the HPA model.}
\centering
\begin{tabular}{ll}
\hline
Variable & Description\\
\hline
$x_1$ & CRH concentration\\
$x_2$ & ACTH concentration\\
$x_3$ & Free GR concentration\\
$x_4$ & Cortisol concentration\\
$d$ & Unknown input modeling\\
{}& stress action\\
$u$ & Control variable modeling\\
{} & treatment action\\
\hline
\end{tabular}
\label{tab:HPA_variables}
\end{table}

We present the HPA model following the development of Gupta et al.~\cite{gupta07}. In this model, stress, modeled by the signal $d$, stimulates the hypothalamus to secrete the CRH denoted $x_1$. The equation modeling the hypothalamus response to stress is
\begin{equation}
\dot{x}_1 = \left(1-\frac{x_4}{k_{i1}}\right)(k_{c}+d)-k_{cd}x_1,
\label{eq:hypthalamus_behav}
\end{equation}
where the terem $-k_{cd}x_1$ models a constant degradation rate of the CRH in the blood. The term $\left(1-\frac{x_4}{k_{i1}}\right)(k_{c}+d)$ models a circadian production\index{circadian production} and stress being reduced by the inhibition term $\left(1-\frac{x_4}{k_{i1}}\right)$. For small values of $x_4/k_{i1}$, we may write
\[
\left(1-\frac{x_4}{k_{i1}}\right)(k_{c}+d) \approx \frac{k_c + d}{1+\frac{x_4}{k_{i1}}}.
\]
Using the above approximation, we obtain the following equation modeling the hypothalamus dynamics\index{hypothalamus dynamics},
\begin{equation}
\dot{x}_1 = \frac{k_{c}+d}{1+\frac{x_4}{k_{i1}}}-k_{cd}x_1 .
\label{eq:hypothalam_dynamics}
\end{equation}
The adrenocorticotropic hormone~(ACTH) produced by the pituitary gland\index{pituitary gland} is modeled as
\begin{equation}
\dot{x}_2=\frac{x_1}{1+\frac{x_1 x_4}{k_{i2}}}-k_{ad}x_2 .
\label{eq:pituitary_dyna}
\end{equation}
The dynamics of the glucocorticoid receptor~(GR) is modeled by Gupta et al.~\cite{gupta07} as
\begin{equation}
\dot{x}_3 = \frac{(x_1 x_4)^2}{k+(x_3 x_4)^2}-k_{rd}x_3 .
\label{eq:gr_concen_dyn}
\end{equation}
Finally, the production of the cortisol\index{cortisol} by the adrenal glands is modeled as
\[
\dot{x}_4=x_2-x_4.
\]
Ben-Zvi et al.~\cite{benzvi09} propose to augment the above equation with the control term $u$ modeling the control action,
\begin{equation}
\dot{x}_4=x_2-x_4 + u.
\label{eq:corti_adrenal}
\end{equation}
Combining equations~\eqref{eq:hypothalam_dynamics},~\eqref{eq:pituitary_dyna},~\eqref{eq:gr_concen_dyn}, and~\eqref{eq:corti_adrenal}, we obtain a fourth-order model of the HPA axis\index{HPA axis!model of},
\begin{equation}
\left[\begin{array}{c}
\dot{x}_1\\
\dot{x}_2\\
\dot{x}_3\\
\dot{x}_4
\end{array}\right]=\left[\begin{array}{c}
\frac{1}{1+\frac{x_4}{k_{i1}}}-k_{cd}x_1\\
\frac{x_1}{1+\frac{x_1 x_4}{k_{i2}}}-k_{ad}x_2\\
\frac{(x_1 x_4)^2}{k+(x_3 x_4)^2}-k_{rd}x_3\\
x_2-x_4
\end{array}\right]+ \left[\begin{array}{c}
0\\
0\\
0\\
1
\end{array}\right]u + \left[\begin{array}{c}
\frac{1}{1+\frac{x_4}{k_{i1}}}\\
0\\
0\\
0
\end{array}\right]d .
\label{eq:hpa_model}
\end{equation}
The normalized parameter values listed in Ben-Zvi et al.~\cite{benzvi09} are given in Table~\ref{tab:HPA_parameters}.
\begin{table}
%\renewcommand{\arraystretch}{1.3}
\caption{Parameter values in the HPA model}
\centering
\begin{tabular}{lll}
\hline
Parameter & Description & Value\\
\hline
$k_{i1}$ & Inhibition constant for CRH synthesis & 0.100\\
$k_{cd}$ & CRH degradation constant & 1.000\\
$k_{i2}$ & Inhibition constant for ACTH synthesis & 0.100\\
$k_{ad}$ & ACTH degradation constant & 10.000\\
$k_c$ & circadian production & 1\\
$k_{cr}$ & GR synthesis constant & 0.050\\
$k_{rd}$ & GR degradation constant & 0.900\\
$k$ & Inhibition constant for GR synthesis & 0.001\\
\hline
\end{tabular}
\label{tab:HPA_parameters}
\end{table}



%Following the approach of Ben-Zvi et al.~\cite{benzvi09}, we obtain steady-state values of the state variables as a function of the external stressor, $d$. The obtained plots are shown in Figures~\ref{fig:x1_and_x3_vs_d} and~\ref{fig:x2_and_x4_vs_d}.
%\begin{figure}
%  \centering
%  \includegraphics[width=3.25in]{x1_and_x3_vs_d}
%  \caption{Plots of the steady-state values of $x_1$ and $x_3$ versus $d$.}
%  \label{fig:x1_and_x3_vs_d}
%\end{figure}
%
%\begin{figure}
%  \centering
%  \includegraphics[width=3.25in]{x2_and_x4_vs_d}
%  \caption{Plots of the steady-state values of $x_2$ and $x_4$ versus $d$.}
%  \label{fig:x2_and_x4_vs_d}
%\end{figure}
%
%We note that in a chronically stressed individual, cortisol concentration, $x_4$, is very low. Thus a healthy individual subjected to a prolonged extreme stress, $d>0.168$, would settle down in a stable equilibrium state corresponding to depressed cortisol concentration, $x_4$, corresponding to the lower branch of the curve in Figure~\ref{fig:x2_and_x4_vs_d}. When the stress subsides, that is, $d=0$, the individual will stay in the new equilibrium state corresponding to a depressed cortisol concentration. This is because the equilibrium corresponding to $d=0$ is asymptotically stable and so states ``close'' to it will be attracted by this asymptotically stable low cortisol equilibrium.  Ben-Zvi et al.~\cite{benzvi09} propose a treatment strategy whereby the chronically stressed individual is moved to the healthy state corresponding to the normal cortisol concentration, by moving the individual into a negative stress region, after which internal regulatory processes translates the individual into the healthy equilibrium which is asymptotically stable and maintains it there. Effective treatment can be accomplished when the states as well as stress levels are available. Thus a method to estimate state variables as well as the stress level are highly desirable.

\section{Notes}

For more examples of dynamical system models from mechanical and electrical engineering, the reader is referred to Truxal~\cite{truxal72}, Banks~\cite{banks86}, Raven~\cite{raven87}, Mohler~\cite{mohler91}, and D{'}Azzo and Houpis~\cite{dazzo95}. Wertz~\cite{wertz97} has many excellent real-life models from space science. A well-written paper on vehicle modeling\index{vehicle modeling} and control for broad audience is by Ackermann~\cite{ackermann97}. In addition to references~\cite{wong93} and~\cite{will97}, vehicle and tire dynamics\index{tire dynamics} modeling are investigated in detail by Nalecz~\cite{nalecz88}, Nalecz and Bindemann~\cite{nalecz88b,nalecz88c}, and Smith and Starkey~\cite{smith94,smith95}. For more on stepper motors, the reader is referred to Acarnley~\cite{acarnley82}, Kenjo~\cite{kenjo84}, Kenjo and Nagamori~\cite{kenjo85}, Del~Toro~\cite[Chapter 9]{deltoro85}, Krause and Wasynczuk~\cite[Chapter 8]{krause89}, and Zribi and Chiasson~\cite{zribi91}. For further analysis of the equations modeling competing species and the predator-prey equations, we recommend Boyce and DiPrima~\cite[Sections~9.4 and 9.5]{boyce97}. Our derivation of the logistic equation follows that of Sandefur~\cite{sandefur90}. The logistic equation as well as equations modeling interacting species can be used to mathematical modeling of tumor dynamics\index{tumor dynamics} and interactions between tumor and immune system\index{immune system}. A survey of models for tumor-immune system dynamics can be found in the well-edited book by Adam and Bellomo~\cite{adam97}.

%Throughout the previous decades much effort has been devoted to the mathematical modeling of the hypothalamic-pituitary-adrenal (HPA) axis to better understand its functioning.
A number of mathematical models of the hypothalamic-pituitary-adrenal~(HPA) axis\index{HPA axis} were proposed in the last six decades. Building on the previous models, Bingzhen, Zhenye, and Lainsong~\cite{bingzheng90} proposed a third-order dynamical model and tested it on clinical data. This model was later improved by Liu et al.~\cite{liu99HPA}. Kyrylov, Severyanova and Vieira~\cite{kyrylov05} modified the model of Liu et al.~\cite{liu99HPA} to include additional properties of the HPA axis. Conrad et al.~\cite{conrad09} proposed and analyzed a second-order non-linear mathematical model of the HPA axis. Lenbury and Pornsawad~\cite{lenbury05} used delay-differential equations  in their HPA axis model to account for the delays associated with the action of the glands in response to the stimulating hormones. Gupta et al.~\cite{gupta07}, on the other hand, proposed a fourth-order non-linear state-space model of the HPA axis in which the effect of stress on the HPA axis dynamical behavior is modeled. Ben-Zvi, Vernon, and Broderick~\cite{benzvi09} modified the model of Gupta et al.~\cite{gupta07}  by adding the control input representing the treatment. The HPA model proposed by Ben-Zvi et al.~\cite{benzvi09} can be viewed as a dynamical system with unknown input\index{unknown input}.



\section{Exercises}
\begin{exercise}\label{ex:kplusuconst}\rm
Prove the law of conservation of energy as stated in the following theorem:
\begin{thm}
Given conservative vector field $\b F$, that is, $\b F=-\nabla U$, where $U$ is the potential energy. Assume that a particle of mass $m$, satisfying Newton's law, moves on a differentiable curve $\b C$. Then, the sum of the potential energy and kinetic energy is constant.
\end{thm}
{\bf Hint}\quad Denote by $\b C(t)$ a differentiable curve along which a particle of mass $m$ moves. Hence
\[
\b F(\b C(t))=m\ddot{\b C}(t).
\]
Then, show that
\[
U(\b C(t))+\frac{1}{2}m\left\|\dot{\b C}(t)\right\|^2=\mb{constant}.
\]
It is enough to differentiate $U(\b C(t))+\frac{1}{2}m\left\|\dot{\b C}(t)\right\|^2$ with respect to time and show that the time derivative is zero, hence the sum of the energies is constant.
\end{exercise}
\begin{exercise}\label{ex:potenergy}\rm
Show that in a conservative vector field the work done by the field on a particle that moves between two points depends only on these points and not on the path followed.\\
{\bf Hint}\quad Note the work done against the conservative force field $\b F$ along a curve $\b C$ between points $\b C(t_1)$ and $\b C(t_2)$ is
\[
\int_{t_1}^{t_2}\b F(\b C(t))^{\top}\dot{\b C}(t)dt.
\]
By assumption $\b F$ is a conservative vector field. Hence,
\begin{equation}
\b F=-\nabla U.
\label{eq:conup}
\end{equation}
Combining the above yields
\[
\int_{t_1}^{t_2}\b F(\b C(t))^{\top}\dot{\b C}(t)dt=-\int_{t_1}^{t_2}\nabla U(\b C(t))^{\top}\dot{\b C}(t)dt.
\]
Use the chain rule to represent the expression inside the integral on the right hand side as
\[
\nabla U(\b C(t))^{\top}\dot{\b C}(t)=\frac{d U(\b C(t))}{dt}.
\]
\end{exercise}

\begin{exercise}\label{ex:twopumps}\rm
The goal of this exercise is to develop a model of a two-tank pumped storage system. Figure~\ref{fig_pump} illustrates such a system.
\begin{figure}
\centering
\includegraphics[scale=0.6]{PUMP.eps}
%%\vspace{3.5in}
%\centerbmp{3.33in}{2.5in}{c:/mybook/pump.bmp}
%%\hskip1in\special{bmp:c:/mybook/pump.bmp x=4.67in y=3.5in}
\caption{A two-tank pumped storage system for Exercise~\ref{ex:twopumps}.}
\label{fig_pump}
\end{figure}
A pumped storage system can be used to supply water to a turbine producing electric power. The two storage tanks are connected together through a valve that may be modeled as a linear resistance $R_1$. The resistance of the turbine supply pump is modeled by a linear resistance $R_2$. Suppose that we wish to design a controller that would control the inflow rate $q_i$ in order to maintain the output flow rate $q_o$ at some desired value $q_d$. The first stage of our project is to develop a design model of the pumped storage system. We assume that linearized flow rates $q(t)$ and $q_o(t)$ are proportional to the corresponding water heights, that is,
\begin{eqnarray*}
q(t) &=& \frac{h_1(t)-h_2(t)}{R_1},\\
q_o(t) &=& \frac{h_2(t)}{R_2}.
\end{eqnarray*}
Let $A_1$ and $A_2$ represent the surface areas in the two tanks. Then, the principle of conservation of mass for an incompressible fluid yields
\begin{eqnarray*}
q_i(t)-q(t) &=& A_1\dot{h}_1(t),\\
q(t)-q_o(t) &=& A_2\dot{h}_2(t).
\end{eqnarray*}
Numerical data of a two-tank pumped storage system is given in Table~\ref{table_pump}.
\begin{table}
\caption{Numerical data for the two-tank pumped storage system model for Exercise~\ref{ex:twopumps}.}
\label{table_pump}
\begin{center}
\begin{tabular}{cl}\hline
Tank 1 diameter & 2\,m\\
Tank 2 diameter & 3\,m\\
$R_1$ & 10\,sec/m$^2$\\
$R_2$ & 20\,sec/m$^2$\\ \hline
\end{tabular}
\end{center}
\end{table}
Find the transfer function ${\cal L}(q_o(t))/{\cal L}(q_i(t))$.
\end{exercise}
\begin{exercise}\label{ex:vehsusex}\rm
The objective of this exercise is to derive a model of an automotive suspension.
The role of an automotive suspension is to support the vehicle body on the axles. A simple linear model of the suspension assembly associated with a single wheel, called quarter-car model, is shown in Figure~\ref{fig:actsus}. The vehicle body is represented by the ``sprung'' mass, $m_s$, while the tire and axle by the ``unsprung'' mass, $m_u$.
\begin{figure}
\centering
\includegraphics[scale=0.6]{ACTSUS.eps}
%%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/sliding/actsus.bmp}
%%\hskip1.25in\special{bmp:c:/sliding/actsus.bmp x=4.67in y=3.5in}
\caption{Linear quarter-car model for Exercise~\ref{ex:vehsusex}.}
\label{fig:actsus}
\end{figure}
Model the tire as a spring with spring coefficient $k_t$, and assume its damping to be negligible. The suspension consists of the spring, shock absorber, and the variable force element. The role of the variable force element is to generate a force that compensates for the uneven road which is the source of the system disturbances. Assume that the variable force element $F_a$ can instantaneously provide any desired force. The above described  model of an automotive suspension is referred to as a linear active suspension. Write down equations modeling the dynamics of the two degree of freedom system depicted in Figure~\ref{fig:actsus} by applying Newton's law to each mass of the system, and taking into account the positive direction of the $z$ coordinates. Define the state variables to be:
\begin{eqnarray*}
x_1 &=& z_s-z_u=\mb{suspension stroke},\\
x_2 &=& \dot{z}_s =\mb{sprung mass velocity},\\
x_3 &=& z_u-z_r =\mb{tire deflection},\\
x_4 &=& \dot{z}_u=\mb{unsprung mass velocity}.
\end{eqnarray*}
Let $u=F_a$ be the system input, and let $d(t)=\dot{z}_r(t)$ be the disturbance input due to the road roughness. Represent the obtained linear model of active suspension in state-space format.
Use data of a commercial vehicle given in Table~\ref{tab:actsus}. The data comes from Sunwoo and Cheok~\cite{sunwoo90}.
\begin{table}
\caption{Parameters of the vehicle suspension of Exercise~\ref{ex:vehsusex}.}
\begin{center}
\begin{tabular}{cc}\hline
$m_s$ & 355 kg\\
$m_u$ & 46 kg\\
$k_s$ & 14384 N/m\\
$b_s$ & 1860 N/m/sec\\
$k_t$ & 300000 N/m\\ \hline
\end{tabular}
\end{center}
\label{tab:actsus}
\end{table}
Substitute the vehicle data into the modeling equations. A design goal of an automotive engineer may be to design a control law to isolate the sprung mass from road disturbances. In other words, he or she may wish to design a control law such that the vertical acceleration of the sprung mass $\ddot{z}_s(t)=\dot{x}_2$ is as close to $0$ as possible in spite of road disturbances. Assume that we can measure the vertical acceleration of the sprung mass. Write the output equation.
\end{exercise}

\begin{exercise}\label{ex:machtoo}\rm
The objective of this exercise is to develop a model of a single axis machine tool drive in which a motor drives a lead screw through a gear. Figure~\ref{fig_mach_dr} illustrates such a system.
\begin{figure}
\centering
\includegraphics[scale=0.7]{MACH_DR.eps}
%%\vspace{4.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/mach_dr.bmp}
%%\special{bmp:c:/mybook/mach_dr.bmp x=6in y=4.5in}
\caption{A schematic of a machine tool drive system for Exercise~\ref{ex:machtoo}.}
\label{fig_mach_dr}
\end{figure}
This model is adapted from Driels~\cite[pages 613--619]{driels96}. The system operates as follows. The lead screw moves the table with a workpiece on it under the cutter to perform the desired operation. The table is positioned by a control system. We wish to place the table in a desirable position as rapidly and accurately as possible. This means, that while modeling the machine tool drive system we should, in particular, take into account the flexibility of the lead screw. The length of the lead screw implies that its rotation at the gear box will not be the same as its rotation at the table. Thus, the angular twist along the length of the lead screw cannot be neglected. Model the lead screw flexibility using a torsional spring as illustrated in Figure~\ref{fig_ma_dr_s}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{MA_DR_S.eps}
%%\vspace{4.5in}
%\centerbmp{4.67in}{3.5in}{c:/mybook/ma_dr_s.bmp}
%%\special{bmp:c:/mybook/ma_dr_s.bmp x=6in y=4.5in}
\caption{A graphical model of the machine tool drive system for Exercise~\ref{ex:machtoo}.}
\label{fig_ma_dr_s}
\end{figure}
The system input is the motor torque $\tau_m$, and the system output is the rotation of the end of the lead screw labeled by $\theta_l$. Assume that the damping torque, due to viscous friction, acting on the motor is proportional to its angular velocity $\dot{\theta}_m=\omega_m$ with the coefficient of viscous friction denoted $c_m$. Denote the motor moment of inertia by $I_m$, and by $\tau_1$ the load torque on the motor gear due to the rest of the gear train. Write the equation modeling the dynamics of the motor shaft.

Let $\tau_2$ be the torque transmitted to the load shaft, and let $N=N_1/N_2$ be the gear ratio, where $N_1$ denotes the number of teeth of the first gear, and respectively $N_2$ is the number of teeth of the second gear. The work done by the first gear is equal to that of the second. Model the load shaft flexibility with a torsional spring. The torsional spring coefficient is $k$. Assume that the moment of inertia of the lead screw is negligible compared to that of the load. Denote the moment of inertia of the load by $I_l$. The torque due to viscous friction of the load is $c_l\omega_l$, where $c_l$ is the coefficient of the viscous friction of the load. The angular displacement of the end of the lead screw is $\theta_l$. Write the equation modeling the dynamics of the load.
Define the state variables to be
\[
x_1=\theta_m,\quad x_2=\dot{\theta},\quad x_3=\theta_l,\quad x_4=\dot{\theta}_l.
\]
Represent the modeling equations in state space format. Typical numerical data for a machine tool drive system are shown in Table~\ref{table_mach_dr}.
\begin{table}
\caption{Data for the machine tool drive system model for Exercise~\ref{ex:machtoo}.}
\label{table_mach_dr}
\begin{center}
\begin{tabular}{ll}\hline
$I_m$ & 0.0001 Nm$\cdot$s$^2$/rad\\
$I_l$ & 0.001 Nm$\cdot$s$^2$/rad\\
$c_m=c_l$ & 0.01 Nm$\cdot$s/rad\\
$k$ & 1 Nm/rad\\
$N$ & 1/5\\ \hline
\end{tabular}
\end{center}
\end{table}
\end{exercise}
\begin{exercise}\label{ex:bicytwoeqmo}\rm
Consider a bicycle model of a front wheel steering vehicle shown in Figure~\ref{fig:fsbikem}. Construct a two-dimensional state-space model of the vehicle, where $x_1=\dot{y}$ and $x_2=\dot{\theta}$. The input $u=\delta_f$ and the output $y=x_2$. In your derivation assume small angles. Vehicle data are given in Table~\ref{tab:fsuvehd}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{FSBIKEM.eps}
%%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/ee680/fsbikem.bmp}
%%\hskip.5in\special{bmp:c:/ee680/fsbikem.bmp x=4.67in y=3.5in}
\caption{Bicycle model of a vehicle for Exercise~\ref{ex:bicytwoeqmo}.}
\label{fig:fsbikem}
\end{figure}
\begin{table}
\caption{Vehicle data for Exercise~\ref{ex:bicytwoeqmo}.}
\begin{center}
\begin{tabular}{cc}\hline
$l_1=l_2$ & 1 m\\
$v_x$ & 10 m/sec\\
$m$ & 1,000 kg\\
$I_z$ & 1,000 kg$\cdot$m$^2$\\
$C_f=C_r$ & 10,000 N/rad\\ \hline
\end{tabular}
\end{center}
\label{tab:fsuvehd}
\end{table}
\end{exercise}
\begin{exercise}\label{ex:twmaspwl}\rm
Consider the two-mass-spring system depicted in Figure~\ref{fig:a_95_f1}.
\begin{figure}
\centering
\includegraphics[scale=0.7]{A_95_F1.eps}
%%\vspace{4.5in}
%\centerbmp{4.67in}{3.5in}{c:/exams/a_95_f1.bmp}
%%\special{bmp:c:/exams/a_95_f1.bmp x=6in y=4.5in}
\caption{The two-mass-spring system for Exercise~\ref{ex:twmaspwl}.}
\label{fig:a_95_f1}
\end{figure}
\begin{itemize}
\item[(i)] Form the Lagrangian of the system.
\item[(ii)] Write the Lagrange equations of motions for the system.
\item[(iii)] Represent the obtained Lagrange equations of motion in state-space format.
\end{itemize}
\end{exercise}
\begin{exercise}\label{ex:maspdalag}\rm
Consider the mechanical system shown in Figure~\ref{fig_test1f1}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{TEST1F1.eps}
%%\vspace{3in}
%\centerbmp{4in}{3in}{c:/ee382/test1f1.bmp}
%%\hskip0.8in\special{bmp:c:/ee382/test1f1.bmp x=4in y=3in}
\caption{Mechanical system for Exercise~\ref{ex:maspdalag}.}
\label{fig_test1f1}
\end{figure}
\begin{itemize}
\item[(i)] Form the Lagrangian for the system.
\item[(ii)] Write the Lagrange equation of motion for the system.
\item[(iii)] Represent the system model in state space format.
\end{itemize}
\end{exercise}
\begin{exercise}\label{ex:lumppenex}\rm
Consider a dynamical system that consists of a cart with an inverted pendulum (point mass on a massless shaft) attached to it as depicted in Figure~\ref{fig_inv_pend}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{INV_PEND.eps}
%%\vspace{4in}
%\centerbmp{4.67in}{3.5in}{c:/ee675/inv_pend.bmp}
%%\hskip.5in \special{bmp:c:/ee675/inv_pend.bmp x=5.33in y=4in}
\caption{An inverted pendulum moving on a cart for Exercise~\ref{ex:lumppenex}.}
\label{fig_inv_pend}
\end{figure}
\begin{itemize}
\item[(i)] Write the Lagrange equations of motion for the system,
\item[(ii)] Represent the obtained model in state-space format using the state variables
\[
x_1=x,\quad x_2=\dot{x},\quad x_3=\theta,\quad x_4=\dot{\theta}.
\]
\end{itemize}
\end{exercise}
\begin{exercise}\label{ex:fieldmoex}\rm
Consider a schematic of a field-controlled DC motor shown in Figure~\ref{fig_dcmot_fc}.
\begin{figure}
\centering
\includegraphics[scale=0.6]{DCMOT_FC.eps}
%%\vspace{3.5in}
%\centerbmp{4in}{3in}{c:/ee675/dcmot_fc.bmp}
%%\hskip1in\special{bmp:c:/ee675/dcmot_fc.bmp x=4.67in y=3.5in}
\caption{A diagram of a field-controlled DC motor for Exercise~\ref{ex:fieldmoex}.}
\label{fig_dcmot_fc}
\end{figure}
Let $x_1=\theta$, $x_2=\dot{\theta}$, $x_3=i_f$, and $y=x_1$. Treat the field inductance $L_f$ as a parasitic element and denote it $\mu$.
\begin{itemize}
\item[(i)] Form the actual plant model $P_{\mu}$ in state space format;
\item[(ii)] Find the reduced-order design model $P_0$.
\end{itemize}
\end{exercise} 