Improving configuration of complex production lines via simulation-based optimization

Computers & Industrial Engineering 109 (2017) 295–312

Contents lists available at ScienceDirect

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

Improving configuration of complex production lines via simulation-based optimization Mustafa Fatih Yegul b,1, Fatih Safa Erenay b,⇑, Soeren Striepe a,2, Mustafa Yavuz a,3 a b

Department of Mechanical and Mechatronics Engineering, University of Waterloo, Waterloo, ON, Canada Department of Management Sciences, University of Waterloo, Waterloo, ON, Canada

a r t i c l e

i n f o

Article history: Received 18 June 2016 Received in revised form 21 March 2017 Accepted 12 April 2017 Available online 17 April 2017 Keywords: Production line configuration Simulation-based optimization Ant-colony optimization Response-surface method Simulated-annealing

a b s t r a c t Optimizing the configuration of a complex production line is an NP-hard problem in various machine settings. Solving real-life-size instances of this problem becomes a more common challenge because the current trend of reshoring induces multi-national firms to transfer manufacturing facilities from workforceintensive to capital-intensive production environments which usually require re-configuration of the transferred manufacturing systems according to the availability of better machinery in a capitalintensive environment. This paper focuses on the problem of optimizing production line configuration and proposes several simulation-based optimization approaches based on myopic search, ant-colony, simulated annealing, and response-surface methodologies. We investigate the relative performances of these proposed algorithms on a real-life manufacturing system transfer case in automotive industry according to solution quality and computation-time metrics under different parameter scenarios. Thus, our numerical results may guide the decision makers in choosing a suitable solution approach for this problem depending on the problem size and time availability. Our results also illustrate that antcolony optimization, a methodology not widely applied in simulation-based optimization, provides high-solution quality for this problem when matched-up with a myopic search to find a good initial solution. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Determining the best manufacturing system configuration (MSC) may significantly improve reliability, product quality, capacity scalability, and costs in production facilities; hence, it is crucial for profitability. For instance, (Koren, Hu, & Weber, 1998) report that improving MSC may lead to up to 22% savings in costs and up to 100% improvement in productivity. However, optimizing MSC is a complex problem which requires specifying many factors including facility layout, production flow, station structures,

Abbreviations: MSC, Manufacturing System Configuration; RTMS, Reverse Transfer of Manufacturing System; GA, Genetic Algorithm; SA, Simulated Annealing; TS, Tabu Search; ANN, Artificial Neural Networks; WIP, Work-In-Process; ACO, Ant Colony Optimization; GS, Greedy Search; FGS, Fast Greedy Search. ⇑ Corresponding author at: 200 University Avenue West, CPH-4323, Waterloo, ON N2L3G1, Canada. E-mail addresses: [email protected] (M.F. Yegul), [email protected] (F.S. Erenay), [email protected] (S. Striepe), [email protected] (M. Yavuz). 1 Address: 200 University Avenue West, E5-3011, Waterloo, ON N2L3G1, Canada. 2 Address: 18579 Centre Street, Mt. Albert, ON L0G1M0, Canada. 3 Address: 200 University Avenue West, E5-3011, Waterloo, ON N2L3G1, Canada. http://dx.doi.org/10.1016/j.cie.2017.04.019 0360-8352/Ó 2017 Elsevier Ltd. All rights reserved.

configuration of the machinery (machinery type and process capacity, buffer size, job allocation mechanism, etc.), number of workers, and skill sets (Wang & Chatwin, 2005). In productionline-based assembly systems (e.g., flow-shops), the problem complexity of optimizing MSC may quickly increase with the number of considered stations (Jeong & Kim, 2000). For example, the optimal buffer allocation in flow-shops, that considers only one of the aspects of MSC listed above, is reported to be an NP-hard combinatorial optimization problem (Huang, Chang, & Chou, 2002). Optimal MSC problems generally arise when a facility is first built or upgraded. Recently, challenging MSC problems have arisen when manufacturing systems are transferred from workforceoriented to automation-oriented environments as a re-shoring practice. In such cases, the transferred manufacturing system is reconfigured based on the machinery availability of the automation-oriented environment. We refer to such a transfer as ‘reverse transfer of manufacturing system’ (RTMS) following the ‘reverse transfer of technology’ concept of (Elshout, 1995). As the practice of relocating manufacturing services back to Europe and North America draws more attention (Ellram, 2013), RTMS is becoming a more common practice especially for

296

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

high-value-added products prototyped in workforce-oriented production cultures. Therefore, RTMS related manufacturing system configuration problems may draw more attention in the future. Optimizing MSC in case of RTMS could be an especially complex problem because (1) it may need to be solved in real-life-size instances as the transferred manufacturing system already has a particular work-flow that cannot be altered; (2) better machinery availability in the automation-oriented environment increases the number of options to be considered (e.g., replacing workers with machinery, upgrading machinery, determining the working pace for the altered components, etc.). Because close form formulations of the objective functions for real-life-size MSC problems are hard to attain if not impossible (Huang et al., 2002), discreteevent simulation could be an appropriate tool to model RTMSrelated MSC problems. Therefore, efficient simulation-based optimization approaches are needed to solve such MSC problems. In this context, this paper proposes eight simulation-based optimization algorithms to derive an efficient configuration (i.e., a good and feasible combination of configuration parameters) for production-line-based manufacturing systems. We test the performances of the proposed metaheuristics on a complex real-life automotive industry case for which an approximately optimal MSC needs to be derived for a prospective production line of a high-value-added product to be built in a capital-intensive country based on a similar existing workforce-oriented manufacturing system. We develop a stochastic discrete-event simulation model mimicking the processes in the transferred production system for any given combination of decision variables assuming that the demand for end-products is constant, and the stations and process flow in the existing system are preserved. The proposed algorithms derive approximately optimal (i.e., a good feasible solution) configuration parameters of the prospective production line including number of machines, number of workers, processing rates, and buffer sizes to maximize the expected profit. Most of the proposed algorithms are based on conventional metaheuristics commonly used in simulation-based optimization such as greedy search, response surface method, and simulated annealing. We also propose an algorithm based on ant-colony optimization, which has not been widely applied to simulation-based optimization, due to its promising performance in some production line optimization problems (Chehade, Amodeo, & Yalaoui, 2011; Mohan & Baskaran, 2012; Rossi & Lanzetta, 2013). In addition, we also examine a couple of hybrid algorithms to measure the benefit in combining some of the proposed algorithms to enhance their performance. Lastly, we compare the performances of the proposed algorithms with an optimizer available in simulation packages (OptQuest) illustrating to what extend an easy-toimplement commercial simulation-based optimization tool can satisfy decision makers in terms of solution quality and speed. The appropriate parameter settings for these proposed algorithms are specified via initial testing on a test-bed problem with different parameter instances. The test-bed problem represents an artificial smaller-scale production line which has similar characteristics (production flow, machinery environment, etc.) with the aforementioned real-life automotive industry case. The production line in the test-bed problem can be simulated much faster and has less decision variables compared to the one in the real-life case, i.e., the best production line configuration parameters of the test-bed problem can be derived via total enumeration of all possible combinations of discrete decision variables. Therefore, we also use the test-bed problem to evaluate the performances of the proposed algorithms compared to the optimal configuration parameters (i.e., the best solution found by total enumeration) in small problem instances. There are two main contributions of this paper: (1) We propose and test several simulation-based optimization approaches based

on metaheuristics for optimizing configuration of complex production lines, and illustrate their relative strengths under different parameter scenarios in terms of solution quality and speed. (2) This is the first study that considers the problem of optimizing production line configuration in an RTMS case. The considered real-life production system (i.e., an asynchronous and inhomogeneous production line with series-parallel flow including reworking, scraping, and product splitting/merging) is very complex compared to those analyzed in the literature; therefore, appropriately represents the challenge associated with the production line configuration problems in case of RTMS. In short, we provide an interesting case study analysis on eightway comparison of conventional simulation-based optimization methods applied to a challenging production line configuration problem in an RTMS case using real data. Given the current trends (e.g., reshoring) and political developments in the North American manufacturing sector, more companies are likely to face the challenge of optimizing configurations of their new or old production lines in an RTMS setting. Therefore, the case study provided in this paper is of interest for both practitioners and academics since it informs about the potential of simulation-based optimization to address such challenges and provide guidance in selection of an appropriate simulation-based optimization method.

