lpvchkr.m

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

function [eigx,eigy,eigxy,fail] = ...
  lpvchkr(vlpv,nmeas,nctrl,gam,xmat,ymat,vnu,fparm,gparm,gradf,gradg)

% [eigx,eigy,eigxy,fail] = 
%  lpvchkr(vlpv,nmeas,nctrl,gam,xmat,ymat,vnu,fparm,gparm,gradf,gradg)
% checks reduced-order PDLF LMI solutions on a grid.
%
% [eigx,eigy,eigxy,fail] = lpvchkr(vlpv,nmeas,nctrl,gam,xmat,ymat)
% checks reduced-order SQLF LMI solutions on a grid.
%
% Inputs: 
% vlpv 		varying matrix of the lpv plant at parm values
% nmeas,nctrl	number of measurements & controls, respectively
% gam 		optimal gamma from solver
% xmat,ymat	varying-matrix LMI solutions 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.
%
% 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		3-column matrix of indices to 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] = ');
 disp('	lpvchkr(vlpv,nmeas,nctrl,gam,xmat,ymat,vnu,fparm,gparm,gradf,gradg)');
 return
end

% If vlpv, xmat, or ymat 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
[typ,xrow,xcol,numx] = minfo(xmat);
if typ == 'cons'
  xmat = vpck(xmat,1);
  [typ,xrow,xcol,numx] = minfo(xmat);
end
[typ,yrow,ycol,numy] = minfo(ymat);
if typ == 'cons'
  ymat = vpck(ymat,1);
  [typ,yrow,ycol,numy] = minfo(ymat);
end
if (xcol ~= yrow) | (xcol ~= ycol) | (xrow > xcol)
  error('Inconsistent dimensions in X and/or Y');
end
viv = getiv(vlpv);

% If constant X and Y
if nargin == 6
 vnu = vpck(zeros(npts1,1),viv);
 nparmx = 0;
 nparmy = 0;
 nvertx = 1;
 nverty = 1;
 fparm = [];
 gparm = [];
 PDLF = 0;
else
 PDLF = 1;
end

[typ,nparm,nbnds,npts6] = minfo(vnu);
if typ ~= 'vary'
 [nparm,nbnds] = size(vnu);
 npts6 = npts1;
 vnu = vpck(kron(ones(npts6,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);
if (ckvx ~= 1) | (ckvy ~= 1)
 error('Error: Basis data should be column vectors.');
end
[dum1,nbasisgx,nparmgx,npts4] = minfo(gradf);
[dum1,nbasisgy,nparmgy,npts5] = minfo(gradg);
if (nparmgx ~= nparmgy)
 error('Error: Number of parameters in basis gradient data inconsistent.');
end
if (nparmgx ~= nparm)
 error('Error: Number of parameters in nu inconsistent with gradient data.');
end
if (nbasisx ~= numx) | (nbasisy ~= numy)
 error('Error: Basis dimension not consistent');
end
if (npts1 ~= npts2) | (npts2 ~= npts3) | (npts3 ~= npts4) | ... 
 	(npts4 ~= npts5) | (npts5 ~= 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 gradf']);
 disp(['There are ' int2str(npts5) ' points in gradg']);
 disp(['There are ' int2str(npts6) ' 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 should vary 
parmx = any(vunpck(gradf)) & any(vunpck(vtp(vnu)));
parmy = any(vunpck(gradg)) & any(vunpck(vtp(vnu)));
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);
nx1 = xrow;
nx2 = nx-xrow;
ny1 = nmeas-nx2;
ny2 = nx2;
ne = no-nmeas;
ne1 = ne-nctrl;
nd = ni-nctrl;
nd1 = nd-ny1;

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

for i = 1:npts

% Form the Lyapunov matrices at parm_i
 fdat = xtracti(fparm,i,1);
 gdat = xtracti(gparm,i,1);
 gfdat = xtracti(gradf,i,1);
 ggdat = xtracti(gradg,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

 X = zeros(nx1,nx);
 Y = zeros(nx,nx);
 for k = 1:nbasisx
  X = X + fdat(k) * xtracti(xmat,k,1);
 end
 for k = 1:nbasisy
  Y = Y + gdat(k) * xtracti(ymat,k,1);
 end
 X11 = X(:,1:nx1);
 X12 = X(:,nx1+1:nx);

% Get state-space data for parm_i.
 sys = xtracti(vlpv,i,1);
 [A,Bd,Bu,Ce,Cy,Ded] = transfr(sys,nmeas,nctrl,nx2);
 [trow,tcol]=size(Ded);
 if min(eig(eye(tcol)-Ded'*Ded)) <= 0
   disp('I - Ded*Ded < 0');
 end
 B11 = Bd(1:nx1,1:nd1);
 B12 = Bd(1:nx1,nd1+1:nd);
 B21 = Bd(nx1+1:nx,1:nd1);
 B22 = Bd(nx1+1:nx,nd1+1:nd);
 Ce1 = Ce(1:ne1,:);
 Ce2 = Ce(ne1+1:ne,:);
 Cy1 = Cy(1:ny1,1:nx1);
 Ded11 = Ded(1:ne1,1:nd1);
 Ded12 = Ded(1:ne1,nd1+1:nd);
 Ded21 = Ded(ne1+1:ne,1:nd1);
 Ded22 = Ded(ne1+1:ne,nd1+1:nd);
 Ahat = A-Bu*Ce2;
 Bdhat = Bd-Bu*[Ded21 Ded22];
 if ny1 == 0 	% (state feedback only)
  Atld11 = A(1:nx1,1:nx1);
  Atld21 = A(nx1+1:nx,1:nx1);
  Cetld = Ce(:,1:nx1);
 else
  Atld11 = A(1:nx1,1:nx1) - B12*Cy1;
  Atld21 = A(nx1+1:nx,1:nx1) - B22*Cy1;
  Cetld = Ce(:,1:nx1) - [Ded12;Ded22]*Cy1;
 end
 
 % Check state feedback LMIs
 for j = 1:nverty
  Ydot = zeros(nx,nx);
  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
   end
  end  
  LMI_Y = [Y*Ahat' + Ahat*Y  - Ydot - gam*Bu*Bu', Y*Ce1', Bdhat;
	   Ce1*Y, -gam*eye(ne1), [Ded11 Ded12];
	   Bdhat', [Ded11 Ded12]', -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
  X11dot = zeros(nx1,nx1);
  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
     X11dot = X11dot + parmvl * gfdat(k,l) * sel(xtracti(xmat,k,1),':',1:nx1);
    end
   end
  end
  if ny1 == 0
   LMI_X11 = [Atld11;Atld21]'*X' + X*[Atld11;Atld21] + X11dot;
  else
   LMI_X11 = [Atld11;Atld21]'*X' + X*[Atld11;Atld21] + X11dot - gam*Cy1'*Cy1;
  end
  LMI_X = [LMI_X11, X*[B11;B21], Cetld';
	   [B11;B21]'*X', -gam*eye(nd1), [Ded11;Ded21]';
	   Cetld, [Ded11;Ded21], -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
 if PDLF | i == 1
  LMI_XY = [Y, [eye(nx1);zeros(nx2,nx1)];
            [eye(nx1) zeros(nx1,nx2)], X11];
  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

end