Hephaestus.Optimisation 0.1.0-beta.1

⚒️ Hephaestus

A purely functional modelling layer for MILP and SMT solvers, for C# 14 / .NET 10.

💡 Overview

"Nothing at all takes place in the universe in which some rule of maximum or minimum does not appear." — Leonhard Euler.

Optimisation models are specified as mathematics. Hephaestus lets the code be that mathematics, so that reviewing a constraint is a direct comparison with its specification:

occ_A ∧ occ_B  ⟹  (dep_A + h ≤ dep_B) ∨ (dep_B + h ≤ dep_A)

var isSeparated    = (departureA + headway <= departureB) | (departureB + headway <= departureA);
var isConflictFree = !(occupiesA & occupiesB) | isSeparated;

There are no auxiliary booleans to declare, no gadget factories, and no big-M to work out by hand: the boolean structure is lowered to linear constraints automatically, and every big-M is derived from the bounds your constraints already state. The solver underneath is swappable.

⬇️ Installation

dotnet add package Hephaestus.Optimisation --prerelease
dotnet add package Hephaestus.Optimisation.OrTools --prerelease

Hephaestus is in beta, so --prerelease is needed for now, and the API may change between versions.

(The bare name Hephaestus was already taken on nuget.org, hence Hephaestus.Optimisation. The namespace is simply Hephaestus.) Operators are C# 14 extension members, so consumers need the .NET 10 SDK.

📖 A complete example

using Hephaestus;
using Hephaestus.OrTools;

var departureA = Variable.Continuous("departureA");
var departureB = Variable.Continuous("departureB");
var occupiesA  = Variable.Binary("occupiesA");
var occupiesB  = Variable.Binary("occupiesB");
const double headway = 120;

var withinTheHour  = departureA.Between(0, 3600) & departureB.Between(0, 3600);
var isSeparated    = (departureA + headway <= departureB) | (departureB + headway <= departureA);
var isConflictFree = !(occupiesA & occupiesB) | isSeparated;

var problem = Problem.Minimise(
    departureA + departureB,
    subjectTo: withinTheHour & isConflictFree & occupiesA & occupiesB);

var summary = OrToolsSolver.Create().Solve(problem).Match(
    optimal:    solution => $"A leaves at {solution.Value(departureA)}, B at {solution.Value(departureB)}",
    feasible:   solution => $"best found: {solution.ObjectiveValue}",
    infeasible: () => "no timetable exists",
    unbounded:  () => "unbounded",
    unknown:    reason => $"the solver gave up: {reason}");

Swap OrToolsSolver.Create() for OrToolsSolver.Create(OrToolsSolverId.Highs) or new Z3Solver() and nothing else changes. samples/TrainHeadway is a slightly larger version in NodaTime types.

🧠 The ideas

Expressions are data; operators only build it

There are two algebraic data types, each an interface with a handful of sealed records:

• ILinearExpression: Constant, Sum, Product, and the variables.
• IBooleanExpression: BooleanConstant, Comparison, Negation, Conjunction, Disjunction, Implication, Equivalence, and BinaryVariable.

Keeping them apart makes illegal compositions unrepresentable: x * y, (x <= 1) + 1 and if (x <= y) do not compile. A BinaryVariable belongs to both types, so occupiesA + occupiesB <= 1 and occupiesA & occupiesB are both fine.

Operators (+ - * /, <= >= < >, & | ! ^, plus EqualTo, NotEqualTo, Between, Implies, Iff) are extension members that do nothing but construct records: a + b is new Sum(a, b). Plain values mix in on either side: x + 5, 5 + x, 0 <= x, and, for constraints that depend on known data, train.IsFreight.Implies(departure >= curfew) or isPeak & (headway >= 180). Nothing is flattened or simplified at construction time. All interpretation happens later, in separate passes over the data:

