Price of robustness optimization through demand forecasting with an application to waste management

Robust optimization can be effectively used to protect production plans against uncertainties. This is particularly important in sectors where variability is inherent the process to be planned. The drawback of robust optimization is the chance of producing over-conservative solutions with respect to the real occurrences of the stochastic parameters. Information can be added in order to better control the extra cost resulting from considering the parameter variability. This work investigates how demand forecasting can be used in conjunction with robust optimization in order to achieve an optimal planning while considering demand uncertainties. In the proposed procedure, forecast is used to update uncertain parameters of the robust model. Moreover, the robustness budget is optimized at each planned stage in a rolling planning horizon. In this way, the parameters of the robust model can be dynamically updated tacking information from the data. The study is applied to a reverse logistics case, where the planning of sorting for material recycling is affected by uncertainties in the demand, consisting of waste material to be sorted and recycled. Results are compared with a standard robust optimization approach, using real case instances, showing potentialities of the proposed method.


Introduction
Robust optimization approaches are commonly applied in solving mathematical programming problems where a certain set of parameters are subject to uncertainty. Considering either production or procurement planning problems, these are dealing with the stochastic nature of demand, and a deterministic approach would definitely fail to provide concrete decision support when modeling those kind of scenarios. Therefore, demand variability across the planning time horizon should be properly addressed. This can be done by introducing a protection function in each constraint Communicated  according to the probabilistic robust approach presented in Bertsimas and Sim (2004). This approach ensures deterministic and probabilistic guarantees on constraints satisfaction, and it does so in a linear framework.
In order to trade-off between optimality and level of robustness, a risk parameter has to be set to formulate the robust counterpart of the nominal model and its level of conservatism. At the same time, the level of conservatism of robust solutions can be such to constitute a significant cost in terms of optimality reduction, also known as price of robustness. This is often the case of over-conservative robust solutions which consider demand values perturbations with low probability of occurrence. In a production planning, setting these solutions leads to extra costs resulting from additional production and storage.
Either risk-averse or risk-taking decision makers will struggle to deal with the challenging trade-off of risk management by setting the robustness control parameter only, besides its odd interpretation. This topic is studied by the literature of robust approaches aiming at producing models better supporting the reality of decision making in uncertain scenarios. Adjustable robust models, for example, are a branch of robust optimization introduced in Ben- Tal et al. (2004) where some of the decision variables can be adjusted after some portion of the uncertain data reveals itself. This approach offers increased flexibility and produces less conservative solutions with respect to static robust optimization. An interesting application of this approach in scheduling problems is presented in Bruni et al. (2017). In the aforementioned paper, the authors present an adjustable robust formulation where sequencing decisions are taken in a first stage, and scheduling decisions are made in a second stage. We refer the reader to the survey in Yanıkoglu et al. (2019) for a substantial knowledge of adjustable robust optimization (ARO) literature. Other approaches are soft robust optimization (Ben-Tal et al. 2010), light robustness (Fischetti and Monaci 2009), scenario-based robust optimization (Goerigk and Schöbel 2011) and the one proposed in (Roos and den Hertog 2020).
Demand is the most critical information input of production planning and the main source of uncertainty as well. Addressing demand forecasting with proper time series analysis and regression models plays an important role in the overall decision processes. Features such as fitness of demand regressions models should be taken into account along with their parameters setting and forecast accuracy metrics.
Supply chain management is one of the most used scenario to prove the potentialities of robust optimization. Some fundamental supply chain and inventory management models are revised using robust optimization theory in Bertsimas and Thiele (2004); Bertsimas and Thiele (2006). The authors develop robust models for the optimization of the inventory in different settings and policies, such as the (s, S) case, allowing to control the level of conservatism of the solution, without assuming a specific demand distribution. The models are applied to a multi-period planning problem for single or networked warehouses, and for the uncapacitated and capacitated cases. The networked capacitated case is similar to the model studied in our work as an application case. Multiperiod inventory management is addressed in detail in See and Sim (2010) where the concept of budget of uncertainty proposed in Bertsimas and Thiele (2004); Bertsimas and Thiele (2006) is extended for controlling the demand. This is done considering demand descriptive information such as standard deviation, seasonality and autoregressive aspects. The inventory policy introduced by the authors is solved by means of Second-Order Cone Programming (SOCP).
In real-case applications, robust optimization could induce over conservatism. The problem is addressed in Daneshvari and Shafaei (2021) where correlation between uncertain parameters is taken into account in the definition of uncertainty sets in order to mitigate over conservatism. Another methodology used to shape the uncertainty sets based on data analysis has been proposed in Qiu et al. (2020). The authors propose a data-driven robust approach of modeling the multi-product inventory problem with demand uncer-tainties. Basically, they construct an uncertainty set using historical demand data which are estimated not via demand descriptors but using a support vector clustering (SVC) model.
In this paper, we aim to provide a mechanism for setting the level of conservatism of robust solutions according to accurate estimates of robustness costs. Thus, a framework integrating both a forecasting model and two extensions of the nominal multi-period planning model is proposed.
The central idea is to iteratively solve planning and forecasting problems with the goal of finding the best configuration of robust parameters that minimizes the costs that result from overestimating or underestimating risk. These costs are the so-called price of robustness (Pr) as described in Bertsimas and Sim (2004) and the potential extra cost resulting from overtime production (Po) whenever demand is underestimated. In the considered setting, the way of dealing with demand underestimation errors is indeed the activation of overtime production, while demand overestimation contributes to Pr costs.
In order to study the performance of the proposed method, a robust optimization problem arising from waste management and inverse production planning with demand uncertain coefficients is solved and variants of the model with different protection strategies are compared with real-case instances.
The remainder of the paper is organized as follows: Sect. 2 is dedicated to the definition of the planning problem formulations; the proposed framework aimed at the optimization of the risk considering parameters is presented in Sect. 3; Sect. 4 presents the experimental results, while Sect. 5 gives some conclusions, research perspectives and reports acknowledgments.

