Control, Root Locus, Nyquist Analysis
=====================================

In this lecture, feedback control is introduced.

Control
-------
 
We consider a plant :math:`P`, with controller :math:`H`,
as shown in :numref:`figControlFeedback`.

.. _figControlFeedback:

.. figure:: ./figControlFeedback.png
   :align: center

   The top picture represents an open loop system.
   The bottom piciture shows a closed loop system,
   obtained with feeding back the output to the input.

Feeding back the output to the input, an :index:`open loop system`
becomes a :index:`closed loop system`.
An example of :index:`feedback` is steering control for an automobile.

In designing a servomechanism, we work with :math:`P(s) = 1/s`
as the transfer function of the plant.
The controller has a transfer function :math:`H(s) = \kappa + \alpha/s`
and obtain a *Proportional Integral* (PI) controller.
The problem is then to find the gains :math:`\kappa` and :math:`\alpha`
that yield a step response with

* at most a :math:`20\%` overshoot, and

* a 1-second rise time, the time to reach one.

For general plant with transfer function :math:`P`
and controller :math:`H`, 
the open loop function has transfer function :math:`H(s) P(s)`.
The output is subtracted from the input in the closed loop:

.. math::

   Y(s) = H(s) P(s) ( U(s) - Y(s) ).

The transfer function of the closed loop system is then

.. math::

   (1 + H(s) P(s)) Y(s) = H(s) P(s) U(s)
   \quad \Rightarrow \quad
   \frac{Y(s)}{U(s)}
   = \frac{H(s)P(s)}{1 + H(s)P(s)}.

For the model of the servomechanism,
:math:`P(s) = 1/s`.  and :math:`H(s) = \kappa + \alpha/s`.
The transfer function of the servomechanism is

.. math::

   G(s) = \frac{Y(s)}{U(s)}
   = \frac{H(s)P(s)}{1 + H(s)P(s)}
   = \frac{\kappa s + \alpha}{s^2 + \kappa s + \alpha}.

The response to the unit step is :math:`U(s) = 1/s`.

.. math::

   Y(s) = G(s) U(s) = 
   \left( \frac{\kappa s + \alpha}{s^2 + \kappa s + \alpha} \right)
   \frac{1}{s}
   = \frac{1}{s} - \frac{s}{s^2 + \kappa s + \alpha}

We consider again :index: `damped harmonic motion`.
The Laplace transform of the output of the unit step is

.. math::

   Y(s) = \frac{1}{s} - \frac{s}{s^2 + \kappa s + \alpha}.

The denominator is quadratic of :math:`Y(s)`,
which corresponds to a second order differential equation.
Substituting :math:`\kappa = 2 \zeta \omega_n`
and :math:`\alpha = \omega_n^2`
introduces the standard notation for :index:`damped harmonic motion`:

.. math::

   Y(s) = \frac{1}{s} - \frac{s}{s^2 + 2 \zeta \omega_n s + \omega_n^2}.

which is equivalent to

.. math::

   Y(s) = \frac{1}{s} - \frac{s + \zeta \omega_n}{D}
   + \left( \frac{\zeta}{\sqrt{1 - \zeta^2}} \right)
   \frac{\omega_d}{D},

with

.. math::

   D = (s + \zeta \omega_n)^2 + \omega_d^2,
   \quad \omega_d = \sqrt{1 - \zeta^2} \omega_n.

Notice the distinction between :math:`\omega_n`
and :math:`\omega_d`.

To solve the ordinary differential equation,
we apply the inverse Laplace transform, as done
in the code below, executed in a Python session.

::

   >>> from sympy import var
   >>> from sympy import inverse_laplace_transform as ilap

   >>> omega = var('omega', real=True)
   >>> zeta = var('zeta', real=True)
   >>> s, t = var('s, t')

   >>> denominator = (s + zeta*omega)**2 + omega**2

   >>> Y1 = - (s + zeta*omega)/denominator
   >>> ilap(Y1, s, t)
   -exp(-omega*t*zeta)*cos(omega*t)*Heaviside(t)

   >>> Y2 = omega/denominator
   >>> ilap(Y2, s, t)
   exp(-omega*t*zeta)*sin(omega*t)*Heaviside(t)

Let us look at the output in the time domain.
The inverse Laplace transform of

.. math::

   Y(s) = \frac{1}{s} - \frac{s + \zeta \omega_n}{D}
        + \left( \frac{\zeta}{\sqrt{1 - \zeta^2}} \right)
          \frac{\omega_d}{D},

with

.. math::

   D = (s + \zeta \omega_n)^2 + \omega_d^2,
   \quad \omega_d = \sqrt{1 - \zeta^2},

is then

