Recursion

Recursion, more properly known as Successive Linear Programming, is a technique whereby LP may be used to solve certain non-linear problems. Some coefficients in an LP problem are defined to be functions of the optimal values of LP variables. When an LP problem has been solved, the coefficients are re-evaluated and the LP re-solved. Under some assumptions this process may converge to a local (though not necessarily a global) optimum.

Example problem

Consider the following financial planning problem: We wish to determine the yearly interest rate x so that for a given set of payments we obtain the final balance of 0. Interest is paid quarterly according to the following formula:

The balance at time t (t=1,...,T) results from the balance of the previous period t-1 and the net of payments and interest:

Model formulation

This problem cannot be modeled just by LP because we have the T products

which are non-linear. To express an approximation of the original problem by LP we replace the interest rate variable x by a (constant) guess X of its value and a deviation variable dx

The formula for the quarterly interest payment i_t therefore becomes

where balance_t is the balance at the beginning of period t.

We now also replace the balance balance_t-1 in the product with dx by a guess B_t-1 and a deviation db_t-1

which can be approximated by dropping the product of the deviation variables

To ensure feasibility we add penalty variables eplus_t and eminus_t for positive and negative deviations in the formulation of the constraint:

The objective of the problem is to get feasible, that is to minimize the deviations:

Implementation

The Mosel model (file recurse.mos) then looks as follows (note the balance variables balance_t as well as the deviation dx and the quarterly nets net_t are defined as free variables, that is, they may take any values between minus and plus infinity):

model Recurse
 uses "mmxprs"

 forward procedure solve_recurse

 declarations
  T=6                                ! Time horizon
  QUARTERS=1..T                      ! Range of time periods
  P,R,V: array(QUARTERS) of real     ! Payments
  B: array(QUARTERS) of real         ! Initial guess as to balances b(t)
  X: real                            ! Initial guess as to interest rate x

  interest: array(QUARTERS) of mpvar ! Interest
  net: array(QUARTERS) of mpvar      ! Net
  balance: array(QUARTERS) of mpvar  ! Balance
  x: mpvar                           ! Interest rate
  dx: mpvar                          ! Change to x
  eplus, eminus: array(QUARTERS) of mpvar  ! + and - deviations
 end-declarations

 X:= 0.00
 B:: [1, 1, 1, 1, 1, 1]
 P:: [-1000, 0, 0, 0, 0, 0]
 R:: [206.6, 206.6, 206.6, 206.6, 206.6, 0]
 V:: [-2.95, 0, 0, 0, 0, 0]
                                     ! net = payments - interest
 forall(t in QUARTERS) net(t) = (P(t)+R(t)+V(t)) - interest(t)	

                                     ! Money balance across periods
 forall(t in QUARTERS) balance(t) = if(t>1, balance(t-1), 0) - net(t)	

 forall(t in 2..T) Interest(t):=     ! Approximation of interest
   -(365/92)*interest(t) + X*balance(t-1) + B(t-1)*dx + eplus(t) - eminus(t) = 0

 Def:= X + dx = x                    ! Define the interest rate: x = X + dx

 Feas:= sum(t in QUARTERS) (eplus(t)+eminus(t))  ! Objective: get feasible

 interest(1) = 0                     ! Initial interest is zero
 forall (t in QUARTERS) net(t) is_free
 forall (t in 1..T-1) balance(t) is_free
 balance(T) = 0                      ! Final balance is zero
 dx is_free

 minimize(Feas)                      ! Solve the LP-problem

 solve_recurse                       ! Recursion loop

                                     ! Print the solution
 writeln("\nThe interest rate is ", getsol(x))
 write(strfmt("t",5), strfmt(" ",4))
 forall(t in QUARTERS) write(strfmt(t,5), strfmt(" ",3))
 write("\nBalances ")
 forall(t in QUARTERS)  write(strfmt(getsol(balance(t)),8,2))
 write("\nInterest ")
 forall(t in QUARTERS)  write(strfmt(getsol(interest(t)),8,2))

end-model

In the model above we have declared the procedure solve_recurse that executes the recursion but it has not yet been defined. The recursion on x and the balance_t (t=1,...,T-1) is implemented by the following steps:

(a) The B_t-1 in constraints Interest_t get the prior solution value of balance_t-1
(b) The X in constraints Interest_t get the prior solution value of x
(c) The X in constraint Def gets the prior solution value of x

We say we have converged when the change in dx (variation) is less than 0.000001 (TOLERANCE). By setting Mosel's comparison tolerance to this value the test variation > 0 checks whether variation is greater than TOLERANCE.

 procedure solve_recurse
  declarations
   TOLERANCE=0.000001                 ! Convergence tolerance
   variation: real                    ! Variation of x
   BC: array(QUARTERS) of real
   bas: basis                         ! LP basis
  end-declarations

  setparam("zerotol", TOLERANCE)      ! Set Mosel comparison tolerance
  variation:=1.0
  ct:=0

  while(variation>0) do
   savebasis(bas)                     ! Save the current basis
   ct+=1
   forall(t in 2..T)
     BC(t-1):= getsol(balance(t-1))   ! Get solution values for balance(t)'s
   XC:= getsol(x)                     ! and x
   write("Round ", ct, " x:", getsol(x), " (variation:", variation,"), ")
   writeln("Simplex iterations: ", getparam("XPRS_SIMPLEXITER"))

   forall(t in 2..T) do               ! Update coefficients
     Interest(t)+= (BC(t-1)-B(t-1))*dx
     B(t-1):=BC(t-1)
     Interest(t)+= (XC-X)*balance(t-1)
   end-do
   Def+= XC-X
   X:=XC
   oldxval:=XC                        ! Store solution value of x

   loadprob(Feas)                     ! Reload the problem into the optimizer
   loadbasis(bas)                     ! Reload previous basis
   minimize(Feas)                     ! Re-solve the LP-problem

   variation:= abs(getsol(x)-oldxval) ! Change in dx
  end-do
 end-procedure 

With the initial guesses 0 for X and 1 for all B_t the model converges to an interest rate of 5.94413% (x = 0.0594413).