## Sunday, March 28, 2021

### Parallel machine scheduling I, two formulations

Scheduling problems with multiple machines are not that easy for MIP solvers. However, modern solvers are capable to solve reasonably sized problems to optimality quite efficiently. Here are some notes.

#### Data and solution

----     29 PARAMETER jobdata

proctime     release     duedate     weights

job1            4          61          70           1
job2            9          59          78           1
job3            7                      67           1
job4            5          21          35           1
job5            5           6          17           1
job6            4                      19           1
job7            5          57          68           1
job8            9          41          59           1
job9            3                      12           1
job10           7          67          82           1
job11          10                      16           1
job12           7           9          25           1
job13          10          40          57           1
job14           9          58          78           1
job15           4                      26           1
job16           8                      13           1
job17           4                      63           1
job18           5                      66           1
job19           8                      45           1
job20           6          25          42           1
job21           5                      32           1
job22           5                      32           1
job23           4           5          21           1
job24           4                      94           1
job25           7          22          44           1
job26           9          60          81           1
job27           4                      37           1
job28           8                      21           1
job29           9          61          78           1
job30           5           4          16           1
job31           3                      28           1
job32           7                      10           1
job33           4          24          34           1
job34           9          39          55           1
job35           5                      23           1
job36           5                      25           1
job37           7          24          40           1
job38           8                      38           1
job39           8                      39           1
job40           6                      97           1
job41           6                     100           1
job42           3          33          43           1
job43           5          28          43           1
job44           3          64          80           1
job45           5          31          46           1
job46           4                      93           1
job47           8           2          20           1
job48           7          59          76           1
job49           9                      15           1
job50           5          51          62           1

----     33 SET precedence  precedence constraints: i must execute before j

job4        job8        job9       job10       job11       job13       job16       job20       job33

job1          YES
job4                      YES
job6                                                          YES
job7                                              YES
job8                                  YES                     YES         YES
job11                                                                                 YES
job15                                                                                 YES
job18                                                                                             YES
job28                                                                                                         YES
job30                                                                                                         YES

+       job36       job37       job44

job33         YES         YES
job40                                 YES


Notes:
• Some jobs have a nonzero release date: they cannot start before a certain time.
• Some jobs also have some precedence relationships: they cannot start until some other job finishes. For instance: job4 needs to wait until job1 is completed.
• The weights on the jobs can be used used to make some jobs more important: they may be used in some objectives.
• We just indicate time here by integer values. In practice, time may be more like a calendar, with days off (weekends, observed holidays). In that case, you have to decide whether to use dates directly as the time index or map them first to integers. In general, I would prefer to use the dates directly, but in your case, that may be different.

Assume we have four identical machines at our disposal. Then a solution that minimizes the amount jobs are violating their due date (tardiness) can look like:

The picture is a bit messy. It tries a bit too hard to display all information. However, we can see there are a few jobs that are tardy:

----    184 VARIABLE start.L  start time of job

job1  61,    job2  67,    job3  40,    job4  65,    job5   8,    job6   1,    job7  57,    job8  70,    job9  79
job10 71,    job11 79,    job12 18,    job13 79,    job14 58,    job15  5,    job16 89,    job17 21,    job18 20
job19 36,    job20 25,    job21 25,    job22 19,    job23 14,    job24  4,    job25 35,    job26 69,    job27  1
job28 13,    job29 62,    job30  9,    job31  1,    job32  1,    job33 24,    job34 44,    job35 14,    job36 28
job37 33,    job38 25,    job39 31,    job40 42,    job41 44,    job42 33,    job43 30,    job44 66,    job45 39
job46  8,    job47 12,    job48 59,    job49  5,    job50 51

----    184 VARIABLE obj.L  individual objectives

completion 2096,    sumtardy    322,    maxtardy     84,    numtardy      7,    makespan     97

----    184 VARIABLE tardy.L  amount a job is tardy

