Abstract
As a practical inventory and transportation problem, it is important to synthesize several objectives for the joint replenishment and delivery (JRD) decision. In this paper, a new multiobjective stochastic JRD (MSJRD) of the onewarehouse and retailer systems considering the balance of service level and total cost simultaneously is proposed. The goal of this problem is to decide the reasonable replenishment interval, safety stock factor, and traveling routing. Secondly, two approaches are designed to handle this complex multiobjective optimization problem. Linear programming (LP) approach converts the multiobjective to single objective, while a multiobjective evolution algorithm (MOEA) solves a multiobjective problem directly. Thirdly, three intelligent optimization algorithms, differential evolution algorithm (DE), hybrid DE (HDE), and genetic algorithm (GA), are utilized in LPbased and MOEAbased approaches. Results of the MSJRD with LPbased and MOEAbased approaches are compared by a contrastive numerical example. To analyses the nondominated solution of MOEA, a metric is also used to measure the distribution of the last generation solution. Results show that HDE outperforms DE and GA whenever LP or MOEA is adopted.
1. Introduction
The joint replenishment problem (JRP) is a practical inventory problem of a group of products that can be jointly ordered from a single supplier (Goyal [1]; Wang et al. [2]), which can help save the ordering costs and inventory holding costs. According to the characteristic of demand, the existing study of JRPs can be divided into two categories: (1) constant demand; (2) stochastic or dynamic demand. An extensive literature review is available in Khouja and Goyal [3] and Narayanan et al. [4]. Many scholars also discussed more realistic JRPs (J.M. Chen and T.H. Chen [5]; Axsäter et al. [6]; Hsu [7]; AbdulJalbar et al. [8]).
Many companies have realized that a joint replenishment and delivery scheduling (JRD) policy can result in considerable cost savings. But the literature on the JRDs under supply chain environment is limited. A stochastic JRD of the onewarehouse, retailer system has been formulated (Qu et al. [9]). Wang et al. [10] studied the same JRD but reduced the decision variables using specific mathematical method and provided a new differential evolution algorithm. Sindhuchao et al. [11] studied the coordinated inventory and transportation decisions with the vehicle capacity limitation in an inbound commodity collection system. Chan et al. [12] addressed issues in scheduling of the multiitem, multibuyer, and single supplier system. Cha et al. [13] handled the JRD of the onewarehouse and retailer system in which the warehouse supplies items from the supplier and delivers them to retailers. Moon et al. [14] modified the model of [13] by utilizing a consolidated freight delivery policies. Wang et al. [2] extended the JRD model of Cha et al. [13] under fuzzy environment and used the widely used signed distance method to ranking fuzzy numbers. Wang et al. [15] studied the JRD with deterministic demand and fuzzy cost using the graded mean integration representation and centroid approaches to defuzzify the total costs. A common limitation in all the literature studies mentioned above is that they only consider a single objective.
However, managers are usually faced with complex multiobjective optimization problems (MOPs) in reality. For the JRD policy, it is necessary to decrease the total cost while improving the service level. Although there are several papers that studied multiobjective inventory models (Roy and Maiti [16]; Rong et al. [17]; Islam [18]; Wee et al. [19]), no study on the multiobjective JRD can be found.
For MOPs, direct comparison among the solutions is very difficult because of the different measurements between each contradicted target. In this study, total cost and service level are obviously two contradictory targets: reducing total cost may result in the decline of service level and vice verse. So we should coordinate two targets. Different from a singleobjective optimization problem which has unique optimal solution, an MOP has a set of optimal solutions called Pareto optimal solutions. Due to the characteristics of MOPs, they are much more complex and the key is to find an effective method to obtain Pareto optimal solutions. Unfortunately, the classical JRPs and JRDs are already NP hard problems (Arkin et al. [20]), and the multiobjective makes the JRDs become much more difficult to handle.
Many linear or nonlinear weighted methods (Rong et al. [17]; Islam [18]; Wee et al. [19]; Roy and Maiti [16]) were used to convert the multiobjective to a single one in the existing studies. These methods undoubtedly provide one easy way to deal with the multiobjective JRD model. However, these approaches do not solve the MOPs intrinsically, since the solutions for MOPs are multiple rather than one. On the other hand, multiobjective optimization methods based on Paretobased MOEAs are widely used, such as multiobjective genetic algorithm (MOGA) (Aiello et al. [21]), nondominated sorting genetic algorithm (NSGA) (Lin [22]), and strength Pareto evolutionary algorithm (SPEA) (Zitzler and Thiele [23]; Sheng et al. [24]).
In recent years, several MOEAs based on Pareto differential evolution (DE) were utilized to solve MOPs. SantanaQuintero and Coello [25] presented a DEbased multiobjective algorithm using a secondary population and the concept of dominance. The performance of the proposed algorithm was also compared with NSGAII and MOEA. Qian et al. [26] proposed a memetic algorithm based on DE (MODEMA) for multiobjective job shop scheduling problems. Qian and li [27] proposed an adaptive DE (ADEA) and the results of five test functions showed that the ADEA was very efficient to find out the true Pareto front. However, the majority of the above studies focused on the effectiveness verification of the algorithms and always analyzed standard testing functions and ignored the practical applications of them in MOPs. Due to the existing unpredictable and uncertain factors of inventory management, it is difficult to convert shortage rate/quantity to shortage cost. Therefore, discussing shortage rate/quantity independently is meaningful.
The aim of this study is to propose a new multiobjective stochastic JRD (MSJRD) model including two minimum objectives, that is, the total cost and shortage quantity. The main difference between this study and [10] is that two objectives are handled simultaneously. Moreover, effective approaches are provided to handle this MSJRD. Having the successful applications in engineering management, as well as the effectiveness of DE in solving MOPs and JRPs/JRDs (Wang et al. [28–30]), linear programming (LP) and MOEAbased approaches using DE are provided. Results of an example show that the proposed hybrid DE is more effective than the original DE and GA whatever LP or MOEA method is used.
The rest of this paper is organized as follows. Section 2 describes the proposed multiobjective stochastic JRD model. Section 3 introduces the hybrid DE. Section 4 presents two approaches to solve the MSJRD. Section 5 contains numerical examples and results. Section 6 discusses conclusions and provides future research directions.
2. Formulation of the Proposed MSJRD Model
Consider an enterprise has a center warehouse in a proper place and several suppliers in decentralized locations. The center warehouse jointly replenishes items from its suppliers according to the market demand or historical data. Then, the center warehouse will collect replenished items from suppliers. In this situation, the center warehouse can jointly determine the replenishment and distribution policy to obtain the optimal decision.
Refering to the model of Qu et al. [9], the differences between the JRD policy and typical JRP can be concluded as follows: (a) the JRP supposes deterministic demand, while our model considers stochastic demand and allows for shortage; (b) the JRP just discusses the replenishment policy, while our model analyses not only replenishment but also transportation decision; (c) the JRP is a single objective model, while our model is a multiobjective model.
For the MSJRD, reducing the related total cost as well as decreasing shortage quantity should be considered simultaneously. Therefore, the first target is to minimize total cost which consists of replenishment cost, inventory holding cost, and distribution cost. Minimizing the shortage quantity is the second target.
2.1. The First Target: Total Cost
The total cost includes inventory holding cost, replenishment cost, and distribution cost. Distribution cost involves the stopover cost at suppliers and cost related to distance. The infinite capacity of vehicle assumption makes the distribution problem a traveling salesman problem (TSP).
Inventory cost can be given as
In (1), the first term is deterministic inventory and the second one is safety stock. The following replenishment cost is the same as the classical JRP:
When are given, taking the least common multiple of we can get integer . It means that the replenishment and distribution behavior will be repeated every period. For example, if , then . That is to say, in period 1, all items should be replenished; in period 2, items 3 and 4 need to be replenished; in period 3, items 2, 3, and 4 should be replenished; in period 4, items 3 and 4 should be replenished. It is obvious that in period 5 the situation is the same as in period 1. The cycle period is . Therefore, a limited horizon with periods is used to calculate annual distribution cost.
Jointly replenishment of items can take the advantage of scale economies not only in the replenishment but also in distribution process. Four items specifically should be replenished in period 1 and it is advisable for the center warehouse to traverse suppliers that supply these items at one time. The shortest path in each period is considered so distribution cost can be obtained by solving a travelling salesman problem (TSP). The distribution cost is where
Therefore, the first target can be summarized as follows:
2.2. The Second Target: StockOut Quantity
In this study, the demand is assumed to follow normal distribution at interval of unit time, which is widely used in the literature (see Qu et al. [9]; Eynan and Kropp [31, 32]). This assumption means, that given and for item , the demand will follow normal distribution over the interval of length with the expectation , variance and probability density function .
With a periodic replenishment policy, stockout could occur any time during replenishment intervals as long as the real demand exceeds maximum inventory level . The total annual stockout quantity is where .
2.3. The MSJRD Model
The whole multiobjective model can be written as where .
The goal of this multiobjective model is to find out the optimal , , and to simultaneously minimize the total cost and stockout quantity and thus to achieve Pareto solutions in which two objectives can be balanced. Two targets have different units of measurements and it is usually difficult to convert the shortage quantity to stockout cost. In addition, they are often in conflict with each other, that is, decreasing shortage quantity may result in cost increasing.
MOPs are much more complex but closer to reality. Several traditional mathematic methods are used for solving multiobjective models, such as linear programming, goal programming, and analytic hierarchy process. However, they are successful only in small scale problems. Mathematic methods are too complex and too time consuming to solve large scale problems. In the following, we provide two common approaches based on an HDE to deal with the proposed MSJRD. Then a numerical example and comparative study between the proposed LP and MOEA are presented.
3. The Hybrid Differential Evolution Algorithm (HDE)
3.1. The Classical DE
DE has been described as an effective and robust method to optimize some wellknown nonlinear, nondifferentiable, and nonconvex functions. 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 (Wang et al. [33]; Cui et al. [34]). DE contains three operations: mutation, crossover, and selection.
3.1.1. Mutation
The mutation operation creates a new vector by adding the weighted difference of two random vectors to a third one. For each target vector , the mutated vector is created as follows:
In (7), , , and , are three serial numbers of vectors, which are randomly generated with different values and none of them equals . Three vectors , , and will be selected from the population for mutation operation when , , and are confirmed, is a scaling factor and is the current number of iteration.
3.1.2. Crossover
A trail vector is created by mixing the mutated vector with the target vector according to the following formula: where represents the th dimension; is randomly generated from 0 to 1; is a randomly selected integer to ensure the effect of mutated vector; CR is the crossover probability and it is very important for DE since the larger CR is, the more contributes to .
3.1.3. Selection
The selection operation is implemented by comparing the trial vector (obtained through mutation and crossover operations) with the corresponding target vector. For example, to minimize the function, the next generation is formed by where the is the fitness function of DE.
Here, an example is given to illustrate the three operations mentioned previously. For the current number of iteration and target vector , suppose that random generated numbers , , and are 23, 40, and NP respectively, and we obtain the following:
Mutation: if , the mutated vector can be obtained by (7) as follows:
Crossover: if , (here ), and vector , the trial vector can be obtained by (8) as follows:
Selection: then target should be compared with . Since , vector should be selected to the next generation as follows:
3.2. The Proposed Hybrid DE (HDE)
The typical DE is simple and easy to be implemented. However, it is likely to be premature too early. Onetoone competing is one of the main reasons. Therefore, improvements including dynamic parameter adjusting, different mutation and crossover strategies, or hybrid algorithms are necessary to be adopted.
Similar to DE, a genetic algorithm (GA) contains crossover, mutation, and selection operations. The crossover operation of GA is quite complicated and its complexity may grow rapidly when the problem scale becomes larger. Fortunately, GA has several efficient selection operations such as roulette wheel selection, tournament selection, and truncation selection. In this study, an HDE that combines the advantages of DE and GA is proposed. The proposed HDE can simplify the evolutionary process, and it can overcome the limitation of onetoone selection of DE and thus prevent premature convergence.
Actually, several scholars also proposed hybrid DEs based on DE and GA (Hrstka and Kucerová [35]; He et al. [36]; Lin [37]), but their mixing modes are quite different from ours. In the proposed HDE, the mutation and crossover operations are the same as in DE while the selection operation is from truncation selection of GA. That is to say, it will be reserved instead of comparing with the target vector when a trial vector is generated. When all trail and target vectors are determined, top NP vectors with better performance are selected to the next generation. The HDEbased procedure is shown in Figure 1.
4. Two Methods for Solving MSJRD
4.1. Linear Programming (LP) Approach for the MSJRD
This method is to summarize the weighted targets and thus converts the multiobjective model to a single one. Take into consideration that two targets have different measurements; it is necessary to standardize two targets beforehand.
4.1.1. Model Analysis Using Linear Programming Approach
Suppose that the weights of two objectives are and , and then the multiobjective problem can be described as
,??, , and are the tolerant maximum total cost, minimum total cost, maximum stockout quantity, and minimum stockout quantity, respectively, which can be given in advance by decision makers (Wee et al. [18]).
Set , . Then, the objective function is changed to where and
Let , that is, .
Note that for standard normal distribution, , so that
That is to say, when and is known, the optimal value of must satisfy (17).
Taking the second derivation of with respect to , we obtain , which means that the optimal is derived from (17) and .
Substituting into (16), the optimal value of for given and is
Finally, the linear programming model can be written as
The goal is to determine the best and to maximize for the given and .
4.1.2. HDEBased Procedures for MSJRD Using LP Approach
When are determined, the optimal delivery cost can be calculated by solving a TSP. can be calculated by when and are known. Change and with the following steps until maximum is obtained.
Step 1. Initialization: set related parameters (CR, F, and NP) for HDE. Set the lower bound () and the upper bound () of respectively. Note that are integers so, ????is obviously 1. According to experience of [2, 10, 15], ??is set sufficiently large to guarantee that the optimal solution does not escape. In this study, it can be set to 100. is randomly generated in the range of 0 and 1. Combining and we get the th individual . Create initial population randomly.
Step 2. For a given , calculate the objective function. When is determined, can be derived accordingly. ??and jointly determine .
Step 3. Differential operations: while stopping criterion is not met, implement mutation and crossover for each individual. After that, the number of population is two times the original one.
Step 4. Genetic operations: select the individuals according to . Those with larger will be chosen to the next generation. Then the number of population is the same as the original one.
Step 5. When the number of iteration reaches a predefined maximum number, output the optimal results; otherwise, repeat Steps 2–4.
4.2. Multiobjective Evolution Algorithm (MOEA) Approach for the MSJRD
In this section, a brief introduction of MOP is given. Then, an HDEbased procedure to handle the MSJRD using noninferior and crowding distance is designed.
4.2.1. Some Definitions of MOP
Definition 1 (multiobjective?optimization?problems?(MOP)).
A general MOP consists of decision variables, objective functions, and constrains. In Definition 1, refers to the decision space and are constrains of MOP.
Definition 2 (Pareto optimal solution). The? optimal solution of MOP is often referred to as the Pareto optimal solution. Let vector belong to and suppose that is a subset of . If there does not exist any vector in that is better than , then is called noninferior solution (or Pareto optimal solution) of . Moreover, if vector is the noninferior solution of , then vector is the Pareto optimal of the MOP.
4.2.2. HDEBased Procedures for MSJRD Using MOEA Approach
There exist many difficulties when applying DE to solve an MOP compared with single objective problem. The main challenges for solving MOP are as follows: how to generate offspring and how to keep Pareto solutions uniformly distributed. The classical DE is not suitable for an MOP since many good solutions may be abandoned due to its onetoone competing mechanism. This will also be confirmed by a numerical example.
Therefore, we also use an HDE which uses truncation selection to choose next generation based on front rank and crowding distance adopted by Qian and li [27]. The steps of calculating crowding distance are presented in Algorithm 1.

