home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / rlocus.m < prev    next >
Text File  |  1999-12-24  |  7KB  |  203 lines
  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 } { outputs =} rlocus ( 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. ## 
  39. ## 
  40. ## @end format
  41. ## @end deftypefn
  42.  
  43. function [rldata, k_break, rlpol, gvec, real_ax_pts] = rlocus (sys, increment, min_k, max_k)
  44.  
  45.   ## Convert the input to a transfer function if necessary
  46.   ## Written by Clem and Tenison
  47.   ## Updated by Kristi McGowan July 1996 for intelligent gain selection
  48.   ## Updated by John Ingram July 1996 for systems
  49.   
  50.   if (nargin < 1) | (nargin > 4)
  51.     usage("rlocus(sys[,inc,mink,maxk])");
  52.   endif
  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=polyderv(num);        derd=polyderv(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 = ax2dlim([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 = ax2dlim(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] = sysgetsg(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
  203.