home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / RICCOND.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.4 KB  |  83 lines
  1. function [tot] = riccond(a,b,qrn,p1,p2);
  2. %
  3. % [TOT] = RICCOND(A,B,QRN,P1,P2) produces a Riccati condition
  4. %    analysis using the concept of Arnold and Laub (1984).
  5. %
  6. %    INPUTS:
  7. %
  8. %      System: (A,B) matrices, and weighting QRN = [Q N;N' R];
  9. %      Riccati Solution: P1,P2 (P = P2/P1)
  10. %
  11. %    OUTPUT:
  12. %
  13. %         TOT = [norAc norQc norRc conr conP1 conArn conBey res]
  14. %     where                                         -1            -1    
  15. %        norAc, norQc, norRc ----- F-Norm of Ac=(A-BR N'), Qc=(Q-NR N'), 
  16. %                                  and Rc = B/R*B'.
  17. %        conr ----- condition number of R
  18. %        conP1 ---- condition number of P1
  19. %        conArn --- Arnold and Laub's Riccati condition number
  20. %        conBey --- Byers's Riccati condition number
  21. %        res ------ residual of Riccati equation
  22. %
  23.  
  24. % R. Y. Chiang & M. G. Safonov 9/18/1988
  25. % Copyright (c) 1988 by the MathWorks, Inc.
  26. % All Rights Reserved.
  27. % -----------------------------------------------------------------
  28. [ns,ns] = size(a);
  29. [q,n,nt,r] = sys2ss(qrn,ns);
  30. %
  31. ac = a - b/r*nt;
  32. rc = b/r*b';
  33. qc = q - n/r*nt;
  34. %
  35. ham = [ac -rc;-qc -ac'];
  36. perr = [p2' -p1']*ham*[p1;p2];
  37. p = p2/p1;
  38. %
  39. clc
  40. format short e
  41. disp(' -----------------------------------');
  42. disp('    Riccati Condition Analysis ');
  43. disp(' -----------------------------------');
  44. disp(' ')
  45. disp(' 1). F-Norm of Ac, Qc and B/R*B`:');
  46. norAc = norm(ac,'fro');
  47. norQc = norm(qc,'fro');
  48. norRc = norm(rc,'fro');
  49. norAQR = [norAc norQc norRc],
  50. disp('                        hit <ret> to continue');
  51. pause
  52. %
  53. disp(' 2). Condition no. of R:');
  54. conr = cond(r),
  55. disp('                        hit <ret> to continue');
  56. pause
  57. %
  58. disp(' 3). Condition no. of P1:');
  59. conP1 = cond(p1),
  60. disp('                        hit <ret> to continue');
  61. pause
  62. %
  63. disp(' 4). Arnold and Laub Riccati Condition no.:');
  64. acl = ac - rc*p;
  65. sepArn = min(svd(kron(eye(acl'),acl')+kron(acl',eye(acl'))));
  66. conArn = norm(qc,'fro')/norm(p,'fro')/sepArn,
  67. disp('                        hit <ret> to continue');
  68. pause
  69. %
  70. disp(' 5). Byers Riccati Condition no.:');
  71. sepBye = sepArn;
  72. conBye = (norm(qc,'fro')/norm(p,'fro')+2*norm(ac,'fro')...
  73.          +norm(rc,'fro')*norm(p,'fro'))/sepBye,
  74. disp('                        hit <ret> to continue');
  75. pause
  76. %
  77. disp(' 6). Residual = norm(perr,1)/norm(p,1)');
  78. res = norm(perr,1)/norm(p,1),
  79. %
  80. tot = [norAc norQc norRc conr conP1 conArn conBye res]';
  81. %
  82. % ------ End of RICCOND.M % RYC/MGS %
  83.