lmi1cd.m

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

% [Dc,ubnd,ubnd_pk,xopt] = lmi1cd(M,w,struct,solv_options,gam_init,x_init)
%
% Compute constant-scaling upper bound for mu, i.e., solve for constant Dc 
% that minimizes ||Dc*M(jw_i)*Dc^(-1)||_inf over all i (frequency grid)
%
% Inputs: 
% M 		SYSTEM to be analyzed (must be square)
% w		vector of frequency points
% struct  	structure of the TV uncertainty
% 		The rows of struct are given as:
%		[n 1] for parameters repeated n times
%		[n 0] for n by n blocks
% NOTE: blk2strt.m can be used to convert a blk description used for mu.m
%        into a struct description used here
%
% OPTIONAL INPUTS: -------------
%
%  solv_options - is a 1x5 vector that specifies options for LMI solvers
%
%     NOTE: solv_options=[] or not defined gives the default values
%           to specify some of the items and leave some at the default
%           use 0 for the elements to be left at default
%
%		the first element is not used
%     max_its=0;         % max # of total iterations
%     feas_rad=0;        % feasability radius of decision vars
%     per_reduct_its=0;  % # of its done w/ less than 1% relative reduction
%     trace_opt=0;       % if 0, show trace, otherwise don't
%  then,
%		  solv_options=[0,max_its,feas_rad,per_reduct_its,trace_opt];
%
%  gam_init - initial guess at gamma (only used if x_init is specified)
%
%  x_init - initial guess at the vector of decision variables that solve the LMI
%           If this is specified without gam_init, gam_init is calculated
%
% Outputs:
% Dc:		constant scaling matrix associated with struct 
% ubnd,ubnd_pk: VARYING mu upper bound at each frequency and its max over w
% xopt:		vector of optimal LMI decision variables 

function [Dc,ubnd,ubnd_pk,xopt] = ...
	lmi1cd(M,w,struct,solv_options,gam_init,x_init);

%
%  Check for input errors
%
if(nargin<3 | nargin>6);
  disp(' [Dc,ubnd,ubnd_pk,xopt] = lmi1cd(M,w,struct,solv_options,gam_init,x_init)');
  return
end

g=M;
omega=w;

if ~exist('solv_options')
  solv_options=[];
end

if isempty(solv_options)
  max_its=0;         % max # of total iterations
  feas_rad=0;        % feasability radius of decision vars
  per_reduct_its=0;  % this is the # of its done w/ less than 1% relative reduct
  trace_opt=0;       % if 0, show trace, otherwise don't
  solv_options=[0,max_its,0,per_reduct_its,trace_opt];
end

if ~exist('x_init')
  gam_init=[];
  gamsq_init=[];
  x_init=[];
end

if isempty(x_init)       % then an IC has been provided
  if isempty(gam_init)   % this is because might have a x_init w/o gam_init
    gam_init=25;         % note: after lmi is defined check to make sure this
  end                    %       is big enough
end

%%% check that the inputs are correct form/type/etc
[type,nin,nout,ns] = minfo(g);
if type ~= 'syst';
  error('Not a system matrix')
end
if nin ~= nout;
  error('Not a square system')
end
if max(real(spoles(g))) >= 0;
  error('Not a stable system')
end
if min(size(omega)) ~= 1;
  error('The frequencies must be given as a vector')
end
if size(struct,2) ~= 2 | min(struct(:,2)) < 0 | max(struct(:,2)) > 1;
  error('Structure of the parameters is wrong')
end

ts_cpu=cputime;

%  Frequency response of the plant and its conjugate transpose
omega = sort(omega);
g_s = frsp(g,omega);
gstar_s = vcjt(g_s);
nfreq = max(size(omega));
nth = sum(struct(:,1));

% Define LMI variables
setlmis([])
qtheta = lmivar(1,struct);
lmis1 = getlmis;

setlmis(lmis1);
q = lmivar(3,daug(decinfo(lmis1,qtheta),decinfo(lmis1,qtheta)));
lmis1 = getlmis;

setlmis(lmis1);
% Add LMI terms
lmiterm([-1 1 1 qtheta],1,1);
for i = 1:nfreq;
 g_i = xtracti(g_s,i,1);
 gs_i = xtracti(gstar_s,i,1);
 bigg = [real(g_i) imag(g_i); -imag(g_i) real(g_i)];
 biggs = [real(gs_i) imag(gs_i); -imag(gs_i) real(gs_i)];
 lmiterm([1+i 1 1 q],biggs,bigg);
 lmiterm([-(1+i) 1 1 q],1,1);
end

lmis1 = getlmis;
% NOTE: there are nfreq+1 LMIs and the last nfreq of them contain a gamma term

% if IC given check gam_init to make sure starting point is feasible
if(isempty(x_init)==0)                % then an IC has been provided
  gamsq_init_min=tinitial(lmis1,x_init,2);
  gam_init_min=sqrt(gamsq_init_min);
  disp(' ');
  disp(['With the given initial values of the decision variables gamma = ',...
        num2str(gam_init_min),';']);
  disp(' ');
  if(gam_init<gam_init_min)
    gam_init=50*gam_init_min;   % NOTE: gam_init_min is feasable but this 
                                %           seems to converge faster
    disp(['gam_init has been adjusted, gam_init = ',num2str(gam_init)]);
    disp(' ');
  end
  gamsq_init=gam_init^2;
end

no_dec_vars=decnbr(lmis1);
disp(['the number of decision variables = ',num2str(no_dec_vars)]);
nlmis=lminbr(lmis1);
disp(['the number of lmis = ',num2str(nlmis)]);
nlmidisc=max(size(lmis1));
disp(['the lmisys vector is ',num2str(nlmidisc),' x 1']);

if(solv_options(5)==0)
  disp(' ** The current value of gamma^2 will be displayed below **');
  disp(' ');
end

if isempty(x_init)
 [gamsq_opt1,xopt] = gevp(lmis1,nfreq,solv_options);
else
 [gamsq_opt1,xopt] = gevp(lmis1,nfreq,solv_options,gamsq_init,x_init);
end
gam1_pk = sqrt(gamsq_opt1);
qthetaopt = dec2mat(lmis1,xopt,qtheta);
dtheta = sqrtm(qthetaopt);

Dc=dtheta;
ubnd = vnorm(frsp(mmult(dtheta,M,minv(dtheta)),w));
ubnd_pk = gam1_pk;