Decompositions for Large Shift Scheduling Problem & Model Lifecycle

Decompositions for Large Shift Scheduling Problem & Model Lifecycle

  • Post author:
  • Post category:Blog

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

Decompositions for shift scheduling

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)

demand in shifts without assigning the shifts to employees

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("")
    }
}