Price of Robustness Optimization through Demand Forecasting and Robust Planning Models Integration in 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 that must be optimally planned. The drawback is that, in real situations, some 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 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 the waste material to be sorted and recycled. Results are compared with the 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 Diego Maria Pinto · Claudio Gentile · Giuseppe Stecca Istituto di Analisi dei Sistemi ed Informatica "A. Ruberti", Consiglio Nazionale delle Ricerche, Via dei Taurini 19, 00185 Rome, Italy Email: diegomaria.pinto,giuseppe.stecca,gentile@iasi.cnr.it 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 decisions supporting 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 according to the probabilistic robust approach presented in [1]. This approach ensures deterministic and probabilistic guarantees on constraints satisfaction and it does so in a linear framework. Furthermore, in order to trade-off between optimality and level of robustness, a risk considering 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 particularly evident for over-conservative robust solutions, these indeed will consider demand values perturbations with low probability of occurrence, leading to extra costs resulting from additional production and storage. Either risk adverse 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 subject is rarely considered by researchers and hence the proposed robust approaches normally lead to models not really supporting the reality of decision making. Demand forecasting 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 play an important role in the overall decision processes. One of the reason behind this work is considering reasonable to set the level of conservatism of robust solutions according to the expected evaluation of its cost. Indeed, this research addresses the aforementioned challenge of robustness parameters tuning with a framework considering both a forecasting model and two additional variants of the nominal planning model. The planning models considered features is made of a selection of the most informative parameters of the robust formulation, the resulting price of robustness and any potential extra costs observed. These extra costs are estimated by the solution of a formulation devoted to this task. Features such as goodness of fit of demand regressions models should be taken into account along with their parameters setting and forecasting accuracy metrics. 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 reminder of the paper is organized as follows: the literature review is given in Section 1.1; Section 2 is dedicated to the definition of the problems and the proposed routines aimed at the optimization of Γ ; Section 3 presents the experimental results while Section 4 gives some conclusions, research perspectives and reports acknowledgments.

Literature review
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 by Bertsimas and Thiele [2,3]. 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 conservativeness of the solution, without assuming a specific demand distribution. The models are applied to a multiperiod 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 application case. In our case we study the reverse logistics case and we consider specific forecasted parameters. Multiperiod inventory management is addressed in detail by See and Sim in [14]. The authors extend the concept of budget of uncertainty proposed by Bertsimas and Sim for controlling the demand, and they consider descriptive information on the demand, 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 [5] where correlation between uncertain parameters is taken into account in the definition of uncertainty sets in order to mitigate over conservatism. In our paper we have a similar aim but we focus on using the forecasted data in order to mitigate over conservatism while preserving the advantages of robust optimization approach. Another methodology used to shape the uncertainty sets based on data analysis has been proposed in [13]. The authors propose a data-driven robust approach of modeling the multi-product inventory problem with demand uncertainties. In practice, they construct an uncertainty set using historical demand data which are estimated not via demand descriptors but using machine learning algorithm Support Vector Clustering (SVC). The experimental data are randomly generated from Bivariate Gamma and Mixed Gaussian distributions. The produced uncertainty sets are compact and convex, preserving the linearity of the model, as proven by [15]. With the respect to [13], we use forecasted demand in conjuction with a procedure to optimize the uncertainty set, and then apply the proposal to a specific waste recycling optimization model.
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 [16] and [18]. The work of [10], develop a robust model where box unertainty and roubustness budget is used to control uncertainty in a closed loop supply chain in the textile industry. The paper considers recycling operations in a networked environment, considers the over conservatism problem but does not consider demand characteristics. In [9], p-robust constraints are used to control disruption events. In [6] a multi-objective multi-period multiproduct multi-site aggregate production planning model is solved with reverse flow considerations. In this case the robust approach is pretty standard. In [12] the same problem addressed in this work is treated with classical robust theory.

