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 run in parallel, then one thread per gurobi. getCurrentTask
0094         % requires the Parallel Computing Toolbox.
0095         inParallel = false;
0096         if license('test','Distrib_Computing_Toolbox')
0097             try
0098                 inParallel = ~isempty(getCurrentTask);
0099             catch
0100             end
0101         end
0102         if inParallel
0103             solverparams.Threads=1;
0104         end
0105         solverparams = structUpdate(solverparams,params);
0106 
0107         % Restructering problem according to gurobi format
0108         if isfield(prob, 'csense')
0109             prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='});
0110             prob = rmfield(prob, {'csense'});
0111         end
0112         if isfield(prob, 'osense')
0113             osense = prob.osense;
0114             prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'});
0115             prob = rmfield(prob, {'osense'});
0116         end
0117         [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype);
0118         prob = rmfield(prob, {'c','b','vartype'});
0119 
0120         resG = gurobi(prob,solverparams);
0121 
0122         try
0123             % Name output fields the same as COBRA does
0124             res.full     = resG.x;
0125             res.obj      = resG.objval;
0126             res.origStat = resG.status;
0127             if isfield(resG,{'pi','rc'})
0128                 res.dual     = -resG.pi*osense;
0129                 res.rcost    = -resG.rc*osense;
0130             end
0131             if milp && strcmp(resG.status, 'TIME_LIMIT')
0132                 % If res has the objval field, it succeeded, regardless of
0133                 % time_limit status
0134                 resG.status = 'OPTIMAL';
0135             end
0136             switch resG.status
0137                 case 'OPTIMAL'
0138                     res.stat = 1;
0139                 case 'UNBOUNDED'
0140                     res.stat = 2;
0141                 otherwise
0142                     res.stat = 0;
0143             end
0144             if ~milp
0145                 res.vbasis = resG.vbasis;
0146                 res.cbasis = resG.cbasis;
0147             else
0148                 res.mipgap = resG.mipgap;
0149             end
0150         catch
0151             res.stat = 0;
0152             res.origStat = resG.status;  % useful information to have
0153         end
0154         %% Use GLPK using RAVEN-provided binary
0155     case 'glpk'
0156         solverparams.scale   = 1; % Auto scaling
0157         %solverparams.tmlim   = defaultparams.timeLimit;
0158         solverparams.tolbnd  = defaultparams.feasTol;
0159         solverparams.toldj   = defaultparams.optTol;
0160         solverparams.tolint  = defaultparams.intTol;
0161         solverparams.tolobj  = defaultparams.objTol;
0162         solverparams.msglev  = 0; % Level of verbosity
0163         solverparams = structUpdate(solverparams,params);
0164 
0165         prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'});
0166 
0167         if milp
0168             solverparams.tmlim   = solverparams.tmlim*10;
0169             solverparams.msglev  = 1; % Level of verbosity
0170             disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.')
0171         end
0172 
0173         % Ensure that RAVEN glpk binary is used, return to original
0174         % directory afterwards
0175         [ravenDir,currDir]=findRAVENroot();
0176         cd(fullfile(ravenDir,'software','GLPKmex'))
0177         [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams);
0178         cd(currDir)
0179 
0180         switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution
0181             case 5
0182                 res.stat = 1; % Optimal
0183             case 2
0184                 res.stat = 2; % Feasible, but not optimal
0185             otherwise
0186                 res.stat = 0;
0187         end
0188         res.origStat = errnum;
0189         res.full     = xopt;
0190         res.obj      = fmin;
0191         res.dual     = -extra.lambda*prob.osense;
0192         res.rcost    = -extra.redcosts*prob.osense;
0193         %% Use scip
0194     case {'soplex','scip'} % Old 'soplex' option also allowed
0195         [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype);
0196 
0197         %   [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts)
0198         %
0199         %   Input arguments*:
0200         %       H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU])
0201         %       f - linear objective vector
0202         %       A - linear constraint matrix (sparse)
0203         %       rl - linear constraint lhs
0204         %       ru - linear constraint rhs
0205         %       lb - decision variable lower bounds
0206         %       ub - decision variable upper bounds
0207         %       xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary)
0208         %       sos - SOS structure with fields type, index and weight (see below)
0209         %       qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below)
0210         %       nl - Nonlinear Objective and Constraints structure (see below)
0211         %       x0 - primal solution
0212         %       opts - solver options (see below)
0213         %
0214         %   Return arguments:
0215         %       x - solution vector
0216         %       fval - objective value at the solution
0217         %       exitflag - exit status (see below)
0218         %       stats - statistics structure
0219         %
0220         %   Option Fields (all optional, see also optiset for a list):
0221         %       solverOpts - specific SCIP options (list of pairs of parameter names and values)
0222         %       maxiter - maximum LP solver iterations
0223         %       maxnodes - maximum nodes to explore
0224         %       maxtime - maximum execution time [s]
0225         %       tolrfun - primal feasibility tolerance
0226         %       display - solver display level [0-5]
0227         %       probfile - write problem to given file
0228         %       presolvedfile - write presolved problem to file
0229         %
0230         %   Return Status:
0231         %       0 - Unknown
0232         %       1 - User Interrupted
0233         %       2 - Node Limit Reached
0234         %       3 - Total Node Limit Reached
0235         %       4 - Stall Node Limit Reached
0236         %       5 - Time Limit Reached
0237         %       6 - Memory Limit Reached
0238         %       7 - Gap Limit Reached
0239         %       8 - Solution Limit Reached
0240         %       9 - Solution Improvement Limit Reached
0241         %      10 - Restart Limit Reached
0242         %      11 - Problem Solved to Optimality
0243         %      12 - Problem is Infeasible
0244         %      13 - Problem is Unbounded
0245         %      14 - Problem is Either Infeasible or Unbounded
0246         
0247         res.origStat = exitflag;
0248         res.full = xopt;
0249         res.obj  = fval;
0250 
0251         switch exitflag
0252             case 11
0253                 res.stat = 1;
0254             case {1, 5, 6, 7, 8, 9, 10, 13} % If user-interrupted, not sure if optimal
0255                 res.stat = 2;
0256             otherwise
0257                 res.stat = 0;
0258         end
0259     otherwise
0260         error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').');
0261 end
0262 if res.stat>0
0263     res.full=res.full(1:size(prob.a,2));
0264 end
0265 end
0266 
0267 function s_merged=structUpdate(s_old,s_new)
0268 %Remove overlapping fields from first struct;
0269 %Obtain all unique names of remaining fields;
0270 %Merge both structs
0271 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new)));
0272 names = [fieldnames(s_merged); fieldnames(s_new)];
0273 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1);
0274 end
0275 
0276 function paramlist = renameparams(paramlist,old,new)
0277 if ~iscell(paramlist)
0278     wasNoCell = true;
0279     paramlist={paramlist};
0280 else
0281     wasNoCell = false;
0282 end
0283 for i=1:numel(old)
0284     paramlist = regexprep(paramlist,old{i},new{i});
0285 end
0286 if wasNoCell
0287     paramlist=paramlist{1};
0288 end
0289 end