Query Solutions

# Querying Solutions

So far we have seen all the elements and constructs related to writing a JuMP optimization model. In this section we reach the point of what to do with a solved problem. Suppose your model is named `model`. Right after the call to `optimize!(model)`, it's natural to ask JuMP questions about the finished optimization step. Typical questions include:

• Why has the optimization process stopped? Did it hit the time limit or run into numerical issues?
• Do I have a solution to my problem?
• Is it optimal?
• Do I have a dual solution?
• How sensitive is the solution to data perturbations?

JuMP follows closely the concepts defined in MathOptInterface (MOI) to answer user questions about a finished call to `optimize!(model)`. There are three main steps in querying a solution:

First, we can query the `termination_status` which will tell us why the optimization stopped. This could be due to a number of reasons. For example, the solver found an optimal solution, the problem was proven to be infeasible, or a user-provided limit such as a time limit was encountered. For more information, see the Termination statuses section below.

Second, we can query the `primal_status` and the `dual_status`, which will tell us what kind of results we have for our primal and dual solutions. This might be an optimal primal-dual pair, a primal solution without a corresponding dual solution, or a certificate of primal or dual infeasibility. For more information, see the Solution statuses section below.

Third, we can query `value` and `dual` to obtain the primal and dual values of the optimization variables and constraints (if there are values to be queried).

## Termination statuses

The reason why the optimization of `model` was finished is given by

``termination_status(model)``

This function will return a `MOI.TerminationStatusCode` `enum`.

``TerminationStatusCode``

An Enum of possible values for the `TerminationStatus` attribute. This attribute is meant to explain the reason why the optimizer stopped executing in the most recent call to `optimize!`.

If no call has been made to `optimize!`, then the `TerminationStatus` is:

• `OPTIMIZE_NOT_CALLED`: The algorithm has not started.

OK

These are generally OK statuses, i.e., the algorithm ran to completion normally.

• `OPTIMAL`: The algorithm found a globally optimal solution.
• `INFEASIBLE`: The algorithm concluded that no feasible solution exists.
• `DUAL_INFEASIBLE`: The algorithm concluded that no dual bound exists for the problem. If, additionally, a feasible (primal) solution is known to exist, this status typically implies that the problem is unbounded, with some technical exceptions.
• `LOCALLY_SOLVED`: The algorithm converged to a stationary point, local optimal solution, could not find directions for improvement, or otherwise completed its search without global guarantees.
• `LOCALLY_INFEASIBLE`: The algorithm converged to an infeasible point or otherwise completed its search without finding a feasible solution, without guarantees that no feasible solution exists.
• `INFEASIBLE_OR_UNBOUNDED`: The algorithm stopped because it decided that the problem is infeasible or unbounded; this occasionally happens during MIP presolve.

Solved to relaxed tolerances

• `ALMOST_OPTIMAL`: The algorithm found a globally optimal solution to relaxed tolerances.
• `ALMOST_INFEASIBLE`: The algorithm concluded that no feasible solution exists within relaxed tolerances.
• `ALMOST_DUAL_INFEASIBLE`: The algorithm concluded that no dual bound exists for the problem within relaxed tolerances.
• `ALMOST_LOCALLY_SOLVED`: The algorithm converged to a stationary point, local optimal solution, or could not find directions for improvement within relaxed tolerances.

Limits

The optimizer stopped because of some user-defined limit.

• `ITERATION_LIMIT`: An iterative algorithm stopped after conducting the maximum number of iterations.
• `TIME_LIMIT`: The algorithm stopped after a user-specified computation time.
• `NODE_LIMIT`: A branch-and-bound algorithm stopped because it explored a maximum number of nodes in the branch-and-bound tree.
• `SOLUTION_LIMIT`: The algorithm stopped because it found the required number of solutions. This is often used in MIPs to get the solver to return the first feasible solution it encounters.
• `MEMORY_LIMIT`: The algorithm stopped because it ran out of memory.
• `OBJECTIVE_LIMIT`: The algorithm stopped because it found a solution better than a minimum limit set by the user.
• `NORM_LIMIT`: The algorithm stopped because the norm of an iterate became too large.
• `OTHER_LIMIT`: The algorithm stopped due to a limit not covered by one of the above.