Problems definition and modeling
This section presents the formulations of the problems used in the proposed framework which is detailed in Sect. 3.
Let τ be a planning period made of T time slots. Each time slot is indexed by t ∈ T = {1, . . . , T } and corresponds to a working shift where production can be activated. Thus, the planning consists of defining the production lot sizes and the schedules in order to meet the demand quantity d t of each time slot t. Therefore, each planning period τ is given with a foreseen demand vector d T ∈ R T + . Considering a deterministic planning problem D, we refer to R as its robust counterpart, while a third model E replicates the same deterministic formulation of D and includes some additional decision variables regarding overtime production. Overtime production is the assumed strategy to counteract, with an extra cost, the uncovered demand.
The formulations of problems D, R and E are presented in Sects. 2.1, 2.2, and 2.3, respectively.

Deterministic model D
Dealing with a production planning setting, D is a mixedinteger linear programming model to schedule and lot-size production operations. To better introduce a general formulation of D, model notation for parameters and indexes is set out in the following.
T : time horizon length; t ∈ T = {1, . . . , T } : index of working shifts across time f t : setup cost of working shift at time t; d t : production demand at time t; k t : production capacity at time t; c t : unitary production cost at time t; h t : unitary inventory holding cost at time t; I 0 : initial inventory level.
The model considers the following variables. y t ∈ {0, 1} : equal to 1 if production is activated at time t, 0 otherwise u t ∈ R + : production quantity at time t I t ∈ R + : inventory level at time t The model minimizes the sum of production, setup and holding costs and is detailed as follows: u t ≤ k t y t ∀t ∈ T (2) The objective function (1) defines the minimization of the sum of the three cost terms, which are production, setup and inventory costs, respectively. Constraints (2) bound the production quantity for each time period, and constraints (3) define the inventory.

