Complementarity constraints

A complementarity constraint enforces that two variables are complementary to each other; i.e., that the following conditions hold for scalar variables x and y:

x \cdot y = 0, \quad x \ge 0, \quad y \ge 0.

The condition above is sometimes expressed more compactly as

0 \le x \perp y \ge 0.

Intuitively, a complementarity constraint is a way to model a constraint that is combinatorial in nature since, for example, the complementary conditions imply that either x or y must be 0 (both may be 0 as well).

Without special care, these types of constraints may cause problems for nonlinear optimization solvers because problems that contain these types of constraints fail to satisfy constraint qualifications that are often assumed in the theory and design of algorithms for nonlinear optimization. For this reason, we provide a special interface in Knitro for specifying complementarity constraints. In this way, Knitro can recognize these constraints and handle them with special care internally.

Note

The complementarity features of Knitro are not available through all interfaces. Currently, they are accessible only to users of the callable library, the MATLAB interface, and some modeling environments such as AMPL.

If a modeling language does not allow you to specifically identify and express complementarity constraints, then these constraints must be formulated as regular constraints and Knitro will not perform any specializations.

Note

There are various ways to express complementarity conditions, but the complementarity features in the Knitro callable library API and MATLAB API require you to specify the complementarity condition as two non-negative variables complementary to each other as shown above. Any complementarity condition can be written in this form.

Example

This problem is taken from J.F. Bard, Convex two-level optimization, Mathematical Programming 40(1), 15-27, 1988.

Assume we want to solve the following MPEC with Knitro.

\min \quad               & f(x) = (x_0 - 5)^2 + (2 x_1 + 1)^2 \\
\mbox{subject to:} & \\
                                        & c_0(x) = 2(x_1 - 1) - 1.5 x_0 + x_2 - 0.5 x_3 + x_4 = 0 \\
                                        & c_1(x) = 3 x_0 - x_1 - 3 \ge 0 \\
                                        & c_2(x) = -x_0 + 0.5 x_1 + 4 \ge 0 \\
                                        & c_3(x) = -x_0 - x_1 + 7 \ge 0 \\
                                        & c_1(x) \cdot x_2 = 0 \\
                                        & c_2(x) \cdot x_3 = 0 \\
                                        & c_3(x) \cdot x_4 = 0 \\
                                        & x_i \ge 0 \quad \forall i = 0, \dots 4.

Observe that complementarity constraints appear. Expressing this in compact notation, we have:

\min \quad               & f(x) = (x_0 - 5)^2 + (2 x_1 + 1)^2 \\
\mbox{subject to:} & \\
                                        & 2(x_1 - 1) - 1.5 x_0 + x_2 - 0.5 x_3 + x_4 = 0 \quad (c_0) \\
                                        & c_1(x) = 3 x_0 - x_1 - 3 \\
                                        & c_2(x) = -x_0 + 0.5 x_1 + 4 \\
                                        & c_3(x) = -x_0 - x_1 + 7 \\
                                        & 0 \leq c_1(x) \perp x_2 \geq 0 \\
                                        & 0 \leq c_2(x) \perp x_3 \geq 0 \\
                                        & 0 \leq c_3(x) \perp x_4 \geq 0 \\
                                        & x_0, x_1 \ge 0.

Since Knitro requires that complementarity constraints be written as two variables complementary to each other, we must introduce slack variables (x_5, x_6, x_7) and re-write the problem as follows:

\min \quad               & f(x) = (x_0 - 5)^2 + (2 x_1 + 1)^2 \\
\mbox{subject to:} & \\
                                        & 2(x_1 - 1) - 1.5 x_0 + x_2 - 0.5 x_3 + x_4 = 0 \quad (c_0) \\
                                        & 3 x_0 - x_1 - 3 - x_5 = 0 \quad(c_1) \\
                                        & -x_0 + 0.5 x_1 + 4 - x_6 = 0 \quad (c_2) \\
                                        & -x_0 - x_1 + 7 - x_7 = 0 \quad (c_3) \\
                                        & 0 \leq x_5 \perp x_2 \geq 0 \\
                                        & 0 \leq x_6 \perp x_3 \geq 0 \\
                                        & 0 \leq x_7 \perp x_4 \geq 0 \\
                                        & x_i \ge 0, \quad \forall i = 0, \dots 7..

The problem is now in a form suitable for Knitro.