job4  35,    job8  20,    job9  70,    job11 73,    job13 32,    job16 84,    job36  8


A different picture emerges when plotting machine occupancy:

The numbers (1 through 50) and colors indicate the jobs. Presented this way the problem looks much smaller! Of course, in practice, it is beneficial to display the results in a way that makes it easy for the user to interpret the results [2]. This may require some interactivity (e.g. clicking on a job to get more information about it).

#### Objectives

In the following models, I will consider a weighted sum of different objectives [1]:
• min sum weighted completion: $\min \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{completion}}_j$ pushes all jobs towards early completion
• min sum tardiness: $\min \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{tardy}}_j$ tries to minimize the total amount each job being tardy (normal people would call this "late", but "late" has a slightly different meaning in scheduling: late can be negative when early, so minimizing lateness is not a good idea).
• min max tardiness: $\min \max \color{darkred}{\mathit{tardy}}_j$ tries to prevent jobs that are very tardy. (Sometimes min max lateness is used. This is mostly the same, except when all jobs are completed before their due date.)
• min num tardy jobs. Minimize the number of tardy jobs. This can help prevent many jobs from being just a little bit tardy.
• min makespan. This is minimizing the completion time of the last job.

These are the most common objectives. They can be used individually, or they can be combined. In some practical models I have been working on, clients want to have rather complex objectives. Some additional considerations I have seen:

• Make the completion time as close a possible to the due date. This is important when the final product will deteriorate in quality (e.g. fresh produce) or when storage is an issue.
• Another objective is to have as much usable empty space in the schedule. I.e. not a lot of small fragments of unoccupied machines, but rather big blocks.
• Different orders can originate from the same client. That may have some impact. E.g. when we want to minimize the total tardiness per client.

#### Model 1: Time indexed model

Discrete time can be implemented using a time index indicating in which time slot a job is starting. In models of this type we use variables of the form: $\color{darkred}x_{j,m,t} = \begin{cases}1 & \text{if job j starts at time t on machine m} \\ 0 & \text{otherwise}\end{cases}$

We need to be a bit precise in our definition of time. When $$\color{darkred}x_{j,m,t}$$, job $$j$$ occupies machine $$m$$ during periods $$t,t+1.\dots,t+\color{darkblue}{\mathit{proctime}}_j-1$$. The release dates $$\color{darkblue}{\mathit{release}}_j=t$$ is measured at the beginning of the period. So if the release date is $$t$$ we can start working on job $$j$$ during period $$t$$. Similarly, the completion time and the due date $$\color{darkblue}{\mathit{due}}_j$$ are measured at the beginning of a period.

You may use a different convention, but it is important to document it and stick to one.

