## Monday, April 19, 2021

### Parallel machine scheduling II, three more formulations

In [1] two formulations were discussed for a scheduling problem with multiple machines. Here we add a few more. Some of them are a bit strange. None of them really works better than the ones in [1].

So this is a collection of formulations that may sound reasonable at first, but are really not suited for this problem. If you want to read about good models, skip this post.

#### Summary of the input data

----     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


We implement a number of objectives, so we combine them and get something that is similar to what is used in practical situations.

In [1] we discussed two models:
• A time-indexed model that used the start of a job as central variable: $\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}$ This model did very well.
• A continuous-time model, where we model no-overlap using standard big-M constraints based on the Alan Manne formulation to sequence jobs such that they don't overlap [2]. This model did not do as well. It found good solutions fairly quickly but finding and proving the global optimum was challenging.
 An optimal solution

In this post I show some results with the following models:
• M3: A version of a time-indexed model that can be used when the run length of a job is not known in advance. That more general approach leads to a slower model.
• M4: A version of a continuous-time model from [3]. It is more complicated but not better than the model I used in [1].
• M5: A positional formulation. It works, but it is not competitive.

#### M3: A slow time-indexed model

Here we try to use as central variable: $\color{darkred}x_{j,m,t} = \begin{cases} 1 & \text{if job j occupies machine m at time t} \\ 0 &\text{otherwise}\end{cases}$ This approach makes the no-overlap constraint very straightforward: $\sum_j \color{darkred}x_{j,m,t} \le 1 \>\>\forall m,t$ To find the start and end time we can do: \begin{align} & \color{darkred}{\mathit{xstart}}_{j,m,t}\ge\color{darkred}x_{j,m,t} - \color{darkred}x_{j,m,t-1} \\ & \color{darkred}{\mathit{xend}}_{j,m,t}\ge\color{darkred}x_{j,m,t} - \color{darkred}x_{j,m,t+1} \\ & \sum_{m,t} \color{darkred}{\mathit{xstart}}_{j,m,t} \le 1&&\forall j \\ &\sum_{m,t}\color{darkred}x_{j,m,t} = \color{darkblue}{\mathit{proctime}}_j && \forall j \\ & \color{darkred}x_{j,m,t} \in \{0,1\} \\ & \color{darkred}{\mathit{xstart}}_{j,m,t}, \color{darkred}{\mathit{xend}}_{j,m,t} \in [0,1] \end{align} The first two inequalities implement the implications \begin{align} & \color{darkred}x_{j,m,t-1}= 0 \> \mathbf{and}\> \color{darkred}x_{j,m,t}=1 \Rightarrow \color{darkred}{\mathit{xstart}}_{j,m,t} = 1 \\ & \color{darkred}x_{j,m,t}= 1 \>\mathbf{and}\>\color{darkred}x_{j,m,t+1}= 0 \Rightarrow \color{darkred}{\mathit{xend}}_{j,m,t} = 1\end{align} The condition that we have maximum one start has two implications: first we have one start and thus one consecutive production run (without holes) and second that the production run is executed on one machine.

 Relation between x, xstart and xend.

It is noted that $$\color{darkred}x_{j,m,t}$$ is a real binary variable. The variables $$\color{darkred}{\mathit{xstart}}_{j,m,t}$$ and $$\color{darkred}{\mathit{xend}}_{j,m,t}$$ can be continuous between 0 and 1.

I also added the constraint $\sum_{m,t} \color{darkred}{\mathit{xend}}_{j,m,t} \le 1$ to enforce that otherwise unrestricted $$\color{darkred}{\mathit{xend}}$$ variables are set to zero. This makes it easier to calculate the $$\color{darkred}{\mathit{obj}}_{\mathit{numtardy}}$$ objective.

The complete model can look like:

Model M3: Slow time indexed model
\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|t\ge\color{darkblue}{\mathit{due}}(j)}\color{darkred}{\mathit{xend}}_{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& \color{darkred}{\mathit{xstart}}_{j,m,t}\ge\color{darkred}x_{j,m,t} - \color{darkred}x_{j,m,t-1} \\ & \color{darkred}{\mathit{xend}}_{j,m,t}\ge\color{darkred}x_{j,m,t} - \color{darkred}x_{j,m,t+1} \\ & \sum_{m,t} \color{darkred}{\mathit{xstart}}_{j,m,t} \le 1&&\forall j \\ & \sum_{m,t} \color{darkred}{\mathit{xend}}_{j,m,t} \le 1&&\forall j \\ &\sum_{m,t}\color{darkred}x_{j,m,t} = \color{darkblue}{\mathit{proctime}}_j && \forall j \\ &\sum_j \color{darkred}x_{j,m,t} \le 1 && \forall m,t\\ & \color{darkred}{\mathit{start}}_j = \sum_{m,t} t \cdot \color{darkred}{\mathit{xstart}}_{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{xstart}}_{j,m,t}, \color{darkred}{\mathit{xend}}_{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}

This model works, but it is slow. This formulation is more often used in situations where the run lengths are not fixed. An example is:  a power generator can only be turned on at most $$n$$ times during a time period to prevent excessive wear and tear. Here the situation is slightly different: we know there is exactly one production run for a given job and its length is fixed.

model m1
time-indexed
model m2
continuous-time
model m3
time-indexed
Equations62610,02640,726
Variables20,1581,63360,158
binary20,0001,47520,000
Nonzero elements191,02449,896250,752
Objective324.096324.102437.303
completion2,0962,1022,303
sumtardy322322435
maxtardy848487
numtardy7714
makespan9797100
Time (seconds)201800 (time limit)1800 (time limit)
Gapoptimal0.13%30%
Nodes5291,262,8685,317
Iterations28,5535,463,87613,311,418

The conclusion is that relatively minor variations in a model can lead to very different performance. Also, a formulation that exploits as much as possible of the problem at hand (in this case: fixed run lengths), is likely the best.

#### M4: an alternative, complex continuous-time model

In [1] a very simple continuous-time model was implemented, based on the implication:

if jobs $$i$$ and $$j$$ execute on the same machine then job $$j$$ ends before job $$i$$ starts or job $$j$$ starts after job $$i$$ ends.
This was implemented 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} Note that all comparisons are only for $$i\lt j$$. Here $$\color{darkred}x$$ indicates the assignment of jobs to machines and $$\color{darkred} \delta$$ implements the or condition. The constant $$\color{darkblue}M$$ should be equal to the length of the planning period.

In [3], a more convoluted version of this approach is proposed. Let's see how it behaves.

We introduce the following variables: \begin{aligned}& \color{darkred}x_{j,m}=\begin{cases} 1 & \text{if job j is assigned to machine m} \\ 0 & \text{otherwise}\end{cases}\\ & \color{darkred}\delta_{i,j} = \begin{cases} 1 & \text{if job i executes before job j on the same machine}\\ 0 & \text{otherwise}\end{cases} \\ &\color{darkred}d_{i,j} = \begin{cases} 1 & \text{if jobs i and j execute on different machines}\\ 0 & \text{otherwise}\end{cases}\end{aligned}

The first constraint is $\color{darkred}\delta_{i,j} + \color{darkred}\delta_{j,i}+\color{darkred}d_{i,j} = 1\>\>\forall i \lt j$ This means: jobs $$i$$ and $$j$$ are on different machines or they don't overlap. The second constraint, the transitivity constraint, is a bit of mystery to me: $\color{darkred}\delta_{i,j} +\color{darkred}\delta_{j,k} + \color{darkred}\delta_{k,i}\le 2 \>\>\forall i \lt j \lt k$ I don't think this constraint is needed. It may have been added to strengthen the formulation. The next constraint $\color{darkred}x_{i,m}+\color{darkred}x_{j,m}+\color{darkred}d_{i,j}\le 2 \>\>\forall i\lt j$ can be considered as an implication $\color{darkred}x_{i,m}=\color{darkred}x_{j,m}=1 \Rightarrow \color{darkred}d_{i,j}=0$ I.e. if jobs $$i$$ and $$j$$ are on same the machine, we need to have $$\color{darkred}d_{i,j}=0$$. This is used in the no-overlap (first) constraint. Obviously we have the assignment constraint: $\sum_m \color{darkred}x_{i,m}=1\>\>\forall i$ Finally, we have the completion times to consider: $\color{darkred}{\mathit{completion}}_j \ge \color{darkred}{\mathit{completion}}_i + \color{darkblue}{\mathit{proctime}}_j(\color{darkred}\delta_{i,j}+\color{darkred}x_{i,m}+\color{darkred}x_{j,m}-2) - \color{darkblue}M(1-\color{darkred}\delta_{i,j})\>\>\forall i,j,m$ Instead, we can use the simpler: $\color{darkred}{\mathit{start}}_j \ge \color{darkred}{\mathit{completion}}_i - \color{darkblue}M(1-\color{darkred}\delta_{i,j})\>\>\forall i\ne j, m$ Note that I explicitly excluded the case $$i=j$$ which is a bit more precise.

