lpvchkpq.m

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

function [eigx,eigy,eigxy,fail] = ...
 lpvchkpq(vlpv,dim,vnu,fgparm,gradfg,psi,gam,xmat,ymat)

% [eigx,eigy,eigxy,fail] = 
%  lpvchkpq(vlpv,dim,vnu,fgparm,gradfg,psi,gam,xmat,ymat)
% checks LMI solutions to the parameter-dependent OSF LPV control problem
% with dynamic parameter measurement.
%
% INPUTS: 
% vlpv: (VARYING) matrix containing L evaluations of the OLIC:
%       vlpv = vpck([olic1;...;olicL],1:L)
%       Any direct state measurements must be ordered last (in A and in Cy).
% dim = [nmeas nctrl NX NY nsf]: # of measurements, controls, X basis functions,
%       Y basis functions, and (optionally) directly measured states.
% vnu:  Np*Nq-point VARYING (2s x nbnds) matrix grid of parameter velocities
%       and pdot-q.
%       nbnds = 1: -vnu(i)   < parmv(i) < vnu(i)
%                 -vnu(s+i) < parmv(i)-q < vnu(s+i)
%       nbnds = 2: -vnu(i,1) < parmv(i) < vnu(i,2)
%               -vnu(s+i,1) < parmv(i)-q < vnu(s+i,2)
%       nbnds > 2: vnu(1:s,j) is the j-th corner of the velocity polytope.
%               vnu(s+(1:s),j) is the j-th corner of the pdot-q polytope.
% fgparm=abv(fparm,gparm): fparm,gparm are basis function data for X(p,q),Y(p,q)
% gradfg=abv(gradf,gradg): gradf,gradg are basis gradient data for X(p,q),Y(p,q)
%       VARYING matrices with Np*Nq points (p = outer loop, q = inner loop)
%       sel(fparm,i,1) = f_i(p,q),
%       sel(gradf,i,j) = df_i(p,q)/dp_j for j=1,...,s
%       sel(gradf,i,j) = df_i(p,q)/dq_j for j=s+1,...,2s
% psi:  (s x 1) vector of bandwidths for parameter measurement filters
% gam 		optimal gamma from solver
% xmat,ymat	varying-matrix LMI solutions from solver
%
% 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('lpvchkpq(vlpv,dim,vnu,fgparm,gradfg,psi,gam,xmat,ymat)');
 return
end

% Get some problem dimensions
nmeas = dim(1); nctrl = dim(2); NX = dim(3); NY = dim(4);
if length(dim) > 4, nsf = dim(5); else nsf = 0; end
gam = vunpck(gam);
ngam = length(gam);

% If vlpv, xmat, or ymat aren't varying matrices, vpck them
[typ,dum2,dum3,npts] = minfo(vlpv);
if typ == 'syst'
  vlpv = vpck(vlpv,1);
  [typ,dum2,dum3,npts] = 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);

% Analyze vnu,fgparm,gradfg dimensions
[dum1,nbasis,ncols,npts2] = minfo(fgparm);
if ncols ~= 1
 error('Basis functions must be vectors.');
end
[dum1,ngbasis,nparmg,npts3] = minfo(gradfg);
Nq = npts2/npts;
Npts = int2str(Nq*npts)
nparmg = nparmg/2;
[type,nparm,nbnds,npts1] = minfo(vnu);
if type ~= 'vary'
 vnu = vpck(kron(ones(npts*Nq,1),vnu),kron(getiv(fgparm)));
 [type,nparm,nbnds,npts1] = minfo(vnu);
end
if nbnds > 2
 GRID_PARMV = 1;
else
 GRID_PARMV = 0;
 if nbnds == 1, vnu = sbs(mscl(vnu,-1),vnu); end
end
nparm = nparm/2;

% Check for dimension inconsistencies
if all(ngam ~= [1 Nq npts*Nq])
 error('Number of gammas inconsistent with grid data')
elseif any([nbasis ngbasis] ~= NX+NY) | (numx ~= NX) | (numy ~= NY)
 error('Number of basis functions inconsistent.');
elseif any(nparmg ~= [nparm length(psi)])
 error('Number of parameters inconsistent.');
elseif any(npts1 ~= [npts2 npts3]) | (round(Nq) ~= Nq)
 disp(['There are ' int2str(npts) ' points in vlpv'])
 disp(['There are ' int2str(npts1) ' points in vnu'])
 disp(['There are ' int2str(npts2) ' points in fgparm'])
 disp(['There are ' int2str(npts3) ' points in gradfg'])
 error('Number of grid points inconsistent.');
else
 fparm = sel(fgparm,1:NX,':');
 gparm = sel(fgparm,NX+(1:NY),':');
 gradf = sel(gradfg,1:NX,':');
 gradg = sel(gradfg,NX+(1:NY),':');
end

