miqp.m

由matlab开发的hybrid系统的描述语言

% handle the different fields of 'Options'
% ----------------------------------------

if isfield(Options,'solver')
   switch Options.solver
      case {'lp','linprog','lpnag','qp','quadprog','qpnag','qp_dantz'}
         solver = Options.solver;
      otherwise
         warning('solver is not implemented: will take default solver')
         solver = default_solver;
   end   
else
   solver = default_solver;   
end

if isfield(Options,'method')
   switch Options.method
      case {'depth','breadth','best','bestdepth'}
         method = Options.method;
      otherwise
         warning('method is not implemented: will take default method')
         method = default_method;
   end   
else
   method = default_method;   
end

if isfield(Options,'branchrule')
   switch Options.branchrule
      case {'first','max','min'}
         branchrule = Options.branchrule;
      otherwise
         warning('branchrule is not implemented: will take default branchrule')
         branchrule = default_branchrule;
   end   
else
   branchrule = default_branchrule;   
end

if isfield(Options,'order')
   switch Options.order
      case {0,1}
         order = Options.order;
      case {'0','1'}
         order = str2num(Options.order); 
      otherwise
         warning('order is not implemented: will take default order')
         order = default_order;
   end   
else
   order = default_order;   
end

if isfield(Options,'verbose')
   switch Options.verbose
      case {0,1,2}
         verbose = Options.verbose;
      case {'0','1','2'}
         verbose = str2num(Options.verbose); 
      otherwise
         warning('verbose is not implemented: will take default verbose')
         verbose = default_verbose;
   end   
else
   verbose = default_verbose;   
end

if isfield(Options,'maxqp')
   if Options.maxqp >= 1
      maxqp = floor(Options.maxqp);
   else
      warning('maxqp is negative or 0: will take default maxqp')
      maxqp = default_maxqp;
   end   
else
   maxqp = default_maxqp;   
end

if isfield(Options,'inftol')
   if Options.inftol > 0
      inftol = Options.inftol;
   else
      warning('inftol is negative or 0: will take default inftol')
      inftol = default_inftol;
   end   
else
   inftol = default_inftol;   
end

if isfield(Options,'matrixtol')
   if Options.matrixtol >= 0
      matrixtol = Options.matrixtol;
   else
      warning('matrixtol is negative: will take default matrixtol')
      matrixtol = default_matrixtol;
   end   
else
   matrixtol = default_matrixtol;   
end

if isfield(Options,'postol')
   if Options.postol >= 0
      postol = Options.postol;
   else
      warning('postol is negative: will take default postol')
      postol = default_postol;
   end   
else
   postol = default_postol;   
end

if isfield(Options,'integtol')
   if Options.integtol >= 0
      integtol = Options.integtol;
   else
      warning('integtol is negative: will take default integtol')
      integtol = default_integtol;
   end   
else
   integtol = default_integtol;   
end

if isfield(Options,'maxQPiter')
   if Options.maxQPiter >= 1
      maxQPiter = floor(Options.maxQPiter);
   else
      warning('maxQPiter is negative or 0: will take default maxQPiter')
      maxQPiter = default_maxQPiter;
   end   
else
   maxQPiter = default_maxQPiter;   
end

if isfield(Options,'deletefile')
   if Options.deletefile == 1
      deletefile = 1;
   elseif Options.deletefile == 0
      deletefile = 0;
   else
      warning('unallowed value for deletefile: will take default')  
      deletefile = default_deletefile;
   end   
else
   deletefile = default_deletefile;  
end

if isfield(Options,'round')
   if Options.round == 1
      rounding = 1;
   elseif Options.round == 0
      rounding = 0;
   else
      warning('unallowed value for round: will take default')  
      rounding = default_round;
   end   
else
   rounding = default_round;  
end

% checking dimensions  
% -------------------

if size(H,1) ~= size(H,2)
   error('H is not square')
end

if ~isfield(Options,'solver') 
   if max(svd(H)) < matrixtol
      if verbose >= 1
         warning('This is a MILP')
      end
      solver = default_solverlp;  % default solver for LPs
   end
elseif strcmp(Options.solver,'qp_dantz')
   if max(svd(H)) < matrixtol
      error('MILPs cannot be solved with qp_dantz')
   end
end      
      
