Skip to content

Commit 6e33598

Browse files
murrayrmmurrayrm
committed
update xlim, ylim handling in pzmap + other supporting changes
1 parent 1390e37 commit 6e33598

File tree

5 files changed

+180
-51
lines changed

5 files changed

+180
-51
lines changed

‎control/grid.py

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,8 @@
88
from matplotlib.projections import PolarAxes
99
from matplotlib.transforms import Affine2D
1010

11+
from .iosys import isdtime
12+
1113

1214
class FormatterDMS(object):
1315
'''Transforms angle ticks to damping ratios'''
@@ -142,17 +144,22 @@ def sgrid():
142144
def _final_setup(ax):
143145
ax.set_xlabel('Real')
144146
ax.set_ylabel('Imaginary')
145-
ax.axhline(y=0, color='black', lw=1)
146-
ax.axvline(x=0, color='black', lw=1)
147+
ax.axhline(y=0, color='black', lw=0.5)
148+
ax.axvline(x=0, color='black', lw=0.5)
147149
plt.axis('equal')
148150

149151

150-
def nogrid():
151-
f = plt.gcf()
152+
def nogrid(dt=None):
153+
fig = plt.gcf()
152154
ax = plt.axes()
153155

156+
# Draw the unit circle for discrete time systems
157+
if isdtime(dt=dt, strict=True):
158+
s = np.linspace(0, 2*pi, 100)
159+
ax.plot(np.cos(s), np.sin(s), 'k--', lw=0.5, dashes=(5, 5))
160+
154161
_final_setup(ax)
155-
return ax, f
162+
return ax, fig
156163

157164

158165
def zgrid(zetas=None, wns=None, ax=None):

‎control/iosys.py

Lines changed: 34 additions & 12 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, dt=None, strict=False):
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
512+
sys : I/O system, optional
513+
System to be checked.
514+
dt : None or number, optional
515+
Timebase to be checked.
514516
strict: bool (default = False)
515-
If strict is True, make sure that timebase is not None
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/pzmap.py

Lines changed: 52 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,9 @@
44
# Date: 7 Sep 2009
55
#
66
# This file contains functions that compute poles, zeros and related
7-
# quantities for a linear system.
7+
# quantities for a linear system, as well as the main functions for
8+
# storing and plotting pole/zero and root locus diagrams. (The actual
9+
# computation of root locus diagrams is in rlocus.py.)
810
#
911

1012
import numpy as np
@@ -41,14 +43,11 @@ def plot(self, *args, **kwargs):
4143

4244
class RootLocusData:
4345
def __init__(
44-
self, poles, zeros, gains=None, loci=None, xlim=None, ylim=None,
45-
dt=None, sysname=None):
46+
self, poles, zeros, gains=None, loci=None, dt=None, sysname=None):
4647
self.poles = poles
4748
self.zeros = zeros
4849
self.gains = gains
4950
self.loci = loci
50-
self.xlim = xlim
51-
self.ylim = ylim
5251
self.dt = dt
5352
self.sysname = sysname
5453

@@ -78,7 +77,8 @@ def pole_zero_map(sysdata):
7877

7978

8079
# Root locus map
81-
def root_locus_map(sysdata, gains=None, xlim=None, ylim=None):
80+
# TODO: use rlocus.py computation instead
81+
def root_locus_map(sysdata, gains=None):
8282
# Convert the first argument to a list
8383
syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]
8484

@@ -94,22 +94,16 @@ def root_locus_map(sysdata, gains=None, xlim=None, ylim=None):
9494
# Convert numerator and denominator to polynomials if they aren't
9595
nump, denp = _systopoly1d(sys[0, 0])
9696

97-
if xlim is None and sys.isdtime(strict=True):
98-
xlim = (-1.2, 1.2)
99-
if ylim is None and sys.isdtime(strict=True):
100-
xlim = (-1.3, 1.3)
101-
10297
if gains is None:
103-
kvect, root_array, xlim, ylim = _default_gains(
104-
nump, denp, xlim, ylim)
98+
kvect, root_array, _, _ = _default_gains(nump, denp, None, None)
10599
else:
106100
kvect = np.atleast_1d(gains)
107101
root_array = _RLFindRoots(nump, denp, kvect)
108102
root_array = _RLSortRoots(root_array)
109103

110104
responses.append(RootLocusData(
111105
sys.poles(), sys.zeros(), kvect, root_array,
112-
dt=sys.dt, sysname=sys.name, xlim=xlim, ylim=ylim))
106+
dt=sys.dt, sysname=sys.name))
113107

