Документация Engee

Отладка

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

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

Это более сложное руководство, чем другие руководства по началу работы. В разделе «Начало работы» содержится предварительный обзор отладки моделей JuMP. Однако если вы только приступаете к работе с JuMP, то можете бегло просмотреть это руководство и вернуться к нему, когда напишете несколько моделей JuMP.

 
julia> using JuMP


julia> import HiGHS

Получение справки

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

Прежде чем задавать новый вопрос, прочитайте запись Как упростить получение помощи с советами по правильной постановке вопроса.

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

Отладка кода Julia

Прочитайте главу, посвященную отладке, в книге ThinkJulia.jl. В ней содержится ряд полезных советов и рекомендаций по отладке кода Julia.

Сбои решения

Когда из-за какой-либо проблемы решатель не может найти оптимальное решение (или доказать, что его не существует), JuMP может вернуть один из статусов termination_status.

Например, если решатель нашел решение, но с численной неточностью, он может вернуть статус ALMOST_OPTIMAL или ALMOST_LOCALLY_SOLVED, указывающий на то, что задача была решена с ослабленным набором допусков. Помимо этого, решатель может вернуть такой статус, как NUMERICAL_ERROR, SLOW_PROGRESS или OTHER_ERROR, указывающий на то, что решение задачи найти не удалось.

Большинство решателей могут столкнуться с численной неточностью, так как используют арифметику с плавающей запятой для выполнения таких операций, как сложение, вычитание и умножение. Эти операции не являются точными, и между теоретическим значением и значением, вычисленным компьютером, может быть небольшая погрешность. Например:

julia> 0.1 * 3 == 0.3
false

Прочитайте раздел Рекомендации в отношении численных проблем в документации Gurobi, а также раздел Отладка численных проблем в документации YALMIP.

Распространенные причины

Распространенные причины сбоев решения:

  • Очень большие и очень маленькие коэффициенты. Что значит «большие», зависит от решателя и задачи, но, как правило, проблемы возникают при значениях больше 1e6 или меньше 1e-6.

  • Нелинейные задачи с функциями, которые не определены на некоторых участках области определения. Например, это может быть задача минимизации log(x), где x >= 0, которая не определена при x = 0 (общее начальное значение).

Стратегии

Стратегии для устранения причин сбоев решения:

  • Измените масштаб переменных в задаче и связанных коэффициентов так, чтобы все коэффициенты находились в диапазоне от 1e-4 до 1e4. Например, это может означать изменение масштаба переменной расстояния с сантиметров на километры.

  • Попробуйте другой решатель. Некоторые решатели могут оказаться более устойчивыми при решении определенной задачи, чем другие.

  • Прочитайте документацию по решателю и попробуйте настройки, которые способствуют численной устойчивости.

  • Задайте границы или добавьте ограничения таким образом, чтобы все нелинейные функции были определены во всей допустимой области. Это особенно касается таких функций, как 1 / x и log(x), которые не определены при x = 0.

Неверные результаты

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

Такую проблему бывает сложно обнаружить и отладить, так как решатель никак не сообщает о ней. Статус termination_status, скорее всего, будет OPTIMAL, и решение будет доступно.

Распространенные причины

Распространенные причины неверных результатов:

  • Ошибка моделирования, из-за которой модель JuMP не соответствует формулировке на бумаге.

  • Не учитываются допуски решателя (например, если x является двоичной переменной, значение x = 1.0000001 все равно может считаться допустимым).

  • Ошибка в JuMP или решателе.

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

Стратегии

Стратегии для устранения причин неверных результатов:

  • Выведите модель JuMP на экран, чтобы проверить, соответствует ли она формулировке на бумаге. Обращайте внимание на правильность знаков, например + вместо -, и ошибки смещения на единицу, например x[t] вместо x[t-1].

  • Убедитесь в том, что не используются точные сравнения, такие как value(x) == 1.0; всегда используйте isapprox(value(x), 1.0; atol = 1e-6), чтобы вручную указать допуск при сравнении.

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

