The Decomposition Algorithm

The most natural formulation for this model is in the form of a mixed integer nonlinear program (MINLP). Let denote the set of ATMs and denote the set of days used in the training data. The predictive model fit is defined by the following data for each ATM on each day : . The model fit parameters define the variables for each ATM that, when applied to the predictive model, give the estimated cash flow need per day, per ATM. In addition, define a surrogate variable for each ATM on each day that defines the cash inventory (replenished from the branch) minus withdrawals given by the fit. The variable also represents the error in the regression model. Let define the budget per day, define the limit on cash-outs per ATM, and define the historical withdrawals at a particular ATM on a particular day. Then the following MINLP models this problem:

Constraint (cash) defines the surrogate variable , which gives the estimated net cash flow. Inequalities (budget) and (count) ensure that the solution satisfies the budget and cash-out constraints, respectively.

To express this model in a more standard form, you can first use some standard model reformulations to linearize the absolute value and the cash-out constraint (count).

A well-known reformulation for linearizing the absolute value of a variable is to introduce one variable for each side of the absolute value. The following systems are equivalent:

is equivalent to

Let and represent the positive and negative parts, respectively, of the net cash flow . Then, you can rewrite the model, removing the absolute value, as the following:

To count the number of times a cash-out occurs, you need to introduce a binary variable to keep track of when this event occurs. Let be an indicator variable that takes value 1 when the net cash flow is negative. You can model the implication , or its contrapositive , by adding the constraint

Now, you can model the cash-out constraint by counting the number of days the net-cash flow is negative for each ATM, as follows:

The MINLP model can now be written as follows:

This MINLP is difficult to solve, in part because the prediction function is not convex. Another approach is to use mixed integer linear programming (MILP) to formulate an approximation of the problem, as described in the next section.

Since the predictive model is a forecast, finding the optimal parameters that are based on nondeterministic data is not of primary importance. Rather, you want to provide as good a solution as possible in a reasonable amount of time. So, using MILP to approximate the MINLP is perfectly acceptable. In the original problem you have products of two continuous variables that are both bounded by (lower bound) and (upper bound). This arrangement enables you to create an approximate linear model using a few standard modeling reformulations.

The first step is to discretize one of the continuous variables . The goal is to transform the product of a continuous variable with another continuous variable instead to a continuous variable with a binary variable. By doing this, you can linearize the product form.

You must assume some level of approximation by defining a binary variable for each possible setting of the continuous variable to be chosen from some discrete set. For example, if you let , then you allow to be chosen from the set . Let represent the possible steps and . Then, you apply the following transformation to variable :

The MINLP model can now be approximated as the following:

You still need to linearize the product terms in the cash flow constraint. Since these terms are products of a bounded continuous variable and a binary variable, you can linearize them by introducing for each product another variable , which serves as a surrogate. In general, you know the following relationship between the original variables and their surrogates:

is equivalent to

Using this relationship to replace each product form, you now can write the problem as an approximate MILP as follows:

Since it is difficult to solve the MINLP model directly, the approximate MILP formulation is attractive. Unfortunately, the size of the approximate MILP is much bigger than the associated MINLP. Direct methods for solving this MILP do not work well. However, the problem is nicely suited for the decomposition algorithm.

When you examine the structure of the MILP model, you see clearly that the constraints can be easily decomposed by ATM. In fact, the only set of constraints that involve decision variables across ATMs is the budget constraint (budget). That is, if you relax constraint (budget), you are left with independent blocks of constraints, one for each ATM.

To show how this is done in PROC OPTMODEL, consider the following data sets which describe an example with 20 ATMs over a period of 100 days. This particular example was submitted to MIPLIB 2010, which is a collection of difficult MILPs in the public domain (Koch et al., 2011).

The first data set, budget_data, provides the cash budget on each particular day:

data budget_data;
   input d $ budget;
   datalines;
DATE0      70079
DATE1      66418
DATE10     52656
DATE11     50439
DATE12     58688
DATE13     45002
DATE14     52369
...
;

The second data set, cashout_data, provides the limit on the number of cash-outs allowed at each ATM:

data cashout_data;
   input a $ cashOutLimit;
   datalines;
ATM0            31
ATM1            24
ATM2            41
ATM3            43
ATM4            29
ATM5            24
ATM6            52
ATM7            44
ATM8            35
ATM9            48
ATM10           31
ATM11           47
ATM12           26
ATM13           34
ATM14           29
ATM15           32
ATM16           33
ATM17           32
ATM18           43
ATM19           28
;

The final data set, polyfit_data, provides the polynomial fit coefficients for each ATM on each date. It also provides the historical cash withdrawals.

data polyfit_data;
   input a $ d $ cx cy cz cu c withdrawal;
   datalines;
