A hybrid heuristic for a stochastic production-inventory-routing problem

We consider a stochastic single item production-inventory-routing problem with a single producer and multiple clients. At the clients, demand is allowed to be backlogged incurring a penalty cost. Demands are considered uncertain. A recourse model is presented where the production and routing decisions are taken before the scenario is known, and the quantities to deliver to the clients and the inventory levels are adjustable to the scenario. Valid inequalities are introduced and a hybrid heuristic that combines ideas from the sample average approximation method and from relax-and-fix approaches is proposed. Preliminary tests based on randomly generated instances are reported showing that the hybrid heuristic performs better than the classical sample approximation algorithm for hard instances.


Introduction
We consider a single item stochastic production-inventory-routing problem (SPIRP) with a single supplier/producer and multiple retailers/clients. A vendor managed inventory approach is followed where the supplier monitors the inventory at the retailers and decides on the replenishment policy for each retailer. Inventory aspects are considered both at the supplier and at the retailers. Demand is allowed to be backlogged at the retailers and in that case a penalty cost is incurred. Backlog is not allowed at the supplier. Demands are considered uncertain, following a uniform distribution. A constant production capacity at the supplier is assumed. The decision maker has to decide the production and the distribution plans for a finite time horizon. The production plan consists in defining the production periods and the amount to produce in each one of those periods. The distribution plan defines the retailers that should be visited in each time period, the quantities to deliver to each visited retailer, and the corresponding route in each time period. For the distribution plan a single vehicle is considered. The goal is to minimize the production and routing cost plus the expected inventory and penalty for backlogged demand costs. We assume the production plan and the choice of which clients to visit in each time period (and consequently the routing) are first stage decisions, that is, such decisions must be taken before the scenario is revealed. The quantities to deliver to each client at each time period, and the inventory levels, can be adjusted to the scenario (known as second stage decisions). Such assumptions may hold for short planning horizons.
Complex problems combining production, inventory and routing decisions have been receiving a great attention in recent years. For a recent survey see [?]. Stochastic approaches for related inventory routing problems have been also considered, e.g. [?,?,?,?,?]. Among these [?,?,?] present heuristic approaches. A classical approach for handling stochastic problems is the sample average approximation (SAA) method, see [?]. In this method the expected inventory and penalty for backlogged demand costs are approximated by a sample average estimate obtained from a random sample. The resulting sample average approximating problem (SAAP) is solved for different samples in order to obtain a set of candidate solutions. Then these candidate solutions are tested on a larger sample and the best solution for that larger sample is chosen.
However, most of the deterministic production-inventory-routing problems reported in the literature are not solved to optimality. Therefore, solving each SAAP to optimality, within reasonable running time, may be impractical for many instances. In [?] it is noticed that when the SAAP is not solved to optimality for each sample, then the use of a commercial solver based on a branch and cut algorithm with a running time limit acts as a heuristic which may lead to lack of stability, since one can obtain solutions whose objective function can significantly vary from sample to sample. Here we propose a hybrid heuristic (HH) that combines different approaches, namely the SAA method and relax-and-fix heuristics. First, following the SAA method, several small size samples are considered. For each sample, a relaxation is used to obtain a partial feasible solution. Then, the obtained partial solutions are explored in order to identify variables that are frequently fixed to zero or to one. Those variables are fixed and the restricted model is solved for a larger sample keeping all the remaining variables free. To improve the efficiency of both methods (SAA and HH), valid inequalities are added to the initial formulation.
To the best of our knowledge, the use of the information from several solutions obtained with the SAA method combined with relax-and-fix heuristics is new. The proposed heuristic has several advantages in relation to the classical SAA method. It does not require solving each sample set to optimality, it does not require the use of large sample sets, and it allows one to adjust the final solution to a larger sample set since many variables are kept free. Further, the number of variables that are fixed can be easily controlled.
The mathematical model is presented in Section 2 as well as the valid inequalities included. The SAA method and the proposed HH are introduced in Section 3. Computational tests that compare the SAA method and the HH and show the advantages of the HH are presented in Section 4. Final conclusions are given in Section 5.

