home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / MIXEDALG.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.5 KB  |  193 lines
  1. function [mult,xqo] = mixedalg(Tt,xqi,K);
  2. % [MULT,XQO] = MIXEDALG(T,XQI,K) finds the existence of a complex,diagonal 
  3. % multiplier, M, with positive real parts such that
  4. %
  5. %           Re ( Tt.M ) > 0                   .....   ($)
  6. %
  7. %  inputs : Tt   -- a square complex matrix.
  8. %           xqi  -- initial guess for any vector of unit norm.
  9. %            K   -- a column vector defining the type of uncertainties 
  10. %                   in the diagonal  uncertainty block
  11. %                -- 1 : denotes real parametric uncertainties normalised
  12. %                       in sector[-1,1]
  13. %                   2 : denotes complex parametric uncertainties normalise 
  14. %                       in unit disk
  15. % outputs : mult -- the multiplier vector obtained such that ($) holds.
  16. %                   If ($) do not holds, mult is the zero vector.
  17. %           xqo  -- the final vector x of unit norm such that ($) holds
  18. %
  19.  
  20. % Authors: Penghin Lee and M. G. Safonov 7/92
  21. %
  22. % R. Y. Chiang & M. G. Safonov 
  23. % Copyright (c) 1988-92 by the MathWorks, Inc.
  24. % All Rights Reserved.
  25.  
  26. n = max(size(Tt)); qc = 1; zc = 1e-12; nzc = 1e-12; zcriteria = 5;j = sqrt(-1);
  27. qc_check = 150;
  28.  
  29. % Scaling
  30. abs_max_Tt_ele = abs ( max(max(Tt)) );
  31. abs_min_Tt_ele = abs ( min(min(Tt)) );
  32.  if abs_max_Tt_ele < 1e-01 | abs_min_Tt_ele < 1e-01;
  33.   scaling_factor = 1e+04;
  34.   else scaling_factor = 1;
  35.  end;
  36.  
  37. Tt = scaling_factor * Tt;
  38.  
  39. % Finding the Q(i) matrix for the given Tt matrix
  40. for i = 1:n;
  41.     for m = 1:n;
  42.         if m == i;
  43.          ti(m,:) = Tt(i,:);
  44.         else 
  45.             for l = 1:n;
  46.              ti(m,l) = 0;
  47.             end;
  48.         end;
  49.     end;
  50.  q(1:n,((((i-1)*n)+1):(i*n))) = j * (ti - ti');
  51. end;
  52.  
  53. for i = (n+1):2*n;
  54.     for m = 1:n;
  55.         if m == i - n;
  56.          ti(m,:) = Tt(m,:);
  57.         else
  58.             for l = 1:n;
  59.              ti(m,l) = 0;
  60.             end;
  61.         end;
  62.     end;
  63.  q(1:n,((((i-1)*n)+1):(i*n))) = ti + ti';
  64. end;
  65.  
  66. for i = 1:2*n;
  67.  zqir(i) = xqi' * q(1:n,((((i-1)*n)+1):(i*n))) * xqi;
  68. end;
  69. zq = real(zqir');
  70.  
  71. % To ensure initial guess is in the convex hull of the set.
  72. for i = n+1:2*n;
  73.   if zq(i,1) < 0;
  74.    zq(i,1) = 0;
  75.   end;
  76. end;
  77.  
  78. % To take care of mixed real and/or complex parametric uncertainty
  79. complx_type = find( K == 2 );
  80. complx_type_size = max( size(complx_type) );
  81.  
  82. if complx_type_size > 0;
  83.   for i = 1 : complx_type_size;
  84.    zq_zero = complx_type(i,1);
  85.    zq(zq_zero,1) = 0;
  86.   end;
  87. end;
  88.  
  89. % Begin iterations
  90. while zcriteria > zc
  91.  
  92. % Finding the Q(z) matrix for a z and the Q(i)
  93. qzm = zeros(n);
  94.   for i = 1:2*n;
  95.    qzi = zq(i,1) * q(1:n,((((i-1)*n)+1):(i*n))); 
  96.    qzm = qzi + qzm;
  97.   end;
  98.  
  99. % To find the smallest eigenvalue and the associated eigenvector
  100. [v,lambda] = eig(qzm,'nobalance');
  101.  for i = 1:n;
  102.   lambx(i,1) = real((lambda(i,i)));
  103.  end;
  104. lambmin = min(lambx);
  105.  
  106.  for i = 1:n;
  107.      if lambmin == lambx(i,1);
  108.       mi=i;
  109.      end;
  110.  end;
  111.  
  112.   for i = 1:n;
  113.    xq(i,1) = v(i,mi);
  114.   end;
  115. lam(qc,1) = lambmin;
  116.  
  117. % For computing f(xq) from the computed xq 
  118.  for i = 1:2*n;
  119.   fxqi(i) = xq' * q(1:n,((((i-1)*n)+1):(i*n))) * xq;
  120.  end;
  121. fxq = real(fxqi');% The new point is fxq.
  122.  
  123. %     To compute the optimum value of alpha in
  124. %     min  ||zq1 = alpha*zq  + (1-alpha)fxq ||;
  125. %     for all alpha between 0 and 1
  126.  
  127. if (zq-fxq) == 0; 
  128.  alphaopt = 1;
  129.   else alphaopt = -((zq-fxq)'*fxq) / ((zq-fxq)'*(zq-fxq));
  130.   if alphaopt < 0;
  131.    alphaopt = 0;
  132.   elseif alphaopt > 1;
  133.    alphaopt = 1;
  134.   end;
  135. end; 
  136.  
  137. % The new point zq1
  138. zq1 = alphaopt * zq + (1-alphaopt) * fxq;
  139.  
  140. % To handle restrictions on  zq
  141. for i = n+1:2*n;
  142.    if zq1(i,1) < 0;
  143.     zq1(i,1) = 0;
  144.    end;
  145. end;
  146.  
  147. % To take care of mixed real and/or complex parametric uncertainty
  148. if complx_type_size > 0;
  149.   for i = 1 : complx_type_size;
  150.    zq_zero = complx_type(i,1);
  151.    zq1(zq_zero,1) = 0;
  152.   end;
  153. end;
  154.  
  155. % For convergence of zq
  156.  
  157. if qc > 3; % In case the initial guess has a very small norm
  158.  
  159. zcriteria = norm(zq1-zq) / norm(zq);
  160.  
  161.  if zcriteria > zc & (norm(zq1)/scaling_factor) < nzc;
  162.   zq1 = zeros(2*n,1);
  163.   zcriteria = zc - 1;
  164.  end;
  165.  
  166. end;
  167.  
  168. if qc > 10;
  169.  
  170.   if lam(qc,1) > 0 & lam(qc-1,1) > 0;
  171.    zcriteria = zc - 1;
  172.   end;
  173.  
  174. if qc == 248 & lam(qc,1) <= 0 & lam(qc,1) >= -30;
  175.     qc_check = 300;
  176. end;
  177.  
  178.         if (lam(qc,1)<=0 & qc > qc_check);
  179.          zq1 = zeros(2*n,1); 
  180.          zcriteria = zc - 1 ;
  181.         end;
  182.  
  183. end; %if qc>10
  184.  
  185. qc = qc + 1;
  186. zq = real(zq1);nzqq(qc,1) = norm(zq);
  187.  
  188. end; %while z
  189.  
  190. xqo = xq;
  191. mult = zq / scaling_factor;
  192. %
  193. % ----------------- End of MIXEDALG.M %PHL/MGS