2.2. More Complex Models

2.2.1. Arrays and Sets

Almost always we are interested in building models where the number of constraints and variables is dependent on the input data. In order to do so we will usually use arrays.

Consider a simple finite element model for modelling temperatures on a rectangular sheet of metal. We approximate the temperatures across the sheet by breaking the sheet into a finite number of elements in a two-dimensional matrix. A model is shown in Listing 2.2.1. It declares the width w and height h of the finite element model. The declaration

set of int: HEIGHT = 0..h;
set of int: CHEIGHT = 1..h-1;
set of int: WIDTH = 0..w;
set of int: CWIDTH = 1..w-1;  
array[HEIGHT,WIDTH] of var float: t; % temperature at point (i,j)

declares four fixed sets of integers describing the dimensions of the finite element model: HEIGHT is the whole height of the model, while CHEIGHT is the centre of the height omitting the top and bottom, WIDTH is the whole width of the model, while CWIDTH is the centre of the width omitting the left and right sides. Finally a two dimensional array of float variables t with rows numbered \(0\) to \(h\) (HEIGHT) and columns \(0\) to \(w\) (WIDTH), to represent the temperatures at each point in the metal plate. We can access the element of the array in the \(i^{th}\) row and \(j^{th}\) column using an expression t[i,j].

Laplace’s equation states that when the plate reaches a steady state the temperature at each internal point is the average of its orthogonal neighbours. The constraint

% Laplace equation: each internal temp. is average of its neighbours
constraint forall(i in CHEIGHT, j in CWIDTH)(
              4.0*t[i,j] = t[i-1,j] + t[i,j-1] + t[i+1,j] + t[i,j+1]);

ensures that each internal point \((i,j)\) is the average of its four orthogonal neighbours. The constraints

% edge constraints
constraint forall(i in CHEIGHT)(t[i,0] = left);
constraint forall(i in CHEIGHT)(t[i,w] = right);
constraint forall(j in CWIDTH)(t[0,j] = top);
constraint forall(j in CWIDTH)(t[h,j] = bottom);

restrict the temperatures on each edge to be equal, and gives these temperatures names: left, right, top and bottom. While the constraints

% corner constraints
constraint t[0,0]=0.0;
constraint t[0,w]=0.0;
constraint t[h,0]=0.0;
constraint t[h,w]=0.0; 

ensure that the corners (which are irrelevant) are set to 0.0. We can determine the temperatures in a plate broken into 5 \(\times\) 5 elements with left, right and bottom temperature 0 and top temperature 100 with the model shown in Listing 2.2.1.

int: w = 4;
int: h = 4;

% arraydec
set of int: HEIGHT = 0..h;
set of int: CHEIGHT = 1..h-1;
set of int: WIDTH = 0..w;
set of int: CWIDTH = 1..w-1;  
array[HEIGHT,WIDTH] of var float: t; % temperature at point (i,j)
var float: left;   % left   edge temperature
var float: right;  % right  edge temperature
var float: top;    % top    edge temperature
var float: bottom; % bottom edge temperature

% equation
% Laplace equation: each internal temp. is average of its neighbours
constraint forall(i in CHEIGHT, j in CWIDTH)(
              4.0*t[i,j] = t[i-1,j] + t[i,j-1] + t[i+1,j] + t[i,j+1]);
% sides
% edge constraints
constraint forall(i in CHEIGHT)(t[i,0] = left);
constraint forall(i in CHEIGHT)(t[i,w] = right);
constraint forall(j in CWIDTH)(t[0,j] = top);
constraint forall(j in CWIDTH)(t[h,j] = bottom);
% corners
% corner constraints
constraint t[0,0]=0.0;
constraint t[0,w]=0.0;
constraint t[h,0]=0.0;
constraint t[h,w]=0.0; 
left = 0.0;
right = 0.0;
top = 100.0;
bottom = 0.0;

solve satisfy;

output [ show_float(6, 2, t[i,j]) ++ 
         if j == w then "\n" else " " endif |
         i in HEIGHT, j in WIDTH
];

Running the command

$ minizinc --solver cbc laplace.mzn

gives the output

 -0.00 100.00 100.00 100.00  -0.00
 -0.00  42.86  52.68  42.86  -0.00
 -0.00  18.75  25.00  18.75  -0.00
 -0.00   7.14   9.82   7.14  -0.00
 -0.00  -0.00  -0.00  -0.00  -0.00
----------

set of <type-inst> : <var-name> ;

{ <expr-1>, ..., <expr-n> }

<expr-1> .. <expr-2>

Our cake baking problem is an example of a very simple kind of production planning problem. In this kind of problem we wish to determine how much of each kind of product to make to maximise the profit where manufacturing a product consumes varying amounts of some fixed resources. We can generalise the MiniZinc model in Listing 2.1.3 to handle this kind of problem with a model that is generic in the kinds of resources and products. The model is shown in Listing 2.2.2 and a sample data file (for the cake baking example) is shown in Listing 2.2.3.

% Products to be produced
enum Products;  
% profit per unit for each product
array[Products] of int: profit;
% Resources to be used
enum Resources;
% amount of each resource available
array[Resources] of int: capacity; 

