0001 function solution = solveCobraLP(LPproblem, 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
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085 global CBTLPSOLVER
0086 if (~isempty(CBTLPSOLVER))
0087 solver = CBTLPSOLVER;
0088 else
0089 error('No solver found. call changeCobraSolver(solverName)');
0090 end
0091 optParamNames = {'minNorm','printLevel','primalOnly','saveInput', ...
0092 'feasTol','optTol','EleNames','EqtNames','VarNames','EleNameFun', ...
0093 'EqtNameFun','VarNameFun','PbName','MPSfilename'};
0094 parameters = '';
0095 if nargin ~=1
0096 if mod(length(varargin),2)==0
0097 for i=1:2:length(varargin)-1
0098 if ismember(varargin{i},optParamNames)
0099 parameters.(varargin{i}) = varargin{i+1};
0100 else
0101 error([varargin{i} ' is not a valid optional parameter']);
0102 end
0103 end
0104 elseif strcmp(varargin{1},'default')
0105 parameters = 'default';
0106 elseif isstruct(varargin{1})
0107 parameters = varargin{1};
0108 else
0109 display('Warning: Invalid number of parameters/values')
0110 solution=[];
0111 return;
0112 end
0113 end
0114 [minNorm, printLevel, primalOnlyFlag, saveInput, feasTol, optTol] = ...
0115 getCobraSolverParams('LP',optParamNames(1:6),parameters);
0116
0117
0118
0119 if ~isempty(saveInput)
0120 fileName = parameters.saveInput;
0121 if ~find(regexp(fileName,'.mat'))
0122 fileName = [fileName '.mat'];
0123 end
0124 display(['Saving LPproblem in ' fileName]);
0125 save(fileName,'LPproblem')
0126 end
0127
0128
0129 [A,b,c,lb,ub,csense,osense] = deal(LPproblem.A,LPproblem.b,LPproblem.c,LPproblem.lb,LPproblem.ub,LPproblem.csense,LPproblem.osense);
0130
0131
0132
0133
0134
0135
0136 f = [];
0137 x = [];
0138 y = [];
0139 w = [];
0140 origStat = -99;
0141 stat = -99;
0142
0143 t_start = clock;
0144 switch solver
0145
0146 case 'glpk'
0147 params.msglev = printLevel;
0148 params.tolbnd = feasTol;
0149 params.toldj = optTol;
0150 if (isempty(csense))
0151 clear csense
0152 csense(1:length(b),1) = 'S';
0153 else
0154 csense(csense == 'L') = 'U';
0155 csense(csense == 'G') = 'L';
0156 csense(csense == 'E') = 'S';
0157 csense = columnVector(csense);
0158 end
0159
0160 b=full(b);
0161 [x,f,y,w,stat,origStat] = solveGlpk(c,A,b,lb,ub,csense,osense,params);
0162
0163 case {'lindo_new','lindo_old'}
0164
0165 if (strcmp(solver,'lindo_new'))
0166
0167 [f,x,y,w,s,origStat] = solveCobraLPLindo(A,b,c,csense,lb,ub,osense,primalOnlyFlag,false);
0168
0169 if (origStat == 1 || origStat == 2)
0170 stat = 1;
0171 elseif (origStat == 4)
0172 stat = 2;
0173 elseif (origStat == 3 || origStat == 6)
0174 stat = 0;
0175 else
0176 stat = -1;
0177 end
0178 else
0179
0180 [f,x,y,w,s,origStat] = solveCobraLPLindo(A,b,c,csense,lb,ub,osense,primalOnlyFlag,true);
0181
0182 if (origStat == 2 || origStat == 3)
0183 stat = 1;
0184 elseif (origStat == 5)
0185 stat = 2;
0186 elseif (origStat == 4 || origStat == 6)
0187 stat = 0;
0188 else
0189 stat = -1;
0190 end
0191 end
0192
0193
0194 case 'lp_solve'
0195
0196 if (isempty(csense))
0197 [f,x,y,origStat] = lp_solve(c*(-osense),A,b,zeros(size(A,1),1),lb,ub);
0198 f = f*(-osense);
0199 else
0200 e(csense == 'E') = 0;
0201 e(csense == 'G') = 1;
0202 e(csense == 'L') = -1;
0203 [f,x,y,origStat] = lp_solve(c*(-osense),A,b,e,lb,ub);
0204 f = f*(-osense);
0205 end
0206
0207 if (origStat == 0)
0208 stat = 1;
0209 elseif (origStat == 3)
0210 stat = 2;
0211 elseif (origStat == 2)
0212 stat = 0;
0213 else
0214 stat = -1;
0215 end
0216 s = [];
0217 w = [];
0218 case 'mosek'
0219
0220
0221
0222
0223 switch printLevel
0224 case 0
0225 options.Display='off';
0226 case 1
0227 options.Display='final';
0228 case 2
0229 options.Display='iter';
0230 otherwise
0231
0232 options = optimset;
0233 end
0234
0235 if (isempty(csense))
0236 [x,f,origStat,output,lambda] = linprog(c*osense,[],[],A,b,lb,ub,[],options);
0237 else
0238 Aeq = A(csense == 'E',:);
0239 beq = b(csense == 'E');
0240 Ag = A(csense == 'G',:);
0241 bg = b(csense == 'G');
0242 Al = A(csense == 'L',:);
0243 bl = b(csense == 'L');
0244 clear A;
0245 A = [Al;-Ag];
0246 clear b;
0247 b = [bl;-bg];
0248 [x,f,origStat,output,lambda] = linprog(c*osense,A,b,Aeq,beq,lb,ub,[],options);
0249 end
0250 y = [];
0251 if (origStat > 0)
0252 stat = 1;
0253 f = f*osense;
0254 y = lambda.eqlin;
0255 elseif (origStat < 0)
0256 stat = 0;
0257 else
0258 stat = -1;
0259 end
0260
0261 case 'gurobi'
0262
0263
0264
0265
0266
0267
0268
0269 clear opts
0270 if printLevel == 0
0271
0272
0273 opts.Display = 0;
0274 opts.DisplayInterval = 0;
0275 else
0276 opts.Display = 1;
0277 end
0278
0279 opts.FeasibilityTol = feasTol;
0280 opts.OptimalityTol = optTol;
0281
0282 if (isempty(csense))
0283 clear csense
0284 csense(1:length(b),1) = '=';
0285 else
0286 csense(csense == 'L') = '<';
0287 csense(csense == 'G') = '>';
0288 csense(csense == 'E') = '=';
0289 csense = csense(:);
0290 end
0291
0292 c = double(c);
0293 [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0294 csense,lb,ub,[],opts);
0295 if origStat==2
0296 stat = 1;
0297 elseif origStat==3
0298 stat = 0;
0299 elseif origStat==5
0300 stat = 2;
0301 elseif origStat==4
0302 stat = 0;
0303 else
0304 stat = -1;
0305 end
0306
0307 case 'gurobi5'
0308
0309
0310
0311 resultgurobi = struct('x',[],'objval',[],'pi',[]);
0312 LPproblem.A = deal(sparse(LPproblem.A));
0313 clear params
0314
0315 if printLevel == 0
0316 params.OutputFlag = 0;
0317 params.DisplayInterval = 1;
0318 else
0319 params.OutputFlag = 1;
0320 params.DisplayInterval = 5;
0321 end
0322
0323 params.FeasibilityTol = feasTol;
0324 params.OptimalityTol = optTol;
0325
0326 if (isempty(LPproblem.csense))
0327 clear LPproblem.csense
0328 LPproblem.csense(1:length(b),1) = '=';
0329 else
0330 LPproblem.csense(LPproblem.csense == 'L') = '<';
0331 LPproblem.csense(LPproblem.csense == 'G') = '>';
0332 LPproblem.csense(LPproblem.csense == 'E') = '=';
0333 LPproblem.csense = LPproblem.csense(:);
0334 end
0335
0336 if LPproblem.osense == -1
0337 LPproblem.osense = 'max';
0338 else
0339 LPproblem.osense = 'min';
0340 end
0341
0342 LPproblem.modelsense = LPproblem.osense;
0343 [LPproblem.rhs,LPproblem.obj,LPproblem.sense] = deal(LPproblem.b,double(LPproblem.c),LPproblem.csense);
0344 resultgurobi = gurobi(LPproblem,params);
0345
0346 if strcmp(resultgurobi.status,'OPTIMAL')
0347 stat = 1;
0348 [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0349 elseif strcmp(resultgurobi.status,'INFEASIBLE')
0350 stat = 0;
0351 elseif strcmp(resultgurobi.status,'UNBOUNDED')
0352 stat = 2;
0353 elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0354 stat = 0;
0355 else
0356 stat = -1;
0357 end
0358
0359 case 'matlab'
0360
0361 if (isempty(csense))
0362 [x,f,origStat,output,lambda] = linprog(c*osense,[],[],A,b,lb,ub);
0363 else
0364 Aeq = A(csense == 'E',:);
0365 beq = b(csense == 'E');
0366 Ag = A(csense == 'G',:);
0367 bg = b(csense == 'G');
0368 Al = A(csense == 'L',:);
0369 bl = b(csense == 'L');
0370 clear A;
0371 A = [Al;-Ag];
0372 clear b;
0373 b = [bl;-bg];
0374 [x,f,origStat,output,lambda] = linprog(c*osense,A,b,Aeq,beq,lb,ub);
0375 end
0376 y = [];
0377 if (origStat > 0)
0378 stat = 1;
0379 f = f*osense;
0380 y = lambda.eqlin;
0381 elseif (origStat < 0)
0382 stat = 0;
0383 else
0384 stat = -1;
0385 end
0386
0387 case 'tomlab_cplex'
0388
0389 if (~isempty(csense))
0390 b_L(csense == 'E') = b(csense == 'E');
0391 b_U(csense == 'E') = b(csense == 'E');
0392 b_L(csense == 'G') = b(csense == 'G');
0393 b_U(csense == 'G') = 1e6;
0394 b_L(csense == 'L') = -1e6;
0395 b_U(csense == 'L') = b(csense == 'L');
0396 else
0397 b_L = b;
0398 b_U = b;
0399 end
0400 tomlabProblem = lpAssign(osense*c,A,b_L,b_U,lb,ub);
0401
0402
0403
0404
0405 tomlabProblem.optParam = optParamDef('cplex',tomlabProblem.probType);
0406 tomlabProblem.QP.F = [];
0407 tomlabProblem.PriLevOpt = printLevel;
0408
0409
0410 if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis)
0411 tomlabProblem.MIP.basis = LPproblem.basis;
0412 end
0413
0414
0415 tomlabProblem.MIP.cpxControl.EPRHS = feasTol;
0416 tomlabProblem.MIP.cpxControl.EPOPT = optTol;
0417
0418
0419 Result = cplexTL(tomlabProblem);
0420
0421
0422 x = Result.x_k;
0423 f = osense*sum(tomlabProblem.QP.c.*Result.x_k);
0424
0425
0426 origStat = Result.Inform;
0427 w = Result.v_k(1:length(lb));
0428 y = Result.v_k((length(lb)+1):end);
0429 basis = Result.MIP.basis;
0430 if (origStat == 1)
0431 stat = 1;
0432 elseif (origStat == 3)
0433 stat = 0;
0434 elseif (origStat == 2 || origStat == 4)
0435 stat = 2;
0436 else
0437 stat = -1;
0438 end
0439 case 'cplex_direct'
0440
0441
0442
0443
0444
0445
0446
0447 if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis)
0448 LPproblem.LPBasis = LPproblem.basis;
0449 end
0450 [solution LPprob] = solveCobraLPCPLEX(LPproblem,printLevel,1,[],[],minNorm);
0451 solution.basis = LPprob.LPBasis;
0452 solution.solver = solver;
0453
0454 case 'lindo'
0455 error('Solver type lindo is obsolete - use lindo_new or lindo_old instead');
0456 case 'pdco'
0457
0458
0459
0460
0461
0462
0463
0464 [nMet,nRxn]=size(LPproblem.A);
0465 x0 = ones(nRxn,1);
0466 y0 = ones(nMet,1);
0467 z0 = ones(nRxn,1);
0468
0469
0470
0471
0472
0473 d1=0;
0474 d2=1e-6;
0475 options = pdcoSet;
0476 options.FeaTol = 1e-12;
0477 options.OptTol = 1e-12;
0478
0479
0480
0481
0482
0483
0484 xsize = 1000;
0485 zsize = 10000;
0486
0487 options.Method=2;
0488 options.MaxIter=100;
0489 options.Print=printLevel;
0490 [x,y,w,inform,PDitns,CGitns,time] = ...
0491 pdco(osense*c*10000,A,b,lb,ub,d1,d2,options,x0,y0,z0,xsize,zsize);
0492 f= c'*x;
0493
0494
0495
0496
0497
0498 if (inform == 0)
0499 stat = 1;
0500 elseif (inform == 1 || inform == 2 || inform == 3)
0501 stat = 0;
0502 else
0503 stat = -1;
0504 end
0505 origStat=inform;
0506 case 'mps'
0507
0508
0509
0510
0511
0512 display('Solver set to MPS. This function will output an MPS matrix string for the LP problem');
0513
0514 [EleNames,EqtNames,VarNames,EleNameFun,EqtNameFun,VarNameFun,PbName,MPSfilename] = ...
0515 getCobraSolverParams('LP',{'EleNames','EqtNames','VarNames','EleNameFun','EqtNameFun','VarNameFun','PbName','MPSfilename'},parameters);
0516
0517 Ale = A(csense=='L',:);
0518 ble = b(csense=='L');
0519 Aeq = A(csense=='E',:);
0520 beq = b(csense=='E');
0521
0522
0523 [neq nvar]=size(Aeq);
0524 nle=size(Ale,1);
0525 if isempty(EleNames)
0526 EleNames=arrayfun(EleNameFun,(1:nle),'UniformOutput', false);
0527 end
0528 if isempty(EqtNames)
0529 EqtNames=arrayfun(EqtNameFun,(1:neq),'UniformOutput', false);
0530 end
0531 if isempty(VarNames)
0532 VarNames=arrayfun(VarNameFun,(1:nvar),'UniformOutput', false);
0533 end
0534
0535
0536 [solution] = BuildMPS(Ale, ble, Aeq, beq, c, lb, ub, PbName,'MPSfilename',MPSfilename,'EleNames',EleNames,'EqtNames',EqtNames,'VarNames',VarNames);
0537
0538
0539 otherwise
0540 error(['Unknown solver: ' solver]);
0541
0542 end
0543 if ~strcmp(solver,'cplex_direct') && ~strcmp(solver,'mps')
0544
0545 t = etime(clock, t_start);
0546 if ~exist('basis','var'), basis=[]; end
0547 [solution.full,solution.obj,solution.rcost,solution.dual,solution.solver,solution.stat,solution.origStat,solution.time,solution.basis] = ...
0548 deal(x,f,w,y,solver,stat,origStat,t,basis);
0549 end
0550
0551
0552 function [x,f,y,w,stat,origStat] = solveGlpk(c,A,b,lb,ub,csense,osense,params)
0553
0554
0555
0556 [x,f,origStat,extra] = glpk(c,A,b,lb,ub,csense,[],osense,params);
0557 y = extra.lambda;
0558 w = extra.redcosts;
0559
0560 if (origStat == 180 || origStat == 5)
0561 stat = 1;
0562 elseif (origStat == 182 || origStat == 183 || origStat == 3 || origStat == 110)
0563 stat = 0;
0564 elseif (origStat == 184 || origStat == 6)
0565 stat = 2;
0566 else
0567 stat = -1;
0568 end