Synchronized scheduling of production and outbound shipping using bilevel-based simulated annealing algorithm

Computers & Industrial Engineering 137 (2019) 106050

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Synchronized scheduling of production and outbound shipping using bilevel-based simulated annealing algorithm

T

Jian Chena,b, , George Q. Huangb, Jun-Qiang Wangc,d ⁎

a

College of Economics and Management, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, PR China HKU-ZIRI Lab for Physical Internet, Department of Industrial and Manufacturing Systems Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong c Performance Analysis Center of Production and Operations Systems (PacPos), Northwestern Polytechnical University, Xi'an, Shaanxi 710072, PR China d Department of Industrial Engineering, School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an, Shaanxi 710072, PR China b

ARTICLE INFO

ABSTRACT

Keywords: Interdivisional synchronization Production scheduling Shipping scheduling Bilevel model Simulating annealing

In general, different divisions even in a same company may make scheduling decisions irrespective of the overall benefit of the company. It is practically impossible to carry out the integrated optimization, and the sequential optimization results in isolated solutions within the divisions often leads to inferior overall solutions. In this study, a novel bilevel approach is proposed for joint synchronized production and shipping scheduling, considering production division as the leader and shipping division as the follower. A bilevel-based simulating annealing (SA) algorithm is proposed as the solution algorithm. Specifically, first, the production division prepares a tentative production schedule obtained by the SA. Based on the tentative schedule, shipping scheduling is addressed by an effective heuristic named Earliest Completion Machine. An inner SA is also presented for shipping scheduling to generate good quality shipping schedules. The tentative schedule is iteratively updated until the stopping criterion is satisfied. The effectiveness and efficiency of the proposed approach were verified via comparison with the traditional sequential and integrated approaches. Furthermore, the sensitivity analysis results indicate that shipping capacity has a considerable effect on the overall makespan. The allocation of the fixed shipping capacity onto blocks was also investigated, indicating that the large more strategy outperforms the average allocation and the large less strategies.

1. Introduction Many manufacturing companies operate in functional structure, in which the long-term planning of production and shipping is jointly carried out by the company, whereas the short-term scheduling is managed by individual divisions of the company. In reality, the synchronization of production and shipping scheduling, described as an interdivisional synchronized scheduling problem (ISP), is a crucial issue in many companies. ISP is of great importance, because it may incur many bad consequences such as high intermediate inventories, delay of order delivery, and low utilization of resources. A lot of efforts have been devoted to solving the ISP and are mainly classified into two types, (a) sequential approaches and (b) integrated approaches. In sequential approach, the upstream division first prepares production division scheduling. Then, the downstream division follows based on the generated production schedule. Sequential approaches are simple, but their performance from the systematic perspective is difficult to guarantee. To hoist the performance, in recent years, the integrated scheduling of production and

⁎

shipping has received extensive attention (Chen, 2010; Condotta, Knust, Meier, & Shakhlevich, 2013). Though the integrated scheduling of production and shipping can offer very competitive overall profits of the whole system, it faces difficulties in the real world. In many companies, each division has its own decision autonomy with individual performance evaluation, and the managers of each division make scheduling decisions irrespective of the convenience of other divisions. This study is originally motivated by the real-world production and shipping scheduling in our collaborating company. As shown in Fig. 1, the production division comprises a number of production lines to process various customer orders. Shipping order is managed by the shipping division in a space limited shipping area. The shipping area is divided into several blocks, each of which owns several loading bays. Each block is operated by its own shipping team, and its loading speed varies from each other. In addition, each block has different space sizes in the form of storage lanes. Shipping scheduling is to allocate orders onto blocks for typical shipping operations, e.g. checking, packaging, and truck loading.

Corresponding author at: Room 1114, 29 Jiangjun Avenue, Nanjing, Jiangsu Province, PR China. E-mail address: [email protected] (J. Chen).

https://doi.org/10.1016/j.cie.2019.106050

Available online 09 September 2019 0360-8352/ © 2019 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Fig. 1. Overview of ISP in the collaborating company.

In this study, the production division is nominated as the leader and the shipping division is nominated as the follower. The production division makes production schedule first and sends the tentative production schedule to the shipping division. The shipping division then prepares its shipping schedule based on the tentative production schedule. As the leader, the production division observes the functions of shipping division and updates the corresponding production schedule. A bilevel scheduling model was built, with each of its levels as a mixed integer programming (MIP) model. The bilevel procedure iterates via solving MIP at each level. This incurs heavy computational burden, which greatly hinders its wide reallife application in practice (Colson, Marcotte, & Savard, 2007). Thus, we developed an efficient bilevel-based simulated annealing algorithm (BSA). In the BSA, SA is employed to generate upper-level production schedules. An effective heuristic named Earliest Completion Machine (ECM) is proposed for the lower-level shipping scheduling. An inner SA is also presented as an alternative to ECM to generate a better shipping schedule when computation time is sufficient. Bilevel approaches are well known to generate worse solutions compared to centralized optimization, which is described as price of anarchy (POA) (Heydenreich, Müller, & Uetz, 2007). The proposed BSA was compared to the integrated approach to observe its gap with the ideal optimal solutions. The benefits of the bilevel approach are also demonstrated by comparing the BSA with the traditional sequential approach. The third question is answered through sensitivity analysis with respect to shipping capacity in two different situations. As aforementioned, the novelty of this study can be summarized. First, an ISP between the production and shipping is considered based on a real-life application case in practice. Second, a corresponding bilevel scheduling model is proposed, considering the respective decision autonomy of the production and the shipping divisions, different from many other scheduling models. Third, a bilevel-based SA algorithm, considering the leader-follower iterative procedure of the bilevel model, is proposed as the solution algorithm. The rest of the study is organized as follows. Section 2 briefly reviews the relevant literature. Section 3 formulates the problem as the bilevel MIP model. Section 4 discusses the proposed BSA. Section 5

At the planning horizon, the two divisions make individual schedule with their own objective based on the master production schedule. Each division aims to minimize its own makespan in order to maximize machine utilization and shipping block utilization (Pinedo, 2012). Currently, the two divisions have individual decision autonomy and operate scheduling in a sequential way. The production division first prepares the production schedule, and then shipping division responds to the production schedule and prepares its own shipping schedule. However, this sequential strategy lacks coordination leading to high inbetween work-in-process. The real situation prompted us to find out another synchronization way for the two divisions considering their selfish decisions with good coordination. In addition, the production scheduling problem seeks to allocate orders to the production lines and determine their production time slots to minimize manufacturing makespan and is regarded as a classic parallel machine scheduling problem. Shipping scheduling is to assign the block and shipping time slots to the finished orders to minimize the shipping makespan. It is somewhat similar to a uniform parallel machine if each block is viewed as a machine and each order as a job. Our shipping scheduling problem is more complex and further considers resource constraint for each block and release time for each order (Edis, Oguz, & Ozkarahan, 2013). The parallel machine scheduling problem is well known to be NP-hard even when there are only two machines (Pinedo, 2012). Thus, efficiently and effectively achieving synchronized scheduling is very challenging, considering the complex scheduling problems faced by the production and shipping divisions. Therefore, the objective of this study was to propose a bilevel approach for synchronized scheduling of production and shipping over a given planning horizon. More specifically, this study aims to answer four important research questions: (1) building a bilevel scheduling model for synchronized production and shipping scheduling; (2) developing efficient solution algorithm to address the complex bilevel scheduling model; (3) Investigating the benefits of the proposed bilevel approach compared to the sequential decisions and the centralized decisions; (4) The effect of the performance of major parameters on the synchronized production and shipping scheduling. 2

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

presents numerical studies and corresponding sensitivity analysis. Finally, Section 6 summarises the contributions and future outlook.

