5. A hybrid simulation-based heuristics solution method The use of simulation modeling can be both cost efficient and time effective for evaluating system behaviors using multiple performance measures. Simulation can be particularly more useful when system complexities intensify by shifting from deterministic to stochastic modeling. In most cases, a set of scenarios are developed representing possible settings of critical system parameters. The scenarios are then used to analyze and compare the system performance under each setting. For example, using simulation modeling, Teunter (1999) showed that the capital tied up in a non-serviceable product depends on the stock holding alternatives. Teunter et al. (2000) utilized simulation analysis to compare five different replenishment policies using system costs and return rate of non-serviceable items as performance measures. Teunter and Vlachos (2002) evaluated a production system with remanufacturing capabilities under different demand and return rates. Kiesmller (2003) addressed a stochastic recovery system with two stock points assuming separate production and remanufacturing lead times. The problem is simulated and analyzed under six scenarios considering different demand mean, variance and return ratios. Aras et al. (2006) studied the behaviors of a remanufacturing facility, as a stock point, under two replenishment policies. OptQuest search engine of Arena 9 is used for finding the best system parameters. Li and Zhu (2009) proposed a simulation optimization method based on an integrated hybrid genetic algorithm to optimize production planning and control policies for dedicated remanufacturing system. The simulation study was conducted in Arena 10, validated using data from an existing case. A simulation method was also used by Hellstrm and Johansson (2010) to investigate different control policies for managing the transportation of returnable products in a closed-loop supply chain. Recently, Alinovi et al. (2011) used a simulation approach to investigate an efficient return policy for a mixed inventory system with manufacturing and remanufacturing channels. They also used real data to validate the model. The classification and analysis of the published simulation optimization models can be found in the comprehensive reviews of Swisher et al. (2004) and Fu et al. (2005). The solution method that we propose in this paper is based on the integration of discrete event simulation with a heuristics search method called Variable Neighborhood Search (VNS). VNS which was first proposed by Glover (1986) and further studied by Hansen and Mladenovic (2001a) uses a set of pre-defined neighborhood structures Nk, (k = 1, kmax) with kmax > 1. What differentiates VNS from other local search heuristics is the systematic change of neighborhood structures during the search procedure. Many extensions have been proposed to the original heuristic enhancing its performance in handling large-size problems. Some of the most notable variants of VNS are Variable Neighborhood Descent (VND), Reduced Variable Neighborhood Search (RVNS), and Basic Variable Neighborhood Search (BVNS). More information about the types and applications of BVNS can be found in Mladenovic and Hansen (1997) and Hansen and Mladenovic (2001b). We propose a simulation-based Hybrid Variable Neighborhood Search (HVNS) method which is a BVNS that uses a tabu list to prevent the algorithm from being trapped at locality. Tabu list is the short-term memory of the search method that 148 H. Zolfagharinia et al. / Transportation Research Part E 67 (2014) 141161 prevents cycling. We use a tabu list with a circular fixed-length. Even though it has been shown that a fixed-length list cannot always avoid trapping at locality (Glover and Konchenberger, 2003), the probabilistic nature of neighbor selection in VNS algorithm enriches the anti-cycling mechanism. Moreover, using tabu-list helps the search method to improve its diversification feature. The proposed simulation-based HVNS algorithm (illustrated in Fig. 5) can be summarized in the following steps: Step 1. Initialization. A set of neighborhood structures is developed (Nk, k = 1, 2, ..., kmax). Step 2. Generate an initial solution x. Step 3. Construct a tabu list by adding the initial solution. Step 4. The solution is evaluated by calling the simulation model. Since the model complexity hinder finding a closed form solution for the proposed problem, a simulation model is used to assess the system total cost for any particular setting. Step 5. The followings are repeated until the termination criteria are satisfied. a. Set k = 1. While k 6 kmax, do the followings: b. Shaking. Generate a random solution at the kth neighborhood of x (x 0 e Nk(x)) . c. Local search. For the initial solution x 0 , a local search is applied by creating two probabilistic neighbors. The simulation model is called to evaluate the two generated neighbors and the best solution is selected as a local optimum, named x00. max 1. Determine Kmax 2. Generate initial Solution x 3.Create a tabu list 5a. K =1 Stopping Criteria met? Move to neighbour Update the tabu list 4. Call the simulation model to evaluate the solution End YES NO K Kmax K =K + 1 5b. Generate a random solution at the kth neighbour of x 5c. Call the simulation model for the local search 5d. Move or not K = 1 5e. Update the tabu list NO YES NO YES Fig. 5. The procedure of the proposed simulation-based HVNS. H. Zolfagharinia et al. / Transportation Research Part E 67 (2014) 141161 149 d. Move or not. Compare the obtained local optimum solution with the incumbent move. If the new solution shows improvement and does not belong to the tabu list, then x x00 and continue the search with Nk (k = 1); otherwise increase k by one unit. e. Update the tabu list if the movement decision is made in the step d. Additional details of the proposed solution method are outlined below. (I) Coding scheme: A simple coding scheme is adopted in which a vector A, S1, S2, S3, S4, S5 of integer values is defined as a feasible solution. An example of a feasible solution is illustrated in Fig. 6 where 5, 10, 40, 80, 100 and 30 are allocated to A, S1, S2, S3, S4 and S5, respectively. A controlling mechanism is also used ensuring A P 0 when generating a neighbor solution. (II) Initial solution generation: An initial solution is generated in a way that the inventory level A at serviceable stock point (at the beginning of planning horizon) is equal to zero (recall at the introduction stage the rate of demand increase is very slow) and maximum inventory level for each of 5 stages (A, S1, S2, S3, S4 and S5) equals to the average demands in corresponding periods. (III) Neighborhood generation: two neighborhood structures (Kmax = 2) are used in the algorithm: the first structure is as follows: A random number between 1 and 3 (with equal probability) is generated for each variable. If the generated random number equals one, a fixed value (called neighbor structure step) is deducted from the current value of that variable. If the generated random number equals two, the variable value remains unchanged. If the generated random number equals three, a fixed value is added to the current value of the variable. The neighbor structure step is a dynamic value in the sense that its magnitude changes during algorithm evolution. In the initial one-third of algorithm iterations, the neighbor structure step equals three. Its value is dropped by one unit in the second one-third and equals to one in the remaining one-third of iterations. This dynamic neighborhood generation approach provides the flexibility to better explore the solution space in the initial iterations, while enforcing less exploitation when fitter solutions found in the later stages of the search process. The second neighborhood structure is similar to the first structure with this difference that the value of only one variable changes at a time. (IV) Fitness evaluation: The fitness value of a solution determines whether to move to a neighbor solution. Simulation is used in the proposed HVNS to evaluate the fitness of each neighbor (Teunter et al., 2000; Teunter and Vlachos, 2002). Based on the problem characteristics, the system starts with A unit of products at each simulation run. With no delay (no warm-up period required), all statistics (including all cost components) are gathered right from the beginning until the end of product life cycle. Despite the complexities of finding analytical solutions for the proposed problem (see Section 4.2), simulating the system behavior is relatively straightforward if chronological events of the system are known (Fig. 4). A simulation model can hence make it to evaluate the policy parameters over the product life cycle within a very short model runtime. (V) Tabu list: Since the search mechanism of BVNS uses a probabilistic neighborhood structure, a simple circular fixed-length tabu list can effectively ensure the anti-cycling behavior and diversify the search. From experience (our observations on multiple test problems), the last ten solutions are stored in the tabu list. (VI) Stopping criteria: Different criteria can be used as the termination condition. The most common criteria include (1) stopping the algorithm after certain passage of time or number of iterations, (2) stopping when no more improvement occurs in the objective value for consecutive iterations, and (3) stopping when the value of objective function is better than a predefined value. It is a common practice to use a combination of these criteria in complex meta-heuristics algorithms (Glover and Konchenberger, 2003). In our experiments, we use a combination of the first two stopping criteria. 6. Implementation and computational results 6.1. The case company and decision scenarios Test scenarios are developed using realistic data obtained from a recognized provider of toner cartridges for laser printers in Australia. The case company is a member of the Australian Cartridge Remanufacturers Association and supplies schools, small business and several home users in South Australia. The purchasing- and remanufacturing-related values that we use to generate our test scenarios are located within the actual ranges provided by the case company. The planning horizon in our analysis is set at 24 months (24 units) comprised of 24 one-month units. The fixed parameters used in all test scenarios A S1 S2 S3 S4 S5 5 10 40 80 100 30 Fig. 6. A vector of integer values in a feasible solution. 150 H. Zolfagharinia et al. / Transportation Research Part E 67 (2014) 141161 include remanufacturing costs and lead time (provided by the regional remanufacturing plant), ordering, shortage and disposal costs (provided by the head-office in Adelaide) and transportation cost and lead time (provided by the strategic 3PL provider). The fixed cost of storage in the remanufacturable stock point is assumed to be negligible. The decision scenarios are developed using various combinations of critical input parameters. Five critical parameters have been identified from the existing literature: product life time (PLT), purchasing lead time (PL), return probability (RP), holding cost ratio (HC) and cost ratio (CR) (refer to Teunter and Vlachos, 2002). HC represents the ratio of holding cost of new items to the returned items and CR shows the ratio of purchasing to remanufacturing costs (Table 1). We developed 108 realistic decision scenarios with PLT and CR taking two values and PL, RP and HC taking three values each (i.e., 2 2 3 3 3 = 108). Obviously, the values of the aforementioned fixed parameters remain unchanged in all scenarios. The demand function presented in Section 3.1 is used to provide an estimate of the mean of demand in period t (kt). The continuous distribution is applied at discrete intervals. For this, 10,000 random numbers are generated and the relative frequency of each interval can be viewed as an appropriate estimation of that particular demand ratio to the total demand of the planning horizon. Here, the relative frequency of observed demand at each interval is used as a proxy of kt in period t. The proposed hybrid algorithm is coded in MATLAB 2010b from where the Arena 13 software is called to evaluate the fitness of each generated neighbor. SAS 9 software is utilized for variance analysis. We first examine the performance of the proposed solution method (HVNS) in the proposed 108 decision scenarios against OptQuest solver in terms of both solution quality and runtime. Then, the two inventory models (i.e., the basic single-stock model and the proposed model with two stock points) are compared in a set of test problems in situation when the product life cycle follows the revised Beta distribution (as depicted in Fig. 1). We next conduct a series of sensitivity analyses. First, the impact of variation in critical parameters is explored on the degree to which the proposed model outperforms a single-stock model. These critical parameters include product life time, purchasing lead time, return probability, holding cost ratio, and cost ratio. Finally, the impacts on the gains are shown when demand follows different distribution patterns during product life cycle. 6.2. Numerical results: HVNS algorithm vs. OptQuest solver OptQuest have been widely adopted in the existing literature as a general-purpose simulation and optimization software (Yun et al., 2011). It utilizes several meta-heuristics such as Genetic Algorithm and Neural Network. We use OptQuest to benchmark our HVNS algorithm. The stopping condition is set at the point when no improvement in the value of objective function is observed in 1000 consecutive iterations, or after 30 min program runtime, whichever occurs first. The number of replications (n) is determined to ensure achieving a narrow confidence interval for the average total system cost. For this, we use Eq. (25) from (Kelton et al., 2004): n n0 h h0 2 ; 25 where n0 is the number of initial replications, and h0 is the half width (i.e., half of confidence interval) resulting from initial replications. We hence need n replications to adjust the half width from h0 to h. In our experiments, a narrow confidence interval (within 2% of the average) can be obtained by setting the replications equal to 50 across the board. Numerical results for solving 108 decision scenarios using both OptQuest solver and HVNS are shown in Table 2. In all scenarios, except for one (scenario 8), HVNS outperforms OptQuest solver in terms of average total cost. The dominance in solution quality varies from zero to close to 3% improvement in average total cost, with an average dominance of 1.34%. HVNS also beats OptQuest in terms of runtime. As mentioned earlier, the maximum runtime for OptQuest was set to be 30 min. The average runtime of OptQuest in 108 test experiments was recorded to be about 22 min which is considerably longer than the average runtime of about 12 min for HVNS (see Table 2). This is primarily due to the fast convergence speed of HVNS. For example, in scenario 13, OptQuest solver reaches a solution with the average total cost of 32,391 in iteration 1575, while the hybrid algorithm hits 32,089 in iteration 70.

What did you learn from this research? (Be as specific as possible please!)