home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / DARESOLV.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  78 lines
  1. function [p1,p2,lamp,perr,wellposed,p] = daresolv(A,B,Q,R,Type)
  2. %
  3. % [P1,P2,LAMP,PERR,WELLPOSED,P] = DARESOLV(A,B,Q,R,TYPE) computes the 
  4. % solution to the discrete algebraic Riccati equation
  5. %                                   -1
  6. %            A'PA - P - A'PB(R+B'PB)  B'PA + 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 matrix
  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. RR = B/R*B';
  43. ham = [A+RR/A'*Q -RR/A';-A'\Q inv(A')];
  44. %
  45. % ------ Eigenstructure approach:
  46. %
  47. if Type == 'eigen'
  48.    [v,d] = reig(ham,2);
  49.    if rank(v) < 2*n
  50.       Type = 'Schur';
  51.    end
  52.    p1 = v(1:n,1:n);
  53.    p2 = v((n+1):(2*n),1:n);
  54. end
  55. %
  56. % ------ Schur vector approach:
  57. %
  58. if Type == 'Schur'
  59.    [v,t,m,swap] = cschur(ham,6);
  60.    vv = [real(v(:,1:m)) imag(v(:,1:m))];
  61.    [vh,rh] = qr(vv);
  62.    p1 = vh(1:n,1:n);
  63.    p2 = vh(n+1:2*n,1:n);
  64. end
  65. %
  66. P = real(p2/p1);
  67. %
  68. perr = A'*P*A-P-A'*P*B/(R+B'*P*B)*B'*P*A+Q;
  69. %
  70. % ------ Check jw-axis roots:
  71. %
  72. lamp = eig(A-B/(R+B'*P*B)*B'*P'*A);
  73. %
  74. if max(abs(lamp)) > 1-1.e-8
  75.      wellposed='FALSE';
  76. end
  77. %
  78. % ------- End of DARESOLV.M ---- % RYC/MGS %