Анализ чувствительности линейной программы

В этом руководстве объясняется, как использовать функцию lp_sensitivity_report для создания отчетов о чувствительности, подобных создаваемым Excel Solver. Чаще всего они используются на вводных занятиях по линейному программированию.

Если говорить коротко, анализ чувствительности линейной программы состоит из двух вопросов.

При данном оптимальном решении насколько могут измениться коэффициенты

целевой функции, прежде чем другое решение станет оптимальным?

При данном оптимальном решении насколько может измениться правая часть линейного

ограничения, прежде чем другое решение станет оптимальным?

В JuMP есть функция lp_sensitivity_report для вычисления этих значений, но в данном руководстве она расширяется для построения более информативных таблиц в виде DataFrame.

Настройка

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

using JuMP
import HiGHS
import DataFrames

а также следующая небольшая линейная программа:

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)

* Solver : HiGHS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "kHighsModelStatusOptimal"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 2.04000e+02
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : 2.04000e+02
  Primal solution :
    x : 1.50000e+01
    y : 1.25000e+00
    z : 1.00000e+00
  Dual solution :
    c1 : 2.50000e-01
    c2 : 1.50000e+00
    c3 : 0.00000e+00

* Work counters
  Solve time (sec)   : 3.29733e-04
  Simplex iterations : 2
  Barrier iterations : 0
  Node count         : -1

Можно ли определить:

целевой коэффициент каждой переменной?
правую часть каждого ограничения?
оптимальные прямое и двойственное решения?

Отчеты о чувствительности

Теперь вызовем метод lp_sensitivity_report:

report = lp_sensitivity_report(model)

SensitivityReport{Float64}(Dict{ConstraintRef, Tuple{Float64, Float64}}(y ≥ 0 => (-Inf, 1.25), z ≤ 1 => (-Inf, Inf), c1 : 6 x + 8 y ≥ 100 => (-4.0, 2.857142857142857), c3 : x + y ≤ 20 => (-3.75, Inf), x ≥ 0 => (-Inf, 15.0), c2 : 7 x + 12 y ≥ 120 => (-3.3333333333333335, 4.666666666666667), y ≤ 3 => (-1.75, Inf)), Dict{VariableRef, Tuple{Float64, Float64}}(y => (-4.0, 0.5714285714285714), x => (-0.3333333333333333, 3.0), z => (-Inf, 1.0)))

Он возвращает объект SensitivityReport, который сопоставляет:

каждую ссылку на переменную с кортежем (d_lo, d_hi)::Tuple{Float64,Float64}, объясняя, насколько может измениться целевой коэффициент соответствующей переменной так, чтобы исходный базис оставался оптимальным;
каждую ссылку на ограничение с кортежем (d_lo, d_hi)::Tuple{Float64,Float64}, объясняя, насколько может измениться правая часть соответствующего ограничения так, чтобы базис оставался оптимальным.

Оба кортежа являются относительными, а не абсолютными. Таким образом, если заданы целевой коэффициент 1.0 и кортеж (-0.5, 0.5), целевой коэффициент может находиться в диапазоне от 1.0 - 0.5 до 1.0 + 0.5.

Например:

report[x]

(-0.3333333333333333, 3.0)

указывает, что целевой коэффициент x, то есть 12, может уменьшиться на -0.333 или увеличиться на 3.0, при этом прямое решение (15, 1.25) останется оптимальным. Кроме того:

report[c1]

(-4.0, 2.857142857142857)

означает, что правая часть ограничения c1 (100) может уменьшиться на 4 единицы или увеличиться на 2,85 единицы, при этом прямое решение (15, 1.25) останется оптимальным.

Чувствительность переменных

Сами по себе кортежи неинформативны. Давайте рассмотрим их в контексте, сопоставив ряд других сведений о переменной:

