/* Name:       pooling-slp.c
** Author:     Leo Liberti
** Purpose:    solving Haverly's pooling problem with SLP
** Source:     C
** History:    060220 work started
*/

// standard include files
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <ilcplex/cplex.h>
#include <assert.h>

// maximum length of variable/constraint names
#define MAXNAME 64
#define INFINITY 1e30
#define EPSILON 1e-6
#define MAXITN 10
#define TRUE 1
#define FALSE 0

/* Solve Haverly's pooling problem by Sequential Linear Programming with CPLEX

   Haverly's pooling problem (bilinear) formulation:
   min  6 x11 + 10 x12 +  16 x21   - 9 y11 - 15 y12 -   9 y21 -  15 y22
   s.t.   x11 +              x21   -   y11 -    y12                      = 0
                   x12                              -     y21 -     y22  = 0
                                       y11          +     y21           <= 100
                                                y12           +     y22 <= 200
        3 x11              + x21   - p y11  - p y12                      = 0
                               (p-2.5) y11          - 0.5 y21           <= 0
                                        (p-1.5) y12           + 0.5 y22 <= 0
        x, y, p >= 0


   Decompose into two linear problems (alternatively fixing p, y):

   - LP1 (fix p)

       x11 x12 x21 y11     y12     y21  y22
     obj fun coeffs:                                       sense      rhs
        6  10  16  -9     -15      -9   -15       
     constraint matrix:
        1   0   1  -1      -1       0     0                 =          0
        0   1   0   0       0      -1    -1                 =          0
        0   0   0   1       0       1     0                <=        100
        0   0   0   0       1       0     1                <=        200
        3   0   1  -p      -p       0     0                 =          0
        0   0   0  (p-2.5)  0      -0.5   0                <=          0
        0   0   0   0      (p-1.5)  0     0.5              <=          0

   - LP2 (fix y11, y12):

       x11 x12 x21                 y21  y22   p
     obj fun coeffs:                                      sense     rhs
        6  10  16                  -9   -15   0                -9 y11 - 15 y12
     constraint matrix:
        1   0   1                   0     0   0             =     y11 + y12
        0   1   0                  -1    -1   0             =         0
        0   0   0                   1     0   0            <=     100 - y11
        0   0   0                   0     1   0            <=     200 - y12
        3   0   1                             -(y11+y12)    =     0
        0   0   0                  -0.5   0   y11          <=     2.5 y11
        0   0   0                   0     0.5 y12          <=     1.5 y12

   Haverly's pooling problem global optimum is at f* = -400. 

*/

// test if x, xsave are equal (w.r.t. LP1 and LP2, respectively)
int is_equal(double* x, double* xsave, int thetype) {
  if (thetype == 1) {
    // lp1
    if (fabs(x[0] - xsave[0]) < EPSILON && 
	fabs(x[1] - xsave[1]) < EPSILON && 
	fabs(x[2] - xsave[2]) < EPSILON && 
	fabs(x[3] - xsave[3]) < EPSILON && 
	fabs(x[4] - xsave[4]) < EPSILON && 
	fabs(x[5] - xsave[5]) < EPSILON && 
	fabs(x[6] - xsave[6]) < EPSILON) {
      return TRUE;
    } else {
      return FALSE;
    }
  } else if (thetype == 2) {
    // lp2
    if (fabs(x[0] - xsave[0]) < EPSILON && 
	fabs(x[1] - xsave[1]) < EPSILON && 
	fabs(x[2] - xsave[2]) < EPSILON && 
	fabs(x[5] - xsave[5]) < EPSILON && 
	fabs(x[6] - xsave[6]) < EPSILON && 
	fabs(x[7] - xsave[7]) < EPSILON) {
      return TRUE;
    } else {
      return FALSE;
    }
  }
}

void AddConstraint(CPXENVptr env, CPXLPptr lp, int nz, 
		   double* coeffs, int* ind, char sense, double rhs) {
  int status = 0;
  // add constraints in row-wise sparse format
  int rmatbeg[2];
  // nonzeroes stored in rmatval starting at position:
  rmatbeg[0] = 0;
  // nonzeroes stored in rmatval ending at position:
  rmatbeg[1] = nz-1;
  // create the constraints
  status = CPXaddrows(env, lp, 0, 1, nz, &rhs, &sense, 
		      rmatbeg, ind, coeffs, NULL, NULL);
  if (status != 0) {
    fprintf(stderr, "pooling-slp: could not create constraints, error %d\n", 
	    status);
    exit(3);
  }
}