Отладка недопустимой модели

Модель является недопустимой, если не существует прямого решения, удовлетворяющего всем ограничениям. Как правило, недопустимость модели означает одно из двух:

  • У задачи нет допустимого решения.

  • В модели есть ошибка.

Пример

Вот простой пример недопустимой модели:

julia> model = Model(HiGHS.Optimizer);


julia> set_silent(model)


julia> @variable(model, x >= 0)
x

julia> @objective(model, Max, 2x + 1)
2 x + 1

julia> @constraint(model, con, 2x - 1 <= -2)
con : 2 x ≤ -1

Задана граница x >= 0, однако ограничение можно переписать в виде x <= -1/2. Если задача недопустима, JuMP может вернуть один из нескольких статусов. Чаще всего возвращается INFEASIBLE:

julia> optimize!(model)


julia> termination_status(model)
INFEASIBLE::TerminationStatusCode = 2

В зависимости от решателя можно также получить статус INFEASIBLE_OR_UNBOUNDED или LOCALLY_INFEASIBLE.

Статус завершения INFEASIBLE_OR_UNBOUNDED означает, что решателю не удалось доказать недопустимость или неограниченность. Удалось доказать лишь то, что у модели нет конечного допустимого оптимального решения.

Нелинейные оптимизаторы, такие как Ipopt, могут возвращать статус LOCALLY_INFEASIBLE. Это означает не то, что решатель доказал отсутствие допустимого решения, а лишь то, что ему не удалось его найти. Если известна прямая допустимая точка, попробуйте указать ее в качестве начальной точки с помощью set_start_value и выполнить оптимизацию повторно.

Распространенные причины

Распространенные причины недопустимости:

  • Неверные единицы измерения, например, нижняя граница задана в мегаваттах, а верхняя — в киловаттах.

  • Использование знака + вместо - в ограничении.

  • Ошибки смещения на единицу и другие схожие ошибки, например использование x[t] вместо x[t-1] в части ограничения.

  • Недопустимость математической формулировки по иным причинам.

Стратегии

Стратегии для устранения причин недопустимости:

  • По очереди закомментируйте ограничения (или блоки ограничений) и еще раз решите задачу. Когда найдете ограничение, добавление которого делает задачу недопустимой, внимательно проверьте его на наличие ошибок.

  • Если задача по-прежнему недопустима, когда закомментированы все ограничения, проверьте границы переменных. Правильные ли данные используются?

  • Если известно допустимое решение, используйте primal_feasibility_report для вычисления ограничений и проверки их нарушения. Скорее всего, в одном из ограничений обнаружится опечатка.

  • Попробуйте другой решатель. Иногда из-за ошибки в решателе он может ошибочно сообщать о недопустимости задачи, хотя это не так. Если в вашем случае один решатель сообщает о недопустимости задачи, а другой находит оптимальное решение, сообщите об этом, открыв проблему в репозитории GitHub решателя, который сообщает о недопустимости.

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

 
julia> compute_conflict!(model)
ERROR: ArgumentError: The optimizer HiGHS.Optimizer does not support `compute_conflict!`

В данном случае HiGHS не поддерживает вычисление конфликтов, но другие решатели, такие как Gurobi и CPLEX, поддерживают такую возможность. Если решатель поддерживает вычисление конфликтов, дополнительные сведения см. в разделе Conflicts.

Ослабление со штрафом

Еще одна стратегия отладки причин недопустимости заключается в использовании функции relax_with_penalty!.

При ослаблении со штрафом ограничения вида преобразуются в форму , где , а затем в целевую функцию вводится штрафной член (при минимизации, в противном случае ), где  — это штраф.

julia> map = relax_with_penalty!(model)
┌ Warning: Skipping PenaltyRelaxation for ConstraintIndex{MathOptInterface.VariableIndex,MathOptInterface.GreaterThan{Float64}}
└ @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/jqDoD/src/Utilities/penalty_relaxation.jl:289
Dict{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}, AffExpr} with 1 entry:
  con : 2 x - _[2] ≤ -1 => _[2]

