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

Нелинейное моделирование (устаревшая информация)

На этой странице описывается устаревший нелинейный интерфейс JuMP. Он имеет ряд особенностей и ограничений, которые побудили разработать новый нелинейный интерфейс. Новый интерфейс описывается в разделе Nonlinear Modeling. Этот устаревший интерфейс сохранится во всех будущих выпусках v1.X JuMP. Два нелинейных интерфейса нельзя использовать совместно.

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

В JuMP внесены три основных изменения в решение нелинейных программ.

Существует ряд ограничений в отношении синтаксиса, который можно использовать в макросах @NLxxx. Обязательно прочтите раздел Syntax notes.

Задание нелинейной целевой функции

Чтобы задать нелинейную целевую функцию, используйте макрос @NLobjective.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @NLobjective(model, Min, exp(x[1]) - sqrt(x[2]))

Чтобы изменить нелинейную целевую функцию, еще раз вызовите @NLobjective.

Добавление нелинейного ограничения

Для добавления нелинейного ограничения используйте макрос @NLconstraint.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @NLconstraint(model, exp(x[1]) <= 1)
exp(x[1]) - 1.0 ≤ 0

julia> @NLconstraint(model, [i = 1:2], x[i]^i >= i)
2-element Vector{NonlinearConstraintRef{ScalarShape}}:
 x[1] ^ 1.0 - 1.0 ≥ 0
 x[2] ^ 2.0 - 2.0 ≥ 0

julia> @NLconstraint(model, con[i = 1:2], prod(x[j] for j = 1:i) == i)
2-element Vector{NonlinearConstraintRef{ScalarShape}}:
 (*)(x[1]) - 1.0 = 0
 x[1] * x[2] - 2.0 = 0

Нелинейные ограничения можно создавать только с использованием <=, >= и ==. Более общие ограничения типа «Nonlinear в Set» не поддерживаются.

Чтобы удалить нелинейное ограничение, используйте метод delete:

julia> delete(model, con[1])

Создание нелинейного выражения

Для создания объектов нелинейных выражений используйте макрос @NLexpression. Синтаксис такой же, как для макроса @expression, за тем исключением, что выражение может содержать нелинейные члены.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> expr = @NLexpression(model, exp(x[1]) + sqrt(x[2]))
subexpression[1]: exp(x[1]) + sqrt(x[2])

julia> my_anon_expr = @NLexpression(model, [i = 1:2], sin(x[i]))
2-element Vector{NonlinearExpression}:
 subexpression[2]: sin(x[1])
 subexpression[3]: sin(x[2])

julia> @NLexpression(model, my_expr[i = 1:2], sin(x[i]))
2-element Vector{NonlinearExpression}:
 subexpression[4]: sin(x[1])
 subexpression[5]: sin(x[2])

Нелинейное выражение можно использовать в @NLobjective, @NLconstraint и даже вкладывать в другие выражения @NLexpression.

julia> @NLobjective(model, Min, expr^2 + 1)

julia> @NLconstraint(model, [i = 1:2], my_expr[i] <= i)
2-element Vector{NonlinearConstraintRef{ScalarShape}}:
 subexpression[4] - 1.0 ≤ 0
 subexpression[5] - 2.0 ≤ 0

julia> @NLexpression(model, nested[i = 1:2], sin(my_expr[i]))
2-element Vector{NonlinearExpression}:
 subexpression[6]: sin(subexpression[4])
 subexpression[7]: sin(subexpression[5])

Чтобы запросить значение нелинейного выражения, используйте метод value:

julia> set_start_value(x[1], 1.0)

julia> value(start_value, nested[1])
0.7456241416655579

julia> sin(sin(1.0))
0.7456241416655579

Создание нелинейного параметра

Для нелинейных моделей в JuMP имеется синтаксис для создания явных объектов «параметров», представляющих собой константы в модели, которые можно эффективно обновлять между решениями.

