From d91e4d0e7a03ef5f5b812ecabf546b62623a171c Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 15 Feb 2024 14:32:33 +1300 Subject: [PATCH 1/5] [docs] add tutorial on querying the basis status --- docs/make.jl | 1 + docs/src/tutorials/linear/basis.jl | 134 +++++++++++++++++++++++++++++ 2 files changed, 135 insertions(+) create mode 100644 docs/src/tutorials/linear/basis.jl diff --git a/docs/make.jl b/docs/make.jl index a2b5014eaa6..aaf473ee8b6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -313,6 +313,7 @@ const _PAGES = [ "tutorials/linear/constraint_programming.md", "tutorials/linear/callbacks.md", "tutorials/linear/lp_sensitivity.md", + "tutorials/linear/basis.md", "tutorials/linear/mip_duality.md", ], "Nonlinear programs" => [ diff --git a/docs/src/tutorials/linear/basis.jl b/docs/src/tutorials/linear/basis.jl new file mode 100644 index 00000000000..68991d31f53 --- /dev/null +++ b/docs/src/tutorials/linear/basis.jl @@ -0,0 +1,134 @@ +# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors #src +# This Source Code Form is subject to the terms of the Mozilla Public License #src +# v.2.0. If a copy of the MPL was not distributed with this file, You can #src +# obtain one at https://mozilla.org/MPL/2.0/. #src + +# # Basis matrices + +# This tutorial explains how to query basis of a linear program. + +# ## Setup + +# This tutorial uses the following packages: + +using JuMP +import HiGHS + +# ## Standard form example + +# Consider the following example, which is from the Wikipedia article on +# [Basic feasible solutions](https://en.wikipedia.org/wiki/Basic_feasible_solution): +# +# ```math +# \begin{aligned} +# \max \; & 0 \\ +# \text{s.t.} & x_1 + 5x_2 + 3x_3 + 4x_4 + 6x_5 = 14 \\ +# & x_2 + 3x_3 + 5x_4 + 6x_5 = 7 \\ +# & x_i \ge 0, \forall i = 1,\ldots,5. +# \end{aligned} +# ``` +# +# The `A` matrix is: + +A = [1 5 3 4 6; 0 1 3 5 6] + +# and the right-hand side `b` vector is: + +b = [14, 7] + +# We can create an optimize the problem in the standard form: + +n = size(A, 2) +model = Model(HiGHS.Optimizer) +@variable(model, x[1:n] >= 0) +@constraint(model, A * x == b) +optimize!(model) +@assert is_solved_and_feasible(model) + +# This has a solution: + +value.(x) + +# Query the basis status using [`MOI.VariableBasisStatus`](@ref): + +MOI.get(model, MOI.VariableBasisStatus(), x[1]) + +# the result is a [`MOI.BasisStatusCode`](@ref). Query all of the basis statuses +# with the broadcast `MOI.get.(`: + +MOI.get.(model, MOI.VariableBasisStatus(), x) + +# The values are either [`MOI.NONBASIC_AT_LOWER`](@ref) or [`MOI.BASIC`](@ref). +# All of the [`MOI.NONBASIC_AT_LOWER`](@ref) have a value at their lower bound. +# The [`MOI.BASIC`](@ref) variables correspond to the columns of the optimal +# basis. + +# We can get the columns using: + +indices = + [MOI.get(model, MOI.VariableBasisStatus(), x[i]) == MOI.BASIC for i in 1:5] + +# Filter the basis matrix from `A`: + +B = A[:, indices] + +# Since the basis matrix is non-singular, solving the system `Bx = b` must yield +# the optimal primal solution: + +B \ b + +#- + +value.(x[indices]) + +# ## A more complicated example + +# Often, you may want to work with the basis of a model that is not in a nice +# standard form. For example: + +model = Model(HiGHS.Optimizer) +@variable(model, x >= 0) +@variable(model, 0 <= y <= 3) +@variable(model, z <= 1) +@objective(model, Min, 12x + 20y - z) +@constraint(model, c1, 6x + 8y >= 100) +@constraint(model, c2, 7x + 12y >= 120) +@constraint(model, c3, x + y <= 20) +optimize!(model) +@assert is_solved_and_feasible(model) +solution_summary(model; verbose = true) + +# A common way to query the basis status of every variable is: + +v_basis = Dict( + xi => MOI.get(model, MOI.VariableBasisStatus(), xi) for + xi in all_variables(model) +) + +# Despite having three constraints, there are only two basis variables. Since +# the basis matrix must be square, where is the other basis variable? + +# The answer is that solvers will reformulate inequality constraints: +# ```math +# \begin{aligned} +# \text{s.t.} & l \le A x \le u +# \end{aligned} +# ``` +# into the system: +# ```math +# \begin{aligned} +# \text{s.t.} & A x - Is = 0 +# & l \le s \le u +# \end{aligned} +# ``` +# Thus, for every inequality constraint there is a slack variable `s`. + +# Query the basis status of the slack variables associated with a constraint +# using [`MOI.ConstraintBasisStatus`](@ref): + +c_basis = Dict( + ci => MOI.get(model, MOI.ConstraintBasisStatus(), ci) for ci in + all_constraints(model; include_variable_in_set_constraints = false) +) + +# Thus, the basis is formed by `x`, `y`, and the slack associated with `c3`. From 2a88c07dabe562b72a85f23ce042eb17aa142d48 Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 15 Feb 2024 15:25:30 +1300 Subject: [PATCH 2/5] Update --- docs/src/tutorials/linear/basis.jl | 46 ++++++++++++++---------------- 1 file changed, 21 insertions(+), 25 deletions(-) diff --git a/docs/src/tutorials/linear/basis.jl b/docs/src/tutorials/linear/basis.jl index 68991d31f53..b0eca8283f7 100644 --- a/docs/src/tutorials/linear/basis.jl +++ b/docs/src/tutorials/linear/basis.jl @@ -22,8 +22,8 @@ import HiGHS # ```math # \begin{aligned} # \max \; & 0 \\ -# \text{s.t.} & x_1 + 5x_2 + 3x_3 + 4x_4 + 6x_5 = 14 \\ -# & x_2 + 3x_3 + 5x_4 + 6x_5 = 7 \\ +# \text{s.t.}\; & 1x_1 + 5x_2 + 3x_3 + 4x_4 + 6x_5 = 14 \\ +# & 0x_1 + 1x_2 + 3x_3 + 5x_4 + 6x_5 = 7 \\ # & x_i \ge 0, \forall i = 1,\ldots,5. # \end{aligned} # ``` @@ -40,6 +40,7 @@ b = [14, 7] n = size(A, 2) model = Model(HiGHS.Optimizer) +set_silent(model) @variable(model, x[1:n] >= 0) @constraint(model, A * x == b) optimize!(model) @@ -49,31 +50,30 @@ optimize!(model) value.(x) -# Query the basis status using [`MOI.VariableBasisStatus`](@ref): +# Query the basis status of a variable using [`MOI.VariableBasisStatus`](@ref): -MOI.get(model, MOI.VariableBasisStatus(), x[1]) +get_attribute(x[1], MOI.VariableBasisStatus()) # the result is a [`MOI.BasisStatusCode`](@ref). Query all of the basis statuses -# with the broadcast `MOI.get.(`: +# with the broadcast `get_attribute.(`: -MOI.get.(model, MOI.VariableBasisStatus(), x) +get_attribute.(x, MOI.VariableBasisStatus()) # The values are either [`MOI.NONBASIC_AT_LOWER`](@ref) or [`MOI.BASIC`](@ref). -# All of the [`MOI.NONBASIC_AT_LOWER`](@ref) have a value at their lower bound. -# The [`MOI.BASIC`](@ref) variables correspond to the columns of the optimal -# basis. +# All of the [`MOI.NONBASIC_AT_LOWER`](@ref) variables have a value at their +# lower bound. The [`MOI.BASIC`](@ref) variables correspond to the columns of +# the optimal basis. -# We can get the columns using: +# Get the columns using: -indices = - [MOI.get(model, MOI.VariableBasisStatus(), x[i]) == MOI.BASIC for i in 1:5] +indices = get_attribute.(x, MOI.VariableBasisStatus()) .== MOI.BASIC # Filter the basis matrix from `A`: B = A[:, indices] # Since the basis matrix is non-singular, solving the system `Bx = b` must yield -# the optimal primal solution: +# the optimal primal solution of the basic variables: B \ b @@ -87,6 +87,7 @@ value.(x[indices]) # standard form. For example: model = Model(HiGHS.Optimizer) +set_silent(model) @variable(model, x >= 0) @variable(model, 0 <= y <= 3) @variable(model, z <= 1) @@ -96,30 +97,25 @@ model = Model(HiGHS.Optimizer) @constraint(model, c3, x + y <= 20) optimize!(model) @assert is_solved_and_feasible(model) -solution_summary(model; verbose = true) # A common way to query the basis status of every variable is: v_basis = Dict( - xi => MOI.get(model, MOI.VariableBasisStatus(), xi) for + xi => get_attribute(xi, MOI.VariableBasisStatus()) for xi in all_variables(model) ) -# Despite having three constraints, there are only two basis variables. Since -# the basis matrix must be square, where is the other basis variable? +# Despite the model having three constraints, there are only two basis +# variables. Since the basis matrix must be square, where is the other basic +# variable? # The answer is that solvers will reformulate inequality constraints: # ```math -# \begin{aligned} -# \text{s.t.} & l \le A x \le u -# \end{aligned} +# l \le A x \le u # ``` # into the system: # ```math -# \begin{aligned} -# \text{s.t.} & A x - Is = 0 -# & l \le s \le u -# \end{aligned} +# A x - Is = 0, \quad l \le s \le u # ``` # Thus, for every inequality constraint there is a slack variable `s`. @@ -127,7 +123,7 @@ v_basis = Dict( # using [`MOI.ConstraintBasisStatus`](@ref): c_basis = Dict( - ci => MOI.get(model, MOI.ConstraintBasisStatus(), ci) for ci in + ci => get_attribute(ci, MOI.ConstraintBasisStatus()) for ci in all_constraints(model; include_variable_in_set_constraints = false) ) From 0655470eb8f2c5ee29a44eab929f057915d0f79f Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Thu, 15 Feb 2024 16:30:36 +1300 Subject: [PATCH 3/5] Apply suggestions from code review --- docs/src/tutorials/linear/basis.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/linear/basis.jl b/docs/src/tutorials/linear/basis.jl index b0eca8283f7..f824e8c96ea 100644 --- a/docs/src/tutorials/linear/basis.jl +++ b/docs/src/tutorials/linear/basis.jl @@ -5,7 +5,7 @@ # # Basis matrices -# This tutorial explains how to query basis of a linear program. +# This tutorial explains how to query the basis of a linear program. # ## Setup @@ -36,7 +36,7 @@ A = [1 5 3 4 6; 0 1 3 5 6] b = [14, 7] -# We can create an optimize the problem in the standard form: +# We can create and optimize the problem in the standard form: n = size(A, 2) model = Model(HiGHS.Optimizer) @@ -59,7 +59,7 @@ get_attribute(x[1], MOI.VariableBasisStatus()) get_attribute.(x, MOI.VariableBasisStatus()) -# The values are either [`MOI.NONBASIC_AT_LOWER`](@ref) or [`MOI.BASIC`](@ref). +# For this problem, the values are either [`MOI.NONBASIC_AT_LOWER`](@ref) or [`MOI.BASIC`](@ref). # All of the [`MOI.NONBASIC_AT_LOWER`](@ref) variables have a value at their # lower bound. The [`MOI.BASIC`](@ref) variables correspond to the columns of # the optimal basis. @@ -105,7 +105,7 @@ v_basis = Dict( xi in all_variables(model) ) -# Despite the model having three constraints, there are only two basis +# Despite the model having three constraints, there are only two basic # variables. Since the basis matrix must be square, where is the other basic # variable? From 3279f86e0eef03f71565590bba6f91efe14d046b Mon Sep 17 00:00:00 2001 From: odow Date: Sat, 17 Feb 2024 13:20:26 +1300 Subject: [PATCH 4/5] Add more docs --- docs/src/tutorials/linear/basis.jl | 74 +++++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/linear/basis.jl b/docs/src/tutorials/linear/basis.jl index f824e8c96ea..b33c29c7a73 100644 --- a/docs/src/tutorials/linear/basis.jl +++ b/docs/src/tutorials/linear/basis.jl @@ -111,11 +111,11 @@ v_basis = Dict( # The answer is that solvers will reformulate inequality constraints: # ```math -# l \le A x \le u +# A x \le b # ``` # into the system: # ```math -# A x - Is = 0, \quad l \le s \le u +# A x + Is = b # ``` # Thus, for every inequality constraint there is a slack variable `s`. @@ -128,3 +128,73 @@ c_basis = Dict( ) # Thus, the basis is formed by `x`, `y`, and the slack associated with `c3`. + +# A simple way to get the `A` matrix of an unstructured linear program is with +# [`lp_matrix_data`](@ref): + +matrix = lp_matrix_data(model) +matrix.A + +# You can check the permutation of the rows and columns using + +matrix.variables + +# and + +matrix.affine_constraints + +# We can construct the slack column associated with `c3` as: + +s_column = zeros(size(matrix.A, 1)) +s_column[3] = 1.0 + +# The full basis matrix is therefore: + +B = hcat(matrix.A[:, [1, 2]], s_column) + +# [`lp_matrix_data`](@ref) returns separate vectors for the lower and upper row +# bounds. Convert to a single right-hand side vector by taking the finite +# elements: + +b = ifelse.(isfinite.(matrix.b_lower), matrix.b_lower, matrix.b_upper) + +# Solving the Basis system as before yields: + +B \ b + +# which is the value of `x`, `y`, and the slack associated with `c3`. + +# ## Identifying degenerate variables + +# Another common task is identifying degenerate variables. A degenerate variable +# is a basic variable that has an optimal value at its lower or upper bound. + +# Here is a function that computes whether a variable is degenerate: + +function is_degenerate(x) + if get_attribute(x, MOI.VariableBasisStatus()) == MOI.BASIC + return (has_lower_bound(x) && ≈(value(x), lower_bound(x))) || + (has_upper_bound(x) && ≈(value(x), upper_bound(x))) + end + return false +end + +# A simple example of a linear program with a degenerate solution is: + +A, b, c = [1 1; 0 1], [1, 1], [1, 1] +model = Model(HiGHS.Optimizer); +set_silent(model) +@variable(model, x[1:2] >= 0) +@objective(model, Min, c' * x) +@constraint(model, A * x == b) +optimize!(model) +degenerate_variables = filter(is_degenerate, all_variables(model)) + +# The solution is degenerate because: + +value(x[1]) + +# and + +get_attribute(x[1], MOI.VariableBasisStatus()) + From 633f94d1e7dbb9b276601b2c84926057699285de Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 17 Feb 2024 15:40:13 +1300 Subject: [PATCH 5/5] Update docs/src/tutorials/linear/basis.jl --- docs/src/tutorials/linear/basis.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/src/tutorials/linear/basis.jl b/docs/src/tutorials/linear/basis.jl index b33c29c7a73..030f64428b7 100644 --- a/docs/src/tutorials/linear/basis.jl +++ b/docs/src/tutorials/linear/basis.jl @@ -197,4 +197,3 @@ value(x[1]) # and get_attribute(x[1], MOI.VariableBasisStatus()) -