the leader managing the retailers’ inventory and retailer as the follower. BLP problems are intrinsically hard. Even the simplest instance has been shown to be NP-Hard problem. When the lower-level problem is convex and regular, it can be replaced by its Karush–Kuhn–Tucker (KKT) conditions, yielding the single-level reformulation of the problem. Then, the BLP can be solved by some exact algorithms such as branch-and-bound algorithm (Shi, Lu, Zhang, & Zhou, 2006). For the nonlinear BLP problems, most algorithms in literature can only guarantee the local optima (Li & Fang, 2012). El-Sobky and Abo-Elnaga (2018) presented a penalty method with trust-region technique for nonlinear BLP problem transform the nonlinear bilevel programming problem to unconstrained optimization problem. This method is simpler and does not need to compute a base of the null space. Recently, evolutionary algorithms are widely used to address complex optimization problems (Wang et al., 2018; Xiang, Xing, Wang, & Zou, 2019; Yi et al., 2018). Some of them are successfully adopted to solve the BLP problems (Angelo & Barbosa, 2015). Marinakis, Migdalas, and Pardalos (2007) proposed a genetic algorithm-based approach as the solution algorithm for solving BLP problem, in which both the upper and lower level problems are addressed by the designed GA. In this study, a bilevel-based simulated annealing is presented to find near-optimal solution to the BLP problem within the acceptable computation time.

2. Literature review The literature related to this study can be divided into two streams: (1) Synchronized scheduling of production and shipping and (2) Bilevel approaches in supply chain management. 2.1. Synchronized scheduling of production and shipping In the past two decades, most of the studies of integrated production and distribution focused on strategic and tactical planning levels. In contrast, the scheduling level research on integrated production and distribution is fairly recent (Chen, 2010; Chen, Huang, Luo, & Wang, 2015). For the production part, the single machine and parallel machine environments are mostly studied (Chen & Vairaktarakis, 2005; Delavar, Hajiaghaei-Keshteli, & Molla-Alizadeh-Zavardehi, 2010). Stecke and Zhao (2007) studied integrated production and delivery scheduling under commit-to-delivery business mode that earlier customer orders produced could choose a slow shipping mode to save shipping costs. They modelled the production as a capacity constrained single machine and considered two scenarios whether partial order delivery is allowed. Both scenarios are provided with a mixed integer programming model. Wang and Cheng (2000) considered a parallel machine scheduling with batch delivery costs. Jobs can be delivered in batches, and delivery cost depends on the number of deliveries. Hall and Potts (2005) considered both single and parallel machine scenarios with each batch delivered to the customer in a single shipping. However, for the shipping part, most studies in the area concentrated on transportation related decisions, e.g., vehicle allocation, routing, or batch delivery. The study on outbound shipping scheduling has hardly received any attention. Regarding the shipping scheduling problem, similar studies on machine scheduling literature are found, related to uniform parallel machine scheduling problem. In uniform parallel machine, each machine has its specific speed. Liao and Lin (2003) designed an optimal algorithm for two uniform machine scheduling to minimize makespan. Koulamas and Kyparisis (2004) investigated the uniform parallel machine scheduling taking account the release time of jobs. They developed a heuristic for this NP-hard problem and derived a tight worsecase ratio bound. Some studies considered the resource constraint for the uniform machine scheduling problem (Edis et al., 2013); however, limited resources are shared by machines, which still differs from our setting in which each machine is imposed with a specific space constraint.

3. Problem formulation and modeling 3.1. Problem formulation The concerned problem is divided into two divisions: the production division and shipping division. The production division consists of m parallel production lines for processing various customer orders. Assuming n customer orders to be processed during the given planning horizon, order j can be processed in any production line with the processing time pj , j = 1,2,. ..,n . The shipping division is responsible for dispatching the finished orders. The shipping division owns a shipping area, which is divided into b blocks with fixed loading bays. Each order is allowed to be dispatched in one shipment. The block is dedicated to dispatch customer orders. Once the space is allocated to an order, it cannot be released until the order is dispatched. The dispatching time (for checking, packaging, truck loading, etc.) of order j is given as t j . Since different blocks are operated by different shipping teams, the dispatching speed of each block varies. Assuming the dispatching speed of block h is sh indicates that if order j is assigned to block h , then it requires t j sh time units to be dispatched. Besides, the limited space of block h is denoted as Qh . Assuming order j requires storage space qj obviously indicates that the space constraint requires qj Qh , if order j is assigned to block h . In addition, customer orders can be available for shipping only after they are finished in the production division. Each customer order has a release time of r j , which equals to the completion time of this customer order Cj . Both production and shipping division are considered as selfish divisions, as they tend to optimize their own objectives when they generate their own schedules. The case company was interested in maximizing their resource utilization to increase throughput. In this study, makespan was used as the performance measure, as minimizing makespan is equivalent to maximizing machine utilization (Pinedo, 2012). Though makespan was used as the objective in this study, our model can also be used for other objectives with some revisions. Furthermore, the definition of the makespan for each division is clarified as follows.

2.2. Bilevel approaches in supply chain management Bilevel programming (BLP) is a branch of optimization models for hierarchical decision making with two non-cooperative decision makers: the leader at the upper level of the hierarchy and the follower at the lower level. BLP was first introduced in the Stackelberg game on market economy (Osborne & Rubinstein, 1994) and has been applied in many other areas including supply chain management. Calvete, Galé, and Oliveros (2011) addressed a hierarchical production–distribution planning problem in which two different decision makers control the production and the distribution processes. A bilevel model was proposed to solve this problem. Lukač, Šorić, and Rosenzweig (2008) investigated a production planning problem confronted by a pharmaceutical company. The senior manager wants to minimize the total setup time, whereas the middle manager focuses to minimize the cost of production, storage, and setup. This problem was solved by BLP where the senior manager acts as the leader and the middle manager as the follower. Yu, Chu, and Chen (2009) proposed a bilevel model for vendor managed inventory (VMI) system that treats manufacturer as

(a) The production makespan is defined as the production duration, which is the difference of the starting time and completion time of the production, Cmax = max (Cj ) . 1 j n

(b) The shipping makespan is defined as the shipping duration, which is 3

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

the difference of the starting time and completion time of shipping, max (C j ) min (r j ) .

1 j n

because min(r j ) is a constant number. Therefore, the shipping scheduling problem is a NP-hard problem, since its simple variant Qm r j Cmax is a NP-hard problem even for m = 2 (Koulamas & Kyparisis, 2004). In our collaborating company, the production division plays a leading role and pushes the downstream shipping operations. In our bilevel scheduling model, the production and shipping divisions are regarded as the leader and follower, respectively. It is reasonable to delegate the production division to control the performance of the entire production and shipping process, since the leader has the opportunity to know the reactions of his follower when he attempts to maximize his objective. Hence, the shipping division’s objective was incorporated in the leader’s own objective to minimize the overall makespan of production and shipping, Cmax . This global measure can also decrease the waiting time between two divisions to improve the interdivisional synchronization. The production division as the leader first determines a production schedule and outputs the production completion time Cj . Then, the shipping division begins to optimize its objective Cmax min(r j ) based on the release time r j = Cj . The obtained shipping schedule updates the leader’s objective, Cmax . With the response of shipping division, the production division has to adjust its decisions correspondingly. Based on the bilevel procedure, the bilevel scheduling model is formulated as follows.

1 j n

(c) The overall makespan of the production and shipping is defined as the duration of the overall production and shipping process and is equal to the difference of the starting time of production and completion time of the shipping, Cmax = max (Cj ) . 1 j n

