# Single-Machine Scheduling with Fixed Periodic Preventive Maintenance to Minimise the Total Weighted Completion Times

Authors: Hongming Zhou; Ya-Chih Tsai; Shenquan Huang; Yarong Chen; Fuh-Der Chou
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The single-machine scheduling problem with fixed periodic preventive maintenance, in which preventive maintenance is implemented periodically to maintain good machine operational status and decrease the cost caused by sudden machine failure, is studied in this paper. The adopted objective function is to minimise the total weighted completion time, which is representative of the minimisation of the global holding/inventory cost in the system. This problem is proven to be NP-hard; a position-based mixed integer programming model and an efficient heuristic algorithm with local improvement strategy are developed for the total weighted completion time problem. To evaluate the performances of the proposed heuristic algorithms, two new lower bounds are further developed. Computational experiments show that the proposed heuristic can rapidly achieve optimal results for small-sized problems and obtain near-optimal solutions with tight average relative percentage deviation for large-sized problems.


---

## Body

## 1. Introduction

Production scheduling is concerned with the allocation of limited resources (machines or operators) to tasks in the manufacturing systems of enterprises to optimise certain objective functions such that their competitive position is maintained or improved in fast-changing markets. In traditional scheduling theory, it is usually assumed that machines can operate continuously during a scheduling planning horizon. However, uncertain machine breakdown or preventive maintenance activities cannot be ignored in the real world, as they cause the machines to be unavailable and production continuity to be slowed or stopped. As a result, the overall performance of the system will be decreased. Therefore, the constraint of machine unavailability taken into account during production scheduling is very important, and it makes scheduling problems more complex and realistic.Preventive maintenance (PM) causes machine unavailability, as mentioned above, but it has gradually become a common practice in many companies as a priori maintenance measure to prevent potential malfunctions and critical nonavailability of a system. For this reason, many researchers and practitioners have gradually regarded production scheduling jointly with PMas a common and significant issue. In the literature, most scheduling models assume that the intervals between maintenance activities are known in advance and that maintenance time (unavailability time) is a constant. This assumption is reasonable in the case that PM is determined by engineering staff in advance. To date, such models have received considerable attention from researchers, and comprehensive reviews of this type of research have also been provided [1–4].Regarding the scheduling problem with PM activities, the single-machine problem with a single interval of machine nonavailability is a basic model, denoted as1,h1‖Obj, where h1 indicates the PM activity occurs once in the planning horizon, and Obj is the objective function adopted in scheduling problems. Adiri et al. [5] and Lee and Liman [6] proved the single-machine problem with the total completion time denoted as1,h1‖∑Ci to be NP-hard. They showed that the SPT rule has a tight worst-case error bound of 2/7. For the same problem, Sadifi et al. [7] showed their proposed approximation algorithm MSPT (modified SPT) had a worst-case error bound of3/17. Breit [8] provided a parametric O (nlogn)-algorithm H and the best error bound of this algorithm was 7.4%. Another related flow-time objective function is the weighted completion time, which also has been discussed in the literature, because weights can measure the importance of processing jobs, or quantify the stocking cost per a unit of time of jobs in a system. Kacem et al. [9] studied the problem of 1,h1‖∑wiCi. They proposed a branch-and-bound (BAB) method, a mixed integer programming model (MIP), and a dynamic programming (DP) method. For the same problem, Kacem and Paschos [10] modified the WSPT rule slightly and provided a polynomial 3−5/2-differential algorithm.All the research cited above addressed cases with a single PM activity in the planning horizon. Without loss of generality, it is necessary to consider multiple PM activities, as maintenance is scheduled periodically in many manufacturing systems. That is, there is usually more than one maintenance period in the planning horizon. Liao and Chen [11] extended a single PM activity to multiple and periodic maintenance activities for a single-machine scheduling problem; they proposed a BAB algorithm and an efficient heuristic for minimising the maximum tardiness. For the objective of minimising the number of tardy jobs, Chen [12] proposed an efficient heuristic based on Moore’s algorithm [13] and the BAB method to solve the problem. Lee and Kim [14] proposed a two-phase heuristic algorithm for the same problem to minimise the number of tardy jobs. Batun and Azizoglu [15] proposed a BAB algorithm to solve the same problem, in terms of total flow time. To minimise the makespan, Perez-Gonzalez and Framinan [16] developed an MIP model and two-phase heuristic algorithm based on bin-packing dispatching rules to obtain optimal or near-optimal solutions. Zhou et al. [17] also proposed different heuristic algorithms with a local improvement mechanism for the same problem.In the related research, a variety of production scheduling problems exist that consider different maintenance types in different manufacturing shop floors. Relative to fixed PM activity with known starting time and duration, flexible maintenance is another attractive research topic where the maintenance start time is not fixed in advance. In some cases of steel processing or semiconductor manufacturing, the maintenance start time depends on a specific threshold value such as maximum continuous working time or maximum dirt. This kind of assumption can be found in Qi et al. [18], Mosheiov and Sarig [19], Cui and Lu [20], Su and Wang [21], Pang et al. [22], Huang et al. [23], Chung et al. [24], etc. Another deterministic case is when the maintenance period [u, v] is arranged in advance and maintenance operations should be undertaken within the period [u, v], where the maintenance time (w) is not longer than the maintenance period, i.e., w≤v−u. Some contributions under this assumption for different scheduling problems can be found in Chen [25], Yang et al. [4], etc. All studies mentioned above assumed that the maintenance time was a constant. More recently, some researchers have taken variable maintenance time into account, because in some realistic applications the maintenance time either increases or decreases depending on the machine conditions or the maintenance skill improvement. Bock et al. [26] considered the case of job-dependent machine deterioration, where the maintenance time is dependent on the amount by which the maintenance level is increased; that is, more improvement of the machine’s state requires more maintenance time. Several recent studies have been conducted to address variable maintenance activity on scheduling problems [27, 28].This study focuses on scheduling jobs on a single machine with periodic unavailability periods due to implementing PM activities; the objective is to minimise the total weighted completion times. In our study, the unavailability periods are fixed and known in advance. In the literature, a recent study by Krim et al. [29] considered the same scheduling problem. Our contribution lies in developing a different MIP model, two lower bounds, and a new heuristic algorithm while examining the effectiveness of the proposed methods in terms of various problems. In the following section, we propose an MIP model and some lower bounds, which are compared with the results by Krim et al. [29]. The results show that our model finds the optimal solution in all problems of the model by Krim et al. [29] at a lower CPU time, and some lower bounds are addressed in Section 3. In Section 4, we propose a heuristic algorithm based on the analytical properties. The proposed heuristic algorithm is examined on large instances, and the results are compared and reported in Section 5. Finally, Section 6 states our conclusions along with future research directions.