ATM0     DATE0       2822     1984    -1984    1045    1373        780
ATM0     DATE1       1337     2530    -2530    1510     174       2351
ATM0     DATE2       2685      -67       67     145    2820       2288
ATM0     DATE3       -595    -3135     3135     581    3319       1357
...
ATM19    DATE96      -734     3392    -3392     162    1648        914
ATM19    DATE97     -1062      969     -969     444    1746       2264
ATM19    DATE98      7676     2308    -2308      59    1388        972
ATM19    DATE99      3062     1308    -1308    1080     654        698
;

The following PROC OPTMODEL statements read in the data and define the necessary sets and parameters:

proc optmodel;
   set<str> DATES;
   set<str> ATMS;

   /* cash budget per date */
   num budget{DATES};
   
   /* maximum amount of cash-outs allowed at each atm */
   num cashOutLimit{ATMS};
   
   /* historical withdrawal amount per atm each date */
   num withdrawal{ATMS, DATES};

   /* polynomial fit coefficients for predicted cash flow needed */
   num c {ATMS, DATES};
   num cx{ATMS, DATES};
   num cy{ATMS, DATES};
   num cz{ATMS, DATES};
   num cu{ATMS, DATES};

   /* number of points used in approximation of continuous range */
   num nSteps = 10; 
   set STEPS = {0..nSteps};

   read data budget_data into DATES=[d] budget;
   read data cashout_data into ATMS=[a] cashOutLimit;
   read data polyfit_data into [a d] cx cy cz cu c withdrawal;      

The following statements declare the variables:

   var x{ATMS,STEPS}                 binary;
   var v{ATMS,DATES}                 binary;
   var z{ATMS,STEPS}                 >= 0 <= 1;
   var y{ATMS}                       >= 0 <= 1;
   var u{ATMS}                       >= 0;
   var fPlus{ATMS,DATES}             >= 0;
   var fMinus{a in ATMS, d in DATES} >= 0 <= withdrawal[a,d];     

The following statements declare the objective and constraints:

   min CashFlowDiff = 
      sum{a in ATMS, d in DATES} (fPlus[a,d] + fMinus[a,d]);

   con BudgetCon{d in DATES}:
      sum{a in ATMS} (fPlus[a,d] - fMinus[a,d] + withdrawal[a,d]) 
          <= budget[d];

   con CashFlowDefCon{a in ATMS, d in DATES}:
      cx[a,d] * sum{t in STEPS} (t/nSteps) * x[a,t] +
      cy[a,d] * y[a]                                +
      cz[a,d] * sum{t in STEPS} (t/nSteps) * z[a,t] +
      cu[a,d] * u[a]                                +
      c[a,d] - withdrawal[a,d] = fPlus[a,d] - fMinus[a,d];

   con PickOneStepCon{a in ATMS}:
      sum{t in STEPS} x[a,t] = 1;

   con CashOutLinkCon{a in ATMS, d in DATES}:
      fMinus[a,d] <= withdrawal[a,d] * v[a,d];

   con CashOutLimitCon{a in ATMS}:
      sum{d in DATES} v[a,d] <= cashOutLimit[a];

   con Linear1Con{a in ATMS, t in STEPS}:
      z[a,t] <= x[a,t];
   
   con Linear2Con{a in ATMS}:
      sum{t in STEPS} z[a,t] = y[a];         

The following statements define the block decomposition by ATM. The .block suffix expects numeric indices while the ATMS specified in the PROC OPTMODEL statements contain strings. You can create a mapping from the string identifier to a numeric identifier as follows:

   /* create numeric block index */
   num blockIndex {ATMS};
   num index init 0;
   for{a in ATMS} do;
      blockIndex[a] = index;
      index = index + 1;
   end;

Then, each constraint can be added to its associated ATM block as follows:

   /* define blocks for each ATM */
   for{a in ATMS} do;
      PickOneStepCon[a].block  = blockIndex[a];
      CashOutLimitCon[a].block = blockIndex[a];
      Linear2Con[a].block      = blockIndex[a];
      for{d in DATES} do;
         CashFlowDefCon[a,d].block = blockIndex[a];
         CashOutLinkCon[a,d].block = blockIndex[a];
      end;
      for{t in STEPS} do;
         Linear1Con[a,t].block = blockIndex[a];
      end;
   end;     

The budget constraint links all ATMs and remains in the master problem. Finally, the following statements use DECOMP to solve the problem:

   /* set the number of threads and get performance details */
   performance details nthreads=12;

   /* solve with the decomposition algorithm */
   solve with milp / relobjgap=1e-5 decomp=();
quit;

The solution summary, performance information, and procedure task timing tables are displayed in Output 13.5.1.

The iteration log, which contains the problem statistics, the progress of the solution, and the optimal objective value, is shown in Output 13.5.2.