Skip to content

Bode margins plot #199

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Dec 28, 2018
Merged

Bode margins plot #199

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d95ecf0
added margin display to the bode plots
icam0 Mar 31, 2018
9c07bb6
fixed smalll typo in plot title
icam0 Mar 31, 2018
36bd4e7
added unittest,doc string and fixed deg plot of margins
icam0 Apr 22, 2018
0cf3989
removed faulty system import
icam0 May 25, 2018
991daf4
Merge branch 'master' into bode_margins_plot
repagh Jul 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 50 additions & 3 deletions control/freqplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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']
Expand All @@ -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

Expand All @@ -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)

Expand Down Expand Up @@ -209,7 +212,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
Expand All @@ -219,6 +222,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.,-1e5],color='k', linestyle=':')
else:
ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':')

if deg:
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], [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), -1e5],color='k', linestyle=':')
ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k')
else:
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], [1e-8, phase_limit],color='k', linestyle=':')
else:
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.radians(pm),'deg' if deg else 'rad',Wcp))

if nyquistfrq_plot:
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())

Expand All @@ -237,7 +284,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)
Expand Down
29 changes: 29 additions & 0 deletions control/tests/freqresp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
from control.xferfcn import TransferFunction
from control.matlab import ss, tf, bode, rss
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):
Expand Down Expand Up @@ -108,6 +109,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 test_discrete(self):
# Test discrete time frequency response

Expand Down