Example: mixed complementarity problems
This tutorial was generated using Literate.jl. Download the source as a .jl
file.
The purpose of this tutorial is to provide a collection of mixed complementarity programs.
Required packages
This tutorial uses the following packages:
using JuMP
import PATHSolver
Background
A mixed complementarity problem has the form:
\[\begin{align} F_i(x) \perp x_i & i = 1 \ldots n \\ l_i \le x_i \le u_i & i = 1 \ldots n. \end{align}\]
where the $\perp$ constraint enforces the following relations:
- If $l_i < x_i < u_i$, then $F_i(x) = 0$
- If $l_i = x_i$, then $F_i(x) \ge 0$
- If $x_i = u_i$, then $F_i(x) \le 0$
You may have seen a complementarity problem written as $0 \le F(x) \perp x \ge 0$. This is a special case of a mixed complementarity problem in which $l_i = 0$ and $u_i = \infty$.
Importantly, a mixed complementarity problem does not have an objective, and no other constraint types are present.
Linear complementarity
Form a mixed complementarity problem using the perp symbol ⟂
(type \perp<tab>
in the REPL).
M = [0 0 -1 -1; 0 0 1 -2; 1 -1 2 -2; 1 2 -2 4]
q = [2, 2, -2, -6]
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, 0 <= x[1:4] <= 10, start = 0)
@constraint(model, M * x + q ⟂ x)
optimize!(model)
@assert is_solved_and_feasible(model)
value.(x)
4-element Vector{Float64}:
2.8
0.0
0.8
1.2
Other ways of writing linear complementarity problems
You do not need to use a single vector of variables, and the complementarity constraints can be given in any order. In addition, you can use the perp symbol, the complements(F, x)
syntax, or the MOI.Complements
set.
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, 0 <= w <= 10, start = 0)
@variable(model, 0 <= x <= 10, start = 0)
@variable(model, 0 <= y <= 10, start = 0)
@variable(model, 0 <= z <= 10, start = 0)
@constraint(model, complements(y - 2z + 2, x))
@constraint(model, [-y - z + 2, w] in MOI.Complements(2))
@constraint(model, w + 2x - 2y + 4z - 6 ⟂ z)
@constraint(model, w - x + 2y - 2z - 2 ⟂ y)
optimize!(model)
@assert is_solved_and_feasible(model)
value.([w, x, y, z])
4-element Vector{Float64}:
2.8
0.0
0.8
1.2
Transportation
This example is a reformulation of the transportation problem from Chapter 3.3 of Dantzig, G.B. (1963). Linear Programming and Extensions. Princeton University Press, Princeton, New Jersey. It is based on the GAMS model gamslib_transmcp
.
capacity = Dict("seattle" => 350, "san-diego" => 600)
demand = Dict("new-york" => 325, "chicago" => 300, "topeka" => 275)
cost = Dict(
("seattle" => "new-york") => 90 * 2.5 / 1_000,
("seattle" => "chicago") => 90 * 1.7 / 1_000,
("seattle" => "topeka") => 90 * 1.8 / 1_000,
("san-diego" => "new-york") => 90 * 2.5 / 1_000,
("san-diego" => "chicago") => 90 * 1.8 / 1_000,
("san-diego" => "topeka") => 90 * 1.4 / 1_000,
)
plants, markets = keys(capacity), keys(demand)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, w[i in plants] >= 0)
@variable(model, p[j in markets] >= 0)
@variable(model, x[i in plants, j in markets] >= 0)
@constraints(
model,
begin
[i in plants, j in markets], w[i] + cost[i=>j] - p[j] ⟂ x[i, j]
[i in plants], capacity[i] - sum(x[i, :]) ⟂ w[i]
[j in markets], sum(x[:, j]) - demand[j] ⟂ p[j]
end
)
optimize!(model)
@assert is_solved_and_feasible(model)
value.(p)
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
Dimension 1, ["new-york", "chicago", "topeka"]
And data, a 3-element Vector{Float64}:
0.22500000000033224
0.15299999994933886
0.126
Expected utility of insurance
This example is taken from a lecture of the course AAE706, given by Thomas F. Rutherford at the University of Wisconsin, Madison. It models the expected coverage of insurance K
that a rational actor would obtain to insure a risk that occurs with probability pi
and results in a loss of L
.
pi = 0.01 # Probability of a bad outcome
L = 0.5 # Loss with a bad outcome
γ = 0.02 # Premium for coverage
σ = 0.5 # Elasticity
ρ = -1 # Risk exponent
U(C) = C^ρ / ρ
MU(C) = C^(ρ - 1)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, EU, start = 1) # Expected utility
@variable(model, EV, start = 1) # Equivalent variation in income
@variable(model, C_G, start = 1) # Consumption on a good day
@variable(model, C_B, start = 1) # Consumption on a bad day
@variable(model, K, start = 1) # Coverage
@constraints(
model,
begin
(1 - pi) * U(C_G) + pi * U(C_B) - EU ⟂ EU
100 * (((1 - pi) * C_G^ρ + pi * C_B^ρ)^(1 / ρ) - 1) - EV ⟂ EV
1 - γ * K - C_G ⟂ C_G
1 - L + (1 - γ) * K - C_B ⟂ C_B
γ * ((1 - pi) * MU(C_G) + pi * MU(C_B)) - pi * MU(C_B) ⟂ K
end
)
optimize!(model)
@assert is_solved_and_feasible(model)
value(K)
0.20474003534537774
Electricity consumption
This example is a mixed complementarity formulation of Example 3.3.1 from (D’Aertrycke et al., 2017).
This example models a risk neutral competitive equilibrium between a producer and a consumer of electricity.
In our example, we assume a producer is looking to invest in a new power plant with capacity $x$ [MW]. This plant has an annualized capital cost of $I$ [€/MW] and an operating cost of $C$ [€/MWh]. There are 8760 hours in a year.
After making the capital investment, there are five possible consumption scenarios, $\omega$, which occur with probability $\theta_\omega$. In each scenario , the producer makes $Y_ω$ MW of electricity.
There is one consumer in the model, who has a quadratic utility function, $U(Q_ω) = A_ω Q_ω + \frac{B_ω Q_ω^2}{2}$.
We now build and solve the mixed complementarity problem with a few brief comments. The economic justification for the model would require a larger tutorial than the space available here. Consult (D’Aertrycke et al., 2017) for details.
I = 90_000 # Annualized capital cost
C = 60 # Operation cost per MWh
τ = 8_760 # Hours per year
θ = [0.2, 0.2, 0.2, 0.2, 0.2] # Scenario probabilities
A = [300, 350, 400, 450, 500] # Utility function coefficients
B = 1 # Utility function coefficients
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, x >= 0, start = 1) # Installed capacity
@variable(model, Q[ω = 1:5] >= 0, start = 1) # Consumption
@variable(model, Y[ω = 1:5] >= 0, start = 1) # Production
@variable(model, P[ω = 1:5], start = 1) # Electricity price
@variable(model, μ[ω = 1:5] >= 0, start = 1) # Capital scarcity margin
# Unit investment cost equals annualized scarcity profit or investment is 0
@constraint(model, I - τ * θ' * μ ⟂ x)
# Difference between price and scarcity margin is equal to operation cost
@constraint(model, [ω = 1:5], C - (P[ω] - μ[ω]) ⟂ Y[ω])
# Price is equal to consumer's marginal utility
@constraint(model, [ω = 1:5], P[ω] - (A[ω] - B * Q[ω]) ⟂ Q[ω])
# Production is equal to consumption
@constraint(model, [ω = 1:5], Y[ω] - Q[ω] ⟂ P[ω])
# Production does not exceed capacity
@constraint(model, [ω = 1:5], x - Y[ω] ⟂ μ[ω])
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)
* Solver : Path 5.0.03
* Status
Result count : 1
Termination status : LOCALLY_SOLVED
Message from the solver:
"The problem was solved"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Objective value : 0.00000e+00
* Work counters
Solve time (sec) : 1.45000e-04
An equilibrium solution is to build 389 MW:
value(x)
389.31506849315065
The production in each scenario is:
value.(Q)
5-element Vector{Float64}:
240.0000000000001
289.9999999999999
340.0
389.31506849315065
389.31506849315065
The price in each scenario is:
value.(P)
5-element Vector{Float64}:
59.999999999999886
60.0
59.99999999999994
60.68493150684928
110.68493150684935