A general variable neighborhood search for single-machine total tardiness scheduling problem with step-deteriorating jobs

In this article, we study a single-machine scheduling problem of minimizing the total tardiness for a set of independent jobs. The processing time of a job is modeled as a step function of its starting time and a specific deteriorating date. A mixed integer programming model was applied to the problem and validated. Since the problem is known to be NP-hard, we proposed a heuristic named simple weighted search procedure (SWSP) and a general variable neighborhood search algorithm (GVNS). A perturbation procedure with 3-opt is embedded within the GVNS process in order to explore broader spaces. Extensive numerical experiments are carried out on some randomly generated test instances so as to investigate the performance of the proposed algorithms. By comparing to the results of the CPLEX optimization solver, the heuristic SWSP and the standard variable neighborhood search, it is shown that the proposed GVNS algorithm can provide better solutions within a reasonable running time.


Introduction
As an important tool for manufacturing and engineering, scheduling has a major impact on the productivity of a process. In the classical model of scheduling theory, it is assumed that the processing times of all jobs are known in advance and remain constant during the entire production process. However, this assumption may not be applicable to model some real manufacturing and service problems. Examples can be found in steel production, equipment maintenance, cleaning tasks allocation and so forth [21,28], where any delay or waiting in starting to process a job may increase the necessary time for its completion. Such kinds of jobs are called deteriorating jobs. This kind of scheduling problems with deteriorating jobs was firstly considered by Browne and Yechiali [2] and Gupta and Gupta [8]. They assumed that the processing time of a job is a single linearly non-decreasing function of its starting time and aimed at minimizing the expected makespan and the variance of the makespan in a single-machine environment. The function that describes the processing time of a deteriorating job shall henceforth be called a deterioration function. Since then, there is a rapidly growing interest in the literature to study various types of scheduling problems involving deteriorating jobs. A recent survey of scheduling problems with deteriorating jobs was provided by Cheng et al. [5].
Most of the research concentrated on obtaining the optimal schedule for scheduling jobs for which their processing times continuously depend on their starting times [1,31,35,17,34,13].
However, if some jobs fail to be processed prior to a pre-specified threshold, their processing times will be extended by adding an extra penalty time in some situations. Job processing times of this type of deterioration may be characterized by a piecewise defined function. Kunnathur and Gupta [21] firstly studied single-machine scheduling with a piecewise deterioration function. They assumed the actual processing time of a job is split into a fixed part and a variable penalty part. The variable part depends linearly upon the starting time of the job. Sundararaghavan and Kunnathur [33] introduced a total weighted completion time scheduling problem with step deterioration. Then, Mosheiov [29] investigated the scheduling problem with step deterioration in which the objective is to minimize the makespan on a single machine or multiple machines. He also introduced a simple heuristics for all these NP-hard problems.
For the scheduling problem with a piecewise linear deterioration function, Kubiak and Van de Velde [20] considered single-machine scheduling with a generalization of the unbounded deterioration model proposed by Browne and Yechiali [2] and presented NP-hardness proofs. They also developed a pseudo-polynomial dynamic programming algorithm and a branch and bound algorithm to obtain the optimal schedule. Alternatively, Cheng et al. [4] studied the problem of scheduling jobs with piecewise linear decreasing processing times on a single or multiple machines. They aimed to minimize the makespan and the flow time, i.e., the total completion time, and proved that the two problems are NP-hard. Subsequently, Ji and Cheng [16] gave a fully polynomial time approximation scheme for the same case with a single-machine. Afterwards, Wu et al. [36] provided two heuristic algorithms to solve the single-machine problem under the piecewise linear deterioration model. Moslehi and Jafari [30] dealt with the single-machine scheduling problem with the assumption of piecewise-linear deterioration. They suggested a branch and bound algorithm and a heuristic algorithm with complexity O(n 2 ) for the objective of minimizing the total number of tardy jobs, where n is the total number of jobs in the problem.
Some articles focused on the scheduling problem with step-deterioration have been published. Cheng and Ding [3] studied some single-machine scheduling problems and showed that the flow time problem is NP-complete. Jeng and Lin [14] introduced a branch and bound algorithm for the single-machine problem of minimizing the makespan. In addition, They also studied the same problem for minimizing the flow time [15]. He et al. [11] proposed an exact algorithm to solve most of the problems with up to 24 jobs and a heuristic algorithm to derive a nearoptimal solution. Layegh et al. [22] minimized the total weighted completion time by using a memetic algorithm. Cheng et al. [6] developed a variable neighborhood search to minimize the flow time on parallel machines. Furthermore, batch scheduling with step-deterioration has got attentions [24,27]. However, the total tardiness as one important objective in practice has not been concerned on the studies of single-machine scheduling problem with step-deterioration.
Minimizing total tardiness has been drawing considerable attentions of researchers in the past decades and it is NP-hard in the strong sense [7]. The latest theoretical developments for the single-machine scheduling problem and the current state-of-the-art algorithms are reviewed by Koulamas [19]. He found that the single-machine total tardiness scheduling problem continues to attract researchers' interests from both theoretical and practical perspectives. In this article, we consider the single-machine total tardiness scheduling problem with step-deteriorating jobs. Since the total tardiness problem without step-deterioration in a single-machine is a NP-hard problem, the problem tackled here is also a NP-hard problem. To the best of our knowledge, there is no literature on minimizing the total tardiness in a single-machine scheduling problem with step-deteriorating jobs.
The remainder of the study is organized as follows. Section 2 contains the problem description and a mixed integer programming model. In Section 3, a heuristic and a general variable neighborhood search algorithm are developed for the proposed problem. The proposed methods are tested and compared to the CPLEX and the variable neighborhood search (VNS) on variance test instances of both small size and large size in Section 4. In the last section, conclusions and future works are mentioned.

