Laplace Transforms

  1. Calculate the Laplace Transform X(s)X(s) of the signal x(t)x(t):
x(t)=(4t3cos(5t))e2t x(t) = (4t-3cos(5t))e^{-2t}
x(t)=4te2t3cos(5t)e2tX(s)=41(s+2)23(s+2)(s+2)2+25 x(t)=4te^{-2t} - 3cos(5t)e^{-2t} \rightarrow X(s) = 4 \frac{1}{(s+2)^2} - 3 \frac{(s+2)}{(s+2)^2 + 25}

Let's verify it with Sympy

import matplotlib.pyplot as plt # import matplotlib - it will be useful later.
import sympy # import sympy
t, s = sympy.symbols('t, s')
f1 = 4*t*sympy.exp(-2*t)
f2 = - 3*sympy.cos(5*t)*sympy.exp(-2*t)
print('{} {}'.format(f1, f2)) 
4*t*exp(-2*t) -3*exp(-2*t)*cos(5*t)
F1 = sympy.laplace_transform(f1, t, s, noconds=True)
F2 = sympy.laplace_transform(f2, t, s, noconds=True)
print('{} {}'.format(F1, F2)) 
4/(s + 2)**2 -(3*s + 6)/((s + 2)**2 + 25)
  1. Calculate the Inverse Laplace Transform g(t)g(t) of the transfer function G(s)G(s):
G(s)=2+32s3+3s2+s G(s)=2 + \frac{3}{2s^3+3s^2+s}
G(s)=2+32s3+3s2+s=2+3s(2s2+3s+1)=2+32s(s+1)(s+12)=2+3s+3s+1+6s+12 G(s)=2 + \frac{3}{2s^3+3s^2+s} = 2 + \frac{3}{s(2s^2+3s+1)} = 2 + \frac{3}{2s(s+1)(s+\frac{1}{2})} = 2 + \frac{3}{s} + \frac{3}{s+1} + \frac{-6}{s+\frac{1}{2}} g(t)=2δ(t)+3+3et6e0.5t g(t) = 2\delta(t) + 3 + 3e^{-t} -6e^{-0.5t}
F1 = 2
F2 = 3/(2*s**3+3*s**2+s)
f1 = sympy.inverse_laplace_transform(F1, s, t)
f2 = sympy.inverse_laplace_transform(F2, s, t)
print('{} + {}'.format(f1, f2)) 
2*DiracDelta(t) + 3*Heaviside(t) + 3*exp(-t)*Heaviside(t) - 6*exp(-t/2)*Heaviside(t)
f3 = sympy.inverse_laplace_transform(F1+F2, s, t)
print('{}'.format(f3)) 
2*DiracDelta(t) + 3*Heaviside(t) + 3*exp(-t)*Heaviside(t) - 6*exp(-t/2)*Heaviside(t)

Block Diagrams

  1. Calculate the equivalent transfer function Gt(s)=Y(s)R1(s)G_t(s)=\frac{Y(s)}{R_1(s)} of the following block diagram:
block-diagram.png
block-diagram.png block-diagram.png block-diagram.png
G1G4(G2+G3)1+G1G4H1+H2G1G4(G2+G3) \frac{G_1G_4(G_2+G_3)}{1+G_1G_4H_1 + H_2G_1G_4(G_2+G_3)}

System Response

  1. Write the transfer function to a step input of a second order system characterised by:
  • Static gain G(0)=5G(0)=5
  • Damping ratio ξ=0.5\xi=0.5
  • Settling time ts=3st_s=3 s
  • No zeros

We would like a system with form: Ks2+2ξωns+ωn2 \frac{K}{s^2+2\xi\omega_ns+\omega_n^2}

we also know that:

ts1ξωnln(0.05)    s t_s \approx -\frac{1}{\xi\omega_n}ln(0.05)\;\; s

from which:

