Good modeling practice with Mosel

Topics covered in this chapter:

The following recommendations for writing Mosel models establish some guidelines as to how to write ``good'' models with Mosel. By ``good'' we mean re-usability, readability, and perhaps most importantly, efficiency: when observing these guidelines you can expect to obtain the best possible performance of Mosel for the compilation and execution of your models.

Using constants and parameters

Many mathematical models start with a set of definitions like the following:

 NT:= 3
 Months:= {'Jan', 'Feb', 'Mar'}
 MAXP:= 8.4
 Filename= "mydata.dat"

If these values do not change later in the model, they should be defined as constants, allowing Mosel to handle them more efficiently:

 declarations
  NT = 3
  Months = {'Jan', 'Feb', 'Mar'}
  MAXP = 8.4
  Filename= "mydata.dat"
 end-declarations

If such constants may change with the model instance that is solved, their definition should be moved into the parameters block (notice that this possibility only applies to simple types, excluding sets or arrays):

 parameters
  NT = 3
  MAXP = 8.4
  Filename = "mydata.dat"
 end-parameters

Mosel interprets these parameters as constants, but their value may be changed at every execution of a model, e.g.

 mosel exec mymodel NT=5 MAXP=7.5 Filename="mynewdata.dat"

Naming sets

It is customary in mathematical models to write index sets as 1,...,N or the like. Instead of translating this directly into Mosel code like the following:

 declarations
  x: array(1..N) of mpvar
 end-declarations

 sum(i in 1..N) x(i) >= 10

it is recommended to name index sets:

 declarations
  RI = 1..N
  x: array(RI) of mpvar
 end-declarations

 sum(i in RI) x(i) >= 10

The same remark holds if several loops or operators use the same intermediate set(s). Instead of

 forall(i in RI | isodd(i)) x(i) is_integer
 forall(i in RI | isodd(i)) x(i) <= 5
 sum(i in RI | isodd(i)) x(i) >= 10

which calculates the same intermediate set of odd numbers three times, it is more efficient to define this set explicitly and calculate it only once:

 ODD:= union(i in RI | isodd(i)) {i}

 forall(i in ODD) x(i) is_integer
 forall(i in ODD) x(i) <= 5
 sum(i in ODD) x(i) >= 10

Alternatively, loops of the same type and with the same index set(s) may be regrouped to reduce the number of times that the sets are calculated:

 forall(i in RI | isodd(i)) do
  x(i) is_integer
  x(i) <= 5
 end-do
 sum(i in RI | isodd(i)) x(i) >= 10

Finalizing sets and dynamic arrays

The declaration of an array in Mosel has one of these two forms

1. Explicit declaration as sparse array by using one of the keywords dynamic or hashmap.
2. `Standard' declaration, resulting in a dense array that is either static (all index sets are known) or not fixed (some or all indexing sets are unknown at the point where the declaration takes place).

If an array is used to represent dense data one should avoid defining it as a sparse array as that uses more memory and is slower than the corresponding dense array.

In many optimization models, dense arrays are created as non-fixed arrays because their contents is initially unknown—but there is no real need to treat them as dynamic structures throughout the whole model as they remain unchanged once they have been initialized.

The automatic finalization mechanism of Mosel therefore transforms such initially dynamic sets/non-fixed arrays as to handle them more efficiently. As an additional advantage, set finalization allows Mosel to check for `out of range' errors that cannot be detected if the sets are allowed to grow dynamically.

• By default, initializations blocks finalize the sets they initialize and also the index sets of initialized dense arrays.
• Data of non-dynamic arrays is read before finalization of the index sets in order to create the arrays static.
• Arrays that are not explicitly declared as sparse arrays are only allocated when they are first accessed: this allows these arrays to be static even if their index sets are finalized after the declaration of the arrays.

So, code like the following example

 declarations
  S: set of string
  A,B: array(S) of real
  x: array(S) of mpvar
 end-declarations

 initializations from "mydata.dat"
  A
 end-initializations

 sum(s in S) B(s)*x(s)

where all arrays are declared as dense arrays that are not fixed (their size is not known at their declaration) but only A that is initialized using a data file really needs to be non-fixed, will be treated by Mosel as if you had written the following

 declarations
  S: set of string
  A: array(S) of real
 end-declarations

 initializations from "mydata.dat"
  A
 end-initializations

 finalize(S)

 declarations
  B: array(S) of real
  x: array(S) of mpvar
 end-declarations

That is, B and x are created as static arrays, making the access to the array entries more efficient.

As a general rule, the following sequence of actions gives better results (in terms of memory consumption and efficiency):

1. Declare data arrays and sets that are to be initialized from external sources.
2. Perform initializations of data.
3. Finalize all related sets.
4. Declare any other arrays indexed by these sets (including decision variable arrays).

Note: there are several possibilities to stop Mosel from applying automatic finalization to model objects:

• Declare arrays explicitly as dynamic or hashmap arrays. (See examples in Sections A transport example and Conditional variable creation and create.)
• Declare sets explicitly as dynamic in which case they cannot be finalized.
• Use control parameter autofinal to enable/disable automatic finalization locally:

   setparam("autofinal", false)
   initializations from "datafile.dat"
     ...
   end-initializations
   setparam("autofinal", true) 

• Use option noautofinal to disable automatic finalization globally for the whole model:

  model "modelname"
   options noautofinal

Ordering indices

Especially when working with sparse arrays, the sequence of their indices in loops should correspond as far as possible to the sequence given in their declaration. For example an array of variables declared by:

 declarations
  A,B,C: range
  x: array(A,B,C) of mpvar
 end-initializations

that is mostly used in expressions like sum(b in B, c in C, a in A) x(a,b,c) should preferrably be declared as

 declarations
  A,B,C: range
  x: array(B,C,A) of mpvar
 end-declarations

or alternatively the indices of the loops adapted to the order of indices of the variables.

Use of exists

The Mosel compiler is able to identify sparse loops and optimizes them automatically, such as in the following example:

 declarations
  I=1..1000
  J=1..500
  A: dynamic array(I,J) of real
  x: array(I,J) of mpvar
 end-declarations

 initializations from "mydata.dat"
  A
 end-initializations

 C:= sum(i in I,j in J | exists(A(i,j))) A(i,j)*x(i,j) = 0

Notice that we obtain the same definition for the constraint C with the following variant of the code, but no loop optimization takes place:

 C:= sum(i in I,j in J) A(i,j)*x(i,j) = 0

Here all index tuples are enumerated and the corresponding entries of A are set to 0. Similarly, if not all entries of x are defined, the missing entries are interpreted as 0 by the sum operator.

The following rules have to be observed for efficient use of the function exists, :

1. The arrays have to be indexed by named sets (here I and J):

    A: dynamic array(I,J) of real              ! can be optimized
    H: hashmap array(I,J) of real              ! can be optimized
    B: dynamic array(1..1000,1..500) of real   ! cannot be optimized

2. The same sets have to be used in the loops:

    forall(i in I,j in J | exists(A(i,j)))              ! fast
    K:=I; forall(i in K,j in 1..500 | exists(A(i,j)))   ! slow

3. The order of the sets has to be respected, particularly for dynamic arrays:

    forall(i in I,j in J | exists(A(i,j)))              ! fast
    forall(j in J,i in I | exists(H(i,j)))              ! slower
    forall(j in J,i in I | exists(A(i,j)))              ! slowest

4. The exists function calls have to be at the beginning of the condition:

    forall(i in I,j in I | exists(A(i,j)) and i+j<>10)  ! fast
    forall(i in J,j in J | i+j<>10 and exists(A(i,j)))  ! slow

5. The optimization does not apply to or conditions:

    forall(i in I,j in J | exists(A(i,j)) and i+j<>10)  ! fast
    forall(i in I,j in J | exists(A(i,j)) or i+j<>10)   ! slow

Structuring a model

Procedures and functions may be introduced to structure a model. For easy readability, the length of a subroutine should not exceed the length of one page (screen).

Large model files could even be split into several files (and combined using the include statement).

Transforming subroutines into user modules

The definitions of subroutines that are expensive in terms of execution time and are called very often (e.g. at every node of the Branch-and-Bound search) may be moved to a user module. Via the Mosel Native Interface it is possible to access and change all information in a Mosel model during its execution. See the Mosel Native Interface User Guide for a detailed description of how to define user modules.

Algorithm choice and parameter settings

The performance of the underlying solution algorithm has, strictly speaking, nothing to do with the efficiency of Mosel. But for completeness' sake the reader may be reminded that the subroutines getparam and setparam can be used to access and modify the current settings of parameters of Mosel and also those provided by modules, such as solvers.

The list of parameters defined by a module can be obtained with the Mosel command

exam -p module_name

With Xpress Optimizer (module mmxprs) you may try re-setting the following control parameters for the algorithm choice:

• LP: XPRS_PRESOLVE
• MIP: XPRS_PREPROBING, XPRS_MIPPRESOLVE, XPRS_CUTSTRATEGY, XPRS_HEUREMPHASIS, XPRS_SBEFFORT, XPRS_NODESELECTION
• Other useful parameters are the criteria for stopping the MIP search: XPRS_MAXNODE, XPRS_MAXMIPSOL, XPRS_TIMELIMIT, XPRS_SOLTIMELIMIT, the cutoff value (XPRS_MIPADDCUTOFF, XPRS_MIPABSCUTOFF), and various tolerance settings (e.g. XPRS_MIPTOL).

Refer to the Xpress Optimizer Reference Manual for more detail.

You may also add priorities or preferred branching directions with the procedure setmipdir (documented in the chapter on mmxprs in the Mosel Reference Manual).