% units of each resource required to produce 1 unit of product
array[Products, Resources] of int: consumption; 
constraint assert(forall (r in Resources, p in Products) 
           (consumption[p,r] >= 0), "Error: negative consumption");

% bound on number of Products
int: mproducts = max (p in Products) 
                     (min (r in Resources where consumption[p,r] > 0) 
                          (capacity[r] div consumption[p,r]));

% Variables: how much should we make of each product
array[Products] of var 0..mproducts: produce;
array[Resources] of var 0..max(capacity): used;

% Production cannot use more than the available Resources:
constraint forall (r in Resources) (     
      used[r] = sum (p in Products)(consumption[p, r] * produce[p]) 
);
constraint forall (r in Resources) (     
      used[r] <= capacity[r]
);    

% Maximize profit
solve maximize sum (p in Products) (profit[p]*produce[p]);

output [ "\(p) = \(produce[p]);\n" | p in Products ] ++
       [ "\(r) = \(used[r]);\n" | r in Resources ];

% Data file for simple production planning model
Products = { BananaCake, ChocolateCake };
profit = [400, 450]; % in cents

Resources = { Flour, Banana, Sugar, Butter, Cocoa };
capacity = [4000, 6, 2000, 500, 500];

consumption= [| 250, 2, 75, 100, 0,
              | 200, 0, 150, 150, 75 |];

The new feature in this model is the use of enumerated types. These allow us to treat the choice of resources and products as parameters to the model. The first item in the model

enum Products;

declares Products as an unknown set of products.

enum <var-name> ;

enum <var-name> = { <var-name-1>, ..., <var-name-n> } ;

The second item declares an array of integers:

array[Products] of int: profit;

The index set of the array profit is Products. This means that only elements of the set Products can be used to index the array.

The elements of an enumerated type of \(n\) elements act very similar to the integers \(1\dots n\). They can be compared, they are ordered, by the order they appear in the enumerated type definition, they can be iterated over, they can appear as indices of arrays, in fact they can appear anywhere an integer can appear.

In the example data file we have initialized the array using a list of integers

Products = { BananaCake, ChocolateCake };
profit = [400,450];

meaning the profit for a banana cake is 400, while for a chocolate cake it is 450. Internally BananaCake will be treated like the integer 1, while ChocolateCake will be treated like the integer 2. The expression [400,500] represents a literal one-dimensional array. In MiniZinc, the index set of literal arrays always starts at 1 (this is similar to other mathematically inspired languages such as MATLAB, Julia or Mathematica). While MiniZinc does not provide an explicit list type, one-dimensional arrays with an index set 1..n behave like lists, and we will sometimes refer to them as lists.

In a similar fashion, in the next two items we declare a set of resources Resources, and an array capacity which gives the amount of each resource that is available.

More interestingly, the item

array[Products, Resources] of int: consumption;

declares a 2-D array consumption. The value of consumption[p,r] is the amount of resource r required to produce one unit of product p. Note that the first index is the row and the second is the column.

The data file contains an example initialization of a 2-D array:

consumption= [| 250, 2, 75,  100, 0,
              | 200, 0, 150, 150, 75 |];

Notice how the delimiter | is used to separate rows. As for one-dimensional array literals, indexing of two-dimensional array literals also starts at 1.

array [ <index-set-1>, ..., <index-set-n> ] of <type-inst>

[ <expr-1>, ..., <expr-n> ]

[| <expr-1-1>, ..., <expr-1-n> |
   ...                         |
   <expr-m-1>, ..., <expr-m-n> |]

array<n>d(<index-set-1>, ..., <index-set-n>, <list>)

The next item in the model defines the parameter mproducts. This is set to an upper-bound on the number of products of any type that can be produced. This is quite a complex example of nested array comprehensions and aggregation operators. We shall introduce these before we try to understand this item and the rest of the model.

First, MiniZinc provides list comprehensions similar to those provided in many functional programming languages, or Python. For example, the list comprehension [i + j | i, j in 1..3 where j < i] evaluates to [2 + 1, 3 + 1, 3 + 2] which is [3, 4, 5]. Of course [3, 4, 5] is simply an array with index set 1..3.

MiniZinc also provides set comprehensions which have a similar syntax: for instance, {i + j | i, j in 1..3 where j < i} evaluates to the set {3, 4, 5}.

[ <expr> | <generator-exp> ]

<generator>
<generator> where <bool-exp>

<identifier>, ..., <identifier> in <array-exp>

Second, MiniZinc provides a number of built-in functions that take a one-dimensional array and aggregate the elements. Probably the most useful of these is forall. This takes an array of Boolean expressions (that is, constraints) and returns a single Boolean expression which is the logical conjunction of the Boolean expressions in the array.

For example, consider the expression

forall( [a[i] != a[j] | i,j in 1..3 where i < j])

where a is an arithmetic array with index set 1..3. This constrains the elements in a to be pairwise different. The list comprehension evaluates to [ a[1] != a[2], a[1] != a[3], a[2] != a[3] ] and so the forall function returns the logical conjunction a[1] != a[2] /\ a[1] != a[3] /\ a[2] != a[3].

The third, and final, piece in the puzzle is that MiniZinc allows a special syntax for aggregation functions when used with an array comprehension. Instead of writing

