vlpvblnd.m

线性时变系统控制器设计的工具包

function [eigx,eigy,eigxy,fail] = vlpvblnd(vlpv,nmeas,nctrl,...
gam,xmat1,ymat1,xmat2,ymat2,vnu,fparm,gparm,gradf,gradg,bparm,gradb)

% [eigx,eigy,eigxy,fail] = vlpvblnd(vlpv,nmeas,nctrl,gam,...
% xmat1,ymat1,xmat2,ymat2,vnu,fparm,gparm,gradf,gradg,bparm,gradb)
% checks "blended" PDLF LMI solutions in the intersection of two adjacent 
% parameter regions.
%
% Inputs: 
% vlpv 		varying matrix of the lpv plant at parm values
% nmeas,nctrl	number of measurements & controls, respectively
% gam 		the largest of the two optimal gammas
% xmat1,xmat2 	VARYING optimal X's from solver
% ymat1,ymat2 	VARYING optimal Y's from solver
% vnu 		(VARYING parameter-dependent) parameter rate bound vector
%		| parm_i | <= nu_i for i = 1,...,s
%		nu can also be two-sided (2 cols) or a rate grid (> 2 cols)
% fparm,gparm   VARYING basis-function values at grid points
% gradf,gradg   VARYING basis-function partial derivatives at grid points
% 		If X or Y is constant, use [] for fparm/gradf or gparm/gradg.
% bparm		VARYING blending function values at grid points.
%		bparm = 1 outside region #2, bparm = 0 outside region #1.
% gradb		VARYING blending function partial derivatives
%
% Outputs: 
% eigx		max eigenvalues of output injection matrices (X)
% eigy		max eigenvalues of state feedback matrices (Y)
% eigxy		min eigenvalues of spectral radius matrices (XY)
% fail		index of LMIs that failed (column 1 = grid point #, 
% 	column 2 = vertex # of rate polytope, column 3 = 0/1/2 for XY/X/Y)

if nargin == 0
 disp('[eigx,eigy,eigxy,fail] = vlpvblnd(vlpv,nmeas,nctrl,gam,xmat1,ymat1,');
 disp(' xmat2,ymat2,vnu,fparm,gparm,gradf,gradg,bparm,gradb)');
 return
end

% If vlpv, xmat1, xmat2, ymat1, or ymat2 aren't varying matrices, vpck them
[typ,dum2,dum3,npts1] = minfo(vlpv);
if typ == 'syst'
  vlpv = vpck(vlpv,1);
  [typ,dum2,dum3,npts1] = minfo(vlpv);
end
viv = getiv(vlpv);

[typ,xrow1,xcol1,numx1] = minfo(xmat1);
if typ == 'cons'
  xmat1 = vpck(xmat1,1);
  [typ,xrow1,xcol1,numx1] = minfo(xmat1);
end
[typ,yrow1,ycol1,numy1] = minfo(ymat1);
if typ == 'cons'
  ymat1 = vpck(ymat1,1);
  [typ,yrow1,ycol1,numy1] = minfo(ymat1);
end
[typ,xrow2,xcol2,numx2] = minfo(xmat2);
if typ == 'cons'
  xmat2 = vpck(xmat2,1);
  [typ,xrow2,xcol2,numx2] = minfo(xmat2);
end
[typ,yrow2,ycol2,numy2] = minfo(ymat2);
if typ == 'cons'
  ymat2 = vpck(ymat2,1);
  [typ,yrow2,ycol2,numy2] = minfo(ymat2);
end
if (xrow1 ~= xcol1) | (yrow1 ~= ycol1) | (xrow2 ~= xcol2) | (yrow2 ~= ycol2)
  error('X and Y must be square');
end

[typ,nparm,nbnds,npts8] = minfo(vnu);
if typ ~= 'vary'
 [nparm,nbnds] = size(vnu);
 npts8 = npts1;
 vnu = vpck(kron(ones(npts8,1),vnu),viv);
end
if nbnds > 2
 GRID_PARMV = 1;
else
 GRID_PARMV = 0;
end

if isempty(fparm)
 fparm = vpck(ones(npts1,1),1:npts1);
 gradf = vpck(zeros(npts1,nparm),1:npts1);
end
if isempty(gparm)
 gparm = vpck(ones(npts1,1),1:npts1);
 gradg = vpck(zeros(npts1,nparm),1:npts1);
end

% Check that the number of grid points and parameters is consistent
[dum1,nbasisx,ckvx,npts2] = minfo(fparm);
[dum1,nbasisy,ckvy,npts3] = minfo(gparm);
[dum1,nblends,ckvb,npts4] = minfo(bparm);
if (ckvx ~= 1) | (ckvy ~= 1)
 error('Error: Basis data should be column vectors.');
end
[dum1,nbasisgx,nparmgx,npts5] = minfo(gradf);
[dum1,nbasisgy,nparmgy,npts6] = minfo(gradg);
[dum1,ngblends,nparmgb,npts7] = minfo(gradb);
if any([ngblends nblends ckvb] ~= 1)
 error('Error: Only one scalar blending function allowed.');
end
if any([nparmgx nparmgy nparmgb] ~= nparm)
 error('Error: Number of parameters in gradient & nu data inconsistent.');
end
if any([numx1 numx2 nbasisgx] ~= nbasisx) | ...
	any([numy1 numy2 nbasisgy] ~= nbasisy) 
 error('Error: Number of basis functions not consistent');
end
if any([npts2 npts3 npts4 npts5 npts6 npts7 npts8] ~= npts1)
 disp(['There are ' int2str(npts1) ' points in vlpv']);
 disp(['There are ' int2str(npts2) ' points in fparm']);
 disp(['There are ' int2str(npts3) ' points in gparm']);
 disp(['There are ' int2str(npts4) ' points in bparm']);
 disp(['There are ' int2str(npts5) ' points in gradf']);
 disp(['There are ' int2str(npts6) ' points in gradg']);
 disp(['There are ' int2str(npts7) ' points in gradb']);
 disp(['There are ' int2str(npts8) ' points in vnu']);
 error('Error: number of grid points inconsistent');
end  
nparm = nparmgx;
npts = npts1;
Npts = int2str(npts);

% Identify the parameters on which X and Y vary 
parmx = any(vunpck(gradf));
parmy = any(vunpck(gradg));
nparmx = sum(parmx);
nparmy = sum(parmy);
if GRID_PARMV & nparmx > 0
 nvertx = nbnds;
else
 nvertx = 2^nparmx;
end
if GRID_PARMV & nparmy > 0
 nverty = nbnds;
else
 nverty = 2^nparmy;
end

% Get matrix containing all combinations of +/- for nparm-dim vector.
combmatx = corners(nparmx);    % combmatx = 1 if nparmx = 0
combmaty = corners(nparmy);    % combmaty = 1 if nparmy = 0

sys = xtracti(vlpv,1,1);
[systype,no,ni,nx] = minfo(sys);
ne = no-nmeas;
ne1 = ne-nctrl;
nd = ni-nctrl;
nd1 = nd-nmeas;

eigx = zeros(npts,nvertx);
eigy = zeros(npts,nverty);
eigxy = zeros(npts,1);
fail = [];

dxmat = msub(xmat1,xmat2);
dymat = msub(ymat1,ymat2);

for i = 1:npts

% Form the Lyapunov matrices at parm_i
 fdat = xtracti(fparm,i,1);
 gdat = xtracti(gparm,i,1);
 bdat = xtracti(bparm,i,1);
 gfdat = xtracti(gradf,i,1);
 ggdat = xtracti(gradg,i,1);
 gbdat = xtracti(gradb,i,1);
 nu = xtracti(vnu,i,1);

% Check for a two-sided rate bound, or just a bound on absolute value
 if nbnds == 1
  nu = [-nu nu];
 end

 X1 = zeros(nx,nx);
 Y1 = zeros(nx,nx);
 X2 = zeros(nx,nx);
 Y2 = zeros(nx,nx);
 for k = 1:nbasisx
  X1 = X1 + fdat(k) * xtracti(xmat1,k,1);
 end
 for k = 1:nbasisy
  Y1 = Y1 + gdat(k) * xtracti(ymat1,k,1);
 end
 for k = 1:nbasisx
  X2 = X2 + fdat(k) * xtracti(xmat2,k,1);
 end
 for k = 1:nbasisy
  Y2 = Y2 + gdat(k) * xtracti(ymat2,k,1);
 end

 X = bdat*X1 + (1-bdat)*X2;
 Y = bdat*Y1 + (1-bdat)*Y2;
 DX = X1 - X2;
 DY = Y1 - Y2;
 xmat = madd(mscl(xmat1,bdat),mscl(xmat2,1-bdat));
 ymat = madd(mscl(ymat1,bdat),mscl(ymat2,1-bdat));

% Get state-space data for parm_i.
 sys = xtracti(vlpv,i,1);
 [a,b1,b2,c1,c2,d11] = transf(sys,nmeas,nctrl);
 [trow,tcol]=size(d11);
 if min(eig(eye(tcol)-d11'*d11)) <= 0
   disp('I - D11*D11 < 0');
 end
 b11 = b1(:,1:nd1);
 b12 = b1(:,nd1+1:nd);
 c11 = c1(1:ne1,:);
 c12 = c1(ne1+1:ne,:);
 d1111 = d11(1:ne1,1:nd1);
 d1112 = d11(1:ne1,nd1+1:nd);
 d1121 = d11(ne1+1:ne,1:nd1);
 d1122 = d11(ne1+1:ne,nd1+1:nd);
 ahat = a-b2*c12;
 atld = a-b12*c2;
 b1hat = b1-b2*[d1121 d1122];
 c1tld = c1-[d1112;d1122]*c2;

 % Check state feedback LMIs
 for j = 1:nverty
  Ydot = zeros(nx,nx);
  bdot = 0;
  ly = 0;
  for l = 1:nparm
   if parmy(l) ~= 0
    ly = ly + 1;
    if GRID_PARMV
     parmvl = nu(l,j);
    else
     parmvl = (nu(l,1)+nu(l,2))/2 + combmaty(j,ly)*(nu(l,1)-nu(l,2))/2;
    end
    for k = 1:nbasisy
     Ydot = Ydot + parmvl * ggdat(k,l) * xtracti(ymat,k,1);
    end
    bdot = bdot + parmvl * gbdat(1,l);
   end
  end  
  Ydot = Ydot + bdot * DY;
  LMI_Y = [Y*ahat' + ahat*Y  - Ydot - gam*b2*b2', Y*c11', b1hat;
	   c11*Y, -gam*eye(ne1), [d1111 d1112];
	   b1hat', [d1111 d1112]', -gam*eye(nd)];
  eigy(i,j) = max(eig(LMI_Y));
  if eigy(i,j) >= 0
   disp(['Grid point ' int2str(i) ' of ' Npts ' failed']);
   fail = [fail;[i j 2]];
  end
 end

 % Check output injection LMIs
 for j = 1:nvertx
  Xdot = zeros(nx,nx);
  bdot = 0;
  lx = 0;
  for l = 1:nparm
   if parmx(l) ~= 0
    lx = lx + 1;
    if GRID_PARMV
     parmvl = nu(l,j);
    else
     parmvl = (nu(l,1)+nu(l,2))/2 + combmatx(j,lx)*(nu(l,1)-nu(l,2))/2;
    end
    for k = 1:nbasisx
     Xdot = Xdot + parmvl * gfdat(k,l) * xtracti(xmat,k,1);
    end
    bdot = bdot + parmvl * gbdat(1,l);
   end
  end
  Xdot = Xdot + bdot * DX;
  LMI_X = [atld'*X + X*atld  + Xdot - gam*c2'*c2, X*b11, c1tld';
	   b11'*X, -gam*eye(nd1), [d1111;d1121]';
	   c1tld, [d1111;d1121], -gam*eye(ne)];
  eigx(i,j) = max(eig(LMI_X));
  if eigx(i,j) >= 0
   disp(['Grid point ' int2str(i) ' of ' Npts ' failed']);
   fail = [fail;[i j 1]];
  end
 end

 % Check spectral radius LMI
 LMI_XY = [Y, eye(nx);
           eye(nx), X];
 eigxy(i) = min(eig(LMI_XY));
 if eigxy(i) <= 0
  disp(['Grid point ' int2str(i) ' of ' Npts ' failed']);
  fail = [fail;[i 1 0]];
 end

end