Problem formulation
The single-machine total tardiness scheduling problem with step deteriorating jobs can be described as follows. Let N := {1, 2, . . . , n} be the set of n jobs to be scheduled on a singlemachine without preemption. Assume that all jobs are ready at time zero and the machine is available at all times. In addition, the machine can handle only one job at a time, and cannot keep idle until the last job assigned to it is processed and finished. For each job j ∈ N , there is a basic processing time a j , a due date d j and a given deteriorating threshold, also called deteriorating date h j . If the starting time s j of job j ∈ N is less than or equal to the given threshold h j , then job j only requires a basic processing time a j . Otherwise, an extra penalty b j is incurred. Thus, the actual processing time p j of job j can be defined as a step-function: p j = a j if s j h j ; p j = a j + b j , otherwise. Without loss of generality, the four parameters a j , b j , d j and h j are assumed to be integers.
Let π = (π 1 , . . . , π n ) be a sequence that arranges the current processing order of jobs in N , where π k , k = 1, . . . , n, indicates the job in position k. The tardiness T j of a job j in a schedule π can be calculated by T j = max {0, s j + p j −d j }. The objective is to find a schedule π * such that the total tardiness T j is minimized. The total tardiness is a non-decreasing criterion of the job completion times. Using the three-field notation, this problem can be denoted by 1|p j = a j or a j + b j , h j | T j .
Based on the above description, we formulate the problem as a 0-1 integer programming model. Firstly, the decision variable y ij , i, j ∈ N is defined such that y ij is 1 if job i is scheduled before job j on the single-machine, and 0 otherwise. For each pair of jobs i and j, y ij + y ji = 1. Then, for a schedule π of the jobs in N , we minimize where M is a large positive constant such that M → ∞ as n → ∞. For example, M may be chosen as M := max j∈N {d j } + j∈N (a j + b j ) .
In the above mathematical model, equation (2.1) represents the objective of minimizing the total tardiness. Constraint (2.2) defines the processing time of each job. Constraint (2.3) determines the starting time s j of job j with respect to the decision variables y ij . Constraint (2.4) defines the tardiness of job j. Finally, constraints (2.5) and (2.6) define the boundary values of variables y ij , s j , T j , for i, j ∈ N .

