Skip to content

Sisotool and dynamic root locus zoom #209

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
murrayrm merged 25 commits into from
Apr 16, 2019
Merged

Sisotool and dynamic root locus zoom #209

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
54583d5
first sisotool chechpoint
icam0 May 23, 2018
2e1bf7d
second sisotool checkpoint
icam0 May 24, 2018
f4545dd
changed matplotlib event handler
icam0 May 25, 2018
e3c1f31
improved sisotool readability
icam0 May 25, 2018
72ffb63
Merge branch 'bode_margins_plot' into sisotool_merge
icam0 May 25, 2018
bde60f5
Added display of GM and PM in bode plot. Fixed Hz bug support. Fixed …
icam0 May 28, 2018
bc74a0a
Added support for both matplotlib 1.x.x. and 2.x.x and fixed redrawin…
icam0 May 29, 2018
1e3a925
Added docstring and unittest
icam0 May 29, 2018
94a0ea1
removed bug where a pole or zero would be deleted if only 1 of them w…
icam0 May 30, 2018
c0e533e
root locus smooth with zoom bug
icam0 Jun 2, 2018
fe13ae1
Fixed zooming on plot. Now displays smoother root locus plot.
icam0 Jun 2, 2018
2ba8277
Fixed paning on plot and clicking on the root locus plot while zoomed
icam0 Jun 3, 2018
c26187f
Added unit test for rootlocus zoom functionality
icam0 Jun 3, 2018
9962cea
Fixed sisotool unittest
icam0 Jun 3, 2018
938ae88
Fixed frequency response unittest
icam0 Jun 3, 2018
9f13a77
small change to sisotool test
icam0 Jun 3, 2018
6e06a9c
small change to unittests
icam0 Jun 3, 2018
cdb9f42
fixed python 2 integer division error in false gain of rootlocus
icam0 Jun 3, 2018
90959bd
changed root locus test to a numerically more significant one.
icam0 Jun 3, 2018
70ec8be
changed root locus test indices.
icam0 Jun 3, 2018
9f50e13
Merge branch 'master' into sisotool_final
repagh Jul 9, 2018
3e7cad6
Merge branch 'master' into sisotool_final
murrayrm Dec 29, 2018
43d8b95
Merge branch 'master' into sisotool_final
murrayrm Apr 11, 2019
03f87a3
Fix small typo from manual merge
murrayrm Apr 12, 2019
17c4a95
fixed issues identified in review
murrayrm Apr 12, 2019
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
1 change: 1 addition & 0 deletions control/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@
from .canonical import *
from .robust import *
from .config import *
from .sisotool import *

# Exceptions
from .exception import *
Expand Down
184 changes: 129 additions & 55 deletions control/freqplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
# SUCH DAMAGE.
#
# $Id$

import matplotlib
import matplotlib.pyplot as plt
import scipy as sp
import numpy as np
Expand Down Expand Up @@ -89,7 +89,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
omega_num: int
number of samples
margins : boolean
if True, plot gain and phase margin
If True, plot gain and phase margin
\*args, \**kwargs:
Additional options to matplotlib (color, linestyle, etc)

Expand Down Expand Up @@ -193,23 +193,35 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
# The code below should work on all cases.

# Get the current figure
fig = plt.gcf()
ax_mag = None
ax_phase = None

# Get the current axes if they already exist
for ax in fig.axes:
if ax.get_label() == 'control-bode-magnitude':
ax_mag = ax
elif ax.get_label() == 'control-bode-phase':
ax_phase = ax

# If no axes present, create them from scratch
if ax_mag is None or ax_phase is None:
plt.clf()
ax_mag = plt.subplot(211, label='control-bode-magnitude')
ax_phase = plt.subplot(212, label='control-bode-phase',
sharex=ax_mag)

if 'sisotool' in kwargs:
fig = kwargs['fig']
ax_mag = fig.axes[0]
ax_phase = fig.axes[2]
sisotool = kwargs['sisotool']
del kwargs['fig']
del kwargs['sisotool']
else:
fig = plt.gcf()
ax_mag = None
ax_phase = None
sisotool = False

# Get the current axes if they already exist
for ax in fig.axes:
if ax.get_label() == 'control-bode-magnitude':
ax_mag = ax
elif ax.get_label() == 'control-bode-phase':
ax_phase = ax

# If no axes present, create them from scratch
if ax_mag is None or ax_phase is None:
plt.clf()
ax_mag = plt.subplot(211,
label='control-bode-magnitude')
ax_phase = plt.subplot(212,
label='control-bode-phase',
sharex=ax_mag)

# Magnitude plot
if dB:
Expand Down Expand Up @@ -237,58 +249,107 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
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:
# TODO: add some documentation describing why this is here
phase_at_cp = phases[0][(np.abs(omegas[0] - Wcp)).argmin()]
if phase_at_cp >= 0.:
phase_limit = 180.
else:
phase_limit = -180.

