Multiple models and parallel
solving with Mosel

FICOTM Xpress Optimization Suite whitepaper

Last update 4 December, 2009

www.fico.com

Make every decision countTM

Mosel: multiple models and parallel solving

Multiple models and parallel solving with Mosel
Y. Colombani and S. Heipcke
Xpress Team, FICO, Leam House, Leamington Spa CV32 5YN, UK
http://www.fico.com/xpress
4 December, 2009

Abstract
This paper describes several examples of sequential and parallel solving of multiple models with Mosel. Without being able to give an exhaustive list of possible configurations, the examples showcase different uses of theMosel module mmjobs, such as concurrent execution of several instances of a model, the (sequential) embedding of a submodel into a master, and the implementation of decomposition algorithms (Dantzig-Wolfe and
Benders decomposition).
From a more technical point of view, topics discussed in this paper include model management, synchronizationof concurrent models, and the use of the shared memory IO driver.

Contents
1
2

3

4

5

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1
Multi-problem vs. multi-model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Basic tasks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1
Executing a submodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
Stopping the submodel execution . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3
Output from the submodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4
Compilation to memory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5
Runtime parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6Running several submodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1
Sequential submodels . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.2
Parallel submodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7
Communication of data between different models . . . . . . . . . . . . . . . .
2.7.1
Using the shared memory driver . . . . . . . . . . . . . . . . . . . . . .
2.7.2
Using the memory pipe driver . . . . . . . . . . . . . . . . . . . . . . .
Column generation: solving different models in sequence . . . . . . . . . . . . . . . .
3.1
Example problem: cutting stock . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1
Master model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2
Knapsack model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Solving several model instances in parallel . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
Example problem: economic lot sizing . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1
Cutting plane algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1
Master model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2
ELS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dantzig-Wolfe decomposition: combining sequential and parallel solving . . . . . . .

c 2009 Fair Isaac Corporation. All rights reserved.

2
3
4
4
5
6
6
7
8
8
8
9
9
10
11
12
13
13
15
15
16
16
17
18
18
19
22
22

page 1

Mosel: multiple models and parallel solving

5.1

Example problem: multi-item, multi-period production planning . . . . . . . .
5.1.1
Original model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2
Problem decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1
Master model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2
Single factory model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3
Master problem subroutines . . . . . . . . . . . . . . . . . . . . . . . .
5.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
Benders decomposition: working with several different submodels . . . . . . . . . . .
6.1
A small example problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1
Master model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.2
Submodel 1: fixed continuous variables . . . . . . . . . . . . . . . . . .
6.2.3
Submodel 2: fixed integer variables . . . . . . . . . . . . . . . . . . . .
6.2.4
Submodel 0: start solution . . . . . . . . . . . . . . . . . . . . . . . . .
6.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

24
24
25
27
27
30
33
35
36
37
38
38
41
42
44
45
46
46

Introduction
Release 1.6 of Mosel introduces the possibility to work with multiple models directly in the Mosel
language. The new functionality, provided by the module mmjobs, includes facilities for model
management, synchronization of concurrent models based on event queues, and a sharedmemory IO driver. The following list gives an overview on the available functionality. For thecomplete documentation of this module the reader is referred to the section mmjobs of the‘Mosel Reference Manual’.
• Model management: compilation of source model files, loading of bim files, model
execution and interruption, retrieval of model information (status, exit code, ID), redirectionof IO streams.
• Synchronization mechanism: sending and retrieving events, waiting for events or eventclasses, retrieval of event information (class, value, sender model).
• Shared memory IO driver: shared memory version of the mem driver for exchanging data
between concurrent models (write access by a single model, read access by several models
simultaneously), usable wherever Mosel expects a (generalized) filename, in particular ininitializations blocks.
• Memory pipe IO driver: memory pipes for exchanging data between concurrent models(write access by several models, read access by a single model), usable wherever Mosel
expects a (generalized) filename, in particular in initializations blocks.
mmjobs introduces two new types, Model and Event. The type Model is used to reference a
Mosel model. Before using the reference to a model it has to be initialized by loading a bim file.
The type Event represents an event in the Mosel language. Events are characterized by a classand a value and may be exchanged between a model and its parent model. An event queue isattached to each model to collect all events sent to this model and is managed with a FIFO policy(First In – First Out).
The first section of this paper introduces the reader to the basic tasks that are typically performed
when working with several models in Mosel, namely:

Introduction

c 2009 Fair Isaac Corporation. All rights reserved.

page 2

Mosel: multiple models and parallel solving• Executing a submodel from a master model: the compile – load – run – wait sequence• Stopping the submodel execution• Output from the submodel: redirection to a file (especially under IVE), making thesubmodel silent• Compilation to memory• Passing runtime parameters• Running several submodels– in sequence– in parallel• Communication of data between different models: using the shared memory and memory
pipe IO drivers
The remainder of the paper gives some more advanced examples of the use of mmjobs with adetailed explanation of their implementation. All examples are available for download from theXpress website.
• Column generation: re-implementation of the column generation example from the Mosel
User Guide with two separate models that are solved sequentially, passing data via sharedmemory.
• Parallel solving: several instances of the same model are run concurrently with different
solution algorithm parameterizations. Improved solution values are sent for bound updates
to all running models and the first model that finishes stops all others.
• Dantzig-Wolfe decomposition: an iterative sequence of concurrent solving of a set of
subproblem instances, followed by solving of the updated master problem; data exchange
via shared memory.
• Benders decomposition: solving iteratively a sequence of several different subproblems;
data exchange via shared memory.

1.1

Multi-problem vs. multi-model
At this place we would like to stress the difference between multiple models and multiple
problems — Mosel releases prior to version 3.0 always associate a single problem with every
model. This means, for instance, if a model contains several calls to a solver such as
Xpress-Optimizer, then the solver will work with a single problem representation, and only thesolution to the last optimization run can be obtained from the solver at any time.
Release 3.0 of Mosel introduces the possibility of defining several problems within a single model.
At any point a single problem is active. It is possible to switch back and forth between problems,e.g., to retrieve solution information for a decision variable from different problems. To
represent a ’problem’, the Mosel language defines the type mpproblem, modules can provide
their own problem types and users may equally create new problem types on the language level.
Here is a comparison of the main characteristics of multi-problem and multi-model
implementations with Mosel.

Introduction

c 2009 Fair Isaac Corporation. All rights reserved.

page 3

Mosel: multiple models and parallel solving

Multi-problem (mpproblem)

Multi-model (mmjobs)

Single model file

Several model files–

problems share data–

communication of data (in memory)–

integrated; no direct access to(sub)problems by other models/applications–

stand-alone execution of submodels
or use of submodels with other master models/applications

Sequential access to problems only

Sequential or parallel execution of models

Problem solving approaches that involve parallel execution of (sub)models can only be
implemented as multiple models, whereas sequential solving can be formulated with either one.
For sequential algorithms (such as in Column generation and Benders decomposition) the
developer may choose among the two design options.

2

Basic tasks
This section introduces the basic tasks that typically need to be performed when working withseveral models in Mosel using the functionality of the module mmjobs.

2.1

Executing a submodel
Assume we are given the following simple model testsub.mos (the submodel) that we wish to
run from a second model (its master model):model "Test submodel"
forall(i in 10..20) write(i^2, " ")writelnend-model

The reader is certainly familiar with the standard compile-load-run sequence that is always
required for the execution of a Mosel model independent of the place from where it is executed(be it the Mosel command line, a host application, or another Mosel model). In the case of a
Mosel model executing a second model, we need to augment this standard sequence tocompile-load-run-wait. The rationale behind this is that the submodel is started as a separate
thread and we thus make sure that the submodel has terminated before the master model ends
(or continues its execution with results from the submodel, see Section 2.7 below).
The following model runtestsub.mos uses the compile-load-run-wait sequence in its basic
form to execute the model testsub.mos printed above. The corresponding Mosel subroutines
are defined by the module mmjobs that needs to be included with a uses statement. The
submodel remains as is (that is, there is no need to include mmjobs if the submodel itself does notuse any functionality of this module). The last statement of the master model, dropnextevent,
may require some further explanation: the termination of the submodel is indicated by an event(of class EVENT_END) that is sent to the master model. The wait statement pauses the executionof the master model until it receives an event. Since in the present case the only event sent by the
submodel is this termination message we simply remove the message from the event queue ofthe master model without checking its nature.
model "Run model testsub"
uses "mmjobs"

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 4

Mosel: multiple models and parallel solving

declarations
modSub: Model
end-declarations! Compile the model file
if compile("testsub.mos")<>0 then exit(1); end-if
load(modSub, "testsub.bim")! Load the bim filerun(modSub)! Start model execution
wait! Wait for model termination
dropnextevent! Ignore termination event message
end-model

The two models are run by executing the master model in any standard way of executing Mosel
models (from the Mosel command line, within Xpress-IVE, or from a host application). With theMosel command line you may use, for instance, the following command:
mosel -c "exec runtestsub.mos"

As a result you should see the following output printed by the submodel (for readers working
under Xpress-IVE: please see Section 2.3):
100 121 144 169 196 225 256 289 324 361 400

If we wish to obtain more precise information about the termination status of the submodel we
could replace the statement dropnextevent by the following lines that retrieve the event sentby the submodel and print out its class (the termination event has the predefined class
EVENT_END) and the value attached to it (here the default value 0). In addition, we display theexit code of the model (value sent by an exit statement terminating the model execution or thedefault value 0).
declarations
ev: Eventend-declarations
ev:=getnextevent
writeln("Event class: ", getclass(ev))
writeln("Event value: ", getvalue(ev))
writeln("Exit code : ", getexitcode(modSub))

2.2

Stopping the submodel execution
If a submodel execution takes a long time it may be desirable to interrupt the submodel without
stopping the master model itself. The following modified version of our master model (filerunsubwait.mos) shows how this can be achieved by adding a duration (in seconds) to the waitstatement. If the submodel has not yet sent the termination event message after executing for
one second it is stopped by the call to stop with the model reference.
model "Run model testsub"
uses "mmjobs"
declarations
modSub: Model
end-declarations! Compile the model file
if compile("testsub.mos")<>0 then exit(1); end-if
load(modSub, "testsub.bim")! Load the bim filerun(modSub)! Start model execution
wait(1)! Wait 1 second for an event

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 5

Mosel: multiple models and parallel solving

if isqueueempty then
! No event has been sent: model still runs
writeln("Stopping the submodel")stop(modSub)! Stop the model
wait! Wait for model termination
end-ifdropnextevent! Ignore termination event message
end-model

2.3

Output from the submodel
Readers who have tried to execute the master models from the previous sections within
Xpress-IVE will not have seen any output from the submodel since IVE only captures the output
from the master model. The best way to obtain access to the submodel output in this case is toredirect this output to a file. Redirecting the submodel output may also become necessary whenseveral (sub)models are run in parallel as the output on screen is likely to mix up output printed
by several of these models.
Output may be redirected directly within the submodel with statements like the following added
to the model before any output is printed:
fopen("testout.txt", F_OUTPUT+F_APPEND) ! Output to file (in append mode)
fopen("tee:testout.txt&", F_OUTPUT)! Output to file and on screen
fopen("null:", F_OUTPUT)! Disable all output

where the first line redirects the output to the file testout.txt, the second statement
maintains the output on screen while writing to the file at the same time, and the third line
makes the model entirely silent. The output file is closed by adding the statement
fclose(F_OUTPUT) after the printing statements in the model.
The same can be achieved from the master model by adding output redirection before the runstatement for the corresponding submodel, such as:
setdefstream(modSub, F_OUTPUT, "testout.txt") ! Output to filesetdefstream(modSub, F_OUTPUT, "tee:testout.txt&")! Output to file and on screen
setdefstream(modSub, F_OUTPUT, "null:")! Disable all output

The output redirection for a submodel may be terminated by resetting its output stream to thedefault output:
setdefstream(modSub, F_OUTPUT, "")2.4

Compilation to memory
The default compilation of a Mosel file filename.mos generates a binary model file
filename.bim. To avoid the generation of physical BIM files for submodels we may compile thesubmodel to memory, as shown in the following example runsubmem.mos. Working in memory
usually is more efficient than accessing physical files. Furthermore, this feature will also be
helpful if you do not have write access at the place where the master model is executed.
model "Run model testsub"
uses "mmjobs", "mmsystem"
declarations
modSub: Model
end-declarations! Compile the model file

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 6

Mosel: multiple models and parallel solving

