# Hybridization of mixed-integer linear program and discrete event systems for robust scheduling on parallel machines

Alexis Aubry, Pascale Marangé, David Lemoine, Sara Himmiche, Sylvie Norre

► To cite this version:

Alexis Aubry, Pascale Marangé, David Lemoine, Sara Himmiche, Sylvie Norre. Hybridization of mixed-integer linear program and discrete event systems for robust scheduling on parallel machines. International Conference on Advances in Production Management Systems, APMS 2021, Sep 2021, Nantes, France. pp.73-80, 10.1007/978-3-030-85874-2_8 . hal-03337509

HAL Id: hal-03337509

https://hal.science/hal-03337509

Submitted on 15 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

ifip

International Federation for Information Processing

This document is the original author manuscript of a paper submitted to an IFIP
conference proceedings or other IFIP publication by Springer Nature. As such, there
may be some differences in the official published version of the paper. Such
differences, if any, are usually due to reformatting during preparation for publication or
minor corrections made by the author(s) during final proofreading of the publication
manuscript.
# Hybridization of mixed-integer linear program and discrete event systems for robust scheduling on parallel machines

A. Aubry¹, P. Marangé¹, D. Lemoine², S. Himmiche¹, and S. Norre³

¹ Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France
{alexis.aubry, pascale.marange, sara.himmiche}@univ-lorraine.fr

² IMT-Atlantique, LS2N UMR 6004, Nantes
*david.lemoine@imt-atlantique.fr*

³ Clermont Auvergne University, LIMOS UMR 6158
*sylvie.norre@uca.fr*

**Abstract.** This paper proposes an approach for robust scheduling on parallel machines. This approach is based on a combination of robust mathematical and discrete event systems models which are iteratively called in order to converge towards a schedule with the required robustness level defined by the decision maker. Experimentations on a small instance (10 jobs and 2 unrelated machines) and a more complex one (30 jobs and 6 uniform machines) show that this approach permits to converge quickly to a robust schedule even if the probability distribution associated to the uncertainties are not symmetrical. The approach achieves a better rate of convergence than those of the literature's methods.

**Keywords:** Robust scheduling · robust mixed integer programming model · discrete event systems · parallel machines

## Introduction

Scheduling under uncertainties is still a present concern in Operation Research and Decision Aiding. Many researchers are interested in determining a robust schedule which is rather insensitive to the data uncertainties and which is able to absorb the perturbations without unreasonably degrading its performances. However, the concept of robustness is differently defined according to the domains. An approach for robust optimization has been proposed by [1], based on a Robust Mixed Integer Programming Model. In order to obtain the robustness level wanted by the decision maker, a vector $\Omega$, which corresponds to the maximum allowed deviation of input parameters, has to be fixed. This vector can be interpreted as a robustness coefficient. The authors have shown that if the uncertainties on parameters are independent and symmetrically distributed, the vector $\Omega$ can be analytically computed. But, such an hypothesis is not often verified in real scheduling problems. In [2], a methodology for iteratively and numerically tuning $\Omega$ has been proposed thanks to a combination of a robust
mathematical programming and Discrete Event Systems models when the probability distribution associated to the uncertainties are not symmetrical. This generic method has been applied on a scheduling problem with parallel machines and the results show that this approach is a good mean for tuning the $\Omega$ parameters. However, the mechanism for updating coefficients $\Omega$ from one iteration to another was simple and naive. Moreover, it did not take into account the characteristics of the current schedule solution (in particular the load of the machines). As a result, the convergence of the method toward a solution with the required robustness level was relatively slow. Our contribution in this paper is to propose a new mechanism for updating robustness coefficients $\Omega$ to increase the rate of convergence of our method in case of scheduling problem on parallel machines. The problem is thus denoted as $RK||C_{max}$ which consists in minimizing the makespan ($C_{max}$) on K parallel machines.

