milp

README

go-milp

A pure-Go mixed-integer linear programming (MILP) solver. Define linear and mixed-integer linear programs — variables with bounds and a kind (continuous, integer, or binary), a linear objective to minimize or maximize, and linear constraints (≤, ≥, =) — and solve them for the optimal objective, the variable values, and a status (optimal, infeasible, or unbounded).

The LP relaxation is solved with a native, pure-Go sparse revised simplex (bounded-variable, with a Gilbert–Peierls sparse-LU basis factorization and product-form updates); integrality is enforced by branch-and-bound with warm-started re-solves and pseudocost branching over the integer and binary variables.

The pure-Go, zero-runtime-dependency guarantee

The default build is pure Go and links zero C, and zero third-party runtime dependencies:

CGO_ENABLED=0 go build ./...   # succeeds; no cgo, no GLPK, no gonum

The solver core (LP relaxation + branch-and-bound) is entirely native Go over the standard library. gonum and testify are test-only: gonum is retained solely as an in-tree differential cross-check (its dense simplex is run against the native core on random LPs); it is not imported by any non-test code, so a consumer importing go-milp does not pull it into their build.

GLPK appears only as a differential-test parity oracle — compiled via cgo behind the glpk_oracle build tag, in the internal/parity_tests package, to check the native solver against a trusted reference. It is never a runtime backend. A default build never compiles or links it.

Install

go get github.com/daniel-sullivan/go-milp

Requires Go 1.26+.

Quickstart

package main

import (
	"log"

	"github.com/daniel-sullivan/go-milp"
)

func main() {
	// maximize 3x + 2y  s.t.  x + y <= 4,  x + 3y <= 6,  x,y >= 0
	m := milp.NewModel()
	x := m.AddContinuous(0, milp.Inf)
	y := m.AddContinuous(0, milp.Inf)
	m.SetObjective(milp.Maximize, []milp.Term{{x, 3}, {y, 2}})
	m.AddConstraint([]milp.Term{{x, 1}, {y, 1}}, milp.LessEq, 4)
	m.AddConstraint([]milp.Term{{x, 1}, {y, 3}}, milp.LessEq, 6)

	sol, err := m.Solve()
	if err != nil {
		log.Fatal(err)
	}
	log.Printf("status=%v objective=%.2f x=%.2f y=%.2f",
		sol.Status, sol.Objective, sol.Value(x), sol.Value(y))
}

A runnable 0/1 knapsack example lives in examples/knapsack.

API surface

• Variables: AddContinuous(lo, hi), AddInteger(lo, hi), AddBinary(), AddVar(lo, hi, kind). Use milp.Inf / milp.NegInf for absent bounds.
• Objective: SetObjective(sense, []Term) with milp.Minimize / milp.Maximize.
• Constraints: AddConstraint([]Term, rel, rhs) with milp.LessEq, milp.GreaterEq, milp.EqualTo.
• Solve: Solve() (defaults) or SolveWith(Options). Implements the Solver interface. SolveWithStats(Options) additionally returns SolveStats (branch-and-bound nodes, LP relaxations solved, simplex pivots, factorizations, warm starts, cuts, cut rounds).
• Options: IntegralityTol, LPTolerance, NodeLimit, MaxDepth, TimeLimit, RelativeGap, AbsoluteGap (anytime / gap-limited solves), IncumbentCallback (called on each new incumbent), DisableCuts, DisablePresolve.
• Result: Solution{Status, Objective, Values, BestBound, Gap} with Value(v). BestBound is the proven dual bound and Gap the relative optimality gap (0 when proven optimal); on a TimeLimit / NodeLimit stop the incumbent is returned with the gap it was proven to.

How it works