forall( [a[i] != a[j] | i,j in 1..3 where i < j])

the modeller can instead write the more mathematical looking

forall (i,j in 1..3 where i < j) (a[i] != a[j])

The two expressions are completely equivalent: the modeller is free to use whichever they feel looks most natural.

<agg-func> ( <generator-exp> ) ( <exp> )

<agg-func> ( [ <exp> | <generator-exp> ] )

We are now in a position to understand the rest of the simple production planning model shown in Listing 2.2.2. For the moment ignore the item defining mproducts. The item afterwards:

array[Products] of var 0..mproducts: produce;

defines a one-dimensional array produce of decision variables. The value of produce[p] will be set to the amount of product p in the optimal solution. The next item

array[Resources] of var 0..max(capacity): used;

defines a set of auxiliary variables that record how much of each resource is used. The next two constraints

constraint forall (r in Resources)
           (used[r] = sum (p in Products) (consumption[p, r] * produce[p]));
constraint forall (r in Resources)(used[r] <= capacity[r] );

compute in used[r] the total consumption of the resource r and ensure it is less than the available amount. Finally, the item

solve maximize sum (p in Products) (profit[p]*produce[p]);

indicates that this is a maximisation problem and that the objective to be maximised is the total profit.

We now return to the definition of mproducts. For each product p the expression

(min (r in Resources where consumption[p,r] > 0)
     (capacity[r] div consumption[p,r]))

determines the maximum amount of p that can be produced taking into account the amount of each resource r and how much of r is required to produce the product. Notice the use of the filter where consumption[p,r] > 0 to ensure that only resources required to make the product are considered so as to avoid a division by zero error. Thus, the complete expression

int: mproducts = max (p in Products)
                     (min (r in Resources where consumption[p,r] > 0)
                          (capacity[r] div consumption[p,r]));

computes the maximum amount of any product that can be produced, and so this can be used as an upper bound on the domain of the decision variables in produce.

Finally notice the output item is more complex, and uses list comprehensions to create an understandable output. Running

$ minizinc --solver gecode prod-planning.mzn prod-planning-data.dzn

results in the output

BananaCake = 2;
ChocolateCake = 2;
Flour = 900;
Banana = 4;
Sugar = 450;
Butter = 500;
Cocoa = 150;
----------
==========

2.2.2. Global Constraints

MiniZinc includes a library of global constraints which can also be used to define models. An example is the alldifferent constraint which requires all the variables appearing in its argument to be pairwise different.

include "alldifferent.mzn";

var 1..9: S;
var 0..9: E;
var 0..9: N;
var 0..9: D;
var 1..9: M;
var 0..9: O;
var 0..9: R;
var 0..9: Y;

constraint           1000 * S + 100 * E + 10 * N + D
                   + 1000 * M + 100 * O + 10 * R + E
       = 10000 * M + 1000 * O + 100 * N + 10 * E + Y;

constraint alldifferent([S,E,N,D,M,O,R,Y]);

solve satisfy;

output ["   \(S)\(E)\(N)\(D)\n",
        "+  \(M)\(O)\(R)\(E)\n",
        "= \(M)\(O)\(N)\(E)\(Y)\n"];

The SEND+MORE=MONEY problem requires assigning a different digit to each letter so that the arithmetic constraint holds. The model shown in Listing 2.2.4 uses the constraint expression alldifferent([S,E,N,D,M,O,R,Y]) to ensure that each letter takes a different digit value. The global constraint is made available in the model using include item

include "alldifferent.mzn";

which makes the global constraint alldifferent usable by the model. One could replace this line by

include "globals.mzn";

which includes all globals.

A list of all the global constraints defined for MiniZinc is included in the release documentation. See Global Constraints for a description of some important global constraints.

2.2.3. Conditional Expressions

MiniZinc provides a conditional if-then-else-endif expression. An example of its use is

int: r = if y != 0 then x div y else 0 endif;

which sets r to x divided by y unless y is zero in which case it sets it to zero.

if <bool-exp> then <exp-1> else <exp-2> endif

include "alldifferent.mzn";

int: S;
int: N = S * S;
int: digs = ceil(log(10.0,int2float(N))); % digits for output

set of int: PuzzleRange = 1..N;
set of int: SubSquareRange = 1..S;

array[1..N,1..N] of 0..N: start; %% initial board 0 = empty
array[1..N,1..N] of var PuzzleRange: puzzle;

% fill initial board
constraint forall(i,j in PuzzleRange)(
    if start[i,j] > 0 then puzzle[i,j] = start[i,j] else true endif );

% All different in rows 
constraint forall (i in PuzzleRange) (
                   alldifferent( [ puzzle[i,j] | j in PuzzleRange ]) ); 

% All different in columns.
constraint forall (j in PuzzleRange) (
                   alldifferent( [ puzzle[i,j] | i in PuzzleRange ]) ); 

% All different in sub-squares:
constraint
        forall (a, o in SubSquareRange)(
                alldifferent( [ puzzle[(a-1) *S + a1, (o-1)*S + o1] |
                                        a1, o1 in SubSquareRange ] ) );

solve satisfy;

