Решения

В этом разделе руководства описывается, как получить доступ к готовому решению задачи. В качестве примера используется следующая модель:

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

Проверка существования оптимального решения

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

julia> is_solved_and_feasible(model)
true

По умолчанию is_solved_and_feasible возвращает значение true как для глобального, так и для локального оптимума. Чтобы проверить, нашел ли решатель глобально оптимальное решение, передайте аргумент allow_local = false:

julia> is_solved_and_feasible(model; allow_local = false)
true

Чтобы проверить, нашел ли решатель оптимальное двойственное решение, помимо оптимального прямого решения, передайте аргумент dual = true:

julia> is_solved_and_feasible(model; dual = true)
true

Если эта функция возвращает значение false, используйте указанные ниже функции, такие как solution_summary, termination_status, primal_status и dual_status, для определения характера найденного решения (если оно есть).

Сводка по решениям

С помощью метода solution_summary можно просмотреть сводку по решениям задачи оптимизации.

julia> solution_summary(model)
* 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.05143e+02
  Objective bound    : -0.00000e+00
  Relative gap       : Inf
  Dual objective value : -2.05143e+02

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

julia> 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.05143e+02
  Objective bound    : -0.00000e+00
  Relative gap       : Inf
  Dual objective value : -2.05143e+02
  Primal solution :
    x : 1.54286e+01
    y[a] : 1.00000e+00
    y[b] : 1.00000e+00
  Dual solution :
    c1 : 1.71429e+00

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

Почему остановился решатель?

Чтобы понять, почему остановился решатель, используйте метод termination_status.

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

В перечислении MOI.TerminationStatusCode содержится полный список возможных состояний.

Обычно возвращаются такие значения, как OPTIMAL, LOCALLY_SOLVED, INFEASIBLE, DUAL_INFEASIBLE и TIME_LIMIT.

Возвращаемое состояние OPTIMAL означает, что решатель нашел (и подтвердил) глобально оптимальное решение. Возвращаемое состояние LOCALLY_SOLVED означает, что решатель нашел локально оптимальное решение (которое также может быть глобально оптимальным, но доказать это не удалось).

Возвращаемое состояние DUAL_INFEASIBLE не гарантирует, что прямое решение является неограниченным. Если двойственное решение является недопустимым, прямое решение является неограниченным при условии существования допустимого прямого решения.

С помощью метода raw_status можно получить характерную для решателя строку с описанием причины остановки оптимизации:

julia> raw_status(model)
"kHighsModelStatusOptimal"

Прямые решения

Состояние прямого решения