114108
if isinstance(sysdata, (list, tuple)):
115109
return RootLocusList(responses)
@@ -203,7 +197,7 @@ def pole_zero_plot(
203197
return poles, zeros
204198

205199
# Initialize the figure
206-
# TODO: turn into standard utility function
200+
# TODO: turn into standard utility function (from plotutil.py?)
207201
fig = plt.gcf()
208202
axs = fig.get_axes()
209203
if len(axs) > 1:
@@ -213,21 +207,22 @@ def pole_zero_plot(
213207
with plt.rc_context(freqplot_rcParams):
214208
if grid:
215209
plt.clf()
216-
if all([response.dt in [0, None] for response in data]):
210+
if all([isctime(dt=response.dt) for response in data]):
217211
ax, fig = sgrid()
218-
elif all([response.dt > 0 for response in data]):
212+
elif all([isdtime(dt=response.dt) for response in data]):
219213
ax, fig = zgrid()
220214
else:
221215
ValueError(
222216
"incompatible time responses; don't know how to grid")
223217
elif len(axs) == 0:
224-
ax, fig = nogrid()
218+
ax, fig = nogrid(data[0].dt) # use first response timebase
225219
else:
226-
# Use the existing axes
220+
# Use the existing axes and any grid that is there
221+
# TODO: allow axis to be overriden via parameter
227222
ax = axs[0]
228223

229-
# Handle color cycle manually as all singular values
230-
# of the same systems are expected to be of the same color
224+
# Handle color cycle manually as all root locus segments
225+
# of the same system are expected to be of the same color
231226
color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
232227
color_offset = 0
233228
if len(ax.lines) > 0:
@@ -240,6 +235,7 @@ def pole_zero_plot(
240235
for i, j in itertools.product(range(out.shape[0]), range(out.shape[1])):
241236
out[i, j] = [] # unique list in each element
242237

238+
# Plot the responses (and keep track of axes limits)
243239
xlim, ylim = ax.get_xlim(), ax.get_ylim()
244240
for idx, response in enumerate(pzmap_responses):
245241
poles = response.poles
@@ -272,13 +268,14 @@ def pole_zero_plot(
272268
real(locus), imag(locus), color=color,
273269
label=response.sysname)
274270

275-
# Compute the axis limits to use
276-
if response.xlim is not None:
277-
xlim = (min(xlim[0], response.xlim[0]),
278-
max(xlim[1], response.xlim[1]))
279-
if response.ylim is not None:
280-
ylim = (min(ylim[0], response.ylim[0]),
281-
max(ylim[1], response.ylim[1]))
271+
# Compute the axis limits to use based on the response
272+
resp_xlim, resp_ylim = _compute_root_locus_limits(response.loci)
273+
274+
# Keep track of the current limits
275+
xlim = [min(xlim[0], resp_xlim[0]), max(xlim[1], resp_xlim[1])]
276+
ylim = [min(ylim[0], resp_ylim[0]), max(ylim[1], resp_ylim[1])]
277+
278+
# TODO: add arrows to root loci (reuse Nyquist arrow code?)
282279

283280
# Set up the limits for the plot
284281
ax.set_xlim(xlim if xlim_user is None else xlim_user)
@@ -329,4 +326,31 @@ def pole_zero_plot(
329326
return out
330327

331328

329+
# Utility function to compute limits for root loci
330+
def _compute_root_locus_limits(loci):
331+
# Go through each locus
332+
xlim, ylim = [0, 0], 0
333+
for locus in loci.transpose():
334+
# Include all starting points
335+
xlim = [min(xlim[0], locus[0].real), max(xlim[1], locus[0].real)]
336+
ylim = max(ylim, locus[0].imag)
337+
338+
# Find the local maxima of root locus curve
339+
xpeaks = np.where(
340+
np.diff(np.abs(locus.real)) < 0, locus.real[0:-1], 0)
341+
xlim = [min(xlim[0], np.min(xpeaks)), max(xlim[1], np.max(xpeaks))]
342+
343+
ypeaks = np.where(
344+
np.diff(np.abs(locus.imag)) < 0, locus.imag[0:-1], 0)
345+
ylim = max(ylim, np.max(ypeaks))
346+
347+
# Adjust the limits to include some space around features
348+
rho = 1.5
349+
xlim[0] = rho * xlim[0] if xlim[0] < 0 else 0
350+
xlim[1] = rho * xlim[1] if xlim[1] > 0 else 0
351+
ylim = rho * ylim if ylim > 0 else np.max(np.abs(xlim))
352+
353+
return xlim, [-ylim, ylim]
354+
355+
332356
pzmap = pole_zero_plot

‎control/rlocus.py

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@
3232
__all__ = ['root_locus', 'rlocus']
3333

3434
# Default values for module parameters
35+
# TODO: merge these with pzmap parameters (?)
3536
_rlocus_defaults = {
3637
'rlocus.grid': True,
3738
'rlocus.plotstr': 'b' if int(mpl.__version__[0]) == 1 else 'C0',
@@ -41,15 +42,16 @@
4142

4243

4344
# Main function: compute a root locus diagram
45+
# TODO: update to use pzmap data structures and plotting
4446
def root_locus(sys, kvect=None, xlim=None, ylim=None,
4547
plotstr=None, plot=True, print_gain=None, grid=None, ax=None,
4648
initial_gain=None, **kwargs):
4749

4850
"""Root locus plot.
4951
50-
Calculate the root locus by finding the roots of 1+k*TF(s)
51-
where TF is self.num(s)/self.den(s) and each k is an element
52-
of kvect.
52+
Calculate the root locus by finding the roots of 1 + k * G(s) where G
53+
is a linear system with transfer function num(s)/den(s) and each k is
54+
an element of kvect.
5355
5456
Parameters
5557
----------
@@ -127,6 +129,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
127129
raise TypeError("unrecognized keywords: ", str(kwargs))
128130

129131
# Create the Plot
132+
# TODO: replace with pole_zero_plot and move additional functionality there
130133
if plot:
131134
if ax is None:
132135
ax = plt.gca()
@@ -139,6 +142,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
139142
else:
140143
start_roots = None
141144

145+
# TODO: don't rely on `start_roots` (sisotool holdover)
142146
if print_gain and start_roots is None:
143147
fig.canvas.mpl_connect(
144148
'button_release_event',
@@ -148,7 +152,8 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
148152
ax.plot(
149153
[root.real for root in start_roots],
150154
[root.imag for root in start_roots],
151-
marker='s', markersize=6, zorder=20, color='k', label='gain_point')
155+
marker='s', markersize=6, zorder=20, color='k',
156+
label='gain_point')
152157
s = start_roots[0][0]
153158
if isdtime(sys, strict=True):
154159
zeta = -np.cos(np.angle(np.log(s)))
@@ -219,16 +224,23 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
219224
Saddle River, NJ : New Delhi: Prentice Hall..
220225
221226
"""
227+
# Compute the break points on the real axis for the root locus plot
222228
k_break, real_break = _break_points(num, den)
229+
230+
# Decide on the maximum gain to use and create the gain vector
223231
kmax = _k_max(num, den, real_break, k_break)
224232
kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break)))
225233
kvect.sort()
226234

235+
# Find the roots for all of the gains and sort them
227236
root_array = _RLFindRoots(num, den, kvect)
228237
root_array = _RLSortRoots(root_array)
238+
239+
# Keep track of the open loop poles and zeros
229240
open_loop_poles = den.roots
230241
open_loop_zeros = num.roots
231242

243+
# ???
232244
if open_loop_zeros.size != 0 and \
233245
open_loop_zeros.size < open_loop_poles.size:
234246
open_loop_zeros_xl = np.append(
@@ -283,7 +295,8 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
283295
tolerance = x_tolerance
284296
else:
285297
tolerance = np.min([x_tolerance, y_tolerance])
286-
indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
298+
indexes_too_far = _indexes_filt(
299+
root_array, tolerance, zoom_xlim, zoom_ylim)
287300

288301
# Add more points into the root locus for points that are too far apart
289302
while len(indexes_too_far) > 0 and kvect.size < 5000:
@@ -295,7 +308,8 @@ def _default_gains(num, den, xlim, ylim, zoom_xlim=None, zoom_ylim=None):
295308
root_array = np.insert(root_array, index + 1, new_points, axis=0)
296309

297310
root_array = _RLSortRoots(root_array)
298-
indexes_too_far = _indexes_filt(root_array, tolerance, zoom_xlim, zoom_ylim)
311+
indexes_too_far = _indexes_filt(
312+
root_array, tolerance, zoom_xlim, zoom_ylim)
299313

300314
new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4)))
301315
new_points = _RLFindRoots(num, den, new_gains[1:4])
@@ -601,6 +615,7 @@ def _removeLine(label, ax):
601615
del line
602616

603617

618+
# TODO: remove and replace with sgid()?
604619
def _sgrid_func(ax, zeta=None, wn=None):
605620
# Get locator function for x-axis, y-axis tick marks
606621
xlocator = ax.get_xaxis().get_major_locator()

0 commit comments

Comments
 (0)