|
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 |
|
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
|
18 |
|
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {@var{outputs} =} rlocus (@var{inputs}) |
|
|
21 ## @format |
|
|
22 ## [rldata, k] = rlocus(sys[,increment,min_k,max_k]) |
|
|
23 ## Displays root locus plot of the specified SISO system. |
|
|
24 ## |
|
|
25 ## ----- --- -------- |
|
|
26 ## --->| + |---|k|---->| SISO |-----------> |
|
|
27 ## ----- --- -------- | |
|
|
28 ## - ^ | |
|
|
29 ## |_____________________________| |
|
|
30 ## |
|
|
31 ## inputs: sys = system data structure |
|
|
32 ## min_k, max_k,increment: minimum, maximum values of k and |
|
|
33 ## the increment used in computing gain values |
|
|
34 ## Outputs: plots the root locus to the screen. |
|
|
35 ## rldata: Data points plotted column 1: real values, column 2: imaginary |
|
|
36 ## values) |
|
|
37 ## k: gains for real axis break points. |
|
|
38 ## @end format |
|
|
39 ## @end deftypefn |
|
|
40 |
|
|
41 ## Author: David Clem |
|
|
42 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
43 ## Updated by Kristi McGowan July 1996 for intelligent gain selection |
|
|
44 ## Updated by John Ingram July 1996 for systems |
|
|
45 |
|
|
46 function [rldata, k_break, rlpol, gvec, real_ax_pts] = rlocus (sys, increment, min_k, max_k) |
|
|
47 |
|
|
48 if (nargin < 1) | (nargin > 4) |
|
|
49 usage("rlocus(sys[,inc,mink,maxk])"); |
|
|
50 endif |
|
|
51 |
|
|
52 ## Convert the input to a transfer function if necessary |
|
|
53 |
|
|
54 [num,den] = sys2tf(sys) # extract numerator/denom polyomials |
|
|
55 lnum = length(num); lden = length(den); |
|
|
56 if(lden < 2) |
|
|
57 error(sprintf("length of derivative=%d, doesn't make sense",lden)); |
|
|
58 elseif(lnum == 1) |
|
|
59 num = [0, num]; # so that derivative is shortened by one |
|
|
60 endif |
|
|
61 |
|
|
62 ## root locus plot axis limits |
|
|
63 |
|
|
64 ## compute real axis locus breakpoints |
|
|
65 ## compute the derivative of the numerator and the denominator |
|
|
66 dern=polyderiv(num); derd=polyderiv(den); |
|
|
67 |
|
|
68 ## compute real axis breakpoints |
|
|
69 real_ax_pol = conv(den,dern) - conv(num,derd); |
|
|
70 real_ax_pts = roots(real_ax_pol); |
|
|
71 if(isempty(real_ax_pts)) |
|
|
72 k_break = []; |
|
|
73 maxk = 0; |
|
|
74 else |
|
|
75 ## compute gains that achieve the breakpoints |
|
|
76 c1 = polyval(num,real_ax_pts); |
|
|
77 c2 = polyval(den,real_ax_pts); |
|
|
78 k_break = -real(c2 ./ c1); |
|
|
79 maxk = max(max(k_break,0)); |
|
|
80 endif |
|
|
81 |
|
|
82 ## compute gain ranges based on computed K values |
|
|
83 if(maxk == 0) maxk = 1; |
|
|
84 else maxk = 1.1*maxk; endif |
|
|
85 mink = 0; |
|
|
86 ngain = 20; |
|
|
87 |
|
|
88 ## check for input arguments: |
|
|
89 if (nargin > 2) mink = min_k; endif |
|
|
90 if (nargin > 3) maxk = max_k; endif |
|
|
91 if (nargin > 1) |
|
|
92 if(increment <= 0) error("increment must be positive"); |
|
|
93 else |
|
|
94 ngain = (maxk-mink)/increment; |
|
|
95 endif |
|
|
96 endif |
|
|
97 |
|
|
98 ## vector of gains |
|
|
99 ngain = max(3,ngain); |
|
|
100 gvec = linspace(mink,maxk,ngain); |
|
|
101 |
|
|
102 ## Find the open loop zeros and the initial poles |
|
|
103 rlzer = roots(num); |
|
|
104 |
|
|
105 ## update num to be the same length as den |
|
|
106 lnum = length(num); if(lnum < lden) num = [zeros(1,lden - lnum),num]; endif |
|
|
107 |
|
|
108 ## compute preliminary pole sets |
|
|
109 nroots = lden-1; |
|
|
110 for ii=1:ngain |
|
|
111 gain = gvec(ii); |
|
|
112 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
113 endfor |
|
|
114 |
|
|
115 ## compute axis limits (isolate asymptotes) |
|
|
116 olpol = roots(den); |
|
|
117 real_axdat = union(real(rlzer), real(union(olpol,real_ax_pts)) ); |
|
|
118 rmin = min(real_axdat); rmax = max(real_axdat); |
|
|
119 |
|
|
120 rlpolv = [vec(rlpol); vec(real_axdat)]; |
|
|
121 idx = find(real(rlpolv) >= rmin & real(rlpolv) <= rmax); |
|
|
122 axlim = axis2dlim([real(rlpolv(idx)),imag(rlpolv(idx))]); |
|
|
123 xmin = axlim(1); |
|
|
124 xmax = axlim(2); |
|
|
125 |
|
|
126 ## set smoothing tolerance per axis limits |
|
|
127 smtol = 0.01*max(abs(axlim)); |
|
|
128 |
|
|
129 ## smooth poles if necessary, up to maximum of 1000 gain points |
|
|
130 ## only smooth points within the axis limit window |
|
|
131 ## smoothing done if max_k not specified as a command argument |
|
|
132 done=(nargin == 4); # perform a smoothness check |
|
|
133 while((!done) & ngain < 1000) |
|
|
134 done = 1 ; # assume done |
|
|
135 dp = abs(diff(rlpol'))'; |
|
|
136 maxd = max(dp); |
|
|
137 ## search for poles in the real axis limits whose neighbors are distant |
|
|
138 idx = find(maxd > smtol); |
|
|
139 for ii=1:length(idx) |
|
|
140 i1 = idx(ii); g1 = gvec(i1); p1 = rlpol(:,i1); |
|
|
141 i2 = idx(ii)+1; g2 = gvec(i2); p2 = rlpol(:,i2); |
|
|
142 |
|
|
143 ## isolate poles in p1, p2 that are inside the real axis limits |
|
|
144 bidx = find( (real(p1) >= xmin & real(p1) <= xmax) ... |
|
|
145 | (real(p2) >= xmin & real(p2) <= xmax) ); |
|
|
146 if(!isempty(bidx)) |
|
|
147 p1 = p1(bidx); |
|
|
148 p2 = p2(bidx); |
|
|
149 if( max(abs(p2-p1)) > smtol) |
|
|
150 newg = linspace(g1,g2,5); |
|
|
151 newg = newg(2:4); |
|
|
152 if(isempty(newg)) |
|
|
153 printf("rlocus: empty newg") |
|
|
154 g1 |
|
|
155 g2 |
|
|
156 i1 |
|
|
157 i2 |
|
|
158 idx_i1 = idx(ii) |
|
|
159 gvec_i1 = gvec(i1:i2) |
|
|
160 delta_vec_i1 = diff(gvec(i1:i2)) |
|
|
161 prompt |
|
|
162 endif |
|
|
163 gvec = [gvec,newg]; |
|
|
164 done = 0; # need to process new gains |
|
|
165 endif |
|
|
166 endif |
|
|
167 endfor |
|
|
168 |
|
|
169 ## process new gain values |
|
|
170 ngain1 = length(gvec); |
|
|
171 for ii=(ngain+1):ngain1 |
|
|
172 gain = gvec(ii); |
|
|
173 rlpol(1:nroots,ii) = vec(sortcom(roots(den + gain*num))); |
|
|
174 endfor |
|
|
175 |
|
|
176 [gvec,idx] = sort(gvec); |
|
|
177 rlpol = rlpol(:,idx); |
|
|
178 ngain = length(gvec); |
|
|
179 endwhile |
|
|
180 |
|
|
181 ## Plot the data |
|
|
182 if(nargout == 0) |
|
|
183 rlpolv = vec(rlpol); |
|
|
184 idx = find(real(rlpolv) >= xmin & real(rlpolv) <= xmax); |
|
|
185 axdata = [real(rlpolv(idx)),imag(rlpolv(idx))]; |
|
|
186 axlim = axis2dlim(axdata); |
|
|
187 axlim(1:2) = [xmin, xmax]; |
|
|
188 gset nologscale xy; |
|
|
189 grid("on"); |
|
|
190 rldata = [real(rlpolv), imag(rlpolv) ]; |
|
|
191 axis(axlim); |
|
|
192 [stn,inname,outname] = sysgetsignals(sys); |
|
|
193 xlabel(sprintf("Root locus from %s to %s, gain=[%f,%f]: Real axis", ... |
|
|
194 nth(inname,1),nth(outname,1),gvec(1),gvec(ngain))); |
|
|
195 ylabel("Imag. axis"); |
|
|
196 |
|
|
197 plot(real(rlpolv),imag(rlpolv),".1;locus points;", ... |
|
|
198 real(olpol),imag(olpol),"x2;open loop poles;", ... |
|
|
199 real(rlzer),imag(rlzer),"o3;zeros;"); |
|
|
200 rldata = []; |
|
|
201 endif |
|
|
202 endfunction |