The paper is built as follows. The first section presents the robust mixed integer programming model for scheduling on parallel machines. The second section presents the methodology and details the mechanism for updating the robustness coefficients $\Omega$. The third section discusses the results on two instances: a simple instance composed of 10 jobs and 2 unrelated machines and a more complex one composed of 30 jobs and 6 uniform machines. Finally, the last section deals with the conclusion and the perspectives.

# 1 Robust Formulation of $RK||C_{max}$

The main assumptions for $RK||C_{max}$ are the following:

- all jobs are available at time 0,

- the K machines are always available (no breakdown etc.),

- processing times for the jobs are independent,

- a machine cannot process more than one job at any time and preemption is not allowed.

Parameters and decision variables of the model are summarized in Table 1.

<table><caption>Table 1: Notations of the model</caption><tr><td>N</td><td>: Number of jobs which have to be scheduled</td></tr><tr><td>K</td><td>: Number of parallel machines</td></tr><tr><td>t<sub>jk</sub></td><td>: Processing time for job j on the machine k, $\forall (j, k) \in \{1, ..., N\} \times \{1, ..., K\}$</td></tr><tr><td>x<sub>jk</sub></td><td>: $\begin{cases} 1 & \text{if j is executed on machine k} \\ 0 & \text{otherwise} \end{cases}$, $\forall (j, k) \in \{1, ..., N\} \times \{1, ..., K\}$</td></tr><tr><td>C<sub>max</sub></td><td>: is the makespan value</td></tr></table>

Since processing times are uncertain, we assume that

$$ \forall (j, k) \in \{1, ..., N\} \times \{1, ..., K\}, \quad t_{jk} \in [t_{jk}^{\min}, t_{jk}^{\max}] $$
Following the methodology of [1], uncertain processing times are modelled as

$$t_{jk} = \bar{t}_{jk} + \zeta_{jk} \hat{t}_{jk}$$

where

$$- \bar{t}_{jk} = \frac{t_{jk}^{\max} + t_{jk}^{\min}}{2} \quad \text{and} \quad \hat{t}_{jk} = \frac{t_{jk}^{\max} - t_{jk}^{\min}}{2}$$

- $\zeta_{jk}$ is a random variable which takes its values in $[-1, 1]$

Thus, the robust model can be formulated as follows:

$$\text{Minimize } C_{\max} \qquad (1)$$

s.t.

$$\sum_{k=1}^{K} x_{jk} = 1 \quad \forall j \in \{1, \dots, N\} \qquad (2)$$

$$\sum_{j=1}^{N} \bar{t}_{jk} x_{jk} + \max_{\sum_{j=1}^{N} |\zeta_{jk}| \le \Omega_k} \left( \sum_{j=1}^{N} \hat{t}_{jk} \zeta_{jk} x_{jk} \right) \le C_{\max} \quad \forall k \in \{1, \dots, K\} \qquad (3)$$

$$x_{jk} \in \{0, 1\} \quad \forall (j, k) \in \{1, \dots, N\} \times \{1, \dots, K\} \qquad (4)$$

where $\Omega = (\Omega_k)_{k \in \{1, \dots, K\}}$ constrains the maximum deviation of processing times allowed on each machine $k$. It can be interpreted as a robustness coefficient: the larger $\Omega_k$ is, the more conservative is the constraint (3). We denote as $X^\Omega$ an optimal schedule provided by the robust model with the input $\Omega$ and as $C_{\max}(X^\Omega)$, the associated optimal makespan.

Let $\tilde{C}_{\max}(X^\Omega)$ be defined as the random variable associated to the makespan, when $X^\Omega$ is executed in the workshop (with the uncertain values $t_{jk}$). It is possible to define a robustness indicator as the probability that $\tilde{C}_{\max}(X^\Omega)$ is lower than $C_{\max}(X^\Omega)$. More formally, this indicator is given by the equation (5):

$$\Gamma(C_{\max}(X^\Omega)) = \mathbb{P}[\tilde{C}_{\max}(X^\Omega) \le C_{\max}(X^\Omega)] \qquad (5)$$