Solution methodologies
Due to the NP-hardness of the considered problem, only small size instances can be solved optimally by the enumeration techniques such as branch and bound or dynamic programming. However the size of some practical problem encountered in industry is usually large. We need to develop feasible heuristics to deal with large-sized instances. In this section, a simple weighted search procedure and a general variable neighborhood search algorithm are proposed.
Before introducing the algorithms, we discuss two properties of the considered problem that will be used later.
Property 3.1. Consider any two jobs j and k that are processed before their given deteriorating dates. If p j ≤ p k and d j ≤ d k , then job j must be scheduled before job k in an optimal sequence.
Proof. Since the two jobs under consideration are not deteriorated, their processing times are only the basic processing times. Note that they remain to be not deteriorated even if the order of processing of the two jobs is switched. Lemma 3.4.1 [32, p. 51] regarding the total tardiness problem in a single-machine environment says if p j p k and d j d k then there exists an optimal sequence in which job j is scheduled before job k. Thus, applying this Lemma based on the basic processing times and the due dates yields the optimal sub-sequence including the job pair (j, k) in an optimal schedule. Property 3.2. Consider any two jobs j and k that are processed after their given deteriorating dates. If p j + b j ≤ p k + b k and d j ≤ d k , then job j should be scheduled before job k in an optimal sequence.
Proof. Since the two jobs j and k are processed after their deteriorating times, the actual processing time of each job is the sum of its basic processing time and a penalty time. Thus, like Property 3.1, Property 3.2 is a direct result of applying Lemma 3.4.1 of [32] to all such job pairs (j, k) based on their actual processing times (a j + b j , a k + b k ) and the due dates (d j , d k ).
According to the two properties, a job should be scheduled in an earlier position if it has a smaller value of a i (or a i + b i if it is deteriorated) and a smaller value of d i . There is a weighted search algorithm that has applied to solve the single-machine flow time scheduling problem [11]. Thus, we adopt this idea and provide a simple weighted search procedure (SWSP) to obtain a heuristic for the 1|p j = a j or a j + b j , h j | T j problem.

Simple weighted search procedure
We shall need an appropriately chosen triple ω := (ω 1 , ω 2 , ω 3 ) of positive weights to compute a weighted value m i := ω 1 d i +ω 2 p i +ω 3 h i , for each i ∈ N . Since it is difficult to determine a priori a triple of weights that can yield good solutions, we adopt the dynamic update strategy to search all possible triples ω. Weight ω 1 is a linearly increasing weight, given by ω 1 = ω 1 min + (ω 1 max − ω 1 min )(l 1 − 1)/(n − 1), where l 1 varies from 1 to n. Weight ω 2 adopts the same updating formula as weight ω 1 . The only difference between the two weights is their possible minimal and maximal values. Suitable values of ω 1 min , ω 1 max , ω 2 min and ω 2 max may be determined by preliminary experiments. For example, for the numerical experiments we conducted in this article, we have chosen the values of the four parameters ω 1 min , ω 1 max , ω 2 min and ω 2 max to be 0.2, 0.9, 0.1 and 0.7, respectively. With ω 1 and ω 2 determined, ω 3 := 1 − ω 1 − ω 2 . But ω 3 may become negative by using this formula. Should a negative ω 3 occurs, it is adjusted to a pre-specified positive value. In this article, this preset value is chosen to be 0.1. Let Ω be the set of all possible triples of weights obtained.
We now describe the SWSP. Firstly, let a triple (ω 1 , ω 2 , ω 3 ) ∈ Ω be fixed. Then the job with the minimal due date is scheduled in the first position. Subsequently the unscheduled jobs are arranged in a nondecreasing order of the weighted sums m i , i ∈ N . Thus each triple (ω 1 , ω 2 , ω 3 ) ∈ Ω is associated with a schedule. Correspondingly, for each schedule, the value of the objective total tardiness is calculated. Among all these different solutions, the best solution is given by the one with the smallest total tardiness. Furthermore, in order to refine the quality of the solution, a local improvement approach based on pairwise swapping of jobs is used in the procedure. The detailed SWSP is described in Algorithm 1.

