[Previous chapter] [Next chapter]

Extensions to Linear Programming


The two examples (recursion and Goal Programming) in this chapter show how Mosel can be used to implement extensions of Linear Programming.

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:

interestt = (92/365) ·balancet ·interest_rate

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:

nett = Paymentst - interestt
balancet = balancet-1 - nett

Model formulation

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

balancet ·interest_rate

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

x = X + dx

The formula for the quarterly interest payment it therefore becomes

interestt = 92/365 ·(balancet-1 ·x)
= 92/365 ·(balancet-1 ·(X + dx))
= 92/365 ·(balancet-1 ·X + balancet-1 ·dx)

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

We now also replace the balance balancet-1 in the product with dx by a guess Bt-1 and a deviation dbt-1

iinterestt = 92/365 ·(balancet-1 ·X + (Bt-1+dbt-1) ·dx)
= 92/365 ·(balancet-1 ·X + Bt-1 ·dx + dbt-1 ·dx)

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

interestt = 92/365 ·(balancet-1 ·X + Bt-1 ·dx)

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

interestt = 92/365 ·(balancet-1 ·X + Bt-1 ·dx + eplust - eminust)

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

minimize
Maths/sum.png
t Maths/insm.png QUARTERS
(eplust + eminust)

Implementation

The Mosel model (file recurse.mos) then looks as follows (note the balance variables balancet as well as the deviation dx and the quarterly nets nett 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 balancet (t=1,...,T-1) is implemented by the following steps:

(a) The Bt-1 in constraints Interestt get the prior solution value of balancet-1
(b) The X in constraints Interestt 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 Bt the model converges to an interest rate of 5.94413% (x = 0.0594413).

Goal Programming

Goal Programming is an extension of Linear Programming in which targets are specified for a set of constraints. In Goal Programming there are two basic models: the pre-emptive (lexicographic) model and the Archimedian model. In the pre-emptive model, goals are ordered according to priorities. The goals at a certain priority level are considered to be infinitely more important than the goals at the next level. With the Archimedian model weights or penalties for not achieving targets must be specified, and we attempt to minimize the sum of the weighted infeasibilities.

If constraints are used to construct the goals, then the goals are to minimize the violation of the constraints. The goals are met when the constraints are satisfied.

The example in this section demonstrates how Mosel can be used for implementing pre-emptive Goal Programming with constraints. We try to meet as many goals as possible, taking them in priority order.

Example problem

The objective is to solve a problem with two variables x and y (x,y Maths/geq.png 0), the constraint

100·x + 60·y Maths/leq.png 600

and the three goal constraints

Goal1: 7·x + 3·y Maths/geq.png 40
Goal2: 10·x + 5·y = 60
Goal3: 5·x + 4·y Maths/geq.png 35

where the order given corresponds to their priorities.

Implementation

To increase readability, the implementation of the Mosel model (file goalctr.mos) is organized into three blocks: the problem is stated in the main part, procedure preemptive implements the solution strategy via preemptive Goal Programming, and procedure print_sol produces a nice solution printout. 

model GoalCtr
 uses "mmxprs"

 forward procedure preemptive
 forward procedure print_sol(i:integer)

 declarations
  NGOALS=3                          ! Number of goals
  x,y: mpvar                        ! Decision variables
  dev: array(1..2*NGOALS) of mpvar  ! Deviation from goals
  MinDev: linctr                    ! Objective function
  Goal: array(1..NGOALS) of linctr  ! Goal constraints
 end-declarations
 
 100*x + 60*y <= 600                ! Define a constraint
 
! Define the goal constraints
 Goal(1):=  7*x + 3*y >= 40
 Goal(2):= 10*x + 5*y = 60
 Goal(3):=  5*x + 4*y >= 35
  
 preemptive                         ! Pre-emptive Goal Programming

At the end of the main part, we call procedure preemptive to solve this problem by pre-emptive Goal Programming. In this procedure, the goals are at first entirely removed from the problem (`hidden'). We then add them successively to the problem and re-solve it until the problem becomes infeasible, that is, the deviation variables forming the objective function are not all 0. Depending on the constraint type (obtained with gettype) of the goals, we need one (for inequalities) or two (for equalities) deviation variables.

Let us have a closer look at the first goal (Goal_1), a Maths/geq.png constraint: the left hand side (all terms with decision variables) must be at least 40 to satisfy the constraint. To ensure this, we add a dev2. The goal constraint becomes 7·x + 3·y + dev2 Maths/geq.png 40 and the objective function is to minimize dev2. If this is feasible (0-valued objective), then we remove the deviation variable from the goal, thus turning it into a hard constraint. It is not required to remove it from the objective since minimization will always force this variable to take the value 0.

We then continue with Goal2: since this is an equation, we need variables for positive and negative deviations, dev3 and dev_4. We add dev_4-dev3 to the constraint, turning it into 10·x + 5·y +dev_4-dev3 = 60 and the objective now is to minimize the absolute deviation dev_4+dev3. And so on. 

 procedure preemptive

! Remove (=hide) goal constraint from the problem
  forall(i in 1..NGOALS) sethidden(Goal(i), true)

  i:=0
  while (i<NGOALS) do
    i+=1
    sethidden(Goal(i), false)       ! Add (=unhide) the next goal
    
    case gettype(Goal(i)) of        ! Add deviation variable(s)
     CT_GEQ: do
              Deviation:= dev(2*i)
              MinDev += Deviation 
             end-do
     CT_LEQ: do
              Deviation:= -dev(2*i-1)
              MinDev += dev(2*i-1) 
             end-do
     CT_EQ : do
              Deviation:= dev(2*i) - dev(2*i-1)
              MinDev += dev(2*i) + dev(2*i-1) 
             end-do
     else    writeln("Wrong constraint type")
             break
    end-case  
    Goal(i)+= Deviation
   
    minimize(MinDev)                ! Solve the LP-problem
    writeln(" Solution(", i,"): x: ", getsol(x), ", y: ", getsol(y))

    if getobjval>0 then
     writeln("Cannot satisfy goal ",i)
     break
    end-if  
    Goal(i)-= Deviation             ! Remove deviation variable(s) from goal
  end-do

  print_sol(i)                      ! Solution printout
 end-procedure

The procedure sethidden(c:linctr, b:boolean) has already been used in the previous chapter (Section Column generation) without giving any further explanation. With this procedure, constraints can be removed (`hidden') from the problem solved by the optimizer without deleting them in the problem definition. So effectively, the optimizer solves a subproblem of the problem originally stated in Mosel.

To complete the model, below is the procedure print_sol for printing the results. 

 procedure print_sol(i:integer)
  declarations
   STypes={CT_GEQ, CT_LEQ, CT_EQ}
   ATypes: array(STypes) of string
  end-declarations
 
  ATypes::([CT_GEQ, CT_LEQ, CT_EQ])[">=", "<=", "="]

  writeln(" Goal", strfmt("Target",11), strfmt("Value",7))
  forall(g in 1..i) 
   writeln(strfmt(g,4), strfmt(ATypes(gettype(Goal(g))),4), 
     strfmt(-getcoeff(Goal(g)),6), 
     strfmt( getact(Goal(g))-getsol(dev(2*g))+getsol(dev(2*g-1)) ,8)) 

  forall(g in 1..NGOALS)
   if (getsol(dev(2*g))>0) then
    writeln(" Goal(",g,") deviation from target: -", getsol(dev(2*g)))
   elif (getsol(dev(2*g-1))>0) then
    writeln(" Goal(",g,") deviation from target: +", getsol(dev(2*g-1))) 
   end-if 
 end-procedure

end-model

When running the program, one finds that the first two goals can be satisfied, but not the third.


[Next] [Index]
If you have any comments or suggestions about these pages, please send mail to docs@dashoptimization.com.