Здесь map — это словарь, в котором индексы ограничений сопоставляются с аффинным выражением, представляющим .

Если оптимизировать ослабленную модель, на этот раз мы получим допустимое решение:

julia> optimize!(model)


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

Переберите содержимое map, чтобы узнать, какие ограничения нарушаются:

julia> for (con, penalty) in map
           violation = value(penalty)
           if violation > 0
               println("Constraint `$(name(con))` is violated by $violation")
           end
       end
Constraint `con` is violated by 1.0

Обнаружив нарушаемое ограничение в ослабленной задаче, проверьте, нет ли в нем опечатки или другой распространенной ошибки.

Сведения о том, как изменить штрафной член a для каждого ограничения в модели или для определенного подмножества ограничений, см. в docstring функции relax_with_penalty!.

При использовании relax_with_penalty! следует учитывать следующее.

  • Границы переменных и ограничения целочисленности не ослабляются. Если задача по-прежнему недопустима после вызова relax_with_penalty!, проверьте границы переменных.

  • Отменить ослабление со штрафом невозможно. Если нужна первоначальная модель, постройте задачу повторно или вызовите copy_model перед вызовом relax_with_penalty!.

Отладка неограниченной модели

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

Пример

Вот простой пример неограниченной модели:

julia> model = Model(HiGHS.Optimizer);


julia> set_silent(model)


julia> @variable(model, x >= 0)
x

julia> @objective(model, Max, 2x + 1)
2 x + 1

x можно увеличивать бесконечно, и значение целевой функции 2x + 1 улучшается по мере увеличения x.

Если задача неограниченная, JuMP может вернуть один из нескольких статусов. Чаще всего возвращается DUAL_INFEASIBLE:

julia> optimize!(model)


julia> termination_status(model)
DUAL_INFEASIBLE::TerminationStatusCode = 3

В зависимости от решателя можно также получить статус INFEASIBLE_OR_UNBOUNDED или код ошибки, например NORM_LIMIT.

Распространенные причины

Распространенные причины неограниченности:

  • Использование Max вместо Min.

  • Опущение границ переменных, например 0 <= x <= 1.

  • Использование знака + вместо - в члене целевой функции.

Стратегии

Стратегии для устранения причин неограниченности:

  • Тщательно проверьте, правильно ли вы выбрали Min или Max в строке @objective.

  • Выведите целевую функцию с помощью print(objective_function(model)) и проверьте правильность значения и знака каждого коэффициента.

  • Добавьте большие границы ко всем свободным или ограниченным с одной стороны переменным, а затем решите задачу заново. Поскольку все переменные теперь ограничены, задача будет иметь конечное оптимальное решение. Проверьте значение каждой переменной в оптимальном решении: находится ли оно на одной из новых границ. Если да, вам либо необходимо указать более подходящую границу для этой переменной, либо в целевой функции может быть связанная с ней ошибка (например, знак + вместо -).

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

julia> model = Model(HiGHS.Optimizer);


julia> set_silent(model)


julia> @variable(model, x >= 0)
       # @objective(model, Max, 2x + 1)
x

julia> @variable(model, objective <= 10_000)
objective

julia> @constraint(model, objective <= 2x + 1)
-2 x + objective ≤ 1

julia> @objective(model, Max, objective)
objective

Эта новая модель имеет конечное оптимальное решение, поэтому ее можно решить, а затем найти в оптимальном решении переменные с большими положительными или отрицательными значениями.

julia> optimize!(model)


julia> @assert is_solved_and_feasible(model)


julia> for var in all_variables(model)
           if var == objective
               continue
           end
           if abs(value(var)) > 1e3
               println("Variable `$(name(var))` may be unbounded")
           end
       end
Variable `x` may be unbounded