Business, Sauder School of
Operations and Logistics (OPLOG), Division of
DSpace
UBCV
Aghababaei, Bahareh
Master of Science in Business Administration - MScB
University of British Columbia
One of the unique characteristics of online retail is the ability to ship items to customers from different fulfillment centers. This causes interdependence between fulfillment centers and creates new challenges in inventory management. Hence, this thesis investigates inventory management decisions in online retail, focusing on the initial stocking of fulfillment centers, their inventory levels, and order fulfillment processes. The problem is formulated as a dynamic integer program, with the objective of minimizing the sum of receiving, holding, and shipping costs while satisfying demand and capacity constraints. The large scale of the problem arises from the size of the fulfillment centers, the number of items, and the number of demand points. We propose a large-scale solution algorithm that aggregates the three sets of items, warehouses, and customers. The algorithm solves the aggregated problem with an exact solution that is tractable in size and then disaggregates the solution through three hierarchical disaggregation stages using a Greedy algorithm to address the NP-hardness of the problem. We develop a theoretical bound on the optimality gap caused by aggregation method. The numerical examples demonstrate that the large-scale algorithm is a computationally efficient approach which produce low-cost solutions for medium-sized clusters compared to the optimal model. Additionally, we also study the impact of key factors, including cost parameters and variations in item sizes on the optimality gap.
Inventory Management Optimization: A Case Study of Amazon Supply Chain by Bahareh Aghababaei A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN BUSINESS ADMINISTRATION in The Faculty of Graduate and Postdoctoral Studies (Transportation and Logistics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2024 © Bahareh Aghababaei, 2024 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled: Inventory Management Optimization: A Case Study of Amazon Supply Chain submitted by Bahareh Aghababaei in partial fulfilment of the requirements for the degree of Master of Science in Business Administration in Transportation and Logistics Examining Committee: Tamar Cohen-Hillel, Assistant Professor, Operations and Logistics, Sauder School of Business, UBC Business S pervi or Harish Krishnan, Professor, Operations and Logistics, Sauder School of Business, UBC Supervisory Committee Member Julia Yan, Assistant Professor, Operations and Logistics, Sauder School of Business, UBC Supervisory Committee Member Ning Nan, Associate Professor, Accounting and Information Systems, Sauder School of Business, UBC Additional Examiner iii Abstract One of the unique characteristics of online retail is the ability to ship items to customers from different fulfillment centers. This causes interdependence between fulfillment centers and creates new challenges in inventory management. Hence, this thesis investigates inventory management decisions in online retail, focusing on the initial stocking of fulfillment centers, their inventory levels, and order fulfillment processes. The problem is formulated as a dynamic integer program, with the objective of minimizing the sum of receiving, holding, and shipping costs while satisfying demand and capacity constraints. The large scale of the problem arises from the size of the fulfillment centers, the number of items, and the number of demand points. We propose a large-scale solution algorithm that aggregates the three sets of items, warehouses, and customers. The algorithm solves the aggregated problem with an exact solution that is tractable in size and then disaggregates the solution through three hierarchical disaggregation stages using a Greedy algorithm to address the NP-hardness of the problem. We develop a theoretical bound on the optimality gap caused by aggregation method. The numerical examples demonstrate that the large-scale algorithm is a computationally efficient approach which produce low-cost solutions for medium-sized clusters compared to the optimal model. Additionally, we also study the impact of key factors, including cost parameters and variations in item sizes on the optimality gap. iv Lay Summary This thesis addresses the challenge of optimizing inventory management in online retail, focusing on initial stocking, inventory replenishment, and planned order fulfillment. We aim to minimize costs associated with receiving, holding, and outbound shipping while meeting customer demand and capacity constraints. We develop a large-scale framework that integrates an aggregation approach with a Greedy algorithm to create a trade off between computation time and optimality gap. The analysis indicates that this approach provides near-optimal solutions quickly and efficiently, particularly for medium-sized problems, making it a robust alternative to traditional methods that struggle with large datasets. Additionally, our study provides theoretical insights into the algorithm's effectiveness and explores how factors like cost parameters and variations in item sizes impact performance. v Preface This thesis is original, unpublished, independent work by the author, Bahareh Aghababaei. It has benefited greatly from the extensive feedback and insightful guidance provided by Professor Tamar Cohen-Hillel. vi Table of Contents Abstract ..................................................................................................................................... iii Lay Summary ............................................................................................................................ iv Preface ........................................................................................................................................ v Table of Contents ...................................................................................................................... vi List of Tables ........................................................................................................................... viii List of Figures ........................................................................................................................... ix Acknowledgements .................................................................................................................... x 1. Introduction ......................................................................................................................... 1 2. Literature Review................................................................................................................ 4 3. Problem description ............................................................................................................ 8 3.1. Complexity of problem ............................................................................................. 10 4. Approximation Algorithm ................................................................................................. 15 4.1. Greedy Algorithm ..................................................................................................... 15 4.3. Greedy performance.................................................................................................. 18 5. Aggregation algorithm ...................................................................................................... 21 5.1. Aggregation formulation ........................................................................................... 23 5.2. Disaggregation steps and formulations ..................................................................... 24 5.2.1. Optimal Disaggregation of demand zones ........................................................ 26 5.2.2. Optimal Disaggregation of warehouses ............................................................ 27 5.2.3. Optimal Disaggregation of items ...................................................................... 29 5.3. Solution approaches .................................................................................................. 30 5.3.1. Exact method for demand zone disaggregation ................................................ 30 5.3.2. Approximation method for warehouse disaggregation ..................................... 31 5.3.3. Approximation method for item disaggregation ............................................... 33 5.4. The algorithm error bound ........................................................................................ 35 5.4.1. Proof of feasibility ............................................................................................ 36 5.4.2. Proof of optimality gap ..................................................................................... 36 6. Computational experiments .............................................................................................. 44 6.1. Simulation description .............................................................................................. 44 vii 6.2. Effect of exclusive cluster size on optimality gap .................................................... 46 6.3. Effects of cluster count on optimality gap and computation time ............................ 50 6.4. Effects of item size variations on optimality gap and running time ......................... 51 6.5. Effects of objective function parameter ratios on algorithm performance ............... 53 6.6. Comparison with the greedy algorithm..................................................................... 54 6.7. Comparison with the linear programming relaxation solution ................................. 56 7. Conclusion ........................................................................................................................ 59 References ................................................................................................................................ 62 viii List of Tables Table 1. The indices and notations used in the optimal model .................................................. 9 Table 2. The summary of greedy algorithm for solving optimal model................................... 16 Table 3: The lists of sets and notations used in the warehouse disaggregation problem ......... 27 Table 4. The summary of greedy algorithm for each cluster (w) ............................................. 32 Table 5. The summary of the greedy algorithm for each cluster (c) ........................................ 33 Table 6. Sensitivity analysis on the ratio of objective function parameters ............................. 54 Table 7. Comparison with LP relaxation solution .................................................................... 58 ix List of Figures Figure 1. The proposed methodology ........................................................................................ 3 Figure 2. The detailed workflow of aggregation and disaggregation stages ............................ 26 Figure 3. Optimality gap versus cluster size |C| ....................................................................... 47 Figure 4. Optimality gap versus cluster size |K| ....................................................................... 49 Figure 5. Optimality gap versus cluster size |W| ...................................................................... 49 Figure 6. Distribution of demand points and fulfilment centres .............................................. 50 Figure 7. Optimality gap and computation time versus cluster size ........................................ 51 Figure 8. The effect of sparsity in item sizes on optimality gap and running time .................. 53 Figure 9. Comparison of the algorithm with greedy performance ........................................... 56 x Acknowledgements The author acknowledges that the work for this project took place on the traditional, ancestral, and unceded territory of the xwməθkwəy̓əm (Musqueam) People. I would like to thank the faculty, staff, and students at UBC who have enriched my two years of learning here. I am profoundly grateful to my supervisor, Professor Tamar Cohen-Hillel, for the opportunity to collaborate on this challenging project. Her unwavering support, both moral and financial, and her faith in my abilities were instrumental in the completion of this work. Professor Cohen-Hillel's insightful guidance, combined with her extensive experience in the Amazon supply chain, provided invaluable perspectives and inspired me to continue pursuing research in this field. I also appreciate the time and constructive feedback from Professor Harish Krishnan and Professor Julia Yan. My sincere thanks to Professor Joseph Paat for his exceptional teaching in combinatorial optimization and for addressing my numerous questions with clarity. Special thanks go to my husband for his constant support throughout my educational journey at UBC. 1. Introduction The increasingly competitive retail business demands greater customer satisfaction and operational efficiencies, including optimal inventory and logistics planning (Hübner et al., 2013). For example, lack of optimal initial stock planning could lead to stock reduction resulting in shortage costs, delays, and eventually missed customers (Seyedan et al., 2022). To avoid such issues and maintain desired service levels, retailers often tend to keep large inventories, incurring substantial storage costs, or, in some cases, paying additional logistics costs to fulfill customer demand from distant warehouses (Salam et al., 2016). Consequently, in order to meet customer demand at the optimal time and with the minimum possible cost, the retail business must optimize initial stocking and fulfillment processes at the operational level, taking into account various factors affecting the planning. Moreover, large retail industries, such as Amazon, face complex planning decisions involving millions of individual items and numerous facilities to satisfy thousands of demand points. Amazon's fulfillment network consists of three key facility types: fulfillment centers, sortation centers, and delivery stations. Fulfillment centers are extensive facilities specializing in handling large volumes of diverse orders shipped as parcels, and categorized by item types (e.g., apparel, electronics) and sizes (small sortable, large sortable, large non-sortable). Sortation centers, large-sized facilities, organize shipments based on regional or local destinations, sorting parcels into smaller batches by postal code. Delivery stations, medium-sized facilities, focus on sorting parcels for specific local delivery routes, loading them onto appropriate vehicles for efficient last-mile delivery (Rodrigue, 2020). According to MWPVL International dataset, there are currently 376 active fulfilment centers, 111 sortation centers, and 566 delivery stations operating only in the United States in 2024 (Amazon Distribution Network Strategy | MWPVL International, n.d.). With warehouse capacity limitations, supply fluctuations, and the fact that different demand points vary in demand patterns for the same items, managing this extensive network, which aims to satisfy thousands of demand points, is a challenging task. The complexity of the problem is further increased by the fact that we are dealing with discrete optimization, and the integrality of the solution necessitates solving a very big combinatorial problem. As a result, the inventory and fulfilment planning for such an extensive network is hard. 2 In this work, our objective is to address the inventory management decisions for an online retailer through a proposed heuristic algorithm. The primary goal is to determine a cost-minimizing plan that simultaneously address three key decisions: initial inventory placement, inventory replenishment, and planned order fulfilment. The model is developed on a large-scale two-tier network, comprising fulfilment centers and customers, over a multi-period planning horizon. Our algorithm consists of a four-stage approach, wherein the first stage employs an auxiliary model utilizing aggregation analysis to reduce the dimension and complexity compared to the original model. Aggregation analysis incorporates cluster analysis techniques and additionally employs the method of combining similar data. The process of identifying which data are sufficiently similar to be grouped together is called cluster analysis, and it involves a large portion of the aggregation techniques that should be considered for the optimization models. These clusters may frequently be predefined in a wide range of optimization issues (Rogers et al., 1991). In this case, for example, clustering is applied to the three sets of demand points, items, and warehouses due to their huge sizes, showing that the entities to be grouped are variables during the aggregation stage. The next sections will provide comprehensive descriptions of the aggregation step. The next three steps include an attempt to recover the original solution through three stages of disaggregation using the outcomes of the reduced model's solution. Disaggregation analysis refers to techniques for taking an aggregated model and using it to derive components of a more refined model (usually with the same dimensions as the original model) (Rogers et al., 1991). In this context, we design a hierarchical disaggregation method in which demand zones are firstly disaggregated, followed by warehouses, and finally items. The first stage of this process is carried out optimally, while the second and third stages are completed approximately with the aid of a greedy algorithm. Our approach is depicted in Figure 1. The main contributions of this work are threefold. First, we develop an aggregation and disaggregation algorithm to address the NP-hard nature of the inventory and fulfillment planning model for the Amazon supply chain. Second, we apply a greedy algorithm which produces near-optimal solutions to solve the warehouse and item disaggregated models given our parameters and solution obtained from the demand zones disaggregated model. The greedy algorithm is capable of solving in polynomial time. Finally, we provide computational 3 experiments on simulated instances of the Amazon distribution network to demonstrate the algorithm's performance. Figure 1. The proposed methodology 4 2. Literature Review In the retail industry, businesses use a variety of sales channels, including online platforms, physical stores, and a combination of both within their supply chains. Online sales have been experiencing rapid annual growth. For instance, while it includes a 15.4 percent share of total US sales in 2023, it experienced a considerable 7.6 percent rise over the previous year (U.S. Census Bureau, 2023). One of the key aspects of online retail is its capability to fulfill customer orders from various fulfillment centers, regardless of their geographical proximity, each with its associated shipping costs (DeValve et al., 2021). While the shipping flexibility in online retail allows for cost savings and improved service quality, it also needs the development of inventory optimization methods to manage the connection between fulfillment centers. To achieve this, a range of decisions must be made depending on the planning horizon, such as network design, sourcing, inventory placement, replenishment, fulfillment, and transshipment (Ge et al., 2019). Given the interconnected nature of these decisions, authors have suggested a number of approaches to optimize them, either independently in a hierarchical manner or jointly. Chen & Graves (2021), for instance, focus on optimizing initial inventory placement in online retail, developing a cost-minimizing plan for each item while considering factors such as fulfillment center capacity constraints and customer satisfaction. In another study, they also explore operational decisions to address challenges arising from sorting and picking multiple-item orders within a warehouse. Their approach involves analyzing the problem from batch scheduling and routing perspectives, leading to a proposed picking policy that optimizes efficiency in both picking and sorting tasks Chen, (2018). Order fulfillment, a bottleneck for online retailers, has gained attention in research. In a recent study, Arlotto et al. (2023) investigate a joint inventory placement-fulfillment problem via regret analysis. At the beginning, the inventory of a single item is dispersed across different warehouses, then the decision maker determines whether to accept or reject the order at each period. They provide heuristic policies and evaluate their effectiveness in terms of minimum-inventory regret to the problem. Torabi et al. (2015) addressed the joint fulfillment and transshipment decisions by employing a mixed integer programming model with Benders decomposition to efficiently find the optimize solutions. Some studies define the order fulfillment process into three decisions: order 5 allocation, order consolidation and synchronization, and order delivery. In this context, Li et al. (2019) integrate these decisions into an online retailer's fulfillment plan and develop an adaptive large neighbourhood search method to handle large-scale instances of this problem. Building upon this work, Jiang & Li (2020) extend it by incorporating time windows to minimize transportation costs, as well as introducing a decomposition-based approach to reduce the computational complexity of solving large instances. Several recent studies address inventory management decisions within an omnichannel setting (Abouelrous et al., 2022; Wang et al., 2022). Bretthauer et al. (2010) determine optimal inventory allocation across different sites to meet in-store and online demand while minimizing overall costs, including holding, backorder, operating, transportation, and handling costs. Govindarajan et al. (2021) develop network-based strategies considering online demand spillover, where inventory decisions are made at the beginning of the selling horizon and fulfillment decisions are made iteratively. They derive inventory heuristics based on a two-stage approximation. In addition, their fulfillment heuristic provides location-specific inventory thresholds that determine how much inventory is rationed between in-store and online demands. The goal of this research is similar to that of Chen & Graves (2021), who study the problem of inventory placement across fulfilment centers considering capacity restrictions. However, we specifically focus on variable costs related to receiving, storage, and shipping rather than the fixed charge related to assigning each item to a fulfillment center. Additionally, we aim to simultaneously optimize the following decisions to address their interdependent nature: 1) inventory placement or initial stocking (determining which fulfillment centers store each item), 2) inventory replenishment (deciding how much inventory to hold), and 3) planned order fulfillment (selecting which fulfillment center to use for an order). If we ignore the fact that items differ in dimensions and each occupies a specific space in a warehouse, the problem reduces to the integral flow problem, which has been proven to be an easy one given the fact that the matrices forming the set of constraints has a total unimodularity property (Conforti et al., 2014a; Hoffman et al., 1956). In other words, it would be a special case of a transportation problem on a bipartite graph, where two sets of nodes correspond to fulfillment centers and customer regions, respectively (Hitchcock, 1941). The problem, however, is a generalized assignment problem (GAP) which falls into the category of multiple knapsack 6 problems (MKP), in which the number of identical copies available for each object type is restricted, each item j that is assigned to knapsack 𝑖 yields profit 𝑃𝑖𝑗 (rather than profit Pj). This problem has been shown to be np-complete (Kellerer et al., 2004) and will be further discussed in the subsequent sections. According to Martello & Toth (1990), there are three types of approaches to solving generalized assignment problems: 1) exact algorithms, such as Variants of branch-and-bound (Sarin et al., 2014; Woodcock & Wilson, 2010), 2) heuristic algorithms (Fathollahi-Fard et al., 2023; Saeedi et al., 2024), and 3) approximation algorithms (Cohen et al., 2006). When GAP is combined with a network flow problem, its complexity increases, particularly with a substantial number of nodes. Hence, in online retail networks, which involve numerous demand points, items, and warehouses, it is critical to employ appropriate techniques to handle the inherent intractability caused by both NP-hardness and the scale of the problem. Various approaches, such as decomposition-based algorithms like those discussed by Dell’Amico et al. (2019), are employed to handle large-scale linear integer programming models. However, the effectiveness of these algorithms depends on the specific characteristics of the problems, including variables and constraint structures. Despite their practical efficiency, they may not necessarily guarantee a solution within polynomial time. In this situation, it could be necessary to use some approximation techniques. The heuristic approach we provide here is the integration of aggregation technique and a greedy algorithm to handle disaggregation stage. Aggregation is an efficient technique to facilitate the analysis of large optimization models and provide a set of tools to transform the model with a high degree of complexity to a reduced manageable form. The simplified model can then be solved, and its optimal solution disaggregated into a feasible solution for the original problem (Å. Hallefjord et al., 1993). Rogers et al. (1991) establish a framework for aggregation and disaggregation methodology and provide a comprehensive review of previous research and applications in various areas of operations research. Later, Litvinchev & Tsurkov (2003) delve into aggregation analysis by reviewing the findings from previous decades, different aggregation techniques, and discussing aggregation error and bounds. Aggregation, as a technique for reducing computational complexity, has been widely applied in many contexts, particularly in those involving a large number of variables and constraints. For example, applications of aggregation theory in optimization include addressing 7 the classical transportation problem by aggregating demand points (Geoffrion, 1976), supply points (Evans, 1979), and both (Zipkin, 1980a). Aggregation has been also employed to reduce the number of items in a multi-commodities network flow problem (Kazemi, 2021), production planning (Leisten, 1998), and online retail inventory management (Chen & Graves, 2021). In this paper, we aim to perform aggregation on nodes and items simultaneously within an integer programming framework. To guarantee feasibility, we develop a problem-specific hierarchal disaggregation procedure which consists of three stages; the first stage is executed optimally, while the second and third ones employ a greedy algorithm to handle their NP-hardness. One main part of the aggregation approach involves quantifying the degree of accuracy lost under specific aggregation. Using duality theory, Zipkin provided bounds on the loss in objective function value that results from column (variables) aggregation (Zipkin, 1980c), and row (constraints) aggregation (Zipkin, 1980b) in linear programming models. A. Hallefjord & Storoey (1986) remark that Zipkin's results can be generalized to apply to discrete optimization problems as well. Here, our theoretical bound is based on input parameters and is specifically obtained from the final disaggregation stage. 8 3. Problem description In this study, we address the joint optimization of initial stocking, inventory replenishment and planned order fulfilment problem for an online retailer across a multi-period planning horizon. The goal is to determine a cost-minimizing plan which meets demand and supply requirements while taking fulfillment center capacity limitations into account. As previously stated, the retailer operates and manages a two-echelon order fulfilment system. The first tier of this system consists of a set of fulfilment centres (FCs), where retailer’s products are stored, while the second tier comprises individual customer nodes. Although the real-world system also involves a full network of sortation centers and delivery stations between these two layers, our study does not consider these intermediate facilities. Given the projected customer demand for each product, the retailer needs to decide on how to best place items in the fulfilment centers and fulfill customer orders, aiming to minimize overall expenses. Our focus lies on the variable costs associated with both inbound and outbound logistics. Inbound costs encompass expenses related to assigning items to fulfillment centers and storage, while outbound costs involve shipping expenses from fulfillment centers to customer demand points. The item assignment costs to fulfillment centers include any associated receiving expenses such as labor and inbound delivery costs. In addition, shipping costs from fulfillment centers to each customer point are computed based on the distance between nodes and a unit rate per distance, along with any additional handling charges incurred by the retailer. If the retailer employs a third-party service like UPS for item delivery, this rate can reflect the fee paid by the retailer to the carrier, or it can be considered a fuel rate if the retailer utilizes its own vehicles and resources for transportation. To ensure a balance between practicality and tractability, we make the following assumptions: • Transshipment between fulfillment centers (FCs) is prohibited to avoid additional costs in order fulfillment. • We exclude location strategic planning considerations; network structure and sourcing plans are predetermined so that the number and the locations of fulfilment centers and demand points, fulfillment center capacities, and supply amounts are given as input. 9 • Every customer demand must be satisfied within each period; there are no shortages in terms of backorders or lost sales. • Demand and supply values are deterministic but dynamically changes over time. • Edge capacities are unrestricted; there are enough vehicles available for transferring products from FCs to demand zones. Given these assumptions, we formulate the problem of initial stocking and order fulfillment for an online retailer as an integer linear program. In the following, we outline the sets, parameters, and variables defining the proposed model. Table 1. The indices and notations used in the optimal model Sets 𝑇 The set of time periods 𝐼 The set of demand points 𝐹 The set of Fulfilment centers 𝑃 The set of items Parameters 𝑑𝑝𝑖𝑡 Demand of demand zone 𝑖 for item 𝑝 in period 𝑡 𝛼𝑝𝑡 Total number of available items 𝑝 in period 𝑡 𝛽𝑓𝑖𝑝𝑡 The unit of shipment cost from the fulfilment centre f to demand zone 𝑖 in period 𝑡 ℎ𝑝𝑓𝑡 The unit of storage cost of item 𝑝 in fulfilment center 𝑓 in period 𝑡 𝜛𝑓𝑡 The capacity of fulfilment center 𝑓 in period 𝑡 𝑟𝑝𝑓𝑡 The processing cost of item 𝑝 in fulfilment center 𝑓 in period 𝑡 𝑣𝑝 The volume of item 𝑝 Decision Variables 𝑥𝑝𝑓𝑡 Quantity of item 𝑝 shipped to fulfilment center 𝑓 in period 𝑡 𝑦𝑓𝑖𝑝𝑡 Quantity of item 𝑝 shipped from fulfilment center 𝑓 to demand zone 𝑖 in period 𝑡 𝑞𝑝𝑓𝑡 The inventory level of fulfilment center 𝑓 for item 𝑝 in period 𝑡 We can formulate the online-retail inventory management problem (1) as follows: 10 𝑍𝑂𝑃𝑇 =𝑀𝑖𝑛𝑥, 𝑞, 𝑦 (∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗ 𝑞𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡) (1𝑎) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: ∑𝑦𝑓𝑖𝑝𝑡𝑓≥ 𝑑𝑝𝑖𝑡 ∀ 𝑝, 𝑖, 𝑡 (1𝑏) 𝑞𝑝𝑓𝑡 = 𝑞𝑝𝑓,𝑡−1 + 𝑥𝑝𝑓𝑡 −∑𝑦𝑓𝑖𝑝𝑡𝑖 ∀ 𝑝, 𝑓, 𝑡 (1𝑐) ∑𝑥𝑝𝑓𝑡𝑓≤ 𝛼𝑝𝑡 ∀ 𝑝, 𝑡 (1𝑑) ∑𝑣𝑝𝑞𝑝𝑓,𝑡−1𝑝+∑𝑣𝑝𝑥𝑝𝑓𝑡𝑝≤ 𝜛𝑓𝑡 ∀ 𝑓, 𝑡 (1𝑒) 𝑥𝑝𝑓𝑡, 𝑞𝑝𝑓𝑡, 𝑦𝑓𝑖𝑝𝑡 ∈ 𝑍+ (1𝑓) The objective function (1𝑎) aims to minimize the overall costs related to processing, storage, and transportation. The first constraint set (1b) guarantees that the total quantity shipped from all fulfilment centres must meet the demand of item 𝑝 for demand zone 𝑖 in each time period. The second equation (1c) captures the inventory balance for each fulfilment centre within each item and time period. It ensures that the inventory level at each centre is equal to the inventory level from the previous period plus the quantity assigned to the centre, subtracting the total quantity shipped from this centre to all demand zones. The third constraint set (1d) guarantees that the total quantity assigned to all fulfilment centres does not exceed the total supply. The constraint set (1e) indicates that the sum of the quantity assigned to fulfillment center 𝑓 and the inventory carried over from the previous period must be less than or equal to the storage capacity of fulfillment center 𝑓. The last constraint set (1f) defines the domain and type of decision variables. Among variables defined above, 𝑥𝑝𝑓𝑡 and 𝑦𝑓𝑖𝑝𝑡 are key variables, while 𝑞𝑝𝑓𝑡 is a secondary variable, as its value can be obtained from 𝑥𝑝𝑓𝑡 and 𝑦𝑓𝑖𝑝𝑡 through flow balance equation. 3.1. Complexity of problem Here, we interpret the original program as a network flow problem, incorporating additional 'non-network' type side constraints (1e). First, the constraint sets (1b) and (1d) represent capacity constraints, which can be expressed as 𝑙 ≤ 𝑥 ≤ 𝑢. To demonstrate that equation set (1c) is a 11 network conservation equation, shift all variables in (1c) to the left-hand side of the equations. Consequently, each variable in (1c) appears at most twice, once with a coefficient of + 1 and once with a coefficient of -1, which can be written as 𝐵𝑥 = 0, where 𝐵 is a total unimodular matrix. In addition, the constraint set (1e) is a multiple knapsack problem with weights 𝜆𝑗 and a set of knapsacks 𝑀 = {1,… ,𝑚} with positive capacities 𝜓𝑖 which can be written in the form of ∑ 𝜆𝑗𝑁𝑗=1 𝑥𝑖𝑗 ≤ 𝜓𝑖, 𝑖 = 1, … ,𝑚 or equivalently 𝜆𝑇𝑥 ≤ 𝜓. Therefore, the original problem can be presented as: min{𝑐𝑇𝑥: 𝐵𝑥 = 0 , 𝜆𝑇𝑥 ≤ Ψ, 𝑙 ≤ 𝑥 ≤ 𝑢 𝑎𝑛𝑑 𝑥 ∈ 𝑧𝑛} Where 𝑀 = 𝐹,𝑁 = 𝑃, 𝜆𝑇 = 𝑣, 𝜓 = 𝜛, 𝑙 = 𝛼, 𝑢 = 𝑑, 𝑐𝑇 = (𝑟 ℎ 𝛽), and 𝐵 = {−1,0,1}|𝐴|∗|𝐸|. Here, |𝐴| and |𝐸| denote the number of nodes and arcs in the corresponding network, respectively. As a result, the problem combines aspects of network flow and bounded multiple knapsack problems. Our objective is to demonstrate the NP-completeness of the aforementioned problem. The complexity of the problem depends on the complexity of the knapsack constraint. If we apply the Lagrangian relaxation method to this constraint, the reaming problem is transformed to a linear programming (LP) model which can be solved in polynomial time. The following paragraphs provide clarification on why we believe the problem is NP-complete. To demonstrate the dynamic online-retail inventory management Problem is NP-hard in the strong sense, we consider its decision version and prove that the decision version of the Dynamic Promotion Planning Optimization Problem is NP-complete. DEFINITION 1.3.1 (decision dynamic online-retail inventory management problem). Given a set of items 𝑝 ∈ 𝑃, warehouses 𝑓 ∈ 𝐹, demand zones 𝑖 ∈ 𝐼, planning period 𝑡 ∈ 𝑇, and parameters demand 𝑑𝑝𝑖𝑡 , supply 𝛼𝑝𝑡, transportation cost 𝛽𝑓𝑖𝑝𝑡, holding cost ℎ𝑝𝑓𝑡, receiving cost 𝑟𝑝𝑓𝑡, warehouse capacity 𝜛𝑓, and item volume 𝑣𝑝, does there exist a feasible solution to the formulation (1) with an objective value less than or equal Ω? We propose a polynomial time reduction from the bounded multiple knapsack problem to the decision dynamic online-retail inventory management problem. We demonstrate that for every instance of the bounded multiple knapsack problem (denoted by 𝐼𝑀𝐾𝑃) and its reduced decision 12 dynamic online-retail inventory management problem instance (denoted by 𝐼𝐷𝐼𝑀), the answer to 𝐼𝑀𝐾𝑃 is “yes” if, and only if, the answer to 𝐼𝐷𝐼𝑀 is “yes”. DEFINITION 2.3.1 (Bounded Multiple Knapsack Problem). Given a finite set 𝑀 of knapsacks, each with a capacity of 𝜓𝑙 for 𝑙 ∈ 𝑀, and a finite set 𝑁 of items, each with a corresponding weight 𝜆𝑗 and cost 𝜍𝑗 in which 𝑜𝑗 identical copies of item type 𝑗 must be used for 𝑗 ∈ 𝑁, does there exist a feasible solution to the following formulation with objective function value at most 𝛿? 𝐵𝑀𝐾𝑃 − 𝐷𝑒𝑐𝑖𝑠𝑖𝑜𝑛 ={ 𝑡ℎ𝑒𝑟𝑒 𝑒𝑥𝑖𝑠𝑡𝑠 𝑎𝑛 𝑥 𝑤𝑖𝑡ℎ ∑𝜍𝑗𝑥𝑗𝑙 ≤ 𝛿𝑙,𝑗 ∑𝜆𝑗𝑥𝑗𝑙 ≤ 𝜓𝑙,𝑛𝑗=1 𝑙 = 1,… ,𝑀∑𝑥𝑗𝑙𝑚𝑙=1≥ 𝑜𝑗 , 𝑗 = 1,… , 𝑁 𝑥𝑗𝑙 ∈ 𝑍+, 𝑙, 𝑗 = 1,… ,𝑁 } The bounded Multiple knapsack problem has been shown to be NP-complete (Kellerer et al., 2004). In what follows, we describe a reduction from the bounded multiple knapsack problem to the decision dynamic online-retail inventory management problem. Polynomial Reduction. Here, we take the instance of the 𝐼𝑀𝐾𝑃 and transform it to the instance of the 𝐼𝐷𝐼𝑀. Let each knapsack 𝑙 ∈ 𝑀 corresponds to a warehouse 𝑓 ∈ 𝐹, and each item 𝑗 ∈ 𝑁 in 𝐼𝑀𝐾𝑃 corresponds to a product 𝑝 ∈ 𝑃 in 𝐼𝐷𝐼𝑀. Let 𝑃 = 𝑁, 𝐹 = M, and |𝑇| = 1. For each item 𝑝 ∈ 𝑃, demand points 𝑖 ∈ 𝐼 and periods 𝑡 ∈ 𝑇, the parameters of objective function are defined as 𝑟𝑝𝑓𝑡 =𝜍𝑗 , ℎ𝑝𝑓𝑡 = 𝜍𝑗 , 𝛽𝑓𝑖𝑝 = 𝜍𝑗 . For each item 𝑝 ∈ 𝑃, the volume is 𝑣𝑝 = 𝜆𝑗, and demand is set to be 𝑑𝑝𝑖𝑡 =𝑜𝑗. The capacity of warehouse 𝑓 is 𝜛𝑓 = 𝜓𝑙 . Finally, the objective function of 𝐼𝐷𝐼𝑀 is equal to ∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡 + ∑ 𝜍𝑗 ∗ 𝑞𝑝𝑓𝑡𝑝,𝑓,𝑡 +∑ 𝜍𝑗 ∗ 𝑦𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 . According to this reduction |𝑇| = 1, so 𝑞𝑝𝑓𝑡 = 𝑥𝑝𝑓𝑡 − ∑ 𝑦𝑓𝑖𝑝𝑡𝑖 , so we have ∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡 + ∑ 𝜍𝑗𝑝,𝑓,𝑡 (𝑥𝑝𝑓𝑡 − ∑ 𝑦𝑓𝑖𝑝𝑡𝑖 ) + ∑ 𝜍𝑗 ∗ 𝑦𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 = 2∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡 . Thus Ω =2∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡 . We can also map the objective function 𝐼𝑀𝐾𝑃 to the objective function of the 13 original problem. Since ∑ 𝜍𝑗𝑥𝑗𝑙 ≤ 𝛿𝑙,𝑗 , the value of objective function 𝐼𝐷𝐼𝑀 will be at most 2𝛿, therefore Ω = 2𝛿. • Proof that “yes” to 𝑰𝑴𝑲𝑷 Implies “yes” to 𝑰𝑫𝑰𝑴 . Given a feasible solution for 𝐼𝑀𝐾𝑃 defined as follows: ∃ ?⃛?𝑗𝑙 ∈ 𝑍+ 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 ∑ ?⃛?𝑗𝑙𝑚𝑙=1 ≥ 𝑜𝑗 and ∑ 𝜆𝑗?⃛?𝑗𝑙 ≤ 𝜓𝑙𝑛𝑗=1 ∀ 𝑙 ∈ 𝑀, 𝑗 ∈ 𝑁, consider the following policy for 𝐼𝐷𝐼𝑀 such that Α(𝐼𝐷𝐼𝑀) ={ 𝑥𝑝𝑓𝑡 = 𝑥𝑗𝑙 𝑦𝑓𝑖𝑝𝑡 = 𝑚𝑖𝑛 {𝑥𝑗𝑙 ,𝑑𝑝𝑖𝑡}𝑞𝑝𝑓𝑡 = 𝑥𝑝𝑓𝑡 −∑𝑦𝑓𝑖𝑝𝑡𝑖 𝑖𝑓 ∃ 𝑥𝑗𝑙 ∈ 𝑍+ 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 ∑𝑥𝑗𝑙𝑚𝑙=1≥ 𝑜𝑗 𝑎𝑛𝑑 ∑𝜆𝑗𝑥𝑗𝑙 ≤ 𝜓𝑙𝑛𝑗=1 ∀ 𝑙 ∈ 𝑀, 𝑗 ∈ 𝑁 Where 𝑑𝑝𝑖𝑡 = 𝑜𝑗, 𝜛𝑓 = 𝜓𝑙, 𝑣𝑝 = 𝜆𝑗, and 𝑑𝑝𝑖𝑡 ≤ 𝛼𝑝𝑡. The above policy Α(𝐼𝐷𝐼𝑀) is feasible for 𝐼𝐷𝐼𝑀 as well. This is because if 𝑥𝑗𝑙 is a feasible solution to 𝐼𝑀𝐾𝑃, the capacity and demand fulfilment constraints (1b) and (1e) are immediately satisfied. Since we have 𝑑𝑝𝑖𝑡 ≤ 𝛼𝑝𝑡, the constraint (1d) is also satisfied. Thus, a feasible solution for 𝐼𝑀𝐾𝑃 guarantees a feasible policy for 𝐼𝐷𝐼𝑀. Moreover, since 𝑥𝑗𝑙 is feasible for 𝐼𝑀𝐾𝑃, it implies: ∑ 𝜍𝑗𝑥𝑗𝑙 ≤ 𝛿𝑙,𝑗 , Consider 𝑟𝑝𝑓𝑡 = ℎ𝑝𝑓𝑡 =𝛽𝑓𝑖𝑝 = 𝜍𝑗 ∀𝑝 ∈ 𝑃, 𝑖 ∈ 𝐼, 𝑓 ∈ 𝐹. Substituting these into the objective function (1-1), we have: ∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ 𝜍𝑗 ∗ (𝑥𝑝𝑓𝑡 −∑𝑦𝑓𝑖𝑝𝑡𝑖)𝑝,𝑓,𝑡+ ∑ 𝜍𝑗 ∗∑𝑦𝑓𝑖𝑝𝑡𝑖𝑝,𝑓,𝑡= 2∑ 𝜍𝑗 ∗ 𝑥𝑝𝑓𝑡𝑝,𝑓,𝑡= 2∑𝜍𝑗 ∗ 𝑥𝑗𝑙𝑗,𝑙≤ 2𝛿 Where 2𝛿 = Ω. Hence, given that there is a feasible solution for 𝐼𝑀𝐾𝑃, there exists an optimal policy to 𝐼𝐷𝐼𝑀 with the cost at most Ω. □ • Proof that “yes” to 𝑰𝑫𝑰𝑴 Implies “yes” to 𝑰𝑴𝑲𝑷. A feasible solution for 𝐼𝐷𝐼𝑀 implies all the constraints of the formulation (1) are satisfied. Given a feasible solution for 𝐼𝐷𝐼𝑀, consider the following policy for 𝐼𝑀𝐾𝑃 such that Α′(𝐼𝑀𝐾𝑃) = {𝑥𝑗𝑙 = ?̃?𝑝𝑓𝑡 𝑖𝑓 ∃ 𝑓𝑒𝑎𝑠𝑖𝑏𝑙𝑒 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑓𝑜𝑟 𝐼𝐷𝐼𝑀 14 Having |𝑇| = 1, the policy Α′(𝐼𝑀𝐾𝑃) will be also feasible for 𝐼𝑀𝐾𝑃 as well. When ?̃?𝑝𝑓𝑡 is feasible for 𝐼𝐷𝐼𝑀, it means the constraints ∑ 𝜆𝑗𝑥𝑗𝑙 ≤ 𝜓𝑙 𝑛𝑗=1 and ∑ 𝑥𝑗𝑙𝑚𝑙=1 ≥ 𝑜𝑗 are met automatically. Additionally, since the feasible solution to 𝐼𝐷𝐼𝑀 has a cost at most Ω, Ω ≥ ∑ 𝜍𝑗 ∗ ?̃?𝑝𝑓𝑡𝑝,𝑓,𝑡 +∑ 𝜍𝑗 ∗ ?̃?𝑝𝑓𝑡𝑝,𝑓,𝑡 + ∑ 𝜍𝑗 ∗ ?̃?𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 . In the optimal solution ?̃?𝑝𝑓𝑡 = ∑ ?̃?𝑓𝑖𝑝𝑡𝑖 , then ?̃?𝑝𝑓𝑡 = (?̃?𝑝𝑓𝑡 −∑ ?̃?𝑓𝑖𝑝𝑡𝑖 ) , so we have Ω ≥ ∑ 𝑟𝑝𝑓𝑡 ∗ ?̃?𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗ (?̃?𝑝𝑓𝑡 −∑ ?̃?𝑓𝑖𝑝𝑡𝑖)𝑝,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ ?̃?𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 = ∑(𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ ?̃?𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ ?̃?𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 = ∑(𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗∑ ?̃?𝑓𝑖𝑝𝑡𝑖𝑝,𝑓,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ ?̃?𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡 = ∑ (𝑟𝑝𝑓𝑡 + 𝛽𝑓𝑖𝑝𝑡) ∗ ?̃?𝑓𝑖𝑝𝑡𝑝,𝑓,𝑖,𝑡 Considering 𝑟𝑝𝑓𝑡 = 𝛽𝑓𝑖𝑝𝑡 = 𝜍𝑗, we have Ω ≥ 2∑ 𝜍𝑗 ∗ ?̃?𝑓𝑖𝑝𝑡𝑝,𝑓,𝑖,𝑡 which leads to ∑ 𝜍𝑗 ∗𝑝,𝑓,𝑡∑ ?̃?𝑓𝑖𝑝𝑡𝑖 = ∑ 𝜍𝑗 ∗ ?̃?𝑝𝑓𝑡𝑝,𝑓,𝑡 ≤Ω2. According to the reduction described earlier Ω = 2𝛿. As a result, the constraint ∑ 𝜍𝑗𝑥𝑗𝑙 ≤ 𝛿𝑙,𝑗 will be satisfied. Therefore, given that there is an optimal policy to 𝐼𝐷𝐼𝑀 with cost at most Ω, there exists a feasible solution for 𝐼𝑀𝐾𝑃. □ Considering the polynomial time reduction from the bounded multiple knapsack problem to the Decision dynamic online-retail inventory management problem and given that every instance of the bounded multiple knapsack problem, 𝐼𝑀𝐾𝑃, can be reduced to the Decision dynamic online-retail inventory management problem instance 𝐼𝐷𝐼𝑀, BMKP-decision has a feasible solution if and only if the corresponding decision dynamic online-retail inventory management problem has a feasible solution. In other words, the answer to 𝐼𝑀𝐾𝑃 is “yes” if and only if, the answer to 𝐼𝐷𝐼𝑀 is “yes”. Therefore, it can be concluded that the Decision dynamic online-retail inventory management problem is equally hard as the problem of bounded multiple knapsack problem. Hence, the dynamic online-retail inventory management Problem is NP-complete. 15 4. Approximation Algorithm Given that the optimal model (1) is NP-hard in the strong sense, finding an exact solution is computationally intractable even for small instances. As a result, heuristic techniques such as the Greedy algorithm offer a practical approach to achieve near-optimal solutions within a reasonable time. The Greedy algorithm is an optimization technique that makes the best local choice at each stage, but it does not guarantee a globally optimal solution. Although it usually provides a sub-optimal solution, its main advantage lies in its ability to solve problems much faster, often within polynomial time. Below, we describe our suggested Greedy algorithm for solving the optimal model (1). 4.1. Greedy Algorithm The greedy algorithm for solving model (1) begins by forming a group of bundles 𝜑 = (𝑝, 𝑖) for each period where 𝑝 and 𝑖 stand for the items and demand points respectively. It then prioritizes these bundles based on the maximum value potential function 𝜌. In the subsequent section, we will discuss the potential function and how it is designed. The algorithm is designed to determine the flow of each warehouse 𝑓 for each bundle 𝜑, denoted as ?̃?𝑓𝜑𝑡 , from which the values of other variables ?̃?𝑓𝑝𝑡, ?̃?𝑓𝑝𝑡 are derived. ?̃?𝑓𝑝𝑡, ?̃?𝑓𝑝𝑡 are the decision variables in the model (1) which are defined as the value of initial stocking, and the inventory level of warehouse 𝑓, respectively. Here, the tilde (~) is used for each variable to indicate that the solution is obtained through heuristic, distinguishing it from the optimal solution. Additionally, the algorithm employs a backward strategy to address the dynamic nature of model. Hence, starting from the last period, the bundles are arranged in decreasing order of potential function 𝜌, the highest priority bundle is picked, and among different warehouses, the one with minimum cost, denoted as 𝑓∗, is selected for this particular bundle 𝜑. Subsequently, the value of ?̃?𝑓𝜑𝑡 is determined by the minimum of 𝑑𝜑𝑡 and the current capacity of warehouse 𝑓∗, i.e., ?̃?𝑓∗𝜑𝑡 = min {𝑑𝜑𝑡, (𝜛𝑓∗𝑡 − ?̅?𝑓∗𝑡) 𝑣𝑝⁄ } where ?̅?𝑓𝑡 is the auxiliary variable representing the utilized capacity of 𝑓∗ . If ?̃?𝑓∗𝜑𝑡 = 𝑑𝜑𝑡, warehouse 𝑓∗ can fully satisfy 𝑑𝜑𝑡, otherwise, the model proceeds to the next cheapest warehouses until 𝑑𝜑𝑡 is fully satisfied. Once ?̃?𝑓∗𝜑𝑡 has been established for every bundle 𝜑, the algorithm now seeks to find out the values of ?̃?𝑓𝑝𝑡 and ?̃?𝑓𝑝𝑡. Up until to this point, it has been identified which warehouse is responsible for 16 fulfilling the bundle 𝜑 = (𝑝, 𝑖). To determine how to stock the warehouses 𝑓 with each item 𝑝, i.e., the value of ?̃?𝑓𝑝𝑡, we define new bundle 𝜑1 = (𝑓∗, 𝑝). This approach is adopted because variables of ?̃?𝑓𝑝𝑡 and ?̃?𝑓𝑝𝑡 depend on the warehouse 𝑓 and the item 𝑝, but not on the demand point 𝑖. By introducing 𝜑1, we simplify the analysis as it allows us to focus on the relationships between warehouses and items without considering individual demand points directly. The bundles 𝜑1 are arranged in decreasing order based on the total value of ?̃?𝜑1𝑖𝑡 for all demand points 𝑖, and inventory level of each warehouse for the current period, i.e., ?̃?𝜑1𝑡 . Given that the algorithm starts from the last period and ?̃?𝜑1𝑡 represents the inventory level at the end of the period, the value of ?̃?𝜑1𝑡 is initially zero, meaning that ?̃?𝜑1𝑇 = 0. The value of ?̃?𝜑1𝑡 is determined based on the minimum of (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡) and total current supply (𝛼𝑝𝑡 − ?̅?𝑝𝑡) for item p, i.e., ?̃?𝜑1𝑡 = 𝑚𝑖𝑛 {(∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡), 𝛼𝑝𝑡 − ?̅?𝑝𝑡}; where 𝛼𝑝𝑡 stands for the supply of item 𝑝, and ?̅?𝑝𝑡 is the utilized portion of 𝛼𝑝𝑡. If the supply is insufficient such that ?̃?𝜑1𝑡 ≠ (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡), the remaining ?̃?𝜑1𝑡 will be fulfilled by inventory from previous period, ?̃?𝜑1,𝑡−1. This ensures that the flow conservation equation ?̃?𝜑1,𝑡−1 = ∑ ?̃?𝜑1𝑖𝑡𝑖 − ?̃?𝜑1𝑡 +?̃?𝜑1𝑡 holds, maintaining the balance between total inflow and total outflow for each bundle 𝜑1. Once these three variables are established, the model proceeds to the next iteration, period 𝑡 − 1, and keeps repeating the same procedure until the variables are determined for all periods. The summary of algorithm is as follows: Table 2. The summary of greedy algorithm for solving optimal model Step 1 Define 𝑁 as the total number of (𝑃, 𝐼) combinations, set 𝑡 = 𝑇, ?̃?𝜑1𝑇 = 0 Step 2 Set ?̅?𝑓𝑡 = 0, ?̅?𝑝𝑡 = 0 Step 3 Define 𝜑 = (p, 𝑖), sort them in decreasing order of potential function 𝜌 Step 4 Set 𝜑 = 1 (represents the bundle (p, 𝑖) with the highest potential function) Step 5 Choose a warehouse with the least cost where 𝑓∗ = 𝐴𝑟𝑔𝑚𝑖𝑛{ (𝑟𝑝𝑓𝑡 + h𝑝𝑓𝑡 +β𝑓𝑖𝑝𝑡)} If ?̅?𝑓∗𝑡 + 𝑣𝑝 ≤ 𝜛𝑓∗𝑡: Assign bundle 𝜑 to the warehouse 𝑓∗ 17 ?̃?𝑓∗𝜑𝑡 = min {𝑑𝜑𝑡, (𝜛𝑓∗𝑡 − ?̅?𝑓∗𝑡) 𝑣𝑝⁄ } ?̅?𝑓∗𝑡 = ?̅?𝑓∗𝑡 + 𝑣𝑝 ∗ ?̃?𝑓∗𝜑𝑡 𝑑𝜑𝑡 = 𝑑𝜑𝑡 − ?̃?𝑓∗𝜑𝑡 Else move to Step 6 Step 6 While 𝑑𝜑𝑡 > 0, then 𝑓∗ = 𝑓∗ + 1, move to step 5 Step 7 Set 𝜑 = 𝜑 + 1, move to step 5 Step 8 Define 𝜑1 = (𝑓∗, p), sort them in decreasing order of (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡) Step 9 If ?̅?𝑝𝑡 ≤ 𝛼𝑝𝑡: ?̃?𝜑1𝑡 = 𝑚𝑖𝑛 {(∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡), 𝛼𝑝𝑡 − ?̅?𝑝𝑡} ?̃?𝜑1,𝑡−1 = ∑ ?̃?𝜑1𝑖𝑡𝑖 − ?̃?𝜑1𝑡 + ?̃?𝜑1𝑡 ?̅?𝑝𝑡 = ?̅?𝑝𝑡 + ?̃?𝜑1𝑡 Else ?̃?𝜑1𝑡 = 0, ?̃?𝜑1,𝑡−1 = ∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡 Step 10 Set 𝜑1 = 𝜑1 + 1, move to step 9 Step 11 ?̅?𝑓,𝑡−1 = ?̅?𝑓,𝑡−1 + 𝑣𝑝 ∗ ?̃?𝜑1,𝑡−1 Step 12 Set 𝑡 = 𝑡 − 1, move to step 2 In the next section, we describe the performance of the above algorithm. 4.2. Potential function Before delving into the performance of the greedy algorithm, let’s discuss the selection criterion used for bundles, referred to as the potential function. This criterion is defined as the weighted average cost per unit volume for a bundle 𝜑, denoted by 𝜌 =(∑ 𝜗𝑝𝑓𝑓 )∗d𝜑𝑣𝑝, where 𝜑 represents the bundle (𝑝, 𝑖), 𝑑𝜑 is the demand for bundle 𝜑, and 𝜗𝑝𝑓𝑡 is the cost of assigning bundle (𝑝, 𝑖) to warehouse 𝑓 ∈ 𝐹 where 𝜗𝑝𝑓𝑡 is calculated as 𝜗𝑝𝑓𝑡 = 𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡 + 𝛽𝑓𝑖𝑝𝑡. In the following discussion, we aim to justify the rationale behind using this criterion. Let’s first consider the potential function excluding ∑ 𝜗𝑝𝑓𝑓 , defined as the ratio of demand to volume, denoted by 𝜌 =d𝜑𝑣𝑝. When employing the greedy algorithm, bundles are arranged in descending order of this efficiency measure. The algorithm selects the bundle with the highest efficiency and assigning it to the most cost-effective warehouses to meet the demand for that 18 bundle, while taking capacity restrictions into account. Utilizing this approach can yield a feasible solution, but in certain special cases, the outcomes may be unsatisfactory. For instance, consider a scenario with two bundles, each with volumes 𝑣1 = 1, 𝑣2 = 2 and equal demands of 1. These bundles need to be allocated to two warehouses with capacities 𝜛1 = 𝜛2 = 2. The associated matrix 𝜗𝑝𝑓𝑡 for a specific period 𝑡 is as follows: 𝜗𝑝𝑓𝑡 = [2 41 𝑀′] Comparing the efficiencies of the bundles, the Greedy algorithm assigns bundle 1 to warehouse 1 and bundle 2 to warehouse 2 with a value of 2+ 𝑀′, while the optimal solution is 5. For a suitably large selection of 𝑀′, any relative performance guarantee for greedy can be forced to be arbitrarily close to infinity. As a result, to avoid this pathological special case, we must modify our sorting criteria. The aforementioned issue can be avoided if we swap out d𝜑𝑣𝑝 with (∑ 𝜗𝑝𝑓𝑓 )∗d𝜑𝑣𝑝, i.e., sort the bundles according to the weighted average cost per unit volume rather than simply the demand. Thus, if we repeat the algorithm for the above instance with the new potential function, the solutions provided by Greedy and the optimal solution would be the same, and the gap would be zero as 𝑀′ approaches a suitably large number. 4.3. Greedy performance The benefit of using the greedy algorithm is twofold: the approach is straightforward to implement. In fact, the algorithm follows a clear set of steps, enabling it to handle an NP-hard dynamic programming problem through an iterative structure. Additionally, the greedy generally has lower computational complexity making it faster and more suitable for large-scale problems where speed is critical. However, this is not the case for other optimal algorithms. For instance, the branch and bound algorithm requires to do branching on non-integer variables and solving a linear programming model at each iteration. In contrast, the greedy algorithm makes a series of locally optimal choices at each step (e.g., selecting the warehouse with the least cost) without backtracking. Thus, the algorithm can quickly provide a feasible solution to the optimal model. However, there are some weaknesses in the algorithm. The solution obtained by the greedy method is sub-optimal. As illustrated in Section 4.2, the algorithm's effectiveness heavily relies on the potential function 𝜌. This function aims to minimize the worst-case scenario by focusing on 19 bundles with the highest cost. Consequently, the algorithm prioritizes assigning these costly bundles, which can lead to suboptimal outcomes. This assignment strategy might result in bundles with lower potential functions being assigned to more expensive warehouses, even though they could have been assigned to cheaper ones. As a result, this can increase the overall cost, leading to a higher optimality gap and a less efficient solution. Consider the following instance for a specific one period. Let 𝜑 = (𝑝, 𝑖) and 𝐹 denote the respective sets of bundles and warehouses. We define 𝑃ℎ, and 𝑃𝑙 as sets of high-potential and low-potential bundles where the unit assignment cost is defined as follows: 𝜗𝜑𝑓 = {1 𝑖𝑓 𝑓 = 1 𝑎𝑛𝑑 𝜑 ∈ 𝑃𝑙𝑀′ ∀𝑓 ∈ 𝐹\{1} 𝑎𝑛𝑑 𝜑 ∈ 𝑃𝑙𝑀′ ∀𝑓 ∈ 𝐹 𝑎𝑛𝑑 𝜑 ∈ 𝑃ℎ The bundles in 𝑃ℎ have higher potential function than those in 𝑃𝑙, i.e., 𝜌ℎ > 𝜌𝑙. The total number of bundles are 𝑁 which includes 𝑁ℎ for the number of 𝑃ℎ, and 𝑁𝑙 for the number of 𝑃𝑙 where 𝑁 =𝑁ℎ + 𝑁𝑙. Bundles in 𝑃ℎ and 𝑃𝑙 might be different in volume and demand; however, the following relationship holds: ∑ 𝑣𝑝 ∗ 𝑑𝜑𝜑∈𝑃𝑙≤ ∑ 𝑣𝑝 ∗ 𝑑𝜑𝜑∈𝑃ℎ There are 𝑀 number of warehouses where the capacity of the first one is equal to 𝜛1 =∑ 𝑣𝑝 ∗ 𝑑𝜑𝜑∈𝑃𝑙 . The greedy algorithm prioritizes bundles from 𝑃ℎ over those from 𝑃𝑙 due to their higher potential functions. As a result, the algorithm assigns the first bundle from 𝑃ℎ to the minimum cost warehouse and continues to assign bundles from 𝑃ℎ to this warehouse until its capacity is exhausted. Given the assumption that the capacity of the first warehouse (𝜛1 = ∑ 𝑣𝑝 ∗ 𝑑𝜑𝜑∈𝑃𝑙 ) is insufficient to accommodate all bundles from 𝑃ℎ, remaining bundles from 𝑃ℎ and all bundles from 𝑃𝑙 are assigned to subsequent warehouses with unit cost 𝑀′. Consequently, the cost of the greedy approach is 𝑍(𝑆𝐺) = 𝑀′∑ 𝑑𝜑𝜑∈𝑃ℎ +𝑀′ ∑ 𝑑𝜑𝜑∈𝑃𝑙 . In contrast, the optimal model assigns all bundles from 𝑃𝑙 to the first warehouse with unit cost 1, which can accommodate them, and assigns bundles from 𝑃ℎ to the remaining warehouses. This results in a cost 𝑍(𝑆∗) = ∑ 𝑑𝜑𝜑∈𝑃𝑙 +𝑀′∑ 𝑑𝜑𝜑∈𝑃ℎ . Therefore, the relative performance of the greedy for this specific example will be 20 greater then 1, illustrating how the algorithm's prioritization of high-potential bundles can lead to suboptimal solutions. Moreover, although the greedy algorithm performs well in a single period, it does not necessarily extend its effectiveness to a multi-period programming model. The algorithm addresses the dynamic nature of the problem by breaking it down into ∣T∣ individual iterations. However, the greedy approach does not account for the interactions between periods when solving a multi-period problem. Each period is treated independently without considering the effect of previous periods, whereas optimal algorithms consider all periods simultaneously to determine the values of decision variables. This lack of integration across periods can limit the effectiveness of the greedy approach for multi-period planning. The algorithm’s strategy to determine the inventory level of warehouses is suboptimal. In step 9, after sorting warehouses in decreasing order of (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡), the algorithm selects the first bundle 𝜑1 = (𝑓∗, 𝑝). The value of ?̃?𝜑1𝑡 is then determined based on the minimum of (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡) and total current supply (𝛼𝑝𝑡 − ?̅?𝑝𝑡) for each item 𝑝. This approach results in an unbalanced inventory distribution across warehouses. Initially selected warehouses end up with ?̃?𝜑1,𝑡−1 = 0, while those chosen later will accumulate a large value ?̃?𝜑1,𝑡−1 = (∑ ?̃?𝜑1𝑖𝑡𝑖 + ?̃?𝜑1𝑡). This imbalance can lead to increased storage costs for some warehouses, thereby raising the overall cost of the greedy algorithm in certain cases. Although the greedy is an efficient method to cope with the NP-hard nature of model (1), it does not resolve the issue of computational intractability for large-scale problems. In fact, the computation time of greedy is 𝑂(𝑇𝑁𝐹) where 𝑇 is the number of periods, 𝑁 is the number of bundles (𝑝, 𝑖) and 𝐹 is the number of warehouses. The algorithm requires enumerating all possible combinations of (𝑝, 𝑖) throughout each iteration. Given that the optimal model (1) is large scale, this results in an extremely large number of bundles 𝜑. Therefore, we must develop a tractable algorithm which creates a trade-off between the running time and approximation ratio. In Section 5. we describe our new strategy to address this issue. 21 5. Aggregation algorithm The greedy algorithm described in Section 4.1 does not completely resolve the issue of computational intractability for large-scale problems. With a computation time of 𝑂(𝑇𝑁𝐹), the algorithm necessitates enumerating all combinations of bundles (𝑝, 𝑖) in each iteration, leading to a significant number of bundles 𝑁 in a large-scale model. Therefore, a more tractable algorithm is needed to create a trade-off between the running time and the approximation ratio. Here, we propose a new methodology focused on an aggregation approach, which utilizes the greedy algorithm as a subroutine. Using this methodology, we can further reduce the computation time required to solve a large-scale model, albeit at the expense of a larger optimality gap. This section will explain the aggregation approach employed in this work. As mentioned previously, we conduct aggregation on three sets: demand zones, fulfillment centers, and items. The rationale behind this decision comes from the original model's nature as a NP-hard problem, wherein computation time exponentially increases with problem dimensions. While the greedy algorithm can somewhat reduce computation time, it is insufficient for large-scale instances due to its dependence on the number of bundles. As the problem size increases, the number of bundles grows significantly, making the greedy approach impractical for handling large instances. Therefore, a more tractable algorithm is necessary to efficiently manage these large-scale problems. Reducing the model's dimensionality is an effective strategy to make the model tractable from both computational and theoretical perspectives. However, it's essential to acknowledge that aggregation and disaggregation method results in errors into the decision-making process and consequently suboptimality of optimization model. We will delve into this issue in Section 5.4. In this section, our focus is on elaborating the aggregation mechanism. Handling a substantial number of items complicates the management of product for inventory control, particularly as transaction volumes increase. This challenge can be addressed through a more effective approach, employing an aggregation method to categorize all products into manageable clusters based on their similarities. The values of these various groups can then be estimated and assessed by appropriate techniques. Therefore, the aggregation approach will facilitate the resource management decision-making process by reducing the dimensionality of the 22 problem (Gustriansyah et al., 2020). Here, we perform aggregation on the three sets of items, demand points, and warehouses, based on the historical data from all demand points for various items. These aggregated groups remain constant throughout the planning horizon. Let the items 𝑃 = {1,… , 𝑙} be a set of 𝑙 items. Suppose 𝑃 is clustered into 𝐶 groups, where the set of indices in the 𝑐-th group is denoted 𝑃𝑐 such that ⋃ 𝑃𝑐𝐶𝑐=1 = {1,… , 𝑙} and 𝑃𝑐⋂𝑃𝑐′ = ∅ ∀𝑐 and 𝑐′, 𝑐 ≠ 𝑐′. We also have clustering on demand points and fulfilment points. Let 𝐼 = {1,… , 𝑛} be a set of 𝑛 demand points, and 𝐹 = {1,… ,𝑚} be a set of 𝑚 fulfilment points. Suppose 𝐼 and 𝐹 are grouped into set of clusters 𝐼𝑘, 𝑘 = 1, . . . , 𝐾 and 𝐹𝑤 , 𝑤 = 1, … ,𝑊, such that 𝐼𝑘⋂𝐼𝑘′ = ∅ ∀𝑘, ⋃ 𝐼𝑘𝐾𝑘=1 = {1,… , 𝑛}, 𝐹𝑤⋂𝐹𝑤′ = ∅ ∀𝑤 and ⋃ 𝐹𝑤𝑊𝑤=1 = {1,… ,𝑚}. We employ a well-known clustering algorithm, k-means, to efficiently perform aggregation. K-means is a centroid-based algorithm which aims to aggregate data into a predetermined number of groups while minimizing the sum of distances between the points and the centroid of their particular cluster. The clustering procedure is performed independently for each set of items, warehouses, and demand points. This process involves the consideration of suitable clustering features, such as the geographical locations of warehouses and demand sites, as well as item volume, cost, and demand, among others. Having completed the clustering process, we will now proceed to determine the parameters of the aggregation model. The following explains how to find out the aggregated model's parameters. To determine the model parameters for the reduced problem, we need to use a method which determines weights that reflect the relative importance of the parameters in the original problem and using these weights to derive parameters for the reduced problem. To do so, we introduce weight vectors 𝑔𝑤𝑐 = {𝑔𝑟𝑓𝑝𝑤𝑐} , 𝑐 = 1, … , 𝐶, 𝑤 = 1,… ,𝑊, such that ∑ ∑ 𝑔𝑟𝑓𝑝𝑤𝑐𝑝∈𝑃𝑐𝑓∈𝐹𝑤 = 1 for 𝑐 = 1,… , 𝐶 and 𝑤 = 1,… ,𝑊 to convert 𝑟𝑝𝑓𝑡 to the centroid aggregation 𝑟𝑐𝑤𝑡. The nonnegative weight vectors 𝑔ℎ𝑓𝑝𝑤𝑐, 𝑔𝛽𝑓𝑖𝑝𝑤𝑘𝑐 are also defined to transfer parameters ℎ𝑝𝑓𝑡 and 𝛽𝑓𝑖𝑝𝑡 to their aggregation format in the reduced model, respectively. These weight vectors must satisfy the normalization condition: ∑ ∑ 𝑔ℎ𝑓𝑝𝑤𝑐𝑝∈𝑃𝑐𝑓∈𝐹𝑤 = 1 for 𝑤 = 1,… ,𝑊, 𝑐 = 1,… , 𝐶 ∑ ∑ 𝑔𝛽𝑓𝑖𝑝𝑤𝑘𝑐𝑓∈𝐹𝑤𝑖∈𝐼𝑘 = 1 for 𝑤 = 1,… ,𝑊, 𝑘 = 1,… , 𝐾, 𝑐 = 1,… , 𝐶 23 One way to compute these vectors can be a fixed weight approach which involves a convex weighting of the elements inside a cluster, and if elements are homogenous, we can replace it with equal weighting. For instance, each element of vector 𝑔𝑤𝑐 can be calculated as 1|𝑃𝑐||𝐹𝑤|. This type of aggregation is similar to fixed-weight combination approach can be seen as an example of column aggregation, when three sets of items, demand and fulfilment points are aggregated simultaneously, we have a combination of both columns (variables) and rows (constraints) aggregation. 5.1. Aggregation formulation The set of clusters 𝐾 = {𝐼𝑘|𝑘 = 1,… , K}, 𝑊 = {𝐹𝑤|𝑤 = 1, … ,W}, 𝐶 = {𝑃𝑐|𝑐 = 1,… , C}, and the collection of weight vectors completely define an aggregated problem (6). In order to obtain a concise formulation of the aggregate problem, we define: ?̅?𝑐𝑘𝑡 = ∑ 𝑑𝑝𝑖𝑡𝑝∈𝑃𝑐,𝑖∈𝐼𝑘 , ?̅?𝑐𝑡 = ∑ 𝛼𝑝𝑡𝑝∈𝑃𝑐 , ?̅?𝑤𝑡 = ∑ 𝜛𝑓𝑡𝑓∈𝐹𝑤 , ?̅?𝑐 = ∑ 𝑔𝑝𝑐𝑣𝑝𝑝∈𝑃𝑐 , ?̅?𝑐𝑤𝑡 =∑ ∑ 𝑔𝑟𝑓𝑝𝑤𝑐𝑟𝑝𝑓𝑡𝑝∈𝑃𝑐𝑓∈𝐹𝑤 , ℎ̅𝑐𝑤𝑡 = ∑ ∑ 𝑔ℎ𝑓𝑝𝑤𝑐ℎ𝑝𝑓𝑡𝑝∈𝑃𝑐𝑓∈𝐹𝑤 , ?̅?𝑤𝑘𝑐 = ∑ ∑ 𝑔𝛽𝑓𝑖𝑝𝑤𝑘𝑐𝛽𝑓𝑖𝑝 𝑖∈𝐼𝑘𝑓∈𝐹𝑤 The following is then an aggregated version of original problem (1): ?̅?𝐴𝐺𝐺 =𝑚𝑖𝑛?̅?, ?̅?, ?̅? (∑ ?̅?𝑐𝑤𝑡?̅?𝑐𝑤𝑡𝑐,𝑤,𝑡+ ∑ ℎ̅𝑐𝑤𝑡 ?̅?𝑐𝑤𝑡𝑐,𝑤,𝑡+ ∑ ?̅?𝑤𝑘𝑐𝑡?̅?𝑤𝑘𝑐𝑡𝑐,𝑤,𝑘,𝑡) (2𝑎) 𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: ∑?̅?𝑤𝑘𝑐𝑡𝑤≥ ?̅?𝑐𝑘𝑡 ∀ 𝑐, 𝑘, 𝑡 (6𝑏) ?̅?𝑐𝑤𝑡 = ?̅?𝑐𝑤,𝑡−1 + ?̅?𝑐𝑤𝑡 −∑?̅?𝑤𝑘𝑐𝑡𝑘 ∀ 𝑐, 𝑤, 𝑡 (6𝑐) ∑?̅?𝑐𝑤𝑡𝑤≤ ?̅?𝑐𝑡 ∀𝑐, 𝑡 (6𝑑) ∑?̅?𝑐 ?̅?𝑐𝑤,𝑡−1𝑐+∑?̅?𝑐?̅?𝑐𝑤𝑡𝑐≤ ?̅?𝑤𝑡 ∀ 𝑤, 𝑡 (6𝑒) ?̅?𝑐𝑤𝑡, ?̅?𝑐𝑤𝑡, ?̅?𝑤𝑘𝑐𝑡 ≥ 0, 𝐼𝑛𝑡𝑒𝑔𝑒𝑟 (6𝑓) 24 This problem represents a reduced form of the original problem (1), whereby the original sets 𝐼, 𝐹 and 𝑃 are reduced to 𝐾,𝑊 and C, respectively. Therefore, in the worst-case scenario, the time required to solve the aggregated problem grows exponentially with respect to the size of the clusters rather than the original input size. The resulting policy from the aggregated model is denoted by ?̅?∗𝑐𝑤𝑡, ?̅?∗𝑐𝑤𝑡, and ?̅?∗𝑤𝑘𝑐𝑡. 5.2. Disaggregation steps and formulations In this section, we detail the process of converting the aggregated solution of problem (6) into a feasible solution for the original problem (1). This process involves breaking down the aggregated decisions into more detailed, and individual decisions, corresponding to the specific variables in the original problem. The simplest way to obtain a feasible solution to the original problem (1) from an aggregated solution (6) is to use the fixed-weight disaggregation. Denoting ?̂?𝑓𝑝𝑡, ?̂?𝑓𝑝𝑡, and ?̂?𝑓𝑝𝑖𝑡 as the solutions directly derived from the aggregated problem (6), under the fixed-weight disaggregation approach, they would be represented as follows: ?̂?𝑓𝑝𝑡 = 𝑔𝑟𝑓𝑝𝑤𝑐?̅?𝑐𝑤𝑡, ?̂?𝑓𝑝𝑡 = 𝑔ℎ𝑓𝑝𝑤𝑐?̅?𝑐𝑤𝑡, and ?̂?𝑓𝑝𝑖𝑡 =𝑔𝛽𝑓𝑖𝑝𝑤𝑘𝑐?̅?𝑤𝑘𝑐𝑡, respectively. However, in the integer programming case, fixed weight disaggregation does not in general result in an integer solution to the original problem. Thus, in our case, problem-specific procedures can be applied which ensures feasibility. In the following, we have defined a three-step hierarchal disaggregation approach as follows: • Demand zones disaggregation The goal of demand zones disaggregation stage is to determine the flow between each aggregated warehouse 𝑤 and each demand point 𝑖 within zone 𝑘, while considering the optimal solution of the aggregated model, represented as (?̅?𝑐𝑤𝑡∗ , ?̅?𝑐𝑤𝑡∗ , ?̅?𝑤𝑘𝑐𝑡∗ ). Since the clustering process ensures no intersection between points of different clusters, each demand point is assigned to only one zone. Additionally, the flow between the group of warehouses 𝑤 and each zone 𝑘 is predetermined; as obtained by solving the aggregated model (6) which is denoted as ?̅?𝑤𝑘𝑐𝑡∗ . As a result, we can perform the process of disaggregation for each zone 𝑘 independently without impacting the final solution. • Warehouses disaggregation 25 Once the flow between each group of warehouse 𝑤 and each demand zone 𝑖 is specified, we proceed to disaggregate the warehouses. To do so, we require the solutions of both the aggregation model (6) and the previous stage, i.e., (?̅?𝑐𝑤𝑡∗ , ?̅?𝑐𝑤𝑡∗ , 𝑦𝑤𝑖𝑐𝑡(1)∗). The amount of inflow, outflow, and inventory level of each fulfillment center 𝑓 within the aggregated warehouse 𝑤 are the decision variables must be determined in this stage. Similar to the previous logic, the solution space of this stage is divided into distinct subspaces, each representing a specific cluster 𝑤. Therefore, this allows us to perform disaggregation for each group 𝑤 individually. • Items disaggregation The decision variables for the original problem have been almost estimated up until this point, albeit with values obtained for the aggregated items. Now, our task is to specify the inventory level, inflow to fulfillment center 𝑓, and the amount of shipment between fulfillment center 𝑓 and demand zone 𝑖 for each specific item 𝑝. Thus, the final stage of disaggregation is conducted in a similar manner, with the input parameters being the solutions obtained from the previous step, represented as ( 𝑥𝑐𝑓𝑡(2)∗, 𝑞𝑐𝑓𝑡(2)∗ , 𝑦𝑓𝑖𝑐𝑡(2)∗). Similar to the preceding stages, this stage can be disaggregated individually, as each item is assigned to a single cluster and the solution of each cluster 𝑐 is predetermined from the previous stage. 26 Figure 2. The detailed workflow of aggregation and disaggregation stages Now, we will delve into the details of the optimization models corresponding to each of the previously described stages. We will provide a comprehensive explanation of the optimization models that utilized to disaggregate demand zones, warehouses, and items, respectively. 5.2.1. Optimal Disaggregation of demand zones Considering the optimal solution of the aggregated model, represented as (?̅?𝑐𝑤𝑡∗ , ?̅?𝑐𝑤𝑡∗ , ?̅?𝑤𝑘𝑐𝑡∗ ), the optimal disaggregation of ?̅?𝑤𝑘𝑐𝑡∗ can be achieved by solving subproblems for each 𝑘 = 1,… , 𝐾. Each subproblem k includes only the portion of the original problem that coincides with the members of each 𝐼𝑘. The objective function value for disaggregation of all demand zones will be equal to 𝑍𝐷𝑖𝑠𝑠𝐷 = ∑ 𝑍𝐷𝑖𝑠𝑠𝐷𝑘𝑘∈𝐾 where 𝑍𝐷𝑖𝑠𝑠𝐷𝑘 stands for the objective value of disaggregating demand zone 𝑘. To determine the disaggregation of demand zones, we define the decision variable 𝑦𝑤𝑖𝑐𝑡(1) for each subproblem 𝑘, representing the quantity of aggregated items 𝑐 transferred from aggregated 27 warehouses 𝑤 to demand point 𝑖 within zone 𝑘. Additionally, the parameter 𝛽𝑤𝑖𝑐𝑡(1) is defined to denote the unit cost of transporting aggregated items 𝑐 from aggregated warehouses 𝑤 to demand point 𝑖 within zone 𝑘. This parameter is calculated as 𝛽𝑤𝑖𝑐𝑡(1)= ∑ ∑ 𝑔𝛽𝑓𝑝𝑤𝑐𝛽𝑓𝑖𝑝𝑡𝑓∈𝐹𝑤𝑝∈𝑃𝑐 . Therefore, we can formulate the fist stage disaggregation model (7) in the following way: 𝑍𝐷𝑖𝑠𝑠𝐷𝑘=𝑚𝑖𝑛𝑦(1)∑ 𝛽𝑤𝑖𝑐𝑡(1)∗ 𝑦𝑤𝑖𝑐𝑡(1)𝑤,𝑖∈𝐼𝑘,𝑐,𝑡 (3𝑎) 𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: ∑𝑦𝑤𝑖𝑐𝑡(1)𝑤≥ ?̅?𝑐𝑖𝑡 ∀ 𝑖 ∈ 𝐼𝑘, 𝑐, 𝑡 (7𝑏) ∑𝑦𝑤𝑖𝑐𝑡(1)𝑖∈𝐼𝑘≤ ?̅?𝑤𝑘𝑐𝑡∗ ∀ 𝑤, 𝑐, 𝑡 (7𝑐) 𝑦𝑤𝑖𝑐𝑡(1)∈ 𝑍+ (7𝑑) Where the objective function (7a) minimizes total transportation costs for each subproblem 𝑘. The constraint set (7b) ensures the fulfillment of the demand for aggregated item 𝑐 at each demand point 𝑖 within each zone 𝑘. The constraint set (7c) guarantees that total quantity assigned to all demand points 𝑖 within zone 𝑘 does not surpass the aggregated value ?̅?𝑤𝑘𝑐𝑡∗ . The last constraint (7d) determines the type of the decision variable 𝑦𝑤𝑖𝑐𝑡(1). 5.2.2. Optimal Disaggregation of warehouses Considering the solution of the first stage disaggregation denoted by 𝑦𝑤𝑖𝑐𝑡(1)∗, and the solution of the aggregated model, i.e., (?̅?𝑐𝑤𝑡∗ , ?̅?𝑐𝑤𝑡∗ ), we can formulate the optimization model (8) for warehouses disaggregation by defining new sets of decision variables and parameters. These are outlined below, and the model can be solved independently for each subproblem 𝑤 = 1, . . ,𝑊. The objective function value of this stage will be equal to 𝑍𝐷𝑖𝑠𝑠𝑤 = ∑ 𝑍𝐷𝑖𝑠𝑠𝑊𝑤𝑤∈𝑊 where 𝑍𝐷𝑖𝑠𝑠𝑊𝑤 represents the objective value of solving the subproblem 𝑤. Table 3: The lists of sets and notations used in the warehouse disaggregation problem Sets 28 𝐼𝑤 The subset of demand points supplied from the aggregated warehouse w defined as: 𝐼𝑤 = {𝑖 ∈ 𝐼|𝑦𝑤𝑖𝑐𝑡(1)∗> 0 ∀ 𝑐 ∈ 𝐶, 𝑡 ∈ 𝑇} Parameters β𝑓𝑖𝑐𝑡(2) The unit of shipment cost from the fulfilment centre f ∈ Fw to demand zone 𝑖 ∈𝐼𝑤 for cluster 𝑐 in period 𝑡 ℎ𝑐𝑓𝑡(2) The unit of storage cost of cluster 𝑐 in fulfilment center f ∈ Fw in period 𝑡 𝑟𝑐𝑓𝑡(2) The receiving cost of cluster 𝑐 in the fulfilment center f ∈ Fw in period 𝑡 Variables 𝑥𝑐𝑓𝑡(2) Quantity of cluster 𝑐 shipped to fulfilment center f ∈ Fw in period t 𝑦𝑓𝑖𝑐𝑡(2) Quantity of cluster 𝑐 shipped from fulfilment center f ∈ Fwto demand zone 𝑖 ∈𝐼𝑤 in period t 𝑞𝑐𝑓𝑡(2) The inventory level of fulfilment center f ∈ Fw for cluster 𝑐 in period t The parameters β𝑓𝑖𝑐𝑡(2), ℎ𝑐𝑓𝑡(2), 𝑎𝑛𝑑 𝑟𝑐𝑓𝑡(2) are obtained as follows: β𝑓𝑖𝑐𝑡(2)= ∑ 𝑔𝑝𝑐𝛽𝑓𝑖𝑝𝑡𝑝∈𝑃𝑐 ∀ 𝑓 ∈ 𝐹𝑤, 𝑖 ∈ 𝐼𝑤 , t ∈ T ℎ𝑐𝑓𝑡(2)= ∑ 𝑔𝑝𝑐ℎ𝑓𝑝𝑡𝑝∈𝑃𝑐 ∀ 𝑓 ∈ 𝐹𝑤, t ∈ T 𝑟𝑐𝑓𝑡(2)= ∑ 𝑔𝑝𝑐ℎ𝑓𝑝𝑡𝑝∈𝑃𝑐 ∀ 𝑓 ∈ 𝐹𝑤, t ∈ T Each sub-problem 𝑤 takes the following form: 𝑍𝐷𝑖𝑠𝑠𝑊𝑤=𝑚𝑖𝑛𝑥(2), 𝑞(2), 𝑦(2) ( ∑ 𝑟𝑐𝑓𝑡(2)∗ 𝑥𝑐𝑓𝑡(2)𝑓∈𝐹𝑤,𝑐,𝑡+ ∑ ℎ𝑐𝑓𝑡(2)∗ 𝑞𝑐𝑓𝑡(2)𝑐,𝑓∈𝐹𝑤,𝑡+ ∑ 𝛽𝑓𝑖𝑐𝑡(2)∗ 𝑦𝑓𝑖𝑐𝑡(2)𝑓∈𝐹𝑤,𝑖∈𝐼𝑤,𝑐,𝑡) (4𝑎) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: 29 ∑ 𝑦𝑓𝑖𝑐𝑡(2)𝑓∈𝐹𝑤= 𝑦𝑤𝑖𝑐𝑡(1)∗ ∀ 𝑐, 𝑖 ∈ Iw, 𝑡 (8𝑏) 𝑞𝑐𝑓𝑡(2)= 𝑞𝑐𝑓,𝑡−1(2)+ 𝑥𝑐𝑓𝑡(2)−∑ 𝑦𝑓𝑖𝑐𝑡(2)𝑖∈𝐼𝑤 ∀ 𝑐, 𝑓 ∈ 𝐹𝑤 , 𝑡 (8𝑐) ∑ 𝑥𝑐𝑓𝑡(2)𝑓∈𝐹𝑤= ?̅?∗𝑐𝑤𝑡 ∀ 𝑐, 𝑡 (8𝑑) ∑?̅?𝑐𝑞𝑐𝑓,𝑡−1(2)𝑐+∑?̅?𝑐𝑥𝑐𝑓𝑡(2)𝑐≤ 𝜛𝑓𝑡 ∀ 𝑓 ∈ 𝐹𝑤, 𝑡 (8𝑒) ∑ 𝑞𝑐𝑓𝑡(2)𝑓∈𝐹𝑤= ?̅?𝑐𝑤𝑡∗ ∀ 𝑐, 𝑡 (8𝑓) 𝑥𝑐𝑓𝑡(2), 𝑞𝑐𝑓𝑡(2), 𝑦𝑓𝑖𝑐𝑡(2)∈ 𝑍+ (8𝑔) The objective function (8a) minimizes the total receiving, storage, and transportation costs for each subproblem 𝑤. The constraint set (8b) guarantees that total flow from all fulfilment canters within the cluster 𝑤 must match the flow leaving the cluster 𝑤 for a specific demand point 𝑖 and aggregated item 𝑐. The constraint (8c) represents the flow conservation. The constraint set (8d) ensures that the total quantity assigned to all fulfilment centres within cluster 𝑤 is equal to the total aggregated value. Additionally, the remaining constraints include warehouse capacity constraint (8e), upper bound for the inventory level of warehouses (8f), and (8g) defines defines the domain and type of decision variables. 5.2.3. Optimal Disaggregation of items Here, we formulate an optimization model similar to warehouse disaggregation for this stage using the solutions from the previous stage, denoted as (𝑥𝑐𝑓𝑡(2)∗, 𝑞𝑐𝑓𝑡(2)∗ , 𝑦𝑓𝑖𝑐𝑡(2)∗). The disaggregated solution can be achieved by solving the subproblems for each 𝑐 = 1, . . . , 𝐶. The overall objective function value after disaggregating all clusters 𝑐 = 1, . . . , 𝐶 will be equal to 𝑍𝐷𝑖𝑠𝑠𝑝 = ∑ 𝑍𝐷𝑖𝑠𝑠𝑝𝑐𝑐∈𝐶 where 𝑍𝐷𝑖𝑠𝑠𝑝𝑐 represents the objective value of solving the subproblem 𝑐. We introduce new sets of decision variables including 𝑥𝑝𝑓𝑡(3), 𝑞𝑝𝑓𝑡(3), and 𝑦𝑓𝑖𝑝𝑡(3) to develop the optimization model (9) for each subproblem 𝑐. Since the set of parameters are identical to those in the original problem, hence won't go into additional detail here. Each sub-problem 𝑐 takes the following form: 30 𝑍𝐷𝑖𝑠𝑠𝑐𝑃 =𝑚𝑖𝑛 𝑥(3), 𝑞(3), 𝑦(3)( ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗ 𝑞𝑝𝑓𝑡(3)𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡(3)𝑝∈𝑃𝑐,𝑖,𝑓,𝑡) (5𝑎) 𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: ∑ 𝑦𝑓𝑖𝑝𝑡(3)𝑝∈𝑃𝑐= 𝑦𝑓𝑖𝑐𝑡(2)∗ ∀ 𝑓, 𝑖, 𝑡 (9𝑏) ∑𝑦𝑓𝑖𝑝𝑡(3)𝑓≥ 𝑑𝑝𝑖𝑡 ∀ 𝑝 ∈ 𝑃𝑐 , 𝑖, 𝑡 (9𝑐) 𝑞𝑝𝑓𝑡(3)= 𝑞𝑝𝑓,𝑡−1(3)+ 𝑥𝑝𝑓𝑡(3)−∑𝑦𝑓𝑖𝑝𝑡(3)𝑖 ∀ 𝑝 ∈ 𝑃𝑐 , 𝑓, 𝑡 (9𝑑) ∑ 𝑣𝑝𝑞𝑝𝑓,𝑡−1(3)𝑝∈𝑃𝑐+ ∑ 𝑣𝑝𝑥𝑝𝑓𝑡(3)𝑝∈𝑃𝑐≤ ?̅?𝑐 (𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗ ) ∀ 𝑓, 𝑡 (9𝑒) ∑ 𝑥𝑝𝑓𝑡(3)𝑝∈𝑃𝑐= 𝑥𝑐𝑓𝑡(2)∗ ∀ 𝑓, 𝑡 (9𝑓) ∑ 𝑞𝑝𝑓𝑡(3)𝑝∈𝑃𝑐= 𝑞𝑐𝑓𝑡(2)∗ ∀ 𝑓, 𝑡 (9𝑔) 𝑥𝑝𝑓𝑡(3), 𝑞𝑝𝑓𝑡(3), 𝑦𝑓𝑖𝑝𝑡(3)∈ 𝑍+ (9ℎ) Although this problem has the same objective and constraints as the original (1), it is significantly more computationally tractable as it only concerns the subset of items 𝑃𝑐 in each sub-problem 𝑐. However, the level of complexity of this model is the same as warehouse disaggregation. 5.3. Solution approaches In this section, we will proceed to describe the solution methods employed to deal with each of the disaggregation stages. 5.3.1. Exact method for demand zone disaggregation Due to the total unimodularity property of the constraint matrix used in formulation (7), the first stage disaggregation can be solved optimally. In the following we aim to prove this claim. Lemma 1. The first stage disaggregation problem (7) can be reduced to Linear programming model. 31 Proof. Formulation (7) is an instance of integral network flow problem. Let (𝑉, 𝐴) be the corresponding directed graph with capacity constraint (0 ≤ 𝑥 ≤ 𝑢) and flow conservation constraint 𝐵𝑥 = 0. 𝐵 is the node-arc incidence matrix and can be defined as 𝐵 = 𝐵𝑣,𝑎 ∈{−1,0,1}|𝑉|∗|𝐴| where 𝐵𝑣,𝑎 = 1 if 𝑎 = (𝑠, 𝑣), 𝐵𝑣,𝑎 = −1 if 𝑎 = (𝑣, 𝑠), and 𝐵𝑣,𝑎 = 0 otherwise. According to the Ghouila-Houri theorem (Conforti et al., 2014b), the incidence matrix 𝐵 of any directed graph is TU. Additionally, based on Hoffman and Kruskal theorem (Hoffman et al., 1956), let 𝐵 ∈ 𝑍𝑚×𝑛, for (𝑑, 𝑏, 𝑙, 𝑢) ∈ (𝑅 ∪ {∞})𝑚 × (𝑅 ∪ {−∞})𝑚 × (𝑅 ∪ {−∞})𝑛 × (𝑅 ∪ {∞})𝑛, define the polyhedron Δ as: Δ(𝑑, 𝑏, 𝑙, 𝑢) ≔ {𝑥 ∈ 𝑅𝑛: 𝑑 ≤ 𝐵𝑥 ≤ 𝑏, 𝑎𝑛𝑑 𝑙 ≤ 𝑥 ≤ 𝑢} 𝐵 is TU if and only if Δ(𝑑, 𝑏, 𝑙, 𝑢) is integral for every integer-valued choice of (𝑑, 𝑏, 𝑙, 𝑢). In other words, if the constraint matrix is TU, the polyhera is integral for every integer-valued of right-hand side. Hence, the linear relaxation of the integral network flow problem has an optimal integer-valued vertex solution. Therefore, we can omit the integrality constraints from the problem, and reduce it to the linear program. □ 5.3.2. Approximation method for warehouse disaggregation The complexity of optimization model (8) is equivalent to that of the model (1), which is classified as NP-hard. Accordingly, we can adjust the greedy algorithm we described in Section 4.1 to solve the model (8). Given our prior discussion on the algorithm, a detailed explanation is omitted here. Instead, a summary of the algorithm is provided below to keep the document concise. As mentioned earlier, the warehouse disaggregation model can be broken down into 𝑤 subproblems and solved separately. The greedy algorithm can be employed for each cluster 𝑤, resulting in 𝑍𝐷𝑖𝑠𝑠𝐺 = ∑ 𝑍𝐷𝑖𝑠𝑠𝐺𝑤𝑤, where 𝑍𝐷𝑖𝑠𝑠𝐺𝑤 represents the cost of greedy algorithm for disaggregating each group of 𝑤. The greedy algorithm for solving the disaggregated model 𝑤 requires a group of bundles 𝜑 =(𝑐, 𝑖) for each period where 𝑐 and 𝑖 stand for the items and demand points respectively. The value of 𝑦𝑓∗𝜑𝑡(2) is determined by the minimum of 𝑦𝑤𝜑𝑡(1)∗ and the current capacity of warehouse 𝑓∗, i.e., 32 𝑦𝑓∗𝜑𝑡(2)= min {𝑦𝑤𝜑𝑡(1)∗, (𝜛𝑓∗𝑡 − ?̅?𝑓∗𝑡) ?̅?𝑐⁄ } where ?̅?𝑓𝑡 is the utilized capacity of 𝑓∗. Subsequently, we define new bundle 𝜑1 = (𝑓∗, 𝑐) to determine the variables 𝑥𝑓𝑐𝑡(2) and 𝑞𝑓𝑐𝑡(2). The value of 𝑥𝜑1𝑡(2) is determined based on the minimum of (∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 + 𝑞𝜑1𝑡(2)) and total current supply (𝜂𝑐𝑡 − 𝜀?̅?𝑡) for item c, i.e., 𝑥𝜑1𝑡(2)= 𝑚𝑖𝑛 {(∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 + 𝑞𝜑1𝑡(2)) , 𝜂𝑐𝑡 − 𝜀?̅?𝑡}; here, 𝜂𝑐𝑡 stands for the total supply for item c, i.e., 𝜂𝑐𝑡 = ?̅?∗𝑐𝑤𝑡 and 𝜀?̅?𝑡 is the utilized portion of supply 𝜂𝑐𝑡 .The summary of algorithm is provided as follows: Table 4. The summary of greedy algorithm for each cluster (w) Step 1 Define 𝑁𝑤 as the total number of (𝐶, 𝐼) combinations in cluster w Set 𝑡 = 𝑇, 𝑞𝜑1𝑇(2)= 0 Step 2 Set ?̅?𝑓𝑡 = 0, 𝜀?̅?𝑡 = 0, 𝜂𝑐𝑡 = ?̅?∗𝑐𝑤𝑡 Step 3 Define 𝜑 = (𝑐, 𝑖), sort them in decreasing order of potential function 𝜌 =∑ (𝑟𝑐𝑓𝑡(2)+ℎ𝑐𝑓𝑡(2)+β𝑓𝜑𝑡(2))𝑓∈Fw 𝑦𝑤𝜑𝑡(1)∗?̅?𝑐 Step 4 Set 𝜑 = 1 (represents the bundle (p, 𝑖) with the highest potential function) Step 5 Choose a warehouse with the least cost where 𝑓∗ = 𝐴𝑟𝑔𝑚𝑖𝑛 { (𝑟𝑐𝑓𝑡(2)+ ℎ𝑐𝑓𝑡(2)+β𝑓𝜑𝑡(2))} If ?̅?𝑓∗𝑡 + 𝑣𝜑 ≤ 𝜛𝑓∗𝑡: Assign bundle 𝜑 to the warehouse 𝑓∗ 𝑦𝑓∗𝜑𝑡(2)= min {𝑦𝑤𝜑𝑡(1)∗, (𝜛𝑓∗𝑡 − ?̅?𝑓∗𝑡) ?̅?𝑐⁄ } ?̅?𝑓∗𝑡 = ?̅?𝑓∗𝑡 + ?̅?𝑐 ∗ 𝑦𝑓∗𝜑𝑡(2) 𝑦𝑤𝜑𝑡(1)∗= 𝑦𝑤𝜑𝑡(1)∗− 𝑦𝑓∗𝜑𝑡(2) Else move to step 6 Step 6 While 𝑦𝑤𝜑𝑡(1)∗> 0, then 𝑓∗ = 𝑓∗ + 1, move to step 5 Step 7 Set 𝜑 = 𝜑 + 1, move to step 5 Step 8 Define 𝜑1 = (𝑓∗, 𝑐), sort them in decreasing order of (∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 + 𝑞𝜑1𝑡(2) ) 33 Step 9 If ε𝑐𝑡 ≤ 𝜂𝑐𝑡: 𝑥𝜑1𝑡(2)= 𝑚𝑖𝑛 {(∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 + 𝑞𝜑1𝑡(2)) , 𝜂𝑐𝑡 − 𝜀𝑐𝑡} 𝑞𝜑1,𝑡−1(2)= ∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 − 𝑥𝜑1𝑡(2)+ 𝑞𝜑1𝑡(2) 𝜀𝑐𝑡 = 𝜀𝑐𝑡 + 𝑥𝜑1𝑡(2) Else 𝑥𝜑1𝑡(2)= 0, 𝑞𝜑1,𝑡−1(2)= ∑ 𝑦𝜑1,𝑖,𝑡(2)𝑖 + 𝑞𝜑1𝑡(2) Step 10 Set 𝜑1 = 𝜑1 + 1, move to step 9 Step 11 ?̅?𝑓,𝑡−1 = ?̅?𝑓,𝑡−1 + ?̅?𝑐 ∗ 𝑞𝜑1,𝑡−1(2) Step 12 Set 𝑡 = 𝑡 − 1, move to step 3 5.3.3. Approximation method for item disaggregation The optimization model (9) is characterized as NP-hard, the same as model (1). Consequently, we can adjust and apply the greedy algorithm described earlier to solve the model (9) in polynomial time. The items disaggregation model can be decomposed into 𝑐 subproblems and solved individually. The greedy algorithm can be utilized for each cluster 𝑐, resulting in 𝑍𝐷𝑖𝑠𝑠𝑝 =∑ 𝑍𝐷𝑖𝑠𝑠𝑝𝑐𝑐∈𝐶 , where 𝑍𝐷𝑖𝑠𝑠𝑝𝑐 represents the cost of greedy algorithm for disaggregating each cluster 𝑐. We employ a similar approach used for solving the warehouse disaggregation, with the exception that certain parameters are modified as outlined in the table below. Since we are solving the algorithm for each cluster 𝑐, we first define bundle 𝜑 = (𝑝, 𝑖) for 𝑝 ∈ 𝑝𝑐. Subsequently, we introduce parameter 𝛾f𝑡 representing the warehouse capacity for each cluster 𝑐, defined as 𝛾f𝑡 =?̅?𝑐(𝑥𝑐𝑓𝑡(2)𝐺 + 𝑞𝑐𝑓,𝑡−1(2)𝐺). In addition, the supply value for each cluster c equals to 𝜂′𝑓𝑡 = 𝑥𝑐𝑓𝑡(2)𝐺. We also incorporate new conditions in step 5, ensuring that the total values assigned from warehouse 𝑓 to demand point 𝑖 for a specific cluster 𝑐 do not exceed the predetermined flow obtained from the previous stage, i.e., 𝑦𝑓∗𝑖𝑐𝑡(2)𝐺. The summary of algorithm is provided as follows: Table 5. The summary of the greedy algorithm for each cluster (c) Step 1 Define 𝑁𝑐 as the total number of (𝑃, 𝐼) combinations in cluster 𝑐 34 Set 𝑡 = 𝑇, 𝑞(3)φ1,T= 0 Step 2 Set ?̅?𝑓𝑡 = 0, 𝜀?̅?𝑡 = 0, ?̅?𝑓𝑖𝑡 = 0, 𝛾f𝑡 = ?̅?𝑐(𝑥𝑐𝑓𝑡(2)𝐺 + 𝑞𝑐𝑓,𝑡−1(2)𝐺), 𝜂′𝑓𝑡 = 𝑥𝑐𝑓𝑡(2)𝐺 Step 3 Define 𝜑 = (p, 𝑖), sort them in decreasing order of potential function 𝜌 =∑ (𝑟𝑝𝑓𝑡+ℎ𝑝,𝑓,𝑡+β𝑓𝜑𝑡)𝑓∈𝐹 𝑑𝜑𝑡𝑣𝑝 Step 4 Set 𝜑 = 1 (represents the bundle (p, 𝑖) with the highest potential function) Step 5 Choose a warehouse with the least cost where 𝑓∗ = 𝐴𝑟𝑔𝑚𝑖𝑛{ (𝑟𝑝𝑓𝑡 + ℎ𝑝,𝑓,𝑡 +β𝑓𝑖𝑝𝑡)} If 𝑦𝑓∗𝑖𝑐𝑡(2)𝐺 − ?̅?𝑓∗𝑖𝑡 > 0: If ?̅?𝑓∗𝑡 + 𝑣𝑝 ≤ 𝛾𝑓∗𝑡: Assign bundle 𝜑 to the warehouse 𝑓∗ 𝑦𝑓∗𝜑𝑡(3)= min {𝑑𝜑𝑡, (𝛾𝑓∗𝑡 − ?̅?𝑓∗𝑡) 𝑣𝑝⁄ } ?̅?𝑓∗𝑡 = ?̅?𝑓∗𝑡 + 𝑣𝑝 ∗ 𝑦𝑓∗𝜑𝑡(3) 𝑑𝜑𝑡 = 𝑑𝜑𝑡 − 𝑦𝑓∗𝜑𝑡(3) ?̅?𝑓∗𝑖𝑡 = ?̅?𝑓∗𝑖𝑡 + 𝑦𝑓∗𝜑𝑡(3) Else move to step 6 Else move to step 6 Step 6 While 𝑑𝜑𝑡 > 0, then 𝑓∗ = 𝑓∗ + 1, move to step 5 Step 7 Set 𝜑 = 𝜑 + 1, move to step 5 Step 8 Define 𝜑1 = (𝑓∗, p), sort them in decreasing order of (∑ 𝑦𝑓∗𝜑𝑡(3)𝑖 + 𝑞𝜑1𝑡(3) ) Step 9 If 𝜀?̅?𝑡 < 𝜂′𝑓𝑡: 𝑥𝜑1𝑡(3)= 𝑚𝑖𝑛 {(∑ 𝑦𝜑1,𝑖,𝑡(3)𝑖 + 𝑞𝜑1𝑡(3)) , 𝜂′𝑓𝑡 − 𝜀?̅?𝑡} 𝑞𝜑1,𝑡−1(3)= ∑ 𝑦𝜑1,𝑖,𝑡(3)𝑖 − 𝑥𝜑1𝑡(3)+ 𝑞𝜑1𝑡(3) 𝜀?̅?𝑡 = 𝜀?̅?𝑡 + 𝑥𝜑1𝑡(3) Else 𝑥𝜑1𝑡(3)= 0, 𝑞𝜑1,𝑡−1(3)= ∑ 𝑦𝜑1,𝑖,𝑡(3)𝑖 + 𝑞𝜑1𝑡(3) 35 Step 10 Set 𝜑1 = 𝜑1 + 1, move to step 9 Step 11 ?̅?𝑓,𝑡−1 = ?̅?𝑓,𝑡−1 + 𝑣𝑝 ∗ 𝑞𝜑1,𝑡−1(3) Step 12 Set 𝑡 = 𝑡 − 1, move to step 3 5.4. The algorithm error bound In the aggregation approach, a crucial aspect involves quantifying the extent of accuracy loss resulting from aggregation. While accurately calculating this difference necessitates solving both the original and aggregated problems, it's possible to establish bounds for the loss of accuracy without solving the original problem. Typically, two types of error bounds are considered: a priori and a posteriori. A priori error bounds are determined after aggregation but before optimizing the aggregated problem, whereas calculating a posteriori error bound necessitates knowledge of an optimal solution for the aggregated problem. A priori error bounds tend to be more conservative than a posteriori ones and are often employed to guide decisions regarding clustering and cluster procedures, despite no evidence suggesting that methods yielding the tightest a priori error bound necessarily yield the tightest a posteriori error bound (Litvinchev & Tsurkov, 2003; Rogers et al., 1991). Using duality theory, several authors measure the error bounds by comparing the actual optimal value of the objective function with the optimal value for the aggregated problem, while ignoring the disaggregation part. For example, Zipkin (1980c) develops a bound on the optimality gap for linear programming models resulting from column aggregations, as well as both column and row aggregations (Zipkin, 1980a). Subsequently, the improvement on the Zipkin’s bound are presented by Mendelssohn (1980) and Shetty & Taylor (1987). A. Hallefjord & Storoey (1986) also modify one version of Zipkin’s bound for the generalized assignment problem. To the best of our knowledge, Chen & Graves (2021) is the only work employing a series of intermediate formulations based on their problem structure to bound the cost between the original and the disaggregated final solutions. In this section, our goal is to establish a posteriori error bound for the algorithm which is derived from the disaggregated solution. Before delving into that, it's essential to demonstrate that the solution derived from the disaggregation remains feasible for the original model. 36 5.4.1. Proof of feasibility To derive a feasible solution to the optimal problem (1) from the aggregated solution, we conducted three stages of hierarchical partial disaggregation, as previously described. Here, our objective is to demonstrate that the solution obtained from the disaggregation stages satisfies the constraints of the original problem. At each stage of disaggregation, we solve a series of subproblems for each 𝑘,𝑤 and 𝑐 subject to their respective constraints. In addition, the hierarchical approach to disaggregation maintains solution consistency, ensuring that the final disaggregated solution from the last step does not violate any of the original constraints. For instance, 𝑦𝑓𝑤𝑖𝑘𝑝𝑐𝑡(3) obtained from the last step of the disaggregation ensures to fulfill the demand 𝑑𝑝𝑖𝑡. This solution depends on the outcomes of previous steps, thereby ensuring overall consistency. ∑𝑦𝑤𝑐𝑖𝑘𝑡(1)𝑤=∑ ∑ 𝑦𝑓𝑤𝑖𝑘𝑐𝑡(2)𝑓∈𝐹𝑊𝑤=∑∑ ∑ 𝑦𝑓𝑤𝑖𝑘𝑝𝑐𝑡(3)𝑓∈𝐹𝑊≥ 𝑑𝑝𝑖𝑡𝑝∈𝑃𝑐𝑤 (6) Similarly, this holds true for 𝑥𝑝𝑓𝑡(3) and 𝑞𝑝𝑓,𝑡(3): ∑?̅?𝑐𝑤𝑡𝑤=∑ ∑ 𝑥𝑐𝑓𝑡(2)𝑓∈𝐹𝑤𝑤= ∑ (∑ ∑ 𝑥𝑝𝑓𝑡(3)𝑓∈𝐹𝑤𝑤≤ 𝛼𝑝𝑡)𝑝∈𝑃𝑐≤ ?̅?𝑐𝑡 (7) ∑ 𝑣𝑝𝑞𝑝𝑓,𝑡−1(3)𝑝∈𝑃𝑐+ ∑ 𝑣𝑝𝑥𝑝𝑓𝑡(3)𝑝∈𝑃𝑐≤ ?̅?𝑐𝑥𝑐𝑓𝑡(2) ≤ 𝜛𝑓𝑡 (8) 5.4.2. Proof of optimality gap Any feasible solution to the original minimization problem provides an upper bound on its optimal objective value. Therefore, to determine the gap between the disaggregated and optimal solutions of the original model, we must compute: 𝒁𝑫𝒊𝒔𝒔𝒑∗− 𝒁∗𝒐𝒑𝒕 =∑𝒁𝑫𝒊𝒔𝒔𝒑∗𝒄𝒄∈?̅?− 𝒁∗𝑶𝒑𝒕 ≤ |𝑪|𝐦𝐚𝐱𝒄 𝒁𝑫𝒊𝒔𝒔𝒑∗𝒄− 𝒁∗𝑹 (9) Where 𝑍∗𝑅 represents the optimal solution of the relaxed original problem (14) formulated as follows: 37 𝑍𝑅 = 𝑚𝑖𝑛 (∑𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡𝑓+ ∑ ℎ𝑝𝑓𝑡 ∗ 𝑞𝑝𝑓𝑡𝑝,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡𝑓,𝑖,𝑝,𝑡) (10𝑎) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜: Constraints (1b) to (1e) 𝑥𝑝𝑓𝑡, 𝑞𝑝𝑓𝑡, 𝑦𝑓𝑖𝑝𝑡 ∈ 𝑅+ (14𝑏) To find the optimal solution of the relaxed problem, 𝑍∗𝑅, we can use the concept of LP duality. The dual form of the problem (14) is formulated as follows: 𝑍𝑅𝐷 =𝑀𝑎𝑥𝑢, ζ, Θ(∑u𝑝𝑖𝑡𝑑𝑝𝑖𝑡𝑝,𝑖,𝑡−∑ζ𝑝𝑡 𝛼𝑝𝑡𝑝,𝑡−∑Θ𝑓𝑡𝜛𝑓𝑡𝑓,𝑡) (11𝑎) Subject to: u𝑝𝑖𝑡 + Γ𝑝𝑓𝑡 ≤ 𝛽𝑓𝑖𝑝𝑡 ∀ 𝑝, 𝑖, 𝑡 (15𝑏) −Γ𝑝𝑓𝑡 − ζ𝑝𝑡 − 𝑣𝑝Θ𝑓𝑡 ≤ 𝑟𝑝𝑓𝑡 ∀ 𝑝, 𝑓, 𝑡 (15𝑐) Γ𝑝𝑓,𝑡 ≤ ℎ𝑝𝑓,𝑡 ∀ 𝑝, 𝑓, 𝑡 (15𝑑) −Γ𝑝𝑓,𝑡+1 − 𝑣𝑝Θ𝑓,𝑡+1 ≤ ℎ𝑝𝑓,𝑡 ∀ 𝑝, 𝑓, 𝑡 = 2,…𝑇 (15𝑒) u𝑝𝑖𝑡, ζ𝑝𝑡, Θ𝑓𝑡 ≥ 0 ∀ 𝑝, 𝑓, 𝑖, 𝑡 (15𝑓) Where (u𝑝𝑖𝑡, ζ𝑝𝑡, Θ𝑓𝑡, Γ𝑝𝑓𝑡) are introduced as decision variables of dual model (14) which are assigned as follows: • u𝑝𝑖𝑡 is a decision variable for constraint (1b) • Γ𝑝𝑓𝑡 is considered for constraint (1c) • ζ𝑝𝑡 is considered for constraint (1d) • Θ𝑓𝑡 is considered for constraint (1e) Using strong duality, the objective function value of the optimal relaxed original model (14) is equal to its corresponding dual objective value (15), i.e., 𝑍𝑅∗𝐷 = 𝑍∗𝑅 , therefore we have: 38 𝑍∗𝑅 =∑𝑢𝑝𝑖𝑡∗ 𝑑𝑝𝑖𝑡𝑝,𝑖,𝑡−∑𝜁𝑝𝑡∗ 𝛼𝑝𝑡𝑝,𝑡−∑Θ𝑓𝑡∗ 𝜛𝑓𝑡𝑓,𝑡 (12) Furthermore, to find the upper-bound for the first term of the optimality gap, i.e., |𝐶|max𝑐 𝑍𝐷𝑖𝑠𝑠𝑝∗𝑐, we can define the following optimization model: |𝐶|𝑚𝑎𝑥𝑐 𝑍𝐷𝑖𝑠𝑠𝑝∗𝑐= |𝐶|𝑚𝑎𝑥𝑐𝑚𝑖𝑛 ( ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗ 𝑞𝑝𝑓𝑡(3)∗𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡)= |𝐶|𝑚𝑎𝑥 𝛿 (13𝑎) Subject to: 𝛿 ≤ ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗ 𝑞𝑝𝑓𝑡(3)∗𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡 ∀𝑐 ∈ 𝐶 (17𝑏) Constraints (9b) to (9h). If we rewrite 𝑞𝑝𝑓𝑡(3)∗= ∑ (𝑥𝑝𝑓𝑗(3)∗− ∑ 𝑦𝑓𝑖𝑝𝑗(3)∗𝑖 )𝑡𝑗=1 ∀ 𝑝 ∈ 𝑃𝑐, 𝑓, 𝑡, we obtain: 𝛿 ≤ ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ ℎ𝑝𝑓𝑡 ∗∑ (𝑥𝑝𝑓𝑗(3)∗ −∑𝑦𝑓𝑖𝑝𝑗(3)∗𝑖)𝑡𝑗=1𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ 𝛽𝑓𝑖𝑝𝑡 ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡 ∀𝑐 ∈ 𝐶 Then, we have: 𝛿 ≤ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡− ∑ ℎ𝑝𝑓𝑡 (∑ 𝑦𝑓𝑖𝑝𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑖,𝑡 ∀𝑐 ∈ 𝐶 (17𝑐) Using Lagrangian relaxation on constraint (17b), the following relations hold for all 𝜆𝑐 ≥ 0: 39 |𝐶|𝛿∗ ≤ |𝐶|𝛿∗ +∑𝜆𝑐𝑐( ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡− ∑ ℎ𝑝𝑓𝑡 (∑ 𝑦𝑓𝑖𝑝𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑖,𝑡− 𝛿∗)≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡] (17𝑑) Since the inventory level is non-negative, 𝑞𝑝𝑓𝑡(3)∗≥ 0, we have 𝑞𝑝𝑓,𝑡−1(3)∗+ 𝑥𝑝𝑓𝑡(3)∗− ∑ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑖 ≥ 0 ∀ 𝑝 ∈𝑃𝑐, 𝑓, 𝑡, then we have ∑ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑖 ≤ 𝑞𝑝𝑓,𝑡−1(3)∗ + 𝑥𝑝𝑓𝑡(3)∗, which in turn ∑ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑖 ≤ ∑ 𝑥𝑝𝑓𝑗(3)∗𝑡𝑗=1 ∀ 𝑝 ∈ 𝑃𝑐 , 𝑓, 𝑡. We can rewrite the above expression as follows: 40 (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ (𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗ 𝑦𝑓𝑖𝑝𝑡(3)∗𝑝∈𝑃𝑐,𝑖,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ 𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗∑𝑦𝑓𝑖𝑝𝑡(3)∗𝑖𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ 𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡 − ℎ𝑝𝑓𝑡) ∗∑ 𝑥𝑝𝑓𝑗(3)∗𝑡𝑗=1𝑝∈𝑃𝑐,𝑓,𝑡+ ∑ ℎ𝑝𝑓𝑡 (∑ 𝑥𝑝𝑓𝑗(3)∗𝑡−1𝑗=1)𝑝∈𝑃𝑐,𝑓,𝑡] 41 = (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ (𝑟𝑝𝑓𝑡 + ℎ𝑝𝑓𝑡) ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ 𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡) ∗∑ 𝑥𝑝𝑓𝑗(3)∗𝑡𝑗=1𝑝∈𝑃𝑐,𝑓,𝑡− ∑ ℎ𝑝𝑓𝑡 (𝑥𝑝𝑓𝑡(3)∗)𝑝∈𝑃𝑐,𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗ +∑𝜆𝑐𝑐[ ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ 𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡) ∗∑ 𝑥𝑝𝑓𝑗(3)∗𝑡𝑗=1𝑝∈𝑃𝑐,𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[ ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ |𝑡|𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡) ∗ 𝑚𝑎𝑥𝑗=1,..𝑡 (𝑥𝑝𝑓𝑗(3)∗)𝑝∈𝑃𝑐,𝑓,𝑡] According to the inequality (9e), we have ∑ 𝑥𝑝𝑓𝑡(3)∗𝑝∈𝑃𝑐 ≤?̅?𝑐(𝑥𝑐𝑓𝑡(2)∗+𝑞𝑐𝑓,𝑡−1(2)∗)𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝} ∀ 𝑓, 𝑡, therefore we can obtain the following results: (|𝐶| −∑𝜆𝑐𝑐)𝛿∗ +∑𝜆𝑐𝑐[ ∑ 𝑟𝑝𝑓𝑡 ∗ 𝑥𝑝𝑓𝑡(3)∗𝑓,𝑝∈𝑃𝑐,𝑡+ ∑ |𝑡|𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡) ∗ 𝑚𝑎𝑥𝑗=1,..𝑡 (𝑥𝑝𝑓𝑗(3)∗)𝑝∈𝑃𝑐,𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[∑𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 ∑ 𝑥𝑝𝑓𝑡(3)∗𝑝∈𝑃𝑐𝑓,𝑡+ ∑ |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡) ∗ ∑ 𝑚𝑎𝑥𝑗=1,..𝑡 (𝑥𝑝𝑓𝑗(3)∗)𝑝∈𝑃𝑐𝑝∈𝑃𝑐,𝑓,𝑡] 42 ≤ (|𝐶| −∑ 𝜆𝑐𝑐) 𝛿∗+∑ 𝜆𝑐𝑐[∑𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡?̅?𝑐 (𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗ )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑓,𝑡+∑|𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐 (𝛽𝑓𝑖𝑝𝑡) ∗ (?̅?𝑐 (𝑥𝑐𝑓𝑗(2)∗ + 𝑞𝑐𝑓,𝑗−1(2)∗ )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝})𝑓,𝑡]= (|𝐶| −∑ 𝜆𝑐𝑐) 𝛿∗+∑ 𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐 (𝛽𝑓𝑖𝑝𝑡))?̅?𝑐 (𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗ )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑓,𝑡] Since 𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗ ≤ ∑ 𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗𝑐 , we have: (|𝐶| −∑𝜆𝑐𝑐)𝛿∗ +∑𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐 (𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗ )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐 (∑ 𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗𝑐 )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑓,𝑡] According to the inequality (8e): ∑𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗𝑐≤𝜛𝑓𝑡𝑚𝑖𝑛𝑐{?̅?𝑐} Then, 43 (|𝐶| −∑𝜆𝑐𝑐)𝛿∗ +∑𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐 (∑ 𝑥𝑐𝑓𝑡(2)∗ + 𝑞𝑐𝑓,𝑡−1(2)∗𝑐 )𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑓,𝑡]≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐(𝜛𝑓𝑡)𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}𝑚𝑖𝑛𝑐{?̅?𝑐}𝑓,𝑡] As a result, for all 𝜆𝑐 ≥ 0: |𝐶|max𝑐 𝑍𝐷𝑖𝑠𝑠𝑝∗𝑐≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[∑(𝑚𝑎𝑥𝑝∈𝑃𝑐 𝑟𝑝𝑓𝑡𝑓,𝑡+ |𝑡|𝑚𝑎𝑥𝑖,𝑝∈𝑃𝑐(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐(𝜛𝑓𝑡)𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}min𝑐{?̅?𝑐}] (17𝑔) Consequently, 𝑍𝐷𝑖𝑠𝑠𝑝∗ − 𝑍∗𝑂𝑅𝐺 ≤ (|𝐶| −∑𝜆𝑐𝑐)𝛿∗+∑𝜆𝑐𝑐[∑𝑚𝑎𝑥𝑝∈𝑃𝑐 (𝑟𝑝𝑓𝑡 + |𝑡|𝑚𝑎𝑥𝑖(𝛽𝑓𝑖𝑝𝑡))?̅?𝑐(𝜛𝑓𝑡)𝑚𝑖𝑛𝑝∈𝑃𝑐{𝑣𝑝}min𝑐{?̅?𝑐}𝑓,𝑡]− (∑𝑢𝑝𝑖𝑡∗ 𝑑𝑝𝑖𝑡𝑝,𝑖,𝑡−∑𝜁𝑝𝑡∗ 𝛼𝑝𝑡𝑝,𝑡−∑Θ𝑓𝑡∗ 𝜛𝑓𝑡𝑓,𝑡) (17ℎ) 44 6. Computational experiments In this section, we aim to implement our proposed algorithm on a specific instance of an online retail inventory system and assess the algorithm performance. First, we create a synthetic network inspired from Amazon fulfillment network in the United States. The network characteristics, which include fulfillment centers and demand zone features, cost parameters, and item attributes are detailed in Section 6.1. Then we provide the following computational experiments: • In Section 6.2, we investigate the impact of exclusive cluster counts on the optimality gap. • In Section 6.3, we study how cluster combinations affect a posteriori optimality gap and computation time. • In Section 6.4, we analyze how variations in item sizes affect both the algorithm’s running time and the optimality gap. • In Section 6.5, we study how changes in the ratio of the objective function parameters affect the algorithm performance. All experiments are implemented in Python 3.10.14 and its interface to Gurobi 9.5.2. We use UBC computational cluster (Sockeye) to run our experiments under the following configuration: Operating system: Enterprise Linux 7 (EL7) CPU: 2 x Intel Xeon Gold 6130 (2.1GHz) Total Cores: 32 Mem per Core: 24GB 6.1. Simulation description Our goal is to create a realistic simulation of the Amazon fulfillment network. Given the computational intensity of considering the actual number of demand points, fulfillment centers, and items, we select a representative subset. Our example involves |𝑃| = 1000 items over |𝑇| =4 planning periods. Most parameters follow a structure similar to that described by Chen & Graves, (2021). Based on the dataset from Amazon (2022), there are a total of |𝐹| = 545 fulfilment centers in the US, from which we select 200 major centers for our case. The total capacities of fulfillment 45 centers are determined by multiplying the total demand of all items by their volume. In other words, we assume that ∑ 𝜛𝑓𝐹𝑓=1 = ∑ ∑ ∑ 𝑑𝑝𝑖𝑡 ∗ 𝑣𝑝𝐼𝑖=1𝑃𝑝=1𝑇𝑡=1 . This total capacity is then divided among the fulfillment centers in proportion to their square footage, which is uniformly generated between the minimum and maximum possible sizes. According to Amazon (2024) the size of fulfilment centers varies from 600,000 to 1 million square feet. The number of counties in the United States is utilized to determine demand points. The US Census data from 2022 indicates there are 3,145 counties. We select a subset of counties with populations exceeding 160,000, resulting in |𝐼| = 410 demand points. Demand for each item at each demand point is determined proportionally to the population. First, we generate the total demand across all demand points, then distribute it based on each demand point's population. We assume total demand of each item ranges within 𝑑𝑝𝑡 ∈ [10000, 20000], and demand at each demand point is given by 𝑑𝑝𝑖𝑡 ∈ 𝑑𝑝𝑡 ∗ 𝜀𝑖 , where 𝜀𝑖 is the population ratio of demand point 𝑖. Shipping costs are calculated using a piecewise-linear function of the distance between fulfillment centers and demand points, as well as item weight, following the cost structure described in Chen & Graves, (2021). 𝛽𝑓𝑖𝑝 = {0.247847𝜚𝑓𝑖𝜏𝑝 − 0.227878𝜚𝑓𝑖 + 0.132063𝜏𝑝 + 7.553833, 𝜚𝑓𝑖 ≤ 1 1.37674𝜚𝑓𝑖𝜏𝑝 + 4.02487𝜚𝑓𝑖 + 0.75604𝜏𝑝 + 11.10890, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 Where 𝜏𝑝 represents the weight of item 𝑝 and 𝜚𝑓𝑖 is the spherical distance (in thousands of miles) between fulfillment center 𝑓 and region 𝑖 . The haversine formula is used to calculate 𝜚𝑓𝑖 which uses the longitude and latitude coordinates of the demand points and warehouses. Other parameter values are defined as follows: 𝑣𝑝 ∈ [1,35] in cubic meters 𝜏𝑝 ∈ [0,36] in pounds 𝑟𝑝𝑓𝑡 ∈ [1,6] in USD ℎ𝑝𝑓𝑡 ∈ [1,5] in USD 46 In this experiment, we use k-means clustering technique for aggregation. The normalized vectors of (𝜏𝑝, 𝑣𝑝, 𝑑𝑝𝑡, ℎ̃𝑝, ?̃?𝑝) are used as the clustering features for items aggregation, with ℎ̃𝑝 and ?̃?𝑝 standing for the average holding and receiving costs across all warehouses and periods. Furthermore, we use the latitude and longitude coordinates for warehouses and demand zone aggregation. The aggregated parameters are determined as follows: ?̅?𝑐𝑤𝑡 = ∑ ∑𝑟𝑝𝑓𝑡|𝑃𝑐||𝐹𝑤|𝑝∈𝑃𝑐𝑓∈𝐹𝑤 ∀ 𝑐, 𝑤, 𝑡 ℎ̅𝑐𝑤 = ∑ ∑ℎ𝑝𝑓|𝑃𝑐||𝐹𝑤|𝑝∈𝑃𝑐𝑓∈𝐹𝑤 ∀ 𝑐, 𝑤, 𝑡 ?̅?𝑤𝑘𝑐 =∑ ∑ ∑𝛽𝑓𝑖𝑝|𝑃𝑐||𝐹𝑤||𝐼𝑘|𝑝∈𝑃𝑐𝑓∈𝐹𝑤𝑖∈𝐼𝑘 ∀ 𝑤, 𝑘, 𝑐 ?̅?𝑐 = 𝑚𝑎𝑥𝑝∈𝑃𝑐{𝑣𝑝} ∀ 𝑐 Considering the instance of our model which includes |𝑃| = 1000, |𝑇| = 4, |𝐼| = 410, |𝐹| =200, the optimal model (1) has an objective function value of Z𝑂𝑃𝑇∗ = 532,936,115.33 and requires 3.68 hours to compute. The solution of the greedy algorithm for model (1) is Z𝑂𝑃𝑇𝐺 =549645540.5 with a running time of 21.5 minutes. 6.2. Effect of exclusive cluster size on optimality gap This experiment investigates how the size of each cluster (|𝐶|, |𝐾|, and |𝑊|) exclusively affects the optimality gap. It aims to understand the relationship between different cluster sizes and their exclusive effects on the algorithm performance. Four different values are assigned to each cluster: |𝐶| = (10, 50, 100, 200), |𝐾| = (4, 20, 40, 100), and |𝑊| = (2, 10, 20, 50), resulting in a total of 43 = 64 cluster combinations. We first investigate how variations in item cluster 𝐶 affect the optimality gap for a fixed combination of 𝐾 and 𝑊, and similarly replicate the same experiment for 𝐾 and 𝑊. Hence, this comprehensive analysis not only enables us to comprehend the relationship among different clusters but also effectively identifies the appropriate cluster combination. Thus, the experiment is 47 conducted for 16 fixed combinations of two clusters, e.g., 𝐾 and 𝑊, to evaluate the impact of variations in one cluster, e.g., 𝐶, on algorithm performance. The results of this experiment are depicted in Figure 3, and Figure 4, with the combinations arranged in the order of |𝐶|, |𝐾|, and |𝑊|, respectively. The cluster being analyzed is denoted by (*). For instance, when we refer to (*)-4-2, it indicates that the analysis focuses on cluster C, with the sizes of |𝐾| and |𝑊| being held constant at 4 and 2. The same logic applies to the remaining experiments. Based on Figure 3, there is a linear relationship between the size of |𝐶| and optimality gap. Regardless of the sizes of |𝐾| and |𝑊|, we consistently observe improved performance as the size of C increases. This improvement is due to the fact that the accuracy of algorithm is contingent upon the aggregated model (6). As the size of |𝐶| increases, the model (6) becomes a more precise approximation of the optimal model (1), leading to an enhancement of the overall algorithm performance. Figure 3. Optimality gap versus cluster size |C| However, when we repeat the experiment for |𝐾| and |𝑊|, we do not observe a similar linear relationship. Generally, the algorithm's performance is enhanced by increasing the size of 𝑊; yet 10%12%14%16%18%20%22%24%26%10 50 100 200Optimality gapItem cluster counts, |C|(*)-4-2 (*)-4-10 (*)-4-20 (*)-4-50 (*)-20-2 (*)-20-10(*)-20-20 (*)-20-50 (*)-40-2 (*)-40-10 (*)-40-20 (*)-40-50(*)-100-2 (*)-100-10 (*)-100-20 (*)-100-50 48 this is only true when |𝐾| is equal to or greater than |𝑊|. For instance, the combination |𝐾| =100, |𝑊| = 50 yields a superior bound (i.e., 11%) than |𝐾| = 100, |𝑊| = 2 (𝑖. 𝑒. , 16%) for a fixed size of |𝐶|. Therefore, as long as |𝑊| ≤ |𝐾|, an improved bound is achieved by increasing |𝑊| for fixed values of |𝐾| and |𝐶|. As illustrated in Figure 4, the majority of cases do not yield favorable results when the size of |𝑊| exceeds the size of |𝐾|. For instance, the maximum gap of 25% occurs at |𝐾| = 4, and |𝑊| =50 for |𝐶| = 10, whereas the gap is reduced to 18% at |𝐾| = 100, |𝑊| = 2. This pattern also confirms that the optimality gap will be reduced if the size of |𝐾| is greater than the size of |𝑊|. Increasing the size of |𝐾| improves the optimality gap in cases when |𝐾| < |𝑊|. For instance, the gap is reduced by nearly 9% when |𝐾| is increased from 4 to 20 for |𝐶| = 10 and |𝑊| = 50. It is worth noting that the optimality gap is not always positively influenced by an increase in |𝐾|, as illustrated in Figure 4. In other words, the gap either remains constant or becomes worse as |𝐾| increases when |𝐾| exceeds |𝑊|. This is evident in Figure 4, where the optimality gap increases by 2% when the value of |𝐾| is increased from 4 to 20 for the blue and gray lines with |𝑊| = 2. In addition, the optimality gap remains constant at approximately 12.8% when the value of |𝐾| is increased from 4 to 40, for a fixed values of |𝑊| and |𝐶|, when |𝐶| = 200 and |𝑊| = 2. It does not improve the approximation; rather, it merely increases computation time. One possible explanation for the improved performance observed when the size of |𝐾| is greater than |𝑊| is the geographical distribution of demand points compared to warehouses, as it is shown in Figure 6. Having more clusters |𝐾| allows for greater similarity of demand points within the same cluster. Consequently, the number of clusters |𝐾| must be at least as large as the number of warehouse clusters |𝑊|. Furthermore, the first stage of disaggregation begins with demand zones, where the value of 𝑦𝑤𝑖𝑐𝑡(1)∗ is determined optimally. The accuracy of this stage significantly impacts subsequent stages, leading to enhanced overall performance. 49 Figure 5. Optimality gap versus cluster size |W| 10%12%14%16%18%20%22%24%26%2 10 20 50Optimality gapWarehouse cluster counts, |W|10-4-(*) 10-20-(*) 10-40-(*) 10-100-(*)50-4-(*) 50-20-(*) 50-40-(*) 50-100-(*)100-4-(*) 100-20-(*) 100-40-(*) 100-100-(*)200-4-(*) 200-20-(*) 200-40-(*) 200-100-(*)10.00%12.00%14.00%16.00%18.00%20.00%22.00%24.00%26.00%4 20 40 100Optimality gapDemand zone cluster counts, |K|10-(*)-2 10-(*)-10 10-(*)-20 10-(*)-50 50-(*)-2 50-(*)-1050-(*)-20 50-(*)-50 100-(*)-2 100-(*)-10 100-(*)-20 100-(*)-50200-(*)-2 200-(*)-10 200-(*)-20 200-(*)-50Figure 4. Optimality gap versus cluster size |K| 50 Figure 6. Distribution of demand points and fulfilment centres 6.3. Effects of cluster count on optimality gap and computation time This section investigates how the sizes of clusters (C, K, W) influence the trade-off between the optimality gap and computation time. Our results are depicted in Figure 7, where cluster sizes are expressed as percentages of the dataset. For example, a cluster size of 1% means that the clusters (𝐶, 𝐾,𝑊) are 1% of the dataset size, resulting in 𝐶 = 10, 𝐾 = 4,𝑊 = 2. The optimality gap is calculated as (𝒁𝑫𝒊𝒔𝒔𝒑−𝒁∗𝒐𝒑𝒕)𝒁∗𝒐𝒑𝒕. Each aggregation stage is computed using parallel processing, i.e., each cluster is distributed on a single processor. Larger cluster sizes result in a smaller optimality gap because they allow the aggregation model (6) to be a better approximation of the original model (1). The findings indicate that very small cluster sizes do not yield satisfactory bounds. However, for medium-sized clusters (above 10%), the cost of the solution from our algorithm approaches the cost of the optimal solution to the original problem (1). The 7% gap with a cluster size of 40% can be considered an appropriate solution compared to other cases and requires 43 minutes of computation. 51 Figure 7. Optimality gap and computation time versus cluster size On the other hand, the computation time increases with cluster size. For instance, computation time is about 5 minutes for cluster sizes up to 10%. Beyond this point, each additional 10% in cluster size roughly doubles the computation time. This is because larger cluster sizes make the computation time dominated by the aggregated model (3), requiring us to solve a larger integer programming problem, which grows exponentially in running time as dimensions increase. However, the disaggregation stages, particularly stages 2 and 3, do not significantly increase with larger cluster sizes. This is because not only we apply the greedy algorithm for the last two stages, which significantly improve our running time but also all disaggregation stages can be executed independently for each cluster 𝑐 ∈ 𝐶, 𝑘 ∈ 𝐾,𝑤 ∈ 𝑊. As a result, the computation time could be reduced with parallel processors. With 𝑛 parallel processors, we can speed up the computation roughly by a factor of 𝑛 by distributing clusters across processors. 6.4. Effects of item size variations on optimality gap and running time We conduct an experiment to see how the variations in the item dimensions affects the solution quality and running time. We change the amount of volume upper-bound within the range of (5,65), with the step size of 10. Therefore, the intervals span from [1, 5] to [1, 65]. For each interval, we 00.050.10.150.20.250.30204060801001201401% 2% 3% 4% 5% 6% 7% 8% 9% 10% 20% 30% 40% 50% 60%Optimality gapRunning time (Min)Cluster sizeAgg Diss_stage1 Diss_stage2Diss_stage3 Total Approximation ratio 52 solve both the optimal model and the algorithm with the configuration of (|𝑃|, |𝐼|, |𝑊|) =(400, 160, 80). The resulting performance of the algorithm is depicted in Figure 8. The findings indicate that a greater dispersion in item sizes corresponds to a slightly higher optimality gap. The optimality gap increases by approximately 0.1%, i.e., from 0.0791 for 𝑣 = [1, 5] to 0.08 for 𝑣 = [1, 65]. The increase in optimality gap comes from the way we calculate the aggregated volume, defined as ?̅?𝑐 = 𝑚𝑎𝑥𝑝∈𝑃𝑐{𝑣𝑝}. Since the aggregated model uses the maximum size volume within each cluster 𝑐 ∈ 𝐶, it becomes overly conservative, causing the solution obtained by the reduced model to deviate significantly from the optimal solution. This demonstrates how the quality of the reduced model’s solution influences the subsequent stages and the final solution. Furthermore, this result aligns with the optimality gap formulation (17h) caused by disaggregation discussed in the previous section. As expected, increasing disparity in item volumes leads to a higher value of ?̅?𝑐 of each cluster 𝑐 ∈ 𝐶, which in turn increases the optimality gap. Furthermore, the computation time of the algorithm remains almost constant while the optimal model’s computation time is more sensitive but without a clear trend, indicating the robustness of our algorithm in terms of running time. This can be seen from the fact that we employ the Greedy method to handle two disaggregation stages which is less sensitive to the volume fluctuations. 53 Figure 8. The effect of sparsity in item sizes on optimality gap and running time 6.5. Effects of objective function parameter ratios on algorithm performance We conducted an experiment to evaluate the sensitivity of the algorithm's performance to the parameters of the objective function, including the unit costs of receiving, holding, and shipment. The algorithm is tested three times, each time normalizing a different parameter. For instance, while normalizing 𝑟𝑝𝑓𝑡 in the first experiment, we examine the effects of change in the ratios (ℎ𝑝𝑓𝑡𝑟𝑝𝑓𝑡) and (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡). Then, the ratios are adjusted by factors of 2 and 1/2. The results are summarized in Table 6. The results indicate the optimality gap has a negative relationship with the ratio of (ℎ𝑝𝑓𝑡𝑟𝑝𝑓𝑡) and (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡). In other words, halving the ratio almost doubles the optimality gap, while doubling the ratio significantly reduces the gap. However, by normalizing ℎ𝑝𝑓𝑡 in the second experiment, varying the ratios has a minimal effect on the optimality gap. The changes are slight, indicating a weak relationship between the ratio of (𝑟𝑝𝑓𝑡ℎ𝑝𝑓𝑡) , (𝛽𝑓𝑖𝑝ℎ𝑝𝑓𝑡) and the optimality gap. For the last experiment, the relationship is opposite to that observed for the first one. The optimality gap 0.07840000.07880000.07920000.07960000.08000000.080400000.511.522.533.544.5v=[1,5] v=[1,15] v=[1,25] v=[1,35] v=[1,45] v=[1,55] v=[1,65]Optimality GapRunning time (Hours)Item Size variationsAlgorithm running time optimal model running time Approximation ratio 54 increases as the ratio increases, showing a positive relationship. Halving the ratios (𝑟𝑝𝑓𝑡𝛽𝑓𝑖𝑝) and (ℎ𝑝𝑓𝑡𝛽𝑓𝑖𝑝) results in the smallest optimality gap, while doubling it leads to the largest gap. The results show that the ratio (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡) has a more significant impact on the optimality gap compared to the ratio (𝑟𝑝𝑓𝑡ℎ𝑝𝑓𝑡) and (𝛽𝑓𝑖𝑝ℎ𝑝𝑓𝑡). This highlights the importance of accurately computing the parameters 𝛽𝑓𝑖𝑝 and 𝑟𝑝𝑓𝑡. If the ratio (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡) is underestimated, the optimality gap may worsen. In addition, efforts to reduce the unit receiving cost, 𝑟𝑝𝑓𝑡, not only lower the total costs, but it also increases the ratio (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡), potentially leading to a smaller optimality gap. This also is consistent with the optimality gap formulation (17ℎ) previously discussed, where the optimality gap is a function of 𝑟𝑝𝑓𝑡 and 𝛽𝑓𝑖𝑝. Therefore, the algorithm is more sensitive to the ratio (𝑟𝑝𝑓𝑡𝛽𝑓𝑖𝑝) than to (𝛽𝑓𝑖𝑝ℎ𝑝𝑓𝑡) or (𝑟𝑝𝑓𝑡ℎ𝑝𝑓𝑡). Table 6. Sensitivity analysis on the ratio of objective function parameters # Unit Cost Ratio Optimality Gap 1 𝑟𝑝𝑓𝑡 (ℎ𝑝𝑓𝑡𝑟𝑝𝑓𝑡) , (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡) *(0.5) 0.140331377 *(1) 0.079754131 *(2) 0.043944369 2 ℎ𝑝𝑓𝑡 (𝑟𝑝𝑓𝑡ℎ𝑝𝑓𝑡) , (𝛽𝑓𝑖𝑝ℎ𝑝𝑓𝑡) *(0.5) 0.082758264 *(1) 0.079754131 *(2) 0.077929631 3 𝛽𝑓𝑖𝑝 (𝑟𝑝𝑓𝑡𝛽𝑓𝑖𝑝) , (ℎ𝑝𝑓𝑡𝛽𝑓𝑖𝑝) *(0.5) 0.042020462 *(1) 0.079754131 *(2) 0.143236938 6.6. Comparison with the greedy algorithm In this section, we compare our proposed framework with the greedy algorithm when used as the primary method to solve model (1). The experiment aims to evaluate the impact of scaling up the problem size on both the greedy algorithm and our proposed framework in terms of running time and performance. 55 To scale up the problem, we increase the item counts from 1000 to 20,000. We evaluate the effect of scaling up on algorithm performance by considering four cluster sizes (10%, 20%, 30%, 40%). The performance is assessed by computing the ratio of the algorithm solution 𝑍𝐷𝑖𝑠𝑠𝑝 to the greedy solution 𝑍𝐺𝑜𝑝𝑡, i.e., (𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡) for each cluster. The average ratio across the four clusters is then calculated for each item size. The results of this experiment are depicted in Figure 9. The results indicate that while the greedy algorithm is appropriate for handling small to medium-sized problems, it struggles with large-scale problems. As the number of items increases, the greedy algorithm's running time grows substantially compared to our algorithm, and this growth rate becomes more pronounced with larger problem sizes. For example, when increasing the number of items from 2,000 to 5,000, the slope of the running time graph for the greedy algorithm is 0.03. However, this slope increases to 0.05 when scaling from 5,000 to 10,000 items, indicating an accelerating growth rate. In contrast, our proposed algorithm shows only a modest increase in running time for each cluster size. Specifically, the running time for a 40% cluster size grows by approximately 89% (from 43 minutes with 1000 items to 81 minutes with 10,000 items), whereas the greedy algorithm's running time increases by 1700% for the same item range. To quantify the relative performance of our algorithm, we can define the ratio 𝑍𝐷𝑖𝑠𝑠𝑝𝑍∗𝑜𝑝𝑡=𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡×𝑍𝐺𝑜𝑝𝑡𝑍∗𝑜𝑝𝑡, where the first part 𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡 reflects the approximation error due to aggregation, and the second part 𝑍𝐺𝑜𝑝𝑡𝑍∗𝑜𝑝𝑡 represents the approximation error of the greedy algorithm itself. This experiment specifically evaluates the value of 𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡. The results show that the average ratio 𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡 increases by 4% as the problem size scales from 1,000 to 10,000 items. However, as we scale up the problem size, we observe a decreasing growth rate in the average ratio 𝑍𝐷𝑖𝑠𝑠𝑝𝑍𝐺𝑜𝑝𝑡, which emphasizes our algorithm’s stable performance for handling large-scale models. This stability implies that our proposed framework maintains a manageable optimality gap, thereby providing a more reliable and scalable approach against greedy for large-scale problems. 56 Figure 9. Comparison of the algorithm with greedy performance 6.7. Comparison with the linear programming relaxation solution In this section, we aim to compare our algorithm with the LP relaxation model. In Section 5.4.1, we demonstrated that the solution obtained by our algorithm ensures feasibility. Now, we want to evaluate whether we can relax the integrity constraint (1f) and solve model (1) by rounding linear programming solution. The first step is to determine if the solution obtained by rounding the linear relaxation model is feasible for the online-retail inventory management problem (1). To test this, we solve model (1) by relaxing the integrality constraint (1f), and then rounding off the optimal value of each decision variable to the nearest integer point. Subsequently, we construct a new linear programming model to verify the feasibility of this rounded solution for the model (1). We define a new set of variables (𝜀, 𝜖+, 𝜖−, 𝜏, 𝜑) corresponding to each constraint (1b) to (1e) to measure the degree of violation. The model (18) is formulated as follows: 1.051.061.071.081.091.11.111.121.130501001502002503003504004501000 2000 5000 10000Ratio Running time (Min)Number of itemsGreedy 100-40-20 200-80-40 300-120-60 400-160-80 Average Ratio 57 𝑍𝐹𝑆 =𝑀𝑖𝑛𝜀, 𝜖+, 𝜖−, 𝜏, 𝜑 (∑Γ𝑝𝑖𝑡 ∗ 𝜀𝑝𝑖𝑡𝑝,𝑖,𝑡+ ∑ Π𝑝𝑓𝑡 ∗ (𝜖𝑝𝑓𝑡+ + 𝜖𝑝𝑓𝑡− )𝑝,𝑓,𝑡+∑Υ𝑝𝑡 ∗ 𝜏𝑝𝑡𝑝,𝑡+∑Ξ𝑓𝑡 ∗ 𝜑𝑓𝑡𝑓,𝑡) (18𝑎) Subject to: ∑?̿?𝑓𝑖𝑝𝑡𝑓+ 𝜀𝑝𝑖𝑡 ≥ 𝑑𝑝𝑖𝑡 ∀ 𝑝, 𝑖, 𝑡 (18𝑏) ?̿?𝑝𝑓𝑡 − ?̿?𝑝𝑓,𝑡−1 − ?̿?𝑝𝑓𝑡 +∑?̿?𝑓𝑖𝑝𝑡𝑖= 𝜖𝑝𝑓𝑡+ − 𝜖𝑝𝑓𝑡− ∀ 𝑝, 𝑓, 𝑡 (18𝑐) ∑?̿?𝑝𝑓𝑡𝑓− 𝜏𝑝𝑡 ≤ 𝛼𝑝𝑡 ∀ 𝑝, 𝑡 (18𝑑) ∑𝑣𝑝?̿?𝑝𝑓,𝑡−1𝑝+∑𝑣𝑝?̿?𝑝𝑓𝑡𝑝− 𝜑𝑓𝑡 ≤ 𝜛𝑓𝑡 ∀ 𝑓, 𝑡 (18𝑒) 𝜀𝑝𝑖𝑡, 𝜖𝑝𝑓𝑡+ , 𝜖𝑝𝑓𝑡− , 𝜏𝑝𝑡, 𝜑𝑓𝑡 ∈ 𝑅+ (18𝑓) Where (Γ𝑝𝑖𝑡, Π𝑝𝑓𝑡, Υ𝑝𝑡, Ξ𝑓𝑡) denote the penalty units for possible violations of constraint (18b) to (18e). The values (?̿?𝑝𝑓𝑡, ?̿?𝑓𝑖𝑝𝑡, ?̿?𝑝𝑓𝑡) represent the rounded solution derived from the linear relaxation model. If the optimal objective function value, 𝑍𝐹𝑆∗ , equals zero, i.e., 𝑍𝐹𝑆∗ = 0, the solution is feasible for model (1); otherwise, it is not. For the particular instance described in Section 6.1, the solution from linear relaxation Z𝑂𝑃𝑇𝑅 is equal to the optimal model, i.e., Z𝑅∗ = Z𝑂𝑃𝑇∗ and 𝑍𝐹𝑆∗ = 0. However, this is not typical for most instances. According to the findings in Table 7, the solution obtained from LP relaxation is not even feasible for most cases, even though it has a lower optimality gap compared to our algorithm. Moreover, even if the solution is feasible, LP relaxation does not resolve the issue of problem’s scale; the problem still remains large-scale. For instance, solving the LP relaxation of model (1) for our specific instance takes 2.11 hours, whereas our algorithm can handle it in a few minutes, demonstrating superior computational efficiency compared to LP relaxation. 58 Table 7. Comparison with LP relaxation solution # Unit Demand Linear relaxation objective function value (𝐙𝑹∗ ) Optimal model objective function value (𝐙𝑶𝑷𝑻∗ ) 𝐙𝑭𝑺∗ 1 ~ 𝑈[0,5] 41,822,762.46 41,823,385.5 536.1143 2 ~ 𝑈[0,10] 94,085,771.75 94,087,700.5 1130.641 3 ~ 𝑈[0,50] 513,028,291.2 513,028,292 671.1964 4 ~ 𝑈[0,100] 1,041,416,914 1,041,416,917 167.799109 59 7. Conclusion This study focuses on optimizing inventory management for an online retailer, aiming to minimize total costs associated with receiving, holding, and shipping. We formulate the problem as a multi-period integer programming model. Given the complexities posed by its large-scale and NP-hard nature, we propose an approximation algorithm. The contributions of this study can be summarized as follows: • The formulation is general enough which can capture different types of variable costs and constraints. One crucial assumption addressed in our problem is the variation in item sizes, enhancing the accuracy of inventory planning while also adding to the problem's complexity. • We propose a large-scale framework for the problem, which integrates three sets of items, warehouses, and customers as a technique for reducing the problem dimension. The proposed algorithm decomposes the problem into different tractable stages and leverages parallel computing. Each disaggregation stage can be computed independently for each cluster, and these clusters can be distributed across processors. To address the NP-hard nature of the last two stages of disaggregation, we develop a greedy algorithm capable of solving the model within polynomial time. The combination of the parallel computing and the Greedy algorithm leads to a huge reduction in running time. The effectiveness of the algorithm is demonstrated through numerical examples. • We develop a theoretical optimality gap bound that address the effect of applying aggregation on handling the large-scale integer program. The bound provides insights into the algorithm’s performance and indicates how far the algorithm can deviate from the optimal model. Our empirical analysis reveals the impact of several factors on the solution quality, including the number of clusters, the item characteristics, and variable costs. These insights can be used to guide the choice of parameters when designing the inventory management for online-retail businesses using the proposed large-scale algorithm. The results indicate that the algorithm can find near optimal solutions within a few minutes for a medium number of clusters. Additionally, a greater dispersion in item sizes corresponds to a slightly higher optimality gap. The sensitivity 60 analysis on the ratio of objective function parameters highlights the crucial role of the ratio of unit receiving and shipment costs, (𝛽𝑓𝑖𝑝𝑟𝑝𝑓𝑡), in controlling the optimality gap. We also analyzed the effectiveness of our algorithm to that of the greedy algorithm when the latter is used as the primary method to solve model (1). Although the greedy algorithm offers a lower optimality gap, it faces significant challenges with large instances. Specifically, the computational time for the greedy approach exhibits an accelerating growth rate that becomes increasingly pronounced with larger problem sizes. In contrast, our algorithm demonstrates superior scalability and maintains stable performance across larger datasets. This is achieved through breaking down the problem into four manageable stages and using greedy as a subroutine, which allows our algorithm to manage large datasets effectively without significant reduction in solution quality. Thus, for large-scale inventory management problems, our algorithm proves to be more robust and efficient than the greedy algorithm. Moreover, our comparison with LP relaxation reveals that recovering the optimal solution from LP relaxation is often impossible. Rounding the optimal values from LP relaxation often fails to ensure feasibility and, in rare cases where feasibility is achieved, the computation time remains significant. Therefore, an efficient algorithm is required to handle large-scale instances within a reasonable time frame. This underscores the importance of our algorithm’s computational efficiency in addressing these challenges. Although the proposed large-scale framework is a powerful tool for online retail inventory management, it has several limitations: • The model does not account for shortages, whether in the form of backorders or lost sales. • Uncertainty in parameters such as demand, which affects the solution quality, is not considered. • The model simplifies the network to include only fulfillment centers and customers, excluding intermediate centers such as sortation centers and delivery stations between warehouses and customers. 61 For future work, one promising direction is to develop a customized clustering method for our specific problem. The findings indicate that standard algorithms, such as K-means, may not be an ideal technique for our case. This is because the algorithm does not give us satisfactory solutions for small size of clusters. We need to develop a customized clustering method that effectively group together items, warehouses and demand points in a way that makes the reduced model be a better approximation of the optimal model. Another approach to improve the aggregation process is to employ the iterative aggregation schemes, similar to those applied in the network design problem (Bärmann et al., 2015). Incorporating uncertainty into parameters such as demand represents another important area for future work. Current model assumes deterministic demand, which may not always produce accurate solutions, particularly when historical data is insufficient or unreliable. Future research could focus on developing data-driven inventory optimization strategies that account for uncertainty, thereby the solution will be better aligned with real-world complexities. 62 References Abouelrous, A., Gabor, A. F., & Zhang, Y. (2022). Optimizing the inventory and fulfillment of an omnichannel retailer: A stochastic approach with scenario clustering. Computers & Industrial Engineering, 173, 108723. https://doi.org/10.1016/j.cie.2022.108723 Amazon. (2022, December 19). List of Amazon Fulfillment Center Locations [2024 Updated]. https://amzprep.com/fba-locations/ Amazon. (2024). Our Facilities. US About Amazon. https://www.aboutamazon.com/workplace/facilities Amazon Distribution Network Strategy | MWPVL International. (n.d.). Retrieved 8 March 2024, from https://mwpvl.com/html/amazon_com.html Arlotto, A., Keskin, I. N., & Wei, Y. (2023). Online Demand Fulfillment Problem with Initial Inventory Placement: A Regret Analysis (SSRN Scholarly Paper 4666493). https://doi.org/10.2139/ssrn.4666493 Bärmann, A., Liers, F., Martin, A., Merkert, M., Thurner, C., & Weninger, D. (2015). Solving network design problems via iterative aggregation. Mathematical Programming Computation, 7(2), 189–217. https://doi.org/10.1007/s12532-015-0079-1 Bretthauer, K. M., Mahar, S., & Venakataramanan, M. A. (2010). Inventory and distribution strategies for retail/e-tail organizations. Computers & Industrial Engineering, 58(1), 119–132. https://doi.org/10.1016/j.cie.2009.09.005 Chen, A. I., & Graves, S. C. (2021). Item Aggregation and Column Generation for Online-Retail Inventory Placement. Manufacturing & Service Operations Management, 23(5), 1062–1076. https://doi.org/10.1287/msom.2020.0867 Chen, C. D. (2018). Operations management in a large online retailer: Inventory, scheduling and picking [Thesis, Massachusetts Institute of Technology]. https://dspace.mit.edu/handle/1721.1/120189 Cohen, R., Katzir, L., & Raz, D. (2006). An efficient approximation for the Generalized Assignment Problem. Information Processing Letters, 100(4), 162–166. https://doi.org/10.1016/j.ipl.2006.06.003 Conforti, M., Cornuéjols, G., & Zambelli, G. (2014a). Integer Programming Models. In M. Conforti, G. Cornuéjols, & G. Zambelli (Eds.), Integer Programming (pp. 45–84). Springer International Publishing. https://doi.org/10.1007/978-3-319-11008-0_2 Conforti, M., Cornuéjols, G., & Zambelli, G. (2014b). Integer Programming Models. In M. Conforti, G. Cornuéjols, & G. Zambelli (Eds.), Integer Programming (pp. 45–84). Springer International Publishing. https://doi.org/10.1007/978-3-319-11008-0_2 Dell’Amico, M., Delorme, M., Iori, M., & Martello, S. (2019). Mathematical models and decomposition methods for the multiple knapsack problem. European Journal of Operational Research, 274(3), 886–899. https://doi.org/10.1016/j.ejor.2018.10.043 DeValve, L., Wei, Y., Wu, D., & Yuan, R. (2021). Understanding the Value of Fulfillment Flexibility in an Online Retailing Environment. Manufacturing & Service Operations Management. https://doi.org/10.1287/msom.2021.0981 Evans, J. R. (1979). Aggregation in the generalized transportation problem. Computers & Operations Research, 6(4), 199–204. https://doi.org/10.1016/0305-0548(79)90003-0 Fathollahi-Fard, A. M., Wong, K. Y., & Aljuaid, M. (2023). An efficient adaptive large neighborhood search algorithm based on heuristics and reformulations for the generalized 63 quadratic assignment problem. Engineering Applications of Artificial Intelligence, 126, 106802. https://doi.org/10.1016/j.engappai.2023.106802 Ge, D., Pan, Y., Shen, Z.-J. (Max), Wu, D., Yuan, R., & Zhang, C. (2019). Retail supply chain management: A review of theories and practices. Journal of Data, Information and Management, 1(1), 45–64. https://doi.org/10.1007/s42488-019-00004-z Geoffrion, A. M. (1976). Customer Aggregation in Distribution Modeling,. https://apps.dtic.mil/sti/citations/ADA033973 Govindarajan, A., Sinha, A., & Uichanco, J. (2021). Joint inventory and fulfillment decisions for omnichannel retail networks. Naval Research Logistics (NRL), 68(6), 779–794. https://doi.org/10.1002/nav.21969 Hallefjord, Å., Jörnsten, K. O., & Värbrand, P. (1993). Solving large scale generalized assignment problems—An aggregation / disaggregation approach. European Journal of Operational Research, 64(1), 103–114. https://doi.org/10.1016/0377-2217(93)90011-B Hallefjord, A., & Storoey, S. (1986). Aggregation and disaggregation in integer programming problems. https://www.osti.gov/etdeweb/biblio/7768114 Hitchcock, F. L. (1941). The Distribution of a Product from Several Sources to Numerous Localities. Journal of Mathematics and Physics, 20(1–4), 224–230. https://doi.org/10.1002/sapm1941201224 Hoffman, A. J., Kruskal, J. B., DANTZIG, G. B., DUFFIN, R. J., FAN, K., FORD, L. R., FULKERSON, D. R., GALE, D., GOLDMAN, A. J., HELLER, I., KRUSKAL, J. B., KUHN, H. W., MILLS, H. D., THOMPSON, G. L., TOMPKINS, C. B., TUCKER, A. W., & WOLFE, P. (1956). INTEGRAL BOUNDARY POINTS OF CONVEX POLYHEDRA. In H. W. Kuhn & A. W. Tucker (Eds.), Linear Inequalities and Related Systems. (AM-38) (pp. 223–246). Princeton University Press; JSTOR. http://www.jstor.org/stable/j.ctt1b9x27g.16 Hübner, A. H., Kuhn, H., & Sternbeck, M. G. (2013). Demand and supply chain planning in grocery retail: An operations planning framework. International Journal of Retail & Distribution Management, 41(7), 512–530. https://doi.org/10.1108/IJRDM-05-2013-0104 Jiang, D.-P., & Li, X. (2020). Order fulfilment problem with time windows and synchronisation arising in the online retailing. International Journal of Production Research, 59, 1–29. https://doi.org/10.1080/00207543.2020.1721589 Kazemi, A. (2021). New Aggregation Methods for Multicommodity Network Flow Problems: Theory and Applications to Locomotive Refueling [Doctoral dissertation]. Monash University, Australia. Kellerer, H., Pferschy, U., & Pisinger, D. (2004). Multidimensional Knapsack Problems. In H. Kellerer, U. Pferschy, & D. Pisinger (Eds.), Knapsack Problems (pp. 235–283). Springer. https://doi.org/10.1007/978-3-540-24777-7_9 Leisten, R. (1998). An LP-aggregation view on aggregationin multi-level production planning. Annals of Operations Research, 82(0), 233–250. https://doi.org/10.1023/A:1018931224060 Li, X., Li, J., Aneja, Y. P., Guo, Z., & Tian, P. (2019). Integrated order allocation and order routing problem for e-order fulfillment. IISE Transactions, 51(10), 1128–1150. https://doi.org/10.1080/24725854.2018.1552820 Litvinchev, I., & Tsurkov, V. (2003). Aggregated Problem and Bounds for Aggregation. In I. Litvinchev & V. Tsurkov (Eds.), Aggregation in Large-Scale Optimization (pp. 1–60). Springer US. https://doi.org/10.1007/978-1-4419-9154-6_1 64 Martello, S., & Toth, P. (1990). Knapsack problems: Algorithms and computer implementations. John Wiley & Sons, Inc. Mendelssohn, R. (1980). Technical Note—Improved Bounds for Aggregated Linear Programs. Operations Research, 28(6), 1450–1453. https://doi.org/10.1287/opre.28.6.1450 Rodrigue, J.-P. (2020). The distribution network of Amazon and the footprint of freight digitalization. Journal of Transport Geography, 88, 102825. https://doi.org/10.1016/j.jtrangeo.2020.102825 Rogers, D. F., Plante, R. D., Wong, R. T., & Evans, J. R. (1991). Aggregation and Disaggregation Techniques and Methodology in Optimization. Operations Research, 39(4), 553–582. https://doi.org/10.1287/opre.39.4.553 Saeedi, M., Parhazeh, S., Tavakkoli-Moghaddam, R., & Khalili-Fard, A. (2024). Designing a two-stage model for a sustainable closed-loop electric vehicle battery supply chain network: A scenario-based stochastic programming approach. Computers & Industrial Engineering, 190, 110036. https://doi.org/10.1016/j.cie.2024.110036 Salam, M., Panahifar, F., & Byrne, P. (2016). Retail supply chain service levels: The role of inventory storage. Journal of Enterprise Information Management, 29, 887–902. https://doi.org/10.1108/JEIM-01-2015-0008 Sarin, S. C., Sherali, H. D., & Kim, S. K. (2014). A branch-and-price approach for the stochastic generalized assignment problem. Naval Research Logistics (NRL), 61(2), 131–143. https://doi.org/10.1002/nav.21571 Seyedan, M., Mafakheri, F., & Wang, C. (2022). Cluster-based demand forecasting using Bayesian model averaging: An ensemble learning approach. Decision Analytics Journal, 3, 100033. https://doi.org/10.1016/j.dajour.2022.100033 Shetty, C. M., & Taylor, R. W. (1987). Solving large-scale linear programs by aggregation. Computers & Operations Research, 14(5), 385–393. https://doi.org/10.1016/0305-0548(87)90035-9 Torabi, S. A., Hassini, E., & Jeihoonian, M. (2015). Fulfillment source allocation, inventory transshipment, and customer order transfer in e-tailing. Transportation Research Part E: Logistics and Transportation Review, 79, 128–144. https://doi.org/10.1016/j.tre.2015.04.004 U.S. Census Bureau. (2023). QUARTERLY RETAIL E-COMMERCE SALES 4th QUARTER 2023. https://www.census.gov/retail/mrts/www/data/pdf/ec_current.pdf Wang, K., Li, Y., & Zhou, Y. (2022). Execution of Omni-Channel Retailing Based on a Practical Order Fulfillment Policy. Journal of Theoretical and Applied Electronic Commerce Research, 17(3), Article 3. https://doi.org/10.3390/jtaer17030060 Woodcock, A. J., & Wilson, J. M. (2010). A hybrid tabu search/branch & bound approach to solving the generalized assignment problem. European Journal of Operational Research, 207(2), 566–578. https://doi.org/10.1016/j.ejor.2010.05.007 Zipkin, P. H. (1980a). Bounds for aggregating nodes in network problems. Mathematical Programming, 19(1), 155–177. https://doi.org/10.1007/BF01581638 Zipkin, P. H. (1980b). Bounds for Row-Aggregation in Linear Programming. Operations Research. https://doi.org/10.1287/opre.28.4.903 Zipkin, P. H. (1980c). Bounds on the Effect of Aggregating Variables in Linear Programs. Operations Research, 28(2), 403–418. https://doi.org/10.1287/opre.28.2.403
Thesis/Dissertation
eng
Business Administration - Transportation and Logistics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Cohen-Hillel, Tamar
Inventory management optimization : a case study of Amazon supply chain
Text