\documentclass[12pt]{article}
\usepackage{latexsym}
\usepackage{graphicx}
\usepackage[hypertex]{hyperref}
\usepackage{epsfig}
\usepackage{subfigure}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{mathrsfs}
\usepackage{cases}
\usepackage{array}
\renewcommand{\baselinestretch}{1.17}
\setlength{\oddsidemargin}{0in}
\setlength{\evensidemargin}{0in}
\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{0in}
\setlength{\textheight}{8.5in}
%\input setbmp
%\input{psfig.sty}
\usepackage{epsfig}
\newtheorem{defn}{Definition}
\newtheorem{thm}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{example}{Example}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}
\newtheorem{algorithm}{Algorithm}
\newtheorem{exercise}{Exercise}
%\usepackage[centertags]{amsmath}
\newcommand{\markex}{\nopagebreak\begin{center}\rule{5in}{.01in}\end{center}}
\newcommand{\markdefn}{\nopagebreak\hfill$\Diamond$}
\newcommand{\mb}{\mbox}
\newcommand{\bm}{\boldmath}
\def\b#1{\mb{\bm$#1$}}
\newfont{\Bb}{msbm10 scaled\magstep1}
\newfont{\Fr}{eufm10 scaled\magstep1}
\newfont{\Sc}{eusm10 scaled\magstep1}
%\pagestyle{empty}
\begin{document}
\noindent{\large ECE 602}\hfill{\Large Fall 2010}
\begin{center}
{\Large \bf Dynamic System Modeling}\\
by\\
{\large Stanislaw H. \.{Z}ak}
\end{center}
\begin{center}
September 07, 2010
\end{center}
\bigskip
\begin{flushright}
\parbox{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}

\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.

%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 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 course 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, 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 depicted 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 dynamic systems.

For the controller design purposes, we first construct a truth model of the given plant and then its design model.


\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 $m$ subjected to a force $\mb{\bm$F$}$. Then, by 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 $\delta W$ done by $\mb{\bm$F$}$ along 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. 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\/}. 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 now review the notion of {\em 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 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}$.


In a conservative vector field 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 now represent Newton's equation in an equivalent format that establishes a connection with the Lagrange equations of motions 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\/} or just the 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 motions, in 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}\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 there are $3p$ degrees 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 will use the letter $q$ as a symbol for a coordinate regardless of its nature. A generalized coordinate may represent, for example, angular displacement or electric charge. This indicates that the results we are going to obtain will be of general nature and will 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,
\]
etc. We assume that the functions $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.

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 will 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} 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 $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 will need the notion of {\em 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 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 is not constant, but changes by the amount of work done by a nonconservative force according to the relation~(\ref{eq:nonconser}). Let us now denote by $E$ the final 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 is transformed into {\em 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}. 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 non conservative forces acting on the particle.
\begin{example} \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 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*}
The Lagrange equation describing the pendulum motion is
\[
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}} \right) - \frac{\partial L}{\partial \theta} = 0.
\]
In our example the above equation evaluates to
\[
12 \ddot{\theta}+ 60 \sin \theta =0,
\]
or equivalently
\[
\ddot{\theta}=-5 \sin \theta.
\]
\end{example}
\markex

\section{Double Pendulum Modeling}
\label{subsection:two_link_robot}
We construct a model of the double pendulum shown in Figure~\ref{fig:dbl_pend}.
\begin{figure}
\centering
\includegraphics[scale=0.5]{dbl_pend.eps}
%\vspace{3.5in}
%\centerbmp{4.67in}{3.5in}{c:/ee675/dbl_pend.bmp}
%\hskip1in\special{bmp:c:/ee675/dbl_pend.bmp x=4.67in y=3.5in}
\caption{Double pendulum.}
\label{fig:dbl_pend}
\end{figure}
This double pendulum can be viewed as a model of a double-link robot manipulator. 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 $L$ is defined as $L=K-U$, where $K$ is the kinetic energy of the double pendulum 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_1\|\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
\begin{eqnarray*}
\|\b v_{m_1}\|^2 &=& \dot{x_1}^2 + \dot{y_1}^2=\left( l_1\dot{\theta}_1\right)^2\\
\|\b v_{m_2}\|^2 &=& \dot{x_1}^2 + \dot{y_1}^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_1 - \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_1\|\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 pendulum system 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_1 - \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_1 - \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 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_1 - \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.
\]
Hence we can solve~(\ref{eq:matftwlags}) for $\ddot{\theta}_1$ and $\ddot{\theta}_2$. We have
\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*}

%\section{Summary and Notes}
%
%
%
%\section{Exercises}
%\begin{exercise}\label{ex:dist_eigenval_pf}\rm Prove Theorem~\ref{thm:dist_eigenval}.
%\end{exercise}

%\begin{exercise}\label{ex:toiletsysex}\rm
%An example of a closed-loop feedback system is a toilet flushing device connected to a water tank. A toilet system is shown in Figure~\ref{fig:toiletd}. The control objective is to maintain the level of water in the tank at a constant level. Draw a block diagram of the system.
%\begin{figure}[h]
%\centering
%\includegraphics[scale=0.6]{TOILETD.eps}
%%\vspace{3.5in}
%%\hskip1in\special{bmp:c:/mybook/toiletd.bmp x=4.67in y=3.5in}
%\caption{Toilet-flushing system for Exercise~\ref{ex:toiletsysex}.}
%\label{fig:toiletd}
%\end{figure}
%\end{exercise}



\bibliography{root2}
\bibliographystyle{plain}
\end{document} 