Robust model R
The mixed-integer linear programming model R defines the robust counterpart of D according to the robust approach presented in Bertsimas and Sim (2004). We first consider the parameters d t , t ∈ T = {1, . . . , T }, which are subject to uncertainty to take values according to a symmetric distribution with mean equal to the nominal value d t in the interval [d t − σ t , d t + σ t ]. Indeed, σ t is the maximum deviation allowed for d t . In order to meet the standard formulation of the nominal problem presented in Bertsimas and Sim (2004), where parameters subject to uncertainties belong to inequality constraints only, the equality constraints (3) of D regarding demand quantities d t are reformulated to turn them into inequality constraints. This is performed considering, for each period t, the sum of all demanded quantities up to that period, as in constraints (7) and (8) of the formulation of R presented in this section. According to the robust approach presented in Bertsimas and Sim (2004), a parameter i is introduced for each constraint i holding one or more uncertainty coefficients. i is not necessarily integer and takes values in the interval [0, |J i |] where J i is the set of the coefficients of constraint i being subject to uncertainty. The nominal problem D presents only one set of T constraints considering the coefficients d t , and these are the ones reformulated as inequality constraints. Therefore, we get ∈ R T + , and because of this reformulation |J t | = t ∀t ∈ {1, . . . , T }. For each period t, t represents the number of coefficients that we consider as allowed to vary within their interval. In practice, we consider nature behaving like only a subset of the coefficients will change with respect to their nominal value. Indeed, as affirmed in Bertsimas and Sim (2004), it is unlikely that all |J t | will change; so the idea of conservative robustness is to be protected against all cases that up to t of these coefficients are allowed to change, and one coefficient d t changes by ( t − t )d t . Note that when t = 0 ∀t ∈ {1, . . . , T }, we get the nominal deterministic scenario, while setting t =|J t | = t ∀t ∈ {1, . . . , T } represents solving the problem of the worst case scenario. It is clear that by varying the level of robustness can be flexibly adjusted against the level of conservatism of the solution. Considering the structure of the constraints including d t is important: Because of the telescopic expansion of each set J t as t goes from 1 to T (i.e., |J t+1 | =|J t | + 1), we consistently constraint t to be bigger than or equal to t−1 . In the following, we present all the additional variables and parameters that are required to introduce the robustness protection functions presented in Bertsimas and Sim (2004) and formulate the robust counterpart of D: p t,k ∈ R + : variables resulting from duality within Bertsimas and Sim's robustness theory; they contribute to the protection function of constraint t with respect to the specific coefficient d k . s t ∈ R + : variables resulting from duality and Bertsimas and Sim's robustness theory; multiplied by σ t , they set the lower bound of the protection function contribution in each constraint t. t : parameter to adjust the level of robustness of each period t.
The robust counterpart R of D is introduced as follows: (2), (4), (5), (6) Constraint (7) defines inventory considering the cumulative demand d t up to period t, the overall produced quantity u t up to period t, and the uncertainties protection function made of the joint contribution of z t t and the sum of p t,k for k ∈ {1, . . . , t}. Constraint (7) sets the lower bound of the inventory for each period. Equality constraint (8) allows the inventory to be considered within the cost function. Constraints (9) and (10) resulted from duality in Bertsimas and Sim (2004) robustness theory, where (9) sets the lower bound of the protection function contribution in constraints (7) and (8).
Solving model R allows us to evaluate the price of robustness (Pr) resulting from a specific protection . This is done by computing the difference between the optimal objective function value of R and the optimal objective function value of D.