Метод primal_status](../api.md#JuMP.primal_status-Tuple{GenericModel}) возвращает перечисление [MOI.ResultStatusCode, описывающее состояние прямого решения.

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

Другие часто возвращаемые значения — NO_SOLUTION и INFEASIBILITY_CERTIFICATE. Первое из них означает, что у решателя нет решения для возврата, а второе — что прямое решение представляет собой свидетельство о двойственной недопустимости (прямой неограниченный луч).

Можно также воспользоваться методом has_values, который возвращает значение true при наличии решения, которое можно запросить, и false в противном случае.

julia> has_values(model)
true

Целевые значения

Целевое значение решенной задачи можно получить с помощью метода objective_value. Лучшую из известных границ оптимального целевого значения можно получить с помощью метода objective_bound. Если решатель поддерживает такую возможность, значение двойственной целевой функции можно получить с помощью метода dual_objective_value.

julia> objective_value(model)
-205.14285714285714

julia> objective_bound(model)  # HiGHS реализует границу целевого значения только для программ MIP
-0.0

julia> dual_objective_value(model)
-205.1428571428571

Значения прямого решения

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

julia> value(x)
15.428571428571429

value транслируется через контейнеры:

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 также работает с выражениями:

julia> value(my_expr)
100.57142857142857

и ограничениями:

julia> value(c1)
120.0

При вызове value для ограничения возвращается функция ограничения, вычисленная в точке решения.

Двойственные решения

Состояние двойственного решения

Метод dual_status](../api.md#JuMP.dual_status-Tuple{GenericModel}) возвращает перечисление [MOI.ResultStatusCode, описывающее состояние двойственного решения.

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

Другие часто возвращаемые значения — NO_SOLUTION и INFEASIBILITY_CERTIFICATE. Первое из них означает, что у решателя нет решения для возврата, а второе — что двойственное решение представляет собой свидетельство о прямой недопустимости (двойственный неограниченный луч).

Можно также воспользоваться методом has_duals, который возвращает значение true при наличии решения, которое можно запросить, и false в противном случае.

julia> has_duals(model)
true

Значения двойственного решения

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

julia> dual(c1)
1.7142857142857142

Запросить двойственные значения границ переменной можно с помощью LowerBoundRef, UpperBoundRef и 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

Определение двойственности в JuMP не зависит от целевого назначения. То есть знак допустимых двойственных решений, связанных с ограничением, зависит от направления ограничения, а не от того, является ли задача задачей максимизации или минимизации. Это соглашение отличается от определения двойственности линейного программирования в некоторых популярных учебниках. Если в линейной программе требуется определение из учебника, вероятно, лучше будет воспользоваться методами shadow_price и reduced_cost.

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

Ошибки OptimizeNotCalled

Из-за различий во внутренних механизмах кэширования решений решателями изменение модели после вызова optimize! приводит к ее сбросу в состояние OPTIMIZE_NOT_CALLED. Если затем попытаться запросить информацию о решениях, будет выдана ошибка OptimizeNotCalled.

Если вы итеративно запрашиваете информацию о решениях и изменяете модель, сначала запросите все результаты, а затем измените задачу.

Например, вместо следующего кода:

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 ~/work/JuMP.jl/JuMP.jl/src/optimizer_interface.jl:712
ERROR: OptimizeNotCalled()
Stacktrace:
[...]

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

выполните такой:

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

Если известно, что определенный решатель поддерживает запрос информации о решениях после изменений, можно воспользоваться методом direct_model, чтобы обойти состояние OPTIMIZE_NOT_CALLED:

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

При этом будьте осторожны. Если решатель не поддерживает запрос информации о решениях после изменения, он может вернуть неправильные решения без предупреждения или выдать ошибку.

Доступ к атрибутам

В интерфейсе MathOptInterface определено множество атрибутов модели, которые можно запрашивать. К некоторым атрибутам можно обращаться напрямую с помощью функций получения значения. К ним относятся следующие:

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

Если дана некоторая задача линейного программирования (LP) и оптимальное решение, соответствующее базису, может возникнуть вопрос, насколько может измениться коэффициент целевой функции или коэффициент в правой части стандартной формы (см. также normalized_rhs) без нарушения прямой или двойственной выполнимости базисного решения.

Обратите внимание, что не все решатели вычисляют базис, а для анализа чувствительности интерфейс решателя должен реализовывать MOI.ConstraintBasisStatus.

Дополнительные сведения об анализе чувствительности см. в разделе Sensitivity analysis of a linear program.

В качестве простого примера можно проанализировать чувствительность оптимального решения для следующей (невырожденной) задачи LP:

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

Для анализа чувствительности задачи мы могли бы проверить допустимые диапазоны возмущений, например коэффициентов стоимости и коэффициента в правой части ограничения c1, следующим образом:

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.

Связанный с переменной диапазон представляет собой диапазон допустимого возмущения соответствующего коэффициента целевой функции. Обратите внимание, что текущее прямое решение остается оптимальным в этом диапазоне; однако соответствующее двойственное решение может измениться из-за возмущения коэффициента стоимости. Аналогичным образом, связанный с ограничением диапазон представляет собой диапазон допустимого возмущения соответствующего коэффициента в правой части. В этом диапазоне текущее двойственное решение остается оптимальным, но оптимальное прямое решение может измениться.

Если задача вырожденная, то существует несколько оптимальных базисов и, следовательно, эти диапазоны могут быть не столь интуитивно понятными и казаться слишком узкими. Например, более значительное возмущение коэффициента стоимости может не делать недействительной оптимальность текущего прямого решения. Более того, если задача является вырожденной из-за конечной точности, то, например, возмущение может приводить к кажущейся недействительности базиса, хотя это не так (опять же из-за слишком узких диапазонов). Чтобы этого не происходило, увеличьте значение именованного аргумента atol метода lp_sensitivity_report. Обратите внимание, что из-за этого диапазоны могут стать слишком широкими для сложных с численной точки зрения случаев. Поэтому этим диапазонам не стоит слепо доверять, особенно в сильно вырожденных или численно нестабильных случаях.

Конфликты

Если входная модель недопустима, некоторые решатели могут помочь установить причину этого, предложив конфликт, то есть подмножество ограничений, вызывающих недопустимость. В зависимости от решателя конфликт также может называться IIS (неприводимой несогласованной подсистемой).

Если такая возможность поддерживается решателем, используйте метод compute_conflict!](../api.md#JuMP.compute_conflict!-Tuple{GenericModel}), чтобы запустить вычисление конфликта. По завершении этого процесса запросите атрибут [MOI.ConflictStatus, чтобы проверить, был ли обнаружен конфликт.

Если подсистема IIS найдена, скопируйте ее в новую модель с помощью метода copy_conflict. Затем ее можно будет вывести на экран или записать в файл для упрощения отладки:

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

Если требуется более полный контроль над списком ограничений, участвующих в конфликте, выполните итерацию по этому списку и запросите атрибут MOI.ConstraintConflictStatus:

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

Несколько решений

Некоторые решатели могут возвращать несколько решений. Проверить, сколько решений доступно для запроса, можно с помощью метода result_count.

Функции для запроса решений, например primal_status, dual_status, value, dual и solution_summary, принимают дополнительный именованный аргумент result, с помощью которого можно указать, какой результат следует вернуть.

Даже если termination_status имеет значение OPTIMAL, некоторые возвращаемые решения могут быть неоптимальными. Однако если решатель нашел хотя бы одно оптимальное решение, то при result = 1 всегда возвращается оптимальное решение. Для оценки качества остальных решений используйте метод objective_value.

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)
* Solver : MOA[algorithm=MultiObjectiveAlgorithms.Dichotomy, optimizer=HiGHS]

