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