Problems definition and modeling
In the following, we introduce a framework of sequential forecasting and planning problems solving in order to mine the best robust parameters tuning that minimizes both the price of robustness and any future potential extra costs resulting from overestimating or underestimating risk.
Let N be consecutive planning periods and in each of them T timeslots, having a total of NT timeslots. A timeslot is a shift where workers must be allocated. Thus, the planning consists of scheduling and lot sizing each phase of a recycling material sorting process while finding the best allocation of operators to each working shift in order to process the recycling material quantity d t arriving at each timeslot t. At the beginning of each planning period τ ∈ {1, . . . N}, d T is forecasted for the T successive timeslots, as each timeslot of the planning period τ is indexed by t ∈ {1, . . . , T }. The main notation that will be used in the models, such as parameters and indexes, is the following: τ ∈ {1 . . . N} : planning time period T : timeslots for each planning period τ Γ (τ) : robustness control parameters vector for planning period τ ρ : number of planning periods used as forecasting model training data D(ρ, τ) : demand time-series considering the planning periods from τ − ρ to τ corresponding to timeslots from (τ − ρ)T to τT F (D(ρ, τ)) : forecasting model d ∈ R T + : generic demand vector d(τ) ∈ R T + : observed demand vector of planning period τd (τ) ∈ R T + : predicted demand vector of planning period τ σ (τ) ∈ R T + : predicted demand standard deviation vector D (d(τ)) : deterministic version of the planning model R (d(τ), σ (τ),Γ (τ)) : robust counterpart of the planning model E (d(τ)) : D (d(τ)) with additional variables regarding overtime production, which are used to cover unfeasibilities with extra costs when observed demand must be met Ob j * {model} : objective function value at the optimizer of a specific model Pr : price of robustness equal to Ob j * {R} minus Ob j * {D} . Po : price of overtime production as detailed in the following P : sum of Pr and Po Dealing with a waste management setting, D is a mixed integer linear programming model to schedule the sorting operations of each phase of a waste sorting process. The model is presented in [11] and supports several strategic decisions that are critical in the considered application. It can be described as a variant of a lot size model with non linear costs (approximated by mean of piece-wise 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 indexes and variables is set out in the following. j = {1, . . . , J} : index of the J sorting stages p = {1, . . . , P} : index of the P time-shifts T : time horizon partitioned in time shifts with t ∈ {1, . . . , T } C : hourly cost of each operator b t : working hours for time t determined by the corresponding shift p C t = C * b t : cost of each operator at time t f j : set-up cost of sorting stage j d t : quantity of material in kg unloaded from trucks at time t α j : percentage of waste processed in stage j − 1, received in input by buffer j S j : maximum inventory capacity of the sorting stage buffer j LC j : critical stock level threshold of buffer j n j : fraction of material allowed to be left at buffer j at the end of time horizon K j : single operator hourly production capacity [kg/h] of sorting stage j SK j,t = K j * b t : operator sorting capacity in sorting stage j, at time t M : maximum number of operators available in each time shift E j : minimum number of operators to be employed in each time shift of stage j ∂ h i j : slope of the i-th part of linearization of the buffer j stock cost curve The model consider 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: stock 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, this binary variables are used to model the piece-wise linear functions of the buffer stock costs.
Considering a case study where J = 2 sorting stages, for the 1 st sorting phase, u 1,t ≥ 0 and y 1,t ∈ {0, 1} represent the quantity (in kg) of material to be selected and decision to activate the process respectively at time t. For the 2 st sorting phase, u 2,t ≥ 0 and y 2,t ∈ {0, 1} represent the quantity (in kg) of material to be selected and decision to activate the process respectively at time t. I 1,t , I ′ 1,t , I ′′ 1,t ≥ 0 are the inventory levels at 1 st phase sorting buffer while I 2,t , I ′ 2,t , I ′′ 2,t ≥ 0 are inventory levels at 2 nd phase sorting buffer. As previously stated, w 1 and w 2 are used to model the piece-wise linear functions of the buffer stock costs. In detail The model minimizes the sum of sorting and holding costs and is detailed as following: The objective function (1) defines the minimization of the sum of the three cost terms, which are sorting, setup, and inventory costs respectively. (2) and (3) bounds the number of workers that can be assigned to each sorting station and to each time shift. Constraints (4) 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 (5) defines the inventory for the first buffer, considering the inbound material d t and the sorted material u 1t , while (6) defines the inventory for the other buffers corresponding to j > 1. Indeed, constraints (6) describe the waste flow across the sorting stages that follow one another: each subsequent inter-operational 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 (7), (8), and (9) define the piece-wise 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 (10) imposes the maximum unsorted material allowed to be left at the end of the planning period for each buffer.
The mixed-integer linear programming model R defines the robust counterpart of D and is newly introduced in [12]. We consider the set of parameters d t , t ∈ T , that are subject to uncertainty taking values according to a symmetric distribution with mean equal to the nominal value d t in the interval Indeed σ t is the maximum deviation of d t . In order to meet the standard formulation of the nominal problem presented in [1], where parameters subject to uncertainties belong to inequality constraints only, the equality constraints of D regarding waste arrivals d t are reformulated to turn them into inequality constraints. This is performed considering, for each period t, the sum of all the received and processed quantities of waste up to that period, as in constraints (14), (15), (16), (17) of the formulation of R presented in this section. According to the robust approach presented in [1], 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, Γ represents the number of coefficients that we consider as allowed to vary within their interval, ergo we consider nature behaving like only a subset of the coefficients will change with respect to their nominal value. Indeed, as affirmed in [1], 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 then that by varying Γ the level of robustness can be flexibly adjusted against the level of conservatism of the solution. Considering the peculiar 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 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 [1] and formulate the robust counterpart of D: ε t ∈ R + : extra variables multiplying d t ∀t ∈ T . These variables are introduced in order to have a variable multiplying the only set of parameters that are affected by uncertainty. These are indeed constrained to be equal to 1 ∀t ∈ T as in constraint (20) of R. z t ∈ R + : variable resulted of duality within Bertsimas and Sim robustness theory; when multiplied by Γ t provides its overall contribution to the protection function of constraint t. p t,k ∈ R + : variable resulted of duality within Bertsimas and Sim robustness theory; provides its contribution to the protection function of constraint t with respect to the specific coefficient d k . s t ∈ R + : variable resulted of duality and Bertsimas and Sim robustness theory; multiplied by σ t sets 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 following: min (1) s.t.
(2), (3), (4), (6), (7), (8), (9), (10), (11), (12), (13) ∀t ∈ T (20) ∀t ∈ T Constraint sets (14)(15)(16)(17) define and limit the inventories: constraint (14) defines the inventory for the first buffer, considering the cumulative inbound material a t up to period t, the overall sorted material u 1t 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 (15) sets the lower bound of the inventory for each period and (16) imposes the maximum unsorted material allowed to be left at the end of the planning period for the first buffer, as constraint (10) does for all other subsequent buffers. Equality constraint (17) allows the inventory of the main buffer (i.e. buffer no. 1) to be considered in the corresponding piecewise linear part of the cost function. Constraints (18) and (19) resulted from duality in [1] robustness theory; where (18) sets the lower bound of the protection function contribution in constraints (14) and (16).
Both D and R are always solved considering the predicted demand vectord. Instead, E replicates the same deterministic formulation of D and introduces some additional decision variables to model overtime production in order to capture the extra costs resulting whenever the robust solution is not feasible when the true demand vectord is observed. Indeed instances of E formulation are solved considering the observed demand vectord while fixing the first set of standard production capacity decision variables (i.e. y j,t ∈ {0, 1} and x j,t ∈ Z + ) to the values of the robust solution of R for the same time scheduling time window. Considering 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 + ) are 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 following: s.t.