2. Details about the real-life RTMS case A capital-intensive North American corporation, ABC,4 plans to build a new manufacturing plant in North America to produce electric car components based on an existing early-stage manufacturing system designed by a relatively workforce-intensive Asian company. ABC prefers building a new plant because: (1) ABC can better manage the potentially large demands and reduce transportation cost by producing these components in the same continent as its major customers rather than outsourcing them from another continent. (2) ABC may benefit from better product quality and more reliable supply chain by producing these components with better automated-machinery. The existing manufacturing system is a flow-type production line which needs to be reconfigured to incorporate automated machinery options to improve efficiency. ABC desires to reach a throughput level equal to an annual demand estimated through discussions with potential customers. Therefore, the objective is set as maximization of annual profit incorporating equipment, raw-material, inventory, and lost-sale costs. ABC established a team of engineers, production planners, and managers to decide which station structures are available for each station. A station structure includes specifying machinery vs. workforce selection and machine type. For each station, ABC engineering team chose one station structure out of all possible combinations by elimination based on performance and cost. The selected station structures and the material flow in the prospective production line are described in Section 4.1. Then, the team collected the necessary input such as costs and processing times. Based on these information and overall station structure setting, we develop the Operations Research approach summarized in Fig. 1 to find the approximately optimal station configuration variables including the number of machines, machine speeds, min/max buffer allocations, and number of workers for particular stations. The approach starts with the creation of a set of decision variable values based on the selected simulation-based optimization algorithm. These values are fed into the simulation module that is run for three replications. The average result is fed back into the optimization module to be evaluated by the algorithm. If the 4

Actual name is withheld for confidentiality.

297

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Simulation Module

Optimization Module

New Configuration of Decision Variables

Create new set of decision variables NO Stopping Criteria met?

Simulation Model of the Manufacturing System Simulation runs/replications

Profit/Cost Calculations

INITIALIZE

Simulation Results

YES

Optimization Output

Optimization Heuristics (Ant Colony, SA, TS etc...)

Objective Function Value

Fig. 1. Simulation optimization process.

stopping criteria are not met, the algorithm determines another set of decision variable values. 3. Background and literature review Because most technology transfers are from capital-intensive to workforce-intensive markets, there are few studies on RTMS (Elshout, 1995). However, more studies are likely to appear on this concept as several capital-intensive manufacturers in countries such as US, Canada, and Japan consider re-opening some of their domestic manufacturing operations (Economist, 2013; Nikkei, 2013; Carmichael, 2012). In addition, the recent trend of reshoring, i.e., bringing the manufacturing part of the supply chain back to or closer to the western countries for profitability, receives increasing interest (Ellram, 2013). There is a richer literature on the optimal production line configuration reviewed by Bergeron, Jamali, and Yamamoto (2010) and Papadopoulos and Heavey (1996). The studies in this literature generally use either mathematical/analytical models or simulation modeling. While most studies focus on optimizing buffer size allocation under different system settings as summarized by Tempelmeier (2003), several studies consider other variables such as the number of machines, processing times, etc. The studies proposing mathematical modeling approaches generally use stochastic models, LP, NLP, MIP, or NMIP to optimize the production line configuration under simplifying assumptions or/ and by decomposing the production line into manageable subsections (Li, 2005; Liu & Li, 2010; Nahas, Nourelfath, & Ait-Kadi, 2009; Sadr & Malhame, 2004). For tractability, these studies generally limit themselves to production lines with few stations and up-to 3 types of configuration variables. Because this research is based on a real-life case, we need to consider a problem of 24 stations and 5 types of configuration variables, for which developing an analytical model is not practical; however, we provide a conceptual model capturing the objective function in Section 4.3. There are several simulation-based optimization techniques proposed in the literature to a wide range of problems including the optimal production line configuration problems (Long-Fei & Le-Yuan, 2013; Tekin & Sabuncuoglu, 2004). Some studies on the optimal buffer size allocation in production lines use simulationbased optimization based on meta-heuristics including genetic algorithm (GA) (Bulgak, 2006; Yamamoto, Qudeiri, & Marui, 2008), simulated annealing (SA) (Spieckermann, Gutenschwager, Heinzel, & Voss, 2000), and tabu search (TS) (Alabas, Altiparmak, & Dengiz, 2002; Lutz, Davis, & Sun, 1998), or meta-models

including artificial neural networks (ANN) (Altiparmak, Dengiz, & Bulgak, 2007; Bulgak, 2006), response surface (Lin, Madu, & Kuei, 1994), and OptQuest (Melouk, Freeman, Miller, & Dunning, 2013). Similarly, there are studies addressing diverse set of production and logistic problems that propose simulation-based optimization approaches using GA, TS, ANN, local search, scatter search, and response surface methodologies (Bassi, Filho, Martins, & Bahiense, 2012; Baykasoglu, 2008; Cordeau, Legato, Mazza, & Trunfio, 2015; Dengiz, Bektas, & Ultanir, 2006; Keskin, Melouk, & Meyer, 2010; Legato, Mazza, & Gullì, 2014; Melouk et al., 2013; Yang, Kuo, & Chang, 2004). Table 1, which summarizes the recent studies in the production line configuration literature, illustrates that no study considers a production line setting as complex as ours and no more than a few search mechanisms are tested in most simulation-based optimization studies (Kleijnen, 2008). Thus, this study is different than those in the literature because: (1) We consider a more complex system and base our analysis on data from a real-life case. (2) Unlike other studies, we propose several different simulationbased optimization approaches and report their relative performances in the solution quality and speed spectrum. Although ABC engineering team reduced the number of station structure combinations to one in our case, a significant number of different station structure combinations need to be considered in a general RTMS case where a simulation-based optimization needs to be conducted for each station structure combination. Therefore, it is desirable to identify the simulation-based optimization methods that are both accurate and fast. 4. Problem description and simulation model 4.1. Process flow in the transferred production line Fig. 2 illustrates the process flow in the prospective production line. The series-parallel flow is asynchronous, inhomogeneous, saturated, and produces a single product. The flow also requires splitting subparts at Station1 and merging subparts from Line 1 and 2 at Station13 after a cool-down delay at the preceding buffer zones (8, 9, 10, and 16). Raw materials enter the system at eight stations (1, 2, 7, 8, 13, 14, 15, and 24). While six stations (2, 4, 5, 8, 10, and 11) produce scrap, three stations (14, 15, and 24) may reject parts during quality control. Parts visit Station17 three times (i.e., 16 ? 17 ? 18 ? 19 ? 20 ? 21 ? 17 ? 22 ? 17 ? 23) to be reprocessed before leaving for Buffer Zone 26. All machines are associated with

298

Table 1 Classification of production line problems. (See below-mentioned references for further information.)

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

299

300

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Fig. 2. Process flow of ABC production line.

Table 2 Decision variables and range of feasible values. Variable

Decision variables

Min

Max

Number of values

Increment

V1 V2 V3 V4 V5 V6 V7 V8 V9

Machine speed at P2 (parts/minute) Machine speed at P8 (parts/minute) Number of machines at P5&P11 Number of machines at P13 Min. buffer at B17 (unit: days of production by P14)a Min. buffer at B18 (unit: days of production by P15) Number of machines at P19 Number of machines at P21 Number of workers/shift at P17-P24

14.5 14.5 5 17 0.01 0.01 23 27 3

16 16 6 20 1.51 0.05 26 30 4

4 4 2 4 4 3 4 4 2

0.5 0.5 1 1 0.5 0.02 1 1 1

Number of combinations a

49,152

For example, daily production of P15 is 4128 parts, so 1.51 equals 6234 parts.

stochastic downtimes and repair times. While nine stations (2, 3, 4, 5, 8, 9, 10, 11, and 13) require stochastic setup times, six stations (14, 15, 16, 17, 20, and 24) require labor input with stochastic processing times. System works 6-days/week, 3-shifts/day, and 8-h/ shift with the exception of Stations 17–24 that run 7-days/week. Labor requirements for all stations are predetermined with the exception of Stations 17–24, which shares a single type of worker.

4.2. Decision variables for configuration Fig. 2 shows the stations and buffer zones associated with decision variables. Nine decision variables leading to a total of 49,152 possible production line configurations are described in Table 2. In the table, Pi and Bj refer to Station i and Buffer Zone j. The minimum, maximum, and increment values of these decision variables

301

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312 Table 3 Explanation of notations. Variable

Explanation

i j i ðjÞ Ri Vj p r inv

