Scheduling in a two-stage supply chain with uncertain parameters

1 Introduction

Production and transportation operations are two key segments in supply chain. To satisfy the customers’ demands, manufacturers endeavor to integrate these two segments and scheduling jointly in a coordinated manner. In the past decade, a lot of research has been done on various supply chain scheduling. Many different intelligent algorithms were introduced to solve supply chain scheduling problems. The first study of supply chain scheduling by Hall and Potts [11] considered a variety of scheduling, batching and delivery problems in an arborescent supply chain. They proposed models which minimize the total scheduling and batch delivery cost. Finally, they demonstrated that cooperation between a supplier and a manufacturer may reduce the total system cost depending upon the scheduling objective. Chang and Lee [2] studied a class of two-stage scheduling problems in which the first stage is job production and the second stage is job delivery. The worst-case performance ratio for their heuristic is proven to be 5/3 and the bound is tight. Dawande et al. [7] studied conflict and cooperation issues arising in a supply chain where a manufacturer made products which are shipped to customers by a distributor. They considered several practical scenarios about the level of cooperation between the manufacturer and the distributor and provided algorithms for each of these scheduling problems. At the last, they demonstrated that the cost saving provided by cooperation between the decision makers is usually significant. Chen [3] investigated the integrating production and transportation scheduling problem in a make-to-order environment. Results from simulation showed that their heuristics outperformed the human expert schedule. You and Hsieh [28] proposed a mixed integer nonlinear programming model to deal with an assembly and transportation scheduling problem. They developed a heuristic algorithm to minimize overall costs. The presented heuristic algorithm is shown to perform well compared with the well-known GAMS/BARON and Lingo commercial software on testing randomly-generated problems. Zegordi et al. [30] considered the scheduling of products and vehicles in two-stage supply chain environment. They set a model and proposed a gendered genetic algorithm to solve the model. Their experimental results showed that the proposed gendered genetic algorithm outperformed standard genetic algorithms. Ullrich [26] studied machine scheduling and vehicle routing with time windows integrated problem. His numerical experimental results show that the genetic algorithm outperforms the classic decomposition approaches in some cases. Cakici et al. [1] studied the problem of minimizing the total weighted tardiness and total distribution costs in an integrated problem and distribution environment. They developed different genetic algorithms and computational tests had shown that both the mathematical modeling and heuristics approaches were competitive and can produce a significant number of non-dominated solutions. Gordon et al. [9] reviewed a lot of results on scheduling with due date assignment under such conditions as given precedence constraints or various scenarios of processing time deterioration and learning. Mazdeh et al. [21] described a model for minimizing total weighted flow times plus delivery cost for a set of jobs that are to be processed on a single machine and delivered in batches to a single customer. They proved that when the maximum number of batches is fixed, the branch-and-bound scheme outperforms the dynamic programming algorithm. Pei et al. [22] modeled a supply chain scheduling problem considering multiple manufacturers and non-identical job sizes in a two-stage supply chain to minimize the makespan. Numerical results showed that the proposed modified gravitational search algorithm better than PSO and GA. Hamidinia et al. [12] proposed an effectivegenetic algorithm for minimizing total tardiness and earliness of weighted jobs in a batched delivery system. The proposed genetic algorithm outperformed thetraditional genetic algorithm under non-batched systems. Lei et al. [14] proposed a two-phase solution approach to deal with an integrated production, inventory and distribution routing problem. They provided a two-phase methodology to effectively solve this problem. Phase I solved a mixed integer programming model and phase II solved an associated consolidation problem.

Almost all of the literature mentioned above assume that parameters are known and precisely defined. In other words, they are deterministic. However, in real supply chain scheduling systems, we will encounter various uncertain factors. So many scholars made great efforts on stochastic case. But in many cases, using probability theory to deal with an indeterminate variable is not appropriate [19]. Thus uncertainty theory, based on the experts’ judgement, was proposed by Liu [15]. From then on, uncertainty theory has become a powerful mathematical tool to deal with various problems under incomplete information. For example, Chen [4] derived an American option pricing formula from uncertain financial market [16, 20] and some mathematical properties of these formulas were studied. Ding and Zhu [8] proposed two types of uncertain models for project scheduling problems. One type of models can be transformed to their crisp forms, which may be solved by classical optimization methods. The other type of models could not be transformed to their crisp forms, an uncertain simulation is employed to approximate uncertain functions. Ke [13] studied an uncertain random time-cost trade-off model. He showed that the proposed model can be solved by uncertain random simulation and GA. Sheng and Gao [24] investigated chance distribution of the maximum flow of an uncertain random network. A new method is derived to calculate chance distribution of the maximum flow for an uncertain random network. Yao and Chen [27] proposed a method for solving uncertain differential equations. They designed a numerical method for solving uncertain differential equations, which essentially solves each α-path and produces an inverse uncertainty distribution of the solution. Zhu [31] studied an optimal control model for an uncertain differential equation. The principle of optimality for uncertain optimal control is obtained by applying Bellman’s Principle of Optimality. And a fundamental result called the equation of optimality in uncertain optimal control is given. Finally, a portfolio selection model can be solved by the equation of optimality effectively.

