|
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 -*- |
|
5460
|
21 ## @deftypefn {Function File} {[@var{rldata}, @var{k_break}, @var{rlpol}, @var{gvec}, @var{real_ax_pts}] =} rlocus (@var{sys}[, @var{increment}, @var{min_k}, @var{max_k}]) |
|
3431
|
22 ## |
|
5460
|
23 ## Display root locus plot of the specified @acronym{SISO} system. |
|
5016
|
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. |
|
5460
|
52 ## @item k_break |
|
5016
|
53 ## Gains for real axis break points. |
|
5460
|
54 ## @item rlpol |
|
|
55 ## Closed-loop roots for each gain value: 1 locus branch per row; 1 pole |
|
|
56 ## set per column |
|
|
57 ## @item gvec |
|
|
58 ## Gains vector |
|
|
59 ## @item real_ax_pts |
|
|
60 ## Real axis breakpoints |
|
5016
|
61 ## @end table |
|
3431
|
62 ## @end deftypefn |
|
|
63 |
|
|
64 ## Author: David Clem |
|
|
65 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
66 ## Updated by Kristi McGowan July 1996 for intelligent gain selection |
|
|
67 ## Updated by John Ingram July 1996 for systems |
|
|
68 |
|
|
69 function [rldata, k_break, rlpol, gvec, real_ax_pts] = rlocus (sys, increment, min_k, max_k) |
|
|
70 |
|
|
71 if (nargin < 1) | (nargin > 4) |
|
6046
|
72 print_usage (); |
|
3431
|
73 endif |
|
|
74 |
|
|
75 ## Convert the input to a transfer function if necessary |
|
|
76 |
|
5605
|
77 [num,den] = sys2tf(sys); # extract numerator/denom polyomials |
|
3431
|
78 lnum = length(num); lden = length(den); |
|
|
79 if(lden < 2) |
|
|
80 error(sprintf("length of derivative=%d, doesn't make sense",lden)); |
|
|
81 elseif(lnum == 1) |
|
|
82 num = [0, num]; # so that derivative is shortened by one |
|
|
83 endif |
|
|
84 |
|
|
85 ## root locus plot axis limits |
|
|
86 |
|
|
87 ## compute real axis locus breakpoints |
|
|
88 ## compute the derivative of the numerator and the denominator |
|
|
89 dern=polyderiv(num); derd=polyderiv(den); |
|
|
90 |
|
|
91 ## compute real axis breakpoints |
|
|
92 real_ax_pol = conv(den,dern) - conv(num,derd); |
|
|
93 real_ax_pts = roots(real_ax_pol); |
|
|
94 if(isempty(real_ax_pts)) |
|
|
95 k_break = []; |
|
|
96 maxk = 0; |
|
|
97 else |
|
|
98 ## compute gains that achieve the breakpoints |
|
|
99 c1 = polyval(num,real_ax_pts); |
|
|
100 c2 = polyval(den,real_ax_pts); |
|
|
101 k_break = -real(c2 ./ c1); |
|
|
102 maxk = max(max(k_break,0)); |
|
|
103 endif |
|
|
104 |
|
|
105 ## compute gain ranges based on computed K values |
|
|
106 if(maxk == 0) maxk = 1; |
|
|
107 else maxk = 1.1*maxk; endif |
|
|
108 mink = 0; |
|
|
109 ngain = 20; |
|
|
110 |
|
|
111 ## check for input arguments: |
|
|
112 if (nargin > 2) mink = min_k; endif |
|
|
113 if (nargin > 3) maxk = max_k; endif |
|
|
114 if (nargin > 1) |
|
|
115 if(increment <= 0) error("increment must be positive"); |
|
|
116 else |
|
|
117 ngain = (maxk-mink)/increment; |
|
|
118 endif |
|
|
119 endif |
|
|
120 |
|
|
121 ## vector of gains |
|
|
122 ngain = max(3,ngain); |
|
|
123 gvec = linspace(mink,maxk,ngain); |
|
|
124 |
|
|
125 ## Find the open loop zeros and the initial poles |
|
|
126 rlzer = roots(num); |
|
|
127 |
|
|
128 ## update num to be the same length as den |
|
|
129 lnum = length(num); if(lnum < lden) num = [zeros(1,lden - lnum),num]; endif |
|
|
130 |
|
|
131 ## compute preliminary pole sets |
|
|
132 nroots = lden-1; |
|
|
133 for ii=1:ngain |
|
|
134 gain = gvec(ii); |
|
|
135 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
136 endfor |
|
|
137 |
|
|
138 ## compute axis limits (isolate asymptotes) |
|
|
139 olpol = roots(den); |
|
|
140 real_axdat = union(real(rlzer), real(union(olpol,real_ax_pts)) ); |
|
|
141 rmin = min(real_axdat); rmax = max(real_axdat); |
|
|
142 |
|
|
143 rlpolv = [vec(rlpol); vec(real_axdat)]; |
|
|
144 idx = find(real(rlpolv) >= rmin & real(rlpolv) <= rmax); |
|
|
145 axlim = axis2dlim([real(rlpolv(idx)),imag(rlpolv(idx))]); |
|
|
146 xmin = axlim(1); |
|
|
147 xmax = axlim(2); |
|
|
148 |
|
|
149 ## set smoothing tolerance per axis limits |
|
|
150 smtol = 0.01*max(abs(axlim)); |
|
|
151 |
|
|
152 ## smooth poles if necessary, up to maximum of 1000 gain points |
|
|
153 ## only smooth points within the axis limit window |
|
|
154 ## smoothing done if max_k not specified as a command argument |
|
|
155 done=(nargin == 4); # perform a smoothness check |
|
|
156 while((!done) & ngain < 1000) |
|
|
157 done = 1 ; # assume done |
|
|
158 dp = abs(diff(rlpol'))'; |
|
|
159 maxd = max(dp); |
|
|
160 ## search for poles in the real axis limits whose neighbors are distant |
|
|
161 idx = find(maxd > smtol); |
|
|
162 for ii=1:length(idx) |
|
|
163 i1 = idx(ii); g1 = gvec(i1); p1 = rlpol(:,i1); |
|
|
164 i2 = idx(ii)+1; g2 = gvec(i2); p2 = rlpol(:,i2); |
|
|
165 |
|
|
166 ## isolate poles in p1, p2 that are inside the real axis limits |
|
|
167 bidx = find( (real(p1) >= xmin & real(p1) <= xmax) ... |
|
|
168 | (real(p2) >= xmin & real(p2) <= xmax) ); |
|
|
169 if(!isempty(bidx)) |
|
|
170 p1 = p1(bidx); |
|
|
171 p2 = p2(bidx); |
|
|
172 if( max(abs(p2-p1)) > smtol) |
|
|
173 newg = linspace(g1,g2,5); |
|
|
174 newg = newg(2:4); |
|
|
175 if(isempty(newg)) |
|
|
176 printf("rlocus: empty newg") |
|
|
177 g1 |
|
|
178 g2 |
|
|
179 i1 |
|
|
180 i2 |
|
|
181 idx_i1 = idx(ii) |
|
|
182 gvec_i1 = gvec(i1:i2) |
|
|
183 delta_vec_i1 = diff(gvec(i1:i2)) |
|
|
184 prompt |
|
|
185 endif |
|
|
186 gvec = [gvec,newg]; |
|
|
187 done = 0; # need to process new gains |
|
|
188 endif |
|
|
189 endif |
|
|
190 endfor |
|
|
191 |
|
|
192 ## process new gain values |
|
|
193 ngain1 = length(gvec); |
|
|
194 for ii=(ngain+1):ngain1 |
|
|
195 gain = gvec(ii); |
|
|
196 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
197 endfor |
|
|
198 |
|
|
199 [gvec,idx] = sort(gvec); |
|
|
200 rlpol = rlpol(:,idx); |
|
|
201 ngain = length(gvec); |
|
|
202 endwhile |
|
|
203 |
|
|
204 ## Plot the data |
|
|
205 if(nargout == 0) |
|
|
206 rlpolv = vec(rlpol); |
|
|
207 idx = find(real(rlpolv) >= xmin & real(rlpolv) <= xmax); |
|
|
208 axdata = [real(rlpolv(idx)),imag(rlpolv(idx))]; |
|
|
209 axlim = axis2dlim(axdata); |
|
|
210 axlim(1:2) = [xmin, xmax]; |
|
5215
|
211 __gnuplot_set__ nologscale xy; |
|
3431
|
212 grid("on"); |
|
|
213 rldata = [real(rlpolv), imag(rlpolv) ]; |
|
|
214 axis(axlim); |
|
|
215 [stn,inname,outname] = sysgetsignals(sys); |
|
|
216 xlabel(sprintf("Root locus from %s to %s, gain=[%f,%f]: Real axis", ... |
|
4771
|
217 inname{1}, outname{1},gvec(1),gvec(ngain))); |
|
3431
|
218 ylabel("Imag. axis"); |
|
|
219 |
|
5605
|
220 if(isempty(rlzer)) |
|
|
221 plot(real(rlpolv),imag(rlpolv),".1;locus points;", ... |
|
|
222 real(olpol),imag(olpol),"x2;open loop poles;"); |
|
|
223 else |
|
|
224 plot(real(rlpolv),imag(rlpolv),".1;locus points;", ... |
|
|
225 real(olpol),imag(olpol),"x2;open loop poles;", ... |
|
|
226 real(rlzer),imag(rlzer),"o3;zeros;"); |
|
|
227 endif |
|
3431
|
228 rldata = []; |
|
|
229 endif |
|
|
230 endfunction |
