home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / HANKMR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.2 KB  |  110 lines
  1. % function [sh,su,hsvu]=hankmr(sys,sig,k,opt)
  2. %
  3. %   Performs an optimal Hankel norm approximation of order K
  4. %   on the state-space system SYS of order n. The input system
  5. %   SYS must be a balanced realization with Hankel singular
  6. %   values SIG (i.e. SYS and SIG must be the output of SYSBAL).
  7. %   OPT may be omitted or
  8. %
  9. %     -set to 'a' to stop with an anti-causal term. The anticausal
  10. %         part will be absorbed into SH, and SU = 0.
  11. %     -set to 'd' to calculate a D-term which satisfies
  12. %         the H infinity error bound the sum of the neglected
  13. %         Hankel singular values SIG(i). That is SH+SU will always
  14. %         be an optimal H infinity norm approximation to SYS with
  15. %         K stable poles. 
  16. %
  17. %   The anticausal will be absorbed into SH, and SU = 0, if OPT = 'a',
  18. %   otherwise SU contains the anticausal term and SH the Kth order
  19. %   causal system. With the 'd' option the Hankel singular values
  20. %   of SU are given in HSVU.
  21. %
  22. %   See also: SFRWTBAL, SFRWTBLD, SNCFBAL, SRELBAL, SYSBAL, 
  23. %             SRESID, and TRUNC.
  24.  
  25. function [sh,su,hsvu]=hankmr(sys,sig,k,opt)
  26. if nargin < 3
  27.   disp(['usage: [sh,su,hsvu] = hankmr(sys,sig,k,opt)']);
  28.   return
  29. end
  30. if nargin==3,
  31.   opt='s'; 
  32. end
  33.  
  34. [mattype,p,m,n]=minfo(sys);
  35. if k<0 | k>=n,error('k must be between 0 and n-1');return;end;
  36. [a,b,c,d]=unpck(sys);
  37. kp=k+1;
  38. skp=sig(kp);skp2=skp^2;
  39. if k>0,
  40.   while sig(k)<1.01*skp&k>0,
  41.     k=k-1;
  42.     disp('sigma(k)<1.01*sigma(k+1), k reduced');
  43.   end;%while sig(k)
  44. end%if k>0
  45. kp=k+1;
  46. skp=sig(kp);skp2=skp^2;
  47. for r=1:n-k,
  48.   if sig(k+r)<0.99*skp,
  49.     kr=k+r-1;break;
  50.   end;
  51.   kr=k+r;
  52. end;
  53. r=kr-k;sig=sig';krp=kr+1;
  54. u=-c(:,kp:kr)*pinv(b(kp:kr,:)');
  55. ds=diag([(sig(1:k).^2-skp2*ones(1,k)).^(-0.5),-(skp2*...
  56.    ones(1,n-kr)-sig(krp:n).^2).^(-0.5)]);
  57. sig1=diag([sig(1:k),sig(krp:n)]);
  58. aw=ds*(sig1*a([1:k,krp:n],[1:k,krp:n])*sig1+skp2*a([1:k,krp:n],...
  59.    [1:k,krp:n])'-skp*c(:,[1:k,krp:n])'*u*b([1:k,krp:n],:)')*ds;
  60. if kp<=(n-r), aw(:,kp:(n-r))=-aw(:,kp:(n-r));
  61.   end % if kp<=(n-r)
  62. bw=ds*(sig1*b([1:k,krp:n],:)+skp*c(:,[1:k,krp:n])'*u);
  63. cw=(c(:,[1:k,krp:n])*sig1+skp*u*b([1:k,krp:n],:)')*ds;
  64. if kp<=(n-r), cw(:,kp:(n-r))=-cw(:,kp:(n-r));
  65.   end %if kp<=(n-r),
  66. dw=d-skp*u;
  67. if opt=='a',
  68.     sh=pck(aw,bw,cw,dw);
  69.     return;
  70.  end; %if opt=='a'
  71. if krp>n, %i.e. dim of SU is zero
  72.   su=zeros(p,m);
  73.   if n>r, sh=pck(aw,bw,cw,dw);
  74.    else, sh=dw; 
  75.    end %if n>r,
  76. else, %if krp>n
  77.   if k==0,
  78.      au=aw;bu=bw;cu=cw;du=zeros(p,m); 
  79.      su=pck(au,bu,cu,du);
  80.      dh=dw; sh=dh;
  81.    end; %if k==0
  82.   if k>0,
  83.     [sh,su]=sdecomp(pck(aw,bw,cw,dw));
  84.     [au,bu,cu,du]=unpck(su);
  85.    end; %if k>0
  86.   if opt=='d',
  87.     if (m==1)&(p==1), du=0.5*(cu/au)*bu;
  88.         [su,hsvu]=sysbal(pck(-au,bu,cu,du));
  89.         [au,bu,cu,du]=unpck(su);
  90.       else
  91.       [su,hsvu]=sysbal(pck(-au,bu,cu,du));
  92.       [au,bu,cu,du]=unpck(su);
  93.       if ~isempty(au),
  94.          du=-dcalc(bu,cu,hsvu,2*(p+m));
  95.         end; %if ~isempty(au)
  96.      end; %if (m==1)
  97.     dh=dw-du;
  98.     if k==0,
  99.       sh=dh;
  100.      else,
  101.       sh(kp:k+p,kp:k+m)=dh;
  102.      end; %if k==0
  103.     if ~isempty(au),
  104.        su=pck(-au,bu,cu,du);
  105.      end; %if ~isempty(au)
  106.   end; %if opt=='d'
  107. end;%if krp>n
  108. %
  109. % Copyright MUSYN INC 1991,  All Rights Reserved
  110.