lpvperf.m

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

% [gam,xmat,xopt] = lpvperf(vlpv,vnu,fparm,gradf,xinit)
% Calculates solution to the parameter-dependent induced-L2 analysis problem.
%
% [gam,xmat,ymat,xyopt] = lpvperf(vlpv)
% Calculates solution to the (SQLF) quadratic performance analysis problem.
%
% INPUTS:  
% vlpv: VARYING matrix containing N evaluations of the open-loop IC's:
%         vlpv = vpck([olic(parm1);...;olic(parmN)],1:N)
% vnu:  (s x 1) [(s x 2)] matrix of [lower and upper] rate bounds on all s
%       parameters, with garbage entries for parameters without rate bounds:
%         -nu(i) <= parm(i) <= nu(i)  OR  nu(i,1) <= parm(i) <= nu(i,2) 
%       NOTE: vnu can be parameter-varying, i.e. a VARYING matrix like vlpv.
%       NOTE: vnu with more than two columns is considered a
%       grid of parameter velocities, e.g. vertices of the velocity polytope. 
% fparm,gradf: contains basis function (and gradient) data for X(parm).
% xinit: initial guess for the LMI variable. (OPTIONAL)
%
% OUTPUTS: 
% gam: The optimal performance level.
% xmat: X_k's packed in a varying matrix.
% xopt: The LMI variables (W,gam) strung out as a vector.

function [gam,xmat,xopt] = lpvperf(vlpv,vnu,fparm,gradf,xinit);

tstart = cputime;
if all(nargin ~= [1 4 5])
 disp('usage: [gam,xmat,xopt] = lpvperf(vlpv,vnu,fparm,gradf,xinit)')
 return
end
disp('  Setting up LMIs...');

% If vlpv is a system matrix, pack it as a varying matrix.
[typ,dum2,dum3,npts1] = minfo(vlpv);
if typ == 'syst'
 vlpv = vpck(vlpv,1);
 [typ,dum2,dum3,npts1] = minfo(vlpv);
end
viv = getiv(vlpv);

if (nargin == 3)     % Fixed QLF
 viv = getiv(vlpv);
 vnu = vpck(zeros(npts1,1),viv);
 fparm = [];
 gradf = [];
 nparmx = 0;
 nvertx = 1;
end
if ~exist('xinit')
 xinit = [];
end

% If vnu is a constant vector, pack it as a varying matrix.
[typ,nparm,nbnds,npts6] = minfo(vnu);
if typ ~= 'vary'
 [nparm,nbnds] = size(vnu);
 npts6 = npts1;
 vnu = vpck(kron(ones(npts4,1),vnu),viv);
end
if nbnds > 2
 GRID_PARMV = 1;
else
 GRID_PARMV = 0;
end

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

% Check that the number of grid points is consistent.
[dum1,nbasisx,ckvx,npts2] = minfo(fparm);
if (ckvx ~= 1) 
 error('Basis data should be column vectors.');
end
[dum1,nbasisgx,nparmgx,npts3] = minfo(gradf);
if (nparmgx ~= nparm) 
 error('Number of parameters in nu inconsistent with gradient data.');
end
if (nbasisgx ~= nbasisx) 
 error('Number of basis not consistent.');
end
if (npts1 ~= npts2) | (npts2 ~= npts3) | (npts3 ~= npts4) 
 disp(['There are ' int2str(npts1) ' points in vlpv']);
 disp(['There are ' int2str(npts2) ' points in fparm']);
 disp(['There are ' int2str(npts3) ' points in gradf']);
 disp(['There are ' int2str(npts4) ' points in vnu']);
 error('Number of grid points inconsistent');
end  
grid_pts = npts1;

% Identify the parameters whose rates are bounded (W depends on them)
% and are nonzero (vnu_i > 0)
parmx = any(vunpck(gradf)) & any(vunpck(vtp(vnu)));
nparmx = sum(parmx);
if GRID_PARMV & nparmx > 0
 nvertx = nbnds;
else
 nvertx = 2^nparmx;
end

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

sys = xtracti(vlpv,1,1);
[systype,no,ni,nx] = minfo(sys);

nvardec = nx*(nx+1)/2*nbasisx+1;
nvarx = nbasisx;

% Setup lmi matrix variables
setlmis([])
for i = 1:nbasisx
 X = lmivar(1,[nx 1]); 
end
gamma = lmivar(1,[1 0]);

