Complementarity.jl with Path

Specific Domains Optimization (Mathematical)
, ,
5

I took another look at this problem, since I thought it’d be a useful test case for https://github.com/jump-dev/JuMP.jl/pull/3106 and some of the complementarity tutorials I’m writing.

The code below won’t run by default, you’ll need to convert it to the Complementarity.jl syntax, or live dangerously and install the development branch via ] add JuMP#od/nlp-expr.

The example in the textbook is a little tricky to formulate explicitly as an MCP, so I see how you got stuck. This is what I came up with. It finds the same solution (rounded to the nearest MW):

A = [300, 350, 400, 450, 500]
B = 1                        
I = 90_000  # Given in the text in terms of $/kWh!
C = 60                       
τ = 8_760                    
θ = [0.2, 0.2, 0.2, 0.2, 0.2]
model = Model(PATHSolver.Optimizer)
@variables(model, begin
    x >= 0, (start = 1)
    Q[1:5] >= 0, (start = 1)
    Y[1:5] >= 0, (start = 1)
    P[1:5], (start = 1)
    μ[1:5] >= 0, (start = 1)
end)
@constraints(model, begin
    I - τ * θ' * μ ⟂ x
    [i=1:5],  C - P[i] + μ[i] ⟂ Y[i]
    [i=1:5], A[i] - B * Q[i] - P[i] ⟂ Q[i]
    [i=1:5], Y[i] - Q[i] ⟂ P[i]
    [i=1:5], x - Y[i] ⟂ μ[i]
end)
optimize!(model)

julia> value(x)
389.31506849315065

julia> hcat(
           value.(Q),
           value.(P),
           (τ * value.(μ) .- I) / 1_000,
       )'
3×5 adjoint(::Matrix{Float64}) with eltype Float64:
 240.0  290.0  340.0  389.315   389.315
  60.0   60.0   60.0   60.6849  110.685
 -90.0  -90.0  -90.0  -84.0     354.0
2 Likes