Solvers#

MATLAB functions in RAVEN/src/solver of the RAVEN toolbox.

Functions#

Function

Summary

checkSolution

Check if a solver result is feasible and optimal.

getMinNrFluxes

Find the minimal set of fluxes that satisfy the model.

optimizeProb

Optimise an LP or MILP formulated in COBRA terms.

setRavenSolver

Set the solver for RAVEN.

solveLP

Solve a linear programming problem.

splitProbForConditioning

Improve LP conditioning by splitting wide columns.

Reference#

checkSolution#

function [isFeasible, isOptimal]=checkSolution(res)

Check if a solver result is feasible and optimal.

Parameters

res : struct

the output structure from optimizeProb.

Returns

isFeasible : logical

true if the solution is feasible.

isOptimal : logical

true if the solution is optimal.

Examples

[isFeasible, isOptimal] = checkSolution(res);

getMinNrFluxes#

function [x,I,exitFlag]=getMinNrFluxes(model, varargin)

Find the minimal set of fluxes that satisfy the model.

Uses mixed integer linear programming to find the minimal set of fluxes that satisfy the model.

Parameters

model : struct

a model structure.

Name-Value Arguments

toMinimize : cell or logical or double

a cell array of reaction IDs, a logical vector with the same number of elements as reactions in the model, or a vector of indexes for the reactions that should be minimised (default model.rxns).

params : struct

solver parameters passed to optimizeProb. With formulation 'reversible', intTol (10^-9), TimeLimit (300) and Seed (26) are set when absent (default empty struct).

scores : double

vector of weights for the reactions. A reaction with a negative score adds the negated score to the minimised objective when it carries flux. Scores of 0 or above are not supported and are changed to max(scores(scores<0)), or to 0 if scores has no negative entry at all. Must have the same dimension as toMinimize (find(toMinimize) if it is a logical vector) (default -1 for all reactions).

formulation : char

which MILP formulation to build (default 'irrev'):

  • 'irrev' : convert the model to irreversible format and use one binary per irreversible reaction, so a reversible reaction that is minimised needs two binaries. The magnitude of the fluxes is estimated with an LP first, which also provides a start solution.

  • 'reversible' : keep the model reversible and use one binary per reaction, with a larger constraint matrix. The formulation has fewer binaries but needs a solver that handles MILPs well. No LP estimate of the flux magnitudes is made, and a fixed flux of 1000 is assumed.

verbose : logical

if true, the MILP progression will be shown (default false).

resolveTies : logical

