GoalObj - Archimedian and pre-emptive goal programming using objective functions

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

  file xbgoalobj.c
  ````````````````
  Archimedian and pre-emptive goal programming
  using objective functions.

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

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

#define NGOALS 3

/***** Data *****/
const char *Type[] = {"perc", "abs", "perc"};
const char *Sense[] = {"max", "min", "max"};
double Weight[] = {100, 1, 0.1};
double Deviation[] = {10, 4, 20};
double Coeff[NGOALS][3] = {{5, 2, 20}, {-3, 15, 48}, {1.5, 21, 3.8}};

int main(int argc, char **argv)
{
 XPRBvar x[2];
 XPRBctr goalCtr[NGOALS], aCtr, wobj;
 double Target[NGOALS];
 XPRBprob prob;
 int i,g;

 if(XPRBinit()) return 1;

 prob = XPRBnewprob("Goal");

 /* Adding the variables */
 x[0] = XPRBnewvar(prob,XPRB_PL,"x",0,XPRB_INFINITY);
 x[1] = XPRBnewvar(prob,XPRB_PL,"y",0,XPRB_INFINITY);

 /* Adding a constraint */
 aCtr = XPRBnewctr(prob,"Limit",XPRB_L);
 XPRBaddterm(aCtr, x[0], 42);
 XPRBaddterm(aCtr, x[1], 13);
 XPRBaddterm(aCtr, NULL, 100);

 /* Goals */
/* goal[0] = 5*x +  2*y - 20;
   goal[1] = -3*x + 15*y - 48;
   goal[2] = 1.5*x + 21*y - 3.8; */
 for(g=0;g<NGOALS;g++)
 {
  goalCtr[g] = XPRBnewctr(prob,"goal",XPRB_N);
  for(i=0;i<2;i++) XPRBaddterm(goalCtr[g], x[i], Coeff[g][i]);
  XPRBaddterm(goalCtr[g], NULL, Coeff[g][2]);
 }

 /***** Archimedian GP *****/
 printf("Archimedian:");
 wobj = XPRBnewctr(prob,"wobj",XPRB_N);
 for(g=0;g<NGOALS;g++)
 {
  if (strcmp(Sense[g],"max")==0)
  {
   for(i=0;i<2;i++) XPRBaddterm(wobj, x[i], - Weight[g]*Coeff[g][i]);
   XPRBaddterm(wobj, NULL, - Weight[g]*Coeff[g][i]);
  }
  else
  {
   for(i=0;i<2;i++) XPRBaddterm(wobj, x[i], Weight[g]*Coeff[g][i]);
   XPRBaddterm(wobj, NULL, Weight[g]*Coeff[g][i]);
  }
 }
 XPRBsetobj(prob,wobj);
 XPRSsetintcontrol(XPRBgetXPRSprob(prob),XPRS_OUTPUTLOG, 0);
 XPRBlpoptimize(prob,"");

 /* Solution printout */
 printf(" Solution: x: %g, y: %g\n", XPRBgetsol(x[0]), XPRBgetsol(x[1]));
 printf(" Goal   Target     Value\n");
 for(g=0;g<NGOALS;g++)
   printf("  %d       %s     %g\n", g+1, Sense[g], XPRBgetact(goalCtr[g]) - XPRBgetrhs(goalCtr[g]));


 /***** Pre-emptive GP *****/
 printf("Pre-emptive:\n");
 i=-1;
 while (i<NGOALS-1)
 {
  i+=1;
  if (strcmp(Sense[i],"max")==0)
  {
   XPRBsetobj(prob,goalCtr[i]);
   XPRBsetsense(prob,XPRB_MAXIM);
   XPRBlpoptimize(prob,"");
   if (XPRBgetlpstat(prob) != XPRB_LP_OPTIMAL)
   {
    printf("Cannot satisfy goal %d\n", i+1);
    break;
   }
   else
   {
    Target[i]=XPRBgetobjval(prob);
    if (strcmp(Type[i],"perc")==0)
     Target[i]-= fabs(Target[i])*Deviation[i]/100;
    else
     Target[i]-= Deviation[i];
    if (i<NGOALS-1) XPRBaddterm(goalCtr[i], NULL, Target[i]);
    XPRBsetctrtype(goalCtr[i],XPRB_G); printf("%g\n",Target[i]);
   }
  }
  else
  {
   XPRBsetobj(prob,goalCtr[i]);
   XPRBsetsense(prob,XPRB_MINIM);
   XPRBlpoptimize(prob,"");
   if (XPRBgetlpstat(prob) != XPRB_LP_OPTIMAL)
   {
    printf("Cannot satisfy goal %d\n", i+1);
    break;
   }
   else
   {
    Target[i]=XPRBgetobjval(prob);
    if (strcmp(Type[i],"perc")==0)
     Target[i]+= fabs(Target[i])*Deviation[i]/100;
    else
     Target[i]+= Deviation[i];
    if (i<NGOALS-1) XPRBaddterm(goalCtr[i], NULL, Target[i]);
    XPRBsetctrtype(goalCtr[i],XPRB_L); printf("%g\n",Target[i]);
   }
  }
  printf("Solution(%d):  x: %g, y: %g\n", i+1, XPRBgetsol(x[0]), XPRBgetsol(x[1]));
 }

 /* Solution printout */
 printf(" Goal       Target         Value\n");
  for(g=0;g<=i;g++)
  {
   printf("  %d     %s%g", g+1, (XPRBgetctrtype(goalCtr[g])==XPRB_G?" >=  ":" <=  "), Target[g]);
   if(g==NGOALS-1)
    printf("     %g\n", XPRBgetobjval(prob));
   else
    printf("     %g\n", (XPRBgetact(goalCtr[g]) - XPRBgetrhs(goalCtr[g])) + Target[g]);
  }

 XPRBfinish();

 return 0;
}