Decompositions for Large Shift Scheduling Problem & Model Lifecycle

When shift scheduling problems become too large to be solved with a single model, an option often chosen by practitioners is to decompose them.

Shift-scheduling MIP models contain a very large amount of boolean variables which makes them hard to solve. The idea behind model decomposition is to divide the model in two phases

a model to select the shifts to cover the demand
a model assign employees to selected shift

The advantage of the decomposition is to reduce the total number of variables in the model. If there are E employees and S possible shifts, the single MIP model would have $E \times S$ boolean variables. While the decomposed models would have S integer and $c \times E$ boolean variables respectively (with $c \ll S$ )

The first phase of the decomposition covers the demand with shifts without assigning the shifts to employees

Reduction of variables by summation

When implementing a decomposition for an optimisation problem, we recommend to always start from a complete model which will serve as guide, benchmark and validity checker.

Let’s start from the shift-scheduling model built in the last post

 $\begin{aligned} \min \mathrm{lex}\, (z, cost, under)\\ \\ \forall t \quad D_t - \left \lbrack \sum_e \sum_{t\, \in\, s} x^e_s \right \rbrack - 1 \leq z & \quad (1)\\ cost = \sum^{}_e \sum_s \mathrm{duration}_s * x^e_s & \quad (2)\\ under = \sum_t u_t & \quad (3)\\ \forall t \quad D_t = \left \lbrack \sum_e \sum_{t\, \in\, s}x^e_s \right \rbrack + u_t - o_t & \quad(4)\\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e & \quad (5)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e & \quad (6)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ z, u_t,\ o_t \in \mathbb{R}^+ &\\ \end{aligned}$

Constraint (4) computes the over and under coverage for each time point t. Constraints (5) and (6) fix the number of workdays and work hours per employee. Constraint (7) says an employee cannot be assigned to shifts that are apart by less than $\underline{d}$ . Constraints (1), (2) and (3) define the objective functions max under-coverage, cost and total under-coverage.

We want to separate the selection of the shifts and the assignment to an employee. In other words we want to break the variable $x^e_s$ into two parts

one that captures the fact that a shift is selected for assignment to any employee
one that selects the specific employee for an already selected shift

This can be achieved by summation : introducing a variable $n_s$ that is the sum of $x^e_s$ over all employees
$\forall s \quad n_s = \sum_e x^e_s$

Now $n_s \in \mathbb{N}$ counts the number of times shift s has been assigned to an employee, without capturing which employees it is. As such it replaces E boolean variables by one integer variable, reducing the total number of variables in the problem from $E \times S$ booleans to a first phase with S integers.

Reference model with double set of variables

We can now rewrite the model. However, we will not replace $x^e_s$ with $n_s$ yet. Instead, every time possible, we will duplicate each constraint with one version using $x^e_s$ and one version using $n_s$ .

Replacing $x^e_s$ with $n_s$ is possible every time a constraint contains the term $\sum_e x^e_s$ . When this is not the case but the constraint starts with $\forall e$ , we can make the term $\sum_e x^e_s$ appear by summation of the constraints over e

 $\begin{aligned} & \forall e \forall a \dots x^e_s\\\\ \Rightarrow \ & \sum_e \forall e \forall a \ \dots \ x^e_s\\ \leadsto &\ \forall a\ \dots\ \sum_e s^e_s\\ \leadsto &\ \forall a\ \dots\ n_s \end{aligned}$

The result on the shift-scheduling model is

 $\begin{aligned} \min \mathrm{lex}\, (z, cost, under)\\ \\ \forall t \quad D_t - \left \lbrack \sum_e \sum_{t\, \in\, s} x^e_s \right \rbrack - 1 \leq z & \quad (1)\\ \forall t \quad D_t - \sum_{t\, \in\, s} n_s - 1 \leq z & \quad (1\ bis)\\ cost = \sum^{}_e \sum_s \mathrm{duration}_s * x^e_s & \quad (2)\\ cost = \sum_s \mathrm{duration}_s * n_s & \quad (2\ bis)\\ under = \sum_t u_t & \quad (3)\\ \forall t \quad D_t = \left \lbrack \sum_e \sum_{t\, \in\, s}x^e_s \right \rbrack + u_t - o_t & \quad(4)\\ \forall t \quad D_t = \sum_{t\, \in\, s} n_s + u_t - o_t & \quad(4\ bis)\\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e & \quad (5\forall)\\ \sum^{}_s n_s = \sum_e \mathrm{Days}^e & \quad (5\Sigma)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e & \quad (6\forall)\\ \sum_s \mathrm{duration}_s * n_s = \sum^{}_e \mathrm{Hours}^e & \quad (6\Sigma)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d} \ \Rightarrow \ n_s + n_u \leq \ E & \quad (7\Sigma)\\ \forall s \quad n_s = \sum^{}_e x^e_s & \quad (8)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ n_s \in \mathbb{N} &\\ z, u_t,\ o_t \in \mathbb{R}^+ &\\ \end{aligned}$