Robustness control optimization
The robustness control vector must be learned by a routine properly build in order to achieve this parameter tuning objective. In the following we present an algorithm that generates a single evaluation of P according to a specific protection Γ . This algorithm, that we call Eval P, is part of any routine aimed at the optimization of the robustness control vector Γ .
1. Set planning time horizon as the number T of timeslots of a planning period τ 2. Set and consider ρ planning periods as the training data part of demand time-series D(ρ, τ) and get observed demand vectord(τ) accordingly 3. Let Ob j(·) the objective function operator and Ob j * {model} the objective function value at the optimizer of a specific planning model of the three presented 4. Solve F (D(ρ, τ)) to forecastd(τ) and σ (τ) 5. Set robustness protection logic and corresponding Γ 6. Solve R (d(τ), σ (τ),Γ (τ)) to obtain y * , x * and Ob j * . Solve E d (τ) where y and x are fix to y * and x * respectively to obtain y ′ * and x ′ * 10. Get Po(Γ ) = f ′ y ′ * +C ′ x ′ * 11. Get P(Γ ) = Pr + Po The optimization of Γ is intended to find such a Γ that minimizes P, that is the sum of the price of robustness and any future potential extra costs resulting from overestimating or underestimating risk. With this objective we propose a routine that seeks for the best configuration of Γ for a considered planning period τ of T timeslots with its corresponding demand forecastingd(τ) and its true observationsd(τ). Its pseudo-code is set out below.
1. Set Γ 0 ∈ R T + : the routine starting configuration of Γ 2. Set Γ step ∈ R T + : a vector of the updating step sizes of each component Γ t of Γ where t ∈ {1, . . . , T } 3. Set ρ and getd(τ) from D(ρ, τ) 4. Solve F (D(ρ, τ)) to forecastd(τ) and σ (τ) 5. Set max k = f loor (1/Γ step) : maximum number of iterations k for a complete search over the space of Γ ∈ R T We refer to this routine as Optimize Γ . Whatever setting of Γ 0 and Γ step is made, this must be done according to the peculiar structure of constraints (14), (15), (16), (17) of the formulation of R. Indeed, as stated in section 2 before introducing model R, the telescopic expansion of each set J t (i.e. |J t+1 | = |J t | + 1) constraints each component Γ t of the protection vector Γ to be bigger or equal to the component Γ t−1 and smaller than t. The setting of both Γ 0 and Γ step must be consistent with this 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 a new proposed framework integrating demand forecasting techniques. The scope of this framework is providing automatic decision support for implementing a proper robustness control vector that minimizes the aforementioned costs resulting from overestimating or underestimating risk in planning problems under uncertainties. Referring to the immediate past of known demand realization and its corresponding time-series as D(ρ, τ), we consider τ planning period as the latest T time slots part of D(ρ, τ) that have just ended. This demand time-series D(ρ, τ) can be used along with the 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(ρ, τ) The proposed framework essence is using D(ρ, τ) and its corresponding best configuration Γ * of the very close past for solving the planning problem of the future planning period τ + 1. Therefore, for each observed demand time-series D(ρ, τ), the corresponding Γ * effect 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 Section 3. Performances in terms of P costs are compared with P resulting from either the deterministic and worste-case scenario approaches that we refer to as P(Γ nominal ) and P(Γ worstecase ) respectively.

