home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / RESCHMR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.6 KB  |  178 lines
  1. function [ahed,bhed,ched,dhed,aug,hsv] = reschmr(A,B,C,D,Type,no,info)
  2. %
  3. % [AHED,BHED,CHED,DHED,AUG,HSV]=RESCHMR(A,B,C,D,TYPE,NO,INFO) performs 
  4. %    relative error Schur model reduction on a SQRARE, STABLE G(s):=
  5. %    (A,B,C,D). The infinity-norm of the relative error is bounded as
  6. %              -1                          __
  7. %            |G (s)(Gm(s)-G(s))|    <= (   ||  [(1+si)/(1-si)] ) - 1,
  8. %                               inf     k+1 to n
  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, the problem is ill-posed.
  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. [rd,cd] = size(D);
  30. if nargin == 6
  31.    info = 'right';
  32. end
  33. %
  34. if (info == 'left ') | (rd < cd)
  35.    A = A'; temp = B; B = C'; C = temp'; D = D';
  36. end
  37. %
  38. if rank(D*D') < min([rd, cd])
  39.   disp('   WARNING: D MATRIX IS NOT FULL RANK - - -');
  40.   disp('            THE PROBLEM IS ILL-CONDITIONED !');
  41.   return
  42. end
  43. %
  44. disp(' ')
  45. disp('        - - Working on Schur BST-REM model reduction - -');
  46. [ma,na] = size(A);
  47. Flag = exist('xxxxx0');
  48. P = gram(A,B);
  49. G = P*C'+B*D';
  50. phi = D*D';
  51. %
  52. qrnQ = [zeros(ma) -C';-C -phi];
  53. [K,Q,Qerr] = lqrc(A,G,qrnQ);
  54. %
  55. lambda = eig(P*Q);
  56. hsv = sqrt(lambda);
  57. [hsv,index] = sort(real(hsv));
  58. hsv = hsv(ma:-1:1,:);
  59. %
  60. % ------ Model reduction based on your choice of TYPE:
  61. %
  62. if Type == 1
  63.    kk = no+1;
  64. end
  65. %
  66. if Type == 2
  67.    toldb = no;    %in db
  68.    tolabs = 10^(toldb/20);
  69.    resig = 1;
  70.    for i = ma:-1:1
  71.       resig = resig*((1+hsv(i))/(1-hsv(i)));
  72.       rem = resig -1; 
  73.       if rem > tolabs
  74.          kk = i+1;
  75.          break
  76.       end
  77.    end
  78. end
  79. %
  80. if Type == 3
  81.    format short e
  82.    format compact
  83.    [mhsv,nhsv] = size(hsv);
  84.    if mhsv < 60
  85.       disp('    Hankel Singular Values:')
  86.       hsv'
  87.       kk = input('Please assign the k-th index for k-th order model reduction:');
  88.    else
  89.       disp('    Hankel Singular Values:')
  90.       hsv(1:60,:)'
  91.       disp('                              (strike a key for more ...)')
  92.       pause
  93.       hsv(61:mhsv,:)'
  94.       kk = input('Please assign the k-th index for k-th order model reduction:');
  95.    end
  96.    format loose
  97.    kk = kk + 1;
  98. end
  99. %
  100. % ------ Save all the states:
  101. %
  102. if kk > na 
  103.    ahed = a; bhed = b; ched = c; dhed = d;
  104.    aug = [0 0];
  105.    return 
  106. end
  107. %
  108. % ------ Disgard all the states:
  109. %
  110. if kk < 1 
  111.    ahed = zeros(ma,na);   bhed = zeros(ma,nd);
  112.    ched = zeros(md,na);   dhed = zeros(md,nd);
  113.    resig = 1;
  114.    for i = 1:na
  115.       resig = resig*(1+hsv(i))/(1-hsv(i));
  116.    end
  117.    rem = resig-1;
  118.    aug = [na rem];
  119.    return
  120. end 
  121. %
  122. % ------ Computing the relative error bound:
  123. %
  124. resig = 1;
  125. for i = kk:na
  126.     resig = resig*(1+hsv(i))/(1-hsv(i));
  127. end
  128. rem = resig-1;
  129. strm = na-kk+1;
  130. aug = [strm rem];
  131. % 
  132. % ---------------------- Start Model Reduction:
  133. %
  134. % ------ Find the left-eigenspace basis:
  135. %
  136. ro = (hsv(kk-1)^2+hsv(kk)^2)/2.;
  137. gammaa = P*Q-ro*eye(na);
  138. if Flag == 2
  139.    [va,ta] = rschur(gammaa,1);            % FORTRAN code
  140. else
  141.    [va,ta,msa,swpa] = hqr10(gammaa);
  142. end
  143. vlbig = va(:,(na-kk+2):na);
  144. %
  145. % ------ Find the right-eigenspace basis:
  146. %
  147. gammad = -gammaa;
  148. if Flag == 2
  149.    [vd,td] = rschur(gammad,1);           % FORTRAN code
  150. else
  151.    [vd,td,msd,swpd] = hqr10(gammad);
  152. end
  153. vrbig = vd(:,1:kk-1);
  154. %
  155. % ------ Find the similarity transformation:
  156. %
  157. ee = vlbig'*vrbig;
  158. [ue,se,ve] = svd(ee);
  159. %
  160. seih = diag(ones(kk-1,1)./sqrt(diag(se)));
  161. slbig = vlbig*ue*seih;
  162. srbig = vrbig*ve*seih;
  163. %
  164. ahed = slbig'*A*srbig;
  165. bhed = slbig'*B;
  166. ched = C*srbig;
  167. dhed = D;
  168. %
  169. if (info == 'left ') | (rd < cd)
  170.   ahed = ahed'; temp = bhed; 
  171.   bhed = ched'; ched = temp'; dhed = dhed';
  172. end
  173. %
  174. disp(' ')
  175. disp(['               ' int2str(aug(1,1)), '    states removed !!'])
  176. %
  177. % ------ End of RESCHMR.M --- RYC/MGS %
  178.