home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS2.DI$ / DCALC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.7 KB  |  63 lines
  1. % function  d = dcalc(b,c,sig,q)
  2. %
  3. % Calculates a D matrix such that the L-infinity norm of
  4. %  (C(sI-A)^(-1)B+D) is <= sig(1)+...+sig(n).
  5. %  B and C are from a balanced realisation with n states,
  6. %  SIG is the vector of hankel singular values.
  7. %  Q is the dimension in which the system is embedded.
  8.  
  9. function d = dcalc(b,c,sig,q)
  10.  
  11. n=sum(sign(sig));if n<1,error('n<1');end;
  12. [nn,m]=size(b);[p,nn]=size(c);
  13. if nargin<4, q=m+p;end
  14. if q<m+p,q=m+p;end;
  15. d=zeros(p,m);
  16. z=[b';zeros(q-m,n)];
  17. y=[c;zeros(q-p,n)];
  18. signd=-1;u=zeros(q,q);
  19. oneq=ones(q,1);
  20. i=0;r=1;
  21. while i+r<=n,
  22.   i=i+r;
  23.  % test for repeated sig(i)
  24.   r=1;
  25.   if i<n,
  26.     while (sig(i+r-1)<1.001*sig(i+r)), 
  27.       r=r+1;
  28.       if i+r>n, break; end;
  29.      end; %while (sig(i)<
  30.    end; %if i<n
  31.   if r==1,
  32.     a1=norm(y(:,i))*sign(y(1,i)+eps);
  33.     y(1,i)=y(1,i)+a1;
  34.     pi1=a1*y(1,i);
  35.     a2=norm(z(:,i))*sign(z(1,i)+eps);
  36.     z(1,i)=z(1,i)+a2;
  37.     pi2=a2*z(1,i);
  38.     bet=-a1/a2;
  39.     u=eye(q)-y(:,i)*(y(:,i))'/pi1;
  40.     u(:,1)=u(:,1)*bet;
  41.     u(:,2:m+p-1)=u(:,[p+1:m+p-1,2:p]);
  42.     u=u-u*z(:,i)*(z(:,i))'/pi2;
  43.   else, %if r=1
  44.     [uc,sc,vc]=svd(y(:,i:i+r-1));
  45.     [ub,sb,vb]=svd(z(:,i:i+r-1)*vc);
  46.     rb=rank(z(:,i:i+r-1));
  47.     u=-uc*[vb(1:rb,1:rb) zeros(rb,q-rb); zeros(q-rb,rb) eye(q-rb)]*ub';
  48.   end;%if r=1
  49. % update z & y
  50.   if i~=n-r+1,
  51.     gam=oneq*sqrt(sig(i)^2*ones(1,n-i-r+1)-(sig(i+r:n))'.^2);
  52.     tmp=-(y(:,i+r:n).*(oneq*(sig(i+r:n))')+u*z(:,i+r:n)*sig(i))./gam;
  53.     z(:,i+r:n)=(z(:,i+r:n).*(oneq*(sig(i+r:n))')+u'*y(:,i+r:n)*sig(i))./gam;
  54.     y(:,i+r:n)=tmp;
  55.   end; %if i~=n-r+1
  56.   signd=-signd;
  57.   d=d+signd*sig(i)*u(1:p,1:m);
  58. end;%while i<n
  59.  
  60.  
  61. %
  62. % Copyright MUSYN INC 1991,  All Rights Reserved
  63.