import control
import matplotlib.pyplot as plt

Gp = control.tf([2],[1,3,2,0])
Gz = control.c2d(Gp,0.05)

mag, phase, omega = control.bode(Gz)

print(control.margin(Gz))
print(control.margin([mag, phase, omega]))

plt.grid()
plt.show()

(31.505480557979414, inf, 3.7736323313549196, nan)
(inf, 177.40890111291884, nan, 0.749329697586152)

import control as ct 
import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import interpolate

Gp = ct.tf([2],[1,3,2,0])
Gz = ct.c2d(Gp, 0.05)

mag, phase, omega = ct.bode(Gz)

p_w = interpolate.interp1d(phase, omega)
p_m = interpolate.interp1d(phase, mag)
print("Gain margin %.3f at %.3f rad/s"%(1/p_m(-np.pi), p_w(-np.pi)))

m_w = interpolate.interp1d(mag, omega)
m_p = interpolate.interp1d(mag, phase)
print("Phase margin %.3f\N{DEGREE SIGN} at %.3f rad/s"%(math.degrees(m_p(1)+np.pi), m_w(1)))

Gain margin 2.771 at 1.367 rad/s
Phase margin 31.394° at 0.753 rad/s