# Mixed complementarity problems

*This tutorial was generated using Literate.jl.* *Download the source as a .jl file*.

This tutorial is a collection of mixed complementarity programs.

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)
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)
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)
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 utilitiy
@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)
value(K)
```

`0.20474003534537774`

## Electricity consumption

This example is a mixed complementarity formulation of example 3.3.1 from D’Aertrycke, G., Ehrenmann, A., Ralph, D., & Smeers, Y. (2017). Risk trading in capacity equilibrium models.

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 the original text 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)
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) : 2.78000e-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.00000000000014
289.99999999999994
340.0
389.31506849315065
389.31506849315065
```

The price in each scenario is:

`value.(P)`

```
5-element Vector{Float64}:
60.00000000000001
60.0
59.99999999999999
60.68493150684928
110.68493150684935
```