changes to rlocus to be compatible with discrete-time systems (#410) · python-control/python-control@d3142ff
@@ -43,6 +43,9 @@
4343# RMM, 2 April 2011: modified to work with new LTI structure (see ChangeLog)
4444# * Not tested: should still work on signal.ltisys objects
4545#
46+# Sawyer B. Fuller (minster@uw.edu) 21 May 2020:
47+# * added compatibility with discrete-time systems.
48+#
4649# $Id$
47504851# Packages used by this module
@@ -52,9 +55,11 @@
5255import matplotlib.pyplot as plt
5356from numpy import array, poly1d, row_stack, zeros_like, real, imag
5457import scipy.signal # signal processing toolbox
58+from .lti import isdtime
5559from .xferfcn import _convert_to_transfer_function
5660from .exception import ControlMIMONotImplemented
5761from .sisotool import _SisotoolUpdate
62+from .grid import sgrid, zgrid
5863from . import config
59646065__all__ = ['root_locus', 'rlocus']
@@ -131,6 +136,13 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
131136# Convert numerator and denominator to polynomials if they aren't
132137 (nump, denp) = _systopoly1d(sys)
133138139+# if discrete-time system and if xlim and ylim are not given,
140+# that we a view of the unit circle
141+if xlim is None and isdtime(sys, strict=True):
142+xlim = (-1.2, 1.2)
143+if ylim is None and isdtime(sys, strict=True):
144+xlim = (-1.3, 1.3)
145+134146if kvect is None:
135147start_mat = _RLFindRoots(nump, denp, [1])
136148kvect, mymat, xlim, ylim = _default_gains(nump, denp, xlim, ylim)
@@ -163,10 +175,14 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
163175 [root.real for root in start_mat],
164176 [root.imag for root in start_mat],
165177'm.', marker='s', markersize=8, zorder=20, label='gain_point')
178+s = start_mat[0][0]
179+if isdtime(sys, strict=True):
180+zeta = -np.cos(np.angle(np.log(s)))
181+else:
182+zeta = -1 * s.real / abs(s)
166183fig.suptitle(
167184"Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" %
168- (start_mat[0][0].real, start_mat[0][0].imag,
169-1, -1 * start_mat[0][0].real / abs(start_mat[0][0])),
185+ (s.real, s.imag, 1, zeta),
170186fontsize=12 if int(mpl.__version__[0]) == 1 else 10)
171187fig.canvas.mpl_connect(
172188'button_release_event',
@@ -199,20 +215,31 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None,
199215ax.plot(real(col), imag(col), plotstr, label='rootlocus')
200216201217# Set up plot axes and labels
202-if xlim:
203-ax.set_xlim(xlim)
204-if ylim:
205-ax.set_ylim(ylim)
206-207218ax.set_xlabel('Real')
208219ax.set_ylabel('Imaginary')
220+209221if grid and sisotool:
210-_sgrid_func(f)
222+if isdtime(sys, strict=True):
223+zgrid(ax=ax)
224+else:
225+_sgrid_func(f)
211226elif grid:
212-_sgrid_func()
227+if isdtime(sys, strict=True):
228+zgrid(ax=ax)
229+else:
230+_sgrid_func()
213231else:
214232ax.axhline(0., linestyle=':', color='k', zorder=-20)
215-ax.axvline(0., linestyle=':', color='k')
233+ax.axvline(0., linestyle=':', color='k', zorder=-20)
234+if isdtime(sys, strict=True):
235+ax.add_patch(plt.Circle((0,0), radius=1.0,
236+linestyle=':', edgecolor='k', linewidth=1.5,
237+fill=False, zorder=-20))
238+239+if xlim:
240+ax.set_xlim(xlim)
241+if ylim:
242+ax.set_ylim(ylim)
216243217244return mymat, kvect
218245@@ -567,12 +594,17 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
567594if abs(K.real) > 1e-8 and abs(K.imag / K.real) < gain_tolerance and \
568595event.inaxes == ax_rlocus.axes and K.real > 0.:
569596597+if isdtime(sys, strict=True):
598+zeta = -np.cos(np.angle(np.log(s)))
599+else:
600+zeta = -1 * s.real / abs(s)
601+570602# Display the parameters in the output window and figure
571603print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
572- (s.real, s.imag, K.real, -1 * s.real / abs(s)))
604+ (s.real, s.imag, K.real, zeta))
573605fig.suptitle(
574606"Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" %
575- (s.real, s.imag, K.real, -1 * s.real / abs(s)),
607+ (s.real, s.imag, K.real, zeta),
576608fontsize=12 if int(mpl.__version__[0]) == 1 else 10)
577609578610# Remove the previous line
@@ -616,13 +648,13 @@ def _sgrid_func(fig=None, zeta=None, wn=None):
616648if zeta is None:
617649zeta = _default_zetas(xlim, ylim)
618650619-angules = []
651+angles = []
620652for z in zeta:
621653if (z >= 1e-4) and (z <= 1):
622-angules.append(np.pi/2 + np.arcsin(z))
654+angles.append(np.pi/2 + np.arcsin(z))
623655else:
624656zeta.remove(z)
625-y_over_x = np.tan(angules)
657+y_over_x = np.tan(angles)
626658627659# zeta-constant lines
628660@@ -647,30 +679,30 @@ def _sgrid_func(fig=None, zeta=None, wn=None):
647679ax.plot([0, 0], [ylim[0], ylim[1]],
648680color='gray', linestyle='dashed', linewidth=0.5)
649681650-angules = np.linspace(-90, 90, 20)*np.pi/180
682+angles = np.linspace(-90, 90, 20)*np.pi/180
651683if wn is None:
652684wn = _default_wn(xlocator(), ylim)
653685654686for om in wn:
655687if om < 0:
656-yp = np.sin(angules)*np.abs(om)
657-xp = -np.cos(angules)*np.abs(om)
688+yp = np.sin(angles)*np.abs(om)
689+xp = -np.cos(angles)*np.abs(om)
658690ax.plot(xp, yp, color='gray',
659691linestyle='dashed', linewidth=0.5)
660692an = "%.2f" % -om
661693ax.annotate(an, textcoords='data', xy=[om, 0], fontsize=8)
662694663695664696def _default_zetas(xlim, ylim):
665-"""Return default list of dumps coefficients"""
697+"""Return default list of damping coefficients"""
666698sep1 = -xlim[0]/4
667699ang1 = [np.arctan((sep1*i)/ylim[1]) for i in np.arange(1, 4, 1)]
668700sep2 = ylim[1] / 3
669701ang2 = [np.arctan(-xlim[0]/(ylim[1]-sep2*i)) for i in np.arange(1, 3, 1)]
670702671-angules = np.concatenate((ang1, ang2))
672-angules = np.insert(angules, len(angules), np.pi/2)
673-zeta = np.sin(angules)
703+angles = np.concatenate((ang1, ang2))
704+angles = np.insert(angles, len(angles), np.pi/2)
705+zeta = np.sin(angles)
674706return zeta.tolist()
675707676708