Algorithm 1 Simple weighted search procedure (SWSP)
1: Input the initial data of a given instance; 2: Set c [0] = 0, N 0 = {1, · · · , n} and π = [ ]; 3: Generate the entire set Ω of possible triples of weights. For each triple (ω 1 , ω 2 , ω 3 ) ∈ Ω, perform the following steps; 4: Choose job i with minimal due date to be scheduled in the first position;   end if 24: until the set N 0 is empty 25: calculate the total tardiness Z(π) of the obtained schedule π; 26: for i ← 1 to n do 27: for j ← 1 to n do 28: if i = j then 29: create a new solution π ′ by interchanging jobs in the ith and jth position from π; 30: replace π by π ′ if the total tardiness of π ′ is smaller than that of π; 31: end if 32: end for 33: end for 34: Output the finial solution π.

General variable neighborhood search heuristic
The variable neighborhood search (VNS) is a simple and effective meta-heuristic proposed by Mladenović and Hansen [25]. It has the advantage of avoiding entrapment at a local optimum by systematically swapping neighborhood structures during the search. The basic VNS combines two search approaches: a stochastic approach in the shaking step that finds a random neighbor of the incumbent, and a deterministic approach which applies any type of local search algorithm. Serval VNS variants have been proposed and successively applied to numerous combinatorial optimization problems [10]. Among which, the variable neighborhood descent (VND) [9] as the best improved deterministic method is the most frequently used variant. The VND initiates from a feasible solution as the incumbent one and then carries out a series of neighborhood searches through operators N k ,(k = 1, · · · , k max ). If a better solution is found in a neighborhood, the incumbent solution is replaced by the better one and the search continues within the current neighborhood; otherwise it will explore the next neighborhood. The obtained solution at the end of the search is a local optimum with respect to all neighborhood structures. If the deterministic local search approach of the basic VNS is replaced by the VND search, the resulted algorithm is called a general VNS (GVNS). The GVNS has received attentions and showed good performance compared to other VNS variants. So far, the GVNS has been applied to many optimization problems, such as the incapacitated single allocation p-hub median problem [12], the one-commodity pickup-and-delivery traveling salesman problem [26], the capacitated vehicle routing problem [23] and the single-machine total weighted tardiness problem with sequence dependent setup times [18].
As far as we know, there is no published work on solving scheduling problems with deteriorating jobs by the GVNS. Therefore, in this section we develop a GVNS heuristic for the single-machine total tardiness scheduling problem with step-deterioration. In what follows, the main components of the designed GVNS algorithm are described: the initialization process, five neighborhood structures provided for the shaking step and local search procedures. In addition, a perturbation phase with a 3-opt operation without change of direction is embedded within the search process in order to decrease the probability of returning to the previous local optimum.

Initialization
Typically, a good initial sequence can be obtained by using the Earliest Due Date (EDD) rule. That is to say, all jobs are arranged in nondecreasing order of the due dates.