int main(int argc, char** argv) {

  // return code
  int ret = 0;
  // status code returned by cplex callable library functions
  int status = 0;
  // number of variables in the problem
  int n = 8;
  // number of constraints in the problem
  int m = 7;
  // cplex environment to pass to/from cplex callable library
  CPXENVptr env = NULL;
  // linear program 1
  CPXLPptr  lp1 = NULL;
  // linear program 2
  CPXLPptr  lp2 = NULL;
  // use for storing error messages returned by CPLEX
  char errmsg[1024];
  // counters
  int i, j;
  // used whilst adding constraints
  double thecoeffs[8];
  int theind[8];
  char thesense;
  double therhs;
  int thenz = 0;
  
  // try to initialize the CPLEX environment 
  env = CPXopenCPLEX(&status);
  if (env == NULL) {
    fprintf(stderr, "%s: could not open CPLEX environment\n", argv[0]);
    CPXgeterrorstring (env, status, errmsg);
    fprintf (stderr, "%s", errmsg);
    exit(2);
  }

  // turn on output to screen
  status = CPXsetintparam(env, CPX_PARAM_SCRIND, CPX_OFF);

  // bilinear variables values
  double y[2] = { 100, 100 };
  double p = 2;

  // first LP
  // create the first LP; beware: x[7] (p) is not present in LP1
  lp1 = CPXcreateprob(env, &status, "lp1");
  if (lp1 == NULL) {
    fprintf(stderr, "%s: failed to create problem lp1\n");
    exit(4);
  }

  // create variables in lp1
  int n1 = 7;
  int m1 = 7;
  // objective function coefficients
  double c1[] = { 6, 10, 16, -9, -15, -9, -15 };
  // variable lower bounds are important!
  double xL1[7];
  // x11
  xL1[0] = 0;
  // x12
  xL1[1] = 0;
  // x21
  xL1[2] = 0;
  // y11, y12: force to EPSILON to start off, or consider constr y11+y12>=eps
#ifdef ADDITIONALCONSTRAINT  
  xL1[3] = 0; 
  xL1[4] = 0; 
#else
  xL1[3] = EPSILON; 
  xL1[4] = EPSILON; 
#endif
  // y21
  xL1[5] = 0;
  // y22
  xL1[6] = 0;
  // variable upper bounds
  double xU1[7];
  for(j = 0; j < n1; j++) {
    xU1[j] = CPX_INFBOUND;
  }

  // create the variables
  status = CPXnewcols(env, lp1, n1, c1, xL1, xU1, NULL, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: could not create variables in lp1, error %d\n", 
	    argv[0], status);
    exit(3);
  }

  // add constraints in lp1
  // constraint 1: x11 + x21 - y11 - y12 = 0
  thenz = 4;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = 1.0;
  thecoeffs[2] = -1.0;
  thecoeffs[3] = -1.0;
  theind[0] = 0;
  theind[1] = 2;
  theind[2] = 3;
  theind[3] = 4;
  thesense = 'E';
  therhs = 0.0;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 2: x12 - y21 - y22 = 0
  thenz = 3;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = -1.0;
  thecoeffs[2] = -1.0;
  theind[0] = 1;
  theind[1] = 5;
  theind[2] = 6;
  thesense = 'E';
  therhs = 0.0;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 3: y11 + y21 <= 100
  thenz = 2;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = 1.0;
  theind[0] = 3;
  theind[1] = 5;
  thesense = 'L';
  therhs = 100;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 4: y12 + y22 <= 200
  thenz = 2;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = 1.0;
  theind[0] = 4;
  theind[1] = 6;
  thesense = 'L';
  therhs = 200;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 5: 3x11 + x21 - p y11 - p y12 = 0
  thenz = 4;
  thecoeffs[0] = 3.0;
  thecoeffs[1] = 1.0;
  thecoeffs[2] = -p;
  thecoeffs[3] = -p;
  theind[0] = 0;
  theind[1] = 2;
  theind[2] = 3;
  theind[3] = 4;
  thesense = 'E';
  therhs = 0;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 6: (p-2.5) y11 - 0.5 y21 <= 0
  thenz = 2;
  thecoeffs[0] = p-2.5;
  thecoeffs[1] = -0.5;
  theind[0] = 3;
  theind[1] = 5;
  thesense = 'L';
  therhs = 0;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 6: (p-1.5) y12 + 0.5 y22 <= 0
  thenz = 2;
  thecoeffs[0] = p-1.5;
  thecoeffs[1] = 0.5;
  theind[0] = 4;
  theind[1] = 6;
  thesense = 'L';
  therhs = 0;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