## 2. Mixed Integer Programming Model

In this section, we first describe the considered problem. There aren independent nonresumable jobs to be processed on a single machine. PM activities are carried out periodically, so the machine is not constantly available in the planning horizon, which is composed of an availability period of T−δ time units followed by a nonavailability period of δ time units. A nonresumable job refers to a job that cannot be finished before a maintenance activity and thus has to restart. Each availability period of T−δ is frozen, and machine idle may occur because of job nonresumablity, as shown in Figure 1. Using the three-field notation of Graham et al. [30], we denote this scheduling problem as 1,hknr∑wiCi, which has been proven to be NP-hard [29]. For the 1,hknr∑wiCi problem, some assumptions and notations are given as follows:(i)
Assumptions:(1)
All jobs and machine are available for timet = 0(2)
The processing times and weights of jobs are deterministic and known in advance(3)
The machine can process only one job at a time without pre-emption(4)
For a machine, the maintenance timeδ is deterministic and known in advance; after maintenance, the status of the machine becomes the initial status(5)
The fixed intervalT is known in advance, and it is not large enough to process all the jobs (i.e., ∑pi>T−δ); however, the availability period of T−δ should be greater than or equal to the maximum processing time of jobs, i.e., T−δ≥max1≤i≤npi(ii)
Notation:Figure 1
A feasible schedule for a single problem with fixed periodic maintenance activities.n: number of jobspi: processing time of jobs i, i=1,2,…,nwi: weight of jobs i, i=1,2,…,nδ: maintenance timeT−δ: an available period, where T is a given fixed time periodM: a large integer numberXik: a binary decision variable used to judge whether job i is processed at position k, defined as follows:Xik=1if  jobiis scheduled at positionk0otherwise where i=1,2,…,n; k=1,2,…,nYk: a binary decision variable that is used to judge whether position k is assigned to execute maintenance, defined as follows:Yk=1if  maintenance is implemented at positionk0otherwisewherek=1,2,…,nQk: sum of processing time at position k, k=1,2,…,nIk: idle time after position k, k=1,2,…,nCi: completion time of job i, i=1,2,…,nSTk: start time at position k,k=1,2,…,nPTk: processing time at position k, k=1,2,…,nCTk: completion time at position k,k=1,2,…,n.Next, an MIP formulation is provided for the problem, and the model is constructed based on the viewpoint of each position. This is different from another MIP model that is built based on the relationship of jobs [29]. The proposed MIP model can be described as follows:(1)minimise∑i=1nwiCi,subject to(2)∑k=1nXik=1,∀i=1,2,…,n,(3)∑i=1nXik=1,∀k=1,2,…,n,(4)PTk=∑i=1npi×Xik,∀k=1,2,…,n,(5)STk≥CTk−1+Ik−1+δ×Yk−1,∀k=2,…,n,(6)CTk≥STk+PTk,∀k=1,2,…,n,(7)Q1=∑i=1npi×Xi1,(8)Qk−1+∑i=1npi×Xik≤Qk+M×Yk−1,∀k=2,…,n,(9)∑i=1npi×Xik≤Qk+M×1−Yk−1,∀k=2,…,n,(10)Qk+M×Yk−1−1≤∑i=1npi×Xik,∀k=2,…,n,(11)Qk≤T−δ,∀k=1,2,…,n,(12)M×Yk−Ik≥0,∀k=1,2,…,n,(13)M×1−Yk−T−δ−Qk−Ik≥0,∀k=1,2,…,n,(14)Cj+M×1−Xik≥CTk,∀j,k=1,2,…,n,(15)ST1=0,(16)Yn=0,(17)In=0,Constraint (1) is the objective function that minimises the total weighted completion time of all jobs. Constraint (2) and constraint (3) ensure that each job can only be placed at one position and that each position can only be occupied by one job. Constraint (4) specifies the processing time of job that is assigned at position k. Constraint (5) ensures that the time for position k to assign a job should start no earlier than the sum of its predecessor’s finish time, idle time, and possible maintenance time. Constraint (6) establishes a link between the start time and finish time for position k. Constraint (7) calculates the accumulated processing time of position 1. Constraints (8) to (10) calculate the accumulated processing time of position k, excluding position 1. Constraint (11) ensures that the accumulated processing time of position k in the sequence does not exceed the predefined available period. Constraints (12) and (13) are used to calculate the idle time immediately after position k. Constraint (14) is the completion time of each job. Constraint (15) specifies the starting time for position 1 when the sequence is equal to zero. Constraints (16) and (17) specify the maintenance time and idle time for the last position in the sequence is equal to zero, respectively.In the proposed model, the number of constraints isn2+11×n−3. The number of integer and binary variables is 6×n and n2+n, respectively, and n is the number of jobs. In the model proposed by Batun et al. [15] and mentioned above, the number of constraints is 2n2+3×n. The number of integer and binary variables is n and 2n2, respectively. For even relatively modest problem sizes, there are too many integer variables in the two models for the problem to be solved; however, the MIP model is still an appropriate way to describe the problem complexity and obtain optimal solutions for a small-size problem. Form a practical viewpoint, heuristic algorithms are required to solve large-sized problem.


## 3. Lower Bounds

