Bases and warm starts¶

A linear program solved by the simplex method ends on a basis: as many basic entities as there are rows, variables or rows, and every other one nonbasic, at one of its bounds as a rule. MIP++ reads that basis as a snapshot that you query with the handles of your own variables and constraints, and writes a basis back so that the next solve starts from it: the basis of a model solved before, of a twin model, or one of your own making.

Reading a basis¶

A workshop makes chairs, tables and desks from 40 boards and 18 hours of labour, plus up to 4 hours of overtime at 3 each, and has orders for at least 3 tables. This program finds the most profitable plan and prints its basis:

#include <concepts>
#include <iostream>
#include <string_view>
#include <variant>

#include "mippp/solvers/highs/all.hpp"

using namespace mippp;
using namespace mippp::operators;

// The most specific tag a status names, by derivation, so that it compiles
// on every backend's status variant.
template <typename Status>
std::string_view describe(const Status & status) {
    return std::visit(
        []<typename Tag>(const Tag &) -> std::string_view {
            if constexpr(std::derived_from<Tag, basis_status::basic>)
                return "basic";
            else if constexpr(std::derived_from<Tag,
                                                basis_status::nonbasic_fixed>)
                return "nonbasic, fixed";
            else if constexpr(std::derived_from<
                                  Tag, basis_status::nonbasic_at_lower_bound>)
                return "nonbasic, at its lower side";
            else if constexpr(std::derived_from<
                                  Tag, basis_status::nonbasic_at_upper_bound>)
                return "nonbasic, at its upper side";
            else if constexpr(std::derived_from<Tag, basis_status::superbasic>)
                return "nonbasic, at no bound";
            else
                return "nonbasic";
        },
        status);
}

int main() {
    highs_lp model;
    auto chairs = model.add_variable();
    auto tables = model.add_variable();
    auto desks = model.add_variable();
    auto overtime = model.add_variable({.lower_bound = 0, .upper_bound = 4});
    model.set_maximization();
    model.set_objective(4 * chairs + 7 * tables + desks - 3 * overtime);
    auto wood = model.add_constraint(2 * chairs + 4 * tables + 2 * desks <= 40);
    auto labour =
        model.add_constraint(chairs + 2 * tables + desks <= 18 + overtime);
    auto orders = model.add_constraint(tables >= 3);
    model.solve();

    auto basis = model.get_basis();
    std::cout << "chairs:   " << describe(basis.get_status(chairs)) << '\n';
    std::cout << "tables:   " << describe(basis.get_status(tables)) << '\n';
    std::cout << "desks:    " << describe(basis.get_status(desks)) << '\n';
    std::cout << "overtime: " << describe(basis.get_status(overtime)) << '\n';
    std::cout << "wood:     " << describe(basis.get_status(wood)) << '\n';
    std::cout << "labour:   " << describe(basis.get_status(labour)) << '\n';
    std::cout << "orders:   " << describe(basis.get_status(orders)) << '\n';
    return 0;
}

It prints:

chairs:   basic
tables:   basic
desks:    nonbasic, at its lower side
overtime: basic
wood:     nonbasic, at its upper side
labour:   nonbasic, at its upper side
orders:   nonbasic, at its lower side

The plan makes 14 chairs and the 3 tables ordered, uses 2 hours of overtime and no desk. The three variables inside their bounds are basic, one per row; the wood and labour rows are used up, at their upper side, and the orders row is at its lower side.

get_basis() returns the snapshot by value: it describes the basis of the last solve, and a later change of the model leaves it as it was. Its get_status(v) and get_status(c) return a std::variant over the tags of namespace basis_status, which form a hierarchy:

The side of a row is the side of its activity, the value of its left-hand side: wood, whose 2 chairs + 4 tables + 2 desks reaches 40, reads nonbasic_at_upper_bound on every backend. Several solvers name the side of a slack variable instead, which is the other side on some rows, or no side at all, and MIP++ translates their answer.

Each backend reports the most specific tag its solver's answer gives. A fixed entity, a variable with equal bounds or an == row, reads nonbasic_fixed where the solver has a code for it, or, for an == row of gurobi_lp and cplex_lp, where the row's sense says so, and otherwise the side its solver names, which follows the sign of its reduced cost or dual; both derive from nonbasic_at_bound. A free nonbasic variable reads nonbasic_free where the solver tells it apart, and superbasic where one code covers it and the entities strictly inside their bounds, which a run stopped before the optimum leaves, and on mosek_lp the certificate of an unbounded LP. Support by model lists which backend reports what.

So test a status at the level your code cares about. is<S>(s) compiles only for a tag the status's variant lists, and is_a<S>(s) for a tag it lists or that one of its tags derives from. is_a on basic, nonbasic, nonbasic_at_bound, nonbasic_at_lower_bound and nonbasic_at_upper_bound compiles on every status of every backend, and is_a<basis_status::nonbasic_at_bound>(s) holds for every entity at a bound, fixed ones included. is_a<basis_status::superbasic> compiles on every column status, but not on the row statuses of gurobi_lp and cplex_lp, whose rows are never at no bound. is on basic, nonbasic, nonbasic_at_lower_bound and nonbasic_at_upper_bound compiles everywhere, while is and is_a on nonbasic_fixed or nonbasic_free, and is on superbasic, compile only where the status lists them. No tag of the library refines basic, so on every status of every backend is<basis_status::basic>(s) holds for every basic entity and !is<basis_status::basic>(s) for every nonbasic one, whereas is<basis_status::nonbasic>(s) holds for bare nonbasic only; is_a<basis_status::basic>(s) differs from is only on a status of your own that lists a tag refining basic. describe, above, visits the tag and tests it by derivation, so it compiles on every status of every backend: generic code over rows and columns alike does the same. Code that must know whether an entity is fixed compares its bounds. The types are model_basis_t<M> for the snapshot, and model_variable_basis_status_t<M> and model_constraint_basis_status_t<M> for the two statuses, which may differ.

