Taguchi Analysis for Improving Optimization of Integrated Forward/Reverse Logistics

The distribution–allocation problem is known as one of the most comprehensive strategic decisions. In real-world cases, it is impossible to solve a distribution–allocation problem completely in acceptable time. This forces the researchers to develop efficient heuristic techniques for the large-term operation of the whole supply chain. These techniques provide near optimal solution and are comparably fast particularly for large-scale test problems. This paper presents an integrated supply chain model which is flexible in the delivery path. As solution methodology, we apply a memetic algorithm with a novelty in population presentation. To identify the optimum operating condition of the proposed memetic algorithm, Taguchi method is adopted. In this study, four factors, namely population size, crossover rate, local search iteration and number of iteration, are considered. Determining the best level of the considered parameters is the outlook of this research.


Introduction
Supply chain networks describe the flow and movement of goods by linking several facilities such as plants, distributions, and retailers in forward flow [1,2] and collec- integer linear programming, which traditional methods fail to solve in acceptable time, particularly when we are facing with large size problems. To overcome this problem, nontraditional solutions such as heuristics and metaheuristics [8,11,[15][16][17] have been proposed in recent years. Within this paper, we present a memetic algorithm with a new chromosome representation as well as updated operators. Each algorithm has some parameters. Some information regarding these parameters would be useful to improve the performance of the results. In this study, Taguchi method is applied to improve operation condition. Determining the best level of parameters is the main contribution of this work.
In the last several decades, growing attention is being paid to the use of Taguchi method to find the best setting of parameters involved in evolutionary algorithms. Chouhan et al. [18] introduced a novel approach that had the ability in collecting endof-life and end-of-use products from the end-users and in reducing the overall cost of the supply chain. To assess the overall cost of the supply chain, a mixed-integer linear programming model has been formulated, while the Taguchi method was used to obtain the best combinations of the algorithm parameters.
Gümüs et al. [19] analyzed and extended the Taguchi method to tune the parameters of a memetic algorithm for cross-domain search. In their work, they investigated the degree to predict the same good parameter setting faster by using a reduced time budget and showed that it was possible to predict good combinations of parameter settings with a much reduced time budget.
Alavidoost et al. [20] proposed a fuzzy mixed-integer linear programming optimization model for an integrated supply, production, and distribution supply chain of multi-product, multi-level, multi-plant. They proposed a model aimed to optimize customer satisfaction as well as chain supply costs, simultaneously, which was solved by some popular multi-objective metaheuristic algorithms. They developed a Taguchi method to calibrate as well as control the parameters in algorithms. Additionally, they proposed a new framework for ranking the algorithms with the use of the fuzzy VIKOR method.
A mathematical model for a fixed-charge solid transportation problem for a twostage supply chain network is introduced in Akbari et al. [21] They presented a fixed-charge solid transportation problem in a two-stage supply chain network, employing conveyances in each stage of the supply chain. They used three metaheuristic algorithms, namely GA, EM, and CSS, in accordance with priority-based encoding, tackling this NP-hard problem, and employed the Taguchi parameter design method to generate the best solutions. They showed that the best result can be obtained by the use of CSS.
Zahraee et al. [22] used dynamic simulation and the Taguchi method to design a robust blood supply chain. To do the Taguchi method, four main controllable factors including arrival rate of donors, maximum inventory level, minimum inventory level, and blood delivery policy, and one noise factor have been chosen. Their results showed that a robust blood supply chain can be achieved, while all main factors are located at a high level except for factor B.
Vinay and Sridharan [23] developed a solution methodology using ant colony optimization (ACO) for a distribution-allocation problem in a two-stage supply chain with a fixed cost for a transportation route, while Taguchi method for robust design was adopted for finding the optimum combination of parameters of the ACO algorithm. Applying the Taguchi method led to significant reduction in the total number of experiments to be carried out for finding an optimum level of parameters, namely from 5 120 to 125, i.e., 25 experiments each with five replications.
The paper is divided into six sections, and the rest of this paper is organized as follows. In Sect. 2, we describe the proposed network and formulate the considered problem into a mixed integer linear programming. Then, after a brief explanation regarding the solution methodology in Sect. 3, parameter analysis methods are presented in Sect. 4. Section 5 provides Taguchi analysis for improving optimization of the model, and results are shown. Section 6 concludes the paper and suggests future plans.

Description for Integrated Forward/Reverse Logistics Network
In this section, we formally define the considered problem and formulate the proposed supply chain network. The model deals with the minimization of transportation and operation cost of the proposed flexible integrated forward/reverse logistics network in order to determine the optimal capacity and number of each node as well as distribution. For this purpose, the problem is stated as follows. Consider G = (N , E) as a digraph, where N is the set of nodes and E the set of edges. The cost involved in the proposed model is denoted in two types: (1) fixed cost for node i ∈ N by c i and (2) unit transportation cost on edge (i, j) ∈ E by c i j . Furthermore, y i and x i j represent decision variables, where y i ∈ {0, 1} and x i j ∈ N 0 and a i and b i , vectors of known coefficients for the involved constraints. These variables indicate whether a stage i ∈ N is used, and which quantity is transported from node i to j, respectively. In accordance to above description, the following mixed integer linear minimization problem arises min Next, we specialize this model to reflect the problem properties.

Assumption of the Proposed Model
The previously described flexible forward/reverse logistics network setting represents an integrated supply chain with seven echelons consisting of suppliers S, plants P, distribution centers Dc, retailers R and customers C in forward flow, as well as collection/inspection centers Co and disposal centers Di in reverse flow, cf. Fig. 2 for a schematic sketch. We like to point out that in accordance to Fig. 2, we consider a hybrid manufacturing-recovery-recycling facility as well as a hybrid collectioninspection facility. Establishing several facilities at the same location can decrease the price of the whole network in comparison with separated design.
To adapt problem (1), we impose the following assumptions: -There are seven echelons: suppliers, plants, distribution centers, retailers, customers, collection/inspection centers, and disposal centers.
-There are no edges between facilities of the same stage, the delivery graph is complete in forward flow and the return graph is simple, i.e., cost. -The un-recyclable returned products will be sent to the disposal center. The remaining products are returned to the same plant. -The required recycled materials are assumed to be of the same quality as the raw materials bought from suppliers and any plant chooses the raw material from the collection/inspection centers over suppliers. -Customers have no special preference. It means, price is the same in all facilities.

Notation
To support the presentation of the proposed mathematical model, we first provide a verbal description of the model as follows: Cost = transportation costs + Fixed opening cost. ( The notation given in Tables 1, 2, 3 and 4 is used in the formulation of the MILP model.

Mathematical Formulation
Utilizing the notation of Tables 1, 2, 3 and 4, the structure of the MILP model can be presented as follows. The objective of this model is to minimize the total cost of the proposed supply chain. The first part of the cost function represents the variable    costs, which are given by the direct multiplication of the transportation cost and the quantity transferred from the origin to the destination. The fixed costs of the logistics chain manufacturing facility is considered as the second part of cost function. This reveals the cost function as follows: To shorten notation, we omit the bounds on the parameters in the following presentation of the model constraints: For one, the capacities in each node is limited. Constraint numbers 4, 5, 6 and 7 show limitation of the capacity for suppliers, plants, distribution centers and retailers in the forward flow. Specifically, constraint number 4 is shown bound of products supplied by supplier i, constraint number 5 bounds of products produced by plant j, constraint number 6 capacity of distribution center k and constrain number 7 capacity of retailer l. Constraint number 8 and 9 shows limitation of the capacity for collection/inspection and disposal centers in the backward flow. Specifically, constraint number 8 is shown capacity of collection/inspection center n and constraint number 9 capacity of disposal center o.
Secondly, this supply chain network is conservative so the demands of customer have to be satisfied.
Additionally, according to assumption, only a fraction p return m is returned by customers, and fraction p disposal n of the returned products has to be disposed off. Apart from these expectations, the supply chain network subjected to the law of the flow conservation in-flow and out-flow in each node must be identical. These conditions reveal: Specifically, constraint number 11 means that the sum of the total flow of recoverable products from collection/ inspection centers and the flow of the products from suppliers to a plant is equal to the flow of the products from that plant to distribution centers, retailers, and customers in the normal delivery, direct delivery, and direct shipment.
O o=1 X CoDi Finally, we require the decision variables to be nonnegative or binary, respectively.
The resulting optimization problem is NP-hard. To approximate a respective solution, we apply a memetic algorithm.
linear programming (MILP), suitable solution in acceptable time cannot be carried out through standard methods. Therefore, we present a memetic algorithm, which belongs to the class of metaheuristic algorithms as the solution methodology [26]. The performance of memetic algorithm is attached to two main issues, "chromosome representation" and "Memetic operators" [27].

Proposed Memetic Algorithm (An Overview)
According the above explanation, "chromosome representation" and "Memetic operators" applied in this study is presented briefly in following.

Chromosome Representation
In this study, we are using the extended random path-based direct encoding method to generate the first population [28]. By using this method, the integrated logistics problem can be integrated into a memetic algorithms by using direct encoding. At the same time, we captured the complexity of a full graph and reduce the size of encoding and thereby computational time by developing a two segment approach. The delivery path for each customer is shown in the first segment. Also, the guide information regarding plant as well as distribution assignment is existing in the second segment. As an illustration, the representation of the extended random path-based direct encoding method in two segments is shown in Fig. 3.
The different color in this sample of chromosome representation is considered as a sample of gene unit that shows the delivery path for customer M. Since there are three different delivery paths, the second segment is added to clarify which way is selected. Each gene in the second segment is assigned to an integer in the range of {0, 1, 2} for plant and {0, 1} for distribution center. The detailed decoding procedure is shown in Fig. 4. The procedure of encoding by extended random path-based direct encoding is shown in Algorithm 1.
It should be noted that this encoding approach may result to infeasible solution. Therefore, if the demand of a depot exceed the capacity of a source, the depot will be connected to another source with sufficient supply which the transportation cost is minimum.

Memetic Operators
In this work, the cost function is implemented as the fitness value, which is related to the objective function obtained by (1). This value is assigned to the decoded chro- The proposed memetic algorithm in this study has three operators including: crossover, local search and selection. For the crossover method, we select the two point crossover method which keeps the feasibility of the chromosomes. Figure 5 illustrates the two point crossover method.
Local search as one of the key factors of the memetic algorithm has an important role to explore the search space deeply. As it is impossible to swap each gene in the proposed chromosome, the applied local search technique in this paper has been updated according to the characteristics of the chromosomes. For local search technique we apply hill climbing-based random mutation method to generate more adapted offspring, reduce the chance of producing infeasible solutions and improve (first segment) Step 2 (second segment, plant delivery path) 10 for i = 0 : J − 1 do 11 ch k [7 * M + i] ← random(0, 2) 12 end for Step 3 (second segment, Dc delivery path) the chance of finding better solutions [28]. This adaptive local search engine is able to adjust the number of local search iterations dynamically during the process according to the size of the test problem. The following four major steps describe this method in details: -Specify a position randomly to apply mutation for a given incumbent solution.
-Create a neighborhood by replacing the selected gene with feasible possibilities using random mutation operator [29]. The procedure of random mutation operator is simply selecting one position in the current individual and exchange the related selected gene by a regenerated randomly new number [30]. -Select the best member of neighborhood according to the fitness value as the offspring. Figure 6 illustrates an example of the local search based-mutation method, while the number of the chosen gene in the whole network is 4. -Repeat this procedure according to the number of local search iterations. This number is considered as an input parameter of MA, and was first proposed by Hart [31]. It refers to the number of iterations of local searches during each generation cycle. For different problems, and according to the size of problems, the values of local search iterations need to be changed [32]. In this work, an adaptive local search engine is applied to adjust the number of local search iterations dynamically during the process according to the size of the considered test problems. As the number of local search iterations for the test problems should be determined by their size, we set the number of local search iterations increasingly according to 542 E. Behmanesh, J. Pannek After applying local search operator, new neighborhood solutions are obtained. If the fitness value of one of the new solutions is better than the old one, the latter is replaced. An example of the applied local search is shown in Fig. 6. The local search operator is repeated according to the number of local search iteration. For the selection part, the well-known roulette wheel method is selected [33]. The strategy applied in this method is a probabilistic selection based on fitness.
The procedure displayed a flowchart in Fig. 7 to show the steps of the solution methodology clearly. We have included different colors that shows the main focus of this study.

Parameter Analysis Methods
There are many methods for studying the process parameters [34]. "One factor at a time" and "Design of Experiment (DOE)" as the most popular approaches are implemented by researchers in this field.
-One factor at a time "One factor at a time" also known as "monothetic analysis" is a method involving one specific condition instead of multiple factors simultaneously. In the following, the experiment are repeated by changing any other one factor. -Design of Experiment "Design of Experiment (DOE)" is a statistical approach capable to determine the relationship between factors involved in a complex and multi-variable process with few trials. Although some researchers have shown that "One factor at a time" can be more accurate some times [35], it requires more runs for the same estimation in most cases. There are two major approaches to DOE [34]: -Full factorial design A "full factorial design" also known as "fully crossed design" is able to consider two or more factors. Each factor with a number of possible levels takes all possible combinations, generated by those levels across all those factors. This method allows to record the effect of each factor on the response variable individually. Using this method is time-consuming, particularly for increased number of factors. By applying this method with k factors and l levels for each factor, l k trials are needed. -Taguchi method According to the statistical theory, not every combination of parameters need to be checked [36]. In this regard, Taguchi [37] came up with a fractional replicated designs called with his name using orthogonal arrays [38].
In this method, a small set from all the possibilities is selected. By a comparison between these two methods, it is found that Taguchi's method is more efficient than "full factorial design" approach due to same results with lesser number of experiments [34]. As Taguchi's method found more efficient [34], we applied in this research.

Improving Optimization Using Taguchi
In order to improve operation condition of the memetic algorithm regarding the selected parameters, we applied the following steps from the Taguchi analysis [34]: 1. Identify the objective function to be optimized. 2. Identify the independence factors and number of levels for each. 3. Select a suitable orthogonal array and construct the matrix of simulation. 4. Conduct the matrix of simulation. 5. Examine and analyze data and predict the best parameter levels. 6. Confirm by simulation. To apply the presented steps, we generated five test problems. Since the framework of the proposed logistics network in this work is not the same as in previous studies, all parameters were generated randomly using uniform distributions as given in Table 5.
Each column in Table 5 shows a number of facilities involved in each test problem separately, and also in total, optimal cost obtained by branch-and-bound algorithm from LINGO15 and number of possible routes to show how complicated they are. The number of decisions variables are 209, 234, 468, 1 006 and 1 780, respectively. Other parameters are generated randomly using uniform distributions shown in Table 6.

Identify the Objective Function to be Optimized
Taguchi suggested use of loss function for measurement. It is recommended to transform value of the loss function into a signal-to-noise ratio S/N = S N R [38]. Since the signals have a very wide dynamic range, they are mostly expressed using the logarithmic decibel. Number of iteration 20 30 40 where f i denotes the fitness value for run number i and n define the number of repeated experiments. Three different categories are classified in order to analyze the S/N ratio: "nominal is better," "larger is better" and "smaller is better", [39]. As the aim of this study is minimizing the total cost, "smaller is better" is considered. The best results are obtained when the total cost is minimum. This condition corresponds to the highest S/N ratio.

Identify the Independence Factors and Number of Level for Each
According to a comprehensive review [40][41][42], there are four independence factors that need to be investigated regarding the proposed algorithm. These factors are population size, crossover rate, local search iteration and number of iteration. As all of these factors have a noticeable effect on the results of the proposed model, it is important to assess all for the efficiency criterion [42]. Three different levels for each factor are considered according to the information presented in Table 7.
For population size, we considered the numbers which are the most common numbers according to the literature by respect to arithmetic progression [40]. For crossover for sure, it should be in the range of [0, 1], and we considered 0.3, 0.5 and 0.7 to cover low rate, medium rate and high rate. For local search iteration and also number of iteration according to a comprehensive experimental research, we realized that 80 percent of progress happen in the bound we considered for the parameters [42].

Select a Suitable Orthogonal Array and Construct the Matrix
Taguchi method is dealing with two main tools including "S/N ratio" and a special set of arrays called "orthogonal array". Compared to many statistical designs, orthogonal array shows a significant success. This method allows to reduce the number of experiments, while a more reliable estimation for parameter optimization can be reached. Previous studies presented many orthogonal arrays according to the problem they are facing. Still, selecting a proper orthogonal array is difficult as the reported orthogonal arrays in the literature are not able to cover all possibilities [43]. To overcome this problem, some standard orthogonal arrays have been designed according to the number of factors and levels. In order to choose an appropriate orthogonal array, the minimum number of experiment needs to be carried out first, based on the total number of freedom given as [36]: where n f denotes the number of factors and L define the number of levels of each factor. For the interested case, there are four different factors with three levels each, and the degree of freedom is obtained 1 + (4 * 2) = 9. It means at least 9 runs need to carry out, and another selection of orthogonal array with more runs may be used to increase accuracy. Here, we utilized Minitab software to find an automatic design of the orthogonal array as well as for analysis of experiments. Minitab software suggested using L 9 orthogonal array, i.e., 9 instead of 3 4 = 81 experiments. The details of the proper orthogonal array is shown in Table 8:

Conduct the Matrix Simulation
In accordance with the presented orthogonal array mentioned in Table 8, simulation were carried out with their factors and levels. The scheme with the specified values of the factors is given in Table 9. Each of these nine experiment is repeated 10 times to generate more reliable averaged results.

Examine and Analyze Data to Predict the Best Parameter Levels
According to the obtained results presented in [42], we realized that all the proposed parameters are important and needed to be investigated further. Since we are dealing with a strategic problem, there is no time limitation. But still having some information regarding the best values of the parameters can help us to improve our results, applicable to other types of proceeding when time limitation is considered in advance. On the other hand, as there is no similar case in the literature, here Taguchi's method is adapted to this aim. As mentioned before, we are dealing with a strategic problem, so there is no time limitation. In the case of finding optimal parameters, we need to consider a specific time, otherwise increasing each parameters can improve  the results. In these regard, the following steps are applied. First, we need to fix an error percent which is acceptable for most of industrial players and logistics problems. This number has fixed as 0.05 ( [27]). Then, we need the specific value of parameters that lead us to an error percent equal or less than 0.05 percent which are available in [42]. In the next step, a high number of iteration is assumed. This number should be more that the required number of iteration that is needed in practice. Using this way, we are able to find the average number of iteration which is required to let the algorithm find results with the preselected error percent. Afterward, the specific time that these results were obtained is recorded. Finally, the goal is discovering the best value of parameters once the algorithm is forced to find out solution in half of the recorded time using Minitab software. Figures 8, 9, 10, 11 and 12 display the effect of the parameters on the objective function value and S/N ratio, respectively.
Since S/N ratio has the "large is better" characteristic, the analysis of the results leads to the population size at level 2, crossover rate at level 3, local search iteration at level 2 and number of iteration at level 2.

Confirmation Simulation
Once the best level of the parameters were chosen according to the mentioned condition, the final step is needed to be implemented. The idea of the confirmation experiment is to validate the results obtained from the experiments. A confirmation experiment was executed 10 times with the obtained combination of parameters. The actual value of cost obtained by MATLAB was calculated based on the experiments. The predicted value obtained by Minitab software is computed according to statistical data. Our implementation was written in MATLAB R2015b and run on the PC with Intelr CoreTM i5 2.40GHz with12 GB RAM. All results are provided in Table 10.
The results show that the predicted value of cost obtained with the best level of parameters is very close to the value of the actual cost. This validates the presented Taguchi method for predicting the optimal levels of the parameters.

Conclusions and Future Studies
This research focuses on an integrated, flexible, multi-stages, single-period and single-product supply chain network that is formulated as a mixed integer linear programming. A memetic algorithm is adapted as the solution methodology using a new chromosome representation. We applied the Taguchi method to design the optimal parameters of the proposed memetic algorithm. From the confirmation experiments, a good agreement between the predicted results and the actual results is observed. Using Taguchi method in determining the optimal level of parameters is a new contribution of this study. In the present work, only four parameters with three different levels were selected. More levels can also be considered in order to obtain better values. Different crossover, local search and selection method can be added in the analysis as well.