home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / SNCFBAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.7 KB  |  119 lines
  1. % function [sysnlcf,sig,sysnrcf] = sncfbal(sys,tol)
  2. %
  3. %   Finds a (truncated)  balanced realization of the normalized 
  4. %   left and right coprime factors of the system state-space 
  5. %   model SYS. The state-space models in SYSNLCF and SYSNRCF
  6. %   are both balanced realizations with Hankel singular values 
  7. %   given by SIG.  
  8. % 
  9. %    SYSNLCF = [Nl Ml] and Nl Nl~ + Ml Ml~ = I, SYS = Ml^(-1) Nl.
  10. %    SYSNRCF = [Nr;Mr] and Nr~ Nr + Mr~ Mr = I, SYS = Nr Mr^(-1).
  11. %
  12. %   The results are truncated to retain all Hankel singular values 
  13. %   greater than TOL. If TOL is omitted then it is set to 
  14. %   MAX(SIG(1)*1.0E-12,1.0E-16). A warning message is printed out 
  15. %   if SYS is close to being undetectable, in which case the output 
  16. %   may be unreliable.
  17. %
  18. %   See also: HANKMR, SFRWTBAL, SFRWTBLD, SRELBAL, SYSBAL,
  19. %             SRESID, and TRUNC.
  20.  
  21. function [sysnlcf,sig,sysnrcf]  = sncfbal(sys,tol)
  22.  
  23.    if nargin<1
  24.      disp('usage: [sysnlcf,sig,sysnrcf] = sncfbal(sys,tol)')
  25.      return
  26.    end
  27.    [A,B,C,D]=unpck(sys);
  28.    [systype,p,m,n]=minfo(sys);
  29.    if systype~='syst'
  30.      disp('SYS must be SYSTEM')
  31.      return
  32.    end
  33.  
  34.  %  Solve the GFARE soln X=S'*S
  35.  
  36.    Ham = [A' zeros(n); -B*B' -A] - [C'; -B*D']*inv(eye(p)+D*D')*[D*B' C];
  37.    [U,Hamm]=schur(Ham);
  38.    if any(imag(sys)), [U,Hamm]=ocsf(U,Hamm,'s');
  39.      else [U,Hamm]=orsf(U,Hamm,'s'); end
  40.    [S,SS]=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);
  41.     S=S*sqrt(abs(SS))*S';
  42.    [qdr,Dr]=qr([D;eye(m)]);Dr=inv(Dr(1:m,1:m));
  43.    [qdl,Dl]=qr([D';eye(p)]);Dl=inv(Dl(1:p,1:p)');
  44.  
  45. % change coordinates using U(1:n,1:n)
  46.  
  47.    nmuinv= norm(inv(U(1:n,1:n)));
  48.    if nmuinv>1e12, 
  49.      disp(['Warning from SNCFBAL: SYS is close to being undetectable',...
  50.            ' results maybe unreliable'])
  51.      return
  52.    end %if nmuinv>1e12,
  53.    B=U(1:n,1:n)'*B; %temp variable
  54.    Cl=Dl*(C/U(1:n,1:n)');
  55.    C=C*U(n+1:2*n,1:n); %temp variable
  56.    Al=Hamm(1:n,1:n)';
  57.    Bl=(B-C'*D)*Dr*Dr';        
  58.    Hl=-(C'+B*D') *Dl'*Dl;       
  59.  
  60. % realization of nlcf [N M] is now pck(Al,[Bl Hl],Cl,[Dl*D;Dl])
  61. % Calculate observability Gramian for nrcf [Nl Ml], R*R'
  62.  
  63.    nm11=n:-1:1;
  64.    R=sjh6(Al(nm11,nm11),Cl(:,nm11));
  65.    R=R(nm11,nm11)';
  66.  
  67. % calculate the Hankel-singular values of [Nl Ml]
  68.  
  69.    [W,T,V] = svd(S*R);
  70.    sig = diag(T);
  71.  
  72. % balancing coordinates
  73.  
  74.    T = W'*S;
  75.    Cl = Cl*T'; Al = Al*T'; 
  76.    T = R*V; 
  77.    Bl = T'*Bl; Al = T'*Al; Hl = T'*Hl;
  78.  
  79. % calculate the truncated dimension nn
  80.  
  81.    if nargin<2, tol=max([sig(1)*1.0E-12,1.0E-16]);end;
  82.    nn = n;
  83.    for i=n:-1:1, if sig(i)<=tol nn=i-1; end; end;
  84.    if nn==0, sysnrcf=[Dl*D Dl];sig=[];
  85.      else
  86.      sig = sig(1:nn);
  87.      % diagonal scaling  by sig(i)^(-0.5)
  88.     irtsig = sig.^(-0.5);
  89.      onn=1:nn;
  90.      Al=Al(onn,onn).*(irtsig*irtsig');
  91.      Bl=(irtsig*ones(1,m)).*Bl(onn,:);
  92.      Cl=Cl(:,onn).*(ones(p,1)*irtsig');
  93.      Hl=Hl(onn,:).*(irtsig*ones(1,p));
  94.      sysnlcf=pck(Al,[Bl Hl],Cl,[Dl*D Dl]);
  95.     end %if nn==0,
  96.  
  97. %now calculate the right coprime factorization if nargout>2.
  98.  
  99.  if nargout>2,
  100.   if nn==0, sysnrcf=[D*Dr;Dr];
  101.    elseif (sig(1)<0.99)&(nmuinv<=1e12),
  102.      roms2=sqrt(1-sig.^2);
  103.      sysnrcf=[zeros(nn+p+m,nn+m) nn*eye(nn+p+m,1); zeros(1,nn+m), -inf];
  104.      B=(Bl-Hl*D)*Dr;
  105.      sysnrcf(nn+p+1:nn+m+p,1:nn)=...
  106.           -(Dr*B').*(ones(m,1)*(sig./roms2)')...
  107.           -(D'*Dl'*Cl).*(ones(m,1)*roms2');
  108.     sysnrcf(nn+1:nn+p,1:nn)=...
  109.           (Dl\Cl).*(ones(p,1)*roms2')+D*sysnrcf(nn+p+1:nn+m+p,1:nn);
  110.     sysnrcf(1:nn+m+p,nn+1:nn+m)=[B./(roms2*ones(1,m));D*Dr;Dr];
  111.     sysnrcf(1:nn,1:nn)=(Al).*(roms2*(roms2.^(-1))');
  112.    else
  113.     sysnrcf = transp(sncfbal(transp(sys)));
  114.   end %if nn==0, sysnrcf
  115.  end, % if nargout>2
  116.  
  117. %
  118. % Copyright MUSYN INC 1991,  All Rights Reserved
  119.