The problem is to find a synchronized production and shipping schedule in the case where each division owns individual decision autonomy. More specifically, the synchronized schedule consists of production and shipping schedules. The production schedule determines which production line to process each order and when. The shipping schedule determines which block to dispatch each order and when. Further assumptions are as follows: (1) Production scheduling starts at time zero. (2) Production line (or block) can process no more than one order at a time. (3) Each order can be processed on one production line (or block). (4) Pre-emption of jobs is not allowed for the two divisions. (5) Setup times of production and shipping are ignored. (6) Each order is available for shipping once it is completed in the production division. (7) Inner transportation time is negligible. The notations used in model formulation are as follows. Parameters – – – – – – – – –

Subject to m ij

j : index of order; j = 1, 2, ...,n i : index of production line at production division, i = 1, 2, ...,m h : index of block at shipping division; h = 1, 2, ...,b pj : processing time of order j t j : dispatching time of order j qj : space requirement of order j Qh : space of block h sh : dispatching speed of block h A : A large number

j = 1, 2, ...,n

(1b)

n ijk

=

ij ,

j

k , i = 1, 2, ...,m , j = 0, 1, ...,n

ijk

=

ik ,

j

k, i = 1, 2, ...,m , k = 0, 1, 2, ...,n

k= 0

(1c)

n j=0

Ck

Cj

A(

1)

ijk

pk ,

j

i0

= 1,

(1e) (1f)

i = 1, 2, ...,m

(1g)

C0 = 0

{0, 1},

ijk

(1d)

k, i = 1, 2, ...,m , j = 0, 1, ...,n,

k = 1, 2, ...,n

– Cj : completion time of processing order j – Cj : completion time of dispatching order j – r j : release time of order j to shipping division, irrespective of inner transportation time, r j = Cj . – Cmax : the latest completion time of production, Cmax = max (Cj )

j

k, i = 1, 2, ...,m , j = 0, 1, ...,n, k = 0, 1, 2, ...,n (1 h)

1 j n

– Cmax : the latest completion time of shipping, Cmax = max (Cj ) ijk = 1, if order otherwise. – ij = 1, if order j – Xhjk = 1, if order – Yhj = 1, if order j

= 1,

i=1

Decision Variables

–

(1a)

Min Cmax

{0, 1},

ij

1 j n

j precedes order k on production line i , and 0

i = 1, 2, ...,m , j = 1, 2, ...,n

(1i)

where Cmax solves by the following model with rj = Cj.

is assigned on production line i , and 0 otherwise. j precedes order k on block h , and 0 otherwise. is assigned on block h , and 0 otherwise.

Min Cmax

(2a)

min(r j )

Subject to b

Yhj = 1, 3.2. Bilevel scheduling model

j = 1, 2, ...,n

(2b)

h= 1 n

If each production line is regarded as a machine and each order as a job, the production scheduling problem can be denoted as Pm Cmax by using the three-field notation a b c conventionally used in the production scheduling literature. This problem is well known and widely studied. It is easy to see that P 2 Cmax is NP-hard in the ordinary sense, as it is equivalent to PARTITION problem (Pinedo, 2012). The shipping scheduling can also be regarded as a machine scheduling problem. The block can be considered as one machine, with each order as a job. The problem is denoted as Qm resi, r j Cmax min(rj ) in this study. The complexity of the problem is identical to Qm resi, r j Cmax

Xhjk = Yhj,

j

k, h = 1, 2, ...,b, j = 0, 1, ...,n

Xhjk = Yhk ,

j

k, h = 1, 2, ...,b, k = 0, 1, 2, ...,n

k= 0

(2c)

n j=0

Ck

4

Cj

A (Xhjk

1)

tk sh,

j

(2d)

k, h = 1, 2, ...,b, j = 0, 1, ...,n,

k = 0, 1, 2, ...,n

(2e)

Cmax

(2f)

Cj ,

j = 1, 2, ...,n

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

greater than the completion time of order j plus processing time of order k if order k is immediately processed after order j on the same block. This constraint also indicates overlap does not exist when dispatching orders. Constraint (2f) imposes shipping makespan to be greater than the completion of all the orders. Constraint (2g) imposes that each order cannot start dispatching before its release time, which is the completion time of production. Constraint (2h) imposes that the order j cannot exceed the space capacity of block h when order j is shipped by block h . Constraint (2i) similarly defines a dummy job placed on the start and end of a sequence of orders on a block. Constraints (2j) and (2k) define the binary variables Xhjk and Yhj , respectively. 4. Bilevel-based SA The bilevel structure, even in the simplest form, with linear upper and lower level programming problems (known as linear BLPP), is NPhard and difficult to solve. In this problem, the production and shipping scheduling problem are both NP-hard. Finding an optimal solution in the polynomial time is even harder. In this circumstance, heuristic methods could serve as alternative solution algorithms, since they heuristically search within a portion of the solution space, generating good quality solutions in acceptable computation time; however, the solution optimality is not guaranteed. More specifically, a bilevel-based simulated annealing algorithm (BSA) was used to address the above bilevel scheduling model. Here, SA was adopted based on the following two reasons. First, SA has been introduced as an efficient heuristic algorithm in solving many scheduling problems (Ishibuchi, Misaki, & Tanaka, 1995; Lin & Ying, 2015). Second, SA works well in neighborhood search, which is applicable to the local exploration of tentative production schedules. BSA is implemented regarding the characteristics of the bilevel scheduling model. In the bilevel scheduling model, the leader makes his decision first. The follower’s decision depends on the leader’s decision, and the leader has full knowledge of the schedule. To achieve the bilevel decision process, the two optimization problems were iteratively linked together. First, SA was used to generate a tentative production schedule. Then, based on the optimized production schedule, a shipping schedule will be obtained by a proposed a heuristic named Earliest Completion Machine (ECM) based on List Scheduling. Then, the tentative production schedule is locally explored by SA for further improvement, considering the shipping schedule. This procedure will be iterated until the leader’s stopping criterion is satisfied. In contrast, an inner SA is proposed as an alternative to ECM. Inner SA can effectively improve the performance of BSA due to better search for shipping scheduling problem. Unfortunately, it may be more suitable for smallscale instances, because it is fairly time-consuming. Fig. 2 shows the flowchart of the proposed BSA.

Fig. 2. Flowchart of the proposed BSA.

Cj

Yhj t j sh

Qh

A (Yhj

Yh0 = 1, Xhjk

Yhj

rj, 1)

h = 1, 2, ...,b, j = 0, 1, ...,n qj ,

h = 1, 2, ...,b, j = 0, 1, ...,n

h = 1, 2, ...,b

{0, 1},

{0, 1},

j

k, h = 1, 2, ...,b, j = 0, 1, ...,n, k = 1, 2, ...,n

h = 1, 2, ...,b, j = 1, 2, ...,n

(2g) (2h) (2i) (2j) (2k)

In the formulation, the leader’s objective (1a) is to minimize the overall makespan of production and shipping. Constraint (1b) ensures that each order is only assigned on one production line. Constraints (1c) and (1d) ensure that if an order j is assigned to a production line, there must be exactly one other order that precedes and one other order that succeeds order j on the production line. Constraint (1e) enforces that the completion time of order k must be greater than completion time of order j plus processing time of order k if order k is immediately processed after order j on the same production line. This constraint also indicates that overlap does not exist when processing orders. Constraint (1f) defines a dummy job placed on the start and the end of a sequence of orders on a production line. Constraint (1g) places dummy job at time 0. Constraint (1h) and (1i) define the binary variables ijk and ik . Eqs. (2a)–(2k) belong to the followers’ model. The follower’s objective (2a) aims to minimize shipping makespan, which is defined as the duration between the completion time of the latest order and release time of the earliest order. Constraint (2b) ensures that each order can only be assigned on one block. Constraints (2c) and (2d) ensure that if an order j is assigned to a block, there must be exactly one other order that precedes and one order that succeeds order j on the block. Constraint (2e) enforces that the completion time of order k must be