ωn=130.5ln(0.05)=2    rad/s \omega_n = -\frac{1}{3\cdot0.5}ln(0.05) = 2 \;\; rad/s

which means:

Ks2+2s+4 \frac{K}{s^2+2s+4}

since we want G(0)=5K=20G(0)=5 \rightarrow K=20

20s2+2s+4 \frac{20}{s^2+2s+4}
  1. Plot the qualitative behaviour of the step response of the system:
G(s)=20(3+0.1s)(s2+10s+160)(2s+10)(0.1s+5)(s2+2s+400) G(s)=\frac{20(3+0.1s)(s^2+10s+160)}{(2s+10)(0.1s+5)(s^2+2s+400)}

Let's re-write the transfer function slightly:

G(s)=200.1(30+s)(s2+10s+160)2(s+102)0.1(s+50.1)(s2+2s+400) G(s)=\frac{20\cdot0.1(30+s)(s^2+10s+160)}{2(s+\frac{10}{2})0.1(s+\frac{5}{0.1})(s^2+2s+400)}

The poles are:

s=5,50,1.+19.97498436j,1.19.97498436j s = -5, -50, -1.+19.97498436j, -1.-19.97498436j

and the zeros are:

s=30,5.+11.61895004j,5.11.61895004j s = -30, -5.+11.61895004j, -5.-11.61895004j

The system is asymptotically stable, we can use the dominant poles approximation:

G(s)=200.130160250.150(s2+2s+400)=192(s2+2s+400) G(s)=\frac{20\cdot0.1\cdot30\cdot160}{2\cdot5\cdot0.1\cdot50(s^2+2s+400)} = \frac{192}{(s^2+2s+400)}

The response to a step input can be obtained as:

Y(s)=G(s)1s=192(s2+2s+400)1s Y(s) = G(s)\frac{1}{s} = \frac{192}{(s^2+2s+400)}\frac{1}{s}

And, using partial fraction decomposition:

Y(s)=Ks+Asp1+Asp1 Y(s) = \frac{K}{s} + \frac{A}{s-p_1} + \frac{A^*}{s-p_1^*}

Where

K=192(s2+2s+400)s=0=192400=0.48 K = \frac{192}{(s^2+2s+400)} \bigg |_{s=0} = \frac{192}{400} = 0.48

And for the pole p1=1.+19.97498436𝑗p_1=−1.+19.97498436𝑗

We know that when we have complex conjugate poles - see notebook 03_Transfer_function:

f(t)=Aeσt[ej(ΦA+ωt)+ej(ΦA+ωt)]=2Aeσtcos(ΦA+ωt) f(t) = |A|e^{\sigma t} \big [ e^{j(\Phi_A+\omega t)} + e^{-j(\Phi_A+\omega t)} \big ] = 2|A|e^{\sigma t}\cos(\Phi_A+\omega t)

where:

A=(sp1)192s(sp1)(sp1)s=1+20j=192s(sp1)s=1+20j=192(1+20j)(40j)=0.24+0.012j A = (s-p_1)\frac{192}{s(s-p_1)(s-p_1^*)} \bigg |_{s=-1+20j} = \frac{192}{s(s-p_1^*)} \bigg |_{s=-1+20j} = \frac{192}{(-1+20j)( 40j)} = -0.24+0.012j
abs(-0.24+0.012)
0.22799999999999998
abs(192/(40j*(-1+20j)))
0.2397005613306827
import numpy as np
np.sqrt(399)
19.974984355438178

and the response is

y(t)=K+Aep1t+Aep1t=K+2Aeσtcos(ωt+A) y(t) = K + Ae^{p_1}t + A^*e^{p_1^*}t = K + 2|A|e^{\sigma t}cos(\omega t + \angle{A})

where p1=σ+ωtp_1=\sigma+\omega t and A=AejAA=|A|e^{j\angle{A}}

