Data and solution
It is always helpful to start with an actual data set. So here we go:
---- 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.
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 | |
|---|---|---|
| Equations | 626 | 10,026 |
| Variables | 20,158 | 1,633 |
| binary | 20,000 | 1,475 |
| Nonzero elements | 191,024 | 49,896 |
| Objective | 324.096 | 324.102 |
| completion | 2,096 | 2,102 |
| sumtardy | 322 | 322 |
| maxtardy | 84 | 84 |
| numtardy | 7 | 7 |
| makespan | 97 | 97 |
| Time (seconds) | 20 | 1800 (time limit) |
| Gap | optimal | 0.13% |
| Nodes | 529 | 1,262,868 |
| Iterations | 28,553 | 5,463,876 |
For this model (and data set) the time-indexed model does way better than the continuous-time model.
Conclusions
- Different formulations can make a large difference. Often more than using solver options.
- 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.
- 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.
- 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.
References
- 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.
- Visualization of a large scheduling solution, https://yetanothermathprogrammingconsultant.blogspot.com/2019/03/visualization-of-large-scheduling.html
- Cplex bug or not? https://yetanothermathprogrammingconsultant.blogspot.com/2016/08/cplex-bug-or-not.html
No comments:
Post a Comment