4.1. SA for production scheduling SA is a randomized local search heuristic that allows for random uphill movements to prevent the algorithm from getting trapped at local optima. The level of randomization is determined by cooling schedule from initial temperature to zero. The logarithmic cooling schedule was used for the SA with the key parameter, namely, cooling factor c (Hajek, 1988). The stopping criterion is set as maximum rounds of iteration NF . The detailed procedure of SA is described as follows. Step 1: Initialization. An initial solution X is generated randomly. The initial solution is encoded as a permutation of customer orders. A preceding customer order in a solution list has a higher priority to be assigned to machines than the subsequent one. Take a solution \{ 2,5,4,1,3\} for example, there are five orders to be processed in this sequence where order 2 is processed first, followed by orders 5, 4, 1, and 3. The initial iteration index is set as zero, i.e., N = 0. 5

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Step 2: Neighbor generation. To generate a neighboring solution, swap operator or insert operator is often used in the literature. The swap neighborhood operator was used in our SA by randomly choosing a pair of orders to swap their positions in the solution list. Let Xg denote the neighboring solution at the gth iteration. Step 3: Evaluation. The leader’s objective function Cmax was used to evaluate the performance of solutions. Denote the performance of the solution Xg by f (Xg ) . In order to calculate f (Xg ) , decoding was conducted to generate production and shipping schedules. To generate a production schedule, List Scheduling (LS) heuristic was used for decoding. Given a solution Xk , orders are assigned one by one according to the list of the solution. Whenever a machine becomes available, the next order on the list is assigned on that machine. It can be shown that the worst case of an arbitrary list scheduling for Pm Cmax is 2 1 m (Pinedo, 2012). By doing so, the solution Xk is decoded and the production makespan, Cmax can be calculated. In contrast, for Cmax , the shipping schedule should be obtained, as introduced in the next section. Step 4: Accept the neighboring solution Xg as the current solution with the following principle. This step enables SA to avoid local optima by probabilistically taking non-locally optimal steps in the search space. a. If the f (X ) f (Xg ) 0, then the current solution X is replaced with the new solution Xg. b. If the f (X ) f (Xg ) < 0, then the current solution X is replaced with the new solution XN by the probability function, exp((f (X ) f (Xg )) Tg ) . In detail, a random number, rand , was generated where 0 < rand < 1. If exp((f (X ) f (Xg )) Tg ) > rand , the new solution Xg is accepted, otherwise the current solution X is retained. Step 5: Update the temperature Tg . The temperature Tg at the gth c iteration is updated by the logarithmic cooling schedule Tg = log(g + 1) and N = N + 1. Step 6: Repeat steps 2–5 until the stopping criterion is satisfied. The stopping criterion is set as the maximum rounds of iteration NF . If the index of the current iteration is N > NF , the repeat is stopped.

Table 1 Basic information of orders and blocks. j

1

2

3

4

5

rj qj

1

2

2

3

4

pj

4

4

4

4

2

h

1

2

3

Qh sh

1 1

2 2

3 2

3

2

1

2

3

Fig. 3. Illustration of ECM.

Example 1. The basic information of example 1 is given in Table 1, and the obtained schedule is shown in Fig. 3. An order list of the five orders was obtained as \{ 1,2,3,4,5\} in the non-decreasing order of r j . Order 2 was placed before order 3 since q2 > q3 , though orders 2 and 3 are simultaneously released. Then, orders are assigned one by one. Order 1 can only be placed on block 3 due to space constraint. Order 2 can be completed earlier on block 2 than block 3. Therefore, order 2 was assigned on block 2. When assigning order 3, though block 3 can complete the order same as block 1 (speed s1 is larger), order 3 was put on block 1 (with smallest space). Similarly, orders 4 and 5 were assigned block 3. The makespan is equal to 7.

4.2. Heuristic for shipping scheduling Shipping scheduling is embedded in the search loop of production scheduling. Efficiency is of great importance in addressing the inner shipping scheduling while considering the computing burden of the overall BSA. Thus, heuristic seems to be a very good choice. Herein, a heuristic named ECM is proposed. Heuristic ECM sorts orders in the non-decreasing order of r j and assign the sorted orders one by one to the least-makespan block meeting space constraint. Let h denotes the order set assigned on vehicle h ; t ( h) is the completion time of the schedule Sh . If several blocks can simultaneously complete the order, the one with the smallest space is chosen. The basic procedure is as follows.

Lemma 1. The shipping schedule can be obtained by ECM in O (n log n + nb 2 log b) time. Proof. First the sorting can be done in O (n log n) time. Then, for each of n orders, comparison and sorting were performed for b blocks; therefore, the running time is O (nb2 log b). Thus, the total running time is O (n log n + nb 2 log b) . □

Heuristic ECM Input: parameters of orders and blocks: r j , t j , qj , Qh , where h = 1, ...,b , j = 1,2,….,n;

4.3. Convergence analysis of the proposed bilevel-based SA

Ouput: the shipping schedule h , h = 1, ...,b ; Begin Sort jobs so that r1 ≤ r1 ≤…rn; For h = 1 to b h= ; End for For j = 1 to n h= min {max (t ( h ), r j ) + t j } ;

In the section, the global convergence property of the proposed bilevel-based SA is discussed. As indicated by Hajek (1988) and Yang (2000), SA algorithm for discrete or combination optimization can be characterized by a discrete-state Markov chain. To ensure the global convergence of SA algorithm, the considered SA algorithm should be irreducible and weak reversible (Hajek, 1988). In the proposed bilevel-based SA, the swap neighborhood operator develops a new neighbor, which is obviously reachable from the current initial solution. Thus, the proposed bilevel-based SA satisfies the irreducible property. Then, if the weak reversibility is satisfied, the global convergence can be achieved. Therefore, the global convergence will be satisfied if and only if

h

h Qh qj

= h {j } ; End for End for End h

6

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

exp( d Tg ) = +

Insertion operator

(3)

g=1

Input: order j and block h : r j , t j , qj , Qh , T (

where d is the maximum of the depths of all solutions, which are local but not global optimum, i.e., the difference between the best local optimum solution with the worse solution. c As mentioned in Step 5 in Section 4.1, Tg = log(g + 1) , then g=1

=

t1s, t2s, ...,tKs and t1c, t2c , ...,tKc

h;

Output: Ch , the makespan after the order j is arranged onto block h ; Begin For k = 1 to K If k = K Put order j after order K Update T ( h) and h Ch = max (t ( h), rj ) + t j )

c

exp( d Tg ) = g=1

h)

(g + 1) d , the condition (3) is ture if and only if

c d , indicating that the convergence of the proposed bilevel-based SA can be satisfied under the condition where c is larger than d . The global optimal solution can be theoretically obtained only if the number of iterations is sufficient. It needs a large amount of computation time, which is definitely not preferred in practice (Cohn & Fielding, 1999). The cooling schedule as well as the maximum number of iterations will be carefully determined with the trade-off of the computation time and performance. In addition, the procedure for determining the key parameters in numerical studies is demonstrated. On the other hand, obviously the proposed bilevel-based SA inherits the good feature of SA algorithm with which it could escape local optima by allowing hill-climbing moves oriented at finding global optimum. The hill-climbing moves accept solutions with worse fitness with certain probability. With temperature decreasing to zero, hillclimbing moves occur less frequently.

Else if tkc+ 1

t j then

max(tkc, r j )

Put order j after order k Ch = tKc Update T ( h) and h End if End for End

