From d95ecf0abff11d312e6385b7effc337bb256b651 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 1 Apr 2018 00:18:43 +0200 Subject: [PATCH 1/4] added margin display to the bode plots --- control/freqplot.py | 51 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 48 insertions(+), 3 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index a71d44cce..226a093b8 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -47,6 +47,7 @@ import math from .ctrlutil import unwrap from .bdalg import feedback +from .margins import stability_margins __all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot', 'bode', 'nyquist', 'gangof4'] @@ -60,7 +61,7 @@ # Bode plot def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, - Plot=True, omega_limits=None, omega_num=None, *args, **kwargs): + Plot=True, omega_limits=None, omega_num=None,margins=None, *args, **kwargs): """Bode plot for a system Plots a Bode plot for the system over a (optional) frequency range. @@ -208,7 +209,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, color=pltline[0].get_color()) # Add a grid to the plot + labeling - ax_mag.grid(True, which='both') + ax_mag.grid(False if margins else True, which='both') ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude") # Phase plot @@ -218,6 +219,50 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, phase_plot = phase ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs) + #show the phase and gain margins in the plot + if margins: + margin = stability_margins(sys) + gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4] + if pm >= 0.: + phase_limit = -180. + else: + phase_limit = 180. + + ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':') + ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':') + mag_ylim = ax_mag.get_ylim() + phase_ylim = ax_phase.get_ylim() + + if pm != float('inf') and Wcp != float('nan'): + if dB: + ax_mag.semilogx([Wcp, Wcp], [0, -(10**5.)], color='k', linestyle=':') + else: + ax_mag.loglog([Wcp,Wcp], [1,10.**-20.],color='k',linestyle=':') + + if deg: + ax_phase.semilogx([Wcp, Wcp], [(10**5.) if phase_limit > 0. else 10.**-20., phase_limit+pm], color='k', linestyle=':') + ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit], color='k', linestyle='-') + else: + ax_phase.semilogx([Wcp, Wcp], [10.**-20., -math.pi+math.radians(pm)], color='k', linestyle=':') + ax_phase.semilogx([Wcp, Wcp], [-math.pi +math.radians(pm), -math.pi], color='k', linestyle='-') + + if gm != float('inf') and Wcg != float('nan'): + if dB: + ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -(10**5.)], color='k', linestyle=':') + ax_mag.semilogx([Wcg, Wcg], [0,- 20 * np.log10(gm)], color='k', linestyle='-') + else: + ax_mag.semilogx([Wcg, Wcg], [1./gm, 10. ** -20.], color='k', linestyle=':') + ax_mag.semilogx([Wcg, Wcg], [1., 1./gm], color='k', linestyle='-') + + if deg: + ax_phase.semilogx([Wcg, Wcg], [10.**-20., phase_limit], color='k', linestyle=':') + else: + ax_phase.semilogx([Wcg, Wcg], [10. ** -20., -math.pi], color='k', linestyle=':') + + ax_mag.set_ylim(mag_ylim) + ax_phase.set_ylim(phase_ylim) + 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.degrees(pm),'deg' if deg else 'rad',Wcp)) + if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) @@ -236,7 +281,7 @@ def genZeroCenteredSeries(val_min, val_max, period): ylim = ax_phase.get_ylim() ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.)) ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True) - ax_phase.grid(True, which='both') + ax_phase.grid(False if margins else True, which='both') # ax_mag.grid(which='minor', alpha=0.3) # ax_mag.grid(which='major', alpha=0.9) # ax_phase.grid(which='minor', alpha=0.3) From 9c07bb6b1731ef7995f3a8826e84a5786010918a Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 1 Apr 2018 00:26:59 +0200 Subject: [PATCH 2/4] fixed smalll typo in plot title --- control/freqplot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/freqplot.py b/control/freqplot.py index 226a093b8..a1147a502 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -261,7 +261,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, ax_mag.set_ylim(mag_ylim) ax_phase.set_ylim(phase_ylim) - 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.degrees(pm),'deg' if deg else 'rad',Wcp)) + 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.degrees(pm),'deg' if deg else 'rad',Wcp)) if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) From 36bd4e7b30a34a747809fec5eb411dbe78fb5fa1 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 22 Apr 2018 13:20:44 +0200 Subject: [PATCH 3/4] added unittest,doc string and fixed deg plot of margins --- control/freqplot.py | 31 +++++++++++++++++-------------- control/tests/freqresp_test.py | 29 +++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index a1147a502..5cc426c22 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -45,6 +45,7 @@ import scipy as sp import numpy as np import math +import sys as system from .ctrlutil import unwrap from .bdalg import feedback from .margins import stability_margins @@ -85,6 +86,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, If Hz=True the limits are in Hz otherwise in rad/s. omega_num: int number of samples + margins : boolean + if True, plot gain and phase margin *args, **kwargs: Additional options to matplotlib (color, linestyle, etc) @@ -219,7 +222,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, phase_plot = phase ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs) - #show the phase and gain margins in the plot + # Show the phase and gain margins in the plot if margins: margin = stability_margins(sys) gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4] @@ -235,33 +238,33 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if pm != float('inf') and Wcp != float('nan'): if dB: - ax_mag.semilogx([Wcp, Wcp], [0, -(10**5.)], color='k', linestyle=':') + ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':') else: - ax_mag.loglog([Wcp,Wcp], [1,10.**-20.],color='k',linestyle=':') + ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':') if deg: - ax_phase.semilogx([Wcp, Wcp], [(10**5.) if phase_limit > 0. else 10.**-20., phase_limit+pm], color='k', linestyle=':') - ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit], color='k', linestyle='-') + ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':') + ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k') else: - ax_phase.semilogx([Wcp, Wcp], [10.**-20., -math.pi+math.radians(pm)], color='k', linestyle=':') - ax_phase.semilogx([Wcp, Wcp], [-math.pi +math.radians(pm), -math.pi], color='k', linestyle='-') + ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':') + ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k') if gm != float('inf') and Wcg != float('nan'): if dB: - ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -(10**5.)], color='k', linestyle=':') - ax_mag.semilogx([Wcg, Wcg], [0,- 20 * np.log10(gm)], color='k', linestyle='-') + ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':') + ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k') else: - ax_mag.semilogx([Wcg, Wcg], [1./gm, 10. ** -20.], color='k', linestyle=':') - ax_mag.semilogx([Wcg, Wcg], [1., 1./gm], color='k', linestyle='-') + ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':') + ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k') if deg: - ax_phase.semilogx([Wcg, Wcg], [10.**-20., phase_limit], color='k', linestyle=':') + ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':') else: - ax_phase.semilogx([Wcg, Wcg], [10. ** -20., -math.pi], color='k', linestyle=':') + ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':') ax_mag.set_ylim(mag_ylim) ax_phase.set_ylim(phase_ylim) - 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.degrees(pm),'deg' if deg else 'rad',Wcp)) + 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)) if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py index 9e8dd44ce..781c40496 100644 --- a/control/tests/freqresp_test.py +++ b/control/tests/freqresp_test.py @@ -12,6 +12,7 @@ from control.statesp import StateSpace from control.matlab import ss, tf, bode from control.exception import slycot_check +from control.tests.margin_test import assert_array_almost_equal import matplotlib.pyplot as plt class TestFreqresp(unittest.TestCase): @@ -107,6 +108,34 @@ def test_mimo(self): #plt.figure(4) #bode(sysMIMO,self.omega) + def test_bode_margin(self): + num = [1000] + den = [1, 25, 100, 0] + sys = ctrl.tf(num, den) + ctrl.bode_plot(sys, margins=True,dB=False,deg = True) + fig = plt.gcf() + allaxes = fig.get_axes() + + mag_to_infinity = (np.array([6.07828691, 6.07828691]), + np.array([1.00000000e+00, 1.00000000e-08])) + assert_array_almost_equal(mag_to_infinity, allaxes[0].lines[2].get_data()) + + gm_to_infinty = (np.array([10., 10.]), np.array([4.00000000e-01, 1.00000000e-08])) + assert_array_almost_equal(gm_to_infinty, allaxes[0].lines[3].get_data()) + + one_to_gm = (np.array([10., 10.]), np.array([1., 0.4])) + assert_array_almost_equal(one_to_gm, allaxes[0].lines[4].get_data()) + + pm_to_infinity = (np.array([6.07828691, 6.07828691]), + np.array([100000., -157.46405841])) + assert_array_almost_equal(pm_to_infinity, allaxes[1].lines[2].get_data()) + + pm_to_phase = (np.array([6.07828691, 6.07828691]), np.array([-157.46405841, -180.])) + assert_array_almost_equal(pm_to_phase, allaxes[1].lines[3].get_data()) + + phase_to_infinity = (np.array([10., 10.]), np.array([1.00000000e-08, -1.80000000e+02])) + assert_array_almost_equal(phase_to_infinity, allaxes[1].lines[4].get_data()) + def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp) From 0cf39895f398fb977cbe384a109e1b412e8861b8 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Fri, 25 May 2018 16:48:09 +0200 Subject: [PATCH 4/4] removed faulty system import --- control/freqplot.py | 1 - 1 file changed, 1 deletion(-) diff --git a/control/freqplot.py b/control/freqplot.py index 5cc426c22..0686b2162 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -45,7 +45,6 @@ import scipy as sp import numpy as np import math -import sys as system from .ctrlutil import unwrap from .bdalg import feedback from .margins import stability_margins