Mathematical formulation and improvements
In this section we introduce the mathematical model for the SPIRP. 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, that is, decisions taken before the demands are known, while 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 routing cost plus the expected value of the inventory cost plus the penalty value 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, ob-tained by the Monte Carlo method. This larger set of s scenarios is regarded as a benchmark scenario set representing the true distribution [?].
Consider a graph G = (N, A) where N = {0, 1, . . . , 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} is the set of clients. T = {1, . . . , nt} is the set of periods.
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 limit at producer/client i, P is the production limit at each time period, Q it and Q it are the lower and upper limits to the delivery quantity in period t ∈ T at client i. For the objective function parameters define: S is the set up cost for producing in a period, P is the production cost of a unit of product, V is the vehicle usage cost, C ij is the travel cost from node i to node j, (i, j) ∈ A, H i is the holding cost of the product at node i ∈ N , B i is the backlog cost of the product at node i ∈ N .
Consider the following non-adjustable variables: binary variables y t indicate whether there is production in period t ∈ T , variables p t give the production level in period t ∈ T , binary variables z it indicate whether there is a visit to node i ∈ N c in period t ∈ T, the routing variables x ijt indicate whether a vehicle travels from node i to node j, (i, j) ∈ A, in period t ∈ T, v t is a binary variable that is one if the vehicle is used in period t ∈ T and zero otherwise; for (i, j) ∈ A, t ∈ T , f ijt is the artificial flow of a single commodity variable used to prevent cycles in the routing problem. Notice this is not the amount transported from node i to node j in period t ∈ T 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 adjustable variables q it (ξ) indicating the quantity delivered at node i ∈ N c in period t ∈ T ; s it (ξ) indicating the stock level of a 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 SPIRP is modeled by the following formulation.
For brevity we omit the explanation of the model. This model can be tightened through the inclusion of valid inequalities. From the underlying network flow structure of the SIRP one can derive several families of valid inequalities. Next we introduce three such families.
where d = i∈Nc t =1 d i (ξ) − i∈N I i + , r = d − ( d P − 1)P , and (x) + = max{x, 0}. Inequalities (19) are adapted from the uncapacitated lotsizing problem with backlogging [?], and are derived for each client and 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 t =1 z i = 0. Inequalities (20) ensure that the total demand from period k 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 k = 1, inequalities (20) can be written in the stronger format y . Inequalities (21) are mixed integer rounding inequalities (see [?]) derived from the constraint P t =1 y + i∈Nc r it (ξ) ≥ d which is a relaxation of the feasible set that imposes the total set up capacity plus the backlogged demand to cover the net demand of all clients until period t.
These inequalities can be further strengthened and/or extended. We present only the ones used in the computational tests, reported in Section 4, and omit further discussion. The formulation (1)-(18), (19), (21) is denoted by SPIRF.

Solution approaches
In this section we discuss two approaches to solve the SPIRF that will be tested in Section 4. A first approach is a classical SAA method described in [?]. In the SAA method m separate sample sets Ω k , k ∈ M = {1, . . . , m}, are considered, each one containing s scenarios. The model is solved for each one of the scenarios set, Ω k , (replacing Ω by Ω k in model SPIRF) with a time limit of α seconds, giving m candidate solutions. 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 by solving a pure linear programming problem on the recourse variables. The solution with the minimum average cost is chosen.
Next we describe the hybrid heuristic HH in four stages. First stage: obtain m solutions. Generate m samples of small dimension 1 . For each sample a partial solution is obtained by solving the SPIRF assuming the routing variables x ijt are continuous (relaxing constraints (15)) with a time limit of α 1 seconds. Variables z it are kept binary. Second stage: fix a set of z it variables. First assign a different weight w k to each one of the m partial solutions obtained in the first stage. Those weights are computed according to the objective function value, denoted by c k , computed on a larger sample of dimension 2 ≥ 1 . Those weights are normal-ized as follows: letc = max k∈M {c k } andd = min k∈M {c − c k |c k =c}. Then for each solution k define its weight w k =c +d−c k i∈M (c+d−c i ) . So m k=1 w k = 1. Let z k it be the value of variable z it in solution k. Then, for each pair (i, t), i ∈ N, t ∈ T, compute the weighted average of the values of z it obtained in the first stage: The values W (i, t) are between 0 and 1. Define two threshold values γ 1 and γ 2 and set z it = 0 if W (i, t) ≤ γ 1 , and z it = 1 if W (i, t) ≥ γ 2 . In order to define the percentage of variables to be fixed a priori, values γ 1 and γ 2 can be obtained from the empirical distribution of the W (i, t) values. Third stage: solve the simplified model for a larger sample. With the variables z it fixed in the previous stage and without the integrality constraints (15), solve the resulting model SPIRF for a new and larger sample with 3 ≥ 2 scenarios. Then, with all the variables z it fixed to their optimal value, the routing variables are determined by solving a TSP problem for each time period. Fourth stage: compute objective function value. Fix the first stage decision variables. Then, compute the objective function value for the very large sample, by computing the value of the recourse variables for each one of the s scenarios. This can be done by solving a pure linear programming model again.

