An adjustable sample average approximation algorithm for the stochastic production‐inventory‐routing problem

We consider a stochastic single item production‐inventory‐routing problem with a single producer, multiple clients, and multiple vehicles. At the clients, demand is allowed to be backlogged incurring a penalty cost. Demands are considered uncertain. A recourse model is presented, and valid inequalities are introduced to enhance the model. A new general approach that explores the sample average approximation (SAA) method is introduced. In the sample average approximation method, several sample sets are generated and solved independently in order to obtain a set of candidate solutions. Then, the candidate solutions are tested on a larger sample, and the best solution is selected among the candidates. In contrast to this approach, called static, we propose an adjustable approach that explores the candidate solutions in order to identify common structures. Using that information, part of the first‐stage decision variables is fixed, and the resulting restricted problem is solved for a larger size sample. Several heuristic algorithms based on the mathematical model are considered within each approach. Computational tests based on randomly generated instances are conducted to test several variants of the two approaches. The results show that the new adjustable SAA heuristic performs better than the static one for most of the instances.

by companies, and to ensure that a two-stage recourse problem has complete recourse (for each first-stage solution there is a feasible assignment of the second-stage solutions). However, as reported in [8], penalizing the backlogged demand makes the instance harder to solve by exact methods based on mixed integer formulations, since the integrality gaps of the linear solutions tend to be high. Hence, these SPIR problems tend to be even harder than the deterministic ones when no backlog is allowed.
When the SAA problems are not solved to optimality for each sample, theoretical results on convergence fail and the SAA method acts as a heuristic. Based on this observation, and on the assumption that one cannot solve efficiently these complex SPIR problems even for small size instances, a natural approach is to follow the SAA method and solve each one of the SAA problems defined for small size samples heuristically. Then, the best solution among the generated solutions is chosen. We call this the static approach. In contrast with this approach, in which the solution to the problem is one of the candidate solutions, we propose a new heuristic approach that explores the knowledge of the several candidate solutions in order to solve the SAA problem for a larger and more representative set of scenarios. Following this approach, a new heuristic named Adjustable Sample Average Approximation (ASAA) algorithm is proposed. In this algorithm, a feasible solution for each sample is generated using a heuristic scheme. Then, these solutions are explored in order to identify first-stage variables that are frequently fixed to zero or to one, and a partial fixed solution is constructed, since only some of those variables are fixed. Finally, the restricted model is solved for a larger sample keeping all the remaining variables free. The proposed adjustable heuristic approach has several advantages in relation to the classical SAA method: (i) it does not require solving each sample set to optimality; (ii) it does not require the use of large sample sets, and (iii) it allows the final solution to be adjusted to a larger sample set since many first-stage variables are kept free. Further, the number of variables that are fixed can be easily controlled.
To the best of our knowledge, the use of the information from several solutions obtained with the SAA method to obtain a partial solution was introduced in a preliminary version of this work [10]. However, algorithms that use information of previous solutions to derive partial solutions are not new. Rochet and Taillard [35] introduce the Adaptive Memory concept to describe a set of good solutions, that are kept in a pool that is dynamically updated taking into account common features of the solutions. Based on this concept, a very general metaheuristic focused on the exploitation of strategic memory components was proposed by Glover [22]. Such a metaheuristic, known as Adaptive Memory Programming (AMP), has been applied over the years to solve several hard-combinatorial optimization problems, such as vehicle routing problems [23,41], stochastic production distribution network design [15], and supplier selection problems [42]. Another closely related concept is the Concentration Set (CS) introduced by Rosing and ReVelle [36] to solve the deterministic p -median problem heuristically. A large number of solutions are generated heuristically. Then, the best υ solutions are selected to the CS and this set is used to set to one (zero) all the variables that take value one (zero) in all the υ solutions. The remaining variables are kept free. Threshold values were also used in [26] for a dynamic facility location problem with capacity levels. Using a Lagrangian relaxation model, some of those levels are eliminated (the corresponding variables are fixed to zero) according to the information provided by the set of solutions.
In this paper, we introduce a general approach that explores the SAA in order to handle large size instances of the productioninventory-routing problem. We discuss a mixed integer model and possible enhancements. Based on the enhanced model with valid inequalities, two main approaches are compared: a static approach and an adjustable SAA approach. In order to improve solutions, an iterated local search heuristic based on local branching is presented. The adjustable approach can be easily extended for other two-stage stochastic problems where the corresponding deterministic problems are already difficult to solve to optimality for realistic size instances.
The rest of this paper is organized as follows. The mathematical model is presented in Section 2. Model improvements, such as the valid inequalities and an extended formulation, are discussed in Section 3. The SAA method, the static heuristic approaches based on the SAA method, and the new ASAA heuristic are introduced in Section 4. Computational tests that compare variants of the static and the adjustable approaches are presented in Section 5. Final conclusions are drawn in Section 6.

