Modeling in Mosel

Topics covered in this chapter:

Basic formulation

Nonlinear capabilities in Mosel are provided by the mmxnlp module. Please refer to the module documentation for more details. This chapter provides a short introduction only. Note that Mosel modeling for FICO Xpress Nonlinear and FICO Xpress Global work exactly the same, s.t. you can easily switch between both solvers.

The model uses the Mosel module mmxnlp which contains the extensions required for modeling general non-linear expressions. This automatically loads the mmxprs module, so there is no need to include this explicitly as well.

model "Polygon"
uses "mmxnlp"

We can design the model to work for any number of sides, so one way to do this is to set the number of sides of the polygon as a parameter.

parameters
  N=5
end-parameters

The meanings of most of these declarations will become apparent as the modeling progresses.

declarations
  area: nlctr
  rho: array(1..N) of mpvar
  theta: array(1..N) of mpvar
  objdef: mpvar
  D: array(1..N,1..N) of nlctr
end-declarations

• The distances are described as "rho", to distinguish them from the default names for the rows in the generated matrix (which are R1, R2, etc).
• The types nlctr (nonlinear constraint) are defined by the mmxnlp module.

area := sum(i in 2..N-1) (rho(i) * rho(i-1) * sin(theta(i)-theta(i-1)))*0.5

This uses the normal Mosel sum function to calculate the area. Notice that the formula is written in essentially the same way as normal, including the use of the sin function. Because the argument to the function is not a constant, Mosel will not try to evaluate the function yet; instead, it will be evaluated as part of the optimization process.

area is a Mosel object of type nlctr.

objdef = area
objdef is_free

What we really want to do is to maximize area. However, although Xpress NonLinear is happy in principle with a non-linear objective function, the Xpress Optimizer is not, unless it is handled in a special way. Xpress NonLinear therefore imposes the requirement that the objective function itself must be linear. This is not really a restriction, because – as in this case – it is easy to reformulate a non-linear objective function as an apparently linear one. Simply replace the function by a new mpvar and then maximize the value of the mpvar. In general, because the objective could have a positive or negative value, we make the variable free, so that it can take any value. In this example, we say:

objdef = area	defining the variable objdef to be equal to the non-linear expression area
objdef is_free	defining objdef to be a free variable
maximize(objdef)	maximizing the linear objective

This is firstly setting the standard bounds on the variables rho and theta. To reduce problems with sides of zero length, we impose a minimum of 0.1 on rho(i) instead of the default minimum of zero.

forall (i in 1..N-1) do
  rho(i) >= 0.1
  rho(i) <= 1
  setinitval(rho(i), 4*i*(N+1-i)/((N+1)^2))
  setinitval(theta(i), M_PI*i/N)
end-do

We also give Xpress NonLinear initial values by using the setinitval procedure. The first argument is the name of the variable, and the second is the initial value to be used. The initial values for theta are divided equally between 0 and π. The initial values for rho are designed to go from 0 (when i=0 or N) to 1 (when i is about half way) and back.

forall (i in 1..N-2, j in i+1..N-1) do
  D(i,j) := rho(i)^2 + rho(j)^2 - rho(i)*rho(j)*2*cos(theta(j)-theta(i)) <= 1
end-do

This is creating the constraints D(i,j) which constrain the other sides of the triangles to be ≤ 1.

These constraints could be made anonymous – that is, the assignment to an object of type nlctr could be omitted – but then it would not be possible to report the values.

forall (i in 2..N-1) do
  theta(i) >= theta(i-1) + 0.01
end-do

These anonymous constraints put the values of the theta variables in non-decreasing order. To avoid problems with triangles which have zero angles, we make each bearing at least 0.01 greater than its predecessor.

This is the boundary condition on the bearing of the final vertex.

theta(N-1) <= M_PI

Setting up and solving the problem

loadprob(objdef)

This procedure loads the currently-defined non-linear problem into the Xpress NonLinear optimization framework. This includes any purely linear part. Where a constraint has a linear expression as its left or right hand side, that linear expression will be retained as linear relationships (constant coefficients) in the matrix. Thus, for example, in the anonymous constraint defining objdef, the objdef coefficient will be identified as a linear term and will appear as a separate item in the problem.

maximise

Optimization is carried out with the maximise or minimise procedures. They can take a string parameter – for example maxmimise("b") – as described in the Xpress NonLinear and Xpress Optimizer reference manuals.