function variable_report(xi)
    return (
        name = name(xi),
        lower_bound = has_lower_bound(xi) ? lower_bound(xi) : -Inf,
        value = value(xi),
        upper_bound = has_upper_bound(xi) ? upper_bound(xi) : Inf,
        reduced_cost = reduced_cost(xi),
        obj_coefficient = coefficient(objective_function(model), xi),
        allowed_decrease = report[xi][1],
        allowed_increase = report[xi][2],
    )
end

variable_report (generic function with 1 method)

Вызовем функцию для x:

x_report = variable_report(x)

(name = "x", lower_bound = 0.0, value = 15.0, upper_bound = Inf, reduced_cost = 0.0, obj_coefficient = 12.0, allowed_decrease = -0.3333333333333333, allowed_increase = 3.0)

Эти данные немного сложны для восприятия, поэтому давайте вызовем функцию для каждой переменной в модели и поместим результаты в DataFrame:

variable_df =
    DataFrames.DataFrame(variable_report(xi) for xi in all_variables(model))

3×8 DataFrame

Row	name	lower_bound	value	upper_bound	reduced_cost	obj_coefficient	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	x	0.0	15.0	Inf	0.0	12.0	-0.333333	3.0
2	y	0.0	1.25	3.0	0.0	20.0	-4.0	0.571429
3	z	-Inf	1.0	1.0	-1.0	-1.0	-Inf	1.0

Чувствительность ограничений

Нечто подобное можно сделать и с ограничениями:

function constraint_report(c::ConstraintRef)
    return (
        name = name(c),
        value = value(c),
        rhs = normalized_rhs(c),
        slack = normalized_rhs(c) - value(c),
        shadow_price = shadow_price(c),
        allowed_decrease = report[c][1],
        allowed_increase = report[c][2],
    )
end

c1_report = constraint_report(c1)

(name = "c1", value = 100.0, rhs = 100.0, slack = 0.0, shadow_price = -0.25, allowed_decrease = -4.0, allowed_increase = 2.857142857142857)

constraint_df = DataFrames.DataFrame(
    constraint_report(ci) for (F, S) in list_of_constraint_types(model) for
    ci in all_constraints(model, F, S) if F == AffExpr
)

3×7 DataFrame

Row	name	value	rhs	slack	shadow_price	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	c1	100.0	100.0	0.0	-0.25	-4.0	2.85714
2	c2	120.0	120.0	0.0	-1.5	-3.33333	4.66667
3	c3	16.25	20.0	3.75	0.0	-3.75	Inf

Вопросы для анализа

Теперь мы можем использовать эти объекты DataFrame, чтобы задавать вопросы о решении.

Например, можно найти базисные переменные, выполнив поиск переменных с приведенной стоимостью 0:

basic = filter(row -> iszero(row.reduced_cost), variable_df)

2×8 DataFrame

Row	name	lower_bound	value	upper_bound	reduced_cost	obj_coefficient	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	x	0.0	15.0	Inf	0.0	12.0	-0.333333	3.0
2	y	0.0	1.25	3.0	0.0	20.0	-4.0	0.571429

и небазисные переменные, выполнив поиск по ненулевой приведенной стоимости:

non_basic = filter(row -> !iszero(row.reduced_cost), variable_df)

1×8 DataFrame

Row	name	lower_bound	value	upper_bound	reduced_cost	obj_coefficient	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	z	-Inf	1.0	1.0	-1.0	-1.0	-Inf	1.0

Можно также найти связывающие ограничения путем поиска нулевых ослабляющих переменных:

binding = filter(row -> iszero(row.slack), constraint_df)

2×7 DataFrame

Row	name	value	rhs	slack	shadow_price	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	c1	100.0	100.0	0.0	-0.25	-4.0	2.85714
2	c2	120.0	120.0	0.0	-1.5	-3.33333	4.66667

или ненулевой теневой стоимости:

binding2 = filter(row -> !iszero(row.shadow_price), constraint_df)

2×7 DataFrame

Row	name	value	rhs	slack	shadow_price	allowed_decrease	allowed_increase
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	c1	100.0	100.0	0.0	-0.25	-4.0	2.85714
2	c2	120.0	120.0	0.0	-1.5	-3.33333	4.66667