A lower bound value is an estimate of the best possible solution and is useful to evaluate the performance of heuristic algorithms when the optimal solutions cannot be found. Krim et al. [29] proposed four different lower bound procedures by polynomial relaxed problems. According to their research, we develop two different lower bound values for the considered problem in this paper. Before introducing our lower bound values, the four different lower bound procedures of Krim et al. [29] are briefly described as follows:(i)
LB1=∑i=1nwi∑j=1ipj, which is based on the relaxed problem of 1//∑i=1nwici(ii)
LB2=p∑i=1niwi+T∑k=1m−1k∑i=1nk+1wn1+n2+⋯+nk+i, which is based on the special case that the processing times of all jobs are the same(iii)
LB3=∑i=1m1wi∑j=1ipj+∑i=gm2wiT+∑j=gipj, which is based on the relaxed problem of 1,h1//∑i=1nwici and applies a dynamic programming algorithm to solve the 1,h1//∑i=1nwici(iv)
LB4=∑i=1npiwi+T×mm−1/2, which is based on the special case that the weight values of all jobs are the sameFor more detail of these lower bounds, please see the work of Krim et al. [29]. Among the four lower bounds, we find the formula of LB2 to be questionable, and we use a counterexample of 8 jobs in which pj=p, wj=w, and T=3p+δ to illustrate the question.For the above counterexample, the optimal sequence is as shown in Figure2, and the optimal solution is 15pw+7Tw. However, according to the LB2 formula where min1≤i≤8pi=p, wi=w, m=∑pi/3p=8p/3p=3 and x denotes the least integer greater than or equal to x, n1=3, n2=3, n3=2; the calculation of LB2 is as follows:(18)LB2=f∗=∑i∈B1win1p+∑i∈B2win2p+T+∑i∈B3win3p+2T,=∑i∈B1w×3×p+∑i∈B2w3p+T+∑i∈B3w2p+2T,=22pw+7Tw≥15pw+7Twoptimal solutionit is conflict.Figure 2
Optimal sequence of the counterexample.Regarding our proposed lower bound values, first we supposed an optimal sequence, as shown in Figure3. Other notations used in Figure 3 are described below.(i)
m∗: number of batches for an optimal solution(ii)
Bk: kth batch, k=1,2,…,m∗(iii)
nk: number of jobs in the kth batch, k=1,2,…,m∗(iv)
∑k=1m∗nk=n(v)
Δk: idle time occurred in the kth batch, k=1,2,…,m∗−1(vi)
Ji: job in the ith position of an optimal sequence(vii)
pi: processing time of the job in the ith position of an optimal sequence(viii)
wi: weight value of the job in the ith position of an optimal sequenceFigure 3
Schematic diagram of an optimal solution.For the optimal sequence as shown in Figure3, we can obtain the total weighted completion time as follows:(19)TWC=TWCB1+TWCB2+⋯+TWCBm∗,TWCB1=∑i∈B1wiCi=p1w1+p1+p2w2+⋯+p1+p2+⋯+pn1wn1,TWCB2=∑i∈B2wiCi=T+pn1+1wn1+1+T+pn1+1+pn2+1wn1+2+⋯+T+pn1+1+⋯+pn1+n2wn1+n2=p1+p2+⋯+pn1+1+Δ1+δwn1+1+p1+p2+⋯+pn1+2+Δ1+δwn1+2+⋯+p1+p2+⋯+pn1+n2+Δ1+δwn1+n2,TWCB3=∑i∈B3wiCi=p1+p2+⋯+pn1+n2+1+Δ1+Δ2+2δwn1+n2+1+⋯+p1+p2+⋯+pn1+n2+n3+Δ1+Δ2+2δwn1+n2+n3,⋮TWCBm∗=∑i∈Bm∗wiCi=m∗−1T+pn1+n2+⋯+nm∗−1+1wn1+n2+⋯+nm∗−1+1+⋯+m∗−1T+pn1+n2+⋯+nm∗−1+1+⋯+pn1+n2+⋯+nm∗wn1+n2+⋯+nm∗=p1+p2+⋯+p∑k=1m∗−1nk+1+∑k=1m∗−1Δk+m∗−1δw∑k=1m∗−1nk+1+⋯+p1+p2+⋯+pn+∑k=1m∗−1Δk+m∗−1δwn,TWC=TWCB1+TWCB2+⋯+TWCBm∗=∑i=1nwi∑j=1ipi+∑k=2m∗Δk−1∑b=1nkw∑x=1k−1nx+b+δ∑k=1m∗k−1∑b=1nkw∑x=1k−1nx+b.According to Figure3 and the formula mentioned above, it is obvious that the number of batches m∗ is an important factor influencing the schedule results, and the obvious estimation for the number of batches is m′=∑pi/T−δ [29]. In this paper, we adopt another lower bound concept provided by Martello and Toth [31] for the bin-packing problem to estimate the number of batches. We letkmax=T−δ/2, where xdenotes the greatest integer less than or equal to x and 1≤a≤kmax. We put the jobs where their processing times satisfy the condition of pi≥T−δ−a into set A. The other jobs where their processing times satisfy the condition of pi≤a are put into set B. We let NA be the number of jobs in set A. Now, for each a, the waste time is calculated based on the formula of Wastea=max0,NAT−δ−∑i∈Api−∑i∈Bpi, and the maximum waste is obtained by Wastemax=Max1≤a≤kmaxWastea. Considering maximum waste, the estimation of the number of batches will be m″=∑pi+Wastemax/T−δ. For m′, m″, and m∗, the relation among them is m′≤m″≤m∗.For a job sequenceJ1′,J2′,…,Jn′ where all jobs are sorted according to WSPT rule, i.e., p1′/w1′≤p2′/w2′≤⋯≤pn′/wn′, we use the estimation m″ as a substitute for the best number of batches m∗ and ignore the impact of the total weighted idle time, i.e., ∑k=2m∗Δk−1∑b=1nkw∑x=1k−1nx+b on the objective function, i.e., ∑k=2m∗Δk−1∑b=1nkw∑x=1k−1nx+b=0. Then, the new lower bound (NewLB) is derived based on equation (18) as below, and the value of NewLB is less than or equal to the optimal solution TWC.NewLB=TWC′=∑i=1nwi′∑j=1ipi′+δ∑k=1m″k−1∑b=1nk′w∑x=1k−1nx′+b′≤TWC, where nk′ is the minimum number of jobs in Bk. The estimation of nk′ is as follows:(20)nk′=T−δpmax,1≤k≤m″−1,max1,∑pj−T−δm″−1pmax,k=m″,where pmax=maxpi. To illustrate the proposed lower bound values, an instance with 5 jobs is given; the information for the instance is shown in Table 1, and the optimal solution of 194 is obtained by the MIP model.Table 1
An instance with 5 jobs,T = 10, and δ=3.
J1J2J3J4J5Processing timespi55554Weightswi12341First,kmax=T−δ/2=10−3/2=3. We let a be equal to 1, and SetA=∅, NA=0, since none of jobs satisfied the condition of pj≥10−3−1=6, SetB=∅. Waste1=max0,NAT−δ−∑i∈Api−∑i∈Bpi=max0,0−0−0=0.We leta be equal to 2, SetA=J1,J2,J3,J4, NA=4 since these jobs satisfied the condition of pj≥10−3−2=5, SetB=∅.Waste2=max0,NAT−δ−∑i∈Api−∑i∈Bpi=max0,4×7−20−0=8. For a=3, SetA=J1,J2,J3,J4,J5, NA=5 since these jobs satisfied the condition of pj≥10−3−3=4. SetB=∅. Waste3=max0,5×7−24−0=11. The maximum waste time ismax0,8,11=11, and the estimation of the number of batches is m″=∑pj+Wastemax/T−δ=24+11/7=5. For each batch, we can obtain the minimum number of jobs according to equation (19), which is n1′=n2′=n3′=n4′=n5′=1. Next, the jobs are sorted based on the WSPT rule, and NewLB is obtained as follows:NewLB=∑i=1nwi′∑j=1ipi′+δ∑k=1m″k−1∑b=1nk′w∑x=1k−1nx′+b′=123+3×4×1+3×1+2×2+1×3=123+3×14=165RecallingLB4, the estimation is composed of ∑i=1npiwi and T×mm−1/2; now replacem with m″; we can derive LB5=∑i=1npiwi+T×m″m″−1/2. For this example, LB5=54+10×5×4/2=154. Table 2 lists the results obtained by the five lower bounds for this five-job example where the performance of lower bound is evaluated by the related percentage gap (RPG) from the optimal solution.Table 2
Comparison of the five lower bounds according to the 5-job example.
LB1LB3LB4LB5NewLBLower bound values123158114154165RPG (%)36.59818.55741.23720.61914.948


