Draft Writeup of April 3, 1998

COMPLEMENTARITY PROBLEMS

A variety of physical and economic phenomena are most naturally modeled by saying that certain pairs of inequality constraints must be complementary, in the sense that at least one must hold with equality. These conditions may in principle be accompanied by an objective function, but are more commonly used to construct complementarity problems for which a feasible solution is sought. Indeed, optimization may be viewed as a special case of complementarity, since the standard optimality conditions for linear and smooth nonlinear optimization are complementarity problems. Other kinds of complementarity problems do not arise from optimization, however, or cannot be conveniently formulated or solved as optimization problems.

A new AMPL operator, complements, permits complementarity conditions to be specified directly in constraint declarations. Complementarity models can thereby be formulated in a natural way, and instances of such models are easily sent to special solvers for complementarity problems.

To motivate the syntax of the new operator, we begin by considering how AMPL might be extended to express the complementary slackness conditions for standard and bounded-variable linear programs. We then give a general definition of the complements operator for pairs of inequalities, and for more general "mixed" complementarity conditions via double inequalities. We also comment on an AMPL interface to Dirkse and Ferris's PATH solver for "square" mixed complementarity problems. Finally, we describe extensions to accommodate complementarity constraints in several of AMPL's existing features: the presolve phase, constraint-name suffixes, and generic synonyms for constraints.

Illustrations are taken from a collection of AMPL complementarity models that use the new syntax. See the complementarity model index for links to the complete models.


Motivation

The requirements for a new complementarity constraint operator can be shown by considering one of the simplest examples: the linear complementarity problem that gives optimality conditions for a linear program. This introduction is likely to be of most value to those who are familiar with the duality theory of linear programming, but not so familiar with the various forms of complementarity problems. Others may want to skip to the next section where the new operator is defined more formally and concisely.

To begin, consider a linear program that has the following elementary primal form: 

minimize PrimalObj: sum {j in J} c[j] * X[j];

subj to PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

subj to PrimalBounds {j in J}: X[j] >= 0;
The corresponding "complementary slackness" conditions for optimality are easily represented like this: 
PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

PrimalBounds {j in J}: X[j] >= 0; 

DualConstr {j in J}: 
   sum {i in I} Y[i] * a[i,j] + Z[j] = c[j]; 

DualBounds {i in I}: Z[i] >= 0; 

Complementarity {j in J}: X[j] = 0 or Z[j] = 0;
Of course, AMPL does not (currently) accept the or operator in constraints, but the Complementarity[j] constraints can be written equivalently in this case as nonlinear equations: 
Complementarity {j in J}: X[j] * Z[j] = 0;
Consider now the case, only slightly more complicated, of a linear program with arbitrary bounds on the variables: 
minimize PrimalObj: sum {j in J} c[j] * X[j];

subj to PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

subj to PrimalBounds {j in J}: l[j] <= X[j] <= u[j];
The corresponding complementary slackness conditions are as follows: 
PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

PrimalBounds {j in J}: l[j] <= X[j] <= u[j]; 

DualConstr {j in J}: 
   sum {i in I} Y[i] * a[i,j] + Z[j] = c[j]; 

Complementarity {j in J}: 

   X[j] = l[j] implies Z[j] >= 0 and 
   X[j] = u[j] implies Z[j] <= 0 and 

   l[j] < X[j] < u[j] implies Z[j] = 0;
The constraints Complementarity[j] are again not recognized by AMPL, but in this case they also lack any concise and straightforward algebraic equivalent.

It would be desirable to be able to express these two similar forms of complementary problem in similar ways. To make this possible, we observe that there are two "elements" comprising each Complementarity[j] constraint: X[j] = 0 and Z[j] = 0 in the first case, and l[j] < X[j] < u[j] and Z[j] in the second. The members of each element pair "complement" each other in a sense appropriate to their form. This suggests defining a new AMPL operator, complements, to explicitly denote the complementarity relationship. The first form of complementary slackness then becomes 

PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

DualConstr {j in J}: 
   sum {i in I} Y[i] * a[i,j] + Z[j] = c[j]; 

Complementarity {j in J}: 
   X[j] >= 0 complements Z[j] >= 0;
and the second form is the same except for the pair of elements that complement each other: 
PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

DualConstr {j in J}: 
   sum {i in I} Y[i] * a[i,j] + Z[j] = c[j]; 

Complementarity {j in J}: 
   l[j] <= X[j] <= u[j] complements Z[j];
To be feasible for this new constraint type, a solution must satisfy the two inequalities that appear, and must cause the two arguments of the complements operator to complement each other. The meaning of "complement each other" for these two cases is exactly as given in our original statements of the complementary slackness conditions above.

