61. Closed loop stability¶
Stability of closed loop control systems appears to be easy to determine. We can just calculate the closed loop transfer function and invert the Laplace transform.
import sympy sympy.init_printing()
s = sympy.Symbol('s')
K_c, t = sympy.symbols('K_c, t', positive=True)
These are the systems from Example 10.4 in Seborg et al
G_c = K_c G_v = 1/(2*s + 1) G_p = G_d = 1/(5*s + 1) G_m = 1/(s + 1)
K_m = sympy.limit(G_m, s, 0)
forward = G_c*G_v*G_p backward = G_m G_CL = K_m*forward/(1 + forward*backward)
y = G_CL/s
Now for the moment of truth. Uncomment the next line if you have a lot of time on your hands…
#yt = sympy.inverse_laplace_transform(y, s, t)
So that didn’t really work as we expected. Can we at least calculate the roots of the denominator?
ce = sympy.denom(G_CL.simplify()) ce.expand()
roots = sympy.solve(ce, s)
OK, that approach works, but is limited in the analytic case to 4th order polynomials
sympy.plot(*[sympy.re(r) for r in roots], (K_c, 1, 20))
<sympy.plotting.plot.Plot at 0x11b5e0978>
We can see that one root gets a positive real part around \(K_c=12.5\)
62. Using the control library¶
We quickly run out of SymPy’s abilities when calculating closed loop transfer functions. Let’s try to do it with the control library instead:
import control import numpy import scipy.signal import matplotlib.pyplot as plt %matplotlib inline
s = control.tf([1, 0], )
G_v = 1/(2*s + 1) G_p = G_d = 1/(5*s + 1) G_m = 1/(s + 1) K_m = 1
from tbcontrol.loops import feedback
def closedloop(K_c): G_c = K_c G_CL = K_m*feedback(G_c*G_v*G_p, G_m) return G_CL
from ipywidgets import interact
smootht = numpy.linspace(0, 20)
loop = G_v*G_p*G_m
_ = control.rlocus(loop)
/Users/alchemyst/anaconda3/lib/python3.7/site-packages/matplotlib/figure.py:98: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance. "Adding an axes using the same arguments as a previous axes "
def response(K_C): G_CL = closedloop(K_C) poles = G_CL.pole() plt.plot(*control.step_response(G_CL, smootht)) _ = control.rlocus(loop) plt.scatter(poles.real, poles.imag, color='red')
interact(response, K_C=(1., 20.))
Now, the step response is calculated quickly enough that we can interact with the graphics using the slider!
62.1. Direct substitution¶
From our exploration above it is clear there is a zero of the characteristic equation at \(K_C \approx 13\). Let’s solve for this numerically:
def chareq(x): Kc, omega = x s = 1j*omega ce = 1 + Kc*loop ce_of_s = ce(s) return ce_of_s.real, ce_of_s.imag
scipy.optimize.fsolve(chareq, [13, 1])
array([12.6 , 0.89442719])
How can we determine stability without calculating the roots? See the next notebook.