Skip to content

Commit 72ffb63

Browse files
icam0icam0
committed
Merge branch 'bode_margins_plot' into sisotool_merge
 i Please enter a commit message to explain why this merge is necessaryskldfj,
2 parents e3c1f31 + 0cf3989 commit 72ffb63

2 files changed

Lines changed: 79 additions & 3 deletions

File tree

‎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
"""Bode plot for a system
6566
6667
Plots a Bode plot for the system over a (optional) frequency range.
@@ -84,6 +85,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
8485
If Hz=True the limits are in Hz otherwise in rad/s.
8586
omega_num: int
8687
number of samples
88+
margins : boolean
89+
if True, plot gain and phase margin
8790
*args, **kwargs:
8891
Additional options to matplotlib (color, linestyle, etc)
8992
@@ -216,7 +219,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
216219
color=pltline[0].get_color())
217220

218221
# Add a grid to the plot + labeling
219-
ax_mag.grid(True, which='both')
222+
ax_mag.grid(False if margins else True, which='both')
220223
ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
221224

222225
# Phase plot
@@ -226,6 +229,50 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
226229
phase_plot = phase
227230
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)
228231

232+
# Show the phase and gain margins in the plot
233+
if margins:
234+
margin = stability_margins(sys)
235+
gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4]
236+
if pm >= 0.:
237+
phase_limit = -180.
238+
else:
239+
phase_limit = 180.
240+
241+
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':')
242+
ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':')
243+
mag_ylim = ax_mag.get_ylim()
244+
phase_ylim = ax_phase.get_ylim()
245+
246+
if pm != float('inf') and Wcp != float('nan'):
247+
if dB:
248+
ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':')
249+
else:
250+
ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':')
251+
252+
if deg:
253+
ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':')
254+
ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k')
255+
else:
256+
ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':')
257+
ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k')
258+
259+
if gm != float('inf') and Wcg != float('nan'):
260+
if dB:
261+
ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':')
262+
ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k')
263+
else:
264+
ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':')
265+
ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k')
266+
267+
if deg:
268+
ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':')
269+
else:
270+
ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':')
271+
272+
ax_mag.set_ylim(mag_ylim)
273+
ax_phase.set_ylim(phase_ylim)
274+
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))
275+
229276
if nyquistfrq_plot:
230277
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())
231278

@@ -244,7 +291,7 @@ def genZeroCenteredSeries(val_min, val_max, period):
244291
ylim = ax_phase.get_ylim()
245292
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.))
246293
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True)
247-
ax_phase.grid(True, which='both')
294+
ax_phase.grid(False if margins else True, which='both')
248295
# ax_mag.grid(which='minor', alpha=0.3)
249296
# ax_mag.grid(which='major', alpha=0.9)
250297
# 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
@@ -12,6 +12,7 @@
1212
from control.statesp import StateSpace
1313
from control.matlab import ss, tf, bode
1414
from control.exception import slycot_check
15+
from control.tests.margin_test import assert_array_almost_equal
1516
import matplotlib.pyplot as plt
1617

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

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

0 commit comments

Comments
 (0)