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