y(t)=0.48+20.23e1tcos(20tπ2) y(t) = 0.48 + 2*0.23*e^{-1t}cos(20t -\frac{\pi}{2})

We can check that the inverse laplace of the system transfer function using Sympy:

def evaluate(f, times):
    res = []
    for time in times:        
        res.append(f.evalf(subs={t:time}).n(chop=1e-5))
    return res
def L(f):
    return sympy.laplace_transform(f, t, s, noconds=True)
def invL(F):
    return sympy.inverse_laplace_transform(F, s, t)
t, s = sympy.symbols('t, s')
F = 192/(s**2+2*s+400)
invL(F)
64399etsin(399t)θ(t)133\displaystyle \frac{64 \sqrt{399} e^{- t} \sin{\left(\sqrt{399} t \right)} \theta\left(t\right)}{133}

Where:

import numpy as np
import cmath
np.sqrt(399)
19.974984355438178

And now we can plot it to verify the step response:

time = np.linspace(0, 7, 1000)
y = 0.48 + 2*0.23*np.exp(-time)
plt.plot(time, y)
plt.xlabel('time');

Finally, let's compare the result we got with what Sympy calculates

fig, ax = plt.subplots(1,1,figsize=(8,5))

time = np.linspace(0,7,1000)

ax.plot(time, evaluate(invL(F*1/s), time), linewidth=3)
ax.plot(time, y, 'g')

ax.set_title(f'y(end)={y[-1]}');
Text(0.5, 1.0, 'y(end)=0.4804111769189878')

Let's look at the response more closely:

plt.plot(time[:100], y[:100])
[<matplotlib.lines.Line2D at 0x7fd5b31b6970>]

We can also use the Python Control Library to verify our dominant pole approximation

import control
import matplotlib.pyplot as plt
s = control.TransferFunction.s
G = (20*(3+0.1*s)*(s**2+10*s+160))/((2*s+10)*(0.1*s+5)*(s**2 + 2*s + 400))
G_approx = 192/(s**2 + 2*s + 400)
T, yout = control.step_response(G, T=np.linspace(0, 7, 1000));
Tap, youtap = control.step_response(G_approx, T=np.linspace(0, 7, 1000));

plt.plot(T, yout, 'b', label='original')
plt.plot(Tap, youtap, 'r', label='approx d.p.')
plt.legend()
plt.grid()
  1. For the system G(s)G(s) defined above, calculate:
  • Its steady state value yy_\infty to a step input (when tt \rightarrow \infty)
  • Its settling time tst_s
  • If the system oscillates, calculate the period TωT_\omega of the oscillation
  • y(t)=0.48y(t\rightarrow\infty) = 0.48
  • ts1ξωnln(0.05)3t_s \approx -\frac{1}{\xi\omega_n}ln(0.05) \approx 3

given that:

  • ωn=39920\omega_n=\sqrt{399}\approx20
  • 2ξωn=220ξ=1ξ=0.052\xi\omega_n=2 \rightarrow 20 \xi =1 \rightarrow \xi = 0.05

And finally:

  • TP=2πwn1ξ2=0.314sT_{P} = \frac{2\pi}{w_n\sqrt{1-\xi^2}} = 0.314s
xi = 0.05
wn = 20
t_s = -1/(xi*wn)*np.log(0.05)
print('t_s = {}'.format(t_s))

Tp = 2*3.14/(wn*np.sqrt(1-xi**2))

print('T_p = {}'.format(Tp))
t_s = 2.995732273553991
T_ = 0.3143932374740646

Routh Criterion

  1. Determine the values of KK for which the following feedback system is asymptotically stable:
block-diagram-routh.png

The characteristic equation of the feedback system is:

1+Ks2+2s+64s2(s+0.5)=0s2(s+0.5)+K(s2+2s+64)=0s3+(0.5+K)s2+2Ks+64K=0 1+K\frac{s^2+2s+64}{s^2(s+0.5)} = 0 \rightarrow s^2(s+0.5)+K(s^2+2s+64) = 0 \rightarrow s^3 + (0.5+K)s^2 + 2Ks+64K=0

