panmip.m

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

%              
%              min(x)
%              s.t.    -x <= 0 
%                    -40x <= 0.04
%
%          With these numbers "linprog" correctly returns a value of 1 for 
%          "exitflag", while quadprog still returns -1. This is also in 
%          contraction with the documentation of "quadprog", which is supposed
%          to call "linprog" in case of a hessian equal to zero
%
%          The matlab code for this experiment follows
%
%                clear
%
%                A = [-1;  -400   ];
%                B = [ 0;     0.4 ];
%                [x1a, f1a, ex1a]  =  linprog(1, A, B);
%                [x2a, l2a, how2a] =       lp(1, A, B);
%                [x3a, f3a, ex3a]  = quadprog(0, 1, A, B);
%
%                A2  = [-1;  -40   ];
%                B2  = [ 0;    0.04 ];
%                [x1b, f1b, ex1b]  =  linprog(1, A2, B2);
%                [x2b, l2b, how2b] =       lp(1, A2, B2);
%                [x3b, f3b, ex3b]  = quadprog(0, 1, A2, B2);
%
%                xvalues   = [x1a, x2a, x3a, x1b, x2b, x3b]  % should all be 0
%                exitflags = [ex1a,ex3a,ex1b,ex3b]           % should all be > 0
%                how2a, how2b    % should both be 'ok'
% 
%          The results are 
%          xvalues =
%          1.0e-05 *
%          -0.1419     0    0.0091   -0.0000   0   -0.0024
%          exitflags =
%          -1    -1     1     -1
%
%          This error was reported to Mathworks and occured with Matlab 5.3 or
%          below. It should not occur for Matlab 6 and optimization toolbox 2.1
%
%          01.02.2001
%          ----------
%          If one tries to solve an LP in CPLEX with the command MIP, CPLEX
%          does not solve the problem reporting that it is not an MIP
%
%          01.02.2001
%          ----------
%          In some degenerate cases, CPLEX writes the solution using a syntax
%          (a keyword sequence) which is different from the one coded in the
%          handling of the solver-option "cplex_file"
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 
% "I get strange results"
%
% Do you obtain strange, unexplicable results while running this routine?
% According to the solver you are currently using, the following list gives you
% some hints about which options to tune first in such cases.
%
% miqp
% ----
%
% get infeasibilities: - Avoid using Options.miqp.solver = 'quadprog';
%                        or Options.miqp.solver = 'linprog';
%                        Beware, these solvers are set as a default. If you did
%                        not specify any solver, add e.g. 
%                        
%                        Options.miqp.solver = 'lpnag'; 
%                         
% cplex_file
% ----------
%
% get infeasibilities: - Turn off presolver adding
%
%                        Options.cplex_file.presolve = 0; 
%                                    
%===============================================================================

tic; % for timing purposes

% ==============================================================================
% ==============================================================================
% Platform-Dependent settings
% ==============================================================================
% ==============================================================================

% Directory where Baron is installed
% currently this is only the UNIX directory

BaronDirectory    = '/home/hybrid/mip/baron';

% CPLEX pertinent directories
% these are the DOS directories, it is assumed, that under UNIX cplex is an
% executable from the console

IlogLicense       = 'C:\ILOG\ILM\ACCESS.ILM';
CplexPath         = 'C:\ILOG\CPLEX70\bin\msvc6\stat_mta';

% XPRESS-MP pertinent directories
% this is the DOS directory, it is assumed, that under UNIX xpress is an
% executable from the console

xpresslocation = 'c:\XPRESSMP\CONSOLE\mp-opt.exe '; % where xpress is located

% platform dependent definitions can be found by searching for the string:
% PLATFORM DEPENDENT   

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

if nargin == 0
   disp('Version 1.35')
   xopt='1.35';
   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
   rangevar = [];
end
if nargin <= 4
   ctype = [];
end      

if ~isfield(Options,'solver')
   Options.solver = 'miqp';   % Specification of default solver
end

if isfield(Options,'bignumber')
   bignumber = Options.bignumber;
else
   bignumber = inf;   
end

% Check on x0ii item in struct Options is after test on x0

Extendedflag = [];

% ==============================================================================
% Checking dimensions and defining defaults
% ==============================================================================

if ~isymmm(H)
   H=0.5*(H+H');
   warning('H is not symmetric: replaced by its symmetric part')
end

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

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

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

nx = size(H,1);
nc = size(C,1);

if size(C,1) ~= size(d,1)
   error('C and d have incompatible dimensions')
end

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

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

if isempty(ctype)
   ctype = char('L'*ones(1,nc));
end

if isempty(vartype)
   vartype = char('C'*ones(1,nx));
end

if isempty(rangevar)
   rangevar = zeros(1,nc);
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

% Define default values for lb, ub, x0

if isempty(lb),
   lb = -bignumber*ones(nx,1);
end
if isempty(ub),
   ub = bignumber*ones(nx,1);
end
if isempty(x0),
   x0  = zeros(nx,1);
   x0_was_provided = 0;
else
   x0_was_provided = 1;
end
if isfield(Options,'x0ii')
   x0ii = Options.x0ii;
   x0ii = x0ii(:);
   if size(x0ii,1)~=nx 
      error('x0ii has wrong dimensions')
   end
   else
      x0ii = ones(nx,1);
end

if ~isempty(find(ub-lb < 0))
   error('infeasible constraints specified in ub and lb')
end   

n = nx; % number of Optimization variables

% Define default values for output arguments xopt and fopt

x_opt_init = inf*ones(nx,1);
f_opt_init = inf;

% ==============================================================================
% Choice of Solver
% ==============================================================================

switch Options.solver

   case {'miqp'}
   %------------   
      if ~isfield(Options,'miqp')
         Options.miqp = [];
      end
      
      [A, b, Aeq, beq] = ctype2mat(C, d, ctype, rangevar);
      
      eval([ '[xopt, fopt, flag, Extendedflag]= ' Options.solver...
             '(H, f, A, b, Aeq, beq, vartype, lb, ub, x0, Options.miqp); '])

      if (flag == 7) | (flag == -1)
         xopt = x_opt_init;
         fopt = f_opt_init;
      end                         
                 
   case {'fletcher'}
   %----------------   
      error('this solver did not pass the quality test yet:contact the authors')
   
      if ~isfield(Options,'fletcher')
         Options.fletcher = [];
      end      
      if isfield(Options.fletcher,'huge')
         huge = Options.fletcher.huge;
      else 
         huge = 1e10;
      end
      if isfield(Options.fletcher,'epsilon')
         epsilon = Options.fletcher.epsilon;
      else 
         epsilon = 1e10;
      end
      if isfield(Options.fletcher,'priority')         
         priority  = Options.fletcher.priority;
      else
         priority  = ones(size(H,1),1);
      end          
      if isfield(Options.fletcher,'iprint')         
         iprint  = Options.fletcher.iprint;
      else
         iprint  = -1;
      end        
      
      vlb = lb(:);
      vub = ub(:);
      
      [A, b] = ctype2matCompact(C, d, ctype, rangevar);
      ivar   = vartype2ivar(vartype, n, 'printerror');
      [lb,ub]= add01(ivar, lb, ub);
              
      a_dense     = [f(:)'; A]; % linear term + constraint matrix
      nvar        = size(H,1);  % total number of vars
      nlin        = size(C,1);  % total number of linear constraints
      x0          = x0(:);
      problemtype = 'MIQP';