% Add data to variable "lmis".
for i = 1:grid_pts
 disp([' Adding data for grid point ' int2str(i) ' of ' int2str(grid_pts)]);

% Get values of basis functions, gradient, and rate bound for parm_i.
 fdat = xtracti(fparm,i,1);
 gfdat = xtracti(gradf,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

% Get state-space data for grid point i
 sys = xtracti(vlpv,i,1);
 [A,B,C,D] = unpck(sys);
 
 if nargin > 1
  indx = (nvertx+1)*(i-1);
  indx1 = (nvertx+1)*i;
 else 
  indx = i-1;
  indx1 = grid_pts+1;
 end

 ezmatx = [eye(nx) zeros(nx,ni)];
 for j = 1:nvertx
  for k = 1:nbasisx
   if abs(fdat(k)) > eps 
%    lmiterm([indx+j 1 1 k],fdat(k)*[A';B'],ezmatx,'s');
    lmiterm([indx+j 1 1 k],fdat(k),A,'s');
    lmiterm([indx+j 2 1 k],fdat(k)*B',1);
   end
   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
     if abs(gfdat(k,l)) > eps & abs(parmvl) > eps
%      lmiterm([indx+j 1 1 k],parmvl*gfdat(k,l)*ezmatx',ezmatx);
      lmiterm([indx+j 1 1 k],parmvl*gfdat(k,l),1);
     end
    end  
   end
  end
%  lmiterm([indx+j 1 1 nvarx+1],-1,daug(zeros(nx,nx),eye(ni)));
%  lmiterm([indx+j 2 1 0],[C, D]);
%  lmiterm([indx+j 2 2 nvarx+1],-1,1);
  lmiterm([indx+j 2 2 nvarx+1],-1,1);
  lmiterm([indx+j 3 1 0],C);
  lmiterm([indx+j 3 2 0],D);
  lmiterm([indx+j 3 3 nvarx+1],-1,1);
 end

% This number can affect the largest pole of the controller: 
delt = 1e-6;

% If PD X 
 if nargin > 1
  for k=1:nbasisx
   if abs(fdat(k)) > eps
    lmiterm([-indx1 1 1 k],fdat(k),1);
   end
  end
 end

end % grid point loop

% If fixed X
if nargin == 3
 lmiterm([-indx2 1 1 1],1,1);
 lmiterm([indx2 1 1 0],delt);
end

lmis=getlmis;  

% Use initial condition if given.
if (xinit ~= [])
 decinit = xinit;
else
 decinit = [];
end
cobj = [zeros(nx*(nx+1)/2*nvarx,1); 1];

% Save data to file and display setup time
datesim = date;
timesim = clock;
ttwo = cputime;
dtime = ttwo-tstart;
if dtime < 60
 disp(['Done.  CPU Time = ' num2str(dtime) ' sec']);
elseif dtime < 3600
 disp(['Done.  CPU Time = ' num2str(dtime/60) ' min']);
else
 disp(['Done.  CPU Time = ' num2str(dtime/3600) ' hrs!!']);
end
 
% Call MINCX
tic
[copt,xyopt] = mincx(lmis,cobj,[1e-2 200 1e7 10 0],decinit,0);
toc

% Calculate solution time and display results.
tthre = cputime;
dtime1 = tthre-ttwo;
dtime2 = tthre-tstart;
timestr = 'Done with optimization. CPU Time (LINOBJ) = ';
if dtime1 < 60
 disp([timestr num2str(dtime1) ' sec']);
elseif dtime2 < 3600
 disp([timestr num2str(dtime1/60) ' min']);
else
 disp([timestr num2str(dtime1/3600) ' hrs!!']);
end
if dtime2 < 60
 disp(['                CPU Time (total) = ' num2str(dtime2) ' sec']);
elseif dtime2 < 3600
 disp(['                CPU Time (total) = ' num2str(dtime2/60) ' min']);
else
 disp(['                CPU Time (total) = ' num2str(dtime2/3600) ' hrs!!']);
end
 
if isempty(copt)
 disp(['Sorry. No feasible solution found.']);
 return;
end
 
% Form output arguments
gam = xyopt(nvardec);
xmat = [];
for i = 1:nbasisx
 xmat = [xmat;dec2mat(lmis,xopt,i)];
end
if nargin > 1
 xmat=vpck(xmat,[1:nbasisx]');
end