The most complicated constraint is the no-overlap constraint: a machine can only work on one job at a time. If we look at period $$t$$, any job $$j$$ occupies machine $$m$$ at $$t$$ if $$\color{darkred}x_{j,m,t'}=1$$ if $$t' \in \{t-\color{darkblue}{\mathit{proctime}}_j+1,\dots,t\}$$. This leads to a constraint: $\sum_j \sum_{t'=t-\color{darkblue}{\mathit{proctime}}_j+1}^{t} \color{darkred}x_{j,m,t'} \le 1 \>\>\forall t,m$

This is a bit of a handful. I like to split this into two. First calculate a set $$\color{darkblue}{\mathit{occupy}}_{j,t,t'}$$ indicating whether a job $$j$$ starting at $$t$$ will occupy a machine at $$t'$$. This set can look like:

With this, we can simplify our constraint to: $\sum_{j,t'|\color{darkblue}{\mathit{occupy}}(j,t',t)} \color{darkred}x_{j,m,t'}\le 1 \>\>\forall m,t$ The advantage of this is that we can develop and debug $$\color{darkblue}{\mathit{occupy}}(j,t,t')$$ in isolation, and before the model is running. Note that we can achieve the same thing with a binary parameter. The constraint would then look like: $\sum_{j,t'} \color{darkblue}{\mathit{occupy}}_{j,t',t} \cdot \color{darkred}x_{j,m,t'}\le 1 \>\>\forall m,t$
As equations are the most difficult to debug, I prefer to keep tham as simple as possible. One way to achieve that is to precalculate things as much as possible.

The complete model can look like:

Model M1: Time indexed
\begin{align}\min&\sum_{\mathit{objs}} \color{darkblue}{\mathit{objw}}_{\mathit{objs}} \cdot \color{darkred}{\mathit{obj}}_{\mathit{objs}} && \\ \hline & \color{darkred}{\mathit{obj}}_{\mathit{completion}} = \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{completion}}_j && \\ & \color{darkred}{\mathit{obj}}_{\mathit{sumtardy}} = \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{tardy}}_j \\ & \color{darkred}{\mathit{obj}}_{\mathit{numtardy}} = \sum_{j,m,t|\color{darkblue}{\mathit{isTardy}}(j,t)}\color{darkred}x_{j,m,t} \\ & \color{darkred}{\mathit{obj}}_{\mathit{maxtardy}} \ge \color{darkred}{\mathit{tardy}}_j && \forall j \\ & \color{darkred}{\mathit{obj}}_{\mathit{makespan}} \ge \color{darkred}{\mathit{completion}}_j && \forall j\\ \hline & \sum_{m,t} \color{darkred}x_{j,m,t}=1 && \forall j \\ & \sum_{j,t|\color{darkblue}{\mathit{occupy}}(j,t',t)} \color{darkred}x_{j,m,t'}\le 1 && \forall m,t \\ & \color{darkred}{\mathit{start}}_j = \sum_{m,t} t \cdot \color{darkred}x_{j,m,t} && \forall j \\ & \color{darkred}{\mathit{completion}}_j = \color{darkred}{\mathit{start}}_j +\color{darkblue}{\mathit{proctime}}_j && \forall j\\ &\color{darkred}{\mathit{tardy}}_j \ge\color{darkred}{\mathit{completion}}_j -\color{darkblue}{\mathit{due}}_j&& \forall j\\ &\color{darkred}{\mathit{completion}}_i \le\color{darkred}{\mathit{start}}_j && \forall i,j|\color{darkblue}{\mathit{precedence}}_{i,j} && \\ \hline & \color{darkred}x_{j,m,t} \in \{0,1\} && \\ & \color{darkred}{\mathit{obj}}_{\mathit{objs}} \ge 0 \\ & \color{darkred}{\mathit{tardy}}_j \ge 0 \\ &\color{darkred}{\mathit{start}}_j \ge \color{darkblue}{\mathit{release}}_j \end{align}

Notes:
• The sections of the model are: (1) the overall objective, (2) the individual objectives, (3) the constraints, and (4) the variable declarations.
• We assume here the data is integer-valued. If you have fractional values, you can use a combination of scaling and rounding to make the problem suited for this model.
• The variables $$\color{darkred}{\mathit{start}}_j$$ and $$\color{darkred}{\mathit{completion}}_j$$ can be declared free or positive.
• $$\color{darkred}{\mathit{obj}}_{\mathit{maxtardy}}$$ and $$\color{darkred}{\mathit{obj}}_{\mathit{makespan}}$$ may not be correct if the corresponding $$\color{darkblue}{\mathit{objw}}$$ is zero. Even if the weight is strictly positive, I have seen cases where the solver delivers solutions that are not tight [3]. For production code, it is recommended to recalculate objectives that are purely defined by a bound. In our case that is $$\color{darkred}{\mathit{obj}}_{\mathit{maxtardy}}$$ and $$\color{darkred}{\mathit{obj}}_{\mathit{makespan}}$$.
• As long as the planning horizon is long enough (number of time periods $$t$$), the model will always be feasible. That means: if the model is infeasible, increase the number of time periods. A good way to establish how many time periods you need is to find a feasible solution using some simple heuristic. Once you have a feasible solution, you can dimension the model accordingly.
• The set $$\color{darkblue}{\mathit{isTardy}}_{j,t}$$ indicates if  a job $$j$$ starting at $$t$$  is tardy (i.e. that $$t+\color{darkblue}{\mathit{proctime}}_j \gt \color{darkblue}{\mathit{due}}_j$$). This can be calculated in advance.
• There is a degree of symmetry in the model. We can rename the machines without affecting the solution. This means we can fix the first job to machine 1. A better approach is to find a measure for each machine (say the number of jobs allocated to that machine) and order the machines accordingly. This is what I did. I.e. machine 1 has the largest number of jobs, then machine 2, etc.
• This model was used to generate pictures in the previous section. It solves very fast, but that is very problem-dependent: slightly different data may lead to longer solution times.

#### Continuous-time model

A continuous time model wouldn't use a time index. Instead, we use a pairwise no-overlap constraint for jobs assigned to the same machine: $\color{darkred}x_{i,m}=\color{darkred}x_{j,m}=1\Rightarrow \color{darkred}{\mathit{completion}}_i \le \color{darkred}{\mathit{start}}_j\>{\mathbf{or}}\> \color{darkred}{\mathit{completion}}_j \le \color{darkred}{\mathit{start}}_i\>\>\forall i\lt j, \forall m$ We only need these for $$i\lt j$$. We can linearize the or as: \begin{align} &\color{darkred}{\mathit{completion}}_i \le \color{darkred}{\mathit{start}}_j + \color{darkblue} M \cdot (1-\color{darkred}x_{i,m}) +\color{darkblue} M \cdot (1-\color{darkred}x_{j,m}) + \color{darkblue} M \cdot \color{darkred} \delta_{i,j} && \forall i\lt j, \forall m \\ & \color{darkred}{\mathit{completion}}_j \le \color{darkred}{\mathit{start}}_i + \color{darkblue} M \cdot (1-\color{darkred}x_{i,m}) +\color{darkblue} M \cdot (1-\color{darkred}x_{j,m}) + \color{darkblue} M \cdot (1- \color{darkred} \delta_{i,j}) && \forall i\lt j, \forall m \\ & \color{darkred} \delta_{i,j} \in \{0,1\}\end{align}

The complete model looks like:

Model M2: Continuous time
\begin{align}\min&\sum_{\mathit{objs}} \color{darkblue}{\mathit{objw}}_{\mathit{objs}} \cdot \color{darkred}{\mathit{obj}}_{\mathit{objs}} && \\ \hline & \color{darkred}{\mathit{obj}}_{\mathit{completion}} = \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{completion}}_j \\ & \color{darkred}{\mathit{obj}}_{\mathit{sumtardy}} = \sum_j \color{darkblue}w_j \cdot \color{darkred}{\mathit{tardy}}_j \\ & \color{darkred}{\mathit{obj}}_{\mathit{numtardy}} = \sum_j \color{darkred}{\mathit{isTardy}}_j \\ & \color{darkred}{\mathit{obj}}_{\mathit{maxtardy}} \ge \color{darkred}{\mathit{tardy}}_j && \forall j\\ & \color{darkred}{\mathit{obj}}_{\mathit{makespan}} \ge \color{darkred}{\mathit{completion}}_j && \forall j\\ \hline & \sum_m \color{darkred}x_{j,m}=1 && \forall j \\ & \color{darkred}{\mathit{completion}}_ j = \color{darkred}{\mathit{start}}_ j + \color{darkblue}{\mathit{proctime}}_ j && \forall j \\ &\color{darkred}{\mathit{completion}}_i \le \color{darkred}{\mathit{start}}_j + \color{darkblue} M \cdot (1-\color{darkred}x_{i,m}) +\color{darkblue} M \cdot (1-\color{darkred}x_{j,m}) + \color{darkblue} M \cdot \color{darkred} \delta_{i,j} && \forall i\lt j, \forall m \\ & \color{darkred}{\mathit{completion}}_j \le \color{darkred}{\mathit{start}}_i + \color{darkblue} M \cdot (1-\color{darkred}x_{i,m}) +\color{darkblue} M \cdot (1-\color{darkred}x_{j,m}) + \color{darkblue} M \cdot (1- \color{darkred} \delta_{i,j}) && \forall i\lt j, \forall m\\ & \color{darkred}{\mathit{tardy}}_j \ge\color{darkred}{\mathit{completion}}_j -\color{darkblue}{\mathit{due}}_j && \forall j\\ &\color{darkred}{\mathit{tardy}}_j \le \color{darkblue}M \cdot\color{darkred}{\mathit{isTardy}}_j&& \forall j \\ &\color{darkred}{\mathit{completion}}_i \le\color{darkred}{\mathit{start}}_j && \forall i,j|\color{darkblue}{\mathit{precedence}}_{i,j} \\ \hline & \color{darkred}x_{j,m}\in \{0,1\} \\ & \color{darkred}\delta_{i,j} \in \{0,1\} \\ &\color{darkred}{\mathit{isTardy}}_j\in \{0,1\} \\ & \color{darkred}{\mathit{obj}}_{\mathit{objs}} \ge 0 \\ & \color{darkred}{\mathit{tardy}}_j \ge 0 \\ & \color{darkred}{\mathit{start}}_j \ge \max(1,\color{darkblue}{\mathit{release}}_j) \end{align}

