MathProgBase.jl

MathProgBase.jl provides high-level one-shot functions for linear and mixed-integer programming, as well as a solver-independent low-level interface for implementing advanced techniques that require efficiently solving a sequence of linear programming problems.

To use MathProgBase, an external solver must be installed. See choosing solvers.

Contents

Linear Programming

linprog(c, A, sense, b, l, u, solver)

Solves the linear programming problem:

\[\begin{split}\min_{x}\, &c^Tx\\ s.t. &a_i^Tx \text{ sense}_i \, b_i \forall\,\, i\\ &l \leq x \leq u\\\end{split}\]

where:

  • c is the objective vector, always in the sense of minimization
  • A is the constraint matrix, with rows \(a_i\) (viewed as column-oriented vectors)
  • sense is a vector of constraint sense characters '<', '=', and '>'
  • b is the right-hand side vector
  • l is the vector of lower bounds on the variables
  • u is the vector of upper bounds on the variables, and
  • solver is an optional parameter specifying the desired solver, see choosing solvers. If this parameter is not provided, the default solver is used.

A scalar is accepted for the b, sense, l, and u arguments, in which case its value is replicated. The values -Inf and Inf are interpreted to mean that there is no corresponding lower or upper bound.

Note

Linear programming solvers extensively exploit the sparsity of the constraint matrix A. While both dense and sparse matrices are accepted, for large-scale problems sparse matrices should be provided if permitted by the problem structure.

A shortened version is defined as:

linprog(c, A, sense, b, solver) = linprog(c, A, sense, b, 0, Inf, solver)

The linprog function returns an instance of the type:

type LinprogSolution
    status
    objval
    sol
    attrs
end

where status is a termination status symbol, one of :Optimal, :Infeasible, :Unbounded, :UserLimit (iteration limit or timeout), :Error (and maybe others).

If status is :Optimal, the other members have the following values:

  • objval – optimal objective value
  • sol – primal solution vector
  • attrs – a dictionary that may contain other relevant attributes such as:
    • redcost – dual multipliers for active variable bounds (zero if inactive)
    • lambda – dual multipliers for active linear constraints (equalities are always active)

If status is :Infeasible, the attrs member will contain an infeasibilityray if available; similarly for :Unbounded problems, attrs will contain an unboundedray if available.

For example, we can solve the two-dimensional problem (see test/linprog.jl):

\[\begin{split}\min_{x,y}\, &-x\\ s.t. &2x + y \leq 1.5\\ & x \geq 0, y \geq 0\end{split}\]

by:

using MathProgBase

sol = linprog([-1,0],[2 1],'<',1.5)
if sol.status == :Optimal
    println("Optimal objective value is $(sol.objval)")
    println("Optimal solution vector is: [$(sol.sol[1]), $(sol.sol[2])]")
else
    println("Error: solution status $(sol.status)")
end
linprog(c, A, lb, ub, l, u, solver)

This variant allows one to specify two-sided linear constraints (also known as range constraints) to solve the linear programming problem:

\[\begin{split}\min_{x}\, &c^Tx\\ s.t. &lb \leq Ax \leq ub\\ &l \leq x \leq u\\\end{split}\]

where:

  • c is the objective vector, always in the sense of minimization
  • A is the constraint matrix
  • lb is the vector of row lower bounds
  • ub is the vector of row upper bounds
  • l is the vector of lower bounds on the variables
  • u is the vector of upper bounds on the variables, and
  • solver is an optional parameter specifying the desired solver, see choosing solvers. If this parameter is not provided, the default solver is used.

A scalar is accepted for the l, u, lb, and ub arguments, in which case its value is replicated. The values -Inf and Inf are interpreted to mean that there is no corresponding lower or upper bound. Equality constraints are specified by setting the row lower and upper bounds to the same value.

A shortened version is defined as:

linprog(c, A, lb, ub, solver) = linprog(c, A, lb, ub, 0, Inf, solver)

Mixed-integer Programming

mixintprog(c, A, sense, b, vartypes, lb, ub, solver)

Solves the same optimization problem as linprog above, except variables are additionally constrained to take only integer values if the corresponding entry in the varypes vector is the symbol :Int. Continuous variables are indicated by the value :Cont, binary variables should be specified by :Bin, semicontinuous by :SemiCont, and semi-integer by :SemiInt.

A scalar is accepted for the sense, b, vartypes, lb, and ub arguments, in which case its value is replicated. The values -Inf and Inf are interpreted to mean that there is no corresponding lower or upper bound.

The mixintprog function returns an instance of the type:

type MixintprogSolution
    status
    objval
    sol
    attrs