The procedure of the revised ECM is given as follows. Revised ECM Input: parameters of orders and blocks: r j , t j , qj , Qh , where h = 1, ...,b , j = 1, 2, ...,n ; Ouput: h and T ( h) , h = 1, ...,b ; Begin Sort jobs so that r1 ≤ r2 ≤ …rn; For h = 1 to b and T ( h) = ; h= End for For j = 1 to n

4.4. Extension of BAS with inner SA

Use insertion operator to calculate Ch , h∈{h|Qh ≥ qj};

Though ECM may quickly find a shipping schedule, it is hard to find an optimal shipping schedule, which may affect the performance of BSA. This section introduces an inner SA for shipping scheduling. The numerical results indicate that inner SA can generate better BAS solutions with more computation time. The procedure of inner SA is similar to the SA for production scheduling except initialization, neighbour generation, and evaluation.

Then h =

min h

{Ch} ;

h Qh qj

Put order j on block h based on insertion operator; Update h and T ( h) ; End for End

Example 1 in Section 4.2 was used to illustrate the active decoding (Fig. 4). Assume an order list \{ 2,5,4,1,3\} to be decoded. According to the revised ECM, order 2 is first assigned to block 2 at its release time. Order 5 is then assigned to block 5, followed by order 4 on block 2. When assigning order 1, block 3 was found as the only one satisfying the space constraint. Then, all the idle slots were searched on block 3 and insert order 1 at its release time, which is before order 5. Obviously, the schedule will not be active if the idle slots were not searched on block 3 and order 1 was just put after order 5.

Step l: Initialization. The initial solution is encoded as a permutation of customer orders. Initial solution X was randomly generated. Cooling factor c and maximum rounds of iteration NF are prescribed. The initial iteration index is set as zero, i.e., N = 0 . Step 2: Neighbor generation. The swap operator is used that randomly chooses a pair of orders and swaps their positions in the solution list to generate a neighboring solution Xg at gth iteration. The distance of the picked two orders is set less than a given control parameter. In our SA, the control parameter is given as the number of the blocks. For example, we assume orders j and k are picked out for swapping. Let j denote the position of order j in the order list. b was set, where b is the number of blocks, j k j, k = 1, 2, ...,n . This restriction is beneficial for taking advantage of the initial solution generated by heuristic ECM. In addition, it reduces the search space to make SA faster. Step 3: Evaluation. The followers’ objective Cmax min(r j ) is used to evaluate the performance of solutions.

Step 4: Accept the neighboring solution as the current solution according to the principle defined in Step 4 of the SA in Section 4.1. Step 5: Update the temperature Tg . The temperature at iteration gth,

Block

h =3

To decode a solution list, heuristic ECM is revised for generating active shipping schedule. When assign an order on a block, an insertion operator was considered to search the idle slots of the current schedule and the order was put within the earliest-capable time slot, otherwise it was put at the end of the current schedule. To implement the insertion operator, let h denotes the job set assigned on vehicle h ; t ( h) is the completion time of the schedule h . Let t s, t s, ...,tKs T ( h) = 1c 2c be a matrix representing the timetable of the t1 , t2 , ...,tKc schedule h , which specifies the starting time and completion time of its orders in the first and second row, respectively. For instance, t1s denotes the starting time of the first order on block h , while t1c .enotes the completion time of this order.

j=1

h =2

j=2

h =1

j=3

1

2

3

j=5

Q3

j=4

Q2

Q1

4

6

Fig. 4. Active decoding using the revised ECM. 7

Time

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Tg is updated by the logarithmic cooling schedule Tg = log(g + 1) , and N = N + 1. Step 6: Repeat steps 2–5 until the stopping criterion is satisfied. The stopping criterion is set as the maximum rounds of iteration NF . If the index of current iteration N > NF , the repeat is stopped. c

Table 3 Three level factorial DOE. Experiment No.

1 2 3 4 5 6 7 8 9

5. Numerical studies The numerical studies consist of four parts. The first part adopts design-of-experiment (DOE) to determine the optimal parameters for BSA. The second part compares the proposed BSA with the traditional sequential approach and integrated approach. The purpose of the comparison is twofold: (1) to verify the effectiveness of BSA, (2) to demonstrate the value of bilevel model. The third part compares ECM and inner SA in different instance sizes in terms of performance and computation time. The fourth part conducts sensitivity analysis to investigate whether shipping capacity has a significant impact on the overall objective. The computational studies were conducted using a Dell OptiPlex PC with Intel i7 CPU 3.4 GHz and 8 GB RAM. The sequential and integrated models were implemented by CPLEX 12.5 for comparison. The BSA was implemented by MATLAB R2012a software. To generate the test instances, the following basic settings were used. The processing time pj and dispatching time t j are generated from the uniform distribution [2, 10] and [10, 20], respectively. The space requirement and space capacity Qh generated from the uniform distribution [2, 5]. The dispatching speed sh is generated from the uniform distribution [1, 3]. Test instances are randomly generated considering the different instance sizes. To generate the small-scale instances, the number of customer orders (n), the number of machines (m), and the number of blocks (b) are taken as n = 7, 8, 9, 10, m = 2, 4, 6, and b = 2, respectively. For the large-scale instances, n = 20, 40, 60, 80 customer orders, m = 10, 15, 20 machines, and b = 10 blocks were considered.

NF c

25,000 30

50,000 50

1 2 3 1 2 3 1 2 3

39.5 39.7 39.5 39.0 39.1 38.9 38.7 38.7 38.9

(4a)

Subject to constraints (1b)

Min Cmax

(4b)

(1i)

(4c)

min(r j )

Subject to constraints (2b)

(4d)

(2k)

The production scheduling model ((4a) and (4b)) generates the completion time of orders, Cj , j = 1, 2, ...,n . The release time of shipping orders are obtained as r j = Cj in the shipping scheduling model ((4c) and (4d)). Then, the overall production and shipping schedules can be obtained by solving the two models sequentially. (Integrated model) (5a)

Min Cmax (1i) and (2b)

(2k)

(5b)

In the following experiments, both sequential model and integrated model are solved by the widely used commercial software CPLEX. In order to provide a comprehensive comparison, for the small instances, the optimal solutions were obtained for both sequential and integrated model solved by CPLEX. For the large instances, the running time of CPLEX was limited for the sequential model, integrated model within one hour. The follower’s objective, the shipping makespan (S_MS) and

Factor levels

5000 10

1 1 1 2 2 2 3 3 3

Min Cmax

Subject to constraints (1b)

3

c

The proposed BSA was compared to the traditional sequential approach and integrated approach. The integrated approach is presented to show the performance of BSA, though it is not practically suitable for real-life application when the production and the shipping divisions have their own decision autonomy. Compared to centralized coordination, selfish behavior will deteriorate the overall system performance and is known as price of anarchy (POA) in literature (Heydenreich et al., 2007). In our setting, the POA is the worst-case ratio of the maximum makespan by the decentralized approach to the optimal makespan by the integrated approach. Since it is difficult to find the identical problem in literature, the sequential as well as the integrated model was built based on the proposed bilevel scheduling model. (Sequential model)

Table 2 BSA parameters and their levels in DOE.

2

NF

5.2. Comparison of BSA to sequential and integrated approaches

The proposed BSA has two key parameters, namely, cooling factor (c) and maximum number of iterations (NF). To investigate the effect of these parameters on the performance of the proposed BSA, the three level full factorial DOE was carried out using a moderate-scaled instance with n = 60, m = 15 and b = 10. Each parameter is regarded as a factor, and three factor levels are considered for each factor. The 32 design, listed in Table 2, is selected with two factors, each at three levels. The nine treatment combinations for this type of design are listed in Table 3. The overall makespan of production and shipping, PS_MS, was used as the output measure. Each instance is run three times, and the mean values of PS_MS are listed in Table 3. The analysis of variance (ANOVA) assumes that the observations are normally and independently distributed with common variance for each treatment. These assumptions were investigated before performing the DOE. First, normal distribution test was performed with 95% confidence interval for the samples. Fig. 5 shows that the P-value is greater than 0.05, demonstrating the normality assumption. Second, Bartlett’s test was performed with 95% confidence interval to test if the samples are from populations with equal variances. Regarding