if compile("", "testsub.mos", "shmem:testsubbim")<>0 thenexit(1)
end-ifload(modSub, "shmem:testsubbim") ! Load the bim file from memory
fdelete("shmem:testsubbim")! ... and release the memory blockrun(modSub)! Start model execution
wait! Wait for model termination
dropnextevent! Ignore termination event message
end-model

The full version of compile takes three arguments: the compilation flags (e.g., use "g" fordebugging), the model file name, and the output file name (here a label prefixed by the name ofthe shared memory driver). Having loaded the model we may free the memory used by thecompiled model with a call to fdelete (this subroutine is provided by the module mmsystemthat needs to be loaded in addition to mmjobs).

2.5

Runtime parameters
A convenient means of modifying data in a Mosel model when running the model (that is,
without having to modify the model itself and recompile it) is to use runtime parameters. Suchparameters are declared at the beginning of the model in a parameters block where every
parameter is given a default value that will be applied if no other value is specified for thisparameter at the model execution.
Consider the following model rtparams.mos that may receive parameters of four different
types—integer, real, string, and Boolean—and prints out their values.
model "Runtime parameters"
parameters
PARAM1 = 0
PARAM2 = 0.5
PARAM3 = ’’
PARAM4 = false
end-parameters
writeln(PARAM1, "", PARAM2, "", PARAM3, "", PARAM4)

end-model

The master model runrtparam.mos executing this (sub)model may look as follows—all runtime
parameters are given new values:
model "Run model rtparams"
uses "mmjobs"
declarations
modPar: Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
load(modPar, "rtparams.bim")! Load the bim file
! Start model execution
run(modPar, "PARAM1=" + 2 + ",PARAM2=" + 3.4 +",PARAM3=’a string’" + ",PARAM4=" + true)
wait! Wait for model termination
dropnextevent! Ignore termination event message
end-model

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 7

Mosel: multiple models and parallel solving

2.6

Running several submodels
Once we have seen how to run a parameterized model from a master model it is only a small stepto the execution of several different submodel instances from a master model. The two following
sections deal with the two cases of sequential and parallel execution of submodels. For
simplicity’s sake, the submodels in our examples are all parameterized versions of a single model.
It is of course equally possible to compile, load, and run different submodel files from a single
master model.
2.6.1

Sequential submodels

Running several instances of a submodel in sequence only requires small modifications to themaster model that we have used for a single model instance as can be seen from the following
example (file runrtparamseq.mos)—to keep things simple, we now only reset a single
parameter at every execution:model "Run model rtparams in sequence"
uses "mmjobs"
declarations
A = 1..10
modPar: Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
load(modPar, "rtparams.bim")! Load the bim file
forall(i in A) dorun(modPar, "PARAM1=" + i)
wait
dropnextevent
end-do! Start model execution! Wait for model termination! Ignore termination event message

end-model

The submodel is compiled and loaded once and after starting the execution of a submodel
instance we wait for its termination before the next instance is started.
2.6.2

Parallel submodels

The parallel execution of submodel (instances) requires slightly more modifications. We still haveto compile the submodel only once, but it now needs to be loaded as many times as we want to
run parallel instances. The wait statements are now moved to a separate loop since we first want
to start all submodels and then wait for their termination.
model "Run model rtparams in parallel"
uses "mmjobs"
declarations
A = 1..10
modPar: array(A) of Model
end-declarations! Compile the model file
if compile("rtparams.mos")<>0 then exit(1); end-if
forall(i in A) doload(modPar(i), "rtparams.bim")run(modPar(i), "PARAM1=" + i)

Basic tasks! Load the bim file
! Start model execution

c 2009 Fair Isaac Corporation. All rights reserved.

page 8

Mosel: multiple models and parallel solving

end-doforall(i in A) do
wait
dropnextevent
end-do! Wait for model termination! Ignore termination event message

end-model

The order in which the submodel output appears on screen is nondeterministic because themodels are run in parallel. However, since the submodel execution is very quick, this may not
become obvious: try adding the line wait(1) to the submodel rtparams.mos immediately
before the writeln statement (you will also need to add the statement uses "mmjobs" at the
beginning of the model) and compare the output of several runs of the master model. You are
now likely to see different output sequences with every run.

2.7

Communication of data between different models
Runtime parameters are a means of communicating single data values to a submodel but they are
not suited, for instance, to pass data tables or sets to a submodel. Also, they cannot be employed
to retrieve any information from a submodel or to exchange data between models during their
execution. All these tasks are addressed by the two IO drivers defined by the module mmjobs: theshmem driver and the mempipe driver. As was already stated earlier, the shmem driver is meant for
one-to-many communication (one model writing, many reading) and the mempipe driver serves
for many-to-one communication. In the case of one model writing and one model reading wemay use either, where shmem is conceptionally probably the easier to use.
2.7.1

Using the shared memory driverWith the shmem shared memory driver we write and read data blocks from/to memory. The use ofthis driver is quite similar to the way we would work with physical files. We have alreadyencountered an example of its use in Section 2.4: the filename is replaced by a label, prefixed by
the name of the driver, such as "mmjobs.shmem:aLabel". If the module mmjobs is loaded by
the model (or another model held in memory at the same time, such as its master model) we may
use the short form "shmem:aLabel".
The exchange of data between different models is carried out through initializations blocks.
In general, the shmem driver will be combined with the raw driver to save data in binary format.
Let us now take a look at a modified version testsubshm.mos of our initial test submodel(Section 2.1). This model reads in the index range from memory and writes back the resulting
array to memory:model "Test submodel"
declarations
A: range
B: array(A) of real
end-declarations
initializations from "raw:"
A as "shmem:A"
end-initializations
forall(i in A) B(i):= i^2
initializations to "raw:"
B as "shmem:B"
end-initializations

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 9Mosel: multiple models and parallel solving

end-model

The master model runsubshm.mos to run this submodel may look as follows:model "Run model testsubshm"
uses "mmjobs"
declarations
modSub: Model
A = 30..40
B: array(A) of real
end-declarations! Compile the model file
if compile("testsubshm.mos")<>0 then exit(1); end-if
load(modSub, "testsubshm.bim")! Load the bim fileinitializations to "raw:"
A as "shmem:A"
end-initializations
run(modSub)
wait
dropnextevent! Start model execution! Wait for model termination! Ignore termination event message

initializations from "raw:"
B as "shmem:B"
end-initializations
writeln(B)end-model

Before the submodel run is started the index range is written to memory and after its end we
retrieve the result array to print it out from the master model.
If memory blocks are no longer used it is recommended to free up these blocks by calling
fdelete (subroutine provided by module mmsystem), especially if the data blocks are large, since
the data blocks survive the termination of the model that has created them, even if the model isunloaded explicitly, until the module mmjobs is unloaded (explicitly or by the termination of theMosel session). At the end of our master model we might thus add the lines
fdelete("shmem:A")
fdelete("shmem:B")

2.7.2

Using the memory pipe driver

The memory pipe IO driver mempipe works in the opposite way to what we have seen for theshared memory driver: a pipe first needs to be opened before it can be written to. That means
we need to call initializations from before initializations to. The submodel (filetestsubpip.mos) now looks as follows:model "Test submodel"
declarations
A: range
B: array(A) of real
end-declarations
initializations from "mempipe:indata"A
end-initializations

Basic tasks

c 2009 Fair Isaac Corporation. All rights reserved.

page 10

Mosel: multiple models and parallel solving

forall(i in A) B(i):= i^2
initializations to "mempipe:resdata"B
end-initializations
end-modelThis is indeed not very different from the previous submodel. However, there are more changes
to the master model (file runsubpip.mos): the input data is now written by the master model
after the submodel run has started. The master model then opens a new pipe to read the result
data before the submodel terminates its execution.
model "Run model testsubpip"
uses "mmjobs"
declarations
modSub: Model
A = 30..40
B: array(A) of real
end-declarations! Compile the model file
if compile("testsubpip.mos")<>0 then exit(1); end-if
load(modSub, "testsubpip.bim")! Load the bim filerun(modSub)! Start model execution

initializations to "mempipe:indata"A
end-initializations
initializations from "mempipe:resdata"B
end-initializations
wait
dropnextevent! Wait for model termination! Ignore termination event message

writeln(B)end-model

Once a file has opened a pipe for reading it remains blocked in this state until it has received therequested data through this pipe. The program control flow is therefore the following in thepresent case:
1. The master model starts the submodel.
2. The submodel opens the input data pipe and waits for the master to write to it.
3. Once the input data has been communicated the submodel continues its execution while
the master model opens the result data pipe and waits for the results.
4. When the result data pipe is open, the submodel writes to the result data pipe and thenterminates.
5. The master model prints out the result.

3

Column generation: solving different models in sequence
The cutting stock example we are working with in this section is taken from the ‘Mosel User

Column generation

c 2009 Fair Isaac Corporation. All rights reserved.

page 11Mosel: multiple models and parallel solving

Guide’. The reader is refered to this manual for further detail on the column generation
algorithm and its implementation with Mosel.
Column generation algorithms are typically used for solving linear problems with a huge numberof variables for which it is not possible to generate explicitly all columns of the problem matrix.
Starting with a very restricted set of columns, after each solution of the problem a column
generation algorithm adds one or several columns that improve the current solution.
Our column generation algorithm for the cutting stock problem requires us to solve a knapsackproblem based on the dual value of the current solution to determine a new column (= cutting
pattern). The difference between the User Guide implementation and the one shown below
consists in the handling of this knapsack (sub)problem. In the User Guide implementation Mosel’sconstraint hiding functionality is used to blend out subsets of constraints; in the version shown
below the subproblem is implemented in a model on its own. Both versions implement exactly
the same algorithm and their performance is comparable. On larger instances, however, thetwo-model version is likely to be slightly more efficient, since every model defines exactly the
problem to be solved, without any selection of (un)hidden constraints.
In this example, the changes to the problems are such that they cause complete re-loading of the
problems for every optimization run. A clearer advantage of the multi-model version would show
up if there were only slight changes (bound updates) to the main (cutting stock) problem so that
this problem did not have to be reloaded into the solver for every new run.

3.1

Example problem: cutting stock
A paper mill produces rolls of paper of a fixed width MAXWIDTH that are subsequently cut into
smaller rolls according to the customer orders. The rolls can be cut into NWIDTHS different sizes.
The orders are given as demands for each width i (DEMANDi ). The objective of the paper mill is
to satisfy the demand with the smallest possible number of paper rolls in order to minimize thelosses.
The objective of minimizing the total number of rolls can be expressed as choosing the best set of
cutting patterns for the current set of demands. Since it may not be obvious how to calculate all
possible cutting patterns by hand, we start off with a basic set of patterns (PATTERNS1 ,...,PATTERNSNWIDTH ), that consists of cutting small rolls all of the same width as many times as
possible out of the large roll.
If we define variables usej to denote the number of times a cutting pattern j(j ∈ WIDTHS = {1, . . . , NWIDTH}) is used, then the objective becomes to minimize the sum of thesevariables, subject to the constraints that the demand for every size has to be met.
X
minimizeusejj∈WIDTHSXPATTERNSij · usej ≥ DEMANDij∈WIDTHS∀j ∈ WIDTHS : usej ≤ dDEMANDj / PATTERNSjj e, usej ∈ IN
The paper mill can satisfy the demand with just the basic set of cutting patterns, but it is likely toincur significant losses through wasting more than necessary of every large roll and by cutting
more small rolls than its customers have ordered. We therefore employ a column generation
heuristic to find more suitable cutting patterns.
Our heuristic performs a column generation loop at the top node, before starting the MIP search.
Every iteration of the column generation loop executes the following steps:
1. solve the LP and save the basis

Column generation

c 2009 Fair Isaac Corporation. All rights reserved.

page 12Mosel: multiple models and parallel solving

2. get the solution values
3. compute a more profitable cutting pattern based on the current solution
4. generate a new column (= cutting pattern): add a term to the objective function and to thecorresponding demand constraints
5. load the modified problem and load the saved basis
Step 3 of this loop requires us to solve an integer knapsack problem of the formXmaximize z =Ci · xjj∈WIDTHSX

Aj · xj ≤ B

j∈WIDTHS∀j ∈ WIDTHS : xj integer
This second optimization problem is independent of the main, cutting stock problem since thetwo have no variables in common.

3.2

Implementation
The implementation is divided into two parts: the master model (file paperp.mos) with thedefinition of the cutting stock problem and the column generation algorithm, and the knapsackmodel (file knapsack.mos) that is run from the master.
3.2.1

Master model

The main part of the cutting stock model looks as follows:model "Papermill (multi-prob)"uses "mmxprs", "mmjobs"
forward procedure column_gen
forward function knapsack(C:array(range) of real,
A:array(range) of real,
B:real,xbest:array(range) of integer): realdeclarations
NWIDTHS = 5
WIDTHS = 1..NWIDTHS
RP: range
MAXWIDTH = 94
EPS = 1e-6

