{'fraction_non_alphanumeric': 0.03419327916738261, 'fraction_numerical': 0.007462137018581703, 'mean_word_length': 4.621636539257624, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rolling horizon policies for multi-stage stochastic assemble-to-order problems

Multi-stage stochastic programming is a powerful and flexible tool to tackle decision problems under uncertainty. Unfortunately, its application is made difficult by two fundamental issues: the need to know the probability distribution of the underlying risk factors and the explosion of the size of the scenario tree in a multistage setting. In order to ease the latter issue, we may shorten the planning horizon, but this can introduce end-of-horizon effects, resulting in a myopic behavior. A possible remedy is to augment the objective function with an additional term accounting for the value of the terminal state. We do so in the context of a multistage production planning problem within an assemble-to-order setting under demand uncertainty. In order to test the idea in a difficult environment, we introduce seasonality, bimodality, and correlations in the distribution of end items demand. We check the effect of adding a piecewise linear approximation of the value of terminal inventory of components to a two-stage model formulation by out-of-sample, rolling horizon simulations. We compare different scenario tree structures, built under limited information about the underlying distribution. Computational experiments suggest the potential benefit of adding a terminal value function, and illustrate interesting patterns arising from demand correlations and different levels of manufacturing capacity.


Introduction
We consider the application of multistage Stochastic Programming (SP) models with recourse to a production planning problem within an Assemble-to-Order (ATO) setting, whereby components (modules) are produced, and possibly stored in inventories, under capacity constraints and subject to uncertainty on the demand of end items. End items are only assembled for the current time period, after observing demand. In this setting, the power and flexibility of multistage SP models may play a significant role, but their application may be hindered by two issues: the lack of reliable information about probability distribution of the underlying risk factors and the exponential growth in the scenario tree.
A considerable research effort has been devoted to proper scenario tree generation in multistage SP. A key issue is the choice of the tree structure, in terms of branching factors and number of periods. In particular, limiting the number of periods may greatly reduce the size of the tree, with a potential adverse effect in terms of myopia. The nasty effects of truncated decision horizons is well-known even in a multiperiod deterministic setting. In [11] the addition of a terminal value is proposed to overcome end-of-horizon effects. In [19] issues in the truncation of infinite-horizon problems are thoroughly analyzed. The addition of a terminal state value is also a key ingredient in stochastic dynamic programming The ATO environment is interesting from both practical and theoretical points of view, because it is characterized by a proactive and dynamic approach to coping with uncertainty. If a single assembly period is considered, the problem can be formalized as a two-stage stochastic linear programming model [9], which can be tackled with a limited computational effort, but may be not quite realistic in some application settings. By considering multiple demand periods, a much more challenging multi-stage decision problem results, whereby different stages are linked by available components inventory. Actually, myopic policies based on two-stage models may perform well in a rolling horizon setting. For example, [13] shows that, if demands are identically distributed and independent (i.i.d.) across different time periods, the myopic policy based on a two-stage problem is optimal. On the contrary, a longer planning horizon is necessary when the i.i.d. assumption does not apply, possibly due to a local increase in demand (e.g., due to seasonality), which cannot be met by the available production capacity and requires careful planning of inventory buffers.
The aim of this paper is to investigate the performance of alternative model formulations, differing in tree structure and number of stages, in an ATO setting featuring complicated demand patterns, characterized by seasonality, bimodality, and possible correlations. We also consider the impact of available capacity for component production. In particular, we consider three types of models: a plain two-stage linear programming problem, a two-stage problem augmented with an estimate of the value of residual inventory, and a set of multi-stage problems characterized by different type of scenario trees. In order to compare the different different model structures in a realistic environment, we simulate the production system within a rolling horizon framework. Simulations are carried out out-of-sample, whereas the in-sample scenarios used in the model rely on limited information about the underlying demand distribution.
For getting the optimal solution of all the models, we use a standard commercial solver. In fact, we focus on a comparison of alternative modeling approaches and the related performance, rather than on ad hoc solution methods. Since the problem is characterized by several parameters and an exhaustive set of experiments would require too much space, we present in this paper only the most interesting results that we have obtained; nevertheless, we share the code 2 to enable interested readers to carry out further experiments.
The remaining of the paper is organized as follows: in Section 2 we review the literature about ATO problems and rolling horizon policies. In Section 3 we describe the mathematical models considered. Then, in Section 4 we introduce the general experimental setting, and in Section 5 we present the results of several computational experiments. Finally, Section 6 is devoted to conclusions and future research directions.

Literature review and paper positioning
In this paper, we consider finite-capacity ATO production planning problems, with single or multiple assembly periods and multiple end items. The bills of materials are flat and include only two levels: end items and components (modules). Components are manufactured under finite capacity constraints. We assume that final assembly is not a bottleneck, so we do not consider assembly capacity and cost. Demand is uncertain and, since backordering is not possible, unsatisfied demand is lost, with an adverse impact on the overall objective, which is to maximize expected profit. We assume continuous decision variables. This could be justified in a setting where demand volume is large enough. Considering integer variables would clearly increase computational requirements, but it would not have a significant impact on the issues we aim at addressing in our research.
There is a relevant amount of literature on several variants of ATO problems, and some early references from the 80s are [31,10]. Since an exhaustive review of the topic is beyond the scope of the present paper (we refer to [2] for a recent literature review), we focus on three possible ways to deal with ATO production planning within a rolling-horizon stochastic framework: two-stage models, two-stage models that quantify the value of the final inventory, and multistage models.
One of the first papers dealing with two-stage ATO problems, proposing a model similar to those we consider here, is [22]. There the authors formulate the problem as a stochastic integer programming model and solve fairly small problem instances by means of decomposition techniques, relying on progressive hedging and scenario aggregation. After almost 30 years of improvements in commercial solvers, it is likely that even some realistic size instances can be solved to near optimality by off-the-shelves software, at least in a two stage case. Nevertheless, as we have already pointed out, we are more concerned with alternative model formulations and realistic rolling horizon simulations, which is why we consider continuous linear programming (LP) models [9]. Another paper using a model similar to those presented here is [26], where the authors investigate the problem of finding optimal inventory levels.
The main weakness of two-stage models is that, by considering just two stages, they disregard what happens after the end of the second stage. This may lead myopic solutions or, in other words, to the so called end-of-horizon effect [19]. To decrease their myopia, two possibilities are either imposing target safety stocks or associating a value the the terminal inventory state. Both alternatives provide an incentive to store components at the end of the second stage, in order to use them in the following time periods to satisfy the demand.
The problem of sizing safety stocks has been considered in several papers and a recent literature review is provided in [17]. In this paper, the authors analyze a set of 95 articles published from 1977 to 2019 and show that 65% of them apply analytical techniques and do not contain references to realistic applications. Considering more recent papers, this characteristic seems to be still true. The few approaches dealing with realistic applications are mostly based on simulation optimization [12].
Instead, the problem of assigning a value to the terminal inventory state in production planning was addressed in [11], where the authors propose to include a valuation term for end-of-horizon inventory in the objective function to mitigate end-of-horizon effects. By explicitly modeling it, an appropriate level of ending inventory can be induced. Accordingly, we propose to integrate two-stage problems with an approximation of the terminal state value function, paving the way for approximate dynamic programming techniques [27]. Clearly, finding an exact expression of the terminal value function is out of the question. Nevertheless, experience with rollout algorithms [3] shows that even a rough approximation may be remarkably effective to improve performance. A useful strategy to this aim is to resort to separable approximations [30].
While two-stage problems have received a lot of attention, fewer papers consider multi-stage ATO problems. In particular, to the best of the authors' knowledge, there is no paper dealing with multi-stage stochastic ATO problems with multiple items. This gap in the literature is probably due to the increase in the problem size when multiple stages are considered [21,2]. In fact, several studies reduce the complexity of decision problem by concentrating on the replenishment decision and making specific assumptions regarding the allocation of the components [20]. We aim at filling this gap by formulating multistage models and by proving that, in the proposed setting, it is possible to solve multi-stage problems by off-the-shelves software in a reasonable amount of time.

Model formulations
In this section, we introduce different mathematical models to solve the ATO problem. Let I = {1, . . . , I} be the set of components, J = {1, . . . , J} the set of end items, and M = {1, . . . , M } the set of production resources (e.g., machines). Moreover, let us define the following parameters: • P j selling price of end item j ∈ J .
• K j lost sale penalty of end item j ∈ J .
• H i inventory holding cost of component i ∈ I.
• L m production availability (time) of machine m ∈ M.
• T im time required to produce one component i ∈ I on machine m ∈ M.
• G ij amount of component i ∈ I required for assembling end item j ∈ J , commonly known as gozinto factor when describing a bill of materials.
In order to formulate a stochastic model, we need a way to represent demand uncertainty (which is the only risk factor we consider). The typical choice in the SP literature is a scenario tree. Following the notation in [6], we define: • N , the set of nodes of the scenario tree and N + = N \ {0}. Node 0 is the root node, where only components production occurs; assembly occurs at nodes in N + , after observing demand. • p(n) the parent of node n ∈ N + .
• d [n] j the demand for end item j at node n ∈ N + .
We define as branching factor the number of children at each node of a given level of the tree. For example, by repeating three times a branching factor 2, the vector [2, 2, 2] identifies a binary tree over four time instants, including the current time corresponding to the root node, comprising eight scenarios. A scenario is a sequence of nodes from the root node to a leaf node, i.e., a realization of the underlying stochastic process. All of the nodes at the same level are associated with the same time instant. For each node n ∈ N , we define: i ≥ 0, the amount of component i ∈ I produced in node n ∈ N . • I Then, for each node n ∈ N + , we also define: j ≥ 0, the amount of end item j ∈ J assembled in node n ∈ N + . • l [n] j ≥ 0, the lost sales of end item j ∈ J in node n ∈ N + .
It is important to notice that variables y j are not defined for the root node n = 0. In other words, we decide on the end items after observing the current demand, carrying out the corresponding assembly. As we have pointed out before, we assume continuous decision variables. Hence, the resulting model is: whereĪ 0 i is the given initial inventory of component i. The objective function (1) is the expected net profit, expressed as the expected revenue, minus lost sales penalties and inventory holding costs, minus the cost of components production. Constraints (2) represent the capacity constraints on component production. Constraints (3) describe inventory balance of components. In particular, the inventory in node n is equal to the inventory in parent node p(n), plus produced amount, minus the quantity used for assembly. Constraints (3), together with non-negativity of inventory, also preclude assembling end items for which there are not enough components (either produced or in inventory). Initial inventory is imposed by constraints (4). Constraints (5) link end item assembly and lost sales with demand in each node n. Finally, constraints (6) and (7) specify the non-negativity of the decision variables. These constraints would be easily adjusted in the case of an integer model.
The size of the model (1) - (7), like any multistage SP model based on a scenario tree, increases exponentially as more time periods and greater branching factors are considered. In the following, we explore reformulation strategies to solve it effectively, without facing the computational complexity of the full-fledged multistage model. In particular, we consider the following strategies: • a plain two-stage formulation; • a two-stage formulation that associates a value to the terminal inventory state; • alternative multistage formulations characterized by different branching factors.
We start with the smallest possible tree, consisting of just two stages, as shown in Figure 1(a). In this setting, for notational convenience, we drop the superscript * [0] for variables associated with the root node and use the notation * s instead of * [n] for second-stage variables. Moreover, we rename N + as S, consistently with the two-stage programming literature. With these changes, model (1) -(7) becomes: In the following, model (8) -(14) will be called two-stage (TS). The main advantage of TS is that it can be solved by off-the-shelves software with a fairly large number of scenarios [9]. Nevertheless, the optimal solutions of this problem are myopic. Specifically, for each component produced but not assembled, the inventory cost H i is paid. Thus, there are no explicit incentives to leave components in inventory after assembly; leftover inventory will just be the result of hedging demand uncertainty in the immediate future (the second stage).
In order to overcome the myopic behaviour of TS, we consider a model that associates a value to the terminal inventory state. In other words, we add an extra term to the objective function in order to approximate the value of the components Figure 2: Example of a concave approximate value function.
that are produced, but are not immediately used for assembly in the second stage. By doing so, we obtain the following model: subject to the same constraints (9) -(14) of the TS model formulation, where V (I s ) = V (I s 1 , . . . , I s I ) is a function measuring the value of the terminal state. In principle, one way of obtaining this value function is by infinite-horizon dynamic programming [29]. Since finding an exact and tractable expression of value functions is generally not feasible [8], we devise a suitable approximation based on the following ideas: 1. The value function can be approximated by decomposition with respect to components, i.e., by setting V (I) ≈ i∈IV i (I i ). Clearly, by doing so we neglect interactions between components. Nevertheless, this decomposition approach is quite common and works often fairly well in practice [27]. 2. In order to preserve linearity, each component of the sum is approximated by a piecewise linear function. 3. In order to avoid an excessive computational burden, the information needed to build the piecewise linear function is collected by solving the two stage model TS for different values of the initial inventory state. In other words, we also assume V (I) ≈ TS(I), where TS(I) denotes the optimal value of the two-stage model obtained by settingĪ 0 = I in Eq (11).
Putting all of these ideas together, we approximate the value of terminal component inventories as a sum of piecewise linear concave functions like the one depicted in Figure 2. Concavity of the value function is a consequence of the fact that it is the optimal value of an LP in maximization form, as a function of the right-hand sides of constraints (11) (see, e.g., [1,Theorem 4.38]). Each individual term in the approximating sum is concave, as concavity is preserved by considering a concave function along its coordinates [5]. As shown in Figure 2, a slope can be negative for excessive inventory levels. It is important to observe that available inventory does not contribute to production costs in the TS model, but there is an inventory leftover cost, which penalizes huge inventories. These are further discouraged by the inventory holding cost in (15). The termV i (I s i ) can also be interpreted as a markdown value of the components, which is included in a term like is positive, leftover components of type i contribute to value, even if they are not sold within an assembled end item. This will be the case if future stages, neglected in model TS, require a quantity of components that is beyond production capacity.
To defineV i (·), we need to select a sequence of K breakpoints (denoted by u k i ) and a sequence of corresponding slopes 6 Algorithm 1 First order stock value approximation where we fix u 0 i = 0, ∀i ∈ I, and the sequence of slopes is non-increasing due to concavity. To choose the set of breakpoints and slopes, we consider a technique similar to Concave Adaptive Value Estimation (CAVE) introduced in [15,16]. In order to avoid incurring the curse of dimensionality, breakpoints are obtained by plain Monte Carlo sampling. At each iteration, we generate a pointÎ = [Î 1 , . . . ,Î I ] by sampling each component inventory level from a uniform distribution between 0 and an upper bound I max . Since we assume no inventory limit, I max must be set empirically. A good rule of thumb is to set I max to a multiple of the average demand; this is a value large enough to ensure that, up to a certain degree of confidence, demand in the next stages could be satisfied.
Given a sampled breakpoint valueÎ, we have to generate information about the slope on its right and on its left, in order to build the piecewise linear approximation of Figure 2. To this aim, we apply a straightforward scheme based on finite differences. The left and right slopes ofV i inÎ i , denoted by π + i and π − i , are obtained as follows: where e i is the i-th element of the orthonormal base of R I , and ε the perturbation of the approximation step. The idea is perturbing the inventory of component i and to check the impact on the optimal value of model TS, in terms of a forward and a backward finite difference. We assume that the number of scenarios in the two-stage model is sufficient to obtain a suitably reliable information. Nevertheless, some degree of smoothing an dynamic adaptation is needed when adding new sampled points. To this aim, we apply a strategy whose pseudo-code is given by Algorithm 1. Based on the proposed approximation, and using K to denote the set of breakpoints of the piecewise linear functionV model, we obtain the following model: ∀i ∈ I ∀s ∈ S ∀k ∈ K (20) m s ik ≥ 0 ∀i ∈ I ∀s ∈ S ∀k ∈ K (21) (9), (10), (11), (12), (13), (14), where the auxiliary variables m s ik are introduced to represent the piecewise linear approximation of (16), as usual in LP models [7]. In the following, we refer to this model formulation as First Order Stock Value Approximation (FOSVA). Both FOSVA and TS have the same type of scenario tree of Figure 1(a), but FOSVA is expected to be less myopic, at the cost of a preliminary learning step. This step requires an offline estimation of the coefficients by to solving K × I different TS problems.
Another possibility to make the TS model (8)-(14) less myopic is to consider a scenario tree as the one in Figure 1(b), which is still two-stage, but multiperiod. The idea is to add assembly periods, corresponding to nodes for which demand is known after the first branching. Actually, under an independence assumption, it can make sense to consider the expected value of demand (or its estimate) for each node. Clearly, the information structure of this model, which we will refer to as MP, is not realistic, but it aims at mitigating end-of-horizon effects by introducing a sort of certainty equivalent steady state, without increasing the number of scenarios.
We also consider the scenario tree structures shown in Figures 1(c) and 1(d). In particular, the scenario tree of Figure 1(c), labeled as MS3, relies on a more detailed description of uncertainty in the near future, as it comprises three stages. The structure of Figure 1(d) adds to this structure a sequence of nodes with a branching factor equal to 1. This is an extension of MS3 along the lines of MP, and it is labeled as MS3+.

Problem instance generation and design of experiments
A problem instance is characterized by a set of technological and economical data, by the probability distribution of demand, and by the structure of the scenario tree. The generation of problem instances aims at supporting the computational experiments, where the proposed model formulations are compared on a set of out-of-sample rolling horizon simulations. As it will be clear from the following, we purposely choose features that may not reflect a realistic setting, but allow us to create settings that are quite critical and stress the potentially myopic behavior of models that do not assign a value to the terminal state.

Technological and economical data
From the technological and economical viewpoint, a problem instance is characterized by the following groups of data (see Table 1): • capacity requirements and availabilities (T im , L m ); • economic parameters, i.e., costs, penalties, and sales prices (C i , H i , K j , and P j ); • product structures, i.e., the gozinto factors (G ij ).
First of all, we set the number of components I, machines M and end items J according to Table 1. Afterwards, we define a tightness factor γ to generate the available capacity for each machine L m , with respect to the average capacity requirement. In formulas, Whered j is the expected demand for end item j; given the complexity of the demand generation mechanism described below, this is estimated by sampling 5000 demand scenarios. In the experiments, we γ ∈ {0.8, 1.2}, corresponding to instances in which we may meet average demand directly from production or not.
8  [12,7,5,3,3] No. common components 2 No. components per family [11,17,12,6,9] Lost sale penalty K j 0.1 · P j Holding cost H i 0.2 · C i  Concerning the economic values, we randomly generate the cost of the components following a uniform distribution in the range [1,50]. Then, we divide the end items in three different classes of profit margin: low, medium, and high. In Table 1, we define the proportion of items in each class and the ranges from which these margins are uniformly sampled. By considering the bill of material of each end item and the cost of required components, we obtain the cost of each end item; the sales price is obtained by adding a profit markup, corresponding to each profit class. By a similar token, we obtain inventory holding costs for components and lost sales penalties for end items as a percentage of costs and sales prices, respectively.
In the literature, several standard structure for the gozinto matrix have been proposed [2]. Nevertheless, these structure are usually considered for their nice theoretical properties rather than their realism. Therefore, we define the gozinto matrix as follows. First, we define families by partitioning the set of items: each family has a given number of end items and, within a family, end items have a number of required components, either specific or common with other members. Moreover, we also introduce degenerate families, consisting of a single item. Such end items are, in a sense, apart from the others and called outcast items. In Table 1, we report the characteristics of the standard families (number of common/total components and number of items per family) and the number of outcast items. Once the total number of components of a (standard) family is known, we select the common components and the specific ones for each end item. Then, we sample the number of each component required by the end item, according to a discrete uniform distribution ranging from 1 to 9. On the contrary, outcast items do not share any specific pattern with others and their components are picked according to a binomial distribution, so that each component is included with probability 0.2. The gozinto matrix, generated by the aforementioned procedure and adopted for the computational experiments, is depicted in Figure 3. The diagonal blocks correspond to the standard families (five in this case), while the bottom rows correspond to the outcast items. Each standard family shares a number of common components (the first columns of each block).

Demand distributions
In order to stress the risk involved in production planning, we do not consider typical distributions like, say, the plain normal. On the contrary, we consider seasonality patterns and bimodal distributions. Furthermore, we also investigate the impact of correlation patterns inside families.
We consider monthly seasonality, modeled by factors that multiply both expected values and standard deviations of the base-case distributions. The multiplicative factors are given in Table 2 and apply to all end items in the same way. As we mentioned before, this need not be realistic, but in some sense it corresponds to a worst-case setting. In fact, representing seasonal demand accurately would require a long planning horizon, with a correspondingly huge scenario tree, and myopic models are going to suffer when all of the end items feature a demand peak at the same time. The capacity tightness factors of 0.8 and 1.2 that we mentioned before refer to the average (base-case) demand during the year.
Another relevant feature of end items demand that we assume is bimodality. Bimodal distributions may be observed in market conditions like the following ones: • New products. It is often difficult to predict whether a newly launched product will become a top seller or a shelf warmer.
• Big customers. Suppliers serving many small customers with independent demands face an aggregate demand that is approximately Gaussian. In the presence of a large customer that accounts for a large share of the sales, however, demand may be multimodal due to irregular bulk orders.
• New market entrants. The entrance (or exit) of a major competitor can have a significant impact on demand.
• Fashion trends. Apparel-type products are heavily influenced by fashion trends, and it may be unclear which style becomes popular. The corresponding demands are thus governed by multimodal distributions, where each mode reflects a particular popularity state.
We assume that distributions of all end items are bimodal for every time period. Again, this is rather questionable in practice. If bimodality is due to uncertainty about market penetration of a new item, after the first demand periods we should be able to figure out whether it is a success or a flop. We stress again that our assumptions aim at creating a very difficult environment on which to compare model formulations.
Another critical issue is correlation among end items demand. End items in different families can be completely unrelated in terms of consumer markets. Hence, assume that there is no correlation between families. However, there can be correlation within a family. If a consumer buys an end item with a certain configuration, she is not likely to buy a similar item with a different configuration, which implies a negative correlation. We explore both fully independent demands and correlation within families. We sample from a bimodal distribution that arises by mixing, with mixing parameter p, two normal distributions with expected values (µ 1 , µ 2 ) and standard deviations (σ 1 , σ 2 ) respectively. Thus, when correlation is disregarded, we have where BN refers to a bimodal arising by mixing two normals. We list the parameters of the employed bimodal distributions in Table 3.  The case of correlation within families must be handled with some care, as arbitrary negative correlations may lead to nonsensical correlation matrix (i.e., a matrix which is not positive semidefinite). Hence, we sample in two steps. Firstly, we independently sample the aggregated demand F r for each family r ∈ {1, . . . , R}, according to a bimodal normal distribution, but with parameters (n r µ 1 , n r µ 2 ), ( √ n r σ 1 , √ n r σ 2 ), and mixing parameter p, where n r is the number of end items in family r. Thus, F r ∼ BN(n r µ 1 , n r µ 2 , √ n r σ 1 , √ n r σ 2 , p), ∀r ∈ {1, . . . , R}.
Then, the overall demand per family r is split among the items belonging to the family (j ∈ J r ), according to weights w jr randomly sampled from a Dirichlet distribution with different parameters ζ r for each family r, i.e., W 1r , . . . , W nrr ∼ Dirch(ζ r ), ∀r ∈ {1, . . . , R}.
This choice ensures that j∈Jr w jr = 1, w jr > 0, ∀j ∈ J r , r ∈ {1, . . . , R}. Clearly, in case of families consisting of just one item (as for the outcast items), the demand for the unique end item is equal to the demand for the family. Therefore, the resulting distribution is: In Figure 4 we provide a graphical representation of the within-family correlation process. We report the weights sampled from the Dirichlet distribution to generate the within-family correlated demands in Table 4.

Scenario trees
We have outlined alternative structure of scenario trees in Figure 1. Now we must be more specific and describe how they are actually generated. In this case, too, our aim is to devise a challenging test. Hence, we assume that limited information is available. Specifically, we consider time periods of one month and assume that ten years of prior data about the demand are available. This means that we have 120 observations of demand for each end item but, given monthly seasonality, we only have 10 observations for each period. As we have noted before, this need not be realistic. Arguably, given the pace of product innovation, a shorter time span but a higher frequency would be employed. However, this allows us to experiment with a long seasonal pattern under limited distributional information, which makes a challenging testbed. In building the scenario trees, we adopt a purely data-driven approach, without trying to estimate the underlying demand distribution. We just bootstrap from the available data, which makes sense if we assume that there is no inter-temporal dependence. The limited in-sample data and the out-of-sample data used in the rolling horizon simulations are obtained from the same distribution. The seasonal pattern would require a minimum time horizon of 12 time periods, which would result in a huge scenario tree. We shorten the time horizon in order to compare models with some limited look-ahead against the two-stage model enhanced with the FOSVA approximation. The horizon covered by MP and MS3+ trees comprises seven periods (following the branching strategy introduced in Section 3), thus entailing (7 × 10 + 1) and (6 × 10 2 + 10 + 1) nodes, respectively. Instead, MS3 considers three periods and (10 2 + 1) nodes, while MP and FOSVA only rely on a second stage that involves ten different observations.

Experimental results
In this section we compare the different scenario structures on the basis of rolling horizon, out-of-sample simulations. We consider a simulation horizon of one year, with monthly production decisions. The impact of the scenario tree structure has been studied, among others, by [4], where it is shown that adding decision stages may actually have a detrimental effect, if this requires a less accurate representation of uncertainty in early stages. Our purpose is to study the effect of the terminal FOSVA term under a difficult demand testbed. Given the limited information we assume, we cannot afford building seasonality into the terminal value. Rather, we pool all of the information into a single function. Finally, we set the I max parameter of FOSVA equal to 4 j∈J G ij 1 S s d s j , i.e., four times the average demand for component i.

Performance metrics
The metric directly linked to the objective function of the models is cumulative profit, computed over the simulation horizon. Nevertheless, this is not informative with respect to the characteristics of the solution. Thus, we further investigate: the level of the inventory (linked to the holding costs), the cost of components production, and the lost sales penalties. For each of these quantities we provide a graphical representation per type of demand (with and without correlation within families) and level of capacity tightness in Figure 5 and 6. The sample mean and a 95% confidence interval for each policy (i.e., the production plans obtained by the alternative scenario tree structures and model formulations) is computed over ten different replications. Following the same experimental breakdown, we also report the mean of the disaggregated monthly profits across the simulations in Table 6.
Solution times of the different models were obtained using Gurobi 9.5 on a 10-Core, 3.2 GHz Apple M1 Pro processor and are reported in Table 5. Since we did not observe any significant difference in terms of computational effort with respect to the correlation on the demand, data Table 5 are classified only with respect to capacity tightness. It is clear that a tighter production capacity is associated with a larger computational effort, especially when the time horizon of the model is increased (MP, MS3+). However, as expected, the overall effort mainly depends on the number of nodes and variables (FOSVA and TS share the number of nodes, but they involve a different number of variables).

Family-correlated demand
In Figure 5, we report results obtained in the case of family-correlated demand. The 0.8 tightness case is on the left, whereas the 1.2 case is on the right. Overall, FOSVA and MS3 stand out with quite similar cumulative profits. We confirm this behaviour through the statistics about the profits in Table 6, where they provide the best performances as well.
In Figure 5(b) there is dip in profit around time period 4, which is actually due to inventory buildup in view of the seasonality peak: production costs are incurred, before revenue is collected later. This effect is more pronounced in the case of a larger production capacity (tightness factor 1.2 is on the right) and is also confirmed by the inventory plots of   On the contrary, MP and MS3+ feature larger inventory buildup with respect to MS3 and FOSVA, which are more conservative and stay between them and TS. However, this does not result in sharp differences in lost sales penalties. The point is that MP and MS3+ have a longer lookahead, but this relies on average demand values, which turns out to be a drawback in this setting, since the hierarchical demand generation mechanism implies a negative correlation between families. As a result, a longer lookahead based on a wrong anticipation of demand implies some overconfidence in production, and results in the production of the wrong components.
Figure 5(g) shows an interesting pattern in terms of production cost. Given the limited capacity, all of the model formulations seem to saturate production volume, with the exception of TS, which shows marked volatility trying to chase the seasonal demand pattern.

Independent demand
The independence in end items demand makes the problem definitely easier. Even though the demand volume is the same as in the case of correlation within families in terms of expected value, we obtain a definite improvement in cumulative profit: Figures 6(a) and (b) suggest a 25% and 50% improvement, depending on capacity tightness. These cumulative profit plots display similar behavior for the models, again with the exception of TS. The plots of cumulative profit display again a dip around time period 4 and at the end of the simulation horizon, but these are related to inventory buildup that results in revenues only later on. This does not happen with TS, which behaves like the previous case: the myopic model does not build enough inventory, suffers from lost sales, and features pronounced volatility in production costs, as it tries to chase the seasonal demand. In fact, the cumulative profit of TS falls behind the competitors around time periods 4 and 5, where the lack of inventory buildup results in a peak in lost sales. The additional terminal value of FOSVA helps in reducing the issue.
A clearer picture can be obtained by looking at Table 6, which shows that in this case the gap between TS and the other policies is reduced with respect to the case of correlations within the families. This is indeed a result of the reduced complexity of the independent demand case.
The seasonal pattern is more noticeable in the inventory levels of every policy (Figure 6(c) and 6(d)) compared to the family-correlated case ( Figure 5(c) and 5(d)). When the capacity is tighter (Figure 6(g)), the production is not saturated as it was on the family-correlated demand ( Figure 5(g)). This occurs because the easier case of independent demand allows to manage inventory better in view of the seasonal peaks.

Conclusions and directions for further research
The main issue we have addressed in this paper is the impact of augmenting the objective function of a stochastic program with a value measuring the value of the terminal state in a challenging ATO production planning problem. The aims of the additional term is to reduce the potential end-of-horizon effects and to obtain a solution of similar quality as that obtained with a longer planning horizon, while reducing the computational effort.
The computational experience reported in [4] within a financial context suggests that a longer planning horizon need not improve results. On the contrary, when this is obtained by reducing the branching structure for early stages or the number of scenarios, a detrimental effect is observed. Our experience is similar. We have tried adding some lookahead by appending a linear graph to leaves of a tree, which is not feasible in a financial context, as it would entail arbitrage opportunities. Experience shows that, in our setting featuring seasonality and bimodality, this deteriorates performance. A two-stage formulation performs badly, especially in presence of a correlation structure introduced by the hierarchical structure of product families, which is a natural occurrence in an ATO setting. The addition of a terminal value improved performance considerably, even though it is a rough-cut approximation based on a decomposition with respect to components, which neglects interaction among components through the structure of bills-of-materials, and ignores seasonality due to the assumed limited availability of data.