Нелинейные параметры объявляются с помощью макроса @NLparameter и могут индексироваться посредством произвольных множеств аналогично переменным и выражениям JuMP.

Начальное значение параметра должно указываться справа от знака ==.

julia> model = Model();

julia> @variable(model, x);

julia> @NLparameter(model, p[i = 1:2] == i)
2-element Vector{NonlinearParameter}:
 parameter[1] == 1.0
 parameter[2] == 2.0

Анонимные параметры создаются с помощью именованного аргумента value:

julia> anon_parameter = @NLparameter(model, value = 1)
parameter[3] == 1.0

Параметр не является переменной оптимизации. Оно должен иметь фиксированное значение, присвоенное посредством ==. Если нужен параметр <= или >=, создайте вместо этого переменную с помощью макроса @variable.

Для запроса или изменения значения параметра используйте методы value и set_value.

julia> value.(p)
2-element Vector{Float64}:
 1.0
 2.0

julia> set_value(p[2], 3.0)
3.0

julia> value.(p)
2-element Vector{Float64}:
 1.0
 3.0

Нелинейные параметры должны использоваться только в нелинейных макросах.

Когда следует использовать параметры

Нелинейные параметры полезны при решении нелинейных моделей в следующем порядке.

using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, z)
@NLparameter(model, x == 1.0)
@NLobjective(model, Min, (z - x)^2)
optimize!(model)
@show value(z) # Равно 1.0.

# Теперь обновим значение x, чтобы решить другую задачу.
set_value(x, 5.0)
optimize!(model)
@show value(z) # Равно 5.0.
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

value(z) = 1.0
value(z) = 5.0

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

Замечания о синтаксисе

Синтаксис, принятый в нелинейных макросах, более ограничен по сравнению с синтаксисом для линейных и квадратичных макросов. Ниже отметим ряд важных моментов.

Только скалярные операции

За исключением рассматриваемого ниже синтаксиса распаковки, все выражения должны быть простыми скалярными операциями. Нельзя использовать dot, произведения матрицы и вектора, векторные срезы и т. д.

julia> model = Model();

julia> @variable(model, x[1:2]);

julia> @variable(model, y);

julia> c = [1, 2];

