Bo Dai

Haoxun Chen
email: haoxun.chen@utt.fr

Yuan Li
email: yuan.lya@alibaba-inc.com

Yidong Zhang

Xiaoqing Wang

Yuming Deng
email: yuming.dym@alibaba-inc.com

An Alternating Direction Method of Multipliers for Optimizing (s, S) Policies in a Distribution System with Joint Replenishment Volume Constraints

Keywords: Inventory management, distribution system, joint replenishment, (s, S) policy, scenario optimization, e-commerce

In this paper, we study a two-echelon distribution system in which multiple products are jointly replenished at each stocking location and the inventory of each product at each location is controlled by an (s, S) policy. The transportation capacity in volume of products for each joint replenishment is limited, and linear rationing policies are used for both on-hand inventory and transportation capacity allocation in the system in case of lack. We propose a novel scenario-based model for the optimization of the (s, S) policies in the system that considers the rationing policies. Because of its high complexity when the number of scenarios is large, an Alternating Direction Method of Multipliers is proposed to solve the model. Based on real data, forty instances were generated and tested to evaluate the model and the solution method. Our numerical experiments show that for these instances this method could find a better solution in a much shorter computation time compared with CPLEX 12.9, whereas the latter often runs out of memory for large-size instances on a personal computer. Moreover, the inventory policies found by this scenario-based optimization approach can reduce costs by 5.1% and improve fill rates by 9.7% on average compared with those currently used in Alibaba.

INTRODUCTION

In the era of e-commerce today, the competition between e-retailers becomes more and more fierce. To increase market shares and gain competitive advantages over their adversaries, e-commerce companies have been making efforts to improve the inventory management in their supply chains. Effective inventory management can help these companies reduce their logistics costs while assuring high service levels to their customers.

In recent years, e-commerce has been developing very fast in China, with the emergence of some e-commerce giants such as Alibaba and JD.com. Under a project of Alibaba Innovative Research (AIR), we investigate the optimization of inventory policies in a distribution system dedicated to on-line Tmall supermarket. This system consists of a central distribution center (CDC) and multiple regional distribution centers called front distribution centers (FDCs) in Alibaba. The CDC is the supplier of all FDCs, which serve final customers directly, while the CDC's suppliers are manufacturers. This system stocks and sells multiple products, which are jointly replenished at each stocking location. The customer demand of each product observed by each FDC in each period (day) is stochastic. An (s, S) policy is adopted for the inventory control of each product at each stocking location in the system. Because of limited transportation capacity in each distribution channel, the volume of all products in each joint replenishment cannot exceed a given maximum volume. Previously, we studied safety stock policies and inventory replenishment planning for a single period in such a system [START_REF] Dai | Inventory replenishment planning in a distribution system with safety stock policy and minimum and maximum joint replenishment quantity constraints[END_REF][START_REF] Dai | Inventory Replenishment Planning of a Distribution System with Storage Capacity Constraints and Multi-Channel Order Fulfilment[END_REF]. In this paper, we focus on the optimization of (s, S) policies in the system. The objective is to minimize the total cost of the distribution system composed of ordering costs, holding costs and lost sales costs, with the maximum joint replenishment volume constraints.

To our best knowledge, the problem considered in this paper was never studied in the literature except for a conference paper [START_REF] Dai | Optimization of (s, S) inventory policies in a distribution system with minimum and maximum joint replenishment quantity constraints[END_REF]. However, the optimization model proposed in that paper does not consider stock allocation/rationing at the CDC and no algorithm was proposed for solving the model except for commercial solver CPLEX.

Because of its high complexity, this problem is unlikely to be solved analytically or by a commercial solver like CPLEX in an acceptable computation time. Because of this, we propose a novel scenario and decomposition-based approach to solve it. In contrast to the work of [START_REF] Noordhoek | A simulation-optimization approach for a service-constrained multi-echelon distribution network[END_REF], we consider a distribution system with joint inventory replenishment subject to transportation volume constraints rather than with inventory replenishment of a single product. Moreover, our optimization approach is mathematical programming model-based rather than simulation-based. Although scenarios in our model correspond to realizations of stochastic demands in simulation, our approach has an advantage over theirs by allowing the evaluation of the quality of its solution based on the relative gap between an upper bound and a lower bound obtained by the corresponding mathematical programming model.

For this scenario optimization, the main challenge for establishing a deterministic equivalence model for our stochastic inventory policy optimization problem is to formulate the (s, S) policy used in each stocking location and the stock allocation policy/rationing rule used at the CDC. We first build a novel scenario optimization model for this system, which takes account of the rationing rule and all side constraints. Because solving the model optimally by a commercial solver like CPLEX is quite time consuming for instances with a large number of scenarios, we propose an Alternating Direction Method of Multipliers (ADMM), which is a decomposition approach, to solve the model near-optimally and quickly. Our numerical experiments on four sets of instances generated based on real data show that for these instances this method could find better solutions in a much shorter computation time compared with CPLEX 12.9 and could avoid the out-of-memory problem of CPLEX. Moreover, the inventory policies found by our method can reduce costs and improve fill rates compared with those currently used in Alibaba.

The rest of this paper is arranged as follows. The literature related to our study is reviewed in Section 2. Section 3 introduces the distribution system considered and the inventory optimization problem we study. A scenario-based optimization model is formulated for the problem in Section 4. Section 5 describes a decomposition-based ADMM approach for solving this model. The performance of this approach is evaluated by numerical experiments and its inventory policies found are compared with those currently used in Alibaba in Section 6. Section 7 summarizes the results of this paper and suggests directions for further study.

LITERATURE REVIEW

In this section, we review the literature related to our work. Inventory management gained many attentions in the past decades, but most of them dealt with a single stocking location [START_REF] Veinott | On the optimality of (s, S) inventory policies: new conditions and a new proof[END_REF][START_REF] Lin | Single-item repairable inventory system with stochastic new and warranty demands[END_REF][START_REF] Teunter | Inventory control with demand substitution: new insights from a two-product Economic Order Quantity analysis[END_REF][START_REF] Qiu | A robust optimization approach for multi-product inventory management in a dual-channel warehouse under demand uncertainties[END_REF].

Previous studies of inventory management of distribution systems were focused on so called 'one warehouse multi-retailer system' with a single product [START_REF] Federgruen | Allocation policies and cost approximations for multilocation inventory systems[END_REF]Axsater, 1990;Axsater, 2003a;[START_REF] Axsater | Supply chain operations: serial and distribution systems, Supply Chain Management: Design, Coordination and Operation[END_REF][START_REF] Wang | Optimising inventory placement in a two-echelon distribution system with fulfillment-time-dependent demand[END_REF][START_REF] Berling | Controlling inventories in omni/multi-channel distribution systems with variable customer order-sizes[END_REF]. Under the balance assumption, [START_REF] Diks | Optimal control of a divergent multi-echelon inventory system[END_REF] proved it is optimal for a distribution system to control each stocking location by a base-stock policy when ordering costs at each stocking location are negligible. However, the balance assumption may not hold and ordering costs may exist in the system. In either of the two cases, the distribution system's optimal inventory policy is not known. Because of this, most studies assume a given type of inventory policy for each stocking location in a distribution system and try to optimize the parameters of these policies. One challenge for this optimization is that, except for determining the optimal inventory policy parameters, a stock allocation (rationing) decision must be made for each stocking location when its on-hand inventory cannot completely meet all replenishment orders of its immediate downstream locations. Various stock allocation policies/rationing rules were proposed for distribution systems [START_REF] Eppen | Centralized ordering policies in a multi-warehouse system with lead times and random demand[END_REF][START_REF] Van Der Heijden | Stock allocation in general multi-echelon distribution systems with (R, S) order-up-to-policies[END_REF][START_REF] Bollapragada | Centralized ordering and allocation policies in a two-echelon system with non-identical warehouses[END_REF]. All these rules belong to the class of linear rationing rules [START_REF] Lagodimos | Service performance of two-echelon supply chains under linear rationing[END_REF]. Each linear rationing rule allocates the shortage of an upstream stocking location to its immediate downstream stocking locations proportionally according to their rationing fractions. One application of rationing rules in a distribution system can be found in [START_REF] Jula | A supply-chain optimization model of the allocation of containerized imports from Asia to the United States[END_REF]. [START_REF] Chen | Stock Allocation in a Two-Echelon Distribution System Controlled by (s, S) Policies[END_REF] reviewed several stock allocation polices/rationing rules studied in the literature and proposed two nonlinear rationing rules for distribution systems controlled by (s, S) policies. They showed that the nonlinear rationing rules could reduce the total cost of a distribution system by 1.5% on average. Although a nonlinear rationing rule may outperform a linear one, the former is usually much more complicated than the latter. As a result, linear rationing rules are usually used in practice because of their simplicity and near-optimality.

