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)
```