home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / ARE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.5 KB  |  50 lines
  1. function X = are(A,B,C)
  2. %ARE    Algebraic Riccati Equation solution.
  3. %    X = ARE(A, B, C) returns the stablizing solution (if it
  4. %       exists) to the continuous-time Riccati equation:
  5. %
  6. %            A'*X + X*A - X*B*X + C = 0
  7. %
  8. %       assuming B is symmetric and nonnegative definite and C is
  9. %       symmetric.
  10. %
  11. %    See also: RIC.
  12.  
  13. % M. Wette, ECE Dept., Univ. of California 5-11-87
  14. % Revised 6-16-87 MW
  15. %         3-16-88 MW
  16. %    Copyright (c) 1986-93 by The MathWorks, Inc.
  17.  
  18. %  -- check for correct input problem --
  19. [nr,nc] = size(A); n = nr;
  20. if (nr ~= nc), error('Nonsquare A matrix'); end;
  21. [nr,nc] = size(B);
  22. if (nr~=n | nc~=n), error('Incorrectly dimensioned B matrix'); end;
  23. [nr,nc] = size(C);
  24. if (nr~=n | nc~=n), error('Incorrectly dimensioned C matrix'); end;
  25.  
  26. [q,t] = schur([A -B; -C -A']*(1.0+eps*eps*sqrt(-1))); 
  27. tol = 10.0*eps*max(abs(diag(t)));    % ad hoc tolerance
  28. ns = 0;
  29. %
  30. %  Prepare an array called index to send message to ordering routine 
  31. %  giving location of eigenvalues with respect to the imaginary axis.
  32. %  -1  denotes open left-half-plane
  33. %   1  denotes open right-half-plane
  34. %   0  denotes within tol of imaginary axis
  35. %  
  36. for i = 1:2*n,
  37.     if (real(t(i,i)) < -tol),
  38.         index = [ index -1 ];
  39.     ns = ns + 1;
  40.     elseif (real(t(i,i)) > tol),
  41.     index = [ index 1 ];
  42.     else,
  43.     index = [ index 0 ];
  44.     end;
  45. end;
  46. if (ns ~= n), error('No solution: (A,B) may be uncontrollable or no solution exists'); end;
  47. [q,t] = schord(q,t,index);
  48. X = real(q(n+1:n+n,1:n)/q(1:n,1:n));
  49.  
  50.