if Hz:
Wcg, Wcp = Wcg/(2*math.pi),Wcp/(2*math.pi)

ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':')
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',
zorder=-20)
ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit),
color='k', linestyle=':')
color='k', linestyle=':', zorder=-20)
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=':')
ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],
color='k', linestyle=':',
zorder=-20)
else:
ax_mag.loglog([Wcp, Wcp], [1., 1e-8], color='k', linestyle=':')
ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',
linestyle=':', zorder=-20)

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')
ax_phase.semilogx([Wcp, Wcp],
[1e5, phase_limit+pm],
color='k', linestyle=':',
zorder=-20)
ax_phase.semilogx([Wcp, Wcp],
[phase_limit + pm, phase_limit],
color='k', zorder=-20)
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')
ax_phase.semilogx([Wcp, Wcp],
[1e5, math.radians(phase_limit) +
math.radians(pm)],
color='k', linestyle=':',
zorder=-20)
ax_phase.semilogx([Wcp, Wcp],
[math.radians(phase_limit) +
math.radians(pm),
math.radians(phase_limit)],
color='k', zorder=-20)

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')
ax_mag.semilogx([Wcg, Wcg],
[-20.*np.log10(gm), -1e5],
color='k', linestyle=':',
zorder=-20)
ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],
color='k', zorder=-20)
else:
ax_mag.loglog([Wcg, Wcg], [1. / gm, 1e-8], color='k', linestyle=':')
ax_mag.loglog([Wcg, Wcg], [1., 1. / gm], color='k')
ax_mag.loglog([Wcg, Wcg],
[1./gm,1e-8],color='k',
linestyle=':', zorder=-20)
ax_mag.loglog([Wcg, Wcg],
[1.,1./gm],color='k', zorder=-20)

if deg:
ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],
color='k', linestyle=':')
ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],
color='k', linestyle=':',
zorder=-20)
else:
ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],
color='k', linestyle=':')
ax_phase.semilogx([Wcg, Wcg],
[1e-8, math.radians(phase_limit)],
color='k', linestyle=':',
zorder=-20)

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 sisotool:
ax_mag.text(0.04, 0.06,
'G.M.: %.2f %s\nFreq: %.2f %s' %
(20*np.log10(gm) if dB else gm,
'dB ' if dB else '',
Wcg, 'Hz' if Hz else 'rad/s'),
horizontalalignment='left',
verticalalignment='bottom',
transform=ax_mag.transAxes,
fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6)
ax_phase.text(0.04, 0.06,
'P.M.: %.2f %s\nFreq: %.2f %s' %
(pm if deg else math.radians(pm),
'deg' if deg else 'rad',
Wcp, 'Hz' if Hz else 'rad/s'),
horizontalalignment='left',
verticalalignment='bottom',
transform=ax_phase.transAxes,
fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6)
else:
plt.suptitle('Gm = %.2f %s(at %.2f %s), Pm = %.2f %s (at %.2f %s)' %
(20*np.log10(gm) if dB else gm,
'dB ' if dB else '\b',
Wcg, 'Hz' if Hz else 'rad/s',
pm if deg else math.radians(pm),
'deg' if deg else 'rad',
Wcp, 'Hz' if Hz else 'rad/s'))

if nyquistfrq_plot:
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())
Expand All @@ -302,12 +363,19 @@ def gen_zero_centered_series(val_min, val_max, period):
return np.arange(v1, v2 + 1) * period
if deg:
ylim = ax_phase.get_ylim()
ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], 45.))
ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], 15.), minor=True)
ax_phase.set_yticks(gen_zero_centered_series(ylim[0],
ylim[1], 45.))
ax_phase.set_yticks(gen_zero_centered_series(ylim[0],
ylim[1], 15.),
minor=True)
else:
ylim = ax_phase.get_ylim()
ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], math.pi / 4.))
ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], math.pi / 12.),
ax_phase.set_yticks(gen_zero_centered_series(ylim[0],
ylim[1],
math.pi / 4.))
ax_phase.set_yticks(gen_zero_centered_series(ylim[0],
ylim[1],
math.pi / 12.),
minor=True)
ax_phase.grid(False if margins else True, which='both')
# ax_mag.grid(which='minor', alpha=0.3)
Expand All @@ -316,15 +384,17 @@ def gen_zero_centered_series(val_min, val_max, period):
# ax_phase.grid(which='major', alpha=0.9)

# Label the frequency axis
ax_phase.set_xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)")
ax_phase.set_xlabel("Frequency (Hz)" if Hz
else "Frequency (rad/sec)")

if len(syslist) == 1:
return mags[0], phases[0], omegas[0]
else:
return mags, phases, omegas


def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, **kwargs):
def nyquist_plot(syslist, omega=None, Plot=True, color=None,
labelFreq=0, *args, **kwargs):
"""
Nyquist plot for a system

Expand Down Expand Up @@ -683,6 +753,10 @@ def gen_prefix(pow1000):
'z', # zepto (10^-21)
'y'][8 - pow1000] # yocto (10^-24)

def find_nearest_omega(omega_list, omega):
omega_list = np.asarray(omega_list)
idx = (np.abs(omega_list - omega)).argmin()
return omega_list[(np.abs(omega_list - omega)).argmin()]

# Function aliases
bode = bode_plot
Expand Down
Loading