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