Joint Optimization of Distribution Network Design and Two-Echelon Inventory Control with Stochastic Demand and CO2 Emission Tax Charges

This study develops an optimization model to integrate facility location and inventory control for a three-level distribution network consisting of a supplier, multiple distribution centers (DCs), and multiple retailers. The integrated model addressed in this study simultaneously determines three types of decisions: (1) facility location (optimal number, location, and size of DCs); (2) allocation (assignment of suppliers to located DCs and retailers to located DCs, and corresponding optimal transport mode choices); and (3) inventory control decisions on order quantities, reorder points, and amount of safety stock at each retailer and opened DC. A mixed-integer programming model is presented, which considers the carbon emission taxes, multiple transport modes, stochastic demand, and replenishment lead time. The goal is to minimize the total cost, which covers the fixed costs of logistics facilities, inventory, transportation, and CO2 emission tax charges. The aforementioned optimal model was solved using commercial software LINGO 11. A numerical example is provided to illustrate the applications of the proposed model. The findings show that carbon emission taxes can significantly affect the supply chain structure, inventory level, and carbon emission reduction levels. The delay rate directly affects the replenishment decision of a retailer.


Introduction
Sustainable supply chain management (SSCM) has become a hot topic for researchers and practitioners who aim to integrate environmental and social aspects with economic considerations [1,2]. Thus, managing and optimizing sustainable supply chains presents multiple challenges for many companies, among which the three important tasks are facility location, inventory control, and transport model decisions with consideration of CO 2 emissions [3][4][5]. Being able to formulate a decision support system, including the aforementioned three decisions, is a major challenge that can result in a tremendous competitive advantage for a company in the market. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Literature Review
SSCM mainly involves three decisions, namely, operational (short-term), tactical (mediumterm), and strategic (long-term) decisions based on the time frame of decision making. One of the long-term strategic decision problems in supply chain management is distribution network design (DND). The DND problem can be classified into single-, two-, and multi-echelon models, depending on graph levels or echelons [6][7][8]. Wang et al. [9] investigated a two-echelon logistics distribution network design optimization model that is solved by a hybrid algorithm embedded with the particle swarm optimization (PSO) and genetic algorithm (GA). Moreover, facility location problem (FLP) is a classical model that solves the DND. For a more detailed review of integrated DND models, the reader is referred to the work of Lemmens et al. [8,10,11].
The location-inventory problem (LIP) is an extension of the classical facility location problem (FLP), which is also NP-hard [12,13]. LIP is basically an integrated discrete location problem that simultaneously determines location, allocation, and inventory decisions, which has received much attention in the past two decades. Chen and Ting [14] described a multiple ant colony system (MACS) to solve the single-source capacitated facility location problem (SSCFLP). Miranda and Garrido [15] addressed a joint location-distribution-inventory model for a three-layered supply chain, which is based on Lagrangian relaxation with consideration of the validity inequalities derived from the finite set of possible combinations of mean demand and variance. Liao et al. [16] presented a multi-objective location-inventory problem (MOLIP) model in vendor-managed inventory using non-dominated sorting genetic algorithm (NSGA-II). Tsao et al. [17] developed a continuous approximation approach to solve an integrated facility location-allocation and inventory management problem. Tancrez et al. [18] proposed a nonlinear continuous formulation that includes transportation, fixed, handling, and holding costs. Diabat et al. [19] considered a closed-loop location-inventory problem that is solved by a two-phase Lagrangian relaxation algorithm.
Compared with the traditional LIPs, the location-inventory model considering the stochastic environment is more practical. Shu and Sun [20] and Snyder et al. [21] addressed a robust location-inventory problem with uncertain demand at retailers. Shahabi et al. [22] studied a location-inventory problem in a three-level supply chain network with uncertainty parameters. Alhaj et al. [23] investigated a carbon-sensitive two-echelon inventory supply chain model with stochastic demand considering environmental impact.
As far as the inventory control decision problem is concerned, many studies have concentrated on the inventory management of DCs. Shen et al. [24] introduced a distribution center location model and proposed a heuristic algorithm based on Lagrangian relaxation. Cabrera et al. [25] solved a novel inventory-location model with a stochastic capacity constraint based on a periodic inventory control policy. Qu et al. [26] investigated a practical and novel location-inventory problem with stochastic demand considering different replenishment policies. Recently, two-echelon inventory control decisions have elicited concern. Teo et al. [27] addressed the two-echelon joint location-inventory model that arises from coordination of replenishment activities between warehouses and retailers in which the inventory control at both the DCs and retailers was considered. Uster et al. [28] proposed a continuous locationinventory problem with a power-of-two replenishment policy. Shu et al. [29] proposed a twostage stochastic model to address the two-echelon location-inventory network design problem with uncertain demand in which the first stage decides the opening of warehouses and the second determines the warehouse-retailer assignments and two-echelon inventory replenishment strategies.
In relation to the inventory control decisions under multi-echelon distribution networks, previous studies mainly focused on replenishment time, order quantity, and safety stock at DCs or retailers, or both. However, limited studies have been conducted on the interaction between the inventory control decisions of DCs and retailers. The present study considers the inventory control decisions of retailers, which are usually relegated to the traditional LIP model.
The problem addressed in this study is also related to a class of the location-inventory problem with capacity restrictions. Romeijn et al. [30] investigated the two-echelon supply network optimization problem, which considers the capacity and congestion of DCs. Ozsen et al. [31] presented a location-inventory model subject to restriction on the maximum possible inventory accumulation at DCs. Miranda and Garrido [32] considered a capacitated location-inventory model with finite capacities of DCs. The aforementioned papers typically considered the capacity restrictions of candidate DCs but ignored the optimal size choices of candidate DCs. The traditional LIP model mainly focuses on total costs or operator efficiency and only slightly considers transport mode choices and the corresponding external environmental costs. Kutanoglu and Mahajan [33] proposed an integrated optimization model to allocate stock levels at warehouses in a service part logistics network and considered two-echelon variable transportation costs, i.e., the costs of transporting the part from the central to the local warehouses and from the local warehouses to the customers. Luathep et al. [22] considered a three-level location-inventory problem in which the transportation cost per unit shipment, both between plants and warehouses and between warehouses and retailers, is assumed to be proportional to the Euclidean distance. Archetti et al. [34] addressed a closed-loop network with capacity limits, single-item management, and uncertainty on product demands and returns by a fuzzy programming method. Diabat et al. [35] addressed a joint location-inventory problem and extended it to account for the reduction of carbon emissions in which a three-level supply chain network is considered, i.e., one plant, multiple DCs, and multiple retailers. Wu et al. [36] proposed a multi-period location model with transportation economies of scale that distributes a single perishable product, as well as considerations of CO 2 emissions on transport modes.
To the best of our knowledge, existing related studies integrating logistics network design and two-echelon inventory with the considerations of environment impacts are still scarce. The present work aims to fill this gap by focusing on the joint optimization of distribution network design (DND) and two-echelon inventory control (2-EIC) with stochastic demand and CO 2 tax charges. First, an optimal programming model is proposed to address the simultaneous optimization problem of the selection of DCs, transport model choice and supply chain inventory control considering the CO 2 emission taxes charges. Second, some management insights are obtained, which are based on the analysis of numerical experiments.