#ifdef ADDITIONALCONSTRAINT
  // constraint 7 (used to prevent premature trivial convergence): y11+y12>=eps
  thenz = 2;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = 1.0;
  theind[0] = 3;
  theind[1] = 4;
  thesense = 'G';
  therhs = EPSILON;
  AddConstraint(env, lp1, thenz, thecoeffs, theind, thesense, therhs);
#endif

  // second LP
  // create the second LP; beware: x[3], x[4] (y11, y12) are not present in LP2
  // this calls for a reindexing of x when dealing with lp2:
  // x[0]:x11, x[1]:x12, x[2]:x21, x[3]:y21, x[4]:y22, x[5]:p
  lp2 = CPXcreateprob(env, &status, "lp2");
  if (lp2 == NULL) {
    fprintf(stderr, "%s: failed to create problem lp2\n");
    exit(4);
  }

  // create variables in lp2
  int n2 = 6;
  int m2 = 7;
  // objective function coefficients
  double c2[] = { 6, 10, 16, -9, -15, 0};
  // variable lower bounds
  double xL2[6];
  // x11
  xL2[0] = xL1[0];
  // x12
  xL2[1] = xL1[1];
  // x21
  xL2[2] = xL1[2];
  // y21 (y11, y12 missing from LP2)
  xL2[3] = xL1[5];
  // y22
  xL2[4] = xL1[6];
  // p
  xL2[5] = 0;
  // variable upper bounds
  double xU2[6];
  for(j = 0; j < n2; j++) {
    xU2[j] = CPX_INFBOUND;
  }
  xU2[5] = 100; // upper bound for p

  // create the variables
  status = CPXnewcols(env, lp2, n2, c2, xL2, xU2, NULL, NULL);
  if (status != 0) {
    fprintf(stderr, "%s: could not create variables in lp2, error %d\n", 
	    argv[0], status);
    exit(3);
  }

  // add constraints in lp2
  // constraint 1: x11 + x21 = y11+y12 
  thenz = 2;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = 1.0;
  theind[0] = 0;
  theind[1] = 2;
  thesense = 'E';
  therhs = y[0] + y[1];
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 2: x12 - y21 - y22 = 0
  thenz = 3;
  thecoeffs[0] = 1.0;
  thecoeffs[1] = -1.0;
  thecoeffs[2] = -1.0;
  theind[0] = 1;
  theind[1] = 3;
  theind[2] = 4;
  thesense = 'E';
  therhs = 0.0;
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 3: y21 <= 100 - y11
  thenz = 1;
  thecoeffs[0] = 1.0;
  theind[0] = 3;
  thesense = 'L';
  therhs = 100 - y[0];
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 4: y22 <= 200 - y12
  thenz = 1;
  thecoeffs[0] = 1.0;
  theind[0] = 4;
  thesense = 'L';
  therhs = 200 - y[1];
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 5: 3x11 + x21 - (y11+y12) p = 0
  thenz = 3;
  thecoeffs[0] = 3.0;
  thecoeffs[1] = 1.0;
  thecoeffs[2] = -(y[0] + y[1]);
  theind[0] = 0;
  theind[1] = 2;
  theind[2] = 5;
  thesense = 'E';
  therhs = 0;
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 6: -0.5 y21 + y11 p <= 2.5 y11
  thenz = 2;
  thecoeffs[0] = -0.5;
  thecoeffs[1] = y[0];
  theind[0] = 3;
  theind[1] = 5;
  thesense = 'L';
  therhs = 2.5 * y[0];
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);
  // constraint 6:  0.5 y22 + y12 p <= 1.5 y12
  thenz = 2;
  thecoeffs[0] = 0.5;
  thecoeffs[1] = y[1];
  theind[0] = 4;
  theind[1] = 5;
  thesense = 'L';
  therhs = 1.5 * y[1];
  AddConstraint(env, lp2, thenz, thecoeffs, theind, thesense, therhs);

  // select the optimization algorithm
  status = CPXsetintparam(env, CPX_PARAM_LPMETHOD, CPX_ALG_AUTOMATIC);
  if (status != 0) {
    fprintf(stderr, "%s: could not select optimization algorithm, error %d\n", 
	    argv[0], status);
    exit(3);
  }

  // current objective value
  double objval = 0;
  // best objective value
  double bestobj = 1e30;
  // current variable values
  double x[8] = { 1111, 1111, 1111, 1111, 1111, 1111, 1111, 1111 };
  // best variable values
  double xast[8];
  // saved variable values
  double xsave[8];
  // dual variables
  double lambda[7];
  // slacks
  double s[7];
  // reduced costs
  double rc[8];

  // do it!
  fprintf(stderr, "************ SLP: Algorithmic progress ****************\n");
  int termination = FALSE;
  int compareflag1, compareflag2;
  int objcompare1, objcompare2;
  int objcompcount = 0;
  int noprogresscount = 0;
  double objvalsave;
  int itncount = 0;
  fprintf(stderr, "\tOF\tx11\tx12\tx21\ty11\ty12\ty21\ty22\tp\n");
  while(!termination) {

    itncount++;
    fprintf(stderr, "iteration %d\n", itncount);
    compareflag1 = FALSE;
    compareflag2 = FALSE;
    
    // optimize lp1
    status = CPXlpopt(env, lp1);
    if (status != 0) {
      fprintf(stderr, "%s: failed to call optimization algorithm on lp1\n", 
	      argv[0]);
      exit(10);
    }
    objvalsave = objval;
    // solve lp1, this sets x[0]...x[6] to correct values
    CPXsolution(env, lp1, &status, &objval, x, lambda, s, rc);
    // update x[7] with value of variable missing from lp1
    x[7] = p;
    // print solution    
    fprintf(stderr, 
	    "LP1\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\n",
	    objval, x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7]);
    // has the solution changed?
    if (fabs(objval - objvalsave) < EPSILON) {
      objcompcount++;
    } else {
      objcompcount = 0;
    }
    if (objcompcount > MAXITN) {
      fprintf(stderr, 
	      "pooling-slp: objfun constant for %d itns, terminating\n", 
	      MAXITN);
      termination = TRUE;
    }
    if (is_equal(x, xsave, 1)) {
      compareflag1 = TRUE;
    } else {
      for(j = 0; j < n; j++) {
	xsave[j] = x[j];
      }
    }
    // improved solution?
    if (objval < bestobj) {
      bestobj = objval;
      for(j = 0; j < n; j++) {
	xast[j] = x[j];
      }
      noprogresscount = 0;
    }
    // values for y to be used in next iteration
    y[0] = x[3];
    y[1] = x[4];

    // update lp2
    // constr 1, rhs = y11+y12
    CPXchgcoef(env, lp2, 0, -1, y[0] + y[1]);
    // constr 3, rhs = 100-y11
    CPXchgcoef(env, lp2, 2, -1, 100 - y[0]);
    // constr 4, rhs = 200-y12
    CPXchgcoef(env, lp2, 3, -1, 200 - y[1]);
    // constr 5, coeff of p = -(y11+y12)
    CPXchgcoef(env, lp2, 4, 5, -y[0] - y[1]);
    // constr 6, coeff of p = y11, rhs = 2.5*y11
    CPXchgcoef(env, lp2, 5, 5, y[0]);
    CPXchgcoef(env, lp2, 5, -1, 2.5*y[0]);
    // constr 7, coeff of p = y12, rhs = 1.5*y12
    CPXchgcoef(env, lp2, 6, 5, y[1]);
    CPXchgcoef(env, lp2, 6, -1, 1.5*y[1]);

    // optimize lp2
    status = CPXlpopt(env, lp2);
    if (status != 0) {
      fprintf(stderr, "%s: failed to call optimization algorithm on lp1\n", 
	      argv[0]);
      exit(10);
    }
    objvalsave = objval;
    // solve lp2, this sets x as follows:
    // x[0]:x11, x[1]:x12, x[2]:x21, x[3]:y21, x[4]:y22, x[5]:p
    CPXsolution(env, lp2, &status, &objval, x, lambda, s, rc);
    // reindex x so that it includes also y11, y12
    x[7] = x[5];
    x[6] = x[4];
    x[5] = x[3];
    x[4] = y[1];
    x[3] = y[0];
    // add missing term to objective
    objval += (-9*y[0] - 15*y[1]);
    // print solution
    fprintf(stderr, 
	    "LP2\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\n",
	    objval, x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7]);
    // has the solution changed?
    if (fabs(objval - objvalsave) < EPSILON) {
      objcompcount++;
    } else {
      objcompcount = 0;
    }
    if (objcompcount > MAXITN) {
      fprintf(stderr, 
	      "pooling-slp: objfun constant for %d itns, terminating\n", 
	      MAXITN);
      termination = TRUE;
    }
    if (is_equal(x, xsave, 2)) {
      compareflag2 = TRUE;
    } else {
      for(j = 0; j < n; j++) {
	xsave[j] = x[j];
      }
    }
    // improved solution?
    if (objval < bestobj) {
      bestobj = objval;
      for(j = 0; j < n; j++) {
	xast[j] = x[j];
      }
      noprogresscount = 0;
    }
    // find values for p to be used in next iteration
    p = x[7];

    // update lp1
    // constr 5, coeff of y11 = -p, coeff of y12 = -p
    CPXchgcoef(env, lp1, 4, 3, -p);
    CPXchgcoef(env, lp1, 4, 4, -p);
    // constr 6, coeff of y11 = (p - 2.5)
    CPXchgcoef(env, lp1, 5, 3, p-2.5);
    // constr 7, coeff of y12 = (p - 1.5)
    CPXchgcoef(env, lp1, 6, 4, p-1.5);

    // termination condition on solution not changing
    if (compareflag1 == TRUE && compareflag2 == TRUE) {
      fprintf
	(stderr, 
	 "pooling-slp: solution has not changed in this itn, terminating\n");
      termination = TRUE;
    }
    
    // termination condition on no progress for too long
    noprogresscount++;
    if (noprogresscount > MAXITN) {
      fprintf
	(stderr, 
	 "pooling-slp: no progress for > %d itns, terminating\n", MAXITN); 
      termination = TRUE;
    }
  }
  fprintf(stderr, "*******************************************************\n");

  // print solution
  printf("%s: solution of Haverly's pooling problem:\n", argv[0]);
  printf("  last solver return status = %d\n", status);
  printf("  best objective function value = %f\n", objval);
  printf("  best primal variables:\n");
  printf("    x11 = %f\n", x[0]);
  printf("    x12 = %f\n", x[1]);
  printf("    x21 = %f\n", x[2]);
  printf("    y11 = %f\n", x[3]);
  printf("    y12 = %f\n", x[4]);
  printf("    y21 = %f\n", x[5]);
  printf("    y22 = %f\n", x[6]);
  printf("    p   = %f\n", x[7]);

  // free storage
  CPXfreeprob(env, &lp1);
  CPXfreeprob(env, &lp2);
  return ret;
}


////////////////// rubbish bin //////////////////

#ifdef OBLIVION
int is_feasible(double* x) {
  if (fabs(x[0] + x[2] - x[3] - x[4]) < EPSILON &&
      fabs(x[1] - x[5] - x[6]) < EPSILON &&
      x[3] + x[5] <= 100 + EPSILON &&
      x[4] + x[6] <= 200 + EPSILON &&
      fabs(x[0] + x[2] - x[7]*x[3] - x[7]*x[4]) < EPSILON &&
      (x[7] - 2.5)*x[3] - 0.5*x[5] <= EPSILON &&
      (x[7] - 1.5)*x[4] + 0.5*x[6] <= EPSILON) {
    return TRUE;
  } else {
    return FALSE;
  }
}
#endif