solveCobraMILP

PURPOSE ^

solveCobraMILP Solve constraint-based MILP problems

SYNOPSIS ^

function solution = solveCobraMILP(MILPproblem,varargin)

DESCRIPTION ^

solveCobraMILP Solve constraint-based MILP problems

 solution = solveCobraMILP(MILPproblem,parameters)

INPUT
 MILPproblem
  A      LHS matrix
  b      RHS vector
  c      Objective coeff vector
  lb     Lower bound vector
  ub     Upper bound vector
  osense Objective sense (-1 max, +1 min)
  csense Constraint senses, a string containting the constraint sense for
         each row in A ('E', equality, 'G' greater than, 'L' less than).
  vartype Variable types ('C' continuous, 'I' integer, 'B' binary)
  x0      Initial solution

OPTIONAL INPUTS
 Optional parameters can be entered using parameters structure or as
 parameter followed by parameter value: i.e. ,'printLevel',3)

 parameters    Structure containing optional parameters.
  timeLimit    Global solver time limit
  intTol       Integrality tolerance
  relMipGapTol Relative MIP gap tolerance
  logFile      Log file (for CPLEX)
  printLevel    Printing level
               = 0    Silent (Default)
               = 1    Warnings and Errors
               = 2    Summary information 
               = 3    More detailed information
  saveInput    Saves LPproblem to filename specified in field. 
               i.e. parameters.saveInput = 'LPproblem.mat';
               Setting parameters = 'default' uses default setting set in
               getCobraSolverParameters.

 The solver is defined in the CBT_MILP_SOLVER global variable
 (set using changeCobraSolver). Solvers currently available are 
 'tomlab_cplex' and 'glpk'

