YUE Fan,SONG Shiji,JIA Peng,WU Guangping,and ZHAO Han

1.Department of Automation,Tsinghua University,Beijing 100084,China;2.Research Institute of Systems Engineering,Beijing 100072,China;3.Jiuquan Satellite Launch Centre,Jiuquan 735000,China;4.Department of Basic Science,Army Logistics Academy,Chongqing 401311,China

Abstract:The single machine scheduling problem which involves uncertain job due dates is one of the most important issues in the real make-to-order environment. To deal with the uncertainty, this paper establishes a robust optimization model by minimizing the maximum tardiness in the worst case scenario over all jobs.Unlike the traditional stochastic programming model which requires exact distributions, our model only needs the information of due date intervals. The worst case scenario for a given sequence that belongs to a set containing only n scenarios is proved,where n is the number of jobs.Then,the model is simplified and reformulated as an equivalent mixed 0-1 integer linear programming (MILP) problem. To solve the MILP problems efficiently, a heuristic approach is proposed based on a robust dominance rule. The experimental results show that the proposed method has the advantages of robustness and high calculating efficiency, and it is feasible for large-scale problems.

Keywords: robust optimization,single machine scheduling,maximum tardiness,uncertain due date.

1.Introduction

Since the job due dates are usually involved in the measure of performance,the single machine scheduling problem which involves uncertain job due dates is one of the most important issues in the real make-to-order environment. This problem is motivated by a real-world application in an iron-making process.Specifically,the blast furnace produces molten iron and taps it into a set of torpedo cars (TPCs) once in a while. Then they are transported to the pretreatment station for dephosphorization,desulfurization, etc. The pretreatment station and TPCs are regarded as a single machine and jobs in this problem,respectively. From blast furnaces to the pretreatment station,TPCs are transported via rail.Since the speed and the route are difficult to be precisely determined,it is hard to predict the release times accurately when TPCs reach the pretreatment station.In general,the job due date is always closely related to the job release time.

The most common relationship between the due date and release time [1] is that dj= rj+pj+c, where rj,djand pjrepresent the release time, the due date and the processing time of the job j,respectively,c denotes a common slack variable.For the iron-making,once the TPC is not processed completely before the due date,it will have a tardiness.This means the temperature of the molten iron starts to drop. We need to reheat the molten iron in the converter if the temperature drops too much[2,3].Thus,it is necessary to find a robust schedule by minimizing the maximum tardiness when the worst case scenario occurs.

In order to hedge the uncertainty in the production environment, the influence of uncertain parameters is considered in the single machine scheduling problem [4–6]. Fuzzy optimization, stochastic programming and robust optimization are three main approaches to obtain robust schedules [6–8]. However, the stochastic optimization requires exact distribution of uncertain parameters,which limits its scope of applications [9,10]. Robust optimization is another way to optimize the system performance in uncertain environments [11–13]. The advantage of robust optimization is that it only requires a small amount of information for uncertain parameters,which are readily available in production. It was first introduced in 1995 in [14] and its aim was to minimize the maximum costs in the worst case scenario.There are three common ways to representing uncertain parameters: the interval data,a discrete set of scenarios,and the random variables[15,16]. In [17], the robust optimization model was established, in which the processing times of jobs were un-certain, and could only take values from closed intervals.The scenarios-based robust optimization models were presented in [18,19].They considered the uncertain processing times as a discrete set of scenarios. Chang et al. [5]proposed a novel distributional robust optimization model for a single machine scheduling problem with the uncertain job processing times.In this model,the uncertain processing time is regarded as a random variable. They use the mean and variance of the processing times, and develop three polynomial-timealgorithms to solve the model.The objective function of the robust optimization has two main forms: absolute robustness and relative robustness.There are two common categories of absolute robust criteria:minimizing the maximum absolute cost[20]and minimizing the maximum regret[21],which is a deviation from the optimal value. Relative robustness is to minimize the maximum relative deviation from the optimal value[22].

