|
3431
|
1 ## Copyright (C) 1996 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by the |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
18 ## 02110-1301 USA. |
|
3431
|
19 |
|
|
20 ## -*- texinfo -*- |
|
5016
|
21 ## @deftypefn {Function File} {[@var{rldata}, @var{k}] =} rlocus (@var{sys}[, @var{increment}, @var{min_k}, @var{max_k}]) |
|
3431
|
22 ## |
|
5016
|
23 ## Displays root locus plot of the specified @acronym{SISO} system. |
|
|
24 ## @example |
|
|
25 ## @group |
|
3431
|
26 ## ----- --- -------- |
|
|
27 ## --->| + |---|k|---->| SISO |-----------> |
|
|
28 ## ----- --- -------- | |
|
|
29 ## - ^ | |
|
|
30 ## |_____________________________| |
|
5016
|
31 ## @end group |
|
|
32 ## @end example |
|
3431
|
33 ## |
|
5016
|
34 ## @strong{Inputs} |
|
|
35 ## @table @var |
|
|
36 ## @item sys |
|
|
37 ## system data structure |
|
|
38 ## @item min_k |
|
|
39 ## Minimum value of @var{k} |
|
|
40 ## @item max_k |
|
|
41 ## Maximum value of @var{k} |
|
|
42 ## @item increment |
|
|
43 ## The increment used in computing gain values |
|
|
44 ## @end table |
|
|
45 ## |
|
|
46 ## @strong{Outputs} |
|
|
47 ## |
|
|
48 ## Plots the root locus to the screen. |
|
|
49 ## @table @var |
|
|
50 ## @item rldata |
|
|
51 ## Data points plotted: in column 1 real values, in column 2 the imaginary values. |
|
|
52 ## @item k |
|
|
53 ## Gains for real axis break points. |
|
|
54 ## @end table |
|
3431
|
55 ## @end deftypefn |
|
|
56 |
|
|
57 ## Author: David Clem |
|
|
58 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
59 ## Updated by Kristi McGowan July 1996 for intelligent gain selection |
|
|
60 ## Updated by John Ingram July 1996 for systems |
|
|
61 |
|
|
62 function [rldata, k_break, rlpol, gvec, real_ax_pts] = rlocus (sys, increment, min_k, max_k) |
|
|
63 |
|
|
64 if (nargin < 1) | (nargin > 4) |
|
|
65 usage("rlocus(sys[,inc,mink,maxk])"); |
|
|
66 endif |
|
|
67 |
|
|
68 ## Convert the input to a transfer function if necessary |
|
|
69 |
|
|
70 [num,den] = sys2tf(sys) # extract numerator/denom polyomials |
|
|
71 lnum = length(num); lden = length(den); |
|
|
72 if(lden < 2) |
|
|
73 error(sprintf("length of derivative=%d, doesn't make sense",lden)); |
|
|
74 elseif(lnum == 1) |
|
|
75 num = [0, num]; # so that derivative is shortened by one |
|
|
76 endif |
|
|
77 |
|
|
78 ## root locus plot axis limits |
|
|
79 |
|
|
80 ## compute real axis locus breakpoints |
|
|
81 ## compute the derivative of the numerator and the denominator |
|
|
82 dern=polyderiv(num); derd=polyderiv(den); |
|
|
83 |
|
|
84 ## compute real axis breakpoints |
|
|
85 real_ax_pol = conv(den,dern) - conv(num,derd); |
|
|
86 real_ax_pts = roots(real_ax_pol); |
|
|
87 if(isempty(real_ax_pts)) |
|
|
88 k_break = []; |
|
|
89 maxk = 0; |
|
|
90 else |
|
|
91 ## compute gains that achieve the breakpoints |
|
|
92 c1 = polyval(num,real_ax_pts); |
|
|
93 c2 = polyval(den,real_ax_pts); |
|
|
94 k_break = -real(c2 ./ c1); |
|
|
95 maxk = max(max(k_break,0)); |
|
|
96 endif |
|
|
97 |
|
|
98 ## compute gain ranges based on computed K values |
|
|
99 if(maxk == 0) maxk = 1; |
|
|
100 else maxk = 1.1*maxk; endif |
|
|
101 mink = 0; |
|
|
102 ngain = 20; |
|
|
103 |
|
|
104 ## check for input arguments: |
|
|
105 if (nargin > 2) mink = min_k; endif |
|
|
106 if (nargin > 3) maxk = max_k; endif |
|
|
107 if (nargin > 1) |
|
|
108 if(increment <= 0) error("increment must be positive"); |
|
|
109 else |
|
|
110 ngain = (maxk-mink)/increment; |
|
|
111 endif |
|
|
112 endif |
|
|
113 |
|
|
114 ## vector of gains |
|
|
115 ngain = max(3,ngain); |
|
|
116 gvec = linspace(mink,maxk,ngain); |
|
|
117 |
|
|
118 ## Find the open loop zeros and the initial poles |
|
|
119 rlzer = roots(num); |
|
|
120 |
|
|
121 ## update num to be the same length as den |
|
|
122 lnum = length(num); if(lnum < lden) num = [zeros(1,lden - lnum),num]; endif |
|
|
123 |
|
|
124 ## compute preliminary pole sets |
|
|
125 nroots = lden-1; |
|
|
126 for ii=1:ngain |
|
|
127 gain = gvec(ii); |
|
|
128 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
129 endfor |
|
|
130 |
|
|
131 ## compute axis limits (isolate asymptotes) |
|
|
132 olpol = roots(den); |
|
|
133 real_axdat = union(real(rlzer), real(union(olpol,real_ax_pts)) ); |
|
|
134 rmin = min(real_axdat); rmax = max(real_axdat); |
|
|
135 |
|
|
136 rlpolv = [vec(rlpol); vec(real_axdat)]; |
|
|
137 idx = find(real(rlpolv) >= rmin & real(rlpolv) <= rmax); |
|
|
138 axlim = axis2dlim([real(rlpolv(idx)),imag(rlpolv(idx))]); |
|
|
139 xmin = axlim(1); |
|
|
140 xmax = axlim(2); |
|
|
141 |
|
|
142 ## set smoothing tolerance per axis limits |
|
|
143 smtol = 0.01*max(abs(axlim)); |
|
|
144 |
|
|
145 ## smooth poles if necessary, up to maximum of 1000 gain points |
|
|
146 ## only smooth points within the axis limit window |
|
|
147 ## smoothing done if max_k not specified as a command argument |
|
|
148 done=(nargin == 4); # perform a smoothness check |
|
|
149 while((!done) & ngain < 1000) |
|
|
150 done = 1 ; # assume done |
|
|
151 dp = abs(diff(rlpol'))'; |
|
|
152 maxd = max(dp); |
|
|
153 ## search for poles in the real axis limits whose neighbors are distant |
|
|
154 idx = find(maxd > smtol); |
|
|
155 for ii=1:length(idx) |
|
|
156 i1 = idx(ii); g1 = gvec(i1); p1 = rlpol(:,i1); |
|
|
157 i2 = idx(ii)+1; g2 = gvec(i2); p2 = rlpol(:,i2); |
|
|
158 |
|
|
159 ## isolate poles in p1, p2 that are inside the real axis limits |
|
|
160 bidx = find( (real(p1) >= xmin & real(p1) <= xmax) ... |
|
|
161 | (real(p2) >= xmin & real(p2) <= xmax) ); |
|
|
162 if(!isempty(bidx)) |
|
|
163 p1 = p1(bidx); |
|
|
164 p2 = p2(bidx); |
|
|
165 if( max(abs(p2-p1)) > smtol) |
|
|
166 newg = linspace(g1,g2,5); |
|
|
167 newg = newg(2:4); |
|
|
168 if(isempty(newg)) |
|
|
169 printf("rlocus: empty newg") |
|
|
170 g1 |
|
|
171 g2 |
|
|
172 i1 |
|
|
173 i2 |
|
|
174 idx_i1 = idx(ii) |
|
|
175 gvec_i1 = gvec(i1:i2) |
|
|
176 delta_vec_i1 = diff(gvec(i1:i2)) |
|
|
177 prompt |
|
|
178 endif |
|
|
179 gvec = [gvec,newg]; |
|
|
180 done = 0; # need to process new gains |
|
|
181 endif |
|
|
182 endif |
|
|
183 endfor |
|
|
184 |
|
|
185 ## process new gain values |
|
|
186 ngain1 = length(gvec); |
|
|
187 for ii=(ngain+1):ngain1 |
|
|
188 gain = gvec(ii); |
|
|
189 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
190 endfor |
|
|
191 |
|
|
192 [gvec,idx] = sort(gvec); |
|
|
193 rlpol = rlpol(:,idx); |
|
|
194 ngain = length(gvec); |
|
|
195 endwhile |
|
|
196 |
|
|
197 ## Plot the data |
|
|
198 if(nargout == 0) |
|
|
199 rlpolv = vec(rlpol); |
|
|
200 idx = find(real(rlpolv) >= xmin & real(rlpolv) <= xmax); |
|
|
201 axdata = [real(rlpolv(idx)),imag(rlpolv(idx))]; |
|
|
202 axlim = axis2dlim(axdata); |
|
|
203 axlim(1:2) = [xmin, xmax]; |
|
5215
|
204 __gnuplot_set__ nologscale xy; |
|
3431
|
205 grid("on"); |
|
|
206 rldata = [real(rlpolv), imag(rlpolv) ]; |
|
|
207 axis(axlim); |
|
|
208 [stn,inname,outname] = sysgetsignals(sys); |
|
|
209 xlabel(sprintf("Root locus from %s to %s, gain=[%f,%f]: Real axis", ... |
|
4771
|
210 inname{1}, outname{1},gvec(1),gvec(ngain))); |
|
3431
|
211 ylabel("Imag. axis"); |
|
|
212 |
|
|
213 plot(real(rlpolv),imag(rlpolv),".1;locus points;", ... |
|
|
214 real(olpol),imag(olpol),"x2;open loop poles;", ... |
|
|
215 real(rlzer),imag(rlzer),"o3;zeros;"); |
|
|
216 rldata = []; |
|
|
217 endif |
|
|
218 endfunction |