It is fundamental to understand that all bis constraints are equivalent to their original version, however the constraints $5\Sigma$ , $6\Sigma$ and $7\Sigma$ are less strong than their counterparts $5\forall$ , $6\forall$ and $7\forall$ .

We have $(\forall5) \Rightarrow (\Sigma5)$ but the opposite is not necessairly true. Which means that once we remove constraints $(\forall5)$ we may find solutions to the model that are not feasible for the whole problem.

Before cutting the model into two phases, let’s remove all constraints that have a bis version and reorganize the constraints

 $\begin{aligned} \min \mathrm{lex}\, (z, cost, under)\\ \\ \forall t \quad D_t - \sum_{t\, \in\, s} n_s - 1 \leq z & \quad (1\ bis)\\ cost = \sum_s \mathrm{duration}_s * n_s & \quad (2\ bis)\\ under = \sum_t u_t & \quad (3)\\ \forall t \quad D_t = \sum_{t\, \in\, s} n_s + u_t - o_t & \quad(4\ bis)\\ \sum^{}_s n_s = \sum_e \mathrm{Days}^e & \quad (5\Sigma)\\ \sum_s \mathrm{duration}_s * n_s = \sum^{}_e \mathrm{Hours}^e & \quad (6\Sigma)\\ \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d} \ \Rightarrow \ n_s + n_u \leq \ E & \quad (7\Sigma)\\ \\ \forall s \quad n_s = \sum_e x^e_s & \quad (8) \\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e & \quad (5\forall)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e & \quad (6\forall)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ n_s \in \mathbb{N} &\\ z, u_t,\ o_t \in \mathbb{R}^+ &\\ \end{aligned}$

This is our reference model. It is typically slower than the original model, but mathematically equivalent : same solutions, same optimal value. And the model is ready to be divided along constraint (8)

Two-phase decomposition

We can now run a first model with S integer variables $n_s$

 $\begin{aligned} \min \mathrm{lex}\, (z, cost, under)\\ \\ \forall t \quad D_t - \sum_{t\, \in\, s} n_s - 1 \leq z & \quad (1\ bis)\\ cost = \sum_s \mathrm{duration}_s * n_s & \quad (2\ bis)\\ under = \sum_t u_t & \quad (3)\\ \forall t \quad D_t = \sum_{t\, \in\, s} n_s + u_t - o_t & \quad(4\ bis)\\ \sum^{}_s n_s = \sum_e \mathrm{Days}^e & \quad (5\Sigma)\\ \sum_s \mathrm{duration}_s * n_s = \sum^{}_e \mathrm{Hours}^e & \quad (6\Sigma)\\ \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d} \ \Rightarrow \ n_s + n_u \leq \ E & \quad (7\Sigma)\\ \\ n_s \in \mathbb{N} &\\ z, u_t,\ o_t \in \mathbb{R}^+ &\\ \end{aligned}$

Once the solution of this model is known, the variables $n_s$ become constants for the next phase. We can now compute the set of selected shifts

 $S' = \left \lbrace\ s \in S \mid n_s \gt 0\ \right \rbrace$

In this case, we know per constraint $(5)$ $\sum_s n_s = \sum_e \mathrm{Days}^e$ an upper bound on the number of shifts in S’

We then solve a second model with $E \times S'$ boolean variables. In this case the model is a feasibility model. Remember $n_s$ are now constants, we will therefore write them $N_s$ instead

 $\begin{aligned} \forall s \quad \sum_e x^e_s = N_s & \quad (8)\\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e & \quad (5\forall)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e & \quad (6\forall)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ \end{aligned}$

Benefits and limitations of a decomposition

Typical benefits of decompositions are

Ability to solve larger problems
Ability to use different technologies in each problem