Workstation index, i ¼ 1; 2; . . . ; n Buffer zone index, j ¼ 1; 2; . . . ; m Workstation that feeds Buffer Zone j Set of labor types at Workstation i Set of buffer zones preceding Buffer Zone j Sales price/product Rate of inventory carrying cost Rate of lost-sale cost

r lst r int

p T d

e t

s D Dt S t ð UÞ Pr t ðUÞ I j s ð UÞ Imt ðUÞ Ij ðUÞ Im ðUÞ

Annual interest rate Lost-sale cost/product Time horizon for the problem (T ¼ 1 year for our case) Time interval between each end-product shipment Smaller time interval used by the simulation to estimate the average inventory in buffer zones Simulation time in terms of d; t ¼ 0; 1; 2; . . . ; T=d Simulation time in terms of e, s ¼ 0; 1; 2; . . . ; T=e Annual demand Demand between ðt  1Þ and t Sales between ðt  1Þ and t Production between ðt  1Þ and t Work-in-process inventory at Buffer Zone j at time s Finished goods inventory at time t Annual average work-in-process inventory at Buffer Zone j

PT i C vj a

Annual average finished goods inventory Investment and annual amortized cost per machine at Workstation i; respectively Number of machines and type l workers employed at Workstation i; respectively Annual cost of a type l worker Inventory carrying cost for having one-unit inventory for the whole year at Buffer Zone j Total material costs/end-product Isolated cost of labor/part processed at the workstation that feeds Buffer Zone j Isolated cost of equipment/part processed at the workstation that feeds Buffer Zone j Isolated raw material cost/part processed at the workstation that feeds Buffer Zone j Processing time at Workstation i (unit: years) Cumulative value-added cost/part waiting at the Buffer Zone j

Lt ðUÞ

Lost sales between ðt  1Þ and t

A CM i ; Ci

M i ; W li Cw l hj C mtr C lb j C eq j C mt j

reasonable maximum buffer capacities because the profit maximizing objective function penalizes high work-in-process (WIP) inventory. 4.3. Objective function The objective function P that inexplicitly defines the profit associated with a particular set of station configuration variables is given below. The objective function is a function of all decision variables denoted by U ¼ ½V 1 ; V 2 ; . . . ; V 9  2 Q and derived using the simulation outputs. Table 3 provides the notation used in the objective function. eq va Note that, Mi , W li , hj , C lb j , C j , and C j are also functions of U for

some i and j values because number of workers and machines are decision variables for some workstations and buffer zones. The objective function, which maximizes the annual profit, is given below where the first term represents expected sales income, while the rest are expected lost-sale, material, equipment, labor, inventory holding, and other (independent of U, e.g., overhead, sales/marketing) costs, respectively: T=d T=d T=d n X X X X maxPðUÞ ¼ maxp St ðUÞ  pLt ðUÞ  C mtr Prt ðUÞ  C iA Mi

U

U

t¼1

t¼1

t¼1

i¼1

n X m X X  Cw hjIj ðUÞ  C other l W li  i¼1 l2Ri

j¼1

Eqs. (1)–(12) provide further details about the components of the objective function.

Dt ¼ D  d=T;

t ¼ 0; 1; 2; . . . ; T=d

ð1Þ

Imt ðUÞ ¼ max½ðImðt1Þ ðUÞ þ Prt ðUÞ  Dt Þ; 0; Im0 ðUÞ ¼ 0; t ¼ 0; 1; 2; . . . ; T=d

ð2Þ

St ðUÞ ¼ min½ðImðt1Þ ðUÞ þ Prt ðUÞÞ; Dt ; " # 1 1  ð1 þ r int Þ A M Ci ¼ Ci = ; r int

t ¼ 0; 1; 2; . . . ; T=d

ð3Þ

i ¼ 1; 2; . . . ; n ðderiv ed from financial annuity formulationÞ are determined based on the machine specifications provided by the producers, and through consultations with ABC engineering team. Note that, by solving instances of this production line configuration problem using wider decision variable ranges, we confirm that wider ranges than those estimated by the engineering team do not necessarily produce better results. V5 and V6 denote the minimum buffers for B17 and B18. The minimum buffer amounts need to be considered for unbalanced lines where there exists a preceding station operating slightly slower on average than the latter one. For example, P13 consists of at least 17 identical machines with deterministic processing times and variable setup times. Assuming that all workstations are set up to work at full capacity, the variability in the setup and processing times of P13 and P14 may create an imbalance in favor of P14 for some configurations. A minimum buffer amount is defined for B17 to diminish the effect of this imbalance on the overall performance. If ever the amount of parts in B17 becomes zero, P14 stops and waits until the amount reaches the minimum level. In this way, P14 works continuously without stopping repeatedly in short time periods. The bulk time during which P14 waits for B17 can be utilized for maintenance. ABC also needs to specify the space allocated for each buffer zone. For the sake of simplicity, we set the maximum buffer capacities to the maximum sub-product in the corresponding buffer zones during the simulation runs following the method proposed by Battini, Persona, and Regattieri (2009). This approach leads to

C lb j ¼

X

 PT i ðjÞ C w l W li ðjÞ ;

ð4Þ

j ¼ 1; 2; . . . ; m

ð5Þ

l2Ri ðjÞ A   C eq j ¼ PT i ðjÞ C i ðjÞ ;

C vv a ¼

X

C lb j þ

j2Vv

X

j ¼ 1; 2; . . . ; m C eq j þ

j2Vv

hj ¼ r inv C jv a ;

X

C mt j ;

ð6Þ

v ¼ 1; 2; . . . ; m

ð7Þ

j2Vv

j ¼ 1; 2; . . . ; m

ð8Þ

p ¼ rlst p

ð9Þ

Lt ðUÞ ¼ max½ðDt  ðImðt1Þ ðUÞ þ Prt ðUÞÞÞ; 0;

t ¼ 0; 1; 2; . . . ; T=d ð10Þ

Im ðUÞ ¼ ð

ððT=dÞ1Þ X

ðImt ðUÞ þ

t¼0

Ij ðUÞ ¼

T=e X Ijs ðUÞ=ðT=eÞ;

Prðtþ1Þ ðUÞ ÞÞ=ðT=dÞ 2

ð11Þ

j ¼ 1; 2; . . . ; m  1

ð12Þ

s¼1

Eq. (7) calculates the production cost of a sub-product at any Buffer Zone v by adding up the labor, equipment, and material cost accrued so far. Eq. (8) computes the cost of holding one unit of sub-product at Buffer Zone j for one year. Eq. (2) derives the amount

302

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

of finished goods inventory at every t used in Eq. (11) to calculate the annual average inventory. Unlike Eq. (11), Eq. (12) does not take the next term’s production into account, because interval e is small enough to ignore its effect. Fig. 3a illustrates possible movements in the finished goods inventory, Imt ðUÞ, while Fig. 3b shows those in WIP inventory of Buffer Zone j, Ijs ðUÞ. The simulation model updates inventory levels discretely at each time period e. In Fig. 3a, sharp declines at every d represent product shipments. 4.4. Simulation model of the real-life case We develop the simulation model using SIMUL8 software package. The expected annual demand, detailed cost values, interest rate, and product and raw material prices were provided by ABC engineering team. Processing times, setup times, and efficiencies of the machinery are derived from the design specifications provided by the suppliers. The production rates are assumed to be deterministic only for the automated machines. We assume that processing times of non-/semi-automated machinery, and

workforce-based setup times of some automated machinery have triangular distributions; a common assumption in the literature (Law, 2007). Furthermore, 85% machine efficiency is ensured in the model using exponential and Erlang inter-arrival times for breakdowns and repairs. We set the warm-up period to two weeks as the production rate becomes stable after about 12 days. Each simulation replication takes 56.4 s (on a PC with IntelÒ CoreTM i7-3770 3.4 GHz CPU) on average. The simulation-based optimization approaches take three replications of the simulation model to evaluate each configuration. Tests conducted on 100 random configurations fail to show significant difference between mean profits from 3 and 30 replications. 4.5. Test-bed problem In order to determine specifications of the proposed simulationbased optimization methods, a smaller-size test-bed problem (representing real-life case on a smaller scale) was generated, because running the simulation model of the real-life case for testing pur-

Fig. 3. (a) Finished product inventory ðImt ðUÞÞ movements. (b) WIP inventory ðIjs ðUÞÞ movements.

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

