A simulation-based genetic algorithm approach for remanufacturing process planning and scheduling

Accepted Manuscript Title: A simulation-based genetic algorithm approach forremanufacturing process planning and scheduling Author: Rui Zhang S.K. Ong A.Y.C. Nee PII: DOI: Reference:

S1568-4946(15)00562-1 http://dx.doi.org/doi:10.1016/j.asoc.2015.08.051 ASOC 3176

To appear in:

Applied Soft Computing

Received date: Revised date: Accepted date:

25-3-2015 2-8-2015 27-8-2015

Please cite this article as: Rui Zhang, S.K. Ong, A.Y.C. Nee, A simulation-based genetic algorithm approach forremanufacturing process planning and scheduling, Applied Soft Computing Journal (2015), http://dx.doi.org/10.1016/j.asoc.2015.08.051 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Highlights (for review)

ip t

cr us an



M



ed



ce pt



We study the integrated process planning and production scheduling problem for remanufacturing processes under uncertain environments. Two potentially conflicting objective functions are considered simultaneously, one related with process route selection and the other related with production cycle time. A simulation-based genetic algorithm approach is developed to provide approximations to Pareto-optimal solutions. Key parameters of the algorithm have been fine-tuned to save computational time and achieve satisfactory results. Extensive computational experiments and evaluations have been performed on a large set of test instances, and the proposed algorithm is shown to be effective and efficient.

Ac



Page 1 of 34

Crossover

Elitist set (E)

M

Temporary population (TP)

an

us

cr

ip t

*Graphical abstract (for review)

Ac ce pt e

Select a suitable number of solutions for TP

d

Mutation

Update E with the solutions in the first rank

Pareto sorting

Page 2 of 34

*Manuscript

A simulation-based genetic algorithm approach for remanufacturing process planning and scheduling

ip t

Rui Zhanga,b,∗, S. K. Ongb , A. Y. C. Neeb

cr

a School of Management, Xiamen University of Technology, Xiamen 361024, PR China b Department of Mechanical Engineering, National University of Singapore, Singapore 117576, Singapore

us

Abstract

Remanufacturing has attracted growing attention in recent years because of its energy-saving and emission-reduction potential. Process planning and scheduling play important roles in the

an

organization of remanufacturing activities and directly affect the overall performance of a remanufacturing system. However, the existing research on remanufacturing process planning

M

and scheduling is very limited due to the difficulty and complexity brought about by various uncertainties in remanufacturing processes. We address the problem by adopting a simulationbased optimization framework. In the proposed genetic algorithm, a solution represents the

ed

selected process routes for the jobs to be remanufactured, and the quality of a solution is evaluated through Monte Carlo simulation, in which a production schedule is generated following the specified process routes. The studied problem includes two objective functions to be optimized

pt

simultaneously (one concerned with process planning and the other concerned with scheduling), and therefore, Pareto-based optimization principles are applied. The proposed solution approach

Ac ce

is comprehensively tested and is shown to outperform a standard multi-objective optimization algorithm.

Keywords: process planning and scheduling, remanufacturing, genetic algorithm, multi-objective optimization

1. Introduction

In recent decades, the negative environmental impact of manufacturing has come to the attention of the government, and consequently, manufacturers are urged to pay more attention to the ∗ Corresponding author.

Email address: [email protected] (Rui Zhang) Preprint submitted to Applied Soft Computing

August 2, 2015

Page 3 of 34

sustainable aspects of production processes (e.g., energy consumption and pollutant emissions). 5

Under this background, remanufacturing has been recognized as a promising way to achieve the sustainability goals [1, 2]. In a remanufacturing process, a returned worn-out product is com-

ip t

pletely disassembled, and then the usable parts are cleaned and refurbished. Finally, these old

parts (together with some new parts where necessary) are reassembled to produce a like-new

10

cr

product. In many countries, remanufacturing is currently a fast growing industry due to its potentially huge environmental and economical benefits [3, 4]. A lot of research efforts have been made regarding pricing and other decision problems in the remanufacturing industry [5, 6, 7].

us

However, the process planning and production scheduling for remanufacturing is considerably more challenging than for traditional manufacturing. As is well known, a remanufacturing

15

an

system is often characterized by a high level of uncertainties, e.g., the uncertain timing and quantity of returned products, and the highly variable processing routes and processing times for handling the different-quality materials. Because of the inherent uncertainty or randomness

M

in remanufacturing systems, the ordinary process planning and scheduling methods will mostly become ineffective, and thus it is required that new approaches should be designed specifically for the process planning and scheduling in remanufacturing environments. The existing research on process planning and scheduling for remanufacturing has been very

ed

20

limited, probably due to the complicating factors stated above. In terms of remanufacturing process planning, a mixed integer program is established in [8] for the optimization of a reman-

pt

ufacturing process and the evaluation of its economic feasibility. The objective is to maximize the profit of remanufacturing the incoming batch of products under the constraints of machine capacity and inventory levels. A case-based reasoning method is proposed in [9] for remanufac-

Ac ce

25

turing process planning. The method attempts to retrieve and reuse the solutions to previously solved process planning cases by analyzing a local similarity between the pending new case and the past cases stored in a history database. In a recent paper [10], the authors suggest an optimal reconditioning process sequence for the core components based on the conditions of these

30

components. In terms of remanufacturing scheduling, a swarm intelligence algorithm called artificial bee colony has been proposed in [11] for solving a flexible job shop scheduling problem occurring in remanufacturing environments. The jobs must be scheduled in a way that complies with the no-wait constraint. A hybrid meta-heuristic algorithm called discrete harmony search is presented in [12] for solving the re-processing scheduling problem in pump remanufacturing. 2

Page 4 of 34

35

The algorithm aims at minimizing the makespan as well as the earliness and tardiness of product delivery. In addition, the research history of remanufacturing scheduling can be found in the review article [13].

ip t

It should be noted that all the relevant publications treat remanufacturing process planning and remanufacturing scheduling as two separate problems. Process planning and scheduling are typically carried out in a separate and sequential manner (scheduling is performed after process

cr

40

planning has been completed), but in fact, the two decision problems are strongly correlated. If we do not take scheduling related factors into account when doing process planning, the ob-

us

tained process plans may result in poor production performance eventually. For example, when a large number of parts follow the same process route in a job shop style remanufacturing system, the makespan performance will deteriorate considerably because of the increased waiting

an

45

times spent before the bottleneck machines. Therefore, process planning and scheduling should be considered in an integrated manner in order to maximize the performance of remanufacturing

M

systems. The integration of the two functions will lead to significant improvement of production efficiency through the elimination of scheduling conflicts, the reduction of flow time and work 50

in process, and the increase of resource utilization.

ed

Although no existing publications have considered integrated process planning and scheduling (IPPS) for remanufacturing, IPPS has been studied for traditional manufacturing systems for several years. It has been found that the integration of process planning and scheduling al-

55

pt

ways outperforms the situation where the two tasks are conducted separately. For example, a hybrid heuristic method combining genetic algorithm and fuzzy neural network is proposed in

Ac ce

[14]. Meanwhile, multi-objective evolutionary algorithms have been used to solve IPPS with more than one objective functions [15, 16]. A hybrid algorithm based on simulated annealing and tabu search is devised to solve IPPS with stochastic processing times [17]. An ant colony optimization algorithm based on AND/OR graph is described in [18] for the IPPS with makespan

60

criterion. A new heuristic is presented for solving IPPS in reconfigurable manufacturing systems, which are based on reconfigurable machine tools [19]. A new approach consisting of four main modules (process plan selection, scheduling, schedule analysis and process plan modification) is introduced for IPPS [20], which is found to perform better than the conventional hierarchical approach.

65

To the best of our knowledge, this paper is the first attempt to study integrated process plan3

Page 5 of 34

ning and scheduling for remanufacturing. We formulate a stochastic and bi-objective optimization model which captures the features of a typical remanufacturing system. Then, we develop a simulation-based multi-objective genetic algorithm to solve the problem effectively. The key

70

ip t

parameters of the algorithm are fine-tuned based on extensive computational tests. Finally, the

performance of the proposed algorithm is evaluated by comparison with a standard algorithm on

cr

a set of benchmark instances. The remainder of the paper is organized as follows. Section 2

defines the problem under study. Section 3 describes the proposed solution approach in detail. Section 4 presents all the computational results. Section 5 concludes the paper and points out

2. Problem definition

an

75

us

future research directions.

2.1. Description of the problem

M

Since the conditions of the used parts (a.k.a. cores) acquired by a remanufacturer can be widely different, there often exists considerable uncertainty in some characteristics (such as surface quality, dimension and precision) of the jobs waiting to be processed. The high uncertainty on the input side suggests that process route selection is an important decision problem for the

ed

80

scheduling of remanufacturing systems.

We assume that the parts to be remanufactured (i.e., the jobs) can be divided into C categories