output  [ show_int(digs,puzzle[i,j]) ++ " " ++ 
         if j mod S == 0 then " " else "" endif ++
         if j == N then
             if i != N then 
                if i mod S == 0 then "\n\n" else "\n" endif 
         else "" endif else "" endif  
         | i,j in PuzzleRange ] ++ ["\n"];

S=3;
start=[|
0, 0, 0, 0, 0, 0, 0, 0, 0|
0, 6, 8, 4, 0, 1, 0, 7, 0|
0, 0, 0, 0, 8, 5, 0, 3, 0|
0, 2, 6, 8, 0, 9, 0, 4, 0|
0, 0, 7, 0, 0, 0, 9, 0, 0|
0, 5, 0, 1, 0, 6, 3, 2, 0|
0, 4, 0, 6, 1, 0, 0, 0, 0|
0, 3, 0, 2, 0, 7, 6, 9, 0|
0, 0, 0, 0, 0, 0, 0, 0, 0|];

Conditional expressions are very useful in building complex models, or complex output. Consider the model of Sudoku problems shown in Listing 2.2.5. The initial board positions are given by the start parameter where 0 represents an empty board position. This is converted to constraints on the decision variables puzzle using the conditional expression

constraint forall(i,j in PuzzleRange)(
     if start[i,j] > 0 then puzzle[i,j] = start[i,j] else true endif );

Conditional expressions are also very useful for defining complex output. In the Sudoku model of Listing 2.2.5 the expression

if j mod S == 0 then " " else "" endif

inserts an extra space between groups of size S. The output expression also uses conditional expressions to add blank lines after each S lines. The resulting output is highly readable.

The remaining constraints ensure that the numbers appearing in each row and column and \(S \times S\) subsquare are all different.

One can use MiniZinc to search for all solutions to a satisfaction problem (solve satisfy). In the MiniZinc IDE, this can be achieved by unticking the “Stop after this many solutions” option in the “User-defined behavior” part of the solver configuration pane (see Section 3.2.2.2). On the command line, one can use the flag -a or --all-solutions. Running

$ minizinc --solver gecode --all-solutions sudoku.mzn sudoku.dzn

results in

 5 9 3  7 6 2  8 1 4
 2 6 8  4 3 1  5 7 9
 7 1 4  9 8 5  2 3 6

 3 2 6  8 5 9  1 4 7
 1 8 7  3 2 4  9 6 5
 4 5 9  1 7 6  3 2 8

 9 4 2  6 1 8  7 5 3
 8 3 5  2 4 7  6 9 1
 6 7 1  5 9 3  4 8 2
----------
==========

The line ========== is output when the system has output all possible solutions, here verifying that there is exactly one.

2.2.4. Enumerated Types

Enumerated types allows us to build models that depend on a set of objects which are part of the data, or are named in the model, and hence make models easier to understand and debug. We have introduce enumerated types or enums briefly, in this subsection we will explore how we can use them more fully, and show some of the built in functions for dealing with enumerated types.

Let’s revisit the problem of coloring the graph of Australia from Basic Modelling in MiniZinc.

enum Color;
var Color: wa; 
var Color: nt; 
var Color: sa;
var Color: q;
var Color: nsw;
var Color: v;
var Color: t;
constraint wa != nt /\ wa != sa /\ nt != sa /\ nt != q /\ sa != q;
constraint sa != nsw /\ sa != v /\ q != nsw /\ nsw != v; 
solve satisfy;

The model shown in Listing 2.2.7 declares an enumerated type Color which must be defined in the data file. Each of the state variables is declared to take a value from this enumerated type. Running this program using

$ minizinc --solver gecode -D"Color = { red, yellow, blue };" aust-enum.mzn

might result in output

wa = blue;
nt = yellow;
sa = red;
q = blue;
nsw = yellow;
v = blue;
t = red;

<enum-name> : <var-name>
<l>..<u> : <var-name>

var <enum-name> : <var-name>
var <l>..<u> : <var-name>

A key behaviour of enumerated types is that they are automatically coerced to integers when they are used in a position expecting an integer. For example, this allows us to use global constraints defined on integers, such as

global_cardinality_low_up([wa,nt,sa,q,nsw,v,t],
                          [red,yellow,blue],[2,2,2],[2,2,3]);

This requires at least two states to be colored each color and three to be colored blue.

It is often useful to define an enumerated type by extending another enumerated type. A good example for this is adding a “don’t care” value to a type. For instance, consider a typical assignment problem, where we have to assign chores to flatmates. The names of the flatmates are given as an enumerated type Flatmates, and the chores as an enumerated type Chores.

In the case where not every flatmate needs to be assigned a job, and not every chore needs to be assigned to a flatmate, we need to introduce a value that signals “no assignment”. We can then use the alldifferent_except constraint to constrain the assignment.

include "alldifferent_except.mzn";

enum Flatmates = { Anne, Bert, Ceci, Dave };
enum Chores = { Cooking, Vacuuming, Bathroom, Kitchen, Rubbish };

enum ChoreOrNothing = C(Chores) ++ { Nothing };

array[Flatmates] of var ChoreOrNothing: assignment;

constraint alldifferent_except(assignment, {Nothing});

output [ show(f)++":\t"++
         if fix(assignment[f])=Nothing then "Nothing\n"
         else show(C^-1(assignment[f]))++"\n" endif
       | f in Flatmates];

