As part of a Non-Permutation Flowshop Scheduling (NPFS) problem, I would like my MIP model to be able to deal with sequence dependent setup times. That is, for each pair of consecutive jobs, a setup time (which may be 0) may be induced in between.
Modeling:
Partially following the notation from Souiden et al., this means I extend the following inequality: $$S_{i,k} + \sum_{j=1}^{n}p_{i,j}\cdot x_{i,j,k} \leq \quad S_{i,k+1} \quad \forall i\in\mathcal{M};\forall k\in\mathcal{K}_i\setminus\{k_i\}.$$
In the above, the binary decision variable $x_{i,j,k}=1$ IFF job $j$ is assigned to machine $i$ in the $k$'th position on that machine. Variable $S_{i,k}$ denotes the start time of the job scheduled $k$'th on machine $i$ and $p_{i,j}$ the processing time of job $j$ on machine $i$. In other words, the above inequality states that two consecutive jobs on the same machine should be non-overlapping.
When using setup times, this means that we should add the setup time between the job scheduled $k$'th and $k+1$'th to the LHS in the above. In case the setup time between job $j$ and $j'$ on machine $i$ is denoted with $O_{i,j',j}$, my idea was to write setup time as follows:
$$\sum_{j'\in\mathcal{J}}\sum_{j\in\mathcal{J}}x_{i,j',k}\cdot x_{i,j,k+1}\cdot O_{i,j',j}.$$
The reason that I am using this sum of products is that we have to know which jobs are scheduled $k$'th and $k+1$'th on machine $i$. Now only for one pair of $(j',j)$ the product $x_{i,j',k}\cdot x_{i,j,k+1}=1$, hence this term.
Implementation
I implemented my model in Python and use Gurobi to solve it. For the line above, I used the code
quicksum(x[i, j1, k] * x[i, j2, k + 1] * O[i, j1,j2] for j1 in range(J)
I tested if model performance would increase if I explicitly linearized the setup term instead of using the * operator in Gurobi, however this made performance worse (Gurobi does it in a more clever way apparently).
Observation & Question
I observed that adding this setup term to my model, significantly decreases the model performance, even if the setup times are always the same value (regardless of what pair of consecutive jobs). If I would add that always-the-same value explicitly to the top inequality, model performance is much faster compared to when modeling it implicitly via the last term.
- What can be the reason for the observation above?
- Does anyone have some suggestions for a better way of formulating and/or implementing this, like cuts, ways to tighten the problem?
Why is the model slower? Hard to say. Part of it is the "cost" of linearizing the products, where the cost is partly additional constraints and maybe variables added by Gurobi, part is loosening of the best bound, and perhaps part is Gurobi being enticed to separate nodes on variables that are perhaps less productive than the ones used by the original model.
I'm not a Gurobi user, but I assume it has a mechanism for setting branching prioritize. It's possible (but far from guaranteed) that prioritizing the $x$ variables in ascending slot order might help. (Knowing job 3 on machine 1 is in slot 17 doesn't say much about start/end times if you don't know what occupies slots 1-16.)
I think an alternative model formulation is to define $x_{ijk}$ to be 1 if job $j$ precedes job $k$ (not necessarily immediately) on machine $i$. I'll assume that all jobs use all machines; if not, you can eliminate certain variables and certain constraints. For each $j\ne k$ and each machine $i$, $x_{ijk} + x_{ikj} = 1$, which effectively cuts the number of binary variables in half. $S_{ij}$ is now the start time for job $j$ on machine $i$, and you have the constraint $$S_{ij} + p_{ij} \le S_{ik} + M_i(1-x_{ijk})\quad \forall i\in\mathcal{M}, \ \forall j\ne k\in \mathcal{J},$$ where $\mathcal{J}$ is the set of jobs and $M_i$ is the dreaded "big M" (something like the sum of all processing times on machine $i$). To add in setup times, if $\sigma_{ijk}$ is the time to set up machine $i$ for job $k$ after it finishes job $j$, the constraint just becomes $$S_{ij} + p_{ij} + \sigma_{ijk} \le S_{ik} + M_i(1-x_{ijk})\quad \forall i\in\mathcal{M}, \ \forall j\ne k\in \mathcal{J}$$ with $M_i$ adjusted to account for both processing and (worst-case) setup times. Whether that is faster or slower than your model is an empirical question.