pt

based on their quality conditions. For each category c ∈ {1, . . . , C}, the set of feasible process routes are represented by Rc , which means any route r ∈ / Rc cannot be applied to jobs belonging to category c. The candidate routes for each category are determined according to technological ∪ requirements and are supposed to be known in advance. Each route r ∈ C c=1 Rc is associated

Ac ce

85

with a performance score φr (a smaller value is assumed better), which can be regarded as a

measurement of the potential influences of choosing that route (e.g., the electricity consumption or the pollutant emissions). Therefore, we will attempt to minimize the sum of scores among all

90

the selected routes (denoted as Φ). Once a process route has been selected for each job under the feasibility constraints, the subsequent decision problem is equivalent to the job shop scheduling problem (JSSP), and makespan (denoted as MS) is often considered as the performance indicator. Based on the above description, the studied problem can be summarized as follows: given n jobs and an estimated distribution of them (relating to the C categories), we need to decide the

95

process routes for these jobs (final decisions are to be expressed in the form: “xl percent of the 4

Page 6 of 34

jobs in category c should be processed in accordance with route rl ∈ Rc ”), so that both Φ and MS are minimized. The problem is difficult to solve because of two complicating factors:

ip t

• This is a stochastic optimization problem, since we do not know the exact number of jobs in each category during the process planning stage. The reason is that the quantity and 100

quality of used parts are actually exogenous factors not controllable by the remanufac-

cr

turer. Therefore, it is assumed that the category of a job will become clear only in the

remanufacturing stage. When doing process planning, however, we need to rely on esti-

us

mated data regarding the distribution of jobs (e.g., we may be informed that, in the current batch of 100 jobs, about 25% of jobs belong to category 1, about 15% belong to category 105

2, and so forth).

an

• This is a bi-objective optimization problem, since Φ and MS are two potentially conflicting objective functions. If we neglect the scheduling-related factors and only aim at the minimum Φ, we will end up selecting identical or similar process routes for different jobs.

110

M

This will result in a long makespan due to severe job congestions on the bottleneck machines. In order to make full use of the machines’ processing capacities (and thereby to improve MS), the process routes of different jobs should be staggered instead of being

ed

aligned, although this will lead to a larger Φ. Clearly, a trade-off must be made between the two objectives to be minimized.

115

pt

Example 1. A very simple example (in the deterministic sense) is given here to illustrate the rationale of considering process planning and scheduling simultaneously. Suppose there

Ac ce

are two jobs waiting to be processed, each consisting of two operations which need to be processed on two machines. According to the classification scheme, job 1 is in category 1 and job 2 is in category 2. The set of feasible routes for category 1 is R1 = {r1 , r2 } = {(M1 (2), M2 (3)), (M2 (5), M1 (4))}, which means a job in category 1 should be processed in

120

either one of the two ways: (i) process the job first on machine 1 for 2 units of time and then on machine 2 for 3 units of time; (ii) process the job first on machine 2 for 5 units of time and then on machine 1 for 4 units of time. Likewise, the set of feasible routes for category 2 is defined as R2 = {r3 , r4 } = {(M1 (1), M2 (5)), (M2 (4), M1 (2))}. The performance score associated with each route is φr1 = 72, φr2 = 33, φr3 = 55, φr4 = 29 (these scores were generated from the

125

uniform distribution U[1, 100] for this example). Obviously, there are 4 solutions to this problem. If route r1 is selected for job 1 and route 5

Page 7 of 34

r3 is selected for job 2 (solution S1 ), the total score is Φ = 127 and the minimum makespan is MS = 9 (Fig. 1(a)). If r1 is selected for job 1 and r4 is selected for job 2 (solution S2 ), the total score is Φ = 101 and the minimum makespan is MS = 7 (Fig. 1(b)). If r2 is selected for job 1 and r3 is selected for job 2 (solution S3 ), the total score is Φ = 88 and the minimum makespan is

ip t

130

MS = 10 (Fig. 1(c)). If r2 is selected for job 1 and r4 is selected for job 2 (solution S4 ), the total

score is Φ = 62 and the minimum makespan is MS = 11 (Fig. 1(d)). The three Pareto optimal

cr

solutions are: S2 with objective values (101, 7), S3 with objective values (88, 10) and S4 with

discarded. M1

2 0

M1

1 1

3

M2

2

M2

1

1

6

9

5

1 0

9

2

5

2

1

4

7

(b)

10

M2

1 5

1 0

2 9

11

2 5

9

(d)

pt

(c)

ed

M2

6

M1

1 1

4

M

2 0

2

2

0

(a)

M1

1 0

an

135

us

objective values (62, 11). The solution S1 is dominated by the solution S2 , and thus S1 could be

Ac ce

Figure 1: The optimal schedule corresponding to each process route decision

The above example shows that it is necessary to optimize process planning under explicit

consideration of scheduling. In this way, we can obtain a series of Pareto optimal solutions which reflect the possible trade-offs between the process planning objective (represented by Φ)

and the scheduling objective (represented by MS).

140

The difference between this simplified example and the original problem discussed earlier is

that no stochastic factors are involved in this example (i.e., we know exactly the category that each job belongs to). In the original problem that we aim to solve, the number of jobs in each category is not exact but estimated information and is thus subject to random errors. In this case, the objective values of a solution (i.e., route selection expressed in percentages of job numbers)

145

cannot be calculated deterministically but must be evaluated by means of simulation. Therefore, 6

Page 8 of 34

a simulation-based optimization algorithm has to be designed to solve the problem effectively. 2.2. Mathematical model for the deterministic case

ip t

For readers’ convenience, we present a mixed-integer linear programming (MILP) model to formulate the studied problem under the assumption of deterministic job categorization (i.e., the category of each job is assumed to be known).

min

cr

150

f 1 = M S, f 2 = Φ

us

s.t. M S ≥ t jl + p jl n ∑ ∑ Φ= y jr φr

∀ j = 1, . . . , n; l = 1, . . . , m

(1b)

an

j=1 r ∈Rc( j)

(1a)

(1d)

∀ j = 1, . . . , n; l = 1, . . . , m − 1

(1e)

t jl + p jl − M(1 − x jl j ′ l ′ ) ≤ t j ′ l ′ + M z jl j ′ l ′

∀ j, j ′ = 1, . . . , n; l, l ′ = 1, . . . , m

(1f)

t j ′ l ′ + p j ′ l ′ − M x jl j ′ l ′ ≤ t jl + M z jl j ′ l ′ ∑ ∑ z jl j ′ l ′ = w jr j ′ r ′ Vrlr ′ l ′

∀ j, j ′ = 1, . . . , n; l, l ′ = 1, . . . , m

(1g)

∀ j, j ′ = 1, . . . , n; l, l ′ = 1, . . . , m

(1h)

∀ j, j ′ = 1, . . . , n; r ∈ Rc( j) ; r ′ ∈ Rc( j ′ )

(1i)

∀ j, j ′ = 1, . . . , n; r ∈ Rc( j) ; r ′ ∈ Rc( j ′ )

(1j)

∀ j, j ′ = 1, . . . , n; r ∈ Rc( j) ; r ′ ∈ Rc( j ′ )

(1k)

∀ j = 1, . . . , n; l = 1, . . . , m

(1l)

∑

y jr = 1

r ∈Rc( j)

∑

p jl =

πrl y jr

ed

t jl + p jl ≤ t j (l+1)

r ∈Rc( j) r ′ ∈Rc( j ′ )

w jr j ′ r ′ ≤ y j ′ r ′

Ac ce

t jl ≥ 0

pt

w jr j ′ r ′ ≥ y jr + y j ′ r ′ − 1 w jr j ′ r ′ ≤ y jr

(1c)

∀ j = 1, . . . , n; l = 1, . . . , m

M

r ∈Rc( j)

∀ j = 1, . . . , n

∀ j, j ′ , l, l ′ , r, r ′

y jr , x jl j ′ l ′ , z jl j ′ l ′ , w jr j ′ r ′ ∈ {0, 1}

(1m) (1)

The notations used in the above formulation are summarized in Table 1. It is assumed that

the number of jobs waiting to be processed is n and the number of available machines is m, and in accordance with the job shop setting, each job consists of m operations which need to be

155

processed sequentially on each of the machines. Therefore, constraint (1a) defines the makespan (M S), which should be no less than the completion time of any operation. Constraint (1b) defines the total performance score (Φ), which is the second objective function to be minimized. Constraint (1c) states that only one process route could be selected for each job. Constraint (1d) 7

Page 9 of 34

Parameters

ip t cr

Table 1: Notations used in the MILP model

performance score of route r

Rc( j)

set of feasible routes for the category that job j belongs to

πrl

processing time of the l-th operation in route r

Vrlr ′ l ′

a predetermined constant value which equals 0 if the l-th operation in route r and the l ′ -th