The model is shown in Listing 2.2.8. It introduces a new enumerated type ChoreOrNothing, which is a combination of the Chores type and a new constructor Nothing. The Chores type is injected into the new type via a constructor function. In this case, we chose the name C for the constructor function, but you can use any valid identifier you like. The constructor function can then be used to coerce a value of type Chores into a value of the new type ChoreOrNothing. For example, we could write a constraint such as

constraint assignment[Bert] = C(Vacuuming);

This is required because the type ChoreOrNothing is considered different from the type Chores by the compiler, in order to avoid any ambiguities or accidental type conversions. The compiler also generates inverses of all constructor functions, to map back from the extended types into the base types. This can be seen in the output item, where C^-1 is used to map from ChoreOrNothing back to Chores. The inverse function can be written as ASCII using the ^-1 notation, or using the unicode symbols ⁻¹, as in C⁻¹.

In the example above, we could have used option types, instead of adding a new value that acts essentially like the absent value <>. However, extended enum types have some advantages in this scenario. For example, we can extend a type by multiple extra values (say, different reasons for not being assigned any chore). We can also combine multiple enumerated types into one. Consider the scenario in Listing 2.2.9 where we can assign chores to flatmates and/or robots. Flatmates have preferences for certain chores, while robots have integer benefits. The model also shows how to test if a value of the extended enum in fact belongs to the underlying enum: In the output item, the test w in F(Flatmates) uses the F constructor to translate the set Flatmates into the Worker type, which enables the use of the in operator.

include "alldifferent_except.mzn";

enum Flatmates = { Anne, Bert, Ceci, Dave };
enum Robots = { R2D2, C3PO, Marvin };

% Higher number = higher preference
array[Flatmates,Chores] of int: preference =
  [| 1, 2, 3, 4, 5 
   | 2, 1, 3, 4, 5 
   | 3, 5, 4, 1, 2 
   | 1, 5, 4, 2, 3 |];

array[Robots,Chores] of int: benefit =
  [| 20, 100, 20, 100, 30 
   | 10, 120, 40, 40,  60  
   | 50, 500, 30, 10,  70 |];

enum Chores = { Cooking, Vacuuming, Bathroom, Kitchen, Rubbish };

enum Workers = F(Flatmates) ++ R(Robots);

enum ChoreOrNothing = C(Chores) ++ { Nothing };

% Create preference and benefit arrays with added values for Nothing

array[Flatmates,ChoreOrNothing] of int: prefWithNothing =
  array2d(Flatmates,ChoreOrNothing,
          [ if c=Nothing then 0 else preference[f,C^-1(c)] endif
          | f in Flatmates, c in ChoreOrNothing]);

array[Robots,ChoreOrNothing] of int: benefitWithNothing =
  array2d(Robots,ChoreOrNothing,
          [ if c=Nothing then 0 else benefit[r,C^-1(c)] endif
          | r in Robots, c in ChoreOrNothing]);

array[Workers] of var ChoreOrNothing: assignment;

constraint alldifferent_except(assignment, {Nothing});

var int: overall_preferences =
  sum(f in Flatmates)(prefWithNothing[f, assignment[F(f)]]);

var int: overall_benefit =
  sum(r in Robots)(benefitWithNothing[r, assignment[R(r)]]);

solve maximize overall_benefit + overall_preferences;

output [ if w in F(Flatmates) then show(F^-1(w)) else show(R^-1(w)) endif++":\t"++
         if fix(assignment[w])=Nothing then "Nothing\n"
         else show(C^-1(assignment[w]))++"\n" endif
       | w in Workers];

2.2.5. Complex Constraints

Constraints are the core of the MiniZinc model. We have seen simple relational expressions but constraints can be considerably more powerful than this. A constraint is allowed to be any Boolean expression. Imagine a scheduling problem in which we have two tasks that cannot overlap in time. If s1 and s2 are the corresponding start times and d1 and d2 are the corresponding durations we can express this as:

constraint s1 + d1 <= s2  \/ s2 + d2 <= s1;

which ensures that the tasks do not overlap.

enum JOB;
enum TASK;
TASK: last = max(TASK);
array [JOB,TASK] of int: d;                   % task durations
int: total = sum(i in JOB, j in TASK)(d[i,j]);% total duration
int: digs = ceil(log(10.0,int2float(total))); % digits for output
array [JOB,TASK] of var 0..total: s;          % start times
var 0..total: end;                            % total end time


constraint %% ensure the tasks occur in sequence
    forall(i in JOB) (
        forall(j in TASK where j < last) 
            (s[i,j] + d[i,j] <= s[i,enum_next(TASK,j)]) /\
        s[i,last] + d[i,last] <= end
    );

constraint %% ensure no overlap of tasks
    forall(j in TASK) (
        forall(i,k in JOB where i < k) (
            s[i,j] + d[i,j] <= s[k,j] \/
            s[k,j] + d[k,j] <= s[i,j]
        )
    );

solve minimize end;

output ["end = \(end)\n"] ++
       [ show_int(digs,s[i,j]) ++ " " ++ 
         if j == last then "\n" else "" endif |
         i in JOB, j in TASK ];

