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