import numpy
import matplotlib.pyplot as plt
%matplotlib inline

import tbcontrol
tbcontrol.expectversion('0.1.3')

73. Numeric simulation

Let’s start with a very simple numeric simulation of a proportional controller acting on a first order process \(G = \frac{y}{u} = \frac{K}{\tau s + 1}\).

\begin{align} y(s)(\tau s + 1) &= K u(s) \\ \tau s y(s) + y(s) &= K u(s) \\ s y(s) &= \frac{1}{\tau}\left(K u(s) - y(s)\right)\\ \frac{dy}{dt} &= \frac{1}{\tau}\left(Ku - y\right) \end{align}

K = 3
tau = 2

Kc = 2

ts = numpy.linspace(0, 5, 1000)
dt = ts[1]

y_continuous = []
u_continuous = []
y = 0
sp = 1
for t in ts:
    e = sp - y
    u = Kc*e
    dydt = 1/tau*(K*u - y)

    y += dydt*dt

    u_continuous.append(u)
    y_continuous.append(y)

def plot_continuous():
    fig, [ax_u, ax_y] = plt.subplots(2, 1, sharex=True)
    ax_u.plot(ts, u_continuous)
    ax_u.set_ylabel('$u(t)$')
    ax_y.plot(ts, y_continuous)
    ax_y.axhline(1)
    ax_y.set_ylabel('$y(t)$')
    ax_y.set_ylabel('time')

    return ax_u, ax_y

plot_continuous()

(<matplotlib.axes._subplots.AxesSubplot at 0x11dfc79e8>,
 <matplotlib.axes._subplots.AxesSubplot at 0x11e022240>)

Now, let’s use a discrete version of the same controller. We will assume a Zero Order Hold between the controller and the system.

The discrete controller will only run at the sampling points. Now, we have an integration timestep and a discrete timestep. We call the integration timestep dt and the sampling time \(\Delta t\). We may think that if we use the above for loop to update t, it will eventually be equal to \(\Delta t\), but this is not true in general. Instead, we set a target for the sampling time and check if we are at a time greater than that time, then set a new time one sampling time in the future.

DeltaT = 0.5 # sampling time

u_discrete = []
y_discrete = []
y = 0
sp = 1
next_sample = 0
for t in ts:
    if t >= next_sample:
        e = sp - y
        u = Kc*e
        next_sample += DeltaT
    dydt = 1/tau*(K*u - y)
    y += dydt*dt

    u_discrete.append(u)
    y_discrete.append(y)

ax_u, ax_y = plot_continuous()
ax_u.plot(ts, u_discrete)
ax_y.plot(ts, y_discrete)

[<matplotlib.lines.Line2D at 0x11e1877f0>]

Notice the difference? Because the discrete controller only calculates its values at the sampling points and because the ZOH keeps its output constant, the discrete controller takes more action later on, in fact introducing some oscillation where the continuous controller could use arbitrarily large gain.

74. Symbolic calculation

Now we will try to replicate that last figure without doing numeric simulation. The continuous controller is trivially done via the Laplace transform:

import sympy
sympy.init_printing()

s = sympy.Symbol('s')
t = sympy.Symbol('t', positive=True)

Gc = Kc  # controller
G = K/(tau*s + 1)  # system

G_cl = Gc*G/(1 + Gc*G)

rs = 1/s  # step input r(s)

ys = rs*G_cl  # system output y(s)
es = rs - ys  # error
us = Gc*es  # controller output

yt = sympy.inverse_laplace_transform(ys, s, t)
ut = sympy.inverse_laplace_transform(us, s, t)

sympy.plot(ut, (t, 0, 5))
sympy.plot(yt, (t, 0, 5))

<sympy.plotting.plot.Plot at 0x11f90c4a8>

Now for the discrete controller. First we need some new symbols.

z, q = sympy.symbols('z, q')

We get the z transform of a sampled step from the table in the datasheet.

rz = 1/(1 - z**-1)

If we rewrite a z-transformed signal as a polynomial \(r(z) = a_0 + a_1 z^{-1} + a_2 z^{-2} \dots\) , the coefficients give the values at the sampling points, so \(a_0 = r(0)\), \(a_1 = r(\Delta t)\), \(a_2 = r(2\Delta t)\) and so on. We can obtain these coefficients easily using a Taylor expansion in sympy.

rz.subs(z, q**-1).series()

$$1 + q + q^{2} + q^{3} + q^{4} + q^{5} + O\left(q^{6}\right)$$

We can see clearly that all the coefficients are 1 for the step.

There is more detail in this notebook if you want to refresh your memory

The tbcontrol.symbolic.sampledvalues function automates this process:

from tbcontrol.symbolic import sampledvalues

def plotdiscrete(fz, N):
    values = sampledvalues(fz, z, N)
    times = [n*DeltaT for n in range(N)]
    plt.stem(times, values)

plotdiscrete(rz, 10)