an

us

φr

operation in route r ′ are to be processed on the same machine and equals 1 otherwise a positive large number

M

M

Decision variables

the makespan, i.e., the maximum completion time of all jobs

Φ

the total performance score

t jl

starting time of the l-th operation of job j

p jl

processing time of the l-th operation of job j

y jr

a 0-1 decision variable which takes 1 if route r is selected for job j and takes 0 otherwise

x jl j ′ l ′

a 0-1 decision variable which takes 1 if the l-th operation of job j precedes the l ′ -th opera-

pt

ed

MS

tion of job j ′ and takes 0 otherwise (x jl j ′ l ′ is relevant only when z jl j ′ l ′ = 0) a 0-1 decision variable which takes 0 if the l-th operation of job j and the l ′ -th operation of

Ac ce

z jl j ′ l ′

job j ′ are processed on the same machine and takes 1 otherwise

w jr j ′ r ′

a 0-1 decision variable which takes 1 if route r is selected for job j and meanwhile route r ′

is selected for job j ′ and takes 0 otherwise

8

Page 10 of 34

specifies the processing time of the l-th operation of job j, which, by definition, should be equal 160

to the processing time of the l-th operation in the selected route for that job. Constraint (1e) reflects the technological constraint, which requires that consecutive operations of the same job

ip t

should be processed in sequence without any overlap in time. Constraints (1f) and (1g) reflect

the machine capacity constraint (a.k.a. the disjunctive constraint), which is due to the fact that

a machine can process at most one job at a time. The pair of constraints become effective only when z jl j ′ l ′ = 0, which means the l-th operation of job j and the l ′ -th operation of job j ′ must

cr

165

be processed by the same machine. In this case, the order of processing these two operations

us

need to be decided. If x jl j ′ l ′ = 1 (constraint (1g) is then inactivated), the decision is that the l-th operation of job j should be processed first; if x jl j ′ l ′ = 0 (constraint (1f) is then inactivated),

170

an

the decision is that the l ′ -th operation of job j ′ should be processed first. Under either situation, the two operations cannot overlap in time. Constraint (1h) defines the way of obtaining z jl j ′ l ′ : if route r is selected for job j and route r ′ is selected for job j ′ , then z jl j ′ l ′ should take the value of

M

Vrlr ′ l ′ which is known in advance. Constraints (1i)–(1k) define the relationship between y jr , y j ′ r ′ and w jr j ′ r ′ to ensure that w jr j ′ r ′ = 1 if and only if y jr = y j ′ r ′ = 1. It should be noted that the MILP model can only formulate the deterministic version of the integrated optimization problem. When the category of each job is not exactly known, which is

ed

175

the situation in typical remanufacturing systems, we need to resort to simulation-based optimization methods rather than traditional mathematical programming. This is also the motivation for

pt

this study, in which a simulation-based multi-objective evolutionary algorithm will be proposed

180

Ac ce

for remanufacturing process planning and scheduling.

3. The proposed optimization algorithm In this paper, we develop a simulation-based multi-objective genetic algorithm (SMOGA)

for solving the process planning and scheduling problem under uncertainty. The details for the implementation of SMOGA will be described in this section. 3.1. Overview

185

The main steps of the proposed SMOGA can be summarized as follows. The simplified flowchart (omitting initialization and termination) is displayed in Fig. 2.

9

Page 11 of 34

Step 1: Initialize the population. Generate the initial population P0 with a size of PS. Let TP0 = P0 , the generation index g = 0, and the set of elite solutions E = ∅. Step 2: Let Pg = TPg ∪ E. Step 3: Perform crossover for [PS/2] times, each time randomly selecting two individuals from

ip t

190

the current population Pg and executing crossover with probability pc (c.f. Section 3.6). The set of offspring individuals produced is denoted by Pg′ .

cr

Step 4: Perform mutation for each individual in Pg′ with probability pm (c.f. Section 3.6). The mutated population is denoted by Pg′′ .

Step 5: Perform Pareto sorting to the individuals in Pg′′ (c.f. Section 3.4).

us

195

Step 6: Update E with the solutions in the first rank (F1 ) of Pg′′ (c.f. Section 3.5). Remove from Pg′′ those solutions which have just entered E (if any). insert them into TPg+1 (c.f. Section 3.6).

M

Step 8: Let g ← g + 1. If g < GN, then return to Step 2. Otherwise, terminate the algorithm

ed

and output E.

pt

Temporary population (TP)

Ac ce

200

an

Step 7: Using the result of Pareto sorting in Step 5, select (PS − |E|) individuals from Pg′′ and

Crossover

Elitist set (E)

Mutation Update E with the solutions in the first rank

Select a suitable number of solutions for TP

Pareto sorting

Figure 2: The main flowchart of SMOGA

10

Page 12 of 34

3.2. Solution representation and initialization Recall that all the jobs fall into C categories, and for each category c ∈ {1, . . . , C}, the set of feasible process routes is denoted by Rc = {rc,1 , rc,2 , . . . , rc,Uc } (where Uc = |Rc | denotes the number of feasible routes for category c). In addition, a performance score φc,u has been

ip t

205

introduced and associated with each process route rc,u ∈ Rc (u = 1, . . . , Uc ). Based on such

settings, a solution to the considered process planning and scheduling problem can be represented

cr

by x = {xc,u } (c = 1, . . . , C, u = 1, . . . , Uc ), meaning that xc,u percent of the jobs in category

us

210

c should be processed by following the route rc,u . Obviously, the decision variables must satisfy ∑Uc u=1 x c,u = 100 (∀c) and x c,u ≥ 0 (∀c, u). It should be noted at this moment that we are dealing with an integrated optimization prob-

an

lem that incorporates both process planning and scheduling. Although the solution representation scheme described above is only concerned with process planning decisions, the scheduling requirement will be considered when the solutions are being evaluated (see the next subsection). In fact, this can be thought of as a two-layered solution methodology. The selection of process

M

215

routes (crude decision) is implemented in the upper layer, while the scheduling task (detailed decision) is left to the lower layer. This policy well suits the optimization under uncertain envi-

ed

ronments: we can only focus on some crude decisions in a priori optimization, while the detailed decisions should be made after the uncertain factors become known. 220

To initialize the population of solutions, we first produce a “balanced” solution by letting

pt

xc,u = 100/Uc (∀c, u), and then generate the other solutions by imposing some random perturbations on this solution.

Ac ce

For example, if we are faced with a problem with 3 categories (i.e., C = 3) and the number

of feasible routes for each category is 3, 4 and 2, respectively (i.e., U1 = 3, U2 = 4, U3 = 2).

225

Then, the first initial solution can be represented as    x1,1 x1,2 x1,3 33 33       x˜ = x2,1 x2,2 x2,3 x2,4  = 25 25    x3,1 x3,2 50 50

 34 25

  25 . 

3.3. Solution evaluation To evaluate the quality of a solution under randomness, we apply the Monte Carlo simulation 230

method. Three major steps are required in each simulation replication: 11

Page 13 of 34

1. Generate the realized value for the number of jobs belonging to each category c, denoted as n c . 2. Evaluate the total performance score Φ corresponding to the current solution.

ip t

3. Schedule the resulting job shop with a fast algorithm and evaluate the corresponding makespan value M S.

235

cr

According to the principle of simulation, if the objective values obtained from the i-th sim-

us