Complementarity constraints in AMPL

Complementarity constraints should be modeled using the AMPL complements command; e.g.,:

0 <= x complements y => 0;

The Knitro callable library API and MATLAB API require that complementarity constraints be formulated as one variable complementary to another variable (both non-negative). However, in AMPL (beginning with Knitro 8.0), you can express the complementarity constraints in any form allowed by AMPL. AMPL will then translate the complementarity constraints automatically to the form required by Knitro.

Be aware that the AMPL presolver sometimes removes complementarity constraints. Check carefully that the problem definition reported by Knitro includes all complementarity constraints, or switch off the AMPL presolver by setting option presolve to 0, if you don’t want the AMPL presolver to modify the problem.

Complementarity constraints in MATLAB

Complementarity constraints can be specified through two fields of the extendedFeatures structure. The fields ccIndexList1 and ccIndexList2 contain the pairs of indices of variables that are complementary to each other.

Note

Variables which are specified as complementary should be specified to have a lower bound of 0 through the variable lower bound array lb.

Complementarity constraints with the callable library

Complementarity constraints can be specified in Knitro through a call to the function KTR_addcompcons() which has the following prototype:

int  KNITRO_API KTR_addcompcons (KTR_context_ptr    kc,
                                 const int          numCompConstraints,
                                 const int * const  indexList1,
                                 const int * const  indexList2);

In addition to kc, which is a pointer to a structure that holds all the relevant information about a particular problem instance, the arguments are:

  • numCompConstraints, the number of complementarity constraints to be added to the problem (i.e., the number of pairs of variables that are complementary to each other).
  • *indexList1 and *indexList2, two arrays of length numCompConstraints specifying the variable indices for the first and second sets of variables in the pairs of complementary variables.

Note

The call to KTR_addcompcons() must occur after the call to KTR_init_problem(), but before the first call to KTR_solve().

Note

Variables which are specified as complementary through the special KTR_addcompcons() functions should be specified to have a lower bound of 0 through the Knitro lower bound array xLoBnds.

Complementarity constraints with the object-oriented interface

Complementarity constraints can be specified in the object-oriented interface by defining the constraints in a class inheriting from KTRIProblem.

The KTRIProblem should implement the functions:

std::vector<int> complementarityIndexList1();
std::vector<int> complementarityIndexList2();

to return the lists of complementary variables. Parameter indexList1 and indexList2, of the same length, specifying the variable indices for the first and second sets of variables in the pairs of complementary variables.

When using the KTRProblem class, the values can be passed to the function:

KTRIProblem::setComplementarity(const std::vector<int>& indexList1,
                          const std::vector<int>& indexList2)

to set the values returned by the complementarityIndexList functions.

Note

Variables which are specified as complementary through KTRIProblem::setComplementarity() functions should have a lower bound of 0. This can be set using KTRProblem::setVarLoBnds().

AMPL example

The AMPL model for our toy problem above is the following.

# Variables
var x{j in 0..7} >= 0;

# Objective function
minimize obj:
                (x[0]-5)^2 + (2*x[1]+1)^2;

# Constraints
s.t. c0: 2*(x[1]-1) - 1.5*x[0] + x[2] - 0.5*x[3] + x[4] = 0;
s.t. c1: 3*x[0] - x[1] - 3 - x[5] = 0;
s.t. c2: -x[0] + 0.5*x[1] + 4 - x[6] = 0;
s.t. c3: -x[0] - x[1] + 7 - x[7] = 0;
s.t. c4: 0 <= x[5] complements x[2] >= 0;
s.t. c5: 0 <= x[6] complements x[3] >= 0;
s.t. c6: 0 <= x[7] complements x[4] >= 0;

Running it through AMPL, we get the following output.

=======================================
          Commercial License
         Artelys Knitro 10.0.0
=======================================

No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables and 0 constraints.

datacheck:            0
hessian_no_f:         1
par_concurrent_evals: 0
The problem is identified as an MPEC.
Knitro changing algorithm from AUTO to 1.
Knitro changing bar_initpt from AUTO to 3.
Knitro changing bar_murule from AUTO to 4.
Knitro changing bar_penaltycons from AUTO to 1.
Knitro changing bar_penaltyrule from AUTO to 2.
Knitro changing bar_switchrule from AUTO to 1.
Knitro changing linsolver from AUTO to 2.

