Complementarity.jl whit Path

Im trying to solve a problema of d’aertrycke, ehrenmann, ralph & smeers (2017), the problem is the 3.3.1 that you can found here Risk trading in capacity equilibrium models
Im trying to do this problem as a Mixed Complementary problem (MCP), whit complementarity and Path.
Here is my code. I put the @ before variable, complementarity and mapping and show result_value but i had to remove it
using PATHSolver, JuMP, Complementarity,

PATHSolver.c_api_License_SetString(“Licence”)

Datos del problema

A = [300, 350, 400, 450, 500]
theta = [0.2, 0.2, 0.2, 0.2, 0.2]
B = 1
I = 90
C = 60
tau = 8760

Modelo de complementariedad

mcp = Model(PATHSolver.Optimizer)

Variables

variable(mcp, x >= 0) # Capacidad
variable(mcp, Q[i=1:5] >= 0)
variable(mcp, alpha[i=1:5] >= 0)

Mapeos de complementariedad

mapping(mcp, KKT1[i=1:5], I - alpha[i])
mapping(mcp, KKT2[i=1:5], (tau * sum((theta[j] * (C - A[j] + B * Q[j])) for j in 1:5)) + alpha[i])

Complementariedad

for i in 1:5
complementarity(mcp, KKT1[i], Q[i])
complementarity(mcp, KKT2[i], x)
end

RESOLUCION

status = solveMCP(mcp, solver=:PATH, minor_iteration_limit=100000)

show result_value.(x)
show result_value.(Q)

and this is the error:
The ‘get_MCP_data’ function is only for MCP models as in ComplementarityType.jl

Stacktrace:
[1] error(s::String)
Base .\error.jl:35
[2] get_MCP_data(m::Model)
Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:38
[3] add_complementarity(m::Model, var::VariableRef, F::NonlinearExpression, F_name::String; fixed_variable::Bool)
Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:396
[4] add_complementarity(m::Model, var::VariableRef, F::NonlinearExpression, F_name::String)
Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:390
[5] macro expansion
C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:430 [inlined]
[6] top-level scope
In[39]:27

Your formulation isn’t quite right. A mixed complementarity problem has exactly one complementarity constraint for each variable.

In your current formulation, you have multiple @complementarity constraints for x, and you are missing complementarity constraints for alpha.

I see what you mean, and I think now it’s supposed to be the correct code.
using PATHSolver, JuMP, Complementarity

PATHSolver.c_api_License_SetString(“Licence”)

Datos del problema

A = [300, 350, 400, 450, 500]
theta = [0.2, 0.2, 0.2, 0.2, 0.2]
B = 1
I = 90
C = 60
tau = 8760

Modelo de complementariedad

mcp = Model(PATHSolver.Optimizer)

Variables

@variable(mcp, x >= 0)
@variable(mcp, Q[i=1:5] >= 0)
@variable(mcp, alpha[i=1:5] >= 0)

Mapeos de complementariedad

@mapping(mcp, KKT1[i=1:5], I - alpha[i])
@mapping(mcp, KKT2[i=1:5], (tau * sum((theta[j] * (C - A[j] + B * Q[j])) for j in 1:5)) + alpha[i])
@mapping(mcp, KKT3[i=1:5], (x-Q[i]))

Complementariedad

for i in 1:5
@complementarity(mcp, KKT1[i], x)
@complementarity(mcp, KKT2[i], Q[i])
@complementarity(mcp, KKT3[i], alpha[i])
end

RESOLUCION

status = solveMCP(mcp, solver=:PATH, minor_iteration_limit=100000)

@show result_value.(x)
@show result_value.(Q)

But now I have this error, and I don’t know how to make it right:
The ‘get_MCP_data’ function is only for MCP models as in ComplementarityType.jl

Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] get_MCP_data(m::Model)
@ Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:38
[3] add_complementarity(m::Model, var::VariableRef, F::NonlinearExpression, F_name::String; fixed_variable::Bool)
@ Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:396
[4] add_complementarity(m::Model, var::VariableRef, F::NonlinearExpression, F_name::String)
@ Complementarity C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:390
[5] macro expansion
@ C:\Users\sjhon.julia\packages\Complementarity\NbsQM\src\mcp.jl:421 [inlined]
[6] top-level scope
@ In[1]:29

This is still not right. You are adding five complementarity constraints to the variable x. It must have exactly one.

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
1 Like