When no samples are available to estimate a probability distribution, invite some domain experts to evaluate the belief degree that each event will happen is necessary. In the process of production and processing, many procedures do not require high accuracy. For example, rough turning, wood processing, heavy boring and so on. Because there exist many man-made factors during processing, a worker must estimated the processing time of job based on his past experience and does not estimate the processing time according to a probability distribution. Due to the variability of road conditions and weather, a driver can not estimate transportation time according to a probability distribution. But he can estimate transportation time by his past experience. So it is reasonable that the processing time and transportation time are assumed to be uncertain variables. Moreover, it is inappropriate to model belief degrees by probability theory because it may lead to counterintuitive results. Consider a counterexample presented by Liu [19]. Assume there is one truck and 50 bridges in an experiment. Also assume the weight of the truck is 90 tons and the 50 bridge strengths are independent and identically distributed uniform random variables on the interval [95, 110] in tons. For simplicity, suppose a bridge collapses whenever its real strength is less than the weight of the truck. Now let us have the truck cross over the 50 bridges one by one. It is easy to verify that Pr { “ Truck can cross over the 50 bridges ″ } = 1 . That is to say, we are 100% sure that the truck can cross over the 50 bridges successfully. However, when there do not exist any observed samples for the bridge strength at the moment, we have to invite some bridge engineers to evaluate the belief degrees about it. As we stated before, human beings usually estimate a much wider range of values than the bridge strength actually takes because of the conservatism. Assume the belief degree function is Φ ( x ) = { 0 , if x ≤ 80 x - 80 40 , if 80 ≤ x ≤ 120 1 , if x ≥ 120 . Let us imagine what will happen if the belief degree function is treated as a probability distribution. At first, we have to regard the 50 bridge strengths as independent and identically distributed uniform random variables on [80, 120] in tons. If we have the truck cross over the 50 bridges one by one, then we immediately have Pr { “ the truck cross over the 50 bridges ″ } = 0 . 75 50 ≈ 0 . Thus it is almost impossible that the truck crosses over the 50 bridges successfully. Unfortunately, the two results are at opposite poles. This example shows that, by inappropriately using probability theory, a sure event becomes an impossible one. The error seems intolerable for us. Hence the belief degrees cannot be treated as subjective probability.

In this paper, we will study a two-stage supply chain scheduling problem. In the first stage there are some parallel machines for processing jobs, and in the second stage there are some vehicles for carrying completed jobs to customers. We assume that the processing time and transportation time are uncertain variables.

The remaining of this paper is organized as follows. In the next section, some basic definitions and properties about uncertainty theory will be introduced. In the third section, a mathematical model will be established about the two-stage supply chain scheduling problem with uncertain processing time and transportation time. In the fourth section, a solution methodology of the model will be proposed. In the fifth section, an numerical example is provided.

2 Preliminary

This section will introduce some concepts in uncertainty theory. Let Γ be a nonempty set, and L be a σ-algebra over Γ. Each element A ∈ L is called an event. A set function M from L to [0, 1] is called an uncertain measure if it satisfies three axioms 15: (normality) M { Γ } = 1; (duality) M { Λ } + M { Λ c } = 1 for any event Λ; and (subadditivity) M { ⋃ i = 1 ∞ Λ i } ≤ ∑ i = 1 ∞ M { Λ i } for every countable sequence of events Λ₁, Λ₂, ⋯. The triplet ( Γ , L , M ) is called an uncertainty space. To obtain an uncertain measure of a compound event, a product uncertain measure was defined by Liu [16]: (Product Axiom) Let ( Γ k , L k , M k ) be uncertainty spaces for k = 1, 2, ⋯, the product uncertain measure M is an uncertain measure satisfying M { Π i = 1 ∞ Λ k } = ∧ i = 1 ∞ M k { Λ k }, where Λ _k are arbitrarily chosen events from L k for k = 1, 2, ⋯, respectively.

An uncertain variable is a function ξ from an uncertainty space ( Γ , L , M ) to the set R of real numbers such that for any Borel set B of real numbers, the set {ξ ∈ B} = {γ ∈ Γ : ξ (γ) ∈ B} is an event. An uncertain vector ξ is (ξ₁, ξ₂, ⋯ , ξ _n ) where ξ₁, ξ₂, ⋯ , ξ _n are uncertain variables. The uncertain distribution Φ of an uncertain variable ξ is defined by Φ ( x ) = M { ξ ≤ x } for any real number x.

Definition 1. [16] The uncertain variables ξ₁, ξ₂, ⋯, ξ _n are said to be independent if M { ⋂ i = 1 m ( ξ i ∈ B i ) } = min 1 ≤ i ≤ m M { ξ i ∈ B i } for any Borel sets B₁, B₂, ⋯ , B _m of real numbers.

Example 1. Let ξ₁ be an uncertain variable and let ξ₂ be a constant c. For any Borel sets B₁ and B₂, if c ∈ B₂, then M { ξ 2 ∈ B 2 } = 1 and M { ( ξ 1 ∈ B 1 ) ∩ ( ξ 2 ∈ B 2 ) } = M { ξ 1 ∈ B 1 } = M { ξ 1 ∈ B 1 } ∧ M { ξ 2 ∈ B 2 } . If c ∉ B₂, then M { ξ 2 ∈ B 2 } = 0 and M { ( ξ 1 ∈ B 1 ) ∩ ( ξ 2 ∈ B 2 ) } = M { ∅ } = 0 = M { ξ 1 ∈ B 1 } ∧ M { ξ 2 ∈ B 2 } .

Given an increasing function Φ (x), Peng and Iwamura [23] introduced an uncertainty space ( R , B , M ) as follows. Let B be the Borel algebra over R. Let C be the collection of all intervals of the form (- ∞ , a], (b, + ∞), ∅ and R. The uncertain measure M is provided in such a way: First, set M { ( - ∞ , a ] } = Φ ( a ) , M { ( b , + ∞ ) } = 1 - Φ ( b ) , M { ∅ } = 0 , M { R } = 1 . Second, for any B ∈ B, there exists a sequence {A _i } in C such that B ⊂ ⋃ i = 1 ∞ A i . Thus M { B } = { inf B ⊂ ∪ A i ∑ i = 1 ∞ M { A i } if inf B ⊂ ∪ A i ∑ i = 1 ∞ M { A i } < 0.5 , 1 - inf B c ⊂ ∪ A i ∑ i = 1 ∞ M { A i } if inf B c ⊂ ∪ A i ∑ i = 1 ∞ M { A i } < 0.5 , 0.5 otherwise .(1) The uncertain variable defined by ξ (γ) = γ from the uncertainty space ( R , B , M ) to R has the uncertainty distribution Φ. This kind of uncertain variables are called common. An uncertain variable with an uncertainty distribution will be regarded as a common one in order to deal with distributions of its functions.

Let us discuss the distribution of f (ξ) for a common uncertain variable ξ or a common uncertain vector. Assume C is the collection of all intervals of the form (- ∞ , a], (b, + ∞), ∅ and R. Each element A _i emerging in the sequel is in C.

Theorem 1. [32] (i) Let ξ be a common uncertain variable with the continuous distribution Φ (x) and f (x) a Borel function. Then the distribution of the uncertain variable f (ξ) is Ψ ( x ) = M { f ( ξ ) ≤ x } = { inf { f ( ξ ) ≤ x } ⊂ ∪ A i ∑ i = 1 ∞ M { A i } if inf { f ( ξ ) ≤ x } ⊂ ∪ A i ∑ i = 1 ∞ M { A i } < 0.5 , 1 - inf { f ( ξ ) > x } ⊂ ∪ A i ∑ i = 1 ∞ M { A i } if inf { f ( ξ ) > x } ⊂ ∪ A i ∑ i = 1 ∞ M { A i } < 0.5 , 0.5 otherwise .

(ii) Let f : R n → R be a Borel function, and ξ  = (ξ₁, ξ₂, ⋯ , ξ _n ) be a common uncertain vector. Then the distribution of the uncertain variable f (ξ) is Ψ ( x ) = M { f ( ξ 1 , ξ 2 , ⋯ , ξ n ) ≤ x } = M { ( ξ 1 , ξ 2 , ⋯ , ξ n ) ∈ f - 1 ( - ∞ , x ) } = { sup Λ 1 × Λ 2 × ⋯ × Λ n ⊂ Λ min 1 ≤ k ≤ n M k { Λ k } if sup Λ 1 × Λ 2 × ⋯ × Λ n ⊂ Λ min 1 ≤ k ≤ n M k { Λ k } > 0.5 , 1 - sup Λ 1 × Λ 2 × ⋯ × Λ n ⊂ Λ c min 1 ≤ k ≤ n M k { Λ k } if sup Λ 1 × Λ 2 × ⋯ × Λ n ⊂ Λ c min 1 ≤ k ≤ n M k { Λ k } > 0.5 , 0.5 otherwise where Λ = f^-1 (- ∞ , x), and each M k { Λ k } is derived from (1).

Definition 2. [15] Let ξ be an uncertain variable, and α ∈ (0, 1]. Then ξ sup ( α ) = sup { r ∣ M { ξ ≥ r } ≥ α } is called the α-optimistic value to ξ, and ξ inf ( α ) = inf { r ∣ M { ξ ≤ r } ≥ α } is called the α-pessimistic value to ξ.

The uncertain variable ξ will reach upwards of the α-optimistic value ξ_sup with uncertain measure α, and will be below the α-pessimistic value ξ_inf (α) with uncertain measure α.

Founded by Liu [17] in 2009 and refined by Liu [18] in 2010, uncertain programming is a type of mathematical programming involving uncertain variables. Assume that x is a decision vector, ξ is an uncertain vector, f (x, ξ) is an objective function, and g _j  (x, ξ) are constraint functions, j = 1, 2, ⋯ , p. In order to obtain a decision with minimum objective value subject to chance constraints, Liu [17] proposed the following uncertain programming model, { min x min f ¯ f ¯ subject to : M { f ( x , ξ ) ≤ f ¯ } ≥ β M { g j ( x , ξ ) ≤ 0 } ≥ α j , j = 1 , 2 , ⋯ , p , where α _j and β are specified confidence levels for j = 1, 2, ⋯ , p, and min f ¯ is the β-pessimistic value. If we want to minimize the pessimistic value of an objective (such as makespan or cost) subject to some chance constraints, we can use the above model.

A key problem in the research area of uncertain programming is how to solve the above model. Under the framework of uncertainty theory, an uncertain programming model may be transformed into a crisp model if the related functions are monotone [18]. If the related functions are not monotone, we may employ uncertain simulation [32] to approximate the uncertain measure of an event and β-pessimistic value of an uncertain variable. Thus, we can solve the model by simplex method, cutting plane method, tabu search, interior point method, and heuristic algorithms.

We will introduce our problem in two stages. The first stage is processing jobs on machines in a manufactory. And the second stage is to delivery the completed jobs from the manufactory to customs in batch.

At first stage, most of the literature considered single machine. The interested readers may refer to [1, 29]. However, multiple machines are more closer to the reality. So we assume that there are n independent non-preemptive jobs to be sequenced for processing on m parallel identical machines and processing time of each job is assumed to be an uncertain variable. In addition, the transportation times are assumed to be uncertain variables. The following assumptions are also assumed in advance.

All of the jobs are available for processing at time zero and can be processed on any machine. And each machine can process one job at a time.

Our objective is to minimize total cost, which consists of inventory holding cost, delivery earliness penalty cost, transportation cost, and delivery tardiness penalty cost. The structure of the problem is shown in Fig. 1.

Our problem considered in this paper is a parallel-machine scheduling problem in which jobs are ultimately delivery to customers in batches. There are n jobs have to be processed by m machines. M vehicles transport jobs in b batches. Each job has a size v. The capacity of each vehicle is Q. Each job has a processing time ξ _i for i = 1, 2, ⋯ , n and each trip has a transportation time η _r . As previously mentioned, we assume that the parameters ξ _i and η _r are uncertain variables. They have uncertainty distribution Φ _{ξ i} and Φ _{η r} , respectively. The cost delivery each batch to customer is denoted by R. The completion time of job i is denoted by C _i . Completion time of a job in kth location on machine j is denoted by T _jk . C_max denotes the makespan. The actual departure time of job i is denoted by G _i . Due date of batch r is denoted by d _r . The job is on-time and thus no penalty is incurred if the job is delivered to the customer at the due date, while it will incur an earliness penalty if the job is arrived before the due date and a tardiness penalty if the job is arrived after the due date. The job will incur a holding penalty if the job is not delivered to the customer at its completion time. Moreover, associate with each job is a unit earliness penalty cost α _r , a unit tardiness penalty cost β _r and a unit inventory holding cost h _i . In addition, decision variables is denoted by x _ijk .

For the sake of convenience, we list the symbols used in the sequel as follows:

Decision variables:

Objective function:

f = ∑ i = 1 n ∑ r = 1 b δ ir · max { 0 , G i - C i } · h i + ∑ r = 1 b max { 0 , max i ∈ B r { G i } + η r - d r } · β r + ∑ r = 1 b max { 0 , d r - ( max i ∈ B r { G i } + η r ) } · α r + R · ( ∑ r = 1 b ⌈ v · δ ir Q ⌉ )(2) where B _r indicates the set of indices of jobs which belong to batch r, δ _ir  = 1 if job i is assigned to batch r and 0 otherwise.

In the right hand side of (2), the first term indicates holding cost; the second term indicates delivery tardiness penalty cost; the third term indicates delivery earliness penalty cost; and the fourth term indicates transportation cost.

The completion time T _jk of the kth job on machine j: T j 1 = ∑ i = 1 n x ij 1 ξ i when job k is assigned at the first position; otherwise, T jk = ∑ i = 1 n x ijk ξ i + T j , k - 1, k = 2, 3, ⋯ , n.

The completion time C _i of job i: C _i  = T _jk where j and k are unique indices that satisfy x _ijk  = 1.

The earliest time of batch r being sent is equal to the completion time of the last job belonging to batch r if there are enough vehicles. Otherwise, jobs must wait for a vehicle back. That is, G i = max i ′ ∈ B r { C i ′ } if there are enough vehicles and G i = max i ′ ∈ B r { C i ′ } + H · η r for some positive integer H otherwise.

Due to the size of jobs and the capacity of vehicles, when there are enough vehicles, we can carry away all jobs at once; otherwise, we can carry away a portion of jobs and wait for the vehicles coming back to carry the remaining portion. A job can be delivered when all other jobs of the batch are completed. Otherwise, the job must wait for other jobs. A job should be stored from accomplishment to delivery. The period will emerge storage costs. Arrival time of the entire batch is equal to arrival time of the last job of the batch. The symbol ⌈· ⌉ indicates rounding up operation.

4 Solution methodology

The complexity of production and transportation problem is NP-complete [10]. Due to the complexity, heuristic approaches are used to generate good approximate solutions. Genetic algorithm (GA) and simulated annealing algorithm (SA) have been shown to be promising technique by many researchers for solving supply chain scheduling problems. GA can reach the optimization solution quickly, but it is liable to be trapped in a local optima. SA can jump out of the local optimization and search for the best solution. The hybrid algorithm (GSA) combining GA and SA may also be an efficient method for optimization problems. In this paper, we use GSA to deal with the first stage of the problem. The second stage of the problem can be handled according to the method described in Section 3.2. Before using GSA, we must use the aforementioned uncertain simulation algorithm to approximate the objection function min f ¯. Most researchers are familiar with GA and SA. Now we list the steps of GSA.

Step 1. Solution representation

Solution is expressed as a string of real numbers. The length of a string of real numbers is equal to the number of jobs. The integral part of each real number is generated randomly from 1 to m (the number of machines). The fractional part is generated from the interval (0, 1), which is used to determine the sequence of jobs. An example of solution representation with 2 machines and 4 jobs is shown in Fig. 2. The decoding result of the solution is shown in Fig. 3. The integral part means the index of machine. If some real numbers have identical integral part, sort these real numbers in ascending order. The order of corresponding jobs is determined by sequence of real numbers. If two real numbers are equal, sequence of jobs will be randomly generated.

Step 2. Initialization

Generate chromosomes randomly in a range according to weight of batch.

Step 3. Calculate the target

Calculate objective function value of each chromosome by uncertain simulation.

Step 4. Selection procedure

Use roulette wheel method to select the parent chromosomes.

Step 5. Crossover and mutation

In the crossover operation, we select a position from parents’ genes randomly and then exchange all genes before this position. An example of crossover operation as shown in Fig. 4. In the mutation operation, we randomly select two position of a chromosome and exchange their values. An example of mutation operation as shown in Fig. 5.

Step 6. Update optimal value

We can obtain a neighbor of current optimal solution by choosing one position at random and alter its current value according to upper and lower bounds. Let Z_cur be optimal value of the objective function for the current solution and Z_new be the value of the objective function for the new candidate. If Z_new ≤ Z_cur, then new solution is accepted as the current solution. If Z_new > Z_cur, the new solution is accepted with probability Prob accept = exp ( - ( Z new - Z cur ) T t ), where t is the iteration number and T _t is the current temperature.

Step 7. Update temperature

GSA starts with an initial temperature T₀ = T _max . And the temperature decreases according to T _t  = T _max 0 . 95 ^t .

Step 8. Stopping rule

If the maximal generation is reached, then stop; otherwise, go to Step 3 to perform another iteration.

5 A Numerical example

Because our problem is different from the traditional deterministic problem, we can’t solve our problem with general heuristic method directly. First, we must put it into deterministic model. But the objective function of our problem is not a monotone function with respect to uncertain variables, so we can’t transform our model to a crisp form by the inverse uncertain distribution method in [18]. We may deal with our objective function by uncertain simulation [32]. Then our problem can be solved by heuristic algorithm after being transformed to a crisp form. In order to test the efficiency of heuristic method, we used GSA, HGA [28] and GA [26] deal with our problem, respectively.

Assume that there are 5 machines and 20 jobs (denote job i by J _i for i = 1, 2, ⋯ , 20). Consider six batches as shown in Table 1.

The weights of batches are from batch 1 to batch 6 in descending order. The time is specified in minutes. The cost is specified in yuan. Processing time of job i is assumed to be a linear uncertain variable: ξ i ∼ L ( i , i + 1 ) for i = 1, 2, ⋯ , 20, transportation time of a trip for jobs in the rth batch is assumed to be a linear uncertain variable: η r ∼ L ( r , ( r + 1 ) 2 ), the size of each job is 10, the number of vehicles is 3, the capacity of each vehicle is 10, a transportation cost of a vehicle for a trip is 100.

Due date d = [20, 25, 70, 100, 130, 200], inventory holding cost h = [20, 15, 20, 25, 30, 10, 35, 30, 25, 20,10, 30, 25, 35, 30, 20, 30, 15, 10, 20], delivery earliness penalty cost α _r  = [30, 20, 35, 40, 25, 30] , delivery tardiness penalty cost β _r  = [40, 30, 20, 25, 35, 15] , confidence levels τ = 0.8.

Population size of GSA is 30. Mutation rate is 0.25. Crossover rate is 0.75. T_max is 100. Maximum iteration number is 200.

The optimal objective value obtained by running GSA is 16562. The processing and transport time of the case are come from linear distribution of uncertain variable:

ξ = [1.06, 2.09, 3.27, 4.40, 5.48, 6.90, 7.60, 8.32,

  9.48, 10.58, 11.16, 12.83, 13.94, 14.60, 15.03,

  16.81, 17.61, 18.70, 19.10, 20.43],

η = [2.13, 3.16, 13.83, 21.61, 19, 47.13] .

Tables 2–4 list the optimal values of the decision variables x _ijk of Model 1 by using GSA.

First stage:

According to the optimal decision variables, we sort jobs on parallel machines. It is assumed that there is no time gap among these jobs on each machine. In other words, the optimal scheduling plan is

Machine 1: 6 → 8 →11 → 18

Machine 2: 4 → 10 → 13 → 17

Machine 3: 2 → 7 →12 → 20

Machine 4: 5 → 1 →14 → 16

Machine 5: 3 → 9 →15 → 19

Second stage:

Having completed scheduling at the first stage, jobs must be sent in batches. Assume that the earliest beginning processing time at the first stage is 8 a.m. The processing time of a job and transportation time are considered at confidence level τ = 0.8.

Batch 1 will be distributed at the first. Initially there are three idle vehicles. Inventory holding cost of {J₁, J₃, J₄} are higher than J₂, so we transport {J₁, J₃, J₄} preferentially. Because makespan of batch 1 is 6′31″, so the earliest transporting time of {J₁, J₃, J₄} is 8 : 06′31″. After three vehicles come back, use one vehicle to transport J₂. Departure time of J₂ is 8 : 14′15″. The arrival time of whole batch is 8 : 16′21″.

Makespan of batch 2 is 9′42″. That is to say that batch 2 can be delivered after 8 : 09′42″. However, the earliest vehicles are usable at 8 : 14′15″ and only have two. Thus the delivery time of {J₅, J₇} is 8 : 14′15″. The vehicle which transports J₂ in batch 1 will free at 8 : 18′27″. So J₆ can be delivered at 8 : 18′27″. The arrival time of whole batch is 8 : 21′36″.

Makespan of batch 3 is 26′24″. Batch 3 can be distributed after 8 : 26′24″. The earliest two idle vehicles which transport {J₅, J₇} at batch 2 are valid until 8 : 17′21″. The departure time of jobs {J₈, J₁₂} is 8 : 26′24″. The vehicle which transports J₆ in batch 2 will free at 8 : 21′33″. Departure time of job J₉ is 8 : 26′24″. The free time of the two vehicles which transport {J₈, J₁₂} is 8 : 40′02″. So the departure time of jobs J₁₀ and J₁₁ is 8 : 40′02″. The arrival time of whole batch is 8 : 53′50″.

Makespan of batch 4 is 29′. Batch 4 can be distributed after 8 : 29′. The earliest two vehicles are valid until 8 : 44′. Then the departure time of job J₁₄ is 8 : 44′. The vehicle which transports {J₁₀, J₁₁} are usable at 9 : 07′38″. So departure time of job J₁₃ is 9 : 07′38″. A vehicle free after 9 : 07′38″. The arrival time of whole batch is 9 : 29′14″.

Makespan of batch 5 is 46′36″. Batch 5 can be distributed after 8 : 46′36″. The earliest idle vehicle is valid until 9 : 27′12″ at batch 4, so departure time of job J₁₅ is 9 : 27′12″. The vehicle which transports J₁₄ at batch 4 is valid until 9 : 27′12″. Departure time of job J₁₇ is 9 : 27′12″. The vehicle which transports J₁₃ at batch 4 is valid until 9 : 50′12′. Departure time of job J₁₆ is 9 : 50′12″. The arrival time of whole batch is 10 : 09′.

Makespan of batch 6 is 46′54″. Batch 6 can be distributed after 8 : 46′54″. The earliest idle vehicle which transports J₁₅ at batch 5 is valid until 10 : 05′12″. Departure time of job J₂₀ is 10 : 05′12″. Departure time of job J₁₈ is 10 : 28′12″. Departure time of job J₁₉ is 10 : 28′12″. The arrival time of whole batch is 11 : 15′18″. The frequency of total schedule vehicle is 20.

Finally, we use Tables 5–7 to schedule the transportation start time in the second stage.

To compare the proposed method with the existing methods of HGA [28] and GA [26], we assume that the parameters of HGA [28] and GA [26] are the same as our algorithm.

The objective value obtained by running uncertain simulation and HGA is 16852. The scheduling result at the first stage is

Machine 1: 7 → 4 →15 → 17

Machine 2: 3 → 10 → 13 → 18

Machine 3: 5 → 6 →12 → 14

Machine 4: 8 → 1 →16 → 20

Machine 5: 2 → 9 →11 → 19

The objective value obtained by running uncertain simulation and GA is 16963. The scheduling result at the first stage is

Machine 1: 5 → 9 →13 → 14

Machine 2: 6 → 8 →10 → 19

Machine 3: 1 → 12 → 15 → 17

Machine 4: 7 → 3 →11 → 20

Machine 5: 4 → 2 →16 → 18

The results show that the three algorithms mix with uncertain simulation respectively are effective to deal with our problem, but GSA would reach a better result.

6 Conclusions

In this paper, we established a scheduling and distribution model for an uncertain two-stage supply chain problem. The model for the first stage problem was proposed in the form of pessimistic value under the assumptions that processing time and transportation time of the jobs are uncertain variables. The objective function includes holding cost, tardiness penalty cost, earliness penalty cost and transportation cost. It is the non-monotone function of processing time and transportation time. A novel method consisting of an uncertain simulation and GSA algorithm is used for solving the model. At last, we give a numerical example and make a comparison between our method and two methods in literatures. Results show that the three heuristic algorithms are effective and GSA is a better one.