poses was time consuming. Although comparison of proposed simulation-based optimization methods on the real-life case constitutes our main numerical results, we also evaluate the performances of these proposed methods compared to the best solution obtained by total enumeration of the decision variables in the test-bed problem (which is not possible to do on the reallife case). We refer to this best solution for the test-bed problem as the optimal solution for evaluating all meaningful configuration variable combinations. The test-bed problem reflects the features of the real-life case on a smaller scale, e.g. having asynchronous and inhomogeneous flow with assembly stations, re-entrant flow, stochastic processing and failure times, etc. The test-bed problem has 8 decision variables with 46,080 possible configurations. Each simulation replication takes 8.5 s. The process flow and input data for the test-bed problem are given in Appendix A. Corresponding SIMUL8 code can be obtained from the authors for comparisons. 5. Proposed simulation-based optimization approaches 5.1. OptQuest OptQuest is a commercial simulation-based optimization tool available for simulation packages such as SIMUL8, Arena, and FlexSim. OptQuest is a meta-heuristic that moves from solution to solution following a combination of TS, ANN, and Scatter Search logics till the desired number of iterations is reached (Kleijnen & Wan, 2007). We let OptQuest to pitch its default initial solution (mid-point for each decision variables) and set the number of iterations as 500, because OptQuest’s performance became stable after 400 iterations in the initial experiments on the test-bed problem. 5.2. Simulated annealing (SA) We develop an SA approach for benchmarking because it is a commonly used method for simulation-based optimization (Alrefaei & Diabat, 2009; Haddock & Mittenthal, 1992; Prudius & Andradottir, 2012; Rosen & Harmonosky, 2005). The proposed SA algorithm is as follows: BEGIN BestSolution = InitialSolution Temp = InitialTemperature UNTIL (No better solution can be found) DO Search immediate neighborhood of BestSolution for NewSolution(s)   DP Apply acceptance test Pi ¼ exp Temp IF NewSolution is accepted THEN BestSolution = NewSolution IF Temperature decrease criteria met THEN Temp ¼ ð1  dÞ  Temp END UNTIL END

where Pi is the probability of accepting solution i. DP is the difference between the profit values of the NewSolution and the current BestSolution. Temp stands for the current temperature, and d is the rate of reduction in the temperature. We take the best among 10 random solutions as the initial solution since tests with more random solutions failed to produce significantly better solution quality. We set the initial temperature as 20,000 with d = 0.10 and decrease it every 10 iterations. Immediate neighborhood of a solution is formed using two adjacent values of each variable (e.g., see Table 4).

303

