Complementarity constraint

[blegat]

We discussed with @frapac how to create a complementarity constraint type.
One idea would be to add a type

struct Complementarity <: MOI.AbstractVectorFunction end
MOI.dimension(::Complementarity) = 2

And the constraint [x, y]-in-Complementarity would constraint x, y to be nonnegative and x * y to be zero.

[odow]

@xhub, @chkwon: you probably have some views on this.

At the JuMP level, it'd be nice to see something like this

@constraint(model, 0 <= F(x) ⟂ x >= 0)

[frapac]

The idea is to load directly complementarity constraints to the solver if it supports them.
For instance:
https://github.com/JuliaOpt/KNITRO.jl/blob/master/src/kn_constraints.jl#L362-L378

As far as I understand, Complementarity.jl reformulates complementarity constraints as NLP constraints. Once NLP would become first class in MOI, we may want to write bridges to reformulate complementarity constraints as NLP ones.

[mlubin]

SGTM. "How do I add complementarity constraints to JuMP?" was one of the original motivators for MOI.

[chkwon]

This will be great. In fact, I'm not so happy with the reformulations strategy currently in place with Complementarity.jl. A solver-native approach would be awesome.

By the way, I have used something like the following in Complementarity.jl

@complements(m, 0 <= F(x),   x >= 0)

A special symbol ⟂ is nice, but also we need a traditional ASCII approach as well. Perhaps:

@constraint(m, 0 <= F(x) complements   x >= 0)

[matbesancon]

it may be reformulated differently in MILP-based solvers, or solvers with better support for integrality than non-linear constraints, using SOS1 or indicator constraints

[xhub]

Sorry for giving my input so late, was too busy with ICCOPT 2019, and then I first got myself more familiar with MOI by porting my solver ReSHOP to MOI.

TLDR I would like to push for supporting more than just "complementarity constraint" (CC) like 0 <= F(x) ⟂ x >= 0, but LCP, MCP, VI (variational Inequalities), QVI (quasi-variational inequalities), MPCC and MPEC.

I started to type a design document 2 years ago for MBP, but right at the time the transition to MOI was started, so I will update it. I get the feeling that more problem classes can be handle through MOI.

I have a lightly different take on the kind of relations to support and the UI, which I will try to motivate by some perceived shortcomings of the proposed syntax. There is quite a few occurrences of the VI lingo, I hope it is not too off-putting.

Discussion on the "natural interface" @constraint(model, 0 <= F(x) ⟂ x >= 0)

With complementarity constraint, a salient point is correctness of the formulation, and the drawbacks with the kind of a "natural interface" like @constraint(model, 0 <= F(x) ⟂ x >= 0) include the following :

1. Nonsensical specification:@constraint(model, 0 == F(x) ⟂ x >= 0) is not a complementarity constraint.

2. Incorrect specification

@constraint(model, 0 <= F(x) ⟂ x >= 0)
@constraint(model, x <= 3)

The feasible set for x defines a rectangular VI (or MCP à la Ferris): F(x) ⟂ 0 <= x <= 3
Which case by case, means that:

• F(x) >= 0 whenever x = 0
• F(x) = 0 whenever 0 < x < 3
• F(x) <= 0 whenever x = 3

3. Inconsistent syntax: For a NCP, we would have a signature like @constraint(model, 0 <= F(x) ⟂ x >= 0) whereas that an MCP, it would be @constraint(m, F(x) ⟂ 0 <= x <= 3)

It's important to realize that in term of complementarity problems, the dual of set in which x belongs constraints is where F(x) belongs:
Hence, there is no need to specify twice the same information, also reducing the possible mistakes.

Here are 2 proposed ways to input complementarity relations

A) Popular CC

• Allow definition of rectangular VI via @constraint(m, F(x) ⟂ lb <= x <= ub), via ub and lb taken as +Inf, -Inf if not defined.
• An error is thrown if the user does something that breaks the type of constraint that he entered. Basically after the definition @constraint(m, 0 F(x) ⟂ lb <= x <= ub) sets x in stone, and no other constraints on x is allowed.
• If x was constrained before the CC was added, an error must be thrown.

Most of the work is in enforcing the "freezing" of x. Would it be done by JuMP, or left for the solver to be done? It being in JuMP would be nice to avoid duplicating this feature. On the MOI side, I am also not sure how to get from a VariableIndex all the ConstraintIndex where this variable is present. Maybe some solvers have this check internally.