Notice how the second phase of the problem became a satisfaction problem instead of an optimization problem. We could solve it with integer programming (MIP), constraint programming (CP), local search (LS), satisfiability (SAT).

Typical limitations of decompositions are

The "optimal" solutions of each phase of the decomposition often do not result in an optimal solution for the whole problem when put together
Implementation is a little tedious due to limited support from modeling languages, and the need to be transferring data between the models
Decomposed models gain less from engine improvements than single models. As engines improve with time, decomposed models will see smaller gains. The most intrusive / complex is the decomposition, the more it negatively impacts the engine improvements
Decomposed models are harder to maintain and evolve, which may be a problem in environments where the person in charge often (external consultants, etc)

Infeasible sub-problems

An important issue in decompositions is infeasible second-phases (sub / slave problem). This results from the first phase (main / master problem) not having enough constraints.

Without constraints $5\Sigma$ , $6\Sigma$ and $7\Sigma$ , a first phase like this one is likely to generate an assignment of $n_s$ that is infeasible (for instance not have $\sum_s n_s = \sum_e \mathrm{Days}^e$ ) for the second phase.

The first model that comes to mind when writing the decomposition "directly" is under-constrained

 $\begin{aligned} \min \mathrm{lex}\, (z, cost, under)\\ \\ \forall t \quad D_t - \sum_{t\, \in\, s} n_s - 1 \leq z & \quad (1\ bis)\\ cost = \sum_s \mathrm{duration}_s * n_s & \quad (2\ bis)\\ under = \sum_t u_t & \quad (3)\\ \forall t \quad D_t = \sum_{t\, \in\, s} n_s + u_t - o_t & \quad(4\ bis)\\ \\ n_s \in \mathbb{N} &\\ z, u_t,\ o_t \in \mathbb{R}^+ &\\ \end{aligned}$

This is why we recommend to always create decompositions by summation from a full model. It is easy to overlook or forget constraints like $5\Sigma$ , $6\Sigma$ and $7\Sigma$ .

Soft sub-problems

Even if the model was obtained by summation, it can still lead to an infeasible second phase (remember $\Sigma$ constraints are weaker than $\forall$ constraints). Therefore tt is advisable to consider the second-phase as soft problem with soft-constraints

Deviation minimization

The following sub-problem minimizes the deviation between $\sum_e x^e_s$ and $n_s$

 $\begin{aligned} \min\ \sum_s \left \Vert \sum_e x^e_s - n_s \right \Vert\\ \\ \forall s \quad \sum_e x^e_s = N_s & \quad (8)\\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e & \quad (5\forall)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e & \quad (6\forall)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ \end{aligned}$

Partial usage of fixed variables

While this sub-problem variant allows not all shifts or employees to be used

 $\begin{aligned} \max \sum_e y^e\\ \\ \forall s \quad \sum_e x^e_s \leq N_s & \quad (8\ soft)\\ \forall e\ \quad \sum_s x^e_s = y^e & \quad (9)\\ \forall e \quad \sum^{}_s x^e_s = \mathrm{Days}^e * y^e & \quad (5\forall\ soft)\\ \forall e \quad \sum_s \mathrm{duration}_s * x^e_s = \mathrm{Hours}^e * y^e & \quad (6\forall\ soft)\\ \forall e\, \forall s, u \quad \mathrm{dist}(s, u) \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \\ x^e_s \in \lbrace 0, 1 \rbrace &\\ y^e \in \lbrace 0, 1 \rbrace &\\ \end{aligned}$

Continuous variables in main problem

A possible approach to reduce infeasibility issues is to not remove variables $x^e_s$ from the first phase completely, but transform them from boolean to floating point variables.

$x^e_s \in \lbrack \, 0 \dots 1 \, \rbrack$

The reduction in memory consumption is lower, but some of the pruning power of $\forall$ constraints is conserved.

Benders decomposition

Benders decomposition is a decomposition technique where when the second phase doesn’t have a solution, new constraints for the first phase are created, and the first phase is solved again. The loops continues till optimality.

Model lifecycle

Working with decompositions requires maintaining the full model and the decomposed model in sync at all times. This may be tedious and one may be tempted to trash the reference model as soon as the decomposition works. But the reference model is crucial to fix bugs and adjust performance of the decomposition.

We will show illustrate why with the minimum distance constraint

