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

  file xbels.c
  ````````````
  Economic lot sizing, ELS, problem, solved by adding
  (l,S)-inequalities) in several rounds looping over
  the root node.

  ELS considers production planning over a horizon
  of T periods. In period t, t=1,...,T, there is a
  given demand DEMAND[t] that must be satisfied by
  production prod[t] in period t and by inventory
  carried over from previous periods. There is a
  set-up up cost SETUPCOST[t] associated with
  production in period t. The unit production cost
  in period t is PRODCOST[t]. There is no inventory
  or stock-holding cost.

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

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

#define EPS    1e-6

#define T 6                             /* Number of time periods */

/****DATA****/
int DEMAND[]    = { 1, 3, 5, 3, 4, 2};  /* Demand per period */
int SETUPCOST[] = {17,16,11, 6, 9, 6};  /* Setup cost per period */
int PRODCOST[]  = { 5, 3, 2, 1, 3, 1};  /* Production cost per period */
int D[T][T];                            /* Total demand in periods t1 - t2 */

XPRBvar prod[T];                        /* Production in period t */
XPRBvar setup[T];                       /* Setup in period t */

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

void mod_els(XPRBprob prob)
{
 int s,t,k;
 XPRBctr ctr;

 for(s=0;s<T;s++)
  for(t=0;t<T;t++)
   for(k=s;k<=t;k++)
    D[s][t] += DEMAND[k];

/****VARIABLES****/
 for(t=0;t<T;t++)
 {
  prod[t]=XPRBnewvar(prob,XPRB_PL, XPRBnewname("prod%d",t+1),0,XPRB_INFINITY);
  setup[t]=XPRBnewvar(prob,XPRB_BV, XPRBnewname("setup%d",t+1),0,1);
 }

/****OBJECTIVE****/
 ctr = XPRBnewctr(prob,"OBJ",XPRB_N);   /* Minimize total cost */
 for(t=0;t<T;t++)
 {
  XPRBaddterm(ctr, setup[t], SETUPCOST[t]);
  XPRBaddterm(ctr, prod[t], PRODCOST[t]);
 }
 XPRBsetobj(prob,ctr);

/****CONSTRAINTS****/
         /* Production in period t must not exceed the total demand for the
            remaining periods; if there is production during t then there
            is a setup in t */
 for(t=0;t<T;t++)
 {
  ctr = XPRBnewctr(prob,"Production",XPRB_L);
  XPRBaddterm(ctr, setup[t], -D[t][T-1]);
  XPRBaddterm(ctr, prod[t], 1);
 }

         /* Production in periods 0 to t must satisfy the total demand
            during this period of time */
 for(t=0;t<T;t++)
 {
  ctr = XPRBnewctr(prob,"Demand",XPRB_G);
  for(s=0;s<=t;s++)
   XPRBaddterm(ctr, prod[s], 1);
  XPRBaddterm(ctr, NULL, D[0][t]);
 }

}

/**************************************************************************/
/*  Cut generation loop at the top node:                                  */
/*    solve the LP and save the basis                                     */
/*    get the solution values                                             */
/*    identify and set up violated constraints                            */
/*    load the modified problem and load the saved basis                  */
/**************************************************************************/
void solve_els(XPRBprob prob)
{
 double objval;                  /* Objective value */
 int t,l;
 int starttime;
 int ncut, npass, npcut;         /* Counters for cuts and passes */
 double solprod[T], solsetup[T]; /* Solution values for var.s prod & setup */
 double ds;
 XPRBbasis basis;
 XPRBctr cut;

 starttime=XPRBgettime();
 XPRSsetintcontrol(XPRBgetXPRSprob(prob),XPRS_CUTSTRATEGY, 0);
                                 /* Disable automatic cuts - we use our own */
 XPRSsetintcontrol(XPRBgetXPRSprob(prob),XPRS_PRESOLVE, 0);
                                 /* Switch presolve off */
 ncut = npass = 0;

 do
 {
  npass++;
  npcut = 0;
  XPRBlpoptimize(prob,"p");      /* Solve the LP */
  basis=XPRBsavebasis(prob);     /* Save the current basis */
  objval = XPRBgetobjval(prob);  /* Get the objective value */

      /* Get the solution values: */
  for(t=0;t<T;t++)
  {
   solprod[t]=XPRBgetsol(prod[t]);
   solsetup[t]=XPRBgetsol(setup[t]);
  }

      /* Search for violated constraints: */
  for(l=0;l<T;l++)
  {
   for (ds=0.0, t=0; t<=l; t++)
   {
    if(solprod[t] < D[t][l]*solsetup[t] + EPS)  ds += solprod[t];
    else  ds += D[t][l]*solsetup[t];
   }

      /* Add the violated inequality: the minimum of the actual production
         prod[t] and the maximum potential production D[t][l]*setup[t]
         in periods 0 to l must at least equal the total demand in periods
         0 to l.
         sum(t=1:l) min(prod[t], D[t][l]*setup[t]) >= D[0][l]
       */
   if(ds < D[0][l] - EPS)
   {
    cut = XPRBnewctr(prob,XPRBnewname("cut%d",ncut+1), XPRB_G);
    XPRBaddterm(cut, NULL, D[0][l]);
    for(t=0;t<=l;t++)
    {
     if (solprod[t] < D[t][l]*solsetup[t] + EPS)
      XPRBaddterm(cut, prod[t], 1);
     else
      XPRBaddterm(cut, setup[t], D[t][l]);
    }
    ncut++;
    npcut++;
   }
  }

  printf("Pass %d (%g sec), objective value %g, cuts added: %d (total %d)\n",
   npass, (XPRBgettime()-starttime)/1000.0, objval, npcut, ncut);

  if(npcut==0)
   printf("Optimal integer solution found:\n");
  else
  {
   XPRBloadmat(prob);             /* Reload the problem */
   XPRBloadbasis(basis);          /* Load the saved basis */
   XPRBdelbasis(basis);           /* No need to keep the basis any longer */
  }
 } while(npcut>0);

      /* Print out the solution: */
 for(t=0;t<T;t++)
  printf("Period %d: prod %g (demand: %d, cost: %d), setup %g (cost: %d)\n",
      t+1, XPRBgetsol(prod[t]), DEMAND[t], PRODCOST[t], XPRBgetsol(setup[t]),
      SETUPCOST[t]);
}

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

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

 prob=XPRBnewprob("Els");         /* Initialize a new problem in BCL */
 mod_els(prob);                   /* Model the problem */
 solve_els(prob);                 /* Solve the problem */

 return 0;
}