SOS1 Constraints in Bonmin Solver - JuMP / AmplNLWriter

I would like to use Bonmin to perform MINLP. The model I have in mind has some SOS constraints of type 1. According to the bonmin documentation page 2, the B-BB algorithm accepts SOS1 constraints.

When I try to solve the model in JuMP I get the following error:

julia> solve(m)
ERROR: Not implemented for SOS constraints
Stacktrace:
 [1] constraintbounds(::JuMP.Model, ::JuMP.ProblemTraits) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:1033
 [2] _buildInternalModel_nlp(::JuMP.Model, ::JuMP.ProblemTraits) at /home/michael/.julia/v0.6/JuMP/src/nlp.jl:1243
 [3] #build#119(::Bool, ::Bool, ::JuMP.ProblemTraits, ::Function, ::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:303
 [4] (::JuMP.#kw##build)(::Array{Any,1}, ::JuMP.#build, ::JuMP.Model) at ./<missing>:0
 [5] #solve#116(::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::Function, ::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:168
 [6] solve(::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:150
 [7] macro expansion at ./REPL.jl:97 [inlined]
 [8] (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:73

I have created a minimal reproducible example that solves the best subset selection problem. I solve this first using Gurobi, which works, and then with Bonmin, which fails.

using JuMP
using Distributions
using Gurobi
n=200
p=20
B=zeros(p)
B[1]=1
B[2]=2
B[3]=3
X=reshape(rand(Normal(0,1),n*p),n,p)
Y=X*B+rand(Normal(0,1),n)


#GUROBI

m = Model(solver=GurobiSolver())
@variable(m, B[1:p])
@variable(m, G[1:p], Bin)
@variable(m, Z[1:p], Bin)
for i=1:p
    addSOS1(m, [B[i],Z[i]])
end
@constraint(m, Z .== 1 .-G )
@constraint(m, sum(G[j] for j=1:p) <= 3)
@objective(m, Min, sum((Y[i]-sum( X[i,j]*B[j] for j=1:p))^2 for i=1:n))
solve(m)

#BONMIN

using AmplNLWriter
m = Model(solver=AmplNLSolver("bonmin", ["bonmin.algorithm=B-BB"]))
@variable(m, B[1:p])
@variable(m, G[1:p], Bin)
@variable(m, Z[1:p], Bin)
for i=1:p
    addSOS1(m, [B[i],Z[i]])
end
@constraint(m, Z .== 1 .-G )
@constraint(m, sum(G[j] for j=1:p) <= 3)
@NLobjective(m, Min, sum((Y[i]-sum( X[i,j]*B[j] for j=1:p))^2 for i=1:n))
solve(m)

I’m just guessing, but maybe it works for Gurobi because it uses the LinearQuadratic interface, while AmplNLWriter uses the Nonlinear interface, where the latter does not know about SOS1 constraints.

Are there any work-arounds or alternatives I can use to solve minlp with sos constraints in julia?

I have tried a different route, using CoinOptServices.jl

using CoinOptservices
m = Model(solver=OsilBonminSolver())
@variable(m, B[1:p])
@variable(m, G[1:p], Bin)
@variable(m, Z[1:p], Bin)
for i=1:p
    addSOS1(m, [B[i],Z[i]])
end
@constraint(m, Z .== 1 .-G )
@constraint(m, sum(G[j] for j=1:p) <= 3)
@NLobjective(m, Min, sum((Y[i]-sum( X[i,j]*B[j] for j=1:p))^2 for i=1:n))

solve(m)

and get the same error



julia> solve(m)
ERROR: Not implemented for SOS constraints
Stacktrace:
 [1] constraintbounds(::JuMP.Model, ::JuMP.ProblemTraits) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:1033
 [2] _buildInternalModel_nlp(::JuMP.Model, ::JuMP.ProblemTraits) at /home/michael/.julia/v0.6/JuMP/src/nlp.jl:1243
 [3] #build#119(::Bool, ::Bool, ::JuMP.ProblemTraits, ::Function, ::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:303
 [4] (::JuMP.#kw##build)(::Array{Any,1}, ::JuMP.#build, ::JuMP.Model) at ./<missing>:0
 [5] #solve#116(::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::Function, ::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:168
 [6] solve(::JuMP.Model) at /home/michael/.julia/v0.6/JuMP/src/solvers.jl:150
 [7] macro expansion at ./REPL.jl:97 [inlined]
 [8] (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:73

I feel like you will always run into the same problem: only the LinearQuadratic interface supports SOS-type constraints.

In case you can derive some finite bounds for your B variables, you could just reformulate it as big-M type linear inequalities, right?

I do have a finite bound, however, the Big-M constraint isn’t as powerful I believe as the SOS1 constraint. Its my understanding that specifying an SOS1 constraint affects the branching procedure and results in a much faster algorithm than the Big-M.

You’ll have to dig deeper and see precisely through what API interface Bonmin accepts SOS1 constraints. This doesn’t seem clear from the manual. After that, we can see how to pipe the constraints through from JuMP.

I believe it is with their B-BB branch and bound algorithm. Let me try and get an example using SOS constraints to run via AMPL using Bonmin as its solver and I will report back.