With the default settings of the parameters, you will see usually nothing from the optimizer. The following parameters affect what is produced:

xnlp_verbose	Normally set to false. If set to true, it produces standard Xpress NonLinear iteration logging.
xprs_verbose	Normally set to false. If set to true, then information from the optimizer will also be output.
xslp_log	Normally set to -1. If set to 0, limited information is output from the SLP iterations. Settings of 1 or greater produce progressively more information for each SLP iteration.
xslp_slplog	If xslp_log is set to 0, this determines the frequency with which SLP progress is reported. The default is 10, which means that it prints every 10 SLP iterations.

Looking at the results

Within Mosel, the values of the variables and named constraints can be obtained using the getsol, getslack and similar functions. A simple report lists just the area and the positions of the vertices:

writeln("Area = ", getobjval)
forall (i in 1..N-1) do
  writeln("V", i, ": r=", getsol(rho(i)), " theta=", getsol(theta(i)))
end-do

This produces the following result for the case N=5:

Area = 0.657166
V1: r=0.616416 theta=0.703301
V2: r=1 theta=1.33111
V3: r=1 theta=1.96079
V4: r=0.620439 theta=2.58648

User functions

If an analytic description of the model is not available, it is possible to use black box functions implemented either directly in Mosel or by any external application.

The area calculation of the example could be implemented in Mosel as

public function MoselArea(I: array(Indices: range, Types: set of string) of real): real
  returned := (sum (i in 2..N-1) (I(i,"rho")*I(i-1,"rho")*sin(I(i,"theta")-I(i-1,"theta")))) * 0.5
 end-function

The user function is linked to the model using a user function object

declarations
  AreaFunction : userfunc
 end-declarations
 AreaFunction := userfuncMosel("MoselArea")

The arguments, which can be any expression, are passed down using an array of expressions

declarations
  FunctionArg: array(RN,{"rho","theta"}) of nlctr
 end-declarations
 forall (i in 1..N) do
  FunctionArg(i, "rho")   := rho(i)
  FunctionArg(i, "theta") := theta(i)
 end-do

Once the user function is declared and the arguments built, the user function is added to the model using F:

Area := F(AreaFunction,FunctionArg)

The function arguments are copied at the point when the F function is used, any later changes to the arrays holding the arguments are ignored.

Parallel evaluation of Mosel user functions

It is possible to use parallel evaluations of simple Mosel functions that return a single real value. These functions may take an arbitrary array of nlctr expressions as input. It is the modeler's responsibility to ensure that the user functions to be called in parallel are thread-safe (i.e., they do not depend upon shared resources). Assuming the name of the user function is MyFunc, the user function before enabling the parallel version is expected to be declared as usefuncMosel('MyFunc').

In order for mmxnlp to be able to utilize parallel user function evaluations, the user function must be implemented as a public function in a Mosel package. Any initialization necessary to enable the evaluation of the user function should be performed as part of the package initialization (which is the code in in the main body of the package).

To enable parallel evaluations, a parallel enabled version of the user function needs to be generated using the mmxnlp procedure generateUFparallel, which takes two arguments: the compiled package .bim name implementing the user function and the name of the user function within the package. It is good practice to use a separate Mosel model to perform this generation, keeping it separate from the main model. Multiple generated parallel user functions may be used within a single model.

The generator will produce a single Mosel file, the Mosel package MyFunc_master. This package also includes the worker model which will be responsible for the user function evaluations and will be resident in memory during the execution. The package also implements the parallel version of the user function, called MyFunc_parallel.

After compiling and including the master package into your model, it is this function that should be used in the actual model as userfuncMosel('MyFunc_parallel',XSLP_DELTAS). In most cases, no other modifications are necessary, as the parallel function will detect the number of threads in the system and will start that many worker threads automatically. These will be shut down when your model finishes. Each worker's initialization code is performed only once, at the time of its first execution.

It may be necessary to explicitly start the worker threads, either to control the number of threads used, or to pass specific parameter settings to the user function package. This can be done by the procedure MyFunc_StartWorkers( ThreadCount : integer, UfPackageParameters : string ). In case it is necessary to stop the workers, the procedure MyFunc_StopWorkers may be used.

In case the user functions are computationally very expensive, by modifying the connection string in the generated module it is possible to utilize distributed/cloud-based computation of the user functions.

The worker model will only be compiled into memory during execution, but may be modified as necessary within the master model. For debugging purposes, it may be practical to redirect the worker to a file.