1

PS_MS

parameter NF and parameter c, the P-values are obtained as 0.652 and 0.882, respectively, demonstrating the equal variance assumption. Besides, each treatment combination is run independently generating independent samples for our experiment, thus satisfying all the required DOE assumptions. Finally, Fig. 6 shows that the effect of NF on PS_MS is more significant than that of c. The larger the value of NF, the better the solution BSA generates. Smaller c was found to be better. Thus, the parameters of BSA were set as follows: maximum rounds of iterations NF = 50,000 and cooling factor c = 10 for the subsequent experiments.

5.1. Parameter tuning

Factors

Factor levels

8

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Fig. 5. Normal probability plot.

Fig. 6. Main effect plot of BSA parameters.

the leader’s objective, the overall makespan of production and shipping (PS_MS), and the computation time (CPT) were used as the measures. Each instance is run three times, and the mean values are listed in Table 4. For the small instances, the fitness values of most BSA solutions were found identical to the optimal ones generated by CPLEX, e.g., the instances n = 8, 9, indicating that BSA is capable of finding the optimal solutions, considering the selfish decision of the shipping division. For these instances, n = 7, 10; however, BSA failed to find the optimal solutions. The average PS_MS of the sequential, integrated, and bilevel approaches for small instances were easily obtained as 31, 28.8, and 29.3, respectively. For large instances, the average PS_MS of the

sequential, integrated, and bilevel approaches were obtained as 51.7, 68.5, and 35.1, respectively. Then, the average POA of BSA was calculated and compared to the integrated approach for small-scale instances as 29.3/28.8 = 1.02. Furthermore, BSA outperforms the sequential approach for all small instances except the instance n = 7, m = 6. The average PS_MS generated by BSA is approximately 5% better than the one obtained by the sequential approach. For the large instances, BSA outperforms sequential and integrated approaches in most instances. The average PS_MS obtained by BSA is 32% and 49% better than the ones obtained by the sequential and integrated approaches, respectively. Note that for large instances, for 9

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Table 4 Comparison to sequential and integrated approaches. N

m

Small instances 7 2 4 6

Sequential

Integrated

BSA with ECM

S_MS

PS_MS

CPT (s)

S_MS

PS_MS

CPT (s)

S_MS

PS_MS

CPT (s)

16.3 15.3 16.0

26.3 27.0 23.7

0.8 0.6 0.2

14.0 13.7 16.0

24.0 23.7 23.7

1.6 1.6 2.4

21.0 21.0 21.0

25.0 25.0 25.0

36.3 39.0 38.1

8

2 4 6

19.2 21.3 19.3

32.5 30.7 28.0

3.2 0.6 0.3

17.7 17.7 17.7

31.7 27.0 27.0

26.0 21.1 949.1

28.7 24.0 24.0

31.7 27.0 27.0

43.2 42.0 43.3

9

2 4 6

24.0 21.5 22.0

29.3 28.5 29.3

31.1 6.1 4.4

22.7 22.7 22.7

28.0 28.0 28.0

212.0 480.7 2518.3

26.0 26.0 26.0

28.0 28.0 28.0

44.6 45.9 45.7

10

2 4 6

40.0 33.0 32.3

46.0 35.0 35.3

439.1 289.9 204.1

34.7 32.0 32.0

37.0 34.0 34.0

14402.2 15032.0 14577.2

33.0 32.0 32.0

39.0 34.0 34.0

50.4 50.7 50.6

13.7 11.7 8.7

26.0 19.7 18.7

1809.6 1142.8 1803.2

14.7 11.3 11.5

22.0 18.7 17.5

3600.0 3600.0 3600.0

19.4 16.0 17.0

22.0 18.3 19.0

126.5 134.9 132.4

Large instances 20 10 15 20 40

10 15 20

32.0 26.0 25.7

41.0 35.7 32.0

3600.0 3600.0 3600.0

39.3 22.7 25.0

51.3 30.7 32.0

3600.0 3600.0 3600.0

29.4 26.3 25.3

31.4 28.3 27.3

271.1 267.9 257.1

60

10 15 20

46.7 47.2 47.0

55.7 53.5 54.3

3600.0 3600.0 3600.0

79.2 65.3 61.0

88.5 79.7 69.7

3600.0 3600.0 3600.0

40.7 36.9 35.9

42.7 38.9 37.9

390.6 393.1 392.7

80

10 15 20

134.3 61.7 60.3

140.3 72.7 70.7

3600.0 3600.0 3600.0

126.3 118.7 128.3

135.0 136.0 140.3

10800.0* 10800.0* 14400.0*

52.4 48.8 48.2

54.4 50.8 50.2

544.9 541.2 538.8

[ ]* means no feasible solution is found until CPLEX run [ ] seconds. The bold values are the best one of PS_MSs obtained by the sequential approach, the integrated approach and BSA with ECM for each instance.

instance n ≥ 80, CPLEX failed to find any feasible solution within one hour. The computation time of BSA was satisfactory and much less than those of the sequential and integrated approaches, especially for the large instances, providing a fair opportunity for real-life application. Overall, the obtained experimental results demonstrate the value of BSA in terms of both solution quality and computation efficiency.

(2) Regarding the fixed shipping capacity, how to allocate them to blocks to maximize the system performance. In this section, the following sensitivity studies were conducted for some useful managerial implications from the two perspectives. 5.4.1. Effect of shipping capacity A two-way ANOVA was conducted to evaluate whether shipping capacity affects the system performance under different instances. The significance level was set as 95%. The value of n was typically considered as 40 and 60 customer orders, m = 10, 20 machines, and b = 10 blocks to generate the instances. Besides the basic settings, the shipping capacity was altered via assigning a dispatching speed factor K. The dispatching speed sh was then generated from the uniform distribution [1, K]. Obviously, the larger the value of K, the larger the shipping capacity. Four levels of K ranging from 1, 2, 3, and 4 were considered for each setting. A total of 16 test instances were generated. Table 6 shows the results of PS_MS, S_MS, and CPT. The ANOVA assumptions were investigated. The normal distribution test was performed within 95% confidence interval, leading to the P-value of 0.210, which satisfies the normality assumption. Then, Bartlett’s test was performed with 95% confidence interval to test if the samples are from populations with equal variances. Regarding instance size factor and dispatching speed factor K, the P-values were obtained as 0.935 and 0.764, demonstrating the equal variance assumption. The detailed results of assumption checking are not listed here due to the length limit of the article. Table 7 shows the ANOVA results in terms of PS_MS. The P-value of K and P-value of instance size are both less than 0.05. Both the shipping capacity and instance size have significant effect on PS_MS. The similar results were obtained for S_MS, but are not listed here. Both S_MS and PS_MS considerably decrease with increasing shipping capacity (K). The decrease is quite sharp initially for K from 1 to 2,

