home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / SRELBAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.3 KB  |  70 lines
  1. % function [sysout,sig,sysfact] = srelbal(sys,tol)
  2. %
  3. %   Find a truncated stochastically balanced realization of
  4. %   the system state-space model SYS.  All the eigenvalues
  5. %   of SYS must have negative real part. The result is
  6. %   truncated to retain all Hankel singular values greater
  7. %   than TOL.  The output parameter SYSFACT is a stable
  8. %   minimum phase system such that SYS~*SYS = SYSFACT*SYSFACT~.
  9. %   If TOL is omitted then it is set to max(SIG(1)*1.0E-12,1.0E-16).
  10. %   SIG imply the achievable relative errors with STRUNC(SYSOUT,K).
  11. %
  12. %   See also: SFRWTBAL, SFRWTBLD, SNCFBAL, SRELBAL, SYSBAL,
  13. %             SRESID, and TRUNC.
  14.  
  15. function [sysout,sig,sysfact] = srelbal(sys,tol)
  16.  
  17.    if nargin<1
  18.      disp('usage: [sysout,sig,sysfact] = srelbal(sys,tol)');
  19.      return;
  20.    end
  21.    [A,B,C,D]=unpck(sys);
  22.    [systype,p,m,n]=minfo(sys);
  23.    if systype~='syst', error('sys must be SYSTEM');end
  24.    [T,A]=schur(A);
  25.    B = T'*B;
  26.    C = C*T;
  27.  
  28.  % find observability Gramian, S'*S (S upper triangular)
  29.  
  30.    S = sjh6(A,C);
  31.    H = B'*S'*S+D'*C;
  32.    Ham = [A' zeros(n); zeros(n) -A] - [H';B]*inv(D'*D)*[B' -H];
  33.    [U,Hamm]=schur(Ham);
  34.    if any(imag(sys)), [U,Hamm]=ocsf(U,Hamm,'s');
  35.      else [U,Hamm]=orsf(U,Hamm,'s'); end
  36.    [R,RR]=schur((U(1:n,1:n)'*U(n+1:2*n,1:n)+U(n+1:2*n,1:n)'*U(1:n,1:n))/2);
  37.     R=(U(n+1:2*n,1:n)*R)./(ones(n,1)*sqrt(abs(diag(RR)')));
  38.    [qdw,Dw]=qr(D);Dw=Dw(1:m,1:m)';
  39.    L=(B-R*R'*H')/(Dw');
  40.    % calculate the Hankel-singular values 
  41.    [U,T,V] = svd(S*R);
  42.    sig = diag(T);
  43.  
  44.  % balancing coordinates
  45.  
  46.    T = U'*S;
  47.    B = T*B; A = T*A; L=T*L;
  48.    T = R*V; 
  49.    C = C*T; A = A*T; H=H*T;
  50.     % calculate the truncated dimension nn
  51.    if nargin<2 tol=max([sig(1)*1.0E-12,1.0E-16]);end;
  52.    nn = n;
  53.    for i=n:-1:1, if sig(i)<=tol nn=i-1; end; end;
  54.    if nn==0, sysout=D; sysfact=Dw;
  55.      else
  56.      sig = sig(1:nn);
  57.      % diagonal scaling  by sig(i)^(-0.5)
  58.      irtsig = sig.^(-0.5);
  59.      onn=1:nn;
  60.      A(onn,onn)=A(onn,onn).*(irtsig*irtsig');
  61.      B(onn,:)=(irtsig*ones(1,m)).*B(onn,:);
  62.      C(:,onn)=C(:,onn).*(ones(p,1)*irtsig');
  63.      L(onn,:)=(irtsig*ones(1,m)).*L(onn,:);
  64.      H(:,onn)=H(:,onn).*(ones(m,1)*irtsig');
  65.      sysout=pck(A(onn,onn),B(onn,:),C(:,onn),D);
  66.      sysfact=pck(A(onn,onn),L(onn,:),H(:,onn),Dw);
  67.     end
  68. %
  69. % Copyright MUSYN INC 1991,  All Rights Reserved
  70.