Problem Characteristics                    ( Presolved)
-----------------------
Objective goal:  Minimize
Number of variables:                    11 (        11)
    bounded below:                      11 (        11)
    bounded above:                       0 (         0)
    bounded below and above:             0 (         0)
    fixed:                               0 (         0)
    free:                                0 (         0)
Number of constraints:                   7 (         7)
    linear equalities:                   7 (         7)
    nonlinear equalities:                0 (         0)
    linear inequalities:                 0 (         0)
    nonlinear inequalities:              0 (         0)
    range:                               0 (         0)
Number of complementarities:             3 (         3)
Number of nonzeros in Jacobian:         20 (        20)
Number of nonzeros in Hessian:           2 (         2)

  Iter      Objective      FeasError   OptError    ||Step||    CGits
--------  --------------  ----------  ----------  ----------  -------
       0    2.811162e+01   1.548e+00
      10    1.700009e+01   2.169e-07   6.513e-05   1.068e-03        0
      11    1.700000e+01   1.413e-11   4.264e-09   1.774e-04        0

EXIT: Locally optimal solution found.

Final Statistics
----------------
Final objective value               =   1.70000000056682e+01
Final feasibility error (abs / rel) =   1.41e-11 / 9.13e-12
Final optimality error  (abs / rel) =   4.26e-09 / 5.33e-10
# of iterations                     =         11
# of CG iterations                  =          1
# of function evaluations           =         16
# of gradient evaluations           =         13
# of Hessian evaluations            =         11
Total program time (secs)           =       0.00246 (     0.002 CPU time)
Time spent in evaluations (secs)    =       0.00005

===============================================================================

Knitro 10.0.0: Locally optimal or satisfactory solution.
objective 17.000000005668205; feasibility error 1.41e-11
11 iterations; 16 function evaluations

Knitro received our three complementarity constraints correctly (“Number of complementarities: 3”) and converged successfully (“Locally optimal solution found”).

MATLAB example

The following functions can be used in MATLAB to solve the same example as is shown for AMPL.

function exampleMPEC1

Jpattern = [];

Hpattern = sparse(zeros(8));
Hpattern(1,1) = 1;
Hpattern(2,2) = 1;

options = optimset('Algorithm', 'interior-point', 'Display','iter', ...
 'GradObj','on','GradConstr','on', ...
 'JacobPattern',Jpattern,'Hessian','user-supplied','HessPattern',Hpattern, ...
 'HessFcn',@hessfun,'MaxIter',1000, ...
 'TolX', 1e-15, 'TolFun', 1e-8, 'TolCon', 1e-8);

A = []; b = [];
Aeq = [-1.5  2   1 -0.5 1  0  0  0;
        3   -1   0  0   0 -1  0  0;
       -1    0.5 0  0   0  0 -1  0;
       -1   -1   0  0   0  0  0 -1];
beq = [2 3 -4 -7];
lb = zeros(8,1);
ub = Inf*ones(8,1);
x0 = zeros(8,1);

extendedFeatures.ccIndexList1 = [6 7 8];
extendedFeatures.ccIndexList2 = [3 4 5];