! Number of different widths! Range of widths! Range of cutting patterns! Maximum roll width! Zero tolerance

WIDTH: array(WIDTHS) of real
! Possible widths
DEMAND: array(WIDTHS) of integer! Demand per width
PATTERNS: array(WIDTHS, WIDTHS) of integer ! (Basic) cutting patterns

Column generation

use: array(RP) of mpvar
soluse: array(RP) of real
Dem: array(WIDTHS) of linctr
MinRolls: linctr! Rolls per pattern! Solution values for variables ‘use’! Demand constraints! Objective function

Knapsack: Model
end-declarations! Reference to the knapsack modelc 2009 Fair Isaac Corporation. All rights reserved.

page 13Mosel: multiple models and parallel solving

WIDTH:: [ 17, 21, 22.5, 24, 29.5]
DEMAND:: [150, 96,
48, 108, 227]

forall(j in WIDTHS)! Make basic patterns
PATTERNS(j,j) := floor(MAXWIDTH/WIDTH(j))

forall(j in WIDTHS) do
create(use(j))! Create NWIDTHS variables ‘use’
use(j) is_integer! Variables are integer and boundeduse(j) <= integer(ceil(DEMAND(j)/PATTERNS(j,j)))end-doMinRolls:= sum(j in WIDTHS) use(j)! Objective: minimize no. of rolls! Satisfy all demands
forall(i in WIDTHS)
Dem(i):= sum(j in WIDTHS) PATTERNS(i,j) * use(j) >= DEMAND(i)res:= compile("knapsack.mos")load(Knapsack, "knapsack.bim")
column_gen! Compile the knapsack model! Load the knapsack model! Column generation at top node

minimize(MinRolls)! Compute the best integer solution! for the current problem (including! the new columns)
writeln("Best integer solution: ", getobjval, " rolls")
write(" Rolls per pattern: ")forall(i in RP) write(getsol(use(i)),", ")writelnend-model

Before starting the column generation heuristic (the definition of procedure column_gen is left
out here since it remains unchanged from the User Guide example) the knapsack model iscompiled and loaded so that at every column generation loop we merely need to run it with new
data. The knapsack model is run from the function knapsack that takes as its parameters thedata for the knapsack problem and its solution values. The function saves all data to sharedmemory, then runs the knapsack model and retrieves the solution from shared memory. Its return
value is the objective value (zbest) of the knapsack problem.
function knapsack(C:array(range) of real,
A:array(range) of real,
B:real,xbest:array(range) of integer):real
initializations to "raw:noindex"
A as "shmem:A" B as "shmem:B" C as "shmem:C"
end-initializations
run(Knapsack, "NWIDTHS="+NWIDTHS)
wait
dropnextevent! Start solving knapsack subproblem! Wait until subproblem finishes! Ignore termination message

initializations from "raw:"
xbest as "shmem:xbest" returned as "shmem:zbest"
end-initializations
end-function

To enforce a sequential execution of the two models (we need to retrieve the results from theknapsack problem before we may continue with the master) we must add a call to the procedure
wait immediately after the run statement. Otherwise the execution of the master model
continues concurrently to the child model. On termination, the child model sends a ‘termination’
event (an event of class EVENT_END). Since our algorithm does not require this event we simply

Column generation

c 2009 Fair Isaac Corporation. All rights reserved.

page 14

Mosel: multiple models and parallel solving

remove it from the model’s event queue with a call to dropnextevent.
3.2.2

Knapsack model

The implementation of the knapsack model is straightforward. All problem data is obtained fromshared memory and after solving the problem its solution is saved into shared memory.
model "Knapsack"uses "mmxprs"
parameters
NWIDTHS=5
end-parameters! Number of different widths

declarations
WIDTHS = 1..NWIDTHS
A,C: array(WIDTHS) of real
B: real
KnapCtr, KnapObj: linctr
x: array(WIDTHS) of mpvar
xbest: array(WIDTHS) of integer
end-declarations! Range of widths! Constraint + obj. coefficients! RHS value of knapsack constraint! Knapsack constraint+objective! Knapsack variables! Solution values

initializations from "raw:noindex"
A as "mmjobs.shmem:A" B as "mmjobs.shmem:B" C as "mmjobs.shmem:C"
end-initializations! Define the knapsack problem
KnapCtr:= sum(j in WIDTHS) A(j)*x(j) <= B
KnapObj:= sum(j in WIDTHS) C(j)*x(j)forall(j in WIDTHS) x(j) is_integer! Solve the problem and retrieve the solution
maximize(KnapObj)
z:=getobjval
forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))
initializations to "raw:"
xbest as "mmjobs.shmem:xbest" z as "mmjobs.shmem:zbest"
end-initializations
end-model

In this model we have prefixed the shared memory driver shmem with the name of the module
mmjobs. Doing so is only required if we want to be able to run the knapsack model separately,that is, without the cutting stock master model that loads the mmjobs module into memory.
Another case where we would have to explicitly add the name of a module to a driver occurs
when we need to distinguish between several shmem drivers defined by different modules.

3.3

Results
With the data in the model above, the column generation algorithm generates 6 new patterns,
taking the value of the LP-relaxation of the cutting stock problem from originally 177.67 down to160.95. The MIP finds a solution with 161 rolls using the following patterns:

Column generation

c 2009 Fair Isaac Corporation. All rights reserved.

page 15

Mosel: multiple models and parallel solving

Widths4

Pattern

17

21

22.5

24

29.5

Usage

3

0

0

4

0

0

1

5

0

0

0

0

3

15

6

0

1

0

3

0

32

8

2

0

0

0

2

75

10

0

2

1

0

1

32

11

0

0

2

2

0

6

Solving several model instances in parallel
In this section we show how to execute several models in parallel and communicate solution
information among these models. This scheme may be particularly interesting when working
with Mosel on a multi-processor machine, e.g. by starting a number of models that corresponds
to the available number of processors.
Our idea is to run several instances (different only by the parameterization of the solution
algorithm) of the same MIP model concurrently and to stop the entire run when the first model
has finished. If the different solution algorithms are complementary in the sense that somequickly produce (good) solutions and others are better at proving optimality once the bestsolution is found then one may reasonably expect an additional synergy effect from exchanging
solution updates during the MIP search.
To implement this scheme, we define a master model that starts the model runs and coordinates
the solution updates, and a parameterizable child model that is loaded and run with the desired
number of versions. The child models all use the same solver (Xpress-Optimizer) but it would
equally be possible to use a different solver for some of the child models, provided it defines thenecessary functionality for interacting with the search.

4.1

Example problem: economic lot sizing
Economic lot sizing (ELS) considers production planning over a given planning horizon, in our
example a range of time periods TIMES = 1, . . . , T. In every period t, there is a given demand
DEMANDpt for every product p (p ∈ PRODUCTS) that must be satisfied by the production in thisperiod and by inventory carried over from previous periods.
A set-up cost SETUPCOSTt is associated with production in a period, and the total production
capacity per period, CAPt , is limited. The unit production cost PRODCOSTpt per product and time
period is also given. There is no inventory or stock-holding cost.
We introduce the decision variables producept for the amount of product p made in period t andthe binary variables setuppt indicating whether a setup takes place for product p in period t(setuppt = 1) or not (setuppt = 0).
We may then formulate the following mathematical model for this problem:


X
X
X
SETUPCOSTt ·minimizesetuppt +PRODCOSTpt · producept t∈TIMES∀p ∈ PRODUCTS, t ∈ TIMES :p∈PRODUCTSt
X
s=1

Solving several model instances in parallel

produceps ≥p∈PRODUCTSt
X

DEMANDps

s=1

c 2009 Fair Isaac Corporation. All rights reserved.

page 16Mosel: multiple models and parallel solving∀p ∈ PRODUCTS, t ∈ TIMES : producept ≤ DptT · setupptX∀t ∈ TIMES :producept ≤ CAPtp∈PRODUCTS∀p ∈ PRODUCTS, t ∈ TIMES : setuppt ∈ {0, 1}, producept ≥ 0
The objective function is to minimize the total cost. The constraints in the second line formulate
the requirement that the production of p in periods 0 to t must satisfy the total demand for thisproduct during this period of time. The next set of constraints establish the implication ‘if there is
production during t then there is a setup in t’ where Dptl stands for the demand of product p inperiods t to l. The production capacity per period t is limited. And finally, the setuppt variables
are binaries.
4.1.1

Cutting plane algorithm

A well-known class of valid inequalities for ELS are the so-called (l, S)-inequalities [PW94]. If Dptl
denotes the total demand of p in periods t to l, then for each period l and each subset of periods
S of 1 to l, the (l, S)-inequality islX

producept +

t=1
t∈S

l
X

Dptl · setuppt ≥ Dp1l

t=1
t6∈S

It says that actual production producept in the periods included in S plus the maximum potential
production Dptl · setuppt in the remaining periods (those not in S) must at least equal the total
demand in periods 1 to l.
It is possible to develop the following cutting plane algorithm based on these (l, S)-inequalities:
1. Solve the LP.
2. Identify violated (l, S)-inequalities by testing violations oflXmin(producept , Dptl · setuppt ) ≥ Dp1l

t=1

3. Add violated inequalities as cuts to the problem.
4. Re-solve the LP problem.
There are numerous options for how to configure this algorithm. For instance:• Generation of cuts only in the root node or also during the search (Cut-and-Branch versus
Branch-and-Cut).
• Number of cut generation passes at a node (e.g. one pass or looping around steps 2.-4. until
no more cuts are generated).
• Search tree depth for cut generation (up to a given depth or at all nodes).
• Exclusive use of (l, S)-cuts or combination with others (e.g. default cuts generated by thesolver).
The implementation of the (l, S)-cut generation algorithm shown below may be configured togenerate cuts at the top node only (TOPONLY = true) and to generate one or several rounds ofcuts (SEVERALROUNDS = true).

Solving several model instances in parallelc 2009 Fair Isaac Corporation. All rights reserved.

page 17Mosel: multiple models and parallel solving

4.2

Implementation
With mmjobs events are always sent between parent – child pairs, not directly from one child toanother. The ’solution found’ message therefore needs to be sent to the parent model that thenpasses on this message to all other child models.
Another point that should be stressed is the fact that we only compile the ELS model file once,but the number of instances loaded into memory needs to correspond to the number of childmodels we wish to run.
4.2.1

Master model

The master model compiles, loads and runs the child models and coordinates the solution
updates. Some care must be taken with the solution updates since new solutions that arereported are not guaranteed to be better than others previously reported by other child models.
For instance, if two models find solutions almost at the same time, the first solution that reaches
the master may be the better one and it must not be overridden by the next.
For a nice solution display at the end, the master model also reads in parts of the problem datafrom file.
model "Els master"uses "mmjobs"
parameters
DATAFILE = "els5.dat"
T = 45
P = 4
end-parameters
declarations
RM = 0..5! Range of models
TIMES = 1..T
! Time periods
PRODUCTS = 1..P! Set of products
solprod: array(PRODUCTS,TIMES) of real
! Sol. values for var.s produce
solsetup: array(TIMES) of real
! Sol. values for var.s setup
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
modELS: array(RM) of Modelmodid: array(set of integer) of integer
NEWSOL = 2
Msg: Event
end-declarations! Models! Model indices! Identifier for "sol. found" event! Messages sent by models! Compile, load, and run models M1 and M2
M1:= 1; M2:=2
res:= compile("elsp.mos")load(modELS(M1), "elsp.bim")load(modELS(M2), "elsp.bim")
forall(m in RM) modid(getid(modELS(m))):= mrun(modELS(M1), "ALG="+M1+",DATAFILE="+DATAFILE+",T="+T+",P="+P)
run(modELS(M2), "ALG="+M2+",DATAFILE="+DATAFILE+",T="+T+",P="+P)
objval:= MAX_REAL
algsol:= -1; algopt:= -1
repeatwait
Msg:= getnextevent
if getclass(Msg)=NEWSOL thenif getvalue(Msg) < objval thenalgsol:= modid(getfromid(Msg))
objval:= getvalue(Msg)

Solving several model instances in parallel! Wait for the next event! Get the event! Get the event class! Value of the event (= obj. value)! ID of model sending the eventc 2009 Fair Isaac Corporation. All rights reserved.

page 18

Mosel: multiple models and parallel solving