Our goal is to find $\Omega$ to guarantee that the scheduling $X^\Omega$ reaches a certain level of robustness $\Gamma^{ref}$. As [1] have shown that in case of each $\zeta_{jk}$ is symmetrically distributed in $[-1, 1]$, determining such $\Omega$ can be done analytically. Next section provides a methodology which allows to reach this goal in a general case.

## 2 Our methodology

In this section, we denote as $X_k^\Omega$ the schedule extracted from $X^\Omega$ by considering only the set of jobs processed on machine $k$, $k \in \{1, \dots, K\}$.

Basically, our methodology is based on the following observations:
1. In a workshop of parallel machines, machines are independent meaning that

$$ \Gamma(C_{\max}(X^{\Omega})) = \prod_{k=1}^{K} \Gamma(C_{\max}(X_k^{\Omega})) $$

2. if $\Gamma(C_{\max}(X^{\Omega})) < \Gamma^{ref}$ then there exists at least one $k \in \{1, \dots, K\}$ such that $\Gamma(C_{\max}(X_k^{\Omega})) < \sqrt[K]{\Gamma^{ref}}$ which is equivalent to

$$ \mathbb{P}[\tilde{C}_{\max}(X_k^{\Omega}) \le C_{\max}(X_k^{\Omega})] < \sqrt[K]{\Gamma^{ref}} $$

Thus, for such a $k$, we estimate a value $\tilde{d}(X_k^{\Omega}) > C_{\max}(X_k^{\Omega})$ for which the following constraint is satisfied.

$$ \mathbb{P}[\tilde{C}_{\max}(X_k^{\Omega}) \le \tilde{d}(X_k^{\Omega})] = \sqrt[K]{\Gamma^{ref}} \quad (6) $$

Therefore, to obtain this robustness score for the scheduling $X_k^{\Omega}$, we have to compute a new robustness coefficient $\Omega_k^{new}$ in order to satisfy the constraint (7).

$$ \max_{\sum_{j \in X_k^{\Omega}} |\zeta_{jk}| \le \Omega_k^{new}} \left( \sum_{j \in X_k^{\Omega}} \hat{t}_{jk} \zeta_{jk} \right) = \tilde{d}(X_k^{\Omega}) - \sum_{j \in X_k^{\Omega}} \bar{t}_{jk} \quad (7) $$

This constraint is derived from (3) in which we allow the deviation of uncertainties constrained by $\Omega_k^{new}$ until we reach $\tilde{d}(X_k^{\Omega})$.

As we want to be the less conservative as possible (that means, in our case, that $\Omega_k^{new}$ should be as small as possible), the following proposition can be stated:

**Proposition 1.** Solving equation (7) is equivalent to solve the following continuous Knapsack problem:

$$ \text{Minimize } \Omega_k^{new} = \sum_{j \in X_k^{\Omega}} \zeta_{jk} \quad (8) $$

s.t.

$$ \sum_{j \in X_k^{\Omega}} \hat{t}_{jk} \zeta_{jk} = \tilde{d}(X_k^{\Omega}) - \sum_{j \in X_k^{\Omega}} \bar{t}_{jk} \quad (9) $$

$$ \zeta_{jk} \in [0, 1] \quad \forall j \in X_k^{\Omega} \quad (10) $$

As the continuous Knapsack problem is known to be polynomial and can be easily solved by a greedy algorithm, we can use such algorithm to find the value of $\Omega_k^{new}$ that satisfies (7).

A similar reasoning can be applied if $\Gamma(C_{\max}(X^{\Omega})) > \Gamma^{ref}$ and leads to the same result.

Based on the previous proposition, our methodology can be summarized into 4 steps:
**Step 1:** For $\Omega = (\Omega_k)_{k \in \{1, \dots, K\}}$, an optimization module based on a linear solver provides $X^\Omega$ an optimal schedule according to the robust model presented in the previous section [1].