5.3. Performance of inner SA The inner SA may generate better shipping schedules than heuristic ECM; but the overall BAS performance needs to be evaluated when embedded with inner SA. The performance of BAS embedded with ECM and inner SA is illustrated in Table 5. S_MS, PS_MS, and CPT are used as the performance measures. The running time of the two BAS algorithms was limited within one hour for a fair comparison. Table 5 shows that inner SA generates better solutions (approximately by 3% on average) than ECM in terms of PS_MS for small instances. ECM generates the same PS_MS with inner SA for 66.7% of small instances. However, inner SA requires much more time than ECM. For large instances, ECM shows obvious advantage, as it offered better solutions for most cases with much less computation time than inner SA. Overall, the results indicate that inner SA is preferred for BAS in the case of sufficient computation time. 5.4. Sensitivity analysis In the shipping division, shipping capacity is related to the shipping resources. Usually, the more shipping resources, e.g., packaging workers, delivery worker, and forklifts, the more the shipping capacity and faster dispatching speed. The following two questions were addressed: (1) Whether shipping capacity affects the system performance; 10

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al.

Table 5 Comparison of ECM and inner SA. n

m

Table 7 The ANOVA table for shipping capacity.

BSA with ECM

Small instances 7 2 4 6

BSA with inner SA

S_MS

PS_MS

CPT (s)

S_MS

PS_MS

CPT (s)

21.0 21.7 21.0

25.0 25.0 25.0

36.3 39.0 38.1

22.0 21.0 20.0

24.0 23.7 23.7

3600.0 3600.0 3600.0

8

2 4 6

28.7 24.0 24.0

31.7 27.0 27.0

43.2 42.0 43.3

25.7 24.0 24.0

31.7 27.0 27.0

3600.0 3600.0 3600.0

9

2 4 6

26.0 26.0 26.0

28.0 28.0 28.0

44.6 45.9 45.7

26.0 26.2 26.0

28.0 28.0 28.0

3600.0 3600.0 3600.0

10

2 4 6

33.0 32.0 32.0

39.0 34.0 34.0

50.4 50.7 50.6

35.0 31.0 32.0

37.0 34.0 34.0

3600.0 3600.0 3600.0

19.4 16.0 17.0

22.0 18.3 19.0

126.5 134.9 132.4

19.9 16.7 16.8

22.5 19.0 18.8

3600.0 3600.0 3600.0

Large instances 20 10 15 20 40

10 15 20

29.4 26.3 25.3

31.4 28.3 27.3

271.1 267.9 257.1

30.6 27.4 26.8

32.9 29.4 28.8

3600.0 3600.0 3600.0

60

10 15 20

40.7 36.9 35.9

42.7 38.9 37.9

390.6 393.1 392.7

43.0 38.4 37.1

45.3 40.4 39.1

3600.0 3600.0 3600.0

80

10 15 20

52.4 48.8 48.2

54.4 50.8 50.2

544.9 541.2 538.8

54.5 49.9 49.2

56.5 51.9 51.2

3600.0 3600.0 3600.0

40

60

K

S_MS

PS_MS

CPT (s)

10

1 2 3 4

43.0 44.0 42.0 41.0

96.0 68.0 44.3 41.5

98.0 71.0 46.3 43.5

399.2 413.6 418.0 405.2

20

1 2 3 4

25.0 24.0 25.0 24.0

94.0 68.0 39.0 31.5

96.0 70.0 41.0 33.5

436.8 421.7 416.5 399.9

10

1 2 3 4

55.0 54.0 56.0 53.0

122.0 90.0 57.7 55.0

124.0 92.0 59.7 57.0

549.5 551.6 556.0 545.2

20

1 2 3 4

30.0 31.0 31.0 29.0

121.0 87.0 51.0 38.3

123.0 89.0 53.0 40.3

541.7 495.4 488.3 486.1

MS

F

P

K Instance size

3 3

11293.8 1337.2

3764.61 445.73

129.46 15.33

0.000 0.001

Error Corrected Total

9 15

261.7 12892.7

29.08

5.4.2. Effect of allocation strategies for fixed shipping capacity This section investigates allocating the fixed shipping capacity onto each block to minimize the PS_MS. Typically n = 40, 60 customer orders, m = 10, 20 machines, and b = 10 blocks were considered. The total shipping capacity was set as SC = 30 . The following three general allocation strategies were compared: (1) Average allocation (AA) strategy: It allocates average shipping capacity to each block. In this experiment, the dispatching speed of each block is set as SC b . (2) Large more (LM) strategy: Large block is allocated with more shipping capacity. The dispatching speed of block h is generated as SC Qh , where Qh is the space of block h . b h=1

Qh

(3) Large less (LL) strategy: Large block is allocated with less shipping Qh capacity. The dispatching speed of block h is generated as SC . b h= 1

1 Qh

Table 8 shows the results of P_MS, S_MS, PS_MS, and CPT. Fig. 7 shows that the LM strategy obtains the minimum PS_MSs for all instances, clearly indicating that allocating more shipping capacity to large block is the best allocation strategy among the three strategies. It is reasonable that large blocks are the most “valuable” resources and can be used for dispatching either large orders or small orders, while small blocks can only be used for small orders. The allocation of shipping capacity depends on the strategic or the planning decision of the shipping manager. It is usually not changed when conducting detailed scheduling. In contrast, reallocation of shipping capacity will not incur extra manual or device cost. Thus, the optimal strategy, i.e., the LM strategy is valuable and easy to be applied in real-life shipping management. Furthermore, the performance of AA strategy is close to that

BSA with ECM P_MS

SS

shipping system, it will certainly incur further cost. Next, the allocation strategy of the fixed shipping capacity was analyzed.

Table 6 Effect of increasing shipping capacity. m

df

SOV denotes the source of variation; df denotes the degrees of freedom; SS denotes the sum of squares; MS denotes the mean square.

The bold values are the best one of PS_MSs obtained by the sequential approach, the integrated approach and BSA with ECM for each instance.

N

SOV

Table 8 Effect of three allocation strategies. n

40

and then becomes slight with increasing K from 3 to 4. The reason is that shipping capacity may be the major bottleneck of the production and shipping system. The marginal benefit of increasing shipping capacity will become less, when it is gradually not the major bottleneck. This finding is useful for setting an optimal shipping capacity, considering the trade-off between the output of the production and shipping system and the cost of increasing shipping capacity. Though increasing shipping capacity is quite a simple and straightforward way to enhance the output of the production and

60

11

m

Strategies

BSA with ECM P_MS

S_MS

PS_MS

CPT (s)

10

AA LM LL

41.0 41.0 40.0

44.0 41.3 72.2

46.0 43.3 74.2

420.7 422.3 410.2

20

AA LM LL

24.0 23.0 23.0

44.0 32.4 72.2

46.0 34.4 74.2

426.2 425.2 424.2

10

AA LM LL

51.0 51.0 54.0

52.3 51.9 75.1

54.3 53.9 77.1

575.4 580.8 562.9

20

AA LM LL

28.0 30.0 29.0

45.0 41.4 75.1

47.0 43.4 77.1

585.2 569.1 559.7

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al. 80

This research can be extended along the following directions. First, for some make-to-order companies, production scheduling may be pulled by shipping scheduling considering due date requirement. In this situation, shipping division may be considered as the leader in the bilevel model and dominates production division. Second, besides the vertical synchronization of one production division and one shipping division, multiple selfish production divisions may exist. This situation will raise another synchronization problem, i.e., horizontal synchronization of multiple production divisions, which is even more complex and deserves more investigation in the future.

n=40, m=10 n=40, m=20 n=60, m=10 n=60, m=20

75 70 65

PS_MS

60 55 50 45

Acknowledgements

40 35

This research was supported by the National Natural Science Foundation of China (No. 51705250, No. 51675442, No. 71701097), Natural Science Foundation of Jiangsu Province (BK20170797), the Fundamental Research Funds for the Central Universities (No. NR2019004). The authors are grateful to the Zhejiang Provincial, Hangzhou Municipal and Lin’an City governments for partial financial supports.

30 AA

LM

LL

Allocation strategies Fig. 7. Effect of three allocation strategies.

