Portfolio - Quadratic Programming with discrete variables

/********************************************************
  BCL Example Problems
  ====================

  file xbportf.c
  ``````````````
  Quadratic portfolio model.

  (c) 2008 Fair Isaac Corporation
      author: S.Heipcke, Jan. 2000, rev. Mar. 2011
********************************************************/

/* In this model, a choice has to be made which values are taken   *
 * into a portfolio in order to minimize the total cost. The costs *
 * for some values are interrelated, introducing a quadratic part  *
 * to the objective function. Upper bounds are given on the total  *
 * number of values and the share of each value that may be taken. */

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

#define NVal 30                       /* Total number of values */
#define LIMIT 20                      /* Maximum number to be chosen */

#define QFILE XPRBDATAPATH "/portf/pfqcost.dat" /* Quadratic cost coeff.s */
#define BFILE XPRBDATAPATH "/portf/pfubds.dat"  /* Upper bounds on percentages */
#define CFILE XPRBDATAPATH "/portf/pflcost.dat" /* Linear cost coefficients */

/**** DATA ****/
double Cost[NVal];                    /* Coeff. of lin. part of the obj. */
double QCost[NVal][NVal];             /* Coeff. of quad. part of the obj. */
double UBnd[NVal];                    /* Upper bound values */

/***********************************************************************/

void modfolio(XPRBprob prob)
{
 int i,j;
 XPRBctr cobj, c;
 XPRBvar x[NVal];                     /* Amount of a value taken into
                                         the portfolio */
 XPRBvar y[NVal];                     /* 1 if value i is chosen, else 0 */

/**** VARIABLES ****/
 for(i=0;i<NVal;i++)
 {
  x[i] = XPRBnewvar(prob,XPRB_PL, XPRBnewname("x_%d",i+1), 0, UBnd[i]);
  y[i] = XPRBnewvar(prob,XPRB_BV, XPRBnewname("y_%d",i+1), 0, 1);
 }

/****OBJECTIVE****/
 cobj = XPRBnewctr(prob,"OBJ",XPRB_N);   /* Define objective: total cost */
 for(i=0;i<NVal;i++)
  XPRBaddterm(cobj,x[i],Cost[i]);
 for(i=0;i<NVal;i++)                  /* Define quadratic part of the obj. */
 {
  XPRBaddqterm(cobj,x[i],x[i],QCost[i][i]);
  for(j=i+1;j<NVal;j++)
   XPRBaddqterm(cobj,x[i],x[j], QCost[i][j]);
 }
 XPRBsetobj(prob,cobj);               /* Set objective function */

/**** CONSTRAINTS ****/
 c = XPRBnewctr(prob,"C1",XPRB_E);    /* Amounts of values chosen must
                                         add up to 100% */
 for(i=0;i<NVal;i++) XPRBaddterm(c, x[i], 1);
 XPRBaddterm(c, NULL, 100);

 for(i=0;i<NVal;i++)                  /* Upper limits */
 {
  c = XPRBnewctr(prob,"UL",XPRB_L);
  XPRBaddterm(c, x[i], 1);
  XPRBaddterm(c, y[i], -UBnd[i]);
 }

 c = XPRBnewctr(prob,"Card", XPRB_L);    /* Limit on total number of values */
 for(i=0;i<NVal;i++) XPRBaddterm(c, y[i], 1);
 XPRBaddterm(c, NULL, LIMIT);

/****SOLVING + OUTPUT****/
/* XPRBprintprob(prob); */               /* Print out the problem definition */
 XPRBexportprob(prob,XPRB_MPS,"Portf");  /* Output matrix in MPS format */
 XPRBexportprob(prob,XPRB_LP,"Portf");   /* Output matrix in LP format */

 XPRBsetsense(prob,XPRB_MINIM);          /* Choose the sense of the optimization */
 XPRSsetintcontrol(XPRBgetXPRSprob(prob),XPRS_CUTSTRATEGY,0);
 XPRBmipoptimize(prob,"");               /* Solve the QP-problem, use
                                            'mipoptimize' to solve MIQP,
                                            'lpoptimize' for QP */

 for(i=0;i<NVal;i++)
  printf("%s:%g (%g), ", XPRBgetvarname(x[i]), XPRBgetsol(x[i]),
    XPRBgetsol(y[i]));
 printf("\n");
}

/***********************************************************************/

    /**** Read data from files ****/
void readdata(void)
{
 int i,j;
 double value;
 FILE *datafile;

        /* Read the quadratic cost data file */
 memset(QCost,0,NVal*NVal*sizeof(double));      /* Initialize Q to 0 */
 datafile=fopen(QFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,i,g", &i, &j,
        &value) == 3)
  QCost[i-1][j-1]=value;
 fclose(datafile);

        /* Read the linear cost data file */
 datafile=fopen(CFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,g", &i, &value) == 2)
  Cost[i-1]=value;
 fclose(datafile);

        /* Read the bounds data file */
 datafile=fopen(BFILE,"r");
 while (XPRBreadlinecb(XPRB_FGETS, datafile, 200, "i,g", &i, &value) == 2)
  UBnd[i-1]=value;
 fclose(datafile);
}

/***********************************************************************/

int main(int argc, char **argv)
{
 XPRBprob prob;

 prob=XPRBnewprob("Portfolio");       /* Initialize a new problem in BCL */
 readdata();                          /* Data input from file */
 modfolio(prob);                      /* Formulate and solve the problem */

 return 0;
}