## 4. Heuristic Algorithms

Despite the relative success of the two MIP models mentioned above for finding optimal solution reporting in computational experiments, the models are still incapable of solving medium and large instances. It is necessary to develop efficient heuristic algorithms. Many existing heuristic algorithms in the literature use some specific optimality properties as a basis for constructing better schedule solutions [18, 32]. Additionally, some local improvement strategies such as insert and swap methods are integrated in heuristic algorithms to gain better solutions, especially when developing metaheuristic algorithms (Li et al. [33], Imran Hossain et al. [34]). Inspired by this idea, we adopt three important optimality properties proposed by Krim et al. [29] in this paper to construct a heuristic algorithm to find optimal or near-optimal solutions for large-sizes problems. Additionally, a local improvement strategy is combined with the heuristic algorithm.First, the three important optimality properties addressed by Krim et al. [29] are briefly described as follows.Property 1.
In any optimal solution, the jobs are scheduled in increasing order ofpi/wi (WSPT rule) in each batch.Property 2.
In any optimal solution,WBk≥WBk+11≤k≤m∗−1 where WBk=∑i∈Bkwi.Property 3.
In any optimal solution, for each idle timeΔk in the batch Bk and Jh with smallest processing time in the batch Bl, k+1≤l≤m∗, the condition of Δk<ph for k∈1,…,m∗−1 is always satisfied.
Krim et al. [29] proposed nine heuristic algorithms, and the basic framework of these algorithms contains the following three stages: (1) a sequence of jobs is generated by dispatching rules; (2) given the sequence of jobs, an initial schedule with PM is produced by batching rules; and (3) Properties 1 and 2 are applied to improve the initial schedule to obtain final schedule. The detailed steps of the heuristic algorithms can be found in the study by Krim et al. [29]. According to their computational results, the WSPTBF and WSPTFF algorithms had better performances among the nine algorithms. In WSPTBF and WSPTFF algorithms, the job sequence is generated by WSPT rule, whereas the initial schedule was, respectively, generated by best fit (BF) and first fit (FF) batching rules, which are well-known batching rules for solving bin-packing problems [32].
In this paper, we develop two-phase heuristic algorithm with framework similar to that of the heuristic algorithms of Krim et al. [29]. Phase 1 consists of a sorting step using a WSPT dispatching rule to determine a sequence of jobs, followed by assigning steps using a full-batching rule that determines the jobs in each batch (i.e., the interval of T−δ). The full-batching rule is used widely to solve the bin-packing problem that tries to minimise the number of batches [32]. As a result, an initial schedule solution is obtained. Next, Phase 2 aims to improve the initial schedule by the local improvement strategy, which includes insert and swap procedures. The insert procedure attempts to satisfy Property 3 by inserting jobs to the forward batches. That is, if a job α, α∈Bk+l and batch Bk exist in which the condition of ∑i∈Bkpi+pα≤T−δ is satisfied, then job α is removed from batch Bk+l and placed into the idle period of batch Bk. The detailed steps of the insert procedure are described as follows:(Algorithm 1)
The swap procedure attempts swapping jobs belonging to two different batches so that the current idle time is decreased or the total weightsWk of batch k are increased. After swapping two jobs, Property 1 still needs to be satisfied. The detailed steps of the swap procedure are described as follows.(Algorithm 2)
Finally, a complete heuristic algorithm combined with local improvement strategy including insert and swap procedures, namedWSPTFB + LIS, is described below.(Algorithm 3)
A preliminary test is conducted using the 12-job instance of Krim et al. [29] to compare our proposed heuristic algorithm with WSPTBF and WSPTFF. Tables 3 and 4 show the information of the 12-job instance and the results obtained by the heuristic algorithms, respectively. In Table 4, the symbol ∗ indicates the solution is also optimal. For this instance, WSPTFF had a better initial solution than the other algorithms, and the initial solutions obtained by the WSPTBF and WSPTFB + LIS algorithms are the same, equal to 11131; however, the WSPTFB + LIS algorithm had better final solution. This result occurred because, without violating the three optimality properties, our heuristic algorithm could reallocate jobs between batches to obtain a better schedule by the local improvement strategy. The other algorithms (WSPTBF and WSPTFF) only tried to re-sort the sequences of jobs (or batches) to improve the solution quality, and the content (jobs) of each batch did not change.Algorithm 1:Insert procedure (Step 4 of Algorithm WSPT + LIS).
(i)
Step 4: letk = 1 and l = 1(ii)
Step 4.1. If the condition of∑i∈Bki∈Bk+pα≤T−δ is satisfied, remove job α from batch Bk+l and place it into the idle period of batch Bk, let improve_flag = true, substitute the current schedule by the new one, go to Step 3(iii)
Otherwise, go to Step 4.2./∗Property 3 check(iv)
Step 4.2. Ifk+l<m is satisfied, let l=l+1 go to Step 4.1; otherwise, go to Step 4.3(v)
Step 4.3. Ifk<m−1 is satisfied, let k=k+1,l=1, go to Step 4.1; otherwise, go to Step 5Algorithm 2:Swap procedure (Step 6 of Algorithm WSPT + LIS).
(i)
Step 6: letk = 1 and l = 1.(ii)
Step 6.1: if there exist two jobs,α∈Bk,β∈Bk+l, the conditions of ∑i∈Bk\αpi+pβ≤T−δ and pα×wα<pβ×wβ are satisfied, duplicate the current schedule π to schedule π′, and go to Step 6.2. Otherwise, go to Step 6.4.(iii)
Step 6.2: jobα in the position is interchangeable with job β in the position of schedule π′ and resorting the jobs in batch Bk′ and batch Bk+l′ based on WSPT rule, respectively. Then, a new schedule πnew is generated; the corresponding objective value is TWCnew./∗Property 1 check.(iv)
Step 6.3: ifTWCnew<TWC is satisfied, let π=πnew, TWC=TWCnew, improve_flag = true, go to Step 3. Otherwise, go to Step 6.4.(v)
Step 6.4: ifk+l<m is satisfied, let l=l+1; go to Step 6.1; otherwise, go to Step 6.5.(vi)
Step 6.5: ifk<m−1 is satisfied, let k=k+1,l=1; go to Step 6.1; otherwise, go to Step 7.Algorithm 3:WSPTFB + LIS.
(i)
Step 1: input data.(ii)
Step 2: sort the jobs according to theWSPT rule, and assign the sorted jobs into each batch based on full-batching method. Let π, TWC, and improve_flag be the obtained schedule, the corresponding objective value, and true, respectively.(iii)
Step 3: if improve_flag = false, output the scheduleπ, TWCandStop. If improve_flag = true, go to Step 4.(iv)
Step 4: execute insert procedure.(v)
Step 5: if improve_flag = true, go to Step 3; otherwise, go to Step 6.(vi)
Step 6: execute swap procedure.(vii)
Step 7: if improve_flag = true, go to Step 3; otherwise, go to Step 8.(viii)
Step 8: computeWBk,∀k=1,2,…,m.(ix)
Step 9: if the scheduleπ satisfies Property 2, go to Step 3; otherwise, sort the batches in decreasing order of WBk and improve_flag = true; go to Step 3./∗ Property 2 check.Table 3
Information of a 12-job instance of Krim et al. [29].
J1J2J3J4J5J6J7J8J9J10J11J12pi61511177131811410185wi101410141911141812191515T = 55, δ=1Table 4
Results obtained by the heuristic algorithms.
Initial solutionFinal solutionWSPTBF11131 (10526)J9,J12,J5,J10,J1,J8J2,J3,J6J11,J4,J710691 (10266)J9,J12,J5,J10,J1,J8J11,J4,J7J2,J3,J6WSPTFF9831J9,J12,J5,J10,J1,J8,J3J2,J6,J11J4,J79831J9,J12,J5,J10,J1,J8,J3J2,J6,J11J4,J7WSPTFB + LIS11131J9,J12,J5,J10,J1,J8J2,J3,J6J11,J4,J79712∗J9,J12,J5,J10,J1,J8,J3J2,J11,J4J6,J7Remark: ( ) indicates that the value in the () proposed by Krim et al. [15] was somewhat questionable.


