Based on the result from the previous step, fix certain
elements of \(S\)  to be 0 and solve the modified problem:

\(\min \ E\left[ {\sum \limits _{t = 1}^T {{C_t}\left( \tilde{S} \right)} } \right]\quad {\rm where}\ \tilde{S} = S \odot {1_{\left\lbrace  {S_1^* > 0} \right\rbrace }},\) 

where \(\odot \)  denotes element-wise product.
In Stage 1, our goal is to select a small fraction of the nodes to keep inventory. In Stage 2, we only consider the selected inventory locations and solve a much smaller-scale smooth inventory optimization problems to determine their base-stock levels. We use the stochastic FISTA in Stage 1 and the SGD in Stage 2.
Notice that solving the regularized problem in two stages is not new. In the field of high dimensional regression, [1]} proposed a two-stage procedure, termed the relaxed Lasso for high-dimensional data where the number of predictor variables \(p\)  is much larger than the number of observations \(n\) . Their theoretical and numerical results demonstrate that the two-stage procedure performs better than the regular Lasso estimator for high-dimensional data.
Numerical Experiments
In this section, we first test the performances of the simulation and gradient computation algorithms on inventory systems of different scales. Then, we conduct a series of experiments to understand the behaviors of the optimization algorithm and to compare it with the GS model of [2]} on small- to medium-scale problems where the GS model applies. In this section, all the computer programs are coded in Python and all the experiments are run on a computer with two Intel Xeon Gold 6248R CPUs (each with 24 cores) and 256GB RAM.
Simulation and Gradient Computation 
To test the performance of our algorithms, We build test inventory examples based on the BOM characteristics that we observe from our consulting experience. We set the average degrees \(\left\langle k \right\rangle \)  of all BOM networks as 10 and the time horizons as 100 periods.
Performance of the Simulation Algorithms.
We test the traditional simulation algorithm, the tensorized simulation algorithms (i.e., Algorithm REF ) with dense matrices or sparse matrices for inventory systems with different number of nodes, ranging from medium-scale system with 1,000 nodes to very large-scale system with 500,000 nodes, and report the simulation run times of different algorithms, averaging over 10 independent replications, in Table REF , where “-” indicates that the algorithm takes over a day so it is terminated before completion and “/” indicates that there is not enough memory in our computer to run the algorithm.
<TABLE>There are several interesting findings from the results. First, the computational complexities of the traditional algorithm and the tensorized algorithm with dense matrices are approximately \(O(n^2)\) , while that of the tensorized algorithm with sparse matrices is approximately \(O(n)\) , which are consistent with Theorem REF . Second, even though tensorization does not improve the computational complexity, it reduces the computational time significantly. Third, the use of sparse matrices reduces significantly both the memory requirements and the run times, further allowing the algorithm to handle very large-scale problems.

Performance of the Gradient Computation Algorithms.
We use the same test problems to test the two gradient computation algorithms proposed in this paper, the IPA algorithm (Algorithm REF ) and the BP Algorithm (Algorithms REF  and REF ). Notice that both algorithms are tensorized and both may use dense or sparse adjacency matrices. The run times of different algorithms are reported in Table REF .
<TABLE>From the results we see that the computational complexities of the IPA algorithms, with dense and sparse adjacency matrices, and the BP algorithms, with dense and sparse adjacency matrices, are approximately \({\rm O}\left( {{n^3}} \right)\) , \({\rm O}\left( {{n^2}} \right)\) , \({\rm O}\left( {{n^2}} \right)\)  and \({\rm O}\left( {{n}} \right)\)  respectively, which are consistent with Theorems REF  and REF . Furthermore, it is clear that the BP algorithms significantly outperform the IPA algorithms in general, and the BP algorithm with sparse adjacency matrices is capable of handling very large-scale problems.