1. Computational (logical-variable) form. Each bounded variable keeps its bounds natively (no standard-form expansion), and each constraint row gets one logical variable whose bounds encode the relation (≤, ≥, =). The system is [A | −I]·x = 0 over n + m bounded columns, with a trivial all-logicals starting basis. This avoids gonum's row-per-upper-bound blow-up and free- variable splitting.
2. Native sparse revised simplex. A bounded-variable primal simplex prices with Dantzig (Bland's rule as an anti-cycling fallback), enters/leaves via a bounded ratio test with bound flips, and reaches feasibility with a composite Phase-1 (sum of infeasibilities). Each basis solve (FTRAN/BTRAN) goes through a Gilbert–Peierls sparse LU with partial pivoting; basis changes are applied as product-form eta updates with periodic refactorization, so the O(m³) dense refactorization of a classic dense simplex is avoided entirely.
3. Branch-and-bound. The root relaxation is tightened by presolve (activity-based bound propagation, singleton/duplicate/redundant rows, binary coefficient tightening) and root Gomory (GMI) cutting planes, then searched depth-first. Child nodes warm-start from the parent's optimal basis, pseudocost branching (Achterberg product rule) picks the branching variable, and a rounding/dive heuristic seeds a strong incumbent. A global dual bound (min over the open frontier) yields a proven optimality gap, so Options.TimeLimit / RelativeGap / AbsoluteGap give an anytime answer. SolveWithStats reports nodes, pivots, refactorizations, and cuts.

The design and rationale are documented in docs/simplex.md.

Verification tiers

Tier	What it checks	How to run
Unit + stress	Known LP/MILP optima, infeasible/unbounded detection, negative and free variables, self-checked feasibility over generated LPs/MILPs and battery-scheduling problems, sparse-LU vs dense and eta-update vs refactor equality, GMI cut validity (no-cut optimum satisfies every cut), and a presolve round-trip fuzz (no feasible point lost). Small MILPs are cross-checked against exhaustive enumeration.	mise run test (or CGO_ENABLED=0 go test ./...)
GLPK parity	The native solver vs GLPK on generated LPs, MILPs, and an battery-scheduling problem. ε-equivalent, not bit-exact: it asserts the objective matches within tolerance and the native solution is feasible — never that the two solution vectors are equal (LP/MILP have multiple valid optima).	mise run oracle then mise run parity
Lint	gofmt, go vet, and llms.txt freshness.	mise run fmt vet llms:check

Why parity is ε-equivalent

Two correct simplex / branch-and-bound implementations routinely reach different optimal vertices with the same optimal cost. The invariant that matters is therefore the objective value (within objAbsTol = 1e-5 + 1e-6·|obj|) plus feasibility of the returned decision, not equality of the solution vectors. The parity tests assert exactly that (see internal/parity_tests/milp/parity_test.go).

Performance

The native sparse revised simplex replaces a dense O(m³)-per-iteration simplex with a sparse-LU basis factorization and product-form updates, so the LP core scales to problems with over a thousand constraints in well under a second, and branch-and-bound warm-starts every node from its parent basis.

Benchmarks on the battery-scheduling model (a battery-scheduling MILP), Apple M-series, CGO_ENABLED=0, best of 3. The earlier dense backend (gonum) is shown for reference; GLPK 5.0 (the cgo parity oracle, not a runtime backend) is shown as a state-of-the-art baseline. Native matches GLPK's optimal objective to ≤1e-13 at every size.

Full MILP solve:

horizon T	rows	vars	dense (gonum)	native sparse	GLPK
6	55	61	65 ms	114 µs	77 µs
12	109	121	5.30 s	1.44 ms	0.44 ms
24	217	241	6.06 s	2.14 ms	0.29 ms
48	434	481	did not finish	22.7 ms	3.11 ms
144	1300	1441	did not finish	445 ms	54 ms

LP relaxation only (the LP core):

horizon T	rows	vars	native sparse	GLPK
48	434	481	7.65 ms	1.11 ms
144	1300	1441	85 ms	9.0 ms

At T=144 both the full MILP (445 ms) and the LP relaxation (85 ms) solve well under a second — roughly 2800–3700× faster than the dense backend at the sizes it could still finish, and within ~10× of GLPK. On integer-heavy MILPs the root Gomory cuts close ~63% of the root gap and cut the branch-and-bound tree ~2.3× (several instances solve at the root); on the already-tight energy shape cuts are a near-no-op. Correctness is proven by the GLPK parity gate at every size including T=144, plus a cut validity-insurance test and a presolve round-trip fuzz.

License

MIT — see LICENSE. GLPK (GPL) is used only as an optional, build-tagged test oracle and is never distributed or linked in a default build.

Documentation

Overview

Package milp is a pure-Go mixed-integer linear programming solver.

It defines and solves linear programs (LP) and mixed-integer linear programs (MILP): decision variables with lower/upper bounds and a kind (continuous, integer, or binary), a single linear objective to minimize or maximize, and linear constraints (<=, >=, =). Solve returns the optimal objective value, the variable values, and a status.

The LP relaxation is solved with a native, pure-Go sparse revised simplex (bounded-variable, with a Gilbert–Peierls sparse-LU basis factorization and product-form updates); integrality is enforced by a branch-and-bound search with warm-started re-solves and pseudocost branching over the integer and binary variables. The default build is pure Go and links zero C and zero third-party runtime dependencies: CGO_ENABLED=0 go build ./... succeeds.

GLPK (via cgo) is used only as a parity oracle in the build-tagged parity tests under internal/parity_tests; it is never a runtime backend. gonum is a test-only differential cross-check and is not imported by any runtime code.

Concurrency

A Model is not safe for concurrent modification. Solve does not mutate the Model and may be called concurrently on the same Model from multiple goroutines; each call allocates its own working state.

Model building

Construct a Model, add variables, set the objective, add constraints, and Solve:

m := milp.NewModel()
x := m.AddInteger(0, 10)
y := m.AddContinuous(0, milp.Inf)
m.SetObjective(milp.Maximize, []milp.Term{{x, 3}, {y, 2}})
m.AddConstraint([]milp.Term{{x, 1}, {y, 1}}, milp.LessEq, 4)
sol, err := m.Solve()

Index

• Variables
• type Constraint
  • func (c Constraint) ID() int
• type Model
  • func NewModel() *Model
  • func (m *Model) AddBinary() Var
  • func (m *Model) AddConstraint(terms []Term, rel Relation, rhs float64) Constraint
  • func (m *Model) AddContinuous(lower, upper float64) Var
  • func (m *Model) AddInteger(lower, upper float64) Var
  • func (m *Model) AddVar(lower, upper float64, kind VarKind) Var
  • func (m *Model) NumConstraints() int
  • func (m *Model) NumVars() int
  • func (m *Model) SetObjective(sense Sense, terms []Term)
  • func (m *Model) Solve() (*Solution, error)
  • func (m *Model) SolveWith(opts Options) (*Solution, error)
  • func (m *Model) SolveWithStats(opts Options) (*Solution, SolveStats, error)
• type Options
  • func DefaultOptions() Options
• type Relation
  • func (r Relation) String() string
• type Sense
  • func (s Sense) String() string
• type Solution
  • func (s *Solution) Value(v Var) float64
• type SolveStats
• type Solver
  • func NewSolver(opts Options) Solver
• type Status
  • func (s Status) String() string
• type Term
• type Var
  • func (v Var) ID() int
• type VarKind
  • func (k VarKind) String() string

Variables

var (
	// ErrNoObjective is returned by Solve when no objective has been set.
	ErrNoObjective = errors.New("milp: no objective set")
	// ErrNoVariables is returned by Solve when the model has no variables.
	ErrNoVariables = errors.New("milp: model has no variables")
	// ErrInvalidBounds is returned by an AddVar call whose lower bound exceeds
	// its upper bound, or whose bound is NaN.
	ErrInvalidBounds = errors.New("milp: variable lower bound exceeds upper bound")
	// ErrUnknownVar is returned when a Term references a variable that does not
	// belong to the model being solved.
	ErrUnknownVar = errors.New("milp: term references a variable from another model")
	// ErrInvalidCoef is returned when a coefficient or right-hand side is NaN or
	// infinite.
	ErrInvalidCoef = errors.New("milp: coefficient or right-hand side is not finite")
	// ErrNodeLimit is returned by Solve when the branch-and-bound node limit was
	// reached before the search completed. Any incumbent found is still returned
	// in the Solution.
	ErrNodeLimit = errors.New("milp: branch-and-bound node limit reached")
	// ErrTimeLimit is returned by Solve when Options.TimeLimit elapsed before the
	// search completed. Any incumbent found is still returned in the Solution,
	// with BestBound and Gap reporting the proven optimality gap.
	ErrTimeLimit = errors.New("milp: branch-and-bound time limit reached")
	// ErrNumeric is returned when the underlying LP solver fails for numeric or
	// structural reasons (singular basis, cycling, ill-conditioning) that are
	// neither infeasibility nor unboundedness.
	ErrNumeric = errors.New("milp: linear relaxation failed numerically")
)

Sentinel errors returned by the milp package. Every message is prefixed with "milp:".

var (
	Inf    = math.Inf(1)
	NegInf = math.Inf(-1)
)

Inf and NegInf are the positive and negative infinities used to express an absent upper or lower bound on a variable.

Types

type Constraint

type Constraint struct {
	// contains filtered or unexported fields
}

Constraint is an opaque handle to a linear constraint within one Model.

func (Constraint) ID

func (c Constraint) ID() int

ID returns the constraint's zero-based index within its Model.

type Model

type Model struct {
	// contains filtered or unexported fields
}

Model builds a linear or mixed-integer linear program. Construct one with NewModel, add variables and constraints, set an objective, then Solve.

A Model is not safe for concurrent modification. Solve is read-only and may be called concurrently on the same Model.

func NewModel

func NewModel() *Model

NewModel returns an empty Model.

func (*Model) AddBinary

func (m *Model) AddBinary() Var

AddBinary adds a binary variable constrained to {0, 1}.

func (*Model) AddConstraint

func (m *Model) AddConstraint(terms []Term, rel Relation, rhs float64) Constraint

AddConstraint adds the linear constraint (sum of terms) rel rhs and returns a handle to it. Terms referencing the same variable are summed.

func (*Model) AddContinuous

func (m *Model) AddContinuous(lower, upper float64) Var

AddContinuous adds a continuous variable in [lower, upper].

func (*Model) AddInteger

func (m *Model) AddInteger(lower, upper float64) Var

AddInteger adds an integer variable in [lower, upper].

func (*Model) AddVar

func (m *Model) AddVar(lower, upper float64, kind VarKind) Var

AddVar adds a variable with the given lower and upper bounds and kind, and returns a handle to it. Use Inf / NegInf for an absent upper / lower bound. For a Binary variable the bounds are ignored and set to [0, 1].

If the bounds are invalid (lower > upper, or NaN) the returned handle is still usable but Solve will fail with ErrInvalidBounds.

func (*Model) NumConstraints

func (m *Model) NumConstraints() int

NumConstraints reports the number of constraints added to the model.

func (*Model) NumVars

func (m *Model) NumVars() int

NumVars reports the number of variables added to the model.

func (*Model) SetObjective

func (m *Model) SetObjective(sense Sense, terms []Term)

SetObjective sets the optimization sense and the linear objective. Terms with the same variable are summed; variables absent from terms have coefficient zero. Calling SetObjective again replaces the objective.

func (*Model) Solve

func (m *Model) Solve() (*Solution, error)

Solve solves the model with DefaultOptions using the built-in branch-and-bound solver and returns the Solution.

Solve does not mutate the Model and is safe for concurrent use on the same Model. It returns a build error (ErrInvalidBounds, ErrUnknownVar, ErrInvalidCoef), ErrNoVariables, or ErrNoObjective when the model is malformed; ErrNumeric when the LP relaxation fails numerically; and ErrNodeLimit when branch-and-bound is truncated by Options.NodeLimit (the best incumbent, if any, is still returned in the Solution).

func (*Model) SolveWith

func (m *Model) SolveWith(opts Options) (*Solution, error)

SolveWith solves the model with the given Options.

func (*Model) SolveWithStats

func (m *Model) SolveWithStats(opts Options) (*Solution, SolveStats, error)

SolveWithStats solves the model with the given Options and additionally returns SolveStats: branch-and-bound node count, LP relaxations solved, total simplex pivots, basis factorizations, and warm-started re-solves. The stats are diagnostic only; the Solution and error are exactly as SolveWith returns.

type Options

type Options struct {
	// IntegralityTol is the tolerance for treating an LP value as integral: a
	// value within IntegralityTol of the nearest integer is considered integer.
	IntegralityTol float64
	// LPTolerance is the optimality tolerance on the reduced costs used by the LP
	// relaxation.
	LPTolerance float64
	// NodeLimit caps the number of LP relaxations solved during
	// branch-and-bound. Zero means no limit. When the limit is reached the best
	// incumbent found so far is returned together with ErrNodeLimit.
	NodeLimit int
	// MaxDepth caps the branch-and-bound tree depth. Zero means no limit.
	MaxDepth int
	// TimeLimit caps the wall-clock time of the search. Zero means no limit. When
	// the limit is reached the best incumbent found so far is returned together
	// with ErrTimeLimit, and the Solution's BestBound and Gap report the proven
	// optimality gap at that point (an anytime answer).
	TimeLimit time.Duration
	// RelativeGap stops the search once the proven relative optimality gap
	// (Objective − BestBound) / (|Objective| + 1e-10) drops to or below this
	// value. Zero means solve to proven optimality (gap 0).
	RelativeGap float64
	// AbsoluteGap stops the search once the proven absolute gap
	// |Objective − BestBound| drops to or below this value. Zero means solve to
	// proven optimality.
	AbsoluteGap float64
	// IncumbentCallback, when non-nil, is invoked each time branch-and-bound
	// finds a new best integer-feasible solution, with the objective (in the
	// model's own sense) and a copy of the variable values. It must not retain or
	// mutate the slice beyond the call and must not call back into the solver.
	IncumbentCallback func(objective float64, values []float64)
	// DisableCuts turns off root Gomory cut generation. Cuts are on by default;
	// they tighten the relaxation of integer-heavy problems (a no-op on already
	// tight models like the energy shape).
	DisableCuts bool
	// DisablePresolve turns off the root presolve pass (activity-based bound
	// propagation, singleton/duplicate/redundant rows, binary coefficient
	// tightening). Presolve is on by default and never eliminates a variable, so
	// it is transparent to the returned solution.
	DisablePresolve bool
}

Options tunes the branch-and-bound search and the LP relaxation.

func DefaultOptions

func DefaultOptions() Options

DefaultOptions returns the options used by Model.Solve.

type Relation

type Relation uint8

Relation is the comparison operator of a linear constraint.

const (
	// LessEq is <= (less than or equal).
	LessEq Relation = iota
	// GreaterEq is >= (greater than or equal).
	GreaterEq
	// EqualTo is = (equal).
	EqualTo
)

Constraint relations: the left-hand side (a linear expression) relates to the right-hand side scalar.

func (Relation) String

func (r Relation) String() string

String returns the relation symbol.

type Sense

type Sense uint8

Sense is the optimization direction of the objective.

const (
	// Minimize minimizes the objective.
	Minimize Sense = iota
	// Maximize maximizes the objective.
	Maximize
)

Optimization senses.

func (Sense) String

func (s Sense) String() string

String returns the sense name.

type Solution

type Solution struct {
	// Status is the solve outcome.
	Status Status
	// Objective is the objective value in the Model's own sense (already negated
	// back for Maximize problems). It is meaningful only when Status is Optimal.
	Objective float64
	// Values holds the optimal value of every variable, indexed by Var.ID.
	// It is nil when no feasible solution was found.
	Values []float64
	// BestBound is the best proven dual bound on the objective in the Model's own
	// sense: a lower bound for Minimize, an upper bound for Maximize. When the
	// search completes it equals Objective (proven optimal); when the search is
	// truncated (time/node limit or a requested gap) it bounds how far the
	// incumbent could still be from optimal.
	BestBound float64
	// Gap is the relative optimality gap |Objective − BestBound| / (|Objective| +
	// 1e-10). It is 0 (within rounding) when the solution is proven optimal.
	Gap float64
}

Solution is the result of solving a Model.

func (*Solution) Value

func (s *Solution) Value(v Var) float64

Value returns the solved value of v. It panics if v does not belong to the Model that produced the Solution or the Solution has no values.

type SolveStats

type SolveStats struct {
	// Nodes is the number of branch-and-bound nodes explored.
	Nodes int
	// LPSolves is the number of LP relaxations solved (root + nodes).
	LPSolves int
	// SimplexPivots is the total number of simplex pivots across all relaxations.
	SimplexPivots int
	// Refactorizations is the total number of basis factorizations.
	Refactorizations int
	// WarmStarts is the number of relaxations started from a parent basis.
	WarmStarts int
	// Cuts is the number of Gomory cutting planes added at the root.
	Cuts int
	// CutRounds is the number of root cut-generation rounds performed.
	CutRounds int
}

SolveStats reports diagnostics from the native LP/branch-and-bound core.

type Solver

type Solver interface {
	Solve(m *Model) (*Solution, error)
}

Solver solves a Model and returns its Solution. NewSolver returns the built-in branch-and-bound solver; the interface allows callers to substitute their own.

func NewSolver

func NewSolver(opts Options) Solver

NewSolver returns the built-in branch-and-bound solver configured with opts.

type Status

type Status uint8

Status is the outcome of a solve.

const (
	// Optimal indicates a proven optimal solution was found.
	Optimal Status = iota
	// Infeasible indicates the problem has no feasible solution.
	Infeasible
	// Unbounded indicates the objective is unbounded on the feasible region.
	Unbounded
)

Solve outcomes.

func (Status) String

func (s Status) String() string

String returns the status name.

type Term

type Term struct {
	Var  Var
	Coef float64
}

Term pairs a variable with a scalar coefficient in a linear expression.

type Var

type Var struct {
	// contains filtered or unexported fields
}

Var is an opaque handle to a decision variable within one Model. The zero Var is not valid; obtain one from a Model's AddVar family.

func (Var) ID

func (v Var) ID() int

ID returns the variable's zero-based index within its Model. It equals the position of the variable's value in Solution.Values.

type VarKind

type VarKind uint8

VarKind is the integrality class of a decision variable.

const (
	// Continuous is a real-valued variable within its bounds.
	Continuous VarKind = iota
	// Integer is a variable constrained to integer values within its bounds.
	Integer
	// Binary is an integer variable constrained to {0, 1}.
	Binary
)

Variable kinds.

func (VarKind) String

func (k VarKind) String() string

String returns the kind name.