OUTPUT
 solution Structure containing the following fields describing a MILP
          solution
  cont     Continuous solution
  int      Integer solution
  full     Full MILP solution vector
  obj      Objective value
  solver   Solver used to solve MILP problem
  stat     Solver status in standardized form (see below)
            1  Optimal solution found
            2  Unbounded solution
            0  Infeasible MILP
           -1  No integer solution exists
            3  Other problem (time limit etc, but integer solution exists)
  origStat Original status returned by the specific solver 
  time     Solve time in seconds


 Markus Herrgard 1/23/07
 Tim Harrington  05/18/12 Added support for the Gurobi 5.0 solver

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function solution = solveCobraMILP(MILPproblem,varargin)
0002 %solveCobraMILP Solve constraint-based MILP problems
0003 %
0004 % solution = solveCobraMILP(MILPproblem,parameters)
0005 %
0006 %INPUT
0007 % MILPproblem
0008 %  A      LHS matrix
0009 %  b      RHS vector
0010 %  c      Objective coeff vector
0011 %  lb     Lower bound vector
0012 %  ub     Upper bound vector
0013 %  osense Objective sense (-1 max, +1 min)
0014 %  csense Constraint senses, a string containting the constraint sense for
0015 %         each row in A ('E', equality, 'G' greater than, 'L' less than).
0016 %  vartype Variable types ('C' continuous, 'I' integer, 'B' binary)
0017 %  x0      Initial solution
0018 %
0019 %OPTIONAL INPUTS
0020 % Optional parameters can be entered using parameters structure or as
0021 % parameter followed by parameter value: i.e. ,'printLevel',3)
0022 %
0023 % parameters    Structure containing optional parameters.
0024 %  timeLimit    Global solver time limit
0025 %  intTol       Integrality tolerance
0026 %  relMipGapTol Relative MIP gap tolerance
0027 %  logFile      Log file (for CPLEX)
0028 %  printLevel    Printing level
0029 %               = 0    Silent (Default)
0030 %               = 1    Warnings and Errors
0031 %               = 2    Summary information
0032 %               = 3    More detailed information
0033 %  saveInput    Saves LPproblem to filename specified in field.
0034 %               i.e. parameters.saveInput = 'LPproblem.mat';
0035 %               Setting parameters = 'default' uses default setting set in
0036 %               getCobraSolverParameters.
0037 %
0038 % The solver is defined in the CBT_MILP_SOLVER global variable
0039 % (set using changeCobraSolver). Solvers currently available are
0040 % 'tomlab_cplex' and 'glpk'
0041 %
0042 %OUTPUT
0043 % solution Structure containing the following fields describing a MILP
0044 %          solution
0045 %  cont     Continuous solution
0046 %  int      Integer solution
0047 %  full     Full MILP solution vector
0048 %  obj      Objective value
0049 %  solver   Solver used to solve MILP problem
0050 %  stat     Solver status in standardized form (see below)
0051 %            1  Optimal solution found
0052 %            2  Unbounded solution
0053 %            0  Infeasible MILP
0054 %           -1  No integer solution exists
0055 %            3  Other problem (time limit etc, but integer solution exists)
0056 %  origStat Original status returned by the specific solver
0057 %  time     Solve time in seconds
0058 %
0059 %
0060 % Markus Herrgard 1/23/07
0061 % Tim Harrington  05/18/12 Added support for the Gurobi 5.0 solver
0062 
0063 %% Process options
0064 
0065 global CBT_MILP_SOLVER
0066 
0067 if (~isempty(CBT_MILP_SOLVER))
0068     solver = CBT_MILP_SOLVER;
0069 else
0070     error('No solver found.  Run changeCobraSolver');
0071 end
0072 
0073 optParamNames = {'intTol', 'relMipGapTol', 'timeLimit', ...
0074     'logFile', 'printLevel', 'saveInput', 'DATACHECK', 'DEPIND', ...
0075     'feasTol', 'optTol', 'absMipGapTol', 'NUMERICALEMPHASIS', 'EleNames', ... 
0076     'EqtNames', 'VarNames', 'EleNameFun', 'EqtNameFun', 'VarNameFun', ...
0077     'PbName', 'MPSfilename'};
0078 parameters = '';
0079 if nargin ~=1
0080     if mod(length(varargin),2)==0
0081         for i=1:2:length(varargin)-1
0082             if ismember(varargin{i},optParamNames)
0083                 parameters.(varargin{i}) = varargin{i+1};
0084             else
0085                 error([varargin{i} ' is not a valid optional parameter']);
0086             end
0087         end
0088     elseif strcmp(varargin{1},'default')
0089         parameters = 'default';
0090     elseif isstruct(varargin{1})
0091         parameters = varargin{1};
0092     else
0093         display('Warning: Invalid number of parameters/values')
0094         solution=[];
0095         return;
0096     end
0097 end
0098 
0099 %optional parameters
0100 [solverParams.intTol, solverParams.relMipGapTol, solverParams.timeLimit, ...
0101     solverParams.logFile, solverParams.printLevel, saveInput, ...
0102     solverParams.DATACHECK, solverParams.DEPIND, solverParams.feasTol, ...
0103     solverParams.optTol, solverParams.absMipGapTol, ...
0104     solverParams.NUMERICALEMPHASIS] = ...
0105     getCobraSolverParams('MILP',optParamNames(1:12), parameters);
0106 
0107 %Save Input if selected
0108 if ~isempty(saveInput)
0109     fileName = parameters.saveInput;
0110     if ~find(regexp(fileName,'.mat'))
0111         fileName = [fileName '.mat'];
0112     end
0113     display(['Saving MILPproblem in ' fileName]);
0114     save(fileName,'MILPproblem')
0115 end
0116 
0117 % Defaults in case the solver does not return anything
0118 %x = [];
0119 xInt = [];
0120 xCont = [];
0121 %stat = -99;
0122 %solStat = -99;
0123 
0124 [A,b,c,lb,ub,csense,osense,vartype,x0] = ...
0125     deal(MILPproblem.A,MILPproblem.b,MILPproblem.c,MILPproblem.lb,MILPproblem.ub,...
0126     MILPproblem.csense,MILPproblem.osense,MILPproblem.vartype,MILPproblem.x0);
0127 
0128 if any(~(vartype == 'C' | vartype == 'B' | vartype == 'I'))
0129     display ('vartype not C or B or I:  Assuming C');
0130     vartype(vartype ~= 'C' & vartype ~= 'I'& vartype ~= 'B') = 'C';
0131 end
0132 
0133 t_start = clock;
0134 switch solver
0135     
0136     case 'glpk'
0137 %% glpk
0138 
0139         % Set up problem
0140         if (isempty(csense))
0141             clear csense
0142             csense(1:length(b),1) = 'S';
0143         else
0144             csense(csense == 'L') = 'U';
0145             csense(csense == 'G') = 'L';
0146             csense(csense == 'E') = 'S';
0147             csense = columnVector(csense);
0148         end
0149         params.msglev = solverParams.printLevel;
0150         params.tmlim = solverParams.timeLimit;
0151         
0152         %whos csense vartype
0153         csense = char(csense);
0154         vartype = char(vartype);
0155         %whos csense vartype
0156         
0157         % Solve problem
0158         [x,f,stat,extra] = glpk(c,A,b,lb,ub,csense,vartype,osense,params);
0159         % Handle solution status reports
0160         if (stat == 5)
0161             solStat = 1; % optimal
0162         elseif (stat == 6)
0163             solStat = 2; % unbounded
0164         elseif (stat == 4) 
0165             solStat = 0; % infeasible
0166         
0167         elseif (stat == 171)
0168             solStat = 1; % Opt integer within tolerance
0169         elseif (stat == 173)
0170             solStat = 0; % Integer infeas
0171         elseif (stat == 184)
0172             solStat = 2; % Unbounded
0173         elseif (stat == 172)
0174             solStat = 3; % Other problem, but integer solution exists
0175         else
0176             solStat = -1; % No integer solution exists
0177         end
0178 
0179    case 'gurobi'
0180         % Free academic licenses for the Gurobi solver can be obtained from
0181         % http://www.gurobi.com/html/academic.html
0182         %
0183         % The code below uses Gurobi Mex to interface with Gurobi. It can be downloaded from
0184         % http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_A_MATLAB_interface_for_Gurobi
0185 
0186         clear opts % Use the default parameter settings
0187         if solverParams.printLevel == 0
0188            % Version v1.10 of Gurobi Mex has a minor bug. For complete silence
0189            % Remove Line 736 of gurobi_mex.c: mexPrintf("\n");
0190            opts.Display = 0;
0191            opts.DisplayInterval = 0;
0192         else
0193            opts.Display = 1;
0194         end
0195         
0196         %minimum intTol for gurobi = 1e-9
0197         if solverParams.intTol<1e-9, solverParams.intTol=1e-9; end
0198         
0199         opts.TimeLimit=solverParams.timeLimit;
0200         opts.MIPGap = solverParams.relMipGapTol;
0201         opts.IntFeasTol = solverParams.intTol;
0202         opts.FeasibilityTol = solverParams.feasTol;
0203         opts.OptimalityTol = solverParams.optTol;
0204 
0205         if (isempty(csense))
0206             clear csense
0207             csense(1:length(b),1) = '=';
0208         else
0209             csense(csense == 'L') = '<';
0210             csense(csense == 'G') = '>';
0211             csense(csense == 'E') = '=';
0212             csense = csense(:);
0213         end
0214         %gurobi_mex doesn't automatically cast logicals to doubles
0215     c = double(c);
0216         [x,f,stat,output] = gurobi_mex(c,osense,sparse(A),b, ...
0217                                              csense,lb,ub,vartype,opts);
0218         if stat == 2
0219            solStat = 1; % Optimal solutuion found
0220         elseif stat == 3
0221            solStat = 0; % Infeasible
0222         elseif stat == 5
0223            solStat = 2; % Unbounded
0224         elseif stat == 4
0225            solStat = 0; % Gurobi reports infeasible *or* unbounded
0226         else
0227            solStat = -1; % Solution not optimal or solver problem
0228         end
0229         
0230  case 'gurobi5'
0231         %% gurobi 5
0232         % Free academic licenses for the Gurobi solver can be obtained from
0233         % http://www.gurobi.com/html/academic.html
0234         resultgurobi = struct('x',[],'objval',[]);
0235         MILPproblem.A = deal(sparse(MILPproblem.A));
0236         clear params            % Use the default parameter settings
0237         
0238         if solverParams.printLevel == 0 
0239            params.OutputFlag = 0;
0240            params.DisplayInterval = 1;
0241         else
0242            params.OutputFlag = 1;
0243            params.DisplayInterval = 5;
0244         end
0245 
0246         params.TimeLimit = solverParams.timeLimit;
0247         params.MIPGap = solverParams.relMipGapTol;
0248         
0249         if solverParams.intTol <= 1e-09
0250             params.IntFeasTol = 1e-09;
0251         else
0252             params.IntFeasTol = solverParams.intTol;
0253         end
0254         
0255         params.FeasibilityTol = solverParams.feasTol;
0256         params.OptimalityTol = solverParams.optTol;
0257         
0258         if (isempty(csense))
0259             clear csense
0260             csense(1:length(b),1) = '=';
0261         else
0262             csense(csense == 'L') = '<';
0263             csense(csense == 'G') = '>';
0264             csense(csense == 'E') = '=';
0265             MILPproblem.csense = csense(:);
0266         end
0267     
0268         if osense == -1
0269             MILPproblem.osense = 'max';
0270         else
0271             MILPproblem.osense = 'min';
0272         end
0273         
0274         MILPproblem.vtype = vartype;
0275         MILPproblem.modelsense = MILPproblem.osense;
0276         [MILPproblem.A,MILPproblem.rhs,MILPproblem.obj,MILPproblem.sense] = deal(sparse(MILPproblem.A),MILPproblem.b,double(MILPproblem.c),MILPproblem.csense);
0277         resultgurobi = gurobi(MILPproblem,params);
0278         
0279         stat = resultgurobi.status;
0280         if strcmp(resultgurobi.status,'OPTIMAL')
0281            solStat = 1; % Optimal solution found
0282            [x,f] = deal(resultgurobi.x,resultgurobi.objval);
0283         elseif strcmp(resultgurobi.status,'INFEASIBLE')
0284            solStat = 0; % Infeasible
0285         elseif strcmp(resultgurobi.status,'UNBOUNDED')
0286            solStat = 2; % Unbounded
0287         elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0288            solStat = 0; % Gurobi reports infeasible *or* unbounded
0289         else
0290            solStat = -1; % Solution not optimal or solver problem
0291         end
0292         
0293     case 'tomlab_cplex'
0294 %% CPLEX through tomlab
0295         if (~isempty(csense))
0296             b_L(csense == 'E') = b(csense == 'E');
0297             b_U(csense == 'E') = b(csense == 'E');
0298             b_L(csense == 'G') = b(csense == 'G');
0299             b_U(csense == 'G') = inf;
0300             b_L(csense == 'L') = -inf;
0301             b_U(csense == 'L') = b(csense == 'L');
0302         elseif isfield(MILPproblem, 'b_L') && isfield(MILPproblem, 'b_U')
0303             b_L = MILPproblem.b_L;
0304             b_U = MILPproblem.b_U;
0305         else
0306             b_L = b;
0307             b_U = b;
0308         end
0309         intVars = (vartype == 'B') | (vartype == 'I');
0310         %intVars
0311         %pause;
0312         tomlabProblem = mipAssign(osense*c,A,b_L,b_U,lb,ub,x0,'CobraMILP',[],[],intVars);
0313         
0314         % Set parameters for CPLEX
0315         tomlabProblem.MIP.cpxControl.EPINT = solverParams.intTol;
0316         tomlabProblem.MIP.cpxControl.EPGAP = solverParams.relMipGapTol;
0317         tomlabProblem.MIP.cpxControl.TILIM = solverParams.timeLimit;
0318         tomlabProblem.CPLEX.LogFile = solverParams.logFile;
0319         tomlabProblem.PriLev = solverParams.printLevel;
0320         tomlabProblem.MIP.cpxControl.THREADS = 1; % by default use only one thread
0321         
0322 
0323         % Strict numerical tolerances
0324         tomlabProblem.MIP.cpxControl.DATACHECK = solverParams.DATACHECK;
0325         tomlabProblem.MIP.cpxControl.DEPIND = solverParams.DEPIND;
0326         tomlabProblem.MIP.cpxControl.EPRHS = solverParams.feasTol;
0327         tomlabProblem.MIP.cpxControl.EPOPT = solverParams.optTol;
0328         tomlabProblem.MIP.cpxControl.EPAGAP = solverParams.absMipGapTol;
0329         tomlabProblem.MIP.cpxControl.NUMERICALEMPHASIS = solverParams.NUMERICALEMPHASIS;
0330         % Set initial solution
0331         tomlabProblem.MIP.xIP = x0;
0332      
0333         % Set up callback to print out intermediate solutions
0334         % only set this up if you know that you actually need these
0335         % results.  Otherwise do not specify intSolInd and contSolInd
0336         global cobraIntSolInd;
0337         global cobraContSolInd;
0338         if(~isfield(MILPproblem, 'intSolInd'))
0339             MILPproblem.intSolInd = [];
0340         else
0341             tomlabProblem.MIP.callback(14) = 1;
0342         end
0343         cobraIntSolInd = MILPproblem.intSolInd;
0344         if(~isfield(MILPproblem, 'contSolInd'))
0345             MILPproblem.contSolInd = [];
0346         end
0347         cobraContSolInd = MILPproblem.contSolInd;
0348         tomlabProblem.MIP.callbacks = [];
0349         tomlabProblem.PriLevOpt = 0;        
0350         
0351         
0352         % Solve problem
0353         Result = tomRun('cplex', tomlabProblem);
0354         
0355         % Get results
0356         x = Result.x_k;
0357         f = osense*Result.f_k;
0358         stat = Result.Inform;
0359         if (stat == 101 || stat == 102)
0360             solStat = 1; % Opt integer within tolerance
0361         elseif (stat == 103)
0362             solStat = 0; % Integer infeas
0363         elseif (stat == 118 || stat == 119)
0364             solStat = 2; % Unbounded
0365         elseif (stat == 106 || stat == 106 || stat == 108 || stat == 110 || stat == 112 || stat == 114 || stat == 117)
0366             solStat = -1; % No integer solution exists
0367         else
0368             solStat = 3; % Other problem, but integer solution exists
0369         end
0370     case 'mps'
0371         %% BuildMPS
0372         % This calls buildMPS and generates a MPS format description of the
0373         % problem as the result
0374         % Build MPS Author: Bruno Luong
0375         % Interfaced with CobraToolbox by Richard Que (12/18/09)
0376         display('Solver set to MPS. This function will output an MPS matrix string for the MILP problem');
0377         %Get optional parameters
0378         [EleNames,EqtNames,VarNames,EleNameFun,EqtNameFun,VarNameFun,PbName,MPSfilename] = ...
0379             getCobraSolverParams('LP',{'EleNames','EqtNames','VarNames','EleNameFun','EqtNameFun','VarNameFun','PbName','MPSfilename'},parameters);
0380         %split A matrix for L and E csense
0381         Ale = A(csense=='L',:);
0382         ble = b(csense=='L');
0383         Aeq = A(csense=='E',:);
0384         beq = b(csense=='E');
0385         %create index of integer and binary variables
0386         intIndex = find(vartype=='I');
0387         binaryIndex = find(vartype=='B');
0388         
0389         %%%%Adapted from BuildMPS%%%%%
0390         [neq nvar]=size(Aeq);
0391         nle=size(Ale,1);
0392         if isempty(EleNames)
0393             EleNames=arrayfun(EleNameFun,(1:nle),'UniformOutput', false);
0394         end
0395         if isempty(EqtNames)
0396             EqtNames=arrayfun(EqtNameFun,(1:neq),'UniformOutput', false);
0397         end
0398         if isempty(VarNames)
0399             VarNames=arrayfun(VarNameFun,(1:nvar),'UniformOutput', false);
0400         end
0401         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0402 
0403         [solution] = BuildMPS(Ale, ble, Aeq, beq, c, lb, ub, PbName,'MPSfilename',MPSfilename,'EqtNames',EqtNames,'VarNameFun',VarNameFun,'Integer',intIndex,'Binary',binaryIndex);
0404   
0405         return
0406     otherwise
0407         error(['Unknown solver: ' solver]);
0408 end
0409 t = etime(clock, t_start);
0410 
0411 %% Store results
0412 if ~strcmp(solver,'mps')
0413     if (~isempty(x))
0414         %xInt = x(MILPproblem.intSolInd);
0415         %xCont = x(MILPproblem.contSolInd);
0416         xInt = x(vartype == 'B' | vartype == 'I');
0417         xCont = x(vartype == 'C');
0418     end
0419     
0420     solution.cont = xCont;
0421     solution.int = xInt;
0422     solution.obj = f;
0423     solution.solver = solver;
0424     solution.stat = solStat;
0425     solution.origStat = stat;
0426     solution.time = t;
0427     solution.full = x;
0428     if(isfield(MILPproblem, 'intSolInd'))
0429         solution.intInd = MILPproblem.intSolInd;
0430     end
0431 end
0432

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003