end

where status takes the same values as with linprog.

If status does not indicate error or infeasiblity, the other members have the following values:

  • objval – optimal objective value
  • sol – primal solution vector
  • attrs – a dictionary that may contain other relevant attributes such as:
    • objbound – Best known lower bound on the objective value

Analogous shortened and range-constraint versions are available as well.

We can solve a binary knapsack problem

\[\begin{split}max\, &5x_1 + 3x_2 + 2x_3 + 7x_4 + 4x_5\\ s.t. &2x_1 + 8x_2 + 4x_3 + 2x_4 + 5x_5 \leq 10\\ & (x_1, x_2, x_3, x_4, x_5) \in \{0,1\}^5\end{split}\]

with the code:

mixintprog(-[5.,3.,2.,7.,4.],[2. 8. 4. 2. 5.],'<',10,:Int,0,1)

Quadratic Programming

quadprog(c, Q, A, sense, b, l, u, solver)

Solves the quadratic programming problem:

\[\begin{split}\min_{x}\, &\frac{1}{2}x^TQx + c^Tx\\ s.t. &a_i^Tx \text{ sense}_i \, b_i \forall\,\, i\\ &l \leq x \leq u\\\end{split}\]

where:

  • c is the objective vector, always in the sense of minimization
  • Q is the Hessian matrix of the objective
  • A is the constraint matrix, with rows \(a_i\) (viewed as column-oriented vectors)
  • sense is a vector of constraint sense characters '<', '=', and '>'
  • b is the right-hand side vector
  • l is the vector of lower bounds on the variables
  • u is the vector of upper bounds on the variables, and
  • solver is an optional parameter specifying the desired solver, see choosing solvers. If this parameter is not provided, the default solver is used.

A scalar is accepted for the b, sense, l, and u arguments, in which case its value is replicated. The values -Inf and Inf are interpreted to mean that there is no corresponding lower or upper bound.

Note

Quadratic programming solvers extensively exploit the sparsity of the Hessian matrix Q and the constraint matrix A. While both dense and sparse matrices are accepted, for large-scale problems sparse matrices should be provided if permitted by the problem structure.

The quadprog function returns an instance of the type:

type QuadprogSolution
    status
    objval
    sol
    attrs
end

where status is a termination status symbol, one of :Optimal, :Infeasible, :Unbounded, :UserLimit (iteration limit or timeout), :Error (and maybe others).

If status is :Optimal, the other members have the following values:

  • objval – optimal objective value
  • sol – primal solution vector
  • attrs – a dictionary that may contain other relevant attributes (not currently used).

Analogous shortened and range-constraint versions are available as well.

We can solve the three-dimensional QP (see test/quadprog.jl):

\[\begin{split}\min_{x,y,z}\, &x^2+y^2+z^2+xy+yz\\ s.t. &x + 2y + 3z \geq 4\\ &x + y \geq 1\end{split}\]

by:

using MathProgBase

sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf)
if sol.status == :Optimal
    println("Optimal objective value is $(sol.objval)")
    println("Optimal solution vector is: [$(sol.sol[1]), $(sol.sol[2]), $(sol.sol[3])]")
else
    println("Error: solution status $(sol.status)")
end

Low-level interface

The linprog and mixintprog functions are written on top of a solver-independent low-level interface called MathProgSolverInterface, which individual solvers implement. The concept is similar to that of OSI, a C++ library which provides a generic virtual base class for interacting with linear programming solvers. Julia, however, does not quite have virtual classes or interfaces. Instead, multiple dispatch is used with abstract types. The API is designed to support problem modification as needed to solve a sequence of linear programming problems efficiently; linear programming solvers are expected to hot-start the solution process after modifications such as additional constraints or variables. For mixed-integer programming, hot-starting is usually impractical.

The MathProgSolverInterface exports two abstract types: AbstractMathProgModel, which represents an instance of an optimization problem, and AbstractMathProgSolver, which represents a solver (with particular solution options), from which an AbstractMathProgModel is generated.

model(s::AbstractMathProgSolver)

Returns an instance of an AbstractMathProgModel using the given solver.

loadproblem!(m::AbstractMathProgModel, filename::String)

Loads problem data from the given file. Supported file types are solver-dependent.

loadproblem!(m::AbstractMathProgModel, A, l, u, c, lb, ub, sense)

Loads the provided problem data to set up the linear programming problem:

\[\begin{split}\min_{x}\, &c^Tx\\ s.t. &lb \leq Ax \leq ub\\ &l \leq x \leq u\\\end{split}\]

sense specifies the direction of the optimization problem, and must be either :Min or :Max.

