Recurse - A successive linear programming model

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

  file xbrecurs.c
  ```````````````
  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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "xprb.h"
#include "xprs.h"

#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;

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

void modfinnlp(XPRBprob prob)
{
 XPRBvar i[T];                  /* Interest */
 XPRBvar n[T];                  /* Net */
 XPRBvar epl[T], emn[T];        /* + and - deviations */
 XPRBctr cobj, bal, net;
 int t;

/****VARIABLES****/
 for(t=0;t<T;t++)
 {
  i[t]=XPRBnewvar(prob,XPRB_PL, "i", 0, XPRB_INFINITY);
  n[t]=XPRBnewvar(prob,XPRB_PL, "n", -XPRB_INFINITY, XPRB_INFINITY);
  b[t]=XPRBnewvar(prob,XPRB_PL, "b", -XPRB_INFINITY, XPRB_INFINITY);
  epl[t]=XPRBnewvar(prob,XPRB_PL, "epl", 0, XPRB_INFINITY);
  emn[t]=XPRBnewvar(prob,XPRB_PL, "emn", 0, XPRB_INFINITY);
 }
 x=XPRBnewvar(prob,XPRB_PL, "x", 0, XPRB_INFINITY);
 dx=XPRBnewvar(prob,XPRB_PL, "dx", -XPRB_INFINITY, XPRB_INFINITY);
 XPRBfixvar(i[0],0);
 XPRBfixvar(b[T-1],0);

/****OBJECTIVE****/
 cobj = XPRBnewctr(prob,"OBJ",XPRB_N);  /* Objective: get feasible */
 for(t=0;t<T;t++)
 {
  XPRBaddterm(cobj,epl[t],1);
  XPRBaddterm(cobj,emn[t],1);
 }
 XPRBsetobj(prob,cobj);                 /* Select objective function */
 XPRBsetsense(prob,XPRB_MINIM);         /* Choose the sense of the optimization */

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

 for(t=1;t<T;t++)
 {                   /* i(t) = (92/365)*( b(t-1)*X + B(t-1)*dx ) approx. */
  interest[t]=XPRBnewctr(prob,"int", XPRB_E);
  XPRBaddterm(interest[t], i[t], -(365/92.0));
  XPRBaddterm(interest[t], b[t-1], X);
  XPRBaddterm(interest[t], dx, B[t-1]);
  XPRBaddterm(interest[t], epl[t], 1);
  XPRBaddterm(interest[t], emn[t], -1);
 }

 ctrd=XPRBnewctr(prob,"def", XPRB_E);    /* X + dx = x */
 XPRBaddterm(ctrd, NULL, X);
 XPRBaddterm(ctrd, x, 1);
 XPRBaddterm(ctrd, dx, -1);
}

/**************************************************************************/
/*  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(XPRBprob prob)
{
 XPRBbasis basis;
 double variation=1.0, oldval, BC[T], XC;
 int t, ct=0;

 XPRSsetintcontrol(XPRBgetXPRSprob(prob), XPRS_CUTSTRATEGY, 0);
                                   /* Switch automatic cut generation off */

 XPRBlpoptimize(prob,"");          /* Solve the LP-problem */

 while(variation>0.000001)
 {
  ct++;
  basis=XPRBsavebasis(prob);       /* Save the current basis */

                                   /* Get the solution values for b and x */
  for(t=1;t<T;t++) BC[t-1]=XPRBgetsol(b[t-1]);
  XC=XPRBgetsol(x);
  printf("Loop %d: %s:%g  (variation:%g)\n", ct, XPRBgetvarname(x),
    XPRBgetsol(x),variation);

  for(t=1;t<T;t++)                 /* Change coefficients in interest[t] */
  {
   XPRBsetterm(interest[t], dx, BC[t-1]);
   B[t-1]=BC[t-1];
   XPRBsetterm(interest[t], b[t-1], XC);
  }
  XPRBsetterm(ctrd, NULL, XC);     /* Change constant term of ctrd */
  X=XC;

  oldval=XC;
  XPRBloadmat(prob);               /* Reload the problem */
  XPRBloadbasis(basis);            /* Load the saved basis */
  XPRBdelbasis(basis);             /* No need to keep the basis any longer */
  XPRBlpoptimize(prob,"");         /* Solve the LP-problem */
  variation=fabs(XPRBgetsol(x)-oldval);
 }

 printf("Objective: %g\n", XPRBgetobjval(prob));  /* Get objective value */
 printf("Interest rate: %g%%\nBalances:  ", XPRBgetsol(x)*100);
 for(t=0;t<T;t++)                  /* Print out the solution values */
  printf("%s:%g ", XPRBgetvarname(b[t]), XPRBgetsol(b[t]));
 printf("\n");
}

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

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

 prob=XPRBnewprob("Fin_nlp");      /* Initialize a new problem in BCL */
 modfinnlp(prob);                  /* Model the problem */
 solvefinnlp(prob);                /* Solve the problem */

 return 0;
}