Consider the Continuous Time System With Transfer Function H S S 1 S 1 S2 4
4.4 Modelling Continuous time systems in Discrete time¶
4.4.1 Converting H(s) to H(z) using \(z=e^{sT}\)¶
The question you should have at this point in the reading, is how do we create discrete time transfer functions?
It was easy to derive transfer functions for continuous time systems
We used relations between power conjugate variables, such as Newton's law, Hooke's law, or Stokes' law,
and substituted them into conservations laws, such as D'Alembert's principle, all in the frequency domain
Instead of defining discrete time physics, we will simply convert continuous time systems to discrete time systems
And you already know how to make the conversion!
\[z=e^{sT}\]
Say, for example, you derived a continuous time transfer function of
\[H(s)=\frac{10}{s+10}\]
solve the relationship between the frequencies s and z, for s
\[s=\frac{ln(z)}{T}\]
and substitue
\[H(z)=\frac{10}{\frac{ln(z)}{T}+10}\]
That's it, we converted the continuous time transfer function to a discrete time transfer function
unfortunatley, the result is non-linear and unsable
In order to use the relationship between frequencies, s and z, we need to make a linear approximation for \(z=e^{sT}\)
In our examples we will assume a sampling frequency of 20 sample/sec, or \(T=0.05\)[s]
4.4.2 Approximating \(z=e^{sT}\)¶
4.4.2.1 The Forward Euler or forward difference approximation¶
The exponential function can be approximated by choosing a finite number of terms
from the infinite Taylor series expansion, specifically the Maclaurin series for exponentials
Taylor series
\[e^{sT}=\sum_{n=0}^{\infty}\frac{(sT)^n}{n!}=1+sT+\frac{(sT)^2}{2!}+\frac{(sT)^3}{3!}+\frac{(sT)^4}{4!}+\cdots\]
A first-order (linear) approximation from the Maclaurin series takes only the first two terms
\[z=e^{sT}\approx1+sT\]
solving for s
\[s=\frac{z-1}{T}\]
substitute this approximation into our example transfer function, \(H(s)=\frac{10}{s+10}\)
\[H(z)=\frac{10}{\left(\frac{z-1}{T}\right)+10}\]
write as a polynomial over a polynomial
\[H(z)=\frac{10T}{z-1+10T}\]
substitute our assumed sampling period, \(T=0.05\)[s]
\[H(z)=\frac{0.5}{z-0.5}\]
And there it is, we have a discrete time linear transfer function (polynomial of \(z\) over a polynomial of \(z\))
that will approximate the continuous time transfer function, \(H(s)=\frac{10}{s+10}\)
4.4.2.2 The Backward Euler or backward difference approximation¶
We can produce another first order approximation from the Maclaurin series expansion by noting that
\[e^{sT}=\frac{1}{e^{-sT}}=\frac{1}{1-sT+\frac{(-sT)^2}{2!}+\frac{(-sT)^3}{3!}+\frac{(-sT)^4}{4!}+\cdots}\]
the first order approximation is
\[z=e^{sT}\approx\frac{1}{1-sT}\]
solving for s
\[s=\frac{z-1}{zT}\]
substitute this approximation into our example transfer function, \(H(s)=\frac{10}{s+10}\)
\[H(z)=\frac{10}{\left(\frac{z-1}{zT}\right)+10}\]
write as a polynomial over a polynomial
\[H(z)=\frac{10zT}{z-1+10zT}\]
substitute our assumed sampling period, \(T=0.05\)[s]
\[H(z)=\frac{0.5z}{1.5z-1}=\frac{0.33z}{z-0.67}\]
4.4.2.3 The Bilinear approximation or Tustin's method¶
The final first order approximation we can make from the Maclaurin series starts by noting
\[e^{sT}=\frac{e^{\frac{sT}{2}}}{e^{-\frac{sT}{2}}}=\frac{1+\frac{sT}{2}+\frac{(\frac{sT}{2})^2}{2!}+\frac{(\frac{sT}{2})^3}{3!}+\frac{(\frac{sT}{2})^4}{4!}+\cdots}{1-\frac{sT}{2}+\frac{(-\frac{sT}{2})^2}{2!}+\frac{(-\frac{sT}{2})^3}{3!}+\frac{(-\frac{sT}{2})^4}{4!}+\cdots}\]
the first order approximation is
\[z=e^{sT}\approx\frac{1+\frac{sT}{2}}{1-\frac{sT}{2}}=\frac{2+sT}{2-sT}\]
solving for s
\[s=\frac{2}{T}\frac{z-1}{z+1}\]
substitute this approximation into our example transfer function, \(H(s)=\frac{10}{s+10}\)
\[H(z)=\frac{10}{\left(\frac{2}{T}\frac{z-1}{z+1}\right)+10}\]
write as a polynomial over a polynomial
\[H(z)=\frac{10zT+10T}{2z-2+10zT+10T}\]
substitute our assumed sampling period, \(T=0.05\)[s]
\[H(z)=\frac{0.5z+0.5}{2z-2+0.5z+0.5}=\frac{0.5z+0.5}{2.5z-1.5}=\frac{0.2z+0.2}{z-0.6}\]
4.4.3 Creating discrete time transfer functions with Scipy¶
In summary, we learned three different ways to write a first order approximation for \(z=e^{sT}\)
The results, when applied to our example continuous time transfer function \(H(s)=\frac{10}{s+10}\) are
\[H(z)=\frac{0.5}{z-0.5}\qquad or\qquad H(z)=\frac{0.33z}{z-0.67}\qquad or\qquad H(z)=\frac{0.2z+0.2}{z-0.6}\]
so which one is right?
They're all right!
They are just different approximations of the same thing
The code cell below shows you how to create the same discrete time transfer functions using Scipy
to_discrete
you must specify a sampling period for this method
import numpy as np import matplotlib.pyplot as plt import scipy.signal as sig # enter the continuous time transfer function and create a continuous time lti system num = [ 10 ] den = [ 1 , 10 ] H_s = sig . lti ( num , den ) # you must specify a sample period dt = 0.05 # sometimes you will get a badly conditioned coefficients warning from Scipy # its just a warning not an error H_z_fe = H_s . to_discrete ( dt = dt , method = 'euler' ) print ( 'Forward Euler \n ' , H_z_fe ) H_z_be = H_s . to_discrete ( dt = dt , method = 'backward_diff' ) print ( 'Backward Euler \n ' , H_z_be ) H_z_b = H_s . to_discrete ( dt = dt , method = 'bilinear' ) print ( 'Bilinear \n ' , H_z_b )
Forward Euler TransferFunctionDiscrete( array([0.5]), array([ 1. , -0.5]), dt: 0.05 ) Backward Euler TransferFunctionDiscrete( array([3.33333333e-01, 5.55111512e-17]), array([ 1. , -0.66666667]), dt: 0.05 ) Bilinear TransferFunctionDiscrete( array([0.2, 0.2]), array([ 1. , -0.6]), dt: 0.05 ) /opt/tljh/user/lib/python3.6/site-packages/scipy/signal/filter_design.py:1622: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless "results may be meaningless", BadCoefficients)
4.4.4 Finding a discrete time output signal by hand¶
\[H(z)\equiv\frac{y(z)}{x(z)}\Rightarrow y(z)=H(z)x(z)\]
To find an output signal using the discrete time transfer function,
we must find the z-transform of an input signal
For our example input, let's choose a step function with a magnitude of 5, rather than the unit step
The Scipy method, step(), only finds the solution for a unit step
\[x[k]=[5, 5, 5, \cdots]\]
this transforms to
\[x(z)=\frac{5z}{z-1}\]
to find an output, we need a transfer function
we will choose the bilinear transfer function we found in the previous example
\[H(z)=\frac{0.2z+0.2}{z-0.6}\]
from the transfer function
\[y(z)=H(z)x(z)=\left(\frac{0.2z+0.2}{z-0.6}\right)\left(\frac{5z}{z-1}\right)\]
\[y(z)=\frac{z^2+z}{(z-0.6)(z-1)}=\frac{z^2+z}{z^2-1.6z+0.6}\]
to transform this signal back to the time (or indexed) domain we just look it up in the z-transform table
unfortunately, this function is not in our table
we need to perform partial fraction expansion
but before we do, we need to divide the output signal, \(y(z)\), by \(z\)
always divide by \(z\) before performing partial fraction expansion
write the output signal in its factored form and find the residues
\[\frac{y(z)}{z}=\frac{z+1}{z^2-1.6z+0.6}=\frac{R_1}{z-0.6}+\frac{R_2}{z-1}\]
as with continuous time systems, there are many ways to find the residues,
but we will use Scipy
Residues
num = [ 1 , 1 ] den = [ 1 , - 1.6 , 0.6 ] r , p , k = sig . residue ( num , den ) print ( 'residues =' , r ) print ( 'poles =' , p ) print ( 'k=' , k )
residues = [-4. 5.] poles = [0.6 1. ] k= []
After we multiply \(z\) back to the other side, we see the output signal becomes
\[y(z)=\frac{-4z}{z-0.6}+\frac{5z}{z-1}\]
Now we can transform this back to the time (or indexed) domain simply by looking up the terms in the z-transform table
\[y[k]=-4(0.6)^k+5u[k]\]
If we want to convert this to an exponential, we will have to consider the sampling rate
use
\[e^{-\sigma T}=0.6\Rightarrow\sigma=10.2\]
such that
\[y[kT]=5-4e^{-10.2kT}\quad or\quad y(t)=5-4e^{-10.2t}\quad for\;t\geq0\]
We know the continuous time system we modelled had a natural frequency of 10[Np/s],
but through the bilinear approximation, we see a natural frequency of 10.2[Np/s]
Also notice, because of the discrete time bilinear approximation, the system starts at a value of 1.0, rather than 0.0
if we used a smaller sample period, the approximation would get closer to the continuous time solution
4.4.5 Finding an output signal with Scipy¶
The code cell below creates a dlti system for the output signal,
and plots it for the first 10 time samples using the impulse method
it then compares the result to the continuous time solution
# create the solution for the bilinear approximation to a magnitude 5 step t1 = np . arange ( 0 , 0.5 , 0.05 ) # space the points at the sampling period # y(z) = (z^2+z)/(z^2-1.6z+0.6) num = [ 1 , 1 , 0 ] # don't forget the zero den = [ 1 , - 1.6 , 0.6 ] y_z = sig . dlti ( num , den , dt = 0.05 ) t1 , y_kT = y_z . impulse ( t = t1 ) # create the solution for the continuous solution to a magnitude 5 step t2 = np . arange ( 0 , 0.5 , 0.001 ) # use more points for the continuous solution # y(s) = 50/(s^2+10s) num = [ 50 ] den = [ 1 , 10 , 0 ] y_s = sig . lti ( num , den ) t2 , y_t = y_s . impulse ( T = t2 ) plt . plot ( t2 , y_t ) plt . scatter ( t1 , y_kT ) plt . xlabel ( 'time, t [s]' ) plt . ylabel ( 'output' ) plt . ylim ([ - 1 , 6 ]) plt . legend ([ 'continuous' , 'bilinear' ], loc = 'lower right' ) plt . title ( 'A magnitude of 5 step response \n bilinear approximation versus continuous solution' ) plt . grid ()
Source: http://lcs-vc-marcy.syr.edu:8080/Chapter44.html