Skip to content

Commit eb0f3f9

Browse files
murrayrmmurrayrm
authored
Merge pull request #953 from murrayrm/pzmap_plots-22Jul2023
Update pole/zero and root locus plots to use _map/_plot pattern
2 parents 6ed2e77 + f500cbe commit eb0f3f9

19 files changed

+1256
-705
lines changed

‎control/freqplot.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -162,7 +162,7 @@ def bode_plot(
162162
values with no plot.
163163
rcParams : dict
164164
Override the default parameters used for generating plots.
165-
Default is set up config.default['freqplot.rcParams'].
165+
Default is set by config.default['freqplot.rcParams'].
166166
wrap_phase : bool or float
167167
If wrap_phase is `False` (default), then the phase will be unwrapped
168168
so that it is continuously increasing or decreasing. If wrap_phase is

‎control/grid.py‎

Lines changed: 52 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,21 @@
1-
import numpy as np
2-
from numpy import cos, sin, sqrt, linspace, pi, exp
1+
# grid.py - code to add gridlines to root locus and pole-zero diagrams
2+
#
3+
# This code generates grids for pole-zero diagrams (including root locus
4+
# diagrams). Rather than just draw a grid in place, it uses the AxisArtist
5+
# package to generate a custom grid that will scale with the figure.
6+
#
7+
38
import matplotlib.pyplot as plt
4-
from mpl_toolkits.axisartist import SubplotHost
5-
from mpl_toolkits.axisartist.grid_helper_curvelinear \
6-
import GridHelperCurveLinear
79
import mpl_toolkits.axisartist.angle_helper as angle_helper
10+
import numpy as np
811
from matplotlib.projections import PolarAxes
912
from matplotlib.transforms import Affine2D
13+
from mpl_toolkits.axisartist import SubplotHost
14+
from mpl_toolkits.axisartist.grid_helper_curvelinear import \
15+
GridHelperCurveLinear
16+
from numpy import cos, exp, linspace, pi, sin, sqrt
17+
18+
from .iosys import isdtime
1019

1120

1221
class FormatterDMS(object):
@@ -65,14 +74,15 @@ def __call__(self, transform_xy, x1, y1, x2, y2):
6574
return lon_min, lon_max, lat_min, lat_max
6675

6776

68-
def sgrid():
77+
def sgrid(scaling=None):
6978
# From matplotlib demos:
7079
# https://matplotlib.org/gallery/axisartist/demo_curvelinear_grid.html
7180
# https://matplotlib.org/gallery/axisartist/demo_floating_axis.html
7281

7382
# PolarAxes.PolarTransform takes radian. However, we want our coordinate
74-
# system in degree
83+
# system in degrees
7584
tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()
85+
7686
# polar projection, which involves cycle, and also has limits in
7787
# its coordinates, needs a special method to find the extremes
7888
# (min, max of the coordinate within the view).
@@ -89,6 +99,7 @@ def sgrid():
8999
tr, extreme_finder=extreme_finder, grid_locator1=grid_locator1,
90100
tick_formatter1=tick_formatter1)
91101

102+
# Set up an axes with a specialized grid helper
92103
fig = plt.gcf()
93104
ax = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)
94105

@@ -97,15 +108,20 @@ def sgrid():
97108
ax.axis[:].major_ticklabels.set_visible(visible)
98109
ax.axis[:].major_ticks.set_visible(False)
99110
ax.axis[:].invert_ticklabel_direction()
111+
ax.axis[:].major_ticklabels.set_color('gray')
100112

113+
# Set up internal tickmarks and labels along the real/imag axes
101114
ax.axis["wnxneg"] = axis = ax.new_floating_axis(0, 180)
102115
axis.set_ticklabel_direction("-")
103116
axis.label.set_visible(False)
117+
104118
ax.axis["wnxpos"] = axis = ax.new_floating_axis(0, 0)
105119
axis.label.set_visible(False)
120+
106121
ax.axis["wnypos"] = axis = ax.new_floating_axis(0, 90)
107122
axis.label.set_visible(False)
108-
axis.set_axis_direction("left")
123+
axis.set_axis_direction("right")
124+
109125
ax.axis["wnyneg"] = axis = ax.new_floating_axis(0, 270)
110126
axis.label.set_visible(False)
111127
axis.set_axis_direction("left")
@@ -119,43 +135,41 @@ def sgrid():
119135
ax.axis["bottom"].get_helper().nth_coord_ticks = 0
120136

121137
fig.add_subplot(ax)
122-
123-
# RECTANGULAR X Y AXES WITH SCALE
124-
# par2 = ax.twiny()
125-
# par2.axis["top"].toggle(all=False)
126-
# par2.axis["right"].toggle(all=False)
127-
# new_fixed_axis = par2.get_grid_helper().new_fixed_axis
128-
# par2.axis["left"] = new_fixed_axis(loc="left",
129-
# axes=par2,
130-
# offset=(0, 0))
131-
# par2.axis["bottom"] = new_fixed_axis(loc="bottom",
132-
# axes=par2,
133-
# offset=(0, 0))
134-
# FINISH RECTANGULAR
135-
136138
ax.grid(True, zorder=0, linestyle='dotted')
137139