Problem description
The proposed problem is described as follows. A three-level distribution network consisting of multiple suppliers, DCs, and retailers is considered, as shown in Fig 1 The multiple suppliers serve a set of DCs, which in turn fulfill the demands at the locations of multiple retailers. Inventories are retained at the opened DCs and retailers. The (R, Q) inventory policy is applied to DCs and retailers. Under the (R, Q) inventory policy, an order of Q units is triggered when the inventory level falls below the reorder point R; thus, the resulting inventory position after ordering is in the interval [R, R + Q]. The replenishment lead times from suppliers to DCs are assumed constant. However, the replenishment lead times from DCs to retailers are stochastic, considering the effects of supply shortages at DCs.
In this study, the integrated inventory and location model based on the principle of systemwide cost minimization considering the candidate DC's capacity is investigated.
The integrated optimization model in this paper simultaneously determines three types of decisions: (1) facility location, i.e., the number, location, and size of DCs; (2) allocation, i.e., assignment of suppliers to located DCs and retailers to located DCs, as well as the corresponding optimal transport mode selection; and (3) inventory control decisions on the order quantities, reorder points, and amount of safety stock at each retailer and opened DC. The goal is to minimize the total costs that comprise the facility location, transportation, inventory, and carbon emission tax charges.

Assumptions
To facilitate the presentation of essential ideas without loss of generality, the following basic assumptions are made:

A1
The demand of each customer is stochastic and independent, which follows a normal distribution with known mean and variance.

A2
The replenishment lead time of DCs is constant. The replenishment time of retailers is uncertain.
A3 A single sourcing strategy is adopted in which each DC is assigned to a single supplier, and each retailer is assigned to a single DC.
A4 DCs and retailers adopt a continuous inventory policy (R, Q).

A5
Shortages are allowed and result in unfulfilled demand cost.
A6 Transportation cost includes two parts. The variable transportation costs are linearly proportional to the product of the quantity of products and traveling distances. The fixed transportation cost is also considered, which is used to identify the fixed operator cost per shift of different transport tools.
A7 At each DC and each retailer, CO 2 emissions are linearly proportional to the quantity of products processed. In transportation, carbon emissions are linearly proportional to the product of the quantity of products and traveling distances.

A8
The replenishment lead time of retailers includes two parts, i.e., transport and delay times caused by absence of stock at DCs. We assume that the delay time is proportional to the amount of shortages that occurred at DCs. Moreover, the replenishment lead time of retailers is mainly affected by the delay time caused by out-of-stock products at DCs. Thus, in this study, the effect of transport time on the replenishment lead time of retailers can be ignored.

Cost analysis
Fixed cost for opening and operating DCs. The fixed cost for opening and operating DCs per unit time is calculated by Eq (1).
Replenishment costs of DCs. The replenishment costs considered in this study comprise ordering, inventory holding, and shortage costs.
(1) Ordering cost The expected ordering cost of DC i per unit time is expressed as follows: (2) Inventory holding cost Given that the demand of each DC is a linear superposition function of customer demand, which is normally distributed and independent, the demand rate of DC i is x ilk u k , and the variance of demand per unit time at DC i is ðs 0 This study assumes that DCs adopt a continuous inventory policy (R, Q). In many traditional studies in the existing literature, the order quantity is obtained by the first-order optimality conditions of the replenishment cost function. Therefore, the order quantity is no longer variable and can be eliminated from the model. By contrast, this procedure would be invalid if any constraint relies on the order quantity. In the model of the present study, a constraint dependent on the order quantity was included. Thus, an alternative procedure to solve the resulting model was presented.
The capacity restriction of DCs is not measured by the mean of each DC demand. The constraint declares a maximum inventory level for each DC based on chance-constrained programming. Given a fixed number δ, the probability of the inventory level of DC at peak times exceeds its capacity as shown in Fig 2 This probability is known as the level of service for the inventory capacity constraint at each DC. Consequently, this inventory level corresponds to the reorder point minus the stochastic demand during the replenishment lead time, plus the replenishment order quantity. Thus, the inventory capacity constraint can be written as follows: The reorder point of DC i with the level of service linked to the safety stock level κ i is (3) can be reformulated as follows: Thus, this constraint can be expressed as follows (details of the derivation process of Eq (5) are provided in Appendix A1 in S1 File): where Z 1-δ is the value of the standard normal distribution, which accumulates a probability of 1δ.
The average stock quantity of DC i can be expressed as follows: where k i is called a ''service level", which means that stockouts occur during the lead time with a probability of k i or less, and z k i is the standard normal deviation such that PðZ z k i Þ ¼ k i [38]. For example, if k i = 0.025, which means the lead time demands are satisfied with a probability of 97.5%, z k i ¼ 1:96. Therefore, the expected inventory holding cost per unit time at DC i can be written as follows: (3) Shortage cost Under a continuous review inventory policy (R, Q), the stock-out phenomenon occurs only when the cumulative demand during the replenishment lead time exceeds the maximum inventory level R at that particular time. Thus, the expected unfulfilled demand at DC i can be computed as follows (details of the derivation process of Eq (8) are provided in Appendix A2 in S1 File): where In addition, the variance of unfulfilled demand at DC i can be computed as follows (details of the derivation process of Eq (9) are provided in Appendix A3 in S1 File): Therefore, the shortage cost during the replenishment lead time at DC i is as follows: The expected replenishment cost per unit time at DC i is formulated as follows: Replenishment costs of retailers. The stochastic demand of each retailer follows a normal demand distribution. The replenishment lead time at retailers is assumed to be the transportation time from DCs to retailers plus the time delays caused by the shortages at DCs. Thus, x ilk LQ i can be obtained. Therefore, the expected replenishment lead time at retailer k can be expressed as follows: The variance of the replenishment lead time at retailer k can be expressed as follows: Thus, the demand at retailer k in each replenishment interval of lt k is a normal distribution with a mean ofũ k ¼ uðlt k Þu k and standard deviation of ;where E(Á) and Var(Á) represent the corresponding mean and variance functions, respectively.
(1) Ordering cost The expected ordering cost per unit time at retailer k is expressed as follows: (2) Inventory holding cost The expected inventory holding cost per unit time at retailer k is expressed as follows: (3) Shortage cost Similar to the derivation process of the expected unfulfilled demand of DCs, the expected unfulfilled demand at retailers can be expressed as follows: where 2ðs k Þ 2 and t k ¼ r k Àũ k s k . Therefore, the shortage cost per unit time at retailer k is calculated as follows: The expected replenishment cost per unit time at retailer k can be expressed as follows: Transportation costs. Transportation costs include the variable and fixed transportation costs. The two-echelon variable transportation costs include two parts, i.e., inbound network transportation cost (cost of shipping the products from the suppliers to the DCs) and outbound network transportation cost (cost of shipping the products from the DCs to retailers). Furthermore, the fixed costs of transportation are independent of the quantity of products and shipping distances. Thus, the transportation cost per unit time can be formulated as follows: Carbon emission costs. The costs of CO 2 emission tax charges into the objectives were incorporated with consideration of the sustainability of the distribution network. Under the carbon emission tax policy, no strict cap is applied to carbon emissions, but the total carbon emissions will be penalized under a carbon tax policy. The tax is a financial penalty, which is assumed to be equal to the product of the total carbon emissions and the carbon tax rate. With regard to the carbon emission sources, two main activities were considered: transportation and inventory management. Furthermore, the average shipment from supplier o to DC i was calculated by mode l using the corresponding expected demand of DC i (i.e., u 0 i ), while the average shipment from DC i to retailer k by mode l is expressed with the expected demand of retailer k (i.e., u k ). Therefore, the expected total carbon emission cost per unit time can be calculated as follows: ev ilk x ilk u k l ik :ð20Þ