The minimizing maximum tardiness problem with unequal release time(1|rj|Tmax)is a non-deterministic polynomial hard problem (NP-hard) without considering the setting of the due date [23]. Some researchers proposed many algorithms early based on branch and bound which were able to give exact solutions to small-scale problems,but for large-scale problems, they were computationally expensive[24–26].In recent years,many heuristics were proposed to solve large-scale problems, such as the selfadaptive differential evolution heuristic approach[27],the genetic algorithm [28], and the approximate algorithm[29]. In general, if the deterministic problem is NP-hard,the corresponding robust optimization problem is also NP-hard[30]in the presence of uncertainty.For solving a robust single machine scheduling problem,some researchers firstly identified the worst case scenario,then developed efficient algorithms from the methods for deterministic problems in[31].Other researchers proposed fast and effective algorithms such as meta-heuristics and genetic algorithms in[32,33].These algorithms can solve robust solutions in a short time,especially for large-scale problems.

In this paper, we establish a min-max robust optimization model with the uncertain job due date,which only requires time intervals of job due dates and can be easily estimated from historical data in practice. We prove that the worst case scenario for a given sequence belongs to a set containing only n scenarios,where n is the number of jobs.The robust model is simplified and reformulated as an equivalent mixed 0-1 integer linear programming(MILP)problem. Then, an efficient heuristic approach robust dominance heuristic (RDH), based on the proven robust dominance rule is proposed to solve the MILP problem.The experimental results show that the proposed method has the advantages of robustness and high calculating efficiency,and it is feasible for large-scale problems.

The rest of this paper is organized as follows.Section 2 provides the establishment of the robust optimization model.And a property is derived to reformulate the model.In Section 3,we propose the RDH to solve the model efficiently.Numerical experiments and results are detailed in Section 4.Section 5 gives a conclusion for the paper.

2.Problem formulation and reformulation

We suppose that n jobs in set J = {1,2,...,n} are processed on a single machine.It assumes that there is no machine breakdown or preemption.

A feasible schedule is a job sequence,which can be represented by a matrix x=[xij](?i,j =1,2,...,n),where

It is obvious that the set of feasible schedules can be given as follows:

The maximum tardiness of a given sequence is thatwhere[y]+=max{y,0},djis the due date of job j,Ciis the completion time of job at the ith position, we assume C0= 0. Let f(x) = Tmax, then the deterministic single machine scheduling problem aims to find a schedule in X with the minimal f(x),therefore,the problem can be formulated as the following mathematical model, we denote it as D-TP:

where pjis the processing time of job j, rjis the release time of job j.Equation(1)is the objective function.Equations(2)and(3)state that f(x)is no less than the tardiness of every job in the sequence x and is no less than zero.Equation (4) ensures that only one job is processed at a time.Equation(5)ensures that the machine only processes one job at a time.Equation(6)means that the lead time of each job is equal to their processing time plus a common slack variable. Equation (7) means that the best schedule belongs to a set of feasible schedules.

In actual production,it is impossible for us to accurately predict the due date of every job due to various disturbances,but a due date interval can be easily estimated.We assume that the due date interval of job j isThat means job j can be released at any value in its due date interval. Then for a given sequence,there are infinite due date scenarios which can be represented as a set S, i.e.,and d ∈ S.Siis the start processing time of job at the ith position.

For a schedule x ∈X,we define Tmaxin the worst case scenario(W-Tmax)as follows:

where vector d ∈S represents a possible due date scenario,and Tiis the tardiness of job in the ith position.

Therefore,this robust single machine scheduling problem(R-TP)takes the following form:

In the R-TP, since the number of due date scenarios is infinite, the exhaustive attack method is unsuitable for solving this problem.We focus on the identification of the worst case scenario firstly.

Let Tmax(x) = max Ci? d[i]+,?i = 1,2,...,n represent the maximum tardiness for a given sequence x,where [i] denotes the job index in the ith position of x.We denote the Tmax(x) under the worst case scenario as W-Tmax(x), and suppose that the job in the k-position is the one with the maximum tardiness in a given sequence x,then

where[k]represents the job index in the kth position.

Because dj= rj+pj+c,and the value of Ckis positively related to the value of r[k],then Ckis also positively related to the value of d[k].

We cannot intuitively identify the worst case scenario of due dates.That is because in(10),if we want to maximize the value of Ck, we need to set the due dates of all jobs before the (k +1)th position take the maximums in their due date intervals, i.e.,If we want to minimize d[k], we need to set the due date of the job in the kth position that takes the minimum in its due date interval, i.e., d[k]= d[k]. Therefore, we cannot directly identify the worst case scenario which can make[Ck?d[k]]+reach the maximum.