Both sparse and dense matrices are accepted for A. Inf and -Inf indicate that there is no corresponding upper or lower bound. Equal lower and upper bounds are used to indicate equality constraints.

writeproblem(m::AbstractMathProgModel, filename::String)

Writes the current problem data to the given file. Supported file types are solver-dependent.

getvarLB(m::AbstractMathProgModel)

Returns a vector containing the lower bounds \(l\) on the variables.

setvarLB!(m::AbstractMathProgModel, l)

Sets the lower bounds on the variables.

getvarUB(m::AbstractMathProgModel)

Returns a vector containing the upper bounds \(u\) on the variables.

setvarUB!(m::AbstractMathProgModel, u)

Sets the upper bounds on the variables.

getconstrLB(m::AbstractMathProgModel)

Returns a vector containing the lower bounds \(lb\) on the linear constraints.

setconstrLB!(m::AbstractMathProgModel, lb)

Sets the lower bounds on the linear constraints.

getconstrUB(m::AbstractMathProgModel)

Returns a vector containing the upper bounds \(ub\) on the linear constraints.

setconstrUB!(m::AbstractMathProgModel, ub)

Sets the upper bounds on the linear constraints.

getobj(m::AbstractMathProgModel)

Returns a vector containing the linear objective coefficients \(c\).

setobj!(m::AbstractMathProgModel, c)

Sets the linear objective coefficients.

getconstrmatrix(m::AbstractMathProgModel)

Returns the full linear constraint matrix \(A\), typically as a SparseMatrixCSC.

addvar!(m::AbstractMathProgModel, constridx, constrcoef, l, u, objcoef)

Adds a new variable to the model, with lower bound l (-Inf if none), upper bound u (Inf if none), and objective coefficient objcoef. Constraint coefficients for this new variable are specified in a sparse format: the constrcoef vector contains the nonzero coefficients, and the constridx vector contains the indices of the corresponding linear constraints.

addvar!(m::AbstractMathProgModel, l, u, objcoef)

Adds a new variable to the model, with lower bound l (-Inf if none), upper bound u (Inf if none), and objective coefficient objcoef. This is equivalent to calling the above method with empty arrays for the constraint coefficients.

addconstr!(m::AbstractMathProgModel, varidx, coef, lb, ub)

Adds a new linear constraint to the model, with lower bound lb (-Inf if none) and upper bound ub (Inf if none). Coefficients for this new constraint are specified in a sparse format: the coef vector contains the nonzero coefficients, and the varidx vector contains the indices of the corresponding variables.

updatemodel!(m::AbstractMathProgModel)

Commits recent changes to the model. Only required by some solvers (e.g. Gurobi).

setsense!(m::AbstractMathProgModel, sense)

Sets the optimization sense of the model. Accepted values are :Min and :Max.

getsense(m::AbstractMathProgModel)

Returns the optimization sense of the model.

numvar(m::AbstractMathProgModel)

Returns the number of variables in the model.

numconstr(m::AbstractMathProgModel)

Returns the total number of constraints in the model.

numlinconstr(m::AbstractMathProgModel)

Returns the number of linear constraints in the model.

optimize!(m::AbstractMathProgModel)

Solves the optimization problem.

status(m::AbstractMathProgModel)

Returns the termination status after solving. Possible values include :Optimal, :Infeasible, :Unbounded, :UserLimit (iteration limit or timeout), and :Error. Solvers may return other statuses, for example, when presolve indicates that the model is either infeasible or unbounded, but did not determine which.

getobjval(m::AbstractMathProgModel)

Returns the objective value of the solution found by the solver. In particular, this may be the objective value for the best feasible solution if optimality is not attained or proven.

getobjbound(m::AbstractMathProgModel)

Returns the best known bound on the optimal objective value. This is used, for example, when a branch-and-bound method is stopped before finishing.

getsolution(m::AbstractMathProgModel)

Returns the solution vector found by the solver.

getconstrsolution(m::AbstractMathProgModel)

Returns a vector containing the values of the linear constraints at the solution. This is the vector \(Ax\).

getreducedcosts(m::AbstractMathProgModel)

Returns the dual solution vector corresponding to the variable bounds, known as the reduced costs. Not available when integer variables are present.

getconstrduals(m::AbstractMathProgModel)

Returns the dual solution vector corresponding to the linear constraints. Not available when integer variables are present.

getinfeasibilityray(m::AbstractMathProgModel)

Returns a “Farkas” proof of infeasibility, i.e., an unbounded ray of the dual, for the linear constraints. Note that for some solvers, one must specify additional options for this ray to be computed.

getbasis(m::AbstractMathProgModel)