Forecasting model
We describe a time series forecasting model presented in [17] and designed to handle the common features of time series seen both in business and industrial scenarios. The implementation is available as open source software, is called Prophet, and it uses a decomposable time series model [7] with three main model components: trend, seasonality, and holidays. They are combined in the following equation: where g(t) is the trend function which models non-periodic changes in the value of the time series, s(t) represents periodic changes such as seasonality, and h(t) represents the effects of holidays which occur on potentially irregular schedules. The error term ε t represents any changes which are not captured by the model and is assumed to be normally distributed. Therefore, seasonality and potential holidays are treated as additive components. The model is identified within this class of models by solving a curve-fitting problem, either using backfitting or L-BFGS [4]. This is a different approach from time series models that explicitly account for the temporal dependence structure in the data.
Regarding the trend component g(t), a piece-wise constant trend provides a parsimonious and often useful model. In addition, eventual trend changes in the model are modeled with changepoints where the trend is allowed to change. The trend model g(t) is indeed considering S changepoints at times s j , ( j ∈ 1, ..., S) and a vector δ ∈ R S , where δ j is the change in rate that occurs at time s j . The rate at any time t is then the base rate k plus all of the adjustments up to that point: k + ∑ j:t>s j δ j . This is better represented by considering a vector a(t) ∈ {0, 1} S such that: The trend rate at time t is then k + a(t) ⊺ δ . The trend component g(t) considers also an offset parameter m that must be adjusted whenever the rate k is adjusted in order to connect the endpoints of the segments. The adjustment at changepoint j is computed as γ j = −s j δ j . Therefore, the trend component g(t) can be defined as following: Many time series have multi-period seasonality producing evident repeated effects that should be properly addressed by a forecasting model. The model relies on a Fourier series to provide a flexible model of periodic effects [8]. Let Q be the expected regular period the the time-series (e.g. Q = 7 for weekly data), then the smooth seasonal effects component can be approximated with the a standard Fourier series: Fitting s(t) requires the estimation of the 2N parameters β = [a 1 , b 1 , ..., a N , b N ] ⊺ . This is done by firstly setting N and the period Q, then constructing a matrix of seasonality vectors X(t) for each value of t in both historical and future data: By doing so the seasonal component is then The authors suggest to impose a prior smoothing on the seasonality by taking β ∼ Normal(0, σ 2 ).
Incorporating holidays effects with the corresponding component h(t) is made by firstly assuming that these effects are independent. Considering each holiday i, D i is the set of past and future dates for that specific holiday and an indicator function ✶(t ∈ D i ) is defined in order to represent whether time t is during holiday i. Then, a parameter k i is assigned to each holiday being the corresponding change in the forecast. As similarly done for seasonality s(t), a matrix of regressors is defined as follows, where L is the number of yearly holidays considered: By doing so the holidays component is As with seasonality, he authors suggest to use a prior k ∼ Normal(0, σ 2 ).