Computational Results
In this section we report preliminary computational experience performed to assess the quality of the proposed HH, and compare it with the SAA method. All tests were run using a computer with an Intel Core i7-4750HQ 2.00GHz processor and 8GB of RAM, and were conducted using the Xpress-Optimizer 28.01.04 solver with the default options.
We present computational results for 7 sets of instances to the SPIRP, A 20 ,. . . ,A 80 . For each set A n with n clients, 5 instances are generated from different samples of demand values. 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 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 number of periods nt considered is 5. For each client and each period, the nominal demand value d it is randomly generate between 40 and 80 units, and we consider uncertain demands varying in [0.7d it , 1.3d it ]. The initial stock at producer I 0 is zero, 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 level S i is 500 for all i ∈ N . The production capacity P is 50% of the average demand. For all i ∈ N , the holding cost H i is 2 and the backlog cost B i is 3. The production set up cost, the unit production cost and the vehicle usage cost are given by S = 100, P = 1 and V = 50, respectively.
Preliminary tests (not reported here) have shown that inequalities (19)-(21) reduce the integrality gap by more than forty percent, on average. In general, when a given inequality (19) -(21), derived for a scenario ξ, cuts off the linear relaxation solution, then the same happens for most of the inequalities derived for the remaining scenarios in Ω. Hence, the number of cuts can be quite large when the number of scenarios is large. For the SPIRF we opted to omit constraints (20) as the number of inequalities added was too large for some instances.
The obtained computational results are displayed in Table 1. For each problem instance s = 1000 scenarios were generated. The number of candidate solutions used is m = 10.
When using the SAA method, each candidate solution is obtained by using a sample of = 5 and 25 scenarios and a time limit of α = 300 seconds. However, when no integer solution is found within α = 300 seconds the model is solved until the first integer solution is obtained. For each strategy, the obtained results, average computational time in seconds (in columns "T ") and average solution cost (in columns "C"), are displayed in columns SAA 1 and SAA 2 in Table 1. Columns SAA 1 and SAA 2 display the results for = 25 and = 5, respectively.
For the HH, in the first stage, each candidate solution is obtained by using a sample of 1 = 5 scenarios and a time limit of α 1 = 300 seconds. In the second stage a sample of 2 = 15 scenarios is used. In the third stage, when solving the simplified model, a sample of 3 = 50 scenarios is used. In the second stage, several strategies for fixing variables z it were compared corresponding to different values of parameters γ 1 and γ 2 . The results for two of these strategies are reported in Table 1. Columns HH 1 and HH 2 display the results for γ 1 = 0.05 and γ 1 = 0.15, respectively. The value of γ 2 is set to 1 − γ 1 . Now we compare strategies SAA 1 and SAA 2 . For instances with a small number of clients (n = 20, 30), the problems used to obtain the candidate solutions can be solved to optimality. Therefore, using the larger sample sets with 25 scenarios (thus more representative samples), better solution values are obtained. For instances with a higher number of clients the problems to obtain the candidate solutions could not be solved to optimality. In these