Returns the basis set for the optimal solution in the form (cbasis,rbasis), where both return values are vectors of symbols. The vector cbasis indexes the columns of the constraint matrix, while rbasis indexes the rows (values indicate whether the constraint is active at a lower/upper bound). The entries take value :Basic if the element is basic, :NonbasicAtLower if it is nonbasic at a lower bound, and :NonbasicAtUpper if it is nonbasic at upper bound. Other values may appear, taking solver-specific values. Note that this function may not work if the optimization algorithm is not able to provide basis information.

getunboundedray(m::AbstractMathProgModel)

Returns an unbounded ray of the problem, i.e., an objective-improving direction in which one may travel an infinite distance without violating any constraints. Note that for some solvers, one must specify additional options for this ray to be computed.

getsolvetime(m::AbstractMathProgModel)

Returns the total elapsed solution time as reported by the solver.

getsimplexiter(m::AbstractMathProgModel)

Returns the cumulative number of simplex iterations during the optimization process. In particular, for a MIP it returns the total simplex iterations for all nodes.

getbarrieriter(m::AbstractMathProgModel)

Returns the cumulative number of barrier iterations during the optimization process.

getobjgap(m::AbstractMathProgModel)

Returns the final relative optimality gap as optimization terminated. That is, it returns \(\frac{|b-f|}{|f|}\), where \(b\) is the best bound and \(f\) is the best feasible objective value.

getrawsolver(m::AbstractMathProgModel)

Returns an object that may be used to access a solver-specific API for this model.

setwarmstart!(m::AbstractMathProgModel, v)

Provide an initial solution v to the solver, as supported. To leave values undefined, set them to NaN. MIP solvers should ignore provided solutions that are infeasible or cannot be completed to a feasible solution. Nonlinear solvers may use provided solutions as starting points even if infeasible.

Integer Programming

getnodecount(m::AbstractMathProgModel)

Returns the total number of branch-and-bound nodes explored during the MIP optimization process.

setvartype!(m::AbstractMathProgModel, v::Vector{Symbol})

Sets the types of the variables to those indicated by the vector v. Valid types are :Int for integer, :Cont for continuous, :Bin for binary, :SemiCont for semicontinuous, and :SemiInt for semi-integer.

getvartype(m::AbstractMathProgModel)

Returns a vector indicating the types of each variable, with values described above.

addsos1!(m::AbstractMathProgModel, idx, weight)

Adds a special ordered set (SOS) constraint of type 1. Of the variables indexed by idx, at most one can be nonzero. The weight argument induces the ordering of the variables; as such, they should be unique values. A typical SOS1 constraint might look like \(y=\sum_i w_i x_i\), where \(x_i \in \{0,1\}\) are binary variables and the \(w_i\) are weights. See here for a description of SOS constraints and their potential uses.

addsos2!(m::AbstractMathProgModel, idx, weight)

Adds a special ordered set (SOS) constraint of type 2. Of the variables indexed by idx, at most two can be nonzero, and if two are nonzero, they must be adjacent in the set. The weight argument induces the ordering of the variables; as such, they should be unique values. A common application for SOS2 constraints is modeling nonconvex piecewise linear functions; see here for details.

Quadratic Programming

numquadconstr(m::AbstractMathProgModel)

Returns the number of quadratic constraints in the model.

setquadobj!(m::AbstractMathProgModel, Q)

Adds a quadratic term \(\frac{1}{2}x^TQx\) to the objective, replacing any existing quadratic terms. Note the implicit \(\frac{1}{2}\) scaling factor. The argument Q must be either a symmetric positive semidefinite matrix or the upper triangular portion of a symmetric positive semidefinite matrix (when minimizing). Sparse (CSC) or dense representations are accepted.

setquadobj!(m::AbstractMathProgModel, rowidx, colidx, quadval)

Adds a quadratic term \(\frac{1}{2}x^TQx\) to the objective, replacing any existing quadratic terms. Note the implicit \(\frac{1}{2}\) scaling factor. Here the entries of \(Q\) should be provided in sparse triplet form; e.g. entry indexed by k will fill quadval[k] in the (rowidx[k],colidx[k]) entry of matrix Q. Duplicate index sets (i,j) are accepted and will be summed together. Off-diagonal entries will be mirrored, so either the upper triangular or lower triangular entries of Q should be provided. If entries for both (i,j) and (j,i) are provided, these are considered duplicate terms. For example, setquadobj!(m, [1,1,2,2], [1,2,1,2], [3,1,1,1]) and setquadobj!(m, [1,1,2], [1,2,2], [3,2,1]) are both are valid descriptions for the matrix \(Q = \begin{pmatrix} 3 & 2 \\ 2 & 1 \end{pmatrix}\).

