Skip to content

Commit 4286329

Browse files
murrayrmmurrayrm
authored
Merge pull request #199 from icam0/bode_margins_plot
Add margin computation/plotting to bode_plot().
2 parents 32b0c2b + 991daf4 commit 4286329

File tree

2 files changed

+79
-3
lines changed

2 files changed

+79
-3
lines changed

‎control/freqplot.py

Lines changed: 50 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@
4747
import math
4848
from .ctrlutil import unwrap
4949
from .bdalg import feedback
50+
from .margins import stability_margins
5051

5152
__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
5253
'bode', 'nyquist', 'gangof4']
@@ -60,7 +61,7 @@
6061

6162
# Bode plot
6263
def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
63-
Plot=True, omega_limits=None, omega_num=None, *args, **kwargs):
64+
Plot=True, omega_limits=None, omega_num=None,margins=None, *args, **kwargs):
6465
"""
6566
Bode plot for a system
6667
@@ -85,6 +86,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
8586
If Hz=True the limits are in Hz otherwise in rad/s.
8687
omega_num: int
8788
number of samples
89+
margins : boolean
90+
if True, plot gain and phase margin
8891
\*args, \**kwargs:
8992
Additional options to matplotlib (color, linestyle, etc)
9093
@@ -209,7 +212,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
209212
color=pltline[0].get_color())
210213

211214
# Add a grid to the plot + labeling
212-
ax_mag.grid(True, which='both')
215+
ax_mag.grid(False if margins else True, which='both')
213216
ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
214217

215218
# Phase plot
@@ -219,6 +222,50 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
219222
phase_plot = phase
220223
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)
221224

225+
# Show the phase and gain margins in the plot
226+
if margins:
227+
margin = stability_margins(sys)
228+
gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4]
229+
if pm >= 0.:
230+
phase_limit = -180.
231+
else:
232+
phase_limit = 180.
233+
234+
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':')
235+
ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':')
236+
mag_ylim = ax_mag.get_ylim()
237+
phase_ylim = ax_phase.get_ylim()
238+
239+
if pm != float('inf') and Wcp != float('nan'):
240+
if dB:
241+
ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':')
242+
else:
243+
ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':')
244+
245+
if deg:
246+
ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':')
247+
ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k')
248+
else:
249+
ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':')
250+
ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k')
251+
252+
if gm != float('inf') and Wcg != float('nan'):
253+
if dB:
254+
ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':')
255+
ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k')
256+
else:
257+
ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':')
258+
ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k')
259+
260+
if deg:
261+
ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':')
262+
else:
263+
ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':')
264+
265+
ax_mag.set_ylim(mag_ylim)
266+
ax_phase.set_ylim(phase_ylim)
267+
plt.suptitle('Gm = %.2f %s(at %.2f rad/s), Pm = %.2f %s (at %.2f rad/s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp))
268+
222269
if nyquistfrq_plot:
223270
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())
224271

@@ -237,7 +284,7 @@ def genZeroCenteredSeries(val_min, val_max, period):
237284
ylim = ax_phase.get_ylim()
238285
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.))
239286
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True)
240-
ax_phase.grid(True, which='both')
287+
ax_phase.grid(False if margins else True, which='both')
241288
# ax_mag.grid(which='minor', alpha=0.3)
242289
# ax_mag.grid(which='major', alpha=0.9)
243290
# ax_phase.grid(which='minor', alpha=0.3)

‎control/tests/freqresp_test.py

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
from control.xferfcn import TransferFunction
1414
from control.matlab import ss, tf, bode, rss
1515
from control.exception import slycot_check
16+
from control.tests.margin_test import assert_array_almost_equal
1617
import matplotlib.pyplot as plt
1718

1819
class TestFreqresp(unittest.TestCase):
@@ -108,6 +109,34 @@ def test_mimo(self):
108109
#plt.figure(4)
109110
#bode(sysMIMO,self.omega)
110111

112+
def test_bode_margin(self):
113+
num = [1000]
114+
den = [1, 25, 100, 0]
115+
sys = ctrl.tf(num, den)
116+
ctrl.bode_plot(sys, margins=True,dB=False,deg = True)
117+
fig = plt.gcf()
118+
allaxes = fig.get_axes()
119+
120+
mag_to_infinity = (np.array([6.07828691, 6.07828691]),
121+
np.array([1.00000000e+00, 1.00000000e-08]))
122+
assert_array_almost_equal(mag_to_infinity, allaxes[0].lines[2].get_data())
123+
124+
gm_to_infinty = (np.array([10., 10.]), np.array([4.00000000e-01, 1.00000000e-08]))
125+
assert_array_almost_equal(gm_to_infinty, allaxes[0].lines[3].get_data())
126+
127+
one_to_gm = (np.array([10., 10.]), np.array([1., 0.4]))
128+
assert_array_almost_equal(one_to_gm, allaxes[0].lines[4].get_data())
129+
130+
pm_to_infinity = (np.array([6.07828691, 6.07828691]),
131+
np.array([100000., -157.46405841]))
132+
assert_array_almost_equal(pm_to_infinity, allaxes[1].lines[2].get_data())
133+
134+
pm_to_phase = (np.array([6.07828691, 6.07828691]), np.array([-157.46405841, -180.]))
135+
assert_array_almost_equal(pm_to_phase, allaxes[1].lines[3].get_data())
136+
137+
phase_to_infinity = (np.array([10., 10.]), np.array([1.00000000e-08, -1.80000000e+02]))
138+
assert_array_almost_equal(phase_to_infinity, allaxes[1].lines[4].get_data())
139+
111140
def test_discrete(self):
112141
# Test discrete time frequency response
113142

0 commit comments

Comments
 (0)