The algorithm stops when a neighborhood fails to produce an accepted solution. 5.3. Ant-colony optimization (ACO) Ant-colony optimization (ACO) has not been widely applied to simulation-based optimization (Mohan & Baskaran, 2012). ACO logic considers combinatorial discrete optimization problems as a network where agents of an artificial ant-colony try to reach a target node (T) from a source node (S) through randomly-selected paths. Each path represents a solution, and agents leave pheromone trails on the paths they used. Because agents are more likely to choose a path that holds more pheromone than the others, better solutions are updated more often leading to fast convergence to an approximately optimal solution. We formulate our problem as a network based on possible values of each decision variable as shown in Fig. 4. Each column represents a decision variable and the shaded path-of-nodes represents a feasible solution. For each path/solution, the simulation model is run and the corresponding expected total profit value is reported as a measure of solution quality. The following pseudoalgorithm explains the steps of the proposed ACO algorithm. BEGIN Define initial CurrentBest = e  D  p UNTIL (max # of iterations reached) DO Construct the solution for Agent k using the path between Nodes i and j with probability 8 sa < P ij a ; if j 2 N ki k s l2Nk ij P ij ¼ i : 0; if j R N ki Update pheromone trail Evaporate pheromone from all arcs sij ¼ ð1  qÞsij 8ði; jÞ 2 A IF Solution found by Agent k is better than CurrentBest then; Update the pheromone levels on the arcs of the solution path

sij ¼ ð1 þ cÞsij 8ði; jÞ 2 Sk CurrentBest = Solution found by Agent k END IF NEXT k END UNTIL END

In this algorithm, CurrentBest represents the best profit value found so far whose initial value is a fraction (e) of the maximum annual sales income (p  D), and a is a constant which leads to faster convergence as its value gets higher. A is the set of nodes and N ki is the neighborhood of Agent k when in Node i. In each iteration, Agent k follows the arc between Nodes i and j with probability Pkij . Sk is the set of nodes in the solution path found by agent k. sij is the amount of pheromone on the arc between Nodes i and j, which evaporates with rate q after each iteration while pheromone level on new-found good solutions is increased with rate c. In order to derive a, e, p, and c parameters for an efficient implementation of the ACO algorithm, we first obtain an initial set of values for these parameters from the existing applications in the literature (Stützle et al., 2011). We set the search range for each parameter by varying them one-by-one and run the simulation model till the performance of ACO algorithm stops improving. Once the ranges are specified, we divide each range into almostequal-size grids and specify the best parameter combination by

304

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Table 4 Sample neighborhood. V01

V02

V03

V04

V05

V06

V07

V08

V09

Neighborhood of

16 15.5 15 16 16 16 16 16 16 16 16 16 16 16 16 16 16

15 15 15 15.5 14.5 15 15 15 15 15 15 15 15 15 15 15 15

6 6 6 6 6 5 6 6 6 6 6 6 6 6 6 6 6

18 18 18 18 18 18 19 17 18 18 18 18 18 18 18 18 18

0.51 0.51 0.51 0.51 0.51 0.51 0.51 0.51 1.01 0.01 0.51 0.51 0.51 0.51 0.51 0.51 0.51

0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.05 0.01 0.03 0.03 0.03 0.03 0.03

25 25 25 25 25 25 25 25 25 25 25 25 25 26 24 25 25

28 28 28 28 28 28 28 28 28 28 28 28 28 28 28 29 27

4 4 4 4 4 4 4 4 4 4 4 4 3 4 4 4 4

Base Solution V01 V01 V02 V02 V03 V04 V04 V05 V05 V06 V06 V07 V08 V08 V09 V09

V01

V02

14.5

14.5

V03

V04

V05

17

0.01

V06

V07

V08

V09

23

27

3

24

28

4

25

29

26

30

0.01 15.0

15.0

5

18

0.51

15.5

15.5

6

19

1.01

S

0.03

T

0.05 16.0

16.0

20

1.51

Fig. 4. Possible values of decision variables structured as a network.

optimal solution. Table 5 shows that, a = 3, e = 0.025, p = 0.1, and c = 0.35 is the best parameter combination based on these initial experiments. The ACO algorithm starts with a random solution where the initial pheromone amount (sij ) is the same (i.e., 1) for each trail. We define the maximum number of iterations as 300

enumerating all combinations on the corresponding grid system. Table 5 illustrates the percentage performances (compared to optimal solution) for a sample of the evaluated a, e, p, and c combinations. Percentage performance refers to the proportion of the total profit of the ACO solution for a given combination to that of the Table 5 ACO parameter selection. Combination no.

a

q

c

e

Converged at iteration

Percentage performance (%)

22 23 24 28 29 30 34 35 36 39 40 41 42 45 46 47 48 51 52 53 54

3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4

0.05 0.05 0.05 0.10 0.10 0.10 0.15 0.15 0.15 0.05 0.05 0.05 0.05 0.10 0.10 0.10 0.10 0.15 0.15 0.15 0.15

0.35 0.40 0.40 0.35 0.40 0.40 0.35 0.40 0.40 0.35 0.35 0.40 0.40 0.35 0.35 0.40 0.40 0.35 0.35 0.40 0.40

0.025 0.01 0.025 0.025 0.01 0.025 0.025 0.01 0.025 0.01 0.025 0.01 0.025 0.01 0.025 0.01 0.025 0.01 0.025 0.01 0.025

50 215 45 50 206 45 51 265 45 102 7 39 7 102 7 39 7 102 7 39 7

86 80 83 86 80 83 83 80 83 80 21 66 21 80 21 66 21 80 21 66 21

a = (2, 3, 4); q = (0.05, 0.10, 0.15); c = (0.25, 0.35, 0.40); e = (0.01, 0.025) Only the results of the combinations for which ACO convergences to a solution within 300 iterations are displayed. Tie between combinations 22 and 28 is broken by additional tests on the real-life problem.

305

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

based on these initial experiments with the test-bed problem (see Appendix B). 5.4. Response surface methodology (RSM) RSM is a method to analyze the association between a number of input variables and outcomes. Although RSM is intended for solving problems with continuous variables, it is also applicable for problems with discrete variables. Simulation-based optimization with RSM is described in detail by Kleijnen (2008). It starts with the same initial solution as employed in SA approach, and calculates the solutions in the immediate neighborhood (as shown in Table 4). A linear regression model is developed based on the neighborhood solution space. The value of the decision variable with the highest regression coefficient (absolute value) is increased (if coefficient is positive) or decreased (if negative) one step to generate the next solution. The algorithm ends when there is no room left for improvement. BEGIN BestSolution = InitialSolution UNTIL (No better solution can be found) DO Calculate all Solutions in the immediate neighborhood of BestSolution. Develop linear regression model using the Solutions. (dependent variable: profit, independent variables: decision variables). Decrease/increase the value of the decision variable with the highest coefficient, Generate NextSolution IF NextSolution is better than BestSolution THEN BestSolution = NextSolution END UNTIL END

Table 4) and continues the search with it. FGS selects the first better solution while scanning the immediate neighborhood of the current solution instead of enumerating the whole neighborhood. Search ends when a neighborhood fails to produce an improved solution. To determine whether combinations of any these methods can produce better results, we propose two hybrid methods: (1) ACO_GS takes the best combination from ACO as the initial solution for GS; (2) ACO_OptQ uses the result from ACO as the initial solution for OptQuest. Among these two hybrid algorithms, ACO_GS performs very well both in the test-bed and real-life problem. In most cases, it escapes the local optima and reaches to global optimum. A detailed convergence analysis of ACO_GS is provided in Appendix B. 6. Results 6.1. Numerical experiment setting Tables 6 and 7 display the performances of the solutions found by each simulation-based optimization algorithm for the test-bed and real-life problems (respectively) under different demand and inventory cost scenarios. Different scenarios are analyzed because a fluctuation between +10% and 15% in the constant annual demand (D) and a ±25% change in the rate of inventory carrying cost (rinv Þ are anticipated. In total, we consider six different demand and inventory cost scenarios. While the base-case represents an annual demand of 1xD and inventory carrying cost rate of rinv ¼ 0:2, the sensitivity scenarios consider the following demand and inventory carrying cost cases: 0.85  D, rinv ¼ 0:2; 0.9  D, r inv ¼ 0:2; 1.1  D, r inv ¼ 0:2; 1  D, rinv ¼ 0:15; 1  D, rinv ¼ 0:25. 6.2. Performance metrics

5.5. Greedy search (GS), fast Greedy search (FGS), and hybrid methods GS uses the same initial solution as RSM. Then it simply selects the solution with the best profit in the neighborhood (as shown in

We consider two metrics for evaluating the performances of the proposed simulation-based optimization methods: solution quality and required computation time. Solution quality is captured by the expected total profit associated with the best solution found

Table 6 Results for the test-bed problem. Methods Problem type

Max. expected profit (Optimal)

OptQuest

SA

ACO

RSM

GS

FGS

ACO_GS

ACO_OptQ

Best expected profit Max @ Iterations

100.00% 343 500

100.00% 73 89

86.18% 60 72

97.62% 66 87

100.00% 89 111

100.00% 51 65

100.00% 97 113

100.00% 310 500

Sensitivity analysis 0.85  D r inv = 0.2 Best expected profit Max @ Iterations

100.00% 189 500

84.74% 38 51

89.34% 212 215

63.19% 134 138

92.28% 83 96

63.19% 50 88

89.34% 215 231

99.48% 323 500

Base case

$2,966,279

(Optimal) $880,120

0.9  D rinv = 0.2

Best expected profit Max @ Iterations

92.23% 364 500

100.00% 103 107

84.22% 37 39

92.47% 36 68

100.00% 65 82

100.00% 73 88

92.98% 78 92

100.00% 281 500

$2,036,146

1.1  D rinv = 0.2

Best expected profit Max @ Iterations

100.00% 377 500

100.00% 77 85

97.71% 71 75

97.19% 66 81

100.00% 90 108

100.00% 51 60

100.00% 79 101

100.00% 194 500

$2,516,279

1  D r inv = 0.15

Best expected profit Max @ Iterations

100.00% 465 500

100.00% 73 87

98.82% 126 134

97.72% 66 94

100.00% 89 107

100.00% 51 65

100.00% 147 165

100.00% 218 500

$3,107,713

1  D r inv = 0.25

Best expected profit Max @ Iterations

100.00% 465 500

100.00% 73 87

98.69% 87 159

97.50% 66 94

100.00% 89 107

100.00% 51 65

100.00% 166 182

100.00% 241 500

$2,824,829

Best expected profit refers to Pers, i.e., the percentage of the profit associated with the best solution found by each method with respect to the maximum expected profit.

306

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Table 7 Results for the real-life case. Methods Problem type Base case

Best expected profit Max @ Iterations

Sensitivity analysis Best expected profit 0.85  D rinv = 0.2 Max @ Iterations

Max. expected profit (approximately optimal)

OptQuest

SA

ACO

RSM

GS

FGS

ACO_GS

ACO_OptQ

99.85% 253 500

99.95% 74 89

95.48% 48 60

99.51% 49 74

100.00% 84 106

100.00% 55 65

100.00% 130 149

99.85% 293 500

$3,839,504

100.00% 238 500

91.34% 42 61

97.98% 66 74

96.47% 72 106

96.15% 50 66

91.34% 39 54

100.00% 117 135

99.97% 471 500

$1,913,895

0.9  D rinv = 0.2

Best expected profit Max @ Iterations

99.92% 490 500

99.84% 117 131

97.09% 144 149

99.24% 182 190

99.94% 108 142

99.84% 126 140

100.00% 216 234

99.94% 378 500

$2,603,223

1.1  D rinv = 0.2

Best expected profit Max @ Iterations

100.00% 338 500

100.00% 97 114

99.14% 80 82

100.00% 256 260

100.00% 119 138

100.00% 60 75

100.00% 118 142

99.84% 278 500

$3,541,028

1  D r inv = 0.15

Best expected profit Max @ Iterations

99.83% 473 500

100.00% 106 121

95.55% 48 60

99.51% 49 74

100.00% 84 106

100.00% 61 76

100.00% 134 153

99.85% 326 500

$3,904,821

1  D r inv = 0.25

Best expected profit Max @ Iterations

99.97% 343 500

99.98% 74 89

95.41% 48 60

99.50% 49 74

100.00% 84 106

100.00% 57 76

100.00% 149 168

99.97% 327 500

$3,779,095

Best expected profit refers to Per s , i.e., the percentage of the profit associated with the best solution found by each method with respect to the maximum expected profit.

by each metaheuristic. The expected total profit is estimated as the average of the total profits from the corresponding simulation replications. Specifically, let s 2 ð1; 2; . . . ; 8Þ represent the considered metaheuristic and q 2 ð1; 2; . . . ; Ns ) denote a particular iteration of the sth metaheuristic. Then, the expected total profit of the solution at the qth iteration of the sth metaheuristic (Usq ) is P estimated as ðUsq Þ ¼ ð Uu¼1 P u ðUsq ÞÞ=U, where Pu ðUsq Þ is the total expected profit at the uth replication of the simulation models. Note that, P u ðUsq Þ is derived from simulation outputs based on the calculation defined in Section 4.3. In Tables 6 and 7, we present the percentage of the expected profit associated with the best solution derived by each proposed metaheuristic with respect to the best solution found amongst all methods for all parameter scenarios to illustrate the optimality gap. We denote the expected total profit and the best solution from ^ s and Us , i.e., P ^ s ¼ maxq PðUs Þ and the sth metaheuristic by P q

Us ¼ Usq whereq ¼ argmaxq PðUsq Þ. The maximum expected profit refers to the best solution found by total enumeration for the test-bed problem (i.e., maxU2Q PðUÞ where Q is the set of all possible solutions), and the solution of the simulation-based optimization method with the highest expected total profit for the real^1 ; P ^2 ; . . . ; P ^ 8 g). The percentages of the life problem (i.e., maxfP expected profit associated with the best solution of the sth metaheuristic with respect to the maximum expected profit are denoted by Pers and Pers for the test-bed and the real-life problems, respectively, where:

^ s =maxU2Q PðUÞ; Per s ¼ 100%  P

s 2 ð1; 2; . . . ; 8Þ

^ s =maxfP^ 1 ; P^ 2 ; . . . ; P ^ 8 g; Per s ¼ 100%  P

s 2 ð1; 2; . . . ; 8Þ

ð13Þ ð14Þ

In order to evaluate the required computation time and convergence rate of the proposed methods, we keep track of the number of necessary iterations (Iterations), and at which iteration the best solution is attained (Max @) for each metaheuristic, respectively. This is a valid approach because estimating the total expected profit via simulation replications is the most time consuming operation in the iterations of the metaheuristics. Therefore, the number

of necessary iterations, which specifies the total number of simulation replications to be taken, is an accurate indicator of required computation time for each metaheuristic.

6.3. Results for the test-bed and the real-life problems Table 6 illustrates the optimality gap associated with the proposed simulation-based optimization approaches by comparing them with total enumeration results. The comparison implies that ACO and RSM significantly deviate from the best solutions while FGS finds the best configuration parameter combinations in most cases and it is significantly faster than any other algorithms. Although it is myopic, FGS achieves the global optima in most cases while ACO and RSM often get stuck at local minima as shown in Fig. 5. However, in one case (0.85  D and 1  rinv ), FGS proposes a significantly inferior solution. As expected, GS is slower than FGS, whereas, associated with better solution quality. OptQuest, the slowest method, also either finds the best solutions in most cases or do not significantly deviate from the best solution. SA and ACO_GS are faster than OptQuest, while their solution qualities are comparable to that of OptQuest on the test-bed problem. Table 6 implies that OptQuest is a reliable solution technique for small size problems in terms of solution quality. Therefore, practitioners may prefer to use the easy-to-implement commercial tool for such problem instances. However, OptQuest is very slow; thus, it may not be an effective method for larger problems, or the general production line configuration problems in case of RTMS for which OptQuest needs to be run to find the approximatelyoptimal configuration parameters for a large number of possible production line reconfigurations. For small size problems, FGS provides a fast and almost reliable alternative. However, the practitioners may prefer to use GS, SA, or ACO_GS instead by sacrificing some computation time in order to avoid a potentially significant deviation from the optimal solution. The comparisons on the real-life problem in Table 7 show how the problem size may affect the performances of the proposed simulation-based optimization approaches. Note that, the profit values for the real-life case are divided by an undeclared real number for confidentiality purposes. In the real-life problem, both FGS

307

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Fig. 5. Search paths followed by ACO, RSM, and FGS on the solution space of the base-case test-bed problem.

and OptQuest deviate from the best solution in most cases; however, the maximum deviation from the best solution is much lower for OptQuest (100%  99.83% = 0.17%) compared to that for FGS (100%  91.34% = 8.66%). ACO_GS finds the best solutions for all problem instances. The solution quality of SA and GS are similar to each other and relatively lower than that of ACO_GS, while the maximum deviation from the best solution is much lower for SA compared to GS. According to these results, ACO_GS is likely to provide the best result even for large size problems. Therefore, decision makers facing similarly challenging optimal production line configuration problems may rely on an ACO_GS-type simulation-based optimization method to derive a good feasible solution much faster than OptQuest with similar (if not better) solution quality. Note that,

in our real-life case, the engineering team reduced the possible station structures to one for each station; therefore, we optimized the configuration variables once. However, in general RTMS cases, there could be a large number of possible reconfigurations (multiple structure types for each station) for the transferred production line, which may require optimizing configuration parameters multiple times. In such cases, using OptQuest (or even ACO_GS) may not be feasible or efficient. In such cases, practitioners may prefer to use SA or GS to trade some of the solution quality with a significant relief from the computation time. Table 7 also illustrates that the solution quality for ACO_GS and OptQuest does not change significantly with the varying levels of demand and inventory holding cost. Similarly, GS and FGS are robust to the changes in inventory holding costs. However, their

Table 8 Values of decision variables for best expected profit. Method

V01

V02

V03

V04

V05

V06

V07

V08

V09

Performance (%)

OptQuest SA ACO RSM GS FGS ACO_GS ACO_OptQ

15.5 16.0 16.0 16.0 16.0 16.0 16.0 15.5

14.5 14.5 15.0 14.5 14.5 14.5 14.5 14.5

6 6 6 6 6 6 6 6

19 18 17 18 18 18 18 19

1.01 1.01 1.01 1.51 1.51 1.51 1.51 1.01

0.03 0.01 0.03 0.03 0.05 0.05 0.05 0.03

25 24 24 26 24 24 24 25

28 27 28 27 27 27 27 28

4 4 4 4 4 4 4 4

99.85 99.95 95.48 99.51 100.00 100.00 100.00 99.85

308

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

solution quality is better when the demand is high (Demand = 1.1  D). In case of low demand, there may be limited number of effective system configurations with high solution quality; thus, a myopic search algorithm may get stuck in a local optimum more easily. Table 8 shows values of the decision variables for the best solution of each simulation-based optimization approach. Our numerical analysis shows that V03 (Number of machines at P5 and P11) and V09 (Number of workers/shift at P17 – P24) are the most important decision variables as a change in their value significantly affects the profit.

7. Conclusion The optimal manufacturing system configuration in case of RTMS is likely to draw more attention in the future due to the recent trend of re-shoring practices. Decision makers facing such challenges should consider the variance and randomness in the complex nature of the transferred manufacturing system to avoid significant efficiency losses due to suboptimal implementations. In this context, practitioners may significantly benefit from developing simulation-based laboratory environments for testing whatif questions regarding the redesigned/transferred manufacturing system. For example, we observed significant performance differences between the approximately optimal production line configurations found by our proposed approaches, and the initial production line design recommended to ABC by the workforce oriented Asian company. The initial simulation results showed that the approximately optimal configurations are able to produce about 35% more end-products than the initial design, which even fails to satisfy the targeted production capacity,5 and causes a significant amount of lost sale cost. We observed that the proposed simulation-based optimization methods lie in a particular solution quality and speed spectrum. The practitioners facing similar problems may prefer to employ easy-to-implement software packages such as OptQuest for small-size problems and large-size ones with limited number of possible reconfigurations to be evaluated. However, the decision makers may need to invest in custom-developed faster approaches for computational feasibility when the number of possible station structure combinations is high or/and extensive sensitivity analyses are needed within a short time-window. In such cases, the proposed ACO_GS may reduce the computational burden and improve the solution quality. Decision makers may even rely on SA and GS to further reduce the computational burden; however, they should be aware of up to 15% and 8% deviations (may be more) in the solution quality compared to the best solution, respectively. Our analysis also shows that ACO logic can be utilized to develop efficient simulation-based optimization approaches for large-scale optimal MSC problems through enhancing its performance with a good initial solution. The simulation-based optimization methods analyzed in this study are all centralized algorithms considering the production line optimization problem from a holistic perspective. This is reasonable since determining the ideal performance of the production line is important at the design stage. Decentralized/distributed algorithms consider large problems/systems as a combination of smaller sub-problems/systems, and analyze them separately to avoid computational burden (Meng, Wang, & Liu, in press). Therefore, distributed algorithm methodology, which is applied to several problems in manufacturing and energy systems in the literature (Meng, Yang, & Sun, 2016; Weigert, Horn, & Werner, 2006), may be a promising approach for the considered problem. 5

Details of this calculation are omitted for confidentiality.

We analyze the potential of a distributed algorithm approach by dividing the production system for the test-bed problem into two, and optimizing each subsystem separately in an iterative fashion. Although it does not outperform ACO_GS and GS, this basic distributed algorithm approach shows a promising performance by solving several problem instances to optimality. However, we left a detailed implementation for future studies.

Acknowledgements We thank Simul8 Corporation for granting the professional edition research license. We also thank ABC Corporation for partially funding the research project.

Appendix A. Test-bed Problem See Table A1. Process flow is given in Fig. A1. Input data is given in Table A1, where framed shaded boxes represent decision variables. B8 is a special storage with limited capacity. Adding one unit of capacity has a cost and its capacity is a decision variable. Each part has to wait there for a specific period of time. Number of machines for P2 and P4 are identical, a single decision variable for both stations. List of decision variables are given in Table A2. Parts In/Out explains whether the processes break down parts into subparts or assemble subparts. For example, B2 routes in 6250 subparts each batch, which means P2 breaks down a part into 6250 subparts. B2 then routes out 11 subparts in each batch to P3 for assembling. There are similar processes in the real-life case, details of which are not explained in the paper due to our confidentiality agreement with ABC. P7, P8 and P9 shares same type of workforce, cost of which is given above. All other labor costs are predetermined and embedded into overhead cost. To simplify the model, scrap and rejections are not considered in test-bed problem. Also equipment amortization costs, raw material costs, and inventory holding costs are assumed to be pre-calculated constants. End products are shipped to buyers weekly.

Appendix B. Convergence of the ACO_GS algorithm We present a convergence analysis of ACO_GS based on its performance on the test-bed problem. In most cases, ACO_GS escapes the local optima and reaches to global optimum. For instance, Fig. A2 illustrates the solution pathway ACO_GS follows through its iterations (until reaching the global optimum) among all possible solutions for test-bed problem with the base case input parameters. In this figure, the performances of the possible solutions (as a percentage of the optimal profit) are denoted by blue lines and the solutions found by ACO_GS are illustrated using diamond-shaped red6 markers. Fig. A3 illustrates the average proportions of the performances of the current-best ACO_GS solutions to the global optimum under the base case and sensitivity scenarios for all iterations of the ACO_GS algorithm. This figure shows how the average percentage optimality gap reduces through iterations. Fig. A4 illustrates the same performance metric separately for all parameter scenarios. Fig. A4 indicates that the optimality gap reduces to less than 20% in the first 40 iterations and solutions stop improving by the 210th iteration. 6 For interpretation of color in Fig. A2, the reader is referred to the web version of this article.

309

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312 Table A1 Test-bed problem input data.

NM: Number of Machines, PD: Probability Distribution, Det.: Deterministic, Trng.: Triangular, Trigger: Setup time occurs after every trigger amount of parts were processed.

1 Raw Material

1

2

2 5

5 3

3

4

6

6

7

7

4

10

9 10

Final Product Fig. A1. Process flow of the test-bed problem.

9

8

8

310

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Table A2 Decision variables and range of feasible values. Variable ID

Decision variables

Min

Max

Number of values

Increment

V1 V2 V3 V4 V5 V6 V7 V8

Speed of Machine at P1 (meters/minute) Speed of Machine at P3 (meters/minute) Number of Machines at P2 & P4 Number of Machines at P5 Min. Buffer Amount at B5 (days of production by P6) Number of Workers at P7, P8 & P9 Capacity of B8 (cartridges, 64 parts each) Number of Machines at P9

23 23 4 25 0.01 2 28 35

26 26 7 28 1.01 4 32 38

4 4 4 4 3 3 5 4

1 1 1 1 0.5 1 1 1

Number of combinations

45,080

Fig. A2. Solution pathway for ACO_GS on the test-bed problem with base case parameters.

Fig. A3. Convergence graph for ACO_GS on the test-bed problem based on average proportion of the current-best solution to the global optimum.

Fig. A4. Convergence graph for ACO_GS on test-bed problem based on proportion of the current-best solution to the global optimum for each parameter scenario.

References Alabas, C., Altiparmak, F., & Dengiz, B. (2002). A comparison of the performance of artificial intelligence techniques for optimizing the number of kanbans. Journal of the Operational Research Society, 53, 907–914. Alrefaei, M. H., & Diabat, A. H. (2009). A simulated annealing technique for multiobjective simulation optimization. Applied Mathematics and Computation, 215, 3029–3035.

Altiparmak, F., Dengiz, B., & Bulgak, A. A. (2007). Buffer allocation and performance modeling in asynchronous assembly system operations: An artificial neural network metamodeling approach. Applied Soft Computing Journal, 7, 946–956. Aziz, A., Jarrahi, F., & Abdul-Kader, W. (2010). Modeling and performance evaluation of a series-parallel flow line system with finite buffers. Information Systems and Operational Research, 48, 103–120. Bassi, H. V., Filho, F., Martins, V. J., & Bahiense, L. (2012). Planning and scheduling a fleet of rigs using simulation–optimization. Computers & Industrial Engineering, 63, 1074–1088.

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312 Battini, D., Persona, A., & Regattieri, A. (2009). Buffer size design linked to reliability performance: A simulative study. Computers & Industrial Engineering, 56, 1633–1641. Baykasoglu, A. (2008). Gene expression programming based meta-modelling approach to production line design. International Journal of Computer Integrated Manufacturing, 21, 657–665. Bergeron, D., Jamali, M. A., & Yamamoto, H. (2010). Modelling and analysis of manufacturing systems: A review of existing models. International Journal of Product Development, 10, 46–61. Blumenfeld, D. E., & Li, J. (2005). An analytical formula for throughput of a production line with identical stations and random failures. Mathematical Problems in Engineering, 2005, 293–308. Bonvik, A. M., Dallery, Y., & Gershwin, S. B. (2000). Approximate analysis of production systems operated by a CONWIP/finite buffer hybrid control policy. International Journal of Production Research, 38, 2845–2869. Bulgak, A. A. (2006). Analysis and design of split and merge unpaced assembly systems by metamodelling and stochastic search. International Journal of Production Research, 44, 4067–4080. Carmichael, K. (2012). Canada’s ‘reshoring’ opportunity. Globe and Mail, June 4. Chan, F. T. S. (2001). Simulation analysis of maintenance policies in a flow line production system. International Journal of Computer Applications in Technology, 14, 78–86. Chan, F. T. S., & Ng, E. Y. H. (2002). Comparative evaluations of buffer allocation strategies in a serial production line. International Journal of Advanced Manufacturing Technology, 19, 789–800. Chehade, H., Amodeo, L., & Yalaoui, F. (2011). Ant colony optimization for multiobjective buffers sizing problems. In Ant colony optimization methods and applications (p. 303). InTech. Chiang, S., Hu, A., & Meerkov, S. M. (2008). Lean buffering in serial production lines with nonidentical exponential machines. IEEE Transactions on Automation Science and Engineering, 5, 298–306. Colledani, M., Gandola, F., Matta, A., & Tolio, T. (2008). Performance evaluation of linear and non-linear multi-product multi-stage lines with unreliable machines and finite homogeneous buffers. IIE Transactions, 40, 612–626. Cordeau, J., Legato, P., Mazza, R. M., & Trunfio, R. (2015). Simulation-based optimization for housekeeping in a container transshipment terminal. Computers & Operations Research, 53, 81–95. De Vericourt, F., & Gershwin, S. B. (2004). Performance evaluation of a make-tostock production line with a two-parameter-per-machine policy: The control point policy. IIE Transactions, 36, 221–236. Demir, L., Tunali, S., & Lkketangen, A. (2011). A tabu search approach for buffer allocation in production lines with unreliable machines. Engineering Optimization, 43, 213–231. Dengiz, B., Bektas, T., & Ultanir, A. E. (2006). Simulation optimization based DSS application: A diamond tool production line in industry. Simulation Modelling Practice and Theory, 14, 296–312. Economist (2013). Coming home; reshoring manufacturing. The Economist, Jan 19, 6–9. Ellram, L. M. (2013). Offshoring, reshoring and the manufacturing location decision. Journal of Supply Chain Management, 49, 3–5. Elshout, J. (1995). Significance of reverse transfer of technology and experiences with its promotion in industrialised countries. Technovation, 15, 423–435. Enginarlar, E., Li, J., & Meerkov, S. M. (2005). Lean buffering in serial production lines with non-exponential machines. OR Spectrum, 27, 195–219. Gershwin, S. B., & Burman, M. H. (2000). A decomposition method for analyzing inhomogeneous assembly/disassembly systems. Annals of Operations Research, 93, 91–115. Gershwin, S. B., & Schor, J. E. (2000). Efficient algorithms for buffer space allocation. Annals of Operations Research, 93, 117–144. Haddock, J., & Mittenthal, J. (1992). Simulation optimization using simulated annealing. Computers & Industrial Engineering, 22, 387–395. Helber, S. (2000). Approximate analysis of unreliable transfer lines with splits in the flow of material. Annals of Operations Research, 93, 217–243. Helber, S., Schimmelpfeng, K., Stolletz, R., & Lagershausen, S. (2011). Using linear programming to analyze and optimize stochastic flow lines. Annals of Operations Research, 182, 193–211. Hemachandra, N., & Eedupuganti, S. K. (2003). Performance analysis and buffer allocations in some open assembly systems. Computers & Operations Research, 30, 695–704. Hu, A. B., & Meerkov, S. M. (2006). Lean buffering in serial production lines with Bernoulli machines. Mathematical Problems in Engineering, 2006, 1–24. Huang, M., Chang, P., & Chou, Y. (2002). Buffer allocation in flow-shop-type production systems with general arrival and service patterns. Computers & Operations Research, 29, 103–121. Jeong, K. C., & Kim, Y. D. (2000). Heuristics for selecting machines and determining buffer capacities in assembly systems. Computers & Industrial Engineering, 38, 341–360. Keskin, B. B., Melouk, S. H., & Meyer, I. L. (2010). A simulation-optimization approach for integrated sourcing and inventory decisions. Computers & Operations Research, 37, 1648–1661. Kim, S., Park, Y., & Jun, C. (2006). Performance evaluation of re-entrant manufacturing system with production loss using mean value analysis. Computers & Operations Research, 33, 1308–1325. Kleijnen, J. P. C. (2008). In Jack P.C. Kleijnen. (Ed.), Design and analysis of simulation experiments. Boston: Springer.

311

Kleijnen, J. P. C., & Wan, J. (2007). Optimization of simulated systems: OptQuest and alternatives. Simulation Modelling Practice and Theory, 15, 354–362. Koren, Y., Hu, S. J., & Weber, T. W. (1998). Impact of manufacturing system configuration on performance. CIRP Annals - Manufacturing Technology, 47, 369–372. Kouikoglou, V. S. (2000). Optimal rate allocation in unreliable, assembly/ disassembly production networks with blocking. IEEE Transactions on Robotics and Automation, 16, 429–434. Kouikoglou, V. S. (2002). An efficient discrete event model of assembly/disassembly production networks. International Journal of Production Research, 40, 4485–4503. Law, A. M. (2007). Simulation modeling and analysis (4th ed.). Boston, Toronto: McGraw-Hill. Legato, P., Mazza, R. M., & Gullì, D. (2014). Integrating tactical and operational berth allocation decisions via Simulation-Optimization. Computers & Industrial Engineering, 78, 84–94. Li, J. (2004). Performance analysis of production systems with rework loops. IIE Transactions, 36, 755–765. Li, J. (2005). Overlapping decomposition: A system-theoretic method for modeling and analysis of complex manufacturing systems. IEEE Transactions on Automation Science and Engineering, 2, 40–53. Lin, C., Madu, C. N., & Kuei, C. (1994). Experimental design and regression analysis in simulation: An automated flowline case study. Microelectronics and Reliability, 34, 845–861. Liu, Y., & Li, J. (2010). Split and merge production systems: Performance analysis and structural properties. IIE Transactions, 42, 422–434. Long-Fei, W., & Le-Yuan, S. (2013). Simulation optimization: A review on theory and applications. Acta Automatica Sinica, 39, 1957–1968. Lutz, C. M., Davis, K. R., & Sun, M. (1998). Determining buffer location and size in production lines using tabu search. European Journal of Operational Research, 106, 301–316. Manitz, M. (2008). Queueing-model based analysis of assembly lines with finite buffers and general service times. Computers & Operations Research, 35, 2520–2536. Melouk, S. H., Freeman, N. K., Miller, D., & Dunning, M. (2013). Simulation optimization-based decision support tool for steel manufacturing. International Journal of Production Economics, 141, 269–276. Meng, W., Wang, X., & Liu, S. (2016). Distributed load sharing of an inverter-based microgrid with reduced communication. IEEE Transactions on Smart Grid. In press Meng, W., Yang, Q., & Sun, Y. (2016). Guaranteed performance control of DFIG variable-speed wind turbines. IEEE Transactions on Control Systems Technology, 24, 2215–2223. Mohan, B. C., & Baskaran, R. (2012). A survey: Ant colony optimization based recent research and implementation on several engineering domain. Expert Systems with Applications, 39, 4618–4627. Nahas, N., Nourelfath, M., & Ait-Kadi, D. (2008). Ant colonies for structure optimisation in a failure prone series-parallel production system. Journal of Quality in Maintenance Engineering, 14, 7–33. Nahas, N., Nourelfath, M., & Ait-Kadi, D. (2009). Selecting machines and buffers in unreliable series-parallel production lines. International Journal of Production Research, 47, 3741–3774. Nikkei (2013). Panasonic may bring white goods production back to Japan. The Nikkei, June 4. Nourelfath, M., Nahas, N., & Ati-Kadi, D. (2005). Optimal design of series production lines with unreliable machines and finite buffers. Journal of Quality in Maintenance Engineering, 11, 121–138. Papadopoulos, H. T., & Heavey, C. (1996). Queueing theory in manufacturing systems analysis and design: A classification of models for production and transfer lines. European Journal of Operational Research, 92, 1–27. Papadopoulos, C. T., O’Kelly, M. E. J., & Tsadiras, A. K. (2013). A DSS for the buffer allocation of production lines based on a comparative evaluation of a set of search algorithms. International Journal of Production Research, 51, 4175–4199. Prudius, A. A., & Andradottir, S. (2012). Averaging frameworks for simulation optimization with applications to simulated annealing. Naval Research Logistics, 59, 411–429. Qudeiri, J. A., Yamamoto, H., & Ramli, R. (2009). Buffer size decision for flexible transfer line with rework paths using genetic algorithm. International Journal of Intelligent Systems Technologies and Applications, 7, 227–240. Qudeiri, J. A., Yamamoto, H., Ramli, R., & Jamali, A. (2008). Genetic algorithm for buffer size and work station capacity in serial-parallel production lines. Artificial Life and Robotics, 12, 102–106. Rosen, S. L., & Harmonosky, C. M. (2005). An improved simulated annealing simulation optimization method for discrete parameter stochastic systems. Computers & Operations Research, 32, 343–358. Rossi, A., & Lanzetta, M. (2013). Scheduling flow lines with buffers by ant colony digraph. Expert Systems with Applications, 40, 3328–3340. Sabuncuoglu, I., Erel, E., & Kok, A. G. (2002). Analysis of assembly systems for interdeparture time variability and throughput. IIE Transactions, 34, 23–40. Sadr, J., & Malhame, R. P. (2004). Unreliable transfer lines: Decomposition/ aggregation and optimization. Annals of Operations Research, 125, 167–190. Shaaban, S., & McNamara, T. (2009). The effects of joint operations time means and variability unbalance on production line performance. International Journal of Manufacturing Technology and Management, 18, 59–78.

312

M.F. Yegul et al. / Computers & Industrial Engineering 109 (2017) 295–312

Shi, C., & Gershwin, S. B. (2009). An efficient buffer design algorithm for production line profit maximization. International Journal of Production Economics, 122, 725–740. Spieckermann, S., Gutenschwager, K., Heinzel, H., & Voss, S. (2000). Simulationbased optimization in the automotive industry - A case study on body shop design. Simulation, 75, 276–286. Stützle, T., López-Ibánez, M., Pellegrini, P., Maur, M., De Oca, M. M., Birattari, M., Dorigo, M. (2011). Parameter adaptation in ant colony optimization. In Autonomous search (pp. 191–215). Springer. Tekin, E., & Sabuncuoglu, I. (2004). Simulation optimization: A comprehensive review on theory and applications. IIE Transactions, 36, 1067–1081. Tempelmeier, H. (2003). Practical considerations in the optimization of flow production systems. International Journal of Production Research, 41, 149–170. Vergara, H. A., & Kim, D. S. (2009). A new method for the placement of buffers in serial production lines. International Journal of Production Research, 47, 4437–4456. Wang, Q., & Chatwin, C. R. (2005). Key issues and developments in modelling and simulation-based methodologies for manufacturing systems analysis, design

and performance evaluation. The International Journal of Advanced Manufacturing Technology, 25, 1254–1265. Weigert, G., Horn, S., & Werner, S. (2006). Optimization of manufacturing processes by distributed simulation. International Journal of Production Research, 44, 3677–3692. Xia, B., Xi, L., & Zhou, B. (2010). A method for performance evaluation of transfer line with unreliable machines and finite transfer-delay buffers. Advanced Materials Research, 97–101, 3135–3138. Xu, M., Abdul-Kader, W., & Ganjavi, O. (2009). Cycle time estimation of a multiple product production line: An approach with fuzzy and chance-constrained programming. Information Systems and Operational Research, 47, 93–103. Yamamoto, H., Qudeiri, J. A., & Marui, E. (2008). Definition of FTL with bypass lines and its simulator for buffer size decision. International Journal of Production Economics, 112, 18–25. Yang, T., Kuo, Y., & Chang, I. (2004). Tabu-search simulation optimization approach for flow-shop scheduling with multiple processors - A case study. International Journal of Production Research, 42, 4015–4030.