Cascading Plants, Circuits, Stability
=====================================

We connect filters with differential equations and then investigate
how filters can be implemented with circuits, 
using differential equations.
A criterion on stability is applied to a traffic problem.

Cascading Plants
----------------

We defined a plant as a linear, time invariant, causal maps,
which maps signals into signals.  
A filter is one obvious example of a plant.
Other examples are mechanisms and circuits.

Let :math:`f` and :math:`g` be functions of :math:`t`.
The  convolution :math:`f \star g` of :math:`f` and :math:`g` is

.. math::

   f \star g = \int_0^t f(\tau) g(t - \tau) d \tau.

By the linearity, time invariance, and causality, 
the output is :math:`y = h \star u`
where :math:`h` is the response to the unit impulse 
and :math:`u` is the input.

The Laplace transform leads to an operational calculus.
The differential equation

.. math::

  \frac{d}{d~\!t} y(t) + y(t) = u(t), \quad y(0) = 0

is transformed into

.. math::

   s Y(s) + Y(s) = U(s) \quad \Rightarrow \quad
   Y(s) = \left( \frac{1}{s+1} \right) U(s),

where :math:`Y(s) = {\cal L}\{ y \}` and
:math:`U(s) = {\cal L}\{ u \}` are the Laplace transforms
for :math:`y` and :math:`u`.

In the time domain, the solution is then

.. math::

   y(t) = e^{-t} \star u(t) = \int_0^t e^{-\tau} u(t - \tau) d \tau.


As a computational experiment, we take as input two
signals, one at low, another at high frequency:

* :math:`u_{\rm low}(t) = \sin(2 \pi t/1000)`
  takes 1000 seconds to complete one cycle.

* :math:`u_{\rm high}(t) = \sin(2 \pi t 1000)`
  does 1000 cycles in one second.

Consider then the solution :math:`y(t)` of

.. math::

   \frac{d}{d~\!t} y(t) + y(t) = u(t), \quad y(0) = 0,

for :math:`u = u_{\rm low}` and :math:`u = u_{\rm high}`,
respectively shown in :numref:`figLowFreqSignal`
and in :numref:`figHighFreqSignal`.

.. _figLowFreqSignal:

.. figure:: ./figLowFreqSignal.png
   :align: center

   A low frequency signal :math:`\sin(2 \pi t/1000)`.

.. _figHighFreqSignal:

.. figure:: ./figHighFreqSignal.png

   A high frequency signal :math:`\sin(2 \pi t 1000)`.

Using ``SymPy.jl``, we start defining the
low frequency and the high frequency signals:

::

    lowfreq(t) = sin(2*pi*t/1000)
    highfreq(t) = sin(2*pi*t*1000)

and then define the left of the ODE and the initial condition:

::

    t = Sym("t")
    y = SymFunction("y")
    leftode = diff(y(t), t) + y(t)
    initcond = Dict(y(0) => 0)

Then the right hand side for the low frequency signal:

::

    filtlow = Eq(leftode, lowfreq(t))
    outlow = dsolve(filtlow, ics=initcond)

The output for :math:`u(t) = u_{\rm low}(t) = \sin(2 \pi t/1000)` is

.. math::

   \begin{array}{rcl}
   y(t) & = & 0.99996052314088 \sin(0.00628318530717959 t) \\
        &   & -~~ 0.00628293726675839 \cos(0.00628318530717959 t) \\
        &   & +~~ 0.00628293726675839 e^{-t}.
   \end{array}

For the high frequency signal we do:

::

    filthigh = Eq(leftode, highfreq(t))
    outhigh = dsolve(filthigh, ics=initcond)

The output for :math:`u(t) = u_{\rm high}(t) = \sin(2 \pi t 1000)` is

.. math::

   \begin{array}{rcl}
   y(t) & = & 2.53302952689605\mbox{E}\!\!-\!\!8 \sin(6283.18530717959 t) \\
        &   & -~~ 0.000159154939060454 \cos(6283.18530717959 t) \\
        &   & +~~ 0.000159154939060454 e^{-t}.
   \end{array}

The output for both signals is shown 
in :numref:`figLowFreqSignalOutput` and
in :numref:`figHighFreqSignalOutput`.
We see that the low frequency signal passes, whereas
the output of the high frequency signal is suppressed.
The differential equation represents a :index:`lowpass filter`.

