# Complementarity.jl with Path
**URL:** https://discourse.julialang.org/t/complementarity-jl-with-path/103081
**Category:** Optimization (Mathematical)
**Tags:** jump, error, complementarity
**Created:** [August 22, 2023, 4:22pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081 "2023-08-22T16:22:53Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@sjhon](https://discourse.julialang.org/u/sjhon)
#### Post date: [August 22, 2023, 4:22pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/1 "2023-08-22T16:22:53Z")
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](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
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 22, 2023, 8:17pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/2 "2023-08-22T20:17:56Z")
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`.
---
### Author:  [@sjhon](https://discourse.julialang.org/u/sjhon)
#### Post date: [August 22, 2023, 8:44pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/3 "2023-08-22T20:44:51Z")
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
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 22, 2023, 8:57pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/4 "2023-08-22T20:57:55Z")
> [@sjhon](#):
>
> @complementarity(mcp, KKT1[i], x)
This is still not right. You are adding five complementarity constraints to the variable `x`. It must have exactly one.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [August 23, 2023, 2:56am UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/5 "2023-08-23T02:56:12Z")
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](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):
```julia
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
```
---
### Author:  [@sjhon](https://discourse.julialang.org/u/sjhon)
#### Post date: [September 28, 2023, 7:23pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/6 "2023-09-28T19:23:21Z")
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [September 28, 2023, 9:04pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/7 "2023-09-28T21:04:49Z")
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:
- [ANN: JuMP v1.15 is released](https://discourse.julialang.org/t/ann-jump-v1-15-is-released/103875)
- [Mixed complementarity problems · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/complementarity/)
Here’s what the translation would look like:
```plaintext
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
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:
```julia
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.
---
### Author:  [@sebacea](https://discourse.julialang.org/u/sebacea)
#### Post date: [October 1, 2023, 11:12pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/8 "2023-10-01T23:12:41Z")
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](https://github.com/jump-dev/JuMP.jl/commit/a573de22a3c52160dcfc2e3bf06fe58b5243a110) revised this issue and proposes what I think is [a correct computation of the example](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/complementarity/#Electricity-consumption). Thanks you for the excellent material.
I just made the point for future readers.
---
### Author:  [@vicentemunozlh](https://discourse.julialang.org/u/vicentemunozlh)
#### Post date: [December 13, 2023, 6:20pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/9 "2023-12-13T18:20:32Z")
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?
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [December 13, 2023, 7:48pm UTC](https://discourse.julialang.org/t/complementarity-jl-with-path/103081/10 "2023-12-13T19:48:34Z")
The constraint was:
```julia
@constraint(m, x[1] == x[i])
```
This is equivalent to
```julia
@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.