Calculate the Symbolic Formulas for the PID Gains

Define Symbols

import numpy as np
from sympy import *
from IPython.display import Math

PI = np.pi
K, alpha = symbols('K alpha')                   # actuator symbols. K is a velocity gain and alpha is the corner frequency
Ki, Kp, Kd = symbols('Ki Kp Kd')                # controller symbols.  The 
lam, s = symbols('lambda s')                    # transfer function symbols

Define the controller and motor transfer functions and display them

Gc = Ki/s + Kp + Kd*s                           # PID controller 
Ga = (K*alpha)/(s*(s+alpha))                    # actuator transfer function
display(Math('\nGc(s) = '+latex(Gc)))           # the symbolic formula for the integrator gain
display(Math('\nGa(s) = '+latex(Ga)))           # the symbolic formula for the integrator gain

$\displaystyle Gc(s) = Kd s + \frac{Ki}{s} + Kp$

$\displaystyle Ga(s) = \frac{K \alpha}{s \left(\alpha + s\right)}$

The closed loop transfer function

cltf = simplify((Gc*Ga)/(1+Gc*Ga))              # closed loop transfer function
num, den = fraction(simplify(cltf))             # the denominator is the characteristic equation
display(Math('\ncltf(s) = '+latex(cltf)))       # the symbolic formula for the integrator gain
display(Math('\nnum(s) = '+latex(num)))         # the symbolic formula for the integrator gain
display(Math('\nden(s) = '+latex(den)))         # the symbolic formula for the integrator gain

$\displaystyle cltf(s) = \frac{K \alpha \left(Ki + s \left(Kd s + Kp\right)\right)}{K \alpha \left(Ki + s \left(Kd s + Kp\right)\right) + s^{2} \left(\alpha + s\right)}$

$\displaystyle num(s) = K \alpha \left(Ki + s \left(Kd s + Kp\right)\right)$

$\displaystyle den(s) = K \alpha \left(Ki + s \left(Kd s + Kp\right)\right) + s^{2} \left(\alpha + s\right)$

The desired closed loop transfer function

dcltf = collect((s+lam)**3,s)                   # desired closed loop transfer function, three poles at -lamda
display(Math('desired\\ closed\\ transfer\\ function = '+latex(dcltf)))  # notice how to create spaces between words in LaTeX

$\displaystyle desired\ closed\ transfer\ function = \left(\lambda + s\right)^{3}$

difce = simplify(expand(den-dcltf))             # difference between actual and desired CE
# difce = collect(difce,s)
# print(difce)
print(collect(difce,lam))
collect(difce,lam)

K*Kd*alpha*s**2 + K*Ki*alpha + K*Kp*alpha*s + alpha*s**2 - lambda**3 - 3*lambda**2*s - 3*lambda*s**2

$\displaystyle K Kd \alpha s^{2} + K Ki \alpha + K Kp \alpha s + \alpha s^{2} - \lambda^{3} - 3 \lambda^{2} s - 3 \lambda s^{2}$

Symbolically solve for the controller PID gains

ki = solve(difce.coeff(s,0), Ki)[0]
kp = solve(difce.coeff(s,1), Kp)[0]
kd = solve(difce.coeff(s,2), Kd)[0]
display(Math('\nKi = '+latex(ki)))
display(Math('\nKp = '+latex(kp)))
display(Math('\nKd = '+latex(kd)))

$\displaystyle Ki = \frac{\lambda^{3}}{K \alpha}$

$\displaystyle Kp = \frac{3 \lambda^{2}}{K \alpha}$

$\displaystyle Kd = \frac{- \alpha + 3 \lambda}{K \alpha}$

Assign values to the open loop symbols

_difce = difce.subs([(K,10),(alpha,2*PI*10),(lam,2*PI*10)])
_ki = solve(_difce.coeff(s,0), Ki)[0]         # returns a list of one element
_kp = solve(_difce.coeff(s,1), Kp)[0]
_kd = solve(_difce.coeff(s,2), Kd)[0]
print(f'\nKi={_ki:8.6f}, Kp={_kp:8.6f}, Kd={_kd:8.6f}')

Ki=394.784176, Kp=18.849556, Kd=0.200000

Display the closed loop transfer function's numerator and denominator

_cltf = cltf.subs([(K,10),(alpha,2*PI*10),(lam,2*PI*10), (Ki,_ki), (Kp,_kp), (Kd, _kd)])
_num, _den =  fraction(simplify(_cltf))         # the denominator is the characteristic equation
display(Math('numerator   = '+latex(_num))) 
display(Math('denominator = '+latex(_den)))

$\displaystyle numerator = 125.663706143592 s^{2} + 11843.5252813072 s + 248050.213442399$

$\displaystyle denominator = 1.0 s^{3} + 188.495559215388 s^{2} + 11843.5252813072 s + 248050.213442399$

Display the Bode and Pole Zero Plots

from bodeplot import bodeplot                   # in lib
from pzplot import pzplot                       # in lib   
if _num.as_poly() == None:
    num_coeffs = [float(_num)]
else:
    num_coeffs = list(map(float,_num.as_poly().all_coeffs())) # convert list of symbols to list of floats
den_coeffs = list(map(float,_den.as_poly().all_coeffs())) # convert list of symbols to list of floats
zeros = np.roots(num_coeffs)
poles = np.roots(den_coeffs)
bodeplot(num_coeffs, den_coeffs)
pzplot(poles, zeros)