optimizeProb

 optimizeProb
   Optimize an LP or MILP formulated in cobra terms.

   prob    cobra style LP/MILP problem struct to be optimised
   params    solver specific parameters (optional)
   verbose if true MILP progress is shown (optional, default true)

   res        the output structure from the selected solver RAVENSOLVER
           (cobra style)

0001 function res = optimizeProb(prob,params,verbose)
0002 % optimizeProb
0003 %   Optimize an LP or MILP formulated in cobra terms.
0004 %
0005 %   prob    cobra style LP/MILP problem struct to be optimised
0006 %   params    solver specific parameters (optional)
0007 %   verbose if true MILP progress is shown (optional, default true)
0008 %
0009 %   res        the output structure from the selected solver RAVENSOLVER
0010 %           (cobra style)
0011 
0012 if nargin<2 || isempty(params)
0013     params=struct();
0014 end
0015 if nargin<3 || isempty(verbose)
0016     verbose = true;
0017 end
0018 
0019 %Set as global variable for speed improvement if optimizeProb is run many times
0020 global RAVENSOLVER;
0021 if isempty(RAVENSOLVER)
0022     if(~ispref('RAVEN','solver'))
0023         dispEM('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').');
0024     else
0025         RAVENSOLVER = getpref('RAVEN','solver');
0026     end
0027 end
0028 solver=RAVENSOLVER;
0029 if ~all(lower(prob.vartype) == 'c')
0030     milp=true;
0031     errorText = 'glpk is not suitable for solving MILPs, ';
0032     switch solver
0033         case 'glpk'
0034             error([errorText 'select a different solver with setRavenSolver().'])
0035         case 'cobra'
0036             if strcmp(CBT_MILP_SOLVER,'glpk')
0037                 error([errorText 'select a different solver with changeCobraSolver() or setRavenSolver().'])
0038             end
0039     end
0040 else
0041     milp=false;
0042 end
0043 
0044 %% Define default parameters, which will then be used to make solver-
0045 % specific solverparams structures
0046 defaultparams.feasTol        = 1e-9;
0047 defaultparams.optTol         = 1e-9;
0048 defaultparams.objTol         = 1e-6;
0049 defaultparams.timeLimit      = 1000;
0050 %defaultparams.iterationLimit = 1000;
0051 defaultparams.intTol         = 1e-12;
0052 defaultparams.relMipGapTol   = 1e-12;
0053 defaultparams.absMipGapTol   = 1e-12;
0054 if milp
0055     defaultparams.MIPGap     = 1e-12;
0056     defaultparams.Seed       = 1;
0057 end
0058 res.obj=[];
0059 switch solver
0060     %% Use whatever solver is set by COBRA Toolbox changeCobraSolver
0061     case 'cobra'
0062         if milp
0063             cparams=struct('timeLimit',1e9,'printLevel',0,'intTol',1e-6,'relMipGapTol',1e-9);
0064             cparams=structUpdate(cparams,params);
0065             res=solveCobraMILP(prob,cparams);
0066         else
0067             res=solveCobraLP(prob);
0068         end
0069         if isfield(res,{'dual','rcost'})
0070             res.dual=res.dual;
0071             res.rcost=res.rcost;
0072         end
0073 
0074         %% Use Gurobi in a MATLAB environment
0075     case 'gurobi'
0076         if milp
0077             if verbose
0078                 solverparams.OutputFlag = 1;
0079             else
0080                 solverparams.OutputFlag = 0;
0081             end
0082             solverparams.IntFeasTol = 10^-9; %min val for gurobi
0083             solverparams.MIPGap = defaultparams.MIPGap;
0084             solverparams.Seed = defaultparams.Seed;
0085         else
0086             solverparams.OutputFlag = 0;
0087         end
0088         solverparams.DisplayInterval= 1; % Level of verbosity
0089         solverparams.TimeLimit      = defaultparams.timeLimit;
0090         solverparams.FeasibilityTol = defaultparams.feasTol;
0091         solverparams.OptimalityTol  = defaultparams.optTol;
0092         solverparams.Presolve       = 2;
0093         if ~isempty(getCurrentTask) % If run in parallel, then one thread per gurobi
0094             solverparams.Threads=1;
0095         end
0096         solverparams = structUpdate(solverparams,params);
0097 
0098         % Restructering problem according to gurobi format
0099         if isfield(prob, 'csense')
0100             prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='});
0101             prob = rmfield(prob, {'csense'});
0102         end
0103         if isfield(prob, 'osense')
0104             osense = prob.osense;
0105             prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'});
0106             prob = rmfield(prob, {'osense'});
0107         end
0108         [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype);
0109         prob = rmfield(prob, {'c','b','vartype'});
0110 
0111         resG = gurobi(prob,solverparams);
0112 
0113         try
0114             % Name output fields the same as COBRA does
0115             res.full     = resG.x;
0116             res.obj      = resG.objval;
0117             res.origStat = resG.status;
0118             if isfield(resG,{'pi','rc'})
0119                 res.dual     = -resG.pi*osense;
0120                 res.rcost    = -resG.rc*osense;
0121             end
0122             if milp && strcmp(resG.status, 'TIME_LIMIT')
0123                 % If res has the objval field, it succeeded, regardless of
0124                 % time_limit status
0125                 resG.status = 'OPTIMAL';
0126             end
0127             switch resG.status
0128                 case 'OPTIMAL'
0129                     res.stat = 1;
0130                 case 'UNBOUNDED'
0131                     res.stat = 2;
0132                 otherwise
0133                     res.stat = 0;
0134             end
0135             if ~milp
0136                 res.vbasis = resG.vbasis;
0137                 res.cbasis = resG.cbasis;
0138             else
0139                 res.mipgap = resG.mipgap;
0140             end
0141         catch
0142             res.stat = 0;
0143             res.origStat = resG.status;  % useful information to have
0144         end
0145         %% Use GLPK using RAVEN-provided binary
0146     case 'glpk'
0147         solverparams.scale   = 1; % Auto scaling
0148         %solverparams.tmlim   = defaultparams.timeLimit;
0149         solverparams.tolbnd  = defaultparams.feasTol;
0150         solverparams.toldj   = defaultparams.optTol;
0151         solverparams.tolint  = defaultparams.intTol;
0152         solverparams.tolobj  = defaultparams.objTol;
0153         solverparams.msglev  = 0; % Level of verbosity
0154         solverparams = structUpdate(solverparams,params);
0155 
0156         prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'});
0157 
0158         if milp
0159             solverparams.tmlim   = solverparams.tmlim*10;
0160             solverparams.msglev  = 1; % Level of verbosity
0161             disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.')
0162         end
0163 
0164         % Ensure that RAVEN glpk binary is used, return to original
0165         % directory afterwards
0166         [ravenDir,currDir]=findRAVENroot();
0167         cd(fullfile(ravenDir,'software','GLPKmex'))
0168         [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams);
0169         cd(currDir)
0170 
0171         switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution
0172             case 5
0173                 res.stat = 1; % Optimal
0174             case 2
0175                 res.stat = 2; % Feasible, but not optimal
0176             otherwise
0177                 res.stat = 0;
0178         end
0179         res.origStat = errnum;
0180         res.full     = xopt;
0181         res.obj      = fmin;
0182         res.dual     = -extra.lambda*prob.osense;
0183         res.rcost    = -extra.redcosts*prob.osense;
0184         %% Use scip
0185     case {'soplex','scip'} % Old 'soplex' option also allowed
0186         [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype);
0187 
0188         %   [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts)
0189         %
0190         %   Input arguments*:
0191         %       H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU])
0192         %       f - linear objective vector
0193         %       A - linear constraint matrix (sparse)
0194         %       rl - linear constraint lhs
0195         %       ru - linear constraint rhs
0196         %       lb - decision variable lower bounds
0197         %       ub - decision variable upper bounds
0198         %       xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary)
0199         %       sos - SOS structure with fields type, index and weight (see below)
0200         %       qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below)
0201         %       nl - Nonlinear Objective and Constraints structure (see below)
0202         %       x0 - primal solution
0203         %       opts - solver options (see below)
0204         %
0205         %   Return arguments:
0206         %       x - solution vector
0207         %       fval - objective value at the solution
0208         %       exitflag - exit status (see below)
0209         %       stats - statistics structure
0210         %
0211         %   Option Fields (all optional, see also optiset for a list):
0212         %       solverOpts - specific SCIP options (list of pairs of parameter names and values)
0213         %       maxiter - maximum LP solver iterations
0214         %       maxnodes - maximum nodes to explore
0215         %       maxtime - maximum execution time [s]
0216         %       tolrfun - primal feasibility tolerance
0217         %       display - solver display level [0-5]
0218         %       probfile - write problem to given file
0219         %       presolvedfile - write presolved problem to file
0220         %
0221         %   Return Status:
0222         %       0 - Unknown
0223         %       1 - User Interrupted
0224         %       2 - Node Limit Reached
0225         %       3 - Total Node Limit Reached
0226         %       4 - Stall Node Limit Reached
0227         %       5 - Time Limit Reached
0228         %       6 - Memory Limit Reached
0229         %       7 - Gap Limit Reached
0230         %       8 - Solution Limit Reached
0231         %       9 - Solution Improvement Limit Reached
0232         %      10 - Restart Limit Reached
0233         %      11 - Problem Solved to Optimality
0234         %      12 - Problem is Infeasible
0235         %      13 - Problem is Unbounded
0236         %      14 - Problem is Either Infeasible or Unbounded
0237         
0238         res.origStat = exitflag;
0239         res.full = xopt;
0240         res.obj  = fval;
0241 
0242         switch exitflag
0243             case 11
0244                 res.stat = 1;
0245             case [5, 6, 7, 8, 9, 10, 13]
0246                 res.stat = 2;
0247             otherwise
0248                 res.stat = 0;
0249         end
0250     otherwise
0251         error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').');
0252 end
0253 if res.stat>0
0254     res.full=res.full(1:size(prob.a,2));
0255 end
0256 end
0257 
0258 function s_merged=structUpdate(s_old,s_new)
0259 %Remove overlapping fields from first struct;
0260 %Obtain all unique names of remaining fields;
0261 %Merge both structs
0262 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new)));
0263 names = [fieldnames(s_merged); fieldnames(s_new)];
0264 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1);
0265 end
0266 
0267 function paramlist = renameparams(paramlist,old,new)
0268 if ~iscell(paramlist)
0269     wasNoCell = true;
0270     paramlist={paramlist};
0271 else
0272     wasNoCell = false;
0273 end
0274 for i=1:numel(old)
0275     paramlist = regexprep(paramlist,old{i},new{i});
0276 end
0277 if wasNoCell
0278     paramlist=paramlist{1};
0279 end
0280 end