[x,fval,exitflag,output,lambda] = ...
    knitromatlab(@objfun,x0,A,b,Aeq,beq,lb,ub,@constfun,extendedFeatures,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [f,g] = objfun(x)

f = (x(1)-5)^2 + (2*x(2)+1)^2;

if nargout > 1
 g = zeros(8,1);
 g(1) = 2*(x(1)-5);
 g(2) = 4*(2*x(2)+1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c,ceq,Gc,Gceq]= constfun(x)

c = [];
ceq=[];
Gc = [];
Gceq=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [H]= hessfun(x,lambda)

H=sparse(zeros(8));

H(1,1) = 2;
H(2,2) = 4;

Running this file will produce the following output from Knitro.

=======================================
          Commercial License
         Artelys Knitro 10.0.0
=======================================

Knitro presolve eliminated 0 variables and 0 constraints.

algorithm:            1
feastol:              1e-08
honorbnds:            1
maxit:                1000
opttol:               1e-08
outlev:               4
par_concurrent_evals: 0
The problem is identified as an MPEC.
Knitro changing bar_initpt from AUTO to 3.
Knitro changing bar_murule from AUTO to 4.
Knitro changing bar_penaltycons from AUTO to 1.
Knitro changing bar_penaltyrule from AUTO to 2.
Knitro changing bar_switchrule from AUTO to 1.
Knitro changing linsolver from AUTO to 2.
Knitro shifted start point to satisfy presolved bounds (8 variables).

Problem Characteristics                    ( Presolved)
-----------------------
Objective goal:  Minimize
Number of variables:                     8 (         8)
    bounded below:                       8 (         8)
    bounded above:                       0 (         0)
    bounded below and above:             0 (         0)
    fixed:                               0 (         0)
    free:                                0 (         0)
Number of constraints:                   4 (         4)
    linear equalities:                   4 (         4)
    nonlinear equalities:                0 (         0)
    linear inequalities:                 0 (         0)
    nonlinear inequalities:              0 (         0)
    range:                               0 (         0)
Number of complementarities:             3 (         3)
Number of nonzeros in Jacobian:         14 (        14)
Number of nonzeros in Hessian:           2 (         2)

  Iter     fCount     Objective      FeasError   OptError    ||Step||    CGits
--------  --------  --------------  ----------  ----------  ----------  -------
       0         2    2.496050e+01   4.030e+00
       1         3    2.847389e+01   1.748e+00   2.160e+00   1.990e+00        1
       2         4    4.226663e+01   3.832e-01   4.643e+00   1.442e+00        0
       3         5    4.667799e+01   1.126e-02   3.638e+00   5.993e-01        0
       4         6    4.213217e+01   4.179e-03   1.258e+01   1.185e+00        0
       5         7    4.074018e+01   3.072e-03   1.265e+01   1.580e-01        1
       6         8    3.810894e+01   1.133e-04   1.259e+01   3.113e-01        0
       7         9    1.701407e+01   1.682e-04   1.542e+00   4.771e+00        0
       8        10    1.699966e+01   1.522e-04   6.416e-02   2.385e-02        0
       9        11    1.700003e+01   1.799e-06   3.532e-05   2.154e-04        0
      10        12    1.700000e+01   6.354e-11   1.530e-09   5.298e-05        0

EXIT: Locally optimal solution found.

Final Statistics
----------------
Final objective value               =   1.70000000010379e+01
Final feasibility error (abs / rel) =   6.35e-11 / 1.58e-11
Final optimality error  (abs / rel) =   1.53e-09 / 1.91e-10
# of iterations                     =         10
# of CG iterations                  =          2
# of function evaluations           =         12
# of gradient evaluations           =         12
# of Hessian evaluations            =         10
Total program time (secs)           =       0.00827 (     0.019 CPU time)
Time spent in evaluations (secs)    =       0.00464

===============================================================================

C example

The same example can be implemented using the callable library. Arrays indexList1 and indexList2 are used to specify the list of complementarities and the KTR_addcompcons() function is called to register the list.

#include <stdio.h>
#include <stdlib.h>
#include "knitro.h"


/* callback function that evaluates the objective
   and constraints */
int  callback (const int             evalRequestCode,
                     const int             n,
                     const int             m,
                     const int             nnzJ,
                     const int             nnzH,
                     const double * const  x,
                     const double * const  lambda,
                           double * const  obj,
                           double * const  c,
                           double * const  objGrad,
                           double * const  jac,
                           double * const  hessian,
                           double * const  hessVector,
                           void   *        userParams)
{
    if (evalRequestCode == KTR_RC_EVALFC) {
        /* objective function */
        *obj = (x[0]-5)*(x[0]-5) + (2*x[1]+1)*(2*x[1]+1);

        /* constraints */
        c[0] = 2*(x[1]-1) - 1.5*x[0] + x[2] - 0.5*x[3] + x[4];
        c[1] = 3*x[0] - x[1] - 3 -x[5];
        c[2] = -x[0] + 0.5*x[1] + 4 -x[6];
        c[3] = -x[0] - x[1] + 7 - x[7];

        return(0);
    } else {
        printf ("Wrong evalRequestCode in callback function.\n");
        return(-1);
    }
}

/* main */
int  main (int  argc, char  *argv[]) {
    int  nStatus;

    /* variables that are passed to Knitro */
    KTR_context    *kc;
    int            n, m, numCompConstraints, nnzJ, nnzH, objGoal, objType;
    int            *cType, *indexList1, *indexList2;
    int            *jacIndexVars, *jacIndexCons;
    double        obj, *x, *lambda;
    double        *xLoBnds, *xUpBnds, *xInitial, *cLoBnds, *cUpBnds;
    int            i, j, k; // convenience variables

    /*problem size and mem allocation */
    n = 8;                  /* number of variables */
    m = 4;                  /* number of regular constraints */
    numCompConstraints = 3; /* number of complementarity constraints */
    nnzJ = n*m;
    nnzH = 0;
    x      = (double *) malloc (n     * sizeof(double));
    lambda = (double *) malloc ((m+n) * sizeof(double));
    xLoBnds      = (double *) malloc (n * sizeof(double));
    xUpBnds      = (double *) malloc (n * sizeof(double));
    xInitial     = (double *) malloc (n * sizeof(double));
    cType        = (int    *) malloc (m * sizeof(int));
    cLoBnds      = (double *) malloc (m * sizeof(double));
    cUpBnds      = (double *) malloc (m * sizeof(double));
    jacIndexVars = (int    *) malloc (nnzJ * sizeof(int));
    jacIndexCons = (int    *) malloc (nnzJ * sizeof(int));
    indexList1   = (int    *) malloc (numCompConstraints * sizeof(int));
    indexList2   = (int    *) malloc (numCompConstraints * sizeof(int));

    /* objective type */
    objType = KTR_OBJTYPE_GENERAL;
    objGoal = KTR_OBJGOAL_MINIMIZE;

    /* bounds and constraints type */
    for (i = 0; i < n; i++) {
        xLoBnds[i] = 0.0;
        xUpBnds[i] = KTR_INFBOUND;
    }
    for (j = 0; j < m; j++) {
        cType[j] = KTR_CONTYPE_GENERAL;
        cLoBnds[j] = 0.0;
        cUpBnds[j] = 0.0;
    }

    /* complementarities */
    indexList1[0] = 2;   indexList2[0] = 5;
    indexList1[1] = 3;   indexList2[1] = 6;
    indexList1[2] = 4;   indexList2[2] = 7;

    /* sparsity pattern (here, of a full matrix) */
    k = 0;
    for (i = 0; i < n; i++)
        for (j = 0; j < m; j++) {
            jacIndexCons[k] = j;
            jacIndexVars[k] = i;
            k++;
        }

    /* create a Knitro instance */
    kc = KTR_new();
    if (kc == NULL)
        exit( -1 ); // probably a license issue

    /* set options: automatic gradient and hessian matrix */
    if (KTR_set_int_param_by_name (kc, "gradopt", 3) != 0)
        exit( -1 );
    if (KTR_set_int_param_by_name (kc, "hessopt", 6) != 0)
        exit( -1 );

    /* register the callback function */
    if (KTR_set_func_callback (kc, &callback) != 0)
        exit( -1 );

    /* pass the problem definition to Knitro */
    nStatus = KTR_init_problem (kc, n, objGoal, objType,
                                xLoBnds, xUpBnds,
                                m, cType, cLoBnds, cUpBnds,
                                nnzJ, jacIndexVars, jacIndexCons,
                                nnzH, NULL, NULL, xInitial, NULL);

    /* declare complementarities */
    KTR_addcompcons (kc, numCompConstraints, indexList1, indexList2);

    /* free memory (Knitro maintains its own copy) */
    free (xLoBnds);
    free (xUpBnds);
    free (xInitial);
    free (cType);
    free (cLoBnds);
    free (cUpBnds);
    free (jacIndexVars);
    free (jacIndexCons);
    free (indexList1);
    free (indexList2);

    /* solver call */
    nStatus = KTR_solve (kc, x, lambda, 0, &obj,
        NULL, NULL, NULL, NULL, NULL, NULL);

    if (nStatus != 0)
        printf ("\nKnitro failed to solve the problem, final status = %d\n",
                nStatus);
    else
        printf ("\nKnitro successful, objective is = %e\n", obj);

    /* delete the Knitro instance and primal/dual solution */
    KTR_free (&kc);
    free (x);
    free (lambda);

    getchar();
    return( 0 );
}

Running this code produces an output similar to what we obtained with AMPL for the same problem. More function evaluations are made since, for simplicity, we did not provide first derivatives and failed to notify Knitro that the constraints are linear. The final objective value is however the same:

Knitro successful, objective is = 1.700000e+001