Базисные матрицы

В этом руководстве объясняется, как запросить базис линейной программы.

Настройка

В этом руководстве используются следующие пакеты.

using JuMP
import HiGHS

Пример стандартной формы

Рассмотрим следующий пример из статьи в Википедии Базисные допустимые решения.

Матрица A имеет вид:

A = [1 5 3 4 6; 0 1 3 5 6]

2×5 Matrix{Int64}:
 1  5  3  4  6
 0  1  3  5  6

а вектор b правой части:

b = [14, 7]

2-element Vector{Int64}:
 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)
@assert is_solved_and_feasible(model)

У нее есть решение:

value.(x)

5-element Vector{Float64}:
 0.0
 2.0
 0.0
 1.0
 0.0

Запросим статус базиса переменной с помощью MOI.VariableBasisStatus:

get_attribute(x[1], MOI.VariableBasisStatus())

NONBASIC_AT_LOWER::BasisStatusCode = 2

Результатом является MOI.BasisStatusCode. Запросим все статусы базиса с помощью трансляции get_attribute.(:

get_attribute.(x, MOI.VariableBasisStatus())

5-element Vector{MathOptInterface.BasisStatusCode}:
 NONBASIC_AT_LOWER::BasisStatusCode = 2
 BASIC::BasisStatusCode = 0
 NONBASIC_AT_LOWER::BasisStatusCode = 2
 BASIC::BasisStatusCode = 0
 NONBASIC_AT_LOWER::BasisStatusCode = 2

Для этой задачи значением будет MOI.NONBASIC_AT_LOWER или MOI.BASIC. У всех переменных MOI.NONBASIC_AT_LOWER есть значение на нижней границе. Переменные MOI.BASIC соответствуют столбцам оптимального базиса.

Получим столбцы:

indices = get_attribute.(x, MOI.VariableBasisStatus()) .== MOI.BASIC

5-element BitVector:
 0
 1
 0
 1
 0

Отфильтруем базисную матрицу из A:

B = A[:, indices]

2×2 Matrix{Int64}:
 5  4
 1  5

Так как базисная матрица невырожденная, решение системы Bx = b должно давать оптимальное прямое решение для базисных переменных:

B \ b

2-element Vector{Float64}:
 2.0
 1.0

value.(x[indices])

2-element Vector{Float64}:
 2.0
 1.0

Более сложный пример

Часто приходится работать с базисом модели, которая не имеет удобной стандартной формы. Например:

model = Model(HiGHS.Optimizer)
set_silent(model)
@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)

Вот распространенный способ запроса статуса базиса каждой переменной:

v_basis = Dict(
    xi => get_attribute(xi, MOI.VariableBasisStatus()) for
    xi in all_variables(model)
)

Dict{VariableRef, MathOptInterface.BasisStatusCode} with 3 entries:
  y => BASIC
  x => BASIC
  z => NONBASIC_AT_UPPER

Хотя у модели три ограничения, в ней всего две базисные переменные. Но ведь базисная матрица должна быть квадратной. Где же еще одна базисная переменная?

Ответ кроется в том, что решатели переформулируют ограничения неравенства:

в систему:

Поэтому для каждого ограничения неравенства есть ослабляющая переменная s.

Запросить статус базиса ослабляющих переменных, связанных с ограничением, можно с помощью MOI.ConstraintBasisStatus:

c_basis = Dict(
    ci => get_attribute(ci, MOI.ConstraintBasisStatus()) for ci in
    all_constraints(model; include_variable_in_set_constraints = false)
)

Dict{ConstraintRef{Model, C, ScalarShape} where C, MathOptInterface.BasisStatusCode} with 3 entries:
  c2 : 7 x + 12 y ≥ 120 => NONBASIC
  c1 : 6 x + 8 y ≥ 100  => NONBASIC
  c3 : x + y ≤ 20       => BASIC

Таким образом, базис образуется из x, y и ослабляющей переменной, связанной с c3.

Матрицу A неструктурированной линейной программы можно легко получить с помощью lp_matrix_data:

matrix = lp_matrix_data(model)
matrix.A

3×3 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 6.0   8.0   ⋅
 7.0  12.0   ⋅
 1.0   1.0   ⋅

Проверить перестановку строк и столбцов можно с помощью

matrix.variables

3-element Vector{VariableRef}:
 x
 y
 z

matrix.affine_constraints

3-element Vector{ConstraintRef}:
 c1 : 6 x + 8 y ≥ 100
 c2 : 7 x + 12 y ≥ 120
 c3 : x + y ≤ 20

Построить столбец ослабляющих переменных, связанных с c3, можно так:

s_column = zeros(size(matrix.A, 1))
s_column[3] = 1.0

1.0

Таким образом, полная базисная матрица имеет следующий вид:

B = hcat(matrix.A[:, [1, 2]], s_column)

3×3 SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 6.0   8.0   ⋅
 7.0  12.0   ⋅
 1.0   1.0  1.0

lp_matrix_data возвращает отдельные векторы для нижней и верхней границ строк. Преобразуем их в один вектор правой части, взяв конечные элементы:

b = ifelse.(isfinite.(matrix.b_lower), matrix.b_lower, matrix.b_upper)

3-element Vector{Float64}:
 100.0
 120.0
  20.0

Если решить базисную систему так же, как и ранее, получится следующее.

B \ b

3-element Vector{Float64}:
 14.999999999999995
  1.250000000000004
  3.75

Это значение x, y и ослабляющей переменной, связанной с c3.

Выявление вырожденных переменных

Еще одна стандартная задача — выявление вырожденных переменных. Вырожденная переменная — это базисная переменная с оптимальным значением на нижней или верхней границе.

Вот функция, которая определяет, является ли переменная вырожденной:

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

is_degenerate (generic function with 1 method)

Простой пример линейной программы с вырожденным решением:

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))

1-element Vector{VariableRef}:
 x[1]

Решение является вырожденным, потому что:

value(x[1])

-0.0

get_attribute(x[1], MOI.VariableBasisStatus())

BASIC::BasisStatusCode = 0