It derives that the number of the worst possible scenarios is finite for any sequence.

Property 1 For any sequence x, the worst case scenario belongs to a finite scenario set U = {di,i =1,2,...,n}.In di,the due date of job j is defined as follows:

Proof We define a set of worst possible scenarios S(x)=arg maxd∈Sf(x,d)for schedule x.Suppose that the worst case scenario d?∈S(x)and d?/∈U.If the job in the kth position is the one with the W-Tmax(x) in the sequence x,the due date of this job can be denoted as d[k].Then we can always increase the W-Tmax(x)by adjusting d[k]toand postpone the due dates of the rest jobs to?l=1,2,...,k ? 1,k+1,...,n.

Through the above adjustment,the worst case scenario d?∈ U.

Using Property 1,we can reduce the set of the worst possible scenarios from S to U with only n elements.Then we can further reformulate the R-TP into the following MILP:

3.RDH

To illustrate the process of the RDH, we first introduce some definitions. In a given sequence x, if there are two jobs i and j, such that rirj, pipjand job j precedes job i, we define these two jobs as an incompatible pair(j,i).Let B be the set of jobs from job j to i in the sequence x.Then we denote the set of jobs before and after B as B?and B+,respectively.

We suppose that there is an incompatible pair (j,i) in x. We want to prove that if we move job i to the front of job j,the Tmaxof x will be reduced or remain unchanged.Firstly, we define that job i dominates job j, if both inequalitieshold.Then we will prove this conclusion in the following dominance rule.

Theorem 1 Dominance rule. Consider two arbitrary jobs i and j in a given sequence x for the D-TP. If both inequalitieshold, the sequence in which i precedes j has no larger W-Tmaxvalue than other sequences.

Proof We design the following adjustment method to make sure that the Tmaxof the new sequence will be no larger than Tmaxof the unadjusted sequence:

(i) For a sequence x, we search from scratch and find the last incompatible pair(j,i).The process start time and the completion time of the jobs in set B are expressed as SBand CB.

(ii)We divide set B into four subsets: {i},{j},K,and L,where K = {k : rkSB+pi, k ∈ B {i,j}},and L = {k : rk>SB+pi, k ∈B {i,j}}.Let |K| = h1,|L| = h2,[l]Kand[l]Lbe the index of job in the lth position in sets K and L,respectively.

(iii) Let i →j represent the job i precedes job j, then we schedule a new sequenceL → B+and keep the order of jobs in set B?,K, L and B+unchanged.Fig.1 shows the order of jobs in x.

Fig.1 Illustration for the dominance rule

The second items in Tmax(x)andcan be expanded as follows:

In x, Ti(x) = CB? pi? ri, Tj(x) = SB? rj.Since CB? piSBand rirj, Ti(x) Tj(x), i.e.,max{Ti(x),Tj(x)} Ti(x).

When we adjust the sequence from x tothe earliest process start time of jobs in set K is advanced from.Thus,the tardiness of these jobs will not be increased,i.e.,

Based on the above analysis,we know that

If there is a sequence that does not satisfy the dominance rule,we can always decrease Tmaxby adjusting job positions based on the above adjustment method. Next, we derive the robust dominance rule based on Theorem 1 for the problem R-TP.

Theorem 2 Robust dominance rule.Consider two arbitrary jobs i and j in a given sequence x for the MILP.If both inequalitiesandhold,the sequence in which job i precedes job j has no larger W-Tmaxvalue than other sequences.In this case,job i dominates job j.

Based on the Theorem 2,we propose our RDH to generate the optimal solution, and the steps of the RDH are listed as follows:

Step 1 Generate initial solution

Step 2 Search for the incompatible pair(j,i).

After k(k =1,2,...)iterations,we search from scratch and find the last incompatible pair(j,i)in xk.If it can be found,go to Step 3,otherwise exit the process and use the current solution as the optimal solution generated by the RDH.

Step 3 Group jobs in xk.

Divide unscheduled jobs into four blocks: {i},{j},K,and L,where K =SB+pi, k ∈ B {i,j}},and L=SB+pi,k ∈B{i,j}}.

Step 4 Adjust the order of the blocks and generate a new sequence.