of the LM strategy, but much better than the LL strategy. AA strategy could be easily applied in practice considering its simplicity.

References

6. Conclusions

Angelo, J. S., & Barbosa, H. J. (2015). A study on the use of heuristics to solve a bilevel programming problem. International Transactions in Operational Research, 22, 861–882. Calvete, H. I., Galé, C., & Oliveros, M.-J. (2011). Bilevel model for production–distribution planning solved by using ant colony optimization. Computers & Operations Research, 38, 320–327. Chen, Z. L. (2010). Integrated production and outbound distribution scheduling: Review and extensions. Operations Research, 58, 130–148. Chen, J., Huang, G. Q., Luo, H., & Wang, J. (2015). Synchronisation of production scheduling and shipment in an assembly flowshop. International Journal of Production Research, 53, 2787–2802. Chen, Z. L., & Vairaktarakis, G. L. (2005). Integrated scheduling of production and distribution operations. Management Science, 51, 614–628. Cohn, H., & Fielding, M. (1999). Simulated annealing: Searching for an optimal temperature schedule. SIAM Journal on Optimization, 9, 779–802. Colson, B., Marcotte, P., & Savard, G. (2007). An overview of bilevel optimization. Annals of Operations Research, 153, 235–256. Condotta, A., Knust, S., Meier, D., & Shakhlevich, N. V. (2013). Tabu search and lower bounds for a combined production–transportation problem. Computers & Operations Research, 40, 886–900. Delavar, M. R., Hajiaghaei-Keshteli, M., & Molla-Alizadeh-Zavardehi, S. (2010). Genetic algorithms for coordinated scheduling of production and air transportation. Expert Systems with Applications, 37, 8255–8266. Edis, E. B., Oguz, C., & Ozkarahan, I. (2013). Parallel machine scheduling with additional resources: Notation, classification, models and solution methods. European Journal of Operational Research, 230, 449–463. El-Sobky, B., & Abo-Elnaga, Y. (2018). A penalty method with trust-region mechanism for nonlinear bilevel optimization problem. Journal of Computational and Applied Mathematics, 340, 360–374. Hajek, B. E. (1988). Cooling schedules for optimal annealing. Mathematics of operations research, 13, 311–329. Hall, N. G., & Potts, C. N. (2005). The coordination of scheduling and batch deliveries. Annals of Operations Research, 135, 41–64. Heydenreich, B., Müller, R., & Uetz, M. (2007). Games and mechanism design in machine scheduling—an introduction. Production and Operations Management, 16, 437–454. Ishibuchi, H., Misaki, S., & Tanaka, H. (1995). Modified simulated annealing algorithms for the flow shop sequencing problem. European Journal of Operational Research, 81, 388–398. Koulamas, C., & Kyparisis, G. J. (2004). Makespan minimization on uniform parallel machines with release times. European Journal of Operational Research, 157, 262–266. Li, H., & Fang, L. (2012). An evolutionary algorithm for solving bilevel programming problems using duality conditions. Mathematical Problems in Engineering. Liao, C. J., & Lin, C. H. (2003). Makespan minimization for two uniform parallel machines. International Journal of Production Economics, 84, 205–213. Lin, S. W., & Ying, K. C. (2015). A multi-point simulated annealing heuristic for solving multiple objective unrelated parallel machine scheduling problems. International Journal of Production Research, 53, 1065–1076. Lukač, Z., Šorić, K., & Rosenzweig, V. V. (2008). Production planning problem with sequence dependent setups as a bilevel programming problem. European Journal of Operational Research, 187, 1504–1512. Marinakis, Y., Migdalas, A., & Pardalos, P. M. (2007). A new bilevel formulation for the vehicle routing problem and a solution method using a genetic algorithm. Journal of Global Optimization, 38, 555–580. Osborne, M. J., & Rubinstein, A. (1994). A course in game theory. MIT Press. Pinedo, M. L. (2012). Scheduling: Theory, algorithms, and systems. Springer Science & Business Media.

In conclusion, an interdivisional synchronized scheduling problem of production and shipping derived from one of our collaborating manufacturing companies is addressed. The two divisions possess individual decision autonomy, which practically fails to carry out the integrated scheduling. Separate and sequential model results in isolated solutions within the divisions but often leads to inferior overall solutions. A novel bilevel approach is proposed to synchronize production and shipping scheduling jointly. A bilevel scheduling model was built, with each of its level as a mixed integer programming model. A BSA was proposed to generate good quality solutions. Specifically, in the BSA, SA was used to generate tentative upper-level production schedules at each iteration. Based on each tentative production schedule, lower upper shipping scheduling is addressed by an effective heuristic, ECM. An inner SA is also presented as an alternative to ECM. The proposed BSA was compared to the traditional sequential and integrated approaches. Compared to the integrated approach, the POA of BSA is only 1.02 on average for small instances, and BSA offered better solutions for large instances in one-hour computation time time. Compared to the sequential approach, BSA can find much better solutions, especially for large instances. It is noticeable that the computation time of BSA is quite satisfactory and thus is suitable for real-life application. This comparison verified the effectiveness of the proposed approach, demonstrating the value of the proposed approach. Two important findings were identified from the sensitivity analysis of shipping capacity. First, shipping capacity was found to have a significant effect on the overall makespan (PS_MS) as well as the shipping makespan (S_MS). PS_MS decreases with increasing shipping capacity. The decrease is quite sharp first and then increases slightly with increasing shipping capacity. This finding may inspire shipping manager to consider the appropriate amount of shipping capacity in the case of the given makespan. Second, the allocation strategy was investigated for assigning shipping capacity onto blocks, indicating that the LM strategy outperforms the AA and LL strategies in terms of PS_MS. The contributions of this study are summarized. First, a bilevel model was proposed to address the synchronized scheduling of production and shipping when each division having its own decision autonomy. Second, a bilevel-based SA algorithm, considering the leader-follower iterative procedure of the bilevel model, is proposed as the solution algorithm. Third, the benefit of the bilevel approach was compared to the sequential and integrated approaches. Fourth, the best shipping capacity allocation strategy was found via the sensitivity analysis. 12

Computers & Industrial Engineering 137 (2019) 106050

J. Chen, et al. Shi, C. G., Lu, J., Zhang, G. Q., & Zhou, H. (2006). An extended branch and bound algorithm for linear bilevel programming. Applied Mathematics and Computation, 180, 529–537. Stecke, K. E., & Zhao, X. (2007). Production and transportation integration for a make-toorder manufacturing company with a commit-to-delivery business mode. Manufacturing & Service Operations Management, 9, 206–224. Wang, G., & Cheng, T. E. (2000). Parallel machine scheduling with batch delivery costs. International Journal of Production Economics, 68, 177–183. Wang, R., Lai, S., Wu, G., Xing, L., Wang, L., & Ishibuchi, H. (2018). Multi-clustering via evolutionary multi-objective optimization. Information Sciences, 450, 128–140.

Xiang, S., Xing, L. N., Wang, L., & Zou, K. (2019). Comprehensive learning pigeon-inspired optimization with tabu list. Science China Information Sciences, 62, 70208. Yang, R. L. (2000). Convergence of the simulated annealing algorithm for continuous global optimization. Journal of Optimization Theory & Applications, 104, 691–716. Yi, J. H., Xing, L. N., Wang, G. G., Dong, J. Y., Vasilakos, A. V., Alavi, A. H., & Wang, L. (2018). Behavior of crossover operators in NSGA-III for large-scale optimization problems. Information Sciences. Yu, Y., Chu, F., & Chen, H. (2009). A Stackelberg game and its improvement in a VMI system with a manufacturing vendor. European Journal of Operational Research, 192, 929–948.

13