The same conditions can be expressed more concisely by substituting the variables Z[j] out of the constraints. We then have for the case of nonnegative variables, 

PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

Complementarity {j in J}: 
   X[j] >= 0 complements 
      c[j] - sum {i in I} Y[i] * a[i,j] >= 0;
and for the case of bounded variables, 
PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

Complementarity {j in J}: 
   l[j] <= X[j] <= u[j] complements 
      c[j] - sum {i in I} Y[i] * a[i,j];
The first form could also be written even more concisely as 
PrimalConstr {i in I}: 
   sum {j in J} a[i,j] * X[j] = b[i]; 

Complementarity {j in J}: 
   X[j] >= 0 complements 
      sum {i in I} Y[i] * a[i,j] <= c[j];
In fact it makes sense to allow the arguments of complements to be any two single-inequality constraints, or any double-inequality constraint and any expression.

It turns out that these few constraint forms are sufficient to conveniently represent most linear complementarity problems of interest, regardless of whether they have any relationship to a linear program. The same is true of nonlinear complementarity problems, when these forms are extended to permit nonlinear constraints and expressions.

We next proceed to define the complements operator more formally and generally, first for the case of two complementary single inequalities, and then for the case of a double inequality complementing an expression.


Complementary pairs of inequalities

A complementarity relation betweeen two inequalities is expressed by placing the complements operator between them in an AMPL constraint declaration. The general form for this constraint type is
single-inequality complements single-inequality ;
where a single-inequality is any valid constraint expression -- linear or nonlinear -- containing one >= or <= operator. A constraint of this type is satisfied if both of the single-inequality constraints are satisfied, and at least one is satisfied with equality.

As an example, pies.mod makes the following statement about coal shipments and prices:

subject to delct {c in CREG, u in USERS}:
   0 <= Ct[c,u] complements ctcost[c,u] + Cv[c] >= P["C",u];
To satisfy this constraint, the variables Ct[c,u], Cv[c], and P["C",u] must satisfy Even though both inequalities are linear in this case, the complementarity condition must be regarded as nonlinear. In fact, the same restrictions are often expressed as the last being explicitly a nonlinear equality. These are all readily written as AMPL arithmetic constraints, without need for the complements operator. They tend to obscure the nature of the complementarity restriction, however, both for people reading the model and for solvers that take advantage of complementarity. They also do not extend well to the more general "mixed" case to be discussed next.


Complementary double inequalities

A "mixed" complementarity condition is specified by placing the complements operator between a double inequality and an expression in an AMPL constraint declaration. The general forms for this constraint type are
double-inequality complements expression ;
expression complements double-inequality ;
where double-inequality is any AMPL constraint expression -- linear or nonlinear -- containing two >= or two <= operators. The conditions for a constraint of this type to be satisfied are as follows: For the special case in which the double-inequality has the form 0 <= body <= Infinity, these conditions reduce to the previously stated conditions for complementarity of a pair of single inequalities.

As an example, pies.mod makes the following statement about coal production:

subject to delc {c in CREG, t in CTYP}:
  0 <= C[c,t] <= cmax[c,t] complements
  ccost[c,t] + (sum {res in R} cruse[res,c,t] * Mu[res]) - Cv[c];
This implies that the following conditions must be satisfied: There is no equally natural way to state these conditions in terms of arithimetic constraints alone; at the least, a transformation involving additional variables would be required.

For completeness, the special case in which the left-hand side equals the right-hand side of the double inequality may be written using one of the forms

equality complements expression ;
expression complements equality ;
A constraint of this kind is equivalent to an ordinary constraint consisting only of the equality; it places no restrictions on the expression.


Using the PATH solver

PATH, a mixed complementarity solver, takes as its input a "square" complementarity system consisting of m double-inequality complementarity constraints in m variables. An AMPL/PATH interface has been developed to accept a broader variety of AMPL constraint systems that can be converted to the form that PATH requres.

Specifically, AMPL/PATH accepts any constraint system such that the number of variables equals the number of complementarity constraints plus the number of equality constraints. There may be any number of additional inequality constraints, but there must not be any objective. The AMPL/PATH interface automatically transforms constraints from this form to the more resricted form required by PATH, and later inverts the transformation so that the results may be viewed in terms of the originally specified model.

AMPL/PATH is invoked by setting AMPL's solver option to path prior to a solve command:

ampl: model pies.mod;
ampl: data pies.dat;
ampl: option solver path;
ampl: solve;
PATH 3.0: Solution found.
14 iterations (1 for crash); 28 pivots.
30 function, 16 gradient evaluations.
ampl:
An error message such as 
Nonsquare complementarity system:
	38 equations and 42 variables.