P R O B L E M S P E C I F I C A T I O N S A N D M A T H E M A T I C A L F O R M U L A T I O N
In this section, we introduce the mathematical model for the SPIR problem. We assume the production plan (production periods and quantities to produce) and the routing (which clients to visit in each period) are first-stage decisions. The quantities to deliver, the amount backlogged at each time period as well as the inventory levels are adjusted to the scenario. The goal of the stochastic approach is to find the solution that minimizes the production and the routing cost plus the expected cost of both the inventory and the penalty costs for backlogged demand. Following the SAA method, the true expected cost value is replaced by the mean value of a large random sample = ξ 1 , . . . , ξ s of scenarios, obtained by the Monte Carlo method. This larger set of s scenarios is regarded as a benchmark scenario set representing the true distribution [27].
Consider a graph G = (N, A) where N = 0, 1, . . . , n N represents the set of nodes and A represents the set of possible links. Set A is a subset of set N × N. Node 0 denotes the producer and N c = 1, . . . , n N is the set of clients. Set T = 1, . . . , n T is the set of periods and K = 1, . . . , n V is the set of homogeneous vehicles available.
Consider the following parameters: d it (ξ ) is the demand of client i ∈ N c in period t ∈ T in scenario ξ ∈ , I 0 /I i is the initial stock at producer/client i,S 0 /S i is the inventory capacity at producer/client i,P is the production capacity at each time period, Q it andQ it are the lower and upper limits to the delivery quantity in period t ∈ T at client i and L is the vehicle capacity. For the objective function, S is the set up cost for producing in a period, P is the production cost of a unit of product, V is the fixed vehicle usage cost, C ij is the travel cost from node i to node j, (i, j) ∈ A, H it is the holding cost of the product at node i ∈ N in period t ∈ T , B it is the backlog cost of the product at node i ∈ N in period t ∈ T .
Consider the following first-stage variables: binary variables y t are the setup variables that indicate whether there is production in period t ∈ T , variables p t give the production level in period t ∈ T , binary variables z itk indicate whether there is a visit to node i ∈ N c in period t ∈ T , by vehicle k ∈ K, the routing variables x ijtk indicate whether vehicle k ∈ K travels from node i to node j, (i, j) ∈ A, in period t ∈ T , v tk is a binary variable that is one if vehicle k ∈ K is used in period t ∈ T and zero otherwise; for (i, j) ∈ A, t ∈ T , k ∈ K, f ijtk is the artificial flow of a single commodity variable used to prevent cycles in the routing problem. Notice that this is not the amount transported from node i to node j in period t ∈ T by vehicle k since that quantity depends on the scenario. Such adjustable variables could be used but would imply the use of an unnecessarily large model.
For each scenario ξ ∈ , we define the second-stage variables q itk (ξ ) indicating the quantity delivered at node i ∈ N c in period t ∈ T by vehicle k; s it (ξ ) indicating the stock level of node i ∈ N at the end of period t ∈ T ; and r it (ξ ) is the quantity backlogged at node i ∈ N c in period t ∈ T .
The SPIR problem is modeled by the following formulation.
The objective function (1) minimizes the total cost which includes the production set up, the production, the vehicle usage, the traveling costs, and the expected holding and backlog penalty costs. Constraints (2) and (3) are the inventory conservation constraints at the producer and at the clients, respectively. Constraints (4) establish the initial inventory level at the producer and at each client. Constraints (5) impose a storage capacity at the producer and at each client. Constraints (6) impose a maximum capacity on the production at each period. Constraints (7) impose limits on the delivery quantity at each client in each period. The lower bound imposed ensures that for each scenario the delivered quantity is positive, otherwise there could be second-stage solutions with null delivery quantity. Constraints (8) establish a capacity on the quantity transported by a vehicle. Constraints (9) and (10) are the routing constraints. Constraints (11) guarantee that a vehicle leaves the producer when there are visits to clients. Constraints (12) are the flow balance constraints at clients ensuring that visits are satisfied and constraints (13) guarantee that there is flow between two nodes only if a vehicle travels between these two nodes. Constraints (14), (15) and (16)

M O D E L T I G H T E N I N G
The mixed integer model is weak in the sense that the integrality gap is usually large [2,6]. In order to tighten the model, we discuss families of valid inequalities and the use of an extended formulation.