.. math::

   \begin{array}{lcr}
   y(t) & = & {\displaystyle 1 - e^{-\zeta \omega_n t} \cos(\omega_d t)
   + \left( \frac{\zeta}{\sqrt{1 - \zeta^2}} \right)
   e^{-\zeta \omega_n t} \sin(\omega_d t)} \\
   & = & {\displaystyle 1 - e^{-\zeta \omega_n t}
   \left[
   \cos(\omega_d t) 
   - \left( \frac{\zeta}{\sqrt{1 - \zeta^2}} \right)
   \sin(\omega_d t)
   \right].}
   \end{array}

To determine the :index:`overshoot`, we run simulations.
The :math:`\zeta` in the above solution for :math:`y(t)`
determines the overshoot, which should be less than :math:`20\%`.
To study the effect of $\zeta$, we go to normalized time

.. math::

   \tau = \omega_n t, \quad \omega_d = \sqrt{1 - \zeta^2},

and graph :math:`y(\tau)`, for various values of :math:`\zeta`,
shown in :numref:`figServoZeta01`, 
:numref:`figServoZetas`, and :numref:`figServoZeta78`.

.. _figServoZeta01:

.. figure:: ./figServoZeta01.png
   :align: center

   A simulation for :math:`\zeta = 0.1`.

.. _figServoZetas:

.. figure:: ./figServoZetas.png
   :align: center

   Simulations for :math:`\zeta = 0.1, 0.2, 0.4, 0.6, 0.7, 0.8`.

.. _figServoZeta78:

.. figure:: ./figServoZeta78.png
   :align: center

   Simulations for :math:`\zeta = 0.7, 0.8`.

Summarizing the observations from the simulations in
:numref:`figServoZeta01`, :numref:`figServoZetas`, and
:numref:`figServoZeta78`, gives a solution to the design problem.

* With :math:`\zeta = 0.7`, we go slightly above 1.2,
  so :math:`\zeta = 0.8` is a safer value to meet the less
  than :math:`20\%` overshoot.

* With :math:`\omega_n`, we can set the rise time,
  as :math:`\tau = \omega_n t`.
  For :math:`\zeta = 0.8`, at :math:`\tau = 1`,
  we reach 0.97, and 1.01 at :math:`\tau = 1.1`.
  If we want one at :math:`t=1`,
  then choose :math:`\omega = \tau/t = 1.1`.

* Use :math:`\kappa = 2 \zeta \omega_n = 1.76`
  and :math:`\alpha = \omega_n^2 = 1.21`
  in the design of the controller.

Root Locus
----------

Consider a plant with transfer function

.. math::

   P(s) = \frac{s^2 + 7 s + 4}{s^3 - 3 s^2 + 2 s + 5}.

With proportional feedback :math:`H(s) = \kappa`,
the closed loop system has transfer function

.. math::

   Q(s) = \frac{P}{1 + \kappa P}.

The denominator of :math:`Q(s)` is

.. math::

   \kappa \left( s^2 + 7 s + 4 \right) + s^3 - 3 s^2 + 2 s + 5.

The closed loop system is stable if and only if all roots
of the denominator all lie in the open left 
plane :math:`\mbox{real}(s) < 0`.

We apply the :index:`root locus method`, using SymPy,
in a Julia function defined below:

::

   """
       function locus(kp::Float64, p)

   returns the roots of the polynomial p in s defined by
   the coefficients where kappa is replaced by kp.
   """
   function locus(kp::Float64, p)
       cffs = [p.coeff(s,k) for k=0:3]
       c1 = [c.subs(Dict(kappa=>kp)) for c in cffs]
       cf1 = [Float64(c1[k]) for k=1:size(c1,1)]
       p1s = sum([cf1[k]*s^(k-1) for k=1:size(cf1,1)])
       r1 = roots(p1s)
       v1 = [Float64(real(r[1]))+Float64(imag(r[1]))*Complex(0,1)
             for r in r1]
       tripletsort!(v1)
       return v1
   end

The function ``tripletsort!(v1)`` sorts the numbers in ``v1``.

The output of runs with this function is shown in
:numref:`figRootLocusKappa2`,
:numref:`figRootLocusKappa3`, and
:numref:`figRootLocusKappa5`.

.. _figRootLocusKappa2:

.. figure:: ./figRootLocusKappa2.png
   :align: center

   The poles for :math:`\kappa = 2`, unstable.

.. _figRootLocusKappa3:

.. figure:: ./figRootLocusKappa3.png
   :align: center

   The poles for :math:`\kappa = 3`, still unstable.

.. _figRootLocusKappa5:

.. figure:: ./figRootLocusKappa5.png
   :align: center

   The poles for :math:`\kappa = 5`, stable.

Nyquist analysis
----------------

The stability of the :index:`closed loop system`
for proportional feedback is equivalent to finding
a :math:`\kappa` so that the roots of

.. math::

   \frac{1}{\kappa} + H(s) P(s) = 0

all lie in the open left half plane.

The method of Nyquist analysis operates as follows:

1. Set :math:`G(s) = H(s) P(s)`.