138-
_final_setup(ax)
140+
_final_setup(ax, scaling=scaling)
139141
return ax, fig
140142

141143

142-
def _final_setup(ax):
144+
# Utility function used by all grid code
145+
def _final_setup(ax, scaling=None):
143146
ax.set_xlabel('Real')
144147
ax.set_ylabel('Imaginary')
145-
ax.axhline(y=0, color='black', lw=1)
146-
ax.axvline(x=0, color='black', lw=1)
147-
plt.axis('equal')
148+
ax.axhline(y=0, color='black', lw=0.25)
149+
ax.axvline(x=0, color='black', lw=0.25)
148150

151+
# Set up the scaling for the axes
152+
scaling = 'equal' if scaling is None else scaling
153+
plt.axis(scaling)
149154

150-
def nogrid():
151-
f = plt.gcf()
152-
ax = plt.axes()
153155

154-
_final_setup(ax)
155-
return ax, f
156+
# If not grid is given, at least separate stable/unstable regions
157+
def nogrid(dt=None, ax=None, scaling=None):
158+
fig = plt.gcf()
159+
if ax is None:
160+
ax = fig.gca()
161+
162+
# Draw the unit circle for discrete time systems
163+
if isdtime(dt=dt, strict=True):
164+
s = np.linspace(0, 2*pi, 100)
165+
ax.plot(np.cos(s), np.sin(s), 'k--', lw=0.5, dashes=(5, 5))
156166

167+
_final_setup(ax, scaling=scaling)
168+
return ax, fig
157169

158-
def zgrid(zetas=None, wns=None, ax=None):
170+
# Grid for discrete time system (drawn, not rendered by AxisArtist)
171+
# TODO (at some point): think about using customized grid generator?
172+
def zgrid(zetas=None, wns=None, ax=None, scaling=None):
159173
"""Draws discrete damping and frequency grid"""
160174

161175
fig = plt.gcf()
@@ -206,5 +220,9 @@ def zgrid(zetas=None, wns=None, ax=None):
206220
ax.annotate(r"$\frac{"+num+r"\pi}{T}$", xy=(an_x, an_y),
207221
xytext=(an_x, an_y), size=9)
208222

209-
_final_setup(ax)
223+
# Set default axes to allow some room around the unit circle
224+
ax.set_xlim([-1.1, 1.1])
225+
ax.set_ylim([-1.1, 1.1])
226+
227+
_final_setup(ax, scaling=scaling)
210228
return ax, fig

‎control/iosys.py‎

Lines changed: 35 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -503,42 +503,64 @@ def common_timebase(dt1, dt2):
503503
raise ValueError("Systems have incompatible timebases")
504504

505505
# Check to see if a system is a discrete time system
506-
def isdtime(sys, strict=False):
506+
def isdtime(sys=None, strict=False, dt=None):
507507
"""
508508
Check to see if a system is a discrete time system.
509509
510510
Parameters
511511
----------
512-
sys : I/O or LTI system
513-
System to be checked
514-
strict: bool (default = False)
515-
If strict is True, make sure that timebase is not None
512+
sys : I/O system, optional
513+
System to be checked.
514+
dt : None or number, optional
515+
Timebase to be checked.
516+
strict: bool, default=False
517+
If strict is True, make sure that timebase is not None.
516518
"""
517519

518-
# Check to see if this is a constant
520+
# See if we were passed a timebase instead of a system
521+
if sys is None:
522+
if dt is None:
523+
return True if not strict else False
524+
else:
525+
return dt > 0
526+
elif dt is not None:
527+
raise TypeError("passing both system and timebase not allowed")
528+
529+
# Check timebase of the system
519530
if isinstance(sys, (int, float, complex, np.number)):
520-
# OK as long as strict checking is off
531+
# Constants OK as long as strict checking is off
521532
return True if not strict else False
522533
else:
523534
return sys.isdtime(strict)
524535

525536

526537
# Check to see if a system is a continuous time system
527-
def isctime(sys, strict=False):
538+
def isctime(sys=None, dt=None, strict=False):
528539
"""
529540
Check to see if a system is a continuous-time system.
530541
531542
Parameters
532543
----------
533-
sys : I/O or LTI system
534-
System to be checked
544+
sys : I/O system, optional
545+
System to be checked.
546+
dt : None or number, optional
547+
Timebase to be checked.
535548
strict: bool (default = False)
536-
If strict is True, make sure that timebase is not None
549+
If strict is True, make sure that timebase is not None.
537550
"""
538551

539-
# Check to see if this is a constant
552+
# See if we were passed a timebase instead of a system
553+
if sys is None:
554+
if dt is None:
555+
return True if not strict else False
556+
else:
557+
return dt == 0
558+
elif dt is not None:
559+
raise TypeError("passing both system and timebase not allowed")
560+
561+
# Check timebase of the system
540562
if isinstance(sys, (int, float, complex, np.number)):
541-
# OK as long as strict checking is off
563+
# Constants OK as long as strict checking is off
542564
return True if not strict else False
543565
else:
544566
return sys.isctime(strict)