PATH 3.0 requires a square system; this one is 38 x 42.
indicates that the specified problem instance cannot be solved by PATH, because the number of variables does not equal the total of equality and complementarity constraints.


Presolve

AMPL incorporates a presolve phase that can substantially simplify some linear programs. By use of an iterative algorithm, presolve deduces tighter bounds on variables and constraint expressions. It may infer from these bounds that certain variables can be fixed and constraints dropped. In the presence of complementarity constraints, several new kinds of inferences become possible.

As an example, given a constraint of the form

expr1 >= 0 complements expr2 >= 0,
if presolve can deduce that expr1 is strictly positive for all feasible points -- in other words, that it has a positive lower bound -- then it can replace the constraint by expr2 = 0. Similarly, in a constraint of the form
expr1 <= expr2 <= expr3 complements expr4,
there are various possibilities, including the following:

If presolve can deduce
for all feasible points that 		Then the constraint can be replaced by
expr2 < expr3 expr1 <= expr2 complements expr4 >= 0
expr1 < expr2 < expr3 expr4 = 0
expr4 < 0 expr2 = expr3

Transformations of these kinds are carried out automatically, unless you use option presolve 0 to turn off the presolve phase. As is currently the case for ordinary constraints, results are reported in terms of the original model, so that the benefits of presolving are achieved without any special attention from the user.

By displaying a few predefined parameters,

_ncons      the number of ordinary constraints before presolve,
_nccons     the number of complementarity conditions before presolve,
_sncons     the number of ordinary constraints after presolve,
_snccons    the number of complementarity conditions after presolve,
or by setting option show_stats 1, you can get some information on the number of simplifying transformations that presolve has made: 
ampl: model pies.mod;
ampl: data pies.dat;
ampl: option solver path;
ampl: option show_stats 1;
ampl: solve;

Presolve eliminates 42 constraints.
Presolve resolves 8 of 42 complementarity conditions.
Adjusted problem:
42 variables:
	6 nonlinear variables
	36 linear variables
42 constraints; 142 linear nonzeros
	34 nonlinear constraints
	8 linear constraints
34 complementarity conditions among the constraints:
	28 linear, 6 nonlinear.
0 objectives.

PATH 3.0: Solution found.
14 iterations (1 for crash); 28 pivots.
30 function, 16 gradient evaluations.

ampl: display _ncons, _nccons;
_ncons = 84
_nccons = 42

ampl: display _sncons, _snccons;
_sncons = 42
_snccons = 34
When first instantiating the problem, AMPL counts each complementarity constraint as two ordinary constraints (the two arguments to complements) and also as a complementarity condition. Thus _nccons equals the number of complementarity constraints before presolve, and _ncons equals twice _nccons plus the number of non-complementarity constraints before presolve. The presolve messages at the beginning of the show_stats output indicate how much presolve was able to reduce these numbers.


Auxiliary solution values

Solution information associated with an AMPL constraint can be represented by expressions of the form cname.suf, where cname is any constraint name and suf is one of several predefined suffixes. For a constraint in the general form
lbound <= body <= ubound,
where lbound and ubound are constants, the recognized suffixes and their values are:
cname.body      body
cname.lb        lbound
cname.ub        ubound
cname.lslack    body - lbound
cname.uslack    ubound - body
cname.slack     min (cname.lslack, cname.uslack)
Any ordinary constraint can be put into this form, possibly with lbound equal to -Infinity or ubound equal to Infinity.

This feature extends to complementarity constraints, but with two collections of suffixes, of the form cname.Lsuf and cname.Rsuf, corresponding to the left and right operands of complements, respectively. Thus for example after pies.mod has been solved, you can use the following display command to look at values associated with the constraint delc: 

ampl: display delc.Llb, delc.Lbody, delc.Lub,
ampl?         delc.Rlb, delc.Rbody, delc.Rub;

:   delc.Llb delc.Lbody delc.Lub   delc.Rlb   delc.Rbody   delc.Rub :=
1 1     0      300         300     -Infinity   -10.7551    Infinity
1 2     0      300         300     -Infinity    -9.51119   Infinity
1 3     0      227.889     400     -Infinity    -8         Infinity
2 1     0      200         200     -Infinity   -10.5051    Infinity
2 2     0      300         300     -Infinity    -8.91133   Infinity
2 3     0      600         600     -Infinity    -8.46915   Infinity
;
Subscripted forms of the constraint name can also be suffixed: 
ampl: display {t in ctyp} (delc[1,t].Lslack, delc[1,t].Rslack);

