Solvers, nonlinear constraints and user defined functions

Specific Domains Optimization (Mathematical)
1

I apologize if someone had raised questions like this. And I understand the topic sounded similar to some previous topics. I did do a search. But didn’t found specific answers to all the questions I have.

I understand that Gurobi is not a nonlinear programming solver. But it’s one of the lead solvers these days. Is it worth to go through all the trouble trying to reformulate the problem somehow into a convex quadratic problem approximately or a nonlinear solver such as Ipopt works as well?

Another question is on user defined functions in nonlinear constraints. I read from a post that the function will need to be registered and due to JuMP syntax, the augments have to be scalars, not vectors. I understand that this is only due to JuMP syntax. Therefore, a wrapper would be ok. So, I assume if I have something like

function myP(kk, n, x, alpha, Q, w, S, T)

and a wrapper function
function outP(kk, x…)

I could do something like
@variable(model, x[1:n]>=0.0)
JuMP.register(model, :outP, 5, outP, autodiff=true)
@variable(model, y>=0.0)
for i=1:K
@NLconstraint(model, y>=outP(i, x[1], x[2], x[3], x[4]))
end

So, I gave this a try. Here are the errors I saw. The first augment in outP and myP somehow became float, not the integer supposedly passed by i. Thus, there was a index error (1.0). If I force the function definition into

function myP(kk::Int, n, x, alpha, Q, w, S, T)
function outP(kk::Int, x…)

I got rid of the error in index. However, when it continues to
@objective(model, Min, y)
status = solve(model)

I got a method error - ERROR: MethodError: no method matching outP(::Float64, ::Float64, ::Float64, ::Float64, ::Float64)

Is there some default in register that make all the augments Float64? If yes, is there a way to define it explicitly?

Of course, this is just work around for a small example. What if x is of high dimension? Any way to get around the augment restriction or some other clever ways to handle nonlinear constraints?

Thanks!

2

I don’t think there’s any general consensus on this. Best is to try and see.

Yes, all function arguments are assumed to be things that you want to differentiate with respect to. If they’re actually parameters you should use a closure:

JuMP.register(model, (x1,x2,x3,x4) -> outP(1,x1,x2,x3,x4), :outP1, 4)
JuMP.register(model, (x1,x2,x3,x4) -> outP(2,x1,x2,x3,x4), :outP2, 4)

Julia+JuMP: variable number of arguments to function - Stack Overflow. Please read my caveat there about autodiff=true in high dimensions.

3

These are very helpful. I think that you meant a different order on the second and third augment, like

JuMP.register(model, :outP1, 4, (x1,x2,x3,x4) -> outP(i,x1,x2,x3,x4), autodiff=true)

And that worked.

1 Like
4

If you use e.g. NLopt directly, you can employ high-dimensional (e.g. I’ve used 10⁵-dimensional x) nonlinear constraints with arbitrary user-defined functions directly.

5

I think I understand registration and the need for closure. Since we need to iterate (the need for outP1, outP2, …), we need to use symbol. We’ve managed to have that work. However, splatting seems to be a problem. It didn’t work together with symbol. Here is the block of code.
for i=1:K
udfSymbol = Symbol(“udf”, i)
udf = (xargs…) → outP(i, xargs…)
JuMP.register(model, udfSymbol, n, udf, autodiff=true)
@NLconstraint(model, y >= $udfSymbol(x…))
end
If we listed out x, such as @NLconstraint(model, y >= $udfSymbol(x[1], x[2], x[3], x[4])), it works fine. But it will not be practical when x becomes high dimensional.
Is there a way to get around this in JuMP?

6

I’ve used this approach before.

Edit: which I see now is basically just what Miles proposed earlier.

7

Sorry, it seems that I only come back to look at this problem once every month this summer. But could someone point out to me what I did wrong with this simple piece of code?

arg_expressions = [:($(x[i])) for i=1:4]
JuMP.register(model, :outP, 4, outP, autodiff=true)
#@NLobjective(model, Min, outP(x[1], x[2], x[3], x[4]))
@NLobjective(model, Min, outP(:($arg_expressions…)))

The commented out code would work. But the expression wouldn’t.

Thanks,
Katherine

8

I’m really weak with expressions and nonlinear modeling in JuMP, so hopefully @miles.lubin can explain the right way to do this. In the mean time, you can try building the expression yourself and using JuMP.setNLobjective like so:

julia> function my_obj(x1, x2, x3, x4)
           return (x1-1)^2 + (x2-2)^2 + (x3-3)^2 + (x4-4)^2
       end
my_obj (generic function with 1 method)

julia> m = Model()
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver is default solver

julia> @variable(m, x[1:4])
4-element Array{JuMP.Variable,1}:
 x[1]
 x[2]
 x[3]
 x[4]

julia> JuMP.register(m, :my_obj, 4, my_obj, autodiff = true)

julia> obj_expr = Expr(:call, :my_obj, x...)
:(my_obj(x[1], x[2], x[3], x[4]))

julia> JuMP.setNLobjective(m, :Min, obj_expr)

julia> setsolver(m, IpoptSolver())

julia> solve(m)
....Ipopt output....
:Optimal

julia> getvalue(x)
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0

Also note that your args_expressions is exactly the same as x:

julia> [:($(x[i])) for i=1:4] == x
true

I’m not 100% sure why this @NLobjective(m, :Min, my_obj(x...)) doesn’t work, but I think it’s because with JuMP.register we told my_obj to expect 4 arguments, and instead we’re giving it one argument of type Expression:

julia> dump(:(x...))
Expr
  head: Symbol tuple
  args: Array{Any}((1,))
    1: Expr
      head: Symbol ...
      args: Array{Any}((1,))
        1: Symbol x
      typ: Any
  typ: Any
1 Like
9

Thanks! This helps a lot. I think your explanation of its being one argument with type expression makes sense. I had similar problem with @NLconstraint as well. I think it’s for the same reason.

10

Can you please tell me how did you install IPOPT solver? I am getting build issues with it. I do have the solver binaries. How do I use them and can that be generalized for other solver binaries too? I am using JuMP in JuliaPro.
Thanks

11

The Ipopt binaries are provided on OS X which is what I’m using. I assume that you are using Linux which will be compiled from source. Sorry, I don’t have experience with that. If you post it to the mailing list, I’m sure someone will have experience with it.