Specifying variables in add_constraint not working

[margaeor]

When calling the add_non_linear_constraint method, specifying the argument vars as seen in the examples doesn't seem to do anything. In order to make it work, we have to specify expr_vars instead. For instance, in the code snipped below, the first line with add_nonlinear_constraint doesn't work. Whereas the second (commented) line works.

N = 5
M = 1

m = JuMP.Model(with_optimizer(CPLEX_SILENT))
@variable(m, x[1:N+1])
@objective(m, Min, x[N+1])

# Initialize GlobalModel
gm = OCTHaGOn.GlobalModel(model = m, name = "qp")

Q = randn(N, N)
c = randn(N)
expr = :((x) -> -x'*$(Q)*x - $(c)'*x)

lbs = push!(-5*ones(N), -800)
ubs = push!(5*ones(N), 0)
lbs_dict = Dict(x .=> lbs)
ubs_dict = Dict(x .=> ubs)

OCTHaGOn.bound!(gm, Dict(var => [lbs_dict[var], ubs_dict[var]] for var in gm.vars))

# Add constraints
OCTHaGOn.add_nonlinear_constraint(gm, expr, vars=x[1:N])
#OCTHaGOn.add_nonlinear_constraint(gm, expr, vars=x[1:N], expr_vars=[x[1:N]])

OCTHaGOn.globalsolve!(gm)

[1ozturkbe]

Huh, interestingly, both of the input methods work for me, as well as the easiest one:

OCTHaGOn.add_nonlinear_constraint(gm, expr)

What do you mean by it not working? Does it throw an error? Please include when you have a chance.

[margaeor]

when you do globalsolve, it throws an error that there is no feasible point. The reason for this is that it tries to evaluate the constraint with all variables (whereas only the first N are needed). Hence, the evaluation returns -Inf and everything is infeasible

[1ozturkbe]

I know this issue... I have actually faced it before. Let me make a pull request about it. The fundamental problem is that the vector lengths need to be specified within the Expr:

expr = :((x) -> -x[1:$(N)]'*$(Q)*x[1:$(N)] - $(c)'*x[1:$(N)])

or else this implementation can be full of bugs. The reason is that, due to the nature of arrays, the function calls need to be made with the FULL x array as an input, not just the relevant indices. So if a function uses x[2], x[3], x[5], then x[1:5] must be its input at a minimum. If the user does not specify that x is not the full length of the variable (1:N+1), then the function call won't know that, using x[1:N+1]and cause an error. Does that make sense? A better way to make a bug-free a linear algebra function work is to flatten it or give it indices, like I showed above. This way, the function itself knows only to use the first N values.

[1ozturkbe]

I hope you can also now understand the design choice behind using Expr instead of functions in nonlinear function inputs. It creates flexibility, but can also cause some issues. In any case, I'm leaving this open before I make the relevant PR, making changes at least to documentation.

[margaeor]

I understand, but the commented line works (although it is vectorized). So I think we need to let the user explicitly state which variables are used through the vars argument, and if the user uses this argument, then we should respect that. Essentially, right now, the vars argument doesn't work but the expr_vars does. Additionally, we should describe the difference between the vars argument and the expr_vars argument (I think one of those has to go)

Last but not least, regarding the choice of using expressions, I think in the future it might much better to use JuMP expressions and work with those. Because this way the user will be able to interface the solver using the actual jump syntax, instead of using clunky Julia expressions that will be converted to JuMP expressions

[1ozturkbe]

Two things. vars and expr_vars do fundamentally different things. vars only contains a vector of JuMP.variables. expr_vars contains the vectorization of those variables. If you have a good way to implement an either/or input method, I would look forward to it! Make sure to check out the utility functions that I've written for going between vars and expr_vars beforehand. Also, the nonlinearize! function may be a good example to show how such vectorized functions may be flattened for a better UI. But you may note that not even JuMP allows for implicit vectorized operations such as those in linear algebra. They recommend using only scalar operations.

I may have told you in the past, but it is not possible to use JuMP.@NLExpressions. I wish it was. I tried and failed many times. The moment you define an @NLExpression, the whole JuMP model turns into a nonlinear model, making it impossible to use MILP solvers (they are considered invalid). And you cannot remove @NLconstraints or @NLexpressions once you have defined them, so you cannot replace the original nonlinear constraints with their MIO approximators... It makes using @NLExpressions impossible. You are welcome to try.

[1ozturkbe]

@margaeor do you have additional thoughts about this / ideas to improve the UI? It would be sweet to see something more streamlined!

[margaeor]

@1ozturkbe I still don't understand the reason why we need to be able to remove constraints. Can't we just copy the original JuMP model (before the MIO constraints are added)? e.g. look at this function https://jump.dev/JuMP.jl/stable/reference/models/#JuMP.copy_model

If we use something like that, can't we use NLExpressions and native JuMP syntax for constraints (instead of expressions)? And whenever we need a clear JuMP model (e.g. in desend), just copy it?

P.S. I don't think this is a priority right now, but it is useful for the shake of discussion.

[1ozturkbe]

Well, firstly, you need to be able to remove constraints if you want to refine/change your approximation at any point (either via a new type of ML model or by improving the one you are using), or if you want to do a localsolve. You could copy in practice, but that would mean that you would have to copy the whole model every time you ever wanted to change any single MIO constraint, which is super clunky and slow. In any case, that doesn't resolve the problem with @NLconstraint, since you'd have to remove those constraints first before making the MIO approximation.

AFAIK, you will can't use NLexpressions and native JuMP syntax, because even declaring NLexpression turns the model into a nonlinear model. Here is a MWE.

using JuMP, CPLEX
m = JuMP.Model(CPLEX.Optimizer)
@variable(m, x >= 0)
@objective(m, Min, x)
@NLexpression(m, sin(x))
optimize!(m)

You should get an error that the "NLPBlock" is not supported by Model, unfortunately.