home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / BSTSCHMR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.5 KB  |  212 lines
  1. function [ahed,bhed,ched,dhed,aug,hsv] = bstschmr(Z1,Z2,Z3,Z4,Z5,Z6,Z7)
  2. % [SS_H,AUG,HSV] = BSTSCHMR(SS_,TYPE,NO,INFO) or
  3. % [AHED,BHED,CHED,DHED,AUG,HSV]=BSTSCHMR(A,B,C,D,TYPE,NO,INFO) performs 
  4. %    relative error Schur model reduction on a SQUARE, STABLE G(s):=
  5. %    (A,B,C,D). The infinity-norm of the relative error is bounded as
  6. %               -1                             n
  7. %            |Gm (s)(Gm(s)-G(s))|    <= 2 * ( SUM  si / (1 - si) )
  8. %                                inf          k+1
  9. %    where si denotes the i-th Hankel singular value of the all-pass 
  10. %    "phase matrix" of G(s). 
  11. %    The algorithm is based on the Balanced Stochastic Truncation (BST) 
  12. %    theory with Relative Error Method (REM).
  13. %    Based on the "TYPE" selected, you have the following options:
  14. %     1). TYPE = 1  --- no: size "k" of the reduced order model.
  15. %     2). TYPE = 2  --- find k-th order reduced model that 
  16. %                       tolerance (db) <= "no".
  17. %     3). TYPE = 3  --- display all the Hankel SV of phase matrix and 
  18. %                   prompt for "k" (in this case, no need to specify "no").
  19. %    Input variable: "info" = 'left '; or 'right' (default value)
  20. %    Output variable "aug": aug(1,1) = no. of state removed
  21. %                           aug(1,2) = relative error bound
  22. %    Note that if D is not full rank, an error will result.
  23.  
  24. % R. Y. Chiang & M. G. Safonov 2/30/88
  25. % Copyright (c) 1988 by the MathWorks, Inc.
  26. % All Rights Reserved.
  27. % ---------------------------------------------------------------------------
  28.  
  29. nag1 = nargin;
  30. inargs = '(A,B,C,D,Type,no,info)';
  31. eval(mkargs(inargs,nag1,'ss'))
  32. nag1 = nargin;     % NARGIN may have changed
  33.  
  34. [rd,cd] = size(D);
  35. if nag1 == 6
  36.    info = 'right';
  37. end
  38. %
  39. if (info == 'left ') | (rd < cd)
  40.    A = A'; temp = B; B = C'; C = temp'; D = D';
  41. end
  42. %
  43. if rank(D*D') < min([rd, cd])
  44.   disp('   WARNING: D MATRIX IS NOT FULL RANK - - -');
  45.   disp('            THE PROBLEM IS ILL-CONDITIONED !');
  46.   return
  47. end
  48. %
  49. disp(' ')
  50. disp('        - - Working on Schur BST-REM model reduction - -');
  51. [ma,na] = size(A);
  52. P = gram(A,B);
  53. G = P*C'+B*D';
  54. phi = D*D';
  55. %
  56. qrnQ = [zeros(ma) -C';-C -phi];
  57. [K,Q,Qerr] = lqrc(A,G,qrnQ);
  58. %
  59. lambda = eig(P*Q);
  60. hsv = sqrt(lambda);
  61. [hsv,index] = sort(real(hsv));
  62. hsv = hsv(ma:-1:1,:);
  63. %
  64. % ------ Model reduction based on your choice of TYPE:
  65. %
  66. if Type == 1
  67.    kk = no+1;
  68. end
  69. %
  70. if Type == 2
  71.    toldb = no;    %in db
  72.    tolabs = 10^(toldb/20);
  73.    resig = 0;
  74.    for i = ma:-1:1
  75.       resig = resig + hsv(i)/(1-hsv(i));
  76.       rem = 2*resig;
  77.       if rem > tolabs
  78.          kk = i+1;
  79.          break
  80.       end
  81.    end
  82. end
  83. %
  84. if Type == 3
  85.    format short e
  86.    format compact
  87.    [mhsv,nhsv] = size(hsv);
  88.    if mhsv < 60
  89.       disp('    Hankel Singular Values:')
  90.       hsv'
  91.       for i = 1:mhsv
  92.         if hsv(i) == 0
  93.            hsvp(i) = eps;
  94.         else
  95.            hsvp(i) = hsv(i);
  96.         end
  97.       end
  98.       disp(' ')
  99.       disp('                              (strike a key to see the plot ...)')
  100.       pause
  101.       subplot
  102.       plot(20*log10(hsvp),'--');hold
  103.       plot(20*log10(hsvp),'*');grid;hold
  104.       ylabel('DB')
  105.       title('Hankel Singular Values')
  106.       pause
  107.       no = input('Please assign the k-th index for k-th order model reduction: ');
  108.    else
  109.       disp('    Hankel Singular Values:')
  110.       hsv(1:60,:)'
  111.       disp('                              (strike a key for more ...)')
  112.       pause
  113.       hsv(61:mhsv,:)'
  114.       for i = 1:mhsv
  115.         if hsv(i) == 0
  116.            hsvp(i) = eps;
  117.         else
  118.            hsvp(i) = hsv(i);
  119.         end
  120.       end
  121.       disp(' ')
  122.       disp('                              (strike a key to see the plot ...)')
  123.       pause
  124.       subplot
  125.       plot(20*log10(hsvp),'--');hold
  126.       plot(20*log10(hsvp),'*');grid;hold
  127.       ylabel('DB')
  128.       title('Hankel Singular Values')
  129.       pause
  130.       no = input('Please assign the k-th index for k-th order model reduction: ');
  131.    end
  132.    format loose
  133.    kk = no + 1;
  134. end
  135. %
  136. % ------ Save all the states:
  137. %
  138. if kk > na 
  139.    ahed = a; bhed = b; ched = c; dhed = d;
  140.    aug = [0 0];
  141.    return 
  142. end
  143. %
  144. % ------ Disgard all the states:
  145. %
  146. if kk < 1 
  147.    ahed = zeros(ma,na);   bhed = zeros(ma,nd);
  148.    ched = zeros(md,na);   dhed = zeros(md,nd);
  149.    resig = 0;
  150.    for i = 1:na
  151.       resig = resig + hsv(i)/(1-hsv(i));
  152.    end
  153.    rem = 2*resig;
  154.    aug = [na rem];
  155.    return
  156. end 
  157. %
  158. % ------ Computing the relative error bound:
  159. %
  160. resig = 0;
  161. for i = kk:na
  162.     resig = resig + hsv(i)/(1-hsv(i));
  163. end
  164. rem = 2*resig;
  165. strm = na-kk+1;
  166. aug = [strm rem];
  167. % 
  168. % ---------------------- Start Model Reduction:
  169. %
  170. % ------ Find the left-eigenspace basis:
  171. %
  172. ro = (hsv(kk-1)^2+hsv(kk)^2)/2.;
  173. gammaa = P*Q-ro*eye(na);
  174. [va,ta,msa] = blkrsch(gammaa,1,na-kk+1);
  175. vlbig = va(:,(na-kk+2):na);
  176. %
  177. % ------ Find the right-eigenspace basis:
  178. %
  179. gammad = -gammaa;
  180. [vd,td,msd] = blkrsch(gammad,1,kk-1);
  181. vrbig = vd(:,1:kk-1);
  182. %
  183. % ------ Find the similarity transformation:
  184. %
  185. ee = vlbig'*vrbig;
  186. [ue,se,ve] = svd(ee);
  187. %
  188. seih = diag(ones(kk-1,1)./sqrt(diag(se)));
  189. slbig = vlbig*ue*seih;
  190. srbig = vrbig*ve*seih;
  191. %
  192. ahed = slbig'*A*srbig;
  193. bhed = slbig'*B;
  194. ched = C*srbig;
  195. dhed = D;
  196. %
  197. if (info == 'left ') | (rd < cd)
  198.   ahed = ahed'; temp = bhed; 
  199.   bhed = ched'; ched = temp'; dhed = dhed';
  200. end
  201. %
  202. if xsflag
  203.    ahed = mksys(ahed,bhed,ched,dhed);
  204.    bhed = aug;
  205.    ched = hsv;
  206. end
  207. %
  208. disp(' ')
  209. disp(['               ' int2str(aug(1,1)), '    states removed !!'])
  210. %
  211. % ------ End of BSTSCHMR.M --- RYC/MGS %
  212.