0001 function solution = solveCobraMIQP(MIQPproblem,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 global CBT_MIQP_SOLVER;
0060 solver = CBT_MIQP_SOLVER;
0061
0062
0063 optParamNames = {'printLevel', 'saveInput', 'timeLimit'};
0064 parameters = '';
0065 if nargin ~=1
0066 if mod(length(varargin),2)==0
0067 for i=1:2:length(varargin)-1
0068 if ismember(varargin{i},optParamNames)
0069 parameters.(varargin{i}) = varargin{i+1};
0070 else
0071 error([varargin{i} ' is not a valid optional parameter']);
0072 end
0073 end
0074 elseif strcmp(varargin{1},'default')
0075 parameters = 'default';
0076 elseif isstruct(varargin{1})
0077 parameters = varargin{1};
0078 else
0079 display('Warning: Invalid number of parameters/values')
0080 solution=[];
0081 return;
0082 end
0083 end
0084 [printLevel, saveInput, timeLimit] = getCobraSolverParams('QP',optParamNames,parameters);
0085
0086 [A,b,F,c,lb,ub,csense,osense, vartype] = ...
0087 deal(MIQPproblem.A,MIQPproblem.b,MIQPproblem.F,MIQPproblem.c,MIQPproblem.lb,MIQPproblem.ub,...
0088 MIQPproblem.csense,MIQPproblem.osense, MIQPproblem.vartype);
0089
0090 t_start = clock;
0091 switch solver
0092
0093 case 'tomlab_cplex'
0094 if (~isempty(csense))
0095 b_L(csense == 'E') = b(csense == 'E');
0096 b_U(csense == 'E') = b(csense == 'E');
0097 b_L(csense == 'G') = b(csense == 'G');
0098 b_U(csense == 'G') = inf;
0099 b_L(csense == 'L') = -inf;
0100 b_U(csense == 'L') = b(csense == 'L');
0101 else
0102 b_L = b;
0103 b_U = b;
0104 end
0105 intVars = find((vartype == 'B') | (vartype == 'I'));
0106
0107 tomlabProblem = miqpAssign(osense*F, osense*c, A, b_L, b_U, lb, ub,[], ...
0108 intVars, [],[],[],'CobraMIQP');
0109 tomlabProblem.CPLEX.LogFile = 'MIQPproblem.log';
0110
0111
0112 PriLvl = printLevel;
0113
0114
0115 if ~isempty(saveInput)
0116 fileName = parameters.saveInput;
0117 if ~find(regexp(fileName,'.mat'))
0118 fileName = [fileName '.mat'];
0119 end
0120 display(['Saving MIQPproblem in ' fileName]);
0121 save(fileName,'MIQPproblem')
0122 end
0123 tomlabProblem.MIP.cpxControl.TILIM = timeLimit;
0124 tomlabProblem.MIP.cpxControl.THREADS = 1;
0125 Result = tomRun('cplex', tomlabProblem, PriLvl);
0126
0127 x = Result.x_k;
0128 f = osense*Result.f_k;
0129 stat = Result.Inform;
0130 if (stat == 1 ||stat == 101 || stat == 102)
0131 solStat = 1;
0132 elseif (stat == 3 || stat == 4)
0133 solStat = 0;
0134 elseif (stat == 103)
0135 solStat = 0;
0136 elseif (stat == 2 || stat == 118 || stat == 119)
0137 solStat = 2;
0138 elseif (stat == 106 || stat == 108 || stat == 110 || stat == 112 || stat == 114 || stat == 117)
0139 solStat = -1;
0140 elseif (stat >= 10)
0141 solStat = -1;
0142 else
0143 solStat = 3;
0144 end
0145
0146 case 'gurobi'
0147
0148
0149
0150
0151
0152
0153 clear opts
0154 if printLevel == 0
0155
0156
0157 opts.Display = 0;
0158 opts.DisplayInterval = 0;
0159 else
0160 opts.Display = 1;
0161 end
0162
0163
0164
0165 if (isempty(csense))
0166 clear csense
0167 csense(1:length(b),1) = '=';
0168 else
0169 csense(csense == 'L') = '<';
0170 csense(csense == 'G') = '>';
0171 csense(csense == 'E') = '=';
0172 csense = csense(:);
0173 end
0174
0175
0176
0177
0178
0179 [qrow,qcol,qval]=find(F);
0180 qrow=qrow'-1;
0181 qcol=qcol'-1;
0182 qval=0.5*qval';
0183
0184 opts.QP.qrow = int32(qrow);
0185 opts.QP.qcol = int32(qcol);
0186 opts.QP.qval = qval;
0187 opts.Method = 0;
0188 opts.Presolve = -1;
0189 opts.FeasibilityTol = 1e-6;
0190 opts.IntFeasTol = 1e-5;
0191 opts.OptimalityTol = 1e-6;
0192
0193
0194
0195 c = double(c);
0196 [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0197 csense,lb,ub,vartype,opts);
0198 if origStat==2
0199 stat = 1;
0200 elseif origStat==3
0201 stat = 0;
0202 elseif origStat==5
0203 stat = 2;
0204 elseif origStat==4
0205 stat = 0;
0206 else
0207 stat = -1;
0208 end
0209 solStat = stat;
0210 case 'gurobi5'
0211
0212
0213
0214 resultgurobi = struct('x',[],'objval',[],'pi',[]);
0215 clear params
0216 if printLevel == 0
0217 params.OutputFlag = 0;
0218 params.DisplayInterval = 1;
0219 else
0220 params.OutputFlag = 1;
0221 params.DisplayInterval = 5;
0222 end
0223
0224 params.Method = 0;
0225 params.Presolve = -1;
0226 params.IntFeasTol = 1e-5;
0227 params.FeasibilityTol = 1e-6;
0228 params.OptimalityTol = 1e-6;
0229
0230 if (isempty(MIQPproblem.csense))
0231 clear MIQPproblem.csense
0232 MIQPproblem.csense(1:length(b),1) = '=';
0233 else
0234 MIQPproblem.csense(MIQPproblem.csense == 'L') = '<';
0235 MIQPproblem.csense(MIQPproblem.csense == 'G') = '>';
0236 MIQPproblem.csense(MIQPproblem.csense == 'E') = '=';
0237 MIQPproblem.csense = MIQPproblem.csense(:);
0238 end
0239
0240 if MIQPproblem.osense == -1
0241 MIQPproblem.osense = 'max';
0242 else
0243 MIQPproblem.osense = 'min';
0244 end
0245
0246 MIQPproblem.vtype = vartype;
0247 MIQPproblem.Q = 0.5*sparse(MIQPproblem.F);
0248 MIQPproblem.modelsense = MIQPproblem.osense;
0249 [MIQPproblem.A,MIQPproblem.rhs,MIQPproblem.obj,MIQPproblem.sense] = deal(sparse(MIQPproblem.A),MIQPproblem.b,MIQPproblem.c,MIQPproblem.csense);
0250 resultgurobi = gurobi(MIQPproblem,params);
0251 solStat = resultgurobi.status;
0252 if strcmp(resultgurobi.status,'OPTIMAL')
0253 stat = 1;
0254
0255 if exist('resultgurobi.pi')
0256 [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0257 else
0258 [x,f] = deal(resultgurobi.x,resultgurobi.objval);
0259 end
0260 elseif strcmp(resultgurobi.status,'INFEASIBLE')
0261 stat = 0;
0262 elseif strcmp(resultgurobi.status,'UNBOUNDED')
0263 stat = 2;
0264 elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0265 stat = 0;
0266 else
0267 stat = -1;
0268 end
0269
0270 otherwise
0271 error(['Unknown solver: ' solver]);
0272 end
0273
0274 t = etime(clock, t_start);
0275
0276 solution.obj = f;
0277 solution.solver = solver;
0278 solution.stat = solStat;
0279 solution.origStat = stat;
0280 solution.time = t;
0281 solution.full = x;