home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / BALSQ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.8 KB  |  156 lines
  1. function [aq,bq,cq,dq,aug,hsv,slbig,srbig] = balsq(Z1,Z2,Z3,Z4,Z5,Z6)
  2. % [SS_Q,AUG,HSV,SLBIG,SRBIG] = BALSQ(SS_,TYPE,NO) or
  3. % [AQ,BQ,CQ,DQ,AUG,HSV,SLBIG,SRBIG] = BALSQ(A,B,C,D,TYPE,NO) performs 
  4. %    balanced-truncation model reduction on G(s):=(a,b,c,d) via the 
  5. %    square-root of PQ such that the infinity-norm of the error 
  6. %    (Ghed(s) - G(s))<=sum of the 2(n-k) smaller Hankel singular values. 
  7. %
  8. %           (aq,bq,cq,dq) = (slbig'*a*srbig,slbig'*b,c*srbig,d)
  9. %
  10. %   Based on the "TYPE" selected, you have the following options:
  11. %
  12. %    1). TYPE = 1  --- no: order "k" of the reduced order model.
  13. %
  14. %    2). TYPE = 2  --- find a k-th order model such that the total error
  15. %                      is less than "no"
  16. %
  17. %    3). TYPE = 3  --- display Hankel SV and prompt for "k".
  18. %
  19. %    AUG = [No. of state(s) removed, Error bound], HSV = Hankel SV's.
  20.  
  21. % R. Y. Chiang & M. G. Safonov 9/11/87
  22. % Copyright (c) 1988 by the MathWorks, Inc.
  23. % All Rights Reserved.
  24.  
  25. inargs = '(a,b,c,d,Type,no)';
  26. eval(mkargs(inargs,nargin,'ss'))
  27.  
  28. [ma,na] = size(a);
  29. [md,nd] = size(d);
  30. [hsv,p,q] = hksv(a,b,c);
  31. %
  32. % ------ Model reduction based on your choice of TYPE:
  33. %
  34. if Type == 1
  35.    kk = no+1;
  36. end
  37. %
  38. if Type == 2
  39.    tails = 0;
  40.    kk = 1;
  41.    for i = ma:-1:1
  42.        tails = tails + hsv(i);
  43.        if 2*tails > no
  44.           kk = i + 1;
  45.           break
  46.        end
  47.    end
  48. end
  49. %
  50. if Type == 3
  51.    format short e
  52.    format compact
  53.    [mhsv,nhsv] = size(hsv);
  54.    if mhsv < 60
  55.       disp('    Hankel Singular Values:')
  56.       hsv'
  57.       for i = 1:mhsv
  58.         if hsv(i) == 0
  59.            hsvp(i) = eps;
  60.         else
  61.            hsvp(i) = hsv(i);
  62.         end
  63.       end
  64.       disp('                              (strike a key to see the plot ...)')
  65.       pause
  66.       subplot
  67.       plot(20*log10(hsvp),'--');hold
  68.       plot(20*log10(hsvp),'*');grid;hold
  69.       ylabel('DB')
  70.       title('Hankel Singular Values')
  71.       pause
  72.       no = input('Please assign the k-th index for k-th order model reduction:');
  73.    else
  74.       disp('    Hankel Singular Values:')
  75.       hsv(1:60,:)'
  76.       disp('                              (strike a key for more ...)')
  77.       pause
  78.       hsv(61:mhsv,:)'
  79.       for i = 1:mhsv
  80.         if hsv(i) == 0
  81.            hsvp(i) = eps;
  82.         else
  83.            hsvp(i) = hsv(i);
  84.         end
  85.       end
  86.       disp('                              (strike a key to see the plot ...)')
  87.       pause
  88.       subplot
  89.       plot(20*log10(hsvp),'--');hold
  90.       plot(20*log10(hsvp),'*');grid;hold
  91.       ylabel('DB')
  92.       title('Hankel Singular Values')
  93.       no = input('Please assign the k-th index for k-th order model reduction:');
  94.    end
  95.    format loose
  96.    kk = no + 1;
  97. end
  98. %
  99. % ------ Save all the states:
  100. %
  101. if kk > na 
  102.    aq= a; bq= b; cq= c; dq= d;
  103.    slbig = eye(na); srbig = eye(na);
  104.    aug = [0 0];
  105.    return 
  106. end
  107. %
  108. % ------ Disgard all the states:
  109. %
  110. if kk == 1
  111.    aq= zeros(ma,na);   bq= zeros(ma,nd);
  112.    cq= zeros(md,na);   dq= zeros(md,nd);
  113.    slbig = []; srbig = [];
  114.    bnd = 2*sum(hsv); 
  115.    aug = [na bnd];
  116.    return
  117. end 
  118. %
  119. % ------ k-th order Schur balanced model reduction:
  120. %
  121. bnd = 2*sum(hsv(kk:na));
  122. strm = na-kk+1;
  123. aug = [strm bnd];
  124. %
  125. % ----- Find the square root for the basis of left/right eigenspace :
  126. %
  127. [up,sp,vp] = svd(p);
  128. lr = up*diag(sqrt(diag(sp)));
  129. [uq,sq,vq] = svd(q);
  130. lo = uq*diag(sqrt(diag(sq)));
  131. %
  132. [uc,sc,vc] = svd(lo'*lr);
  133. %
  134. % ------ Find the similarity transformation :
  135. %
  136. s1 = sc(1:(kk-1),1:(kk-1));
  137. scih = diag(ones(kk-1,1)./sqrt(diag(s1)));
  138. uc = uc(:,1:(kk-1));
  139. vc = vc(:,1:(kk-1));
  140. slbig = lo*uc*scih;
  141. srbig = lr*vc*scih;
  142. %
  143. aq = slbig'*a*srbig;
  144. bq = slbig'*b;
  145. cq = c*srbig;
  146. dq = d;
  147. %
  148. if xsflag
  149.   aq = mksys(aq,bq,cq,dq);
  150.   bq = totbnd;
  151.   cq = hsv;
  152.   dq = slbig;
  153.   aug = srbig;
  154. end
  155. %
  156. % ------ End of BALSQ.M --- RYC/MGS 9/11/87 %