Skip to content

Instantly share code, notes, and snippets.

@mlubin
Last active December 14, 2015 21:59
Show Gist options
  • Save mlubin/5154771 to your computer and use it in GitHub Desktop.
Save mlubin/5154771 to your computer and use it in GitHub Desktop.
Prototype design for linear(/mixed-integer/quadradic/+) programming interface in Julia

There are four different levels at which users may interact with LP solvers (top to bottom):

  • Algebraic modeling languages such as MathProg.jl
  • Matlab-style linprog functions, provide problem as input and retrieve solution in one function call. Maybe put this in Optim.jl.
  • Standard low-level interface. This is a state-based interface with a solver object where the user may dynamically add and remove variables and constraints, repeatedly resolve the problem with hotstarts (important for decomposition algorithms), etc. Inspired by COIN-OR's OSI project, but hopefully we can improve on some of their mistakes.
  • Solver-specific low-level interface. Typically a lightweight wrapper around the solver's C interface. Names should follow the solver's convention.

Implementers of solver interfaces (or solvers written in Julia) are responsible for providing the standard low-level interface (not linprog). The linprog interface will be built on top of the standard low-level interface.

High-level user-facing linprog function:

solution = linprog(c, A, sense, b, lb, ub[, kwargs])

where

  • c is the objective vector, always in the sense of minimization
  • A is the constraint matrix
  • sense is the vector of constraint sense characters: '<', '=', and '>'.
  • b is the right-hand side vector.
  • lb is the vector of lower bounds on the variables
  • ub is the vector of upper bounds on the variables

A scalar is accepted for the b, sense, lb, and ub arguments, in which case its value is replicated. -Inf and Inf should be treated correctly when given for lb or ub.

A shortened version should be available as

linprog(c, A, b, sense[, kwargs]) = linprog(c, A, b, sense, 0, Inf[, kwargs])

Second version based on range constraints:

solution = linprog(c, A, rowlb, rowub, lb, ub[, kwargs])

where

  • c is the objective vector, always in the sense of minimization
  • A is the constraint matrix
  • rowlb is the vector of row lower bounds.
  • rowub is the vector of row upper bounds.
  • lb is the vector of lower bounds on the variables
  • ub is the vector of upper bounds on the variables

kwargs contains solution parameters as keyword arguments. We'll have to wait and see what the syntax will look like, but this will be in Julia soon enough.

solution should be an instance of some type with the following members

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, maybe others.

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

objval - optimal objective value

sol - primal solution vector

attrs - a dictionary to 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)

colbasis - optimal simplex basis statuses for the variables (columns) if available. Possible values are "NonbasicAtLower", "NonbasicAtUpper", "Basic", "Superbasic"

rowbasis - optimal simplex basis statuses for the constraints (rows) if available. Same statuses.


By convention the dual multipliers should have the sign following the interpretation of marginal change in the objective when the corresponding active right-hand side is increased. This corresponds to the standard that reduced costs should be nonnegative when a variable is at a lower bound and nonpositive when a variable is at an upper bound. Different solvers might have different conventions for the lambda vector, so transformations might be needed.


@lindahua proposes extending linprog with a more convenient approach to specify different types of constraints:

for example, we have this set of linear constraints:

a <= A x <= b
Cx <= c
Aeq x == beq

It is possible to merge these into one constraint as

[a; -inf; beq] <= [A; C; Aeq] x <= [b; c; beq] 

But let the user/caller do this may not be the optimal interface.

On the caller side, an interface can be

linprog(f, (a, '<', A, '<', b), (C, '<', c), (Aeq, '=', aeq), ...)

Proposal for standardized low-level LP interface:

Conventions:

  • All indices are one-based.
  • It should be easy to access the solver-specific model object in order to call solver-specific functions.

model([params])

Create a new empty model, optionally passing a dictionary of parameters.