Notes:
• The constant $$\color{darkblue} M$$ can be set to the length of the planning period. So although we have no time index and associated dimensioning problem, we still need to have some idea about the planning horizon.
• We can use an ordering on the machines to reduce symmetry.
• There are no assumptions about the integrality of the data.
• This model has serious problems in finding proven optimal solutions. It finds good solutions pretty quickly, but if we want a proven optimal solution, we need to wait a long time.

#### Results

I ran the two models with the following objective weights:

----     69 PARAMETER objw  objective weights

completion 0.001,    sumtardy   1.000


This means: first try to reduce the sum of the tardiness in the solution, and secondly, minimize the (weighted) sum of the completion times, so we move things as much to the left as possible.

model m1
time-indexed
model m2
continuous time
Equations62610,026
Variables20,1581,633
binary20,0001,475
Nonzero elements191,02449,896
Objective324.096324.102
completion2,0962,102
sumtardy322322
maxtardy8484
numtardy77
makespan9797
Time (seconds)201800 (time limit)
Gapoptimal0.13%
Nodes5291,262,868
Iterations28,5535,463,876

For this model (and data set) the time-indexed model does way better than the continuous-time model.

#### Conclusions

1. Different formulations can make a large difference. Often more than using solver options.
2. In this case, we have two very different formulations: one that uses time slots (and has time as an index of the main decision variable), and the other that allows continuous time. The difference in performance is striking.
3. The two solutions are almost identical. That gives some confidence that the models are correct. For complex models, it is not a bad idea to verify solutions this way.
4. To interpret solutions of scheduling models, it is paramount to use some GUI or plotting tool. Just printing the solution as a table is not that meaningful.

In a follow-up post, I will discuss some other, and stranger, formulations.

#### References

1. Yasin Unlu, Scott J. Mason, Evaluation of mixed integer programming formulations for non-preemptive parallel machine scheduling problems, Computers & Industrial Engineering 58 (2010) 785–800.
2. Visualization of a large scheduling solution, https://yetanothermathprogrammingconsultant.blogspot.com/2019/03/visualization-of-large-scheduling.html
3. Cplex bug or not? https://yetanothermathprogrammingconsultant.blogspot.com/2016/08/cplex-bug-or-not.html