## 5. Computational Experiments

This section consists of three experiments. Experiment1 focuses on comparing the computational efficiency of the two MIP models provided by Krim et al. [29] and ours. Experiment 2 aims to evaluate different lower bounds. Finally, the performances of the proposed heuristic algorithms are compared in Experiment 3. Since the benchmark data of Krim et al. [29] is not available at the following URL: https://sakai.uphf.fr/wiki/site/1b6900a5-38f1-4003-94e9-317898cbf168/instancesmaintenancewct.html, we generate a set of problems to carry out the computational experiments. The test-bed contains small-sized and large-sized problems. Small-sized problems consist of 5, 6, 7, …, 14, 15 jobs, and large-sized problems consider 20, 25, 30, 35, 40, …, 90, 100, 200, 300, 400, 500, 1000 jobs. The processing times of jobs from a discrete uniform distribution in the range [1, 5, 15, 20] are generated, respectively. The weight of jobs is generated from a uniform distribution in the range [1, 10]. T and δare set to [30, 50] and [3, 10], respectively. 10 instances were generated for each combination of problem size (n), processing time pi, T and δ; thus, a total of 2240 28×2×2×2×10 instances were randomly generated and tested. To make the generated instances be reasonable, the following two restricted conditions should be involved: (1) maxpi≤T−δ; (2) ∑pi>T−δ. All experiments were conducted on a PC with an Intel Xeon E-2124 3.4 GHz CPU and 32 GB RAM.Experiment 1.
Comparison of the two MIP models
In this Experiment1, we compare our model with the one provided by Krim et al. [29]. The two models are executed by IBM ILOG CPLEXOptimization Studio Version 12.7.1, and all problems are solved on the same machine for fair comparison. Additionally, the time limit of ILOG OPL is set to 7200 seconds for solving each instance. Table 5 reports the results for the two models, where MIPk and MIPc are, respectively, denoted as the model provided by Krim et al. [29] and our model. The “AveSol” and “AveLB” columns give the average the weighted completion time and average lower bound value over 80 instances, respectively. Another column, “nOPT,” denotes the number of times each model finds the optimal solution. Since the MIP model may not find optimal solutions within the time limit of 7200 seconds, we record the lower bound values obtained by each MIP model to make sure the identified solutions are optimal, that is, SolMIP=LBMIP. For the MIPk model, nOPT for 13-job instances is 79 (53), where 79 indicates the number of times each solution of MIPk model is equal to lower bound values ofMIPcmodel over 80 instances; furthermore, (53) means the number of times each solution of MIPk model is equal to lower bound values of MIPk model. From this table, we can observe that AveLB obtained by the proposed model MIPC improves significantly as the number of problem size increases, which allows the MIPC model to find better solutions within the time limit of 7200 seconds. Additionally, because the values of LBMIPcare tighter than those of LBMIPk, using the condition of SolMIPk=LBMIPc, feasible solutions obtained by MIPk model within the time limit of 7200 seconds can be proved to be optimal. For 15-job problem, Nopt value of MIPk is increased from 1 to 76, as shown in Table 5.
Regarding computation times, the computational effort required by theMIPk model is significantly higher than that of the MIPc model, especially when the problem sizes increases. Table 5 and Figure 4 reveal that the efficiency of the MIPc model is better than results obtained by the MIPk (Krim et al. [29]). However, the two models cannot solve the NP-hard problems considered in this paper in a reasonable time, and the performances of the two models become worse as the number of jobs increases. The comparison results showed that the two models are computationally inefficient when the number of jobs is beyond 20. Therefore, it is necessary to develop heuristics to obtain near optimal solutions in reasonable time.Table 5
Comparison results of the Krim model and our model.
nMIPkMIPcAveSolAveLBNoptAveSolAveLBNopt5637.763637.76380 (80)637.763637.763806916.075916.07580 (80)916.075916.0758071123.3881123.38880 (80)1123.3881123.3888081497.7251497.72580 (80)1497.7251497.7258091808.9131808.91380 (80)1808.9131808.91380102210.2002210.20080 (80)2210.2002210.20080112623.4502623.45080 (80)2623.4502623.45080123091.4753091.47580 (80)3091.4753091.47580133528.4253446.01379 (53)3528.4133528.41380144155.1753469.96377 (9)4153.6254153.62580154729.3253369.71376 (1)4724.0134724.01380207875.8632942.80020 (0)7862.4387622.038332512471.2252464.22510 (0)12368.56311775.16315Figure 4
Average computational times required by the two models.Experiment 2.
Comparison of the lower bounds
This experiment aims to evaluate the performance of the lower bounds, so we compare the lower bounds values with the optimal solutions for small-size problems. For large-size problems, the different lower bound values are compared with the value ofLBbest, where LBbest=maxLB1,LB3,LB4,LB5,NewLB; that is, the related percentage gap (RPG) is adopted as follows:
RPG = RPG=Best−LB/Best∗100% where Best indicates the best solutions obtained by the MIPc model or LBbest.
Regarding the computational time, excludingLB3, the average computational time required by the lower bound methods, less than 1 second, is negligible. The LB3 method applied the dynamic programming method to obtain a lower bound value for the special case in which only one PM activity exists; therefore, LB3 requires much more computational time as the number of jobs increases. In our computational experiment, the average computational time required by LB3 is approximately 1700 seconds for 1000 jobs. The computation time required by all lower bound methods excluding LB3 is very small; therefore, we do not report it here.
Tables6–8 show the minimum, maximum, and average RPG, respectively, of the lower bounds relative to the optimal solutions over 80 instances in each set of small-sized problems. From Tables 6–8, we can find that the minimum RPG, maximum RPG, and average RPG of NewLB are considerably smaller than others. Additionally, the results also supported the result of Krim et al. [29] in which LB1 and LB3 are quite tight compared to LB4. The overall average RPG of NewLB is at 5.403% on the optimality gap. For large-sized problems, Table 9 shows that NewLB remains the best lower bound. From the results shown in Tables 6–9, we observe that the NewLB is quite superior to others. To that end, NewLB is used to compute the performances of the heuristic algorithms in Experiment 3.Table 6
Best performance for lower bounds compared to the optimal solutions.
nMinimum RPG (%)LB1LB3LB4LB5NewLB50.5880.00028.01624.3070.00060.6120.00038.78226.0990.00070.9920.00044.01344.0130.00080.5390.00049.65336.0250.00090.6760.00051.38549.4970.000101.5850.00049.19538.4520.000111.3740.52753.66350.8150.000122.2600.88860.48255.1030.000130.4630.46350.21042.3440.000141.7180.00061.76055.8800.000153.0641.00064.29351.6600.180Table 7
Worst performance for lower bounds compared to the optimal solutions.
nMaximum RPG (%)LB1LB3LB4LB5NewLB542.42434.29468.23562.1990.000642.62034.30466.68364.97622.162741.72334.46771.65568.75017.949848.66041.91076.70371.56420.907946.65939.08275.13274.59230.5861043.39837.13178.23777.44122.6131145.37039.55777.74777.74724.8201240.32932.90978.52378.52320.9861341.05035.12778.99378.99324.6031438.36829.14980.15680.15622.0071542.54637.53480.65580.65526.603Table 8
Average performance for lower bounds compared to the optimal solutions.
nAverage RPG (%)LB1LB3LB4LB5NewLB511.9479.51348.14846.9102.226612.4649.28855.26354.2842.759713.3169.49358.75858.0733.933814.2189.81862.60761.6354.425914.5549.84864.50963.8545.2681014.7999.84567.26066.1125.7741115.01610.09968.98368.2055.8671215.47910.33570.49169.9656.8181315.49610.24971.31070.6767.1301415.76210.52073.08072.1287.3891516.03911.08573.56072.5837.845Overall average14.46310.00864.90664.0395.403Table 9
Average performance for lower bounds compared to the best lower bound.
nLB1LB3LB4LB5NewLB208.4643.21874.61273.5510.000258.0043.20877.45076.5950.000307.7273.40778.79877.8630.000357.3843.35679.71278.7380.000407.3603.69381.07780.1010.000457.0163.61881.75181.1740.000507.4144.36782.48381.7120.000607.0494.32083.13082.4460.000707.0904.71583.55282.9430.000807.0584.94584.30683.6220.000907.0405.07884.26783.5500.0001007.0925.32584.88284.2750.0002006.9966.03386.23285.6690.0003007.0936.42686.78286.1250.0004007.0196.50786.79786.1850.0005006.9336.52387.10186.5710.00010007.0886.87987.31586.6660.000Experiment 3.
Comparison of the heuristic algorithms
In this experiment, our proposed heuristic algorithm is compared with the two heuristic algorithms provided by Krim et al. [29], and we adopt the RPD (Relative Percentage Deviation) to measure the solution quality of each instance for the heuristic algorithms. RPD is defined as follows:RPD=TWCmethod−Min/Min×100%, where Min is the optimal or lower bound value provided by MIPc model or NewLB method, and TWCmethod is the total weighted completion time obtained by each heuristic algorithm. Tables 10 and 11 provide the results of the WSPTFB and WSPTFB + LIS algorithms, where WSPTFB represents the proposed heuristic algorithm without local improvement strategy. From Tables 10 and 11, algorithm WSPTFB + LIS give notable results by local improvement strategy in which WSPTFB + LIS has found the smaller average RPD, and the number of optimal solutions obtained is more than WSPTFB. Additionally, we also applied ANOVA analysis to examine the differences of the average RPD between the two heuristic algorithms. The results in Table 12 show that a statistically significant difference exists between the two algorithms, because p-value is less than 0.05. Based on the preliminary results, we focus only on comparing WSPTFB + LIS with other heuristic algorithms later.
To the best of our knowledge, the most recent and best algorithms proposed by Krim et al. [29] for solving the 1,hknr∑wiCi problem are the WSPTBF and WSPTFF algorithms. Thus, we adopt these two algorithms as benchmarks and compare the results with those of our WSPTFB + LIS algorithm. For fair comparison, WSPTBF and WSPTFF are recoded in C++ based on their pseudo codes (Krim et al. [29]), and the same test-bed problems are examined under the same computer configuration. In Table 13, we present the minimum, maximum, and average RPD values and the number of optimal solutions found by each heuristic algorithm over the 80 instances of each small-size problem. From this table, we can observe that WSPTFB + LIS is able to obtain better results in all cases (552 optimal solutions found over 880 instances). Additionally, the average RPD values of WSPTFB + LIS are considerably smaller than those of other heuristic algorithms, which is also seen in the results presented in Table 14 in which WSPTFB + LIS has smaller average RPD values for large-size problems. Table 15 also shows that there are significant differences on the average RPD values among the three algorithms based on the result of ANOVA analysis, where p-value is less than 0.05.
In terms of computation times, Table16 shows the proposed WSPTFB + LIS for which the required computational effort is slightly higher than that for the others as the number of jobs increases. This is because it performs the local improvement strategy to find better solutions. Nonetheless, the proposed WSPTFB + LIS can obtain very good solutions within extremely small computational times as shown in the above results, where in Table 16 the average computational effort required by our proposed algorithm is less than 2 seconds for solving 1000 jobs. This indicates that our proposed WSPTFB + LIS algorithm outperforms the WSPTBF and WSPTFF algorithms, and the proposed algorithm can apply to real-world problems.Table 10
Results ofWSPTFB and WSPTFB + LIS for small-sized problems.
nWSPTFBWSPTFB + LISMinMaxAveNoptMinMaxAveNopt50.0039.493.83530.0010.870.227760.0039.415.11430.0011.180.467270.0037.775.94340.0020.811.186580.0040.747.07280.0013.371.015890.0033.538.31190.0019.811.4860100.0035.078.98140.0022.052.0154110.0036.8110.33130.0022.931.9144120.0039.0610.8570.0018.071.9235130.0029.949.7770.0019.131.5143140.0035.5710.6260.0013.942.3826150.0039.7810.8920.0025.272.3018Table 11
Results ofWSPTFB and WSPTFB + LIS for large-sized problems.
nWSPFFBWSPTFB + LISMinMaxAveMinMaxAve200.5275.0725.410.4634.9712.20252.9579.3226.301.3537.3013.73304.2168.9528.012.5939.9614.60355.3669.8429.412.8136.9815.22406.8172.2029.483.1743.0015.24456.0371.6030.903.3741.4815.60506.5177.6030.573.1442.5515.79607.1771.4230.693.7740.7215.61707.7073.4130.904.0438.9915.52808.5967.6831.683.6237.8915.71908.7573.1031.223.9240.7415.621009.9968.3232.423.9339.1215.672009.3270.8733.164.0040.5215.663009.7569.9533.074.1240.2615.4040010.4471.5833.844.1539.0315.4350011.6869.4234.154.2338.9115.57100012.1070.6634.534.2438.0315.25Table 12
ANOVA results based on the averages of the proposed heuristic algorithms.
SourceDfSum of squaresMean squareFp valueFactor12103.8162103.81623.436<0.001Error544847.51989.769Total55Table 13
Results of the heuristic algorithms compared to the optimal solutions.
nWSPTBFWSPTFFWSPTFB + LISMinMaxAveNoptMinMaxAveNoptMinMaxAveNopt50.0039.492.89580.0024.921.81590.0010.870.227760.0035.883.70460.0039.412.30490.0011.180.467270.0035.245.42370.0031.352.81480.0020.811.186580.0029.115.98300.0027.122.59410.0013.371.015890.0028.797.09210.0024.203.30380.0019.811.4860100.0034.008.18170.0022.054.21260.0022.052.0154110.0033.808.65140.0024.084.53190.0022.931.9144120.0027.039.4670.0022.893.90190.0018.071.9235130.0026.327.90110.0019.663.46200.0019.131.5143140.0027.679.1360.0016.894.57140.0013.942.3826150.0038.939.6930.0026.654.4360.0025.722.3018Table 14
Results of heuristic algorithms compared to the lower bounds.
nWSPTBFWSPTFFWSPTFB + LISMinMaxAveMinMaxAveMinMaxAve200.5255.9722.900.5240.4916.090.4634.9712.20252.9566.4724.501.5849.5916.741.3537.3013.73304.2163.5924.853.0746.0517.422.5939.9614.60355.3657.7926.523.3745.9717.862.8136.9815.22406.5667.7227.113.5046.6317.973.1743.0015.24456.0367.2427.823.5144.2918.273.3741.4815.60506.5173.4328.243.5344.6518.883.1442.5515.79607.1760.1027.774.3444.6918.423.7740.7215.61707.7065.0928.114.9143.9418.194.0438.9915.52807.1863.8228.834.5542.3318.553.6237.8915.71908.7565.9528.514.2843.1518.463.9240.7415.621009.9561.8429.314.4942.3618.513.9339.1215.672009.1665.6329.864.7542.0118.314.0040.5215.663008.9963.8329.914.7041.4017.924.1240.2615.404009.8764.6330.484.6441.0818.044.1539.0315.4350010.8564.1930.904.8040.1018.074.2338.9115.57100011.4465.0431.064.7240.2117.794.2438.0315.25Table 15
ANOVA results based on the averages among the three heuristic algorithms.
SourceDfSum of squaresMean squareFp valueFactor21526.017763.00910.728<0.001Error815761.03871.124Total83Table 16
Average computation time of the heuristic algorithms.
nWSPTBFWSPTFFWSPTFB + LIS50.000000.000000.0000060.000000.000000.0000070.000000.000000.0000080.000000.000000.0000090.000000.000000.00000100.000000.000000.00000110.000000.000000.00000120.000000.000000.00000130.000000.000000.00000140.000000.000000.00000150.000000.000000.00000200.000000.000000.00001250.000000.000000.00001300.000000.000010.00002350.000010.000010.00003400.000010.000010.00004450.000010.000010.00006500.000010.000010.00013600.000010.000020.00014700.000020.000020.00024800.000020.000030.00031900.000020.000030.000461000.000030.000040.000662000.000100.000130.005133000.000200.000260.020554000.000320.000410.067355000.000460.000600.1307710000.001590.002111.54310