Extended formulation
In this subsection, we present an extended formulation which results from adapting ideas from the extended formulations for related problems (see [32] for deterministic lot-sizing problems and [4,37] for deterministic production-inventory-routing problems without backlog) to the production-inventory-routing problem with multiple scenarios. One weakness of model SPIRF are the linking constraints (7), especially when the delivered quantities q ikt (ξ ) can be much smaller than the upper boundsQ it .
By decomposing the amount delivered at each client in a given time period into the net demand of the time period destination, tighter linking constraints can be developed.
Let nd it (ξ ) denote the net demand for client i ∈ N c at time period t ∈ T when scenario ξ ∈ occurs. Such net demand represents the minimum amount of product that must be delivered to client i at period t in order to satisfy the client demand taking into account the initial stock level, and can be computed as follows: Consider the new variables u it (ξ ) indicating the fraction of the net demand of client i ∈ N c at time period t that is delivered in time period when scenario ξ occurs, for t ∈ T , ∈ T ∪ n N + 1 , ξ ∈ . Notice that when = n T + 1, variables u it (ξ ) represent an amount that is not delivered during the considered time horizon.
In order to tighten model SPIRF, the following set of constraints, that we call EF, can be added.
Constraints (17) relate variables u it (ξ ) to the delivery quantity variables q itk (ξ ) and compute the quantity that is delivered to client i at time period . Constraints (18) ensure that if some fraction of net demand for client i is delivered at time period , then client i must be visited in that time period. Constraints (19) compute the demand not satisfied. Constraints (20) guarantee that demand is satisfied by production and/or is unmet during the time horizon (when u i,t,n T +1 (ξ ) > 0). Constraints (21) define the new variables as nonnegative.
The linking constraints (18) of this extended formulation are tighter than the linking constraints (7) of the original model, leading to a tighter model.
Even considering a small number of scenarios, the number of constraints and variables in (17)-(21) becomes too large, which means that it is computationally impractical to use such sets of variables and constraints with all scenarios. In Section 5.1, some strategies are proposed to overcome this problem.

Valid inequalities
From the underlying fixed-charge network flow structure of the SPIRF one can derive several families of valid inequalities [31]. Since the possible set of inequalities is large, we select several such families of inequalities based on preliminary computational tests and on the knowledge of its relevance to related problems. Some of these inequalities are adapted from deterministic lot-sizing problems [32].
The first family relates the delivered quantity k∈K i∈Nc q itk (ξ ) with the setups. If there is no set up in time period t, then the delivered amount is bounded by the stock in the producer at the end of time period t − 1. Otherwise, k∈K i∈Nc q itk (ξ ) is bounded by the total vehicles capacities n V L. k∈K i∈Nc To the best of our knowledge, inequalities (22) are new, even for the deterministic case (when only one scenario is considered). These inequalities are useful to cut linear solutions with a fractional value of y t . Their impact is greater when n V L is small, which happens, for instance, when a single small capacity vehicle is used, see Section 5.1.
The following two families of inequalities link the inventory levels (stock and backlog) with the binary variables representing either setups or the visits to clients. We notice that the structure of these families of inequalities is quite different from the ones discussed in [4] for a related problem since backlog is not considered there.
The first family is useful to cut fractional solutions when the aggregated value of the binary variables in the RHS of the inequalities is less than 1. This family of inequalities is derived from the uncapacitated lot-sizing with backlogging feasible set (see [32]) defined by inequalities (3), (4), and (7), for a client i, period t, and scenario ξ : where(x) + = max {x, 0} . These inequalities force the net demand (demand that cannot be satisfied by the initial stock) at client i until period t to be satisfied with backlog in case there is no delivery until period t, that is, if k∈K t =1 z i k = 0. Note that these inequalities are a particular case of the more general family The second family of inequalities is derived from the uncapacitated lot-sizing with backlogging feasible set obtained by considering the aggregated demands, stocks, and backlogs.
Inequalities (24) ensure that the total demand from period h to period t must be satisfied either from stock (at the producer or at the clients) or from backlog at the clients if there is no production in that period.
For the particular case of h = 1, inequalities (24) can be written in the stronger format , is the net demand until period t for scenario ξ . Next, we introduce a family of inequalities that can be derived by Mixed Integer Rounding (MIR) [30]. Given a simple mixed integer set where a, b ∈ R + are arbitrary constants, the MIR inequality is given by is a relaxation of the feasible solutions set and imposes that the total production capacity plus the backlogged demand must cover the net demand of all clients until period t. (25) state that the total setup capacity plus the backlogged demand cover the net demand of all clients until period t.
From preliminary computational experience, we observed that the number of inequalities cutting off the linear relaxation solutions is too large. Hence some strategies need to be taken in order to choose which cuts should be included. Two main strategies, Option (W) and Option (A), that differ in the way the scenarios are considered, were tested.

(W) Worst case scenario
In this option, for each family of inequalities and each set of time period and/or client indices (accordingly to the family) at most one inequality is added. The inequality added is the one that, of all the scenarios in , corresponds to the scenario leading to the largest violated inequality.