Let’s move on to the other transfer functions. The controller is simple:

Gcz = Kc

The controller is connected to a hold element (\(H\)) which is connected to the system itself (\(G\)). The z-transform of this combination can be written \(\mathcal{Z}\{H(s)G(s)\}\). Remember, \(H(s) = \frac{1}{s}(1 - e^{-\Delta t s})\). Now we can show

\begin{align} \mathcal{Z}\left\{{H(s)G(s)}\right\} &= \mathcal{Z}\left\{\frac{1}{s}(1 - e^{-Ts})G(s)\right\} \\ &= \mathcal{Z}\left\{\underbrace{\frac{G(s)}{s}}_{F(s)}(1 - e^{-Ts})\right\} \\ &= \mathcal{Z}\left\{F(s) - F(s)e^{-Ts}\right\} \\ &= \mathcal{Z}\left\{F(s)\right\} - \mathcal{Z}\left\{F(s)e^{-Ts}\right\} \\ &= F(z) - F(z)z^{-1} \\ &= F(z)(1 - z^{-1}) \end{align}

So the z transform we’re looking for will be \(F(z)(1 - z^{-1})\) with \(F(z)\) being the transform on the right of the table of \(\frac{1}{s}G(s)\).

For \(G(s) = \frac{K}{\tau s + 1}\), \(F(s) = \frac{K}{s(\tau s + 1)}\), \(F(z) = \frac{K(1 - b)z^{-1}} {(1 - z^{-1})(1 - bz^{-1})}\)

a = 1/tau
b = sympy.exp(-a*DeltaT)
Fz = K*(1 - b)*z**-1/((1 - z**-1)*(1 - b*z**-1))
HGz = Fz - z**-1*Fz

Let’s verify that this is correct by plotting the continuous response along with the discrete values:

plotdiscrete(rz*HGz, 10)
plt.plot(ts, K*(1 - numpy.exp(-ts/tau)))

[<matplotlib.lines.Line2D at 0x11df60c88>]

Now we have the discrete transfer functions, we can repeat the same calculation as before.

yz = rz*Gcz*HGz/(1 + Gcz*HGz)

plt.plot(ts, y_discrete)
plotdiscrete(yz, 10)
plt.legend(['Numeric simulation', 'Analytical at sampling points'])

<matplotlib.legend.Legend at 0x11fefd9e8>

So now we have recovered the response we calculated numerically before analytically. Let’s see if we can reproduce the numeric values between the sampling points for the continuous system output. Our first step is to construct the Laplace transform of the controller output. Let’s look again what that looked like:

ez = rz - yz
uz = Gcz*ez
plotdiscrete(uz, 10)
plt.plot(ts, u_discrete)

[<matplotlib.lines.Line2D at 0x11ff33f28>]

We can interpret this as shifted pulse signals added together.

If we were trying to calculate the response of the system to a single pulse input, it would be simple.

Since

\begin{align} y(s) &= G(s)u(s) \\ y(t) &= \mathcal{L}^{-1}\{G(s)u(s)\} \end{align}

and a pulse signal of width \(\Delta t\) and height \(v\) has Laplace transform \(\frac{v}{s}\left(1 - e^{-\Delta t s}\right) = vH(s)\), the single output would be easy to obtain:

Hs = 1/s*(1 - sympy.exp(-DeltaT*s))

u_single_pulse = 2*Hs
y_single_pulse = sympy.inverse_laplace_transform(G*u_single_pulse, s, t)

sympy.plot(y_single_pulse, (t, 0, 5))

<sympy.plotting.plot.Plot at 0x1201e82b0>

For subsequent pulses, we simply shift the pulse one time step up. So, given

\[u(z) = a_0 + a_1 z^{-1} + a_2 z^{-2} + \cdots = \sum_{i=0}^{\infty} a_i z^{-i}\]

If the output of the controller is zero-order held so that the held version of the signal is \(u_h\), we can write

\[u_h(s) = a_0H(s) + a_1 H(s) e^{-\Delta t s} + a_2 H(s) e^{-2 \Delta t s} = \sum_{i=0}^{\infty} a_i H(s) e^{-i\Delta t s}\]

This maps quite elegantly to the following generator expression:

uhs = sum(ai*Hs*sympy.exp(-i*DeltaT*s)
          for i, ai in enumerate(sampledvalues(uz, z, 10)))

Equivalently, we can write

uhs = 0
a = sampledvalues(uz, z, 10)
for i in range(10):
    uhs += a[i]*Hs*sympy.exp(-i*DeltaT*s)

Now we can construct the continuous response like this:

ys = uhs*G

yt = sympy.inverse_laplace_transform(ys, s, t)

plt.plot(ts, y_discrete)
plt.plot(ts, tbcontrol.symbolic.evaluate_at_times(yt, t, ts), '--')

[<matplotlib.lines.Line2D at 0x11f9ee160>]

Notice that the analytical solution and numeric solution agree, but make sure you understand what the difference is in approach.