Solutions

This section of the manual describes how to access a solved solution to a problem. It uses the following model as an example:

julia> begin
           model = Model(HiGHS.Optimizer)
           set_silent(model)
           @variable(model, x >= 0)
           @variable(model, y[[:a, :b]] <= 1)
           @objective(model, Max, -12x - 20y[:a])
           @expression(model, my_expr, 6x + 8y[:a])
           @constraint(model, my_expr >= 100)
           @constraint(model, c1, 7x + 12y[:a] >= 120)
           optimize!(model)
           print(model)
       end
Max -12 x - 20 y[a]
Subject to
 6 x + 8 y[a] ≥ 100
 c1 : 7 x + 12 y[a] ≥ 120
 x ≥ 0
 y[a] ≤ 1
 y[b] ≤ 1

Check if an optimal solution exists

Use is_solved_and_feasible to check if the solver found an optimal solution:

julia> is_solved_and_feasible(model)
true

By default, is_solved_and_feasible returns true for both global and local optima. Pass allow_local = false to check if the solver found a globally optimal solution:

julia> is_solved_and_feasible(model; allow_local = false)
true

Pass dual = true to check if the solver found an optimal dual solution in addition to an optimal primal solution:

julia> is_solved_and_feasible(model; dual = true)
true

If this function returns false, use the functions mentioned below like solution_summary, termination_status, primal_status, and dual_status to understand what solution (if any) the solver found.

Solutions summary

solution_summary can be used for checking the summary of the optimization solutions.

julia> solution_summary(model)
solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : -2.05143e+02
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : -2.05143e+02
│ ├ dual_objective_value : -2.05143e+02
│ └ relative_gap         : 1.38546e-16
└ Work counters
  ├ solve_time (sec)   : 5.93345e-04
  ├ simplex_iterations : 2
  ├ barrier_iterations : 0
  └ node_count         : -1

julia> solution_summary(model; verbose = true)
solution_summary(; result = 1, verbose = true)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusOptimal
│ └ objective_bound    : -2.05143e+02
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : -2.05143e+02
│ ├ dual_objective_value : -2.05143e+02
│ ├ relative_gap         : 1.38546e-16
│ ├ value
│ │ ├ x : 1.54286e+01
│ │ ├ y[a] : 1.00000e+00
│ │ └ y[b] : 1.00000e+00
│ └ dual
│   └ c1 : 1.71429e+00
└ Work counters
  ├ solve_time (sec)   : 5.93345e-04
  ├ simplex_iterations : 2
  ├ barrier_iterations : 0
  └ node_count         : -1

Why did the solver stop?

Usetermination_status to understand why the solver stopped.

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

The MOI.TerminationStatusCode enum describes the full list of statuses that could be returned.

Common return values include OPTIMAL, LOCALLY_SOLVED, INFEASIBLE, DUAL_INFEASIBLE, and TIME_LIMIT.

Use raw_status to get a solver-specific string explaining why the optimization stopped:

julia> raw_status(model)
"kHighsModelStatusOptimal"

Primal solutions

Primal solution status

Use primal_status to return an MOI.ResultStatusCode enum describing the status of the primal solution.

julia> primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

Other common returns are NO_SOLUTION, and INFEASIBILITY_CERTIFICATE. The first means that the solver doesn't have a solution to return, and the second means that the primal solution is a certificate of dual infeasibility (a primal unbounded ray).

Objective values

The objective value of a solved problem can be obtained via objective_value. The best known bound on the optimal objective value can be obtained via objective_bound. If the solver supports it, the value of the dual objective can be obtained via dual_objective_value.

julia> objective_value(model)
-205.14285714285714

julia> objective_bound(model)
-205.1428571428571

julia> dual_objective_value(model)
-205.1428571428571

Primal solution values

If the solver has a primal solution to return, use value to access it:

julia> value(x)
15.428571428571429

Broadcast value over containers:

julia> value.(y)
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, [:a, :b]
And data, a 2-element Vector{Float64}:
 1.0
 1.0

value also works on expressions:

julia> value(my_expr)
100.57142857142857

and constraints:

julia> value(c1)
120.0

Unbounded rays

If the termination_status is DUAL_INFEASIBLE and primal_status is INFEASIBILITY_CERTIFICATE, then the value is a certificate of dual infeasibility. If a feasible primal solution exists, the certificate is an unbounded ray of the primal problem. The objective_value is the value of the objective, evaluated using the ray, excluding any constant term. For more details, see Infeasibility certificates.