The Routh table is:

33 11 2K2K
22 0.5+K0.5+K 64K64K
11 (0.5+K)2K64K(0.5+K)2K-64K 00
00 64K64K

from which it is possible to find the following constraints on KK:

0.5+K>0,2K2+1K64K>0,K>0 0.5 + K > 0, \hspace{0.5cm} 2K^2 +1K-64K > 0, \hspace{0.5cm}K>0
K>0.5,K>31.5,K>0 K>0.5, \hspace{0.5cm} K>31.5, \hspace{0.5cm} K>0

And the system is asymptotically stable for K>31.5K>31.5

Bode Plots

  1. Plot the Bode amplitude and phase plots for the system discussed in question 7.
Gd(s)=s2+2s+64s2(s+0.5)=64(s264+2s64+1)0.5s2(s0.5+1) G_d(s) = \frac{s^2+2s+64}{s^2(s+0.5)}= \frac{64(\frac{s^2}{64}+\frac{2s}{64}+1)}{0.5s^2(\frac{s}{0.5}+1)}

KdB=20log(64/0.5)42dBK_{dB}=20\log(64/0.5) \approx 42 dB

bode-1.png

DC Motor Transfer Functions

The figures below represents a DC motor attached to an inertial load.

dc-motor-wiring dc-motor-sketch
dc-motor-pancake

The input voltage can be applied to either the field or the armature terminals. The voltages applied to the field and armature sides of the motor are represented by VfV_f and VaV_a.

The resistances and inductances of the field and armature sides of the motor are represented by RfR_f , LfL_f, RaR_a, and LaL_a. The torque generated by the motor can be assumed to be related linarly to ifi_f and iai_a, the currents in the field and armature sides of the motor, as follows:

Tm=Kifia        (1) T_m = K i_f i_a \;\;\;\;(1)

From the equation above, clearly to retain the linearity, one current must be kept constant, while the other becomes the input current.

This makes it possible to have field-current or armature-current controlled motors. We will focus on the field-current motor for our analysis.

This means that in a field-current controlled motor, the armature current is kept constant, while the field current iai_a is controlled through the field voltage VfV_f.

Tm=Kifia=Kmif        (2) T_m = K i_f i_a = K_m i_f \;\;\;\;(2)

where KmK_m is defined as the motor constant. The motor torque increases linearly with the field current.

1. Write the Transfer Function from the input current to the resulting torque of the previous equation

Tm(s)If(s)=Km \frac{T_m(s)}{I_f(s)} = K_m

For the field side of the motor the voltage/current relationship is Vf=VR+VL=Rfif+Lf(difdt)        (3) V_f =V_R +V_L = R_f i_f +L_f (di_f dt) \;\;\;\;(3)

2. Write the Transfer Function from the input voltage to the resulting current

Vf = RfI+LsI

I = Vf/Lf/(Rf/Lf+s)

If(s)Vf(s)=1Lfs+RfLf        (4) \frac{I_f(s)}{V_f(s)} = \frac{\frac{1}{L_f}}{s+\frac{R_f}{L_f}} \;\;\;\;(4)

3. Write the Transfer Function from the input voltage to the resulting motor torque

Tm(s)Vf(s)=KmLfs+RfLf        (5) \frac{T_m(s)}{V_f(s)} = \frac{\frac{K_m}{L_f}}{s+\frac{R_f}{L_f}} \;\;\;\;(5)

4. Discuss how the motor torque behaves with respect to different input signals (e.g., step input in field voltage, etc.)

This is a first-order system (type 0), and a step input in field voltage VfV_f results in an exponential rise in the motor torque.


We can now add load to the motor and verify how its behaviour changes.

