Quadratic Programming with BCL

As an extension to LP and MIP, BCL also provides support for formulating and solving Quadratic Programming (QP) and Mixed Integer Quadratic Programming (MIQP) problems, that is, problems with linear constraints with a quadratic objective function of the form

where x is the vector of decision variables, c is the cost vector, and Q is the quadratic cost coefficient matrix. The matrix Q must be symmetric. It should also be positive semi-definite if the problem is to be minimized, and negative semi-definite if it is to be maximized, because the Xpress Optimizer solves convex QP problems. If the problem is not convex, the solution algorithms may not converge at all, or may only converge to a locally optimal solution.

Release 4.0 of BCL extends this functionality to Quadratically Constrained Quadratic Programming (QCQP) problems, that is, problems that in addition to a quadratic objevctive function have constraints of the form

where a is the coefficient vector for the linear terms, b the constant RHS value, and the same conditions as in objective functions apply to the quadratic coefficient matrix Q (positive semi-definite in ≤ constraints, and negative semi-definite in ≥ constraints). Quadratic constraints in QCQP problems must be inequalities.

Any other quadratic form supported by the Xpress Optimizer can also be used, e.g. Second Order Cone constraints in Second Order Cone problems (SOCPs).

Add quadratic term

XPRBctr c;
XPRBvar x1;
XPRBaddqterm(c,x1,x1,3);

Set quadratic term

XPRBvar x2;
XPRBsetqterm(c,x1,x2,-7.2);

Delete a quadratic term

XPRBdelqterm(c,x2,x1);

Enumerate quadratic terms

double coeff;
const void *ref;
ref = XPRBgetnextqterm(c,ref,&x1,&x1,&coeff);

Figure 3.4: Defining and accessing quadratic terms in BCL.

In BCL, the quadratic part of constraints is defined termwise, much like what we have seen for the definition of linear constraints in Section Constraints. The coefficient of a quadratic term is either set to a given value (XPRBsetqterm) or its value is augmented by the given value (XPRBaddqterm). Quadratic objective functions are set in the same way as linear ones with a call to XPRBsetobj. Note that the definition of the quadratic constraint terms should always be preceded by the definition of the corresponding variables.

The coefficient of a quadratic constraint term can be retrieved with XPRBgetqcoeff. The quadratic terms of a constraint can be enumerated with XPRBgetnextqterm.

Unless BCL is used in Student Mode, functions XPRBprintprob, XPRBprintobj, XPRBexportprob, and XPRBprintctr will print or output to a file the complete problem / constraint definition, including the quadratic terms.

Example

We wish to distribute a set of points represented by tuples of x-/y-coordinates on a plane minimizing the total squared distance between all pairs of points. For each point i we are given a target location (CX_i,CY_i) and the (square of the) maximum allowable distance R_i to this location.

In mathematical terms, we have two decision variables x_i and y_i for the coordinates of every point i. The objective to minimize the total squared distance between all points is expressed by the following sum.

For every point i we have the following quadratic inequality.

The following BCL program (xbairport.c) implements and solves this problem.

#include <stdio.h>
#include "xprb.h"

#define N 42
double CX[N], CY[N], R[N];

...                                   /* Initialize the data arrays */

int main(int argc, char **argv)
{
 int i,j;
 XPRBprob prob;
 XPRBvar x[N],y[N];                /* x-/y-coordinates to determine */
 XPRBctr cobj, c;

 prob=XPRBnewprob("airport");    /* Initialize a new problem in BCL */

/**** VARIABLES ****/
 for(i=0;i<N;i++)
  x[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("x(%d)",i), -10, 10);
 for(i=0;i<N;i++)
  y[i] = XPRBnewvar(prob, XPRB_PL, XPRBnewname("y(%d)",i), -10, 10);

/****OBJECTIVE****/
/* Minimize the total distance between all points */
 cobj = XPRBnewctr(prob, "TotDist", XPRB_N);
  for(i=0;i<N-1;i++)
   for(j=i+1;j<N;j++)
   {
    XPRBaddqterm(cobj, x[i], x[i], 1);
    XPRBaddqterm(cobj, x[i], x[j], -2);
    XPRBaddqterm(cobj, x[j], x[j], 1);
    XPRBaddqterm(cobj, y[i], y[i], 1);
    XPRBaddqterm(cobj, y[i], y[j], -2);
    XPRBaddqterm(cobj, y[j], y[j], 1);
   }
 XPRBsetobj(prob, cobj);             /* Set the objective function */

/**** CONSTRAINTS ****/
/* All points within given distance of their target location */
 for(i=0;i<N;i++)
 {
  c = XPRBnewctr(prob, XPRBnewname("LimDist_%d",i), XPRB_L);
  XPRBaddqterm(c, x[i], x[i], 1);
  XPRBaddterm(c, x[i], -2*CX[i]);
  XPRBaddterm(c, NULL, -CX[i]*CX[i]);
  XPRBaddqterm(c, y[i], y[i], 1);
  XPRBaddterm(c, y[i], -2*CY[i]);
  XPRBaddterm(c, NULL, -CY[i]*CY[i]);
  XPRBaddterm(c, NULL, R[i]);
 }

/****SOLVING + OUTPUT****/
 XPRBsetsense(prob, XPRB_MINIM);          /* Sense of optimization */
 XPRBlpoptimize(prob,"");                     /* Solve the problem */

 printf("Solution: %g\n", XPRBgetobjval(prob));
 for(i=0;i<N;i++)
  printf(" %d: %g, %g\n", i, XPRBgetsol(x[i]), XPRBgetsol(y[i]));

 return 0;
}