lpvsyn.m

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

function [vctrl,vclp] = lpvsyn(vlpv,ny,nu,gam,xmat,ymat,fparm,gparm)

% function [vctrl,vclp] = lpvsyn(vlpv,ny,nu,gam,xmat,ymat,fparm,gparm)
% or
% function [vctrl,vclp] = lpvsyn(vlpv,ny,nu,gam,xmat,ymat)
%
% Evaluate full/reduced-order LPV controller and closed-loop system 
% at grid points of the open-loop LPV system.
% If X and Y are both varying, then these controllers apply to CONSTANT
% parameter trajectories only.
%
% Inputs: vlpv		varying matrix of the lpv plant at parm values
%	  ny,nu		number of measurements & controls, respectively
%	  gam		optimal gamma from lpvsolr (can be varying)
%	  xmat,ymat 	LMI solutions from lpvsolr in varying matrices.
%			The no. of measured states is implied by xmat.
%	  fparm,gparm 	basis-function values at parm values
% 
%			If X or Y is constant, enter [] for fparm or gparm.
%			If both are constant, exclude fparm and gparm.
%
% Output: vctrl 	varying matrix of packed controllers
% 	  vclp 		varying matrix of packed closed-loop systems

if nargin == 0
 disp('usage: [vctrl,vclp] = lpvsyn(vlpv,ny,nu,gam,xmat,ymat,fparm,gparm)');
 disp('       or');
 disp('       [vctrl,vclp] = lpvsyn(vlpv,ny,nu,gam,xmat,ymat)'); 
 return;
end

% If vlpv, gam, 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
viv = getiv(vlpv);
[typ,dum2,dum3,npts0] = minfo(gam);
if typ ~= 'vary'
 gam = vpck(gam*ones(npts1,1),viv);
 npts0 = npts1;
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 constant X and/or Y
if ~exist('fparm') 
 fparm = vpck(ones(npts1,1),1:npts1);
end
if ~exist('gparm')
 gparm = vpck(ones(npts1,1),1:npts1);
end
if isempty(fparm)
 fparm = vpck(ones(npts1,1),1:npts1);
