0001 function solution = solveCobraQP(QPproblem,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 global CBT_QP_SOLVER;
0062
0063 if (~isempty(CBT_QP_SOLVER))
0064 solver = CBT_QP_SOLVER;
0065 else
0066 error('No solver found');
0067 end
0068
0069 optParamNames = {'printLevel','saveInput'};
0070 parameters = '';
0071 if nargin ~=1
0072 if mod(length(varargin),2)==0
0073 for i=1:2:length(varargin)-1
0074 if ismember(varargin{i},optParamNames)
0075 parameters.(varargin{i}) = varargin{i+1};
0076 else
0077 error([varargin{i} ' is not a valid optional parameter']);
0078 end
0079 end
0080 elseif strcmp(varargin{1},'default')
0081 parameters = 'default';
0082 elseif isstruct(varargin{1})
0083 parameters = varargin{1};
0084 else
0085 display('Warning: Invalid number of parameters/values')
0086 solution=[];
0087 return;
0088 end
0089 end
0090
0091
0092 x = [];
0093 y = [];
0094 w = [];
0095 f = [];
0096 xInt = [];
0097 xCont = [];
0098 stat = -99;
0099 solStat = -99;
0100
0101
0102 [printLevel, saveInput] = getCobraSolverParams('QP',optParamNames,parameters);
0103
0104 [A,b,F,c,lb,ub,csense,osense] = ...
0105 deal(QPproblem.A,QPproblem.b,QPproblem.F,QPproblem.c,QPproblem.lb,QPproblem.ub,...
0106 QPproblem.csense,QPproblem.osense);
0107
0108 t_start = clock;
0109 switch solver
0110
0111 case 'tomlab_cplex'
0112 if (~isempty(csense))
0113 b_L(csense == 'E') = b(csense == 'E');
0114 b_U(csense == 'E') = b(csense == 'E');
0115 b_L(csense == 'G') = b(csense == 'G');
0116 b_U(csense == 'G') = inf;
0117 b_L(csense == 'L') = -inf;
0118 b_U(csense == 'L') = b(csense == 'L');
0119 else
0120 b_L = b;
0121 b_U = b;
0122 end
0123 tomlabProblem = qpAssign(F,osense*c,A,b_L,b_U,lb,ub,[],'CobraQP');
0124
0125
0126 tomlabProblem.PriLvl=printLevel;
0127 tomlabProblem.MIP.cpxControl.QPMETHOD = 1;
0128 tomlabProblem.MIP.cpxControl.THREADS = 1;
0129
0130
0131 if ~isempty(saveInput)
0132 fileName = saveInput;
0133 if ~find(regexp(fileName,'.mat'))
0134 fileName = [fileName '.mat'];
0135 end
0136 display(['Saving QPproblem in ' fileName]);
0137 save(fileName,'QPproblem')
0138 end
0139
0140 Result = tomRun('cplex', tomlabProblem);
0141 x = Result.x_k;
0142 f = osense*Result.f_k;
0143 origStat = Result.Inform;
0144 if (origStat == 1)
0145 stat = 1;
0146 elseif (origStat == 3 || origStat == 4)
0147 stat = 0;
0148 elseif (origStat == 2)
0149 stat = 2;
0150 elseif (origStat >= 10)
0151 stat = -1;
0152 else
0153 stat = 3;
0154 end
0155
0156 case 'cplex_direct'
0157
0158
0159
0160
0161
0162
0163
0164
0165 solution=solveCobraLPCPLEX(QPproblem,printLevel,[],[],[],minNorm);
0166
0167 case 'qpng'
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177 ctype=csense;
0178 ctype(('G'==csense))='L';
0179 ctype(('E'==csense))='E';
0180 ctype(('L'==csense))='U';
0181
0182 x0=ones(size(QPproblem.A,2),1);
0183
0184 [x, f, y, info] = qpng (QPproblem.F, QPproblem.c*QPproblem.osense, full(QPproblem.A), QPproblem.b, ctype, QPproblem.lb, QPproblem.ub, x0);
0185
0186 f = 0.5*x'*QPproblem.F*x + c'*x;
0187
0188 w=[];
0189
0190 if (info.status == 0)
0191 stat = 1;
0192 elseif (info.status == 1)
0193 stat = 0;
0194 else
0195 stat = -1;
0196 end
0197 origStat=info.status;
0198
0199 case 'mosek'
0200 if (~isempty(csense))
0201 b_L(csense == 'E') = b(csense == 'E');
0202 b_U(csense == 'E') = b(csense == 'E');
0203 b_L(csense == 'G') = b(csense == 'G');
0204 b_U(csense == 'G') = inf;
0205 b_L(csense == 'L') = -inf;
0206 b_U(csense == 'L') = b(csense == 'L');
0207 else
0208 b_L = b;
0209 b_U = b;
0210 end
0211
0212 if printLevel>0
0213 cmd='minimize';
0214 else
0215 cmd='minimize echo(0)';
0216 end
0217
0218
0219
0220
0221
0222 [res] = mskqpopt(F,osense*c,A,b_L,b_U,lb,ub,[],cmd);
0223
0224 if isempty(res)
0225 stat=3;
0226 else
0227 if isfield(res,'sol')
0228 origStat=res.sol.itr.solsta;
0229 if strcmp(res.sol.itr.prosta,'PRIMAL_AND_DUAL_FEASIBLE') && (strcmp(res.sol.itr.solsta,'OPTIMAL') || strcmp(res.sol.itr.solsta,'NEAR_OPTIMAL'))
0230 stat=1;
0231
0232 x = res.sol.itr.xx;
0233 f = 0.5*x'*F*x + c'*x;
0234
0235
0236 y=res.sol.itr.y;
0237
0238
0239 w=res.sol.itr.slx - res.sol.itr.sux;
0240 else
0241 stat=3;
0242 end
0243 else
0244 stat=3;
0245 origStat=res.rmsg;
0246 end
0247 end
0248
0249
0250
0251
0252
0253
0254 case 'pdco'
0255
0256
0257
0258
0259
0260
0261
0262 [nMet,nRxn]=size(A);
0263 d1=ones(nRxn,1)*1e-4;
0264
0265 d1(c~=0)=0;
0266 d2=1e-5;
0267 options = pdcoSet;
0268
0269 x0 = ones(nRxn,1);
0270 y0 = ones(nMet,1);
0271 z0 = ones(nRxn,1);
0272 xsize = 1000;
0273 zsize = 1000;
0274 options.Method=2;
0275 options.MaxIter=100;
0276 options.Print=printLevel;
0277
0278 pdObjHandle = @(x) QPObj(x);
0279
0280 [x,y,w,inform,PDitns,CGitns,time] = ...
0281 pdco(pdObjHandle,A,b,lb,ub,d1,d2,options,x0,y0,z0,xsize,zsize);
0282 f= c'*x + 0.5*x'*F*x;
0283
0284
0285
0286
0287
0288 if (inform == 0)
0289 stat = 1;
0290 elseif (inform == 1 || inform == 2 || inform == 3)
0291 stat = 0;
0292 else
0293 stat = -1;
0294 end
0295 origStat=inform;
0296
0297 case 'gurobi'
0298
0299
0300
0301
0302
0303
0304 clear opts
0305 if printLevel == 0
0306
0307
0308 opts.Display = 0;
0309 opts.DisplayInterval = 0;
0310 else
0311 opts.Display = 1;
0312 end
0313
0314
0315
0316 if (isempty(csense))
0317 clear csense
0318 csense(1:length(b),1) = '=';
0319 else
0320 csense(csense == 'L') = '<';
0321 csense(csense == 'G') = '>';
0322 csense(csense == 'E') = '=';
0323 csense = csense(:);
0324 end
0325
0326
0327
0328
0329
0330 [qrow,qcol,qval]=find(F);
0331 qrow=qrow'-1;
0332 qcol=qcol'-1;
0333 qval=0.5*qval';
0334
0335 opts.QP.qrow = int32(qrow);
0336 opts.QP.qcol = int32(qcol);
0337 opts.QP.qval = qval;
0338 opts.Method = 0;
0339 opts.Presolve = -1;
0340 opts.FeasibilityTol = 1e-6;
0341 opts.IntFeasTol = 1e-5;
0342 opts.OptimalityTol = 1e-6;
0343
0344
0345
0346 c = double(c);
0347 [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0348 csense,lb,ub,[],opts);
0349 if origStat==2
0350 stat = 1;
0351 elseif origStat==3
0352 stat = 0;
0353 elseif origStat==5
0354 stat = 2;
0355 elseif origStat==4
0356 stat = 0;
0357 else
0358 stat = -1;
0359 end
0360
0361 case 'gurobi5'
0362
0363
0364
0365 resultgurobi = struct('x',[],'objval',[],'pi',[]);
0366 clear params
0367 if printLevel == 0
0368 params.OutputFlag = 0;
0369 params.DisplayInterval = 1;
0370 else
0371 params.OutputFlag = 1;
0372 params.DisplayInterval = 5;
0373 end
0374
0375 params.Method = 0;
0376 params.Presolve = -1;
0377 params.IntFeasTol = 1e-5;
0378 params.FeasibilityTol = 1e-6;
0379 params.OptimalityTol = 1e-6;
0380
0381
0382 if (isempty(QPproblem.csense))
0383 clear QPproblem.csense
0384 QPproblem.csense(1:length(b),1) = '=';
0385 else
0386 QPproblem.csense(QPproblem.csense == 'L') = '<';
0387 QPproblem.csense(QPproblem.csense == 'G') = '>';
0388 QPproblem.csense(QPproblem.csense == 'E') = '=';
0389 QPproblem.csense = QPproblem.csense(:);
0390 end
0391
0392 if QPproblem.osense == -1
0393 QPproblem.osense = 'max';
0394 else
0395 QPproblem.osense = 'min';
0396 end
0397
0398 QPproblem.Q = 0.5*sparse(QPproblem.F);
0399 QPproblem.modelsense = QPproblem.osense;
0400 [QPproblem.A,QPproblem.rhs,QPproblem.obj,QPproblem.sense] = deal(sparse(QPproblem.A),QPproblem.b,double(QPproblem.c),QPproblem.csense);
0401 resultgurobi = gurobi(QPproblem,params);
0402 origStat = resultgurobi.status;
0403 if strcmp(resultgurobi.status,'OPTIMAL')
0404 stat = 1;
0405 [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0406 elseif strcmp(resultgurobi.status,'INFEASIBLE')
0407 stat = 0;
0408 elseif strcmp(resultgurobi.status,'UNBOUNDED')
0409 stat = 2;
0410 elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0411 stat = 0;
0412 else
0413 stat = -1;
0414 end
0415
0416
0417 otherwise
0418 error(['Unknown solver: ' solver]);
0419 end
0420
0421 t = etime(clock, t_start);
0422 solution.obj = f;
0423 solution.solver = solver;
0424 solution.stat = stat;
0425 solution.origStat = origStat;
0426 solution.time = t;
0427 solution.full = x;
0428 solution.dual = y;
0429 solution.rcost = w;
0430
0431
0432
0433 function [obj,grad,hess] = QPObj(x)
0434 obj = c'*x + 0.5*x'*F*x;
0435 grad = c + F*x;
0436 hess = F;
0437 end
0438 end