# 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()

:

%matplotlib inline

:

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)

:

sympy.simplify(G_CL)

:

$$\frac{K_{c} \left(s + 1\right)}{K_{c} + \left(s + 1\right) \left(2 s + 1\right) \left(5 s + 1\right)}$$
:

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()

:

$$K_{c} + 10 s^{3} + 17 s^{2} + 8 s + 1$$
:

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

:

closedloop(2)

:

$$\frac{20 s^3 + 34 s^2 + 16 s + 2}{100 s^5 + 240 s^4 + 209 s^3 + 103 s^2 + 29 s + 3}$$
:

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.))

:

<function __main__.response(K_C)>


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

:

import scipy.optimize

:

scipy.optimize.fsolve(chareq, [13, 1])

:

array([12.6       ,  0.89442719])


How can we determine stability without calculating the roots? See the next notebook.