2. Graph the image :math:`\Gamma` of :math:`z = G(s)`, 
   for :math:`s = \omega i`, for increasing :math:`\omega`.

All :math:`\kappa` so that :math:`-1/\kappa` lies in the unbounded component
of the *z*-plane disjointed by :math:`\Gamma` 
will DEstabilize the closed loop system.

Consider the example:

.. math::

   P(s) = 1/(s-1), \quad H(s) = 1/(s-2), \quad G(s) = H(s) P(s). 

The code to make the :index:`complex plot` is below:

::

   I = Complex(0, 1)
   r(omega) = G(omega*I)
   omegarange = -100:0.01:100
   rt = [(real(r(t)), imag(r(t)))
          for r in omegarange]
   plot(rt, xlabel="real part",
            ylabel="imag part", label="G(w*I)")

In :numref:`figNyquistPlot`, we see that
the entire negative real axis is outside the curve.
Thus, there is no :math:`\kappa > 0` that stabilizes.

.. _figNyquistPlot:

.. figure:: ./figNyquistPlot.png
   :align: center

   The complex plot of :math:`1/((s-1)(s-2))`.

In :numref:`figNyquistPlot2`, we see that, for another example,
the entire negative real axis is outside the curve,
except for :math:`[-0.25, 0]`,
so :math:`-1/\kappa \in [-0.25, 0]` or :math:`\kappa > 4`,
which agrees with the root locus method.

.. _figNyquistPlot2:

.. figure:: ./figNyquistPlot2.png
   :align: center

   The complex plot of
   :math:`G(s) = (s^2 + 7s + 4)/(s^3 - 3 s^2 + 2 s + 5)`.

We end with the argument principle.

.. topic:: Theorem on the argument principle

   * Let :math:`\Omega` be a domain in the complex plane.

   * Suppose :math:`f(s)` is meromorphic (analytic everywhere,
     except for poles) on :math:`\Omega`.

   * Let :math:`C` be the closed curve in :math:`\Omega`
     that does not pass through any zeros or poles of :math:`f(s)`.

   Then:

   .. math::

      \frac{1}{2 \pi i} \int_C \frac{f'(s)}{f(s)} ds
      = \#\mbox{zeros} - \#\mbox{poles}

   in :math:`\Omega` enclosed by :math:`C`, counted with multiplicities.

.. index:: winding number

As a corollary, *the winding number* :math:`n(z_0, C)` counts
the number of roots minus the number of poles
of :math:`f(s) - z_0 = 0` within :math:`C`.

.. index:: Nyquist contour

.. topic:: Definition of the Nyquist contour

   For :math:`R>0`, *the Nyquist contour* is 

   1. the imaginary axis from :math:`-R i` to :math:`+R i`, and

   2. the semi-circle :math:`R e^{i \theta}`,
      for :math:`\theta` from :math:`+\pi/2` to :math:`-\pi/2`.

To check the stability of a rational plant :math:`P(s) = N(s)/D(s)`,
check that no poles of :math:`P` lie within the Nyquist contour
for large :math:`R`.

Proposals for a Topic of a Project
----------------------------------

1. design an automobile cruise control system

   Our text book has the design of an automobile cruise control system.
   Study the solution to this design problem in our text book.

   Do the following:

   * Describe the solution using your own words.

   * Do the computational experiments to examine the stability.

   * Write a report on your computational experiments.

2. root finding with the winding number

   The winding number leads to a numerical method to approximate
   all roots of a general nonlinear equation in a domain of
   the complex plane.

   Study the proof of the argument principle theorem in our text book
   and make sure you understand the definition of the winding number.

   Do the following:

   * Describe the algorithm to compute all complex roots in a domain
     of a general nonlinear function.

   * Experiment with the implementation in Pari/GP.

   * Write a report on your computational experiments.

Exercises
---------

1. Apply the partial fraction decomposition to justify

   .. math::

      \left( \frac{\kappa s + \alpha}{s^2 + \kappa s + \alpha} \right)
      \frac{1}{s} = \frac{1}{s} - \frac{s}{s^2 + \kappa s + \alpha}.

2. For the root locus example,
   find the critical threshold for :math:`\kappa`,
   so the closed loop system is stable 
   when :math:`\kappa` crosses this threshold.

3. The function ``locus`` is inefficient
   and ``tripletsort!`` a hack.

   A continuation method operates in two stages:

   * Update :math:`\kappa` as :math:`\kappa + \Delta \kappa`,
     with step size :math:`\Delta \kappa`.

   * After each update, apply Newton's method to correct the root.

   Apply this continuation method on the example
   of the root locus method.

Bibliography
------------

1. Charles R. MacCluer, end of Chapter 10 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.

2. Richard C. Dorf and Robert H. Bishop:
   *Modern Control Systems.*  Nineth Edition, Prentice Hall, 2001.