(A) Aggregation of scenarios
In this option, for each family of inequalities and each set of time period and/or client indices (accordingly to the family), the set of inequalities is aggregated for all the scenarios in . It adds the inequality aggregating all the scenarios leading to a violated inequality.
Option (W) tends to add too many cuts since, for a given fractional first-stage solution, the value of the second-stage solution variables, variables s it (ξ ) and r it (ξ ), tends to be null. Moreover, inequalities for a given scenario do not necessarily affect the solution variables of the remaining scenarios (unless they force the values of the first-stage solutions to change). Option (A) tends to add fewer cuts but these cuts are weaker than the individual inequalities for each scenario. This option resembles the scenario-group cuts approach proposed by Adulyasak et al. [4] where groups of scenarios are considered accordingly to the total demand, and Benders cuts are aggregated for each group. Option (A) corresponds to the case where only one group is created which, in the general framework we are proposing, is reasonable since the number of scenarios considered is small.
In Section 5, we report results derived from the tests conducted to assess the performance of these two options.

S O L U T I O N A P P R O A C H E S
In this section, we discuss several algorithms to solve the SPIR problem. These algorithms are based on the SPIRF model which is a two-stage stochastic model. A common approach to solve stochastic problems is the Sample Average Approximation (SAA) method [28,43], where the true expected value of the objective function is approximated by the average value obtained for a very large sample of s scenarios. The SAA method generates M separate sample sets m , m ∈ {1, . . . , M}, each one containing s scenarios. For each one of the scenarios set, m , the resulting SAA problem (where is replaced by m in model SPIRF) is solved generating a candidate solution. Then, for each candidate solution, the first-stage solution is fixed, and the value of the objective function for a very large sample with s scenarios is computed. In the case of the two-stage model SPIRF, this value is computed by solving a pure linear programming problem on the second-stage variables.
Solving the SPIR problem to optimality, even for a very small size scenarios sample, may not be practical, given the difficulty of the problem. Therefore the M subproblems generated by the SAA method can hardly be solved to optimality. Hence, the candidate solutions will be feasible solutions obtained with heuristics.
Next, we describe two approaches to solve the SPIR problem. The first one follows the SAA method and uses heuristics to obtain the M candidate solutions. Then the best one is chosen. This will be named the static approach. The second approach uses the candidate solutions to fix part of the first-stage solution and then solves (heuristically) the resulting restricted problem for a larger sample. Since the final solution can be adjusted to the larger sample we name it the adjustable approach.

Static approach
A first heuristic consists in solving the SPIRF model by branch-and-bound with a commercial solver and imposing a running time limit. In fact, as we can observe in the computational section, this is a natural option that works well on the easiest instances where the solver can find near optimal solutions within reasonable time. For those cases, imposing a time limit does not have great impact on the quality of the solutions. Since it is well known, a large amount of running time is used to prove optimality, which is not a major relevant issue here given that the solution will be evaluated with another sample.
When the instances are harder to solve, the previous approach can perform quite badly and the quality of the solutions varies from sample to sample, leading to a large variability in the solution values. In order to circumvent this fault, we used two heuristic approaches that can be taken separately or combined. One approach is to use an iterative local search method to improve the quality of the candidate solutions. The other approach is to simplify the model by predefining a route, making the restricted model easier to solve.

Iterative local search method
The quality of any feasible solution can be improved using a Local Search procedure. The local search scheme applies a local branching method based on the method proposed by Fischetti and Lodi [21] for reducing the solution space. For a given positive integer parameter z , define the neighborhood N (z, z ) ofz as the set of feasible solutions of the SPIRF satisfying the additional local branching constraint: Hence, the neighborhood ofz, N (z, z ), is the set of solutions that differ by a maximum of z values from the z itk variables of the current solutionz. In fact, the linear constraint (26) limits to z the total number of binary variables z itk flipping their value with respect to the solutionz, either from 1 to 0 or from 0 to 1. The Local Branching heuristic amounts to adding constraint (26) to SPIRF and running the solver for a time limit of α seconds. If a better solution is found, thenz is updated and the Local Branching is run again. The complete heuristic is denoted by Iterated Local Search (ILS) and is described in Algorithm 1.

Using routing precedences from a defined route
The SPIRF includes the TSP as a subproblem. For each time period, the minimum cost hamiltonian cycle that includes all the visited clients must be determined. The TSP is a NP-hard problem and is often very time-consuming. Combining the TSP formulation with production and inventory makes the model SPIRF quite large. In order to simplify the model, a TSP is initially solved for all clients. The order in which clients are visited on this route is used to define a precedence relation between the clients. These precedences are used to eliminate (set to zero) those variables that violate the precedence relation. Similar heuristic approaches have been used before for deterministic problems [38,40]. The TSP involving all the clients is solved using the Concorde software [18].
We name this simplification procedure as the Precedence order relation (P) technique.

