www.gusucode.com > rctutil 工具箱 matlab源码程序 > rctutil/private/mufastub.m

    function [ub,Dr,Dc,Gcr] = mufastub(M,index,bnds)
% function [ub,Dr,Dc,Gcr] = mufastub(M,index,bnds);
%
% 12/23/2005 PJS
% Based on deebal, geestep and rub code by Young/Newlin

% Copyright 2011-2014 The MathWorks, Inc.

% Get problem data
[Mr,Mc]=size(M);
realr = index.allreal.allrows;
realc = index.allreal.allcols;
%lreal = length(realr);
G_mask = index.masks.G_mask;

% Initialize bisection bounds
if nargin<3
    ub = max( [max(svd(M)), 10*eps] );
    lb = 0;
else
    ub = bnds(1);
    lb = bnds(2);
end
ubhigh = ub;
ublow = max(lb,1e-3);  % XXX Nov 21 2013: 10*eps -> 1e-3

%-------------------------------------------------------------------
% Choose G to cancel skew-Hermitian part of M in the balanced form
% Diagonalize this choice by wrapping orthog. transformation back
% into D-scales.
%-------------------------------------------------------------------
G = M(realr,realc).*G_mask;
G = (G - G')/2i; 
[gu,gdiag]=schur(G);  
gdiag = diag(gdiag);

dMd = M;
dMd(realr,:) = gu'*dMd(realr,:);
dMd(:,realc) = dMd(:,realc)*gu;

% Bisect in balanced form to find best (i.e. smallest)
% upper bound which satisfies
%         norm(Lterm*(Cterm/ub)*Rterm)<1
Lterm = ones(Mr,1);
Rterm = ones(Mc,1);
Cterm = dMd;
Cterm(realr,realc) = Cterm(realr,realc)-diag(1i*gdiag);
gdiag2 = gdiag.^2;
while (ubhigh-ublow)> 0.001*ublow  
    ubtry = (ublow+ubhigh)/2;
    gtry2 = gdiag2/(ubtry.^2);
    Lterm(realr) = (1+gtry2).^(-0.25);
    Rterm(realc) = Lterm(realr);
    if norm((Cterm/ubtry).*(Lterm*Rterm')) < 1        
        ubhigh=ubtry;
    else
        ublow =ubtry;
    end
end

% Store best g-scale obtained in balanced form
if ubhigh<ub
    ub = ubhigh;
    g	= gdiag/ub;
else
    g	= 0*gdiag;
end

%-------------------------------------------------------------------
% Do a single gradient step on g in balanced form.  
%-------------------------------------------------------------------
Lterm(realr) = (1+g.^2).^(-0.25);
Rterm(realc) = Lterm(realr);
Cterm = dMd/ub;
Cterm(realr,realc) = Cterm(realr,realc)-diag(1i*g);
[uu,ss,vv]=svd(Cterm.*(Lterm*Rterm'));
smax	= max([max(diag(ss)) 10*eps]);
sidx	= find((diag(ss)/smax)>0.95); % find s.v. clustered near max
if isempty(sidx)
    % This can happen if ss=0, e.g. for mussv(5*i,[-1 0])
    sidx = 1;
end    
lidx = length(sidx); 
uu    = uu(realr,sidx);
ss   = ss(sidx,sidx);
vv   = vv(realc,sidx);

% Compute ds_dg, the gradient of largest singular vals w.r.t. g
gtmp1 = repmat( g./(1+g.^2)/2 , [1 lidx]);
gtmp2 = repmat( (1+g.^2).^(-0.5) , [1 lidx]);
ds_dg = -real( ((abs(uu).^2 + abs(vv).^2).*gtmp1)*ss ...
    + 1i*conj(uu).*vv.*gtmp2 );

% Let g(t)=g0+t*gdir and compute the descent direction, gdir, to give a
% desired derivative of the largest singular values w.r.t to t, ds_dt.
ds_dt_des = -10*(diag(ss)-0.9);
gdir = pinv(ds_dg')*ds_dt_des;
ds_dt = ds_dg(:,1)'*gdir;

% Compute ds_dub, the negative of the gradient of largest 
% singular val w.r.t. ub. Then compute dub_dt, the gradient of the 
% upper bound w.r.t. t. Compute initial step, t0, to reduce upper
% bound by dub (based on first order Taylor series).
ds_dub = real( smax + 1i*uu(:,1)'*(gtmp2(:,1).*vv(:,1).*g) )/ub;
dub_dt = ds_dt/ds_dub;
dub = (lb - ub)/2;
t0 = dub/dub_dt;
t0 = min( max(t0,-1e6) , 1e6 ); % ensure t0 is bounded

% Take initial step and then backtrack, if necessary. 
gt = g+t0*gdir;  
Lterm(realr) = (1+gt.^2).^(-0.25);
Rterm(realc) = Lterm(realr);
Cterm = dMd/ub;
Cterm(realr,realc) = Cterm(realr,realc)-diag(1i*gt);
st=norm(Cterm.*(Lterm*Rterm'));
cnt = 0;
while st>smax && (cnt<10)
    t0 = t0/2;
    gt = g+t0*gdir;  
    Lterm(realr) = (1+gt.^2).^(-0.25);
    Rterm(realc) = Lterm(realr);
    Cterm = dMd/ub;
    Cterm(realr,realc) = Cterm(realr,realc)-diag(1i*gt);
    st=norm(Cterm.*(Lterm*Rterm'));
    cnt = cnt+1;
end

% Update g if smax is reduced.
if st<smax
    g = gt;
end

%-------------------------------------------------------------------
% Convert from Balanced to LMI form
%-------------------------------------------------------------------
gbal = (1+g.*g).^(-0.5);
Dtmp = gu*diag(gbal)*gu';

Dr = eye(Mr);
Dr(realr,realr) = Dtmp;

Dc = eye(Mc);
Dc(realc,realc) = Dtmp;

Gcr= zeros(Mc,Mr);
Gcr(realc,realr)=ub*(gu*diag(g.*gbal)*gu');

% FV certification
if ~isempty(index.FVidx.fixedBlkIdx)
   FixedCols = index.FVidx.FixedCols;
   VaryCols = index.FVidx.VaryCols;
   DcF = zeros(Mc);
   DcV = zeros(Mc);
   DcF(FixedCols,FixedCols) = Dc(FixedCols,FixedCols);
   DcV(VaryCols,VaryCols) = Dc(VaryCols,VaryCols);
   GcrM = Gcr*M;
   A = 1i*(GcrM - GcrM')+M'*Dr*M-DcF;
   B = DcV;
   evals = eig(A, B);
   finiteIdx = ~isinf(evals);
   ubsq = 1.0001*max(real(evals(finiteIdx)));
   if max(real(eig(A-ubsq*B)))<=0
      ub = sqrt( max(0,ubsq) );
   else
      ub = inf;
   end
else
   % Re-compute upper bound based on updated g-scale
   ubsq = max(real(eig( M'*Dr*M+1i*(Gcr*M-M'*Gcr') , Dc )));
   ub = sqrt( max(0,ubsq) );
end