ulation replication are recorded as Φi and M Si , and a total of N simulation replications have / / ∑N ∑N M Si N ) been performed, then we may rely on the average values (i.e., i=1 Φi N and i=1 as an estimate for the evaluation of the current solution. Also apparently, the precision of such an 240

estimation will improve as N increases. In the following, we will elaborate on each of the three

an

steps. 3.3.1. Realizing the number of jobs in each category

245

M

As mentioned before, in our problem setting, we are only aware of an estimated distribution of the n jobs relating to the C categories. Suppose we have been informed that about n¯ c jobs ∑ belong to category c, with c = 1, . . . , C (keeping C c=1 n¯ c = n). However, the actual number of

ed

jobs in each category is unknown for the moment and may deviate from this a priori information. Therefore, in each simulation replication, we need to generate a realized value for the number of jobs in each category. The procedure is detailed as follows.

250

pt

Step 1: Produce n c (c = 1, . . . , C) by sampling a certain random distribution with mean at n¯ c , e.g., the normal distribution N (n¯ c , σc2 ). The produced value should be rounded to the

255

Ac ce

nearest nonnegative integer. ∑ Step 2: If C c=1 n c = n, then terminate the procedure. ∑C ∑ Step 3: If c=1 n c > n, let δ = C c=1 n c − n. Repeat the following task for δ times: randomly

select c′ such that n c′ > 0, and then let n c′ ← n c′ − 1. ∑C ∑ Step 4: If C c=1 n c . Repeat the following task for δ times: randomly c=1 n c < n, let δ = n − select c′ , and then let n c′ ← n c′ + 1.

In Step 1, the type of distribution should be decided according to specific features of the remanufacturing system under study. After the above procedure is terminated, it is guaranteed ∑ that C c=1 n c is equal to n, and thus n c will be regarded as the number of jobs belonging to 260

category c in the subsequent simulation steps. 12

Page 14 of 34

3.3.2. Calculating Φ First, we need to obtain the exact number of jobs that are to be processed via each process route. Given the current solution x = {xc,u }, the number of jobs required to follow route rc,u

ip t

265

can be calculated by n c,u = ⌊n c × xc,u %⌋ + ηc,u , where ηc,u takes either 0 or 1. Since it must ∑ c be ensured that the rounding procedures result in no violation of U u=1 n c,u = n c , the values of ∑ c ηc,u can be determined in the following way. Within each category c, let ∆c = n c − U u=1 ⌊n c ×

cr

xc,u %⌋, and then sort εc,u = n c × xc,u % − ⌊n c × xc,u %⌋ in decreasing order, and finally assign the ηc,u corresponding to the ∆c largest εc,u to 1 (the remaining ηc,u should be assigned with 0).

us

270

The total performance score (Φ) for the current solution can be computed by Φ = ∑C ∑Uc c=1 u=1 n c,u φc,u .

an

3.3.3. Evaluating M S

Now that the process route for each job has been specified by the current solution, the result-

M

ing decision problem is equivalent to a job shop scheduling problem. Ideally, we should find the optimal schedule under the specified routes and use its makespan as the evaluation (M S) for the 275

current solution x. However, obtaining the optimal schedule is unrealistic both in theory and in

ed

practice. Theoretically, JSSP is N P-hard in the strong sense and probably there exists no polynomial algorithm which can guarantee optimality. Practically, it usually happens that the exact categorization of each job is determined only a moment before the production stage starts1 . Once

280

pt

the job categories have been determined, the extra time that can be used for generating a detailed production schedule is usually very limited (so as to avoid production delays). In this case, a fast

Ac ce

algorithm should be used to obtain a feasible (not necessarily optimal) production schedule. To simulate such a practical need, we use a priority rule based scheduling policy for the

evaluation of M S in our optimization algorithm. The schedule generation method is based on the well-known Giffler-Thompson algorithm [21]. Given a job shop consisting of n jobs and m

285

machines, the evaluation algorithm can be described as follows. Input: A priority rule ρ.

1 Usually, an inspection procedure is performed immediately before the actual remanufacturing operations, and the

aim of the inspection is to determine the conditions of the returned items (and thereby to determine the category that each job should belong to).

13

Page 15 of 34

Step 1: Let Q (1) = O (the set of all operations), R (1) = F (O) = { f 1 , . . . , f n } (the set of the first operation of each job). Set t = 1. Step 2: Find the operation i ∗ = arg mini∈R(t) {ri + pi }, and let k ∗ be the index of the machine

ip t

on which this operation should be processed. Use B (t) to denote all the operations from

290

R (t) which should be processed on machine k ∗ . Step 3: Delete from B (t) the operations that satisfy ri ≥ ri ∗ + pi ∗ .

cr

Step 4: Use the priority rule ρ to identify an operation oˆ from B (t) if currently there are more than one candidates. Schedule operation oˆ on machine k ∗ at the earliest possible time. Step 5: Let Q (t + 1) = Q (t) \{o}, ˆ R (t + 1) = R (t) \{o} ˆ ∪ {suc(o)}, ˆ where suc(o) ˆ denotes the immediate job successor of operation oˆ (if any).

us

295

an

Step 6: If Q (t + 1) ̸= ∅, set t ← t + 1 and go to Step 2. Otherwise, the schedule generation procedure is terminated and the resulting makespan is returned as M S. In the above description, ri represents the earliest possible starting time of operation i (determined by the already scheduled operations), so (ri + pi ) is the earliest possible completion

M

300

time of operation i. Q (t) represents the set of operations yet to be scheduled at iteration t, while R (t) represents the set of ready operations (whose job predecessors have all been scheduled) at

ed

iteration t. In Step 4, the operation set B (t) is also called a conflict set. According to many previous studies, the priority rule ρ is selected as MWR (most work 305

remaining) and SPT (shortest processing time) since these rules are beneficial for achieving a

pt

small makespan in job shops [22]. The MWR rule always picks out the operation associated with the job that has the most work remaining to be processed, while SPT favors the operation

Ac ce

with the minimum processing time. Here, MWR is applied first and SPT is used to break ties if needed.

310

3.4. Pareto sorting

Since we are dealing with a bi-objective optimization problem, Pareto sorting is an important

step in the SMOGA. The aim of Pareto sorting is to differentiate the quality of different individuals in the population. The sorting procedure includes two steps: first partition the solutions into a number of Pareto ranks, and then determine the sorting of solutions within each rank.

315

In the first step, the population is divided into K subsets F1 , F2 , . . . , F K such that: (1) If

14

Page 16 of 34

1 ≤ i < j ≤ K , then for any solution x in F j , there must exist a solution in Fi that dominates2 x; (2) For each i, there exists no dominance relation between any pair of solutions in Fi (neither one dominates the other). Clearly, F1 is the set of Pareto non-dominated solutions in the current

320

ip t

population. The Pareto ranking procedure is detailed in [23].

In the second step, the sorting of solutions within the same Pareto rank Fi depends on the

cr

“density”. The solutions that reside in low-density regions (in the objective space) should come

before the solutions in crowded areas. The aim is to ensure that the solutions ultimately found are diverse and distributed evenly. This mechanism is especially important for simulation-based

325

us

optimization. Because the objective values of each solution are not accurate but subject to noises, a group of solutions that are closely crowded together may have no statistical differences in

an

between.

In order to define a density measure αi for each solution x i in the non-dominated solution set F, we must first define the distance between two solutions x 1 , x 2 ∈ F. The distance between x 1

M

(2)

ed

330

and x 2 in the objective space can be calculated as follows: v )2 u Z ( u∑ f i (x 1 ) − f i (x 2 ) t , D (x 1 , x 2 ) = f imax (F) − f imin (F) i=1

where Z is the number of objectives (Z = 2 for our problem); f imax (F) = max{ f i (x) | x ∈ F} and f imin (F) = min{ f i (x) | x ∈ F} respectively denote the maximum and minimum value of

pt

the i-th objective within the scope of F. The normalization ensures that the objectives under different dimensions can be added.

Then, we sort the solutions in F\{x i } according to their distances to x i . Let D˜ i(1) denote the { } distance between x i and the nearest solution to x i , i.e., D˜ i(1) = min D(x i , x j )| x j ∈ F\{x i } .

Ac ce

335

(2) Likewise, let D˜ i denote the distance between x i and the second nearest solution to x i , · · · ,

and let D˜ i(k) denote the distance between x i and its k-th nearest neighbor solution. Finally, the density of x i is defined as

340

( αi =

λ

1 ∑ ˜ (k) Di λ

)−1 ,

(3)

k=1

2 Referring to our problem, a solution x (with objective values (Φ , M S )) is said to dominate another solution x 1 1 1 2

(with objective values (Φ2 , M S2 )) if either of the following conditions is satisfied: (1) Φ1 < Φ2 and M S1 ≤ M S2 ; (2) M S1 < M S2 and Φ1 ≤ Φ2 .

15

Page 17 of 34

where the truncation length λ is assumed to be 5. If |F| > 5, only the five solutions that are closest to x i will be used to calculate the density of x i . The solutions with smaller αi values are

ip t

given higher priority in the final sorting. 3.5. Elitism 345

Elitism is a strategy of preserving the elite solutions found during the evolutionary process

cr

to ensure that they are not lost. This mechanism can usually help to accelerate the convergence of GAs. In the proposed SMOGA, an elitist archive E is maintained during the evolutionary

us

process. In each generation, the individuals in E join the current population and participate in the crossover and mutation procedures.

The maximum number of solutions in E is set as [PS/4], where PS is the population size.

an

350

The steps of updating E with a solution set X are given as follows: Step 1: Perform the following steps for each solution x i in X :

M

(1.1) If x i is dominated by any solution in E, try the next solution x i+1 . Otherwise, proceed to Step 1.2.

(1.2) Remove from E all the solutions that are dominated by x i , and then let E ← E ∪ {x i }.

ed

355

Step 2: If |E| > [PS/4], sort all the solutions in E in an increasing order of αi , and remove the

pt

(|E| − [PS/4]) solutions with the highest density. 3.6. Genetic operators

The genetic operators in a GA typically include crossover, mutation and selection.

Ac ce

360

For crossover, we first choose cˆ randomly from 1, . . . , C, and then exchange the relevant

parts {xc,u ˆ : u = 1, . . . , Ucˆ } (as a whole) between the two parent solutions. For example, when we consider the same instance as in Section 3.2, the crossover will be performed in the following way (x 1 ⊗ x 2 ⇒ x ′1 , x ′2 ) if cˆ happens to be 2:     30

365

21 70

29

41

42

5

30

32

53

⊗ 10 52

37

10

25

55

10



30

⇒ 10 70

48

29

41

25

55

30



53

10 21 52

37

10

42

5

 32

48

For mutation, we modify the values of two randomly selected variables, say xc,u ˆ 1 and x c,u ˆ 2 (u 1 , u 2 ∈ {1, . . . , Ucˆ }), by letting xc,u ˆ 1 ← x c,u ˆ 1 + 1x and x c,u ˆ 2 ← x c,u ˆ 2 − 1x, where 1x 16

Page 18 of 34

is a random number in [dmin , dmax ] and satisfies 1x ≤ xc,u ˆ 2 . Using the previous example, the mutation will be performed in the following way (x 1 ⇒ x ′′1 ) if it happens that cˆ = 2, u 1 = 3, u 2 = 2 and 1x = 20:



30

21 70

29

41

42

5

 32



30

⇒ 21

30

70

29

41

22

25

 32

30

ip t

370

375

cr

In the SMOGA framework, selection is required after Pareto sorting. (PS − |E|) individuals should be selected for the temporary population TP. In this process, the selection probabilities

ranks i-th is set as

2 (τ + 1 − i) , τ (τ + 1)

(4)

an

P[i] (τ ) =

us

are assigned in a linearly decreasing style, i.e., the probability of selecting the individual that

where τ is the current size of the population. For example, if currently there are 5 individuals 380

in the population, the probability of selecting each individual is 1/3, 4/15, 1/5, 2/15, 1/15,

M

respectively, in the Pareto ranking order. After an individual has been selected, we shall remove it from the population, let τ ← τ − 1 and re-evaluate P[i] (τ ). This process is repeated until a

3.7. Summary 385

ed

desired number of solutions have been selected.

The advantages of the proposed SMOGA as compared to the well-known NSGA-II [23] can

pt

be summarized as follows.

(1) The SMOGA incorporates a specifically-designed solution evaluation method, which is effi-

Ac ce

cient for large-scale simulation optimization.

(2) The SMOGA includes an improved elitism strategy, which dynamically maintains a fixed

390

number of elite solutions throughout the evolutionary process.

(3) The SMOGA introduces an improved Pareto sorting mechanism, which relies on accurately defined distance and density measures.

(4) The genetic operators used in the SMOGA have been modified to better suit the properties of the studied problem.

17

Page 19 of 34

395

4. Computational experiments 4.1. Experimental settings

ip t

The test instances are generated in the following way. The number of jobs and the number of categories (n, C) will be considered at six different levels, i.e., (100, 3), (100, 4), (200, 5), (200, 8), (300, 6) and (300, 10). The estimated number of jobs in each category n¯ c is given by

n¯ c = [n/C] + ξc , where ξc is an integer drawn from the discrete uniform distribution U[−15, 15] ∑ and further adjusted to ensure C c=1 n¯ c = n. The uncertainty in n c is represented by a normal

cr

400

us

distribution N (n¯ c , σc2 ) with σc = 2 + n/100. The number of feasible routes for each category (Uc ) is randomly decided to be one of the four values, i.e., 3, 4, 5 or 6. The number of machines m will be considered at three different levels, i.e., 10, 15 and 20. Each process route is obtained by a random permutation of the m machines and the corresponding processing times are generated

an

405

from the uniform distribution U[10, 90]. The performance score (φ) related with each process

M

route is generated from the uniform distribution U[1, 100].

We have compared our proposed SMOGA to the well-known NSGA-II [23], for two reasons. First, the fundamental structure of the NSGA-II is similar and comparable to our algorithm: both of them rely on the GA framework. Second, the NSGA-II is commonly used as a benchmark

ed

410

algorithm to provide baseline results for the evaluation of newly proposed multi-objective optimization algorithms. Both the proposed SMOGA and the NSGA-II have been implemented

pt

using Visual C++ 2010 and tested on a PC of Intel Core i5-750 2.67GHz CPU, 4GB RAM and Windows 7 operating system. For a fair comparison, the maximum running time is limited to the same level for both algorithms, i.e., n × m/10 seconds for solving an instance with n jobs and m

Ac ce

415

machines.

The recommended parameter values for the SMOGA are as follows: the population size

PS = 40; the maximum number of generations GN = 500 (or determined by the exogenous limit on computational time); the crossover probability pc = 0.8; the mutation probability pm = 0.4.

420

For the NSGA-II, the relevant parameters have been set identically as stated above. In addition, the step size applied in the mutation procedure is defined as dmin = 5 and dmax = 30. The number of simulation replications devoted to evaluating a solution (N ) will be switched from 10 to 50 as the optimization process progresses (the most suitable change point will be discussed later). Note that the final solutions output by each algorithm will be evaluated with 100 simulation

425

replications in order to reflect their true quality with satisfactory accuracy. 18

Page 20 of 34

4.2. The adopted performance measures Different from the single-objective case, we needed to consider the following three aspects

ip t

simultaneously when comparing the quality of multi-objective optimization algorithms: (1) The number of Pareto solutions found. Usually, a larger number of output solutions suggests 430

that the algorithm is of higher quality.

cr

(2) The degree of Pareto optimality of the obtained solutions.

(3) The evenness of the distribution of solutions along the Pareto front. It is desired that the

us

solutions be spread evenly so that the end-user can make a satisfactory choice.

In accordance with the above requirements, we adopted three performance measures: NS, 435

RN and TS. Suppose two different algorithms produce two non-dominated solution sets for the

an

same problem, i.e., SS1 and SS2 . To evaluate the performance of algorithm a (a ∈ {1, 2}), the performance measures are calculated in the following way:

440

M

(1) NS (a) = |SSa |, which is the number of Pareto non-dominated solutions obtained by algorithm a. This is exactly the ONVG metric defined in [24]. (2) RN (a) = SSa∗ /NS (a), where SSa∗ denotes the solutions in SSa that are not dominated by

ed

any solution of the other algorithm. Note that we usually have RN (1) + RN (2) > 1. This performance measure is equivalent to the metric C defined in [25], for example, RN(1) =

pt

445

1 − C(SS2 , SS1 ). √ ∑ ¯ 2 (Tan’s Spacing metric [26]), where D¯ = (3) TS(a) = D1¯ |SS1a | x i ∈SSa (Di − D) ∑ 1 x i ∈SSa Di and Di is the Euclidean distance (in the objective space) between x i and its |SSa |

Ac ce

nearest neighbor solution in SSa . Clearly, a smaller value of TS suggests that the solutions in SSa are distributed more evenly.

4.3. Experiment on the number of simulation replications Under the simulation-based optimization framework, the number of simulation replications

450

used for the evaluation of solutions is a key parameter of the optimization algorithm. Naturally, increasing the number of replications results in more accurate evaluation of candidate solutions, which is beneficial for the search process but will inevitably slow down the entire algorithm. On the contrary, adopting a small number of replications will accelerate the optimization algorithm but meanwhile can mislead the search process because of the inaccuracy in the reported objective 19

Page 21 of 34

455

values. Therefore, a trade-off must be made between accuracy and efficiency, which is always an important perspective in the design of simulation optimization algorithms. In this study, we implement a simple yet effective strategy to control the number of simulation

ip t

replications (i.e., the parameter N in our algorithm). The fundamental idea is quite straightforward. In the beginning period of the optimization process, the algorithm is mainly performing an

“exploration” of the solution space in order to locate promising areas. At this stage, the solutions

cr

460

visited by the algorithm can be widely apart from each other, so a crude estimation is already enough to differentiate the quality difference among the solutions. However, as the optimization

us

progresses, the algorithm will gradually move on to the “exploitation” stage, in which the aim is to find the best solutions from a few small regions of the solution space. In this period, the candidate solutions are often closely situated, and the quality differences between the candidate

an

465

solutions can be very subtle, so the algorithm must rely on more accurate evaluation of the solutions in order to correctly pick out the true best. Based on the above reasoning, the number of

M

simulation replications needs to be increased during the optimization process. In the proposed SMOGA, the parameter N is controlled by a simple strategy, i.e., N is ini470

tially set to 10 and will have a chance to switch to 50 at a certain time point in the optimization

ed

process. So the key question is: when is the most suitable time point to switch N ? To provide some insights into this issue, some preliminary experiments have been conducted. We choose 3 different-scaled instances, and study the influence of the change point on the optimization per-

475

pt

formance of the algorithm. Suppose the whole optimization horizon is evenly divided into 5 intervals (in terms of CPU time), and the switch of N could be implemented at any of the 6

Ac ce

end points. Of course, switching N to 50 at the first point means that the value of N will be 50 throughout the optimization process, and switching N at the final point means that the value of N is always 10 during the optimization process. To reflect the influence of change point position, we focus on the R N performance index

480

since it is most relevant to the quality of optimization. For each of the 3 instances and under each setting of change point, 10 independent runs of the SMOGA are executed under the time limit specified above. To evaluate R N , firstly we run the comparative algorithm NSGA-II 10 times (under a constant setting of N = 50 and no additional limit on computational time), and combine the obtained solutions to produce a single reference set of Pareto non-dominated solutions. Then,

485

the result from each execution of the SMOGA is compared with the reference set and the value 20

Page 22 of 34

of R N is calculated. The averaged R N values are reported in Fig. 3, where the legend (n, C, m) represents the size of each instance under consideration.

ip t

1.0

0.8

RN

cr

0.6

(100, 3, 10) (200, 5, 15)

0.4

us

(300, 6, 20)

0.0 0%

20%

40% 60% Change point

an

0.2

80%

100%

M

Figure 3: The impact of the change point of N on the solution quality

Based on the results, we can observe the variation pattern of R N with respect to the position

490

ed

of change point. When solving the small-scale instance, the best time point to switch N is when about 40% of the overall computational time has elapsed. When faced with the medium- and large-scale instances, however, the best choice is to switch N after about 60%∼80% of the total

pt

time has been consumed. The general trend is that an early switch is preferred for solving smaller instances while a late switch is favored for solving larger instances. The underlying reasons can

495

Ac ce

be intuitively explained as follows. For smaller instances, it does not take long for the algorithm to converge to promising areas in the solution space, and thus in the subsequent search process, the algorithm will mainly focus on exploitation. In this case, it is necessary to increase N at a relatively early time point in order to promote the optimization precision for the upcoming exploitation phase. On the contrary, when the size of the instances increases, the algorithm must spend much longer time on the exploration of the solution space, leaving relatively less

500

time for exploitation. In the exploration phase, it is beneficial to keep N at a low value so that more generations can be evolved (and more solutions can be visited) within the limited computational time. Therefore, when the algorithm is solving large instances, the solution quality will deteriorate considerably if N is switched too early. Another intuitive perspective that may 21

Page 23 of 34

explain the fact mentioned above is the law of large numbers. As the number of jobs increases, 505

the final makespan becomes less sensitive to the variations in {n c }, and therefore there is no need to adopt a large number of simulation replications since fewer replications can already produce

ip t

reliable evaluations of the solutions. This may also explain why the best time point to increase N is postponed when faced with large instances.

510

cr

As suggested by the result of this experiment, we set the time point to switch N at 40% of

total computational time for the instances with 100 jobs, and set the time point at 60% of total

4.4. Experiment on the mutation related parameters

us

computational time for the instances with 200 or 300 jobs.

an

Since mutation is important for increasing the diversity of the population, the related parameters ( pm , dmin and dmax ) may be critical to the algorithm’s performance. We have therefore 515

designed a computational experiment to study the influence of these parameters and find the best

M

setting for them.

The 3 instances used in the previous experiment are still adopted here. We have considered 5 options for the parameter pm (i.e., 0.1, 0.2, 0.3, 0.4 and 0.5), and 6 options for the combina-

configurations. Under each configuration, 10 independent runs of the algorithm are executed for each of the 3 instances. To evaluate R N , the same method as described before is applied. The

\ (dmin , dmax ) pm

pt

averaged R N values are shown in Table 2.

Table 2: Influence of the mutation related parameters

Ac ce

520

ed

tion of (dmin , dmax ) (i.e., (1, 10), (1, 20), (5, 20), (5, 30), (10, 30) and (10, 40)), resulting in 30

0.1

0.2

0.3

0.4

0.5

(1, 10)

0.647

0.674

0.698

0.712

0.708

(1, 20)

0.671

0.695

0.715

0.722

0.709

(5, 20)

0.667

0.715

0.732

0.733

0.725

(5, 30)

0.658

0.708

0.734

0.735

0.722

(10, 30)

0.670

0.720

0.730

0.732

0.725

(10, 40)

0.669

0.710

0.715

0.729

0.729

As the results suggest, the best combination of the parameters is (dmin , dmax ) = (5, 30) and 22

Page 24 of 34

pm = 0.4. 525

Meanwhile, we can make some remarks regarding the influence of these parameters. A comparison of the results obtained under different pm suggests that the mutation procedure plays

ip t

a significant role in the proposed algorithm. The best value for pm in the SMOGA is somewhat

larger than the common settings in a typical genetic algorithm. The emphasis on mutation is

530

cr

caused by two factors. First, in the multi-objective optimization background, population diversity is always important since we aim at a group of output solutions rather than a single one. Mutation is encouraged since it increases diversity. Second, in the proposed SMOGA, the function of

us

crossover in modifying a solution is limited since it exchanges decisions for a certain category simply as a whole and does not go into sufficient details. In this respect, the mutation operator

535

an

provides a complementary mechanism by modifying the basic elements of each solution and thus helps to diversify the search. However, the experiment results also show a degradation of solution quality when pm exceeds the optimal value. This is because excessive diversification

M

may introduce misleading information and impair the convergence of the algorithm. In addition to a proper setting of pm , the amplitude of mutation (characterized by dmin and dmax ) must also be controlled into a reasonable range. Note that dmin and dmax correspond to changes in percentage points, so too small values may result in ineffective mutation because

ed

540

it may so happen that, after decoding the solution, the number of jobs following each route does not change at all. On the other hand, the values of the two parameters should not be too

pt

large either. Otherwise, the mutation may bring about major changes to the current solution and

545

Ac ce

thereby destroy some good patterns existing in the solution. 4.5. Main computational results We systematically compare the two algorithms using randomly-generated test instances. For

each problem specification (n, C, m), 5 different instances have been generated. Since we have considered 6 levels for (n, C) and 3 levels for m, a total of 6 × 3 × 5 = 90 instances are involved

in our experiment. These instances are grouped according to the number of jobs (n), so each

550

group contains 30 instances. For each instance, both the SMOGA and NSGA-II are executed 10 times independently so that the RN index can be calculated on a one-to-one basis. NS, RN and TS respectively denote the average values of the three indices resulted from the 10 runs of each algorithm. Tables 3–5 show the computational results in terms of the three performance measures. 23

Page 25 of 34

Table 3: Computational results for the instances with n = 100

NS

(100, 4, 10)

Average

SMOGA

9.06

0.92

0.16

1.34

11.60

10.02

0.73

0.34

0.99

3

8.70

8.46

0.88

0.23

0.98

4

12.41

10.80

0.91

0.20

5

12.64

10.80

0.96

0.16

6

10.53

9.57

0.75

7

8.64

8.59

0.75

8

11.57

8.67

0.79

9

8.30

6.27

0.72

10

8.04

6.46

0.79

1.80 1.34

1.04

1.01

1.07

1.23

1.24

1.26

1.56

0.37

1.26

1.59

0.31

1.10

1.11

0.40

1.33

1.79

0.27

1.12

1.42

an

us

0.32

NSGA-II

cr

10.18

2

11

8.95

8.55

0.76

0.34

1.40

1.76

12

7.56

5.44

0.73

0.37

1.11

1.39

13

7.90

6.90

0.72

0.41

1.22

1.49

14

11.29

8.66

0.90

0.15

1.14

1.41

15

8.54

8.74

0.96

0.15

1.22

1.22

16

10.82

8.54

0.72

0.38

1.24

1.34

17

10.82

7.01

0.95

0.15

1.40

1.40

18

9.46

7.04

0.91

0.20

1.37

1.77

19

12.35

10.12

0.83

0.23

1.11

1.42

13.26

8.70

0.79

0.26

1.27

1.36

12.96

12.68

0.85

0.24

1.19

1.33

22

10.19

7.98

0.75

0.37

1.00

1.35

23

10.05

7.16

0.76

0.32

1.24

1.53

24

12.42

12.56

0.79

0.31

1.32

1.70

25

10.22

9.43

0.92

0.19

1.26

1.62

26

10.27

8.09

0.91

0.17

0.95

1.12

27

8.54

5.23

0.87

0.22

1.09

1.12

28

13.46

11.93

0.88

0.22

1.41

1.55

29

13.20

11.83

0.78

0.34

1.41

1.65

30

11.02

7.90

0.93

0.15

1.10

1.26

10.53

8.77

0.83

0.26

1.20

1.42

21

Ac ce (100, 4, 20)

NSGA-II

1

20 (100, 4, 15)

SMOGA

ed

(100, 3, 20)

NSGA-II

pt

(100, 3, 15)

TS

ip t

SMOGA (100, 3, 10)

RN

No.

M

Instance size

24

Page 26 of 34

Table 4: Computational results for the instances with n = 200

NS

(200, 8, 10)

Average

SMOGA

8.03

0.93

0.15

1.25

0.90

0.21

1.23

33

11.31

10.71

0.88

0.19

1.05

34

10.36

7.12

0.90

0.17

35

12.52

12.82

0.91

0.25

36

11.67

8.09

0.86

37

11.76

11.02

0.83

38

11.29

9.37

0.94

39

13.97

10.92

0.92

40

9.73

8.40

0.93

41

11.86

10.88

42

12.10

8.43

43

12.85

8.80

44

13.80

11.51

45

9.92

8.36

0.27

NSGA-II

ip t

8.12

11.68

1.33 1.51 1.09

1.23

1.51

1.04

1.08

1.08

1.27

0.27

0.92

0.98

0.16

1.20

1.20

0.24

0.99

1.02

0.20

1.18

1.33

0.85

0.26

0.91

1.05

0.93

0.22

0.93

1.09

0.94

0.17

0.92

1.13

0.81

0.30

1.07

1.19

0.91

0.20

1.14

1.16

46

11.73

8.15

0.80

0.31

0.93

0.91

47

12.85

12.40

0.83

0.31

0.94

1.20

48

11.29

7.78

0.91

0.18

1.24

1.25

49

12.24

12.51

0.97

0.14

1.08

1.14

13.61

10.72

0.92

0.19

1.16

1.48

14.94

14.77

0.84

0.24

1.17

1.44

52

12.56

9.28

0.81

0.31

1.01

1.10

53

14.68

10.07

0.79

0.30

0.97

1.17

54

13.35

10.45

0.96

0.18

1.17

1.44

55

11.80

9.63

0.93

0.21

1.23

1.31

56

12.57

10.86

0.80

0.34

1.11

1.35

57

12.09

8.32

0.84

0.25

1.13

1.19

58

13.37

11.25

0.88

0.22

0.97

1.02

59

12.19

9.99

0.97

0.17

1.00

1.07

60

14.22

10.93

0.91

0.19

1.20

1.38

12.32

9.99

0.89

0.23

1.08

1.21

51

Ac ce (200, 8, 20)

NSGA-II

cr

11.44

32

SMOGA

us

31

50 (200, 8, 15)

TS

an

(200, 5, 20)

NSGA-II

M

(200, 5, 15)

SMOGA

ed

(200, 5, 10)

RN

No.

pt

Instance size

25

Page 27 of 34

Table 5: Computational results for the instances with n = 300

NS

(300, 10, 10)

NSGA-II

SMOGA

61

16.76

15.35

0.95

0.13

0.84

62

16.64

14.29

0.85

0.22

1.12

63

12.69

10.48

0.93

0.12

1.18

64

12.61

11.10

0.83

0.23

65

17.08

12.12

0.84

0.22

66

17.27

15.85

0.80

0.26

67

12.68

9.30

0.88

68

12.81

10.22

0.89

69

16.85

13.16

0.79

70

13.96

11.98

71

17.09

72

13.34

73

12.78

11.31

74

15.75

13.01

75

13.17

12.19

76

18.05

77

16.03

78 79

Average

1.13

cr

1.60

0.92

1.15

0.98

1.03

1.46

0.91

1.05

0.19

0.92

1.16

0.26

0.90

1.07

0.80

0.27

0.98

1.14

14.32

0.81

0.24

1.13

1.15

11.74

0.85

0.22

0.90

0.91

0.89

0.17

1.09

1.37

0.93

0.12

0.85

1.09

0.90

0.16

1.16

1.59

16.04

0.90

0.16

1.24

1.37

12.06

0.86

0.22

0.87

1.03

15.78

12.79

0.85

0.23

1.03

1.06

17.61

13.78

0.80

0.26

0.99

1.23

15.13

11.66

0.95

0.11

1.04

1.30

18.72

15.74

0.87

0.19

1.02

1.35

82

14.59

12.50

0.83

0.23

1.11

1.38

83

14.19

10.14

0.89

0.17

1.18

1.40

84

18.28

13.97

0.94

0.13

1.12

1.37

85

15.74

12.62

0.84

0.24

1.04

1.29

86

17.49

14.32

0.83

0.24

1.06

1.41

87

18.32

16.76

0.93

0.15

1.22

1.44

88

15.70

13.59

0.91

0.15

1.14

1.33

89

14.43

11.11

0.93

0.15

1.19

1.23

90

16.74

13.29

0.93

0.16

0.94

1.14

15.61

12.89

0.87

0.19

1.04

1.24

81

us

1.18

Ac ce (300, 10, 20)

1.08

0.18

80 (300, 10, 15)

NSGA-II

ip t

SMOGA

an

(300, 6, 20)

NSGA-II

M

(300, 6, 15)

TS

SMOGA

ed

(300, 6, 10)

RN

No.

pt

Instance size

26

Page 28 of 34

555

Based on the computational results, the following remarks could be made: (1) Focusing on N S, we can conclude that, on average, the SMOGA is able to find more Pareto

ip t

solutions than the NSGA-II. Meanwhile, we notice that the number of Pareto solutions found by each algorithm increases considerably as the problem size grows. When the number of jobs and the number of categories increase, the ways of utilizing the process routes become

more diverse, and thus there will potentially be more opportunities to compromise between the two objectives, resulting in a larger number of Pareto solutions.

cr

560

us

(2) Focusing on R N , we can see that the solutions found by the SMOGA are of much higher quality than the solutions obtained by the NSGA-II. For example, with respect to the 200-job instances, 89% of the SMOGA solutions are not dominated by any of the NSGA-II solutions, while 77% of the NSGA-II solutions are dominated by at least one of the SMOGA solutions

an

565

(in the average sense). The advantage of SMOGA can be attributed to the improved search mechanisms, especially the elitism strategy which preserves the highest-quality solutions and

M

at the same time maintains diversity of the population. By contrast, the NSGA-II does not impose a limit on the number of elite solutions and thus may lead to premature convergence. 570

(3) Focusing on T S, we can judge that the solutions obtained by the SMOGA are distributed

ed

more evenly than the solutions found by the NSGA-II. This certainly verifies the effectiveness of the Pareto sorting and the selection procedure, which give priority to the solutions in

pt

less crowded areas while not neglecting population diversity. 4.6. Experiment on the computational time performance Since the algorithm’s execution speed is crucial for avoiding production delays in practical

Ac ce

575

application, we have designed the following experiments to examine the time performance of the proposed algorithm.

First, we remove the computational time limit and allow both algorithms to complete the

desired number of generations (i.e., G N = 500). Nine instances with different sizes have been

580

chosen for this experiment, and the average computational time consumed by each algorithm is recorded in Fig. 4. As suggested by the figure, the time complexity of SMOGA is generally comparable with that of NSGA-II, mainly because both algorithms are built upon the same evolutionary optimization framework. Meanwhile, SMOGA is more time-efficient for solving larger-scale instances since it adopts an adaptive solution evaluation method which helps to bring down the 27

Page 29 of 34

required number of simulation replications.

ip t

900 800 700 600 500 400 300 200 100 0

cr

SMOGA

us

NSGA-II

an

Computational time (sec)

585

Figure 4: Comparison of computational time under G N = 500

M

Second, we will test the robustness of the proposed algorithm with respect to computational time limit. As stated earlier, in the previous experiments the computational time limit has been set as n × m/10 seconds for solving an n-job, m-machine instance. Accordingly, for example,

590

ed

the available time for solving a (300, 10, 20) instance is 10 minutes. However, if we consider a dynamic production environment in which the algorithm must respond very quickly to possible changes of job information, it is likely that even 10 minutes’ time is not affordable. In this experi-

pt

ment, we focus on the five largest instances with size (300, 10, 20) and compare the performance of the two algorithms under modified computational time limits. In particular, the time limit is

595

Ac ce

varied from 2.5 minutes to 15 minutes, and the average values of RN are reported in Fig. 5. Here the RN for each algorithm is evaluated based on the same reference set, which is obtained by combining the solutions from 10 runs of NSGA-II with no additional time limit. For each algorithm, according to the figure, the solution quality will naturally deteriorate

as the time budget becomes tighter. Comparatively, the proposed SMOGA has quite consistent performance under different time limits. The solution quality achieved by SMOGA under such

600

a tight limit as 2.5 minutes is still good enough (on average, as many as 67.04% of the solutions found by SMOGA are not dominated by any solution in the reference set). This is certainly an advantage of the proposed algorithm especially for real application. By contrast, NSGA-II is more sensitive to the execution time constraint, and the solution quality decreases considerably 28

Page 30 of 34

1.0 0.8

RN

ip t

0.6

SMOGA

0.4

cr

NSGA-II

0.2

2.5

5

us

0.0

7.5 10 12.5 Computational time limit (min)

15

an

Figure 5: Performance of the algorithms under different computational time limits

when the algorithm is not run for sufficiently long time. Based on such a comparison, we can judge that SMOGA not only finds better solutions but also converges much faster than NSGA-II.

M

605

ed

5. Conclusion

In this paper, we consider the integrated process planning and scheduling problem for a typical remanufacturing system. Since the quality conditions of the returned items are unknown

610

pt

in advance, the classification of jobs is uncertain. To cope with the uncertainty, a simulationbased optimization framework is adopted. In the proposed algorithm, the solution representation

Ac ce

takes care of the process route decisions, while the detailed production schedule is generated in the simulation procedure. Such a two-layered optimization scheme is beneficial for maintaining the flexibility of final decisions. In addition, the studied problem involves two objectives to be minimized simultaneously, and therefore Pareto-based optimization principles are applied to

615

identify a group of potential tradeoff solutions. To evaluate the effectiveness of the proposed algorithm, a series of computational experiments have been conducted. The impact of some key parameters have been investigated. The computational results reveal that the proposed algorithm outperforms the well-known NSGA-II, a benchmark algorithm for multi-objective optimization. The proposed approach is potentially applicable to the process planning and scheduling of

620

real-world remanufacturing systems. Before applying the proposed algorithm, the production manager must first obtain some information about the condition of used parts in the current batch. 29

Page 31 of 34

This could be achieved by a crude inspection procedure or a regression analysis on the selected features using historical data. Based on the obtained knowledge, the distribution of jobs could be estimated. The proposed SMOGA will then be used to determine the process route for each job. Finally, the production manager may resort to an independent job shop scheduling software

ip t

625

to further optimize the ultimate production schedule under the determined process routes.

cr

The future research will be carried out in the following aspects:

• Efforts will be made to further promote the computational efficiency of the simulation-

630

us

based optimization framework. In the current work, each solution is evaluated with an identical number of simulation replications, which may lead to inefficiency because some apparently inferior solutions can be excluded from consideration after only a few simula-

an

tion replications. To utilize the limited number of simulation replications more efficiently, a computing budget allocation mechanism similar to the one described in [27] might be designed under the multi-objective optimization environment.

• Additional types of uncertainties in remanufacturing systems will be modelled and con-

M

635

sidered. In the current work, we mainly focus on the uncertainty on the input side, i.e., the categorization of the input jobs. Meanwhile, uncertain factors can also exist within the re-

ed

manufacturing processes. For example, the processing times are usually uncertain because the operation time is partly determined by the damage condition of the part being pro640

cessed, which is unknown in advance. Taking such kind of uncertainty into consideration

pt

means that we need to handle a stochastic scheduling subproblem inside the simulation procedure, which certainly increases the complexity of the entire problem.

Ac ce

• The algorithm may be extended to solve other decision problems related with the remanufacturing process. For example, given the desired working time, the algorithm can be

645

utilized by an outer-layer optimization procedure to find the optimal number of jobs (i.e., batch size) for the current remanufacturing system. In addition, the algorithm can be adapted for re-scheduling purpose if we design a strategy that is able to quickly refocus on feasible schedules when unexpected perturbations occur in the remanufacturing system.

Acknowledgement 650

The paper is supported by the Natural Science Foundation of China under Grant Nos. 61473141 and 61273233. The authors are grateful to the three anonymous reviewers for 30

Page 32 of 34

their helpful comments.

ip t

References [1] J. Li, W. Du, F. Yang, G. Hua, The carbon subsidy analysis in remanufacturing closed-loop supply chain, Sustain655

ability 6 (6) (2014) 3861–3877.

[2] W. Shi, K. J. Min, Product remanufacturing: A real options approach, IEEE Transactions on Engineering Manage-

cr

ment 61 (2) (2014) 237–250.

[3] V. V. Agrawal, A. Atasu, K. Van Ittersum, Remanufacturing, third-party competition, and consumers’ perceived

us

value of new products, Management Science 61 (1) (2015) 60–72. 660

[4] Q. Zhu, J. Sarkis, K.-h. Lai, Supply chain-based barriers for truck-engine remanufacturing in China, Transportation Research Part E: Logistics and Transportation Review 68 (2014) 103–117.

Operations Management 36 (2015) 130–146.

an

[5] J. D. Abbey, J. D. Blackburn, V. D. R. Guide, Optimal pricing for new and remanufactured products, Journal of

[6] Z.-Z. Zhang, Z.-J. Wang, L.-W. Liu, Retail services and pricing decisions in a closed-loop supply chain with 665

remanufacturing, Sustainability 7 (3) (2015) 2373–2396.

M

[7] T. Schulz, G. Voigt, A flexibly structured lot sizing heuristic for a static remanufacturing system, Omega - The International Journal of Management Science 44 (2014) 21–31.

[8] S. Kernbaum, S. Heyer, S. Chiotellis, G. Seliger, Process planning for IT-equipment remanufacturing, CIRP Journal of Manufacturing Science and Technology 2 (1) (2009) 13–20.

[9] F. Zhou, Z. Jiang, H. Zhang, Y. Wang, A case-based reasoning method for remanufacturing process planning,

ed

670

Discrete Dynamics in Nature and Societydoi:10.1155/2014/168631. [10] S. Tsang Mang Kin, S. K. Ong, A. Y. C. Nee, Remanufacturing process planning, Procedia CIRP 15 (2014) 189– 194.

675

pt

[11] S. Sundar, P. N. Suganthan, T. J. Chua, A swarm intelligence approach to flexible job-shop scheduling problem with no-wait constraint in remanufacturing, Lecture Notes in Computer Science 7895 (2013) 593–602.

Ac ce

[12] K. Gao, P. N. Suganthan, T. J. Chua, T. X. Cai, C. S. Chong, Hybrid discrete harmony search algorithm for scheduling re-processing problem in remanufacturing, in: Proceedings of the 15th annual conference on genetic and evolutionary computation, ACM, 2013, pp. 1261–1268.

[13] S. D. Morgan, R. J. Gagnon, A systematic literature review of remanufacturing scheduling, International Journal

680

of Production Research 51 (16) (2013) 4853–4879.

[14] A. Seker, S. Erol, R. Botsali, A neuro-fuzzy model for a new hybrid integrated process planning and scheduling system, Expert Systems with Applications 40 (13) (2013) 5341–5351.

[15] L. Zhang, L. Gao, X. Li, A hybrid genetic algorithm and tabu search for a multi-objective dynamic job shop scheduling problem, International Journal of Production Research 51 (12) (2013) 3516–3531.

685

[16] P. Mohapatra, A. Nayak, S. K. Kumar, M. K. Tiwari, Multi-objective process planning and scheduling using controlled elitist non-dominated sorting genetic algorithm, International Journal of Production Research 53 (6) (2015) 1712–1735. [17] M. Haddadzade, M. R. Razfar, M. H. F. Zarandi, Integration of process planning and job shop scheduling with

31

Page 33 of 34

stochastic processing time, The International Journal of Advanced Manufacturing Technology 71 (1-4) (2014) 690

241–252. [18] J. Wang, X. Fan, C. Zhang, S. Wan, A graph-based ant colony optimization approach for integrated process planning

ip t

and scheduling, Chinese Journal of Chemical Engineering 22 (7) (2014) 748–753. [19] A. Bensmaine, M. Dahane, L. Benyoucef, A new heuristic for integrated process planning and scheduling in reconfigurable manufacturing systems, International Journal of Production Research 52 (12) (2014) 3583–3594. 695

[20] R. K. Phanden, A. Jain, R. Verma, An approach for integration of process planning and scheduling, International

cr

Journal of Computer Integrated Manufacturing 26 (4) (2013) 284–302.

[21] B. Giffler, G. L. Thompson, Algorithms for solving production-scheduling problems, Operations Research 8 (4) (1960) 487–503.

us

[22] K. R. Baker, D. Trietsch, Principles of Sequencing and Scheduling, John Wiley & Sons, 2009. 700

[23] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6 (2) (2002) 182–197.

an

[24] D. A. Van Veldhuizen, G. B. Lamont, On measuring multiobjective evolutionary algorithm performance, in: Proceedings of the IEEE Congress on Evolutionary Computation, Vol. 1, La Jolla, CA, 2000, pp. 204–211. [25] E. Zitzler, L. Thiele, Multiobjective evolutionary algorithms: A comparative case study and the strength Pareto approach, IEEE Transactions on Evolutionary Computation 3 (4) (1999) 257–271.

M

[26] K. C. Tan, C. K. Goh, Y. J. Yang, T. H. Lee, Evolving better population distribution and exploration in evolutionary multi-objective optimization, European Journal of Operational Research 171 (2) (2006) 463–495. [27] C.-H. Chen, E. Y¨ucesan, L. Dai, H.-C. Chen, Optimal budget allocation for discrete-event simulation experiments,

pt

ed

IIE Transactions 42 (1) (2009) 60–70.

Ac ce

705

32

Page 34 of 34