home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / CTRBF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  78 lines
  1. function [abar,bbar,cbar,t,k] = ctrbf(a, b, c, tol)
  2. %CTRBF    Controllability staircase form.
  3. %    [ABAR,BBAR,CBAR,T,K] = CTRBF(A,B,C) returns a decomposition
  4. %    into the controllable/uncontrollable subspaces.
  5. %    [ABAR,BBAR,CBAR,T,K] = CTRBF(A,B,C,TOL) uses tolerance TOL.
  6. %
  7. %    If Co=CTRB(A,B) has rank r <= n, then there is a similarity
  8. %    transformation T such that
  9. %
  10. %    Abar = T * A * T' ,  Bbar = T * B ,  Cbar = C * T'
  11. %
  12. %    and the transformed system has the form
  13. %
  14. %           | Anc    0 |           | 0 |
  15. %    Abar =  ----------  ,  Bbar =  ---  ,  Cbar = [Cnc| Cc].
  16. %           | A21   Ac |           |Bc |
  17. %                                               -1          -1
  18. %    where (Ac,Bc) is controllable, and Cc(sI-Ac)Bc = C(sI-A)B.
  19. %
  20. %    See also: CTRB.
  21.  
  22. %    Author : R.Y. Chiang  3-21-86
  23. %    Revised 5-27-86 JNL
  24. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  25.  
  26. % This M-file implements the Staircase Algorithm of Rosenbrock, 1968.
  27. [ra,ca] = size(a);
  28. [rb,cb] = size(b);
  29. %
  30. % ------ Assign Initial Conditions :
  31. %
  32. ptjn1 = eye(ra);
  33. ajn1 = a;
  34. bjn1 = b;
  35. rojn1 = cb;
  36. deltajn1 = 0;
  37. sigmajn1 = ra;
  38. k = zeros(1,ra);
  39. if nargin == 3
  40.     tol = ra*norm(a,1)*eps;
  41. end
  42. %
  43. % ------ Begin Major Loop :
  44. %
  45. for jj = 1 : ra
  46.     [uj,sj,vj] = svd(bjn1);
  47.     [rsj,csj] = size(sj);
  48.     p = rot90(eye(rsj),1);
  49.     uj = uj*p;
  50.     bb = uj'*bjn1;
  51.     roj = rank(bb,tol);
  52.     [rbb,cbb] = size(bb);
  53.     sigmaj = rbb - roj;
  54.     sigmajn1 = sigmaj;
  55.     k(jj) = roj;
  56.     if roj == 0, break, end
  57.     if sigmaj == 0, break, end
  58.     abxy = uj' * ajn1 * uj;
  59.     aj   = abxy(1:sigmaj,1:sigmaj);
  60.     bj   = abxy(1:sigmaj,sigmaj+1:sigmaj+roj);
  61.     ajn1 = aj;
  62.     bjn1 = bj;
  63.     [ruj,cuj] = size(uj);
  64.     ptj = ptjn1 * ...
  65.           [uj zeros(ruj,deltajn1); ...
  66.            zeros(deltajn1,cuj) eye(deltajn1)];
  67.     ptjn1 = ptj;
  68.     deltaj = deltajn1 + roj;
  69.     deltajn1 = deltaj;
  70. end
  71. %
  72. % ------ Final Transformation :
  73. %
  74. t = ptjn1';
  75. abar = t * a * t';
  76. bbar = t * b;
  77. cbar = c * t';
  78.