%---TeX-start-of-header---
\documentclass{article}
\usepackage{epsfig,mathematica,figi}

\mmanobreak

\begin{comment}
\begin{mathematica}
ohide = DisplayFunction->Identity;
oshow = DisplayFunction -> $DisplayFunction;
ojoin = PlotJoined -> True;
oall = PlotRange->All;
osurf = HiddenSurface->False;
\end{mathematica}
\begin{mathematica}
<<tek.m;
.
 -- Tektronix graphics initialized --
\end{mathematica}

\end{comment}

\begin{document}
\psfiginit
%---TeX-end-of-header---

\begin{titlepage}

\author{Tomasz J.\ Cholewo}
\title{EE\,665 Computer Project}
\date{\today}
\maketitle
\vfil
\tableofcontents
\vfil
\thispagestyle{empty}
\end{titlepage}
%====================================================================
\section{Controller design}
%--------------------------------------------------------------------
\subsection{Linearization}
The given nonlinear system has the following form:
\begin{eqnarray}
\nonumber  \dot{x_1} &=& x_1 - x_2 + u \\
\nonumber  \dot{x_2} &=& -x_1 - x_2^3 + u \\
\nonumber  y &=& x_1 - x_1 x_2^2
\end{eqnarray}
After linearization this system can be written using matrices:
\begin{eqnarray}
\nonumber  \mathbf{\dot{x}} &=& \mathbf{A x + B u + f(x, u)}  \\
\nonumber  \mathbf{y} &=& \mathbf{C x + D u + g(x, u)}
\end{eqnarray}
In my notation \texttt{S} and \texttt{T} denote system equations and
output vectors, accordingly:
\begin{mathematica}
x = {x1[t], x2[t]};
S = {x1[t] - x2[t] + u[t], -x1[t] - x2[t]^3 + u[t] };
T = x1[t] - x1[t] x2[t]^2;
\end{mathematica}
  Matrix \texttt{A} is equal to the jacobian evaluated at the point $(0,0)$.
\begin{mathematica}
(A = Outer[D, S, x] /. {x1[t] -> 0, x2[t] -> 0, u[t] -> 0}) // MatrixForm
.
//MatrixForm= 1    -1

              -1   0
\end{mathematica}
  Vector \texttt{B} can be calculated as a simple derivative:
\begin{mathematica}
B = D[S, u[t]] /. {x1[t] -> 0, x2[t] -> 0, u[t] -> 0}
.
= {1, 1}
\end{mathematica}
  Vector \texttt{f} represents higher order terms left after linearization.
\begin{mathematica}
f = S - A . x /. u[t] -> 0
.
            3
= {0, -x2[t] }
\end{mathematica}
  Now let's calculate a vector \texttt{c}:
\begin{mathematica}
c = D[T, #]& /@ {x1[t], x2[t]} /. {x1[t] -> 0, x2[t] -> 0}
.
= {1, 0}
\end{mathematica}
  a scalar \texttt{d}:
\begin{mathematica}
d = D[T, u[t]]
.
= 0
\end{mathematica}
and a nonlinear term \texttt{g}:
\begin{mathematica}
g = T - c . x /. u[t] -> 0
.
               2
= -(x1[t] x2[t] )
\end{mathematica}
  Let us check if the given system is really unstable because maybe
  there is no need for additional work :-).
\begin{mathematica}
Eigenvalues[A] // N
.
= {-0.618034, 1.61803}
\end{mathematica}
  Unfortunately, one of the eigenvalues is positive and it makes this system
  unstable.


%--------------------------------------------------------------------
\subsection{Feedback}

  The system is controllable if matrix $(\mathbf{B, AB})$ has full rank.
\begin{mathematica}
Det[ {B, A . B} ] != 0
.
= True
\end{mathematica}
  It follows that it is possible to set the eigenvalues of the closed
  loop system matrix $\mathbf{A - BK}$ to the desired $-1$ and $-1$.
  The linear feedback matrix, assuming that we can observe $x_1$ and
  $x_2$, can be determined as:
\begin{mathematica}
K = {k1, k2} /. Solve[
    Eigenvalues[ A - Outer[Times, B, {k1, k2}] ] == {-1, -1}] [[1]]
.
= {8, -5}
\end{mathematica}
%--------------------------------------------------------------------
\subsection{Luenberger Observer}
  Let us check whether the pair $\mathbf{(A, C)}$ is observable.
\begin{mathematica}
Det[{c, c . A}] != 0
.
= True
\end{mathematica}
  So, it is observable and we can move both eigenvalues to the point $-2$
  properly choosing vector \texttt{F}:
\begin{mathematica}
F = {f1, f2} /. Solve[
    Eigenvalues[A - Outer[Times, {f1, f2}, c]] == {-2, -2}] [[1]]
.
= {5, -5}
\end{mathematica}

  We can ``manually'' check that the eigenvalues of the
  linearized system are like expected and that this matrix is Hurwitz.
\begin{mathematica}
Eigenvalues[{ {-7,4,8,-5}, {-9,5,8,-5}, {0,0,-4,-1}, {0,0,4,0}}]
.
= {-2, -2, -1, -1}
\end{mathematica}

