home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / place.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  3.7 KB  |  127 lines
  1. ## Copyright (C) 1997 Jose Daniel Munoz Frias
  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{K} =} place (@var{sys}, @var{P})
  21. ## Computes the matrix  K such that if the state
  22. ## is feedback with gain K, then the eigenvalues  of the closed loop
  23. ## system (i.e. A-BK) are those specified in the vector P.
  24. ## 
  25. ## Version: Beta (May-1997): If you have any comments, please let me know.
  26. ##                 (see the file place.m for my address)
  27. ## 
  28. ## Written by: Jose Daniel Munoz Frias.
  29. ## @end deftypefn
  30.  
  31. function K = place (sys, P) 
  32.  
  33.   ## Universidad Pontificia Comillas
  34.   ## ICAIdea
  35.   ## Alberto Aguilera, 23
  36.   ## 28015 Madrid, Spain
  37.   ##
  38.   ## E-Mail: daniel@dea.icai.upco.es
  39.   ##
  40.   ## Phone: 34-1-5422800   Fax: 34-1-5596569
  41.   ##
  42.   ## Algorithm taken from "The Control Handbook", IEEE press pp. 209-212
  43.   ##
  44.   ## code adaped by A.S.Hodel (a.s.hodel@eng.auburn.edu) for use in controls
  45.   ## toolbox
  46.  
  47.   sav_val = empty_list_elements_ok;
  48.   empty_list_elements_ok = 1;
  49.  
  50.   ## check arguments
  51.  
  52.   if(!is_struct(sys))
  53.     error("sys must be in system data structure format (see ss2sys)");
  54.   endif
  55.   sys = sysupdate(sys,"ss");    # make sure it has state space form up to date
  56.   if(!is_controllable(sys))
  57.     error("sys is not controllable.");
  58.   elseif( min(size(P)) != 1)
  59.     error("P must be a vector")
  60.   else
  61.     P = reshape(P,length(P),1);    # make P a column vector
  62.   endif
  63.   ## system must be purely continuous or discrete
  64.   is_digital(sys);
  65.   [n,nz,m,p] = sysdimensions(sys);
  66.   nx = n+nz;    # already checked that it's not a mixed system.
  67.   if(m != 1)
  68.     error(["sys has ", num2str(m)," inputs; need only 1"]);
  69.   endif
  70.  
  71.   ## takes the A and B matrix from the system representation
  72.   [A,B]=sys2ss(sys);
  73.   sp = length(P);
  74.   if(nx == 0)
  75.     error("place: A matrix is empty (0x0)");
  76.   elseif(nx != length(P))
  77.     error(["A=(",num2str(nx),"x",num2str(nx),", P has ", num2str(length(P)), ...
  78.     "entries."])
  79.   endif
  80.  
  81.   ## arguments appear to be compatible; let's give it a try!
  82.   ## The second step is the calculation of the characteristic polynomial ofA
  83.   PC=poly(A);
  84.  
  85.   ## Third step: Calculate the transformation matrix T that transforms the state
  86.   ## equation in the controllable canonical form.
  87.  
  88.   ## first we must calculate the controllability matrix M:
  89.   M=B;
  90.   AA=A;
  91.   for n = 2:nx
  92.     M(:,n)=AA*B;
  93.     AA=AA*A;
  94.   endfor
  95.  
  96.   ## second, construct the matrix W
  97.   PCO=PC(nx:-1:1);
  98.   PC1=PCO;     # Matrix to shift and create W row by row
  99.  
  100.   for n = 1:nx
  101.     W(n,:) = PC1;
  102.     PC1=[PCO(n+1:nx),zeros(1,n)];
  103.   endfor
  104.  
  105.   T=M*W;
  106.  
  107.   ## finaly the matrix K is calculated 
  108.   PD = poly(P); # The desired characteristic polynomial
  109.   PD = PD(nx+1:-1:2);
  110.   PC = PC(nx+1:-1:2);
  111.   
  112.   K = (PD-PC)/T;
  113.  
  114.   ## Check if the eigenvalues of (A-BK) are the same specified in P
  115.   Pcalc = eig(A-B*K);
  116.  
  117.   Pcalc = sortcom(Pcalc);
  118.   P = sortcom(P);
  119.  
  120.   if(max( (abs(Pcalc)-abs(P))./abs(P) ) > 0.1)
  121.     disp("Place: Pole placed at more than 10% relative error from specified");
  122.   endif
  123.  
  124.   empty_list_elements_ok = sav_val;
  125. endfunction
  126.  
  127.