Problematic

This group of statuses means that something unexpected or problematic happened.

• `SLOW_PROGRESS`: The algorithm stopped because it was unable to continue making progress towards the solution.
• `NUMERICAL_ERROR`: The algorithm stopped because it encountered unrecoverable numerical error.
• `INVALID_MODEL`: The algorithm stopped because the model is invalid.
• `INVALID_OPTION`: The algorithm stopped because it was provided an invalid option.
• `INTERRUPTED`: The algorithm stopped because of an interrupt signal.
• `OTHER_ERROR`: The algorithm stopped because of an error not covered by one of the statuses defined above.
source

Additionally, we can receive a solver specific string explaining why the optimization stopped with `raw_status`.

## Solution statuses

These statuses indicate what kind of result is available to be queried with `value` and `dual`. It's possible that no result is available to be queried.

We can obtain these statuses by calling `primal_status` for the primal status, and `dual_status` for the dual status. Both will return a `MOI.ResultStatusCode` `enum`.

``ResultStatusCode``

An Enum of possible values for the `PrimalStatus` and `DualStatus` attributes. The values indicate how to interpret the result vector.

• `NO_SOLUTION`: the result vector is empty.
• `FEASIBLE_POINT`: the result vector is a feasible point.
• `NEARLY_FEASIBLE_POINT`: the result vector is feasible if some constraint tolerances are relaxed.
• `INFEASIBLE_POINT`: the result vector is an infeasible point.
• `INFEASIBILITY_CERTIFICATE`: the result vector is an infeasibility certificate. If the `PrimalStatus` is `INFEASIBILITY_CERTIFICATE`, then the primal result vector is a certificate of dual infeasibility. If the `DualStatus` is `INFEASIBILITY_CERTIFICATE`, then the dual result vector is a proof of primal infeasibility.
• `NEARLY_INFEASIBILITY_CERTIFICATE`: the result satisfies a relaxed criterion for a certificate of infeasibility.
• `REDUCTION_CERTIFICATE`: the result vector is an ill-posed certificate; see this article for details. If the `PrimalStatus` is `REDUCTION_CERTIFICATE`, then the primal result vector is a proof that the dual problem is ill-posed. If the `DualStatus` is `REDUCTION_CERTIFICATE`, then the dual result vector is a proof that the primal is ill-posed.
• `NEARLY_REDUCTION_CERTIFICATE`: the result satisfies a relaxed criterion for an ill-posed certificate.
• `UNKNOWN_RESULT_STATUS`: the result vector contains a solution with an unknown interpretation.
• `OTHER_RESULT_STATUS`: the result vector contains a solution with an interpretation not covered by one of the statuses defined above.
source

Common status situations are described in the MOI docs.

## Obtaining solutions

Provided the primal status is not `MOI.NO_SOLUTION`, the primal solution can be obtained by calling `value`. For the dual solution, the function is `dual`. Calling `has_values` for the primal status and `has_duals` for the dual solution is an equivalent way to check whether the status is `MOI.NO_SOLUTION`.

It is important to note that if `has_values` or `has_duals` return false, calls to `value` and `dual` might throw an error or return arbitrary values.

The container type (e.g., scalar, vector, or matrix) of the returned solution (primal or dual) depends on the type of the variable or constraint. See `AbstractShape` and `dual_shape` for details.

Info

To call `value` or `dual` on containers of `VariableRef` or `ConstraintRef`, use the broadcast syntax, e.g., `value.(x)`.

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`.

The following is a recommended workflow for solving a model and querying the solution:

``````using JuMP
model = Model()
@variable(model, x[1:10] >= 0)
# ... other constraints ...
optimize!(model)

if termination_status(model) == MOI.OPTIMAL
optimal_solution = value.(x)
optimal_objective = objective_value(model)
elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
suboptimal_solution = value.(x)
suboptimal_objective = objective_value(model)
else
error("The model was not solved correctly.")
end``````
Querying after modification

If a solved model is modified, then querying the solution is undefined behavior. Adding, deleting, or modifying a constraint (or variable) may invalidate any part of the solution.

## Accessing MathOptInterface attributes

MathOptInterface defines a large number of model attributes that can be queried. Examples include `MOI.RelativeGap` and `MOI.SimplexIterations`.

Some attributes can be directly accessed by getter functions. These include

To query these attributes, use:

``````using JuMP
model = Model()
# ...
optimize!(model)

@show relative_gap(model)
# or
@show MOI.get(model, MOI.RelativeGap())

@show simplex_iterations(model)
# or
@show MOI.get(model, MOI.SimplexIterations())``````

## 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 rhs 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`.