Reschedule a new sequencex(k+1)asB? → i →K → j → L → B+and keep the order of jobs in setB?,K,LandB+unchanged.Then returnx(k+1)to Step 2.

Theoretically, since the RDH is an algorithm based on Theorem 2, the optimization goal can be continuously improved when searching for a better solution by the RDH.The algorithm can converge quickly and be efficient.Therefore, for the large-scale problem that the business solver IBM ILOG CPLEX optimization studio (CPLEX)based on the branch and bound method cannot solve in an acceptable time,the RDH can quickly find the suboptimal solution.

The pseudocode of the RDH is given in Algorithm 1.In the RDH, we useito denote the position of the job in ?xwith W-Tmax,and usePto keep candidate sequences.

Algorithm 1

Find the job with W-Tmaxin ?xand denote its position asi,its job index is[i].

whilei>1 do

Calculate the set of position indicesφi={j:ji,job[j]dominates job[i]}.

ifφi= ? then

else

k=1.

end if

whilei>kdo

i=i ?1.

end while

Leti=k.

end while

Select the best sequence with the minimal W-TmaxfromPand denote it asx?.

Ensure:x?andfR(x?).

Then, we use the following example with six jobs to describe the process of the RDH. We suppose that in the initial sequence,the order of jobs is that 1→2→3→4→5→6,wherei →jrepresents jobiprecedes jobj,while job 6 is the one with W-Tmax, job 1 dominates job 3,and job 3 dominates job 6.

Fig.2 shows the search process of the RDHon the initial sequence 1→2→3→4→5→6.xkexpresses the sequence obtained by thekth search forx?.

Fig. 2 Search process of RDH based on the initial sequence 1 →2 → 3 → 4 → 5 → 6

Since job 5 does not dominate job 6,we exchange them to generate a new sequencex1. Since job 4 also does not dominate job 6,we exchange them and generatex2.Job 6 cannot be exchanged with job 3 in terms of the fact job 3 dominates job 6,then we fix the position of job 6 and move job 3 forward.We exchange job 3 and job 2 to generatex3.Job 3 cannot be moved forward since job 1 dominates it.Then the RDH terminates. The sequence obtained by the RDH is the robust sequence of the R-TP.

4.Numerical experiments and results

In this section, we implement computational experiments to validate the effectiveness and efficiency of the proposed robust model and the proposed RDH based on artificial instances.The experiments are coded in MATLAB and executed on a personal computer with an Intel Core CPU i7-4710K and 8 GB RAM at 2.5 GHz.Experimental settings and results are presented and analyzed in following subsections.

4.1 Experimental setup

4.1.1 Experimental design

(i)Experiments on the performance of RDH

The common method to solve the single machine scheduling problem is to give an initial solution, improve the initial solution by local search(LS),variable neighborhood search(VNS),etc.[34–37],and get a better solution.Among these methods,the VNS has the best performance.The main idea of the VNS is to explore multiple neighborhood structures systematically instead of a single neighborhood, and escape local minima (in the case of minimization).

To evaluate the effectiveness and efficiency of the proposed RDH,we compare it with the VNS[34]and CPLEX solver which calculates the MILP with the branch and bound method.

(ii)Experiments on robustness of solutions

To evaluate the robustness of solutions,we compare robust sequencex?obtained by the RDH with nominal sequenceobtained by FCFS according toThe average W-Tmaxofandx?can be calculated.Let DEV express the average reduction of the averageW-Tmaxby adopting robust sequencex?instead of nominal sequencewhile R-DEV is the ratio of the reduction.These four evaluation indexes are defined as follows:

From these definitions, it is clear that the DEV expresses the average reduction of the average W-Tmaxby adopting robust sequencex?instead of nominal sequencewhile R-DEV is the ratio of the reduction.(iii) Robustness of solutions under different distributions of uncertain release times Note that the proposed approach does not require any information about the distribution of release times, and these experiments are designed to evaluate the robustness of solutions under some other distribution assumptions. Four classical probability distributions are examined, namely the uniform distribution, normal distribution, poisson distribution, and geometric distribution.We use 10 jobs to evaluate the robustness of the solutions. Average values of release time intervals are set aswhile processing time values are set aspj= [12,15,10,11,14,10,15,

12,14,11], and the deviation values are set asδj=[11,13,8,5,15,9,12,15,5,13]. Letc= 5. By setting these parameters,the release time intervals are fixed.Five thousand scenarios are generated based on the intervalsaccording to the four distributions above. The mean objective values and standard deviations ofandx?are calculated, denoted asW-Tmax(x?), SD()and SD(x?), respectively. We introduce another two indexes to evaluate the robustness ofx?,defined as follows:

4.1.2 Data setting

The above first two types of experiments are based on the same data setting, i.e., our experiments consist of the test problems involving|J|={10,50,100,150,200,250,300,400,500}jobs and a single machine. Letbe the average value of the release time interval,andδjrepresents the maximum deviation of the release times fromthat is,Then the real release time interval of jobjcan be denoted as

To simulate a real single machine scheduling problem of steel-making, we set the average release time gap of two adjacent jobs to be the same with the average processing time.The processing timepjis a random integer generated from a uniform distribution of interval[8,12].Thenis randomly generated from a uniform distribution of interval [0,10(n ?1)].δjis randomly generated from a uniform distribution of interval[5,15].We setc= 3,5,8(representing short,medium,long lead time,respectively),thendj=rj+pj+c.

The maximum computation time for all experiments is set to 1 800 s(30 min).For problem instances of different sizes,100 trials are conducted to obtain the average results.

4.2 Experimental results and analysis

4.2.1 Evaluation of the proposed RDH

To evaluate the performance of the proposed RDH, we compare it with variable neighborhood search(VNS) and CPLEX in terms of average W-Tmax, gap (%) which means the relative deviations of obtained solutions from the best solutions,and CPU time consumption in Table 1.The optimal solutions (or the best solution found by the methods)are highlighted in bold.

Table 1 Computational results of RDH,VNS and CPLEX

Continued

The results show that for anyc, CPLEX can calculate the accurate optimal solutions for small-scale problems(n150)within an acceptable period,while the RDH and VNS can give approximate solutions with small gaps. In these cases,CPLEX outperforms RDH in terms of the solution quality.However,for large-scale problem instances(150

As for the computational efficiency of the tested solution approaches,in Table 1,we can see that the RDH and VNS require fewer computation efforts to calculate solutions.For the same scale problems,RDH has a shorter computation time than VNS.The computation time of CPLEX increases dramatically as problem scales increase.In particular,its time consumption exceeds 1 800 s whenn150.In view of this, the RDH holds absolute advantages compared to CPLEX and VNS.

4.2.2 Evaluation of the robust model and solutions

To compare the performance of robust sequencex?with nominal sequencewhich is obtained by FCFS based on.Results in Table 2 show that robust sequencex?,which is theoretically the optimal solution of the MILP problem,certainly possesses smaller W-Tmax. These results verify the robustness of the proposed robust scheduling model.

Table 2 Comparison between robust and nominal sequences

Table 3 shows the robustness of solutions with different probability distributions.It reveals that solutions based on different distributions are generally in line with solutions based on average values and deviations obtained by the RDH,even though higher mean values and lower standard deviations of the AWT can be obtained by using robust sequencex?. This means the performance of the robust sequence is more stable. The results here again confirm the distributional robustness of the proposed heuristic approach.

Table 3 Simulation results with different distributions(n=10)

5.Conclusions

The problem with uncertain job due dates is one of the most important issues in single machine scheduling,which is a realistic difficult problem in iron-making.To deal with the uncertainty,we establish a robust model by minimizing the maximum tardiness in the worst case scenario over all jobs.Unlike the traditional stochastic programming model which requires exact distributions, our model needs only the information of due date intervals. The optimization objective is designed as minimizing the maximum tardiness in the worst case scenario.Based on the proven property in which the number of the worst possible scenarios is finite for any sequence, the R-TP can be reformed as a solvable MILP. Then, the RDH is proposed to solve the linear programming problem.In this algorithm,we firstly generate an initial sequence by the FCFS rule based on the average release times of all jobs.Then,the initial sequence is further improved by a variable neighborhood search according to the robust dominance rule.With three extensive experiments,we validate and verify that the RDH has clear advantages over CPLEX in terms of solution quality, robustness and computation time especially for large-scale problems.

There are many possible directions in future research.For example,we can consider other uncertain parameters,such as the processing time and delivery time.We can also design more efficient meta-heuristics to improve the effectiveness and efficiency of the proposed problem.