Constraints $(\forall 7)$ and $(\Sigma7)$ generate a large number of rows in the MIP

 $\begin{aligned} \forall e\, \forall s, u \quad t_s < t_u \wedge t_u + - t_s - d_s \lt \underline{d}\ \Rightarrow \ x^e_s + x^e_u \leq 1 & \quad (7\forall)\\ \forall s, u \quad t_s < t_u \wedge t_u - t_s - d_s \lt \underline{d}\ \Rightarrow \ n_s + n_u \leq \ E & \quad (7\Sigma)\\ \end{aligned}$

If we take $\underline{d} = 12$ hours between shifts, shifts ranging from 6h to 10h and a time granularity of 15 minutes over a 1 week horizon we have

(7 24 4) = 672 choices for the starting time of s
(12 * 4) = 48 choices for the starting time of t
less than ((10 – 6) * 4) choices for the duration of s
((10 – 6) * 4) = 16 choices for the duration of t

In total $E * 8M$ constraints for $(7\forall )$ and $8M$ constraints for $(7\Sigma)$ . The constraint needs to be rewritten in the original model and re-aggregated.

Because an employee cannot have two shifts in the same day d, we can separate the same day inequality and the next day minimum distance. And then simplify the constraint to

 $\begin{aligned} \forall e\, \forall d \quad \sum_{s \in d} x^e_s \leq 1 & \quad (7a) \\ \forall e\, \forall d\, \forall s \in d \quad x^e_s + \sum_{u \in d + 1 \wedge d_{us} \leq \underline{d}} s^e_u \leq 1 & \quad (7b) \end{aligned}$

Where the sum in (7b) is restricted to the shifts u that violate the minimum distance to s.

This reduces the number of constraints to one (7a) per day and one (7b) per shift and employee ; with the same data as before, that’s 10K constraints.

After summation we get

 $\begin{aligned} \forall d \quad \sum_{s \in d} n_s \leq E & \quad (7a \Sigma) \\ \forall d\, \forall s \in d \quad n_s + \sum_{u \in d + 1 \wedge d_{us} \leq \underline{d}} n_u \leq E & \quad (7b \Sigma)\\ \end{aligned}$

The process would be the same for any change to the model like adding contract selection.

OPL code

Writing decompositions in OPL is not very practical because it requires writing control flow in OPL script which is a little cumbersome. Instead we are going to do a the following

the first model writes a .dat file with the selected shifts
the second model reads that .dat

Random demand is generated in the models to make them self-contained

Results

When ran on an Intel core i5 @1.60Ghz with 8MB of ram

the reference model ran out of memory and never returned any result
the decomposed model ran in 1 minute overall

When ran on instances small enough for both models to run without problems, the reference model and the two-phase decompositions seem to find the same quality of solution.

Reference model