Given an LP optimal solution (and both `has_values` and `has_duals` returns `true`) `lp_objective_perturbation_range` returns a range of the allowed perturbation of the cost coefficient corresponding to the input variable. 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, `lp_rhs_perturbation_range` returns a range of the allowed perturbation of the rhs coefficient corresponding to the input constraint. And in this range the current dual solution remains optimal but the primal solution might change since a rhs coefficient is perturbed.

However, if the problem is degenerate, there are multiple optimal bases and hence these ranges might not be as intuitive and seem too narrow. E.g., 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, e.g., a perturbation seems to invalidate a basis even though it doesn't (again providing too narrow ranges). To prevent this `feasibility_tolerance` and `optimality_tolerance` is introduced, which in turn, 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.

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

``````julia> model = Model();
julia> @variable(model, x[1:2]);
julia> @constraint(model, c1, x + x <= 1);
julia> @constraint(model, c2, x - x <= 1);
julia> @constraint(model, c3, -0.5 <= x <= 0.5);
julia> @objective(model, Max, x);``````

To analyze the sensitivity of the problem we could check the allowed perturbation ranges of, e.g., the cost coefficients and the rhs coefficient of constraint `c1` as follows:

``````julia> optimize!(model);

julia> value.(x)
2-element Array{Float64,1}:
1.0
0.0
julia> lp_objective_perturbation_range(x)
(-1.0, Inf)
julia> lp_objective_perturbation_range(x)
(-1.0, 1.0)
julia> lp_rhs_perturbation_range(c1)
(-1.0, 1.0)``````
``````lp_objective_perturbation_range(var::VariableRef;
optimality_tolerance::Float64)
::Tuple{Float64, Float64}``````

Gives the range by which the cost coefficient can change and the current LP basis remains optimal, i.e., the reduced costs remain valid.

Notes

• The range denotes valid changes, Δ ∈ [l, u], for which cost[var] += Δ do not violate the current optimality conditions.
• `optimality_tolerance` is the dual feasibility tolerance, this should preferably match the tolerance used by the solver. The defualt tolerance should however apply in most situations (c.f. "Computational Techniques of the Simplex Method" by István Maros, section 9.3.4).
source
``````lp_rhs_perturbation_range(constraint::ConstraintRef;
feasibility_tolerance::Float64)
::Tuple{Float64, Float64}``````

Gives the range by which the rhs coefficient can change and the current LP basis remains feasible, i.e., where the shadow prices apply.

Notes

• The rhs coefficient is the value right of the relation, i.e., b for the constraint when of the form a*x □ b, where □ is ≤, =, or ≥.
• The range denotes valid changes, e.g., for a*x <= b + Δ, the LP basis remains feasible for all Δ ∈ [l, u].
• `feasibility_tolerance` is the primal feasibility tolerance, this should preferably match the tolerance used by the solver. The default tolerance should however apply in most situations (c.f. "Computational Techniques of the Simplex Method" by István Maros, section 9.3.4).
source

## Conflicts

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

The function `compute_conflict!` is used to trigger the computation of a conflict. Once this process is finished, the attribute `MOI.ConflictStatus` returns a `MOI.ConflictStatusCode`.

If there is a conflict, you can query from each constraint whether it participates in the conflict or not using the attribute `MOI.ConstraintConflictStatus`, which returns a `MOI.ConflictParticipationStatusCode`.

For instance, this is how you can use this functionality:

``````using JuMP
model = Model() # You must use a solver that supports conflict refining/IIS computation, like CPLEX or Gurobi
@variable(model, x >= 0)
@constraint(model, c1, x >= 2)
@constraint(model, c2, x <= 1)
optimize!(model)