Neighborhood structures
In a local search algorithm, a neighborhood structure is designed by introducing moves from one solution to another. In order to conduct a local search in the proposed GVNS, we next develop five neighborhood structures based on swap, insertion and k-opt for the problem 1|p j = a j or a j + b j , h j | T j . These neighborhood structures are defined by their corresponding operators. We remark that for simplicity a neighborhood structure shall be referred to by the name of its corresponding generating operator. For example, if a neighborhood structure is generated by an operator N 1 , the neighborhood structure shall be called N 1 neighborhood structure or simply N 1 neighborhood.
Swap Operator (N 1 ): The swap operator (Figure 3.1) selects a pair of jobs π i and π j in the current sequence π of jobs, exchanges their positions, and computes the total tardiness of the resulting solution. If an improvement is found, the new sequence is accepted as the incumbent one. The operator repeats this process for all indices i,j ∈ N , with i = j until all the neighborhoods have been searched. The size of the swap operator is n(n − 1)/2.
Insertion Operator (N 2 ): For a given incumbent arrangement of jobs, the insertion neighborhood (Figure 3.2) can be obtained by removing a job from its position and inserting it into  The pairwise exchange operator is similar to the swap operator except that the swap is applied to a pair of jobs rather than a single job. In our implementation, these pairwise jobs are selected from all the combinations of two out of n jobs. It exchanges their positions and computes the total tardiness of the resulting sequence. Once an improvement is obtained, the incumbent solution is updated by the new solution.
Couple Insertion Operator (N 4 ): This procedure is similar to the insertion operator N 2 . But the operation is for a pair of successive jobs. For each couple of jobs π i and π i+1 , 1 ≤ i < n − 1, the operator extracts these two jobs and inserts them in another pair of positions j and j + 1, 1 ≤ j ≤ n − 1. Note that i = j.
2-opt Operator (N 5 ): The 2-opt is the most classical heuristic for the traveling salesman problem in which it removes two edges from the tour and reconnects the two paths created. In our implementation (Algorithm 2), the 2-opt operator selects two jobs π i and π j in the current sequence. It then deletes the edge connecting job π i and its successor and the edge connecting job π j with its successor. Afterwards, it constructs a new connection of π i to π j and a new connection between their respective successors. Furthermore, the partial sequence between the successors of π i and π j is reversed. Figure 3.3 illustrates an example of the 2-opt procedure. When an improvement is found in terms of the objective value, the incumbent is updated with the improved sequence. The search continues with applying the 2-opt operator to all possible pairs of jobs that are at least three positions away from each other. The solution obtained by this operator may not significantly differ from the incumbent in terms of the objective values. But some random jumps in the objective value may be achieved to escape from a current local optimum due to the path reversion.

Shaking and local search
A shaking operation is performed before the local search in each iteration. The shaking procedure plays a role in diversifying the search and facilitating escape from a local optimum. To preserve the simplicity of the principles of the VNS, both the shaking procedure and the local search of the GVNS make use of the same set N k (k=1,. . . , 5) of neighborhood structures. The apply the 2-opt procedure to the jobs π i and π j to generate a new sequence π ′ ; 7: if the new sequence π ′ is better than the current sequence in terms of objective value then 8: update the incumbent;  shaking operation is implemented by generating randomly a neighboring solution π ′ of the current one π using a given neighborhood operator N k . The neighborhood operator N k will be chosen to cycle from N 1 through N 5 . If the given neighborhood structure is N 1 or N 5 , then the shaking procedure randomly selects two jobs. If the given neighborhood structure is N 2 , then the shaking procedure randomly selects a job and an insertion position. If the given neighborhood structure is N 3 , then the shaking procedure randomly selects two pairs of jobs. If the given neighborhood structure is N 4 , then the shaking procedure randomly selects a couple of two consecutive jobs and a couple of two consecutive positions. A complete local search is organized as a VND using the proposed neighborhood structures. To efficiently explore possible solutions, a sequential order K of applying these neighborhoods are randomly generated. For instance, assuming that the sequence is K = (3, 1, 2, 5, 4), the search starts from N 3 and ends at N 4 . For each neighborhood, a new local optimum π ′′ is obtained by carrying out the corresponding local search operation. If π ′′ is better than π ′ , the new solution π ′′ will be accepted as a descent so that π ′ is updated with π ′′ , and the search continues for a new π ′′ within the current neighborhood; otherwise the search turns to the next neighborhood in K. The search stops until the last neighborhood in K is explored. It is worthwhile to point out that after one iteration of a shaking and local search, the solution π ′ generated by the shaking procedure may not be improved.
The shaking and local search steps are summarized below.
1. Begin the shaking procedure. Use the current solution π to randomly generate a neighbor π ′ and set i = 1. 2. Obtain the sequential order K of neighborhood search at random.  3. Apply the neighborhood operator N K(i) to π ′ to achieve a local optimum π ′′ , where K(i) means the i-th entry of the sequence K. 4. If π ′′ is better than π ′ , update π ′ and keep the value i so the search will continue within the current neighborhood; otherwise, i ← (i + 1). 5. If i > k max = 5, complete one iteration of the shaking and local search; otherwise, go to Step 3.

