miqp.m

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

function [xmin, fmin, flag, Extendedflag] = ...
          miqp(H, f, A, b, Aeq, beq, vartype, lb, ub, x0, Options)
%===============================================================================
% function [xmin, fmin, flag, Extendedflag] = ...
%           miqp(H, f, A, b, Aeq, beq, vartype, lb, ub, x0, Options)
%
% Title:   miqp.m
%
% Version: 1.06
% 
% Project: Development of a Mixed Integer Quadratic Program (MIQP) solver
%          for Matlab 
%
% Purpose: Solves the following MIQP
% 
%             min         0.5*x'H x + f' x
%             subject to  A   x <= b
%                         Aeq x  = beq
%              
%          if lb, ub are supplied, then the additional constraint
%
%             lb <= x <= ub
%
%          is imposed. If x0 is supplied, then x0 is taken as initial condition
%          for the search of the optimum xmin.
%          The variables indexed by vartype are binary, or more precisely
%
%             x(vartype) in {0,1} 
%
%          where {0,1} is understood as subset of the real numbers
%
% Requires:This function requires the availability of solvers for linear
%          programs (LP) and quadratic programs (QP). The following solvers are
%          supported:
%          
%          QP
%          --
%          quadprog.m : QP solver from matlab's optimization toolbox
%          qp.m       : Old version of quadprog.m
%          e04naf.m   : QP solver from NAG foundation toolbox
%          dantzgmp.m : QP solver from MPC toolbox
%
%          LP
%          --
%          linprog.m  : LP solver from matlab's optimization toolbox 
%          lp.m       : Old version of linprog.m
%          e04mbf.m   : LP solver from NAG foundation toolbox
%
%          Using e04naf.m requires also the presence of an m-function 
%          perfomring the simple matrix vector product H*x. Use for instance
%          the function qphess.m contained in the package.
%
%          If one or more solvers for LP and QP are not available, it is
%          recommended to force the use of the available solvers with the
%          parameter Options explained below.
%
% Authors: Alberto Bemporad   
%          Domenico Mignone  
%              
% History: 
% version  date         subject    
%              
% 0.1      1998.05.03   Preliminary Version 0.1, codename: miqp3_naf.m
% 1.02     2000.09.28   Public Beta Version 1.0
% 1.03     2000.10.27   Default initialization has been set to -inf instead of 0
% 1.04     2000.10.30   The options parameter for linprog and quadprog has been
%                       set as a field in the structure Options, namely
%                       Options.optimset. This variable can be set using the
%                       optimization toolbox routine "optimset"
% 1.05     2001.05.08   Systematic rounding of binary variables introduced
%                       according to the Flag Options.round 
%                       Routine qphess.m is not required anymore since it is
%                       created and deleted on the fly. To avoid deletion, see
%                           Options.deletefile
% 1.06     2001.05.09   Solver option qp_dantz added
%               
% Inputs:  mandatory arguments:
%
%          H, f    : parameters of the cost function     
%          A, b    : parameters defining the inequality constraints    
%
%          optional arguments:
%
%          Aeq, beq: parameters defining the equality constraints
%
%          vartype : vector specifying indices of binary variables 
%                    miqp.m supports 2 notations for this parameter:
%                    1.) vartype is a vector of positive integers
%                        the entries in vartype are the indices of those 
%                        variables constrained to be binary 
%                        the length is from 0 to size(H,1) (default [])
%                    2.) vartype is a vector of characters of length size(H,1)
%                        vartype(i) = 'C'  i-th variable is continuous
%                        vartype(i) = 'B'  i-th variable is binary: 0,1
%                        vartype(i) = 'I'  i-th variable is integer
%                    the option 'I' is not supported and is included for
%                    compatibility
%          lb      : Lower bounds on x (default -infinity)
%          ub      : Upper bounds on x (default +infinity) 
%          x0      : Initial condition (default 0)
%          Options : Variable in matlab's structure format, allowing to set
%                    various options for the solution algorithm. The following
%                    fields are considered. For detailed explanation of the
%                    various options see the documentation mentioned below.
%                    The mentioned default values are assumed either if Options
%                    is not defined, or if the fields below are not defined.
%          
%          Options.solver {'lp'|'linprog'|'lpnag'|'qp'|'quadprog'|'qpnag'|
%                          'qp_dantz'}
%
%                    'lp'       : LP solver from matlab's optimization toolbox
%                    'linprog'  : LP solver from matlab's optimization toolbox
%                    'lpnag'    : LP solver from matlab's NAG foundation toolbox
%                    'qp'       : QP solver from matlab's optimization toolbox
%                    'quadprog' : QP solver from matlab's optimization toolbox
%                    'qpnag'    : QP solver from matlab's NAG foundation toolbox
%                    'qp_dantz' : QP solver from matlab's MPC toolbox
%                    (default   : 'quadprog')
%                    see also Options.matrixtol below
%
%          Options.method {'depth'|'breadth'|'best'|'bestdepth'}
% 
%                    'depth'    : depth first tree exploring strategy
%                    'breadth'  : breadth first tree exploring strategy
%                    'best'     : best first tree exploring strategy
%                    'bestdepth': best first tree exploring strategy, norm. cost
%                    (default   : 'depth')
%
%          Options.branchrule {'first'|'max'|'min'}
%   
%                    'first'    : first free variable
%                    'max'      : variable, where the relaxed solution has the
%                                 largest distance to the next integer in {0,1}
%                    'min'      : variable, where the relaxed solution has the
%                                 smallest distance to the next integer in {0,1}
%                    (default   : 'first')
%          
%          Options.order {0|1}
%
%                    0          : problem, where the binary variable has been
%                                 set to 0 is solved first
%                    1          : problem, where the binary variable has been
%                                 set to 1 is solved first
%                    (default   : 0)
%                     
%          Options.verbose {0|1|2}
%
%                    0          : quiet
%                    1          : medium number of messages
%                    2          : high number of messages
%                    (default   : 0)
%
%          Options.maxqp {positive integer}
%
%                    maximum number of relaxed QPs (LPs) allowed to be solved
%                    (default   : inf)
%
%          Options.inftol {positive real}
%                    
%                    large number to be considered as infinity in constraints
%                    this is only used with the solvers from the NAG toolbox
%                    (default   : 1e8)
%
%          Options.matrixtol {nonnegative real}
%
%                    tolerance for recognizing that the MIQP is actually an MILP
%                    this option is only utilizable, if Options.solver is un-
%                    defined, the default solver for LPs is 'linprog'
%                    a problem is considered an MILP if
%                        max(svd(H)) <= matrixtol
%                    Set this parameter to inf in order to always ignore the 
%                    matrix H and always treat the problem as an MILP
%                    This parameter is also used to determine, whether the 
%                    problem can be solved with qp_dantz.m
%                    (default   : 1e-6)
%
%          Options.postol {nonnegative real}
%
%                    tolerance for recognizing H > 0, if any relaxed QP has
%                        cond(H) <= Options.postol
%                    then a warning message is produced for verbose >= 1
%                    to avoid the computation of the condition number for each
%                    relaxed QP, leave this field undefined
%                    (default   : no check, otherwise 1e-6 for values out of
%                                 range)
%
%          Options.integtol {nonnegative real}
%   
%                    tolerance to recognize integers
%                    (default   : 1e-4)
%
%          Options.maxQPiter {positive integer}
%          
%                    maximum number of iterations within each QP or LP 
%                    (default   : 1000)
%
%          Options.optimset {structure}
%                    
%                    optional arguments to be passed to linprog or quadprog
%                    (default   : optlinprog, where
%                                 optlinprog =optimset('LargeScale','off', ...
%                                                      'Display','off', ...
%                                                      'MaxIter',maxQPiter);   )
%          
%          Options.deletefile {0|1}
%
%                    if the auxiliary routine qphess.m had to be generated by
%                    miqp.m (only for the solver 'qpnag'), then this flag
%                    determines, whether the file is deleted at the end of the
%                    execution of miqp.m
%
%                    0: qphess.m is left in the current directory
%                    1: qphess.m is deleted
%                    (default : 1)
%
%          Options.round {0|1}
%
%                    0: relaxed variables, that fulfill the integer tolerances
%                       are not rounded to integer values
%                    1: as soon as a variable satisfies the integer tolerance
%                       it is rounded to the next integer value
%                    (default : 1)
%
%          Calling miqp.m without input arguments reports the version number
%
% Outputs: xmin:  minimizer of the MIQP
%          fmin:  minimum value of the cost function
%          flag:  characterization of solution according to the following list
%                   1 : integer feasible
%                   5 : feasible, but not integer feasible
%                   7 : infeasible, i.e. relaxed problem is infeasible
%                  11 : integer feasible, the limit maxqp of QPs has been
%                       reached, i.e. the solution might be suboptimal
%                  15 : feasible, but not integer feasible, the limit maxqp of
%                       QPs has been reached, i.e. did not search long enough
%                  -1 : the solution is unbounded
%
%          Extendedflag: Variable in matlab's structure format, which contains
%                        other informations about the MIQP or MILP. The
%                        following fields can be accessed in Extendedflag:
%
%          Extendedflag.QPiter: Number of relaxed QPs or LPs solved in the
%                               branch and bound tree
%          Extendedflag.time  : Time elapsed for running miqp.m
%          Extendedflag.optQP : Number of relaxed QP at which the optimal
%                               has been found
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
% Legal note:
%          This library is free software; you can redistribute it and/or
%          modify it under the terms of the GNU Lesser General Public
%          License as published by the Free Software Foundation; either
%          version 2.1 of the License, or (at your option) any later version.
%
%          This library is distributed in the hope that it will be useful,
%          but WITHOUT ANY WARRANTY; without even the implied warranty of
%          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%          Lesser General Public License for more details.
% 
%          You should have received a copy of the GNU Lesser General Public
%          License along with this library; if not, write to the 
%          Free Software Foundation, Inc., 
%          59 Temple Place, Suite 330, 
%          Boston, MA  02111-1307  USA
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%                               
% Documentation: A technical report presenting all the features and the main
%                implementation ideas of this solver can be found at:
%
%                http://www.aut.ee.ethz.ch/~hybrid/miqp/ 
%
% Contact:       Alberto Bemporad
%                Automatic Control Laboratory
%                ETH Zentrum
%                Zurich, Switzerland
%                [e-mail]
%
%                Domenico Mignone
%                Automatic Control Laboratory
%                ETH Zentrum
%                Zurich, Switzerland
%                [e-mail]
%
%                Comments and bug reports are highly appreciated
%
% (C) 1998-2000  Alberto Bemporad, Domenico Mignone
%                                    
%===============================================================================