Adjustable approach
In this subsection, we describe the Adjustable Sample Average Approximation (ASAA) heuristic. Since this heuristic can easily be extended to other two-stage stochastic problems we introduce it as a general framework procedure.
In contrast with the SAA method, where a set of M candidate solutions are generated by solving M SAA problems and the corresponding first-stage solutions are fixed, we propose a new approach that uses the M candidate solutions to generate a partial solution. The idea is to identify within the first-stage variables those that have the same value in all (or almost all) the M generated solutions. Thus, it is likely that this value coincides with the variable value in the optimal solution. By fixing the value of those first-stage variables, the problem is simplified and the solution can be completed by solving a problem with a new larger sample.
It is important that the M candidate solutions result from different solution structures, so that the variable fixing corresponds to common structures under different sample scenarios. Hence, the number of scenarios considered in each sample can be small which facilitates the generation of the M candidate solutions. Our experience (not reported in the computational tests) show that for the SPIR problem working with only one scenario leads to solutions that are quite different and few first-stage variables are fixed. Thus, the number of scenarios should not be too small.
Another important ingredient is to decide which variables to fix. If a binary first-stage variable has the same value in all the candidate solutions it is natural to fix it to that value. However, we may wish to fix variables that have the same value in at least a given percentage of solutions or satisfy another criterion. Moreover, we may wish to assign different weights to each candidate solution either depending on the probability (if the solution is obtained for a scenario with higher probability its weight could be higher) or on its objective function value (a solution with a better objective function value than the others would be given a greater weight).
Hence Vectors 0 M and 1 M represent the null vector and the vector with all components equal to 1 of dimension M, respectively. Condition 1 ensures that if a variable ω is zero (one) in all the candidate solutions, then the value of the function will be zero (one).
Let be a subset of binary first-stage decision variables that makes the problem easier to solve when their values are fixed. Thus, this new method can be described in four steps.
First step: solve M independent problems. Generate M separate samples m , m ∈ {1, . . . , M}, of dimension 1 . Then, solve the two-stage stochastic model for each sample set, m , with a time limit of α seconds (either a feasible solution is found within this time limit or it stops after the first feasible solution is found). This gives M candidate solutions, ϕ 1 , . . . , ϕ M .
Second step: obtainM partial solutions. Use the candidate solutions ϕ 1 , . . . , ϕ M to obtain partial candidate solutions by fixing the value of some first-stage variables in set and keeping all the remaining variables free.
Let ϕ P m denote the projection of the solution vector ϕ m into the space of variables in . Define two threshold values, γ 1 and γ 2 and generate a new partial solution in the following way: denote by (ϕ P m ) τ the value of a generic first-stage decision variable τ from set in its partial solution ϕ P m , with the value of this variable in the generated partial solution, ψ P , computed as follows Either consideringM different functions f I or varying the parameters γ 1 and γ 2 ,M different partial solutions ψ P1 , . . . , ψ PM can be obtained.
Third step: complete theM partial solutions. For each partial candidate solution ψ P1 , . . . , ψ PM constructed and having some of the first-stage variables fixed, solve the two-stage stochastic model using a larger sample set with 2 ≥ 1 scenarios. Use the optimal solution obtained and fix to its optimal value the first-stage variables not fixed yet. Repeating the process for each partial candidate solution yieldsM candidate solutions with all the first-stage variables fixed.
Fourth step: solveM simplified problems. With all the first-stage decision variables fixed, the two-stage stochastic model is solved for a very large sample set with 3 2 scenarios. The value of the recourse variables (the second-stage variables) is computed for each scenario and, consequently, the objective function value of each solution is computed. The solution with the lowest average cost is chosen.
The second and the third steps correspond to the adjustable part of the procedure. If we let M =M and remove the second and the third steps, the adjustable approach ASAA coincides with the SAA method.
In the computational tests we considered the following aspects: • Only one partial solution is derived, that is,M = 1.
• The binary first-stage variables considered in set are the z ijk variables. • In the third step, the model for a large sample is solved by relaxing the routing variables, and all the first-stage variables are fixed, except the routing variables. Then, defining the routing variables as binary, a TSP is solved for each vehicle and each time period.

C O M P U T A T I O N A L E X P E R I M E N T S
This section reports the computational experiments carried out to illustrate the performance of the methods described to solve the SPIR problem. All tests were run using a computer with an Intel Core i7-4750HQ 2.00 GHz processor and 8 GB of RAM, and were conducted using the Xpress-Optimizer 28.01.04 solver with the default options. The performed experiments are based on the instances introduced in [6] which are generated as follows.
The coordinates of the n clients are randomly generated in a 100 by 100 square grid, and the producer is located in the center of the grid. For each value of n a complete graph G = (N, A) is obtained and symmetric traveling costs are associated to the set of arcs. The traveling costs C ij are the Euclidean distance between nodes i and j in the grid.
The computational experiments used to test the improvements of both the valid inequalities and the extended formulation are based on instances with n N = 5, 15 and 25 clients. For the computational experiments conducted to assess the performance of the ASAA, the number of clients used is n N = 10, 20, 30, 40, 50, 60, 70, 80. In both cases, n T = 5, 10 time periods are considered.
For each client and each period, the nominal demand value d it is randomly generated between 40 and 80 units, and the uncertain demands vary in [0.7d it , 1.3d it ].
The initial stock I 0 at the producer is randomly generated in the interval [0, 30], and the initial stock I i at client i is randomly generated between 0 and three times the average demand of client i. The maximum inventory levelS i is 500 for all i ∈ N. The production capacityP is 80% of the average net demand. The number n V of available homogeneous vehicles is one and two and their capacity L is set to 80% and to 40% of the average net demand, respectively. The lower Q it and upperQ it delivery limits are 1 and L (the vehicles capacity), respectively. The production set up cost, the unit production cost and the fixed vehicle usage cost are given by S = 100, P = 1, and V = 50, respectively. For all i ∈ N, t ∈ T , the holding cost H it is 0.2 and the backlog cost B it is 0.5, except for t = n T where B i,n T is 5, since this cost penalizes the demand of client i that is not satisfied during the planning horizon.