writeln("Improved solution ", objval, " found by model ", algsol)
forall(m in RM | m <> algsol) send(modELS(m), NEWSOL, objval)
else
writeln("Solution ", getvalue(Msg), " found by model ",modid(getfromid(Msg)))
end-ifend-ifuntil getclass(Msg)=EVENT_END! A model has finishedalgopt:= modid(getfromid(Msg))
forall(m in RM) stop(modELS(m))! Retrieve ID of terminated model! Stop all running models! Retrieve the best solution from shared memory
initializations from "raw:noindex"
solprod as "shmem:solprod"+algsol
solsetup as "shmem:solsetup"+algsol
end-initializations
initializations from DATAFILE
DEMAND
end-initializations! Solution printing
writeln("Best solution found by model ", algsol)
writeln("Optimality proven by model ", algopt)
writeln("Objective value: ", objval)
write("Period setup")
forall(p in PRODUCTS) write(strfmt(p,-7))
forall(t in TIMES) do
write("\n ", strfmt(t,2), strfmt(solsetup(t),8), "")forall(p in PRODUCTS) write(strfmt(solprod(p,t),3), " (",DEMAND(p,t),")")
end-dowritelnend-model

In this implementation we define an array modid that establishes the correspondence between
the model index used in our model and Mosel’s internal ID of the model. Whenever a child model
sends an event to the master, we retrieve its ID (with function getfromid) and store thecorresponding model index, to be able to use it for solution printing later on.
4.2.2

ELS model

The ELS child model is written in such a way that the model can be executed separately. In
particular, every model performs the complete initialization of its data from file, a task that for
greater efficiency could be reserved to the master model, communicating data via sharedmemory to the child models (however, in our example data handling time is negligible comparedto the running time of the solution algorithms).
The main part of the ELS model contains the definition of the model itself and the selection of
the solution algorithm:
model Els
uses "mmxprs","mmjobs"
parameters
ALG = 0
DATAFILE = "els4.dat"
T = 60
P = 4
end-parameters

Solving several model instances in parallel! Default algorithm: no user cuts

c 2009 Fair Isaac Corporation. All rights reserved.

page 19

Mosel: multiple models and parallel solving

forward procedure tree_cut_gen
forward public function cb_node: boolean
forward public function cb_updatebnd(node:integer): integer
forward public procedure cb_intsol
declarations
NEWSOL = 2
EPS = 1e-6
TIMES = 1..TPRODUCTS = 1..P! "New solution" event class! Zero tolerance! Time periods! Set of products

DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
SETUPCOST: array(TIMES) of integer! Setup cost per period
PRODCOST: array(PRODUCTS,TIMES) of real ! Production cost per period
CAP: array(TIMES) of integer! Production capacity per period
D: array(PRODUCTS,TIMES,TIMES) of integer ! Total demand in periods t1 - t2
produce: array(PRODUCTS,TIMES) of mpvar
setup: array(TIMES) of mpvar! Production in period t! Setup in period t

solprod: array(PRODUCTS,TIMES) of real
solsetup: array(TIMES) of real

! Sol. values for var.s produce! Sol. values for var.s setup

Msg: Event
end-declarations! An event

initializations from DATAFILE
DEMAND SETUPCOST PRODCOST CAP
end-initializations
forall(p in PRODUCTS,s,t in TIMES) D(p,s,t):= sum(k in s..t) DEMAND(p,k)! Objective: minimize total cost
MinCost:= sum(t in TIMES) (SETUPCOST(t) * setup(t) +sum(p in PRODUCTS) PRODCOST(p,t) * produce(p,t) )! Production in period t must not exceed the total demand for the! remaining periods; if there is production during t then there! is a setup in t
forall(t in TIMES)sum(p in PRODUCTS) produce(p,t) <= sum(p in PRODUCTS) D(p,t,T) * setup(t)! Production in periods 0 to t must satisfy the total demand! during this period of time
forall(p in PRODUCTS,t in TIMES)sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)! Capacity limits
forall(t in TIMES) sum(p in PRODUCTS) produce(p,t) <= CAP(t)forall(t in TIMES) setup(t) is_binary! Variables setup are 0/1

setparam("zerotol", EPS/100)
SEVERALROUNDS:=false; TOPONLY:=false! Set Mosel comparison tolerance

case ALG of
1: do
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_HEURSTRATEGY", 0)
end-do2: do
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_HEURSTRATEGY", 0)
setparam("XPRS_PRESOLVE", 0)
end-do3: tree_cut_gen
4: doSolving several model instances in parallel! No cuts! No heuristics! No cuts! No heuristics! No presolve! User branch&cut (single round)! User branch&cut (several rounds)

c 2009 Fair Isaac Corporation. All rights reserved.

page 20

Mosel: multiple models and parallel solving

tree_cut_gen
SEVERALROUNDS:=true
end-do5: dotree_cut_gen
SEVERALROUNDS:=true
TOPONLY:=true
end-doend-case! User cut&branch (several rounds)! Parallel setup
setcallback(XPRS_CB_PRENODE, "cb_updatebnd") ! Node pre-treatment callback
setcallback(XPRS_CB_INTSOL, "cb_intsol") ! Integer solution callback
setparam("XPRS_SOLUTIONFILE",0)! Do not save solutions to file! Solve the problem
minimize(MinCost)
end-model

The procedure tree_cut_gen sets up a user cut generation routine, configurable to generatecuts only at the top node of the branch-and-bound search (TOPONLY) or to execute one or several
cut generation iterations per node (SEVERALROUNDS). The definition of the cut generation
routine cb_node itself is left out here.
procedure tree_cut_gen
setparam("XPRS_HEURSTRATEGY", 0)
setparam("XPRS_CUTSTRATEGY", 0)
setparam("XPRS_PRESOLVE", 0)
setparam("XPRS_EXTRAROWS", 5000)
setcallback(XPRS_CB_CUTMGR, "cb_node")
end-procedure! Switch heuristics off! Switch automatic cuts off! Switch presolve off! Reserve extra rows in matrix! Define the cut manager callback

The communication between concurrently running child models has two parts: (a) any integer
solution found must be saved and comunicated to the master model and (b) bound updates sentby the master problem must be incorporated into the search. Xpress-Optimizer provides a specificinteger solution callback for saving solutions into user structures. An obvious place for boundupdates in nodes is the cut-manager callback function. However, this function being already inuse for other purposes with certain settings of the algorithm, we employ a different callback
function that also gets called at every node, the node pre-treatment callback.
! Update cutoff value
public function cb_updatebnd(node:integer): integer
if not isqueueempty thenrepeatMsg:= getnextevent
until isqueueempty
newcutoff:= getvalue(Msg)
setparam("XPRS_MIPABSCUTOFF", newcutoff)
if (newcutoff < getparam("XPRS_LPOBJVAL")) thenreturned:= 1! Node becomes infeasible
end-ifend-ifend-function! Store and communicate new solution
public procedure cb_intsol
objval:= getparam("XPRS_LPOBJVAL")
cutoff:= getparam("XPRS_MIPABSCUTOFF")if(cutoff > objval) thensetparam("XPRS_MIPABSCUTOFF", objval)
end-ifSolving several model instances in parallel! Retrieve current objective valuec 2009 Fair Isaac Corporation. All rights reserved.

page 21

Mosel: multiple models and parallel solving! Get the solution values
forall(t in TIMES) do
forall(p in PRODUCTS) solprod(p,t):=getsol(produce(p,t))solsetup(t):=getsol(setup(t))
end-do! Store the solution in shared memory
initializations to "raw:noindex"
solprod as "shmem:solprod"+ALG
solsetup as "shmem:solsetup"+ALG
end-initializations! Send "solution found" signal
send(NEWSOL, objval)
end-procedure

The bound update callback function checks whether the event queue contains any events, if thisis the case, it takes all events from the queue and sets the value of the last event as the new
cutoff value. The rationale behind the loop for emptying the event queue is that the mastermodel may have sent several improved solution values since the last check, the best value isalways the one sent last, that is, the last in the queue.
The integer solution callback writes the solution values to shared memory, adding the identifierof the model (= value of ALG). The latter ensures that two child models that possibly write out
their solution at the same time do not use the same memory area.

4.3

Results
A run with two models may generate a log similar to the following one (note that the model thatterminates the search is not the same that has found the optimal solution).
Improved solution 1283 found by model 2
Improved solution 1250 found by model 2
Improved solution 1242 found by model 1
Improved solution 1236 found by model 2
Improved solution 1234 found by model 2
Best solution found by model 2
Optimality proven by model 1
Objective value: 1234

5Dantzig-Wolfe decomposition: combining sequential and parallel solving
Dantzig-Wolfe decomposition (see [Teb01] for further detail) is a solution method for problems
where, if a relatively small number of constraints were removed, the problem would fall apartinto a number of independent problems. This means, it is possible to re-order the rows andcolumns of the constraint matrix as shown in Figure 1, where non-zero coefficients only occur
within the gray shaded areas. Such a primal block angular structure may become immediately
apparent by visualizing a problem matrix with Xpress-IVE. However, in most cases it will benecessary to re-organize the constraint definitions, grouping them by common index (sub)sets
such as time periods, products, plant locations, and so on.
The constraints (including the objective function) involving variables of several or all subproblems
are referred to as global constraints (also: common, linking, or master constraints). These
constraints are used to form the master problem. The individual blocks on the diagonal of thecoefficient matrix are solved as pricing subproblems, coordinated by the master problem. By

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 22

Mosel: multiple models and parallel solving..

.

Figure 1: Coeffcient matrix with primal block angular structuresolving the master problem we obtain a solution to the original problem. Since the masterproblem has a large number of variables (defined by the set of basic feasible solutions andunbounded directions of the pricing problems), we work with a restricted master problem over asmall subset of the variables. The variables to enter the active set of variables of the restricted
master problem are determined by solving the pricing subproblems. With objective functions forthe pricing problems that are based on the dual values of the restricted master problem we canobtain that the objective function value at each extreme point is the reduced cost (or price) of
the master problem variable associated with the extreme point.
For maximization problems, solving the modified pricing problems generates basic feasible
solutions of maximum reduced cost. If the objective value at an extreme point is positive, thenthe associated master problem variable is added to the master problem; if the minimum objective
value over all extreme points is negative, then no master problem variables exists to improve thecurrent master problem solution.
The computational advantage of Dantzig-Wolfe decomposition arises from performing asignificant amount of the computational work in the pricing problems that are roughly an orderof magnitude smaller than the original problem and thus easier to solve. An aspect of themethod that is of interest in the context of this paper is that the subproblems are independent ofeach other, and may therefore be solved concurrently.
A potential drawback of the decomposition approach is the huge size of the master problem, it
has many more variables—though fewer constraints—than the original problem. In general it is
not required to generate all variables explicitly but since the feasible region of the masterproblem is more complex than that of the original problem, the solution path may be longer.
Furthermore, numerical problems may occur through the dynamic generation of variables of themaster problem.
Many factors may influence the performance of a decomposition approach, so for a particularapplication computational experiments will be required to find out whether this solution method
is suitable. Such tests may include different ways of decomposing a given problem. As a generalrule for the definition of a problem decomposition, one should aim for few global constraints
since it is important to be able to (re)solve the master problem quickly. In addition, the pricing
problems should be constructed in such a way that they are well formed problems in their own
right to avoid computational problems due to degeneracy.

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 23

Mosel: multiple models and parallel solving

5.1

Example problem: multi-item, multi-period production planning
The company Coco has two plants that can produce two types of cocoa powder. The first factory
has a total capacity of 400 tons per month and the second of 500 tons per month. The marketing
department has provided estimations for the maximum amount of every product that may besold in each of the next four months, and also the expected sales prices. Several raw materials areused in the production with known raw material requirements per ton of finished products.
Finished products and raw material may be stored at the factories from one tim period to thenext, incurring a given cost per ton and per time period. The raw material storage capacity at thefactories is limited to 300 tons. How should the two plants be operated during the planning
period to maximize the total profit?
5.1.1

Original model

Let PRODS be the set of finished products, FACT the set of factories, RAW the set of raw
materials, and TIMES = {1, . . . , NT} the set of time periods under consideration.
We define decision variables makepft representing the quantity of product p made at factory fduring time period t. Furthermore, to model the transition from one time period to the next andto account for the different types of cost incurred, we need several other sets of variables: sellpft ,the amount of product p sold at factory f in period t; buyrft the amount of raw material r bought
at f in t; and finally pstockpft and rstockrft (both defined for t = 1, . . . , NT + 1) respectively theamount of product and raw material held in stock at factory f at the beginning of period t.
Let further MXSELLpt be the maximum sales quantity of product p in period t, MXMAKEf thecapacity limit of factory f, and MXRSTOCK the raw material storage capacity.
Let IPSTOCKpf be the quantity of product p initially held in stock at factory f and IRSTOCKrf theinitial stock level of raw material r. We denote by CPSTOCK and CRSTOCK respectively the unit
cost for storing finished product and raw material.
The objective function of maximizing the total profit is to maximize the sales revenues (REVp ),minus the cost of production (CMAKEpf ), buying raw material (CBUYrt ), and storing finishedproducts and raw material.

