Recurse - A successive linear programming model

/********************************************************
  Xpress-BCL C++ Example Problems
  ===============================

  file xbrecurs.cxx
  `````````````````
  Recursion solving a non-linear financial planning problem
  The problem is to solve
  	net(t) = Payments(t)  - interest(t)
  	balance(t) = balance(t-1) - net(t)
  	interest(t) = (92/365) * balance(t) * interest_rate
  where
        balance(0) = 0
        balance[T] = 0
  for interest_rate

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

#include <iostream>
#include <cmath>
#include "xprb_cpp.h"
#include "xprs.h"

using namespace std;
using namespace ::dashoptimization;

#define T 6

/****DATA****/
double X = 0.00;                /* An INITIAL GUESS as to interest rate x */
double B[] = /* {796.35, 589.8918, 398.1351, 201.5451, 0.0, 0.0}; */
             {1,1,1,1,1,1};
                                /* An INITIAL GUESS as to balances b(t) */
double P[] = {-1000, 0, 0, 0, 0, 0};                 /* Payments */
double R[] = {206.6, 206.6, 206.6, 206.6, 206.6, 0}; /*  "       */
double V[] = {-2.95, 0, 0, 0, 0, 0};                 /*  "       */

XPRBvar b[T];                   /* Balance */
XPRBvar x;                      /* Interest rate */
XPRBvar dx;                     /* Change to x */
XPRBctr interest[T], ctrd;

XPRBprob p("Fin_nlp");          /* Initialize a new problem in BCL */
 
/***********************************************************************/

void modFinNLP()
{
 XPRBvar i[T];                  /* Interest */
 XPRBvar n[T];                  /* Net */
 XPRBvar epl[T], emn[T];        /* + and - deviations */
 XPRBexpr cobj, le;
 int t;
 
/****VARIABLES****/
 for(t=0;t<T;t++)
 {
  i[t]=p.newVar("i");
  n[t]=p.newVar("n", XPRB_PL, -XPRB_INFINITY, XPRB_INFINITY);
  b[t]=p.newVar("b", XPRB_PL, -XPRB_INFINITY, XPRB_INFINITY);
  epl[t]=p.newVar("epl");
  emn[t]=p.newVar("emn");
 }
 x=p.newVar("x");
 dx=p.newVar("dx", XPRB_PL, -XPRB_INFINITY, XPRB_INFINITY);
 i[0].fix(0);
 b[T-1].fix(0);
  
/****OBJECTIVE****/
 for(t=0;t<T;t++)               /* Objective: get feasible */
  cobj += epl[t]+emn[t];
 p.setObj(cobj);                /* Select objective function */ 
 p.setSense(XPRB_MINIM);        /* Choose the sense of the optimization */   

/****CONSTRAINTS****/
 for(t=0;t<T;t++)
 {                              /* net = payments - interest */
  p.newCtr("net", n[t] == P[t]+R[t]+V[t] - i[t]);
                                /* Money balance across periods */
  if(t>0) p.newCtr("bal", b[t] == b[t-1]-n[t]);
  else  p.newCtr("bal", b[t] == -n[t]);
 }

 for(t=1;t<T;t++)    /* i(t) = (92/365)*( b(t-1)*X + B(t-1)*dx ) approx. */
  interest[t] =
   p.newCtr("int", X*b[t-1] + B[t-1]*dx + epl[t] - emn[t] == 365/92.0 * i[t]);
 
 ctrd = p.newCtr("def", x == X + dx);
}

/**************************************************************************/
/*  Recursion loop (repeat until variation of x converges to 0):          */
/*    save the current basis and the solutions for variables b[t] and x   */
/*    set the balance estimates B[t] to the value of b[t]                 */
/*    set the interest rate estimate X to the value of x                  */
/*    reload the problem and the saved basis                              */
/*    solve the LP and calculate the variation of x                       */
/**************************************************************************/
void solveFinNLP()
{
 XPRBbasis basis;
 double variation=1.0, oldval, BC[T], XC;
 int t, ct=0;

 XPRSsetintcontrol(p.getXPRSprob(), XPRS_CUTSTRATEGY, 0); 
                                 /* Switch automatic cut generation off */

 p.lpOptimize("");               /* Solve the LP-problem */
 
 while(variation>0.000001)
 {
  ct++;
  basis=p.saveBasis();           /* Save the current basis */

                                 /* Get the solution values for b and x */
  for(t=1;t<T;t++) BC[t-1]=b[t-1].getSol();
  XC=x.getSol();
  cout << "Loop " << ct << ": " << x.getName() << ":" << x.getSol();
  cout << "  (variation: " << variation << ")" << endl;
                   
  for(t=1;t<T;t++)               /* Change coefficients in interest[t] */
  {   
   interest[t].setTerm(dx, BC[t-1]);
   B[t-1]=BC[t-1];
   interest[t].setTerm(b[t-1], XC);
  }
  ctrd.setTerm(XC);              /* Change constant term of ctrd */ 
  X=XC;

  oldval=XC;
  p.loadMat();                   /* Reload the problem */
  p.loadBasis(basis);            /* Load the saved basis */
  basis.reset();                 /* No need to keep the basis any longer */

  p.lpOptimize("");              /* Solve the LP-problem */
  variation=fabs(x.getSol()-oldval);
 }

 cout << "Objective: " << p.getObjVal() << endl; /* Get objective value */
 cout << "Interest rate: " << x.getSol()*100 << "%" << endl << "Balances:  ";
 for(t=0;t<T;t++)                /* Print out the solution values */
  cout << b[t].getName() << ":" << b[t].getSol() << " ";
 cout << endl;
} 

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

int main(int argc, char **argv)
{
 modFinNLP();                    /* Model the problem */
 solveFinNLP();                  /* Solve the problem */
  
 return 0;
}