In this algorithm, the low front rank corresponds to the high quality of a solution. As to the those individuals with the same front rank, the larger crowding distance means better distribution. Therefore, individuals with lower front rank and larger crowding distance are selected to the next generation.
The first target can be divided into an inventory problem and a delivery problem. When all are determined, the optimal delivery cost can be calculated by solving a TSP. In addition, for a stochastic JRP with normal distributed demand, when , , and are known, the stochastic JRP can then be solved. With the same value of , , and in the second target, we can obtain the corresponding value of the second target. Then change , , and with the following steps until the termination condition is satisfied. The steps of HDEbased approach are described as follows.
Step 1. Initialization: set related parameters (CR, F, and NP) for the HDE. Set the lower bound and the upper bound of , respectively; that is, and . is randomly generated in the range of 0 and 3, which can cover 99.7% of the demand. is randomly generated in the range of 0 and 1. Combining , , and we get the th individual . Create initial population randomly.
Step 2. Calculate the objective function, that is, the total cost and the total shortage quantity of all items.
Step 3. Calculate the Pareto front and crowding distance of each individual.
Step 4. Differential operations: while stopping criterion is not met, implement mutation and crossover for each individual. After that, the number of population is two times the original one.
Step 5. Genetic operations: select the individuals according to the front rank and crowding distance. Then the number of population is the same as the original one.
Step 6. When the number of iteration reaches a predefined maximum number, output the optimal results; otherwise, repeat Steps 2–5.
5. Contrastive Example and Results Analysis
5.1. Basic Data of Numerical Example
The data come from Qu et al. [9]. Table 1 describes the items to be replenished and the center warehouse correspondingly. Tables 2 and 3 are the related parameters of items and distances between suppliers and warehouse, respectively.
In the following, two approaches named LP and MOEA are compared. The comparison contains two aspects: the Pareto solutions and some specific solutions obtained by each method. In the meanwhile, three algorithms used in each method are compared with each other. Table 4 reports related parameters of HDE, DE, and GA.
For LPbased approach, we directly set , , , and according to the advice of the decision makers. This approach is also widely used by other scholars (Wee et al. [18]).
5.2. Comparisons for LPBased and MOEABased Solutions
In this section, the above numerical example is handled using LP and MOEA. For LP, the weight of each objective must be assigned firstly. In order to compare with MOEA, the objectives can be converted to single index by setting the total cost and total shortage quantity with the same weights for MOEA when the Pareto solutions are obtained. The best results for LP when are presented in Table 5. As to MOEA, the highest index after converting is shown in Table 6.
Table 5 shows that HDE and DE are better than GA for LP; Table 6 implies that HDE is better than GA and GA is better than DE for MOEA. In order to further verify the conclusion, we obtained for different , is set from 0.1 to 0.9 and the results are reported in Table 7.
Set the total cost and total shortage quantity with the same weights and , respectively, for MOEA. Then, solutions with the highest weighted objective from the obtained Pareto solutions are shown in Table 8.
Table 7 shows that values of each target and change correspondently when varies, which means that different settings of the weight will result in different decisions. Moreover, when the weight of the first objective (total cost) equals 0.1, the solution is the best. When the difference between and becomes smaller, the solutions become worse. Table 8 implies a similar conclusion for HDE, DE, and GA.
From the comparisons for specific solutions, the following conclusion can be easily drawn: (1) HDE is better than DE or GA no matter whether LP or MOEA is adopted: HDE and DE are more suitable than GA when LP is used; HDE and GA are better than DE when MOEA is used. (2) Different weights for objectives will influence the solutions, for the conflicted objectives, and the assigned weights with large ratios (i.e., ) may result in better solutions.
5.3. Nondominated Solution Analysis of MOEA
For the MOP, there have several metrics to evaluate the quality of the nondominated solutions (Robic and Filipic [38]). However, the implementation of most metrics needs a prerequisite; that is, the true Pareto front must be known. In this study, it is impossible to find the true Pareto front because the MSJRD is a practical problem. So we adopt the metric (Spacing, SP) used by EsparciaAlcázar et al. [39, 40] to measure the distribution of solutions on the Pareto front by evaluating the variance of neighboring solutions. The lower value of SP means that better nondominated solution is obtained.
SP measures the relative distances between the members of Pareto front as where is the number of the first nondominated solutions found. The distance is given by where is the fitness of point on objective and is the mean of all . Table 9 shows the mean and variance of SP by 10 runs using MOEA.
Table 9 shows that the mean and variance of SP obtained by HDE is the lowest, and corresponding values obtained by DE are biggest. That is to say, HDE is better than DE and GA for the MOEA method. The conclusion is consistent with Section 5.2. At the same time, it verifies that the conversion using weights for the MSJRD is feasible.
In order to have a better understanding of the solutions’ distribution of the last generation, the entire nondominated fronts found by HDE, DE, and GA are presented from Figure 2 to Figure 4.
Figures 2, 3, and 4 show that HDE and GA are capable of obtaining Pareto solutions, while the effectiveness and distribution of solutions are much worse for DE. The main reason is the onetoone competing of DE. When trial individual dominates target individual, the trial vector (otherwise the target vector) remains to the next generation. This mechanism is different from the selection operation of the HDE, where all the target and trial individuals are kept to be chosen according to the front rank and crowding distance. So, DE is not suitable for solving this MSJRD.
6. Conclusions and Future Research
In this study, a new multiobjective JRD model with stochastic demand is proposed which takes into account the service level while making the replenishment and delivery decisions. Then, two approaches to solve this complex optimization problem are designed using an improved HDE. The main contributions are as follows.(1)Considering the difficulty to estimate the shortage cost in reality, the shortage quantity is utilized as another standard to evaluate the rationality of decisions besides the total cost. To our best knowledge, this is the first time to propose a practical multiobjective stochastic JRD model. (2)The MOEA is adopted to solve the proposed multiobjective JRD model. The results of the numerical example and Pareto solution analysis show the feasibility of the MOEA to handle the proposed MSJRD. It enriches the application field of the MOEA.(3)The comparison of two approaches for the MSJRD verifies that LP and MOEA are suitable for solving this MSJRD problem. Furthermore, results show that the proposed HDE is more effective than DE and GA whatever LP or MOEA method is used. This illustrates that DE is likely to combine with other algorithms so as to provide a more effective way to solve complex problems.
The future research on the multiobjective JRD problem should consider more realistic assumptions such as uncertain costs, freight consolidation, and budget constraint.
Abbreviations
:  Annual demand rate of item 
:  Major ordering cost 
:  Minor ordering cost of item 
:  Annual inventory holding cost of item 
:  Basic cycle time (decision variable) 
:  Multiplier of (decision variable) 
:  Cycle time of item (decision variable) 
:  Maximum inventory level of item during replenishment interval 
:  Safety stock factor of item (decision variable) 
:  Lead time of items 
:  Standard deviation of item 
:  Least common multiple of 
:  Stockout cost per unit for item 
:  Distribution cost per unit distance 
:  Stopover cost at supplier 
:  Shortest path needed for distribution in th basic cycle. 
Acknowledgments
This research is partially supported by the National Natural Science Foundation of China (70801030; 71101060; 71131004; 71371080) and Humanities and Social Sciences Foundation of the Chinese Ministry of Education (no. 11YJC630275).