setquadobjterms!(m::AbstractMathProgModel, rowidx, colidx, quadval)

Provides an alternative “terms”-based interface to setquadobj!. A list of quadratic terms is specified instead of the matrix Q. For example, the objective \(x_1^2 + 2x_1x_2\) is specified by setquadobjterms!(m,[1,1],[1,2],[1.0,2.0]). Duplicate terms are summed together. Note: this method does not need to be implemented by solvers.

addquadconstr!(m::AbstractMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)

Adds the quadratic constraint \(s^Tx + \sum_{i,j} q_{i,j}x_ix_j \,\, sense \, rhs\) to the model. The linearidx and linearval arrays specify the sparse vector s. The quadratic terms are specified as in setquadobjterms! in the “terms” format. Sense must be '<', '>', or '='. If supported by the solver, addquadconstr! may also be used to specify second-order cone (SOCP) and rotated second-order cone constraints. These should be of the form \(x^Tx -y^2 \le 0\) or \(x^Tx -yz \le 0\), where \(y\) and \(z\) are restricted to be non-negative (in particular, \(Q\) can have at most one off-diagonal term).

getquadconstrsolution(m::AbstractMathProgModel)

Returns a vector containing the values of the quadratic constraints at the solution.

getquadconstrduals(m::AbstractMathProgModel)

Returns the Lagrangian dual solution vector corresponding to the quadratic constraints. Some solvers do not compute these values by default. Not available when integer variables are present.

getquadinfeasibilityray(m::AbstractMathProgModel)

Returns a “Farkas” proof of infeasibility, i.e., an unbounded ray of the dual, for the quadratic constraints. Note that for some solvers, one must specify additional options for this ray to be computed.

getquadconstrRHS(m::AbstractMathProgModel)

Returns a vector containing the right-hand side values on the quadratic constraints.

setquadconstrRHS!(m::AbstractMathProgModel, lb)

Sets the right-hand side values on the quadratic constraints. If the constraint was provided in the special second-order conic format, the solver may reject changing the right-hand side from zero.

Choosing solvers

Solvers and solver-specific parameters are specified by AbstractMathProgSolver objects, which are provided by particular solver packages. For example, the Clp package exports a ClpSolver object, which can be passed to linprog as follows:

using Clp
linprog([-1,0],[2 1],'<',1.5, ClpSolver())

Options are passed as keyword arguments, for example, ClpSolver(LogLevel=1). See the Clp, Cbc, GLPKMathProgInterface, and Gurobi packages for more information.

If no solver is specified, a default is chosen. See src/defaultsolvers.jl for the list of default solvers.

MIP Callbacks

MathProgBase supports a standardized and abstracted way to implement common MIP callbacks on the model. Currently there is support for adding:

  • Lazy constraints (only added to model if violated by integer-feasible solution)
  • Cut callbacks (only cuts off non-integer feasible solutions)
  • Heuristic callbacks (proposes heuristically constructed integer-feasible solutions at MIP nodes)

A more detailed description of the three types of supported callbacks can be found in the JuMP documentation here.

The MathProgSolverInterface exports an abstract type MathProgCallbackData which represents the solver-specific data needed to implement the callback.

If a callback function returns :Exit, the solver is expected to handle this and terminate the optimization.

setlazycallback!(m::AbstractMathProgModel, f)

Adds lazy constraint callback f to the model. Function f takes as argument only a MathProgCallbackData object.

setcutcallback!(m::AbstractMathProgModel, f)

Adds cut callback f to the model. Function f takes as argument only a MathProgCallbackData object.

setheuristiccallback!(m::AbstractMathProgModel, f)

Adds heuristic callback f to the model. Function f takes as argument only a MathProgCallbackData object.

setinfocallback!(m::AbstractMathProgModel, f)

Adds informational callback f to the model. Function f takes as argument only a MathProgCallbackData object.

cbgetmipsolution(d::MathProgCallbackData[, output])

Grabs current best integer-feasible solution to the model. The optional second argument specifies an output vector.

cbgetlpsolution(d::MathProgCallbackData[, output])

Grabs current best linear relaxation solution to the model. The optional second argument specifies an output vector.

cbgetobj(d::MathProgCallbackData)

Grabs objective value for current best integer-feasible solution.

cbgetbestbound(d::MathProgCallbackData)

Grabs best bound for objective function found so far (lower bound when minimizing, upper bound when maximizing).

cbgetexplorednodes(d::MathProgCallbackData)