Notes:

  • Solver-specific parameters are allowed. Eventually we will also have a standard set of parameters.
  • Maybe use keyword arguments instead of a dictionary? We'll have to wait and see what keyword arguments look like when implemented in Julia.
  • In order to load a problem from a file in Gurobi, the file name needs to be passed to the GRBModel constructor. I suggest that this function create a GRBEnv object and wait until later to create a GRBModel.

loadproblem(m, filename::String)

Load a problem from a file into the model m. File extension should be used to determine the format.

Notes:

  • Different solvers support different input formats. An error should be thrown if an unsupported format is used.
loadproblem(m, A, collb, colub, obj, rowlb, rowub)

Load a problem from the given input data. A is the constraint matrix, collb are lower bounds on variables, colub are upper bounds on variables, rowlb are lower bounds on constraints, rowub are upper bounds on constraints.

Notes:

  • The solver may choose what type of matrix to accept for A. As a minimum, this should include Array{Float64,2} and SparseMatrixCSC{Float64,Int}. Exact solvers may want to allow matrices with rational entries.
  • Inf and -Inf indicate no corresponding upper or lower bound, respectively. Equal lower and upper bounds are used to indicate equalities.

writeproblem(m, filename::String)

Write problem to a file. File extension should be used to determine the format.

Notes:

  • Different solvers support different output formats. An error should be thrown if an unsupported format is used.

getvarLB(m)
setvarLB(m, collb)
getvarUB(m)
setvarUB(m, colub)
getconstrLB(m)
setconstrLB(m, rowlb)
getconstrUB(m)
setconstrUB(m, rowub)
getobj(m)
setobj(m, obj)

Setters and getters for vector properties.


addvar(m, rowidx, rowcoef, collb, colub, objcoef)

Add a column to the problem. rowcoef specifies non-zero coefficients in the new column in corresponding indices specified in rowcoef. collb, colub, and objcoef are scalars specifying the lower bound, upper bound, and objective coefficient for the new variable.


addconstr(m, colidx, colcoef, rowlb, rowub)

Add a row to the problem.


updatemodel(m)

Update the model after making changes. Only some solvers need this.


setsense(m,sense)
getsense(m)

Set/get problem sense. Either :Min or :Max.


numvar(m)
numconstr(m)

Get number of variables and constraints in the model.


optimize(m)

Solve (or resolve) the current model.


status(m)

Get termination status symbol, one of :Optimal, :PrimalInfeasible, :DualInfeasible, :UserLimit (iteration limit or timeout), :Error, maybe others.


getobjval(m)

Optimal objective value.

getsolution(m)

Primal solution.

getconstrsolution(m)

Constraint matrix times solution vector, also known as "row activity".

getreducedcosts(m)

Dual solution corresponding to (active) variable bounds.

getconstrduals(m)

Dual solution corresponding to (active) row constraints.

Notes:

  • By convention the dual values should have the sign following the interpretation of marginal change in the objective when the corresponding active right-hand side is increased. This corresponds to the standard that reduced costs should be nonnegative when a variable is at a lower bound and nonpositive when a variable is at an upper bound. Different solvers might have different conventions for the row duals, so transformations might be needed.

getrawsolver(m)

Return a solver-specific object that can be used directly with the solver's low-level interface.


setvartype(m,v)
getvartype(m)

Get/set vector of variable types. v should be a Char vector where 'C' indicates continuous and 'I' indicates integer. (Binaries should be specified as 'I' with lower bound 0 and upper bound 1).


What's missing:

  • Setting/getting solution method
  • Add multiple rows/columns. Remove rows/columns.
  • Setting/getting simplex basis (if applicable)
  • Standard subset of parameters (Timeout, output level, iteration limit, etc.)
  • SOS constraints for integer problems
  • Quadratic objectives
@lindahua
Copy link

what about calling variables as "var" instead of "col", and calling "constraints" and "constraints" instead of "row"?

This sounds more like terminologies of optimization ("col" and "row" are more like linear algebraic things). Then the function names would be "getvarlb", "getvarub", (or more descriptive, "get_variable_lbound", "get_variable_ubound".

Except for the naming, everything else looks good.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment