Constraint seems linear but the compiler determined it was quadratic. Not sure why or how to fix

Specific Domains Optimization (Mathematical)
1

Hi all,

I’m doing some optimization on weighted graphs, and one of my constraints is showing up as quadratic for some reason, and JuMP says it’s can’t convert it to Affine, even though the docs say this is possible. (“Because solvers can take advantage of the knowledge that a constraint is quadratic, prefer adding quadratic constraints using constraint, rather than NLconstraint.” Linked here: ) Additionally, I’m not sure why JuMP thinks it’s quadratic; I’m not multiplying or dividing my two variables, just adding an inequality constraint between them.

The model in question is as follows, with some missing constraints for brevity. The only one that fails is the one surrounded by asterisks; the rest have been tested and work fine from what I can tell.

model=Model(HiGHS.Optimizer)
flows=collect(Graphs.weights(graph))
@variable(model, ACR[1:nv(graph),1:nv(graph)], Int);
@variable(model, an[1:nv(graph)], Bin);
@constraint(model,  ACR .>= 0)
**@constraint(model,  (ACR .<= Diagonal(1 .- an)*flows))** 
@constraint(model,  sum(ACR)+dot(cn,an) <= B);
etc...

What I don’t understand is that the constraint is simply a decision matrix of integers (ACR) - a binary decision vector (an) embedded in another matrix along the diagonal* a constant matrix (flows). How is this quadratic, or any different than an x.>=y constraint? Also, is there any way to fix it?

Thanks in advance for any help anyone is able to provide here, I’m completely stumped.

2

I can’t reproduce this.

julia> using JuMP, LinearAlgebra

julia> model = Model();

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

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

julia> flows = [1.0, 2.0];

julia> @constraint(model, ACR .<= Diagonal(1 .- an) * flows)
2×2 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
 an[1] + ACR[1,1] ≤ 1.0    an[1] + ACR[1,2] ≤ 1.0
 2 an[2] + ACR[2,1] ≤ 2.0  2 an[2] + ACR[2,2] ≤ 2.0

Can you provide a reproducible example? Note that your formulation is a bit weird. On the left-hand side of the constraint you have a matrix, and on the right-hand side you have a vector?

3

Flows should be a 2x2 matrix. The result should be one matrix .<= another matrix, since Diagonal embeds the an vector along the diagonal of a matrix of dimension length(an),length(an).

4

And yeah, I’ll post a more reproducible example shortly, I’m also stuck with the issue that sometimes when I try to run the model, it just freezes with no output, and I’m not sure what’s going on there either.

5

Please see my two other replies, also, may I ask what julia version and jump version you’re using, just to rule that out?

6

Ah. That explains it:

julia> using JuMP, LinearAlgebra

julia> model = Model();

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

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

julia> flows = ones(2, 2);

julia> @constraint(model, ACR .<= Diagonal(1 .- an) * flows)
ERROR: MethodError: Cannot `convert` an object of type QuadExpr to an object of type AffExpr
Closest candidates are:
  convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{T, V}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:521
  convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{S, V}) where {S, T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:528
  convert(::Type{GenericAffExpr{T, V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:517
  ...
Stacktrace:
  [1] setindex!
    @ ./array.jl:845 [inlined]
  [2] setindex!
    @ ./multidimensional.jl:645 [inlined]
  [3] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [4] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [5] copyto!
    @ ./broadcast.jl:983 [inlined]
  [6] copyto!
    @ ./broadcast.jl:936 [inlined]
  [7] materialize!
    @ ./broadcast.jl:894 [inlined]
  [8] materialize!
    @ ./broadcast.jl:891 [inlined]
  [9] lmul!(D::Diagonal{AffExpr, Vector{AffExpr}}, B::Matrix{AffExpr})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:212
 [10] *(D::Diagonal{AffExpr, Vector{AffExpr}}, A::Matrix{Float64})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:201
 [11] macro expansion
    @ ~/.julia/packages/MutableArithmetics/Lnlkl/src/rewrite.jl:294 [inlined]
 [12] macro expansion
    @ ~/.julia/packages/JuMP/Y4piv/src/macros.jl:823 [inlined]
 [13] top-level scope
    @ REPL[8]:1

This is another case of LinearAlgebra matrix types and expressions · Issue #66 · jump-dev/MutableArithmetics.jl · GitHub.

The work-around is to wrap Diagonal in Matrix:

julia> using JuMP, LinearAlgebra

julia> model = Model();

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

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

julia> flows = ones(2, 2);

julia> @constraint(model, ACR .<= Matrix(Diagonal(1 .- an)) * flows)
2×2 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}, ScalarShape}}:
 an[1] + ACR[1,1] ≤ 1.0  an[1] + ACR[1,2] ≤ 1.0
 an[2] + ACR[2,1] ≤ 1.0  an[2] + ACR[2,2] ≤ 1.0
7

Okay, thanks! Hopefully, if you have a minute, you can help me with one more related issue?

I noticed an linear algebra error in the constraint and had to adjust it to the following based on your suggested constraint (it’s supposed to set a row and column to 0s and I didn’t catch it sooner because the guy I was collaborating with suggested this solution and I went with it):

 @constraint(model, ACR .<= Matrix(Diagonal(1 .- an)) * flows*Matrix(Diagonal(1 .- an)))

I see now this ends up being a quadratic expression, since I’m multiplying from both sides by my decision variables. And then, I get the following error:

“Constraints of type MathOptInterface.ScalarQuadraticFunction{Float64}-in-MathOptInterface.LessThan{Float64} are not supported by the solver.”

Is there any hope of salvaging this with a linear solver, or am I going to have to figure out Ipopt next?

Thank you so much for your help so far; I’m very new to this, but I was roped in on a grant and the problem just keeps getting bigger, lol.

8

A couple of things:

  • This is a mixed-integer quadratically constrained program, so you’ll need a MISOCP solver for it. Gurobi is a good choice if you have access.
  • Don’t use Ipopt, because it doesn’t support binary and integer variables

But…

  • So long as an is binary, there is a linear reformulation.
model = Model()
N = 2
@variable(model, x[1:N], Bin)
# x_ij[i, j] = x[i] * x[j]
@variable(model, x_ij[1:N, 1:N], Bin)
@constraints(model, begin
    # x_ij is 0 if x_i is 0
    [i=1:N, j=1:N], x_ij[i, j] <= x[i]
    # x_ij is 0 if x_j is 0
    [i=1:N, j=1:N], x_ij[i, j] <= x[j]
    # x_ij is 1 if x_i and x_j are 1
    [i=1:N, j=1:N], x_ij[i, j] >= x[i] + x[j] - 1
end)

I think that means your model is the following, but I didn’t run so I might have made a mistake. It should point you in the right direction though.

using JuMP, LinearAlgebra
N = 2
model = Model()
@variable(model, an[1:N], Bin)
@variable(model, an_ij[1:N, 1:N], Bin)
@variable(model, ACR[1:N, 1:N], Int)
@constraint(model, [i=1:N, j=1:N], an_ij[i, j] <= an[i])
@constraint(model, [i=1:N, j=1:N], an_ij[i, j] <= an[j])
@constraint(model, [i=1:N, j=1:N], an_ij[i, j] >= an[i] + an[j] - 1)
flows = ones(N, N)
@constraint(
    model, 
    [i=1:N, j=1:N], 
    # ACR[i, j] <= (1 - an[i]) * flows[i, j] * (1 - an[j])
    # ACR[i, j] <= flows[i, j] * (1 - an[i]) * (1 - an[j])
    # ACR[i, j] <= flows[i, j] * (1 - an[i] - an[j] + an[i] * an[j])
    ACR[i, j] <= flows[i, j] * (1 - an[i] - an[j] + an_ij[i, j]),
)
1 Like
9

Thank you so much for the clarification and help with the reformulation! I’ll see if I can get my hands on Gurobi and if not, I’ll try out this reformulation instead, since an is indeed binary. Thanks again!

1 Like
10

Hi odow,

I have one more question related to reformulation but on a different objective. Should I continue this thread, or make a new one?

P.S. Gurobi worked a treat for that other constraint! I just have one more tricky one for my objective that seems like it’s quadratic, but it’s because it’s multiple quadratic expressions being multiplied. If that’s the case, I may need to use your (or a similar) reformulation anyway.

2 Likes
11

Start a new thread with a reproducible example

1 Like