Can we support SCS solver?

[pqrz]

Though, I know Juniper is a MINLP solver mainly - meaning IPOPT it is. I wish we can also support SCS solver for the following reasons:

1. This way we might better justify Juniper's support for MISOCP problems.
2. Testing SCS+Juniper on a sample problem (below) says - it does works for SCS, though not very perfectly (gave dual infeasible for a feasible problem). So we aren't very far from SCS.

using Juniper, Ipopt, JuMP, LinearAlgebra, SCS
u0 = [3,5]
# nl_solver = optimizer_with_attributes(Ipopt.Optimizer)
nl_solver = optimizer_with_attributes(SCS.Optimizer)
m = Model(optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>nl_solver))
@variable(m, x[1:2], Int)
@objective(m, Min, dot(u0,x))
@constraint(m, [10, x[1], x[2]] in SecondOrderCone())          # line works after corrections in models.jl
optimize!(m)
println(JuMP.objective_value(m))

3. Given that our usual MISOCP solver: Pajarito is discontinued, we can now have atleast one alternative in Julia.

[ccoffrin]

I think in theory this should work, the root cause of the issue you are having may be numerical issues in SCS.

Thoughts @blegat and @Wikunia?

[blegat]

Making this work was one of the objective of https://github.com/lanl-ansi/Juniper.jl/pull/207. I haven't tested using conic solvers yet so there might still be some work testing it and handling corner cases but we shouldn't be far.

[Mastomaki]

This is a very relevant question! Here's what I get with SCS:

using JuMP
using SCS
using Juniper

optimizer = Juniper.Optimizer
nl_solver = SCS.Optimizer

model = Model(optimizer_with_attributes(optimizer, "nl_solver"=>nl_solver))

@variable(model, 0 <= y[1:2] <= 25)
@variable(model, 0 <= yaux)
@variable(model, x)
@variable(model, ksi[1:1], Bin)

@constraint(model, q_1, [y[1] , y[2] , x ] in SecondOrderCone())

@constraint(model, q_2, yaux == y[1] + (1 - ksi[1] ))

@objective(model, Max, x )
optimize!(model)

The error message:

ERROR: LoadError: UndefVarError: value not defined Stacktrace: [1] set_subsolver_option! at /home/jube/.julia/packages/Juniper/8wso7/src/util.jl:289 [inlined] [2] set_subsolver_option!(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{SCS.Optimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}, ::String, ::String, ::Pair{Float64,Float64}) at /home/jube/.julia/packages/Juniper/8wso7/src/util.jl:292 [3] process_node!(::Juniper.JuniperProblem, ::Juniper.StepObj, ::Juniper.BnBNode, ::Array{Int64,1}, ::Bool; restarts::Int64) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:67 [4] process_node! at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:56 [inlined] [5] branch!(::Juniper.JuniperProblem, ::Juniper.SolverOptions, ::Juniper.StepObj, ::Int64, ::Array{Int64,1}; temp::Bool) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:161 [6] branch_strong_on!(::Juniper.JuniperProblem, ::Juniper.SolverOptions, ::Juniper.StepObj, ::Array{Int64,1}, ::Array{Int64,1}, ::Bool, ::Int64) at /home/jube/.julia/packages/Juniper/8wso7/src/bb_strategies.jl:164 [7] branch_strong! at /home/jube/.julia/packages/Juniper/8wso7/src/bb_strategies.jl:275 [inlined] [8] get_branching_disc_idx!(::Juniper.JuniperProblem, ::Juniper.StepObj, ::Juniper.SolverOptions, ::Array{Int64,1}, ::Juniper.GainObj, ::Int64) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:29 [9] one_branch_step!(::Juniper.JuniperProblem, ::Nothing, ::Juniper.SolverOptions, ::Juniper.StepObj, ::Array{Int64,1}, ::Juniper.GainObj, ::Int64) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:311 [10] solve_sequential(::Juniper.BnBTreeObj, ::Array{Any,1}, ::Float64, ::Array{Symbol,1}, ::Array{Int64,1}, ::Juniper.TimeObj) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:408 [11] solvemip(::Juniper.BnBTreeObj) at /home/jube/.julia/packages/Juniper/8wso7/src/BnBTree.jl:608 [12] optimize!(::Juniper.Optimizer) at /home/jube/.julia/packages/Juniper/8wso7/src/MOI_wrapper/MOI_wrapper.jl:297 [13] optimize!(::MathOptInterface.Bridges.LazyBridgeOptimizer{Juniper.Optimizer}) at /home/jube/.julia/packages/MathOptInterface/5WwpK/src/Bridges/bridge_optimizer.jl:293 [14] optimize!(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at /home/jube/.julia/packages/MathOptInterface/5WwpK/src/Utilities/cachingoptimizer.jl:237 [15] optimize!(::Model, ::Nothing; bridge_constraints::Bool, ignore_optimize_hook::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/jube/.julia/packages/JuMP/y5vgk/src/optimizer_interface.jl:139 [16] optimize! at /home/jube/.julia/packages/JuMP/y5vgk/src/optimizer_interface.jl:115 [inlined] (repeats 2 times)

[blegat]

This looks like an easy fix. The value at https://github.com/lanl-ansi/Juniper.jl/blob/f88cfbeb016947fd8da1f71671363691f14b8e95/src/util.jl#L289 should be removed

[Mastomaki]

Yes but the return value is used in process_node! (BnBTree.jl#L54) to store the old value of the parameter "mu_init". The function further assumes that subsolver is Ipopt. Probably there is no such parameter for SCS.

[Wikunia]

That isn't a problem though as nothing will be returned then and it will not be set we don't use Ipopt again when resetting. There currently is another problem though which I try to debug right now. Will hopefully make a PR soon for this.

[blegat]

When the return value is used, a Pair is given so it calls this method instead: https://github.com/lanl-ansi/Juniper.jl/blob/f88cfbeb016947fd8da1f71671363691f14b8e95/src/util.jl#L272

[Wikunia]

The current problem is that the default tolerance has some problems. We use a tolerance of 1e-6 by default but the tolerance of SCS might be different such that Juniper fails due to the fact that the lower and upper bound of a variable are the same but the solution value is still not discrete. When you (@Mastomaki ) run your problem as:

using JuMP
using SCS
using Juniper

optimizer = Juniper.Optimizer
nl_solver = SCS.Optimizer

model = Model(optimizer_with_attributes(optimizer, "nl_solver"=>optimizer_with_attributes(nl_solver, "verbose"=>0), "atol"=>1e-4))

@variable(model, 0 <= y[1:2] <= 25)
@variable(model, 0 <= yaux)
@variable(model, x)
@variable(model, ksi[1:1], Bin)

@constraint(model, q_1, [y[1] , y[2] , x ] in SecondOrderCone())

@constraint(model, q_2, yaux == y[1] + (1 - ksi[1] ))

@objective(model, Max, x )
optimize!(model)

it finishes with LOCALLY_SOLVED but the value for y[1] is a bit higher than 25.

[Mastomaki]

Cool! Thanks for the tip to adjust the tolerance. I suppose one could also tune SCS tolerance like below:

nl_solver = optimizer_with_attributes(SCS.Optimizer, "verbose"=>0, "eps"=>1e-6)

[odow]

Can this be closed? It seems like an SCS tolerance issue, not a problem in Juniper.

[ccoffrin]

If this is working, probably should add a test to that effect before closing this issue. Maybe the example shown in this issue is sufficient?