(!*******************************************************
  * Mosel Example Problems                              *
  * ======================                              *
  *                                                     *
  * file recurse.mos                                    *
  * ````````````````                                    *
  * Example for the use of the Mosel language           *
  * (Financial application of Recursion, that is        *
  *  Non-linear programming)                            *  
  *                                                     *
  *  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                               *
  *                                                     *
  *  The problem is that we have the T products:        *
  *     balance(t) * interest_rate                      *
  *  which cannot be modelled just by LP.               *
  *                                                     *
  * (c) 2008 Fair Isaac Corporation                     *
  *     author: S. Heipcke, 2001                        *
  *******************************************************!)

model Fin_nlp                       ! Start a new model

uses "mmxprs"                       ! Load the optimizer library

forward procedure solverec          ! Declare a procedure defined later

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

 i: array(RT) of mpvar              ! Interest
 n: array(RT) of mpvar              ! Net
 b: array(RT) of mpvar              ! Balance
 x: mpvar                           ! Interest rate
 dx: mpvar                          ! Change to x
 epl, emn: array(RT) 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 RT) net(t):= n(t) = (P(t)+R(t)+V(t)) - i(t)

                                    ! Money balance across periods
 forall(t in RT) bal(t):= b(t) = if(t>1, b(t-1), 0) - n(t)

!       interest = (92/365)*( balance * interest_rate)
! i.e.  i(t)     = (92/365)*( b(t-1)  * x )
!       i(t)     = (92/365)*( b(t-1)  * ( X + dx ) )
!       i(t)     = (92/365)*( b(t-1)*X + (B(t-1)+db(t-1))*dx )
!       i(t)     = (92/365)*( b(t-1)*X + B(t-1)*dx ) (approx)
!       Use penalty variables (epl and emn) to ensure feasibility:
 forall(t in 2..T) 
  interest(t):= -(365/92)*i(t) + X*b(t-1) + B(t-1)*dx + epl(t) - emn(t) = 0
                          
 def:= X + dx = x                    ! Define the interest rate: x = X + dx
                                     
 feas:= sum(t in RT) (epl(t)+emn(t)) ! Objective: get feasible

 i(1) = 0                            ! Initial interest is zero
 forall (t in RT) n(t) is_free
 forall (t in 1..T-1) b(t) is_free
 b(T) = 0                            ! Final balance is zero
 dx is_free

 minimize(feas)                      ! Solve the LP-problem

 solverec                            ! Recursion loop
 
                                     ! Print the solution
 writeln("\nThe interest rate is ", x.sol)
 write(strfmt("t",5), strfmt(" ",4))
 forall(t in RT) write(strfmt(t,5), strfmt(" ",3))
 write("\nBalances ")
 forall(t in RT)  write(strfmt(b(t).sol,8,2))
 write("\nInterest ")
 forall(t in RT)  write(strfmt(i(t).sol,8,2))
 writeln

 
!************************************************************************
!  Recursion loop: we recurse on X and the b(t)'s.
!    The 'B(t-1)' in rows interest(t) get the prior value of b(t-1)
!    The 'X' in rows interest(t) get the prior value of x
!    The 'X' in row def gets the prior value of x
!    We say we have converged when the change in dx is less than 1.0E-6
!************************************************************************
procedure solverec

 declarations
  TOLERANCE=0.000001                  ! Convergence tolerance
  variation: real                     ! Variation of x
  BC: array(RT) of real
  bas: 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):= b(t-1).sol            ! Get the solution values for b(t)'s
    XC:= x.sol                        ! and x
    write("Round ", ct, " x:", x.sol, " (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)*b(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(x.sol-oldxval)    ! Change in dx
  end-do
 
end-procedure 

end-model