## 6. Conclusions

Scheduling problems with PM activity have received increasing attention in the last decade, due to the enthusiastic adoption of the PM philosophy in many real manufacturing systems. However, not much work has been performed for problems with the objective of total weighted completion time, denoted as1,hknr∑wiCi. For the 1,hknr∑wiCi problem, we proposed a new MIP model, two new lower bounds, and a heuristic algorithm. The proposed heuristic algorithm adopted a bin-packing scheme that considers batch formation and batch sequence and applied the optimality properties in the local improvement strategy to efficiently obtain optimal or near optimal solutions. To the best of our knowledge, a recent new study on the same problem was conducted by Krim et al. [29]. Thus, we took their proposed methods including the MIP model, lower bound methods, and better heuristic algorithms as the benchmark approaches.Based on the results of a comprehensive comparison experiment with the benchmark approaches, in terms of MIP models, the number of optimal solutions found by our MIP model was more than that of the benchmark model within the same time limit (7200 seconds). Regarding lower bound methods,NewLB was significantly more effective than other methods, as the average RPG of NewLB was smaller (closer to optimal solutions/best lower bounds) than that of other methods. For solving small-sized problems, approximately 63% of our proposed heuristic solutions were found to be optimal, a ratio that is higher than that of the other algorithms (28% and 39% solutions found to be optimal for the WSPTBF and the WSPTFF, respectively). Furthermore, the results for WSPTFF + LIS were significantly better than the benchmark heuristic results for the average RPD values (reaching lower bound values) for larger instances. The average computation time of our heuristic algorithm was less than 2 seconds. These results show that our proposed heuristic algorithm is very efficient and has better performance with respect to problems with different dimensions.In summary, the results presented in this paper are very encouraging regarding the use of the proposed MIP model and theWSPTFF + LIS algorithm for the single-machine total weighted completion time problem with periodic preventive maintenance. In future work, it will be worth considering different types of maintenance activities, such as variable maintenance duration or flexible maintenance activities, or extending this problem to parallel machines or, flow or job shop environments, as well as including different operation constraints. Additionally, another extension may be a consideration of other criterion such as total tardiness-earliness cost, and total weighted tardiness, which are also important objective functions for system performance measures.Finally, benchmark data sets and solutions for our proposed methods including the MIP model,NewLB, and WSPTFF + LIS are made available at the following URL: https://drive.google.com/file/d/1Fh5S7j6pHGHhKuy2O6WcO0yfl2eMno-t/view?usp=sharing for researchers to compare their solutions.



---

*Source: 8891322-2021-04-17.xml*