Joint inventory replenishment is an industrial practice. It has been proved that joint replenishment can reduce ordering costs when multiple products are involved in an inventory system [START_REF] Goyal | Joint replenishment inventory control: Deterministic and stochastic models[END_REF][START_REF] Abouee-Mehrizi | Optimal joint replenishment and transshipment policies in a multi-period inventory system with lost sales[END_REF]. Readers can refer to a recent review paper for the studies on joint replenishment [START_REF] Bastos | A systematic literature review on the joint replenishment problem solutions: 2006-2015[END_REF]. In most joint replenishment models, ordering costs of multiple products consist of product-independent major ordering costs and product-dependent minor ordering costs [START_REF] Khouja | A review of the joint replenishment problem literature: 1989-2005[END_REF]. Various joint replenishment policies were proposed in the literature, such as can-order policy, periodic (s, S) policy, (Q, S), and (T, S) policy. The inventory policies considered in this paper can be considered periodic (s, S) policies [START_REF] Viswanathan | Periodic Review (s, S) Policies for Joint Replenishment Inventory Systems[END_REF] where the review period (or review interval) of each stocking location is given. Most studies on joint replenishment only consider a single stocking location with only few exceptions. [START_REF] Axsäter | A joint replenishment policy for multi-echelon inventory control[END_REF] studied a continuous-review two-echelon distribution system composed of a central warehouse and multiple identical retailers, where the retailer who has the lowest inventory position orders a batch when the aggregate inventory position of all retailers in the system reaches a joint reorder point. [START_REF] Zhou | A multi-product multi-echelon inventory control model with joint replenishment strategy[END_REF] considered a multi-product multi-echelon inventory system composed of a supply network and a distribution network linked through a producer, where the inventory of each stocking location is controlled by a (T, S) policy. A mathematical programming model is established for optimizing the order cycles and the order up to levels of all products at all stocking locations. However, this model is only an approximate model, because the holding cost and the shortage cost of each stocking location are only approximately formulated. Moreover, the model is only approximately solved by using a genetic algorithm, so our modeling and solution method are both different from theirs. Çapar (2013) studied joint shipment consolidation in a two echelon distribution system, where the distribution center uses a time-based shipment consolidation policy to dispatch to retailers their accumulated orders at the end of each consolidation cycle. An exact method was developed to determine the optimal replenishment quantity for the distribution center and the base-stock level for each retailer. However, only one product is involved in the system and the joint shipment consolidation of a single product is different from the joint replenishment of multiple products we study. [START_REF] Carvajal | Heuristic approaches for a two-echelon constrained joint replenishment and delivery problem[END_REF] studied joint inventory replenishment in a two-echelon distribution system under constant final demands, budgetary and storage capacity constraints. They proposed a hybrid meta-heuristic algorithm and a memetic algorithm for optimizing the system. However, they considered a distribution system with constant demands, whereas we consider a distribution system with stochastic demands.