The motor torque Tm(s)T_m(s) is equal to the torque delivered to the load, and this relation can be expressed as:

Tm(s)=TL(s)+Td(s)        (6) T_m(s) = T_L(s) + T_d(s) \;\;\;\;(6)

where TL(s)T_L(s), and Td(s)T_d(s) is the disturbance torque (e.g., external forces acting on the load).

We can calculate the rotational motion of the inertial load summing moments (see Figures above):

M=TLfω=Jω˙        (7) \sum{M} = T_L - f\omega = J \dot{\omega} \;\;\;\;(7)

where ff is due to friction, and JJ is the load inertia.

5. Write the Laplace transform of the load torque TL(s)T_L(s) using using equation (7) above

TL(s)=JsΩ(s)+fΩ(s) T_L(s) = Js\Omega(s)+f\Omega(s)

6. Given that the relationship between position and angular velocity ω\omega, refine the equation you have just calculated to explicit the rotor position θ\theta (hint: θ˙=ω\dot\theta = \omega)

TL(s)=Js2θ(s)+fsθ(s) T_L(s) = Js^2\theta(s) + fs\theta(s)

where Ω(s)=L(ω(t))\Omega(s) = \mathcal{L}(\omega(t))

7. Put together equations (2,4,6) and calculate the transfer function of the motor-load combination

θ(s)Vf(s)=Kms(Js+f)(Lfs+Rf)=KmJLfs(s+fJ)(s+RfLf).        (8) \frac{\theta(s)}{V_f(s)} = \frac{K_m}{s(Js+f)(L_fs+R_f)} = \frac{\frac{K_m}{JL_f}}{s(s+\frac{f}{J})(s+\frac{R_f}{L_f})}. \;\;\;\;(8)

8. Discuss the resuting system (order, dominant pole approximations, what happens when parameters change, select specific values for the parameters and draw its step response.)

θ(s)Vf(s)=KmJLfs(s+fJ)(s+RfLf)=KmfRfs(τfs+1)(τLs+1) \frac{\theta(s)}{V_f(s)} = \frac{\frac{K_m}{JL_f}}{s(s+\frac{f}{J})(s+\frac{R_f}{L_f})} = \frac{\frac{K_m}{fR_f}}{s(\tau_f s+1)(\tau_L s+1)}

where τf=LfRf\tau_f=\frac{L_f}{R_f} and τL=Jf\tau_L=\frac{J}{f}.

When τL>τf\tau_L > \tau_f, the field time constant τf\tau_f can be neglected.

Let's choose the following values:

  • (J) moment of inertia of the rotor 0.01 kg.m^2
  • (f) motor viscous friction constant 0.1 N.m.s
  • (Ke) electromotive force constant 0.01 V/rad/sec
  • (Kt) motor torque constant 0.01 N.m/Amp
  • (R) electric resistance 1 Ohm
  • (L) electric inductance 0.5 H

_Note that in SI units, Ke=Kt=KK_e = K_t = K_

θ(s)Vf(s)=0.010.11s(0.5/1s+1)(0.01/0.1s+1)=0.1s(0.5s+1)(0.1s+1) \frac{\theta(s)}{V_f(s)} = \frac{\frac{0.01}{0.1\cdot1}}{s(0.5/1s+1)(0.01/0.1s+1)} = \frac{0.1}{s(0.5s+1)(0.1s+1)}
s = control.TransferFunction.s
G_motor = 0.1/(s*(0.5*s+1)*(0.1*s+1))
T, yout = control.step_response(G_motor, T=np.linspace(0, 7, 1000));

plt.plot(T, yout, 'b', label='dc-motor-step-response')
plt.legend()
plt.xlabel('time')
plt.grid()

9. Draw a block diagram of the field controlled DC motor from the field voltage to the position output using all the relevant equations calculated so far. Make sure you include the disturbance Td(s)T_d(s) as per equation (6).

dc-motor-field-current-bd