From the corollary on page~91 of notes it follows that the closed-loop
system equilibrium $(0, 0)$ is locally asymptotically stable.

%====================================================================
\section{Simulation}
  Function \texttt{Simulation} solves the system of differential
  equations for given initial conditions and time interval.
\begin{mathematica}
e = {e1[t], e2[t]};
BK = Outer[Times, B, K];
Fc = Outer[Times, F, c];
Simulation[x0_List, e0_List, tmax_] := NDSolve[
    Flatten[{
        MapThread[Equal, {{x1'[t], x2'[t]}, (A - BK) . x + BK . e + f}],
        x1[0] == x0[[1]], x2[0] == x0[[2]],
        MapThread[Equal, {{e1'[t], e2'[t]}, (A - Fc) . e + f - F g}],
        e1[0] == e0[[1]], e2[0] == e0[[2]]
    }],
    {x1, x2, e1, e2},
    {t, 0, tmax}, MaxSteps->Infinity
];
\end{mathematica}

%====================================================================
\section{Results}
  The following function plots one of the states of the system versus time
  for initial conditions $x_{10} = x_0$, and $x_{20} = e_{10} = e_{20} =0$.
\begin{mathematica}
plotx[x_, x0_, tmax_] := Plot[
    Evaluate[ x[t] /. Simulation[ {x0, 0.0}, {0, 0}, tmax] ],
    {t, 0, tmax}, PlotRange->All
];
\end{mathematica}
  Function \texttt{plotp} creates a parametric plot in the phase plane
  $x_1$--$x_2$.
\begin{mathematica}
plotp[x0_, tmax_] := ParametricPlot[
    Evaluate[{x1[t], x2[t]} /. Simulation[{x0, 0}, {0,0}, tmax]],
    {t, 0, tmax}, PlotRange->All
];
\end{mathematica}
  Plots for different initial conditions are generated below
  and shown in following figures.
  First let's plot a stable case.
\begin{mathematica}
Display["x1s.eps", plotx[x1, 0.103, 10], "EPS"];
Display["x2s.eps", plotx[x2, 0.103, 10], "EPS"];
Display["e1s.eps", plotx[e1, 0.103, 10], "EPS"];
Display["e2s.eps", plotx[e2, 0.103, 10], "EPS"];
Display["ps1.eps", plotp[0.103, 100], "EPS"];
Display["ps2.eps", plotp[0.01, 100], "EPS"];
\end{mathematica}
\pagevi{x1s}{x2s}{e1s}{e2s}{ps1}{ps2}{Plots for the initial point within the
  domain of attraction: (a)--(d)~plots of $x_1$, $x_2$, $e_1$ and $e_2$
  versus time for time interval $[0,10]$ and $x_{10} = 0.103$;
  (e)~analogous phase-plane plot; (f)~phase-plane plot for $x_{10} =
  0.01$. ($x_{20} = e_{10} = e_{20} = 0$.)}
  Below plots for unstable solutions are calculated.
\begin{mathematica}
Display["x1u.eps", plotx[x1, 0.105, 6], "EPS"];
Display["x2u.eps", plotx[x2, 0.105, 6], "EPS"];
Display["e1u.eps", plotx[e1, 0.105, 6], "EPS"];
Display["e2u.eps", plotx[e2, 0.105, 6], "EPS"];
Display["pu1.eps", plotp[0.105, 2], "EPS"];
Display["pu2.eps", plotp[0.105, 6], "EPS"];
\end{mathematica}
\pagevi{x1u}{x2u}{e1u}{e2u}{pu1}{pu2}{Plots for the initial point not
  within the domain of attraction: (a)--(d)~plots of $x_1$, $x_2$, $e_1$
  and $e_2$  versus time for time interval $[0,10]$ and $x_{10} =
  0.105$; (e)~analogous phase-plane plot for interval $[0,20]$;
  (f)~phase-plane plot for interval $[0,50]$. ($x_{20} = e_{10} =
  e_{20} = 0$.) }
  The above plots are not too representative because they are close to the
  stability boundary. We can conclude that the boundary of the domain of
  attraction goes somewhere between points $(0.103, 0)$ and
  $(0.105, 0)$.

 Let's check for another equilibria of the stabilized
  system:
\begin{mathematica}
{x1[t], x2[t], e1[t], e2[t]} /. Solve[
    {{0, 0} == (A - BK) . x + BK . e + f,
    {0, 0} == (A - Fc) . e + f - F g},
    {x1[t], x2[t], e1[t], e2[t]}
] // N
.
Out[20]= {{0, 0, 0, 0}, {-1. I, -1. I, -1. I, -1. I},

>    {1. I, 1. I, 1. I, 1. I}, {-0.173925 I, -0.316228 I, -0.013835 I,

>     -0.0316228 I}, {0.173925 I, 0.316228 I, 0.013835 I, 0.0316228 I}}
\end{mathematica}
  As we can see there are four more equlibria but they are all complex
  so they rather don't influence the behavior of plotted trajectories.

\end{document}