Simulation and Gradient Computation Using TensorFlow.
One major impediment to the use of our algorithms in general inventory optimization is the construction of the computational graph and the derivation of the sample-path gradient. Small changes of the logic of the simulation process may result in quite different computational graphs and thus very different gradient estimators. Therefore, the use of these algorithms may require a significant amount of analyst's effort, which may be very expensive in practice. This motivates us to program the tensorized simulation algorithm using the machine-learning tool TensorFlow, which may construct the computational graph and compute the sample-path gradient automatically. Moreover, the use of TensorFlow also makes it possible to utilize GPUs for matrix operations used in the algorithms, which may provide further speedups.
We implement the tensorized simulation algorithm using TensorFlow, and also use its automatic BP tool to compute the sample-path gradient. We test both algorithms with dense or sparse adjacency matrices. We further test the algorithms using an additional NVIDIA GeForce RTX 3090 GPU with 24GB memory. The results are reported in Tables REF  and REF .
<TABLE><TABLE>There are several interesting findings from these results. First, when the inventory problems are of medium to large scales (e.g., 1,000 to 50,000 nodes), the TensorFlow implementations of the simulation algorithm and the automatic BP algorithm, both with sparse adjacency matrices, provide acceptable run times, making TensorFlow a very promising tool for inventory optimization. Second, it appears that the use of GPU does not provide any advantage when the sparse adjacency matrices are used, in terms of the computational time. Third, compared to Tables REF  and REF , it is clear that our original method with sparse adjacency matrices is significantly faster than the TensorFlow implementations, making it an ideal tool to solve large- to very large-scale problems, though it requires a significant amount of analyst's effort to construct the computational graph and to derive the BP gradient.

Performance of Multiple Replications.
To solve the inventory optimization problems, the objective value and the gradient need to be evaluated based on multiple replications (i.e., a mini-batch) in each epoch (i.e., iteration of the optimization algorithms). Therefore, it is important to understand the run times of the algorithms when multiple replications are run in parallel. We test the algorithms with a mini-batch of 10 and report the results in Tables REF  and REF .
<TABLE><TABLE>We employ the multiprocessing parallel computing technique and each replication is executed on a separate CPU core. As the communication and data transmission between the parent and children processes takes some time, the run time of multiple replications is higher than that of a single replication. However, as can be seen in the tables, the run time of a mini-batch of 10 is much lower than 10 times that of a single replication, suggesting that we should take advantage of the multi-core capability in the optimization process.