.. _figLowFreqSignalOutput:

.. figure:: ./figLowFreqSignalOutput.png
   :align: center

   The output of the low frequency signal :math:`\sin(2 \pi t/1000)`.

.. _figHighFreqSignalOutput:

.. figure:: ./figHighFreqSignalOutput.png
   :align: center

   The output of the high frequency signal :math:`\sin(2 \pi t 1000)`.

The first exercise considers a highpass filter.

A :index:`bandpass filter` can be defined as a lowpass filter
followed by highpass filter.
Consider :math:`u(t)` as input to

.. math::

   \frac{d}{d~\!t} y(t) + y(t) = u(t), \quad y(0) = 0,

and then give :math:`y(t)` as input to

.. math::

   \frac{d}{d~\!t} z(t) + z(t) = \frac{d}{d~\!t} y(t), \quad z(0) = 0,

which gives a bandpass filter.
The frequency domain representation is

.. math::

   Z(s) = \left( \frac{s}{s+1} \right)  \left( \frac{1}{s+1} \right) U(s),
   \mbox{ as }
   Y(s) = \left( \frac{1}{s+1} \right) U(s),

with :math:`Y(s) = {\cal L}\{ y \}`, :math:`U(s) = {\cal L}\{ u \}`,
and :math:`Z(s) = {\cal L}\{ z \}`.
This frequency domain representation is shown
in :numref:`figBandPassFilter`.

.. _figBandPassFilter:

.. figure:: ./figBandPassFilter.png
   :align: center

   Frequency domain representation of a bandpass filter.

This is an example of two *cascading plants*.
We have that
*The transfer function of plants in cascade is the product
of the transfer functions of the components in the cascade.*

Electric Circuits
-----------------

We consider the realization of lowpass and highpass filters
as electric circuits.

The relation of voltage to current is shown in :numref:`figVoltageCurrent`.

.. _figVoltageCurrent:

.. figure:: ./figVoltageCurrent.png
   :align: center

   The relation between the voltage :math:`v(t)` across, and
   the current :math:`i(t)` through a resistor, capacitor, and inductor.

.. index:: voltage, current, resistor, capacitor, inductor, surge impedance

.. topic:: Definition of surge impedance of a device

   The ratio :math:`Z = V/I` is the *surge impedance* of a device.

Circuits are combined into networks, either by arranging them
in series or in parallel, as shown in :numref:`figSeriesParallelCircuits`.
If :math:`Z_1` and :math:`Z_2` are the transfer functions of two circuits,
then the network defined by arranging the circuits in series has
transfer function

.. math::

   Z = Z_1 + Z_2.

If the circuits are arranged in parallel, then the transfer function
of the network is

.. math::

   Z = \frac{1}{1/Z_1 + 1/Z_2}.

.. _figSeriesParallelCircuits:

.. figure:: ./figSeriesParallelCircuits.png
   :align: center

   Series and parallel circuits, with their transfer functions.

A circuit realization of a :index:`lowpass filter`
is shown in :numref:`figCircuitLowPassFilter`.
For the resistor :math:`Z_1 = R`, and 
the capacitor :math:`Z_2 = \frac{1}{Cs}`,
the transfer function of 
the realization in :numref:`figCircuitLowPassFilter` is

.. math::

   \frac{Y}{U} = \frac{\frac{1}{Cs}}{R + \frac{1}{Cs}} = \frac{1}{RC s + 1}.

.. _figCircuitLowPassFilter:

.. figure:: ./figCircuitLowPassFilter.png
   :align: center

   A realization of a lowpass filter.

In :numref:`figCircuitHighPassFilter` is a circuit realization
of a :index:`highpass filter`.
For the capacitor :math:`Z_1 = \frac{1}{Cs}`, 
and the resistor :math:`Z_2 = R`,
the transfer function of 
the realization in :numref:`figCircuitHighPassFilter` is

.. math::

   \frac{Y}{U} = \frac{R}{R + \frac{1}{Cs}} = \frac{R C s}{RC s + 1}.

.. _figCircuitHighPassFilter:

.. figure:: ./figCircuitHighPassFilter.png
   :align: center

   A realization of a highpass filter.

Stability
---------

In this section we develop a criterion to decide
whether a plant is stable or not.