Because expressions are immutable values with structural equality, a constraint can be composed from smaller named pieces, logged, compared and unit-tested without a solver or a "problem" object in sight. (== is reserved by C# for that structural equality, which is why the equality constraint is spelt EqualTo.)

A problem is one constraint, not a list

Problem.Satisfy(constraint)
Problem.Minimise(objective, subjectTo: constraint)
Problem.Maximise(objective, subjectTo: constraint)

A list of constraints would only duplicate what & already means. Conjoin a collection with constraints.AllOf() (and AnyOf(), and terms.Sum()); these build balanced trees, so a hundred thousand constraints are no deeper than seventeen levels. Problems are records too: "rollback" is keeping the old value, and extending is problem with { Constraint = problem.Constraint & extra }.

Variables have no bounds; bounds are constraints

A variable is a name and a kind (Continuous, Integer, Binary). 0 <= x & x <= 3600 is a constraint like any other. The encoder recognises unconditional single-variable constraints and hands them to the solver as native column bounds, so this costs nothing.

Big-M is derived, never supplied

The MILP encoding runs in four pure steps:

1. Normalise to negation normal form, with comparisons as affine form <= 0 or == 0.
2. Encode the logic as guarded rows, "if these literals all hold then e <= 0", introducing as few auxiliary binaries as possible. In a disjunction the existing literals become guards, and only the compound disjuncts other than the last need an auxiliary. (a + h <= b) | (b + h <= a) comes out as the textbook either-or with one binary; flag.Iff(x >= 5) needs none. Equal subformulas share one auxiliary.
3. Propagate bounds (feasibility-based bound tightening) over the unconditional rows. x <= 3600 bounds x directly; arrival == departure + run bounds arrival once departure and run are bounded; and so on down the line until nothing changes.
4. Relax each guard with e <= M · (number of guards that are off), where M is the largest value e can take inside the propagated bounds. That is the tightest valid M, computed separately for every row. Rows that can never bind are dropped; rows that can never hold forbid their guards outright.

If some M is infinite, encoding fails with an error naming the variables that lack bounds. That is deliberate: a guessed big-M that is too small silently cuts off solutions, and the whole point is that wrong constraints should not fail silently. If you really want a guess, opt in with new EncodingOptions(FallbackBigM: 1e6).

The first two steps are problem.EncodeLogic(), giving an IndicatorProblem; the last two are .RelaxGuards(), giving a MilpProblem; problem.Encode() is both. Either result can be printed with .Format(). Solvers with indicator or half-reified constraints (Gurobi, CP-SAT) stop after the first half and never see a big-M, and the Z3 backend skips all of it, because an SMT solver takes the boolean structure as it stands.

Strictness

Over whole-valued expressions, n < 5 is exactly n <= 4. Over the reals a MILP cannot express strictness, so x < 5 becomes x + ε <= 5 with EncodingOptions.StrictnessEpsilon (default 1e-4, comfortably above solver tolerances). Negation is where this usually bites: !(x <= 5) is x > 5. Z3 has true strict inequalities and no epsilon.

Reified truth values

To use the truth of a constraint as a number (say, to count violated soft constraints), tie it to a binary variable yourself: isLate.Iff(arrival >= deadline), then use isLate in the objective. Iff binds in both directions.

Typed expressions

A Quantity<T> is a linear expression read as an amount of type T (a duration); a Point<T, TDelta> is one read as a position (a date-time) whose differences are amounts of TDelta. Each pairs an ordinary ILinearExpression with an IProjection<T>, an affine map between T and the solver's number line ("seconds since 08:00").

The two types carry the right algebra, once, generically:

Quantity<Duration>               runTime  = arrival - departure;          // point - point
Point<LocalDateTime, Duration>   earliest = departure + dwell + minimum;  // point + quantity
IBooleanExpression               onTime   = arrival <= deadline;          // compare with plain values
Point<LocalDateTime, Duration>   release  = start + dwell;                // plain value + quantity
//                                          arrival + departure           // does not compile
//                                          2 * arrival                   // does not compile
//                                          arrival <= dwell              // does not compile

Plain values of the right type are welcome wherever an expression is, on either side (Duration.FromMinutes(2) + departure, start <= departure, dwell.Between(minimum, arrival - departure)), and delay.In(Duration.FromMinutes(1)) turns a quantity back into a plain linear expression ("delay in minutes") for a cost function.

Solutions are read back in the same types: solution.Value(arrival) is a LocalDateTime, and so is solution.Value(departure + dwell); any expression can be read, not just variables. Expressions under different projections (minutes against seconds, different origins) are reconciled automatically.

The core ships projections for TimeSpan, DateTime and DateTimeOffset; Hephaestus.Optimisation.NodaTime adds Duration, Instant, LocalDateTime, LocalDate (in Periods of whole days), LocalTime, OffsetDateTime and ZonedDateTime:

var departure = Variable.LocalDateTime("departure", origin: start);
var dwell     = Variable.Duration("dwell", unit: Duration.FromSeconds(30), inWholeUnits: true);  // quantised

Supporting another type means writing one small record that implements IProjection<T> (or IPointProjection<T, TDelta>), plus, for a point type, the three one-line T ± Quantity<TDelta> operators that C# will not let the core declare generically (they delegate to quantity.Beyond(origin, projection)). Typed reads round the underlying number to five decimal places of the unit by default, so that solver noise does not turn 08:04:00 into 08:03:59.99999999.

Swapping the solver

ISolver is the seam: ISolveResult Solve(IProblem problem, CancellationToken cancellationToken = default). A backend joins at whichever level suits the solver:

IndicatorSolver and MilpSolver wrap the latter two into an ISolver: they encode, solve, hide the auxiliaries and snap whole-number variables. A backend is about a hundred lines.

GurobiSolver.Create()                       // indicator constraints, no big-M
GurobiSolver.Create(useIndicators: false)   // classic big-M, derived per row
CpSatSolver.Create()                        // whole-number problems only
HighsSolver.Create()
OrToolsSolver.Create(OrToolsSolverId.Scip)
new Z3Solver()

The test suite runs one contract against SCIP, CBC, HiGHS (standalone and through OR-Tools), Gurobi (both formulations) and Z3, and a second, whole-number contract against those and CP-SAT. The Gurobi tests are skipped where no licence is found.

⚠️ Things worth knowing

• Sum and conjoin collections with Sum(), AllOf() and AnyOf() rather than folding + or & yourself. (Folding still works: deep trees are handled without overflowing the stack. It is merely slower, and the records' built-in ToString/Equals do recurse.)
• An equivalence duplicates its operands when negations are pushed inwards, so deeply nested Iffs grow exponentially. Name the inner ones with binary variables.
• Variables are identified by name and kind. Two variables of different kinds sharing a name is an error.
• Do not reference Hephaestus.Optimisation.Highs and Hephaestus.Optimisation.OrTools from the same application. Both upstream packages ship a native highs.dll, of different versions, and whichever is copied last breaks the other.
• CP-SAT works in whole numbers over finite domains. It refuses continuous variables by name, needs every variable bounded (directly or by implication), and scales fractional coefficients by a power of ten, refusing those that none makes whole (a third, say).
• HiGHS, as driven by OR-Tools, prints a one-line banner to standard output on every solve. That is upstream behaviour; the standalone HiGHS backend and the other solvers are silent.

🛠️ Building

dotnet build
dotnet test
dotnet run --project samples/TrainHeadway
dotnet pack --configuration Release --output artefacts

The test projects use xunit.v3 on Microsoft.Testing.Platform (see global.json). tests/Hephaestus.Tests/EncodingEquivalenceTests.cs is the one to read first: it checks, exhaustively over a grid and for hundreds of random formulas, that an assignment satisfies a formula exactly when it extends to a solution of the encoded rows.

Publishing a GitHub release tagged v1.2.3 tests, packs and publishes version 1.2.3 of every package to nuget.org, and publishes the API documentation. See CONTRIBUTING.md.

👮‍♂️ Licence

MIT; see LICENSE.md.

👁️ See also

• This library was made using C# Library Template.