# termination_status(model) will likely be MOI.INFEASIBLE,
# depending on the solver

compute_conflict!(model)
if MOI.get(model, MOI.ConflictStatus()) != MOI.CONFLICT_FOUND
error("No conflict could be found for an infeasible model.")
end

# Both constraints should participate in the conflict.
MOI.get(model, MOI.ConstraintConflictStatus(), c1)
MOI.get(model, MOI.ConstraintConflictStatus(), c2)``````

## 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, e.g., `primal_status` and `value`, all take an additional keyword argument `result` which can be used to specify which result to return.

``````using JuMP
model = Model()
@variable(model, x[1:10] >= 0)
# ... other constraints ...
optimize!(model)

if termination_status(model) != MOI.OPTIMAL
error("The model was not solved correctly.")
end

num_results = result_count(model)
@assert has_values(model; result = num_results)
an_optimal_solution = value.(x; result = num_results)
an_optimal_objective = objective_value(model; result = num_results)``````

## Reference

``termination_status(model::Model)``

Return the reason why the solver stopped (i.e., the MathOptInterface model attribute `TerminationStatus`).

source
``raw_status(model::Model)``

Return the reason why the solver stopped in its own words (i.e., the MathOptInterface model attribute `RawStatusString`).

source
``primal_status(model::Model; result::Int = 1)``

Return the status of the most recent primal solution of the solver (i.e., the MathOptInterface model attribute `PrimalStatus`) associated with the result index `result`.

See also: `result_count`.

source
``has_values(model::Model; result::Int = 1)``

Return `true` if the solver has a primal solution in result index `result` available to query, otherwise return `false`.

See also `value` and `result_count`.

source
``value(con_ref::ConstraintRef; result::Int = 1)``

Return the primal value of constraint `con_ref` associated with result index `result` of the most-recent solution returned by the solver.

That is, if `con_ref` is the reference of a constraint `func`-in-`set`, it returns the value of `func` evaluated at the value of the variables (given by `value(::VariableRef)`).

Use `has_values` to check if a result exists before asking for values.

See also: `result_count`.

Note