**Step 2:** Thanks to Discrete Event Systems models and tools, $\Gamma(C_{\max}(X^\Omega))$ is then evaluated. If $\Gamma(C_{\max}(X^\Omega)) \in [\Gamma^{\text{ref}} - \epsilon, \Gamma^{\text{ref}} + \epsilon]$ then the process stops and $X^\Omega$ is the required schedule solution, else the third step is engaged. $\epsilon$ is a parameter that helps to fix the required level of accuracy.

**Step 3:** For each $k \in \{1, \dots, K\}$, the value of $\tilde{d}(X_k^\Omega)$ that satisfies constraint (6) is then determined thanks to another Discrete Event Systems module. $\Omega$ is then updated in the Step 4.

**Step 4:** Thanks to Equation (7), $\forall k \in \{1, \dots, K\}$, $\Omega_k^{\text{new}}$ is then computed using the proposition 1 and we loop back to Step 1 with $\Omega = (\Omega_k^{\text{new}})_{k \in \{1, \dots, K\}}$.

Steps 2 and 3 use models and tools from discrete event systems (DES) [3]. Discrete event systems (DES) allow to model the behavior of a system by considering its possible states and the possible events (allowing the evolution from one state to another). The event is an instantaneous occurrence of an action or a phenomenon in the system environment. The evolution on the event occurrence can be deterministic when the behavior is known with certainty or stochastic when the behavior is uncertain and the evolution can lead to different states. The works of [4], [5], [6] have shown that DES are particularly relevant for scheduling evaluation due to their ability to model the behavior of industrial systems and perturbations. Indeed, DES allow to represent many dynamic features such as the communication between the elements of the workshop (jobs, resources), the time and the probabilistic behavior of perturbations. Many stochastic discrete event system languages allow to model these characteristics: Stochastic Petri Nets [7], Stochastic automata [8], Stochastic Automata Networks [9]. The method for evaluating the impact of a set of perturbations on a given schedule uses the approach proposed by [4]. The approach models the characteristics of the workshop (operations, resources), and the probability distribution associated to the perturbation by a discretization of this one and allows the evaluation of different elements: the robustness indicator $\Gamma(C_{\max}(X^\Omega))$ as in step 2, the minimal duration $\tilde{d}(X_k^\Omega)$ as in step 3, ...

# 3 Case Studies

## 3.1 A simple instance

In this first application, 10 jobs are considered. These jobs can be executed on two unrelated parallel machines (the problem is $R2||C_{\max}$). $\Gamma^{\text{ref}}$ is fixed to 90% and $\epsilon$ is fixed to 1% such that $\Gamma(C_{\max}(X^\Omega))$ has to be in [89%, 91%] for concluding to the acceptability of the solution $X^\Omega$. We applied the two approaches defined (from [2] and the presented one) by starting with $\Omega = [0, 0]$ (meaning that no uncertainty is considered).
<table><caption>Table 2: Obtained iterations for $R2||C_{max}$ with the approach from [2]</caption><thead><tr><th>Iteration i</th><th>Ω</th><th colspan="2">Solution</th><th>Ω<sup>new</sup></th></tr><tr><th></th><th></th><th>X<sup>Ω</sup></th><th>C<sub>max</sub>(X<sup>Ω</sup>)</th><th></th></tr></thead><tbody><tr><td>0</td><td>[0, 0]</td><td>X<sub>1</sub></td><td>12</td><td>[0.10, 0.18]</td></tr><tr><td>⋮</td><td>⋮</td><td>⋮</td><td>⋮</td><td>⋮</td></tr><tr><td>10</td><td>[0.53, 0.66]</td><td>X<sub>1</sub></td><td>13.98</td><td>∅</td></tr></tbody></table>

