Adding complementarity constraints

Specific Domains Optimization (Mathematical)
1

I’m new to both AMPL and JuMP. I’m trying to translate constraints below to jump but don’t understand how to do it.

compl_1{k in K}: 0 <= ll[1,k]  complements    y[1,k] - x[1] <= 0;
compl_2{k in K}: 0 <= ll[2,k]  complements    y[2,k] - x[2] <= 0;
compl_3{k in K}: 0 <= ll[3,k]  complements  - y[1,k] + x[3] <= 0;
compl_4{k in K}: 0 <= ll[4,k]  complements  - y[2,k] + x[4] <= 0;

One way I tried is to

@variable(model, ll[L, K] >= 0) #specifically constraint to be positive?
@constraint(model, [k=K], ll[1,k] ⟂ -y[1,k] + x[1])
@constraint(model, [k=K], ll[2,k] ⟂ -y[2,k] + x[2])
@constraint(model, [k=K], ll[3,k] ⟂  y[1,k] - x[3])
@constraint(model, [k=K], ll[4,k] ⟂  y[2,k] - x[4])

More generally I want to ask how to model constraints of type G_i(x)H_i(x) = 0 \;, G_i(x) \geq 0,\, H_i(x) \geq 0\,\forall x; i=1\ldots n in Julia.

link to AMPL model : http://www.mcs.anl.gov/~leyffer/MacMPEC/problems/design-cent-4.mod

2

Hi @anand.kumar, welcome to the forum :smile:

JuMP supports a particularly form of mixed complementarity problems in which a variable (or vector of variables) x complements a function F(x).

See the documentation: Mixed complementarity problems · JuMP

Those particular constraints would need to be written as:

using JuMP
I = 1:4
J = 1;2
K = 1:3
L = 1:4
model = Model()
@variable(model, x[i in I], start = x0[i])
@variable(model, y[j in J, k in K], start = y0[j, k])
@variable(model, ll[l in L, k in K] >= 0, start = ll0[l, k])
@constraint(model, compl_1[k in K], x[1] - y[1,k] ⟂ ll[1, k])
@constraint(model, compl_2[k in K], x[2] - y[2,k] ⟂ ll[2, k])
@constraint(model, compl_3[k in K], y[1,k] - x[3] ⟂ ll[3, k])
@constraint(model, compl_4[k in K], y[2,k] - x[4] ⟂ ll[4, k])

Note the signs. Since ll >= 0, then the complementing function F(x) is 0 if ll > 0 and F(x) >= 0 if ll == 0.

However, your larger example is an MPEC. Currently, JuMP is missing first-class support for converting MPECs into NLPs, . so you’ll need to use a library like Complementarity.jl:

using JuMP, Ipopt, Complementarity
I = 1:4
J = 1:2
K = 1:3
L = 1:4
x0 = [1, 1, -1, -1]
y0 = zeros(length(J), length(K))
ll0 = ones(length(L), length(K))
model = Model(Ipopt.Optimizer)
@variable(model, x[i in I], start = x0[i])
@variable(model, y[j in J, k in K], start = y0[j, k])
@variable(model, ll[l in L, k in K], start = ll0[l, k])
@objective(model, Max, (x[1] - x[3]) * (x[2] - x[4]))
@constraints(model, begin
    -y[1,1] - y[2,1]^2  <= 0
    y[1,2]/4 + y[2,2] - 3/4 <= 0
    -y[2,3] - 1 <= 0
    1        + ll[1,1] - ll[3,1] == 0
    2*y[2,1] + ll[2,1] - ll[4,1] == 0
    -1/4     + ll[1,2] - ll[3,2] == 0
     -1      + ll[2,2] - ll[4,2] == 0
    0        + ll[1,3] - ll[3,3] == 0
    1        + ll[2,3] - ll[4,3] == 0
end)
for k in K
    @complements(model, 0 >=   y[1,k] - x[1], ll[1,k] >= 0)
    @complements(model, 0 >=   y[2,k] - x[2], ll[2,k] >= 0)
    @complements(model, 0 >= - y[1,k] + x[3], ll[3,k] >= 0)
    @complements(model, 0 >= - y[2,k] + x[4], ll[4,k] >= 0)
end
optimize!(model)
solution_summary(model)
value.(x)
value.(y)
value.(ll)

This gets a trivial solution with zero volume though, so I don’t know if I made a typo or something?

3

Ah. I was missing the correct starting values for y0 and ll0: https://www.mcs.anl.gov/~leyffer/MacMPEC/problems/design-cent-4.dat

y0 = [0 3 0; 0 0 -1]
ll0 = [0 0 0; 0 0 0; 0 0 0; 0 0 1]

Now I get the more reasonable:

julia> solution_summary(model)
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 3.07920e+00
  Dual objective value : 6.17863e+00

* Work counters
  Solve time (sec)   : 3.42009e-02
  Barrier iterations : 47


julia> value.(x)
1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, 1:4
And data, a 4-element Vector{Float64}:
  3.6188022179102477
 -0.15470053521845614
 -0.02393227263092365
 -1.0000000143111953
4

Thanks a ton @odow for quick and elaborate response. There’s lot to learn from your solution.

1 Like
5

Let me know if you have any questions.

I hope to one day add proper MPEC support to JuMP so you dont have to use the special @complements macro from Complementarity.jl. It’s annoying that MCP works nicely but MPEC doesn’t!

1 Like