Returns number of nodes that have been explored so far in the solve process.

cbgetstate(d::MathProgCallbackData)

Returns current location in solve process: :MIPNode if at node in branch-and-cut tree, :MIPSol at an integer-feasible solution, and :Other otherwise.

cbaddcut!(d::MathProgCallbackData, varidx, varcoef, sense, rhs)

Adds cut to model. The coefficient values are represented sparsely, with (one-indexed) indices in varidx and values in varcoef. The constraint sense sense is a character taking value <, >, or =, and the right-hand side value is rhs.

cbaddlazy!(d::MathProgCallbackData, varidx, varcoef, sense, rhs)

Adds lazy constraint to model. The coefficient values are represented sparsely, with (one-indexed) indices in varidx and values in varcoef. The constraint sense sense is a character taking value <, >, or =, and the right-hand side value is rhs.

cbaddsolution!(d::MathProgCallbackData)

Submit a (possibly partially defined) heuristic solution for the model. Should reset the solution stored in d to the original state at the start of callback.

cbsetsolutionvalue!(d::MathProgCallbackData, varidx, value)

Sets the value of a variable with (one-based) index varidx to value in the current partial solution being constructed by a user heuristic.

Conic Programming

Conic programming is an important class of convex optimization problems for which there exist specialized efficient solvers. We describe extensions to the MathProgSolverInterface for conic problems.

The design of this interface is inspired by the CBLIB format and the MOSEK modeling manual.

We consider the following primal problem to be in canonical conic form:

\[\begin{split}\min_{x}\, &c^Tx\\ s.t.\, &b - Ax \in K_1\\ &x \in K_2\\\end{split}\]

where \(K_1\) and \(K_2\) are cones (likely a product of a number of cones), with corresponding dual

\[\begin{split}\max_y\, &-b^Ty\\ s.t.\, &c + A^Ty \in K_2^*\\ &y \in K_1^*\end{split}\]

where \(K_1^*\) and \(K_2^*\) are the dual cones of \(K_1\) and \(K_2\), respectively.

The recognized cones are:

  • :Free, no restrictions (equal to \(\mathbb{R}^n\))
  • :Zero, all components must be zero
  • :NonNeg, the nonnegative orthant \(\{ x \in \mathbb{R}^n : x_i \geq 0, i = 1,\ldots,n \}\)
  • :NonPos, the nonpositive orthant \(\{ x \in \mathbb{R}^n : x_i \leq 0, i = 1,\ldots,n \}\)
  • :SOC, the second-order (Lorentz) cone \(\{(p,x) \in \mathbb{R} \times \mathbb{R}^{n-1} : ||x||_2^2 \leq p^2, p \geq 0\}\)
  • :SOCRotated, the rotated second-order cone \(\{(p,q,x) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{n-2} : ||x||_2^2 \leq 2pq, p \geq 0, q \geq 0\}\)
  • :SDP, the cone of symmetric positive semidefinite matrices \(\{ X \in \mathbb{R}^{n\times n} : X \succeq 0\}\)
  • :ExpPrimal, the exponential cone \(\operatorname{cl}\{ (x,y,z) \in \mathbb{R}^3 : y > 0, y e^{x/y} \leq z \}\)
  • :ExpDual, the dual of the exponential cone \(\{ (u,v,w) \in \mathbb{R}^3 : u < 0, -ue^{v/q} \leq ew\} \cup \{(0,v,w) : v \geq 0, w \geq 0\}\)

Not all solvers are expected to support all types of cones. However, when a simple transformation to a supported cone is available, for example, from :NonPos to :NonNeg or from :SOCRotated to :SOC, solvers should perform this transformation in order to allow users the extra flexibility in modeling.

loadconicproblem!(m::AbstractMathProgModel, c, A, b, constr_cones, var_cones)

Load the conic problem in primal form into the model. The parameter c is the objective vector, the parameter A is the constraint matrix (typically sparse), and the parameter b is the vector of “right-hand side” values. The parameters constr_cones and var_cones, which specify \(K_1\) and \(K_2\), are lists of (Symbol,indices) tuples, where Symbol is one of the above recognized cones and indices is a list of indices of constraints or variables (respectively) which belong to this cone (may be given as a Range). All variables and constraints must be listed in exactly one cone, and the indices given must correspond to the order of the columns and rows in the constraint matrix A. Cones may be listed in any order, and cones of the same class may appear multiple times. For the semidefinite cone, the number of variables or constraints present correspond to the lower (or upper) triangular elements in column-major (resp., row-major) order. Since an \(n \times n\) matrix has \(\frac{n(n+1)}{2}\) lower-triangular elements, by inverting this formula, when \(y\) elements are specified in indices, the corresponding matrix has \(\left(\sqrt{\frac{1}{4}+2y}-\frac{1}{2}\right) \times \left(\sqrt{\frac{1}{4}+2y}-\frac{1}{2}\right)\) elements.