/*********************************************
 * OPL 20.1.0.0 Model
 * Author: diegoolivierfernandezpons
 * Creation Date: Jan 15, 2021 at 4:51:52 PM
 *********************************************/

 int ticks_per_hour = 4;
 int undercoverage_tolerance = 1;
 int inter_shift = 12;

 range days = 0 .. 6;
 range hours = 6 * ticks_per_hour .. 24 * ticks_per_hour;

 int demand [d in days][h in hours] = rand(10);

 {string} employees = { "Alice", "Bob", "Charlie", "Diana", "Eleonor", "Frida", "Gilles", "Helena", "Irina", "Jasmin", "Keith", "Lily", "Mathias", "Nathalie", "Olga", "Patricia", "Quinn",  "Ralf", "Samatha", "Taylor", "Ulrich", "Vanessa", "William", "Xavier", "Yannis", "Zoe" };

 int contract_days = 5;
 int contract_hours = 40 * ticks_per_hour;
 int min_hours = 5 * ticks_per_hour;
 int max_hours = 10 * ticks_per_hour;
 range duration = min_hours .. max_hours;

 tuple shift_t { int start; int duration; }
 {shift_t} shifts = { < s, d > | s in hours, d in duration : s + d - 1 in hours };

 {int} md [s in shifts] = { u.start | u in shifts : s.start + s.duration + inter_shift * ticks_per_hour > 24 * ticks_per_hour + u.start };

 dvar boolean use [employees][days][shifts];
 dvar int+ open [days][shifts];

 dvar float+ over  [days][hours];
 dvar float+ under [days][hours];
 dvar float+ z;

 dexpr float undercoverage = sum (d in days, h in hours) under[d][h];
 dexpr float overcoverage = sum (d in days, h in hours) over[d][h];
 dexpr int headcount = card(employees);

 minimize staticLex (z, undercoverage, overcoverage);

 constraints {

    forall (d in days, s in shifts) open[d][s] == sum (e in employees) use[e][d][s];

    // Undercoverage threshold
    forall (d in days, h in hours)
        under[d][h] - undercoverage_tolerance <= z;

    // Compute under and over coverage
    forall (d in days, h in hours) 
        sum (s in shifts : s.start <= h && h  < s.start + s.duration) open[d][s] 
        == demand[d][h] + over[d][h] - under[d][h];

    // Number of days in the contract
    sum (d in days, s in shifts) open[d][s] == headcount * contract_days;
    forall (e in employees) sum (d in days, s in shifts) use[e][d][s] == contract_days; 

    // Number of hours in the contract
    sum (d in days, s in shifts) s.duration * open[d][s] == headcount * contract_hours;
    forall (e in employees) sum (d in days, s in shifts) s.duration * use[e][d][s] == contract_hours; 

    // At most one shift per day per employee
    forall (d in days) sum (s in shifts) open[d][s] <= headcount;
    forall (e in employees, d in days) sum (s in shifts) use [e][d][s] <= 1;

    // Minimum distance between shifts
    //forall (e in employees, d in days : d + 1 in days, s in shifts : card(md[s]) > 0)
    //  use[e][d][s] + sum (t in md[s], u in shifts : u.start == t) use[e][d + 1][u] <= 1;

    // Minimum distance between shifts
    forall (d in days : d + 1 in days, s in shifts : card(md[s]) > 0)
          open[d][s] 
        + sum (t in md[s], u in shifts : u.start == t) open[d + 1][u] 
        <= headcount;
 }

 tuple sol_shift { int start; int duration; }
 sorted {sol_shift} solution [e in employees] = { <24 * d * ticks_per_hour + s.start, s.duration> | s in shifts, d in days : use[e][d][s] == 1 };
 range horizon = 0 .. max (d in days) (d + 1) * 24 * ticks_per_hour ;
 int present [h in horizon][e in employees] = sum (<s, u> in solution[e] : s <= h && h < s + u) 1;

 execute { 
    for (var e in employees) {
        write (e.charAt(0), " ")
        for (var d in days) {
            for (var h in hours)
                if (present[d * 24 * ticks_per_hour + h][e]) write ("X"); else write (" ")
            write ("|")
        }
        writeln("")
    }
}

First phase

/*********************************************
 * OPL 20.1.0.0 Model
 * Author: diegoolivierfernandezpons
 * Creation Date: Jan 15, 2021 at 4:51:52 PM
 *********************************************/

 int ticks_per_hour = 4;
 int undercoverage_tolerance = 1;
 int inter_shift = 12;

 range days = 0 .. 6;
 range hours = 6 * ticks_per_hour .. 24 * ticks_per_hour;

 int demand [d in days][h in hours] = rand(10);

 {string} employees = { "Alice", "Bob", "Charlie", "Diana", "Eleonor", "Frida", "Gilles", "Helena", "Irina", "Jasmin", "Keith", "Lily", "Mathias", "Nathalie", "Olga", "Patricia", "Quinn",  "Ralf", "Samatha", "Taylor", "Ulrich", "Vanessa", "William", "Xavier", "Yannis", "Zoe" };

 int contract_days = 5;
 int contract_hours = 40 * ticks_per_hour;
 int min_hours = 5 * ticks_per_hour;
 int max_hours = 10 * ticks_per_hour;
 range duration = min_hours .. max_hours;

 tuple shift_t { int start; int duration; }
 {shift_t} shifts = { < s, d > | s in hours, d in duration : s + d - 1 in hours };

 {int} incompatible [s in shifts] = { u.start | u in shifts : s.start + s.duration + inter_shift * ticks_per_hour > 24 * ticks_per_hour + u.start };

 dvar int+ open [days][shifts];

 dvar float+ over  [days][hours];
 dvar float+ under [days][hours];
 dvar float+ z;

 dexpr float undercoverage = sum (d in days, h in hours) under[d][h];
 dexpr float overcoverage = sum (d in days, h in hours) over[d][h];
 dexpr int headcount = card(employees);

 minimize staticLex (z, undercoverage, overcoverage);

 constraints {

    // Undercoverage threshold
    forall (d in days, h in hours)
        under[d][h] - undercoverage_tolerance <= z;

    // Compute under and over coverage
    forall (d in days, h in hours) 
        sum (s in shifts : s.start <= h && h  < s.start + s.duration) open[d][s] 
        == demand[d][h] + over[d][h] - under[d][h];

    // Number of days in the contract
    sum (d in days, s in shifts) open[d][s] == headcount * contract_days;

    // Number of hours in the contract
    sum (d in days, s in shifts) s.duration * open[d][s] == headcount * contract_hours;

    // At most one shift per day per employee
    forall (d in days) sum (s in shifts) open[d][s] <= headcount;

    // Minimum distance between shifts
    forall (d in days : d + 1 in days, s in shifts : card(incompatible[s]) > 0)
          open[d][s] 
        + sum (t in incompatible[s], u in shifts : u.start == t) open[d + 1][u] 
        <= headcount;
 }

 tuple shift { int day; int start; int duration; int quantity; }
 sorted {shift} open_shifts = { <d, s.start, s.duration, open[d][s]> | s in shifts, d in days : open[d][s] > 0 };
 int q = sum (d in days, s in shifts) open[d][s];

 execute {
    var o = new IloOplOutputFile("shifts.dat");
    o.writeln ("// ", open_shifts.size, " shifts types")
    o.writeln ("// ", q, " open shifts")
    o.writeln("shifts = ", open_shifts, ";")    
    o.close()
  }