Model enhancements
In this subsection, we report the computational experiments conducted to test the model enhancements, namely the extended formulation and the valid inequalities. Based on these results the final model is chosen.

Inclusion of the extended formulation
First, we evaluate the improvements observed when the extended formulation EF, described in Subsection 3.1, is used. Several strategies can be used to deal with the scenarios in the extended formulation. On the one hand, when the number of scenarios used is large the size of the formulation becomes too large. On the other hand, by considering few scenarios the model becomes weaker. To evaluate advantages and disadvantages of the use of the extended formulation, four strategies were studied: (i) SPIRF, only SPIRF is used; (ii) SPIRF+EF, SPIRF is used together with EF considering all the scenarios; (iii) SPIRF+EF l , SPIRF is used together with EF for the scenario with the lowest value of total net demand, ξ l = argmin ξ ∈ i∈Nc t∈T nd it (ξ ) ; (iv) SPIRF+EF g , SPIRF is used together with EF for the scenario with the greatest value of total net demand, ξ g = argmax ξ ∈ i∈Nc t∈T nd it (ξ ) .
We consider instances, as described above, with n N = 5, 15 and 25 clients, n T = 5, 10 periods and a single vehicle. For each instance, 10 scenarios were considered and the time limit was set to 600 seconds. Table 1 summarizes the computational experience. For each one of the four strategies and for each instance, identified by the number n T of periods and the number n N of clients, Table 1 displays the values of the linear programming relaxation (LR), the lower bound (LB), the upper bound (UB), and two gaps, where BFS is the value of the best feasible solution obtained among the four strategies mentioned. The values of LB and UB are obtained at the end of the running time.
As expected, the formulation SPIRF + EF provides the best quality linear relaxation solutions, however the computational time is higher, and the lower and upper bounds (and consequently the gaps Gap 1 and Gap 2 ), are the worst. The differences are very small between strategies SPIRF, SPIRF + EF l , and SPIRF + EF g . Therefore, no clear conclusions can be drawn. Although SPIRF + EF l provides better results for small size instances, the results for n N = 25, n T = 10 seem to indicate that for the larger instances model SPIRF is better. Also, the gains in terms of the linear relaxation value indicate that the improvements obtained by including the complete model EF are minor when the size of the model increases.
Hence, since there is no clear benefit in the use of any one of the strategies SPIRF + EF, SPIRF + EF l , and SPIRF + EF g , mainly for the largest instances, and since the number of constraints and variables tend to increase rapidly as the size of the instances increase, in what follows the extended formulation is not included.