Perturbation phase
Since a local search applied to optimization problems often suffers from getting trapped in a local optimum, the well-know approach for this deficiency is to adopt a multi-start method when no improvement is observed. Multi-start heuristics can usually be characterized as iterative procedures consisting of two phases: the first phase in which a starting solution is generated and the second one in which the starting solution is typically (but not necessarily) improved by local search methods. However, a far more effective approach is to generate the new starting solution from an obtained local optimum by a suitable perturbation method.
Inspired by the multi-start strategy, we provide a perturbation phase with 3-opt operator without change of direction for producing a new starting solution in case a local search is trapped in a local optimum. In our implementation, if the current solution π cannot be further improved through a predetermined number γ of iterations of the shaking and local search procedures, the GVNS algorithm assumes that there is no hope to continue the local search based on the current best solution π. Then, the 3-opt operator (Figure 3.4) is implemented on π. Three jobs are randomly selected from the incumbent sequence. The 3-opt operator removes the three edges connecting the selected jobs with their successors. Then it reconnects the four subsequences created. The direction of these subsequences remains unchanged because reversing the direction of these subsequences may produce a much worse solution after a number of local searches.

Proposed GVNS algorithm
The GVNS algorithm starts with an initialization phase in which an initial solution π 0 is generated by the EDD rule and is set as the current solution π. To evaluate the quality of solutions, the solution value is considered as the total tardiness.
At each iteration, the proposed algorithm needs to carry out mainly two steps. First, a shaking procedure is performed to generate a random neighboring solution π ′ of the incumbent solution π. Subsequently, a local search based on the VND with π ′ as the input solution is performed in an attempt to improve π ′ . The local solution π ′ , possibly improved by the local search is then compared to the current best solution π in terms of the objective function value. If π ′ is better than π, the current solution π is updated by π ′ . This completes one iteration of the shaking and local search. The GVNS algorithm then continues with next iteration of shaking and local search. In addition, if no improvement of the current best solution π can be found through a predetermined number γ of iterations of shaking and local search, a 3-opt perturbation procedure is used to generate a newly starting solution π.
We next describe the stopping criterion. A pre-specified maximum number Iter max of iterations of shaking and local search or a maximum number Iter nip between two iterations without improvement are used as the stopping criterion of the algorithm. In our implementation, we choose Iter max = 500 and Iter nip = 150, respectively. Moreover, the predetermined number γ of iterations to perform a 3-opt perturbation is chosen to be γ = 0.5 · Iter nip . This strategy is determined by our preliminary experiments. A Pseudo-code of the proposed GVNS algorithm is summarized in Algorithm 3.

Algorithm 3 General variable neighborhood search heuristic
1: Initialize the parameters of the algorithm; 2: Define a set of neighborhood structures N k (k = 1, . . . , 5), that will be used in the shaking phase and the local search phase; 3: Generate the initial solution π 0 by the EDD rule; 4: Calculate the objective value f (π 0 ) for the solution π 0 ; 5: Set the current best solution π = π 0 ; 6: Choose the stopping criterion; initialize the counter: iter 1 = 0, iter 2 = 0 and iter 3 = 0; 7: Set the first neighborhood structure for the shaking procedure to be k ← 1; 8: repeat 9: Shaking: Generate a point π ′ at random in the N k neighborhood of π;