Second phase

/*********************************************
 * OPL 20.1.0.0 Model
 * Author: diegoolivierfernandezpons
 * Creation Date: Jan 15, 2021 at 4:51:52 PM
 *********************************************/

 int ticks_per_hour = 4;
 int undercoverage_tolerance = 1;
 int inter_shift = 12;

 range days = 0 .. 6;
 range hours = 6 * ticks_per_hour .. 24 * ticks_per_hour;

 int demand [d in days][h in hours] = rand(10);

 {string} employees = { "Alice", "Bob", "Charlie", "Diana", "Eleonor", "Frida", "Gilles", "Helena", "Irina", "Jasmin", "Keith", "Lily", "Mathias", "Nathalie", "Olga", "Patricia", "Quinn",  "Ralf", "Samatha", "Taylor", "Ulrich", "Vanessa", "William", "Xavier", "Yannis", "Zoe" };

 int contract_days = 5;
 int contract_hours = 40 * ticks_per_hour;
 int min_hours = 5 * ticks_per_hour;
 int max_hours = 10 * ticks_per_hour;
 range duration = min_hours .. max_hours;

 tuple shift_t { int day; int start; int duration; int quantity; }
 {shift_t} shifts = ...;

 {int} md [s in shifts] = { u.start | u in shifts : s.start + s.duration + inter_shift * ticks_per_hour > 24 * ticks_per_hour + u.start };

 dvar boolean use [employees][shifts];

 constraints {

    forall (d in days, s in shifts) sum (e in employees) use[e][s] == s.quantity;

    // Number of days in the contract
    forall (e in employees) sum (s in shifts) use[e][s] == contract_days; 

    // Number of hours in the contract
    forall (e in employees) sum (s in shifts) s.duration * use[e][s] == contract_hours; 

    // At most one shift per day per employee
    forall (e in employees, d in days) sum (s in shifts : s.day == d) use[e][s] <= 1;

    // Minimum distance between shifts
    forall (e in employees, s in shifts : card(md[s]) > 0)
        use[e][s] + sum (t in md[s], u in shifts : u.start == t && u.day == s.day + 1) use[e][u] <= 1;
 }

 tuple sol_shift { int start; int duration; }
 sorted {sol_shift} solution [e in employees] = { <24 * s.day * ticks_per_hour + s.start, s.duration> | s in shifts : use[e][s] == 1 };
 range horizon = 0 .. max (d in days) (d + 1) * 24 * ticks_per_hour ;
 int present [h in horizon][e in employees] = sum (<s, u> in solution[e] : s <= h && h < s + u) 1;

 execute { 
    for (var e in employees) {
        write (e.charAt(0), " ")
        for (var d in days) {
            for (var h in hours)
                if (present[d * 24 * ticks_per_hour + h][e]) write ("X"); else write (" ")
            write ("|")
        }
        writeln("")
    }
}