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.