julia> using HiGHS

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> set_attribute(model, "presolve", "off")

julia> @variable(model, x[1:2]);

julia> @objective(model, Max, sum(x) + 2.0);

julia> @constraint(model, 2 * x[1] + x[2] <= 0);

julia> optimize!(model)

julia> termination_status(model)
DUAL_INFEASIBLE::TerminationStatusCode = 3

julia> primal_status(model)
INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4

julia> d = value.(x)
2-element Vector{Float64}:
 -0.5
  1.0

julia> objective_value(model)
0.5

Dual solutions

Dual solution status

Use dual_status to return an MOI.ResultStatusCode enum describing the status of the dual solution.

julia> dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

Other common returns are NO_SOLUTION, and INFEASIBILITY_CERTIFICATE. The first means that the solver doesn't have a solution to return, and the second means that the dual solution is a certificate of primal infeasibility (a dual unbounded ray).

Dual solution values

If the solver has a dual solution to return, use dual to access it:

julia> dual(c1)
1.7142857142857142

Query the duals of variable bounds using LowerBoundRef, UpperBoundRef, and FixRef:

julia> dual(LowerBoundRef(x))
0.0

julia> dual.(UpperBoundRef.(y))
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, [:a, :b]
And data, a 2-element Vector{Float64}:
 -0.5714285714285694
  0.0

julia> shadow_price(c1)
1.7142857142857142

julia> reduced_cost(x)
-0.0

julia> reduced_cost.(y)
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, [:a, :b]
And data, a 2-element Vector{Float64}:
  0.5714285714285694
 -0.0

Unbounded rays

If the termination_status is INFEASIBLE and dual_status is INFEASIBILITY_CERTIFICATE, then the dual is a certificate of primal infeasibility. If a feasible dual solution exists, the certificate is an unbounded ray of the dual problem. The dual_objective_value is the value of the dual objective, evaluated using the ray, excluding any constant term. For more details, see Infeasibility certificates.

Fore more details, see Infeasibility certificates.

julia> using HiGHS

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> set_attribute(model, "presolve", "off")

julia> @variable(model, x[1:2] >= 0);

julia> @objective(model, Max, sum(x) + 0.5);

julia> @constraint(model, c, 2 * x[1] + x[2] <= -1);

julia> optimize!(model)

julia> termination_status(model)
INFEASIBLE::TerminationStatusCode = 2

julia> dual_status(model)
INFEASIBILITY_CERTIFICATE::ResultStatusCode = 4

julia> d_c = dual(c)
-1.0

julia> d_x = dual.(LowerBoundRef.(x))
2-element Vector{Float64}:
 2.0
 1.0

julia> dual_objective_value(model)
-1.0

Recommended workflow

You should always check whether the solver found a solution before calling solution functions like value or objective_value.

A simple approach for small scripts and notebooks is to use is_solved_and_feasible:

julia> function solve_and_print_solution(model)
           optimize!(model)
           if !is_solved_and_feasible(model; dual = true)
               error(
                   """
                   The model was not solved correctly:
                   termination_status : $(termination_status(model))
                   primal_status      : $(primal_status(model))
                   dual_status        : $(dual_status(model))
                   raw_status         : $(raw_status(model))
                   """,
               )
           end
           println("Solution is optimal")
           println("  objective value = ", objective_value(model))
           println("  primal solution: x = ", value(x))
           println("  dual solution: c1 = ", dual(c1))
           return
       end
solve_and_print_solution (generic function with 1 method)

julia> solve_and_print_solution(model)
Solution is optimal
  objective value = -205.14285714285714
  primal solution: x = 15.428571428571429
  dual solution: c1 = 1.7142857142857142

For code like libraries that should be more robust to the range of possible termination and result statuses, do some variation of the following:

julia> function solve_and_print_solution(model)
           optimize!(model)
           status = termination_status(model)
           if status in (OPTIMAL, LOCALLY_SOLVED)
               println("Solution is optimal")
           elseif status in (ALMOST_OPTIMAL, ALMOST_LOCALLY_SOLVED)
               println("Solution is optimal to a relaxed tolerance")
           elseif status == TIME_LIMIT
               println(
                   "Solver stopped due to a time limit. If a solution is available, " *
                   "it may be suboptimal."
               )
           elseif status in (
               ITERATION_LIMIT, NODE_LIMIT, SOLUTION_LIMIT, MEMORY_LIMIT,
               OBJECTIVE_LIMIT, NORM_LIMIT, OTHER_LIMIT,
           )
               println(
                   "Solver stopped due to a limit. If a solution is available, it " *
                   "may be suboptimal."
               )
           elseif status in (INFEASIBLE, LOCALLY_INFEASIBLE)
               println("The problem is primal infeasible")
           elseif status == DUAL_INFEASIBLE
               println(
                   "The problem is dual infeasible. If a primal feasible solution " *
                   "exists, the problem is unbounded. To check, set the objective " *
                   "to `@objective(model, Min, 0)` and re-solve. If the problem is " *
                   "feasible, the primal is unbounded. If the problem is " *
                   "infeasible, both the primal and dual are infeasible.",
               )
           elseif status == INFEASIBLE_OR_UNBOUNDED
               println(
                   "The model is either infeasible or unbounded. Set the objective " *
                   "to `@objective(model, Min, 0)` and re-solve to disambiguate. If " *
                   "the problem was infeasible, it will still be infeasible. If the " *
                   "problem was unbounded, it will now have a finite optimal solution.",
               )
           else
               println(
                   "The model was not solved correctly. The termination status is $status",
               )
           end
           if primal_status(model) in (FEASIBLE_POINT, NEARLY_FEASIBLE_POINT)
               println("  objective value = ", objective_value(model))
               println("  primal solution: x = ", value(x))
           elseif primal_status(model) == INFEASIBILITY_CERTIFICATE
               println("  primal certificate: x = ", value(x))
           end
           if dual_status(model) in (FEASIBLE_POINT, NEARLY_FEASIBLE_POINT)
               println("  dual solution: c1 = ", dual(c1))
           elseif dual_status(model) == INFEASIBILITY_CERTIFICATE
               println("  dual certificate: c1 = ", dual(c1))
           end
           return
       end
solve_and_print_solution (generic function with 1 method)

julia> solve_and_print_solution(model)
Solution is optimal
  objective value = -205.14285714285714
  primal solution: x = 15.428571428571429
  dual solution: c1 = 1.7142857142857142

OptimizeNotCalled errors

Due to differences in how solvers cache solutions internally, modifying a model after calling optimize! will reset the model into the OPTIMIZE_NOT_CALLED state. If you then attempt to query solution information, an OptimizeNotCalled error will be thrown.

If you are iteratively querying solution information and modifying a model, query all the results first, then modify the problem.

For example, instead of:

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x >= 0);

julia> optimize!(model)

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> set_upper_bound(x, 1)

julia> x_val = value(x)
┌ Warning: The model has been modified since the last call to `optimize!` (or `optimize!` has not been called yet). If you are iteratively querying solution information and modifying a model, query all the results first, then modify the model.
└ @ JuMP ~/.julia/dev/JuMP/src/optimizer_interface.jl:1085
ERROR: OptimizeNotCalled()
Stacktrace:
[...]

julia> termination_status(model)
OPTIMIZE_NOT_CALLED::TerminationStatusCode = 0

do

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x >= 0);

julia> optimize!(model);

julia> x_val = value(x)
0.0

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> set_upper_bound(x, 1)

julia> set_lower_bound(x, x_val)

julia> termination_status(model)
OPTIMIZE_NOT_CALLED::TerminationStatusCode = 0

If you know that your particular solver supports querying solution information after modifications, you can use direct_model to bypass the OPTIMIZE_NOT_CALLED state:

julia> model = direct_model(HiGHS.Optimizer());

julia> set_silent(model)

julia> @variable(model, x >= 0);

julia> optimize!(model)

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

julia> set_upper_bound(x, 1)

julia> x_val = value(x)
0.0

julia> set_lower_bound(x, x_val)

julia> termination_status(model)
OPTIMAL::TerminationStatusCode = 1

Accessing attributes

MathOptInterface defines many model attributes that can be queried. Some attributes can be directly accessed by getter functions. These include:

• solve_time
• relative_gap
• simplex_iterations
• barrier_iterations
• node_count

Sensitivity analysis for LP

Given an LP problem and an optimal solution corresponding to a basis, we can question how much an objective coefficient or standard form right-hand side coefficient (c.f., normalized_rhs) can change without violating primal or dual feasibility of the basic solution.