.. topic:: Definition of stability of a plant

   A plant is *stable* if bounded inputs give bounded outputs.

The inpulse response of a stable plant is decaying.

.. index:: pole

.. topic:: Definition of the poles of a rational function

   Consider the quotient :math:`q = f(x)/g(x)`
   of two polynomials :math:`f` and :math:`g`.
   The *poles* of :math:`q` are the roots
   of the denominator :math:`g(x) = 0`.

A criterion on the :index:`stability of a plant` is provided
by the theorem below.

.. topic:: Theorem on the stability criterion of a plant

   Let the plant :math:`P` have transfer
   function :math:`P(s) = f(s)/g(s)`.
   :math:`P` is stable if and only if all poles 
   lie in the open left complex plane:
   for any pole :math:`z`, :math:`g(z) = 0`:
   :math:`\mbox{Re}(z) < 0`.

For example, consider

.. math::

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

The poles are :math:`-7/2 \pm \sqrt{33}/2`,
so :math:`P` is stable.

.. index:: two vehicle traffic

Two vehicle traffic is a situation where stability matters.
Let :math:`u(t)` be the position of the leading vehicle,
and :math:`y(t)` be the position of the trailing vehicle.

Modeling with an ordinary differential equation,
with delay :math:`\tau` and parameter :math:`\alpha`:

.. math::

   \frac{d^2}{d~\!t^2} y(t + \tau)
   = \alpha \left( \frac{d}{d~\!t} u(t) - \frac{d}{d~\!t} y(t) \right),

subject to the initial 
conditions :math:`y(0) = y'(0) = u(0) = u'(0) = 0`.

Integration leads to the differential equation:

.. math::

   \frac{d}{d~\!t} y(t+\tau) + \alpha y(t) = \alpha u(t).

In the frequency domain, we have:

.. math::

   Y(s) = \frac{\alpha e^{-\tau s}}{ s + \alpha e^{-\tau s}} U(s).

This extends to the stability of following cars.

Proposal for a Topic of a Project
---------------------------------

Consider the model for two vehicle traffic

.. math::

   \frac{d}{d~\!t} y(t+\tau) + \alpha y(t) = \alpha u(t).

* Obtain the stability condition for this problem.

  In particular, what is the relation
  between :math:`\tau` and :math:`\alpha`
  for this problem to be stable?

* Show that even when stability is satisfied,
  a long platoon of vehicles is inherently unsafe because changes 
  in speed of the lead car produce larger and larger deviations
  for cars farther down the platoon.

Exercises
---------

1. The differential equation

   .. math::

      \frac{d}{d~\!t} y(t) + y(t) = \frac{d}{d~\!t} u(t), \quad y(0) = 0

   is transformed into

   .. math::

      s Y(s) + Y(s) = s U(s) \quad \Rightarrow \quad
      Y(s) = \left( \frac{s}{s+1} \right) U(s),
 
   where :math:`Y(s) = {\cal L}\{ y \}` and :math:`U(s) = {\cal L}\{ u \}`.

   * Solve the above ODE
     for :math:`u = u_{\rm low}` and :math:`u = u_{\rm high}`,
     that is: for a low and a high frequency signal, as defined earlier.

   * Verify that the ODE represents a :index:`highpass filter`.

2. Consider a plant composed of

   * an inverting linear amplifier of gain :math:`-5`,

   * followed by an integrator :math:`\int_0^t`,

   * followed by a convolution with :math:`e^{-2t}`,

   * followed by a 3-second delay,

   shown in :numref:`figCascadingPlants`.

   Compute the time domain version of this cascading plant.

.. _figCascadingPlants:

.. figure:: ./figCascadingPlants.png
   :align: center

   An example of cascading plants.

3. Consider the network in :numref:`figCircuitsNetwork`.

   Prove that the network has transfer function
   :math:`\displaystyle \frac{Y}{U} = \frac{Z_2}{Z_1 + Z_2}`.

.. _figCircuitsNetwork:

.. figure:: ./figCircuitsNetwork.png
   :align: center

   A network with input :math:`U` and output :math:`Y`.

4. Let :math:`g(s) = a s^2 + b s + c`, with :math:`a, b, c > 0`.

   What is the sign of the real parts of the roots of :math:`g`?

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

1. Charles R. MacCluer, middle 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.