Inclusion of valid inequalities
Now, we assess the improvements observed when SPIRF is tightened with the inclusion of the valid inequalities described in Subsection 3.2.
Since inequalities (22) are new, an independent experiment was conducted on these inequalities for instances with n T = 5 periods, n N = 10 clients, and n V = 1, 2 vehicles. The linear gap of SPIRF (Gap L ) and the linear gap of SPIRF tightened with inequalities (22) (Gap T ) were computed for L, the vehicles capacity, which was set to 20%, 40%, 60%, and 80% of the average net demand. The average gap reduction obtained by the inclusion of the inequalities, computed by Gap L −Gap T Gap T × 100, was 27.5, 20.6, 10.6, and 11.1, respectively, for one vehicle, and 26.2, 17.8, 1.4, and 0, respectively, for two vehicles. Thus, the impact of these inequalities is greater for small values of L and for one vehicle.
The valid inequalities (22)- (25) are dynamically added to the SPIRF and three strategies are compared: (i) SPIRF, only SPIRF is used, no inequalities are included.
(ii) SPIRF-W, corresponds to Option (W) for the inclusion of valid inequalities discussed in Section 3.2. (iii) SPIRF-A, this strategy corresponds to Option (A) applied to inequalities (22), (23), and (24) and Option (W) applied to inequalities (25). Notice that Option (A), aggregation of the inequalities for all the scenarios, applied to inequalities (25) leads to very weak inequalities.
Twelve sets of instances were used, six considering one vehicle and six considering two vehicles. For each set, five instances were generated considering different samples of demand values. A time limit of 600 seconds was imposed. The obtained results are displayed in Table 2, Figures 1 and 2. Table 2 displays the results for the three strategies: SPIRF, SPIRF-A, and SPIRF-W. Each line has the results for an instance set identified by its number n V of vehicles, its number n T of periods, and its number n N of clients. Columns identified with #vi display the average number of valid inequalities added to the SPIRF. Columns identified with #nodes display the average number  Table 2 shows that the average number of inequalities added to the SPIRF in strategy SPIRF-A is lower than the average number of inequalities added to the SPIRF in strategy SPIRF-W. Also, the number of nodes solved is lower. Columns Gap show that the average gap is reduced by the inclusion of valid inequalities and smaller average gaps are always obtained with strategy SPIRF-A. It is also possible to observe that instances with more than one vehicle are more difficult to solve, even for small size instances. Figures 1 and 2 display the relation between the average linear relaxation, average lower bound, and average upper bound values for the three strategies, for each instance set. Since these values can vary considerably from instance to instance, for better comparison each value was divided by the best (lower) upper bound of the corresponding instance set. Figure 1 depicts results for the single-vehicle case, while Figure 2 depicts results for the case where two vehicles are available. In both figures, the values obtained by using approaches SPIRF-A and SPIRF-W are associated with the solid and dashed lines, respectively. The dotted lines are associated with SPIRF, the case where no inequalities are added. The circle points (in the first/down set of three lines) represent the ratio between the average linear relaxation and the best upper bound value. Solid points (in the second/middle set of three lines) represent the ratio between the average lower bound and the best upper bound value. Square points (in the third/upper set of three lines) represent the ratio between the average upper bound and the best upper bound value. Figures 1 and 2 help to explain the gaps computed in Table 2, since the largest average lower bound and smallest average upper bound values are obtained with strategy SPIRF-A, except in one case. Note that, for the three approaches, the average lower bound values are very close.
Since the best results are obtained with strategy SPIRF-A, in what follows, the SPIRF will be used together with the inclusion of inequalities following strategy SPIRF-A.