getconicdual(m::AbstractMathProgModel)

If the solve was successful, returns the optimal dual solution vector \(y\).

supportedcones(m::AbstractMathProgSolver)

Returns a list of cones supported by the solver.

The solution vector, optimal objective value, termination status, etc. should be accessible from the standard methods, e.g., getsolution, getobjval, status, respectively.

Nonlinear Programming

MathProgBase provides an interface for nonlinear programming which is independent of both the solver and the user’s representation of the problem, whether using an algebraic modeling language or customized low-level code.

The diagram below illustrates MathProgBase as the connection between typical NLP solvers IPOPT, MOSEK, and KNITRO, and modeling languages such as JuMP and AMPL (via ampl.jl).

graph foo {
node [shape="box"];

subgraph clusterA {

 "IPOPT" -- "MOSEK" -- "KNITRO" -- "..." [style="invis", constraint="false"];
 label="Solvers";
 penwidth=0;

 }
 "IPOPT" -- "JuMP" [style="invis"];
"IPOPT" -- "MathProgBase";
"MOSEK" -- "MathProgBase";
"KNITRO" -- "MathProgBase";
"..." -- "MathProgBase";
"MathProgBase" -- "JuMP";
"MathProgBase" -- "AMPL";
"MathProgBase" -- "User";
"MathProgBase" -- x;

 subgraph clusterB {
 x [label="..."];
 "JuMP" -- "AMPL" -- "User" -- x [style="invis", constraint="false"];
 label="Modeling";
 penwidth=0;
 }
 rankdir=LR;
}

This structure also makes it easy to connect solvers written in Julia itself with user-provided instances in a variety of formats.

We take the prototypical format for a nonlinear problem as

\[\begin{split}\min_{x}\, &f(x)\\ s.t. &lb \leq g(x) \leq ub\\ &l \leq x \leq u\\\end{split}\]

Where \(x \in \mathbb{R}^n, f: \mathbb{R}^n \to \mathbb{R}, g: \mathbb{R}^n \to \mathbb{R}^m\), and vectors \(lb \in \mathbb{R}^m \cup \{-\infty\}, ub \in \mathbb{R}^m \cup \{\infty\},l \in \mathbb{R}^n \cup \{-\infty\}, u \in \mathbb{R}^n \cup \{\infty\}\).

The objective function \(f\) and constraint function \(g\) may be nonlinear and nonconvex, but are typically expected to be twice differentiable.

Below we describe extensions to the MathProgSolverInterface for these nonlinear programming problems.

loadnonlinearproblem!(m::AbstractMathProgModel, numVar, numConstr, l, u, lb, ub, sense, d::AbstractNLPEvaluator)

Loads the nonlinear programming problem into the model. The parameter numVar is the number of variables in the problem, numConstr is the number of constraints, l contains the variable lower bounds, u contains the variable upper bounds, lb contains the constraint lower bounds, and ub contains the constraint upper bounds. Sense contains the symbol :Max or :Min, indicating the direction of optimization. The final parameter d is an instance of an AbstractNLPEvaluator, described below, which may be queried for evaluating \(f\) and \(g\) and their corresponding derivatives.

The abstract type AbstractNLPEvaluator is used by solvers for accessing the objective function \(f\) and constraints \(g\). Solvers may query the value, gradients, Hessian-vector products, and the Hessian of the Lagrangian.

initialize(d::AbstractNLPEvaluator, requested_features::Vector{Symbol})

Must be called before any other methods. The vector requested_features lists features requested by the solver. These may include :Grad for gradients of \(f\), :Jac for explicit Jacobians of \(g\), :JacVec for Jacobian-vector products, :HessVec for Hessian-vector and Hessian-of-Lagrangian-vector products, :Hess for explicit Hessians and Hessian-of-Lagrangians, and :ExprGraph for expression graphs.

features_available(d::AbstractNLPEvaluator)

Returns the subset of features available for this problem instance, as a list of symbols in the same format as in initialize.

eval_f(d::AbstractNLPEvaluator, x)

Evaluate \(f(x)\), returning a scalar value.

eval_g(d::AbstractNLPEvaluator, g, x)

Evaluate \(g(x)\), storing the result in the vector g which must be of the appropriate size.

eval_grad_f(d::AbstractNLPEvaluator, g, x)