* Status
  Result count       : 3
  Termination status : OPTIMAL
  Message from the solver:
  "Solve complete. Found 3 solution(s)"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : [0.00000e+00,0.00000e+00]
  Objective bound    : [0.00000e+00,-9.00000e+00]
  Relative gap       : Inf
  Dual objective value : -9.00000e+00

* Work counters
  Solve time (sec)   : 1.82720e-03
  Simplex iterations : 0
  Barrier iterations : 0
  Node count         : -1

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]

В руководстве Multi-objective knapsack приводятся дополнительные примеры запроса нескольких решений.

Проверка допустимости решений

Чтобы проверить допустимость прямого решения, используйте функцию primal_feasibility_report, которая принимает model, словарь, в котором каждая переменная сопоставляется со значением прямого решения (по умолчанию — последнего полученного решения), и допуск atol (по умолчанию 0.0).

Эта функция возвращает словарь, в котором ссылки на недопустимые ограничения сопоставляются с расстоянием между прямым значением ограничения и ближайшей точкой в соответствующем множестве. Точка классифицируется как недопустимая, если расстояние превышает указанный допуск 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

Если точка допустима, возвращается пустой словарь:

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

Чтобы использовать прямое решение из решателя, опустите аргумент point:

julia> optimize!(model)

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

Вызов primal_feasibility_report без аргумента point полезен в случае, когда primal_status имеет значение FEASIBLE_POINT или NEARLY_FEASIBLE_POINT и нужно оценить качество решения.

Чтобы применить primal_feasibility_report](../api.md#JuMP.primal_feasibility_report-Union{Tuple{GenericModel{T}}, Tuple{T}, Tuple{GenericModel{T}, AbstractDict{GenericVariableRef{T}, T}}} where T) к недопустимой модели, необходимо также предоставить потенциальную точку (решатели обычно ее не предоставляют). Сведения о диагностике причины недопустимости см. в разделе [Conflicts.

Чтобы пропустить ограничения, содержащие переменные, отсутствующие в point, передайте аргумент skip_mising = true.

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

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

julia> optimize!(model)

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

Решения

Проверка существования оптимального решения

Сводка по решениям

Почему остановился решатель?

Прямые решения

Состояние прямого решения

Целевые значения

Значения прямого решения

Двойственные решения

Состояние двойственного решения

Значения двойственного решения

Рекомендованный рабочий процесс

Ошибки OptimizeNotCalled

Доступ к атрибутам

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

Конфликты

Несколько решений

Проверка допустимости решений