Overtime production model E
In our approach, overtime production is to be activated whenever standard scheduled production is not enough to meet the entire demand across the time horizon. For this reason, the third model E replicates the same deterministic formulation of D and includes some additional decision variables regard-ing overtime production. This additional production has of course a higher unitary cost with respect to standard scheduled production. Solving model E is intended to evaluate the need and the cost of overtime production (Po).
Considering the main decision variables of D, such as y t ∈ {0, 1} and u t ∈ R + , the formulation of E incorporates also y t ∈ {0, 1} and u t ∈ R + which correspond to the overtime production equivalent of the previous ones. These additional variables allow E to capture the extra costs resulting whenever the robust solution of R is not feasible when the true demand vectord is observed. Indeed, while both D and R are solved considering the predicted demand vector d, instances of E are solved w.r.t. the true demand vector d and its set of standard production variables (i.e., y t and u t ) are fixed to the values of the robust solution of R for the same corresponding scheduling time window. The corresponding unitary costs of these overtime production variables are f t > f t and c t > c t , respectively. Therefore, the formulation of E is the following: Solving E allows the evaluation of the cost of overtime production (Po) as the difference between the optimal objective function value of E and the optimal objective function value of D. Thus, the sum P of the price of robustness (Pr) and overtime production (Po) can be computed according to a specific protection (i.e., P( ) = Pr( ) + Po( )). We introduce a procedure, which we call Eval_P, that generates a single evaluation of P( ). This procedure is part of the routine, detailed later in Sect. 3, aimed at the optimization of the robustness control vector . The pseudo-code of Eval_P is set out below.
Eval_P ,d,d, σ : . Solve E d where y and u are fixed to y * and u * , respectively, to obtain Obj *

Robustness control optimization
The main notation used in the proposed framework is defined in Table 1. The forecasting model F employed in the framework to evaluate the expected demand is presented in Taylor (2017); its implementation is available as open source software and is called Prophet, and it uses a decomposable time series model Harvey and Peters (1990) with three main model components: trend, seasonality and holidays.
The robustness control vector is learned by a routine properly built in order to achieve this parameter tuning objective. The optimization of is to find the value that minimizes the sum P( ) of the price of robustness Pr( ) and the potential extra costs resulting from overtime production Po( ). We present a routine that seeks for the best configuration of for a considered planning period τ with its corresponding demand forecastd(τ ) and its true observationsd(τ ). Its pseudo-code is set out below. We refer to this routine as Opti-mize_ 0 , step ,d(τ ),d(τ ), σ (τ ) : 1. Set 0 ∈ R T + : the starting configuration of 2. Set step ∈ R T + : a vector of the updating step sizes of each component t of 3. Set max k = f loor (1/ step) : maximum number of iterations k for a complete search over the space of ∈ R T Whatever setting of 0 and step is made, this must be done according to the peculiar structure of constraints (7), (8) of the formulation of R. Indeed, as stated in Sect. 2, before introducing model R, the telescopic expansion of each set J t (i.e., |J t+1 | =|J t | + 1) constrains each component t of the protection vector to be bigger than or equal to the component t−1 and smaller than t. The setting of both 0 and step must be consistent with these constraints. Given a scalar γ step ∈ {0 . . . 1}, we propose a rule for setting each component step (t) of step such that the previous constraints are satisfied: This Optimize_ routine can be embedded in the proposed framework for solving the multi-period planning problem under uncertainties. The purpose of this framework is to provide automated decision support to implement an appropriate robustness control vector that minimizes the aforementioned Deterministic version of the planning model Robust counterpart of the planning model E (d(τ )) D (d(τ )) with additional variables modelling overtime production Referring to the most recent demand realizations, and to its corresponding time series D(ρ, τ ), we consider τ planning period as the latest T time slots part of D(ρ, τ ). This demand time series D(ρ, τ ) can be used along with Optimize_ to mine what would have been the best robustness control vector * for the planning period τ . This is indeed described in the following: 1. Set ρ and τ , then consider D(ρ, τ ) and getd (τ ) 2. Solve F (D(ρ, τ )) to forecastd(τ ) andσ (τ ) 3. Set 0 and step The rationale of the proposed framework is the use of D(ρ, τ ) and its corresponding best configuration * to solve the planning problem of the future planning period τ + 1. Therefore, for each observed demand time series D(ρ, τ ), the corresponding * impact in term of P( * ) is evaluated when * is used to protect against the uncertainties of the very immediate future period τ + 1. Considering a sequence of framework implementations, at each of these it is reasonable to set 0 starting configuration to * τ −1 . The experimental results of both Optimize_ and its implementation framework are presented in Sect. 4. Performances in terms of P costs are compared with P resulting from either the deterministic and worst-case scenario approaches that we refer to as P( nominal ) and P( worstcase ), respectively.

