Complementarity.jl with Path

Specific Domains Optimization (Mathematical)
, ,
1

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 https://www.repository.cam.ac.uk/items/6b997afc-3412-452e-af38-75ef1c00e29a
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

2

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.

3

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

4

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

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
6

I am giving this problem another try by converting it to an MCP. I think I have done it correctly, but I am encountering a problem when passing it to Julia with complementarity and Path:

julia

Problem Data

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;
technologies = [1, 2, 3, 4, 5]

Complementarity Model

using Complementarity, JuMP
using PATHSolver
PATHSolver.c_api_License_SetString(“License”)

m = MCPModel()

Variables

@variable(m, x[i in technologies] >= 0) # Capacity
@variable(m, Q[i in technologies] >= 0) # Supply and Demand
@variable(m, alpha[i in technologies] >= 0) # Dual associated with Q = yesto e

Ensure that all instances of x[i] are equal

for i in 2:length(technologies)
@constraint(m, x[1] == x[i])
end

Complementarity 1

@mapping(m, max_capacity[i in technologies], x[i] - Q[i])
@complementarity(m, max_capacity, alpha)

Complementarity 2

@mapping(m, kkt_Q[i in technologies], Tau * sum((theta[j] * C - theta[j] * A[j] + theta[j] * B * Q[j]) for j in technologies) + alpha[i])
@complementarity(m, kkt_Q, Q)

Complementarity 3

@mapping(m, kkt_x[i in technologies], I*1000 - alpha[i])
@complementarity(m, kkt_x, x)

Solving the model

status = solveMCP(m)

println("x = ", result_value.(x))
println("Q = ", result_value.(Q))

And I get these results:
x = 1-dimensional DenseAxisArray{Float64,1,…} with index sets:
Dimension 1, [1, 2, 3, 4, 5]
And data, a 5-element Vector{Float64}:
329.72602740364846
329.72602738700664
329.72602742164287
329.7260273870066
329.72602738700647
Q = 1-dimensional DenseAxisArray{Float64,1,…} with index sets:
Dimension 1, [1, 2, 3, 4, 5]
And data, a 5-element Vector{Float64}:
329.72602740364886
329.726027387007
329.72602742164327
329.726027387007
329.72602738700687

However, these are not the same as the real problem. Do you see any issues with my code? As you can see, the variable x has to be the same for all technologies.

7

Some good news is that you can stop using Complementarity.jl. We’ve recently added first-class support for mixed-complementarity problems to JuMP. See:

Here’s what the translation would look like:

using JuMP
import PATHSolver
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
technologies = [1, 2, 3, 4, 5]
m = Model(PATHSolver.Optimizer)
@variables(m, begin
    x[technologies] >= 0
    Q[technologies] >= 0
    alpha[technologies] >= 0
end)
for i in 2:length(technologies)
    @constraint(m, x[1] == x[i])
end
@constraints(m, begin
    max_capacity[i in technologies], x[i] - Q[i] ⟂ alpha[i]
    kkt_Q[i in technologies], Tau * sum(theta[j] * (C - A[j] + B * Q[j]) for j in technologies) + alpha[i] ⟂ Q[i]
    kkt_x[i in technologies], 1000 * I - alpha[i] ⟂ x[i]
end)
optimize!(m)
println("x = ", value.(x))
println("Q = ", value.(Q))

But if you run the code, you’ll get the error:

julia> for i in 2:length(technologies)
           @constraint(m, x[1] == x[i])
       end
ERROR: Constraints of type MathOptInterface.ScalarAffineFunction{Float64}-in-MathOptInterface.EqualTo{Float64} are not supported by the solver.

If you expected the solver to support your problem, you may have an error in your formulation. Otherwise, consider using a different solver.

The list of available solvers, along with the problem types they support, is available at https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] _moi_add_constraint(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.GenericOptimizer{Float64, MathOptInterface.Utilities.ObjectiveContainer{Float64}, MathOptInterface.Utilities.VariablesContainer{Float64}, PATHSolver.OptimizerFunctionConstraints{Float64}}}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, f::MathOptInterface.ScalarAffineFunction{Float64}, s::MathOptInterface.EqualTo{Float64})
   @ JuMP ~/.julia/packages/JuMP/ToPd2/src/constraints.jl:679
 [3] add_constraint(model::Model, con::ScalarConstraint{AffExpr, MathOptInterface.EqualTo{Float64}}, name::String)
   @ JuMP ~/.julia/packages/JuMP/ToPd2/src/constraints.jl:706
 [4] macro expansion
   @ ~/.julia/packages/JuMP/ToPd2/src/macros.jl:1345 [inlined]
 [5] top-level scope
   @ ./REPL[14]:2

That’s because you cannot add equality constraints to a mixed-complementarity problem.

You need to do instead:

using JuMP
import PATHSolver
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
technologies = [1, 2, 3, 4, 5]
m = Model(PATHSolver.Optimizer)
@variables(m, begin
    x[technologies] >= 0
    Q[technologies] >= 0
    alpha[technologies] >= 0
    μ[technologies[2:end]]   # !!! NEW !!!
end)
@constraints(m, begin
    [i=technologies[2:end]], x[1] - x[i] ⟂ μ[i]   # !!! NEW !!!
    max_capacity[i in technologies], x[i] - Q[i] ⟂ alpha[i]
    kkt_Q[i in technologies], Tau * sum(theta[j] * (C - A[j] + B * Q[j]) for j in technologies) + alpha[i] ⟂ Q[i]
    kkt_x[i in technologies], 1000 * I - alpha[i] ⟂ x[i]
end)
optimize!(m)
println("x = ", value.(x))
println("Q = ", value.(Q))

and now I get the same solution as Complementarity.jl.

If this solution isn’t what you expect, then you have a mistake somewhere in your equations.

1 Like
8

Thanks @odow, I am following this discussion with @sjhon and the translation you propose understands the example with the possibility of multiple technologies. But Example 3.3.1 has only one.

Nevertheles, I see that the tutorial you prepared revised this issue and proposes what I think is a correct computation of the example. Thanks you for the excellent material.

I just made the point for future readers.

1 Like
9

Hi @odow ! I am finding JuMP super good for MCP programming.
This is not really about the package itself but about how you got to the MCP formulation on this problem. I am finding it really hard to do so in a problem that initially appears to be super easy.
Do you have somewhere stored the iterations to get to it?

10

The constraint was:

@constraint(m, x[1] == x[i])

This is equivalent to

@constraint(m, x[1] - x[i] == 0)

By the definition of a complementarity constraint, f(x) = 0 if the complementing variable has free bounds. So we add a new variable μ[i] and the complementarity condition x[1] - x[i] ⟂ μ[i].

You can also think of μ[i] as being the dual variable associated with the equality constraint. In this model, μ[i] shows up only in the one complementarity constraint, but in others, it might be needed elsewhere.