Model formulation
Based on the preceding analysis, the formulation of the integrated inventory and location models was presented, which determines the optimal location and size of DCs, as well as the optimal order quantity and reorder point for the opened DCs and retailers. The model is formulated as follows: s.t.
x ilk y lk 8i; l; k; ð32Þ x 0 oil ; x ijk ; u ia 2 f0; 1g 8o; a; i; j; k; ð33Þ The objective function of Eq (21) represents the total distribution network costs. The first term is the fixed cost of DCs. The second represents the fixed order, inventory holding, and penalty costs for unfilled demand at DCs, while the third term is that of retailers. The fourth term includes the total transport costs for inbound and outbound network transportation. The fifth and sixth terms are CO 2 emission tax charging costs of inventory and transport process, respectively.
Constraints (22) and (23) represent the mean and variance of each DC demand. Constraint (24) ensures that each retailer is served by exactly one distribution center by only one transport mode. Constraint (25) implies that the candidate DCs select an option size from the alternative size set. Constraints (26) and (27) ensure that if two points on one echelon are related to each other, then one mode of transportation must be shared between them. Constraint (28) ensures that the order quantity of retailers is less than the corresponding upper bound. Constraint (29) states the capacity restriction for each DC. Constraint (30) implies that if DC i is opened, it must be serviced by a certain supplier. Constraint (31) ensures that only when DC i is opened can it serve retailers using various transport modes. Constraint (32) represents the transport mode selection constraints between inbound and outbound networks. Constraint (33) enforces that location decision and transport mode selection decision variables are binary. Finally, constraint (34) specifies that other decision variables are non-negative.