Note that not all solvers compute the basis, and for sensitivity analysis, the solver interface must implement MOI.ConstraintBasisStatus.

To give a simple example, we could analyze the sensitivity of the optimal solution to the following (non-degenerate) LP problem:

julia> begin
           model = Model(HiGHS.Optimizer)
           set_silent(model)
           @variable(model, x[1:2])
           set_lower_bound(x[2], -0.5)
           set_upper_bound(x[2], 0.5)
           @constraint(model, c1, x[1] + x[2] <= 1)
           @constraint(model, c2, x[1] - x[2] <= 1)
           @objective(model, Max, x[1])
           print(model)
       end
Max x[1]
Subject to
 c1 : x[1] + x[2] ≤ 1
 c2 : x[1] - x[2] ≤ 1
 x[2] ≥ -0.5
 x[2] ≤ 0.5

To analyze the sensitivity of the problem we could check the allowed perturbation ranges of, for example, the cost coefficients and the right-hand side coefficient of the constraint c1 as follows:

julia> optimize!(model)

julia> value.(x)
2-element Vector{Float64}:
  1.0
 -0.0

julia> report = lp_sensitivity_report(model);

julia> x1_lo, x1_hi = report[x[1]]
(-1.0, Inf)

julia> println("The objective coefficient of x[1] could decrease by $(x1_lo) or increase by $(x1_hi).")
The objective coefficient of x[1] could decrease by -1.0 or increase by Inf.

julia> x2_lo, x2_hi = report[x[2]]
(-1.0, 1.0)

julia> println("The objective coefficient of x[2] could decrease by $(x2_lo) or increase by $(x2_hi).")
The objective coefficient of x[2] could decrease by -1.0 or increase by 1.0.

julia> c_lo, c_hi = report[c1]
(-1.0, 1.0)

julia> println("The RHS of c1 could decrease by $(c_lo) or increase by $(c_hi).")
The RHS of c1 could decrease by -1.0 or increase by 1.0.

The range associated with a variable is the range of the allowed perturbation of the corresponding objective coefficient. Note that the current primal solution remains optimal within this range; however the corresponding dual solution might change since a cost coefficient is perturbed. Similarly, the range associated with a constraint is the range of the allowed perturbation of the corresponding right-hand side coefficient. In this range the current dual solution remains optimal, but the optimal primal solution might change.

If the problem is degenerate, there are multiple optimal bases and hence these ranges might not be as intuitive and seem too narrow, for example, a larger cost coefficient perturbation might not invalidate the optimality of the current primal solution. Moreover, if a problem is degenerate, due to finite precision, it can happen that, for example, a perturbation seems to invalidate a basis even though it doesn't (again providing too narrow ranges). To prevent this, increase the atol keyword argument to lp_sensitivity_report. Note that this might make the ranges too wide for numerically challenging instances. Thus, do not blindly trust these ranges, especially not for highly degenerate or numerically unstable instances.

Conflicts

When the model you input is infeasible, some solvers can help you find the cause of this infeasibility by offering a conflict, that is, a subset of the constraints that create this infeasibility. Depending on the solver, this can also be called an IIS (irreducible inconsistent subsystem).

If supported by the solver, use compute_conflict! to trigger the computation of a conflict. Once this process is finished, query the MOI.ConflictStatus attribute to check if a conflict was found.

If found, copy the IIS to a new model using copy_conflict, which you can then print or write to a file for easier debugging:

julia> using JuMP

julia> import Gurobi

julia> model = Model(Gurobi.Optimizer);

julia> set_silent(model)

julia> @variable(model, x >= 0)
x

julia> @constraint(model, c1, x >= 2)
c1 : x ≥ 2.0

julia> @constraint(model, c2, x <= 1)
c2 : x ≤ 1.0

julia> optimize!(model)

julia> compute_conflict!(model)