Application to waste management
This section holds the experimental results obtained by applying the proposed methodology to a waste management (WM) setting.
With respect to the problem addressed, while it is possible to find several applications of robust optimization to supply chain and production problems, fewer research works address circular economy and reverse logistics related cases. Stochastic-based approaches in waste flow optimization have been addressed in Sun and Huang (2015) and Zhu and Huang (2011). The work of Kim et al. (2018) develops a robust model where box uncertainty and robustness budget are used to control uncertainty in a closed loop supply chain in the textile industry. The paper considers recycling operations in a networked environment and considers the over conservatism problem but does not consider demand characteristics. In Hatefi and Jolai (2014), p-robust constraints are used to control disruption events. In Entezaminia et al. (2017), a multi-objective multi-period multi-product multisite aggregate production planning model is solved with reverse flow considerations. In this case, the robust approach is pretty standard. In Pinto and Stecca (2020), the same problem addressed in this work is treated with classical robust theory.
In WM, the planning consists of scheduling and lot sizing each phase of a recycling material sorting process. The planning model presented in Pinto et al. (2021) is considered as the deterministic model D for testing the proposed framework. The formulation of D is reported in Sect. A of the appendix. This model supports several strategic decisions that are critical in the considered WM application. It can be described as a variant of a lot sizing model with nonlinear costs (approximated by mean of piecewise linear functions) with the additional features of scheduling the operations and allocating the appropriate workforce dimension. The aim of the model is finding the best allocation of operators to each working shift in order to process the recycling material quantity d t arriving at each time slot t. In this scenario, a time slot is indeed a working shift where teams of operators may be allocated. The robust counterpart of D is newly introduced in Pinto and Stecca (2020) and is considered as R for what concerns the experimental results. The formulation of R is reported in Sect. B of the appendix. Model E is derived from model D in Pinto et al. (2021) as illustrated in Sect. 2.3, and its formulation is reported in Sect. C of the appendix.
Optimize_ routine is supposed to densely explore the search space of for an exhaustive search aimed at the optimization of P. In order to validate its performances, some preliminary experiments are performed over a set of subsequent planning periods: For each planning period τ , performances in terms of P costs resulting from applying * found by Optimize_ are compared with P resulting from the deterministic (i.e., nominal) and worst-case scenario approaches. The results are illustrated in Table 2 and discussed below. The accuracy of the forecast model is reported in terms of symmetric mean absolute percentage error (i.e., SMAPE), an accuracy measure based on percentage errors which is usually defined as follows: where A t is the actual value and F t is the forecast value. The absolute difference between these values is divided by half the sum of absolute values of the actual value A t and the forecast value F t . Then, the value of this calculation is summed for every fitted point t and divided by the number of fitted points n.
The table's columns present the reduction of robustness and overtime production costs P achieved by P( * ) with respect to the deterministic P( nominal ) and worstcase P( worstcase ) approaches. The presented values of costs reduction result from the following expressions: The following columns of Table 2 present the percentage of protection guaranteed by each * considered, obtained as the ratio between the sum of the components of * and the sum of the components of the vector corresponding to the worst case, the SMAPE accuracy metrics of the forecasting model, while the last column presents the sum of all the forecasting residuals of the planning period to consider forecast bias, that is overestimating (a positive value of bias) or underestimating (a negative value of bias) the future demand. Indeed, the forecast model is implicitly introducing robust-ness to the planning solution whenever is overestimating the demand realization. At the same time, the model F is incautiously fostering the importance of robust planning solutions whenever is underestimating the demand realization. Indeed, considering systematical biased behavior of forecast models is crucial for a proper planning. Results of Table 2 show how * moderate its robustness contribution with respect to forecast bias over the planning period: its percentage protection remains minimal (or close to zero) for positively biased prediction in order to produce solutions close or equal to a deterministic (nominal) solution, while it provides a stronger protection to uncertainties when dealing with underestimating demand forecasts. However, for both types of biased prediction scenarios, Optimize_ routine produces configurations of minimizing risk costs with respect to either deterministic and worst-case approaches. Figure 1 shows an example of Optimize_ instance solving.