Numerical Experiments
The proposed model and solution algorithm are applied to a simple three-echelon distribution network, which comprises a single supplier, three candidate DCs with three alternative sizes, and ten retailers. Three alternative transport modes can be chosen between supplier and DCs, as well as between DCs and retailers.

Data input
The corresponding parameters are the following:  Table 1.

Effect of carbon emission tax rate
This section studies the effect of carbon emission tax rate on distribution network design, inventory control decisions, and total cost invested. Six scenarios were considered with different carbon emission taxes to analyze the corresponding results. From scenarios 1 to 6, the carbon emission tax charges are set to 10, 20, 30, 40, 50, and 60 $/ton, respectively.
We first analyze the effects of carbon emission tax rate on location and inventory decisions of the distribution network. The results in Table 2 show that the total CO 2 emissions of distribution network gradually decreases from 154.68 tons to 116.40 tons when the carbon emission tax rate increases from $10/ton to $60/ton. However, the total cost of the entire distribution network increases, thereby implying that the decrease range of carbon emissions is insufficiently large to eliminate the effect of the carbon tax rate increase on the carbon emission cost. Tables 3 and 4 indicate the effects of the carbon emission tax rate on the inventory control parameters of the opened DCs and retailers, respectively. Fig 3 shows that the order point of DCs remains the same and that the order quantity continuously decreases. The results presented in Table 3 indicate that the replenishment lead time and order point of retailers remain the same. Moreover, the order quantity in DCs continuously decreases from 127.82 tons to 92.61 tons with the increase of the carbon tax rate from 10$/ton to 60$/ton. Moreover, the

Parameters Values
Mean daily demand of retailer k (ton) u k = U [5,20] Daily demand standard deviation of retailer k (ton) σ k = U [1,5] Three Ordering cost at DC i ($) OR i = U [3,200,3,800] Ordering cost at retailer k ($) or k = U[1,000, 1,400] Holding cost per ton per time unit at DC i ($) H i = U [18,20] Holding cost per ton per time unit at retailer k ($) h k = U [8,13] Unit penalty of shortage per time unit at DC i($) P i = U [8,9] Unit penalty of shortage per time unit at retailer k ($) p k = U [8,10] Carbon emission per unit shipment to deal with within DC (ton/ton) order quantity of retailers continuously decreases with the increase in tax rate as shown in Table 4. Decreasing the average stock level is reasonable to decrease the total carbon emission of the entire distribution network, especially with the increase of the carbon tax rates.
The effect of various carbon emission rates on the distribution network design was also analyzed. The results presented in Fig 4. confirm that all opened DCs are consistently served by transport mode 1 and that all retailers are consistently served by transport mode 4. The number of opened DCs is equal to 1; it is stable, i.e., its size remains the same. However, the location of opened DCs varies with the increase in the carbon emission tax rate. For example, when the carbon emission tax rate is $10/ton, DC2 is selected to open, whereas for the carbon emission tax rate of $20/ton, the DC3 is selected.

