Linear and integer programming

Solve a linear programme exactly over the rationals, with a certificate.

minimise (or, with sense="max", maximise) c @ x subject to A_ub @ x <= b_ub, A_eq @ x == b_eq and bounds — the same calling convention as scipy.optimize.linprog, but every coefficient must be exact (int, Fraction, "p/q" or a rational Expr) and every result is exact. bounds is one (lower, upper) pair for all variables or one per variable; None (or ±inf) is no bound; the default is (0, None).

method picks the simplex: "dense" — a two-phase tableau with Bland's rule; "revised" — a revised simplex keeping A sparse with an exact product-form factorisation of the basis (Dantzig pricing with a Bland fallback, so it cannot cycle either), much faster on medium sparse programmes; "auto" (default) — dense for at most 12 rows, revised otherwise. basis warm-starts the revised method from a previous result's LPResult.basis (or the result itself): a basis still primal feasible goes straight to phase 2, one still dual feasible (after a right-hand-side change or added rows) runs the dual simplex, anything else falls back to a cold start (warm_started says which). Both honour the active alkahest.context(budget=...) at every pivot.

LPResult — status with its certificate (dual multipliers,

Farkas vector, or ray), already verified exactly, plus basis,

iterations, method and warm_started.

LinprogError

E-LP-001 malformed shapes, E-LP-002 a float or other inexact coefficient, E-LP-006 a certificate that failed its own check (a bug; the result is withheld).

BudgetExceededError

the active budget tripped (E-BUDGET-*).

>>> r = alkahest.linprog([1, 1], A_ub=[[1, 2], [3, 1]], b_ub=[4, 6], sense="max")
>>> r.fun, r.x
(Fraction(14, 5), [Fraction(8, 5), Fraction(6, 5)])

Solve a mixed-integer linear programme exactly by branch and bound.