end
if isempty(gparm)
 gparm = vpck(ones(npts1,1),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
 if (nbasisx ~= numx) | (nbasisy ~= numy)
   error('Error: Basis dimension not consistent');
 end
 if (npts1 ~= npts2) | (npts2 ~= npts3) | (npts3 ~= npts0) 
  disp(['There are ' int2str(npts0) ' points in gam']);
  disp(['There are ' int2str(npts1) ' points in vlpv']);
  disp(['There are ' int2str(npts2) ' points in fparm']);
  disp(['There are ' int2str(npts3) ' points in gparm']);
  error('Error: number of grid points inconsistent');
 end  
grid_pts = npts1;

sys = xtracti(vlpv,1,1);
[systype,no,ni,nx] = minfo(sys);
if any(nx ~= [xcol ycol yrow]) | (xrow > xcol)
 error('Size of state vector, [X11 X12], and/or Y inconsistent.');
end
nx1 = xrow;
nx2 = nx-nx1;
ny1 = ny-nx2;
ny2 = nx2;
ne = no-ny;
ne1 = ne-nu;
nd = ni-nu;
nd1 = nd-ny;

if min(numx,numy) > 1
 disp('Warning: Both X and Y are varying. Output applies to LTI case only.')
end

sysrows = nx1+nu+1;
syscols = nx1+ny+1;
clprows = (nx+nx1)+ne+1;
clpcols = (nx+nx1)+nd+1;
vctrl = zeros(sysrows*grid_pts,syscols);
vclp = zeros(clprows*grid_pts,clpcols);

[a,b,c,d] = veval('unpck',vlpv);
d11 = sel(d,1:ne,1:nd);
DED_FLAG = pkvnorm(veval('*',d11,minv(gam)));
if DED_FLAG > 1 - eps
 disp('Strictly proper controller infeasible, using controller with Dk1 term');
 DK1_FLAG = 1;
else
 disp('Using strictly proper controller formulae');
 DK1_FLAG = 0;
end
DED_FLAG = veval('sign',DED_FLAG);

for i=1:grid_pts

 % Form the lyapunov matrices for parm_i
 X = zeros(nx1,nx);
 Y = zeros(nx);
 gamma = xtracti(gam,i,1);
 fvec = xtracti(fparm,i,1);
 gvec = xtracti(gparm,i,1);
 for k=1:nbasisx
  X = X + fvec(k) * xtracti(xmat,k,1);
 end
 for k=1:nbasisy
  Y = Y + gvec(k) * xtracti(ymat,k,1);
 end
 Y11 = Y(1:nx1,1:nx1);
 Y12 = Y(1:nx1,nx1+1:nx);
 Y22 = Y(nx1+1:nx,nx1+1:nx);
 Yi = inv(Y);
 Yi11 = Yi(1:nx1,1:nx1);
 Yi12 = Yi(1:nx1,nx1+1:nx);
 Yi22 = Yi(nx1+1:nx,nx1+1:nx);
 X11 = X(:,1:nx1);
 X12 = X(:,nx1+1:nx);
 Z11 = X11 - Yi11;
 Z12 = X12 - Yi12;
 X22 = Z12'/Z11*Z12 + Yi22;
 X = [X11 X12; X12' X22];

 N = [eye(nx1); Z12'/Z11];
 M = -Y*[Z11;Z12'];
 N1 = N(1:nx1,:);
 N2 = N(nx1+1:nx,:);
 M1 = M(1:nx1,:);
 M2 = M(nx1+1:nx,:);
  
% Get state-space data for parm_i.
 sys = xtracti(vlpv,i,1);
 [A,Bd,Bu,Ce,Cy,Ded,Deu,Dyd,Dyu,r12inv,r21inv,q12,q21] = ...
	transfr(sys,ny,nu,nx2);
 Cy1 = Cy(1:ny1,:);
 Dy1d = Dyd(1:ny1,:);

 if DK1_FLAG
  D11 = Ded(1:ne1,1:nd1);
  D12 = Ded(1:ne1,nd1+1:nd);
  D21 = Ded(ne1+1:ne,1:nd1);
  D22 = Ded(ne1+1:ne,nd1+1:nd);
  dk1 = -D22-D21/(gamma^2 - D11'*D11)*D11'*D12;
  A = A + Bu*dk1*Cy1;
  Bd = Bd + Bu*dk1*Dy1d;
  Ce = Ce + Deu*dk1*Cy1;
  Ded = Ded + Deu*dk1*Dy1d;
 else
  dk1 = zeros(nu,ny1);
 end

 if DED_FLAG
  Dh = inv(1 - Ded*Ded'/(gamma^2));
  Dt = inv(1 - Ded*Ded'/(gamma^2));
  A = A + Bd*Ded'*Dh*Ce/(gamma^2);
  Bu = Bu + Bd*Ded'*Dh*Deu/(gamma^2);
  Cy = Cy + Dyd*Ded'*Dh*Ce/(gamma^2);
  F = -(Deu'*Dh*Deu)\(gamma*Bu'/Y + Deu'*Ce);
  L = -(X\Cy'*gamma + Bd*Dyd')*daug(Dy1d*Dt*Dy1d',eye(ny2));
 else
  F = -(gamma*Bu'/Y + Deu'*Ce);
  L = -(X\Cy'*gamma + Bd*Dyd');
 end

 akhat = -X*(Bd+L*Dyd)*Bd'/gamma - Ce'*(Ce+Deu*F)*Y/gamma;
 aktld = sel(A'+X*(A+Bu*F+L*Cy)*Y-akhat,1:nx1,':'); % + \dot{X}*Y
 dcbak = [F*Y; -N1\aktld]/([Y12 M1; Y22 M2]');
 ak = dcbak(nu+(1:nx1),ny2+(1:nx1));
 bk = N1\X(1:nx1,:)*L + [zeros(nx1,ny1) dcbak(nu+(1:nx1),1:ny2)];
 ck = dcbak(1:nu,ny2+(1:nx1));
 dk = [dk1, dcbak(1:nu,1:ny2)];

 dk = r12inv*dk*r21inv';
 ck = r12inv*ck;
 bk = bk*r21inv';

 if any(any(Dyu))
  ck = (eye(nu) + dk*Dyu)\ck;
  dk = (eye(nu) + dk*Dyu)\dk;
  ak = ak - bk*Dyu*ck;
  bk = bk - bk*Dyu*dk;
 end
	      
 ctrl = pck(ak,bk,ck,dk);
 clp = starp(sys,ctrl);
 vctrl((i-1)*sysrows+1:i*sysrows,:) = ctrl;
 vclp((i-1)*clprows+1:i*clprows,:) = clp;
   
end

if grid_pts > 1
 vctrl = vpck(vctrl,viv);
 vclp = vpck(vclp,viv);
end