An Effective Hybrid Self-Adapting Differential Evolution Algorithm for the Joint Replenishment and Location-Inventory Problem in a Three-Level Supply Chain

The integration with different decisions in the supply chain is a trend, since it can avoid the suboptimal decisions. In this paper, we provide an effective intelligent algorithm for a modified joint replenishment and location-inventory problem (JR-LIP). The problem of the JR-LIP is to determine the reasonable number and location of distribution centers (DCs), the assignment policy of customers, and the replenishment policy of DCs such that the overall cost is minimized. However, due to the JR-LIP's difficult mathematical properties, simple and effective solutions for this NP-hard problem have eluded researchers. To find an effective approach for the JR-LIP, a hybrid self-adapting differential evolution algorithm (HSDE) is designed. To verify the effectiveness of the HSDE, two intelligent algorithms that have been proven to be effective algorithms for the similar problems named genetic algorithm (GA) and hybrid DE (HDE) are chosen to compare with it. Comparative results of benchmark functions and randomly generated JR-LIPs show that HSDE outperforms GA and HDE. Moreover, a sensitive analysis of cost parameters reveals the useful managerial insight. All comparative results show that HSDE is more stable and robust in handling this complex problem especially for the large-scale problem.


Introduction
As a multiitem inventory replenishment policy, the joint replenishment (JR) which can save the total costs by grouping multiitems in the same order had received numerous attentions [1]. While research on the classic joint replenishment problem (JRP) reached a saturation point, many extension versions had been proposed. These extensions can be divided into two branches: (a) relax the assumptions (such as storage, demand, and capacity limitation) of the classic JRP to simulate more realistic problems. The extensive literature review is available in [1]; (b) integrate more supply chain decisions (strategic or operational).
Silva and Gao [2] pointed out that three important decisions, facility location decisions, inventory management decisions, and distribution decisions, are highly related within a supply chain. Although independent policy always leads to a degree of suboptimization [3], these decisions usually are undertaken separately due to the complexity of integration.
The existing works of literature were mainly the any two of three decisions integration. Wang et al. [4] proposed a two-level joint replenishment and delivery problem (JRD) by extending the model of Qu et al. [5]. The comparative results verified that the hybrid differential evolution (HDE) algorithm is effective. Cha et al. [6] considered the JRD of a three-level supply chain. A heuristic algorithm and an intelligent algorithm were proposed to solve this JRD. An extension of this JRD by adopting the freight consolidation policy was proposed by Moon et al. [7]. All the papers mentioned above were the integration of inventory and distribution decisions within one warehouse or distribution center (DC). The cost-saving can be achieved by sharing fixed ordering cost and transportation cost and economies of scale. With the development of global purchasing and supply chain management, there had been a strong move towards integration of location and inventory decisions between multi-DC (Berman et al. [8]). 2 The Scientific World Journal Except Berman et al. [8] and Silva and Gao [2], the early studies on location-inventory problems (LIPs) had almost exclusively assumed the continuous-review ( , ) inventory policy. Teo et al. [9] formulated a LIP ignoring the shipping costs from the DCs to retailers. Daskin et al. [10] studied a LIP including location cost, inventory related costs, transportation costs from the supplier to the DCs, and local delivery costs from the DCs to the retailers. They proposed a Lagrangian relaxation solution algorithm for the proposed model. Shen et al. [11] considered the similar problem of Daskin et al. [10] by restructuring the model as a set-covering model and used column generation to solve the model. Based on the work of Daskin et al. [10] and Shen et al. [11], several works of literature had been published on this topic [3,[12][13][14].
The periodic-review ( , ) policy was introduced into the LIP by Berman et al. [8]. They pointed out that the ( , ) policy is easier to coordinate than the ( , ) policy and integrating the JR policy is a possible direction for future research. Silva and Gao [2] firstly proposed a LIP by considering the JR policy simultaneously. They solved this model using a two-stage method, the first stage, a Greedy Randomize Adaptive Search Procedure (GRASP) was used to determine the location decision; the second stage solved the JRP corresponding to the locations defined in the first stage. The proposed method is more suitable for the cases of specified number of opened DCs. This is the only one paper about the joint replenishment and location-inventory problem (JR-LIP). Thus it is meaningful to extend the JR-LIP and provide a new effective solving approach.
JRPs and uncapacitated facility location problems had been proven to be the NP-hard problems and they were rather hard to find effective algorithms [15,16]. Current solution approaches for LIPs included different heuristics based on the complex mathematical analysis of the models, that is, Lagrangian relaxation based solution [10,12], column generation algorithm [3,11], conic integer programming approach [17,18], and greedy randomized adaptive search procedure [2]. However the problem becomes much more complex than traditional LIPs because of the introduction of the JR policy (the ordering frequency of each DC and the basic ordering time needed to determine). It is rather hard to solve this JR-LIP effectively by traditional approaches. Firstly, available heuristics are too problem-specific and rather difficult to design. There exists no versatility approaches. Secondly, the enumeration is inefficient. It may take years to find a solution for a large problem size. On the other hand, due to the advantage of ignoring the mathematical property to search the optimal solution by starting from a feasible solution, intelligent algorithms had grown quickly in handling JRPs and JRDs. Results of similar studies illustrated that genetic algorithms (GAs) and differential evolution algorithms (DEs) are comparable suitable approaches [4,[19][20][21][22]. Although the existing works of literature hinted the superior of ant colony optimization (ACO) in solving combinatorial optimization, Wang et al. [22] pointed out that ACO was inferior in convergence rate for the similar JRD and more difficult to be expanded to complex problems.
The aim of this study is to propose a new and effective approach to handle the modified JR-LIP model based on the study of Silva and Gao [2]. The differences between our model and Silva and Gao [2] are as follows: (I) we assume that the maximum number of DC is known (additional variables were added in our model), the objective is to determine the number of DC to be opened, the locations of these DC, the assignment of customers, the replenishment frequency of each DC, and the replenishment cycle time to minimize the total system cost; (II) in Silva and Gao [2], the optimal number of DCs is obtained through contrasting the total cost by ADDING (or DROPPING) DC one by one. This method is not suitable for large-scale problems. The proposed hybrid self-adapting DE (HSDE) can directly find the optimal number of DCs through an intelligent search in a given space and is relatively easier to extend for large-scale problems. In fact, the HSDE can overcome the disadvantage of one-to-one competition used by classic DE and can avoid the manual parameter testing of mutation factor and crossover factor.

Results of benchmark functions tests and numerical examples
show the robustness of the HSDE.
The rest of this paper is organized as follows. Section 2 proposes the mathematical model. In Section 3, a DE-based algorithm is proposed to solve the JR-LIP. Section 4 is numerical examples and results discussion. Section 5 contains the conclusions and future research directions.

Assumptions and Notations.
A three-level supply chain consisting of multidistribution centers (DCs), an outside supplier, and multicustomers is considered. The item is ordered and collected by the DCs and then is distributed to customers. The objective is to decide the following policy: (1) how many DCs should be opened, and where to locate them; (2) how the customers are assigned to appropriate DCs; (3) how many and when to order, to minimize the total cost. The JR-LIP model is studied based on the following assumptions.
(i) Demand rates and costs are known and constant.
(ii) Shortages are not allowed.
(iii) Replenishment lead time is constant.
(iv) Each customer is only assigned to one DC, while other DCs cannot serve it.
(v) There is no limitation for the capacity of storage and shipment.
Only one item is considered in our model, DCs replenish their demand jointly. DC replenishes its item at every integer multiple of the basic cycle time and then delivers them to customer . Figure 1 gives a simple description about our model.
Notions used in the model are as follows: : index of DCs, 1 ≤ ≤ ; : index of customers, 1 ≤ ≤ ; : index of potential sites for DC, 1 ≤ ≤ ; : annual demand rate for the item per unit time at customer ; The Scientific World Journal : the variable to decide whether DC to be opened (decision variable), it is defined as : the location variable (decision variable), it is defined as : the assignment variable (decision variable), it is defined as

Mathematical Model and Analysis.
The total cost is composed of the fixed location costs of DCs, assignment costs of DCs, and replenishment costs of DCs. For a given problem, we assume that the maximum number ( ) of DCs is known and each site of customer is a potential site for a DC.
(a) Location Costs. The distance between customers is not considered in our model. Thus, besides the effect of replenishment policy, the main concern for making the location decision is the fixed opening cost and the distances from DCs to customers which are usually used to measure transportation costs in a supply chain network [2,4]. So the total location costs ( ) can be written as In (4), the first term is the fixed cost of locating DCs. The second term is the transportation cost of shipping from DCs to customers.
(b) Replenishment Costs. Similar to the classic JRP, the replenishment costs of DCs include ordering costs which consist of major ordering cost and minor ordering cost and inventory holding cost. The difference is the frequency of orders and the basic ordering cycle time at each DC is determined by the demand served by the DC. In turn, the demand served by the DC is a function of the assignment of customers to the DC. Thus the annual replenishment costs ( ) are formulated as where = ∑ =1 is the demand for the item by unit of time allocated to DC . The first term is the major ordering cost of DCs; the second term is the total minor ordering cost of DCs; the third term is the annual inventory holding cost of DCs in (5).
(c) The Objective. According to the above analysis, the annual total cost (TC) is The objective of JR-LIP model is min TC ( , , , , ) The goal of this problem is to find the optimal , , , , and to minimize the TC.

Problem Solving Methodology
Since all involved costs in the model are associated with the location of the DC, the thinking of solving methodology is converting (7) to a function with firstly. Assuming that the index of opened DC ( ), the customers assigned to the opened DC ( ), and the replenishment policy of the opened DC ( , ) are known, the problem we faced is where to locate these opened DCs ( ) to minimize the total cost. For simplification, we denote that A is the set of opened DCs and A is the customer set of DCs . Then the annual total cost can be rewritten as follows: For a given set { , , , }, the total cost of each opened DC at each potential site (TC * ) can be calculated easily. Sorting the total cost of each opened DC of all potential sites by ascending, the optimal * and TC * of DC is the first site in the list. If there are DCs with the same best location index, the optimal * and TC * is obtained by comparing at least combinations and at most ! combinations from the first site to the th site in the list. Finally, sum TC * and the optimal total cost TC * ( ) are obtained.

The Proposed Hybrid Self-Adapting DE Algorithm (HSDE).
Due to its easy implementation, quick convergence, and robustness, DE has turned to be one of the best evolutionary algorithms in a variety of fields [23][24][25][26]. However, the limitations on DE structure had inspired many scholars to improve upon DE by proposing modifications to the original DE fields [27].

The works of Literature on DE's Modification. Neri and
Tirronen [28] made a comprehensive study on DE's recent advances. They used twenty-four benchmark functions to survey the performance of eight DE-based algorithms. The comparative results showed that self-adapting parameters of DE (jDE) proposed by Brest et al. [29] were the effective and most simple improvement. Furthermore this modification also can avoid the manual parameter setting of and CR.
The DE-based algorithms referred in Neri and Tirronen [28] were the modification for DE structure and many scholars also devoted to integrate other algorithm's superior scheme into DE. The most common integration is the combination of DE and GA operations. He et al. [30] combined GA and SQP operation into DE which is rather complex to be carried out for complex application cases. To avoid the limitation of one-to-one competition in DE, Lin [31] integrated the roulette wheel selection into DE. Further, Wang et al. [4] proposed a hybrid DE (HDE) which adopted a simpler selection scheme of GA named truncation selection. The numerical results showed that HDE is effective and can easily be applied in practical problems. The newest integration algorithm of DE that appeared recently is quantum-inspired differential evolution algorithms (QDEs) The dimension of the specific problem

GenM
The maximum generation for evolution The target vector of individuals t in G generation

V
The mutated vector of individuals t in G generation The trial vector of individuals t in G generation Mutation factor of individuals t in G generation CR Crossover factor of individuals t in G generation [32][33][34][35]. Since a quantum system with qubits can represent 2 states simultaneously, the population size of the algorithm can be smaller, even one individual [36]. However, due to the complex encoding and decoding scheme, it was not easy for complex practical problems. Based on the above analysis, we propose a DE-based algorithm named HSDE by integrating the GA and selfadapting parameters of DE.

The Operations of HSDE.
The difference between HSDE and original DE is the structure of individual and selection operation. Table 1 lists the notations used in the intelligent algorithm. The details are discussed as follows.
(1) Individual Structure. In HSDE, the control parameters mutation factor and crossover factor are related to the individual and evolution generation, not constant for all individuals in the whole evolution process. Based on the study of Brest et al. [29], the individual structure of HSDE is where and CR are updated by the following formulation: where rand 1 , . . . , rand 4 are randomly generated number with uniform distributed between 0 and 1; 1 and 2 are constant values which represent the probabilities of parameters' update; min and max are the minimum and maximum , respectively.
The Scientific World Journal

The Procedures of HSDE for the JR-LIP.
Since the common operations of HSDE have been described in Section 3.1, we only focus on the procedures for the specific problem in this section, that is, initialization and decoding scheme, which are bolded in Figure 2.

Initialization.
In the stage of initialization, we should determine the parameters of HSDE, the representation of chromosome, and the encoding and decoding scheme for the JR-LIP.
The Parameters of HSDE. Since the mutation factor and crossover factor for HSDE are encoded in the individual, the parameters that need to be determined are population size ( ), the number of iterations (GenM), the lower bound ( min ) and upper bound ( max ) of mutation factor, and the probabilities of parameters' update ( 1 , and 2 ). As suggested in Brest et al. [29], we set min = 0.1, max = 0.9, and 1 = 2 = 0.1. The values of and GenM should be confirmed according to the practical problem as discussed in Section 4.
The Representation of the Decoded Individual. The proper representation of a solution plays a key role in the development of the HSDE. There are four parts included in the solution, the first part is the assignment information of customers; the second part is the replenishment frequency of DCs; the third part is the basic replenishment cycle time of DCs; the fourth part is the control parameters of HSDE. The dimension of the specific problem is = + + 1 ( customers' assignment information; DCs' ordering frequency; a basic ordering cycle time), and the total length of chromosome  equals to + 2. Figure 3 shows a decoded chromosome of twelve customers, three DCs.
The first part of the decoded individual represents that customers 1, 2, 4, 6, and 9 are assigned to DC 1 ; customers 3, 5, 8, 11, and 12 are assigned to DC 2 ; customers 7, 10 are assigned to DC 3 . Thus, the decision variables and can be confirmed by assignment information. The rest of solution  represents that DC 1 replenishes its item at every 2 interval; DC 2 and DC 3 replenish their item at each interval.
The Initialization. According to the representation of the decoded individual, the dimension of initial individual can be confirmed as + + 3. Then the values of individuals with + + 3 dimension are generated between 0 and 1 randomly and then are mapped into practical values (see Figure 3) through decoding scheme.

Decoding
The practical values of the last two parts of chromosome are between 0 and 1, since the value of gene can directly map into practical value.
From (14), we can see that the lower and upper bound is the key for decoding. It is obvious that the lower bound of opened DC and basic order cycle time are equal to 1. As to the upper bound, if the value is set too small, the optimal solution may be excluded; if the value is set too big, the search scope is enlarged which has directly impact on performance of algorithm. The maximum number of opened DCs ( ) is usually decided according to the number and scope of customers; the upper bound of replenishment frequency ( ) is set to 15 which is almost 4 times of the maximum optimal replenishment frequency obtained by Silva and Gao [2].

Numerical Examples and Results Discussion
To verify the performance of the proposed HSDE, three numerical examples were designed. Another two intelligent algorithms, GAs and HDE which had been proven to be effective approaches for solving JRPs and JRDs [4,19,20], were chosen to compare with it. Example 1 is numerical optimization problems with benchmark functions; Example 2 is used to compare three algorithms by different sizes of the JR-LIPs; Example 3 is designed to observe the impact of different cost parameters on the optimal decision.
The decoding scheme of GA and HDE is the same with HSDE. Besides the number of populations ( ) and maximum number of iterations (GenM) and the rest of common parameters for three numerical examples are set as follows: (a) for HSDE, the parameters had been given in Section 3.2.1; (b) for HDE, crossover rate CR is set to 0.3, is set to 0.6; (c) for GA, the probabilities of crossover and mutation are set to 0.9 and 0.1, respectively. The decisions for using these values are based on the experience from the works of literature [4,29,38].

Seven Benchmark Functions.
In order to verify the performance of HSDE, seven well-known benchmark functions are used in the following experiments. To assure a relatively fair comparison, the functions were selected according to their different properties [29]. Functions 1 -3 are unimodal functions, 4 is a step function, and 5 -7 are multimodal functions which appear to be the most difficult class of problems for many intelligent algorithms. The details about seven test functions are listed in Table 2 ( denotes the dimension and denotes the decision space of the problem).

Comparative Results and Analysis.
Three algorithms including HSDE, HDE, and GA are compared with two different dimensions, that is, the speed of convergence and the ability to obtain optimal solution. In all cases, the number of populations ( ) is 100; the maximum number of iterations (GenM) is 300 and 500 for = 10 and = 30, respectively. The average min ± standard deviations of each algorithm for different test functions over 10 independent runs are shown in Tables 3 and 4. The best results are highlighted in bold face. Furthermore, in order to have an intuitive understanding about the convergence speed performance of each algorithm, Figures 4, 5, 6, and 7 present the average fitness curve of two unimodal functions ( 1 and 3 ), step function 4 and multimodal function 6 , with thirty dimensions.

Results in Tables 3 and 4 and Figures 4-7 show that
(1) the ability to obtain optimal solutions of HSDE is better than HDE and GA, especially for functions with the large dimension (when = 10, the cases to obtain best performance of HSDE, HDE, and GA are 6, 5, 0; for = 30, the cases are 6, 2, 0, resp.); (2) although HSDE and HDE can both obtain the global optimum for 4 with thirty dimensions, the convergence curves presented in Figure 6 illustrate that HSDE has better convergence speed. In summary, the proposed HSDE has better performance for the largescale optimization problems.  The Scientific World Journal   Table 5 comes from Silva and Gao [2]. In all examples, the demand points were randomly located in a 50 × 50 square.

Results and Analysis.
We denote as the scale of problems. As for of the HSDE and HDE, Storn and Price [38] proposed that the population size ∈ [4 , 10 ]; Nobakhti and Wang [39] and Wang et al. [40] suggested that ∈ [2 , 20 ] is rational. According to their experiments, we set = 200, 300, and 450 for 5 30, 10 50, and 20 100, respectively. Table 6 reported the results of the average CPU times (Avg CPU times), the optimal total cost, the average minimum total cost for 20 times (Avg TC min ), and the ratio of finding the optimal total cost. Figures 8 and 9 showed the convergence trend of three algorithms (based on average fitness) for 5 30 and 10 50.
Results in Table 6 show that (1) both HSDE and HDE can find the optimal results with 100% for the smaller scale problem; (2) with the increasing of problem size, the performance of HSDE is better than HDE and GA whatever in CPU times and total cost. Furthermore, the running time of HSDE for 10 50 (99 seconds) is superior to Silva and Gao [2] in which the running time is more than 270 seconds when = 50. Convergence results (5,30) The average best value of TC  The Scientific World Journal 9   Section 4.2 mainly focuses on the performance of algorithm in handling different scales of JR-LIPs. The more detailed analysis about the impacts of cost parameters on the optimal decision is presented in Section 4.3.

Example 3: Impacts of Cost Parameters on Optimal Decision.
The size of problem in this section is 10 50 which is the middle scale of three problems in Section 4.2 and the optimal TC can be found at each running by HSDE and HDE. GA is not used in this section due to its inferior performance (the ratio to find optimal TC is 0%). Tables 7  and 8 show the computational results by varying the values of inventory holding cost and major ordering cost. From the tables we can see clearly the impacts of cost parameters on DC selection, replenishment frequency ( ), location costs ( ), replenishment costs ( ), and total running time of algorithms.
The comparative results of Tables 7 and 8 show the following useful conclusions.
(1) The robustness of HSDE is better than HDE.
(2) When inventory holding cost and major ordering cost vary from −40% to 20%, the optimal number of DCs is always three. Moreover, the location sites of DCs have little change (Δ = 20%), and the varied direction of replenishment cost is the same with the cost parameters.
(3) DCs have more motivation to order item jointly when the major ordering cost S increases ( 1 = 2 = 3 = 1 indicates that DCs replenish jointly in each order). This observation is consistent with the principle of JR policy.

Conclusions and Future Research
In this paper, we proposed an effective intelligent algorithm for a modified joint replenishment and location-inventory (JR-LIP) model. The objective of the JR-LIP is to determine the number and locations of DCs, the assignment decision, and replenishment policy to minimize the total system cost.
To handle this NP-hard problem effectively, an intelligent algorithm named HSDE is designed to solve the proposed model. To verify the effectiveness of HSDE, GA and HDE were chosen to be compared with it by benchmark functions tests and numerical examples. We can easily come to useful conclusions and managerial insight as follows.
(1) Results of benchmark functions tests show the good ability of HSDE in handling the large-scale problems.
When the dimension of test function is 10, HSDE and HDE have the same precision, while the dimension is 30, HSDE has the higher precision and faster speed than HDE. (2) The similar conclusion can be obtained from example 2. The rate of convergence obtained by HSDE and HDE is both 100% when = 30 and = 50, whereas when = 100, the rate of HDE is 0% and HSDE is 45%. (3) Example 3 illustrates the impacts of cost parameters on the optimal decision and reveals that when the major ordering cost is bigger, DCs have more incentive to replenish jointly for sharing related costs.
All numerical examples verify that the HSDE is an easy and effective algorithm to handle the JR-LIP. To our best knowledge, this is the first time to use the improved DE-based algorithm to solve this NP-hard problem.
In our model, the demand is constant which generated from a uniform distribution, and there is no capacity limitation. These assumptions are not required. The future direction about this problem includes relaxing the assumption to match real-world scenario and looking for more quickly and efficiently solving method.