Evaluate \(\nabla f(x)\) as a dense vector, storing the result in the vector g which must be of the appropriate size.

jac_structure(d::AbstractNLPEvaluator)

Returns the sparsity structure of the Jacobian matrix \(J_g(x) = \left[ \begin{array}{c} \nabla g_1(x) \\ \nabla g_2(x) \\ \vdots \\ \nabla g_m(x) \end{array}\right]\) where \(g_i\) is the \(i\text{th}\) component of \(g\). The sparsity structure is assumed to be independent of the point \(x\). Returns a tuple (I,J) where I contains the row indices and J contains the column indices of each structurally nonzero element. These indices are not required to be sorted and can contain duplicates, in which case the solver should combine the corresponding elements by adding them together.

hesslag_structure(d::AbstractNLPEvaluator)

Returns the sparsity structure of the Hessian-of-the-Lagrangian matrix \(\nabla^2 f + \sum_{i=1}^m \nabla^2 g_i\) as a tuple (I,J) where I contains the row indices and J contains the column indices of each structurally nonzero element. These indices are not required to be sorted and can contain duplicates, in which case the solver should combine the corresponding elements by adding them together. Any mix of lower and upper-triangular indices is valid. Elements (i,j) and (j,i), if both present, should be treated as duplicates.

eval_jac_g(d::AbstractNLPEvaluator, J, x)

Evaluates the sparse Jacobian matrix \(J_g(x) = \left[ \begin{array}{c} \nabla g_1(x) \\ \nabla g_2(x) \\ \vdots \\ \nabla g_m(x) \end{array}\right]\). The result is stored in the vector J in the same order as the indices returned by jac_structure.

eval_jac_prod(d::AbstractNLPEvaluator, y, x, w)

Computes the Jacobian-vector product \(J_g(x)w\), storing the result in the vector y.

eval_jac_prod_t(d::AbstractNLPEvaluator, y, x, w)

Computes the Jacobian-transpose-vector product \(J_g(x)^Tw\), storing the result in the vector y.

eval_hesslag_prod(d::AbstractNLPEvaluator, h, x, v, σ, μ)

Given scalar weight σ and vector of constraint weights μ, computes the Hessian-of-the-Lagrangian-vector product \(\left(\sigma\nabla^2 f(x) + \sum_{i=1}^m \mu_i \nabla^2 g_i(x)\right)v\), storing the result in the vector h.

eval_hesslag(d::AbstractNLPEvaluator, H, x, σ, μ)

Given scalar weight σ and vector of constraint weights μ, computes the sparse Hessian-of-the-Lagrangian matrix \(\sigma\nabla^2 f(x) + \sum_{i=1}^m \mu_i \nabla^2 g_i(x)\), storing the result in the vector H in the same order as the indices returned by hesslag_structure.

isobjlinear(d::AbstractNLPEvaluator)

true if the objective function is known to be linear, false otherwise.

isobjquadratic(d::AbstractNLPEvaluator)

true if the objective function is known to be quadratic (convex or nonconvex), false otherwise.

isconstrlinear(d::AbstractNLPEvaluator, i)

true if the \(i\text{th}\) constraint is known to be linear, false otherwise.

obj_expr(d::AbstractNLPEvaluator)

Returns an expression graph for the objective function as a standard Julia Expr object. All sums and products are flattened out as simple Expr(:+,...) and Expr(:*,...) objects. The symbol x is used as a placeholder for the vector of decision variables. No other undefined symbols are permitted; coefficients are embedded as explicit values. For example, the expression \(x_1+\sin(x_2/\exp(x_3))\) would be represented as the Julia object :(x[1] + sin(x[2]/exp(x[3]))). See the Julia manual for more information on the structure of Expr objects. There are currently no restrictions on recognized functions; typically these will be built-in Julia functions like ^, exp, log, cos, tan, sqrt, etc., but modeling interfaces may choose to extend these basic functions.

constr_expr(d::AbstractNLPEvaluator, i)

Returns an expression graph for the \(i\text{th}\) constraint in the same format as described above. The head of the expression is :comparison, indicating the sense of the constraint. The right-hand side of the comparison must be a constant; that is, :(x[1]^3 <= 1) is allowed, while :(1 <= x[1]^3) is not valid. Double-sided constraints are allowed, in which case both the lower bound and upper bounds should be constants; for example, :(-1 <= cos(x[1]) + sin(x[2]) <= 1) is valid.

The solution vector, optimal objective value, termination status, etc. should be accessible from the standard methods, e.g., getsolution, getobjval, status, respectively. Nonlinear solvers may also provide optimal Lagrange multipliers if available through getreducedcosts and getconstrduals.