Nonlinear optimization using JuMP, The solver does not support nonlinear problems (i.e., NLobjective and NLconstraint)

Specific Domains Optimization (Mathematical)
, ,
1

I am solving a nonlinear optimization problem using Mosek, I am running the following script

model = Model(Mosek.Optimizer)
@variable(model, y[1:n])
@NLconstraint(model, sum((Ay-b)[i]^4 for i=1:m) <= 1.0) # The constraint is norm(Ay-b,4) <=1
@objective(model, Min, 0.5*y’Qy + c’*y-1)
optimize!(model)

Then I got

The solver does not support nonlinear problems (i.e., NLobjective and NLconstraint).

When I do the L-2 norm constraint norm(Ay-b,2)<=1, I use
@constraint(model, [1; Ay-b] in SecondOrderCone()),

everything works well. But now for L-4 norm, am I doing something wrong?

1 Like
2

Mosek does not support arbitrary nonlinear constraints. You need to formulate your problem as a conic optimization problem.

See

In this particular case, the p-norm can be modeled using PowerCones: 4 The power cones — MOSEK Modeling Cookbook 3.3.0

julia> using JuMP

julia> import LinearAlgebra

julia> import SCS

julia> function pnorm_model(x0, p)
           model = Model(SCS.Optimizer)
           set_silent(model)
           @variable(model, x[i = 1:4] == x0[i])
           @variable(model, r[1:4])
           @variable(model, t)
           @constraint(model, [i=1:4], [r[i], t, x[i]] in MOI.PowerCone(1 / p))
           @constraint(model, sum(r) == t)
           @objective(model, Min, t)
           optimize!(model)
           return value(t)
       end
pnorm_model (generic function with 1 method)

julia> x0 = rand(4)
4-element Vector{Float64}:
 0.5033115464997369
 0.748802243142874
 0.0906063454411572
 0.44528993571248954

julia> pnorm_model(x0, 4)
0.8040190643234455

julia> LinearAlgebra.norm(x0, 4)
0.8040443072979447
1 Like
3

I’ve opened a PR to add this example to the documentation: [docs] add p-norm example by odow · Pull Request #3282 · jump-dev/JuMP.jl · GitHub

1 Like
4

Thank you so much! It works well for my problem!

1 Like
5

For future readers, the full problem with the norm constraint is something like:

julia> using JuMP

julia> import LinearAlgebra

julia> import SCS

julia> begin
           n = 3
           A = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
           b = [13, 14, 15, 16]
           p = 4
           model = Model(SCS.Optimizer)
           set_silent(model)
           @variable(model, y[1:n])
           @expression(model, f, A * y - b)
           @variable(model, r[1:length(f)])
           @variable(model, t <= 1)
           @constraint(model, [i=1:4], [r[i], t, f[i]] in MOI.PowerCone(1 / p))
           @constraint(model, sum(r) == t)
           @objective(model, Max, sum(y))
           optimize!(model)
           println("t = ", value(t))
           println("pnorm(Ay - b, $p) = ", LinearAlgebra.norm(A * value.(y) - b, p))
       end
t = 1.0002567650485095
pnorm(Ay - b, 4) = 1.0004631209883683

So t is the p-norm, give or take a little tolerance.

3 Likes
6

One more question, if I use Gurobi instead of Mosek, I got

Constraints of type MathOptInterface.VectorAffineFunction{Float64}-in-MathOptInterface.PowerCone{Float64} are not supported by the solver.

If you expected the solver to support your problem, you may have an error in your formulation. Otherwise, consider using a different solver.

Does MOI.PowerCone not work with Gurobi?

7

Nope. Gurobi is limited to linear or quadratic programs.

1 Like
8

But with Convex.jl, Gurobi is able to solve this 4-norm constraint problem. With JuMP, it isn’t. Is there a reason why Convex.jl works?

9

Convex must use a different reformulation that doesn’t use the PowerCone, but I don’t know the details.

I think this is really asking for an explicit PNormCone in JuMP. Thatd let us write proper reformulations for each solver.

10

For convex I just use the following code

y = Variable(n)
Prob = minimize(0.5quadform(y,Q; assume_psd=true)+dot(c,y)-1, norm(Ay - b, 4) <= 1)
Convex.solve!(Prob, Gurobi.Optimizer)

I don’t know why this works.

11

I don’t know why this works.

Because Convex includes a reformulation that JuMP has not implemented yet.

I’ve opened an issue to get this implemented: Add NormPCone · Issue #2118 · jump-dev/MathOptInterface.jl · GitHub

1 Like