<table><caption>Table 3: Obtained iterations for the problem $R2||C_{max}$ with our approach</caption><thead><tr><th>Iteration i</th><th>Input Ω</th><th colspan="2">Step 1</th><th>Step 2</th><th>Step 3</th><th>Step 4</th></tr><tr><th></th><th></th><th>X<sup>Ω</sup></th><th>C<sub>max</sub>(X<sup>Ω</sup>)</th><th>Γ(C<sub>max</sub>(X<sup>Ω</sup>))</th><th>[d̃(X<sub>k</sub><sup>Ω</sup>)]<sub>k∈{1,2}</sub></th><th>Ω<sup>new</sup></th></tr></thead><tbody><tr><td>0</td><td>[0, 0]</td><td>X<sub>1</sub></td><td>12</td><td>0.65</td><td>[13, 14]</td><td>[0.50, 0.66]</td></tr><tr><td>1</td><td>[0.50, 0.66]</td><td>X<sub>1</sub></td><td>14</td><td>0.90</td><td>∅</td><td>∅</td></tr></tbody></table>

The table 2 summarizes the results by giving the extremal iterations of the combined approach presented in [2]. Two solutions are explored during the different iterations. The solution $X1$ allocates the first machine to jobs 1, 4, 6, 7, 8 and the second machine to jobs 2, 3, 5, 9, 10. The solution $X2$ allocates the first machine to jobs 1, 4, 6, 7, 8, 10 and the second machine to jobs 2, 3, 5, 9.

We can conclude that the improved approach finds the right solution in the first iteration (as $Γ(C_{max}(X^Ω)) = 90\%$, it is not necessary to update $Ω$ and to run a new iteration). Thus, the improved method decreases drastically the number of necessary iterations for converging to a comparable solution.

## 3.2 A more complex instance

We consider 30 jobs which can be executed on 6 uniform machines (the problems is $Q6||C_{max}$). The machines are such that the first machine is the fastest and the sixth is the slowest. $Γ^{ref}$ is fixed to 90% and $ε$ is fixed to 1%.

We applied the two approaches by starting with $Ω = [0, 0, 0, 0, 0, 0]$.

The figure 1 presents the evolution of the average value of $Ω$ and $Γ(C_{max}(X_k^Ω))$ according to the iteration when applying the two approaches. First, it can be observed that even after 10 iterations, the targeted performance $Γ^{ref}$ is not reached when applying the approach of [2] and $Γ(C_{max}(X_k^Ω))$ remains low. Moreover, $Ω$ increases very slowly when applying the approach of [2] in comparison with our approach. We can postulate that a lot of iterations is still necessary for reaching the target when applying the approach of [2].

The table 4 detailed the iterations when applying our improved approach. After the first iteration, the obtained value of $Γ(C_{max}(X^Ω))$ is lightly too high such that a second iteration for adjusting the values of $Ω$ is necessary. The second
<figure><img src="image_3.png" /><figcaption>Fig. 1: Comparison of the results when applying the two approaches</figcaption></figure>

<table><caption>Table 4: Obtained iterations for the problem $Q6||C_{max}$ with our approach</caption><thead><tr><th>Iteration i</th><th>Input Ω</th><th colspan="2">Step 1</th><th colspan="2">Step 2</th><th colspan="2">Step 3</th><th colspan="2">Step 4</th></tr><tr><td></td><td></td><td colspan="2">X<sup>Ω</sup> C<sub>max</sub>(X<sup>Ω</sup>)</td><td colspan="2">Γ(C<sub>max</sub>(X<sup>Ω</sup>))</td><td colspan="2">[d̃(X<sub>k</sub><sup>Ω</sup>)]<sub>k∈{1,2}</sub></td><td colspan="2">Ω<sup>new</sup></td></tr></thead><tbody><tr><td>0</td><td>[0, 0, 0, 0, 0]</td><td>X<sub>0</sub></td><td>926</td><td colspan="2">3.6%</td><td colspan="2">[968, 969, 974, 971, 976, 979]</td><td colspan="2">[2.67, 2.80, 2.50, 2.50, 2.23, 2.18]</td></tr><tr><td>1</td><td>[2.67, 2.80, 2.50, 2.50, 2.23, 2.18]</td><td>X<sub>1</sub></td><td>973</td><td colspan="2">91.1%</td><td colspan="2">[972, 970, 971, 973, 973, 976]</td><td colspan="2">[2.67, 2.67, 2.39, 2.50, 2.26, 2.32]</td></tr><tr><td>2</td><td>[2.67, 2.67, 2.39, 2.50, 2.26, 2.32]</td><td>X<sub>2</sub></td><td>972.72</td><td colspan="2">90.3%</td><td colspan="2">Ø</td><td colspan="2">Ø</td></tr></tbody></table>

