Initializing help system before first use

# Extensions to Linear Programming

Topics covered in this chapter:

- Recursion
  - Example problem
  - Model formulation
  - Implementation
- Goal Programming
  - Example problem
  - Implementation
  - Multi-objective optimization with Xpress Optimizer

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:

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

∑
t ∈ 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 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 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 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 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 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

Goal1: 7·x + 3·y ≥ 40
Goal2: 10·x + 5·y = 60
Goal3: 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 (Goal1), 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 dev2. The goal constraint becomes 7·x + 3·y + dev2 ≥ 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 dev4. We add dev4-dev3 to the constraint, turning it into 10·x + 5·y +dev4-dev3 = 60 and the objective now is to minimize the absolute deviation dev4+dev3. 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

Parent Topic Advanced language features

© 2001-2026 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.