B) General purpose complementarity specifications

• Provide an interface that support a much larger class of problem types: @complements(m, F, x)
  This would mean that the user doesn't specify directly a CC like A), but rather that there is a specific relation (named matching or complements) between a variable and a mapping.
• The actual type of constraint or problem, is not known until the optimize phase. The set where x belongs, is determined by looking at the the constraints involving only x. Then the proper problem class (NCP, MCP, VI, MPEC or MPCC) can be inferred. This is akin to an implicit splitting of the constraints in the model.

I favor B over A, but I have a limited view of the way other solvers (beyond ReSHOP and PATH) do support complementarity in models. I understand that both departs what most people are (unfortunately) taught about CC in the way that it is a relation where 0 <= F(x), x>=0, and F(x)^T x = 0, and that a majority of users may want to specify that kind of constraints. I hope that a good documentation/tutorial would clarify the point that there is no need for an explicit constraint on F.

MOI/JuMP rough implementation sketch for B

For B, at the MOI level, I would define a constraint [F, x]-in-generalized-equation, or [F, x]-in-VI that represent the generalized equation (another way of representing Variational Inequalities

where N_C is the normal cone from convex analysis,

Whenever C is closed convex cone, this is a (generalized) complementarity problem:

When C is the positive orthant or the whole space, one recovers "classical" complementarity condition.

At the JuMP level, I like the @complements(m, F, x) or @constraints(m, F ⟂ x), even @constraints(m, F ⟂ lb <= x <= ub) as a syntactic sugar that would translate into adding a constraint to x and setting the relationship between F and x. If the user specifies

@constraints(m, F ⟂ lb <= x)
@constraints(m, x <= ub)

This is understood as an MCP F ⟂ lb <= x <= ub

Quasi-Variational Inequalities (QVI)

For the Quasi-Variational Inequalities (QVI), one way to model them it to introduce a 3rd type of variable, y, such that the problem definition reads

At the solution, we require y = x. This would be covered at the MOI level by [F, x, y]-in-QVI. At the JuMP level, an interface of the kind @complements(F, x, y) or @QVI(F, x, y) could be considered

There are also opportunity for problem transformations (not an equivalence), MCP => NCP; VI => MCP; in addition to the aforementioned NCP=> binary/integer variables. In general, the first two should be avoided if the solver supports the original problem types, since the reformulated problem may not be processable. For instance, any VI over a compact set with a continuous functional has a solution (simple application of Browder FP theorem). The reformulation as a NCP may destroy that structure. This has been observed in practice over non-monotone problems.

And also relaxation MPCC => NCP, ..., have seen a lot of attentions in the past (NLPEC, work by S. Leyffer, ...).

In terms of solvers, AFAIK Baron does automatic detection. Knitro needs the mapping to be transformed into a defining equation for an auxiliary variable, which are then used in the CC definition.

[odow]

Thanks for the input @xhub!

I believe Michaels' preference is exactly your option B): complements(model, F, x), where F is a scalar (nonlinear) function and x is a variable, and that the bounds on x define the VI. (And for exactly the reasons you describe... (unsurprisingly).)

So we would define

struct Complements <: MOI.AbstractVectorSet end

Solvers (e.g., PATH.jl) would support

function MOI.add_constraint(model, ::MOI.VectorAffineFunction, ::MOI.Complements)

and throw an error is dimension(::MOI.VectorAffineFunction) != 2 or if the second dimension was not a straight variable.

Then, bridges would map people writing

@constraint(model, [y, x] in Complements())

into

@constraint(model, [1.0 y, x] in Complements())

JuMP can't (currently) handle

@constraint(model, [sin(y), x] in Complements())

because it doesn't have an AbstractScalarNonlinearFunction type. But this would be worked around by

@variable(model, z)
@NLconstraint(model, z == sin(y)
@constraint(model, [y, x] in Complements())

Handling quasi-variational inequalities seems harder, and a stretch goal. You could define

@constraint(model, [F(x), x, y] in QuasiVariationalInequality())

But the number of solvers supporting this seems minimal. It might be a good function/set combo for Path-VI to support until we work out the details, rather than add to MOI from the outset.

[xhub]

Would the following also work (inspired by SOS(1|2)) ?

struct Complements <: AbstractVectorSet
    variable::Union{MOI.VariableIndex,MOI.VectorOfVariables}
end

This would also allow for the use of vectors in the definition of CC.
I also like the semantics of having the variable and the set begin tied together.

The canonical signature would then be

function MOI.add_constraint(model, S, ::MOI.Complements) where S <: Union{MOI.VariableIndex, VOV, SF, VAF, VQF}

On the JuMP side, we would have

@constraint(model, y in Complements(x))
@constraint(model, x*x - 1 in Complements(x))

For the NLP workaround, it's a bit sad, I have read #846 . In EMP.jl I used some macro @NLvipair. Maybe we could have a nicer but hackisch @CCconstraint(m, func, x) which would translate into

@variable(m, z)
@NLconstraint(m, func == z)
@constraint(m, z in ComplementsNL(x))

where ComplementsNL <: Complements. I would push for a different signature for the VI/CP case.
Indeed, in that case if z were to stay as a variable, it would have to be matched with the constraint z = sin(x). For VI/CP, all variables must be defined has matching some function.

Furthermore, the variable z should be completely eliminated from the VI/CP since it most likely makes the feasible set nonconvex. Having a dedicated Set would be nice. I understand that one may not want to

For MPCC, the situation may be different, since for instance KNITRO requires the CC constraints to be between 2 variables. It seems that they would not need to differentiate.

Finally for QVI, how about a

struct QuasiVariationalInequalitySet <: AbstractVectorSet
    x::Union{MOI.VariableIndex,MOI.VectorOfVariables}
    y::Union{MOI.VariableIndex,MOI.VectorOfVariables}
end

That would allow for

@constraint(model, F(x) in QuasiVariationalInequality(x, y))

[odow]

I started an interface to PATH here: https://github.com/odow/PATH.jl.

(It's intended to replace https://github.com/chkwon/PATHSolver.jl, which uses an additional C wrapper instead of going directly into the libpath shared object, and has some funky ways of dealing with callbacks which are an artifact of old Julia code.)

My next step is to write the MOI interface. Talking with Michael, we will support VectorAffineFunction -in- Complements and VectorQuadraticFunction-in-Complements. At the JuMP level, you. will write

@constraint(model, [i=1:2], [i * x[1] + x[2], x[i]] in Complements())

You could also make F(x) a quadratic.

This is a little ungainly, so I might try to implement

@constraint(model, [i = 1:2], i * x[1] + x[2] ⟂ x[i])

Alternatively, we might take advantage of jump-dev/JuMP.jl#2051, and write

@constraint(model, [i=1:2], complements(i * x[1] + x[2], x[i]))

In any case, the JuMP syntax is purely cosmetic. The MOI syntax is the one required by the solvers.

Would the following also work (inspired by SOS(1|2)) ?

No. One limitation with the current MOI implementation is that we really don't want the sets to contain variables because this violates a lot of the assumptions made by bridges etc.

For the NLP workaround, it's a bit sad

I agree. In regard to nonlinear: I retract my previous work-around because you're right, it might destroy convexity/monotonicity.

We won't be targeting nonlinear equilibrium problems until MOI has a ScalarNonlinearFunction and VectorNonlinearFunction function. Then it should be seamless to write nonlinear F(x).

For MPCC, the situation may be different, since for instance KNITRO requires the CC constraints to be between 2 variables.

Does KNITRO allow

@constraint(model, z == x + y)
@constraint(model, [z, x] in Complements())

If so, it can just support VectorOfVariables-in-Complements and throw an error if the dimension is not 2.

[chkwon]

@odow A direct Julia interface to PATH is great! When PATHSolver.jl was written, installing PATH libraries was a bit of headache. Hope it will work out smoothly this time. By the way, they will not allow the name PATH.jl to be registered in METADATA.jl.

[odow]

installing PATH libraries was a bit of headache.

Yes, hopefully this will just require installing libpathXX, rather than compiling the Standalone_C interface. I believe we can use https://github.com/ampl/pathlib/tree/master/lib directly, although I have been developing for the latest Path 5.0.0. Michael was going to look into the API differences.

By the way, they will not allow the name PATH.jl to be registered in METADATA.jl.

Ah. Good to know. I guess PATH is a little generic.

[frapac]

@odow I confirm that Knitro would allow your formulation.

[xhub]

@frapac Based on the comments in https://github.com/JuliaOpt/KNITRO.jl/blob/master/examples/mpec1.jl I infered that it even requires that formulation. Am I missing something?

[xhub]

No. One limitation with the current MOI implementation is that we really don't want the sets to contain variables because this violates a lot of the assumptions made by bridges etc.

Ok that's unfortunate. Is it something that is planned on being worked on?

Sorry if I sound pushy, it would be great to be able to write

@constraint(model, A*x + b ⟂ x)

Could we have a

struct ComplementsPairing{F,V} <: AbstractFunction
    func::F
    variable::V
end

Then we can allow for scalar and vector use of Complements. The MOI signature would then be

function MOI.add_constraint(model, ::MOI.ComplementsPairing{F,V}, ::MOI.Complements)

[odow]

Sorry if I sound pushy, it would be great to be able to write
@constraint(model, A*x + b ⟂ x)

Not pushy at all. But, we should distinguish between the JuMP syntax, and the MOI syntax.

The immediate question is what the MOI syntax should be, and I favor VectorAffineFunction-in-Complements, VectorQuadraticFunction-in-Complements, and the hypothetical VectorNonlinearFunction-in-Complements, where each constraint gives the i'th element F_i(x) complements x_i.

The JuMP syntax you wrote looks great, and is what we should aim for. JuMP would lower that syntax into the vector of the complements constraints.

[odow]

Progress:

julia> using JuMP, PATH

julia> model = Model(with_optimizer(PATH.Optimizer))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SolverName() attribute not implemented by the optimizer.

julia> @variable(model, x >= 0)
x

julia> @constraint(model, [x + 2, x] in PATH.Complements())
[x + 2, x] ∈ PATH.Complements()

julia> optimize!(model)
Path 5.0.00 (Mon Aug 19 10:57:18 2019)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Crash Log
major  func  diff  size  residual    step       prox   (label)
    0     0             1.5279e+00             0.0e+00 (f[    1])
    1     1     1     1 0.0000e+00  1.0e+00    0.0e+00 (f[    1])

Major Iteration Log
major minor  func  grad  residual    step  type prox    inorm  (label)
    0     0     2     2 0.0000e+00           I 0.0e+00 0.0e+00 (f[    1])

Major Iterations. . . . 0
Minor Iterations. . . . 0
Restarts. . . . . . . . 0
Crash Iterations. . . . 1
Gradient Steps. . . . . 0
Function Evaluations. . 2
Gradient Evaluations. . 2
Basis Time. . . . . . . 0.000015
Total Time. . . . . . . 0.132081
Residual. . . . . . . . 0.000000e+00

julia> value(x)
0.0

julia> termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4

[blegat]

Could we have a

struct ComplementsPairing{F,V} <: AbstractFunction
    func::F
    variable::V
end

Then we can allow for scalar and vector use of Complements.

This function would not only be useful for complementarity constraints but also for indicator constraints. In fact, it was already proposed in #709 (comment).
However, we decided not to pursue this as we would like to keep a small number of different MOI functions (that depends on MOI variables) and a large number of sets (that do not depend on any MOI indices).
So unless it really allows to do something more, we prefer not adding new function type.
Here to allow for vectors, we could say that in func-in-Complement(set), we should have MOI.output_dimension(func) == 2MOI.dimension(set) and F is given by the first half of func and x is given by the second half of func (error thrown if it's not convertible to a VectorOfVariables). This would constraint F(x) to belong to the normal cone of set at x. [F(x), x]-in-Compements(MOI.Nonnegatives(1)) would be the classical complementarity constraint that @odow just implemented for PATH.

[xhub]

The immediate question is what the MOI syntax should be, and I favor VectorAffineFunction-in-Complements, VectorQuadraticFunction-in-Complements, and the hypothetical VectorNonlinearFunction-in-Complements, where each constraint gives the i'th element F_i(x) complements x_i.

I would favour a tailored data structure where the functional part and the variable part are clearly separated. Just by looking at the definition of ComplementsPairing it is clear how to get the data. Typed structure are clearer than structure by convention. Also in the nonlinear case, it would be harder to extract the variable part (need to ensure that a Expr is just a variable).

Furthermore, it would require the type VectorNonlinearFunction to exists, which does not currently exists. I may have misunderstood the vector NLP discussion recently, but I got the feeling that vector nonlinear functions support in JuMP would need a major rewrite. I am concerned that the complementarity would be the only consumer of a VectorNonlinearFunction in MOI. Would that impact the likelihood of having the type defined, and hence, having NCP support in MOI?

I understand a drawback is to have to support another type, with the risk of having to depreciate it.

The JuMP syntax you wrote looks great, and is what we should aim for. JuMP would lower that syntax into the vector of the complements constraints.

What would be cost of having to lower it vs just building a data structure? I'm just wondering :)

[xhub]

@blegat Thanks for your input, I guessed right about the resistance of adding new types, which totally makes sense.

I will try to further rest my case for a data type of a tuple (variable, function). One may wish to write some utilities to extract from a VAF a tuple (variable, function). Convention are usually based on the semantics of the applications, which can be inconsistent.
Here, there is an example of this phenomenon, since Complementarity and indicator constraints could use the same data structure. However one could not write a utility to extract the variables and functions from a (VAF|VQF) for both type of constraint, since the semantics of the constraints favour opposite conventions: for indicator function, it is more natural to have [x, F], and for complementarity it is [F, x]. So any kind of common utility functions to get the variable and function out of a (VAF|VQF) would have to handle the 2 orders, or the code is duplicated.

It's true that such a data structure is not necessary in "classical" mathematical programming. It would also be useful in disjunctive programming, where it is also necessary to associate constraints with a variable. It would also work to encode those relationship by convention.

I will try to think about this data strcture allowing something more. I doubt I will come up with anything, since a vector function are always capture a tuple (var, func).

[mlubin]

@xhub To be a bit more concrete, defining a new MOI function is about 50x more effort than defining a new set given how the infrastructure has been designed.

So any kind of common utility functions to get the variable and function out of a (VAF|VQF) would have to handle the 2 orders, or the code is duplicated.

Relatively speaking, this is a really tiny amount of code.

[odow]

So any kind of common utility functions to get the variable and function out of a (VAF|VQF) would have to handle the 2 orders, or the code is duplicated.

This is exactly what I do in PATH. A rough implementation without errors that need improving is here: https://github.com/odow/PATH.jl/blob/f15f931c64bfd4a81cfb3aeaac51214ca8885a2d/src/MOI_wrapper.jl#L82-L138

What would be cost of having to lower it vs just building a data structure? I'm just wondering :)

Not sure. But probably << than the cost of a PATH solve. I'd be surprised if it was the bottleneck.

[odow]

So I think we should define:

struct Complements <: MOI.AbstractVectorFunction
    dimension::Int
end

Then we require that the vector valued function in VectorFunction-in-Complements(dim) has dimension 2 * dim, where the ith component of the function complements the i+dim/2 component.

Then you could go:

model = Model()
@variable(model, x[1:2] >= 0)
@variable(model, y[1:2])
@constraint(model, [M * x .+ q; x] in Complements(2))
@constraint(model, [y; x] in Complements(2))

Edit: this syntax is now implemented in PATH.jl: https://github.com/odow/PATH.jl#example-usage
I think it's quite nice.

[odow]

And here's the transmcp example from GAMS in JuMP: https://github.com/odow/PATH.jl/blob/master/examples/transmcp.jl

Current syntax forces us to write:

@constraints(model, begin
    vec([profit x]) in PATH.Complements(length(x))
    vec([supply w]) in PATH.Complements(length(w))
    vec([fxdemand p]) in PATH.Complements(length(p))
end)

If we added the Complements set to JuMP, it could automatically detect size. It would also be nice if we could automatically vec matrices into a vector. Although dealing with [profit x] and [profit; x] and [profit, x] might be annoying.

[xhub]

Sounds good, you guys all know better how much work it is to introduce new functions, and how MOI was intended to work.

struct Complements <: MOI.AbstractVectorFunction
    dimension::Unsigned
end

sounds like a very good idea.

From a modeler standpoint, we should support [profit x] and [profit; x] and [profit, x].

Or we would only expose @constraint(model, A*x + b ⟂ x) and @constraint(model, A*x + b perp x) or @constraint(model, A*x + b complements x) as ascii variant. This way we avoid dealing with such vector/matrix issues. Again, no idea of the cost to support such things.

[blegat]

Again, no idea of the cost to support such things.

This will be at parse time in JuMP.parse_constraint so it's equivalent to doing [profit; x] as far as runtime is concerned.