Starting a solve from a basis¶

A new product comes up: a stool, made from one board and one hour, sold at 5. Column generation adds it to the solved model as a column, and starts the next solve from the basis read before:

auto stools = model.add_column({{wood, 1.}, {labour, 1.}},
                               {.obj_coef = 5, .lower_bound = 0});
// stools did not exist when basis was taken: it reads nonbasic
std::cout << "stools:   " << describe(basis.get_status(stools)) << '\n';
model.set_basis(basis);
model.solve();
std::cout << "profit:   " << model.get_solution_value() << '\n';

It prints:

stools:   nonbasic
profit:   89

set_basis(b) hands the basis to the solver, and the next solve() starts from it. A column the snapshot does not cover, created after it was taken, reads bare nonbasic, and a row it does not cover reads basic, which is how every solver extends its own basis: the basis stays square, one basic entity per row, and the solve starts one column away from the optimum. A re-solve of the model you changed in place usually warm-starts without set_basis, since the solver keeps its own basis; set_basis matters when it does not, see Support by model, and to start from a basis the solver no longer holds:

• Another model. A snapshot is keyed by handle, so it seeds a second model built in the same order, of the same backend or another one: other.set_basis(model.get_basis()). Models are move-only, so this is how a second model, built per thread or per node, starts where the first ended. Removals break the correspondence: a target that did not remove the same variables needs a basis of its own making, mapping its handles to the source's.
• An earlier basis. After a dive that changed bounds, set_basis(saved) takes the solver back to where it was.
• A basis of your own making. set_basis accepts any type whose const member functions get_status(v) and get_status(c) return a std::variant holding a basis_status::basic tag and a basis_status::nonbasic tag. A default-constructed model_basis_t<M> is the slack basis, every column nonbasic and every row basic, and set_status(e, tag) edits a snapshot, growing it for handles it does not cover: this is how you warm-start a model never solved, here the workshop's, whose class is Model:

model_basis_t<Model> start;  // the slack basis
start.set_status(chairs, basis_status::basic{});
start.set_status(tables, basis_status::basic{});
start.set_status(wood, basis_status::nonbasic_at_upper_bound{});
start.set_status(orders, basis_status::nonbasic_at_lower_bound{});
model.set_basis(start);
model.solve();

set_status stores the most specific tag the snapshot's status lists: nonbasic_fixed written on a backend that reports no fixed entity is stored as bare nonbasic, which the solver resolves against the bounds. set_basis writes any basis as the model's own snapshot would hold it: a tag of your own as the tag it refines, a superbasic that another backend read as bare nonbasic where the model lists none, and an alternative derived from neither basic nor nonbasic as the status of an entity no basis covers.

A basis is never an error. Every solver snaps a side that the bounds do not have, a variable at an infinite upper bound say, and on an entity with two finite bounds a wrong side costs a bound flip, free under the dual simplex. A basis with the wrong number of basic entries costs more: Gurobi, COPT and GLPK start the solve cold, and the other solvers repair it. Removals are where that happens. A snapshot taken before a basic column was removed is one basic entity short over the model that remains: on Gurobi the next solve starts cold, and the eviction strategies of column generation evict nonbasic columns only, at an optimum, under the conditions given there. A handle created after a removal may reuse the removed variable's id, and then reads that variable's entry in an older snapshot: a snapshot describes the model as it was.

When a basis is available¶

get_basis() is defined right after a solve() that ran a simplex method, before any other call that changes the model or what the solver holds: it then returns the basis that solve ended with. Where the solver holds none, it throws mippp::solver_error, but soplex_lp returns a slack basis built from the bounds: after presolve proved the model infeasible or unbounded, after a time limit that stopped the solve before the simplex method started, and after an interior-point run that skipped its crossover, which only native settings can ask for. Before any solve, glpk_lp and soplex_lp, and copt_lp on a model with constraints, return the solver's starting basis, every row basic, and the other backends throw. Which of these each backend reports, and how, is in Support by model. A model without variables has one basis, every row basic, which get_basis() returns on every backend without a solve.

At any other time its result is unspecified: between a change of the model and the next solve(), between set_basis() and the next solve(), after a solve() that MIP++ answers without running the solver, as on a model whose bounds cross, and after refine_lp_status(), compute_iis() or compute_iis_by_deletion. Some backends then throw, others return the basis their next solve would start from.

The status afterwards¶

set_basis() is a change of the model like any other: until the next solve(), get_status() reports unknown and the result getters throw mippp::solve_required_error, see When results can be read. Gurobi, CPLEX and COPT discard their solution when a basis is written, so read what you need of the last solve before writing one.

Support by model¶

has_lp_basis<M> holds where get_basis() exists, and has_lp_basis_warm_start<M> where set_basis() does. Both are LP features: no *_milp model has them, and neither has highs_qp, whose quadratic solver neither reports a simplex basis nor starts from one.

Before the first solve and after any change of the model, get_basis() throws mippp::solve_required_error on every model, see When results can be read. The last column tells when a solved model has no basis to report, where get_basis() throws mippp::solver_error.

The backends differ in a few more ways, see LP bases:

• soplex_lp reads a basis but cannot write one: SoPlex's C interface has no setter.
• mosek_lp solves an LP by the interior point, which cannot start from a basis: set_basis() makes the next solve() run MOSEK's free simplex, unless you chose an optimizer through native_api(), and later solves return to the interior point.
• glpk_lp replaces a basis it cannot factorize by its standard basis, so the next solve starts cold.

Next¶

Re-solving and model updates — changing a model between solves without rebuilding it.