///|
/// Public entry points of `moonopt`.
///
/// This package is orchestration only: it validates a `@model.Model`, reduces it
/// with `@presolve` unless the caller turns that off, adapts it into the form the
/// kernel accepts, calls the kernel, and reports the outcome. Algorithm details
/// live in the kernel packages — see `docs/design.md`.
///
/// Current capability: continuous models with `x >= 0`, finite upper bounds and any
/// mix of `≤`, `≥`, `=` rows, solved by the sparse revised simplex with the presolve
/// reduction on by default; and models with integer variables, handed to the branch
/// and bound package (`mip`) which proves every relaxation it branches on. Anything
/// outside that class returns `NotSolved` together with the reason — never a
/// plausible but unverified answer.
///
/// A model with integer variables is the one class the presolve reduction is not
/// offered. The kernel answers a relaxation when it is asked for one, but this entry
/// point never reports a relaxation as if it were an answer to an integer model: the
/// model goes to the branch and bound, whose statuses map onto `Optimal`,
/// `Infeasible`, and — when the node budget runs out with work left — `NodeLimit`,
/// where `values` is a feasible integer point that is *not* proven optimal and
/// `bound`/`gap` say how much room is left; an unbounded relaxation whose point and
/// ray are integral on every integer variable maps onto `Unbounded`, and any other
/// unbounded relaxation stays `NotSolved`.
///
/// With presolve on — the default — the model is reduced first and the kernel's
/// solution is mapped back to the original variables. The way back is *verified*,
/// not assumed: a reconstructed solution is reported as `Optimal` only after the
/// original model agrees that it satisfies its rows and bounds and that the
/// objective, recomputed from the original cost vector, is the number being
/// reported. A reduction that fails that check becomes `NotSolved` with the
/// measured violations, because a wrong answer carrying a confident status is the
/// one outcome this package exists to prevent.
///|
/// Row and bound violation a reconstructed solution may carry before it is
/// refused. The same scale the kernel uses for its own feasibility decision.
let reconstruction_tolerance : Double = 1.0e-7
///|
/// Outcome of a solve call.
pub(all) enum SolveStatus {
/// The answer is proven: for a linear model the kernel's own residual and dual checks
/// stand behind it, and for a model with integer variables every relaxation the
/// branch and bound branched on was checked by `verify` and the tree was exhausted (or
/// every open bound was closed against the incumbent).
Optimal
/// No feasible point exists, and that is proven — a Farkas ray for a linear model, an
/// exhausted tree with no integer point for one with integer variables.
Infeasible
/// The objective has no bound on the feasible set. For a linear model that is the
/// kernel's own verdict; for an integer model it is reported only when the relaxation's
/// point and improving ray are integral on every integer variable — otherwise the
/// relaxation is `NotSolved`, because an improving ray need not carry any lattice
/// points.
Unbounded
/// The search stopped at the node budget with work left. `values` is a feasible
/// integer point that is **not** proven optimal, and `bound` and `gap` are the two
/// numbers a caller decides on: the best objective still reachable among the nodes
/// left open, and how far the incumbent can still be beaten.
NodeLimit
/// Nothing was proven, and `message` says why: the iteration or node budget ran out
/// without an answer, the model is outside the class this solver accepts, an
/// unbounded relaxation says nothing about the integer model, or a certificate was
/// refused by the independent checker.
NotSolved
} derive(Debug, Eq)
///|
/// What a solve call produced.
///
/// `values` is indexed by variable index. It is empty unless the status is `Optimal` or
/// `NodeLimit`; for `NodeLimit` it is a feasible integer point that is not proven
/// optimal. `message` carries the reason otherwise, and also the kernel's or the
/// search's own note when a run needed its recovery attempt or reached no verdict on a
/// relaxation.
pub struct Solution {
status : SolveStatus
values : Array[Double]
/// The objective at `values`, in the model's own sense. Meaningful only when `values`
/// is not empty: a run that never stood on a feasible point reports zero here, and it
/// is `values` — not this number — that says whether there is a point at all.
objective : Double
/// Pivots the kernel ran. A branch and bound run counts relaxations instead, in
/// `nodes`, because that is what its budget is spent on.
iterations : Int32
/// Relaxations a branch and bound run solved, including the ones its primal heuristic
/// spent walking a fractional node down to a whole point. Zero for a linear solve.
nodes : Int32
/// The best objective still reachable among the work left open, in the model's sense.
/// Equal to `objective` when nothing is left to prove — a linear solve is its own
/// bound, so it is equal there too — and the number `gap` is measured against
/// otherwise.
bound : Double
/// How far the incumbent can still be beaten: the difference between `objective` and
/// `bound` in the model's own scale, so zero means proven optimal. Carried separately
/// because it is the number a caller decides on, and recomputing it from two numbers
/// whose sense the caller has to know is how that decision goes wrong.
gap : Double
message : String
} derive(Debug)
///|
/// Tunables for a solve call.
///
/// `pub(all)` so callers can write `{ max_iterations: 1000, presolve: false }`
/// directly; `SolveOptions::new` is the readable alternative.
pub(all) struct SolveOptions {
max_iterations : Int32
/// Reduce the model before the kernel sees it and reconstruct the solution
/// afterwards. On by default: every reduction is proved before it fires and every
/// reconstruction is checked against the original model. Turn it off to observe
/// the kernel's own behaviour on the model exactly as written.
///
/// A model with integer variables skips the reduction either way: see the package
/// doc for why that promise must not depend on what a reduction happened to fix.
presolve : Bool
/// Relaxations a model with integer variables may solve, including the ones the primal
/// heuristic spends. Reaching it ends the run as `NodeLimit` with the incumbent and the
/// bound still open. Ignored by a linear solve, which spends `max_iterations` instead.
max_nodes : Int32
/// Rounds of cuts the root relaxation may be tightened with before the tree is built, or
/// `0` to search the model as written. Each round re-solves the root, and every cut that
/// lands has been re-derived and checked by the independent verifier first — an
/// unchecked cut would remove integer points from every node below it without anything
/// downstream noticing. Ignored by a linear solve, which has nothing to round.
cut_rounds : Int32
} derive(Debug)
///|
/// Default solve options: a generous iteration budget and node budget, presolve on, cuts
/// on.
pub fn SolveOptions::default() -> SolveOptions {
{ max_iterations: 10000, presolve: true, max_nodes: 20000, cut_rounds: 2, }
}
///|
/// Solve options with an explicit iteration cap, keeping presolve on.
pub fn SolveOptions::new(max_iterations : Int32) -> SolveOptions {
{ max_iterations, presolve: true, max_nodes: 20000, cut_rounds: 2, }
}
///|
/// Solve options with an explicit node budget, for a model with integer variables.
pub fn SolveOptions::with_max_nodes(max_nodes : Int32) -> SolveOptions {
{ max_iterations: 10000, presolve: true, max_nodes, cut_rounds: 2, }
}
///|
/// Maps a kernel result onto the public one.
///
/// A question the kernel cannot answer comes back as `NotSolved` with the kernel's
/// own words; an optimal run keeps its message too, because that is where a run
/// says it had to restart with Bland's rule to get there.
///
/// A linear solve is its own bound: an optimal one has nothing left open, so `bound` is
/// the objective and `gap` is zero, which is what makes the two fields mean the same
/// thing on every path that fills them in.
fn from_simplex(result : @simplex.SimplexResult) -> Solution {
match result.status {
@simplex.SimplexStatus::Optimal =>
{
status: Optimal,
values: result.values,
objective: result.objective,
iterations: result.iterations,
nodes: 0,
bound: result.objective,
gap: 0.0,
message: result.message,
}
@simplex.SimplexStatus::Infeasible =>
{
status: Infeasible,
values: [],
objective: 0.0,
iterations: result.iterations,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
@simplex.SimplexStatus::Unbounded =>
{
status: Unbounded,
values: [],
objective: 0.0,
iterations: result.iterations,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
@simplex.SimplexStatus::IterationLimit =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: result.iterations,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
@simplex.SimplexStatus::NumericalFailure =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: result.iterations,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
@simplex.SimplexStatus::TooLarge =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
}
}
///|
/// Maps a branch and bound result onto the public one.
///
/// Only the search's own `Optimal` becomes `Optimal` here, and `NodeLimit` becomes
/// `NodeLimit` rather than `NotSolved`: a budget that ran out with an incumbent and an
/// open bound is a state a caller can act on, and flattening it into "nothing was
/// proven" would throw away the two numbers that say how much is left.
///
/// The three statuses that are *not* an answer to the question all come back as
/// `NotSolved`, each with the search's own words: an unbounded relaxation whose ray
/// carries no integrality (an improving ray of the relaxation need not carry lattice
/// points), a certificate the independent checker refused (the kernel is wrong
/// somewhere, which is not a statement about this model), and a model that failed
/// validation. None of them may be reported as `Unbounded`, `Optimal` or `Infeasible`,
/// because none of them proves any of those — only the search's own `Unbounded`, an
/// integral point and ray, maps onto the public `Unbounded`.
fn from_mip(result : @mip.MipResult) -> Solution {
match result.status {
@mip.MipStatus::Optimal =>
{
status: Optimal,
values: result.values,
objective: result.objective,
iterations: 0,
nodes: result.nodes,
bound: result.bound,
gap: result.gap,
message: result.message,
}
@mip.MipStatus::Infeasible =>
{
status: Infeasible,
values: [],
objective: 0.0,
iterations: 0,
nodes: result.nodes,
bound: result.bound,
gap: result.gap,
message: result.message,
}
@mip.MipStatus::NodeLimit =>
{
status: NodeLimit,
values: result.values,
objective: result.objective,
iterations: 0,
nodes: result.nodes,
bound: result.bound,
gap: result.gap,
message: result.message,
}
@mip.MipStatus::UnboundedRelaxation =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: result.nodes,
bound: result.bound,
gap: result.gap,
message: result.message,
}
@mip.MipStatus::Unbounded =>
{
// The search proved the integer model unbounded with an integral point and an
// integral ray (both carried on `result`), so this is the same verdict as a linear
// solve: `Unbounded`, with no values in the public result — the point and ray are
// the proof, and they live on `MipResult` the way the kernel's do on `SimplexResult`.
status: Unbounded,
values: [],
objective: 0.0,
iterations: 0,
nodes: result.nodes,
bound: 0.0,
gap: 0.0,
message: result.message,
}
@mip.MipStatus::Unverified =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: result.nodes,
bound: result.bound,
gap: result.gap,
message: result.message,
}
@mip.MipStatus::Invalid =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: result.message,
}
}
}
///|
/// Solves `model` with default options.
pub fn solve(model : @model.Model) -> Solution {
solve_with(model, SolveOptions::default())
}
///|
/// The same model with every integrality flag cleared — its **LP relaxation**.
///
/// This package answers a model with integer variables with the branch and bound (`mip`), on
/// purpose: a relaxation is not an answer to an integer model. A caller who wants the *bound*
/// instead — the number a plan they already have is measured against — asks for that
/// explicitly, and this is the question: the same variables in the same order, the same rows,
/// the same objective, no integrality. What comes back is an ordinary linear model, so it goes
/// through `Model::solve` like any other, reduction and independent verification included, and
/// `bound`/`gap` mean what they mean there (a linear solve is its own bound).
///
/// The variables keep their indices, so a point the caller already holds is still a point of
/// the same space and can be read back on the relaxation: compare its `objective_value` with
/// the relaxation's optimum to get the gap the plan is worth, and check
/// `max_row_violation`/`max_bound_violation` on it if the plan is not known to be feasible.
/// `examples/schedule_bound` is that whole workflow in one file.
pub fn relaxation(model : @model.Model) -> @model.Model {
let relaxed = @model.Model::new(model.sense())
for j = 0; j < model.num_vars(); j = j + 1 {
let variable = model.get_var(j)
relaxed.add_var(variable.name, lb=variable.lb, ub=variable.ub, is_int=false)
|> ignore
}
relaxed.set_objective(@model.normalize_terms(model.objective()))
for i = 0; i < model.num_constraints(); i = i + 1 {
let constraint = model.constraint(i)
relaxed.add_named_constraint(
@model.normalize_terms(constraint.terms),
constraint.rel,
constraint.rhs,
constraint.name,
)
}
relaxed
}
///|
/// `true` when the model has at least one integer variable.
fn has_integer_variables(model : @model.Model) -> Bool {
for j = 0; j < model.num_vars(); j = j + 1 {
if model.get_var(j).is_int {
return true
}
}
false
}
///|
/// Solves `model` and returns a `Solution` describing what happened.
///
/// The work happens in the `simplex` kernel; this function validates the model,
/// reduces it, maps the kernel's status onto the public one, and makes sure that a
/// question the kernel cannot answer comes back as `NotSolved` with a reason rather
/// than as a plausible looking answer.
///
/// The reduction is never trusted on its own account. When presolve proves the
/// model infeasible or unbounded, that verdict is returned without a kernel run;
/// when it returns a reduced model, the kernel's solution is reconstructed and then
/// measured against the **original** rows, bounds and objective, and a
/// reconstruction that does not verify turns into `NotSolved` rather than an
/// answer.
///
/// Integer models are the one class the reduction is not offered: the reduction must not
/// be the reason an integer model gets an answer, because that answer would depend on
/// whether the reduction happened to fix every integer variable — and that is exactly
/// where a fractional fixed value would slip through. Such a model goes straight to the
/// branch and bound, on the model as written.
pub fn solve_with(model : @model.Model, options : SolveOptions) -> Solution {
let issues = model.validate()
if issues.length() > 0 {
return {
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: "model is not well formed: " + issues[0],
}
}
if has_integer_variables(model) {
// Built from the branch and bound's own defaults so this layer cannot drift from
// them: only the two budgets the caller owns are overridden.
let base = @mip.MipOptions::default()
let mip_options : @mip.MipOptions = {
max_nodes: options.max_nodes,
max_iterations: options.max_iterations,
tolerance: base.tolerance,
verify_nodes: base.verify_nodes,
dive_steps: base.dive_steps,
dive_attempts: base.dive_attempts,
cut_rounds: options.cut_rounds,
primal_rounding: base.primal_rounding,
neighbourhood_passes: base.neighbourhood_passes,
row_strengthening: base.row_strengthening,
column_dominance: base.column_dominance,
node_solver: base.node_solver,
cut_hook: base.cut_hook,
}
return from_mip(@mip.solve_mip(model, mip_options))
}
let kernel_options = @simplex.SimplexOptions::with_max_iterations(
options.max_iterations,
)
if !options.presolve {
return match @simplex.solve_model_with(model, kernel_options) {
Err(reason) =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: reason,
}
Ok(result) => from_simplex(result)
}
}
let presolved = @presolve.presolve(model)
match presolved.verdict() {
@presolve.PresolveVerdict::Infeasible =>
{
status: Infeasible,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: presolved.reason(),
}
@presolve.PresolveVerdict::Unbounded =>
{
status: Unbounded,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: presolved.reason(),
}
@presolve.PresolveVerdict::Reduced => {
let reduced = presolved.reduced()
if reduced.num_vars() == 0 {
// Every variable was decided by the reduction, so its own answer is the
// whole answer — and it still goes through the same reconstruction check
// as any other solution before it leaves this function.
let values = presolved.reconstruct([])
let objective = presolved.objective(0.0)
let rows = @presolve.max_row_violation(model, values)
let bounds = @presolve.max_bound_violation(model, values)
let recomputed = @presolve.objective_value(model, values)
if rows > reconstruction_tolerance ||
bounds > reconstruction_tolerance ||
!@core.approx_eq(objective, recomputed, rel=1.0e-9, abs=1.0e-6) {
return {
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: "the reduction decided every variable but its solution did not verify against the original model (worst row violation " +
rows.to_string() +
", worst bound violation " +
bounds.to_string() +
", objective " +
objective.to_string() +
" vs " +
recomputed.to_string() +
" recomputed from the original model)",
}
}
{
status: Optimal,
values,
objective,
iterations: 0,
nodes: 0,
bound: objective,
gap: 0.0,
message: "",
}
} else {
match @simplex.solve_model_with(reduced, kernel_options) {
Err(reason) =>
{
status: NotSolved,
values: [],
objective: 0.0,
iterations: 0,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: reason,
}
Ok(result) =>
match result.status {
@simplex.SimplexStatus::Optimal => {
let values = presolved.reconstruct(result.values)
let objective = presolved.objective(result.objective)
let rows = @presolve.max_row_violation(model, values)
let bounds = @presolve.max_bound_violation(model, values)
let recomputed = @presolve.objective_value(model, values)
if rows > reconstruction_tolerance ||
bounds > reconstruction_tolerance ||
!@core.approx_eq(
objective,
recomputed,
rel=1.0e-9,
abs=1.0e-6,
) {
return {
status: NotSolved,
values: [],
objective: 0.0,
iterations: result.iterations,
nodes: 0,
bound: 0.0,
gap: 0.0,
message: "the reduced model solved but its reconstruction did not verify against the original model (worst row violation " +
rows.to_string() +
", worst bound violation " +
bounds.to_string() +
", objective " +
objective.to_string() +
" vs " +
recomputed.to_string() +
" recomputed from the original model)",
}
}
{
status: Optimal,
values,
objective,
iterations: result.iterations,
nodes: 0,
bound: objective,
gap: 0.0,
message: result.message,
}
}
_ => from_simplex(result)
}
}
}
}
}
}