Xmaximize

X

X

REVpt · sellpftp∈PRODS f∈FACT t∈TIME

−

X

X

X

CMAKEpf · makepft −p∈PRODS f∈FACT t∈TIME−X

X NT+1X

X

X

CBUYrt · buyrftr∈RAW f∈FACT t∈TIME

CPSTOCK · pstockpft −p∈PRODS f∈FACT t=2

XX

X NT+1X

CRSTOCK · rstockrftr∈RAW f∈FACT t=2

The possibility to store products between time periods gives rise to three sets of constraints: theinventory balance constraints for finished products (PBalpft ) and for raw material (RBalrft ), and thelimit on the raw material storage capacity (MxRStockft ). The stock pstockp,f,t+1 of product p at the
beginning of t + 1 is given by the stock at the beginning of t plus the production in t reduced by
the amount sold on t. The stock rstockr,f,t+1 of raw material r at the beginning of t + 1 is given by
the corresponding stock at the beginning of t plus the amount bought in t reduced by thequantity used in production during t.
∀p ∈ PRODS, ∀f ∈ FACT, ∀t ∈ TIME : PBalpft := pstockp,f,t+1 = pstockpft + makepft − sellpft∀r ∈ RAW, ∀f ∈ FACT, ∀t ∈ TIME :Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 24

Mosel: multiple models and parallel solving

X

RBalrft := rstockr,f,t+1 = rstockrft + buyrft −REQpr · makepftp∈PRODSX∀f ∈ FACT, ∀t ∈ {2, . . . , NT + 1} : MxRStockft :=rstockrft ≤ MXRSTOCKr∈RAW

We further have two sets of capacity constraints: the production capacity of the factories islimited (constraints MxMakeft ) and we can only sell up to a given maximum amount of finishedproducts per time period (constraints MxSellft ).
Below the complete mathematical model is given. The initial stock levels (t = 1) of finishedproducts and raw material are fixed to the given values.

X

maximize

X

X

REVpt · sellpftp∈PRODS f∈FACT t∈TIME

−

X

X

X

CMAKEpf · makepft −p∈PRODS f∈FACT t∈TIME−X

X NT+1X

X

X

X

CBUYrt · buyrftr∈RAW f∈FACT t∈TIME

CPSTOCK · pstockpft −X NT+1X

X

CRSTOCK · rstockrftr∈RAW f∈FACT t=2

p∈PRODS f∈FACT t=2∀p ∈ PRODS, ∀f ∈ FACT, ∀t ∈ TIME : PBalpft := pstockp,f,t+1 = pstockpft + makepft − sellpft∀r ∈ RAW, ∀f ∈ FACT, ∀t ∈ TIME :X

RBalrft := rstockr,f,t+1 = rstockrft + buyrft −REQpr · makepftp∈PRODS∀p ∈ PRODS, ∀t ∈ TIME : MxSellpt :=X

sellpft ≤ MXSELLp

f∈FACT∀f ∈ FACT, ∀t ∈ TIME : MxMakeft :=X

makepft ≤ MXMAKEfp∈PRODSX∀f ∈ FACT, ∀t ∈ {2, . . . , NT + 1} : MxRStockft :=rstockrft ≤ MXRSTOCKr∈RAW∀p ∈ PRODS, ∀f ∈ FACT : pstockpf1 = IPSTOCKpf∀r ∈ RAW, ∀f ∈ FACT : rstockrf1 = IRSTOCKrf∀p ∈ PRODS, ∀f ∈ FACT, ∀t ∈ TIME : makepft ≥ 0, sellpft ≥ 0∀r ∈ RAW, ∀f ∈ FACT, ∀t ∈ TIME : buyrft ≥ 0∀p ∈ PRODS, ∀f ∈ FACT, ∀t ∈ {1, . . . , NT + 1} : pstockpft ≥ 0∀r ∈ RAW, ∀f ∈ FACT, ∀t ∈ {1, . . . , NT + 1} : rstockrft ≥ 0
5.1.2

Problem decomposition

We now decompose the problem described above according to production locations. Notice that
this is not the only way of decomposing this problem: we may just as well choose adecomposition by products or by time periods. However, both of these choices result in a larger
number of global constraints than the decomposition by factories, meaning that the masterproblem may become more difficult to solve.
For every factory f we obtain the following subproblem (including the sales limit constraints
MxSell in the form of bounds in the submodels is not required, but may lead to better, that is,more exact or faster solutions).

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 25

Mosel: multiple models and parallel solvingX

maximize

X

REVpt · sellptp∈PRODS t∈TIME

X

−

X

p∈PRODS t∈TIME
NT+1
X

X

−

X

CMAKEp · makept −X

CBUYrt · buyrtr∈RAW t∈TIME

X NT+1X

CPSTOCK · pstockpt −CRSTOCK · rstockrtr∈RAW t=2p∈PRODS t=2∀p ∈ PRODS, ∀t ∈ TIME : PBalpt := pstockp,t+1 = pstockpt + makept − sellptX∀r ∈ RAW, ∀t ∈ TIME : RBalrt := rstockr,t+1 = rstockrt + buyrt −REQpr · makeptp∈PRODSX∀t ∈ TIME : MxMaket :=makept ≤ MXMAKEp∈PRODS∀t ∈ {2, . . . , NT + 1} : MxRStockt :=X

rstockrt ≤ MXRSTOCKr∈RAW∀p ∈ PRODS, ∀t ∈ TIME : MxSellpt := sellpt ≤ MXSELLp∀p ∈ PRODS : pstockp1 = IPSTOCKp∀r ∈ RAW : rstockr1 = IRSTOCKr∀p ∈ PRODS, ∀t ∈ TIME : makept ≥ 0, sellpt ≥ 0∀r ∈ RAW, ∀t ∈ TIME : buyrt ≥ 0∀p ∈ PRODS, ∀t ∈ {1, . . . , NT + 1} : pstockpt ≥ 0∀r ∈ RAW, ∀t ∈ {1, . . . , NT + 1} : rstockrt ≥ 0
The master problem only contains a single set of global constraints, namely the sales limit
constraints MxSell (for clarity’s sake, the non-negativity constraints are left out here).

Xmaximize

X

X

REVpt · sellpftp∈PRODS f∈FACT t∈TIME

−

X

X

X

CMAKEpf · makepft −p∈PRODS f∈FACT t∈TIME−X

X NT+1X

X

X

CBUYrt · buyrftr∈RAW f∈FACT t∈TIME

CPSTOCK · pstockpft −X

X NT+1X

CRSTOCK · rstockrftr∈RAW f∈FACT t=2

p∈PRODS f∈FACT t=2∀p ∈ PRODS, ∀t ∈ TIME : MxSellpt :=XX

sellpft ≤ MXSELLp

f∈FACT

In the decomposition algorithm, the decision variables in the master problem are expressed via
the solutions (proposals) generated by the subproblems, such as
nPROPf∀p ∈ PRODS, ∀f ∈ FACT, ∀t ∈ TIME : sellpft =X

Prop_sellpftk · weightfk

k=1

where Prop_sellpftk stands for the solution value of variable sellpt in the kth proposal generated bysubproblem f and Prop_costfk is the objective value this proposal. For every subproblem f weneed to add a convexity constraint Convexf on the weightfk variables.

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 26

Mosel: multiple models and parallel solving

nPROPf

X

X

f∈FACT

k=1

maximize

Prop_costfk · weightfk
nPROPf∀p ∈ PRODS, ∀t ∈ TIME : MxSellpt :=X

X

f∈FACT

k=1

Prop_sellpftk · weightfk ≤ MXSELLp

nPROPf∀f ∈ FACT : Convexf :=Xweightfk = 1

k=1∀f ∈ FACT, ∀k ∈ 1, . . . , nPROPf : weightfk ≥ 0
We shall refer to this problem as the modified master problem. Without going any further into
technical detail we simply remark that a correspondence exists between the solution of theoriginal problem and those of the modified master problem.

5.2

Implementation
The decomposition algorithm has several phases:• Phase 1: generation of a set of proposals to obtain a feasible solution to the modified
master problem.
• Phase 2: optimization of the modified master problem.
• Phase 3: calculate the solution to the original problem.
The suproblems solved in phases 1 and 2 take as their objective functions sums of the dual values
from the modfied master problem. To avoid starting off with an empty master problem andhence random dual information that may be misleading we add a crash Phase 0 to ourimplementation that generates one proposal for each subproblem with the original objectivefunction.
5.2.1

Master model

Below follows the body of the master model file. The definitions of the function bodies will beshown later in Section 5.2.3.
model "Coco3 Master"uses "mmxprs","mmjobs","mmsystem"
parameters
DATAFILE = "coco2.dat"
ALG = 0! 0: stop phase with 1st failed subpb.
! 1: stop when all subprob.s failend-parameters
forward procedure process_sub_result
forward procedure solve_master(phase:integer)
forward procedure process_master_result
forward function calc_solution:real
forward procedure print_solution
declarations
PHASE_0=2
PHASE_1=3Dantzig-Wolfe decomposition! Event codes sent to submodels

c 2009 Fair Isaac Corporation. All rights reserved.

page 27

Mosel: multiple models and parallel solving

PHASE_2=4
PHASE_3=5
EVENT_SOLVED=6
EVENT_FAILED=7
EVENT_READY=8
NPROD, NFACT, NRAW, NT: integer
end-declarations! Event codes sent by submodels

initializations from DATAFILE
NPROD NFACT NRAW NT
end-initializations
declarations
PRODS = 1..NPRODFACT = 1..NFACTRAW = 1..NRAWTIME = 1..NTnIter: integer
nPROP: array(FACT) of integer
end-declarations! Range of products (p)!
factories (f)!
raw materials (r)!
time periods (t)! Iteration counter! Counters of proposals from subprob.s!**** Master problem ****
declarations
MXSELL: array(PRODS,TIME) of real
! Max. amount of p that can be soldexcessS: mpvar! Violation of sales/buying limits
weight: array(FACT,range) of mpvar ! weights for proposals
MxSell: array(PRODS,TIME) of linctr ! Sales limit constraints
Convex: array(FACT) of linctr! Convexity constraints
Price_convex: array(FACT) of real
! Dual price on convexity constraints
Price_sell: array(PRODS,TIME) of real ! Dual price on sales limits
end-declarations
initializations from DATAFILEMXSELL
end-initializations!**** Submodels ****
declarations
submod: array(FACT) of Model! One subproblem per factory
Stopped: set of integer
modid: array(set of integer) of integer ! Model indicesend-declarations
res:= compile("g","cocoSubF.mos")! Compile the submodel file
forall(f in FACT) do! Load & run one submodel per product
Price_convex(f):= 1
load(submod(f), "cocoSubF.bim")
modid(getid(submod(f))):= frun(submod(f), "Factory=" + f + ",DATAFILE=" + DATAFILE)wait! Wait for child model to be readydropnextevent
end-do!**** Phase 0: Crash ****
nIter:=1; finished:=false
writeln("\nPHASE 0 -- Iteration ", nIter); fflush
forall(f in FACT)
send(submod(f), PHASE_0, 0)! Start solving all submodels (Phase 1)

forall(f in FACT) do
wait! Wait for child (termination) events
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result! Add new proposal to master problem
elif getclass(ev)=EVENT_FAILED thenDantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 28

Mosel: multiple models and parallel solving

finished:= true
end-ifend-doif finished thenwriteln("Problem is infeasible")exit(1)
end-ifsolve_master(1)process_master_result! Solve the updated Ph. 1 master problem! Store initial pricing data for submodels!**** Phase 1: proposal generation (feasibility) ****
repeatnoimprove:= 0
nIter+=1
writeln("\nPHASE 1 -- Iteration ", nIter); fflush
forall(f in FACT)! Start solving all submodels (Phase 1)
send(submod(f), PHASE_1, Price_convex(f))
forall(f in FACT) do
wait! Wait for child (termination) events
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result! Add new proposal to master problem
elif getclass(ev)=EVENT_FAILED thennoimprove += 1
end-ifend-doif noimprove = NFACT thenwriteln("Problem is infeasible")exit(2)
end-ifif ALG=0 and noimprove > 0 thenwriteln("No improvement by some subproblem(s)")
break
end-ifsolve_master(1)if getobjval>0.00001 thenprocess_master_result
end-ifuntil getobjval <= 0.00001! Solve the updated Ph. 1 master problem! Store new pricing data for submodels!**** Phase 2: proposal generation (optimization) ****
writeln("\n**** PHASE 2 ****")
finished:=false
repeatsolve_master(2)! Solve Phase 2 master problem
process_master_result! Store new pricing data for submodels
nIter+=1
writeln("\nPHASE 2 -- Iteration ", nIter); fflush
forall(f in FACT)! Start solving all submodels (Phase 2)
send(submod(f), PHASE_2, Price_convex(f))
forall(f in FACT) do
wait! Wait for child (termination) events
ev:= getnextevent
if getclass(ev)=EVENT_SOLVED thenprocess_sub_result! Add new proposal to master problem
elif getclass(ev)=EVENT_FAILED thenDantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 29

