0001 function solution = solveCobraMILP(MILPproblem,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
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
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
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
0118
0119 xInt = [];
0120 xCont = [];
0121
0122
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
0138
0139
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
0153 csense = char(csense);
0154 vartype = char(vartype);
0155
0156
0157
0158 [x,f,stat,extra] = glpk(c,A,b,lb,ub,csense,vartype,osense,params);
0159
0160 if (stat == 5)
0161 solStat = 1;
0162 elseif (stat == 6)
0163 solStat = 2;
0164 elseif (stat == 4)
0165 solStat = 0;
0166
0167 elseif (stat == 171)
0168 solStat = 1;
0169 elseif (stat == 173)
0170 solStat = 0;
0171 elseif (stat == 184)
0172 solStat = 2;
0173 elseif (stat == 172)
0174 solStat = 3;
0175 else
0176 solStat = -1;
0177 end
0178
0179 case 'gurobi'
0180
0181
0182
0183
0184
0185
0186 clear opts
0187 if solverParams.printLevel == 0
0188
0189
0190 opts.Display = 0;
0191 opts.DisplayInterval = 0;
0192 else
0193 opts.Display = 1;
0194 end
0195
0196
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
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;
0220 elseif stat == 3
0221 solStat = 0;
0222 elseif stat == 5
0223 solStat = 2;
0224 elseif stat == 4
0225 solStat = 0;
0226 else
0227 solStat = -1;
0228 end
0229
0230 case 'gurobi5'
0231
0232
0233
0234 resultgurobi = struct('x',[],'objval',[]);
0235 MILPproblem.A = deal(sparse(MILPproblem.A));
0236 clear params
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;
0282 [x,f] = deal(resultgurobi.x,resultgurobi.objval);
0283 elseif strcmp(resultgurobi.status,'INFEASIBLE')
0284 solStat = 0;
0285 elseif strcmp(resultgurobi.status,'UNBOUNDED')
0286 solStat = 2;
0287 elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0288 solStat = 0;
0289 else
0290 solStat = -1;
0291 end
0292
0293 case 'tomlab_cplex'
0294
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
0311
0312 tomlabProblem = mipAssign(osense*c,A,b_L,b_U,lb,ub,x0,'CobraMILP',[],[],intVars);
0313
0314
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;
0321
0322
0323
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
0331 tomlabProblem.MIP.xIP = x0;
0332
0333
0334
0335
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
0353 Result = tomRun('cplex', tomlabProblem);
0354
0355
0356 x = Result.x_k;
0357 f = osense*Result.f_k;
0358 stat = Result.Inform;
0359 if (stat == 101 || stat == 102)
0360 solStat = 1;
0361 elseif (stat == 103)
0362 solStat = 0;
0363 elseif (stat == 118 || stat == 119)
0364 solStat = 2;
0365 elseif (stat == 106 || stat == 106 || stat == 108 || stat == 110 || stat == 112 || stat == 114 || stat == 117)
0366 solStat = -1;
0367 else
0368 solStat = 3;
0369 end
0370 case 'mps'
0371
0372
0373
0374
0375
0376 display('Solver set to MPS. This function will output an MPS matrix string for the MILP problem');
0377
0378 [EleNames,EqtNames,VarNames,EleNameFun,EqtNameFun,VarNameFun,PbName,MPSfilename] = ...
0379 getCobraSolverParams('LP',{'EleNames','EqtNames','VarNames','EleNameFun','EqtNameFun','VarNameFun','PbName','MPSfilename'},parameters);
0380
0381 Ale = A(csense=='L',:);
0382 ble = b(csense=='L');
0383 Aeq = A(csense=='E',:);
0384 beq = b(csense=='E');
0385
0386 intIndex = find(vartype=='I');
0387 binaryIndex = find(vartype=='B');
0388
0389
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
0412 if ~strcmp(solver,'mps')
0413 if (~isempty(x))
0414
0415
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