10:
Produce a random sequence K of applying the 5 neighborhood structures;

20:
if π has not improved for a given number of iterations γ, that is, when iter 3 > γ, and if iter 2 < iter N IP then 21: use the 3-opt perturbation procedure to generate a new starting solution π; and reset iter 3 = 0; 22: end if 23: Set k = k mod 5 + 1 to cycle through the neighborhood structures for shaking. 24: until the stopping criterion is met, that is, iter 1 > Iter max or iter 2 > Iter nip ; 25: Output the current best solution π.

Computational experiments
In this section, in order to evaluate the performance of the GVNS algorithm and the proposed heuristic SWSP, a set of testing instances were generated and solved on a PC with Intel Core i3 3.20 GHz CPU and 4 GB of memory in the environment of Windows 7 OS. The procedure of generating the testing instances and analysis of the results are described below.

Experiments design
In this article, the basic processing times (a j ) are randomly generated from a discrete uniform distribution on the interval (0, 100]. The deteriorating dates (h j ) are picked from the uniform distributions over three different intervals (0, A/2], [A/2, A] and (0, A], where, the value of A is obtained from the equation A = j∈N a j . The deterioration penalties b j , j ∈ N , are drawn from the uniform distribution (0, 100×τ ], where we choose τ =0.5. In addition, the due dates are randomly generated according to a discrete uniform distribution from two different intervals (0, 0.5×C max ], and (0,C max ], whereC max is the value of the makespan obtained by arranging the jobs in the non-decreasing order of the ratios a j /b j , j ∈ N [1]. The algorithms were tested over five different job sizes for small-sized instances. Those sizes are respectively n= 8, 10, 15, 20 and 25. For all possible combinations of deteriorating dates and due dates, to account for the three intervals in which deteriorating dates are generated and the two intervals in which due dates are drawn, six types of problems are needed to solve for a specific job size. For convenience, these types of problems are denoted by symbols S_k 1 k 2 . For example, S_12 represents a type of problems with deteriorating dates drawn from the interval (0, A/2] and due dates drawn from the interval (0,C max ]. In addition, a set of large sized instances with n={50, 60, 70, 80, 90, 100} are being considered. As a consequence, for both small-sized and large-sized problems, 66 ((5 + 6) × 6) sample instances were randomly generated.

Experimental results
Because there are no comparative data and no competing heuristic for our problem, comparisons with best know solutions are not possible. Thus, we have designed and coded the standard VNS heuristic (Algorithm 4) as a comparison to assess our approach. The VNS heuristic does not include the perturbation phase. The implementation sequence of neighborhoods is chosen by the deterministic order (N 1 , N 2 , N 3 , N 4 and N 5 ). The stopping criteria are chosen to be identical to the GVNS algorithm.
The described problem can be solved optimally by a commercial solver CPLEX 12.5 for some small-sized instances. However, because of its NP-hardness, it is impossible to obtain optimal solution by the CPLEX for medium or large instances. The preset run time for CPLEX is one hour. If CPLEX fails to either converge to the optimum or prove the optimality of the incumbent within the time limit, a GAP is provided by the software. The GAP shows the quality of a solution given by the CPLEX within the run time limit to some extent. A large GAP implies that the CPLEX needs more time to converge to an optimal solution.
All proposed algorithms are implemented in Matlab 7.11. The CPLEX and the SWSP are deterministic hence only one run is necessary for each problem instance. While the VNS and the GVNS are stochastic so we have to run some replications in order to better assess the experimental results. In our experiment, for each instance we run 10 times when using the VNS and the GVNS.
Each algorithm's performance is measured by computing a relative percentage deviation (RPD) defined by the equation where Z alg is the solution value obtained for a given algorithm and instance, Z best is the best solution of all approaches. For the VNS and the GVNS, the comparison results are in terms of average RPD and average computational time.

