Extensions to Linear Programming

Topics covered in this chapter:

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 (SLP), 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.

This section describes the complete implementation of a recursion algorithm with Mosel using an LP solver. Section Xpress NonLinear and Xpress Global shows how the example problem discussed in this section can be implemented and solved in a significantly simpler manner with Xpress NonLinear using the SLP solver.

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:

interest_t = (92/365) ·balance_t ·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:

net_t = Payments_t - interest_t
balance_t = balance_t-1 - net_t

Model formulation

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

balance_t ·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 i_t therefore becomes

interest_t = 92/365 ·(balance_t-1 ·x)
= 92/365 ·(balance_t-1 ·(X + dx))
= 92/365 ·(balance_t-1 ·X + balance_t-1 ·dx)

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

iinterest_t = 92/365 ·(balance_t-1 ·X + (B_t-1+db_t-1) ·dx)
= 92/365 ·(balance_t-1 ·X + B_t-1 ·dx + db_t-1 ·dx)

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:

interest_t = 92/365 ·(balance_t-1 ·X + B_t-1 ·dx + eplus_t - eminus_t)

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

minimize

∑

t ∈ QUARTERS

(eplus_t + eminus_t)

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 solverecurse

 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

 solverecurse                          ! Recursion loop

                                       ! Print the solution
 writeln("\nThe interest rate is ", x.sol)
 write(strfmt("t",5), strfmt(" ",4))
 forall(t in QUARTERS) write(strfmt(t,5), strfmt(" ",3))
 write("\nBalances ")
 setparam("realfmt", "%8.2f")
 forall(t in QUARTERS) write(balance(t).sol)
 write("\nInterest ")
 forall(t in QUARTERS) write(interest(t).sol)

end-model

In the model above we have declared the procedure solverecurse 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 solverecurse
   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).

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 first 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.

When the goals are specified as objective functions the resulting problem is often refered to as a multi-objective optimization problem. Xpress Optimizer supports the specification of multiple objective functions along with some configuration options, within Mosel this functionality is accessible through the mmxprs module as shown in the second example in section.

Example problem

The objective is to solve a problem with two variables x and y (x,y ≥ 0), the constraint

100·x + 60·y ≤ 600

and the three goal constraints

Goal₁: 7·x + 3·y ≥ 40
Goal₂: 10·x + 5·y = 60
Goal₃: 5·x + 4·y ≥ 35

where the order given corresponds to their priorities.

Implementation

To increase readability, the implementation of the Mosel model (file goalctrp.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 printsol produces a nice solution printout.

model GoalCtrPreempt
 uses "mmxprs", "mmsystem"

 forward procedure preemptive
 forward procedure printsol(i:integer)

 declarations
   NGOALS=3                            ! Number of goals
   GOALS=1..NGOALS                     ! Set of goals
   x,y: mpvar                          ! Decision variables
   dev: array(1..2*NGOALS) of mpvar    ! Deviation from goals
   MinDev: linctr                      ! Objective function
   Goal: array(GOALS) 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₁), a ≥ 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 dev₂. The goal constraint becomes 7·x + 3·y + dev₂ ≥ 40 and the objective function is to minimize dev₂. 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 Goal₂: since this is an equation, we need variables for positive and negative deviations, dev₃ and dev₄. We add dev₄-dev₃ to the constraint, turning it into 10·x + 5·y +dev₄-dev₃ = 60 and the objective now is to minimize the absolute deviation dev₄+dev₃. And so on.

 procedure preemptive

! Remove (=hide) goal constraint from the problem
   forall(i in GOALS) 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: ", x.sol, ", y: ", y.sol)

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

   printsol(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 printsol for printing the results.

 procedure printsol(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", textfmt("Target",11), textfmt("Value",7))
   forall(g in 1..i)
     writeln(formattext("%4d %3s %5g %7g", g, ATypes(Goal(g).type),
       -Goal(g).coeff, Goal(g).act-dev(2*g).sol+dev(2*g-1).sol))

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

end-model

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

Multi-objective optimization with Xpress Optimizer

The minimize/maximize routines of mmxprs have additional forms that allow the user to specify multiple objective functions along with configuration settings for each objective. The objconfig constructor used in the example below takes as its arguments the priority (all equal values to indicate the Archimedian case, otherwise the pre-emptive case applies where objectives with higher priority values will be considered first), a weight value (negative values inverse the sense of the objective), followed by absolute or relative tolerance values to be applied in the pre-emptive case when transforming an objective into a constraint after a solution has been determined for it.

Specific solution reporting functionality for multi-objective solves includes the number of objectives for which the problem has been solved XPRS_SOLVEDOBJS (this value may be lower than the total number of specified objectives if the final problems status is 'infeasible') and the possibility to access summary information about the iterations of the pre-emptive solving via getobjrealattr or getobjintattr.

model GoalObjPreempt
 uses "mmxprs", "mmsystem"

 declarations
   NGOALS=3                            ! Number of goals
   GOALS=1..NGOALS                     ! Set of goals
   x,y: mpvar                          ! Variables
   Goal: array(GOALS) of linctr        ! Goal objective expressions
   ObjCfg: array(GOALS) of objconfig   ! Goal configuration
 end-declarations

 Limit:= 42*x + 13*y <= 100            ! Define a constraint

! Define the goal objectives
 Goal(1):=   5*x +  2*y - 20
 Goal(2):=  -3*x + 15*y - 48
 Goal(3):= 1.5*x + 21*y - 3.8

! Configuration of the goal objectives
 ObjCfg(1):= objconfig(3, -1, 0, 0.1)  ! 1st, maximize, 10% tolerance
 ObjCfg(2):= objconfig(2,  1, 4, 0)    ! 2nd, minimize, absolute tolerance of 4
 ObjCfg(3):= objconfig(1, -1, 0, 0.2)  ! 3rd, maximize, 20% tolerance

 minimize(Goal,ObjCfg)                 ! Solve the multi-objective problem

! Solution reporting
 writeln("Solved for ", getparam("XPRS_SOLVEDOBJS"), " goals")

! Solution values in iterations of pre-emptive solving
 forall(i in GOALS)
   writeln(" Iteration ", i, ": val=", getobjrealattr(i, "XPRS_LPOBJVAL"))

! Final solution
 writeln(" Solution: x: ", x.sol, ", y: ", y.sol)
 writeln(" Goal", textfmt("Target",9), textfmt("Value",12))
 forall(g in GOALS)
   writeln(formattext("%4d %8s %12.6f", g, if(ObjCfg(g).weight=1,"min","max"),
      Goal(g).sol))
end-model