- Email: [email protected]
A mixed integer programming approach for the single machine problem with unequal release dates
Computers & Operations Research 51 (2014) 323–330
Contents lists available at ScienceDirect
Computers & Operations Research journal homepage: www.elsevier.com/locate/caor
A mixed integer programming approach for the single machine problem with unequal release dates Anis Kooli a,n, Mehdi Serairi b a b
Ecole Supérieure des Sciences Economiques et Commerciales de Tunis, Université de Tunis. 4, Rue Abou Zakaria El Hafsi, 1089 Montfleury, Tunisie Université de Technologie de Compiègne, Heudiasyc, CNRS UMR 7253. Rue Roger Couttolenc, CS 60319, 60203 Compiègne Cedex, France
art ic l e i nf o
a b s t r a c t
Available online 18 June 2014
In this paper, we consider the problem of scheduling on a one-machine, a set of operations subject to unequal release dates with respect to the total completion time. This problem is known to be NP-hard in the strong sense. We propose an algorithm based on a Mixed Integer Linear Programming. This algorithm includes the implementation of a preprocessing procedure together with the consideration of valid inequalities. A computer simulation to measure the performance of the algorithm shows that our proposed method outperforms state-of-the-art branch-and-bound algorithms. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Mixed integer programming Preprocessing Valid inequalities Release dates
1. Introduction We investigate the single machine total completion time subject to unequal release dates. Let J ¼ f1; …; ng be a set of jobs, where each job j A J is available at its release date rj and has to be processed, without preemption, during its processing time pj. We seek for a schedule of the set J where only one job is processed at a time on the machine and the total completion time is minimized. Following the scheduling notation introduced by Graham et al. [22], this problem is denoted by 1jr j j∑C j where Cj is the completion time of job j; 8 j ¼ 1; …; n. Studying the 1jr j j∑C j problem is motivated by both practical and theoretical reasons. From a practical point of view, Conway et al. [11] stated that the single machine problem can represent many situations involving complex equipments as in the chemical, processing and mechanical manufacturing industries. Other realworld applications can be cited as in systems involving queuing in order to improve the service quality through minimizing the total flow time and the total waiting time. Actually, solving the foregoing problem is equivalent to minimize the total flow time and the total waiting time subject to release dates. Furthermore, many situations where job setup time is required, such request in database servers, can be modeled by the problem under consideration. This is an immediate consequence of the equivalence result, established by Gharbi et al. [17], between the 1jr j j∑C j and 1jr j ; sj j∑C j where sj denotes the setup time of job j. From a theoretical point of view, the simplicity and the ease of
n
Corresponding author. E-mail addresses: [email protected] (A. Kooli), [email protected] (M. Serairi). http://dx.doi.org/10.1016/j.cor.2014.06.013 0305-0548/& 2014 Elsevier Ltd. All rights reserved.
understanding of the problem are useful in suggesting future research directions. Moreover, it is frequently used as a subproblem in more complex problems, for example as a relaxation in more general flow shop problems (see [17]). This problem is known to be NP-hard in the strong sense (see [29]). However, well solvable cases exist. On one hand, when release dates are identical, the problem can be solved by the Shortest Processing Time (SPT) rule priority [32]. The SPT rule consists in processing jobs in nondecreasing order of their processing time. On the other hand, when preemption is allowed, the problem can be solved by the Shortest Remaining Processing Time priority rule (SRPT) [2]. The solution is obtained by scheduling, at a given time, the job having the shortest remaining processing time. The one machine total completion time with release dates has been the scope of numerous investigations. Thorough studies have addressed lower bounding strategies. Among the lower bounds proposed in the literature, those based on the SRPT solution are the best known ones. Ahmadi and Bagchi [1] proved that the SRPTbased lower bound outperforms several lower bounding strategies both in terms of efficiency and effectiveness. Later, Della Croce and T'Kindt [14] exploited the properties of the SRPT solution in order to propose an enhancement procedure computed by a dynamic programming. So far, several authors have investigated heuristic and metaheuristic approaches [4–6,9,23,30,15,21,31]. In particular, we mention the excellent Recovering Beam Search (RBS) heuristic of Della Croce and T'Kindt [15] which is based on a truncated branch-andbound algorithm without a backtracking procedure. Moreover, the authors considered a recovering procedure mainly based on interchange operations. At each level of the beam search, only one node is selected and then, in order to recover less interesting decisions that were taken in the higher level of the beam search,
324
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
variables indicating in which time instant a job will start are used. Another formulation proposed by Dyer and Wolsey [16] is based on the ordering between jobs, i.e., binary variables are used in order to indicate whether a job precedes another one or not. A survey of different mathematical models for a class of single machine problems and their performance in terms of exact resolution can be found in [26]. In our approach, we used the assignment-based formulation originally introduced by Lasserre and Queyranne [28]. The formulation is based on two types of variables:
the recovering procedure is applied on the corresponding partial solution. Jouglet et al. [21] introduced new dominance rules that were used to improve heuristics from the literature. Moreover, they proposed a Tabu Search method. The Tabu list is composed of the values of the total completion times of the last visited solutions. Thus, a solution is considered as “Tabu” if its corresponding total completion time belongs to the Tabu list. Furthermore, the authors considered a neighborhood based on the proposed dominance rules. Recently, Rakrouki et al. [31] proposed two new metaheuristics. The first one, named Genetic Local Search (GLS) Algorithm, consists of a combination between the genetic algorithm and the local search procedures that are based on constructive heuristics. The second metaheuristic, called the Genetic Recovering Beam Search (GRBS), consists of a hybridization between the GLS and the RBS heuristic of Della Croce and T'Kindt [15]. More precisely, a light version of GLS algorithm is invoked on the selected nodes of the RBS heuristic. The authors provide strong evidence that GRBS outperforms the state of the art heuristics. Exact methods based on branch-and-bound procedures were the scope of numerous papers [7,12,13,3,25,8,20,34]. Among these contributions, the most important ones are probably those of [8,20,34]. Indeed, while the first works [7,12,13,25] failed to solve medium sized instances, Chu [8] proposed a branch-and-bound algorithm which solved instances with up to 100 jobs. Later, T'Kindt et al. [34] revisited the latter algorithm by considering a dynamic programming dominance property introduced for a class of machine scheduling problems. Nearly at the same time, Jouglet et al. [20] proposed a different method which includes the nogood recording technique and new dominance rules based on the sets of scheduled and unscheduled jobs. These two branch-andbound algorithms showed almost the same performance and have improved the results of Chu [8]. In addition to the branch-and-bound based methods, several mixed integer programming formulations were proposed for single machine problems. Recently, Keha et al. [26] studied the performance of four different formulations on various single machine problems. The authors remarked that these formulations cannot solve large scale instances and are not competitive with specific designed procedures like branch-and-bound algorithms. The scope of this paper is to propose an easy to implement method, based on a previously proposed mathematical model. Our approach demonstrates good performance for solving large instances. To our aim, we introduce a preprocessing procedure in order to reduce the number of decision variables of the formulation. In addition, a set of valid inequalities are proposed to reduce the search space. Computational results show that our method is competitive with best known branch-and-bound procedures of the literature. The remainder of this paper is organized as follows. Section 2 recalls the assignment-based formulation that will be considered in this paper. In Section 3, we describe the preprocessing procedure. Two types of valid inequalities, namely position-based valid inequalities and SRPT-based valid inequalities, are detailed in Section 4. Computational experiments on a set of randomly generated instances are reported in Section 5. In Section 6, we present potential extensions of our approach to general scheduling problems. Finally, concluding remarks are given in Section 7.
The objective function (1) minimizes the total completion time. Constraints (2) and (3) state that each job is scheduled at only one position and at each position only one job is processed, respectively. In a feasible schedule, the job scheduled at position k cannot start before its release date and the completion time of the job at position k 1. Constraints (4) and (5) ensure these conditions. Finally, constraints (6) and (7) are the continuous and integrality constraints, respectively. As stated in [26], the considered formulation fails to solve large scale instances. One major key to enhance the performance of ðMIP0Þ is to perform a preprocessing procedure. Preprocessing could enhance the performance of mathematical formulation by reducing the number of decision variables. Actually, the number of the binary variables xjk is equal to n2. This can be considered as one of the reasons why this model cannot solve large scale instances. Indeed, the resolution of a MIP is based on a branch-and-cut procedure. Thus, the number of variables could have an impact on the performance of the resolution. In the next section, we describe the preprocessing procedure which will drastically reduce the number of such variables.
2. Mixed integer programming formulation
3. Preprocessing procedure
Several mathematical formulations were proposed for single machine problems. One of the first models is the time indexed formulation introduced by Sousa and Wolsey [33]. In this formulation, the planning horizon is subdivided into periods and binary
The aim of the preprocessing procedure is to fix some decision variables in order to reduce the computational burden of the proposed Mixed Integer Programming (MIP). Our preprocessing is based on computing the total completion time by taking into
Binary variables xjk determine the position of the job j in the schedule, i.e., xjk is equal to 1 if job j is scheduled at position k and 0 otherwise. Continuous variables ςk denote the completion time of the job scheduled at position k.
The mathematical model, hereafter denoted by ðMIP0Þ, is written as the following: n
∑ ςk
min
ð1Þ
k¼1 n
∑ xjk ¼ 1;
8 j ¼ 1; …; n
ð2Þ
8 k ¼ 1; …; n
ð3Þ
8 k ¼ 1; …; n
ð4Þ
ςk Zςk 1 þ ∑ pj xjk ;
8 k ¼ 2; …; n
ð5Þ
ςk Z0;
8 k ¼ 1; …; n
ð6Þ
xjk A f0; 1g;
8 j; k ¼ 1; …; n
ð7Þ
s:t:
k¼1 n
∑ xjk ¼ 1;
j¼1
n
ςk Z ∑ ðr j þ pj Þxjk ; j¼1
n
j¼1
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
account that a given job j is scheduled at position k, say LBnjk . If this latter value is greater than an upper bound on the total completion time then the decision variable xjk should be removed. It is noteworthy that calculating LBnjk can be time consuming. Therefore, in the sequel, we shall use preemptive relaxations instead of exact resolutions. Hence, we compute a lower bound on LBnjk , hereafter denoted by LBjk. This computation consists in computing a lower bound for the completion time of job j that is assumed to be scheduled at position k, denoted hereafter λjk, and lower bounds for the completion times of the jobs processed at position 0 k a k. In the sequel of this section, we start by estimating the cost λjk. After that, we show how to compute LBjk. Before proceeding further, let introduce the following notation:
Given a job j A J, S and Sj denote the SRPT sequence obtained on
the set J and J\fjg, respectively. Given a sequence σ, C k ðσÞ denotes the completion time of the job processed at position k in σ. Furthermore, we set C 0 ðσÞ to 0.
We begin by recalling the following lemma initially introduced by Chu [8]. Lemma 1 (Chu [8]). Ck(S) is a lower bound on the completion time of the job processed at position k in any sequence. From this Lemma, it is easy to see that maxðC k 1 ðSÞ; r j Þ þpj is a valid lower bound on the completion time of job j if it is processed at position k. However, a better estimation of λjk can be obtained by omitting job j when computing C k 1 ðSÞ. Thus, the following Lemma holds: Lemma 2. The quantity λjk ¼ maxðC k 1 ðSj Þ; r j Þ þ pj is a valid lower bound on the completion time of job j if it is processed at position k. Proof. Since job j is supposed to be processed at position k, it can be omitted from the computation of the SRPT rule in order to determine the completion time of the job at position k 1. □ We can now compute an estimation of the total completion time knowing that job j is scheduled at position k. First, a lower bound on the completion time of the job scheduled at position 0 k A f1; …; kg can be obtained by omitting job j from the computa0 tion of the SRPT sequence, namely Sj. Second, for k A fk þ1; …; ng, in addition to the sequence S, we explore an SRPT subsequence Skj of jobs J\fjg where the release dates are updated to maxðr i ; λjk Þ; 8 i A J\fjg. The following lemma resumes the computation of LBjk.
325
before λjk. Thus, a lower bound on the completion time of the job 0 processed at position k A fk þ 1; …; ng is C k0 k ðSkj Þ. □ As a consequence, we have the following straightforward corollary. Corollary 1. Let UB be an upper bound on the total completion time. If LBjk 4 UB, then in an optimal solution job j cannot be scheduled at position k. Therefore, the decision variable xjk should be removed. It should be noticed that we invoke Corollary 1 as a preprocessing procedure in order to detect non-allowed positions for each job jA J. A pseudocode of the procedure is described in Algorithm 1. Algorithm 1. Preprocessing procedure (Input UB). LBjk ’0 Compute the SRPT sequence S on the set J. for j¼1 to n do Compute the SRPT sequence Sj on the set J\fjg. for k ¼1 to n do 0 for k ¼ 1 to k 1 do LBjk ’LBjk þ C k0 ðSj Þ end for LBjk ’LBjk þ λjk Compute the SRPT sequence Skj on the set J\fjg starting from time λjk. 0 for k ¼ k þ 1 to n do LBjk ’LBjk þ maxðC k0 ðSÞ; C k0 k ðSkj ÞÞ end for if LBjk 4 UB then Remove xjk end if end for end for
Note that in our implementation we used the Recovering Beam Search heuristic of Della Croce and T'Kindt [15] with a beam width equal to 10. Before closing this section, let us now examine the time complexity associated with Algorithm 1. First, note that we invoke the SRPT procedure n2 times to determine all the Skj sequences. Since the SRPT procedure is computed on Oðn log nÞ-time, a running time of Oðn3 log nÞ follows in order to complete the preprocessing procedure.
0
Lemma 3. The value LBjk ¼ ∑nk0 ¼ 1 LBjk ðk Þ where 8 0 if k o k ðiÞ C 0 ðSj Þ > > < k 0 0 if k ¼ k ðiiÞ LBjk ðk Þ ¼ λjk > > : maxðC 0 ðSÞ; C 0 ðSk ÞÞ if k0 4 k ðiiiÞ k
k k
j
is a valid lower bound on the total completion time if job j is processed at position k. Proof. Since job j is supposed to be processed at position k, it can be omitted from the computation of the SRPT rule in order to determine the completion time of the job at position 0 0 k ; 8 k A f1; …; k 1g. It follows the result in (i). The second statement (ii) is an immediate result of Lemma 2. With regard to the third statement of the Lemma. On one hand, from Lemma 1, C k0 ðSÞ is a valid lower bound on the completion 0 time of a job processed at position k . On the other hand, after the processing of job j we seek for a sub-sequence of jobs that have to be processed after j. Since λjk is a valid lower bound on the completion time of job j, the jobs processed after j cannot start
4. Valid inequalities The preprocessing procedure described in the previous section allows us to eliminate a huge number of decision variables. Indeed, this procedure is able to remove about 78% of the decision variables (see Section 5.1). Another way in order to enhance the performance of the resolution of a MIP is to add cuts that will reduce the search space. In this section, we propose valid inequalities based on two different approaches. The first one is based on the lower bound on the completion time of a job j if it is processed at position k (λjk) described above. The second type of valid inequalities is based on the well-known improved preemptive lower bound of Della Croce and T'Kindt [14]. 4.1. Position-based valid inequalities The first type of valid inequalities is based on the fact that the completion time of the job at position k must be greater than a
326
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
○ Case (a): In an optimal solution, position ℓ þ 1 is occupied by a job scheduled after j in the SRPT sequence. Thus, improvements on the completion times of the jobs processed at positions ℓ and ℓ þ 1 are given by γ 1ℓ ¼ Ψ C ℓ ðSÞ and γ 2ℓ ¼ maxð0; minh A SuccðjÞ fmaxðr h ; Ψ Þ þ ph g C ℓ þ 1 ðSÞÞ, respectively. ○ Case (b): in an optimal solution, position ℓ þ 1 is occupied by a job scheduled before j in the SRPT sequence. Thus, improvements on the completion times of the jobs processed at position ℓ 1, ℓ and ℓ þ 1 are given by δ1ℓ ¼ maxð0; X C ℓ 1 ðSÞÞ; δ2ℓ ¼ maxð0; Y C ℓ ðSÞÞ and δ3ℓ ¼ max ð0; Z C ℓ þ 1 ðSÞÞ, respectively.
lower bound. Thus, we can use the estimation calculated on the preprocessing phase, namely λjk. The following proposition introduces the position-based valid inequalities. Proposition 1. n ςk Z ∑ λjk xjk ; 8 k ¼ 1; …; n
ð8Þ
j¼1
is a valid constraint for the model (1)–(7) Proof. Recall that ςk is the completion time of the job processed at position k. From the definition of λjk, ∑nj¼ 1 λjk xjk represents a lower bound on the completion time of the job processed at position k. □ It is note worthy to indicate that Haouari and Kharbeche [24] proposed similar valid inequalities for generic two machine flow shop problems denoted by F2 J ∑f j ðC j Þ, where f j ðÞ is a regular performance criterion. 4.2. SRPT-based valid inequalities The second type of valid inequalities is based on the SRPT properties exploited by Della Croce and T'Kindt [14] with the aim to improve the SRPT-based lower bound. Before introducing this type of valid inequalities, we shall recall the SRPT properties involved in the improved preemptive bound. In the sequel, we begin by developing the necessary properties. Then, we introduce the valid constraints. 4.2.1. SRPT properties For the sake of clarity, we adopt a slightly modified notation of Della Croce and T'Kindt [14]. Assume that job j A J is completed at ℓth position in the SRPT sequence, we denote by the following:
Pred(j) and Succ(j) the subset of predecessors and successors of job j in the SRPT sequence, respectively.
Ψ ¼ maxðC ℓ ðSÞ; C ℓ 1 ðSÞ þ pj Þ a lower bound with respect to C ℓ ðSÞ.
X ¼ C ℓ ðSj Þ ph ; Y ¼ maxfrj ; X; C ℓ 1 ðSÞg þ pj and Z ¼ maxfY; C ℓ
ðSÞg þ ph (where h is a job such that hA PredðjÞ and ph 4 pj ) a valid lower bound with respect to the completion times of the jobs at ℓ 1, ℓ and ℓ þ 1 in an optimal solution, respectively.
The improved SRPT bound [14] is based on determining a lower bound on the completion time increase of the job processed at a given position with respect to its SRPT completion time. The main idea of the procedure is based on the fact that, in an optimal solution, a job is scheduled either before, in the same position, or after its SRPT position. Thus, the following three cases may be distinguished:
Case 1: In an optimal solution, job j occupies a position k after
its SRPT position, i.e., k Zℓ þ 1. Then, α1ℓ ¼ C ℓ ðSj Þ C ℓ ðSÞ is the corresponding increase on the completion time of position ℓ. Furthermore, α2ℓ ¼ minfC ℓ ðSj Þ þ pj ; C ℓ þ 1 ðSj Þg C ℓ þ 1 ðSÞ is the increase in position ℓ þ 1. Case 2: In an optimal solution, job j occupies a position k before its SRPT position, i.e., k r ℓ 1. Then, β1ℓ ¼ maxf0; r j þ pj C ℓ þ 1 ðSÞ; C ℓ ðSÞ maxh A SuccðjÞ ph C ℓ 1 ðSÞg and β2ℓ ¼ max ð0; minh a j fmaxfr h ; C ℓ 1 ðSÞ; r j þpj ; C ℓ ðSÞ max i A PredðjÞ pi g þph g C ℓ ðSÞÞ are valid lower bounds on the increase of the completion times of the jobs processed at positions ℓ 1 and ℓ, respectively. Case 3: In an optimal solution, job j occupies a position k that corresponds to its SRPT position, i.e., k ¼ ℓ. Then, two subcases should be distinguished:
In order to simplify, hereafter, the presentation of the valid constraints, we introduce the fundamental Theorem of Della Croce and T'Kindt [14] in a slightly modified manner. Theorem 1 (Della Croce and T'Kindt [14]). Consider a job j completing in position 2 rℓ rn 1 in the SRPT sequence S. Let k be the position of job j in an optimal solution. We have the following results:
β1ℓ and δ1ℓ are valid lower bounds on the completion time increase
of the job processed at position ℓ 1 with respect to C ℓ 1 ðSÞ if k o ℓ and k ¼ ℓ (if the job scheduled at position ℓ þ 1 belongs to Pred(j)), respectively. β2ℓ ; minðγ 1ℓ ; δ2ℓ Þ and α1ℓ are valid lower bounds on the completion time increase of the job processed at position ℓ with respect to C ℓ ðSÞ if k o ℓ; k ¼ ℓ and k 4 ℓ, respectively. minðγ 2ℓ ; δ3ℓ Þ and α2ℓ are valid lower bounds on the completion time increase of the job processed at position ℓ þ 1 with respect to C ℓ þ 1 ðSÞ if k ¼ ℓ and k 4 ℓ, respectively.
Remark 1. As stated in Theorem 1, the corresponding increase in the completion times is valid for 2 rℓ r n 1. Therefore, Della Croce and T'Kindt [14] generalized the obtained results by taking into account the special cases where ℓ ¼1 and ℓ ¼ n. The authors reported that when ℓ ¼ 1 (respectively ℓ ¼ n) Cases 2 and 3b (respectively Cases 1 and 3b) do not hold. Moreover, when ℓ ¼ n, the Case 3a reduces to compute only γ1n. 4.2.2. Valid inequalities After having established the necessary background, we can now present the valid inequalities. The following proposition introduces valid inequalities related to position ℓ, ℓ þ1 and ℓ 1. Proposition 2. The following inequalities ςℓ Z C ℓ ðSÞ þ α1ℓ
ℓ1
n
∑
k ¼ ℓþ1
xjk þ β2ℓ ∑ xjk þ minðγ 1ℓ ; δ2ℓ Þxjℓ ; k¼1
8 ℓ ¼ 1; …; n ð9Þ
ςℓ þ 1 ZC ℓ þ 1 ðSÞ þ α2ℓ
n
∑
k ¼ ℓþ1
xjk þ minðγ 2ℓ ; δ3ℓ Þxjℓ ;
8 ℓ ¼ 1; …; n 1 ð10Þ
ℓ1
ςℓ 1 Z C ℓ 1 ðSÞ þ β1ℓ ∑ xjk þ δ1ℓ ðxjℓ þ k¼1
∑
h A PredðjÞ
xh;ℓ þ 1 1Þ;
8 ℓ ¼ 2; …; n
ð11Þ are valid for the model (MIP0) given by (1)–(7). Proof. Let k be the position where job j completes its execution in an optimal solution. From Theorem 1, we know that α1ℓ ; β2ℓ and minðγ 1ℓ ; δ2ℓ Þ are valid lower bounds on the completion time increase of the job in position ℓ with respect to C ℓ ðSÞ if k o ℓ; k ¼ ℓ and k 4 ℓ, respectively. Thus, the RHS of (9) constitutes a valid lower bound on the completion time of the job processed at position ℓ. Applying the same argument, constraints (10) and (11) hold. □
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
Remark 2. Taking into account the cases described in Remark 1, we shall initialize certain quantities as follows. On one hand, if ℓ ¼ 1, then β11 ¼ β21 ¼ β31 ¼ 0 and δ11 ¼ δ21 ¼ δ31 ¼ 1 . On the other hand, if ℓ ¼ n, then α1n ¼ α2n ¼ γ 2n ¼ 0 and δ1n ¼ δ2n ¼ δ3n ¼ 1: Thus, we must consider these special cases in the implementation of constraints (9)–(11). Remark 3. In inequality (11), ðxjℓ þ ∑h A PredðjÞ xh;ℓ þ 1 1Þ could be less than 0, it can be interesting to use inequality (11) with the following inequality: ℓ1
ςℓ 1 Z C ℓ 1 ðSÞ þ β1ℓ ∑ xjk ; k¼1
8 ℓ ¼ 2; …; n
ð12Þ
5. Computational results In this section, we report the results of a computational study that aims at assessing the performance of the proposed approach. To that aim, we carried out three sets of experiments. In the first set, we analyze the impact of the preprocessing on the original formulation. In the second set, we investigate the best combination of valid constraints that offers a good trade-off between tightness and efficiency. Finally, in the third set, we compare our best method with state-of-the-art exact methods. The tests were conducted on a set of generated instances following the scheme described in [25]. The number of jobs n was taken to be equal to 20, 40, 60, 80, 100, 120, 140, 160, 180 and 200 jobs. The processing times were randomly generated using the discrete uniform distribution on ½1; 100. The release dates were uniformly drawn from ½0; 50:5nR, where R controls the range of the distribution and belongs to the set f0:2; 0:4; 0:6; 0:8; 1; 1:25; 1:5; 1:75; 2:0; 3:0g. For each combination of n and R, five instances were generated. Hence, 50 instances are generated for each value of n. All the experiments were carried out on an Intel Xeon E5620 2.4 GHz processor, with 8 GB RAM memory, running Windows 7. The tested algorithms were coded in C and compiled under Microsoft Visual Studio 2010 using CPLEX 12.2. For all the experiments, the time limit for executing the procedures was set to 3600 s.
5.1. Impact of the preprocessing In this first experiment, we analyzed the impact of the preprocessing on the formulation. To this end, we tested the following two versions:
F0: The model (MIP0). F1: The model (MIP0) with the preprocessing procedure (see Section 3).
The results are displayed in Table 1. For each version, we provide the average required CPU time in seconds. The numbers between parenthesis are the number of unsolved instances. We note that we stopped the resolution of F0 to 60-jobs instances since it is clear that this version fails to solve larger instances. We see from Table 1 that the original formulation F0 cannot solve instances with up to 40 jobs. We remark that formulation F1 greatly improves the results since it is able to solve 78% of the 80-jobs instances. This is due to the fact that the preprocessing is able to eliminate a large number of decision variables xjk. In fact, the tests show that this procedure is capable of removing 78.24% of the decision variables on average.
327
Table 1 Impact of the preprocessing. n
F0
F1
20 40 60 80 100
4.53 2578.38 (33) 3276.47 (44) – –
0.04 0.88 189.82 (1) 946.49 (11) 2047.42 (26)
Table 2 Impact of the valid constraints. n
F1
F2
F3
F4
20 40 60 80 100
0.04 0.88 189.82 (1) 946.49 (11) 2047.42 (26)
0.02 0.24 1.16 4.49 81.91
0.05 0.39 1.18 4.44 54.12
0.02 0.26 0.76 2.69 17.62
5.2. Impact of the valid constraints In this second phase of experiments, we investigated the impact of the constraints on the formulation. We tested the following versions:
F2: Formulation F1 along with position-based constraints (8). F3: Formulation F1 along with SRPT-based constraints (9)–(11). F4: Formulation F1 after appending all valid constraints. We first reported in Table 2 the performance of F1, F2, F3 and F4 and then observed the following points:
While F1 fails to solve 52% of 100-jobs instances, the three
versions F2, F3 and F4 were able to solve all the instances from 20 to 100 jobs within the time limit of 1 h. The best performance is given by appending all the new cuts. Actually, F4 is less time consuming than F2 and F3.
We stopped the simulation at n ¼100 as we are confident that the proposed constraints improve the performance of the resolution. Therefore, we assumed in our implementation that all the new cuts are appended to MIP0. Thus, F4 will constitute our reference approach in order to compare it with methods of the literature. 5.3. Comparison with state-of-the-art exact methods In this third set of experiments, we compare the performance of F4 with respect to state-of-the-art exact methods. It is worth noting that the best known exact methods for this problem are those of Jouglet et al. [20] and T'Kindt et al. [34]. The two methods are based on a branch-and-bound strategy. In order to compare our formulation with these methods, we developed a branch-andbound algorithm based on the strategies and dominance rules found in the latter papers. More precisely, the branch-and-bound algorithm that we implemented is based on these concepts:
Best active node search strategy: T'Kindt et al. [34] compared
different search strategies for this problem and found that this strategy is the best one. SRPT rule for the computation of lower bounds on nodes. Dynamic programming dominance rule as described in [34]. Dominance rules relying on scheduled and unscheduled jobs introduced by Jouglet et al. [20]. It is worth noting that these rules increase the average resolution time but in return they
328
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
are able to prune more nodes. Thus, the algorithm is able to solve larger instances. In the following, this branch-and-bound procedure is denoted by BB. Moreover, in order to get a better idea between the developed method in this paper and the branch-and-bound algorithm, we tested a version of F4 with only one thread hereafter denoted by F41Th . Indeed, CPLEX solves the formulation using as many threads as the number of available cores in the testing machine, i.e., in our case, eight threads were launched. We added this version in order to have a “fair” comparison between the mathematical formulation and the branch-and-bound method. For each procedure, we provide
Time: The average CPU time in seconds. #Unsolved: The number of unsolved instances. #Best: The number of times the procedure achieves the best value among all other procedures.
F4 exhibits an excellent performance both in terms of efficiency
n
ζ k Z ςk ∑ dj xjk ; j¼1
In this phase, we focused our attention to larger instances from 120 to 200 jobs. The results are presented in Table 3. From Table 3, we made the following observations:
scheduling problems. For the sake of clarity, we discuss this foregoing question on two single machine scheduling problems subject to release dates where the objective is to minimize the total tardiness and the total weighted tardiness. Before proceeding further, let us denote by wj, dj and T j ¼ maxð0; C j dj Þ the weight, the due date and the tardiness of job jA J, respectively. Following these notations, the former problems are denoted by 1jr j j∑T j and 1jr j j∑wj T j and will be referred to, in the sequel, as the P problems. The assignment-based model (1)–(7) cannot be used as it is in the P problems since a valid instance of these problems requires some decision variables and constraints to be added. For 1jr j j∑T j , we introduce the semi-continuous decision variables ζk that represent the tardiness of the job processed at position k; 8 k A f1; …; ng. The objective function (1) is replaced by min ∑nk ¼ 1 ζ k . In addition, the following set of constraints must be added:
and effectiveness. Indeed, it consistently exhibits the best average time of resolution and it solves more instances. Interestingly, we remark that the version with one thread F41Th is competitive with the branch-and-bound procedure and even manages to solve larger instances. Surprisingly, for n ¼200, F41Th yields the best value more times than F4 and BB. This observation can be explained by “the variability performance” (see [27]). Actually, on average, using several threads may lead to enhancing the solver performance. However, on a single instance, the performance may arbitrarily depend on the number of threads.
Pushing our analysis a step further, we performed a detailed analysis of the results according to the parameter R. Table 4 presents the mean CPU time on all instances (from 120 to 200 jobs) for a fixed value of R. The total number of unsolved instances is shown between parenthesis. We remark from Table 4 that the hard instances come for the range of values 0.4–1.25. For these instances, the two versions of our method F4 and F41Th dominate the branch-and-bound procedure BB. For the other easier instances, we noticed that BB is competitive with the proposed methods. Specially for the instances with the value 1.75 and higher, the branch-and-bound is more efficient than our methods.
6. Extension for general scheduling problems In this section, we discuss how the different components of our approach should be adapted to be used for other general
8 k ¼ 1; …; n
On the contrary, to represent the tardiness of a job in the 1jr j j∑wj T j problem, one needs to introduce new decision variables. Indeed, for the 1jr j j∑T j problem, we can consider the tardiness of the job processed at position k instead of the tardiness of the job itself. However, this is not possible for the 1jr j j∑wj T j problem since both the job's weight and the tardiness are considered. Let us denote by θj the tardiness of the job j; 8 j A f1; …; ng. The objective function (1) is consequently replaced by min ∑nj¼ 1 wj θj and the following constraints are added: θj Z ςk dj Mð1 xjk Þ;
8 j; k ¼ 1; …; n
where M is an upper bound on how large ςk can be. Unfortunately, we face a big-M model. It is worth noting that this type of models is known to be the hardest to solve (see [10]). In order to remove the big-M dependency, a possible solution consists in proposing alternative inequalities. We refer the reader to Keha et al. [26] who proposed alternative inequalities for a set of single machine problems without release dates. Another solution could be the use of the automatic problem reformulation [10]. These solutions, however, need both a theoretical and a computational thorough investigation. The lower bound on the completion time of job j if it is processed at position k (λjk) is one of the most important parts of our algorithm. Indeed, this lower bound allows us to derive the preprocessing procedure and valid inequalities (8). Recall that, the lower bounds λjk is a valid lower bound for any sequence. This result is a straightforward consequence of Lemmas 1 and 2. Thus, inequalities (8) remain valid for the P problems. Concerning inequalities (9)–(11), recall also that Cl(S) is a valid lower bound on the completion time of the job processed at position k (see Lemma 1). Therefore, from Theorem 1 and Proposition 2, we can conclude that these inequalities are also valid for the P problems. In order to adapt the preprocessing procedure, we must propose a lower bound on the total (weighted) tardiness while
Table 3 Comparison of methods. n
120 140 160 180 200
F4
BB
F41Th
Time (s)
#Unsolved
#Best
Time (s)
#Unsolved
#Best
Time (s)
#Unsolved
#Best
51.09 38.71 309.89 437.02 712.00
0 0 3 4 7
50 50 49 49 48
76.26 89.53 520.40 893.90 1310.22
0 0 4 8 16
50 50 47 42 36
163.60 116.59 454.24 572.12 876.65
1 0 4 5 9
50 50 48 45 49
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
329
Table 4 Detailed results. R
0.2
0.4
0.6
0.8
1
1.25
1.5
1.75
2
3
F4 BB F41Th
52.82 60.33 85.46
478.15 (2) 569.35 (3) 769.22 (2)
1039.09 (5) 1289.96 (5) 1383.48 (7)
1160.03 (6) 2023.95 (12) 1529.46 (8)
255.93 (1) 1559.13 (7) 421.17 (1)
87.36 255.91 (1) 156.30 (1)
8.24 15.22 6.87
6.44 3.99 5.66
5.28 1.74 4.83
4.08 1.04 3.92
taking into account that a given job j is scheduled at position k. To this end, we use the lower bound proposed by Haouari and Kharbeche [24] that consists in solving the linear assignment problem. In our problem, the assignment costs are κ j0 k0 ¼ 0 wj0 maxðλj0 k0 dj0 ; 0Þ; 8 j0 ; k A f1; …; ng (it is obvious that the weights wj are equal to 1 for 1jr j j∑T j ). In order to ensure that the solution 0 will take job j at position k, we set to infinity the costs κ jk0 ; 8 k a k 0 0 and the costs κ j k ; 8 j a j. As described in Section 3, if this value is greater than an upper bound on the total (weighted) tardiness, then a decision variable can be eliminated from the model. As the linear assignment problem is solved in Oðn3 Þtime and will be repeated for each job j A f1; …; ng and each position k A f1; …; ng, the total complexity of the procedure is Oðn5 Þ. This leads to investigate more appropriate lower bounds on the (weighted) tardiness of job j scheduled at position k. Finally, we note that the approach developed in this paper can be extended to permutation flowshop problems where the objective function is to minimize total completion time or total tardiness problem. For these problems we think that the assignmentbased formulation could be very useful. The weighted version of these problems needs a deeper analysis to overcome the big-M problem. Moreover many components of our approach could be used for other general scheduling problems. However, to achieve a good performance, we believe that it is very important to exploit the characteristics of the problem when designing a method. For example, for the two-machine total completion time flowshop problem, we noticed that Haouari and Kharbeche [24] used the Johnson's rule [19] to calculate the λjk values and derived a valid inequality to strengthen the continuous relaxation of a mathematical model proposed by Hoogeveen et al. [18]. These values can be used in our preprocessing procedure to reduce the binary variables of the mathematical model. The valid inequality along with the preprocessing procedure should enhance the performance of the exact resolution of the mathematical model. In a nutshell, two important workarounds would prove to be useful to achieve good performances of our approach for general scheduling problems. The first is to use suitable lower bounds for each studied problem to be able to eliminate decision variables. The second issue consists in deriving more tailored valid inequalities.
7. Conclusion In this paper, we presented a new mathematical approach to solve the one machine scheduling problem subject to release dates. First, we presented a preprocessing procedure which aims to reduce the number of binary variables from the mathematical model. Then, we developed two types of valid inequalities. The first set of inequalities is deduced from the computation of lower bounds of the completion times of the jobs scheduled at these positions. The second set is based on improvements on completion times of the Shortest Remaining Processing Time sequence. Computational results showed that our simple method through the use of a commercial solver is competitive with state-of-the-art branch-and-bound algorithms and is able to solve large instances up to 140 jobs.
Acknowledgments Many thanks to Professor Vincent T'Kindt for valuable discussion about the branch-and-bound procedure. Also, the authors would like to thank the anonymous referees for their insightful comments that substantially improved the quality of the paper. References [1] Ahmadi RH, Bagchi U. Lower bounds for single-machine scheduling problems. Naval Res Logist 1990;37:967–79. [2] Baker KR. Introduction to sequencing and scheduling. New York: Wiley; 1974. [3] Bianco L, Ricciardelli S. Scheduling of a single machine to minimize total weighted completion time subject to release dates. Naval Res Logist Q 1982;29:151–67. [4] Chand S, Traub R, Uzsoy R. Single-machine scheduling with dynamic arrivals: decomposition results and an improved algorithm. Naval Res Logist 1996;43:709–19. [5] Chand S, Traub R, Uzsoy R. An iterative heuristic for the single-machine dynamic total completion time scheduling problem. Comput Oper Res 1996;23:641–51. [6] Chand S, Traub R, Uzsoy R. Rolling horizon procedures for the single machine deterministic total completion time scheduling problem with release dates. Ann Oper Res 1997;70:115–25. [7] Chandra R. On n=1=F dynamic deterministic systems. Naval Res Logist Q 1979;26:537–44. [8] Chu C. A branch-and-bound algorithm to minimize total tardiness with different release dates. Naval Res Logist 1992;39:265–83. [9] Chu C. Efficient heuristics to minimize total flow time with release dates. Oper Res Lett 1992;12:321–30. [10] Codato G, Fischetti M. Combinatorial benders' cuts for mixed-integer linear programming. Oper Res 2006;54:756–66. [11] Conway RW, Maxwell WL, Miller LW. Theory of scheduling. MA: AddisonWesley Publishing Company; 1967. [12] Dessouky MI, Deogun JS. Sequencing jobs with unequal ready times to minimize mean flow time. SIAM J Comput 1981;10:192–202. [13] Deogun JS. On scheduling with ready times to minimize mean flow time. Comput J 1983;26:320–8. [14] Della Croce F, T'Kindt V. Improving the preemptive bound for the onemachine dynamic total completion time scheduling problem. Oper Res Lett 2003;31:142–8. [15] Della Croce F, T'Kindt V. A recovering beam search algorithm for the onemachine dynamic total completion time scheduling problem. J Oper Res Soc 2002;53:1275–80. [16] Dyer ME, Wolsey LA. Formulating the single machine sequencing problem with release dates as a mixed integer program. Discr Appl Math 1990;26: 255–70. [17] Gharbi A, Ladhari T, Msakni MK, Serairi M. The two-machine flowshop scheduling problem with sequence-independent setup times: new lower bounding strategies. Eur J Oper Res 2013;231:69–78. [18] Hoogeveen H, Van Norden L, Van de Velde S. Lower bounds for minimizing total completion time in a two-machine flow shop. J Sched 2006;9:559–68. [19] Johnson S. Optimal two- and three-stage production schedules with set up times included. Naval Res Logist Q 1954;1:61–8. [20] Jouglet A, Baptiste P, Carlier J. Branch-and-bound algorithms for total weighted tardiness. In: Leung JYT, editor. Handbook of scheduling: algorithms, models and performance analysis; 2004. [21] Jouglet A, Savourey D, Carlier J, Baptiste P. Dominance-based heuristics for one-machine total cost scheduling problems. Eur J Oper Res 2008;184: 879–99. [22] Graham RL, Lawler EL, Lenstra JK, Rinnooy Kan AHG. Optimization and approximation in deterministic sequencing and scheduling: a survey. Ann Discr Math 1979;5:287–326. [23] Guo Y, Lim A, Rodriguez B, Yu S. Minimizing total flow time in single machine environment with release time: an experimental analysis. Comput Ind Eng 2004;47:123–40. [24] Haouari M, Kharbeche M. An assignment-based lower bound for a class of two-machine flow shop problems. Comput Oper Res 2013;40:1693–9. [25] Hariri AMA, Potts CN. An algorithm for single machine sequencing with release dates to minimize total weighted completion time. Discr Appl Math 1983;5:99–109.
330
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
[26] Keha AB, Khowala K, Fowler JW. Mixed integer programming formulations for single machine scheduling problems. Comput Ind Eng 2009;56:357–67. [27] Koch T, Achterberg T, Andersen E, Bastert O, Berthold T, Bixby RE, et al. MIPLIB 2010. Math Program Comput 2011;3:103–63. [28] Lasserre JB, Queyranne M. Generic scheduling polyhedra and a new mixedinteger formulation for single-machine scheduling. In: Proceedings of the second IPCO conference. Pittsburgh: Carnegie-Mellon University; 1992. p. 136–49. [29] Lawler EL, Lenstra JK, Rinnooy Kan AHG, Shmoys DB. Sequencing and scheduling: algorithms and complexity. In: Graves SC, Rinnooy Kan AHG, Zipkin P, editors. Handbooks in operations research and management science, Logistics of production and inventory, vol. 4. Amsterdam: North-Holland; 1993.
[30] Liu J, MacCarthy BL. Effective heuristics for the single machine sequencing problem with ready times. Int J Prod Res 1991;29:1521–33. [31] Rakrouki MA, Ladhari T, T'Kindt V. Coupling genetic local search and recovering beam search algorithms for minimizing the total completion time in the single machine scheduling problem subject to release dates. Comput Oper Res 2012;39:1257–64. [32] Smith WE. Various optimizers for single-state production. Naval Res Logist Q 1956;3:59–66. [33] Sousa JP, Wolsey LA. A time indexed formulation of non-preemptive single machine scheduling problems. Math Program 1992;54:353–67. [34] T'Kindt V, Della Croce F, Esswein C. Revisiting branch and bound search strategies for machine scheduling problems. J Sched 2004;7:429–40.
Contents lists available at ScienceDirect
Computers & Operations Research journal homepage: www.elsevier.com/locate/caor
A mixed integer programming approach for the single machine problem with unequal release dates Anis Kooli a,n, Mehdi Serairi b a b
Ecole Supérieure des Sciences Economiques et Commerciales de Tunis, Université de Tunis. 4, Rue Abou Zakaria El Hafsi, 1089 Montfleury, Tunisie Université de Technologie de Compiègne, Heudiasyc, CNRS UMR 7253. Rue Roger Couttolenc, CS 60319, 60203 Compiègne Cedex, France
art ic l e i nf o
a b s t r a c t
Available online 18 June 2014
In this paper, we consider the problem of scheduling on a one-machine, a set of operations subject to unequal release dates with respect to the total completion time. This problem is known to be NP-hard in the strong sense. We propose an algorithm based on a Mixed Integer Linear Programming. This algorithm includes the implementation of a preprocessing procedure together with the consideration of valid inequalities. A computer simulation to measure the performance of the algorithm shows that our proposed method outperforms state-of-the-art branch-and-bound algorithms. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Mixed integer programming Preprocessing Valid inequalities Release dates
1. Introduction We investigate the single machine total completion time subject to unequal release dates. Let J ¼ f1; …; ng be a set of jobs, where each job j A J is available at its release date rj and has to be processed, without preemption, during its processing time pj. We seek for a schedule of the set J where only one job is processed at a time on the machine and the total completion time is minimized. Following the scheduling notation introduced by Graham et al. [22], this problem is denoted by 1jr j j∑C j where Cj is the completion time of job j; 8 j ¼ 1; …; n. Studying the 1jr j j∑C j problem is motivated by both practical and theoretical reasons. From a practical point of view, Conway et al. [11] stated that the single machine problem can represent many situations involving complex equipments as in the chemical, processing and mechanical manufacturing industries. Other realworld applications can be cited as in systems involving queuing in order to improve the service quality through minimizing the total flow time and the total waiting time. Actually, solving the foregoing problem is equivalent to minimize the total flow time and the total waiting time subject to release dates. Furthermore, many situations where job setup time is required, such request in database servers, can be modeled by the problem under consideration. This is an immediate consequence of the equivalence result, established by Gharbi et al. [17], between the 1jr j j∑C j and 1jr j ; sj j∑C j where sj denotes the setup time of job j. From a theoretical point of view, the simplicity and the ease of
n
Corresponding author. E-mail addresses: [email protected] (A. Kooli), [email protected] (M. Serairi). http://dx.doi.org/10.1016/j.cor.2014.06.013 0305-0548/& 2014 Elsevier Ltd. All rights reserved.
understanding of the problem are useful in suggesting future research directions. Moreover, it is frequently used as a subproblem in more complex problems, for example as a relaxation in more general flow shop problems (see [17]). This problem is known to be NP-hard in the strong sense (see [29]). However, well solvable cases exist. On one hand, when release dates are identical, the problem can be solved by the Shortest Processing Time (SPT) rule priority [32]. The SPT rule consists in processing jobs in nondecreasing order of their processing time. On the other hand, when preemption is allowed, the problem can be solved by the Shortest Remaining Processing Time priority rule (SRPT) [2]. The solution is obtained by scheduling, at a given time, the job having the shortest remaining processing time. The one machine total completion time with release dates has been the scope of numerous investigations. Thorough studies have addressed lower bounding strategies. Among the lower bounds proposed in the literature, those based on the SRPT solution are the best known ones. Ahmadi and Bagchi [1] proved that the SRPTbased lower bound outperforms several lower bounding strategies both in terms of efficiency and effectiveness. Later, Della Croce and T'Kindt [14] exploited the properties of the SRPT solution in order to propose an enhancement procedure computed by a dynamic programming. So far, several authors have investigated heuristic and metaheuristic approaches [4–6,9,23,30,15,21,31]. In particular, we mention the excellent Recovering Beam Search (RBS) heuristic of Della Croce and T'Kindt [15] which is based on a truncated branch-andbound algorithm without a backtracking procedure. Moreover, the authors considered a recovering procedure mainly based on interchange operations. At each level of the beam search, only one node is selected and then, in order to recover less interesting decisions that were taken in the higher level of the beam search,
324
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
variables indicating in which time instant a job will start are used. Another formulation proposed by Dyer and Wolsey [16] is based on the ordering between jobs, i.e., binary variables are used in order to indicate whether a job precedes another one or not. A survey of different mathematical models for a class of single machine problems and their performance in terms of exact resolution can be found in [26]. In our approach, we used the assignment-based formulation originally introduced by Lasserre and Queyranne [28]. The formulation is based on two types of variables:
the recovering procedure is applied on the corresponding partial solution. Jouglet et al. [21] introduced new dominance rules that were used to improve heuristics from the literature. Moreover, they proposed a Tabu Search method. The Tabu list is composed of the values of the total completion times of the last visited solutions. Thus, a solution is considered as “Tabu” if its corresponding total completion time belongs to the Tabu list. Furthermore, the authors considered a neighborhood based on the proposed dominance rules. Recently, Rakrouki et al. [31] proposed two new metaheuristics. The first one, named Genetic Local Search (GLS) Algorithm, consists of a combination between the genetic algorithm and the local search procedures that are based on constructive heuristics. The second metaheuristic, called the Genetic Recovering Beam Search (GRBS), consists of a hybridization between the GLS and the RBS heuristic of Della Croce and T'Kindt [15]. More precisely, a light version of GLS algorithm is invoked on the selected nodes of the RBS heuristic. The authors provide strong evidence that GRBS outperforms the state of the art heuristics. Exact methods based on branch-and-bound procedures were the scope of numerous papers [7,12,13,3,25,8,20,34]. Among these contributions, the most important ones are probably those of [8,20,34]. Indeed, while the first works [7,12,13,25] failed to solve medium sized instances, Chu [8] proposed a branch-and-bound algorithm which solved instances with up to 100 jobs. Later, T'Kindt et al. [34] revisited the latter algorithm by considering a dynamic programming dominance property introduced for a class of machine scheduling problems. Nearly at the same time, Jouglet et al. [20] proposed a different method which includes the nogood recording technique and new dominance rules based on the sets of scheduled and unscheduled jobs. These two branch-andbound algorithms showed almost the same performance and have improved the results of Chu [8]. In addition to the branch-and-bound based methods, several mixed integer programming formulations were proposed for single machine problems. Recently, Keha et al. [26] studied the performance of four different formulations on various single machine problems. The authors remarked that these formulations cannot solve large scale instances and are not competitive with specific designed procedures like branch-and-bound algorithms. The scope of this paper is to propose an easy to implement method, based on a previously proposed mathematical model. Our approach demonstrates good performance for solving large instances. To our aim, we introduce a preprocessing procedure in order to reduce the number of decision variables of the formulation. In addition, a set of valid inequalities are proposed to reduce the search space. Computational results show that our method is competitive with best known branch-and-bound procedures of the literature. The remainder of this paper is organized as follows. Section 2 recalls the assignment-based formulation that will be considered in this paper. In Section 3, we describe the preprocessing procedure. Two types of valid inequalities, namely position-based valid inequalities and SRPT-based valid inequalities, are detailed in Section 4. Computational experiments on a set of randomly generated instances are reported in Section 5. In Section 6, we present potential extensions of our approach to general scheduling problems. Finally, concluding remarks are given in Section 7.
The objective function (1) minimizes the total completion time. Constraints (2) and (3) state that each job is scheduled at only one position and at each position only one job is processed, respectively. In a feasible schedule, the job scheduled at position k cannot start before its release date and the completion time of the job at position k 1. Constraints (4) and (5) ensure these conditions. Finally, constraints (6) and (7) are the continuous and integrality constraints, respectively. As stated in [26], the considered formulation fails to solve large scale instances. One major key to enhance the performance of ðMIP0Þ is to perform a preprocessing procedure. Preprocessing could enhance the performance of mathematical formulation by reducing the number of decision variables. Actually, the number of the binary variables xjk is equal to n2. This can be considered as one of the reasons why this model cannot solve large scale instances. Indeed, the resolution of a MIP is based on a branch-and-cut procedure. Thus, the number of variables could have an impact on the performance of the resolution. In the next section, we describe the preprocessing procedure which will drastically reduce the number of such variables.
2. Mixed integer programming formulation
3. Preprocessing procedure
Several mathematical formulations were proposed for single machine problems. One of the first models is the time indexed formulation introduced by Sousa and Wolsey [33]. In this formulation, the planning horizon is subdivided into periods and binary
The aim of the preprocessing procedure is to fix some decision variables in order to reduce the computational burden of the proposed Mixed Integer Programming (MIP). Our preprocessing is based on computing the total completion time by taking into
Binary variables xjk determine the position of the job j in the schedule, i.e., xjk is equal to 1 if job j is scheduled at position k and 0 otherwise. Continuous variables ςk denote the completion time of the job scheduled at position k.
The mathematical model, hereafter denoted by ðMIP0Þ, is written as the following: n
∑ ςk
min
ð1Þ
k¼1 n
∑ xjk ¼ 1;
8 j ¼ 1; …; n
ð2Þ
8 k ¼ 1; …; n
ð3Þ
8 k ¼ 1; …; n
ð4Þ
ςk Zςk 1 þ ∑ pj xjk ;
8 k ¼ 2; …; n
ð5Þ
ςk Z0;
8 k ¼ 1; …; n
ð6Þ
xjk A f0; 1g;
8 j; k ¼ 1; …; n
ð7Þ
s:t:
k¼1 n
∑ xjk ¼ 1;
j¼1
n
ςk Z ∑ ðr j þ pj Þxjk ; j¼1
n
j¼1
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
account that a given job j is scheduled at position k, say LBnjk . If this latter value is greater than an upper bound on the total completion time then the decision variable xjk should be removed. It is noteworthy that calculating LBnjk can be time consuming. Therefore, in the sequel, we shall use preemptive relaxations instead of exact resolutions. Hence, we compute a lower bound on LBnjk , hereafter denoted by LBjk. This computation consists in computing a lower bound for the completion time of job j that is assumed to be scheduled at position k, denoted hereafter λjk, and lower bounds for the completion times of the jobs processed at position 0 k a k. In the sequel of this section, we start by estimating the cost λjk. After that, we show how to compute LBjk. Before proceeding further, let introduce the following notation:
Given a job j A J, S and Sj denote the SRPT sequence obtained on
the set J and J\fjg, respectively. Given a sequence σ, C k ðσÞ denotes the completion time of the job processed at position k in σ. Furthermore, we set C 0 ðσÞ to 0.
We begin by recalling the following lemma initially introduced by Chu [8]. Lemma 1 (Chu [8]). Ck(S) is a lower bound on the completion time of the job processed at position k in any sequence. From this Lemma, it is easy to see that maxðC k 1 ðSÞ; r j Þ þpj is a valid lower bound on the completion time of job j if it is processed at position k. However, a better estimation of λjk can be obtained by omitting job j when computing C k 1 ðSÞ. Thus, the following Lemma holds: Lemma 2. The quantity λjk ¼ maxðC k 1 ðSj Þ; r j Þ þ pj is a valid lower bound on the completion time of job j if it is processed at position k. Proof. Since job j is supposed to be processed at position k, it can be omitted from the computation of the SRPT rule in order to determine the completion time of the job at position k 1. □ We can now compute an estimation of the total completion time knowing that job j is scheduled at position k. First, a lower bound on the completion time of the job scheduled at position 0 k A f1; …; kg can be obtained by omitting job j from the computa0 tion of the SRPT sequence, namely Sj. Second, for k A fk þ1; …; ng, in addition to the sequence S, we explore an SRPT subsequence Skj of jobs J\fjg where the release dates are updated to maxðr i ; λjk Þ; 8 i A J\fjg. The following lemma resumes the computation of LBjk.
325
before λjk. Thus, a lower bound on the completion time of the job 0 processed at position k A fk þ 1; …; ng is C k0 k ðSkj Þ. □ As a consequence, we have the following straightforward corollary. Corollary 1. Let UB be an upper bound on the total completion time. If LBjk 4 UB, then in an optimal solution job j cannot be scheduled at position k. Therefore, the decision variable xjk should be removed. It should be noticed that we invoke Corollary 1 as a preprocessing procedure in order to detect non-allowed positions for each job jA J. A pseudocode of the procedure is described in Algorithm 1. Algorithm 1. Preprocessing procedure (Input UB). LBjk ’0 Compute the SRPT sequence S on the set J. for j¼1 to n do Compute the SRPT sequence Sj on the set J\fjg. for k ¼1 to n do 0 for k ¼ 1 to k 1 do LBjk ’LBjk þ C k0 ðSj Þ end for LBjk ’LBjk þ λjk Compute the SRPT sequence Skj on the set J\fjg starting from time λjk. 0 for k ¼ k þ 1 to n do LBjk ’LBjk þ maxðC k0 ðSÞ; C k0 k ðSkj ÞÞ end for if LBjk 4 UB then Remove xjk end if end for end for
Note that in our implementation we used the Recovering Beam Search heuristic of Della Croce and T'Kindt [15] with a beam width equal to 10. Before closing this section, let us now examine the time complexity associated with Algorithm 1. First, note that we invoke the SRPT procedure n2 times to determine all the Skj sequences. Since the SRPT procedure is computed on Oðn log nÞ-time, a running time of Oðn3 log nÞ follows in order to complete the preprocessing procedure.
0
Lemma 3. The value LBjk ¼ ∑nk0 ¼ 1 LBjk ðk Þ where 8 0 if k o k ðiÞ C 0 ðSj Þ > > < k 0 0 if k ¼ k ðiiÞ LBjk ðk Þ ¼ λjk > > : maxðC 0 ðSÞ; C 0 ðSk ÞÞ if k0 4 k ðiiiÞ k
k k
j
is a valid lower bound on the total completion time if job j is processed at position k. Proof. Since job j is supposed to be processed at position k, it can be omitted from the computation of the SRPT rule in order to determine the completion time of the job at position 0 0 k ; 8 k A f1; …; k 1g. It follows the result in (i). The second statement (ii) is an immediate result of Lemma 2. With regard to the third statement of the Lemma. On one hand, from Lemma 1, C k0 ðSÞ is a valid lower bound on the completion 0 time of a job processed at position k . On the other hand, after the processing of job j we seek for a sub-sequence of jobs that have to be processed after j. Since λjk is a valid lower bound on the completion time of job j, the jobs processed after j cannot start
4. Valid inequalities The preprocessing procedure described in the previous section allows us to eliminate a huge number of decision variables. Indeed, this procedure is able to remove about 78% of the decision variables (see Section 5.1). Another way in order to enhance the performance of the resolution of a MIP is to add cuts that will reduce the search space. In this section, we propose valid inequalities based on two different approaches. The first one is based on the lower bound on the completion time of a job j if it is processed at position k (λjk) described above. The second type of valid inequalities is based on the well-known improved preemptive lower bound of Della Croce and T'Kindt [14]. 4.1. Position-based valid inequalities The first type of valid inequalities is based on the fact that the completion time of the job at position k must be greater than a
326
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
○ Case (a): In an optimal solution, position ℓ þ 1 is occupied by a job scheduled after j in the SRPT sequence. Thus, improvements on the completion times of the jobs processed at positions ℓ and ℓ þ 1 are given by γ 1ℓ ¼ Ψ C ℓ ðSÞ and γ 2ℓ ¼ maxð0; minh A SuccðjÞ fmaxðr h ; Ψ Þ þ ph g C ℓ þ 1 ðSÞÞ, respectively. ○ Case (b): in an optimal solution, position ℓ þ 1 is occupied by a job scheduled before j in the SRPT sequence. Thus, improvements on the completion times of the jobs processed at position ℓ 1, ℓ and ℓ þ 1 are given by δ1ℓ ¼ maxð0; X C ℓ 1 ðSÞÞ; δ2ℓ ¼ maxð0; Y C ℓ ðSÞÞ and δ3ℓ ¼ max ð0; Z C ℓ þ 1 ðSÞÞ, respectively.
lower bound. Thus, we can use the estimation calculated on the preprocessing phase, namely λjk. The following proposition introduces the position-based valid inequalities. Proposition 1. n ςk Z ∑ λjk xjk ; 8 k ¼ 1; …; n
ð8Þ
j¼1
is a valid constraint for the model (1)–(7) Proof. Recall that ςk is the completion time of the job processed at position k. From the definition of λjk, ∑nj¼ 1 λjk xjk represents a lower bound on the completion time of the job processed at position k. □ It is note worthy to indicate that Haouari and Kharbeche [24] proposed similar valid inequalities for generic two machine flow shop problems denoted by F2 J ∑f j ðC j Þ, where f j ðÞ is a regular performance criterion. 4.2. SRPT-based valid inequalities The second type of valid inequalities is based on the SRPT properties exploited by Della Croce and T'Kindt [14] with the aim to improve the SRPT-based lower bound. Before introducing this type of valid inequalities, we shall recall the SRPT properties involved in the improved preemptive bound. In the sequel, we begin by developing the necessary properties. Then, we introduce the valid constraints. 4.2.1. SRPT properties For the sake of clarity, we adopt a slightly modified notation of Della Croce and T'Kindt [14]. Assume that job j A J is completed at ℓth position in the SRPT sequence, we denote by the following:
Pred(j) and Succ(j) the subset of predecessors and successors of job j in the SRPT sequence, respectively.
Ψ ¼ maxðC ℓ ðSÞ; C ℓ 1 ðSÞ þ pj Þ a lower bound with respect to C ℓ ðSÞ.
X ¼ C ℓ ðSj Þ ph ; Y ¼ maxfrj ; X; C ℓ 1 ðSÞg þ pj and Z ¼ maxfY; C ℓ
ðSÞg þ ph (where h is a job such that hA PredðjÞ and ph 4 pj ) a valid lower bound with respect to the completion times of the jobs at ℓ 1, ℓ and ℓ þ 1 in an optimal solution, respectively.
The improved SRPT bound [14] is based on determining a lower bound on the completion time increase of the job processed at a given position with respect to its SRPT completion time. The main idea of the procedure is based on the fact that, in an optimal solution, a job is scheduled either before, in the same position, or after its SRPT position. Thus, the following three cases may be distinguished:
Case 1: In an optimal solution, job j occupies a position k after
its SRPT position, i.e., k Zℓ þ 1. Then, α1ℓ ¼ C ℓ ðSj Þ C ℓ ðSÞ is the corresponding increase on the completion time of position ℓ. Furthermore, α2ℓ ¼ minfC ℓ ðSj Þ þ pj ; C ℓ þ 1 ðSj Þg C ℓ þ 1 ðSÞ is the increase in position ℓ þ 1. Case 2: In an optimal solution, job j occupies a position k before its SRPT position, i.e., k r ℓ 1. Then, β1ℓ ¼ maxf0; r j þ pj C ℓ þ 1 ðSÞ; C ℓ ðSÞ maxh A SuccðjÞ ph C ℓ 1 ðSÞg and β2ℓ ¼ max ð0; minh a j fmaxfr h ; C ℓ 1 ðSÞ; r j þpj ; C ℓ ðSÞ max i A PredðjÞ pi g þph g C ℓ ðSÞÞ are valid lower bounds on the increase of the completion times of the jobs processed at positions ℓ 1 and ℓ, respectively. Case 3: In an optimal solution, job j occupies a position k that corresponds to its SRPT position, i.e., k ¼ ℓ. Then, two subcases should be distinguished:
In order to simplify, hereafter, the presentation of the valid constraints, we introduce the fundamental Theorem of Della Croce and T'Kindt [14] in a slightly modified manner. Theorem 1 (Della Croce and T'Kindt [14]). Consider a job j completing in position 2 rℓ rn 1 in the SRPT sequence S. Let k be the position of job j in an optimal solution. We have the following results:
β1ℓ and δ1ℓ are valid lower bounds on the completion time increase
of the job processed at position ℓ 1 with respect to C ℓ 1 ðSÞ if k o ℓ and k ¼ ℓ (if the job scheduled at position ℓ þ 1 belongs to Pred(j)), respectively. β2ℓ ; minðγ 1ℓ ; δ2ℓ Þ and α1ℓ are valid lower bounds on the completion time increase of the job processed at position ℓ with respect to C ℓ ðSÞ if k o ℓ; k ¼ ℓ and k 4 ℓ, respectively. minðγ 2ℓ ; δ3ℓ Þ and α2ℓ are valid lower bounds on the completion time increase of the job processed at position ℓ þ 1 with respect to C ℓ þ 1 ðSÞ if k ¼ ℓ and k 4 ℓ, respectively.
Remark 1. As stated in Theorem 1, the corresponding increase in the completion times is valid for 2 rℓ r n 1. Therefore, Della Croce and T'Kindt [14] generalized the obtained results by taking into account the special cases where ℓ ¼1 and ℓ ¼ n. The authors reported that when ℓ ¼ 1 (respectively ℓ ¼ n) Cases 2 and 3b (respectively Cases 1 and 3b) do not hold. Moreover, when ℓ ¼ n, the Case 3a reduces to compute only γ1n. 4.2.2. Valid inequalities After having established the necessary background, we can now present the valid inequalities. The following proposition introduces valid inequalities related to position ℓ, ℓ þ1 and ℓ 1. Proposition 2. The following inequalities ςℓ Z C ℓ ðSÞ þ α1ℓ
ℓ1
n
∑
k ¼ ℓþ1
xjk þ β2ℓ ∑ xjk þ minðγ 1ℓ ; δ2ℓ Þxjℓ ; k¼1
8 ℓ ¼ 1; …; n ð9Þ
ςℓ þ 1 ZC ℓ þ 1 ðSÞ þ α2ℓ
n
∑
k ¼ ℓþ1
xjk þ minðγ 2ℓ ; δ3ℓ Þxjℓ ;
8 ℓ ¼ 1; …; n 1 ð10Þ
ℓ1
ςℓ 1 Z C ℓ 1 ðSÞ þ β1ℓ ∑ xjk þ δ1ℓ ðxjℓ þ k¼1
∑
h A PredðjÞ
xh;ℓ þ 1 1Þ;
8 ℓ ¼ 2; …; n
ð11Þ are valid for the model (MIP0) given by (1)–(7). Proof. Let k be the position where job j completes its execution in an optimal solution. From Theorem 1, we know that α1ℓ ; β2ℓ and minðγ 1ℓ ; δ2ℓ Þ are valid lower bounds on the completion time increase of the job in position ℓ with respect to C ℓ ðSÞ if k o ℓ; k ¼ ℓ and k 4 ℓ, respectively. Thus, the RHS of (9) constitutes a valid lower bound on the completion time of the job processed at position ℓ. Applying the same argument, constraints (10) and (11) hold. □
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
Remark 2. Taking into account the cases described in Remark 1, we shall initialize certain quantities as follows. On one hand, if ℓ ¼ 1, then β11 ¼ β21 ¼ β31 ¼ 0 and δ11 ¼ δ21 ¼ δ31 ¼ 1 . On the other hand, if ℓ ¼ n, then α1n ¼ α2n ¼ γ 2n ¼ 0 and δ1n ¼ δ2n ¼ δ3n ¼ 1: Thus, we must consider these special cases in the implementation of constraints (9)–(11). Remark 3. In inequality (11), ðxjℓ þ ∑h A PredðjÞ xh;ℓ þ 1 1Þ could be less than 0, it can be interesting to use inequality (11) with the following inequality: ℓ1
ςℓ 1 Z C ℓ 1 ðSÞ þ β1ℓ ∑ xjk ; k¼1
8 ℓ ¼ 2; …; n
ð12Þ
5. Computational results In this section, we report the results of a computational study that aims at assessing the performance of the proposed approach. To that aim, we carried out three sets of experiments. In the first set, we analyze the impact of the preprocessing on the original formulation. In the second set, we investigate the best combination of valid constraints that offers a good trade-off between tightness and efficiency. Finally, in the third set, we compare our best method with state-of-the-art exact methods. The tests were conducted on a set of generated instances following the scheme described in [25]. The number of jobs n was taken to be equal to 20, 40, 60, 80, 100, 120, 140, 160, 180 and 200 jobs. The processing times were randomly generated using the discrete uniform distribution on ½1; 100. The release dates were uniformly drawn from ½0; 50:5nR, where R controls the range of the distribution and belongs to the set f0:2; 0:4; 0:6; 0:8; 1; 1:25; 1:5; 1:75; 2:0; 3:0g. For each combination of n and R, five instances were generated. Hence, 50 instances are generated for each value of n. All the experiments were carried out on an Intel Xeon E5620 2.4 GHz processor, with 8 GB RAM memory, running Windows 7. The tested algorithms were coded in C and compiled under Microsoft Visual Studio 2010 using CPLEX 12.2. For all the experiments, the time limit for executing the procedures was set to 3600 s.
5.1. Impact of the preprocessing In this first experiment, we analyzed the impact of the preprocessing on the formulation. To this end, we tested the following two versions:
F0: The model (MIP0). F1: The model (MIP0) with the preprocessing procedure (see Section 3).
The results are displayed in Table 1. For each version, we provide the average required CPU time in seconds. The numbers between parenthesis are the number of unsolved instances. We note that we stopped the resolution of F0 to 60-jobs instances since it is clear that this version fails to solve larger instances. We see from Table 1 that the original formulation F0 cannot solve instances with up to 40 jobs. We remark that formulation F1 greatly improves the results since it is able to solve 78% of the 80-jobs instances. This is due to the fact that the preprocessing is able to eliminate a large number of decision variables xjk. In fact, the tests show that this procedure is capable of removing 78.24% of the decision variables on average.
327
Table 1 Impact of the preprocessing. n
F0
F1
20 40 60 80 100
4.53 2578.38 (33) 3276.47 (44) – –
0.04 0.88 189.82 (1) 946.49 (11) 2047.42 (26)
Table 2 Impact of the valid constraints. n
F1
F2
F3
F4
20 40 60 80 100
0.04 0.88 189.82 (1) 946.49 (11) 2047.42 (26)
0.02 0.24 1.16 4.49 81.91
0.05 0.39 1.18 4.44 54.12
0.02 0.26 0.76 2.69 17.62
5.2. Impact of the valid constraints In this second phase of experiments, we investigated the impact of the constraints on the formulation. We tested the following versions:
F2: Formulation F1 along with position-based constraints (8). F3: Formulation F1 along with SRPT-based constraints (9)–(11). F4: Formulation F1 after appending all valid constraints. We first reported in Table 2 the performance of F1, F2, F3 and F4 and then observed the following points:
While F1 fails to solve 52% of 100-jobs instances, the three
versions F2, F3 and F4 were able to solve all the instances from 20 to 100 jobs within the time limit of 1 h. The best performance is given by appending all the new cuts. Actually, F4 is less time consuming than F2 and F3.
We stopped the simulation at n ¼100 as we are confident that the proposed constraints improve the performance of the resolution. Therefore, we assumed in our implementation that all the new cuts are appended to MIP0. Thus, F4 will constitute our reference approach in order to compare it with methods of the literature. 5.3. Comparison with state-of-the-art exact methods In this third set of experiments, we compare the performance of F4 with respect to state-of-the-art exact methods. It is worth noting that the best known exact methods for this problem are those of Jouglet et al. [20] and T'Kindt et al. [34]. The two methods are based on a branch-and-bound strategy. In order to compare our formulation with these methods, we developed a branch-andbound algorithm based on the strategies and dominance rules found in the latter papers. More precisely, the branch-and-bound algorithm that we implemented is based on these concepts:
Best active node search strategy: T'Kindt et al. [34] compared
different search strategies for this problem and found that this strategy is the best one. SRPT rule for the computation of lower bounds on nodes. Dynamic programming dominance rule as described in [34]. Dominance rules relying on scheduled and unscheduled jobs introduced by Jouglet et al. [20]. It is worth noting that these rules increase the average resolution time but in return they
328
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
are able to prune more nodes. Thus, the algorithm is able to solve larger instances. In the following, this branch-and-bound procedure is denoted by BB. Moreover, in order to get a better idea between the developed method in this paper and the branch-and-bound algorithm, we tested a version of F4 with only one thread hereafter denoted by F41Th . Indeed, CPLEX solves the formulation using as many threads as the number of available cores in the testing machine, i.e., in our case, eight threads were launched. We added this version in order to have a “fair” comparison between the mathematical formulation and the branch-and-bound method. For each procedure, we provide
Time: The average CPU time in seconds. #Unsolved: The number of unsolved instances. #Best: The number of times the procedure achieves the best value among all other procedures.
F4 exhibits an excellent performance both in terms of efficiency
n
ζ k Z ςk ∑ dj xjk ; j¼1
In this phase, we focused our attention to larger instances from 120 to 200 jobs. The results are presented in Table 3. From Table 3, we made the following observations:
scheduling problems. For the sake of clarity, we discuss this foregoing question on two single machine scheduling problems subject to release dates where the objective is to minimize the total tardiness and the total weighted tardiness. Before proceeding further, let us denote by wj, dj and T j ¼ maxð0; C j dj Þ the weight, the due date and the tardiness of job jA J, respectively. Following these notations, the former problems are denoted by 1jr j j∑T j and 1jr j j∑wj T j and will be referred to, in the sequel, as the P problems. The assignment-based model (1)–(7) cannot be used as it is in the P problems since a valid instance of these problems requires some decision variables and constraints to be added. For 1jr j j∑T j , we introduce the semi-continuous decision variables ζk that represent the tardiness of the job processed at position k; 8 k A f1; …; ng. The objective function (1) is replaced by min ∑nk ¼ 1 ζ k . In addition, the following set of constraints must be added:
and effectiveness. Indeed, it consistently exhibits the best average time of resolution and it solves more instances. Interestingly, we remark that the version with one thread F41Th is competitive with the branch-and-bound procedure and even manages to solve larger instances. Surprisingly, for n ¼200, F41Th yields the best value more times than F4 and BB. This observation can be explained by “the variability performance” (see [27]). Actually, on average, using several threads may lead to enhancing the solver performance. However, on a single instance, the performance may arbitrarily depend on the number of threads.
Pushing our analysis a step further, we performed a detailed analysis of the results according to the parameter R. Table 4 presents the mean CPU time on all instances (from 120 to 200 jobs) for a fixed value of R. The total number of unsolved instances is shown between parenthesis. We remark from Table 4 that the hard instances come for the range of values 0.4–1.25. For these instances, the two versions of our method F4 and F41Th dominate the branch-and-bound procedure BB. For the other easier instances, we noticed that BB is competitive with the proposed methods. Specially for the instances with the value 1.75 and higher, the branch-and-bound is more efficient than our methods.
6. Extension for general scheduling problems In this section, we discuss how the different components of our approach should be adapted to be used for other general
8 k ¼ 1; …; n
On the contrary, to represent the tardiness of a job in the 1jr j j∑wj T j problem, one needs to introduce new decision variables. Indeed, for the 1jr j j∑T j problem, we can consider the tardiness of the job processed at position k instead of the tardiness of the job itself. However, this is not possible for the 1jr j j∑wj T j problem since both the job's weight and the tardiness are considered. Let us denote by θj the tardiness of the job j; 8 j A f1; …; ng. The objective function (1) is consequently replaced by min ∑nj¼ 1 wj θj and the following constraints are added: θj Z ςk dj Mð1 xjk Þ;
8 j; k ¼ 1; …; n
where M is an upper bound on how large ςk can be. Unfortunately, we face a big-M model. It is worth noting that this type of models is known to be the hardest to solve (see [10]). In order to remove the big-M dependency, a possible solution consists in proposing alternative inequalities. We refer the reader to Keha et al. [26] who proposed alternative inequalities for a set of single machine problems without release dates. Another solution could be the use of the automatic problem reformulation [10]. These solutions, however, need both a theoretical and a computational thorough investigation. The lower bound on the completion time of job j if it is processed at position k (λjk) is one of the most important parts of our algorithm. Indeed, this lower bound allows us to derive the preprocessing procedure and valid inequalities (8). Recall that, the lower bounds λjk is a valid lower bound for any sequence. This result is a straightforward consequence of Lemmas 1 and 2. Thus, inequalities (8) remain valid for the P problems. Concerning inequalities (9)–(11), recall also that Cl(S) is a valid lower bound on the completion time of the job processed at position k (see Lemma 1). Therefore, from Theorem 1 and Proposition 2, we can conclude that these inequalities are also valid for the P problems. In order to adapt the preprocessing procedure, we must propose a lower bound on the total (weighted) tardiness while
Table 3 Comparison of methods. n
120 140 160 180 200
F4
BB
F41Th
Time (s)
#Unsolved
#Best
Time (s)
#Unsolved
#Best
Time (s)
#Unsolved
#Best
51.09 38.71 309.89 437.02 712.00
0 0 3 4 7
50 50 49 49 48
76.26 89.53 520.40 893.90 1310.22
0 0 4 8 16
50 50 47 42 36
163.60 116.59 454.24 572.12 876.65
1 0 4 5 9
50 50 48 45 49
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
329
Table 4 Detailed results. R
0.2
0.4
0.6
0.8
1
1.25
1.5
1.75
2
3
F4 BB F41Th
52.82 60.33 85.46
478.15 (2) 569.35 (3) 769.22 (2)
1039.09 (5) 1289.96 (5) 1383.48 (7)
1160.03 (6) 2023.95 (12) 1529.46 (8)
255.93 (1) 1559.13 (7) 421.17 (1)
87.36 255.91 (1) 156.30 (1)
8.24 15.22 6.87
6.44 3.99 5.66
5.28 1.74 4.83
4.08 1.04 3.92
taking into account that a given job j is scheduled at position k. To this end, we use the lower bound proposed by Haouari and Kharbeche [24] that consists in solving the linear assignment problem. In our problem, the assignment costs are κ j0 k0 ¼ 0 wj0 maxðλj0 k0 dj0 ; 0Þ; 8 j0 ; k A f1; …; ng (it is obvious that the weights wj are equal to 1 for 1jr j j∑T j ). In order to ensure that the solution 0 will take job j at position k, we set to infinity the costs κ jk0 ; 8 k a k 0 0 and the costs κ j k ; 8 j a j. As described in Section 3, if this value is greater than an upper bound on the total (weighted) tardiness, then a decision variable can be eliminated from the model. As the linear assignment problem is solved in Oðn3 Þtime and will be repeated for each job j A f1; …; ng and each position k A f1; …; ng, the total complexity of the procedure is Oðn5 Þ. This leads to investigate more appropriate lower bounds on the (weighted) tardiness of job j scheduled at position k. Finally, we note that the approach developed in this paper can be extended to permutation flowshop problems where the objective function is to minimize total completion time or total tardiness problem. For these problems we think that the assignmentbased formulation could be very useful. The weighted version of these problems needs a deeper analysis to overcome the big-M problem. Moreover many components of our approach could be used for other general scheduling problems. However, to achieve a good performance, we believe that it is very important to exploit the characteristics of the problem when designing a method. For example, for the two-machine total completion time flowshop problem, we noticed that Haouari and Kharbeche [24] used the Johnson's rule [19] to calculate the λjk values and derived a valid inequality to strengthen the continuous relaxation of a mathematical model proposed by Hoogeveen et al. [18]. These values can be used in our preprocessing procedure to reduce the binary variables of the mathematical model. The valid inequality along with the preprocessing procedure should enhance the performance of the exact resolution of the mathematical model. In a nutshell, two important workarounds would prove to be useful to achieve good performances of our approach for general scheduling problems. The first is to use suitable lower bounds for each studied problem to be able to eliminate decision variables. The second issue consists in deriving more tailored valid inequalities.
7. Conclusion In this paper, we presented a new mathematical approach to solve the one machine scheduling problem subject to release dates. First, we presented a preprocessing procedure which aims to reduce the number of binary variables from the mathematical model. Then, we developed two types of valid inequalities. The first set of inequalities is deduced from the computation of lower bounds of the completion times of the jobs scheduled at these positions. The second set is based on improvements on completion times of the Shortest Remaining Processing Time sequence. Computational results showed that our simple method through the use of a commercial solver is competitive with state-of-the-art branch-and-bound algorithms and is able to solve large instances up to 140 jobs.
Acknowledgments Many thanks to Professor Vincent T'Kindt for valuable discussion about the branch-and-bound procedure. Also, the authors would like to thank the anonymous referees for their insightful comments that substantially improved the quality of the paper. References [1] Ahmadi RH, Bagchi U. Lower bounds for single-machine scheduling problems. Naval Res Logist 1990;37:967–79. [2] Baker KR. Introduction to sequencing and scheduling. New York: Wiley; 1974. [3] Bianco L, Ricciardelli S. Scheduling of a single machine to minimize total weighted completion time subject to release dates. Naval Res Logist Q 1982;29:151–67. [4] Chand S, Traub R, Uzsoy R. Single-machine scheduling with dynamic arrivals: decomposition results and an improved algorithm. Naval Res Logist 1996;43:709–19. [5] Chand S, Traub R, Uzsoy R. An iterative heuristic for the single-machine dynamic total completion time scheduling problem. Comput Oper Res 1996;23:641–51. [6] Chand S, Traub R, Uzsoy R. Rolling horizon procedures for the single machine deterministic total completion time scheduling problem with release dates. Ann Oper Res 1997;70:115–25. [7] Chandra R. On n=1=F dynamic deterministic systems. Naval Res Logist Q 1979;26:537–44. [8] Chu C. A branch-and-bound algorithm to minimize total tardiness with different release dates. Naval Res Logist 1992;39:265–83. [9] Chu C. Efficient heuristics to minimize total flow time with release dates. Oper Res Lett 1992;12:321–30. [10] Codato G, Fischetti M. Combinatorial benders' cuts for mixed-integer linear programming. Oper Res 2006;54:756–66. [11] Conway RW, Maxwell WL, Miller LW. Theory of scheduling. MA: AddisonWesley Publishing Company; 1967. [12] Dessouky MI, Deogun JS. Sequencing jobs with unequal ready times to minimize mean flow time. SIAM J Comput 1981;10:192–202. [13] Deogun JS. On scheduling with ready times to minimize mean flow time. Comput J 1983;26:320–8. [14] Della Croce F, T'Kindt V. Improving the preemptive bound for the onemachine dynamic total completion time scheduling problem. Oper Res Lett 2003;31:142–8. [15] Della Croce F, T'Kindt V. A recovering beam search algorithm for the onemachine dynamic total completion time scheduling problem. J Oper Res Soc 2002;53:1275–80. [16] Dyer ME, Wolsey LA. Formulating the single machine sequencing problem with release dates as a mixed integer program. Discr Appl Math 1990;26: 255–70. [17] Gharbi A, Ladhari T, Msakni MK, Serairi M. The two-machine flowshop scheduling problem with sequence-independent setup times: new lower bounding strategies. Eur J Oper Res 2013;231:69–78. [18] Hoogeveen H, Van Norden L, Van de Velde S. Lower bounds for minimizing total completion time in a two-machine flow shop. J Sched 2006;9:559–68. [19] Johnson S. Optimal two- and three-stage production schedules with set up times included. Naval Res Logist Q 1954;1:61–8. [20] Jouglet A, Baptiste P, Carlier J. Branch-and-bound algorithms for total weighted tardiness. In: Leung JYT, editor. Handbook of scheduling: algorithms, models and performance analysis; 2004. [21] Jouglet A, Savourey D, Carlier J, Baptiste P. Dominance-based heuristics for one-machine total cost scheduling problems. Eur J Oper Res 2008;184: 879–99. [22] Graham RL, Lawler EL, Lenstra JK, Rinnooy Kan AHG. Optimization and approximation in deterministic sequencing and scheduling: a survey. Ann Discr Math 1979;5:287–326. [23] Guo Y, Lim A, Rodriguez B, Yu S. Minimizing total flow time in single machine environment with release time: an experimental analysis. Comput Ind Eng 2004;47:123–40. [24] Haouari M, Kharbeche M. An assignment-based lower bound for a class of two-machine flow shop problems. Comput Oper Res 2013;40:1693–9. [25] Hariri AMA, Potts CN. An algorithm for single machine sequencing with release dates to minimize total weighted completion time. Discr Appl Math 1983;5:99–109.
330
A. Kooli, M. Serairi / Computers & Operations Research 51 (2014) 323–330
[26] Keha AB, Khowala K, Fowler JW. Mixed integer programming formulations for single machine scheduling problems. Comput Ind Eng 2009;56:357–67. [27] Koch T, Achterberg T, Andersen E, Bastert O, Berthold T, Bixby RE, et al. MIPLIB 2010. Math Program Comput 2011;3:103–63. [28] Lasserre JB, Queyranne M. Generic scheduling polyhedra and a new mixedinteger formulation for single-machine scheduling. In: Proceedings of the second IPCO conference. Pittsburgh: Carnegie-Mellon University; 1992. p. 136–49. [29] Lawler EL, Lenstra JK, Rinnooy Kan AHG, Shmoys DB. Sequencing and scheduling: algorithms and complexity. In: Graves SC, Rinnooy Kan AHG, Zipkin P, editors. Handbooks in operations research and management science, Logistics of production and inventory, vol. 4. Amsterdam: North-Holland; 1993.
[30] Liu J, MacCarthy BL. Effective heuristics for the single machine sequencing problem with ready times. Int J Prod Res 1991;29:1521–33. [31] Rakrouki MA, Ladhari T, T'Kindt V. Coupling genetic local search and recovering beam search algorithms for minimizing the total completion time in the single machine scheduling problem subject to release dates. Comput Oper Res 2012;39:1257–64. [32] Smith WE. Various optimizers for single-state production. Naval Res Logist Q 1956;3:59–66. [33] Sousa JP, Wolsey LA. A time indexed formulation of non-preemptive single machine scheduling problems. Math Program 1992;54:353–67. [34] T'Kindt V, Della Croce F, Esswein C. Revisiting branch and bound search strategies for machine scheduling problems. J Sched 2004;7:429–40.