Experimental results
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 worste-case scenario approaches. Results are shown in Table 1 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 worste-case P(Γ worstecase ) approaches. The presented values of costs reduction result from the following expressions: The follwing columns of Table 1 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 gamma 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 robustness 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 behaviour of forecast models is crucial for a proper planning. Results of Table 1 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 is providing 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 worste-case approaches. Figure 1 shows ad example of Optimize Γ instance solving. As described at the end of Section 2.1, 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 related performances concerning 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 2. As in Table 1, performances in terms of P(Γ * ) costs are compared with P resulting from either the deterministic P(Γ nominal ) and worst-case P(Γ worstecase ) 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 timeslots considering 6 weekly working days of 2 working shifts each). Table 2 highlight 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 worste-case approaches. Indeed, out of the 30 instances where F underestimates the 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 overestimate 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(Γ worstecase ).

Conclusions
In this paper we analyzed the problem of how to considering extra costs resulting from uncertain demand in waste management. Robust optimization theory is applied to a mathematical programming model optimizing the planning of recycling operations. The demand forecasts and 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 workers. The framework is tested in a real case environment where demand does not follow a specific probability distribution. The results shown that the most of the time, the classical robust approach is over conservative. Moreover, when the forecasts over estimates the observed demand, the robustness budget Gamma can be effectively optimized. Fore-  casting model features and their overall repercussion over the protection cost are not addressed by this study and they will in a future work. Moreover, the proposed approach can be applied to different planning problems and tested against different scenario in future works.