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`	Optimize 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

either 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 minimized (default model.rxns).

params : struct

obsolete option.

scores : double

vector of weights for the reactions. Negative scores should not have flux. Positive scores are not possible in this implementation, and they 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 minimized costs 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, at the cost of a larger constraint matrix. Fewer binaries makes this the faster formulation on large problems, but it needs a solver that handles MILPs well and no LP estimate of the flux magnitudes is made, so an arbitrary large flux is assumed.

verbose : logical

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

resolveTies : logical

if true, pin the MILP's degenerate optimum to a canonical answer instead of relying on Seed alone: hold the primary (scores) objective at its optimum, then minimize the count of "on" reactions, then, among the sparsest, minimize their summed reaction-id rank (lowest ids preferred). Two extra MILP solves; best-effort, so a phase that does not converge within TimeLimit is dropped with a warning rather than changing the result. Only implemented for formulation='reversible' -- errors if combined with 'irrev'. Matches raven-toolbox's resolve_ties (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, but it may use more fluxes than the minimum

Examples

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

Notes

Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly cause problems if the fluxes in the model are larger than that.

optimizeProb#

function res = optimizeProb(prob,varargin)

Optimize an LP or MILP formulated in COBRA terms.

Parameters

prob : struct: COBRA-style LP/MILP problem struct to be optimised.

Name-Value Arguments

params : struct: solver-specific parameters. In addition to solver parameters, the field "maxRatio" can be set to a number > 1 to improve numerical conditioning: before solving (LP only) any column whose coefficients span more than maxRatio is split via auxiliary metabolites/variables, preserving the feasible region. See splitProbForConditioning. The field is consumed here and not forwarded to the solver.
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. Whatever solution is 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). A RAVEN:glpkMILP warning then recommends Gurobi, which is much faster and more reliable on larger MILPs (such as ftINIT or getMinNrFluxes on genome-scale models).

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 a RAVEN:glpkMILP warning recommends Gurobi, which is much faster and more reliable on 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. MILPs are refused 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 optimization should be performed in order to get rid of loops in the flux distribution (default 0):

0 : no such optimization is performed
1 : the sum of abs(fluxes) is minimized. This is the fastest way of getting rid of loops
2 : option is deprecated
3 : the number of fluxes is minimized. This can result in the flux distributions that are the easiest to interpret. Note that this optimization can be very slow

params : struct

solver parameters, forwarded to optimizeProb. Mostly obsolete, but the field "maxRatio" (a number > 1) can be set to improve numerical conditioning of ill-scaled models: any reaction whose stoichiometric coefficients span more than maxRatio is transiently split via 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 by default).

hsSol : struct

hot-start solution for the LP solver. This can significantly speed up the process if many similar optimization problems are solved iteratively. Only used if minFlux is 0 or 1.

Returns

solution : struct

LP solution with fields:

f : objective value
x : primal (flux distribution)
stat : exit flag (1 = optimization terminated successfully, 0 = solution is feasible but not necessarily optimal, -1 = no feasible solution found, -2 = solution found but flux minimization failed)
msg : string describing the status of the optimization
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). Only used if minFlux is 0 or 1.

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 (e.g. a biomass reaction with both major substrates and trace cofactors) cannot be removed by the diagonal row/column scaling that LP solvers apply, and may lead to unreliable solutions or spurious infeasibility.

Splitting is done via auxiliary metabolites and variables that preserve the feasible region of the original problem exactly, so the values of all original variables are unchanged.

The function is meant to be called inside optimizeProb, on the COBRA-style problem struct, immediately before the problem is passed to the solver. Auxiliary rows/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 (possibly) 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 auxiliary metabolite/conversion pairs added
nOrigRows : number of constraint rows before augmentation
nOrigCols : number of variables (columns) before augmentation

Examples

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