home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / ARESOLV.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.1 KB  |  90 lines
  1. function [p1,p2,lamp,perr,wellposed,p] = aresolv(a,q,R,Type)
  2. %
  3. % [P1,P2,LAMP,PERR,WELLPOSED,P] = ARESOLV(A,Q,R,TYPE) computes the 
  4. % solution to the algebraic Riccati equation
  5. % 
  6. %                     A'P + PA - PRP + Q = 0
  7. %
  8. %   based on the "Type" one selects:
  9. % 
  10. %                Type = 'eigen' ---- Eigenstructure approach
  11. %
  12. %                Type = 'Schur' ---- Schur vector approach
  13. %
  14. %  Outputs:   P = P2/P1 riccati solution
  15. %             [P1;P2] stable eigenspace of hamiltonian [A,-R;-Q,-A']
  16. %             LAMP closed loop eigenvalues
  17. %             PERR residual error matrix 
  18. %             WELLPOSED = 'TRUE ' if hamiltonian has no imaginary-axis
  19. %                  eigenvalues, otherwise 'FALSE'
  20.  
  21. % R. Y. Chiang & M. G. Safonov 2/88
  22. % Copyright (c) 1988 by the MathWorks, Inc.
  23. % All Rights Reserved.
  24.  
  25. wellposed='TRUE ';  %initialize 
  26. if nargin == 3
  27.    Type = 'eigen';
  28. end
  29. %
  30. [n,n] = size(a);
  31. %
  32. % ------ If A is stable and Q is zero, p2 = zero(n), p1 = eye(n)
  33. %
  34. lamA = eig(a);
  35. ind = find(lamA>0);
  36. if (sum(ind) == 0) & (norm(q,'fro') == 0)
  37.    p2 = zeros(n); p1 = eye(n);
  38.    p = p2; perr = p2;
  39.    return
  40. end
  41. %   
  42. ham = real([a -R;-q -a']);
  43. %
  44. % ------ Eigenstructure approach:
  45. %
  46. if Type == 'eigen'
  47.    [v,d] = reig(ham);
  48.    if rank(v) < 2*n
  49.       Type = 'Schur';
  50.    end
  51.    p1 = v(1:n,1:n);
  52.    p2 = v((n+1):(2*n),1:n);
  53. end
  54. %
  55. % ------ Schur vector approach:
  56. %
  57. if Type == 'Schur'
  58.    [v,t,m,swap] = cschur(ham,1);
  59.    vv = [real(v(:,1:m)) imag(v(:,1:m))];
  60.    [vh,rh] = qr(vv);
  61.    p1 = vh(1:n,1:n);
  62.    p2 = vh(n+1:2*n,1:n);
  63. end
  64. %
  65. if nargout == 6
  66.    p = real(p2/p1);
  67. end
  68. %
  69. % ------ Residual of the Riccati solution:
  70. %
  71. perr = [p2' -p1']*[a -R;-q -a']*[p1;p2];
  72. %
  73. % ------ Check jw-axis roots:
  74. %
  75. [gvv,dd] = eig(a*p1-R*p2,p1);
  76. %
  77. % ------ if the generalized e-value is NaN or Inf, assign to -1.e15:
  78. %
  79. for i = 1:n
  80.    if abs(dd(i,i)) == NaN | abs(dd(i,i)) == inf
  81.       dd(i,i) = -1.e15;
  82.    end
  83. end
  84. lamp = diag(dd);
  85. %
  86. if max(real(lamp)+1.e-6*abs(imag(lamp))) > -1.e-13
  87.      wellposed='FALSE';
  88. end
  89. %
  90. % ------- End of ARESOLV.M ---- % RYC/MGS %