Performance of the Optimization Algorithms
Remember that in the Stage 1 of the optimization procedure, we solve a \(L_1\)  regularized stochastic optimization problem (REF ). Both the SSGD and FISTA algorithms may be used. We first compare the performance of these two algorithms on the inventory optimization problem with different scales. Notice that the subgradient of the \(L_1\)  norm (\({{{\left\Vert  S \right\Vert }_1}}\) ) is needed in the implementation of the SSGD algorithm and we use the sign function (\({\mathop {\rm sign}} (S)\) ). To provide a fair comparison, for both algorithms, we solve the same test problems used in Section REF , we apply the same mini-batch of 10 replications to evaluate the cost and gradient in each epoch, and we use the same regularization parameter (\(\lambda \) ) and the step-size (\({t_k})\) . Furthermore, because the two algorithms may stop at different epochs, we report the results based on fixed numbers of epochs.
The results are shown in Figure REF . The four rows represent four different problem scales, i.e., \(10,000\) , \(50,000\) , \(100,000\)  and \(500,000\)  nodes, and the two columns shows the objective values (i.e., the average total cost) and the sparsity of the solutions (i.e., the numbers of non-zero base-stock levels). From these plots, we see that, in terms of the objective values, both algorithms have similar performance and the FISTA algorithm finds slightly better solutions. In terms of the sparsity of the solutions, the FISTA algorithm tends to find solutions that are more sparse. Therefore, we conclude that in the Stage 1, the FISTA algorithm is a preferred one among the two.
<FIGURE>Notice that the solution of Stage 1 provides the information not only on where to keep inventory but also on the base-stock levels. However, in our optimization procedure, we only keep the location information and suggest to use Stage 2 optimization to optimize the base-stock levels. To test the advantage of a second stage, we conduct additional simulation experiments to estimate the objective values of the solutions of Stage 1 and Stage 2 and report the results in Table REF . From the results it is clear that the Stage 2 solution is clearly better than the Stage 1 solution for all four problems, validating the use of the Stage 2 to further improve the quality of the solution from Stage 1.
<TABLE>In Table REF  we report the computational times of a typical run of the two-stage optimization procedure, where the FISTA is used in Stage 1 and the SGD is used in Stage 2 and both algorithms are run until the termination conditions are satisfied. From the table, we see that a large-scale inventory optimization problem with 50,000 nodes may be solved in about 2 hours and a very large-scale problem with 500,000 nodes may be solved in about a day. Since inventory decisions are not real-time decisions and they are updated infrequently (for instance, once every six months), these run times are in general acceptable.
<TABLE>
Comparison with the Guaranteed Service Model
The simulation optimization algorithm appears to work well in Section REF . However, without a benchmark algorithm to compare with, it is difficult to know the actual performance of the algorithm. To partly solve this problem, we consider the GS model of [2]}. Notice that the GS model requires the BOM to be a spanning tree with bounded demands and deterministic lead times, which are difficult to satisfy for complex production systems. Nevertheless, it provides an optimization algorithm that can be used as a benchmark. In this subsection we consider three different types of test problems. The first one satisfies all requirements of the GS model, the second one has a more complex network structure, and the third one has random lead times. We compare our algorithm with the algorithm of [2]}. From these three types of test problems, we conclude that when the assumptions of the GS model are satisfied, both algorithms produce solutions with similar objectives and solution sparsity. However, when the assumptions are not satisfied, our algorithm performs significantly better than the GS algorithm, as expected.
Throughout this subsection we assume that the outside demands of all items are independent and normally distributed with different means and variances. To implement the GS model, we need an upper bound for each demand and we set it as an upper quantile of the normal distribution as suggested by [2]}. Moreover, the GS model quotes a committed service time to each external demand. We revise the logic of our simulation model slightly to adopt the same logic to ensure a more fair comparison.
Spanning Tree.
We start with a simple example of Kodak digital camera supply chain pictured in Figure REF  [6]}. The lead times and holding costs are listed in Table REF .
<FIGURE><TABLE>We apply both algorithms and the optimal base-stock levels found by both algorithms are shown in Table REF . As can be seen from the table, both algorithms produce solutions with similar costs and sparsity, with the GS algorithm performs slightly better than our algorithm.
<TABLE>To understand the performances of both algorithms for larger-scale problems, we built some larger inventory systems with spanning tree topology. The results of both algorithms are reported in Table REF . Again, we see that both algorithms produce solutions with similar costs and sparsity. Furthermore, we notice that, when the number of nodes is large, e.g., 50,000, our algorithm is significantly faster.
<TABLE>
Network with Cycles.
The GS model proposed by [2]} is developed for supply chains that can be modeled as spanning trees. For general systems which may include cycles, this model is no longer applicable. Next, we modify the Kodak network (see Figure REF ) to explore the influence of the cycles.
<FIGURE>To apply the GS model to the network in Figure REF , a simple idea is to cut off an edge and regard it as a spanning tree. Here, we cut off the arc between “D” and “J” so that the optimal base-stock level is the same as that of Table REF . We also obtain the optimal base-stock level by our algorithm, and the results are shown in Table REF . From the results we can see that our algorithm finds the solution with two more non-zero base-stock levels and better cost, compared to the solution of the GS model.
<TABLE>
Random Lead Times.
The GS model assumes a deterministic lead time for each node in the supply chain. Nevertheless, this is typically not true in practice. For instance, procurement lead times are often quite volatile. We consider a case that the procurement lead times of the raw material nodes in Figure REF  are random, i.e., \({l_{t,i}} = \max \left( {1,round\left( {N\left( {{\mu _{{l_i}}},0.05{\mu _{{l_i}}}} \right)} \right)} \right)\) . For the GS model, we use \(\mu _{l_i}\)  as the (deterministic) lead time. The optimal solutions found by both algorithms are evaluated based on the simulation model, and the total cost of 100 periods for the GS model and our algorithm are \(1.97 \times {10^6}\)  and \(9.38 \times {10^5}\) , respectively. It is clear that our algorithm performs significantly better than the GS model.