julia> if get_attribute(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
           iis_model, _ = copy_conflict(model)
           print(iis_model)
       end
Feasibility
Subject to
 c1 : x ≥ 2.0
 c2 : x ≤ 1.0

If you need more control over the list of constraints that appear in the conflict, iterate over the list of constraints and query the MOI.ConstraintConflictStatus attribute:

julia> list_of_conflicting_constraints = ConstraintRef[]
ConstraintRef[]

julia> for (F, S) in list_of_constraint_types(model)
           for con in all_constraints(model, F, S)
               if get_attribute(con, MOI.ConstraintConflictStatus()) == MOI.IN_CONFLICT
                   push!(list_of_conflicting_constraints, con)
               end
           end
       end

julia> list_of_conflicting_constraints
2-element Vector{ConstraintRef}:
 c1 : x ≥ 2.0
 c2 : x ≤ 1.0

Multiple solutions

Some solvers support returning multiple solutions. You can check how many solutions are available to query using result_count.

Functions for querying the solutions, for example, primal_status, dual_status, value, dual, and solution_summary all take an additional keyword argument result which can be used to specify which result to return.

julia> using JuMP

julia> import MultiObjectiveAlgorithms as MOA

julia> import HiGHS

julia> model = Model(() -> MOA.Optimizer(HiGHS.Optimizer));

julia> set_attribute(model, MOA.Algorithm(), MOA.Dichotomy())

julia> set_silent(model)

julia> @variable(model, x1 >= 0)
x1

julia> @variable(model, 0 <= x2 <= 3)
x2

julia> @objective(model, Min, [3x1 + x2, -x1 - 2x2])
2-element Vector{AffExpr}:
 3 x1 + x2
 -x1 - 2 x2

julia> @constraint(model, 3x1 - x2 <= 6)
3 x1 - x2 ≤ 6

julia> optimize!(model)

julia> solution_summary(model; result = 1)
solution_summary(; result = 1, verbose = false)
├ solver_name          : MOA[algorithm=MultiObjectiveAlgorithms.Dichotomy, optimizer=HiGHS]
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 3
│ ├ raw_status         : Solve complete. Found 3 solution(s)
│ └ objective_bound    : [0.00000e+00,-9.00000e+00]
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ └ objective_value      : [0.00000e+00,0.00000e+00]
└ Work counters
  └ solve_time (sec)   : 1.88589e-03

julia> for i in 1:result_count(model)
           println("Solution $i")
           println("   x = ", value.([x1, x2]; result = i))
           println(" obj = ", objective_value(model; result = i))
       end
Solution 1
   x = [0.0, 0.0]
 obj = [0.0, 0.0]
Solution 2
   x = [0.0, 3.0]
 obj = [3.0, -6.0]
Solution 3
   x = [3.0, 3.0]
 obj = [12.0, -9.0]

Checking feasibility of solutions

To check the feasibility of a primal solution, use primal_feasibility_report, which takes a model, a dictionary mapping each variable to a primal solution value (defaults to the last solved solution), and a tolerance atol (defaults to 0.0).

The function returns a dictionary which maps the infeasible constraint references to the distance between the primal value of the constraint and the nearest point in the corresponding set. A point is classed as infeasible if the distance is greater than the supplied tolerance atol.

julia> model = Model(HiGHS.Optimizer);

julia> set_silent(model)

julia> @variable(model, x >= 1, Int);

julia> @variable(model, y);

julia> @constraint(model, c1, x + y <= 1.95);

julia> point = Dict(x => 1.9, y => 0.06);

julia> primal_feasibility_report(model, point)
Dict{Any, Float64} with 2 entries:
  x integer         => 0.1
  c1 : x + y ≤ 1.95 => 0.01

julia> primal_feasibility_report(model, point; atol = 0.02)
Dict{Any, Float64} with 1 entry:
  x integer => 0.1

If the point is feasible, an empty dictionary is returned:

julia> primal_feasibility_report(model, Dict(x => 1.0, y => 0.0))
Dict{Any, Float64}()

To use the primal solution from a solve, omit the point argument:

julia> optimize!(model)

julia> primal_feasibility_report(model; atol = 0.0)
Dict{Any, Float64}()

Calling primal_feasibility_report without the point argument is useful when primal_status is FEASIBLE_POINT or NEARLY_FEASIBLE_POINT, and you want to assess the solution quality.

Pass skip_mising = true to skip constraints which contain variables that are not in point:

julia> primal_feasibility_report(model, Dict(x => 2.1); skip_missing = true)
Dict{Any, Float64} with 1 entry:
  x integer => 0.1

You can also use the functional form, where the first argument is a function that maps variables to their primal values:

julia> optimize!(model)

julia> primal_feasibility_report(v -> value(v), model)
Dict{Any, Float64}()