iteration allows to reach the targeted performance. If we accept to degrade the optimal deterministic makespan ($C_{max} = 926$ with the solution $X_0$) of only 5%, then it is possible to guarantee this makespan despite the uncertainties with a probability of 90% giving a good compromise between optimality ($C_{max}$) and robustness ($\Gamma$). This confirms globally the results obtained with the simple instance: the improved approach permits to converge faster to a robust solution.

## 4 Conclusion and perspectives

We have proposed an approach combining robust mathematical programming and Discrete Event Systems models for the building of a robust scheduling on parallel machines. This allows to reach the level of robustness desired by the decision-maker by finely assessing the degree of robustness of the solutions provided by the optimization module, regardless of the probability distributions that follow the uncertainties on the model input data. The probability distribution associated to the uncertainties are supposed independent but not necessarily
symmetrical. Experiments on two instances show that our approach permits to converge quickly to a robust schedule and improves the rate of convergence of literature's methods. Several perspectives to this work can be considered. First, at short term, a more specific distribution of levels of robustness can be considered, taking into account, for example, the configuration of the production system, the criticality of certain machines (for instance, requiring greater robustness for bottleneck machines, ...). It would also be interesting to consider dependent probability distributions associated to the uncertainties. Long term perspectives will concern the extension of this approach to more complex shop scheduling problems as flow shop, job shop or hybrid flow shop.

## References

1. Bertsimas, D., Sim, M.: The Price of Robustness. Operations Research **2**(1), 35–53 (2004)

2. Marangé, P., Lemoine, D., Aubry, A., Himmiche, S., Norre, S., Bloch, C. Pétin, J.-F.: Coupling Robust Optimization and Model-Checking Techniques for Robust Scheduling in the Context of Industry 4.0. In: Sokolov B., Ivanov D., Dolgui A. (eds) Scheduling in Industry 4.0 and Cloud Manufacturing. International Series in Operations Research & Management Science, vol 289. Springer, Cham (2020)

3. Cassandras, C. G., Lafortune, S.: Introduction to discrete event systems. 3rd edn. Springer International Publishing, Cham (2021)

4. Himmiche, S., Aubry, A., Marangé, P., Pétin, J.-F., Duflot, M.: Using statistical-model-checking-based simulation for evaluating the robustness of a production schedule. In: 7th Workshop on Service Orientation in Holonic and Multi-Agent Manufacturing, SOHOMA’17, Oct 2017, Nantes, France (2017)

5. Lefebvre, D. and Mejia, G.: Robust scheduling in uncertain environment with petri nets and beam search. In: IFAC-PapersOnLine, **51**(11), 1077–1082 (2018)

6. Cherif, G., Leclercq, E., Lefebvre, D.: Scheduling Problems for a Class of Hybrid FMS Using T-TPN and Beam Search. Journal of Control, Automation and Electrical Systems, (2021)

7. Chiola, G., Marsan, M. A., Balbo, G., Conte, G.: Generalized stochastic petri nets: a definition at the net level and its implications. IEEE Transactions on Software Engineering, **19**(2), 89–107 (1993)

8. Alur, R., Dill, D. L.: A theory of timed automata. Theoretical Computer Science, **126**(2), 183–235 (1994)

9. Plateau, B., Atif, K.: Stochastic automata network of modeling parallel systems. IEEE Transactions on Software Engineering, **17**(10), 1093–1108 (1991)