The complete model can be summarized as follows:

Model M4: Alternative continuous-time model
\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 &\color{darkred}\delta_{i,j} + \color{darkred}\delta_{j,i}+\color{darkred}d_{i,j} = 1&& \forall i \lt j \\ &\color{darkred}\delta_{i,j} +\color{darkred}\delta_{j,k} + \color{darkred}\delta_{k,i}\le 2 &&\forall i \lt j \lt k\\ & \color{darkred}x_{i,m}+\color{darkred}x_{j,m}+\color{darkred}d_{i,j}\le 2 && \forall i\lt j \\ & \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{start}}_j \ge \color{darkred}{\mathit{completion}}_i - \color{darkblue}M(1-\color{darkred}\delta_{i,j})&&\forall i\ne j, 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}d_{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}

The performance is as follows:

model m1
time-indexed
model m2
continuous-time
model m4
continuous-time
Equations62610,02628,451
Variables20,1581,6334,133
binary20,0001,4752,750
Nonzero elements191,02449,89685,571
Objective324.096324.102324.139
completion2,0962,1022,139
sumtardy322322322
maxtardy848484
numtardy777
makespan979797
Time (seconds)201800 (time limit)1800 (time limit)
Gapoptimal0.13%0.15%
Nodes5291,262,868114,732
Iterations28,5535,463,8763,455,261

Notes:
• The performance is very close to model M2. This formulation does not add much value above the more intuitive and simpler model M2. The M2 model is easier to remember. More complex models are not necessarily better.
• I don't know what the purpose is of the transitivity constraint.
• Some of the constraints should have $$\forall i\ne j$$ instead of $$\forall i,j$$.
• Of course, we can add ordering constraints to reduce symmetry. I used: machine $$m$$ has more jobs (or the same number) than machine $$m+1$$. As the machines in my model are identical, we can renumber machines.

#### M5: positional formulation.

This is also from [3]. Here we use as central variable: $\color{darkred}x_{j,m,p}=\begin{cases} 1 & \text{if job j has position p on machine m}\\ 0 & \text{otherwise}\end{cases}$ The assignment constraints related to this variable are: \begin{align}&\sum_{m,p}\color{darkred}x_{j,m,p}=1 &&\forall j \\ &\sum_j\color{darkred}x_{j,m,p}\le 1 &&\forall m,p \end{align} The corresponding completion time is denoted by $$\color{darkred}f_{m,p}$$.  We have: $\color{darkred}f_{m,p} \ge \color{darkred}f_{m,p-1} + \sum_j \color{darkblue}{\mathit{proctime}}_j\cdot \color{darkred}x_{j,m,p}$ We assume here the first position is 1 and $$\color{darkred}f_{m,0}=1$$ (this is because my first time period is 1 and completion time is at the beginning of a time period). Finally we link the job completion time through: $\color{darkred}{\mathit{completion}}_j \ge \color{darkred}f_{m,p} - \color{darkblue}M(1-\color{darkred}x_{j,m,p}) \>\>\forall j,m,p$