JOB = _(1..5);
TASK = _(1..5);
d = [| 1, 4, 5, 3, 6 
     | 3, 2, 7, 1, 2
     | 4, 4, 4, 4, 4  
     | 1, 1, 1, 6, 8
     | 7, 3, 2, 2, 1 |];

The job shop scheduling model given in Listing 2.2.10 gives a realistic example of the use of this disjunctive modelling capability. In job shop scheduling we have a set of jobs, each consisting of a sequence of tasks on separate machines: so task [i,j] is the task in the \(i^{th}\) job performed on the \(j^{th}\) machine. Each sequence of tasks must be completed in order, and no two tasks on the same machine can overlap in time. Even for small instances of this problem it can be quite challenging to find optimal solutions.

The command

$ minizinc --solver gecode --intermediate jobshop.mzn jdata.dzn

solves a small job shop scheduling problem, and illustrates the behaviour of --intermediate for optimisation problems (note that when running this in the MiniZinc IDE, --intermediate is the default behaviour for optimisation problems). Here the solver outputs increasingly better solutions as it finds them. This flag has no effect for satisfaction problems. The output from this command is:

end = 39
 5  9 13 22 30
 6 13 18 25 36
 0  4  8 12 16
 4  8 12 16 22
 9 16 25 27 38
----------
end = 37
 4  8 12 17 20
 5 13 18 26 34
 0  4  8 12 16
 8 12 17 20 26
 9 16 25 27 36
----------
end = 34
 0  1  5 10 13
 6 10 15 23 31
 2  6 11 19 27
 1  5 10 13 19
 9 16 22 24 33
----------
end = 30
 5  9 13 18 21
 6 13 18 25 27
 1  5  9 13 17
 0  1  2  3  9
 9 16 25 27 29
----------
==========

indicating an optimal solution with end time 30 is finally found, and proved optimal. We can generate all optimal solutions by adding a constraint that end = 30 and changing the solve item to solve satisfy and then executing

$ minizinc --solver gecode --all-solutions jobshop.mzn jdata.dzn

For this problem there are 3,444,375 optimal solutions. In the MiniZinc IDE, you can select “User-defined behavior” in the configuration editor and uncheck the option “Stop after this many solutions” in order to display all solutions of this satisfaction problem.

2.2.6. Array Access Constraints

Another powerful modelling feature in MiniZinc is that decision variables can be used for array access. As an example, consider the (old-fashioned) stable marriage problem. We have n (straight) women and n (straight) men. Each man has a ranked list of women and vice versa. We want to find a husband/wife for each women/man so that all marriages are stable in the sense that:

• whenever m prefers another women o to his wife w, o prefers her husband to m, and

• whenever w prefers another man o to her husband m, o prefers his wife to w.

This can be elegantly modelled in MiniZinc. The model and sample data is shown in Listing 2.2.12 and Listing 2.2.13.

int: n;

enum Men = _(1..n);
enum Women = _(1..n);

array[Women, Men] of int: rankWomen;
array[Men, Women] of int: rankMen;

array[Men] of var Women: wife;
array[Women] of var Men: husband;

% assignment
constraint forall (m in Men) (husband[wife[m]]=m);
constraint forall (w in Women) (wife[husband[w]]=w);
% ranking
constraint forall (m in Men, o in Women) (
     rankMen[m,o] < rankMen[m,wife[m]] -> 
         rankWomen[o,husband[o]] < rankWomen[o,m] );

constraint forall (w in Women, o in Men) (
     rankWomen[w,o] < rankWomen[w,husband[w]] -> 
         rankMen[o,wife[o]] < rankMen[o,w] );
solve satisfy;

output ["wives= \(wife)\nhusbands= \(husband)\n"];

n = 5;
rankWomen = 
 [| 1, 2, 4, 3, 5,
  | 3, 5, 1, 2, 4,
  | 5, 4, 2, 1, 3,
  | 1, 3, 5, 4, 2,
  | 4, 2, 3, 5, 1 |];
  
rankMen = 
 [| 5, 1, 2, 4, 3,
  | 4, 1, 3, 2, 5,
  | 5, 3, 2, 4, 1,
  | 1, 5, 4, 3, 2,
  | 4, 3, 2, 1, 5 |];

The first three items in the model declare the number of men/women and the set of men and women. Here we introduce the use of anonymous enumerated types. Both Men and Women are sets of size n, but we do not wish to mix them up so we use an anonymous enumerated type. This allows MiniZinc to detect modelling errors where we use Men for Women or vice versa.

The matrices rankWomen and rankMen, respectively, give the women’s ranking of the men and the men’s ranking of the women. Thus, the entry rankWomen[w,m] gives the ranking by woman w of man m. The lower the number in the ranking, the more the man or women is preferred.

There are two arrays of decision variables: wife and husband. These, respectively, contain the wife of each man and the husband of each women.

The first two constraints

constraint forall (m in Men) (husband[wife[m]]=m);
constraint forall (w in Women) (wife[husband[w]]=w);

ensure that the assignment of husbands and wives is consistent: w is the wife of m implies m is the husband of w and vice versa. Notice how in husband[wife[m]] the index expression wife[m] is a decision variable, not a parameter.

The next two constraints are a direct encoding of the stability condition:

constraint forall (m in Men, o in Women) (
     rankMen[m,o] < rankMen[m,wife[m]] -> 
         rankWomen[o,husband[o]] < rankWomen[o,m] );

constraint forall (w in Women, o in Men) (
     rankWomen[w,o] < rankWomen[w,husband[w]] -> 
         rankMen[o,wife[o]] < rankMen[o,w] );

This natural modelling of the stable marriage problem is made possible by the ability to use decision variables as array indices and to construct constraints using the standard Boolean connectives. The alert reader may be wondering at this stage, what happens if the array index variable takes a value that is outside the index set of the array. MiniZinc treats this as failure: an array access a[e] implicitly adds the constraint e in index_set(a) to the closest surrounding Boolean context where index_set(a) gives the index set of a.

Thus for example, consider the variable declarations

array[1..2] of int: a= [2,3];
var 0..2: x;
var 2..3: y;

The constraint a[x] = y will succeed with \(x=1 \wedge y=2\) and \(x=2 \wedge y=3\). And the constraint not a[x] = y will succeed with \(x=0 \wedge y=2\), \(x=0 \wedge y=3\), \(x=1 \wedge y=3\) and \(x=2 \wedge y=2\).

In the case of invalid array accesses by a parameter, the formal semantics of MiniZinc treats this as failure so as to ensure that the treatment of parameters and decision variables is consistent, but a warning is issued since it is almost always an error.

int: n;
array[0..n-1] of var 0..n: s;

constraint forall(i in 0..n-1) (
   s[i] = (sum(j in 0..n-1)(bool2int(s[j]=i))));

solve satisfy;   

output [ "s = \(s);\n" ] ;

The coercion function bool2int can be called with any Boolean expression. This allows the MiniZinc modeller to use so called higher order constraints. As a simple example consider the magic series problem: find a list of numbers \(s= [s_0,\ldots,s_{n-1}]\) such that \(s_i\) is the number of occurrences of \(i\) in \(s\). An example is \(s = [1,2,1,0]\).

A MiniZinc model for this problem is shown in Listing 2.2.14. The use of bool2int allows us to sum up the number of times the constraint s[j]=i is satisfied. Executing the command

$ minizinc --solver gecode --all-solutions magic-series.mzn -D "n=4;"

leads to the output

s = [1, 2, 1, 0];
----------
s = [2, 0, 2, 0];
----------
==========

indicating exactly two solutions to the problem (the effect of --all-solutions can be achieved in the MiniZinc IDE using the “User-defined behavior” option in the solver configuration pane).

Note that MiniZinc will automatically coerce Booleans to integers and integers to floats when required. We could replace the the constraint item in Listing 2.2.14 with

constraint forall(i in 0..n-1) (
   s[i] = (sum(j in 0..n-1)(s[j]=i)));

and get identical results, since the Boolean expression s[j] = i will be automatically coerced to an integer, effectively by the MiniZinc system automatically adding the missing bool2int.

2.2.7. Partiality and Relational Semantics

Since MiniZinc supports complex constraints, it is important to understand how partial functions are handled. For example, we would expect the following model to have the solution b=true; x=1, but also b=false; x=0 and b=false; x=1:

var bool: b;
var 0..10: x;
constraint b -> 10 div x = 10;

However, division is a partial function, and 10 div 0 is undefined. So how exactly does MiniZinc handle partiality?

MiniZinc implements the relational semantics, which means that any undefinedness “bubbles up” to the closest Boolean context and becomes false there. In the case above, the closest Boolean context of the division is the = operator, so the right hand side of the implication becomes false. This forces the left hand side to be false as well.

MiniZinc has a number of built-in partial functions. The most common one is probably array access, since an expression x[i] is undefined if i is not in the index set of x.

In order to “guard” against partiality, we can use conditional statements like if y=0 then 0 else x div y endif, or if i in index_set(x) then x[i] else 0 endif, which define a default value for the cases where the expression is undefined.

A shorthand notation for this is the default operator, which allows you to write (x div y) default 0 or x[i] default 0 instead.

2.2.8. Set Constraints

Another powerful modelling feature of MiniZinc is that it allows sets containing integers to be decision variables: this means that when the model is evaluated the solver will find which elements are in the set.

As a simple example, consider the 0/1 knapsack problem. This is a restricted form of the knapsack problem in which we can either choose to place the item in the knapsack or not. Each item has a weight and a profit and we want to find which choice of items leads to the maximum profit subject to the knapsack not being too full.

It is natural to model this in MiniZinc with a single decision variable: var set of ITEM: knapsack where ITEM is the set of possible items. If the arrays weight[i] and profit[i] respectively give the weight and profit of item i, and the maximum weight the knapsack can carry is given by capacity then a natural model is given in Listing 2.2.15.

enum ITEM;
int: capacity;

array[ITEM] of int: profits;
array[ITEM] of int: weights;

var set of ITEM: knapsack;

constraint sum (i in knapsack) (weights[i]) <= capacity;

solve maximize sum (i in knapsack) (profits[i]) ;

output ["knapsack = \(knapsack)\n"];