Arguments as for linprog(), plus integrality — one flag per variable (1/True integer, 0/False continuous, SciPy's convention), or a single bool for all; the default makes every variable an integer. Binary variables are integers with bounds=(0, 1).

Before branching, up to cut_rounds rounds of Gomory mixed-integer cuts (at most max_cuts in all) strengthen the root relaxation (cut-and-branch). The cuts are read off the exact optimal tableau, so each one is valid: it never removes an integer-feasible point (they are listed in MILPResult.cuts). Pass cut_rounds=0 to branch on the plain relaxation. The search is then depth-first until a first integer point is found, then best-bound; every node LP is solved exactly, warm-started from its parent's basis (warm_start), and its certificate checked. max_nodes caps the number of node LPs, and the active budget is honoured at every node and pivot.

MILPResult — "optimal" with gap == 0,

"infeasible", or "unbounded" (integer point plus integral ray).

LinprogError

E-LP-005 when max_nodes is reached — the error carries incumbent (best integer point, unproven, or None), incumbent_fun, bound (best open bound) and nodes; and the codes of linprog().

BudgetExceededError

the budget tripped; carries incumbent, incumbent_fun and bound likewise.

>>> r = alkahest.milp([10, 13, 7, 8], A_ub=[[4, 6, 3, 5]], b_ub=[12],
...                   bounds=(0, 1), sense="max")
>>> r.fun, r.x
(25, [1, 0, 1, 1])

Solve a linear (or mixed-integer) programme stated with expressions.

objective is an affine Expr in vars; constraints is a list of relations built with pool.le, pool.ge or pool.pred_eq (or an And of them) between affine expressions with rational coefficients. The variables are free unless a constraint bounds them (write pool.ge(x, 0) for x ≥ 0). Each <=/>= becomes one A_ub row and each == one A_eq row, in the order given, so ineqlin and eqlin line up with constraints (result.program() shows the rows).

integer — None for an LP; True to make every variable an integer; or a list of the variables that must be integers (then the result is a MILPResult, and max_nodes applies).

parameter — a symbol θ that may appear affinely either in the objective's coefficients ((2 + θ)*x) or in the constraints' constant terms (x + y <= 4 + θ); the result is then a ParametricLPResult over parameter_range (a (lo, hi) pair, None for an infinite end), as from linprog_parametric().

LPResult, MILPResult or ParametricLPResult;

result.solution (solution_at(θ) for a parametric result) maps

each variable to its value.

LinprogError

E-LP-003 non-affine objective/constraint or a symbol that is not a decision variable; E-LP-002 a float or irrational coefficient; E-LP-004 a strict inequality or a non-relation; E-LP-007 an unsupported parametric form (the parameter multiplying a variable inside a constraint, entering non-linearly, moving the objective and the right-hand side at once, or combined with integer); and the codes of linprog() / milp().

>>> p = alkahest.ExprPool(); x, y = p.symbol("x"), p.symbol("y")
>>> r = alkahest.linprog_expr(3*x + 2*y,
...     [p.le(x + y, p.integer(4)), p.le(x + 3*y, p.integer(6)),
...      p.ge(x, p.integer(0)), p.ge(y, p.integer(0)), p.le(x, p.integer(3))],
...     [x, y], sense="max")
>>> r.fun
11

Solve a linear programme whose right-hand side or objective is affine in one parameter θ, exactly, for every θ at once.

Arguments as for linprog(), plus the θ-direction: either c_theta (the objective is (c + θ·c_theta) @ x) or b_ub_theta / b_eq_theta (the right-hand sides are b_ub + θ·b_ub_theta and b_eq + θ·b_eq_theta; a missing one is zero) — not both. theta_range is a (lo, hi) pair, None (or ±inf) for an infinite end; the default is the whole line.

The optimal value is piecewise affine in θ (convex or concave); the result lists the exact breakpoints and, for each critical region, the optimal basis, value, point and multipliers as affine functions of θ, each region's certificate verified exactly for every θ in it. Regions where the programme is infeasible (RHS case) or unbounded (objective case) are reported as such.

ParametricLPResult.

LinprogError

E-LP-007 when both c_theta and a right-hand-side direction are given (a parameter in the constraint matrix cannot be stated here; linprog_expr() refuses it with the same code); E-LP-001 a direction of the wrong length or an empty theta_range; and the codes of linprog().

BudgetExceededError

the active budget tripped.

>>> r = alkahest.linprog_parametric([1, 1], A_ub=[[1, 2], [3, 1]], b_ub=[4, 6],
...                                 b_ub_theta=[1, 0], sense="max")
>>> r.breakpoints
[-4, -2, 8]
>>> r.value_at(0), r.value_at(10)
(Fraction(14, 5), 6)

Certified result of linprog() / linprog_expr().

status is "optimal", "infeasible" or "unbounded"; every number is an int or fractions.Fraction. The certificate for the status is attached and verify() re-checks it exactly:

• optimal — x, fun and the dual multipliers ineqlin, eqlin, lower, upper (SciPy's marginals convention) with c = A_ub.T @ ineqlin + A_eq.T @ eqlin + lower + upper, complementary slackness, and dual objective equal to fun;
• infeasible — farkas: a dict of non-negative weights (ineqlin, lower, upper) and free weights (eqlin) combining the constraints into 0 <= negative;
• unbounded — x a feasible point and ray an improving direction.

status

success

x

fun

ineqlin

eqlin

lower

upper

farkas

ray

solution

basis

iterations

method

warm_started

program()

verify()

Result of milp().

status is "optimal" (the branch-and-bound tree closed: gap is exactly 0 and bound == fun), "infeasible" (every leaf was LP-infeasible with a verified Farkas vector), or "unbounded" (x is an integer-feasible point and ray a direction, integral on the integer variables, along which the objective improves without bound). A search that hit max_nodes or the active budget raises instead, carrying the incumbent.

status

success

x

fun

bound

gap

nodes

ray

cuts

solution

program()

verify()

Result of linprog_parametric() (or linprog_expr() with a parameter): the critical regions of θ, left to right.

Each region is a dict: theta — (lo, hi) (None for an infinite end); status — "optimal" (closed region, one basis optimal throughout), "infeasible" or "unbounded" (open at the end shared with an optimal region); and, when optimal, value, x, ineqlin, eqlin, lower, upper as affine pairs (a, b) meaning a + θ·b (value a pair of numbers, the rest pairs of lists), and the region's basis. Every optimal region's certificate was verified exactly for all θ in it (at its ends, and through the certificate's slope for an infinite end).

kind

regions

breakpoints

value_at(theta)

status_at(theta)

x_at(theta)

solution_at(theta)

program()

verify()