To allow for release dates and precedence constraints, we need to add: \begin{align}& \color{darkred}f_{m,p} \ge \sum_j (\color{darkblue}{\mathit{release}}_j+\color{darkblue}{\mathit{proctime}}_j)\cdot \color{darkred}x_{j,m,p}\\ & \color{darkred}f_{m,p} \ge \color{darkred}{\mathit{completion}}_i + \color{darkblue}{\mathit{proctime}}_j\cdot \color{darkred}x_{j,m,p} - \color{darkblue}M\cdot (1-\color{darkred}x_{j,m,p})&& \forall i,j|\color{darkblue}{\mathit{precedence}}_{i,j} \end{align} The last constraint is not in [3]: they only consider release dates but no precedence constraints.

Model M5: Positional model
\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,p}\color{darkred}x_{j,m,p}=1 &&\forall j \\ & \sum_j\color{darkred}x_{j,m,p}\le 1 &&\forall m,p \\ & \color{darkred}f_{m,p} \ge \color{darkred}f_{m,p-1} + \sum_j \color{darkblue}{\mathit{proctime}}_j\cdot \color{darkred}x_{j,m,p} &&\forall m,p\\ & \color{darkred}f_{m,0}=1 \\ & \color{darkred}f_{m,p} \ge \sum_j (\color{darkblue}{\mathit{release}}_j+\color{darkblue}{\mathit{proctime}}_j)\cdot \color{darkred}x_{j,m,p} && \forall m,p \\ & \color{darkred}f_{m,p} \ge \color{darkred}{\mathit{completion}}_i + \color{darkblue}{\mathit{proctime}}_j\cdot \color{darkred}x_{j,m,p} - \color{darkblue}M\cdot (1-\color{darkred}x_{j,m,p})&& \forall i,j|\color{darkblue}{\mathit{precedence}}_{i,j},m,p \\ & \color{darkred}{\mathit{completion}}_j \ge \color{darkred}f_{m,p} - \color{darkblue}M(1-\color{darkred}x_{j,m,p}) && \forall j,m,p \\ & \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{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,p}\in \{0,1\} \\ & \color{darkred}f_{m,p} \ge 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}

We can add a few refinements. First, we can renumber the machines as they are identical. I reduce symmetry by ordering the machines by how many jobs they execute. The second source of symmetry is that we can have "empty" positions in this model. I added a constraint to make sure the empty, unused positions are at the end and not in the middle. This also adds a symmetry-breaking effect to the model. After this we see:

model m1
time-indexed
model m2
continuous-time
model m5
positional formulation
Equations62610,0265,742
Variables20,1581,6334,368
binary20,0001,4754,050
Nonzero elements191,02449,89636,564
Objective324.096324.102324.135
completion2,0962,1022,135
sumtardy322322322
maxtardy848484
numtardy777
makespan979797
Time (seconds)201800 (time limit)1800 (time limit)
Gapoptimal0.13%0.14%
Nodes5291,262,86849,186
Iterations28,5535,463,8762,368,537

Notes:
• This model is not as effective as the models in [1].
• There are some interesting possibilities to reduce symmetry in this model.

#### Conclusion

There are many formulations for this problem, choosing the right one can pay off.

Having access to different formulations can be an excellent debugging tool: they all should deliver the same optimal objective. We can plug in a solution from one model into another (using mipstart or by fixing): it should be accepted without hesitancy. Complex models are not so easy to get right immediately, so an extra tool to verify solutions is valuable.

#### References

1. Parallel machine scheduling I, two formulations, https://yetanothermathprogrammingconsultant.blogspot.com/2021/03/parallel-machine-scheduling-i-two.html
2. Manne, Alan S. “On the Job-Shop Scheduling Problem.” Operations Research, vol. 8, no. 2, 1960, pp. 219–223.
3. 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.

#### 1 comment:

1. Interesting analysis.
I'm pleased to see that the simplest model performed best - often modellers have a tendency to complicate models unnecessarily, just because they can.
Your point about model validation is very important. Too many times I've seen models that get an optimal solution to the wrong problem, due to errors in either the situation specification or the model implementation. Having multple models can give us confidence that we have a sensible result.
FYI, I've written a brief summary at https://www.solvermax.com/blog/parallel-machine-scheduling