17:
Update the counter of total number of iterations iter 1 = iter 1 + 1; 18: until the stopping criterion is reached, that is, iter 1 > Iter max or iter 2 > Iter nip ; 19: Output the current best solution π. Table 1 summaries the results obtained from the CPLEX, the SWSP, the VNS and the GVNS. As shown in the table, the GVNS finds the best solutions for all small-sized instances. Since the computational time of the CPLEX was preset to one hour, the CPLEX could not find the optimal solutions for all instances with 20 or 25 jobs. The computational time of the CPLEX exponentially increases as the number of jobs increases due to the NP-hardness of the problem. Note that 4 solutions obtained by the CPLEX are worse than those obtained by the GVNS. The average RPD given by the SWSP is 11.63%. The run times of the SWSP are very short and are ignored in table 1. Compared to the CPLEX, the run times of the GVNS and the VNS do not significantly increase as the number of jobs increases. The computational time of the GVNS is longer than that of the VNS because the GVNS includes the VND scheme as the local search and additionally the perturbation phase. However, the VNS delivers the best solutions for only 13 out of the 30 small-sized test instances. The RPD delivered by VNS is only 2.32%. This suggests that these designed neighborhood structures are well suitable for solving the problem under consideration.
For large-sized instances, it is impossible to find the optimal solution by using the CPLEX within the preset one-hour time limit. Thus we tested the performance of the proposed algorithms with respect to the standard VNS. The results are shown in Table 2. The average RP D of the GVNS is about 0.78% which is significantly smaller than that of the VNS. Except for instance S_12 with 60 jobs, the RPDs of other instances given by the GVNS are very small. The relative higher RPD of instance S_12 may be because its optimal solution value is relatively small and there are 2 local optima found in 10 replications. According to figure 4.1, it is observed that the quality of solutions delivered by the GVNS is very good without respect to the distribution of the deteriorating dates. On the other hand, the quality of solutions obtained by the SWSP is inferior compared to the GVNS and the VNS.
To measure the robustness of the algorithm, a mean absolute deviation (MAD) of 10 runs when applying the VNS or the GVNS was calculated for each of large-sized instances according to the equation where R is the number of replications,Z alg denotes the average solution value obtained from the R runs for a given algorithm, and Z alg (r) indicates the solution value obtained for a given algorithm in the r-th run. The average M AD of the GVNS is merely 0.81%. Figure 4.2 reports the M AD obtained by the GVNS for large-sized instances. Overall, the GVNS demonstrates to be a good alternative to solve single-machine total tardiness scheduling problems with stepdeteriorating jobs, because the algorithm finds good solutions in a reasonable computational time.

Conclusions
In this article, we considered a single-machine total tardiness scheduling problem with stepdeteriorating jobs. The objective of this problem is to determine the sequence policy of the jobs under consideration so as to minimize the total tardiness. In order to solve the problem, an integer programming model was developed. Solutions were obtained by the CPLEX up to the size of 25 jobs. Due to the intractability of the problem, it is impossible to solve large instances to optimality by using the CPLEX. Based on two properties introduced, a heuristic called the SWSP and the GVNS algorithm were proposed to obtain near-optimal solutions of the problem. Meanwhile, the standard VNS algorithm was presented as reference for evaluating the performance of our algorithms. Experimental results showed that the proposed GVNS outperforms other heuristics in terms of relative percentage deviation from the best solution for both smalland large-sized instances. The computational times are in general short for the GVNS and the VNS. The GVNS performed better than the VNS in both cases. Furthermore, the GVNS is more robust than the VNS with regards to the choice of different intervals of deteriorating dates. For further study, a suggestion is to consider the setup times between different jobs for the single-machine scheduling problem with step-deteriorating jobs. In addition, multi-machine or multi-criteria cases with step-deterioration are also encouraged.