Analysis of results
As described at the end of Sect. 3, considering an observed demand time series D(ρ, τ ), the proposed framework essence is using the very close past best configuration * to protect against the uncertainties of the very immediate future, that is, solving the robust planning problem R d , σ, * of the subsequent future planning problem τ + 1. Therefore, a new set of experiments is performed to test the performance relative to the observed costs P of robustness and overtime production of the very next future. All instances are created by a real-world case study from a waste sorting plant located next to Rome, Italy. Results are shown in Table 3. As in Table 2, performances in terms of P( * ) costs are compared with P resulting from either the deterministic P( nominal ) and worst-case P( worstcase ) scenario approaches. Demand arising from the considered real-case study is particularly unstable, and it seems reasonable to consider this adverse feature of the time series D particularly challenging for the proposed framework. Indeed, no autocorrelations are present either for small or larger seasonal lagged values of D. Therefore, the larger the planning period τ + 1 towards the future, the smaller the chances of a good performance of * learned over the immediate past planning period τ . Experiments presented in this section consider 50 instances of subsequent planning periods of two weeks (i.e., 12 time slots considering 6 weekly working days of 2 working shifts each). Table 3 highlights two main important results: The framework provides stronger protections to uncertainties (as featured in * column) for most of the occasions of underestimation of demand realization, and it does so while guaranteeing a lower cost P of robustness and eventual overtime production with respect to the deterministic and worst-case approaches. Indeed, out of the 30 instances where F underestimates the Fig. 1 Example of an Optimize_ instance solving  demand, 21 of these reveal a valuable P cost reduction, that is about 70% of the occasions. At the same time, costs reduction are also observed for 7 of the remaining 20 occasions (35%) where F overestimates the demand. Moreover, 42 of all 50 instances (84%) prove how the framework is providing protection vectors * with smaller associated P( * ) with respect to P( worstcase ).

Conclusions
In this paper, we addressed the problem of controlling the extra costs resulting from considering demand uncertainties in planning problems. Robust optimization theory is applied to both a general production planning model in the general discussion and to a waste management use case regarding recycling operations planning. An optimization procedure determining the best possible robustness budget to use in successive planning stages is proposed. In this way, it is possible to dynamically update both the estimated demand and the robustness budget, helping to balance the extra costs of robustness and extra costs due to overtime production. The framework is tested in a real-case environment where demand does not follow a specific probability distribution. The results show that, most of the time, the classical robust approach is over-conservative. Moreover, when the forecasts overestimate the observed demand, the robustness budget can be effectively optimized. Forecasting model features and their overall effect over the protection cost are not addressed by this study, and they will be done in a future work. Moreover, the proposed approach can be applied to different planning problems and tested against different scenarios in future works.

