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

Языки алгебраического моделирования

JuMP — это язык алгебраического моделирования для математической оптимизации, написанный на языке Julia. На этой странице объясняется, что на самом деле представляет собой язык алгебраического моделирования.

Что такое язык алгебраического моделирования

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

где c, A и b — это векторы и матрицы данных соответствующего размера, а обозначает набор целочисленных переменных.

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

Решатель — это программный пакет, который вычисляет решения для одного или нескольких классов задач.

 
Например, https://github.com/ERGO-Code/HiGHS[HiGHS] — это решатель задач

линейного программирования (LP) и частично целочисленного программирования (MIP). Он включает в себя такие алгоритмы, как симплексный метод и метод внутренней точки.

В настоящее время JuMP поддерживает ряд решателей с открытым исходным кодом и коммерческих решателей,

которые можно просмотреть в таблице Supported-solvers.

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

Выглядит знакомо? Это задача о рюкзаке.

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

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

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

Язык алгебраического моделирования — это инструмент, который упрощает преобразование между алгебраической формой конструктора моделей и стандартной формой решателя.

Каждый язык алгебраического моделирования состоит из двух основных частей:

  1. предметно-ориентированный язык, позволяющий пользователю записывать задачи в алгебраической

форме;

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

решателем (и обратно).

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

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

Что такое MathOptInterface

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

MathOptInterface состоит из трех основных частей:

  1. Независимый от решателя API, который абстрагирует такие концепции, как добавление и удаление

переменных и ограничений, задание и получение параметров, а также запрос результатов. Дополнительные сведения об API MathOptInterface см. в документации.

  1. Автоматическая система преобразования, основанная на эквивалентных формулировках

ограничения. Дополнительные сведения об этой системе преобразования см. в разделе LazyBridgeOptimizer руководства и в нашей статье на сайте arXiv.

  1. Вспомогательные средства для управления тем, как и когда модели копируются в решатели. Дополнительные

сведения об этом см. в разделе CachingOptimizer руководства.

От пользователя к решателю

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

Шаг I. Запись в алгебраической форме

JuMP предоставляет первую часть языка алгебраического моделирования посредством макросов @variable, @objective и @constraint.

Например, задача о рюкзаке записывается в JuMP так:

julia> using JuMP, HiGHS

julia> function algebraic_knapsack(c, w, b)
           n = length(c)
           model = Model(HiGHS.Optimizer)
           set_silent(model)
           @variable(model, x[1:n] >= 0, Int)
           @objective(model, Max, sum(c[i] * x[i] for i = 1:n))
           @constraint(model, sum(w[i] * x[i] for i = 1:n) <= b)
           optimize!(model)
           if termination_status(model) != OPTIMAL
               error("Not solved correctly")
           end
           return value.(x)
       end
algebraic_knapsack (generic function with 1 method)

julia> algebraic_knapsack([1, 2], [0.5, 0.5], 1.25)
2-element Vector{Float64}:
 0.0
 2.0

Эта формулировка лаконична и близка к алгебраической формулировке модели, написанной выше.

Шаг II. Преобразование алгебраической формы в функциональную

На следующем этапе макросы JuMP переписывают переменные и ограничения в функциональной форме. Вот как выглядит код JuMP после этого шага.

julia> using JuMP, HiGHS

julia> function nonalgebraic_knapsack(c, w, b)
           n = length(c)
           model = Model(HiGHS.Optimizer)
           set_silent(model)
           x = [VariableRef(model) for i = 1:n]
           for i = 1:n
               set_lower_bound(x[i], 0)
               set_integer(x[i])
               set_name(x[i], "x[$i]")
           end
           obj = AffExpr(0.0)
           for i = 1:n
               add_to_expression!(obj, c[i], x[i])
           end
           set_objective(model, MAX_SENSE, obj)
           lhs = AffExpr(0.0)
           for i = 1:n
               add_to_expression!(lhs, w[i], x[i])
           end
           con = build_constraint(error, lhs, MOI.LessThan(b))
           add_constraint(model, con)
           optimize!(model)
           if termination_status(model) != OPTIMAL
               error("Not solved correctly")
           end
           return value.(x)
       end
nonalgebraic_knapsack (generic function with 1 method)

julia> nonalgebraic_knapsack([1, 2], [0.5, 0.5], 1.25)
2-element Vector{Float64}:
 0.0
 2.0

Думаем, вы согласитесь, что версию с макросами гораздо легче читать.

Часть III. Из JuMP в MathOptInterface

На третьем этапе JuMP преобразует функциональную форму задачи, то есть nonalgebraic_knapsack, в API MathOptInterface:

julia> import MathOptInterface as MOI

julia> import HiGHS

julia> function mathoptinterface_knapsack(optimizer, c, w, b)
           n = length(c)
           model = MOI.instantiate(optimizer)
           MOI.set(model, MOI.Silent(), true)
           x = MOI.add_variables(model, n)
           for i in 1:n
               MOI.add_constraint(model, x[i], MOI.GreaterThan(0.0))
               MOI.add_constraint(model, x[i], MOI.Integer())
               MOI.set(model, MOI.VariableName(), x[i], "x[$i]")
           end
           MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
           obj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0)
           MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)
           MOI.add_constraint(
               model,
               MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(w, x), 0.0),
               MOI.LessThan(b),
           )
           MOI.optimize!(model)
           if MOI.get(model, MOI.TerminationStatus()) != MOI.OPTIMAL
               error("Not solved correctly")
           end
           return MOI.get.(model, MOI.VariablePrimal(), x)
       end
mathoptinterface_knapsack (generic function with 1 method)

julia> mathoptinterface_knapsack(HiGHS.Optimizer, [1.0, 2.0], [0.5, 0.5], 1.25)
2-element Vector{Float64}:
 0.0
 2.0

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

Шаг IV. Из MathOptInterface в HiGHS

На последнем этапе пакет HiGHS.jl преобразует форму MathOptInterface, то есть mathoptinterface_knapsack, в API HiGHS:

julia> using HiGHS

julia> function highs_knapsack(c, w, b)
           n = length(c)
           model = Highs_create()
           Highs_setBoolOptionValue(model, "output_flag", false)
           for i in 1:n
               Highs_addCol(model, c[i], 0.0, Inf, 0, C_NULL, C_NULL)
               Highs_changeColIntegrality(model, i-1, 1)
           end
           Highs_changeObjectiveSense(model, -1)
           Highs_addRow(
               model,
               -Inf,
               b,
               Cint(length(w)),
               collect(Cint(0):Cint(n-1)),
               w,
           )
           Highs_run(model)
           if Highs_getModelStatus(model) != kHighsModelStatusOptimal
               error("Not solved correctly")
           end
           x = fill(NaN, 2)
           Highs_getSolution(model, x, C_NULL, C_NULL, C_NULL)
           Highs_destroy(model)
           return x
       end
highs_knapsack (generic function with 1 method)

julia> highs_knapsack([1.0, 2.0], [0.5, 0.5], 1.25)
2-element Vector{Float64}:
 0.0
 2.0

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

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