if ~isymm(H)
   H=0.5*(H+H');
   if verbose >= 1
      warning('H is not symmetric: replaced by its symmetric part')
   end   
end

if (size(f,1) ~= 1) & (size(f,2) ~= 1)
   error('f must be a vector')
end

if (size(b,1) ~= 1) & (size(b,2) ~= 1)
   error('b must be a vector')
end

f = f(:);   % f and b are column vectors
b = b(:);

nx = size(H,1);

if size(A,1) ~= size(b,1)
   error('A and b have incompatible dimensions')
end

if size(A,2) ~= nx
   error('A and H have incompatible dimensions')
end

if size(Aeq,1) ~= size(beq,1)
   error('Aeq and beq have incompatible dimensions')
end

if (size(Aeq,2) ~= nx)&(~isempty(Aeq))
   error('Aeq and H have incompatible dimensions')
end

lb      = lb(:);
ub      = ub(:);
x0      = x0(:);      
vartype = vartype(:);

if ischar(vartype)
   if length(vartype) ~= nx
      error('wrong dimension of vartype as character array')
   end      
   Ccons = find(vartype=='C');
   Bcons = find(vartype=='B');
   Icons = find(vartype=='I');
   if length([Ccons(:); Bcons(:); Icons(:)]) ~=  nx;
      error('wrong entries in vartype')
   end   
   if ~isempty(Icons)
      error('specifications for integer variables are not supported')
   end
   % deletes the variable in character syntax and uses the syntax with vector
   % of indices instead
   vartype = Bcons(:);
end   

if size(lb,1)~=nx & ~isempty(lb)
   error('lb has wrong dimensions')
end
if size(ub,1)~=nx & ~isempty(ub)
   error('ub has wrong dimensions')
end
if size(x0,1)~=nx & ~isempty(x0)
   error('x0 has wrong dimensions')
end
if max(vartype) > nx
   error('largest index in vartype is out of range')
end
if min(vartype) < 1
   error('smallest index in vartype is out of range')
end
if find(diff(sort(vartype)) == 0)
   error('binary variables are multiply defined')
end   
if length(vartype) > nx
   error('too many entries in vartype') 
end 
if floor(vartype) ~= vartype
   error('fractional number in vartype not allowed')
end     

% Define default values for lb,ub,x0
% ------------------------------------

if isempty(lb),
   lb =-inf*ones(nx,1);
end
if isempty(ub),
   ub = inf*ones(nx,1);
end
if isempty(x0),
   x0  = zeros(nx,1);
end
if find(ub-lb < 0)
   error('infeasible constraints specified in ub and lb')
end   

cont          = (1:nx)';
cont(vartype) = [];       % Indices of continuous variables

% zstar    denotes the best value for the cost function so far          
% QPiter   counts how many times the QP or LP algorithm is invoked            
% nivar    is the total number of integer variables                     
% xstar    denotes the current optimal vector

zstar   = inf;        
xstar   = inf*ones(nx,1); 
QPiter  = 0;
nivar   = length(vartype);
flag    = 7;		  % by default it is infeasible
optQP   = 0;		  % number of QP where optimum is found

% Initialize parameters for the QP routines
% -----------------------------------------

if strcmp(solver,'qpnag')
   cold   = 1;
   wu     = sqrt(eps);
   orthog = 1;
   lpf    = 0;
   if ~exist('qphess.m')
      fid = fopen('qphess.m','w');
      fprintf(fid,'function [hx] = qphess(n,nrowh,ncolh,jthcol,hess,x)\n' );
      fprintf(fid,'hx=hess*x;\n' );
      term = fclose(fid);
      justcreated = 1;
   else
      justcreated = 0;
   end      
end   

if strcmp(solver,'qp_dantz')
   if (min(lb) > -inf) & (max(ub) < inf)
      xxmax = abs(max(min(lb),max(ub)))*ones(nx,1);
   else
      xxmax = default_xxmax*ones(nx,1);  
   end 
   
   if ~exist('dantzgmp')
      error('the option qp_dantz requires the MPC toolbox')
   end
else
   xxmax = default_xxmax*ones(nx,1);
end   

if strcmp(solver,'lpnag')
   % no parameters
end

if strcmp(solver,'linprog')
   if isfield(Options,'optimset')
      optlinprog = Options.optimset;
   else
      optlinprog = optimset('LargeScale','off','Display','off', ...
                            'MaxIter',maxQPiter);
   end			    
end

if strcmp(solver,'quadprog')
   if isfield(Options,'optimset')
      optquadprog = Options.optimset;
   else      
      optquadprog = optimset('LargeScale','off','Display','off', ...
                             'MaxIter',maxQPiter);