julia> @NLobjective(model, Min, c' * x + 3y)
ERROR: Unexpected array [1 2] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.
[...]

Преобразуйте векторные операции в явные операции sum():

julia> @NLobjective(model, Min, sum(c[i] * x[i] for i = 1:2) + 3y)

Либо используйте макрос @expression:

julia> @expression(model, expr, c' * x)
x[1] + 2 x[2]

julia> @NLobjective(model, Min, expr + 3y)

Распаковка

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

julia> model = Model();

julia> @variable(model, x[1:3]);

julia> @NLconstraint(model, *(x...) <= 1.0)
x[1] * x[2] * x[3] - 1.0 ≤ 0

julia> @NLconstraint(model, *((x / 2)...) <= 0.0)
ERROR: Unsupported use of the splatting operator. JuMP supports splatting only symbols. For example, `x...` is ok, but `(x + 1)...`, `[x; y]...` and `g(f(y)...)` are not.

Пользовательские функции

JuMP изначально поддерживает набор одномерных и многомерных функций, распознаваемых подмодулем MOI.Nonlinear. В дополнение к функциям из этого списка можно регистрировать собственные пользовательские нелинейные функции. Пользовательские функции можно использовать в любом месте в макросах @NLobjective, @NLconstraint и @NLexpression.

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

Пользовательские функции должны возвращать скалярный результат. Обходное решение см. в разделе User-defined operators with vector outputs.

Автоматическое дифференцирование

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

Автоматическое дифференцирование не является конечным. Автоматически вычисляемые JuMP производные не подвержены погрешности аппроксимации.

Для автоматического дифференцирования в JuMP применяется пакет ForwardDiff.jl. Инструкции по написанию функций, пригодных для автоматического дифференцирования, см. в документации по ForwardDiff.jl.

Распространенные ошибки при написании пользовательских функций

Получаете ошибку наподобие No method matching Float64(::ForwardDiff.Dual)? Прочитайте этот раздел, а также рекомендации в документации по ForwardDiff.jl.

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

f(x::Float64) = 2 * x  # Это не сработает.
f(x::Real)    = 2 * x  # Так правильно.
f(x)          = 2 * x  # Так тоже правильно.

Еще одна причина этой ошибки — создание внутри функции массивов типа Float64.

function bad_f(x...)
    y = zeros(length(x))  # Создается массив типа `Float64`!
    for i = 1:length(x)
        y[i] = x[i]^i
    end
    return sum(y)
end

function good_f(x::T...) where {T<:Real}
    y = zeros(T, length(x))  # Создайте вместо этого массив типа `T`!
    for i = 1:length(x)
        y[i] = x[i]^i
    end
    return sum(y)
end

Регистрация функции

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

square(x) = x^2
f(x, y) = (x - 1)^2 + (y - 2)^2

model = Model()

register(model, :square, 1, square; autodiff = true)
register(model, :my_f, 2, f; autodiff = true)

@variable(model, x[1:2] >= 0.5)
@NLobjective(model, Min, my_f(x[1], square(x[2])))

В приведенном выше коде создается модель JuMP с целевой функцией (x[1] - 1)^2 + (x[2]^2 - 2)^2. Метод register имеет следующие аргументы:

  1. Модель, для которой регистрируются функции.

  2. Объект символа Julia, который служит именем пользовательской функции

в выражениях JuMP.

  1. Количество принимаемых функцией входных аргументов.

  2. Метод Julia, вычисляющий функцию.

  3. Флаг, предписывающий JuMP автоматически вычислить точные градиенты.

Символ :my_f не обязательно должен совпадать с именем функции f. Однако так код будет более удобочитаемым. В макросах обязательно используйте my_f, а не f.

Пользовательские функции нельзя регистрировать повторно, и они не обновляются при изменении базовых функций Julia. Если необходимо изменить пользовательскую функцию между решениями, создайте модель повторно или используйте другое имя. Сведения о том, как изменить имя программным путем, см. в разделе Raw expression input.

Регистрация функции и градиента

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

Одномерные функции

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

f(x) = x^2
∇f(x) = 2x
model = Model()
register(model, :my_square, 1, f, ∇f; autodiff = true)
@variable(model, x >= 0)
@NLobjective(model, Min, my_square(x))

При autodiff = true JuMP использует автоматическое дифференцирование для вычисления гессиана.

Многомерные функции

Для многомерных функций функция градиента ∇f должна принимать градиентный вектор в качестве первого аргумента, который заполняется на месте:

f(x, y) = (x - 1)^2 + (y - 2)^2
function ∇f(g::AbstractVector{T}, x::T, y::T) where {T}
    g[1] = 2 * (x - 1)
    g[2] = 2 * (y - 2)
    return
end

model = Model()
register(model, :my_square, 2, f, ∇f)
@variable(model, x[1:2] >= 0)
@NLobjective(model, Min, my_square(x[1], x[2]))

Убедитесь в том, что первый аргумент функции ∇f поддерживает тип AbstractVector; не следует предполагать, что входные данные имеют тип Float64.

Регистрация функции, градиента и гессиана

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

Одномерные функции

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

f(x) = x^2
∇f(x) = 2x
∇²f(x) = 2
model = Model()
register(model, :my_square, 1, f, ∇f, ∇²f)
@variable(model, x >= 0)
@NLobjective(model, Min, my_square(x))

Многомерные функции

Для многомерных функций функция гессиана ∇²f должна принимать в качестве первого аргумента матрицу AbstractMatrix, нижний треугольник которой заполняется на месте:

f(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
function ∇f(g, x...)
    g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2
    g[2] = 200 * (x[2] - x[1]^2)
    return
end
function ∇²f(H, x...)
    H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2
    # H[1, 2] = -400 * x[1]  <-- Не требуется. Заполняется только нижний треугольник.
    H[2, 1] = -400 * x[1]
    H[2, 2] = 200.0
    return
end

model = Model()
register(model, :rosenbrock, 2, f, ∇f, ∇²f)
@variable(model, x[1:2])
@NLobjective(model, Min, rosenbrock(x[1], x[2]))

Предполагается, что гессианова матрица H инициализирована нулями, и поскольку матрица H симметрична, достаточно заполнить ненулевыми значениями элементы нижнего треугольника. Тип передаваемой матрицы H зависит от системы автоматического дифференцирования, поэтому первый аргумент функции гессиана должен поддерживать тип AbstractMatrix (он может быть отличен от Matrix{Float64}). Однако изначально предполагается только то, что H поддерживает size(H) и setindex!. Наконец, матрица считается плотной, поэтому для функций с входными данными высокой размерности производительность будет низкой.

Пользовательские функции с векторными входными данными

Пользовательские функции, которые принимают векторы в качестве входных аргументов (например, f(x::Vector)), не поддерживаются. Вместо этого используйте синтаксис распаковки Julia для создания функции со скалярными аргументами. Например, вместо

f(x::Vector) = sum(x[i]^i for i in 1:length(x))

определите:

f(x...) = sum(x[i]^i for i in 1:length(x))

Эту функцию f можно использовать в модели JuMP следующим образом:

model = Model()
@variable(model, x[1:5] >= 0)
f(x...) = sum(x[i]^i for i in 1:length(x))
register(model, :f, 5, f; autodiff = true)
@NLobjective(model, Min, f(x...))

Обязательно ознакомьтесь с ограничениями синтаксиса в разделе Splatting.

Факторы, влияющие на время решения

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

Total CPU secs in IPOPT (w/o function evaluations)   =      7.412
Total CPU secs in NLP function evaluations           =      2.083

Для Ipopt производительность можно повысить, установив дополнительные пакеты разреженной линейной алгебры; см. раздел Installation Guide. Советы по повышению производительности других решателей см. в соответствующей документации.

В свою очередь, время вычисления функций зависит от языка моделирования. JuMP вычисляет производные, применяя автоматическое дифференцирование в обратном режиме и методы раскраски графа для эффективного использования разреженности гессиановой матрицы. Соответствующую производительность JuMP можно осторожно оценить как в пять раз превышающую показатели AMPL. Дополнительные сведения см. в нашей статье в издании SIAM Review.

Запрос производных из модели JuMP

В некоторых сложных случаях может потребоваться напрямую запросить производные модели JuMP, а не поручать эту задачу решателю. В JuMP реализован внутренний интерфейс MOI.AbstractNLPEvaluator. Чтобы получить объект средства оценки NLP из модели JuMP, используйте NLPEvaluator](../api.md#JuMP.NLPEvaluator-Tuple{Model}). index возвращает индекс [MOI.VariableIndex, соответствующий переменной JuMP. Объект MOI.VariableIndex сам по себе является типобезопасной оболочкой для значения типа Int64 (хранящегося в поле .value).

Например:

julia> raw_index(v::MOI.VariableIndex) = v.value
raw_index (generic function with 1 method)

julia> model = Model();

julia> @variable(model, x)
x

julia> @variable(model, y)
y

julia> @NLobjective(model, Min, sin(x) + sin(y))

julia> values = zeros(2)
2-element Vector{Float64}:
 0.0
 0.0

julia> x_index = raw_index(JuMP.index(x))
1

julia> y_index = raw_index(JuMP.index(y))
2

julia> values[x_index] = 2.0
2.0

julia> values[y_index] = 3.0
3.0

julia> d = NLPEvaluator(model)
Nonlinear.Evaluator with available features:

  * :Grad

  * :Jac

  * :JacVec

  * :Hess

  * :HessVec

  * :ExprGraph

julia> MOI.initialize(d, [:Grad])

julia> MOI.eval_objective(d, values)
1.0504174348855488

julia> sin(2.0) + sin(3.0)
1.0504174348855488

julia> ∇f = zeros(2)
2-element Vector{Float64}:
 0.0
 0.0

julia> MOI.eval_objective_gradient(d, ∇f, values)

julia> ∇f[x_index], ∇f[y_index]
(-0.4161468365471424, -0.9899924966004454)

julia> cos(2.0), cos(3.0)
(-0.4161468365471424, -0.9899924966004454)

В области NLPEvaluator существуют только нелинейные ограничения (добавленные с помощью @NLconstraint) и нелинейные целевые функции (добавленные с помощью @NLobjective).

NLPEvaluator не вычисляет производные линейных и квадратичных ограничений и целевых функций.

При применении к объекту ссылки на нелинейное ограничение метод index](../api.md#JuMP.index-Tuple{ConstraintRef}) возвращает его индекс в виде [MOI.Nonlinear.ConstraintIndex. Например:

julia> model = Model();

julia> @variable(model, x);

julia> @NLconstraint(model, cons1, sin(x) <= 1);

julia> @NLconstraint(model, cons2, x + 5 == 10);

julia> typeof(cons1)
NonlinearConstraintRef{ScalarShape} (alias for ConstraintRef{GenericModel{Float64}, MathOptInterface.Nonlinear.ConstraintIndex, ScalarShape})

julia> index(cons1)
MathOptInterface.Nonlinear.ConstraintIndex(1)

julia> index(cons2)
MathOptInterface.Nonlinear.ConstraintIndex(2)

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

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

Ввод необработанного выражения

Этот раздел требует углубленного знания типа Expr в Julia. Сначала следует ознакомиться с разделом Выражения и вычисление документации Julia.

Помимо использования макросов @NLexpression, @NLobjective и @NLconstraint, можно также предоставлять объекты Julia Expr напрямую с помощью методов add_nonlinear_expression, set_nonlinear_objective и add_nonlinear_constraint.

Такой способ ввода может быть полезен, если выражения генерируются программным образом или возникают проблемы с компиляцией входных данных макроса (дополнительные сведения см. в разделе Known performance issues).

Добавление нелинейного выражения

Чтобы добавить нелинейное выражение в модель, используйте метод add_nonlinear_expression.

julia> model = Model();

julia> @variable(model, x)
x

julia> expr = :($(x) + sin($(x)^2))
:(x + sin(x ^ 2))

julia> expr_ref = add_nonlinear_expression(model, expr)
subexpression[1]: x + sin(x ^ 2.0)

Это равносильно следующему коду:

julia> model = Model();

julia> @variable(model, x);

julia> expr_ref = @NLexpression(model, x + sin(x^2))
subexpression[1]: x + sin(x ^ 2.0)

Переменные следует интерполировать напрямую в выражение expr.

Задание целевой функции

Чтобы задать нелинейную целевую функцию, используйте метод set_nonlinear_objective.

julia> model = Model();

julia> @variable(model, x);

julia> expr = :($(x) + $(x)^2)
:(x + x ^ 2)

julia> set_nonlinear_objective(model, MIN_SENSE, expr)

Это равносильно следующему коду:

julia> model = Model();

julia> @variable(model, x);

julia> @NLobjective(model, Min, x + x^2)

Вместо Min и Max следует использовать MIN_SENSE или MAX_SENSE.

Добавление ограничения

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

julia> model = Model();

julia> @variable(model, x);

julia> expr = :($(x) + $(x)^2)
:(x + x ^ 2)

julia> add_nonlinear_constraint(model, :($(expr) <= 1))
(x + x ^ 2.0) - 1.0 ≤ 0

Это равносильно следующему коду:

julia> model = Model();

julia> @variable(model, x);

julia> @NLconstraint(model, Min, x + x^2 <= 1)
(x + x ^ 2.0) - 1.0 ≤ 0

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

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

В качестве примера создадим модель с нелинейными ограничениями f(x) <= 1, где f(x) = x^2 и f(x) = sin(x)^2:

julia> function main(functions::Vector{Function})
           model = Model()
           @variable(model, x)
           for (i, f) in enumerate(functions)
               f_sym = Symbol("f_$(i)")
               register(model, f_sym, 1, f; autodiff = true)
               add_nonlinear_constraint(model, :($(f_sym)($(x)) <= 1))
           end
           print(model)
           return
       end
main (generic function with 1 method)

julia> main([x -> x^2, x -> sin(x)^2])
Feasibility
Subject to
 f_1(x) - 1.0 ≤ 0
 f_2(x) - 1.0 ≤ 0

В качестве еще одного примера создадим модель с ограничением x^2 + sin(x)^2 <= 1:

julia> function main(functions::Vector{Function})
           model = Model()
           @variable(model, x)
           expr = Expr(:call, :+)
           for (i, f) in enumerate(functions)
               f_sym = Symbol("f_$(i)")
               register(model, f_sym, 1, f; autodiff = true)
               push!(expr.args, :($(f_sym)($(x))))
           end
           add_nonlinear_constraint(model, :($(expr) <= 1))
           print(model)
           return
       end
main (generic function with 1 method)

julia> main([x -> x^2, x -> sin(x)^2])
Feasibility
Subject to
 (f_1(x) + f_2(x)) - 1.0 ≤ 0

Зарегистрированные функции с переменным числом аргументов

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

julia> f(x...) = sum(exp(x[i]^2) for i in 1:length(x));

с разным количеством аргументов.

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

julia> A = [[1], [1, 2], [2, 3, 4], [1, 3, 4, 5]];

julia> model = Model();

julia> @variable(model, x[1:5]);

julia> funcs = Set{Symbol}();

julia> for a in A
           key = Symbol("f$(length(a))")
           if !(key in funcs)
               push!(funcs, key)
               register(model, key, length(a), f; autodiff = true)
           end
           add_nonlinear_constraint(model, :($key($(x[a]...)) <= 1))
       end

julia> print(model)
Feasibility
Subject to
 f1(x[1]) - 1.0 ≤ 0
 f2(x[1], x[2]) - 1.0 ≤ 0
 f3(x[2], x[3], x[4]) - 1.0 ≤ 0
 f4(x[1], x[3], x[4], x[5]) - 1.0 ≤ 0

Известные проблемы с производительностью

Ввод данных на основе макросов в нелинейный интерфейс JuMP может вызвать проблемы с производительностью в следующих случаях:

  1. Макрос содержит большое количество (сотни) членов.

  2. Макрос вызывается из функции, а не из верхнего уровня

глобальной области.

Первая проблема обусловлена не итоговым количеством членов в математическом выражении, а количеством членов в представлении Expr этого выражения в Julia. Например, выражение sum(x[i] for i in 1:1_000_000) содержит миллион математических членов, однако представление Expr — это всего лишь одна сумма.

Наиболее распространенной причиной, помимо утомительного набора текста, является написание программы, которая автоматически записывает модель JuMP в текстовый файл, который вы затем запускаете. Одним из примеров является пакет MINLPlib.jl, который автоматически компилирует модели в скалярном формате GAMS в примеры на языке JuMP.

Согласно общему правилу, если вы пишете программы для автоматического генерирования выражений для макросов JuMP, следует придерживаться подхода, описанного в разделе Raw expression input. Дополнительные сведения см. в описании проблемы № 1997 в репозитории MathOptInterface.