Modelling arbitrary non-linear expression (returned by a function) in JuMP

Specific Domains Optimization (Mathematical)
1

Hi, I’m new to Julia and JuMP so please forgive me if this is a trivial problem.

I have a constraint function of the form:

function g(x)
  [ 
      exp(x[1]) + exp(x[2]),
      x[1]^4 + x[2]^4
  ]
end

I’m trying to model g(x) .>= 0 using JuMP with no luck. Which is to say that this does not work:

@variable(m, x[1:2])
for i = 1:2
  @NLconstraint(m, g(x)[i] >= 0)
end

Whereas this does:

@NLconstraint(m, x[1]^3 + x[2]^3 >= 0)
@NLconstraint(m, exp(x[1]) >= 0)

Is there a way I could do this? I understand that JuMP has some restrictions that may present a challenge but I’m open to any reasonable implementations!

2

Unfortunately, you cannot do this very easily.

You can register a function which returns a scalar: Nonlinear Modeling · JuMP

There is also a hack for user-defined functions which return a vector:

Longer-term, we’re in the middle of a rewrite to JuMP’s nonlinear interface, https://github.com/jump-dev/JuMP.jl/pull/3106, which will make this work with no hacks or registered functions.

(p.s. I was having a conversation with @ccoffrin a few hours ago that this is the most common question I get about JuMP. It comes up a few times every week… So hopefully we’ll have a better answer once that PR is merged.)

3 Likes
3

For a similar problem and the suggested work-around, see Adding NL expressions to an AffExpr type Vector - #5 by odow

4

I ended up doing a very hackish implementation using Symbolics

using JuMP
using Ipopt
import Symbolics

@Symbolics.variables x[1:2]
x = Symbolics.scalarize(x)

function g(x)
    [
        x[1]^4 + x[2]^5,
        exp(x[1])
    ]
end

function o(x)
    exp(x[1]) + exp(x[2])
end

C = g(x)
O = o(x)

m = Model(Ipopt.Optimizer)

# replace variable
@variable(m, x[1:2])

for i in eachindex(C)
    eval("@NLconstraint(m, $(repr(C[i])) <= 0")
end

eval("@NLobjective(m, Min, $(repr(O)))")
optimize!(m)
println(value.(x))
1 Like
5

The use of eval isn’t very safe, and you’ll likely run into a bunch of scoping issues. A better approach is to use the raw expression input: Nonlinear Modeling · JuMP

using JuMP
import Symbolics

to_expression(f, x::Vector{VariableRef}) = f

function to_expression(f::Expr, x::Vector{VariableRef})
    if Meta.isexpr(f, :ref)
        return x[f.args[2]]
    end
    for i in 1:length(f.args)
        f.args[i] = to_expression(f.args[i], x)
    end
    return f
end

function to_expression(f::Function, x::Vector{VariableRef})
    Symbolics.@variables(y[1:length(x)])
    f_y = f(Symbolics.scalarize(y))
    if f_y isa Vector
        return [to_expression(Meta.parse("$f_i"), x) for f_i in f_y]
    else
        return to_expression(Meta.parse("$f_y"), x)
    end
end


g(x) = [x[1]^4 + x[2]^5, exp(x[1])]
o(x) = exp(x[1]) + exp(x[2])
model = Model()
@variable(model, x[1:2])
set_nonlinear_objective(model, MOI.MIN_SENSE, to_expression(o, x))
for gi in to_expression(g, x)
    add_nonlinear_constraint(model, :($gi <= 0))
end
print(model)
2 Likes