|
37 | 37 | # * Added BSD copyright info to file (per Ryan) |
38 | 38 | # * Added code to convert (num, den) to poly1d's if they aren't already. |
39 | 39 | # This allows Ryan's code to run on a standard signal.ltisys object |
40 | | -# or a controls.TransferFunction object. |
| 40 | +# or a control.TransferFunction object. |
41 | 41 | # * Added some comments to make sure I understand the code |
| 42 | +# |
| 43 | +# RMM, 2 April 2011: modified to work with new Lti structure (see ChangeLog) |
| 44 | +# * Not tested: should still work on signal.ltisys objects |
42 | 45 | # |
43 | 46 | # $Id$ |
44 | 47 |
|
45 | 48 | # Packages used by this module |
46 | 49 | from scipy import * |
| 50 | +import scipy.signal # signal processing toolbox |
| 51 | +import pylab # plotting routines |
| 52 | +import xferfcn # transfer function manipulation |
47 | 53 |
|
48 | 54 | # Main function: compute a root locus diagram |
49 | | -def RootLocus(sys, kvect, fig=None, fignum=1, \ |
50 | | - clear=True, xlim=None, ylim=None, plotstr='-'): |
| 55 | +def RootLocus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True): |
51 | 56 | """Calculate the root locus by finding the roots of 1+k*TF(s) |
52 | 57 | where TF is self.num(s)/self.den(s) and each k is an element |
53 | | - of kvect.""" |
| 58 | + of kvect. |
| 59 | +
|
| 60 | + Parameters |
| 61 | + ---------- |
| 62 | + sys : linsys |
| 63 | + Linear input/output systems (SISO only, for now) |
| 64 | + klist : gain_range (default = None) |
| 65 | + List of gains to use in computing diagram |
| 66 | + Plot : boolean (default = True) |
| 67 | + If True, plot magnitude and phase |
| 68 | +
|
| 69 | + Return values |
| 70 | + ------------- |
| 71 | + rlist : list of computed root locations |
| 72 | + klist : list of gains |
| 73 | + """ |
54 | 74 |
|
55 | 75 | # Convert numerator and denominator to polynomials if they aren't |
56 | 76 | (nump, denp) = _systopoly1d(sys); |
57 | 77 |
|
58 | | - # Set up the figure |
59 | | - if fig is None: |
60 | | - import pylab |
61 | | - fig = pylab.figure(fignum) |
62 | | - if clear: |
63 | | - fig.clf() |
64 | | - ax = fig.add_subplot(111) |
65 | | - |
66 | 78 | # Compute out the loci |
67 | 79 | mymat = _RLFindRoots(sys, kvect) |
68 | 80 | mymat = _RLSortRoots(sys, mymat) |
69 | 81 |
|
70 | | - # plot open loop poles |
71 | | - poles = array(denp.r) |
72 | | - ax.plot(real(poles), imag(poles), 'x') |
73 | | - |
74 | | - # plot open loop zeros |
75 | | - zeros = array(nump.r) |
76 | | - if zeros.any(): |
77 | | - ax.plot(real(zeros), imag(zeros), 'o') |
78 | | - |
79 | | - # Now plot the loci |
80 | | - for col in mymat.T: |
81 | | - ax.plot(real(col), imag(col), plotstr) |
82 | | - |
83 | | - # Set up plot axes and labels |
84 | | - if xlim: |
85 | | - ax.set_xlim(xlim) |
86 | | - if ylim: |
87 | | - ax.set_ylim(ylim) |
88 | | - ax.set_xlabel('Real') |
89 | | - ax.set_ylabel('Imaginary') |
| 82 | + # Create the plot |
| 83 | + if (Plot): |
| 84 | + ax = pylab.axes(); |
| 85 | + |
| 86 | + # plot open loop poles |
| 87 | + poles = array(denp.r) |
| 88 | + ax.plot(real(poles), imag(poles), 'x') |
| 89 | + |
| 90 | + # plot open loop zeros |
| 91 | + zeros = array(nump.r) |
| 92 | + if zeros.any(): |
| 93 | + ax.plot(real(zeros), imag(zeros), 'o') |
| 94 | + |
| 95 | + # Now plot the loci |
| 96 | + for col in mymat.T: |
| 97 | + ax.plot(real(col), imag(col), plotstr) |
| 98 | + |
| 99 | + # Set up plot axes and labels |
| 100 | + if xlim: |
| 101 | + ax.set_xlim(xlim) |
| 102 | + if ylim: |
| 103 | + ax.set_ylim(ylim) |
| 104 | + ax.set_xlabel('Real') |
| 105 | + ax.set_ylabel('Imaginary') |
| 106 | + |
90 | 107 | return mymat |
91 | 108 |
|
92 | 109 | # Utility function to extract numerator and denominator polynomials |
93 | 110 | def _systopoly1d(sys): |
94 | 111 | """Extract numerator and denominator polynomails for a system""" |
| 112 | + # Allow inputs from the signal processing toolbox |
| 113 | + if (isinstance(sys, scipy.signal.lti)): |
| 114 | + nump = sys.num; denp = sys.den; |
| 115 | + |
| 116 | + else: |
| 117 | + # Convert to a transfer function, if needed |
| 118 | + sys = xferfcn._convertToTransferFunction(sys) |
| 119 | + |
| 120 | + # Make sure we have a SISO system |
| 121 | + if (sys.inputs > 1 or sys.outputs > 1): |
| 122 | + raise ControlMIMONotImplemented() |
95 | 123 |
|
96 | | - # Start by extracting the numerator and denominator from system object |
97 | | - nump = sys.num; denp = sys.den; |
| 124 | + # Start by extracting the numerator and denominator from system object |
| 125 | + nump = sys.num[0][0]; denp = sys.den[0][0]; |
98 | 126 |
|
99 | 127 | # Check to see if num, den are already polynomials; otherwise convert |
100 | 128 | if (not isinstance(nump, poly1d)): nump = poly1d(nump) |
|
0 commit comments