For scalar contraints, the constant is moved to the `set` so it is not taken into account in the primal value of the constraint. For instance, the constraint `@constraint(model, 2x + 3y + 1 == 5)` is transformed into `2x + 3y`-in-`MOI.EqualTo(4)` so the value returned by this function is the evaluation of `2x + 3y`. ```

source
``value(con_ref::ConstraintRef, var_value::Function)``

Evaluate the primal value of the constraint `con_ref` using `var_value(v)` as the value for each variable `v`.

source
``value(v::VariableRef; result = 1)``

Return the value of variable `v` associated with result index `result` of the most-recent returned by the solver.

Use `has_values` to check if a result exists before asking for values.

See also: `result_count`.

source
``value(v::VariableRef, var_value::Function)``

Evaluate the value of the variable `v` as `var_value(v)`.

source
``value(ex::GenericAffExpr, var_value::Function)``

Evaluate `ex` using `var_value(v)` as the value for each variable `v`.

source
``value(v::GenericAffExpr; result::Int = 1)``

Return the value of the `GenericAffExpr` `v` associated with result index `result` of the most-recent solution returned by the solver.

Replaces `getvalue` for most use cases.

See also: `result_count`.

source
``value(v::GenericQuadExpr; result::Int = 1)``

Return the value of the `GenericQuadExpr` `v` associated with result index `result` of the most-recent solution returned by the solver.

Replaces `getvalue` for most use cases.

See also: `result_count`.

source
``value(p::NonlinearParameter)``

Return the current value stored in the nonlinear parameter `p`.

Example

``````model = Model()
@NLparameter(model, p == 10)
value(p)

# output
10.0``````
source
``value(ex::NonlinearExpression, var_value::Function)``

Evaluate `ex` using `var_value(v)` as the value for each variable `v`.

source
``value(ex::NonlinearExpression; result::Int = 1)``

Return the value of the `NonlinearExpression` `ex` associated with result index `result` of the most-recent solution returned by the solver.

Replaces `getvalue` for most use cases.

See also: `result_count`.

source
``dual_status(model::Model; result::Int = 1)``

Return the status of the most recent dual solution of the solver (i.e., the MathOptInterface model attribute `DualStatus`) associated with the result index `result`.

See also: `result_count`.

source
``has_duals(model::Model; result::Int = 1)``

Return `true` if the solver has a dual solution in result index `result` available to query, otherwise return `false`.

See also `dual`, `shadow_price`, and `result_count`.

source
``dual(con_ref::ConstraintRef; result::Int = 1)``

Return the dual value of constraint `con_ref` associated with result index `result` of the most-recent solution returned by the solver.

Use `has_dual` to check if a result exists before asking for values.

See also: `result_count`, `shadow_price`.

source
``solve_time(model::Model)``

If available, returns the solve time reported by the solver. Returns "ArgumentError: ModelLike of type `Solver.Optimizer` does not support accessing the attribute MathOptInterface.SolveTime()" if the attribute is not implemented.

source
``struct OptimizeNotCalled <: Exception end``

A result attribute cannot be queried before `optimize!` is called.

source
``optimize!(optimizer::AbstractOptimizer)``

Start the solution procedure.

``result_count(model::Model)``

Return the number of results available to query after a call to `optimize!`.

source
``relative_gap(model::Model)``

Return the final relative optimality gap after a call to `optimize!(model)`. Exact value depends upon implementation of MathOptInterface.RelativeGap() by the particular solver used for optimization.

source
``simplex_iterations(model::Model)``

Gets the cumulative number of simplex iterations during the most-recent optimization.

Solvers must implement `MOI.SimplexIterations()` to use this function.

source
``barrier_iterations(model::Model)``

Gets the cumulative number of barrier iterations during the most recent optimization.

Solvers must implement `MOI.BarrierIterations()` to use this function.

source
``node_count(model::Model)``

Gets the total number of branch-and-bound nodes explored during the most recent optimization in a Mixed Integer Program.

Solvers must implement `MOI.NodeCount()` to use this function.

source
``compute_conflict!(model::Model)``

Compute a conflict if the model is infeasible. If an optimizer has not been set yet (see `set_optimizer`), a `NoOptimizer` error is thrown.

The status of the conflict can be checked with the `MOI.ConflictStatus` model attribute. Then, the status for each constraint can be queried with the `MOI.ConstraintConflictStatus` attribute.

source
``compute_conflict!(optimizer::AbstractOptimizer)``

Computes a minimal subset of constraints such that the model with the other constraint removed is still infeasible.

Some solvers call a set of conflicting constraints an Irreducible Inconsistent Subsystem (IIS).

Note

If the model is modified after a call to `compute_conflict!`, the implementor is not obliged to purge the conflict. Any calls to the above attributes may return values for the original conflict without a warning. Similarly, when modifying the model, the conflict can be discarded.

``ConflictStatus()``

A model attribute for the `ConflictStatusCode` explaining why the conflict refiner stopped when computing the conflict.

``ConflictStatusCode``

An Enum of possible values for the `ConflictStatus` attribute. This attribute is meant to explain the reason why the conflict finder stopped executing in the most recent call to `compute_conflict!`.

Possible values are:

• `COMPUTE_CONFLICT_NOT_CALLED`: the function `compute_conflict!` has not yet been called
• `NO_CONFLICT_EXISTS`: there is no conflict because the problem is feasible
• `NO_CONFLICT_FOUND`: the solver could not find a conflict
• `CONFLICT_FOUND`: at least one conflict could be found
source
``ConstraintConflictStatus()``

A constraint attribute indicating whether the constraint participates in the conflict. Its type is `ConflictParticipationStatusCode`.

``ConflictParticipationStatusCode``

An Enum of possible values for the `ConstraintConflictStatus` attribute. This attribute is meant to indicate whether a given constraint participates or not in the last computed conflict.

Possible values are:

• `NOT_IN_CONFLICT`: the constraint does not participate in the conflict
• `IN_CONFLICT`: the constraint participates in the conflict
• `MAYBE_IN_CONFLICT`: the constraint may participate in the conflict, the solver was not able to prove that the constraint can be excluded from the conflict
source