% Identify parameters having finite, nonzero rate bounds
parmx = any(vunpck(gradf)) & any(vunpck(vtp(vnu)));
parmx = parmx(1:nparm) | parmx(nparm+(1:nparm));
parmy = any(vunpck(gradg)) & any(vunpck(vtp(vnu)));
parmy = parmy(1:nparm) | parmy(nparm+(1:nparm));
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

% Get more problem dimensions
[type,no,ni,nx] = minfo(xtracti(vlpv,1,1));
nx1 = nx-nsf; nx2 = nsf;
ny1 = nmeas-nsf; ny2 = nsf;
nd = ni-nctrl; nd1 = nd-ny1;
ne = no-nmeas; ne1 = ne-nctrl;

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

for i = 1:npts

% 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
 
% Cycle through q grid
for ii = 1:Nq
 iii = Nq*(i-1)+ii;
 if (ngam == Nq)
  igam = ii;
 elseif (ngam == npts*Nq)
  igam = iii;
 else
  igam = 1;
 end
 gamma = gam(igam);

% Form the Lyapunov matrices at parm_i
 nu = xtracti(vnu,iii,1);
 fdat = xtracti(fparm,iii,1);
 gdat = xtracti(gparm,iii,1);
 gfdat = xtracti(gradf,iii,1);
 ggdat = xtracti(gradg,iii,1);
 X = zeros(nx1,nx);
 Y = zeros(nx,nx);
 for k = 1:NX
  X = X + fdat(k) * xtracti(xmat,k,1);
 end
 for k = 1:NY
  Y = Y + gdat(k) * xtracti(ymat,k,1);
 end
 X11 = X(:,1:nx1);
 X12 = X(:,nx1+1:nx);

 % Check state feedback LMIs
 for j = 1:nverty
  ly = 0;
  parmv = zeros(2*nparm,1);
  for l = 1:nparm
   ll = nparm + l;
   if parmy(l) ~= 0
    ly = ly + 1;
    if GRID_PARMV
     parmv(l) = nu(l,j);
     parmv(ll) = nu(ll,j) * psi(l);
    else
     parmv(l) = (nu(l,1)+nu(l,2))/2 + combmaty(j,ly)*(nu(l,1)-nu(l,2))/2;
     parmv(ll) = (nu(ll,1)+nu(ll,2))/2 + combmaty(j,ly)*(nu(ll,1)-nu(ll,2))/2;
    end
   end
  end  
  Ydot = zeros(nx,nx);
  for k = 1:NY
   Ydot = Ydot + (ggdat(k,:) * parmv) * xtracti(ymat,k,1);
  end
  LMI_Y = [Y*Ahat' + Ahat*Y  - Ydot - gamma*Bu*Bu', Y*Ce1', Bdhat;
	   Ce1*Y, -gamma*eye(ne1), [Ded11 Ded12];
	   Bdhat', [Ded11 Ded12]', -gamma*eye(nd)];
  eigy(iii,j) = max(eig(LMI_Y));
  if eigy(iii,j) >= 0
   disp(['Grid point ' int2str(iii) ' of ' Npts ' failed']);
   fail = [fail;[iii j 2]];
  end
 end

 % Check output injection LMIs
 for j = 1:nvertx
  lx = 0;
  parmv = zeros(2*nparm,1);
  for l = 1:nparm
   ll = nparm + l;
   if parmx(l) ~= 0
    lx = lx + 1;
    if GRID_PARMV
     parmv(l) = nu(l,j);
     parmv(ll) = nu(ll,j) * psi(l);
    else
     parmv(l) = (nu(l,1)+nu(l,2))/2 + combmatx(j,lx)*(nu(l,1)-nu(l,2))/2;
     parmv(ll) = (nu(ll,1)+nu(ll,2))/2 + combmaty(j,ly)*(nu(ll,1)-nu(ll,2))/2;
    end
   end
  end
  X11dot = zeros(nx1,nx1);
  for k = 1:NX
   X11dot = X11dot + (gfdat(k,:) * parmv) * sel(xtracti(xmat,k,1),':',1:nx1);
  end
  if ny1 == 0
   LMI_X11 = [Atld11;Atld21]'*X' + X*[Atld11;Atld21] + X11dot;
  else
   LMI_X11 = [Atld11;Atld21]'*X' + X*[Atld11;Atld21] + X11dot - gamma*Cy1'*Cy1;
  end
  LMI_X = [LMI_X11, X*[B11;B21], Cetld';
	   [B11;B21]'*X', -gamma*eye(nd1), [Ded11;Ded21]';
	   Cetld, [Ded11;Ded21], -gamma*eye(ne)];
  eigx(iii,j) = max(eig(LMI_X));
  if eigx(iii,j) >= 0
   disp(['Grid point ' int2str(iii) ' of ' Npts ' failed']);
   fail = [fail;[iii j 1]];
  end
 end

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

end % q loop
end % p loop