tic; % for timing purposes

% ==============================================================================
% Argument verifications                                           
% ==============================================================================

if nargin == 0
   disp('Version 1.06')
   return
end   
                                     
error(nargchk(4,11,nargin));

% define optional arguments
% -------------------------

if nargin <= 10
   Options = [];
end      
if nargin <= 9
   x0 = [];
end
if nargin <= 8
   ub = [];
end
if nargin <= 7
   lb = [];
end
if nargin <= 6
   vartype = [];
end
if nargin <= 5
   beq = [];
end
if nargin <= 4
   Aeq = [];
end         

% specify defaults for all undefined parameters in Options
% --------------------------------------------------------

default_solver     = 'quadprog'; % specify here the default solver for QPs
default_method     = 'depth';    % specify here the default method
default_branchrule = 'first';    % specify here the default branchrule
default_order      = 0;          % specify here the default order
default_verbose    = 0;          % specify here the default verbose
default_maxqp      = inf;        % specify here the default maxqp
default_inftol     = 1e8;        % specify here the default inftol
default_matrixtol  = 1e-6;       % specify here the default matrixtol
default_postol     = 1e-6;       % specify here the default postol
default_integtol   = 1e-4;       % specify here the default integtol
default_maxQPiter  = 1000;       % specify here the default maxQPiter
default_solverlp   = 'linprog';  % specify here the default solver for LPs
default_deletefile = 1;          % specify here the default deletefile-flag
default_round      = 1;          % specify here the default round-flag
default_xxmax      = 1e7;        % specify here the default xmax for qp_dantz