Notice that the var keyword comes before the set declaration indicating that the set itself is the decision variable. This contrasts with an array in which the var keyword qualifies the elements in the array rather than the array itself since the basic structure of the array is fixed, i.e. its index set.

include "partition_set.mzn";
int: weeks;    set of int: WEEK = 1..weeks; 
int: groups;   set of int: GROUP = 1..groups;
int: size;     set of int: SIZE = 1..size;
int: ngolfers = groups*size; 
set of int: GOLFER = 1..ngolfers;

array[WEEK,GROUP] of var set of GOLFER: Sched;

% constraints
constraint 
    forall (i in WEEK, j in GROUP) (
           card(Sched[i,j]) = size
        /\ forall (k in j+1..groups) (
                Sched[i,j] intersect Sched[i,k] = {}
           )
    ) /\
    forall (i in WEEK) (
          partition_set([Sched[i,j] | j in GROUP], GOLFER)
    ) /\
    forall (i in 1..weeks-1, j in i+1..weeks) (
        forall (x,y in GROUP) (
            card(Sched[i,x] intersect Sched[j,y]) <= 1
        )
    );
% symmetry
  constraint 
    % Fix the first week %
    forall (i in GROUP, j in SIZE) (
        ((i-1)*size + j) in Sched[1,i]
    ) /\
    % Fix first group of second week %
    forall (i in SIZE) (
        ((i-1)*size + 1) in Sched[2,1]
    ) /\
    % Fix first 'size' players
    forall (w in 2..weeks, p in SIZE) (
        p in Sched[w,p]
    );
  
solve satisfy;

output [ show(Sched[i,j]) ++ " " ++
         if j == groups then "\n" else "" endif |
         i in WEEK, j in GROUP ];

weeks = 4;
groups = 4;
size =3 ;

As a more complex example of set constraint consider the social golfers problem shown in Listing 2.2.16. The aim is to schedule a golf tournament over weeks using groups \(\times\) size golfers. Each week we have to schedule groups different groups each of size size. No two pairs of golfers should ever play in two groups.

The variables in the model are sets of golfers Sched[i,j] for the \(i^{th}\) week and j^{th} group.

The constraints shown in lines 11-32 first enforces an ordering on the first set in each week to remove symmetry in swapping weeks. Next they enforce an ordering on the sets in each week, and make each set have a cardinality of size. They then ensure that each week is a partition of the set of golfers using the global constraint partition_set. Finally the last constraint ensures that no two players play in two groups together (since the cardinality of the intersection of any two groups is at most 1).

There are also symmetry breaking initialisation constraints shown in lines 34-46: the first week is fixed to have all players in order; the second week is made up of the first players of each of the first groups in the first week; finally the model forces the first size players to appear in their corresponding group number for the remaining weeks.

Executing the command

$ minizinc --solver gecode social-golfers.mzn social-golfers.dzn

where the data file defines a problem with 4 weeks, with 4 groups of size 3 leads to the output

1..3 4..6 7..9 10..12
{1,4,7} {2,9,12} {3,5,10} {6,8,11}
{1,6,9} {2,8,10} {3,4,11} {5,7,12}
{1,5,8} {2,7,11} {3,6,12} {4,9,10}
----------

Notice hows sets which are ranges may be output in range format.

2.2.9. Putting it all together

We finish this section with a complex example illustrating most of the features introduced in this chapter including enumerated types, complex constraints, global constraints, and complex output.

enum Guests = { bride, groom, bestman, bridesmaid, bob, carol, 
  ted, alice, ron, rona, ed, clara}; 
set of int: Seats = 1..12;
set of int: Hatreds = 1..5;
array[Hatreds] of Guests: h1 = [groom, carol, ed, bride, ted];
array[Hatreds] of Guests: h2 = [clara, bestman, ted, alice, ron];
set of Guests: Males = {groom, bestman, bob, ted, ron,ed};
set of Guests: Females = {bride,bridesmaid,carol,alice,rona,clara}; 

array[Guests] of var Seats: pos; % seat of guest
array[Hatreds] of var Seats: p1;  % seat of guest 1 in hatred
array[Hatreds] of var Seats: p2;  % seat of guest 2 in hatred 
array[Hatreds] of var 0..1: sameside;  % seats of hatred on same side 
array[Hatreds] of var Seats: cost;  % penalty of hatred  

include "alldifferent.mzn";
constraint alldifferent(pos);
constraint forall(g in Males)( pos[g] mod 2 == 1 );
constraint forall(g in Females)( pos[g] mod 2 == 0 );
constraint not (pos[ed] in {1,6,7,12});
constraint abs(pos[bride] - pos[groom]) <= 1 /\ 
           (pos[bride] <= 6 <-> pos[groom] <= 6);
constraint forall(h in Hatreds)(
               p1[h] = pos[h1[h]] /\
               p2[h] = pos[h2[h]] /\
               sameside[h] = bool2int(p1[h] <= 6 <-> p2[h] <= 6) /\
               cost[h] = sameside[h] * abs(p1[h] - p2[h]) +
                  (1 - sameside[h]) * (abs(13 - p1[h] - p2[h]) + 1) );

solve maximize sum(h in Hatreds)(cost[h]);

output [ show(g)++" " | s in Seats,g in Guests where fix(pos[g]) == s]
       ++ ["\n"];