if true, select one optimum among the degenerate optima of the MILP instead of relying on Seed: the primary (scores) objective is held at its optimum, the count of "on" reactions is minimised, and among the sparsest solutions the summed reaction-id rank is minimised (lowest ids preferred). This takes two extra MILP solves. A phase that is not proven optimal within TimeLimit gives a warning, and the solution of the last phase is used when it is feasible. Only implemented for formulation='reversible'; combined with 'irrev' it gives an error. Matches resolve_ties in raven-toolbox (raven-gecko-parity#104) (default false).

Returns

x : double

the corresponding fluxes for the full model.

I : double

the indexes of the reactions in toMinimize that were used in the solution.

exitFlag : double

exit status:

  • 1 : optimal solution found

  • -1 : no solution found, either because the problem is infeasible or because the solver reached its time limit before finding one

  • -2 : a solution was found but is not proven optimal, because the solver stopped before reaching optimality. It is returned and may use more fluxes than the minimum

Examples

[x, I, exitFlag] = getMinNrFluxes(model, toMinimize, params, scores);

Notes

The binary of a reaction is linked to its flux through a large flux value. With formulation 'irrev' the value is the larger of 1000 and 5 times the largest flux of an LP solution; with 'reversible' it is 1000 (mmol/gDW/h). A flux above this value is cut off by the binary constraint of the reaction, so the MILP can become infeasible.

optimizeProb#

function res = optimizeProb(prob,varargin)

Optimise an LP or MILP formulated in COBRA terms.

Parameters

prob : struct

COBRA-style problem struct of an LP or MILP to be optimised.

Name-Value Arguments

params : struct

solver-specific parameters. In addition to solver parameters, the field "maxRatio" can be set to a finite number > 1 for numerical conditioning: before solving an LP, any column whose coefficients span more than maxRatio is split with auxiliary metabolites and variables, which preserves the feasible region. See splitProbForConditioning. The field is removed before the parameters are passed to the solver, and is ignored for MILPs.

verbose : logical

if true, MILP progress is shown (default true).

Returns

res : struct

the output structure from the selected solver RAVENSOLVER (COBRA style). In addition to the fields set by the solver, the logical field "hitTimeLimit" is true if the solver stopped because it reached its time limit. The solution in res is then either absent or not proven optimal.

Notes

MILPs can be solved with GLPK, either RAVEN's bundled GLPK ('glpk') or the COBRA Toolbox's GLPK interface ('cobra' with GLPK as its MILP solver). The RAVEN:glpkMILP warning then states that GLPK works for small problems and that Gurobi is much faster and more reliable for larger ones (such as ftINIT or getMinNrFluxes on genome-scale models).

See also

splitProbForConditioning

setRavenSolver#

function setRavenSolver(solver)

Set the solver for RAVEN.

Sets the solver for RAVEN and optionally saves a default value to MATLAB startup.

Parameters

solver : char

string with solver name. The following options are available:

  • glpk : uses RAVEN-provided binaries. Solves LPs and MILPs; for MILPs the RAVEN:glpkMILP warning states that Gurobi is much faster and more reliable for larger problems.

  • gurobi : requires a working Gurobi installation.

  • scip : requires a SCIP MEX interface. RAVEN provides a precompiled binary for Windows only.

  • cobra : requires a working COBRA Toolbox installation. RAVEN parses the problem to COBRA and uses the solver that is set by changeCobraSolver. For MILPs, the RAVEN:glpkMILP warning is issued when COBRA's MILP solver is GLPK.

Examples

setRavenSolver(solver);

solveLP#

function [solution, hsSolOut]=solveLP(model,varargin)

Solve a linear programming problem.

Parameters

model : struct

a model structure.

Name-Value Arguments

minFlux : double

determines if a second optimisation is performed to remove loops from the flux distribution (default 0):

  • 0 : no second optimisation is performed

  • 1 : the sum of abs(fluxes) is minimised with an LP, which does not need the MILP that the minimisation of the number of fluxes requires

  • 2 : deprecated; gives an error

  • 3 : the number of fluxes is minimised with getMinNrFluxes, which solves a MILP and can take long. The result has the fewest active reactions, which can make the flux distribution easier to interpret

params : struct

solver parameters, forwarded to optimizeProb. The field "maxRatio" (a number > 1) can be set for numerical conditioning of ill-scaled models: any reaction whose stoichiometric coefficients span more than maxRatio is split with auxiliary metabolites before solving, which preserves the feasible region (and thus the flux distribution of the original reactions). A common value is 1e6. No splitting is done when maxRatio is absent (default struct with relGap=0.4).

hsSol : struct

hot-start solution for the LP solver, as hsSolOut from a previous call. It speeds up the solution when many similar optimisation problems are solved iteratively. A hsSol with stat<1 is ignored.

Returns

solution : struct

LP solution with fields:

  • f : objective value

  • x : primal (flux distribution)

  • stat : exit flag (1 = optimisation terminated successfully, 0 = solution is feasible but not necessarily optimal, -1 = no feasible solution found, -2 = solution found but flux minimisation failed)

  • msg : string describing the status of the optimisation

  • sPrice : shadow price (only reported if minFlux was 0)

  • rCost : reduced cost (only reported if minFlux was 0)

hsSolOut : struct

solution to be used as hot-start solution (see the input parameters). Set only if minFlux is 0 or 1, and only when the solver returns a basis.

Examples

[solution, hsSolOut] = solveLP(model, minFlux, params, hsSol);

splitProbForConditioning#

function [prob, condInfo] = splitProbForConditioning(prob, maxRatio)

Improve LP conditioning by splitting wide columns.

Reformulates an LP problem to improve its numerical conditioning by splitting columns whose coefficients span more than a given ratio. In a metabolic model a column corresponds to a reaction, and a wide spread of stoichiometric coefficients within one reaction (for example a biomass reaction with both large and small coefficients) cannot be removed by the diagonal row and column scaling that LP solvers apply, and may lead to unreliable solutions or spurious infeasibility.

Splitting is done with auxiliary metabolites and variables that preserve the feasible region of the original problem exactly, so the values of all original variables are unchanged. A column is split when its ratio exceeds maxRatio by more than 1 %. Columns with a single nonzero coefficient are never split.

optimizeProb calls this function on the COBRA-style problem struct of an LP, immediately before the problem is passed to the solver. Auxiliary rows and columns are appended after the original ones, so the augmentation can be stripped from the solver output afterwards (see condInfo).

Parameters

prob : struct

COBRA-style LP problem struct, with at least the fields A, b, c, lb, ub, csense and vartype.

maxRatio : double

maximum allowed ratio between the largest and smallest absolute coefficient within any single column.

Returns

prob : struct

the augmented problem struct. Identical to the input if no column needed splitting.

condInfo : struct

struct describing the augmentation, with fields:

  • applied : true if any column was split

  • nSplits : number of pairs of an auxiliary metabolite and a conversion variable added

  • nOrigRows : number of constraint rows before augmentation

  • nOrigCols : number of variables (columns) before augmentation

Examples

[prob, condInfo] = splitProbForConditioning(prob, maxRatio);

See also

optimizeProb