: delc[1,t].Lslack  delc[1,t].Rslack  :=
1        0               Infinity
2        0               Infinity
3      172.111           Infinity
;
Notice that since the right operand of delc is an expression, it is treated as a "constraint" with infinite lower and upper bounds, and hence infinite slack.

The suffix cname.slack is also defined for complementarity constraints. For complementary pairs of single inequalities, it is equal to min(cname.Lslack, cname.Rslack). Hence it is nonnegative if and only if both inequalities are satisfied. For complementary double inequalities of the form

expr complements lbound <= body <= ubound
lbound <= body <= ubound complements expr
cname.slack is defined to be

min(expr, body - lbound) if body <= lbound
min(-expr, ubound - body) if body >= ubound
-abs(expr) otherwise

Hence in this case it is always nonpositive, and is zero when one part of the double inequality is satisfied exactly.

If cname for a complementarity constraint appears unsuffixed in an expression, it is interpreted as representing cname.slack.


Generic synonyms

AMPL's generic synonyms for constraints have been extended to encompass complementarity. The generic synonyms for ordinary (equality and inequality) constraints remain the same, while a parallel collection of synonyms for complementarity constraints is constructed mainly by substituting ccon for con where appropriate.

From the modeler's view (before presolve), the ordinary constraint synonyms remain

_ncons     the number of ordinary constraints before presolve,
_conname   names of the ordinary constraints before presolve,
_con       synonyms for the ordinary constraints before presolve,
The new complementarity constraint synonyms are
_nccons     the number of complementarity constraints before presolve,
_cconname   names of the complementarity constraints before presolve,
_ccon       synonyms for the complementarity constraints before presolve,
Because each complementarity constraint also gives rise to two ordinary constraints, as explained in the preceding discussion of presolve, there are two entries in _conname corresponding to each entry in _cconname: 
ampl: model josephy.mod;
ampl: data josephy.dat;

ampl: display _conname,_cconname;

:   _conname _cconname    :=
1   'f[1].L'   'f[1]'
2   'f[1].R'   'f[2]'
3   'f[2].L'   'f[3]'
4   'f[2].R'   'f[4]'
5   'f[3].L'      .
6   'f[3].R'      .
7   'f[4].L'      .
8   'f[4].R'      .
;

For each complementarity constraint cname, the left and right arguments to the complements operator are the ordinary constraints named cname.L and cname.R. You can see this by using the synonym terminology to expand complementarity constraint f[1] and the corresponding two ordinary constraints from the example above: 

ampl: expand f[1];
s.t. f[1]:
	-x[1] <= 0
     complements
	3*x[1]*x[1] + 2*x[1]*x[2] + 2*x[2]*x[2] + x[3] + 3*x[4] >= 6;

ampl: expand {i in 1..2} _con[i];
s.t. f[1].L:
	-x[1] <= 0;

s.t. f[1].R:
	3*x[1]*x[1] + 2*x[1]*x[2] + 2*x[2]*x[2] + x[3] + 3*x[4] >= 6;
From the solver's view (after presolve), a more limited collection of synonyms is defined:
_sncons     the number of all constraints after presolve,
_snccons    the number of complementarity constraints after presolve,
_sconname   names of all constraints after presolve,
_scon       synonyms for all constraints after presolve,
Necessarily _snccons is less than or equal to _sncons, with equality only when all constraints are complementarity constraints. By using solexpand in place of expand, you can see the form in which AMPL sent complementarity constraint f[1] to the solver: 
ampl: solexpand f[1];
s.t. f[1]:
	3*x[1]*x[1] + 2*x[1]*x[2] + 2*x[2]*x[2] + -6 + x[3] + 3*x[4] >= 0
     complements
	0 <= x[1];
To simplify the problem description that is sent to the solver, AMPL converts every complementarity constraint into one of the following canonical forms,
expr complements lbound <= var <= ubound,
expr <= 0 complements var <= ubound,
expr >= 0 complements lbound <= var,
where var is the name of a different variable for each constraint. A predefined array of integers, _scvar, gives the indices of these canonical complementing variables in the generic variable arrays _var and _varname. This terminology can be used to display a list of names of such variables: 
ampl: display {i in 1.._sncons} _varname[_scvar[i]];

_varname[_scvar[i]] [*] :=
1  'x[1]'
2  'x[2]'
3  'x[3]'
4  'x[4]'
;
When constraint i is an ordinary equality or inequality, _scvar[i] is 0.


Comments or questions?
Write to info@ampl.com or use our comment form.

Return to the AMPL update page.

Return to the AMPL home page.

LAST MODIFIED 3 APRIL 1998 BY 4er.