Concluding Remarks
Inspired by the similarities between inventory optimization and recurrent neural networks, in this paper we develop a new framework for inventory simulation and inventory optimization for large-scale inventory optimization problems on complex production networks. The framework utilizes modern computational tools to reduce the computational time drastically, so that large-scale problems with tens to hundreds of thousands nodes may be solved in a reasonable amount of time.
There are a few directions that one may further extend our work. First, many practical inventory problems have constraints, e.g., capacity constraints or fill-rate constraints. It is interesting to consider how to incorporate these constraints in the framework. Second, the similarity between RNN and dynamic simulation holds not only for inventory simulation, but also for general periodically reviewed dynamic systems. It is interesting to further study how the framework developed in this paper may be extended to simulation and optimization of general dynamic systems.
Derivation of the Inventory Position (Equation REF )
By definition, the inventory position equals the on-hand inventory plus the on-order quantity minus the backorders (demand that have occurred but have not been satisfied), i.e.,
\(I{P_{t,i}} = I_{t-1,i} + O{O_{t,i}} - \left( {B_{t-1,i}^{out} + B_{t - 1,i}^{in}} \right),\) 
where \(O{O_{t,i}}\)  denotes the on-order quantity of item \(i\)  at period \(t\) , \({B_{t - 1,i}^{out}}\)  and \({B_{t - 1,i}^{in}}\)  denote the backlogged outside and internal demand of item \(i\)  at period \(t-1\) , respectively. The on-order quantity refers to the amount that is ordered but not yet received through either procurement or production, i.e.,
\(O{O_{t,i}} = \sum \limits _{u = 1}^{t - 1} {\left( {{O_{u,i}} - {P_{u,i}}} \right)}.\) 
The backlogged internal demand can be calculated by
\(B_{t,i}^{in} = \sum \limits _j {{a_{ij}}O_{t,j}^b}.\) 
Notice that, unlike outside demands where at least one of the on-hand inventory and outside backlog have to be zero (see Equations REF  and REF ), it is possible to have positive on-hand inventory and positive internal backlog at the same time, because the backlog of the downstream item may be caused by the shortage of some other components.
According to Equation (REF ), we have
\(I{P_{t,i}} = I_{t - 1,i}^0 - \sum \limits _j {{a_{ij}}{M_{t - 1,j}}} + O{O_{t,i}} - B_{t - 1,i}^{out} - B_{t - 1,i}^{in}.\) 
Based on Equations (REF ) and (REF ), we have
\(I_{t - 1,i}^0 - B_{t - 1,i}^{out} = {I_{t - 2,i}} + {P_{t - 1,i}} - B_{t - 2,i}^{out} - D_{t - 1,i}^{out}.\) 
Thus,
\(I{P_{t,i}} &=& {I_{t - 2,i}} + {P_{t - 1,i}} - B_{t - 2,i}^{out} - D_{t - 1,i}^{out} - \sum \limits _j {{a_{ij}}{M_{t - 1,j}}} + O{O_{t - 1,i}} + {O_{t - 1,i}} - {P_{t - 1,i}} - \sum \limits _j {{a_{ij}}O_{t - 1,j}^b} \\&=& {I_{t - 2,i}} + O{O_{t - 1,i}} - \left( {B_{t - 2,i}^{out} + B_{t - 2,i}^{in}} \right) + B_{t - 2,i}^{in} - D_{t - 1,i}^{out} + {O_{t - 1,i}} - \sum \limits _j {{a_{ij}}\left( {{M_{t - 1,j}} + O_{t - 1,j}^b} \right)} \\&=& I{P_{t - 1,i}} + {O_{t - 1,i}} - D_{t - 1,i}^{out} + \sum \limits _j {{a_{ij}}O_{t - 2,j}^b} - \sum \limits _j {{a_{ij}}\left( {{O_{t - 1,j}} + O_{t - 2,j}^b} \right)} \\&=& I{P_{t - 1,i}} + {O_{t - 1,i}} - D_{t - 1,i}^{out} - D_{t - 1,i}^{in}.\) 
Therefore, Equation (REF ) holds.
