Calculating duals with SOS2 constraints

Specific Domains Optimization (Mathematical)
,
1

Hi,

I am working with a model in JuMP and Gurobi which is usually a LP or MILP. I have added a set of SOS2 constraints to create a piecewise linear function. Gurobi successfully reformulates the constraints to avoid needing integer or binary variables (if the original problem is an LP) and the model solves correctly. However, the model no longer has duals.

When faced with this problem with MILPs, we solved the MILP and then used the function below to convert it to an LP and call optimize!() again. This is not a perfect approach but worked fine for us.

Does anyone have suggestions for how to achieve the same results with the SOS2 constraints? I also tried to fix my lambda variables but the SOS2 constraints are still there and still no duals.

Thanks

function fix_model(jump_model)
        values = Dict(v => value(v) for v in all_variables(jump_model))
	for v in all_variables(jump_model)
		if is_integer(v)
            fix(v,values[v],force=true)
			unset_integer(v)
        elseif is_binary(v)
            fix(v,values[v],force=true)
			unset_binary(v)
        end
	end
2 Likes
2

The answer to this depends on “what do you expect the dual to return and why do you need it?”

Note that an SOS2 constraint is implemented either by the solver branching, or by adding new binary variables. In both cases, it does not have a “natural” dual solution.

In the special case that you are using SOS2 to form a piecewise linear function, then it might make sense to talk about a dual.

Since the SOS2 is “at most two elements are non-zero, and if there are two non-zero elements, they must be adjacent in the ordering,” then one possible relaxation is to fix all the lambda variables to their optimal solution. The other is to fix the zero lambda variables to 0, leaving the two non-zero variables free. In both cases you can then delete the SOS2 constraint.

As one example:

julia> using JuMP, Gurobi

julia> begin
           x̂ = -0.5:0.6:2.5
           ŷ = x̂ .^ 2
           N = length(x̂)
           model = Model(Gurobi.Optimizer)
           set_silent(model)
           @variable(model, -1 <= x <= 2)
           @variable(model, y)
           @variable(model, 0 <= λ[1:N] <= 1)
           @objective(model, Max, y)
           @constraints(model, begin
               c_dual, x == sum(x̂[i] * λ[i] for i in 1:N)
               y == sum(ŷ[i] * λ[i] for i in 1:N)
               sum(λ) == 1
               c_sos, λ in SOS2()
           end)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           dual_status(model)
       end
NO_SOLUTION::ResultStatusCode = 0

julia> λ_value = value.(λ)
6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.8333333333333334
 0.16666666666666666

julia> fix.(λ[λ_value .<= 1e-6], 0.0; force = true);

julia> delete(model, c_sos)

julia> optimize!(model)

julia> @assert is_solved_and_feasible(model)

julia> dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

julia> dual(c_dual)
4.3999999999999995

The other option is:

julia> using JuMP, Gurobi

julia> begin
           x̂ = -0.5:0.6:2.5
           ŷ = x̂ .^ 2
           N = length(x̂)
           model = Model(Gurobi.Optimizer)
           set_silent(model)
           @variable(model, -1 <= x <= 2)
           @variable(model, y)
           @variable(model, 0 <= λ[1:N] <= 1)
           @objective(model, Max, y)
           @constraints(model, begin
               c_dual, x == sum(x̂[i] * λ[i] for i in 1:N)
               y == sum(ŷ[i] * λ[i] for i in 1:N)
               sum(λ) == 1
               c_sos, λ in SOS2()
           end)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           dual_status(model)
       end
NO_SOLUTION::ResultStatusCode = 0

julia> λ_value = value.(λ)
6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.8333333333333334
 0.16666666666666666

julia> fix.(λ, λ_value; force = true);

julia> delete(model, c_sos)

julia> optimize!(model)

julia> @assert is_solved_and_feasible(model)

julia> dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1

julia> dual(c_dual)
0.0

Note that the resulting “dual” depends on how you choose to relax the problem.

You’d need to integrate this into your own fix_model function. There is nothing built-in to JuMP.

3 Likes
3

Thanks very much @odow! I was almost there but for some reason I thought setting force=true when fixing the lambdas would automatically delete the SOS2 constraints. Clearly not.

Both approaches worked for me and I’m going to do some testing to check which outcome r.e. the duals makes most sense. I think it will be the former from some initial runs.

For anyone looking in the future, the resolved code is in the Dolphyn.jl repository. It will be saved in src/core/piecewise_linear and src/solve_model . I hope to push it this week.

I’m using anonymous variables and constraints as I don’t know the number of piecewise variables ahead of time. As a quick solution, I’m storing them in the Model.obj_dict. I’m happy to take suggestions on better approaches.

2 Likes