Mosel: multiple models and parallel solving

if ALG=0 thenfinished:=true! 1st submodel w/o prop. stops phase 2
else
Stopped += {modid(getfromid(ev))} ! Stop phase 2 only when no submodel! generates a new proposal
end-ifend-ifend-doif getsize(Stopped) = NFACT then finished:= true; end-if
until finishedsolve_master(2)! Re-solve master to integrate! proposal(s) from last ph. 2 iteration!**** Phase 3: solution to the original problem ****
writeln("\n**** PHASE 3 ****")
forall(f in FACT) dosend(submod(f), PHASE_3, 0)! Stop all submodels
wait
dropnextevent
end-dowriteln("Total Profit=", calc_solution)
print_solution
end-model

The initial declarations block of this model defines a certain number of event codes that areused to identify the messages sent between this master model and its child (sub)models. The
same declarations need to be repeated in the child models.
5.2.2

Single factory model

The model file cocoSubF.mos with the definition of the subproblems has the following contents.
model "Coco Subproblem (factory based decomp.)"uses "mmxprs", "mmjobs"
parameters
Factory = 0
TOL = 0.00001
DATAFILE = "coco3.dat"
end-parameters
forward procedure process_solution
declarations
PHASE_0=2
PHASE_1=3
PHASE_2=4
PHASE_3=5
EVENT_SOLVED=6
EVENT_FAILED=7
EVENT_READY=8
NPROD, NFACT, NRAW, NT: integer
end-declarations
send(EVENT_READY,0)! Event codes sent to submodels! Event codes sent by submodels! Model is ready (= running)

initializations from DATAFILE
NPROD NFACT NRAW NT
end-initializations
declarations

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 30

Mosel: multiple models and parallel solving

PRODS = 1..NPRODFACT = 1..NFACTRAW = 1..NRAWTIME = 1..NT! Range of products (p)!
factories (f)!
raw materials (r)!
time periods (t)REV: array(PRODS,TIME) of real
CMAKE: array(PRODS,FACT) of real

! Unit selling price of products! Unit cost to make product p! at factory f! Unit cost to buy raw materials! Requirement by unit of product p! for raw material r! Max. amount of p that can be sold! Max. amount factory f can make! over all products! Initial product stock levels! Initial raw material stock levels! Unit cost to store any product p! Unit cost to store any raw mat. r! Raw material storage capacity

CBUY: array(RAW,TIME) of real
REQ: array(PRODS,RAW) of real
MXSELL: array(PRODS,TIME) of real
MXMAKE: array(FACT) of real
IPSTOCK: array(PRODS,FACT) of real
IRSTOCK: array(RAW,FACT) of real
CPSTOCK: real
CRSTOCK: real
MXRSTOCK: real

make: array(PRODS,TIME) of mpvar! Amount of products made at factorysell: array(PRODS,TIME) of mpvar! Amount of product sold from factorybuy: array(RAW,TIME) of mpvar! Amount of raw material bought
pstock: array(PRODS,1..NT+1) of mpvar ! Product stock levels at start! of period t
rstock: array(RAW,1..NT+1) of mpvar! Raw material stock levels! at start of period t
sol_make: array(PRODS,TIME) of real
! Amount of products madesol_sell: array(PRODS,TIME) of real
! Amount of product sold
sol_buy: array(RAW,TIME) of real
! Amount of raw mat. bought
sol_pstock: array(PRODS,1..NT+1) of real ! Product stock levels
sol_rstock: array(RAW,1..NT+1) of real
! Raw mat. stock levels
Profit: linctr! Profit of proposal
Price_sell: array(PRODS,TIME) of real ! Dual price on sales limits
end-declarations
initializations from DATAFILECMAKE REV CBUY REQ MXSELL MXMAKEIPSTOCK IRSTOCK MXRSTOCK CPSTOCK CRSTOCK
end-initializations! Product stock balance
forall(p in PRODS,t in TIME)
PBal(p,t):= pstock(p,t+1) = pstock(p,t) + make(p,t) - sell(p,t)! Raw material stock balance
forall(r in RAW,t in TIME)
RBal(r,t):= rstock(r,t+1) =rstock(r,t) + buy(r,t) - sum(p in PRODS) REQ(p,r)*make(p,t)! Capacity limit
forall(t in TIME)
MxMake(t):= sum(p in PRODS) make(p,t) <= MXMAKE(Factory)! Limit on the amount of prod. p to be soldforall(p in PRODS,t in TIME) sell(p,t) <= MXSELL(p,t)! Raw material stock limit
forall(t in 2..NT+1)
MxRStock(t):= sum(r in RAW) rstock(r,t) <= MXRSTOCK! Initial product and raw material stock levels
forall(p in PRODS) pstock(p,1) = IPSTOCK(p,Factory)
forall(r in RAW) rstock(r,1) = IRSTOCK(r,Factory)Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 31

Mosel: multiple models and parallel solving! Total profit
Profit:=sum(p in PRODS,t in TIME) REV(p,t) * sell(p,t) sum(p in PRODS,t in TIME) CMAKE(p,Factory) * make(p,t) sum(r in RAW,t in TIME) CBUY(r,t) * buy(r,t) sum(p in PRODS,t in 2..NT+1) CPSTOCK * pstock(p,t) sum(r in RAW,t in 2..NT+1) CRSTOCK * rstock(r,t)! revenue! prod. cost
! raw mat.
! p storage! r storage! (Re)solve this model until it is stopped by event "PHASE_3"
repeatwait
ev:= getnextevent
Phase:= getclass(ev)if Phase=PHASE_3 then! Stop the execution of this model
break
end-ifPrice_convex:= getvalue(ev)! Get new pricing data
if Phase<>PHASE_0 theninitializations from "raw:noindex"
Price_sell as "shmem:Price_sell"
end-initializations
end-if! (Re)solve this model
if Phase=PHASE_0 thenmaximize(Profit)
elif Phase=PHASE_1 thenmaximize(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) + Price_convex)
else! PHASE 2
maximize(
Profit - sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) Price_convex)
end-ifwriteln("Factory ", Factory, " - Obj: ", getobjval," Profit: ", getsol(Profit), " Price_sell: ",getsol(sum(p in PRODS,t in TIME) Price_sell(p,t)*sell(p,t) )," Price_convex: ", Price_convex)
fflushif getobjval > TOL thenprocess_solution
elif getobjval <= TOL thensend(EVENT_FAILED,0)
else
send(EVENT_READY,0)
end-ifuntil false! Solution found: send values to master! Problem is infeasible (Phase 0/1) or! no improved solution found (Phase 2)!----------------------------------------------------------! Process solution dataprocedure process_solution
forall(p in PRODS,t in TIME) do
sol_make(p,t):= getsol(make(p,t))sol_sell(p,t):= getsol(sell(p,t))
end-doforall(r in RAW,t in TIME) sol_buy(r,t):= getsol(buy(r,t))
forall(p in PRODS,t in 1..NT+1) sol_pstock(p,t):= getsol(pstock(p,t))
forall(r in RAW,t in 1..NT+1) sol_rstock(r,t):= getsol(rstock(r,t))
Prop_cost:= getsol(Profit)
send(EVENT_SOLVED,0)
initializations to "mempipe:noindex,sol"
Factorysol_make sol_sell sol_buy sol_pstock sol_rstock
Prop_cost

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 32

Mosel: multiple models and parallel solving

end-initializations
end-procedure
end-model

The child models are re-solved until they receive the PHASE_3 code. At every iteration they write
their solution values to memory so that these can be processed by the master model.
5.2.3

Master problem subroutines

The following three subroutines of the master model recover the solutions produced by the
subproblems (process_sub_result), re-solve the master problem (solve_master), andcommunicate the solution of the master problem to its child models (process_master_result).
declarations
Prop_make: array(PRODS,FACT,TIME,range) of real ! Amount of products made
Prop_sell: array(PRODS,FACT,TIME,range) of real ! Amount of product sold
Prop_buy: array(RAW,FACT,TIME,range) of real
! Amount of raw mat. bought
Prop_pstock: array(PRODS,FACT,1..NT+1,range) of real ! Product stock levels
Prop_rstock: array(RAW,FACT,1..NT+1,range) of real
! Raw mat. stock levels
Prop_cost: array(FACT,range) of real ! Cost/profit of each proposal
end-declarations
procedure process_sub_result
declarations
f: integer! Factory index! Solution values of the proposal:
sol_make: array(PRODS,TIME) of real
! Amount of products madesol_sell: array(PRODS,TIME) of real
! Amount of product sold
sol_buy: array(RAW,TIME) of real
! Amount of raw mat. bought
sol_pstock: array(PRODS,1..NT+1) of real ! Product stock levels
sol_rstock: array(RAW,1..NT+1) of real
! Raw mat. stock levels
pc: real! Cost of the proposal
end-declarations! Read proposal data from memory
initializations from "mempipe:noindex,sol"fsol_make sol_sell sol_buy sol_pstock sol_rstockpc
end-initializations! Add the new proposal to the master problem
nPROP(f)+=1
create(weight(f,nPROP(f)))
forall(p in PRODS,t in TIME) do
Prop_make(p,f,t,nPROP(f)):= sol_make(p,t)
Prop_sell(p,f,t,nPROP(f)):= sol_sell(p,t)
end-doforall(r in RAW,t in TIME) Prop_buy(r,f,t,nPROP(f)):= sol_buy(r,t)
forall(p in PRODS,t in 1..NT+1) Prop_pstock(p,f,t,nPROP(f)):= sol_pstock(p,t)
forall(r in RAW,t in 1..NT+1) Prop_rstock(r,f,t,nPROP(f)):= sol_rstock(r,t)
Prop_cost(f,nPROP(f)):= pc
writeln("Sol. for factory ", f, ":\n make:", sol_make, "\n sell:",
sol_sell, "\n buy:", sol_buy, "\n pstock: ", sol_pstock,"\n rstock: ", sol_rstock)end-procedure!----------------------------------------------------------procedure solve_master(phase: integer)
forall(f in FACT)
Convex(f):= sum (k in 1..nPROP(f)) weight(f,k) = 1
if phase=1 thenforall(p in PRODS,t in TIME)

Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 33

Mosel: multiple models and parallel solving

MxSell(p,t):=sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) excessS <= MXSELL(p,t)
minimize(excessS)
else
forall(p in PRODS,t in TIME)
MxSell(p,t):=sum(f in FACT,k in 1..nPROP(f)) Prop_sell(p,f,t,k)*weight(f,k) <=MXSELL(p,t)
maximize(sum(f in FACT, k in 1..nPROP(f)) Prop_cost(f,k) * weight(f,k))
end-ifwriteln("Master problem objective: ", getobjval)
write(" Weights:")
forall(f in FACT,k in 1..nPROP(f)) write(" ", getsol(weight(f,k)))writeln
end-procedure!----------------------------------------------------------procedure process_master_result
forall(p in PRODS,t in TIME) Price_sell(p,t):=getdual(MxSell(p,t))
forall(f in FACT) Price_convex(f):=getdual(Convex(f))
initializations to "raw:noindex"
Price_sell as "shmem:Price_sell"
end-initializations
end-procedure

Finally, the master model is completed by two subroutines for calculating the solution to theoriginal problem (calc_solution), Phase 3 of the decomposition algorithm, and printing out
the solution (print_solution). The solution to the original problem is obtained from thesolution values of the modfied master problem and the proposals generated by the subproblems.
declarations
REV: array(PRODS,TIME) of real
CMAKE: array(PRODS,FACT) of real
CBUY: array(RAW,TIME) of real
COPEN: array(FACT) of real

! Unit selling price of products! Unit cost to make product p! at factory f! Unit cost to buy raw materials! Fixed cost of factory f being! open for one period! Unit cost to store any product p! Unit cost to store any raw mat. r

