How can I plot the margins in a python bode plot?

Question

I want to plot a bode plot of a system with the python control systems library. This is fairly easy. The problem is the plot of the margins. It is no problem to plot the phase margin. But how can I plot the gain margin?

So far, this is a part of my code:

import control as cn
%matplotlib notebook
import matplotlib.pyplot as plt

Ks=2
T1=5
T2=0.3
T3=0.1
Gs=cn.tf(Ks,[T1*T2*T3, T1*T2+T1*T3+T2*T3, T1+T2+T3, 1])

Vr=1
Tn=1
plt.close()
R=cn.tf([Vr*Tn, Vr],[1, 0])
L=Gs*R
gm, pm, wg, wp = cn.margin(L)

_,_,_ = cn.bode(L,dB=True)
plt.axvline(x = wp,color='r')

Answer 1

Not the most elegant solution but hey it works for me.

###Import modules
import numpy as np
import control as ctl
import matplotlib.pyplot as plt

##Functions
def plot_margins(sys):
    mag,phase,omega = ctl.bode(sys,dB=True,Plot=False)
    magdB = 20*np.log10(mag)
    phase_deg = phase*180.0/np.pi
    Gm,Pm,Wcg,Wcp = ctl.margin(sys)
    GmdB = 20*np.log10(Gm)
    ##Plot Gain and Phase
    f,(ax1,ax2) = plt.subplots(2,1)
    ax1.semilogx(omega,magdB)
    ax1.grid(which="both")
    ax1.set_xlabel('Frequency (rad/s)')
    ax1.set_ylabel('Magnitude (dB)')
    ax2.semilogx(omega,phase_deg)
    ax2.grid(which="both")
    ax2.set_xlabel('Frequency (rad/s)')
    ax2.set_ylabel('Phase (deg)')
    ax1.set_title('Gm = '+str(np.round(GmdB,2))+' dB (at '+str(np.round(Wcg,2))+' rad/s), Pm = '+str(np.round(Pm,2))+' deg (at '+str(np.round(Wcp,2))+' rad/s)')
    ###Plot the zero dB line
    ax1.plot(omega,0*omega,'k--',lineWidth=2)
    ###Plot the -180 deg lin
    ax2.plot(omega,-180+0*omega,'k--',lineWidth=2)
    ##Plot the vertical line from -180 to 0 at Wcg
    ax2.plot([Wcg,Wcg],[-180,0],'r--',lineWidth=2)
    ##Plot the vertical line from -180+Pm to 0 at Wcp
    ax2.plot([Wcp,Wcp],[-180+Pm,0],'g--',lineWidth=2)
    ##Plot the vertical line from min(magdB) to 0-GmdB at Wcg
    ax1.plot([Wcg,Wcg],[np.min(magdB),0-GmdB],'r--',lineWidth=2)
    ##Plot the vertical line from min(magdB) to 0db at Wcp
    ax1.plot([Wcp,Wcp],[np.min(magdB),0],'g--',lineWidth=2)
    return Gm,Pm,Wcg,Wcp

#%%%Actuator Dynamics
G = ctl.tf([1],[1,2,1,0])
Gm,Pm,Wcg,Wcp=plot_margins(G)
plt.show()