Only few papers studied the optimization of an inventory system with multiple stocking locations and controlled by (s, S) policies [START_REF] De Kok | A typology and literature review on stochastic multi-echelon inventory models[END_REF]. [START_REF] Kukreja | A model for lumpy demand parts in a multi-location inventory system with transhipments[END_REF] considered a single echelon continuous review inventory system with multiple stocking locations controlled by (s, S) policies, where lateral transshipments among stocking locations are allowed. [START_REF] Fattahi | Investigating replenishment policies for centralised and decentralised supply chains using stochastic programming approach[END_REF] investigated a serial system with a manufacturer and a retailer controlled by (s, S) policies. They considered maximum production capacity in the manufacturer and maximum ordering capacity in the retailer. The systems studied in these papers are simpler than distribution systems and neither of them considers joint inventory replenishment. For distribution systems, the optimization of (R, Q) policies or (R, nQ) policies was also studied in the literature (Axsäter, 2003a[START_REF] Al-Rifai | An efficient heuristic optimization algorithm for a two-echelon (R, Q) inventory system[END_REF][START_REF] Berling | Multi-echelon inventory control: an adjusted normal demand model for implementation in practice[END_REF]. All these papers only consider one product without joint replenishment and assume identical review period (the time between two consecutive reviews) for all stocking locations. Moreover, they did not explicitly consider the stock allocation policy/rationing rule used in each upstream stocking location of a distribution system. Three kinds of methods have been used to optimize inventory polices of a stochastic inventory system in the literature, i.e., analytical, simulation-based, and scenario-based methods. Among them, most analytical approaches were used for the optimization of base-stock policies [START_REF] Clark | Optimal policies for a multi-echelon inventory problem[END_REF], or (R, Q) policies in one warehouse multi-retailer systems (Forsberg, 1996, Axsater, 2000). However, using analytical methods to optimize (s, S) policies for a multi-echelon system was rare because of its high complexity. Simulation-based optimization approaches use a metaheuristic as optimization engine and simulation for cost evaluation of an inventory system. [START_REF] Jalali | Simulation optimization in inventory replenishment: a classification[END_REF] provided an overview of simulation-based approaches for inventory optimization. [START_REF] Tsai | A simulation-based multi-objective optimization framework: A case study on inventory management[END_REF] considered multiple objectives and used a simulation-based approach to optimize the reorder points and the order quantities in a two-echelon serial system. The three objectives considered are the expected inventory cost, inventory level, and frequency of shortage. [START_REF] Noordhoek | A simulation-optimization approach for a service-constrained multi-echelon distribution network[END_REF] developed an approach that combines scatter search and simulation for the optimization of (s, S) policies in a multi-echelon distribution system of Dutch food retail industry. Numerical results show this approach outperforms a nested bisection search method. [START_REF] Avci | A multi-objective simulation-based optimization approach for inventory replenishment problem with premium freights in convergent supply chains[END_REF] proposed a simulation-based approach for the optimization of an inventory system in an automotive supply chain. They considered multiple objectives and designed a heuristic algorithm to optimize the objectives simultaneously. Scenario-based methods were also used in inventory optimization. [START_REF] Fattahi | Investigating replenishment policies for centralised and decentralised supply chains using stochastic programming approach[END_REF] formulated the problem of optimizing inventory policies in a single product serial inventory system with a retailer and a manufacturer as a mixed-integer linear program using scenario approximation. The model was solved by applying either an evolutionary or imperialist competitive algorithm. [START_REF] Nguyen | Supplier selection and operation planning in biomass supply chains with supply uncertainty[END_REF] formulated a stochastic programming model for supplier selection in a biomass supply chain with uncertain supply. The deterministic equivalence model of the stochastic program was solved by applying a regularized L-shaped algorithm.

Because only a meta-heuristic is used as its optimization engine, the quality of the solution found by a simulation-based approach for the optimization of a stochastic inventory system cannot be well evaluated. This is one major weakness of such approach. Although each scenario corresponds to a realization of random variables in a stochastic inventory system as in simulation, the quality of the solution found by a scenario-based approach can be well evaluated, because for a given number of scenarios, the stochastic optimization problem can be equivalently transformed into a deterministic mathematical program. This is one advantage of scenario-based approach versus simulation-based approach.

In the literature, both inventory systems with backorders and lost-sales were widely studied. A recent literature review on inventory systems with lost-sales was given by [START_REF] Bijvank | Lost-sales inventory theory: A review[END_REF]. They concluded that the optimization of an inventory system with lost-sales is usually much more difficult than the optimization of its counterpart with backorders. For this reason, most studies on inventory systems with lost-sales only considered one stocking location with one product. Only few papers studied distribution systems with lost sales [START_REF] References Andersson | A two-echelon inventory model with lost sales[END_REF][START_REF] Seifbarghy | Cost evaluation of a two-echelon inventory system with lost sales and approximately Poisson demand[END_REF][START_REF] Hill | A two-echelon inventory model with lost sales[END_REF]. Because lost sales occur in our considered distribution system, we adopt an inventory model with lost sales for this system. In addition, in this system, the inventory review period of the CDC is different from that of each FDC.

With respect to the previous works reviewed above, our work makes the following contributions to the literature.

1. A two-echelon joint replenishment distribution system with lost sales and managed by (s, S) policies with transportation volume constraints and different review periods at the CDC and the FDCs is considered. Most of previous works on distribution systems consider backorders and no previous work studied a distribution system with all the features we consider.

2. By applying scenario approximation, a novel mixed-integer linear programming model is formulated for the optimization of the inventory policies. This model well considers stock allocation/rationing in the system and all side constraints.

3. An alternating direction method of multipliers is proposed to solve the model. This method can find a near-optimal solution for the instances generated based on real data in a reasonable computation time.

PROBLEM DESCRIPTION

A two-echelon distribution system dedicated to on-line Tmall supermarket of Alibaba is studied, which is illustrated in Figure 1. This distribution system consists of a Central Distribution Center (CDC) and multiple regional distribution centers called Front Distribution Centers (FDCs), where the CDC is the supplier of each FDC which serves final customers. Each stocking location in the system stores multiple fast moving products. It is assumed that each product's demand observed by each FDC in each period (day) is normally distributed, and the demands of different products at different stocking locations in different periods are independent. Note that our scenario-based optimization model and solution method are also applicable in case when the demand of each product at each stocking location in each period is subject to any other distribution such as Poisson or compound Poisson distribution.

Figure 1: A two-echelon distribution system for Tmall supermarket

The inventory of each product at each stocking location is periodically reviewed, where each period corresponds to one day. All FDCs have the same review period (the time between two consecutive reviews), whereas the CDC's review period is a multiple of each FDC's. Suppose that the review period of the CDC and each FDC are R and 1, respectively. To distinguish between the two review periods, the CDC's review period is also called review cycle hereafter. Each stocking location's replenishment lead time is given, and the CDC may have a lead time different from that of the FDCs. Each stocking location jointly replenishes multiple products involved. Because of limited transportation capacity for each distribution channel, the volume of each joint replenishment of a stocking location cannot exceed a given maximum volume. The inventory of each product at each stocking location in this distribution system is controlled by an (s, S) policy, where s and S are reorder point and order-up-to level, respectively. In addition, when the on-hand inventory of the CDC cannot satisfy all replenishment orders of the FDCs, a stock allocation/rationing rule is used to allocate the on-hand inventory among the orders.

We want to optimize the (s, S) policies of the distribution system in a planning horizon starting from period r and ending with period T, where 1 ≤ ≤ r R and T = 1 + ×r N R for some positive integer N. The number of periods considered in the planning horizon is × N R . In this planning horizon, each FDC reviews its inventory in each period, and the periods in which the CDC reviews its inventory can be written as 1 × + n R , for integer n with 0 1 ≤ ≤n N and 1 × + n R ≥ r. Note that considering a finite horizon in the optimization of a dynamic system is also adopted by lot-sizing models for production planning [START_REF] Balkhi | On a finite horizon production lot size inventory model for deteriorating items: An optimal solution[END_REF][START_REF] Zhou | A finite horizon lot-sizing problem with time-varying deterministic demand and waiting-time-dependent partial backlogging[END_REF].

Holding costs, lost sales costs, and two types of ordering costs are considered in the system. Minor ordering cost of each stocking location is determined by transportation costs of its inventory replenishments. Major ordering cost of the CDC is charged every time when it places a joint replenishment order to external suppliers, whereas that cost is negligible for each FDC. The negligibility of major ordering cost for each FDC is because this cost is much smaller than transportation costs for its order delivery. Our objective of optimizing the (s, S) policies of the distribution system is to minimize its expected total cost in the planning horizon.

The problem studied is a very complex stochastic optimization problem since it deals with a distribution system controlled by (s, S) policies with stochastic demands of multiple products jointly replenished at each stocking location with maximum joint replenishment volume constraints. Because of the complexity, it is unlikely to find an analytical approach for this optimization. For this reason, we propose a scenario-based approach for this optimization, where the random demands are approximated by a finite number of their possible realizations (scenarios) with respective occurrence probabilities. For a given number of scenarios, this stochastic optimization problem can be transformed into its deterministic equivalence model. In this equivalence model, the expected total cost of the system is replaced by its total cost averaged over all demand scenarios.

MATHEMATICAL MODELLING OF THE PROBLEM

In this section, a mixed integer linear program model is established for optimizing (s, S) policies in the distribution system considered. Since the demands of the distribution system are relatively stationary, we assume s and S for each product at each stocking does not change over time, but the optimization model and the solution method proposed in this paper can be easily extended to ones with time-dependent s and S.

Since most stock allocation/rationing rules proposed in the literature and that used in the distribution system of Alibaba are linear ones, we restrict rationing rules used in this system to linear ones. That is, a linear rationing rule is applied to both the allocation of on-hand inventory of the CDC among orders of the FDCs for each product and the allocation of transportation capacity of each distribution channel among multiple products for each stocking location in the system. This restriction will not have a significant impact on the cost of the system [START_REF] Chen | Stock Allocation in a Two-Echelon Distribution System Controlled by (s, S) Policies[END_REF] while making its stock allocation decisions simpler. Like the on-hand inventory allocation, when a linear rationing rule is applied to the transportation capacity allocation, the shortage of transportation capacity is allocated proportionally among the products involved in each joint replenishment according to their rationing fractions pre-specified. We assume all rationing fractions (coefficient) in each linear rationing rule used are given.

Before presenting the model, we first introduce the following notations used in it.

Indices

i: index of a product, i ∈ Ι = {1, 2, …, I}, where I is the number of products considered. j: index of a stocking location, j ∈ ϑ = {0, 1, 2, …, J}, where j = 0 represents the CDC and j > 0 represents a FDC, and J is the number of FDCs in the distribution system considered.

t, τ: index of a period, t,τ ∈ Τ = {r, r+1, …, 1 + ×r N R }, where r with 1 ≤ r ≤ R is the starting period of the planning horizon, R is the inventory review period of the CDC, N is a positive integer, and × N R is the number of periods considered in the planning horizon.

k: index of a scenario index, k ∈ Κ = {1, 2, …, K}, where K is the number of demand scenarios considered.

Parameters

= = ∑ I ij i f , ' 0 ij f ≥ k p : occurrence probability of scenario k Decision Variables , t k ij
IP : product i's inventory position at stocking location j in period t prior to receiving shipments from its supplier in this period under scenario k , t k ij I : product i's on-hand inventory at stocking location j in period t prior to receiving shipments from its supplier in this period under scenario k , t k ij LS : lost sales of product i at stocking location j in period t under scenario k ij s : the reorder point of product i at stocking location j ij S : the order-up-to level of product i at stocking location j , t k ij x : product i's real replenishment quantity for stocking location j in period t under scenario k after both on-hand inventory allocation/rationing of the CDC for this product and the transportation capacity allocation/rationing among all products involved in the corresponding joint replenishment , t k ij q : product i's replenishment quantity for FDC j in period t under scenario k, determined by its (s, S) policy , ' t k ij q : product i's replenishment quantity for FDC j in period t under scenario k after on-hand inventory allocation/rationing of the CDC for this product but before the transportation capacity allocation/rationing among all products involved in the corresponding joint replenishment Model SBO:

( ) , , , , 0 0 \{0} ' ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ =     ⋅ ⋅ + ⋅ ⋅ + ⋅ + ⋅           ∑ ∑∑ ∑ ∑ ∑∑ ∑ SBO t k t k t k t k k ij ij ij i ij ij ij k t i j t t i j Z Min p h I t v x O y l LS K T I J T T I J (1)
subject to:

( ), ( 1), , ,
, , , \ {0}, ,

- + = + - + ∈ ∈ ∈ ∈ ij t L k t k t k t k t k ij ij ij ij ij i I I x d LS j t k I J T K (2) 0 J ( ), ( 1), , , 0 0 0 1 , , , - + = = + - ∈ ∈ ∈ ∑ i t L k t k t k t k i i i ij j i I I x x t k I T K (3) 1 , , , , , \{0}, , τ τ -
= - = + ∈ ∈ ∈ ∈ ∑ ij t t k t k k ij ij ij t L i IP I x j t k I J T K (4) 0 1 , , , 0 0 0 , , , τ τ -
= - = + ∈ ∈ ∈ ∑ i t t k t k k i i i t L i IP I x t k I T K (5) If , ≤ t k ij ij IP s , , = 1, , \ {0}, , ∈ ∈ ∈ ∈ t k ij i y j t k I J T K (6) If , > t k ij ij IP s , , = 0, , \ {0}, , ∈ ∈ ∈ ∈ t k ij i y j t k I J T K (7) , , , , \{0}, , ≤ - ∈ ∈ ∈ ∈ t k t k ij ij ij i q S IP j t k I J T K (8) If , 0 0 ≤ t k i i IP s , , 0 = 1, , 1 , , 0 1, ∈ = + × ≥ ≤ ≤ - ∈ t k i i y t n R t r n N k I K (9) If , 0 0 > t k i i IP s , , 0 = 0, , 1 , , 0 1, ∈ = + × ≥ ≤ ≤ - ∈ t k i i y t n R t r n N k I K (10) , , 0 0 0 , , 1 , , 0 1, ≤ - ∈ = + × ≥ ≤ ≤ - ∈ t k t k i i i i q S IP t n R t r n N k I K (11) , , 0 0 ' , , , ≤ ∈ ∈ ∈ t k t k i i y y t k I T K (12) , , max , , , , ⋅ ≤ ⋅ ∈ ∈ ∈ ∈ t k t k i ij ij j i v x y C j t k I J T K (13) If , 1 , 0 = > ∑ J t k ij j t k i q I , , 1 , , , 0 ' , , \{0}, , =   = - - ∈ ∈ ∈ ∈     ∑ J t k ij j t k t k t k ij ij ij i i q q f q I j t k I J T K (14) If , 1 , 0 = ≤ ∑ J t k ij j t k i q I , , , ' , , \{0}, , = ∈ ∈ ∈ ∈ t k t k ij ij i q q j t k I J T K (15) If , 1 max ' = ⋅ > ∑ I t k i ij i j v q C , , , 1 , max '
' ' , = ⋅ ⋅ ⋅   = - -     ∑ I t k t k i ij i i ij i t k ij ij j v x v v q f q C , \{0}, , ∈ ∈ ∈ ∈ i j t k I J T K (16) If , 1 max ' = ⋅ ≤ ∑ I t k i ij i j v q C , , , ' , , \{0}, , = ∈ ∈ ∈ ∈ t k ij t k ij x i q j t k I J T K (17) If , 0 1 max 0 = ⋅ > ∑ I t k i i i v q C , , , 0 0 1 , max 0 0 0 ' , = ⋅ ⋅ ⋅   = - -     ∑ I t k t k i i i i i i t k i i v x v v q f q C , 1 , , 0 1, ∈ = + × ≥ ≤ ≤ -∈ i t n R t r n N k I K (18) If , 0 1 max 0 = ⋅ ≤ ∑ I t k i i i v q C , , 0 , 0 , , 1 , , 0 1, = ∈ = + × ≥ ≤ ≤ - ∈ t k i t k i x i q t n R t r n N k I K (19) , , ≥ ∈ ∈ ij ij i S s j I J (20) 0, , ≥ ∈ ∈ ij i s j I J (21) , 0 0, , , 1 , 0 1, = ∈ ∈ ≠ + × ≤ ≤ - ∈ t k i i x t t n R n N k I T K (22) { } , , , , , , , , , , , , , ' 0;
, ' 0,1 , , , ,

t k t k t k t k t k t k t k t k ij ij ij ij ij ij ij ij i IP I LS x q q y y j t k ≥ ∈ ∈ ∈ ∈ ∈ I J T K (23)
The objective function (1) represents the total cost composed of inventory holding costs at all stocking locations, transportation costs (minor ordering costs) for inventory replenishments of all stocking locations, major ordering costs of the CDC, and lost sales costs at all FDCs for the distribution system considered. Equations ( 2) and ( 3) are inventory balance equations of the FDCs and the CDC, respectively. Equations ( 4) and ( 5) establish the relationship among inventory position, inventory level, and outstanding orders at each stocking location.

Constraints ( 6) to (11) ensure that an (s, S) inventory policy is used at each stocking location. Among them, logical constraints ( 6) and ( 7) imply that inventory replenishment of each product at each FDC only occurs when its inventory position drops below its reorder point, constraints (8) imply that each product's replenishment quantity at each FDC does not exceed its order-up-to level minus its inventory position. Similarly, logical constraints ( 9) and ( 10) imply that the inventory replenishment of each product at the CDC only occurs when its inventory position drops below its order-up-to level, constraints (11) ensure that the inventory replenishment quantity of each product at the CDC is no more than its order-up-to level minus its inventory position, where constraints ( 9) to ( 11) are imposed only for periods t

= 1 , + × ≥ n R t r with 0 1 ≤ ≤ -
n N , this is because the CDC replenishes its inventory only in its review periods.

Note that, in this model, the real replenishment quantity of a product at some stocking location in some period may be smaller than the order quantity determined by its (s, S) policy because of on-hand inventory rationing at the CDC and transportation capacity rationing among products jointly replenished.

Constraints (12) imply that any product is replenished at the CDC only if it is one of the products jointly replenished. Constraints (13) link real variable t ij x and binary variable t ij y .

Logical constraints ( 14)-( 19) describe the inventory rationing at the CDC in each period and the transportation capacity rationing at each distribution channel for each joint replenishment. For each product i, the inventory rationing constraints ( 14) and ( 15) indicate the replenishment quantity , ' t k ij q allocated to each FDC j determined by a linear rationing rule and according to on-hand inventory of the CDC without considering the maximum joint replenishment volume constraints. This quantity , ' t k ij q is a function of the replenishment quantity , t k ij q determined by an (s, S) policy, the shortage , 1 , 0

J t k ij j t k i q I = -

∑

of the CDC, and the rationing fraction ij f for each FDC j, where , 0 t k i I is product i's on-hand inventory at the CDC at the moment of rationing.

Note that only the linear inventory rationing equations ( 14) themselves may generate negative allocation , ' t k ij q , when the balance assumption does not hold for the distribution system considered [START_REF] Diks | Optimal control of a divergent multi-echelon inventory system[END_REF]. This problem is overcome by imposing the non-negativity of , ' t k ij q in (23) and allowing , , ≤ -

t k t k ij ij ij q S IP even if , ≤ t k ij ij IP s .

= -

t k t k ij ij ij q S IP to , , ≤ - t k t k ij ij ij q S IP in case of , ≤ t k ij ij
IP s . This relaxation is reasonable since the real replenishment quantity of a product at some FDC in some period may be smaller than the order quantity determined by its (s, S) policy as explained above.

The transportation capacity rationing constraints ( 16) and ( 17

I t k i ij i j v q C = ⋅ - ∑
, and the rationing fraction ' ij f for each product i.

Similarly, the rationing constraints ( 18) and ( 19) indicate the real replenishment quantity , 0 t k i x of product i for the CDC determined by a linear rationing rule of transportation capacity that considers the maximum joint replenishment volume max 0 C for each replenishment of the CDC, where , 0 t k i q is the quantity of product i replenished by the CDC determined by an (s, S) policy.

Similarly, only the linear transportation capacity rationing equations ( 16) or (18) themselves may generate negative allocation

≤ -

t k t k ij ij ij q S IP generally.
Constraints (20) imply that the order-up-to level of product i at stocking location j is always larger or equal than its corresponding reorder point. Since in this inventory model with lost sales, the inventory position , t k ij IP is never negative, constraints (21) are added for all products i and all stocking locations j without loss of generality.

Finally, ( 22) and ( 23) indicate the ranges of all decision variables. As a replenishment from its supplier to the CDC in each scenario happens only once every N periods, and the distribution system starts in the r-th period in the review cycle of N periods of the CDC, the replenishment quantity , 0 0 ≥

t k i x in the periods 1, = × + ≥ t n R t r with 0 1 ≤ ≤ - n N and , 0 0 = t k i x
in the other periods.

All logical constraints in the above model SBO can be equivalently transformed into linear constraints by introducing auxiliary variables and big numbers. After the transformation, this model becomes a MILP model. In the remainder of this section, we only linearize logical constraints ( 6), ( 7), [START_REF] References Andersson | A two-echelon inventory model with lost sales[END_REF], and (10) which will be cited and reformulated in the next section when describing a solution approach for the model. For the linearization of other logical constraints and the setting of big numbers introduced in the linearization of all logical constraints in model SBO, please see the appendix at the end of this paper.

By introducing big numbers

11 ij M and 12 ij M for each i and j, logical constraints ( 6) and ( 7) can be transformed into the following linear constraints ( 24) and ( 25).

( )

11 , , 1 , , \ {0}, , ⋅ -≤ - ∈ ∈ ∈ ∈ t k t k ij ij ij ij i M y s IP j t k I J T K (24) , 12 , , , \ {0}, , - ≤ ⋅ ∈ ∈ ∈ ∈ t k t k ij ij ij ij i s IP M y j t k I J T K (25)
Similarly, logical constraints ( 9) and ( 10) can be transformed into the following linear constraints ( 26) and ( 27) by introducing big numbers

21 i M and 22 i M for each i. ( ) 21 , , 0 0 0 1 , , 1 , , 0 1, ⋅ -≤ - ∈ = + × ≥ ≤ ≤ -∈ t k t k i i i i i M y s IP t n R t r n N k I K (26) , 22 , 0 0 0 , , 1 , , 0 1, - ≤ ⋅ ∈ = + × ≥ ≤ ≤ - ∈ t k t k i i i i i s IP M y t n R t r n N k I K (27)

SOLUTION APPROACH

After the linearization of all logical constraints, model SBO becomes a mixed integer linear program (MILP). This MILP can be solved by using a commercial MILP solver such as CPLEX. However, the time required by CPLEX to solve the model increases exponentially as the number of its integer variables increases, so when the number of scenarios is large, a more efficient method is required to solve this model. Because when sij and Sij are given for all products and all stocking locations, this model can be decomposed into K independent sub-models, where K is the number of scenarios, we try to design a scenario-based decomposition approach to solve it. In order to do so, we first reformulate SBO by introducing scenario-dependant variables k ij s and k ij S to replace sij and Sij and adding constraints

' '∈ = ∑ k k ij ij k s s K K and ' ' ∈ = ∑ k k ij ij k S S K K
, for any k ∈ Κ, where K is the number of scenarios in Κ.

The reformulated model is given as follows:

Model SBO2:

( ) 2 , , , , 0 0 \{0} ' ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ =     ⋅ ⋅ + ⋅ ⋅ + ⋅ + ⋅           ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ SBO t k t k t k t k k ij ij j i ij ij ij k t i j t t i j Z Min p h I t v x O y l LS K T I J T T I J (28)
subject to constraints (2)-( 5), ( 12)-( 13), ( 22)-( 23),

( ) 11 , , 1 , , \ {0}, , ⋅ -≤ - ∈ ∈ ∈ ∈ t k k t k ij ij ij ij i M y s IP j t k I J T K (29) , 12 , , , \{0}, , - ≤ ⋅ ∈ ∈ ∈ ∈ k t k t k ij ij ij ij i s IP M y j t k I J T K (30) , , , , \{0}, , ≤ - ∈ ∈ ∈ ∈ t k k t k ij ij ij i q S IP j t k I J T K (31) ( ) 21 , , 0 0 0 1 , , 1 , , 0 1, ⋅ -≤ - ∈ = + × ≥ ≤ ≤ -∈ t k k t k i i i i i M y s IP t n R t r n N k I K (32) , 22 , 0 0 0 , , 1 , , 0 1, - ≤ ⋅ ∈ = + × ≥ ≤ ≤ - ∈ k t k t k i i i i i s IP M y t n R t r n N k I K (33) , , 0 0 0 , , 1 , , 0 1, ≤ - ∈ = + × ≥ ≤ ≤ - ∈ t k k t k i i i i q S IP t n R t r n N k I K (34) , , , ≥ ∈ ∈ ∈ k k ij ij i S s j k I J K (35) ' ' , , , ∈ = ∈ ∈ ∈ ∑ k k ij ij k i s s K j k K I J K (36) ' ' , , , ∈ = ∈ ∈ ∈ ∑ k k ij ij k i S S K j k K I J K (37) 0, , , ≥ ∈ ∈ ∈ k ij i s j k I J K (38)
and constraints ( 46)-( 69) in the Appendix.

In model SBO2, the objective function ( 28) is the same as that of model SBO; constraints (29) to (34) are the reformulation of constraints ( 6) to (11), or equivalently the reformulation of constraints ( 24), ( 25), ( 8), ( 26 We first attempted to solve model SBO2 by Lagrangian relaxation that relaxes constraints (36) and (37) by introducing Lagrange multipliers. However, although such approach can obtain a good lower bound, we cannot construct a high quality feasible solution of the model from the solution of its Lagrangian relaxed problem/model, because for the relaxed model, its optimal values of some variables k ij s and k ij S may be infinite during the Lagrangian dual optimization process as the coefficients associated with those variables in its objective function may become negative. Of course, we can use the augmented Lagrangian relaxation method [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] to solve the infinite problem, but it requires solving all subproblems optimally in each iteration, which are quite time consuming, especially for large-size instances. Because of this, we finally chose Alternating Direction Method of Multipliers (ADMM, [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF] to solve the reformulated model. This method, which can avoid the above mentioned infinite problem, is computationally more efficient and can generate more stable and better solutions than the other two methods mentioned. ADMM has been successfully used in distributed convex optimization and machine learning [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF].

ADMM Procedure

By introducing Lagrange multipliers λ k ij and µ k ij and relaxing constraints ( 36) and (37) respectively, the ADMM relaxed model of SBO2 can be formulated as RSBO2 in the following, where the objective function (28) in model SBO2 is replaced by an augmented Lagrangian function (39) with penalty terms. Here, ρ is a small positive real number.

Model RSBO2: ( ) , , , , 0 0 \{0} 2 ' ' ' ' ' ' ' 2 ρ µ λ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ =   ⋅ ⋅ + ⋅ ⋅ + ⋅ + ⋅           + ⋅ - + ⋅ -                + ⋅ -   ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ t k t k t k t k k ij ij j i ij ij ij k t i j t t i j RSBO2 k k k k k ij ij ij ij ij k i j k k k k k ij ij ij k Z Min p h I t v x O y l LS s s K s s K S S K K T I J T T I J K I J K K K 2 ' ' 2 ρ ∈ ∈ ∈ ∈                        + ⋅ -                   ∑ ∑ ∑ ∑ k k ij ij k i j k S S K K I J K (39)
subject to constraints (2)-( 5), ( 12)-( 13), ( 22)-( 23), ( 29)-( 35), ( 38), and constraints ( 46)-( 69) in the Appendix.

ADMM is an iterative method. In each iteration, the relaxed model is approximately solved by applying a Gauss-Seidel like approach. That is, when ADMM optimizes the decision variables related to each scenario k in model RSBO2, all other decision variables' values are fixed, and all decision variables of the relaxed model are optimized scenario by scenario and sequentially. In this way, the relaxed model RSBO2 can be decomposed into K sub-models, denoted by RSBO2k, k ∈ Κ, where '

( ) 2 , , , , 0 0 \{0} 2 ' ' ' , ' ' , ' ' 2 ρ µ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ≠ ∈ ≠ =   ⋅ ⋅ + ⋅ ⋅ + ⋅ + ⋅                 + ⋅ - + + ⋅ - +                         ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ RSBO t k t k t k t k k ij ij j i ij ij ij k t i j t t i j k k k k k k k k ij ij ij ij ij ij ij k k k k k k Z Min p h I t v x O y l LS s s s K s s s K K T I J T T I J K K 2 ' ' ' , ' ' , ' 2 ρ λ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ≠ ∈ ≠                                 + ⋅ - + + ⋅ - +                               ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ k i j k k k k k k k ij ij ij ij ij ij ij k i j k k k k k k S S S K S S S K K I J K I J K K (40)
subject to constraints (2)-( 5), ( 12)-( 13), ( 22)-( 23), ( 29)-( 35

∈   = + ⋅ - ∈ ∈ ∈     ∑ k m k m k m k m ij ij ij ij k i s s K j k K I J K (41) , 1 , , ', ' , , , µ µ ρ +
∈   = + ⋅ - ∈ ∈ ∈     ∑ k m k m k m k m ij ij ij ij k i S S K j k K I J K (42) ', ' , , , λ λ ρ +
where , In summary, the iterative procedure of ADMM can be described as the following.

ADMM:

Step 0. Set Step 1. Update the values of Lagrangian multipliers k ij µ and k ij λ for all i, j and k according to equations ( 41) and ( 42), and set m + 1 → m.

Step 2. Solve all sub-models one by one with Lagrange multipliers Otherwise, go to Step 1.

The stopping criterion of ADMM can be defined in different ways, such as a given computation time is achieved, a predefined gap between an upper bound and a lower bound is achieved, or a given number of iterations is reached. Note that in each iteration of ADMM, an upper bound of the model SBO can be obtained by solving it with sij and Sij replaced by

* ij s = , ∈ ∑ k m ij k s K K and * ij S = , ∈ ∑ k m ij k S K
K respectively for all i, j, and a lower bound of the model can be obtained by solving the ADMM relaxed model RSBO2 with the quadratic penalty terms in its objective function (39) ignored.

Solving ADMM Submodels by Linear Approximation

Each sub-model RSBO2k to be solved in each iteration of ADMM is a mixed-integer quadratic program. Although all the sub-models can be solved by a solver of quadratic programs, it is quite time consuming if the number of stocking locations, the number of products, and/or the number of periods considered are large. For this reason, we propose a linear approximation approach to solve the sub-models more quickly. ( )

∈ ≠     = = - +           ∑ k k k k k ij s ij ij ij ij k k k u G s s s s K K and 2 ' ' , ' ( ) ∈ ≠     = = - +           ∑ k k k k k ij ij S ij ij ij k k k w G S S S S K K .
The quadratic terms in the objective function of RSBO2k can be written as

( ) ( ) 2 2 ρ ρ ∈ ∈ ∈ ∈ ∈ ∈ + ∑ ∑ ∑ ∑ ∑ ∑ k k s ij S ij k i j k i j G s G S K I J K I J . Since all ( ) k s ij G s and ( ) k S ij
G S are convex functions with a single variable and RSBO2k is a minimization problem, these functions can be well approximated from below by a set of linear inequalities.

Let

' ( )

k s ij G s and ' ( ) k S ij G S
be the first derivative of function ( )

k s ij G s and ( ) k S ij G S , respectively, we have ' ( ) k s ij G s = 2 ' 2 ' , ' 1 1 2 2 ∈ ≠ - -   ⋅ - ⋅     ∑ k k ij ij k k k K K s s K K K , and ' ( ) k S ij G S = 2 ' 2 ' , ' 1 1 2 2 ∈ ≠ - -   ⋅ - ⋅     ∑ k k ij ij k k k K K S S K K K . Function ( ) k s ij G s and ( ) k S ij
G S can be approximated from below by the following two sets of linear inequalities ( 43) and ( 44), respectively. ( ) 43) and ( 44) to this model at the beginning. This MILP model is solved once to obtain the optimal solution of RSBO2k. In our numerical experiments in the next section, the second way is adopted.

' ( ) ( )( )

k k ij s ij s ij ij ij u G s G s s s ≥ + - , ij s ∈ s ij Π , , ∈ ∈ i j I J (43) ' ( ) ( )( ) ≥ + - k k ij S ij S ij ij ij w G S G S S S , ij S ∈ S ij Π , , ∈ ∈ i j I J ( 44
2 , , , , 0 0 \{0} ' ' , ' ' ' ' 2 ρ µ λ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ∈ ≠ ∈ =   ⋅ ⋅ + ⋅ ⋅ + ⋅ + ⋅           + ⋅ - + + ⋅                 + ⋅ - + ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ RSBO t k t k t k t k k ij ij j i ij ij ij k t i j t t i j k k k k k k ij ij ij ij ij k i j k k k k k k k ij ij ij ij k Z Min p h I t v x O y l LS s s s K u S S S K T I J T T I J K I J K K , ' 2 ρ ∈ ∈ ∈ ≠                         + ⋅                       ∑ ∑ ∑ ∑ k ij k i j k k K v K I J (

NUMERICAL RESULTS

In this section, the performances of model SBO and its ADMM algorithm are evaluated on four sets of instances generated according to real data from Alibaba, and the (s, S) policies found by the algorithm are compared with those used in Alibaba for these instances.

The parameters used to generate the instances are set in the following ways: the number of products I is set to 5 or 10, the number of FDCs J is set to 9, and the number of periods in the planning horizon is set to 30 periods (days). In addition, the review period of each FDC and that of the CDC are 1 period and 3 periods respectively. The replenishment lead time of each stocking location is taken as an integer in [0, 2]. This lead time may be 0 because the locations of some suppliers and FDCs are close to the CDC, so that goods can be delivered from those suppliers to the CDC or from the CDC to those FDCs on the same day. The parameter r, i.e., the starting period of the planning horizon, is an integer randomly taken from [0, 3], and the ending period of the planning horizon is taken as r+30-1.

Since only fast moving goods are involved in the distribution system, we assume that each product's demand observed by each FDC in each period (day) is stationary and normally distributed, with mean and standard deviation generated from historical demand data of Alibaba. For the demand of each product, its mean is in the interval [18,239], and its coefficient of variation is in the interval [0.1, 0.4]. Each product's volume is ranged from 0.001 to 0.01 m 3 . For the case with I = 5, the joint replenishment volume of each FDC and that of the CDC are limited to 5 m 3 and 50 m 3 , respectively, whereas for the case with I = 10, this volume is limited to 10 m 3 and 100 m 3 for each FDC and the CDC respectively. In addition, the unit holding cost of each product in each day at each stocking location is set based on the price of the product, and the unit lost sales cost of each product at each stocking location is set as the unit holding cost × 0.95 / (1-0.95) = 19 × the unit holding cost, where 0.95 = 95% is the expected service level. For each stocking location, its transportation cost for each replenishment is set based on its volume in m 3 and its transportation distance in km. For this distribution system, all FDCs and the CDC are located in a region with the distance between any two stocking locations ranged from 80 to 150 km, and the transportation cost per m 3 per km is estimated as 1 RMB Yuan. Accordingly, the transportation cost for the replenishment of each m 3 by each stocking location is in the interval [80,150]. The major ordering cost for each joint replenishment of the CDC is set between 100 and 2000 RMB Yuan. For the reason of confidentiality, price data of the products are not provided in detail here.

Let ij µ denote the mean and ij σ the standard deviation of product i's demand observed by stocking location j (j = 0 for the CDC and j = 1, 2, …, J for the FDCs), where 0

1 J i ij j µ µ = = ∑ and 2 0 1 J i ij j σ σ = = ∑ .
The fair share rationing rule introduced in Eppen and Schrage (1981) is applied for both the allocation of on-hand inventory of the CDC among the FDCs and the allocation of transportation capacity among multiple products for each joint replenishment, with the rationing fractions ij f and ' ij f given as In the planning horizon, each product's initial inventory position and final inventory position at each stocking location is set to its order-up-to level. For each scenario, each product i's demand observed by each FDC j in each period t was randomly generated according to its normal distribution by using the Monte Carlo Sampling method.

1 J ij ij ij j ij ij j j f L R L R σ σ = = + + ∑ for each product i and 1 ' σ σ = = + + ∑ I ij i ij ij j i ij ij j i f v L R v L R
µ σ + + + ij j ij ij ij ij j L R z L R and max ij S = max ( ) ( ) µ σ + + + + ij j ij ij ij ij j j i L R z L R C v ,
t ij F -is the inverse function of t ij F .
With this scenario generation method, all scenarios generated are assumed to have the same probability of occurrence, i.e., the probability of each scenario is 1/K, that is 1

= k p K for any k = 1, 2, …, K,
where K is the number of scenarios generated.

The number of scenarios K is set by applying the statistical analysis method introduced in [START_REF] Shapiro | Stochastic programming by monte carlo simulation methods[END_REF] and [START_REF] Nguyen | Supplier selection and operation planning in biomass supply chains with supply uncertainty[END_REF]. This method calculates the number of scenarios K required according to a given level of confidence α and a given length of confidence interval

H, i.e., (
)

2 /2 0 0 max( ( ) , ) K z K H K α σ = ⋅
, where 0 ( ) x z α ≤ ) =1-α/2, where x∼N(0, 1). We set K0 = 100 and tested 10 randomly generated instances with 5 and 10 products, respectively. After conducting the statistical analysis of their results, we observed that our scenario-based MILP model requires 55 and 339 scenarios on average for the instances with 5 products and 35 and 211 scenarios on average for the instances with 10 products to ensure their solution within the deviation 5% and 2% respectively from their true optimal solution in terms of cost at a probability no less than 99.9%. Based on this observation, we set the number of scenarios to 100 and 500 respectively for each instance in our numerical experiments.

σ K = ( ) 0 2 0 1 ( ) ( 1) = - - ∑ K SBO SBO k Z Z k K ,
Four sets of 10 instances were generated. Each set is indicated by the number of products and the number of scenarios. For example, the set 5x100 contains 10 instances with 5 products and 100 scenarios, and the set 10x500 contains 10 instances with 10 products and 500 scenarios.

In the numerical experiments, all MILP models were solved by calling the solver of CPLEX 12.9 on a personal PC with a CPU of i7-8565U CPU and a RAM of 16GB. The parameter ρ in ADMM is set to 10 -6 for all instances tested. Moreover, ADMM is terminated when a predefined relative gap 10% between an upper bound and a lower bound is achieved. The two bounds are obtained by ADMM and the Lagrangian method described at the end of section 5, respectively. Table 1 to Table 8 show the computational results of the tested instances.

To evaluate the efficiency and effectiveness of our ADMM, we compare it with CPLEX MIP solver for the instances sets 5x100, 10x100 and 5x500. For the instances set 10x500, this comparison is impossible because CPLEX MIP solver run out of memory when it solved these instances. The reason is that for each instance in set 10x500 with 10 products, 10 stocking locations, 30 periods and 500 scenarios, its model SBO contains 1500000 integer variables , t k ij y . However, our ADMM could be used to solve these instances because it requires much less memory due to its decomposition nature. Therefore, ADMM is more adapted to solve large-size instances. For each large-size instance tested, the quality of its solution found by ADMM is evaluated by the relative gap between its lower bound and upper bound obtained by ADMM (see section 5). These results are presented in Table 1 to Table 4. In these tables, row 'Costcplex' gives the best total cost of model SBO, i.e., the smallest objective value (1), found by CPLEX MIP solver with the computation time limited to 2 hours (7200s) for each instance in set 5x100 and 4 hours (14400s) for each instance in set 10x100 and 5x500. Row 'CostADMM' provides the best total cost of model SBO found by ADMM for each instance, CostADMM represents an upper bound of the optimal total cost of model SBO. 'CostLB' gives the lower bound of the total cost found by the Lagrangian relaxation method mentioned at the end of section 5. Row 'ReduceCost' gives the percentage cost reduction of the solution found by ADMM with respect to the solution found by CPLEX, i.e., (Costcplex -CostADMM) / Costcplex. 'Fillratecplex' indicates the average fill rate of all products at all FDCs of the solution found by CPLEX, and 'FillrateADMM' gives the average fill rate of all products at all FDCs of the solution found by ADMM. Row 'IncreaseFR1' represents the percentage increase of the average fill rate of the solution found by ADMM compared with the solution found by CPLEX, i.e., (FillrateADMM -Fillratecplex) / Fillratecplex. 'TimeADMM' provides the computation time (in seconds) of ADMM. Row 'ReduceTime' indicates the percentage computation time reduction of ADMM with respect to CPLEX, i.e., (Timecplex -TimeADMM) / Timecplex. Row 'GapADMM' represents the relative gap between the upper bound CostADMM and the lower bound CostLB.

Note that Table 4 only provides the results of ADMM since CPLEX run out of memory for all instances in set 10x500. In addition, in our numerical experiments, the fill rate of product i at FDC j indicated by ij fillrate is calculated as the average of , t k ij fillrate 's for all periods t under all scenarios k considered, where , t k ij fillrate is defined as the percentage of product i's demand (in units) fulfilled (satisfied) by FDC j in period t under scenario k. Because of the space limitation, we only report the average fill rate of all products at all FDCs for each instance, this average fill rate is obtained by averaging ij fillrate over products i and all FDC j; individual fill rate for each of the products at each FDC is not reported in the following tables. From the results of instance sets 5x100, 10x100 and 5x500 in Table 1 to 3, we can observe that ADMM could find a better solution in a much shorter time compared with the computation time of CPLEX MIP solver with relative cost reduction between 1.6% and 5.5% for set 5x100, between 1.1% and 58.3% for set 10x100, and between 46.6% to 66.5% for set 5x500. At the same time, ADMM could find a solution with a higher average customer fill rate than that of the solution found by CPLEX for all instances, with relative fill rate increase more than 1% for 38 of the 40 instances and maximum relative fill rate increase 226.3%. Note that 226.3% is not an absolute fill rate increase for an instance but the relative increase in percentage of the average fill rate of this instance obtained by ADMM compared with that obtained by CPLEX, because for this instance, the solution found by CPLEX is very poor. In addition, the computation time of ADMM is significantly less than that of CPLEX with relative time reduction between 60.3% and 95.8% for all instances. Moreover, CPLEX could not solve any instance in set 5x500 because of out of memory, but for all instances in this set ADMM could find a solution with the relative gap ranged from 3.8% to 8.2% as shown in Table 4.

To prove the superiority of the (s, S) policies found by our scenario-based optimization (SBO) approach for the distribution system considered, we compare them with those currently used in Alibaba. For ease of exposition, our approach and the approach used by Alibaba to define (s, S) policies in its distribution system are referred to as SBO and IPA respectively hereafter.

In IPA, for each product at each stocking location, its inventory policy parameters s and S are determined based on the mean μ and standard deviation σ of its demand, its replenishment lead time L, and its review period R. More precisely, the two parameters are defined as ( )

µ σ = ⋅ + + ⋅ ⋅ + s L R z L R
, and S = s + Q, respectively, where z is the z-value (safety stock factor) determined by the expected service level of the product at this stocking location and Q is determined by the mean μ and the maximum between the replenishment lead time L and the review period R, i.e., max(L, R). For the confidentiality reason, the formula for calculating Q is not detailed here.

We compare the two approaches by two criterions: the expected total cost and the average fill rate of the distribution system considered, where the expected total cost is defined as the objective value of model SBO and the average fill rate of the system is defined as the average of the fill rates of all products at all FDCs. Since it is very hard to derive analytical formulas for calculating the two criterions for the considered distribution system, we evaluate/estimate the two criterions by simulation of the system with its (s, S) policies set by either SBO or IPA. To keep consistence, we use the same demand scenarios (realizations of random demands) in the simulation of the system in both cases of SBO and IPA. These demand scenarios are the same as those used in the optimization model SBO.

In the following Table 5 to 8, row 'FillrateSBO' and 'CostSBO' provide respectively the average fill rate of all products at all FDCs of the distribution system and its expected total cost if it is controlled by the (s, S) policies found by SBO. Row 'FillrateIPA' and 'CostIPA' provide respectively the average fill rate of all products at all FDCs of the distribution system and its expected total cost if it is controlled by the (s, S) policies found by IPA. Row 'IncreaseFR2' indicates the percentage increase of the average fill rate obtained by SBO compared with that obtained by IPA, i.e., (FillrateSBO -FillrateIPA) / FillrateIPA. Row 'ReduceCost' indicates the percentage reduction of the expected total cost obtained by SBO with respect to that obtained by IPA, i.e., (CostIPA -CostSBO) / CostIPA. From Table 5 to 8, we can observe that SBO can find the (s, S) policies with a higher average customer fill rate than IPA for all instances, with relative fill rate increase ranged from 2.6% to 27.7% and 9.7% on average. At the same time, the inventory polices found by SBO can result in a lower total cost with respect to IPA for all instances, with relative cost reduction ranged from 1.8% to 9.6% and 5.1% on average. That is, the (s, S) policies found by SBO can achieve a larger fill rate at a smaller total cost with respect to those currently used in Alibaba. Therefore, SBO is more effective than IPA.

It should be noted that all numerical experiments were conducted on a personal PC, and the results in Table 1 to 8 are only used to verify the feasibility of SBO and the effectiveness of ADMM. Since ADMM is a decomposition algorithm, when it is applied to solve real instances in Alibaba, we can use parallel computing techniques and run it on multiple distributed high-performance cloud servers of Alibaba, since all sub-models RSBOk in each iteration of ADMM can be solved in parallel. This can dramatically reduce the running time of ADMM. By contrast, CPLEX MIP solver cannot realize this and has the out-of-memory problem for large-size instances.

Furthermore, we analyse in more detail the fill rates of each product at FDCs obtained by ADMM and IPA respectively on instance set 5x100. Table 9 and 10 provide the maximum and minimum fill rate of each product at all FDCs obtained by the two methods respectively. For example, the rows 'P1Max' and 'P1Min' give the maximal and minimal fill rate of product 1 at all FDCs obtained by ADMM or IPA. The row 'Fillrateave' gives the average fill rate of each product at all FDCs. From the two tables, we can observe that the fill rates of some products at all FDCs are very close, such as product 2 for instance 1, 4, 5, 6 and 9, and product 3 for instance 1, 2, 3, 5 and 8 (see Table 9). However, there are products whose fill rates at different FDCs are different, which may be because of the joint replenishment constraints. Compared with the method IPA currently used in Alibaba, SBO solved by ADMM can achieve a maximal fill rate very close to the best one and a much better minimal fill rate (in 35 of all 50 fill rates) for each product at each FDC. The minimal fill rates of some products obtained by IPA are very small, such as 29.5% for one product in instance 6, 16.31%, 27.66% and 28.1% for three products in instance 7, 15.99% for one product in instance 8, and 16.55% and 29.22% for two products in instance 10. Thus, there is a much larger variation of the fill rates of each product at all FDCs obtained by IPA. In practice, e-commerce companies as Alibaba usually seek for similar service levels for each product at all FDCs to achieve the best customer satisfaction. Therefore, SBO is more effective than IPA in achieving this goal. To further evaluate the performance of our ADMM, we compare the solution obtained by it with the optimal solution of model SBO for 10 instances of a distribution system with one warehouse (CDC), one retailer (FDC), and one product. The ten instances were randomly generated with the parameters set in the same ways as described above. The comparison results are given in Table 11, where the rows Costcplex, CostADMM, Timecplex, TimeADMM, Fillratecplex and FillrateADMM have the same meanings as the corresponding rows in Table 1-3 except that the optimal solution of model SBO for each instance is found by CPLEX. GapCost, GapTime and GapFR denote the cost gap, the computation time gap and the fillrate gap between the solution of ADMM and the optimal solution. From Table 11, we can see that our ADMM can find a solution very close to the optimal solution for most instances. The computation time of ADMM is much shorter than that of CPLEX for most instances. In the above numerical experiments, the planning horizon is set to 30 days based on the practice of Alibaba. To evaluate the effect of the planning horizon on the solution of model SBO, we compare the total cost of the distribution system obtained with the planning horizon of 30 days (case PH1) and that obtained with the planning horizon of 60 days (case PH2) respectively. More precisely, we consider 60 days of operation of the distribution system. In case PH1, the inventory policy optimization of the system is done twice with the planning horizon of 30 days for each SBO model, and the first and the second 30-day SBO model are solved independently, under the condition that the final state of the first 30-day model is turned into the initial state of the second 30-day model. The sum of the total costs of the two models are regarded as the total cost of the system in the 60 days of operation in case PH1, which is compared with the total cost of model SBO (the system) with the planning horizon of 60 days in case PH2. Table 12 shows the comparison results, from which we can see that the relative gap between the two total costs is no larger than 3.28% for all instances. Thus, it is reasonable to set the planning horizon to 30 days by considering both the effect of the planning horizon and the increase of computational complexity of model SBO if we take a longer planning horizon. Finally, we evaluate the effect of the end-of-horizon conditions on the solution of model SBO. We tested 10 instances of the distribution system with 5 products and 100 demand scenarios and compared its total cost obtained by ADMM with the end-of-horizon conditions for model SBO mentioned above (EC1) and its total cost obtained by the same solution method but with no end-of-horizon conditions for the model (EC2: the inventory position of each product at each stocking location at the end of the last period is not fixed). Table 13 presents the comparison results. From this table, we can find that the relative gap between the two total costs is very small for all instances. Therefore, the effect of the end-of-horizon conditions is minor.

CONCLUSION

In this paper, we consider a two-echelon distribution system in Alibaba with joint inventory replenishment constraints and have studied the optimization of its (s, S) policies. A novel scenario-based optimization model and an Alternating Direction Method of Multipliers (ADMM) method are proposed to solve this model. ADMM allows solving all subproblems in parallel in each iteration on multiple distributed personal computers or servers. We evaluated the model and the solution method on the instances generated according to real data of Alibaba. Our numerical results show the effectiveness of the model and the method. Compared with the MILP solver of CPLEX 12.9, our ADMM algorithm could find a better solution with a much shorter computation time and could avoid the out-of-memory for large-size instances. The inventory policies found by the algorithm can result in a better average fill rate with a lower expected total cost with respect to those used in Alibaba.

One direction for further study is to extend our proposed model and solution method to more complex distribution systems such as two-echelon distribution systems with non-stationary demands and three-echelon distribution systems. Besides, inventory replenishment planning for Alibaba's supply chain in some special seasons such as the season of double 11 promotion is worthy to be investigated. It is a great challenge for Alibaba to make cost-effective replenishment plans to achieve high customer fill rates in such seasons.

ACKNOWLEDGEMENT

APPENDIX I Linearization of the Model

In the following, we transform logical constraints ( 14) to ( 19) in model SBO into linear constraints so that it becomes a mixed integer linear program.

Logical constraints ( 14) and ( 15) can be transformed into linear constraints in the following way. Firstly, we introduce binary variables , 0

t k i z to indicate whether , 1 , 0 = > ∑ J t k ij j t k i q I , such that , 0 t k i z = 1 if , 1 , 0 = > ∑ J t k ij j t k i q I and , 0 t k i z = 0 if , 1 , 0 = ≤ ∑ J t k ij j t k i q I .
, , , = - ≤ ⋅ ∈ ∈ ∈ ∑ J t k ij j t k t k i i i i q I M z t k I T K (46) , 1 32 , , 0 0 ( 1) , , , = ⋅ -≤ - ∈ ∈ ∈ ∑ J t k ij j t k t k i i i i M z q I t k I T K (47)
Two cases may happen for the logical constraints. Firstly, in case without inventory rationing, i.e., , From constraints (54), ( 55), ( 59) and ( 64

≤ ∑ J t k ij j t k i q I , we have , , ' = t k t k ij ij q q . This 1 , 0 =

  i's demand observed by stocking location j in period t under scenario k ij L : product i's replenishment lead time at stocking location j j R : inventory review period of stocking location j with 1 = j R if j ≥ 1 and 0 = R R . That is, each FDC is reviewed once every period, and the CDC is reviewed once every R periods, R is a positive integer i v : volume of product i max j C : maximum joint replenishment volume of stocking location j ij h : unit holding cost of product i in each period at stocking location j ij l : unit lost sales cost of product i at stocking location j ij t : unit transportation cost of product i in any joint replenishment of stocking location j 0 O : major ordering cost for each joint replenishment of the CDC ij f : rationing fraction of FDC j for on-hand inventory allocation of the CDC among orders of the FDCs for each product i, fraction of product i for transportation capacity allocation of each joint replenishment of stocking location j among all products, 1 ' 1

  j replenishes its inventory of product i in period t under scenario k, and , the CDC makes a joint replenishment in period t under scenario k, and , , we can establish the following mixed-integer linear programming model for scenario-based optimization of the inventory policies of the distribution system with lost sales considered.

  x of product i for each FDC j determined by a linear rationing rule of transportation capacity that considers the maximum joint replenishment volume max j C for each replenishment of the FDC. This quantity is a function of ,

  ), (27), and (11); constraints (35) to (37) are the reformulation of constraints (20); constraints (38) are the reformulation of constraints (21). Since constraints (36) and (37) imply k ij s = sij and k ij S = Sij for any k for some sij and Sij, model SBO2 is equivalent to model SBO.

  of the Lagrange multipliers in each iteration of ADMM are updated according to the following equations (41) and (42).

  by solving sub-model RSBO2k in m-th iteration.

  all i, j and k. Solve all sub-models RSBO2k with the objective function (35) ignoring the penalty terms (i.e., the two terms with coefficient 2 . Set m = 0.

S

  obtained in the last or current iteration of ADMM.Step 3. If a given stop stopping criterion is achieved, output * i, j, calculate an upper bound of the optimal objective value of model SBO by solving it with sij and Sij replaced by * ij s and * ij S respectively for all i, j, and stop.

  approximation, sub-model RSBO2k can be reformulated as a mixed integer linear program by transforming its quadratic objective function into the following linear objective function (45) and adding two sets of linear constraints (43) and (44) to the set of original (existing) constraints of RSBO2k.

Π

  and their corresponding linear inequalities (43) and (44).One is to generate the two sets dynamically. That is, we first take an integer value between min ij s and max ij s for ij s and an integer value between min ij S and max ij S for ij S , add the corresponding constraints (43) and (44) to model RSBO2k, and replace its original quadratic objective function by linear objective function (45). The obtained MILP model is then solved by a MILP solver. If the optimal solution of the MILP model, denoted by for all i ∈ Ι and j ∈ ϑ, then s ij Π and S ij Π have been found, and this solution is the optimal solution of RSBO2k. Otherwise, , and the above procedure is repeated until an optimal solution of RSBO2k is found. The other way is to generate s ij the quadratic objective function of model RSBO2k by linear objective function (45), and add all linear constraints (

  at each stocking location j.The minimum possible value of ij s and ij S are both set to zero, i.e.

  used in the setting of some big numbers in model SBO as explained in subsection 3.3 and in the linear approximation approach used to solve all sub-models RSBO2k in ADMM.

  That is, for each scenario k, to generate ,

		t k ij d , a real number , t k ij u was first generated from the uniform distribution on [0, 1],
	, t k ij d is then generated as	1 ( ) ( ) , t t k ij ij F u -	, where	t ij F	(.)	is the probability distribution function of
	t ij d and	( )	1

Table 1 : Comparison of ADMM and CPLEX on instances 5x100

 1

	Instance	1	2	3	4	5	6	7	8	9	10
	Costcplex	54148.5 45927.3 72308.4 71311.6 45073.5 61723.4 78200.9 64957.3 49526.9 77990.7
	CostADMM	53123.4 44698.2 69741.9 70140.6 43832.6 59505.2 73924.8 63218.9	48125	75927.4
	ReduceCost	1.9%	2.7%	3.5%	1.6%	2.8%	3.6%	5.5%	2.7%	2.8%	2.6%
	Fillratecplex	88.7%	88.6%	87.4%	85.4%	92.1%	87.6%	80.7%	91.8%	89.9%	86.2%
	FillrateADMM	90.9%	90.3%	89.0%	87.9%	93.3%	88.8%	85.2%	92.6%	90.9%	87.4%
	IncreaseFR1	2.4%	1.9%	1.9%	2.9%	1.3%	1.4%	5.6%	0.9%	1.1%	1.5%
	TimeADMM	415	343	306	1361	350	435	427	365	396	342
	ReduceTime	94.2%	95.2%	95.8%	81.1%	95.1%	94.0%	94.1%	94.9%	94.5%	95.3%
	CostLB	50818.1 43690.6 67825.7 67988.1 41760.4 57779.3 72893.6 61693.2 46525.3	74875
	GapADMM	4.5%	2.3%	2.8%	3.2%	5%	3%	1.4%	2.5%	3.4%	1.4%

Table 2 : Comparison of ADMM and CPLEX on instances 10x100

 2

	Instance	1	2	3	4	5	6	7	8	9	10
	Costcplex	216714	235497	265637	253580	261575	99569.7 94146.3 96343.9	239737	91116.7
	CostADMM	110263	101235	130895	105792	142650	98101.7 92929.9 94704.2	109825	90126
	ReduceCost	49.1%	57.0%	50.7%	58.3%	45.5%	1.5%	1.3%	1.7%	54.2%	1.1%
	Fillratecplex	47.9%	50.6%	45.5%	50.5%	44.2%	94.9%	95.0%	94.4%	48.4%	92.2%
	FillrateADMM	90.9%	94.2%	91.6%	95.1%	86.8%	96.2%	96.2%	95.8%	93.5%	93.1%
	IncreaseFR1	89.9%	86.2%	101.2%	88.1%	96.6%	1.3%	1.3%	1.5%	93.3%	0.9%
	TimeADMM	1472	1121	1041	1240	1218	955	1043	1386	1463	1391

Table 3 : Comparison of ADMM and CPLEX on instances 5x500

 3

	Instance	1	2	3	4	5	6	7	8	9	10
	Costcplex	142471	128374	108164	113174	135716	120875	138493	137840	126383	138000
	CostADMM	68712.9 49656.2 39304.3 45540.8	72423	48179.9 46437.8 67227.4	59191	63151.8
	ReduceCost	51.8%	61.3%	63.7%	59.8%	46.6%	60.1%	66.5%	51.2%	53.2%	54.2%
	Fillratecplex	29.5%	40.2%	39.0%	36.7%	28.2%	39.9%	43.1%	30.6%	31.6%	30.3%
	FillrateADMM	96.2%	92.8%	96.4%	95.3%	89.5%	95.9%	94.8%	90.3%	90.0%	94.4%
	IncreaseFR1	226.3%	130.9%	147%	159.7%	217.3%	140%	120.1%	195%	185.1%	211.6%
	TimeADMM	2257	3767	1545	3851	4746	5719	2380	3451	4682	2353
	ReduceTime	84.3%	73.8%	89.3%	73.3%	67%	60.3%	83.5%	76%	67.5%	83.7%
	CostLB	65256.1 46492.2 36514.9 41974.6 67488.4 44649.1 43292.5 64465.4 56298.8 57611.4
	GapADMM	5.3%	6.8%	7.6%	8.5%	7.3%	7.9%	7.3%	4.3%	5.1%	9.6%

Table 4 : Computational results of ADMM on instances 10x500

 4

	Instance	31	32	33	34	35	36	37	38	39	40
	CostLB	90567.5	100740	125097	100226	112858	97113.2	134968	118209	112987	117087
	CostADMM	95520.9	106416	133027	105902	122107	104413	141952	124308	117320	124081
	GapADMM	5.5%	5.6%	6.3%	5.7%	8.2%	7.5%	5.2%	5.2%	3.8%	6%
	TimeADMM	11173	15097	18893	23851	13580	23922	16413	17749	16458	14856

Table 5 : Comparison of SBO and IPA on instances 5x100

 5

	Instance	1	2	3	4	5	6	7	8	9	10
	FillrateIPA	84.2%	85.8%	79.3%	80.0%	89.3%	82.7%	66.7%	79.7%	86.5%	73.8%
	FillrateSBO	90.9%	90.3%	89.0%	87.9%	93.3%	88.8%	85.2%	92.6%	90.9%	87.4%
	IncreaseFR2	7.9%	5.3%	12.3%	9.8%	4.5%	7.4%	27.7%	16.1%	5.1%	18.4%

Table 6 : Comparison of SBO and IPA on instances 10x100

 6

	Instance	11	12	13	14	15	16	17	18	19	20
	FillrateIPA	86.6%	90.9%	81.3%	91.3%	73.9%	93.1%	93.3%	92.8%	88.6%	89.4%
	FillrateSBO	90.9%	94.2%	91.6%	95.1%	86.8%	96.2%	96.2%	95.8%	93.5%	93.1%
	IncreaseFR2	5%	3.6%	12.7%	4.1%	17.6%	3.3%	3%	3.3%	5.5%	4.1%
	CostIPA	117286	107270	138506	114752	154252	101247	98783.8	100029	118629 95810.3
	CostSBO	110263	101235	130895	105792	142650 98101.7 92929.9 94704.2	109825	90126
	ReduceCost	6%	5.6%	5.5%	7.8%	7.5%	3.1%	5.9%	5.3%	7.4%	5.9%

Table 7 : Comparison of SBO and IPA on instances 5x500

 7

	Instance	21	22	23	24	25	26	27	28	29	30
	FillrateIPA	79.6%	86.5%	92.0%	89.0%	80.0%	92.0%	91.4%	81.6%	75.1%	80.3%
	FillrateSBO	96.2%	92.8%	96.4%	95.3%	89.5%	95.9%	94.8%	90.3%	90.0%	94.4%
	IncreaseFR2	20.9%	7.3%	4.8%	7.1%	11.8%	4.2%	3.8%	10.6%	19.9%	17.5%
	CostIPA	73510.6 51754.9	40153	47250.5 75942.9 49512.2 48662.8 68932.9 62010.6 64316.6
	CostSBO	68712.9 49656.2 39304.3 45540.8	72423	48179.9 46437.8 67227.4	59191	63151.8
	ReduceCost	6.5%	4.1%	2.1%	3.6%	4.6%	2.7%	4.6%	2.5%	4.5%	1.8%

Table 8 : Comparison of SBO and IPA on instances 10x500

 8

	Instance	31	32	33	34	35	36	37	38	39	40
	FillrateIPA	93.0%	84.8%	83.9%	84.7%	88.5%	91.1%	75.6%	0.822	81.4%	84.8%
	FillrateSBO	95.4%	93.5%	94.6%	93.1%	95.5%	95.1%	94.0%	0.921	90.0%	95.4%
	IncreaseFR2	2.6%	10.2%	12.7%	9.9%	7.9%	4.4%	24.3%	12.1%	10.6%	12.4%
	CostIPA	100098	111053	137970	110116	125884	109469	148445	128599	121266	128619

Table 9 : Fill rate of each product at FDCs obtained by ADMM on instances 5x100

 9

	Instance	1	2	3	4	5	6	7	8	9	10
	P1Max	99.95% 99.94% 99.77%	99.88%	99.01%	96.27%	100%	99.91%	99.93%	97.04%
	P1Min	61.71% 70.42% 70.25%	67.15%	59.79%	65.58%	82.55%	77.43%	99.70%	61.89%
	P2Max	100%	99.93% 99.91%	100%	99.97%	100%	100%	100%	99.93%	99.57%

Table 10 : Fill rates of each product at FDCs obtained by IPA on instances 5x100

 10

	Instance	1	2	3	4	5	6	7	8	9	10
	P1Max	100%	100%	99.77%	100%	100%	92.07%	100%	100%	100%	83.99%
	P1Min	34.06%	57.32%	36.44%	47.8%	39.36%	34.25% 57.25% 55.90% 99.97% 29.22%
	P2Max	100%	100%	100%	100%	100%	100%	100%	100%	100%	87.21%
	P2Min	99.998% 58.134% 35.133% 99.993% 99.996% 99.99% 37.48% 46.47% 99.98% 16.55%
	P3Max	100%	100%	100%	100%	100%	100%	100%	100%	100%	99.99%
	P3Min	97.39%	99.64%	99.99%	53.97%	99.99%	56.84% 27.66% 99.996% 33.74% 33.39%
	P4Max	97.23%	60.34%	100%	80.30%	100%	100%	99.21%	100%	100%	100%
	P4Min	39.18%	30.94%	37.10%	34.41%	44.14%	29.5%	28.1%	46.13% 99.99% 56.69%
	P5Max	100%	100%	100%	100%	100%	100%	99.99% 57.69%	100%	100%
	P5Min	60.93%	87.85% 99.9995% 44.48% 99.998% 99.996% 16.31% 15.99% 42.73% 99.99%
	Fillrateave 84.24%	85.77%	79.28%	80.04%	89.27%	82.74% 66.70% 79.72% 86.50% 73.84%

Table 11 : Comparison of the solution of ADMM and the optimal solution

 11

	Instance	1	2	3	4	5	6	7	8	9	10
	Costcplex	741.732 1280.22 1751.77 1996.48 902.055	449.7	2308.28 3282.67 803.649 1059.98
	CostADMM	743.693 1286.67	1802.9	2093.92 924.718 495.551 2379.96 3339.95 810.718 1066.72
	GapCost	0.3%	0.5%	2.9%	4.9%	2.5%	10.2%	3.1%	1.7%	0.9%	0.6%
	Timecplex	309	7204	7210	38	136	101	204	43	815	207
	TimeADMM	43	54	43	33	42	43	52	44	47	45
	GapTime	86.1%	99.3%	99.4%	13.2%	69.1%	57.4%	74.5%	-2.3%	94.2%	78.3%
	Fillratecplex	80.66% 76.48% 70.45% 63.06% 82.28% 98.71% 74.01% 59.89% 97.46% 71.32%
	FillrateADMM 82.86% 78.53%	75.7%	66.18% 86.93% 83.86% 79.58% 66.46%	87.8%	76.4%
	GapFR	2.7%	2.7%	7.5%	4.9%	5.6%	-15.0%	7.5%	11.0%	-9.9%	7.1%

Table 12 : Comparison of the results with different planning horizons

 12

	Instance	1	2	3	4	5	6	7	8	9	10
	PH1	2423.42 4187.6 6458.03 7794.42	3032.5	1064.453 7985.84 12890.57 1759.431 3621.15
	PH2	2505.53 4278.3	6480.1	7889.53 3000.37	1043.67	8178.63	12949.3	1703.52	3702.76
	Gap	3.28%	2.12%	0.34%	1.21%	1.07%	1.99%	2.36%	0.45%	3.28%	2.20%

Table 13 : Comparison of the results with different end-of-horizon conditions

 13

	Instance	1	2	3	4	5	6	7	8	9	10
	EC1	53123.4 44698.2	69741.9	70140.6 43832.6 59505.2 73924.8 63218.9	48125	75927.4
	EC2	53129.7 44692.9	69741.6	70136.7 43825.8 59499.6 73932.3 63229.6 48125.5	75927.2
	Gap	0.01%	-0.01% -0.0004% -0.01%	-0.02%	-0.01%	0.01%	0.02%	0.001% -0.0003%

  This work is supported financially by the Alibaba Innovative Research Project 'Optimization of Safety Stock Placement in Supply Chains with Demand and Lead Time Uncertainty' and the Natural Science Foundation of Hunan Province (No. 2022JJ30204).

Secondly, in case with transportation capacity rationing, i.e., if

. By introducing auxiliary variables

II Parameter setting of the big numbers

In the linear constraints ( 24) to ( 27) and ( 46) to (69) transformed from the logical constraints of model SBO, all big numbers must be defined.