A : Model D
Dealing with a waste management setting, D is a mixedinteger linear programming model scheduling the sorting operations of each phase of a waste sorting process. The model can be described as a variant of a lot sizing model with nonlinear costs (approximated by mean of piecewise linear functions) with the additional features of scheduling the operations and allocating the appropriate workforce dimension.
To better introduce the formulation of D, model notation for parameters and indexes is set out in the following.  The model considers the following variables.
x j,t ∈ Z + : operators employed in the sorting stage j at time t u j,t ∈ R + : processed quantity at stage j at time t y j,t ∈ {0, 1} : equal to 1 if stage j is activated at time t , 0 otherwise I j,t = I j,t + I j,t ≥ 0: inventory level of material in buffer j at time t; for each stage j the corresponding I j,t and I j,t represent the inventory level before and after reaching the critical threshold, respectively. w j,t ∈ {0, 1} : equal to 1 if I j,t > 0, 0 otherwise. Indeed, these binary variables are used to model the piecewise linear functions of the buffer inventory costs.
The model minimizes the sum of sorting and holding costs and is detailed as follows: min Z = j∈J t∈T C t x j,t + j∈J t∈T f j y j,t + j∈J t∈T ∂h 1 j I j,t + ∂h 2 j I j,t u j,t ≤ SK j,t x j,t ∀ j ∈ J , ∈ T (21) I 1,t = I 1,t−1 + d t − u 1,t ∀t ∈ T I j,t = I j,t−1 − u j,t + α j u j−1,t ∀t ∈ T , j ∈ J \ 1 (23) I j,t = I j,t + I j,t ∀ j ∈ J , t ∈ T LC j w j,t ≤ I j,t ≤ LC j ∀ j ∈ J , t ∈ T 0 ≤ I j,t ≤ (S j − LC j ) w j,t ∀ j ∈ J , t ∈ T (26) I j,T ≤ n j LC j ∀ j ∈ J x j,t ∈ Z + ∀ j ∈ J , t ∈ T (28) u j,t ∈ R + ∀ j ∈ J , t ∈ T (29) y j,t ∈ {0, 1} ∀j ∈ J , t ∈ T The objective function (18) defines the minimization of the sum of the three cost terms, which are sorting, setup and inventory costs, respectively. Constraints (19) and (20) bound the number of workers that can be assigned to each sorting station and to each time shift. Constraints (21) limit the quantity sorted u j,t to the sorting capacity dependent on the number of workers x j,t . The remaining constraint sets define and limit the inventories: Constraint set (22) defines the inventory for the first buffer, considering the inbound material d t and the sorted material u 1t , while (23) defines the inventory for the other buffers corresponding to j > 1. Indeed, constraints (23) describe the waste flow across the sorting stages that follow one another: Each subsequent interoperational buffer j receives by the previous sorting stage j − 1 a quantity of waste equal to a α j percentage of the waste processed in stage j −1. Constraint sets (24), (25), and (26) define the piecewise linear functions for inventories; in these constraints, level S j and maximum capacity LC j are connected with the inventory levels through the variable w j,t . The last constraint set (27) imposes the maximum unsorted material allowed to be left at the end of the planning period for each buffer.

C : Model E
Model E replicates the same deterministic formulation of D and introduces some additional decision variables to model overtime production. Consider the following set of main decision variables of D: y j,t ∈ {0, 1} : equal to 1 if production stage j is activated at time t, 0 otherwise x j,t ∈ Z + : operators employed in the production stage j at time t u j,t ∈ R + : processed quantity at production stage j at time t The formulation of E incorporates also the variables y j,t , x j,t and u j,t which correspond to the overtime production equivalent of the previous ones. The corresponding unitary cost of these overtime production capacity variables (i.e., y j,t ∈ {0, 1} and x j,t ∈ Z + ) is f j and C t , respectively, while M is the maximum number of operators available for overtime production. Therefore, the formulation of E is the I 1,t = I 1,t−1 + d t − u 1,t − u 1,t ∀t ∈ T (41) I j,t = I j,t−1 − u j,t + α j u j−1,t +α j u j−1,t ∀t ∈ T , j ∈ J \ 1 ( 4 2 ) x j,t ∈ Z + ∀ j ∈ J , t ∈ T (43) u j,t ∈ R + ∀ j ∈ J , t ∈ T (44)