Computational comparison of the static and adjustable strategies
Extensive computational experiments to assess the performance of the proposed Adjustable Sample Average Approximation method are presented. Three variants of the ASAA, that differ from each other in the strategy used to obtain the initial candidate solutions, are compared. These three variants of the ASAA method are: (i) ASAA ILS , the ILS procedure is used; (ii) ASAA ILS+P , the ILS procedure is used and precedence relations in the routing are considered; (iii) ASAA R , the routing variables x ijtk are relaxed.
These three variants were tested against the following four static approaches: (i) SAA, the SAA method; (ii) SAA ILS , the SAA with the ILS procedure; (iii) SAA ILS+P , the SAA with the ILS procedure in which precedence relations in the routing are considered; (iv) EVP, the Expected Value Problem (EVP), which corresponds to solving the deterministic problem with all the uncertainty parameters replaced by their expected values.
Whenever the ILS procedure is used, the initial scenario considered is the one that corresponds to the EVP. However, since we are dealing with difficult instances, usually the corresponding EVP problem cannot be solved to optimality within a reasonable amount of running time. Thus, a two-stage heuristic procedure to obtain a feasible solution is followed. First, the EVP with the routing variables x ijtk relaxed is solved with a time limit of 300 seconds. Then, all the visit variables z itk are fixed and n T pure TSP problems are solved by using the software Concorde. Even when the number of visit variables that are allowed to change their value in the local branching constraint (26) is low, the time used to solve each subproblem is large, since all the routing variables remain free. Hence, in order to solve such subproblems faster, we used a local branching constraint similar to (26) for the routing variables x ijtk , in which the number of variables that are allowed to flip their value is three times the value used in the local branching constraint associated to the visit variables.
To evaluate the performance of these seven strategies to solve the SPIRF, 32 instances, randomly generated as described before, are used. For each instance and each strategy, except for strategy named EVP (corresponding to solve the EVP), the number of candidate solutions used is M = 10 and each candidate solution is obtained by using a sample of 1 = 10 scenarios and a time limit of 300 seconds. However, when no integer solution is found within this time limit, the procedure only stops when the first integer solution is obtained. For strategy EVP, only one candidate solution is obtained using a single scenario, with a time limit of 3600 seconds (one hour). For all strategies, the dimension of the final sample used is 3 = 1000.
For the ASAA variants, a sample of * = 25 scenarios is used to define the weight of each partial candidate solution in the second step and several strategies for fixing variables z itk are used, corresponding to different choices of the parameters γ 1 and γ 2 . The threshold values 0, 0.05, 0.15, and 0.25 were tested for γ 1 , while γ 2 is set to 1 − γ 1 . In the third step, when solving the simplified model, a sample of 2 = 50 scenarios is used.
The obtained results are displayed in Table 3 and Table 4. Each strategy is identified in the first line of the table and each instance is identified by its number n V of vehicles, its number n T of periods and its number n N of clients in the first three columns.
In Table 3, the computational time (in seconds) required by each approach to obtain the final solution is displayed in columns named Time. For the adjustable approaches, the computational time values displayed correspond to the value obtained for γ 1 , the one which required more computational time. In columns VSS, an approximation of the Value of the Stochastic Solution is displayed for each approach considered. Notice that, as explained before, it is not possible to obtain the exact solution of the Expected Value Problem, thus all the values are computed according to the heuristic solution obtained for the EVP. For the adjustable approaches, the approximation for the value of the optimal solution is computed based on the best solution obtained for the different threshold values of γ 1 .
The computational times presented in Table 3 for all the approaches are very close and, in general, the lowest ones correspond to the approaches in which precedence relations are used. Furthermore, the additional computational time required to improve the solution in the adjustable approaches is lower than the time required to evaluate all the 10 candidate solutions in the final step of the corresponding static approaches.
The cost of the solutions obtained by each one of the strategies is reported in Table 4. In each set of four columns associated with each variant of the ASAA, the results for the different threshold values of γ 1 are reported.
The results displayed in Table 4 show that the cost of the solutions obtained by SAA tend to be excessively high for medium and large size instances. However, for the smallest instances, the cost of the solutions obtained by strategy ASAA R , where the structure of the problem changes due to the relaxation of the routing variables, is even higher than the one obtained by SAA.
The solutions obtained by the static strategies that use the ILS procedure, strategies SAA ILS and SAA ILS+P , are in general better than the ones obtained by solving the EVP. Furthermore, such solutions can be improved by using the corresponding adjustable strategies, ASAA ILS and ASAA ILS+P , mainly when the size of the instances increases. This means that the second and  Time  Cost  Time  Cost  VSS  Time  Cost  VSS   60  3710  27 285  2361  26 802  483  1531  26 590  695  5  70  3748  31 667  2874  30 442  1225  1740  30 204  1463  80  3789  35 110  3207  32 962  2148  1745  32 729  2381  1  60  4014  84 826  6718  67 447  17 379  3456  67 170  17 656  10  70  10 191  81 235  7215  73 633  7602  2809  73 290  7945  80  3973  96 860  6250  85 318  11 542  3316  84 981  11 879  60  3708  28 017  2534  27 388  629  1720  27 237  780  5  70  3745  33 008  3875  31 139  1869  2765  31  considering threshold values different from zero, which corresponds to fixing more variables than the ones that take the same value in all the candidate solutions. However, the threshold value γ 1 used should not be close to 0.5, since the quality of the solution can be lost when many variables are fixed a priori. From a general point of view, the computational results presented in Tables 3 and 4 allow us to conclude that in almost all the cases, the solutions obtained by the static approaches can be improved by the second and third steps of the adjustable approaches using a short running time. Hence, the proposed adjustable strategies prove to be effective and efficient, mainly for medium and large instances, since good solutions can be obtained in a short time. Moreover, among the three adjustable strategies tested, strategy ASAA ILS+P is the one that, in general, leads to the best results, solution and time.
Results not reported here show that each subproblem associated to the determination of each candidate solution is highly influenced by the variation of the demand values. Hence, in order to validate our conclusions, further tests were conducted for the large size instances in which the demand samples are derived from a normal distribution with parameters d it (mean) and 0.15d it (standard deviation), and for a gamma distribution with parameters (1/0.15) 2 (shape) and 1/(0.15 2 d it ) (rate). For these results we compare the EVP, the best static approach SAA ILS+P , and the best adjustable approach ASAA ILS+P using the threshold value 0.15. The results are reported in Tables 5 and 6, respectively.
Both tables show that the static and the adjustable approaches lead to gains when compared to the feasible solution obtained by solving the EVP. Although the running times are similar (which is expected since the first step in both approaches is the same and corresponds to computing the set of candidate solutions, which is the most time-consuming step), the adjustable approach is always slightly faster. In relation to the quality of the feasible solution, the adjustable approach is always better than the static one.

C O N C L U S I O N
A stochastic production-inventory-routing problem with recourse is considered. A sample average approximation method is applied where the stochastic problem is solved for several samples and the best solution among the solutions obtained is selected.
Since the stochastic problem is very complex and problem instances can hardly be solved to optimality within reasonable amount of running time, two different approaches are proposed: a static approach that consists of using heuristics to generate a solution for the stochastic problem defined for each sample; and a new adjustable heuristic procedure where the generated solutions are used to fix part of the first-stage solution (corresponding to those variables which have common values in the most of the candidate solutions). The partial solution is completed by solving a new restricted problem defined for a larger sample. Several algorithms were tested for each approach. Computational results have shown that the static procedure performs well on small size instances using as heuristic the (time) truncated branch-and-bound algorithm. For medium size instances, the static procedure performs well if an iterated local search heuristic is used to improve the initial candidate solutions. However, such solutions can be improved by using the adjustable approach with a short computational time.
Future research would be to extend the adjustable approach introduced here to other complex two-stage stochastic problems, such as network design problems.