CPSTOCK: real
CRSTOCK: real

Sol_make: array(PRODS,FACT,TIME) of real ! Solution value (products made)Sol_sell: array(PRODS,FACT,TIME) of real ! Solution value (product sold)Sol_buy: array(RAW,FACT,TIME) of real
! Solution value (raw mat. bought)
Sol_pstock: array(PRODS,FACT,1..NT+1) of real ! Sol. value (prod. stock)
Sol_rstock: array(RAW,FACT,1..NT+1) of real
! Sol. value (raw mat. stock)end-declarations
initializations from DATAFILECMAKE REV CBUY CPSTOCK CRSTOCK COPEN
end-initializations
function calc_solution: real
forall(p in PRODS,f in FACT,t in TIME) do
Sol_sell(p,f,t):=sum(k in 1..nPROP(f)) Prop_sell(p,f,t,k) * getsol(weight(f,k))
Sol_make(p,f,t):=sum(k in 1..nPROP(f)) Prop_make(p,f,t,k) * getsol(weight(f,k))
end-doforall(r in RAW,f in FACT,t in TIME) Sol_buy(r,f,t):=sum(k in 1..nPROP(f)) Prop_buy(r,f,t,k) * getsol(weight(f,k))
forall(p in PRODS,f in FACT,t in 1..NT+1) Sol_pstock(p,f,t):=sum(k in 1..nPROP(f)) Prop_pstock(p,f,t,k) * getsol(weight(f,k))
forall(r in RAW,f in FACT,t in 1..NT+1) Sol_rstock(r,f,t):=Dantzig-Wolfe decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 34

Mosel: multiple models and parallel solving

sum(k in 1..nPROP(f)) Prop_rstock(r,f,t,k) * getsol(weight(f,k))
returned:=sum(p in PRODS,f in FACT,t in TIME) REV(p,t) * Sol_sell(p,f,t) sum(p in PRODS,f in FACT,t in TIME) CMAKE(p,f) * Sol_make(p,f,t) sum(r in RAW,f in FACT,t in TIME) CBUY(r,t) * Sol_buy(r,f,t) sum(p in PRODS,f in FACT,t in 2..NT+1) CPSTOCK * Sol_pstock(p,f,t) sum(r in RAW,f in FACT,t in 2..NT+1) CRSTOCK * Sol_rstock(r,f,t)
end-function!----------------------------------------------------------procedure print_solution
writeln("Finished products:")
forall(f in FACT) do
writeln("Factory ", f, ":")
forall(p in PRODS) do
write(" ", p, ":")forall(t in TIME) write(strfmt(Sol_make(p,f,t),6,1), "(",
strfmt(Sol_pstock(p,f,t+1),5,1), ")")
writeln
end-doend-dowriteln("Raw material:")forall(f in FACT) do
writeln("Factory ", f, ":")
forall(r in RAW) do
write(" ", r, ": ")forall(t in TIME) write(strfmt(Sol_buy(r,f,t),6,1), "(",
strfmt(Sol_rstock(r,f,t+1),5,1), ")")
writeln
end-doend-dowriteln("Sales:")
forall(f in FACT) do
writeln("Factory ", f, ":")
forall(p in PRODS) do
write(" ", p, ": ")forall(t in TIME) write(strfmt(Sol_sell(p,f,t),4))
writeln
end-doend-dowriteln("\nComputation time: ", gettime)
end-procedure

5.3

Results
For the test cases that have been tried the solutions produced by our decomposition algorithm
are close to the optimal solution, but the latter is not always reached. The reason behind this isthat the decomposition algorithm is a sequence of iterations that may accumulate errors at
different levels—lowering the tolerances used as stopping criterion in the submodels most of thetime does not improve the solution. However, the configuration of the decomposition algorithm
itself shows some impact on the solution: in phases 1 and 2 one may choose, for instance, to stop
once the first submodel returns no improvement or continue until no more proposals are
generated. Generating more proposals sometimes helps finding a better solution, but it also
increases the number of times (sub)problems are solved and hence prolongates the solving time.

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 35

Mosel: multiple models and parallel solving

6

Benders decomposition: working with several different submodels
Benders decomposition is a decomposition method for solving large Mixed Integer Programming
problems. Instead of solving a MIP problem that may be too large for standard solution methods
all-in-one, we work with a sequence of linear and pure integer subproblems (the latter with asmaller number of constraints than the original problem). The description of the decomposition
algorithm below is taken from [Hu69] where the interested reader will also find proofs of its
finiteness and of the statement that it always results in the optimal solution.
Consider the following standard form of a mixed-integer programming problem (problem I).
minimize z = Cx + Dy
Ax + By ≥ bx, y ≥ 0, y integer
In the above, and all through this section, we use bold letters for vectors and matrices. For
PNCTVARPNINTVAR
Ci · xi + i=1
Di · yi , where NCTVAR denotes the number of
instance, Cx + Dy stands for i=1
continuous variables and NINTVAR the number of integer variables.
For given values y0 of y the problem above is reduced to a linear program (problem II)—we leaveout the constant term in the objective:
minimize Cx
Ax ≥ b − By0x≥0
The dual program of problem II is given by problem IId:
maximize u(b − By0 )uA ≤ Cu≥0
An interesting feature of the dual problem IId is that the feasible region of u is independent of y.
Furthermore, from duality theory it follows that if problem IId is infeasible or has no finite
optimum solution, then the original problem I has no finite optimum solution. Again fromduality theory we know that if problem IId has a finite optimum solution u∗ then this solution has
the same value as the optimum solution to the primal problem (that is, u∗ (b − By0 ) = Cx∗ ), and fora solution up (at a vertex p of the feasible region) we have up (b − By0 ) ≤ Cx∗ . Substituting thelatter into the objective of the original MIP problem we obtain a constraint of the formz ≥ up (b − By) + Dy
To obtain the optimum solution z∗ of the original MIP problem (problem I) we may use thefollowing partitioning algorithm that is known as Benders decomposition algorithm:
Step 0Benders decomposition

Find a u0 that satisfies u’A ≤ C
if no such u0 exists
then the original problem (I) has no feasible solution
else continue with Step 1
end-ifc 2009 Fair Isaac Corporation. All rights reserved.

page 36

Mosel: multiple models and parallel solving

Step 1

Solve the pure integer programminimize zz ≥ u0 (b − By) + Dyy ≥ 0, y integer
If z is unbounded from below, take a z to be any small value z0 .

Step 2

With the solution y0 of Step 1, solve the linear programmaximize u(b − By0 )uA ≤ Cu≥0
If u goes to inifinity withPu(b − By0 )then add the constraint i ui ≤ M, where M is a large positive constant, and resolve
the problem
end-if
Let the solution of this program be u00 .
if z0 − Dy0 ≤ u00 (b − By0 )then continue with Step 3
else return to Step 1 and add the constraint z0 ≥ Dy + u00 (b − By)
end-ifStep 3

With the solution y0 of Step 1, solve the linear programminimize Cx
Ax ≥ b − By0x≥0x0 and y0 are the optimum solution z∗ = Cx0 + Dy0

This algorithm is provably finite. It results in the optimum solution and at any time during its
execution lower and upper bounds on the optimum solution z∗ can be obtained.

6.1

A small example problem
Our implementation of Benders decomposition will solve the following example problem withNCTVAR = 3 continuous variables xi , NINTVAR = 3 integer variables yi , and NC = 4 inequality
constraints.
maximize 5 · x1 + 4 · x2 + 3 · x3 + 2 · y1 + 2 · y2 + 3 · y3x1 −2 · x1x2 + 5 · x3 + 3 · y1

+ 2 · y3 ≤ 5

4 · x2 + 3 · x3 −

y1 +

y2

≤7

− 2 · x3 +

y1 −

y2 +

y3 ≤ 43 · x1+ 5 · y1 + 5 · y2 + 5 · y3 ≤ −2x, y ≥ 0, y integer

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 37

Mosel: multiple models and parallel solving

6.2

Implementation
6.2.1

Master model

The master model reads in the data, defines the solution algorithm, coordinates thecommunication between the submodels, and prints out the solution at the end. For step 2 of thealgorithm (solving the dual problem with fixed integer variables) we have the choice to solve
either the primal problem and retrieve the dual solution values from the Optimizer or to define
the dual problem ourselves and solve it. Model parameter ALG lets the user choose between
these two options.
The main part of the master model looks as follows. Prior to the start of the solution algorithm
itself all submodels are compiled, loaded, and started so that in each step of the algorithm we
simply need to trigger the (re)solving of the corresponding submodel.
model "Benders (master model)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
ALG = 1

! 1: Use Benders decomposition (dual)! 2: Use Benders decomposition (primal)

DATAFILE = "bprob33.dat"
end-parameters
forward procedure start_solution
forward procedure solve_primal_int(ct: integer)
forward procedure solve_contforward function eval_solution: boolean
forward procedure print_solution
declarations
STEP_0=2
STEP_1=3
STEP_2=4
EVENT_SOLVED=6
EVENT_INFEAS=7
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

CtVars = 1..NCTVAR
IntVars = 1..NINTVAR
Ctrs = 1..NC! Continuous variables! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer
D: array(IntVars) of integer
Ctr: array(Ctrs) of linctr
CtrD: array(CtVars) of linctr
MC: array(range) of linctr! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables! Obj. coeff.s of discrete variables! Constraints of orig. problem! Constraints of dual problem! Constraints generated by alg.

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real
sol_obj: real

! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)! Objective function value (primal)

RM: range
stepmod: array(RM) of Model
end-declarations! Model indices! Submodels

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 38

Mosel: multiple models and parallel solving

initializations from DATAFILE
A B b C D
end-initializations! **** Submodels ****
initializations to "raw:noindex"! Save data for submodels
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C" D as "shmem:D"
end-initializations! Compile + load all submodels
if compile("benders_int.mos")<>0 then exit(1); end-if
create(stepmod(1)); load(stepmod(1), "benders_int.bim")
if compile("benders_dual.mos")<>0 then exit(2); end-if
if ALG=1 thencreate(stepmod(2)); load(stepmod(2), "benders_dual.bim")elsecreate(stepmod(0)); load(stepmod(0), "benders_dual.bim")
if compile("benders_cont.mos")<>0 then exit(3); end-if
create(stepmod(2)); load(stepmod(2), "benders_cont.bim")run(stepmod(0), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
end-if! Start the execution of the submodels
run(stepmod(1), "NINTVAR=" + NINTVAR + ",NC=" + NC)run(stepmod(2), "NCTVAR=" + NCTVAR + ",NINTVAR=" + NINTVAR + ",NC=" + NC)
forall(m in RM) dowait! Wait for "Ready" messages
ev:= getnextevent
if getclass(ev) <> EVENT_READY thenwriteln("Error occurred in a subproblem")exit(4)
end-ifend-do! **** Solution algorithm ****
start_solution
ct:= 1
repeatwriteln("\n**** Iteration: ", ct)
solve_primal_int(ct)solve_cont
ct+=1
until eval_solution
print_solution! 0. Initial solution for cont. var.s! 1. Solve problem with fixed cont.
! 2. Solve problem with fixed int.
!
Test for optimality! 3. Retrieve and display the solution

The subroutines starting the different submodels send a ‘start solving’ event and retrieve thesolution once the submodel solving has finished. For the generation of the start solution we needto choose the right submodel, according to the settings of the parameter ALG. If this problem isfound to be infeasible, then the whole problem is infeasible and we stop the execution of themodel.
! Produce an initial solution for the dual problem using a dummy objective
procedure start_solution
if ALG=1 then
! Start the problem solving
send(stepmod(2), STEP_0, 0)
else
send(stepmod(0), STEP_0, 0)
end-ifwait! Wait for the solution
ev:=getnextevent

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 39

Mosel: multiple models and parallel solving

if getclass(ev)=EVENT_INFEAS thenwriteln("Problem is infeasible")exit(6)
end-ifend-procedure!----------------------------------------------------------! Solve a modified version of the primal problem, replacing continuous! variables by the solution of the dual
procedure solve_primal_int(ct: integer)
send(stepmod(1), STEP_1, ct)! Start the problem solving
wait! Wait for the solution
ev:=getnextevent
sol_obj:= getvalue(ev)! Store objective function value
initializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
end-procedure! Retrieve the solution!----------------------------------------------------------! Solve the Step 2 problem (dual or primal depending on value of ALG)
! for given solution values of yprocedure solve_contsend(stepmod(2), STEP_2, 0)! Start the problem solving
wait! Wait for the solution
dropnextevent
initializations from "raw:noindex"
sol_u as "shmem:u"
end-initializations
end-procedure! Retrieve the solution

The master model also tests whether the termination criterion is fulfilled (function
eval_solution) and prints out the final solution (procedure print_solution). The latterprocedure also stops all submodels:function eval_solution: boolean
write("Test optimality: ", sol_obj - sum(i in IntVars) D(i)*sol_y(i)," = ", sum(j in Ctrs) sol_u(j)* (b(j) sum(i in IntVars) B(j,i)*sol_y(i)) )returned:= ( sol_obj - sum(i in IntVars) D(i)*sol_y(i) =sum(j in Ctrs) sol_u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i)) )writeln(if(returned, " : true", " : false"))end-function!----------------------------------------------------------procedure print_solution! Retrieve results
initializations from "raw:noindex"
sol_x as "shmem:x"
end-initializations
forall(m in RM) stop(stepmod(m))! Stop all submodels

