76. Dahlin controller

We will replicate the controller output in figure 17.11a. This notebook has a video commentary.

Note I am replicating these results using analytic methods to show that the artefacts are not numerical but rather fundamental to the calculations. If you simply want to simulate the action of a discrete controller on a continuous system, have a look at the Simple discrete controller simulation notebook.

import sympy
sympy.init_printing()
import tbcontrol
tbcontrol.expectversion('0.1.3')

77. Discretise the system

We need to find the corresponding z transform of the hold element and the system. Since \(H=1/s(1 - e^{-Ts})\), we can find \(F=G/s\), from there \(f(t)\) and then work out the \(z\) transform

s, t = sympy.symbols('s, t')
Gs = 1/(15*s**2 + 8*s + 1)
Gs

$$\frac{1}{15 s^{2} + 8 s + 1}$$

f = sympy.inverse_laplace_transform(Gs/s, s, t).simplify()
sympy.nsimplify(sympy.N(f)).simplify()

$$\theta\left(t\right) + \frac{3 e^{- \frac{t}{3}} \theta\left(t\right)}{2} - \frac{5 e^{- \frac{t}{5}} \theta\left(t\right)}{2}$$

We can see that \(f(t)\) is the sum of 1 and two exponentials. It is easy to determine the corresponding z transforms from the table

Time domain	Laplace-transform	z-transform (\(b=e^{-aT}\))
\(e^{-at}\)	\(\frac{1}{s+a}\)	\(\frac{1}{1-bz^{-1}}\)

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

Now the sampling interval

T = 1 # Sampling interval

def expz(a):
    b = sympy.exp(-a*T)
    return 1/(1 - b*z**-1)

Fz = -5/2*expz(1/5)  + 3/2*expz(1/3) + 1/(1 - z**-1)

Fz

$$\frac{1.5}{1 - \frac{0.716531310573789}{z}} - \frac{2.5}{1 - \frac{0.818730753077982}{z}} + \frac{1}{1 - \frac{1}{z}}$$

Let’s see if we did that right

import tbcontrol.symbolic

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

import matplotlib.pyplot as plt
%matplotlib inline

import numpy

ts = numpy.linspace(0, 10, 300)

plotdiscrete(Fz, 10)
plt.plot(ts, tbcontrol.symbolic.evaluate_at_times(f, t, ts))

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

Here is the transform of the system and the hold element. See the Discrete control notebook for the derivation.

HG_z = Fz*(1 - z**-1)

78. Dahlin Controller

The desired closed loop response is first order

θ = 0 # system dead time
λ = 1 # Dahlin's lambda
h = θ # Dahlin's h

First order response in eq 17-63

N = θ/T
A = sympy.exp(-T/λ)
yclz = (1 - A)*z**(-N-1)/(1 - A*z**-1)

K = (1/HG_z*yclz/(1 - yclz)).simplify()

K

$$\frac{1.0 \left(0.632120558828558 z^{2} - 0.970470713623842 z + 0.370831136111342\right)}{0.0279700831657292 z^{2} - 0.00455601047319787 z - 0.0234140726925314}$$

#K = 0.632/(1 - z*-1) * (1 - 1.5353*z**-1 + 0.5866*z**-2)/(0.028 + 0.0234*z**-1)

We will model the control response to a unit step in reference signal

rz = 1/(1 - z**-1) # unit step in z

Now we can calculate the z-domain version of the controller output

uz = K/(1 + K*HG_z)*rz

N = 10

plotdiscrete(uz, N)

This oscillating controller output is known as “ringing”. It is an undesirable effect.

79. Continuous response

By design the controlled variable follows an exponential at the sampling points.

yz = K*HG_z/(1 + K*HG_z)*rz

plotdiscrete(yz, N)
plt.plot(ts, 1 - numpy.exp(-ts/λ))

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

But what does it really look like between data points? First, let’s construct the response of the system to a single sampling time length pulse input

p = f - f.subs(t, t-T)

Note that there is a slight issue with finding the value of this function at zero, so we avoid that point in plotting by starting at some small value which is not zero.

smallvalue = 0.001

sympy.plot(p, (t, smallvalue, 10))

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

Let’s get the values of the controller output as a list:

u = tbcontrol.symbolic.sampledvalues(uz, z, N)

Now, we calculate the output of the system as the sum of the various pulse inputs.

yt = numpy.zeros_like(ts)
for i in range(0, N):
    yt += [float(sympy.N(u[i]*p.subs(t, ti-i + smallvalue))) for ti in ts]

Finally, we present the discrete system response, the designed response and the analytical continuous response on the same graph.

plotdiscrete(yz, 10)
plt.plot(ts, 1 - numpy.exp(-ts/λ))
plt.plot(ts, yt)

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