‎control/matlab/wrappers.py‎

Lines changed: 103 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,14 @@
1212
from ..lti import LTI
1313
from ..exception import ControlArgument
1414

15-
__all__ = ['bode', 'nyquist', 'ngrid', 'dcgain', 'connect']
15+
__all__ = ['bode', 'nyquist', 'ngrid', 'rlocus', 'pzmap', 'dcgain', 'connect']
1616

1717
def bode(*args, **kwargs):
1818
"""bode(syslist[, omega, dB, Hz, deg, ...])
1919
2020
Bode plot of the frequency response.
2121
22-
Plots a bode gain and phase diagram
22+
Plots a bode gain and phase diagram.
2323
2424
Parameters
2525
----------
@@ -195,6 +195,106 @@ def _parse_freqplot_args(*args):
195195
return syslist, omega, plotstyle, other
196196

197197

198+
# TODO: rewrite to call root_locus_map, without using legacy plot keyword
199+
def rlocus(*args, **kwargs):
200+
"""rlocus(sys[, klist, xlim, ylim, ...])
201+
202+
Root locus diagram.
203+
204+
Calculate the root locus by finding the roots of 1 + k * G(s) where G
205+
is a linear system with transfer function num(s)/den(s) and each k is
206+
an element of kvect.
207+
208+
Parameters
209+
----------
210+
sys : LTI object
211+
Linear input/output systems (SISO only, for now).
212+
kvect : array_like, optional
213+
Gains to use in computing plot of closed-loop poles.
214+
xlim : tuple or list, optional
215+
Set limits of x axis, normally with tuple
216+
(see :doc:`matplotlib:api/axes_api`).
217+
ylim : tuple or list, optional
218+
Set limits of y axis, normally with tuple
219+
(see :doc:`matplotlib:api/axes_api`).
220+
221+
Returns
222+
-------
223+
roots : ndarray
224+
Closed-loop root locations, arranged in which each row corresponds
225+
to a gain in gains.
226+
gains : ndarray
227+
Gains used. Same as kvect keyword argument if provided.
228+
229+
Notes
230+
-----
231+
This function is a wrapper for :func:`~control.root_locus_plot`,
232+
with legacy return arguments.
233+
234+
"""
235+
from ..rlocus import root_locus_plot
236+
237+
# Use the plot keyword to get legacy behavior
238+
kwargs = dict(kwargs) # make a copy since we modify this
239+
if 'plot' not in kwargs:
240+
kwargs['plot'] = True
241+
242+
# Turn off deprecation warning
243+
with warnings.catch_warnings():
244+
warnings.filterwarnings(
245+
'ignore', message='.* return values of .* is deprecated',
246+
category=DeprecationWarning)
247+
retval = root_locus_plot(*args, **kwargs)
248+
249+
return retval
250+
251+
252+
# TODO: rewrite to call pole_zero_map, without using legacy plot keyword
253+
def pzmap(*args, **kwargs):
254+
"""pzmap(sys[, grid, plot])
255+
256+
Plot a pole/zero map for a linear system.
257+
258+
Parameters
259+
----------
260+
sys: LTI (StateSpace or TransferFunction)
261+
Linear system for which poles and zeros are computed.
262+
plot: bool, optional
263+
If ``True`` a graph is generated with Matplotlib,
264+
otherwise the poles and zeros are only computed and returned.
265+
grid: boolean (default = False)
266+
If True plot omega-damping grid.
267+
268+
Returns
269+
-------
270+
poles: array
271+
The system's poles.
272+
zeros: array
273+
The system's zeros.
274+
275+
Notes
276+
-----
277+
This function is a wrapper for :func:`~control.pole_zero_plot`,
278+
with legacy return arguments.
279+
280+
"""
281+
from ..pzmap import pole_zero_plot
282+
283+
# Use the plot keyword to get legacy behavior
284+
kwargs = dict(kwargs) # make a copy since we modify this
285+
if 'plot' not in kwargs:
286+
kwargs['plot'] = True
287+
288+
# Turn off deprecation warning
289+
with warnings.catch_warnings():
290+
warnings.filterwarnings(
291+
'ignore', message='.* return values of .* is deprecated',
292+
category=DeprecationWarning)
293+
retval = pole_zero_plot(*args, **kwargs)
294+
295+
return retval
296+
297+
198298
from ..nichols import nichols_grid
199299
def ngrid():
200300
return nichols_grid()
@@ -254,6 +354,7 @@ def dcgain(*args):
254354

255355
from ..bdalg import connect as ct_connect
256356
def connect(*args):
357+
257358
"""Index-based interconnection of an LTI system.
258359
259360
The system `sys` is a system typically constructed with `append`, with

0 commit comments

Comments
 (0)