write("\n**** Solution (Benders): ", sol_obj, "\nforall(i in CtVars) write(sol_x(i), " ")write(" y: ")forall(i in IntVars) write(sol_y(i), " ")writeln
end-procedure

Benders decomposition

x: ")c 2009 Fair Isaac Corporation. All rights reserved.

page 40

Mosel: multiple models and parallel solving

6.2.2

Submodel 1: fixed continuous variables

In the first step of the decomposition algorithm we need to solve a pure integer problem. Whenthe execution of this model is started it reads in the invariant data and sets up the variables. It
then halts at the wait statement (first line of the repeat-until loop) until the master model
sends it a (solving) event. At each invocation of solving this problem, the solution values of theprevious run of the continuous problem—read from memory—are used to define a new
constraint MC(k) for the integer problem. The whole model, and with it the endless loop into
which the solving is embedded, will be terminated only by the ‘stop model’ command from themaster model. The complete source of this submodel (file benders_int.mos) is listed below.
model "Benders (integer problem)"uses "mmxprs", "mmjobs"
parameters
NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
declarations
STEP_0=2
STEP_1=3
EVENT_SOLVED=6
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

IntVars = 1..NINTVAR
Ctrs = 1..NC! Discrete variables! Set of constraints (orig. problem)

B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
D: array(IntVars) of integer
MC: array(range) of linctr! Coeff.s of discrete variables! RHS values! Obj. coeff.s of discrete variables! Constraints generated by alg.

sol_u: array(Ctrs) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob.

y: array(IntVars) of mpvar
z: mpvar
end-declarations! Discrete variables! Objective function variable

initializations from "raw:noindex"
B as "shmem:B" b as "shmem:b" D as "shmem:D"
end-initializations
z is_free! Objective is a free variable
forall(i in IntVars) y(i) is_integer ! Integrality condition
forall(i in IntVars) y(i) <= BIGM
! Set (large) upper boundsend(EVENT_READY,0)! Model is ready (= running)

repeatwait
ev:= getnextevent
ct:= integer(getvalue(ev))
initializations from "raw:noindex"
sol_u as "shmem:u"
end-initializations! Add a new constraint
MC(ct):= z >= sum(i in IntVars) D(i)*y(i) +sum(j in Ctrs) sol_u(j)*(b(j) - sum(i in IntVars) B(j,i)*y(i))

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 41

Mosel: multiple models and parallel solving

minimize(z)! Store solution values of yforall(i in IntVars) sol_y(i):= getsol(y(i))
initializations to "raw:noindex"
sol_y as "shmem:y"end-initializations
send(EVENT_SOLVED, getobjval)
write("Step 1: ", getobjval, "\n y: ")forall(i in IntVars) write(sol_y(i), " ")write("\n Slack: ")forall(j in 1..ct) write(getslack(MC(j)), " ")writelnfflushuntil false
end-model

Since the problems we are solving differ only by a single constraint from one iteration to the next,
it may be worthwhile to save the basis of the solution to the root LP-relaxation (not the basis to
the MIP solution) and reload it for the next optimization run. However, for our small test case we
did not observe any improvements in terms of execution speed. For saving and re-reading thebasis, the call to minimize needs to be replaced by the following sequence of statements:
declarations
bas: basis
end-declarations
loadprob(z)loadbasis(bas)minimize(XPRS_TOP, z)
savebasis(bas)minimize(XPRS_GLB, z)

6.2.3

Submodel 2: fixed integer variables

The second step of our decomposition algorithm consists in solving a subproblem where all
integer variables are fixed to their solution values found in the first step. The structure of themodel implementing this step is quite similar to the previous submodel. When the model is run, it
reads the invariant data from memory and sets up the objective function. It then halts at the line
wait at the beginning of the loop to wait for a step 2 solving event sent by the master model. At
every solving iteration the constraints CTR are redefined using the coefficients read from memory
and the solution is written back to memory. Below follows the source of this model (filebenders_cont.mos).
model "Benders (continuous problem)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
declarations

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 42Mosel: multiple models and parallel solving

STEP_0=2
STEP_2=4
STEP_3=5
EVENT_SOLVED=6
EVENT_READY=8! Event codes sent to submodels

CtVars = 1..NCTVAR
IntVars = 1..NINTVAR
Ctrs = 1..NC! Continuous variables! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer
Ctr: array(Ctrs) of linctr! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables! Constraints of orig. problem

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)

x: array(CtVars) of mpvar
end-declarations! Continuous variables! Event codes sent by submodels

initializations from "raw:noindex"
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C"
end-initializations
Obj:= sum(i in CtVars) C(i)*x(i)send(EVENT_READY,0)! Model is ready (= running)! (Re)solve this model until it is stopped by event "STEP_3"
repeatwait
dropnextevent
initializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
forall(j in Ctrs)
Ctr(j):= sum(i in CtVars) A(j,i)*x(i) +sum(i in IntVars) B(j,i)*sol_y(i) >= b(j)minimize(Obj)! Solve the problem! Store values of u and x
forall(j in Ctrs) sol_u(j):= getdual(Ctr(j))
forall(i in CtVars) sol_x(i):= getsol(x(i))
initializations to "raw:noindex"
sol_u as "shmem:u" sol_x as "shmem:x"
end-initializations
send(EVENT_SOLVED, getobjval)
write("Step 2: ", getobjval, "\n u: ")forall(j in Ctrs) write(sol_u(j), " ")write("\n x: ")forall(i in CtVars) write(getsol(x(i)), " ")writelnfflushuntil false
end-model

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 43

Mosel: multiple models and parallel solving

6.2.4

Submodel 0: start solution

To start the decomposition algorithm we need to generate an initial set of values for thecontinuous variables. This can be done by solving the dual problem in the continuous variables
with a dummy objective function. A second use of the dual problem is for Step 2 of thealgorithm, replacing the primal model we have seen in the previous section. The implementationof this submodel takes into account these two cases: within the solving loop we test for the type(class) of event that has been sent by the master problem and choose the problem to be solved
accordingly.
The main part of this model is implemented by the following Mosel code (filebenders_dual.mos).
model "Benders (dual problem)"uses "mmxprs", "mmjobs"
parameters
NCTVAR = 3
NINTVAR = 3
NC = 4
BIGM = 1000
end-parameters
forward procedure save_solution
declarations
STEP_0=2
STEP_2=4
EVENT_SOLVED=6
EVENT_INFEAS=7
EVENT_READY=8! Event codes sent to submodels! Event codes sent by submodels

CtVars = 1..NCTVAR
IntVars = 1..NINTVAR
Ctrs = 1..NC! Continuous variables! Discrete variables! Set of constraints (orig. problem)

A: array(Ctrs,CtVars) of integer
B: array(Ctrs,IntVars) of integerb: array(Ctrs) of integer
C: array(CtVars) of integer! Coeff.s of continuous variables! Coeff.s of discrete variables! RHS values! Obj. coeff.s of continuous variables

sol_u: array(Ctrs) of real
sol_x: array(CtVars) of real
sol_y: array(IntVars) of real

! Solution of dual problem! Solution of primal prob. (cont.)! Solution of primal prob. (integers)

u: array(Ctrs) of mpvar
end-declarations! Dual variables

initializations from "raw:noindex"
A as "shmem:A" B as "shmem:B"b as "shmem:b" C as "shmem:C"
end-initializations
forall(i in CtVars) CtrD(i):= sum(j in Ctrs) u(j)*A(j,i) <= C(i)send(EVENT_READY,0)! Model is ready (= running)! (Re)solve this model until it is stopped by event "STEP_3"
repeatwait
ev:= getnextevent
Alg:= getclass(ev)if Alg=STEP_0 thenBenders decomposition! Produce an initial solution for thec 2009 Fair Isaac Corporation. All rights reserved.

page 44

Mosel: multiple models and parallel solving! dual problem using a dummy objective
maximize(sum(j in Ctrs) u(j))
if(getprobstat = XPRS_INF) thenwriteln("Problem is infeasible")send(EVENT_INFEAS,0)
else
write("**** Start solution: ")save_solution
end-if! Problem is infeasible

else! STEP 2: Solve the dual problem for! given solution values of yinitializations from "raw:noindex"
sol_y as "shmem:y"end-initializations
Obj:= sum(j in Ctrs) u(j)* (b(j) - sum(i in IntVars) B(j,i)*sol_y(i))
maximize(XPRS_PRI, Obj)if(getprobstat=XPRS_UNB) then
BigM:= sum(j in Ctrs) u(j) <= BIGM
maximize(XPRS_PRI, Obj)
end-ifwrite("Step 2: ")save_solution
BigM:= 0
end-if! Write solution to memory! Reset the ’BigM’ constraint

until false

This model is completed by the definition of the subroutine save_solution that writes thesolution to memory and informs the master model of it being available by sending the
EVENT_SOLVED message.
! Process solution dataprocedure save_solution! Store values of u and x
forall(j in Ctrs) sol_u(j):= getsol(u(j))
forall(i in CtVars) sol_x(i):= getdual(CtrD(i))
initializations to "raw:noindex"
sol_u as "shmem:u" sol_x as "shmem:x"
end-initializations
send(EVENT_SOLVED, getobjval)
write(getobjval, "\n u: ")forall(j in Ctrs) write(sol_u(j), " ")write("\n x: ")forall(i in CtVars) write(getdual(CtrD(i)), " ")writelnfflush
end-procedure
end-model

6.3

Results
The optimal solution to our small test problem has the objective function value 18.1852. Our
program produces the following output, showing that the problem is solved to optimality with 3
iterations (looping around steps 1 and 2) of the decomposition algorithm:**** Start solution: 4.05556

Benders decomposition

c 2009 Fair Isaac Corporation. All rights reserved.

page 45

Mosel: multiple models and parallel solving

u: 0.740741 1.18519 2.12963 0
x: 0.611111 0.166667 0.111111**** Iteration: 1
Step 1: -1146.15
y: 1000 0 0
Slack: 0
Step 2: 1007
u: 0 1 0 0
x: 0 251.75 0
Test optimality: -3146.15 = 1007 : false**** Iteration: 2
Step 1: 17.0185
y: 3 0 0
Slack: 0 -1.01852
Step 2: 12.5
u: 0 1 2.5 0
x: 0.5 2.5 0
Test optimality: 11.0185 = 12.5 : false**** Iteration: 3
Step 1: 18.1852
y: 2 0 0
Slack: 0 -5.18519 -0.185185
Step 2: 14.1852
u: 0.740741 1.18519 2.12963 0
x: 1.03704 2.22222 0.037037
Test optimality: 14.1852 = 14.1852 : true**** Solution (Benders): 18.1852
x: 1.03704 2.22222 0.037037
y: 2 0 0

7

Summary
The examples in this white paper show how to use the functionality provided by the mmjobs
module. Without giving an exhaustive overview of the technical possibilities they provide starting
points for implementation of, and experimentation with, parallel solving and othermulti-problem solution approaches.
Any solver module available for Mosel may be used in conjunction with mmjobs. However,parallel solving of multiple problems with the same solver is only possible if the underlying solver
can work in a multi-threaded environment. This is the case for Xpress-Optimizer and its derivates(QP, SLP, SP) as well as for Xpress-Kalis.

Bibliography[Hu69] T. C. Hu. Integer programming and network flows. Addison-Wesley, Reading, MA, 1969.
[PW94] Y. Pochet and L. A. Wolsey. Algorithms and Reformulations for Lot Sizing Problems. Technical Report Discussion
Paper 9427, CORE, Université Catholique de Louvain, 1994.
[Teb01] J. R. Tebboth. A computational study of Dantzig-Wolfe decomposition. PhD thesis, University of Buckingham,
2001.

Summaryc 2009 Fair Isaac Corporation. All rights reserved.

page 46