Effect of delay rate
The impacts of the delay rate parameter λ on the supply chain distribution network configuration and corresponding inventory control decisions are examined. Six different scenarios, i.e., λ = 0.04, 0.08, 0.12, 0.16, 0.20, and 0.24, are considered. The remaining parameters are set to be the same as those of the previous test. The analysis of the effects delay rate parameter on system performance of distribution networks is given in Table 5, clearly indicating that when the delay rate increases from 0.04 to 0.24, the replenishment cost of retailers gradually increases from $8,146.09 to $8,287.82. The inventory control parameters of the opened DCs and retailers corresponding to different delay rates are shown in Tables 6 and 7, respectively. According to Table 6 and Fig 5, the order quantity and reorder point of opened DCs remain the same during the increase of the delay rate. For example, the order quantity and reorder point at DC2 are 109.52 tons and 99.00 tons, respectively, when the delay rate is varied. Table 7 indicates that the replenishment lead time, reorder point, and order quantity of retailers continuously increase with the growth of delay rate. Based on Assumption A8, the delay rate of the replenishment lead time in retailers show significant impacts on inventory decisions at DCs and retailers. When the value of the delay rate parameter λ increases from 0.04 to 0.24, the replenishment lead time of retailers lengthens because of the shortage of DCs. This condition implies that the replenishment lead time increases with the increase of the delay rate. Consequently, the retailers must set a higher reorder point and larger order quantity to reduce the risk of shortage.
In this study, six different delay rate scenarios were considered as shown in Fig 6 Obviously, the number of opened DCs is stable and equal to 1 when the delay rate increases from 0.04 to 0.24. Moreover, the sizes of opened DCs are always kept the same. The transport mode between the DCs and retailers remains the same, whereas the transport mode between the supplier and DCs changes when the delay rate increases from 0.09 to 0.24. The reason may be that the available discrete size series of DCs have different costs; the best choice of the size is a result of the balance between the transportation and inventory costs.

Discussion and analysis
The aforementioned results indicate that an increase in the CO 2 emission tax rate leads not only to a reduction in the total distribution network CO 2 emissions but also to an increase in the total distribution network costs. The benefits gained from CO 2 emission decrease by a charging policy that does not compensate for the extra costs of CO 2 emission charges. Therefore, CO 2 charging policy may contribute to a reduction of the CO 2 emissions especially when CO 2 is charged at a high rate. Moreover, the order quantities of opened DCs and retailers decrease with the increase of the CO 2 emission tax rate. The reason is that CO 2 emissions of inventory management are related to the average stock quantity of retailers and DCs. The preceding analysis implies that the policy of CO 2 emission taxes will decrease the total CO 2  Joint Optimization of DND and 2-EIC with Stochastic Demand and CO 2 Emission Tax Charges emission of distribution network and induce supply chain managers to choose more green transport modes and corresponding inventory control policies. The delay parameter λ at DCs shows a significant effect on the corresponding inventory control parameters of retailers such as replenishment lead time, reorder point, and order quantity of retailers. Specifically, the replenishment lead time, reorder point, and order quantity of  retailers will increase with the increase of the delay rate at DCs, while the corresponding inventory control parameters of opened DCs are kept the same. An increase in the replenishment lead time of retailers will lead to higher-order quantity and reorder point. Moreover, the location decision of DCs is insensitive to the changes of the delay rate. Specifically, the locations and sizes of opened DCs remain the same with the delay rate increase while the location of opened DCs varies with the increase of the CO 2 emission tax rate.

Conclusions
A new mathematical model to design a three-level distribution network design problem that considers the capacity constraints of DCs, transportation mode choices, and carbon emission tax policy was studied. Location, transportation, and inventory issues are included. The objective is to determine the distribution network decisions, such as the number of opened DCs, their locations, their sizes, and the assignments of each retailer and each opened DC, as well as the transportation mode choices. Furthermore, the inventory control decisions on the order quantity, reorder point, and service level at each retailer and opened DC are also determined. This study differs from previous related literature because the interaction effects of the twoechelon inventory control decisions have been analyzed. These effects include the shortage at DCs that possibly prolongs the replenishment lead time of retailers and the stock level of DCs being influenced by the demands of assigned retailers. In addition, the order quantity of products and reorder point for opened DCs and retailers were simultaneously determined.
An integrated inventory-location model was formulated as a mixed integer nonlinear programming problem with chance constraints corresponding to the capacity constraints of DCs and nonlinear objective function. LINGO software was used to solve the resulting MINLP problem. The effects of carbon emission tax rate and delay rate on the distribution network design, inventory control decisions, and total costs were discussed.
Future research directions are as follows: (1) investigating inventory sharing and lateral trans-shipments, which are becoming increasingly common in third-party logistics systems; (2) considering the extension of allowing the retailers to be sourced from multiple DCs; and (3) solving the presented model through heuristic or meta-heuristic algorithms.