Christian Artigues
email: artigues@laas.fr

Jean-Charles Billaut
email: jean-charles.billaut@univ-tours.fr

Online recoverable robustness based on groups of permutable jobs for integrated production scheduling and delivery routing Azeddine Cheref *

Keywords: Integrated production and distribution scheduling, Recoverable robustness, Groups of permutable jobs

It is well known that the magnitude and the impact on schedule quality of uncertainty in problem data both increase when the integration of scheduling decisions takes place. This paper aims at comparing the standard robust optimization and the recoverable robust optimization frameworks to tackle scenario-based uncertainty in integrated production and distribution scheduling. While the standard robust optimization aims at providing an integrated baseline schedule that optimizes a performance criterion for the worst scenario, the recoverable robustness framework adds a repairing/recovering component, which is able to react to the realized scenario by modifying in a limited way the issued baseline schedule. Hence, the recoverable robustness framework has two merits: it is more realistic, as such repairing elements are often present in practice, and it is able in theory to reach better worst-case performances on the scenario set. However, the question remains whether under limited CPU time, recoverable robustness algorithms are able to actually reach better performance than robust algorithms for reasonable problem instance size. This paper proposes new MILP models and heuristics for a recoverable robust optimization framework based on the concept of groups of permutable jobs and carries out such a comparison on a one-machine production facility coupled with a single vehicle distribution system. For this problem, the experimental results exhibit conditions under which a recoverable robust approach is likely to yield better results than the standard robust approach.

Introduction

In a supply chain, two essential elements are the production and the distribution. These two problems used to be solved separately, but it it well known that considering them jointly may lead to a global optimization of the supply chain performances. In the literature, more and more research deal with integrated approaches, involving production and distribution decisions, at strategic, tactical and operational levels. The relation between production and distribution may involve an intermediate storage area and, in this case, inventory levels are considered as part of the production and distribution problems. In the more general terms and depending on the decision level, the production problem is a planning problem and/or a scheduling problem, where a manufacturer has to decide in which location the orders will be processed, at which time period and in which quantity. Then, the delivery problem is an assignment and/or a routing problem, where the manufacturer has to decide how many vehicles will be used at each location, how orders are grouped into delivery batches and when they will be delivered. We refer to [Chen, 2004, Chen, 2010a, Wang et al., 2014] for a global overview of integrated production and distribution problems.

In this paper, we incorporate the delivery plan of a vehicle into a single machine scheduling problem, representing a single manufacturing facility. We consider the simplified case where there is only one vehicle available to serve the customers, with infinite capacity. However, we assume that the data are known with uncertainty. The objective is to find a schedule and a delivery plan so that a robustness criterion is minimized, under a scenario-based uncertainty modeling. In contrast with standard robust optimization approaches for scheduling [Kouvelis and Yu, 1997], we do not propose a single complete solution to the problem that has to be feasible w.r.t. all scenarios and minimizes a worst-case criterion over the scenarios. Instead, we adopt the recoverable robustness framework that considers first-stage decisions and second-stage recovery options, as in stochastic programming except that no probabilities are attached to scenarios and that the possibilities for recovery are limited by a predefined set of recovery algorithms [Liebchen et al., 2009]. Following this scheme, we propose for the first-stage a set of solutions by using the concept of groups of permutable jobs presented in [Billaut andRoubellat, 1996, Artigues et al., 2005]. During second stage, a recovery algorithm uses the information from the revealed scenario to compute a complete schedule from the set of solutions issued from the first-stage. However, we consider that the scenario is unrevealed progressively during schedule execution, so the recovery algorithm has to be both greedy and online, exploiting only known information about the jobs available to be scheduled or to be delivered at decision time. This appear to be a practical characteristic of many applications both on production scheduling (where for example the exact release date of a job is known only when the job arrives) and on vehicle routing (where the traveling time from a location to another is only accurate in a short term basis due to traffic fluctuation). We will refer to this framework as online recoverable robustness.

Despite its practical appeal, the online recoverable robustness framework raises two major issues, one linked to modeling difficulties and another one linked to tractability in solution methods. For the first issue, a basic question is to know whether is it possible model a non trivial online recoverable robust integrated scheduling and vehicle routing problem such as the one considered in this paper in a standard optimization paradigm such as the mixed-integer linear programming framework. In this paper, we answer positively to this question by providing new MILP models that extends the one proposed in [Kouvelis and Yu, 1997] for the standard robust one machine scheduling problem. The second issue is to design, under limited CPU time, recoverable robustness methods that are able to actually reach better performance than standard robust methods for reasonable problem instance size. In this work, a comparison of a robust heuristic approach with an online recoverable robust heuristic is carried out on a one-machine production facility coupled with a single vehicle distribution system. The experimental results allow to determine conditions under which the proposed recoverable robust approach is likely to yield better results than the standard robust approach.

The paper is organized as follows. Section 2 presents a literature review on integrated production scheduling and vehicle routing problems on the one hand and robust scheduling and vehicle routing on the other hand. In Section 3, we give a formal description of the problem, introducing the proposed online recoverable robustness scheme. In Section 4, we provide a (single stage) robust mixed-integer linear programming formulation of the problem. Section 5 presents a two-stage robust mixed-integer linear programming formulation, instantiating the concept of online recoverable robustness and making use of the group-of-permutable-jobs structure. As both the single-stage and two-stage robust MILP quickly become intractable as instance size grows, we present in Section 6 a heuristic based on tabu search to build a robust solution on realistic instances. In Section 7, we provide computational experiments to validate our proposal. Concluding remarks and directions for further research are pointed out in Section 8.

Related work

To our knowledge, no approach for integrated production scheduling and vehicle routing under uncertainty can be found in the literature. Hence we divide this state-of-the-art section into deterministic integrated production and distribution scheduling (Sect. 2.1), robust scheduling (Sect. 2.2) and robust vehicle routing (Sect. 2.3).

Deterministic integrated production and distribution scheduling

Few papers deal with an integrated model of production and distribution at an operational level, i.e. with vehicle routing as a part of the decision. In [Chang and Lee, 2004], the authors consider a single machine or a two-parallel machines environment, a single vehicle and one or two customer areas. Each job has a size and the capacity of the vehicle is limited. The authors give NP-hardness proofs and algorithms with worst-case analysis. The same problem is considered by [He et al., 2006] and [Zhong et al., 2007], in which some improvement are given. In [Lu et al., 2008], the authors consider a single machine problem with release dates and a single vehicle to deliver the products to one customer area. The vehicle has a bounded capacity and the objective is to minimize the maximum delivery time of the jobs. The authors prove the NP-hardness of the problem and propose an algorithm with worst-case analysis. They consider also the case when preemption is allowed and propose a polynomial time algorithm for this problem. In [Condotta et al., 2013], the authors treat a similar problem where due dates for delivery are integrated. Multiple vehicles are available for the delivery and the objective is to minimize the maximum lateness. The authors propose a mixed linear program, lower bounds and a tabu search algorithm for solving the problem. In order to minimize the maximum delivery time of the jobs, [Dong et al., 2013] consider the two-machines open shop problem and a single vehicle with bounded capacity to deliver the products to one customer area. The authors consider the special case when the vehicle can take only one job and a general capacity case. The authors propose a polynomial time algorithm with a worst case ratio for the problem. In [Li et al., 2005a], the authors consider a single machine and single vehicle problem. The vehicle has a bounded or unbounded capacity and the objective is to minimize the sum of delivery completion times. The authors prove the NP-hardness of the problem. They propose an efficient dynamic programming algorithm for the single-customer case and, for the general case, they propose a polynomial time dynamic programming algorithm when the number of customers is fixed. A similar problem is considered in [Levin and Penn, 2008], with an arbitrary number of customers and an uncapacitated delivery vehicle. The objective function is to minimize the total latency. The authors propose a ∼ =16.3-approximation algorithm. In [Li and Ou, 2005b], the authors consider a factory (single machine) and a warehouse with two different locations. The raw materials (jobs) are located at the warehouse before processing, have to be delivered to the factory for their processing and then delivered back to the warehouse when completed (there is no vehicle routing problem). Only one vehicle is available with a limited capacity (in terms of number of jobs) for transporting unprocessed jobs, and another, also of limited capacity, for transporting processed jobs. The authors prove the NP-completeness of the problem and propose a polynomial time dynamic programming algorithm when the sequence of processing jobs is fixed. A heuristic algorithm with performance guarantee is proposed for the general case. This model is extended in [START_REF] Wang | Production scheduling with supply and delivery considerations to minimize the makespan[END_REF] by considering one warehouse for material supply and another warehouse for finished products. There are two vehicles, one for the travels between the first warehouse and the factory and one for the travels between the factory and the second warehouse. The authors derive some optimal properties, study some polynomially solvable cases and propose heuristic algorithms with performance guarantee for the general case and some particular cases. In [Li and Vairaktarakis, 2007], the authors consider a machine processing environment called "bundling operations", where each job is composed by two tasks that are processed on dedicated machines. A job is ready when both tasks are completed and can then be delivered to the customer location. There is an unlimited number of vehicles with limited capacity. The objective is to schedule the jobs and to define the routes of the vehicles in order to minimize a delivery cost plus a weighted sum of the delivery completion times. The authors propose polynomial-time heuristic algorithms and dynamic programming-based approximation schemes. In [Scholz-Reiter et al., 2010], the authors consider a global supply chain context, and focus on the interface between manufacturing and logistic systems along the supply chain. Orders are produced in an hybrid flowshop (flowshop with assignment problem) and then assigned to tours. The objective function is a cost function, composed by a production cost, a storage cost (for intermediate and finished products), a tardiness cost, the cost of conducting a tour and the cost of a tour. The authors propose a mathematical programming formulation, tested on a small instance. In [Chen, 2010b], the author considers a problem with a single machine and multiple vehicles. To each job is associated with a processing time and a volume. All the jobs have to be delivered at only one customer location. Vehicles have different limited capacities and different routing times. The objective is to minimize the makespan. The author proves NP-completeness and propose a MILP formulation. For the problem with a single machine, an on-line version is considered by [Ng and Lu, 2011], in which the processing time of each job becomes known at its arrival time. The authors consider the case where preemption is allowed and algorithms with worst case analysis are proposed for both of them. In [Wang and Liu, 2013], the authors consider a problem with a single machine and a single vehicle with limited capacity. The customers are located at the nodes of a star-shaped network (disjoint branches with a common endpoint corresponding to the production center). Each job has a physical space and the vehicle capacity is bounded. The objective is to find a processing sequence and a delivery rule in order to minimize the makespan (maximum delivery completion time). The problem is NP-hard and the authors propose approximation algorithms in the identical-job-size case and in the general case. In [Ullrich, 2013], the author considers a parallel machine scheduling environment with machine ready times. The jobs have to be assigned to certain tours while the routes and the schedule of the tours have to be decided, so that the total tardiness is minimized. A fleet of several vehicles is available and a ready time is associated with each vehicle. Jobs have a size and the vehicles capacity is limited. The author proposes MILP formulations for the scheduling problem and the routing problem and for the integrated problem. Two decomposition heuristics and a genetic algorithm are proposed and compared on specific benchmark instances. Many work on integration of production and distribution scheduling concern perishable products. [Armstrong et al., 2008] consider a zero-inventory production and distribution problem with a single transporter and a fixed sequence of customers. They propose a branch-and-bound procedure for the problem. [Geismar et al., 2008] develop lower bounds and heuristics for an integrated production and transportation scheduling problem with capacity constraints and also no inventory, which means that a completed production job must be transported immediately to the customer. In [Gao et al., 2015], the authors propose a heuristic for the problem of minimizing the makespan on a single machine scheduling problem with a unique capacitated vehicle and a no wait constraint. In [START_REF] Lee | [END_REF], the authors consider the problem of nuclear medicine production and delivery. Each job has to be assigned to a production machine and has to be delivered in a given time window. Several vehicles with different capacities are available. The objective is to minimize a cost function composed of a production cost, a fixed cost for the use of a vehicle and a travel cost. The authors propose an MILP formulation and a large neighborhood search heuristic, tested on benchmark instances issued from the VRP literature. [Viergutz and Knust, 2014] pointed out and corrected an error in the branch-and-bound method of [Armstrong et al., 2008]. They also extend the models to variable production and delivery sequences and propose several heuristics.

Robust scheduling

"Uncertainty affects a wide range of decisions managers, engineers, and other decision makers have to make" [Kouvelis and Yu, 1997]. This is specially true for production managers, where data such as processing times or release dates are generally not known with an absolute certitude, and for delivery problems, where transportation durations are time dependent and subject to fluctuations due to traffic. In this context, it is generally important to find a solution to the problem which is not only a good solution, but also a solution that keeps a good performance in presence of uncertainties. This is the field of robust approaches. There are several ways to understand and to structure uncertainty and several ways to define a robust approach [Billaut et al., 2008]. The aim of such an approach is to produce a solution that will have a "reasonable objective value" under any input data scenario [Kouvelis and Yu, 1997].

Robustness considerations in scheduling are frequently treated in the literature and it is not possible to make here an exhaustive review of the state-of-the-art. We voluntarily restrain our work to robust scheduling approaches where no probability distribution is available for the uncertain data. Instead there is a set of possible scenarios for the input data. Under this framework, based on uncertainty scenarios, the concept of robust scheduling appeared in [Daniels and Kouvelis, 1995] (in 1992 for a working paper version). The authors consider a single machine environment, and processing time uncertainty is represented by interval data, which yields an infinite number of scenarios. The authors consider the total flow time as the deterministic objective function and the robustness is based on absolute or relative deviations from an optimal performance. Since this seminal paper, a lot of papers have addressed uncertainty in scheduling. Two survey papers deal with the generation of proactive (robust) schedules for the resource constrained project scheduling problem: [Herroelen andLeus, 2004, Herroelen andLeus, 2005]. Aytug et al. [Aytug et al., 2005] review the literature on production scheduling in the presence of unexpected events/unforeseen disruptions and propose a taxonomy for uncertainty in production manufacturing. Another recent survey on general approaches for scheduling under uncertainty can be found in [Verderame et al., 2010].

In this paper we integrate three families of approaches for scheduling under uncertainty. The first family is based on implicit consideration of uncertainty via the concept of computing a set of solutions rather than a single solution for a scheduling problem in which data is subject to unmodeled disruptions. In production scheduling, a family of methods described notably in [Billaut andRoubellat, 1996, Artigues et al., 2005] aimed at computing on each machine a sequence of groups of permutable jobs such that any permutation of the job inside each group yields a feasible job sequence with a guaranteed value on the objective function as soon as the job durations do not exceed a predefined upper bound. In other words, a sequence of groups of permutable jobs is a structured partial solution such that each complete solution issued from it has a guaranteed performance. The drawback of these approaches is that there is no indication on how to select the job order inside each group in response to a particular disruption. However, in [Wu et al., 1999], the author have shown through simulation of various processing time disturbances that computing a static way such a solution structure 1 and then allowing the remaining decisions to be taken dynamically via dispatching rules yields superior results.

The second family of approaches is based on robust discrete optimization for scheduling as stated in [Kouvelis and Yu, 1997]. In this framework (see also several other surveys or seminal works on robust discrete optimization [START_REF] Bertsimas | Robust discrete optimization and network flows[END_REF], Bertsimas and Sim, 2003, Aissi et al. 2009, Gabrel et al. 2014]), a complete sequence of jobs on each machine is the output of the robust scheduling method. The uncertainty is represented by a continuous or discrete set of scenarios, each scenario corresponding to a deterministic value of the problem parameters. On a particular scenario, the performance of the sequence is generally measured by computing the earliest start schedule compatible with the prescribed sequence and the realized scenario. Hence, at least for a finite scenario set, the worst case performance of a sequence on the scenario set can be computed a priori. Note that most robust scheduling problems (aiming at finding the optimal sequence, having the best worst-case performance) are NP-hard even for single-machine scheduling problems [Aloulou and Della Croce, 2008].

The third family of approaches is based on recoverable robustness. This framework, defined in [Liebchen et al., 2009], aims at reducing the conservativeness of robust optimization as described in [START_REF] Bertsimas | Robust discrete optimization and network flows[END_REF]). Recoverable robustness considers two sets of decision variables and consists in providing an assignment for the first-stage decision variables and a family of recovery algorithms such that, for each scenario, there exists one recovery algorithm in the family that computes feasible assignments for the second-stage decision variables. Thanks to this algorithm, the baseline solution adapts to the realized scenario. Note that we can consider that the classical robust scheduling framework (second family of approaches) falls into this scheme, with a simplified recovery algorithm. Indeed, its outcome is made of both the prescribed sequence (which can be defined as an assignment of first-stage decision variables) and the earliest scheduling algorithm (which can be defined as the prescribed algorithm for assigning the start times considered as second-stage decision variables). Recently, recoverable robustness has been successfully applied to scheduling problems [START_REF] Caprara | [END_REF].

Robust vehicle routing

Robustness considerations in vehicle routing have been addressed since 1990s. Bertsimas and Simchi-Levi in [Bertsimas and Simchi-Levi, 1996] present a survey on papers dealing with uncertainty issues in vehicle routing. Most of the problems consider that the demand is stochastic (demands arrive randomly, demands have a random size, etc.), and the authors propose robust algorithms, i.e. algorithms "independent of the specific environment and the variability of the data". The same year, in [Gendreau et al., 1996], the authors propose a summary of the literature on stochastic vehicle routing problems. In their survey, stochastic demands, stochastic travel times and stochastic customers are considered. Besides this literature on vehicle routing based on stochastic programming approaches, robust optimization methods have been more recently addressed. A nice state-of-the-art review and a robust optimization approach can be found in [START_REF] Agra | The robust vehicle routing problem with time windows A[END_REF] for the vehicle routing problem with time windows. In this paper, the authors consider an uncertainty polytope for the travel times. Along with standard robust optimization, they also consider the adjustable robustness framework [Ben-Tal et al., 2004], which is close to the recoverable robustness framework, except that the second-stage variables can be adjusted without any algorithmic restriction to the realized scenario. Notably, they allow arrival times to be adjustable to the realized scenario. Other authors considered robust optimization for the capacitated vehicle routing with uncertain demand [Gounaris et al., 2013].

permutable jobs" structure defined in [Billaut and Roubellat, 1996] 3 Problem description and online recoverable robustness Following the classical three-field notation of scheduling problems [Graham et al., 1979], Z-L. Chen introduces in [Chen, 2010a] a five-field notation α|β|π|δ|γ for the problems where scheduling and routing are integrated at an operational level. In this notation, α, β and γ correspond to the machine configuration, scheduling constraints and objective function respectively, similarly to the notation of scheduling problems. Fields π and δ specify the characteristics of the delivery process and the number of customers, respectively. In the following, we focus on a single machine environment (α = 1), with routing as a delivery problem (i.e. solving a vehicle routing problem is part of the decision), we have a single vehicle with unbounded capacity and n customers to deliver

(π = {routing, V (1, ∞)}, δ = n).
We consider a set J = {J 1 , ..., J n } of n jobs to schedule on a single machine representing a manufacturing facility (plant). We denote by p j the processing time of J j , r j the release date of J j that corresponds to the earliest possible start time of production and d j the delivery due date of J j , i.e. the date at which J j is supposed to be delivered to the customer. It is assumed that the storage capacity of jobs after production is not limited. However, we also consider that the jobs have strict due dates, which is equivalent to say that the maximum lateness is bounded. Each job J j is associated with a customer location j, and 0 denotes the location of the plant. t 0,j denotes the travel time between the plant and the customer location j and t i,j is the travel time between customer locations i and j (∀i, j, 1 ≤ i, j, ≤ n). Decision variable C j denotes the production completion time of J j at the plant (∀j, 1 ≤ j ≤ n). Decision variable D j denotes the delivery completion time of J j , i.e. the date of delivery at the customer location (∀j, 1 ≤ j ≤ n). Remember that jobs are delivered by a single vehicle that may return to the plant several times (as in the multi-trip traveling salesman problem) with unlimited capacity (capacity K = n). Without uncertainty considerations, the problem can be denoted by 1|r j |routing, V (1, ∞)|n|γ where γ is the objective function. In the following, we have γ ∈ {D max , L max } with D max the maximum delivery completion time, defined by D max = max 1≤j≤n D j and L max the maximum delivery lateness, defined by L max = max 1≤j≤n (D j -d j ). The problem is to find a production sequence of jobs (equivalently, to determine jobs production completion times), a batching of jobs and a route for each batch, so that the objective function is minimized. A batch can only start when all the jobs of the batch are completed.

We assume in addition that uncertainty affects the input data and that our objective is to find a robust solution. To represent uncertainty, we use a scenario based approach similar to the one introduced in [Kouvelis and Yu, 1997]. The data associated to the jobs and to the transportation times between sites vary according to predefined scenarios. We denote by S the set of possible scenarios over the planning horizon. Now, we say that to each job J j is associated a release date r s j , a processing time p s j and a due date of delivery d s j for each scenario s ∈ S. The matrix of travel times T s = (t s i,j ) 0≤i≤n,0≤j≤n gives the transportation time from any site i to any site j under scenario s ∈ S.

We denote by C s j the production completion time of job J j under scenario s and by D s j its delivery completion time under scenario s. The maximum lateness of J j under scenario s is denoted by L s j = D s j -d s j . Under uncertainty, the following objective functions are defined:

the delivery makespan D max = max s∈S (max 1≤j≤n D s j ), the maximum lateness L max = max s∈S (max 1≤j≤n L s j ).

We now formally define the concept of online recoverable robustness. Let us consider the standard robust optimization problem that we can denote as min

x∈∪ s∈S Xs max s∈S f (x, s)
where S is the set of scenarios, X s is the set of solutions feasible for scenario s and f (x, s) is some performance indicator of solution x ∈ X s under scenario s ∈ S. This definition implies that x has to be feasible for any scenario, inducing a high degree of conservativeness [START_REF] Bertsimas | Robust discrete optimization and network flows[END_REF].

The concept of recoverable robustness [Liebchen et al., 2009] is often considered in the literature as a two-stage extension to the robust optimization problem, in which the (first-stage) solution x can be modified in a limited way in the second-stage to fit with the realized scenario to obtain solution y = A(x, s) where A(x, s) is the application-dependent recovery algorithm of solution x under scenario s. Hence the recoverable robust optimization problem becomes:

min x∈X max s∈S f (A(x, s), s)
where X is the scenario-independent set of first-stage solutions, which is a more realistic definition in practical applications. Note that x does not necessarily assign all the decision variables of the problem. As in stochastic programming x can be only the subset of first-stage decision variables, while algorithms A(x, s) returns a complete solution, giving values to both first-stage and second stage decision variables.

Instantiating this concept to our integrated scheduling and vehicle routing problem, we denote by π a complete or partial sequence of the jobs on the machines, B a batching (i.e. a partition of the jobs) and Σ = {σ(B), B ∈ B} the set of complete or partial routing σ(B), i.e. the complete or partial job delivery sequence inside each batch B ∈ B). Let C s (D s ) denote the vector of second-stage job completion times (delivery times, respectively). We assume that given a first stage solution x = (π, B, σ) the recovery algorithm A(x, s) output a feasible solution y = A(x, s) = (C s , D s ). Function f (y, s) is either D max or L max as defined above. However, as both the scheduling and routing problems are time-constrained problems, we assume that algorithm A(x, s) must be an online algorithm, i.e. called during schedule execution. At a given decision time t the algorithm must take quickly a basic decision (i.e. what is the next job to be scheduled at time t or what is the next job to be delivered at time t) with the only knowledge on the scenario that has been revealed up to time t (i.e. what are the jobs available for scheduling at time t and what are the estimated travel times at time t). Sections 4 and 5 presents two different definitions of the first-stage decision variables x = (π, B, σ), the first one, inspired by [Kouvelis and Yu, 1997], defining complete sequences of jobs both in production and delivery and the second one defining partial sequences based on the concept of groups of permutable jobs [Billaut andRoubellat, 1996, Artigues et al., 2005] under the recoverable robustness framework. Associated with these two different first-stage decisions, two different online scheduling algorithms are defined to output the two-stage decision variables y = (C s , D s ).

In Sections 4 and 5, we extend the robust scheduling mixed-integer linear programming models proposed by [Kouvelis and Yu, 1997] to the fixed-sequence robust scheduling problem and to the online recoverable robust scheduling problem, respectively. Section 6 considers a particular case of the second model to devise a tabu search heuristic able to deal with large-size instances.

The fixed sequence robust model

This Section presents mixed integer linear programs (MILPs) for the robust integrated scheduling and vehicle routing problem using the concept proposed by [Kouvelis and Yu, 1997]. At first stage a complete production sequence π, a set of batch B and, for each batch B ∈ B, a complete routing σ(B) are computed. The online recovery algorithm A(x, s) aims at strictly following the prescribed sequences and routings, assigning the production and delivery completion times as early as possible according to the prescribed sequence on the realized scenario. In section 4.1, we first consider the particular case in which the delivery part of the problem is ignored. Three MILPs are proposed to model the optimal assignment of the first-stage variables and second-stage variables according to the on-line algorithm. In Section 4.2, a MILP is proposed for the integrated scheduling and vehicle routing problem. In both cases, an illustrative example is given.

MILP for the fixed sequence robust scheduling problem

In this section, we consider the scheduling problem with release dates and maximum lateness minimization under several scenarios based on fixed sequences and earliest start online recovery algorithm. Note that the problem is NP-hard with only one scenario. The objective is thus to find a single machine sequence π that is the best possible for all scenarios, given that the online algorithm will adjust at the earliest the completion times C s i to the realized scenario. The following model (SM) is based on the one proposed by [Kouvelis and Yu, 1997] for the 1|| C j problem. We define positional variables x j,k equal to 1 if job J j is in position k in sequence π and 0 otherwise.

(SM) min L max s. t. n j=1 x j,k = 1, ∀k ∈ {1, ..., n} (1)
n k=1 x j,k = 1, ∀j ∈ {1, ..., n} (2)
L max ≥ n j=1 r s j x j,k + k q=k n j=1 p s j x j,q - n j=1 d s j x j,k , ∀k, k ∈ {1, ..., n}, k ≥ k , ∀s ∈ S (3) x j,k ∈ {0, 1}, ∀j, k ∈ {1, ..., n} (4)
L max ≥ 0, (5)
Constraints ( 1) and ( 2) guarantee the feasibility of the sequence, assigning exactly one job at each position. Constraints (3) express the fact that the worst-case maximum lateness (L max ) has to be larger than the maximum lateness under scenario s, given by the right hand side. Indeed n j=1 r s j x j,k gives the release date of the job in position k under scenario s, and k q=k n j=1 p s j x j,q is equal to the total duration of the jobs between positions k and k under scenario s, while n j=1 d s j x j,k gives the due date of the job at position k in scenario s. These constraints are based on the fact that, according to the earliest time online recovery algorithm, there exists a block of activities starting by a job scheduled at its release date, followed by jobs scheduled without idle times and ending by the job that fixes the L max . Hence the explicit usage of the second-stage completion time variables is not even necessary. This model contains n 2 binary variables, 1 continuous variable and O(n 2 |S|) constraints.

Example

Let us consider an instance with n = 5 jobs and the following data.

s = 1 J 1 J 2 J 3 J 4 J 5 r 1 j 0 7 3 4 3 p 1 j 3 4 1 2 4 d 1 j 3 14 4 6 10 s = 2 J 1 J 2 J 3 J 4 J 5 r 2 j 3 3 0 1 7 p 2 j 2 5 1 3 3 d 2 j 6 11 2 5 14
The optimal robust sequence for this example is (J 3 , J 1 , J 4 , J 5 , J 2 ) with a maximum lateness L max = 5. Fig. 1 provides the second-stage completion time variables and maximum lateness of this fixed sequence for each scenario.

-0 2 4 6 8 10 12 14 16 18 20

J 1 J 2 J 3 J 4 J 5 Scenario s = 2 L max = max{5 -6, 16 -11, 1 -2, 8 -5, 11 -14} = 5 - 0 2 4 6 8 10 12 14 16 18 20 J 1 J 2 J 3 J 4 J 5 L max = max{7 -3, 17 -14, 4 -4, 9 -6, 13 -10} = 4 Scenario s = 1
Figure 1: Sequence (J 3 , J 1 , J 4 , J 5 , J 2 ) under the two possible scenarios

MILP for the fixed sequence robust integrated scheduling and vehicle routing problem

In this section, we a present a new MILP model in which the robust scheduling and the vehicle routing problems are integrated. For these models, the scheduling part is based on the previous model. Here the first stage variables are the complete production sequence π, the batching B and the routings σ(B) for each batch B ∈ B.

For the following model, the binary variables associated to the vehicle routing part of the problem are the following. Let B r denote the r th batch. r varies from 1 to n and some batches can be empty. A binary first-stage variable z i,j,r is equal to 1 if tour σ(B r ) visits the site of job i and then the site of job j consecutively, and 0 otherwise, while another binary first-stage variable y j,r is equal to 1 if tour σ(B r ) visits site j (in other words if j belongs to batch B r ). We denote by D s j,r ≥ 0 the second-stage delivery completion time of job J j (D s 0,r is the start time of tour σ(B r )) if it is in tour σ(B r ), for scenario s. The dummy site (n + 1) is introduced as a finishing site of a tour (D s n+1,r is the finishing time of tour σ(B r ) in scenario s). The model called SRM comes from model SM but introduces the second stage completion times of the kth job in scenaro s ( Cs k ) with constraints (6) that is still based on the concept of blocks, and uses constraints ( 7)-( 13) to set the worst-case L max relatively to the worst case delivery times.

Note that variable Cs

k is necessary to serve as an earliest start time in scenario s of any tour that delivers the site of a job scheduled at position k.

(SRM) min L max s. t. n j=1 r s j x j,k + k q=k n j=1 p s j x j,q ≤ Cs k , ∀k, k ∈ {1, ..., n}, k ≥ k , ∀s ∈ S (6) n j=0,j =i z i,j,r = y i,r , ∀i ∈ {0, ..., n}, ∀r ∈ {1, ..., n} (7)
n i=0,i =j z i,j,r = y j,r , ∀j ∈ {0, ..., n}, ∀r ∈ {1, ..., n} (8)
n r=1 y j,r = 1, ∀j ∈ {1, ..., n} (9)
D s 0,r ≥ D s n+1,r-1 , ∀r ∈ {2, ..., n}, ∀s ∈ S (10) D s j,r ≥ D s i,r + t s i,j -M (1 -z i,j,r ), ∀i ∈ {0, ..., n}, ∀j, r ∈ {1, ..., n}, i = j, ∀s ∈ S (11) D s 0,r ≥ Cs k -M (2 -x j,k -y j,r ), ∀j, k, r ∈ {1, ..., n}, ∀s ∈ S (12) L max ≥ D s j,r -d s j , ∀j, r ∈ {1, ..., n}, ∀s ∈ S (13) z i,j,r ∈ {0, 1} ∀i, j ∈ {0, ..., n}, i = j, ∀r ∈ {1, ..., n} (14)
y i,r ∈ {0, 1} ∀i, j ∈ {0, ..., n}, ∀r ∈ {1, ..., n} (15)
D s i,r ≥ 0 ∀i, j ∈ {0, ..., n}, i = j, ∀r ∈ {1, ..., n} (16)
C s k ≥ 0 ∀k ∈ {1, ..., n}, ∀s ∈ S (17)
and Constraints (1), ( 2), ( 4), ( 5) Constraints ( 7) and ( 8) indicate that if a site i belongs to a tour σ(B r ), then there exists a unique ingoing arc and a unique outgoing arc selected for i in tour σ(B r ). Constraints (9) ensure that a site j can belong only to one tour. Constraints (10)-( 13) are the scenario-dependent constraints. Constraints (10) ensure that the tour σ(B r ) starts after the end of previous tour σ(B r-1 ) in scenario s. Constraints (11) sets for each scenario s the delivery completion time of job J j as the completion time of the delivery of a job J i plus the traveling time from site i to site j, provided that arc (i, j) has been selected in a tour σ(B r ). These constraints also automatically eliminate subtours. Constraints (12) ensure that tour σ(B r ) starts after the execution of all the jobs which are delivered by σ(B r ). Constraints (13) gives the expression of the maximum lateness for each scenario s. This model contains O(n 3 ) binary variables, O(n 2 |S|) continuous variables and O(n 3 |S|) constraints.

Example

For this example, we change the delivery due dates d j and data routing are added to the previous example. Release dates and processing times are similar to those of the previous example.

s = 1 J 1 J 2 J 3 J 4 J 5 r 1 j 0 7 3 4 3 p 1 j 3 4 1 2 4 d 1 j 11 18 9 10 17 s = 2 J 1 J 2 J 3 J 4 J 5 r 2 j 3 3 0 1 7 p 2 j 2 5 1 3 3 d 2 j 9 17 10 11 18
The transportation time matrices (t 1 i,j ) and (t 2 i,j ) are the following:

(t 1 i,j ) =         0 2 4 3 2 2 2 0 3 2 1 3 3 3 0 3 2 2 3 2 3 0 1 3 2 1 2 1 0 2 2 3 2 3 2 0         (t 2 i,j ) =         0 3 2 2 2 3 3 0 4 1 2 3 2 4 0 4 3 2 2 1 4 0 2 4 2 2 3 2 0 4 3 3 2 4 4 0        
Let us consider the schedule (J 3 , J 4 , J 1 , J 5 , J 2 ) for the production phase and then two batches composed by jobs B 1 = {J 1 , J 3 , J 4 } and B 2 = {J 2 , J 5 } with routings (J 4 → J 3 → J 1 ) and (J 5 → J 2 ), respectively. The production and transportation phases are illustrated in Fig. 2 for each scenario. In the figure, the scheduling part is displayed with a standard Gantt chart while the routing part is symbolized via arrows whose length gives the transportation time. Visited sites are displayed in parenthesis, where (j) denotes the site of job J j for j > 1 and (0) denotes the depot. The start time of each batch is marked with a double bar (||). While the visit of each site is marked with a simple bar (|). In scenario 1, site J 3 is delivered at time 12 while the deadline of J 3 for scenario 1 is d 1 3 = 9. In scenario 2, site J 2 is delivered at time 20 and d 2 2 = 20. For both scenario we have

L max = 3. - (4) - (0) B 1 (1) - - (0) (3) (2) B 2 -- (5)
L max = 3 - Scenario s = 2 J 1 J 2 J 3 J 4 J 5 0 2
(2)

B 2 (4) - (0) B 1 (1)
--

(0)(0) -- (5)
- In this section we present a more flexible online recoverable robustness model for the integrated scheduling and vehicle routing problem. We now suppose that the first-stage solution x = (π, B, Σ) is made of a set of batches B that partitions the jobs (just as in the previous model) but the production and routing sequences π and σ(B) for each batch B ∈ B are now partial sequences. Furthermore, we consider only partial sequences that can be expressed as a sequence of groups of permutable jobs. Hence a partial sequence involving a set E of jobs is of the form g 1 , g 2 , . . . , g q where ∪ q h=1 g h = E and ∩ q h=1 g h = ∅ and represent the precedence relation J i ≺ J j , ∀i ∈ g x , ∀j ∈ g y , ∀x < y. Consequently π is a sequence of groups of permutable jobs on J , while σ(B) is a sequence of groups of permutable jobs on B for each batch B ∈ B. Hence the relative ordering of the jobs inside each group is left undetermined during the first-stage. It follows that the recovery algorithm A(x, s) has to determine the total order besides assigning the completion and delivery times. We will specify A(x, s) and provide an online recovery robust MILP for scheduling in Section 5.1. We will extend the approach to the integrated scheduling and routing problem in Section 5.2.

L max = 3 J 1 J 2 J 3 J 4 J 5 Scenario s = 1 - (0)

MILP for the online recoverable robust scheduling problem

In this subsection, we only consider production scheduling without delivery routing. Delivery routing will be considered in the next subsection. As already mentioned, we assume that the recovery algorithm is a greedy online algorithm. When a job J j completes at a given time t, the online recovery algorithm A(x, s) has only the knowledge on the operations available at time t for scenario s (i.e. such that r s i ≤ t) and no other information. If there are operations available in the group of J j , the algorithm selects one of them. Otherwise, the A(x, s) scans the available operations of the next group. To summarize, A(x, s) follows the sequence of groups and schedule the jobs inside a group following the earliest release date (ERD) rule according to the realized scenario.

To simulate this rule in a MILP formulation, let us introduce the list of predecessors and successors of each job for each scenario according to the ERD rule. We denote by Succ s j and by P rec s j the list of successors and of predecessors of J j under scenario s according to tule ERD, defined as follows.

P rec s j = {J i ∈ J /(r s i < r s j ) or (r s i = r s j and i < j)} Succ s j = {J i ∈ J /(r s i > r s j ) or (r s i = r s j and i > j)} These lists are used to ensure that the ERD rule is respected in each group for each scenario. To model the group sequence, we also define variables χ j,k equal to 1 if job J j is in group G J k , and 0 otherwise, 1 ≤ j ≤ n, 1 ≤ k ≤ n. Note that the number of groups cannot exceed n (There cannot be more than one group per job). The online recoverable robust MILP for the one machine scheduling problem (OSM) can now be proposed.

(OSM) min L max (18) s. t. n k=1 χ j,k = 1, ∀j ∈ {1, ..., n} (19)
r s i + p s i + l∈(Succ s i ∩P rec s j ) p s l χ l,k + p s j -M (2 -χ i,k -χ j,k ) ≤ C s j , ∀i, j, k ∈ {1, ..., n}, j ∈ Succ s i , ∀s ∈ S (20)
r s i + p s i + l∈Succ s i p s l χ l,k + k -1 q=k+1 n l=1 p s l χ l,q + l∈P rec s j p s l χ l,k + p s j -M (2 -χ i,k -χ j,k ) ≤ C s j , ∀i, j, k, k ∈ {1, ..., n}, k > k, ∀s ∈ S (21) L max ≥ C s j -d s j , ∀j ∈ {1..n}, ∀s ∈ S (22)
Constraints ( 19) assign each job to exactly one group. Constraints ( 20) and ( 21) are used to compute the minimum value of the job completion time under each scenario as set by the online algorithm. For a given sequence, at the earliest, the completion time of a job J j is equal to the release date of some job J i plus the duration of the jobs that are between J i and J j in the sequence, these two jobs included. Constraints (20) considers the case were J i and J j are in the same group. In this case l∈(Succ s i ∩P rec s j ) p s l χ l,k is the total duration of the jobs that are scheduled between J i and J j under scenario s according to the ERD rule. Constraints (21) considers the case where J i and J j are in different groups. The expression l∈Succ s i p s l χ l,k + k -1 q=k+1 n l=1 p s l χ l,q + l∈P rec s j p s l χ l,k is the total duration of the jobs between jobs J i and J j . Constraints ( 22) bound the worst-case L max . This model contains n 2 binary variables, n|S| continuous variables and O(n 4 |S|) constraints.

Example

We consider the same example as in Section 4.1. An optimal solution for this example gives two groups {J 1 , J 3 , J 4 } and {J 2 , J 5 } with a maximum lateness L max = 0 instead of L max = 5 in the version without groups. Fig. 3 illustrates the solution for each scenario (in each group the jobs are presented in ERD order). 0 2 4 6 8 10 12 14 16 18 20

- J 3 Scenario s = 2 J 4 J 1 J 2 J 5 0 2 4 6 8 10 12 14 - Scenario s = 1 J 1 J 3 J 4 J 5 J 2 16 18 20 G J 1 G J 2 G J 1 G J 2 Figure 3: Groups G J 1 = {J 1 , J 3 , J 4
} and G J 2 = {J 2 , J 5 } under the two possible scenarios

MILP for the online recoverable robust scheduling and vehicle routing problems

In the first-stage solution x = (π, B, Σ) of the integrated robust scheduling and vehicle routing problem, π is a sequence of groups of permutable delivery jobs while σ(B) for each batch B ∈ B is a sequence of groups of permutable delivery jobs. The online recovery algorithm A(x, s) as the same behavior on production partial sequence π as in the previous section. For the delivery part, A(x, s) follows the order of the groups and when a delivery is completed inside a group it scans the remaining deliveries and select the nearest one according to the traveling time matrix of the realized scenario, as this travel time is assumed to be known with a sufficient precision at time of departure. This amounts to perform the routing inside a group according to the nearest neighbor rule. The following model is called OSRM and is based on model OSM for the scheduling part.

The binary variables associated to the vehicle routing part of the problem are the following: λ s j,k is equal to 1 if the site J j is in position k (i.e. the job J j is the k th job delivered) for scenario s, and 0 otherwise. δ k,r is equal to 1 if the tour r (i.e. the r th tour) visits the site in position k, and 0 otherwise, γ k,h is equal to 1 if the site in position k is in the group of permutable deliveries G D h , and a binary variable θ h,r equal to 1 if the group G D h is in the route r. For a scenario s, continuous variables D s j and α s r represent respectively the arrival date to the site J j and the departure date of the tour r.

According to these decision variables, we detail model ORSM below.

Objective (18) and constraints ( 19) -( 21) are taken from OSM. As already mentioned, Constraints ( 19) -( 21) set the production job completion times w.r.t. the production group orderings and the ERD rule inside each production group. The constraints related to the routing part are as follows. Constraints ( 23) and ( 24) bound from below the departure dates of the tours. Constraints (23) state that a tour r cannot start in scenario s before the completion of all the jobs which are to be delivered by r.

α s r ≥ C s j + M (2 -δ kr -λ s jk ), ∀j, r, k ∈ {1, ..., n}, ∀s ∈ S (23)
Constraints ( 24) prevent, for a given scenario s, tour r to start before the end of previous tour r -1, which is equal to the delivery date of the last visited site J j plus the return time to the depot.

α s r ≥ D s j + t s j0 + M (2 -δ k,r-1 -λ s jk ), ∀j, k ∈ {1, ..., n}, ∀r ∈ {2, ..., n}, ∀s ∈ S (24)
Constraints ( 25) to (31) below set the delivery dates of the jobs based on the computed batches, delivery group ordering inside each batch via the decision variables and on the nearest neighbor rule to schedule deliveries inside each groups. Constraints (25) prevent a job J i at the position k + 1 in scenario s to be delivered before the delivery of the job J j at the position k in scenario s.

D s j ≥ D s i + t s ij + M (2 -λ s i,k -λ s j,k+1
), ∀i, j ∈ {1, ..., n}, j = i, ∀k ∈ {1, ..., n -1}, ∀s ∈ S (25)

Constraints ( 26) prevent in turn to deliver a job J j at position k for scenario s before the departure time of the tour r that delivers it.

D s j ≥ α s r + ω s 0j -M (2 -δ kr -λ s j,k ), ∀j, r, k ∈ {1, ..., n}, ∀s ∈ S (26)
Constraints ( 27) ensure that the nearest neighbor rule is respected into each group, for each scenario. This is achieved by enforcing that the job J j delivered at position k in scenario s in a group h is the nearest in this scenario from the previous job J i (delivered at position k -1 in the same group h for the same scenario s) among all jobs J l delivered in group h.

D s j ≤ D s i + t s il + M (6 -γ k-1,h -λ s i,k-1 -γ k ,h -λ s l,k -γ k,h -λ s j,k
), ∀i, j, l, h ∈ {1, ..., n}, j = i, l = i, ∀k ∈ {2, ..., n}, ∀k ∈ {k, ..., n}, ∀s ∈ S (27) Constraints (28) considers a job J j at position k in scenario s which is the first delivered job of its group h for tour r. Then, if there is a preceding delivery group h -1 for the same tour r, J j must be the nearest among all other delivery jobs J l in group h in scenario s from the last job J i of group h -1.

D s j ≤ D s i + t s il + M (8 -γ k-1,h-1 -λ s i,k-1 -γ k ,h -λ s l,k -γ k,h -λ s j,k -θ h-1,r -θ h,r ), ∀i, j, l, r ∈ {1, ..., n}, j = i, l = i, ∀h, k ∈ {2, ..., n}, ∀k ∈ {k, ..., n}, ∀s ∈ S (28) (29)
Constraints (30) considers again a job J j at position k in scenario s which is the first delivered job of its group h for tour r, but in the case that group h is the first delivery group of tour r. In that case, if there exists a preceding tour r -1, which is assessed by the fact that some job J i is delivered at position k -1, then job J i must be the nearest from the depot.

D s j ≤ α s r + t s 0l + M (8 -γ k-1,h-1 -λ s i,k-1 -γ k ,h -λ s l,k -γ k,h -λ s j,k -θ h-1,r-1 -θ h,r
), ∀i, j, l ∈ {1, ..., n}, j = i, l = i, ∀r, h, k ∈ {2, ..., n}, ∀k ∈ {k, ..., n}, ∀s ∈ S (30)

Constraints (31) considers the last case for a job J j at position k in scenario s which is the first delivered job of its group h for tour r. In the case that group h is the first delivery group of tour r but when there is no preceding tour, which is assesed by the fact that job J j is delivered at position 1, then job J j must also be the nearest from the depot among all other delivery jobs l of its group.

D s j ≤ α s 1 + t s 0l + M (4 -γ 1,1 -λ s j,1 -γ k,1 -λ s l,k ) ∀j, l, k ∈ {1, ..., n}, ∀s ∈ S (31)
Since the variables γ k,h depend to the positions of the jobs, Constraints (32) ensure that a group contains the same jobs at each scenario.

γ k,h ≥ 1 -M (3 -λ s j,k -λ s j,k -γ k ,h ), ∀j, h, k, k ∈ {1, ..., n}, k = k , ∀s , s ∈ S, s > s (32)
Constraints ( 33)-( 41) ensure the feasibility of the tours. For a scenario s, each site is assigned to exactly one position (33) and each position to one site (34

n r=1 θ h,r = 1, ∀h ∈ {1, ..., n} (36)
δ kr + δ k+1,r ≤ 1, ∀k ∈ {1, ..., n -1}, ∀r, r ∈ {1, ..., n}, r = r, r = r + 1 (38)

γ kh + γ k+1,h ≤ 1, ∀k ∈ {1, ..., n -1}, ∀h, h ∈ {1, ..., n}, h = h, h = h + 1 (39) δ kr ≥ γ kh + θ hr , ∀k ∈ {1, ..., n -1}, ∀h, h ∈ {1, ..., n}, h = h, h = h + 1 (40) δ 11 = 1, γ 11 = 1, θ 11 = 1, ( 41
)
L max ≥ D s j -d s j , ∀j ∈ {1, ...,

Example

We consider here the same example as in Section 4.2 . We define the sequence of groups of permutable jobs G J 1 = {J 1 , J 3 , J 4 } and G J 2 = {J 2 , J 5 }. For the transportation phase, one batch contains jobs {J 1 , J 3 , J 4 } and the other batch contains the jobs {J 2 , J 5 }. In the first batch, there are two groups of permutable deliveries G D 1 = {J 3 } and G D 2 = {J 1 , J 4 } and in the second batch there is only one group G D 3 = {J 2 , J 5 }. Delivery phases are illustrated in Fig. 4 for each scenario and the global solutions that include the production and delivery phases are illustrated in Fig. 5 for each scenario. For this example, the maximum lateness over all scenarios is L max = 0.

0) (3) (0) - J 3 J 4 J 1 J 2 J 5 L max = 0 Scenario s = 2 (5) r 2 - (2) - (0) L max = 0 Scenario s = 1 J 1 J 3 J 4 J 5 J 2 0 2 (
4 6 8 10 12 14 16 18 20 22 24 26

- (0) r 1 (4) (3) -- (0) (2) r 2 -- (5) - (0) - r 1 (4) - - (1) (0)
(1) Although we were able to propose MILPs for both the standard robust and the online recovery robust integrated scheduling and routing problems, they can only be used to solve very small problem instances. To overcome this difficulty, we propose in this section two tabu search methods. The first one is called TSR (Tabu search for robust Scheduling & Routing) and solves the fixed sequence robust problem while the second one, called TOSR (Tabu search for Online recoverable robust Scheduling & Routing), deals with the online recovery robust problem. To make a fair comparison of the two robust solution paradigms, the two proposed algorithms are based on the same principles and can be presented jointly.

- ( 0
For TSR, a solution is a production job sequence on the machine, a sequence of delivery batches and, for each batch a delivery sequence (a routing). A solution is evaluated by computing for each scenario, the earliest start time of each production and delivery jobs following the prescribed production and delivery sequences, yielding a maximum delivery lateness over all scenarios.

For TOSR, a solution is a sequence of groups of permutable production jobs, a sequence of delivery batches and, for each batch a sequence of groups of permutable deliveries. However, to simplify the method, we restrict the solution space to the case where there is exactly one group of deliveries per batch, which means that no routing decisions have to be made in advance apart from batching. A solution is evaluated by computing for each scenario an earliest start time for each production job (following the prescribed production group sequence and the earliest release date first rule for scheduling the jobs inside each production group) as well as an earliest delivery time for each delivery job (following the nearest neighbor rule on the realized scenario for delivering the jobs inside each delivery group/batch). The maximum lateness over all scenarios can also be computed from these earliest start times.

In both cases, a quite standard version of the the Tabu Search algorithm has been implemented [Glover, 1990]. This metaheuristic starts from an initial complete solution and then searches the most appropriate solution among neighbors. A solution can be chosen if it is not tabu (i.e. not in the tabu list) or if it satisfies an aspiration criterion. The chosen neighbor becomes tabu and the search process is then repeated using this chosen neighbor as a new basic solution.

We propose a decomposition of the neighborhood search into two neighborhoods for TOSR and three neighborhoods for TSR that we search sequentially.

The two algorithms differ in the scheduling sub-problem where the neighbors are different as TOSR includes group-based neighborhoods and in the routing part where there is no decision to take in TOSR as the scenario-aware greedy nearest neighbor heuristic is used, while TSR ultimately optimizes the routes via the scenario-unaware 2-opt heuristic. We selected this policy (i.e. having an optimization process for the routing in TSR and not in TOSR) to have a fair comparison between the two approaches that have different strengths and weaknesses. On the one hand, the standard robust approach has the advantage of optimizing the routing part, while the recovery robust approach can not optimize any routing decision apart from the batching decision as only one group is allowed inside each batch. On the other hand the recovery robust approach has the possibility to adapt to the realized scenario inside each production group and each batch with the above-described online heuristics while the standard robust approach can only use worst-case evaluations. In the following sub-sections we details some key elements of both algorithms.

TOSR and TSR are sketched in Algorithms 6 and 6, respectively. Algorithm TSR uses the concept of reference scenario, which is one of the scenarios as explained later in Section 7.1. The following sections present different key components of the methods.

Algorithm 1 Tabu search algorithm (Online recoverable robust version -TOSR) COMPUTE INITIAL SOLUTION

• Generate a feasible solution (sequence of groups of permutable jobs and batching) using a greedy algorithm (details in Section 6.1.2).

WHILE THE GLOBAL STOP CONDITION IS NOT SATISFIED DO

• Production sequencing neighborhood -For a fixed batching and starting with the current production group sequence, search for the best non-tabu production group sequence move (details in Section 6.1.4).

-Store the newly selected solution in the scheduling tabu list.

• Batching neighborhood -For a fixed production group sequence and starting with the current batching, search the best batching move to get the best non-tabu neighbor solution (details in Section 6.1.4).

-Store the newly selected solution in the batching tabu list.

Algorithm 2 Tabu search algorithm (Robust version -TSR) COMPUTE INITIAL SOLUTION

• Generate a feasible solution (sequence of production jobs, batching, routing) using a greedy algorithm (details in Section 6.1.1).

WHILE THE GLOBAL STOP CONDITION IS NOT SATISFIED DO

• Production sequencing neighborhood -For a fixed batching and starting with the current production sequence, search the best production sequencing move to get the best non-tabu neighbor solution (routing inside each batch is computed with the NN heuristic on the reference scenario, details in Section 6.1.4).

-Store the newly selected solution in the scheduling tabu list.

• Batching neighborhood -For a fixed production sequence and starting with the current batching, search the best batching move to get the best non-tabu neighbor solution (routing inside each batch is computed with the NN heuristic on the reference scenario, details in Section 6.1.4).

-Store the newly selected solution in the batching tabu list.

• Routing sub-neighborhood -For a fixed production sequencing and batching and starting with the current routing, search the best permutation inside each route using the 2-opt heuristic (descent method only, details in Section 6.1.4).

Initial solutions

To initialize the methods, two constructive approaches are used. We first present the constructive approach for initializing the TSR algorithm called GSR (Greedy, Scheduling, Routing) then a second approach adapted to the TOSR algorithm called GOSR (Greedy, Online, Scheduling, Routing) is presented.

GSR: initial solution for TSR

In order to build an initial solution for the TSR algorithm, we first find a sequence of the jobs on the machine for the scheduling part. Then, from this sequence, we build batches. First, we sort the jobs according to their average due dates over the scenarios (Earliest Due Date rule -EDD) and from this sequence, batches are built by inserting the job in the current batch or by creating a new batch for this job. More precisely, Algorithm GSR is as follows:

step 1: Sort the jobs according to their average due dates over the scenarios.

step 2: There are two options for each job, the best option being selected at each iteration:

-Insert the job in the current batch.

-Create a new batch and insert the job.

GOSR: Initial solution for TOSR

A second constructive approach called GOSR is used to obtain an initial solution for TOSR. This time, we build simultaneously the groups for the scheduling part and the batches for the routing part. First we sort the jobs according to EDD and, from this sequence, production groups and batches are built at the same time by inserting the job in the current group or by by creating a new group for this job. The same procedure is done for batching: the job is inserted in the current batch or in a new batch. Finally, for each job, we keep the best solution found and we proceed to the next job. Heuristic GOSR is as follows:

step 1: Sort the jobs according to their average due dates over the scenarios.

step 2: There are four options for each job, the best option being selected at each iteration:

-Insert the job in the current production group; insert the job in the current routing group/batch.

-Insert the job in the current production group; create a new routing group/batch with the job.

-Create a new production group with the job; insert the job in the current routing group/batch.

-Create a new production group and a new current routing group/batch with the job.

Encoding

Encoding a solution for TOSR: In the TOSR algorithm, the problem of routing the jobs of a batch is systematically solved by the nearest neighbor algorithm (NN) on the realized scenario. In fact, since we enforce in addition that there is only one group of permutable delivery for each batch, it suffices to have the scheduling part and the composition of the batches to encode a solution. A solution of the problem S is consequently encoded directly by a vector of size 2n. The first part of size n contains the scheduling part, whereas the second part of size n contains the batching part. Let v denote the used vector for encoding a solution of the TOSR method. For the first part, we consider that v j for 1 ≤ j ≤ n is the index of the production group to which job J j is assigned, where the groups are sequenced in the order of their numbering. In the second part of the vector, we find the assignment to delivery groups/batches, where v j , for n + 1 ≤ j ≤ 2n, is the index of the delivery group/batch to which job J j-n is assigned. The delivery groups/batches are sequenced in the order of their numbering. The example in Fig. 6 shows a solution encoding corresponding to production groups

G J 1 = {J 1 , J 3 , J 4 } and G J 2 = {J 2 , J 5 } and delivery batches/groups G D 1 = {J 3 , J 4 }, G D 2 = {J 1 , J 5 }, G D 3 = {J 3 }.
Encoding a solution for TSR: the encoding for TSR is more subtile as for the first two neighborhoods, the routing is generated with NN on the reference scenario from a fixed production sequence and batching and, consequently, does not require an encoding. For the third neighborhood the routing requires an encoding for the 2-opt neighborhood. Hence for the first two neighborhoods, the solution is encoded with a vector v of size 2n where the first part of v corresponds to the position of the production jobs in the sequence where ( v i is the position of job J i in the schedule) and the second part represents the batching, as in the TOSR algorithm. For the third neighborhood, the vector has a third part of size n representing the routing, where v [2n + k] corresponds the the job delivered at position k in the concatenation of all batch routings. The encoding of the solution displayed on top of Figure 6 in TSR would be vector (3, 5, 1, 2, 4|2, 3, 1, 1, 2) for the first two neighborhood, with the adjunction of vector (3, 4, 5, 1, 2) to represent the routing in the third neighborhood.

-

J 1 J 3 J 4 - (3) -tour1 - (4) J 2 J 5 - (1) tour2 - (5) - (2) tour3 --
1 2 1 1 2 2 3 1 1 2 Encoding Figure 6: A solution for a scenario and an encoding for the production group sequence/batching

Neighborhood definition

We present in this section the neighborhoods used in each algorithm. The production scheduling neighborhoods are first presented for TOSR and TSR. Then, we present the batching neighborhood, which is the same for both algorithms (although evaluations of a neighbor differ, as already mentioned). Finally we present the 2-opt neighborhood used only in the TSR algorithm. We use a perimeter δ as a limit for the search area from a solution. For each neighborhood, the best non-tabu move is retained over all possible changes. Production group sequence neighborhood for TOSR:. We denote by |G| the number of groups. The following four neighborhoods are considered:

-Groups exchange: this neighborhood performs the exchange of two values in the first part of the encoding, i.e., the swap of the assignment of two production jobs, taking perimeter δ into consideration. Let v i denote the group of job J i while v j is the group of job J j . The exchange is possible if and only if

v i = v j and v i -δ|G| ≤ v j ≤ v i + δ|G|.
-Insertion of a job in a different group: this neighborhood allows a job J i of group v i to be inserted in an existing group

v k if v i -δ|G| ≤ v k ≤ v i + δ|G|.
-Group split: this neighborhood allows one group to be split into two groups, which postpone the subsequent groups. Since the sequence inside the groups depends on the realized scenario, the number of possible cases to split a group is exponential. Therefore, we order the jobs of a group according to the average due dates of the jobs over the scenarios, which reduces the number of cases to the size of the group.

-Group fusion: this neighborhood performs the merging of two consecutive groups into one, which advances the second group.

Production job sequence neighborhood for TSR: Two types of neighborhoods are used:

-Position sequence exchange: this neighborhood performs the exchange of two values v i and v j in the first part of the encoding. The perimeter δ is taken into consideration in the sense that we enforce that |v i -v j | ≤ δn.

-Insertion at a different position in the sequence: this neighborhood performs the exchange of the position v i of the job J i by inserting the job J i at a position v k , where, |v i -v k | ≤ δn. This yields that jobs at positions v i+1 , .., v k are advanced in the sequence if v i > v k or that jobs at positions v k , .., v i-1 are postponed if v i < v k .

Delivery batching neighborhood for TSR and Delivery grouping/batching neighborhood for TOSR The neighborhood generation is then the same for determining the delivery batches in TSR and the delivery groups/batches in TOSR. In fact, we have for a delivery batch (TSR) and for a delivery group/batch (TOSR) the same modifications as for a production group for TOSR. This yields four types of neighborhoods that can be derived from the TOSR production group sequence neighborhoods: batches exchange; insertion of a job in a different batch, batch split and batch fusion).

Routing neighborhoods for TSR This neighborhood is the standard 2-opt operator that swaps two jobs at different positions in the third part of the encoding, i.e. for each position j, k with 1 ≤ j < k ≤ n, the exchange of jobs v [2n + j] and v [2n + k] is evaluated.

Tabu list

To discourage the tabu algorithm to select the solutions that have been visited recently, the tabu list contains solutions that cannot be selected. These solutions are called Tabu. In our case, since the scheduling part and the batching part are treated separately, we opted for two tabu lists with a fixed size equal to the number of jobs. The first list contains the part of the vector corresponding to scheduling while the second list contains the part corresponding to the batches. We have no tabu list for the routing part in TSR as the 2-opt algorithm is applied as a descent method.

Aspiration criterion

A solution is considered tabu if either the scheduling part or the routing part is tabu. For this reason, a solution is selected if it is better than the best solution ever found, even if it is tabu as in standard tabu search.

Global stopping condition

The general algorithm stops when a global time limit is reached. However, for the scheduling subproblem and the batching sub-problem, there is no specific time limit but the algorithm stops after a fixed number of iterations without improvment.

Experimentations

Data generation

To validate our models, we have generated a random data set as follows. We call scenario s = 1 the "reference scenario" and we first generate instance (r 1 j , p 1 j , d 1 j ). For this, we first generate the p 1 j for all the jobs in the interval [1, 50] and we denote by P = n j=1 p 1 j . The release dates r 1 j are generated in the interval [0, γP ] where γ is a given parameter. Due dates d 1 j are generated in the interval [α -β 2 P, α + β 2 P ], where α and β are given parameters. For the generation of the distance matrix between sites, we generate randomly the coordinates (X 1 j , Y 1 j ) of each site j in the interval [0, 50] and the distance between two sites is the classical euclidian distance.

t 1 i,j = t 1 j,i = (X i -X j ) 2 + (Y i -Y j ) 2
For modeling the uncertainty on the data, each scenario s (s ≥ 2) is generated as follows:

• r s j is generated in [(1 -ω r )r 1 j , (1 + ω r r 1 j )] with a given coefficient ω r

• p s j is generated in [(1 -ω p )p 1 j , (1 + ω p p 1 j )] with a given coefficient ω p

• d s j is generated in [(1 -ω d )d 1 j , (1 + ω d d 1 j )] with a given coefficient ω d • X s j is generated in [(1 -ω X )X 1 j , (1 + ω X X 1 j )] with a given coefficient ω X • Y s j is generated in [(1 -ω Y )Y 1 j , (1 + ω Y Y 1 j )] with a given coefficient ω Y
The changes of the coordinates of the sites can be explained, if we consider that a site is a city and the change of position corresponds to different possible locations in the same city but it can also be seen as a simple way to vary the traveling time between two sites. For practical reasons, all the values of the coefficients ω are identical.

Computational results

In this section, we evaluate our tabu search algorithms on the instances presented above. The number of jobs varies from 10 to 100, the parameters α, β and γ are fixed to to 1. There are different classes of problems depending on the number of jobs n and the coefficient of variation ω between the scenarios. An instance is then represented by a tuple (n, ω), and, for each n and ω, 10 instances are generated with a number of scenarios S = 50. The practical limitations of the proposed MILP formulations has led us to omit the presentation of the results. However, for the solved instances with a number of jobs n = 10, the obtained results are equal to that obtained by the tabu search heuristic for both cases (TSR and TOSR). Note that, the fact that in the TOSR algorithm it is assumed that there are exactly one group of permutable delivery in each batch led us to adapt our formulation for the comparison.

The comparison takes place in two phases. In a "training" phase, we define a subset of scenarios (the training set) S ⊂ S . We compare on S the GSR and TSR algorithms for the fixed sequence robust approach with the GOSR and TOSR algorithms for the online robust approach. In a second phase, we keep the best production job sequence, batching and routing found by GSR and TSR on set S, and the best production group sequence and batching found by GOSR and TOSR on the same set. Then, we evaluate and compare the performance of these solutions on the scenario set S . This allows to evaluate these solutions on a scenario that was not anticipated and hence to provide additional informal information on the robustness of the methods.

The instances are solved on a 2.4 GHz Intel i3 computer with 3 GB of RAM.

Comparisons on the training scenario set S

In this section, the tabu search algorithms TSR and TOSR as well as the greedy algorithms GSR and GOSR presented above are tested on instances having a number of jobs from 10 to 100, a deviation ω ∈ {20, 40, 60} and a number of scenarios |S| ∈ {2, 5, 10}. The size of the tabu lists is equal to the number of jobs n and the perimeter δ is set to 0.4. For each (n, S, ω), the algorithms are tested for the ten instances and the results displayed in Table 1 give the average worst-case maximum lateness on the ten instances. The presented results are divided into three blocks. The first block provides the results of the fixed sequence robust approaches (GSR and TSR). The second one concerns the online recoverable robust approaches (GOST and TOSR) and the last one compares the fixed sequence and the online recoverable robust approaches. Let z(A) denote the objective function value (i.e. the maximum lateness over scenario set S) returned by algorithm A ∈ {TSR, TOSR, GSR, GOSR} on a particular instance. Columns ttb give the average CPU time needed. ∆ is the average gap between the initial solution found the the greedy approach and the tabu search solution and it is equal to (z(GSR) -z(T SR))/z(GSR) in the first block and (z(GOSR) -z(T OSR))/z(GOSR) for the second. In the last block, ∆ 1 = (z(GSR) -z(GOSR))/z(GSR) is the average gap between the initial solutions found by GSR and GOSR and ∆ 2 = (z(T SR) -z(T OSR))/z(T SR) is the average gap between the solutions found by TSR and TOSR. First, Column ∆ 1 reveals that the greedy online recoverable robust algorithm GOSR is far more better than the greedy fixed sequence robust algorithm GSR. This could be expected as the optimization level is very low in GSR whereas GOSR can react to the realized scenario thanks to the group structure. Second, the first and the second block show that for both approaches the tabu search algorithm improves significantly the performances of the greedy algorithms. Interestingly, TSR and TOSR have comparable CPU times on comparable instances, which means that considering an additional decision level for grouping does not yield a computational overhead. Finally, column ∆ 2 shows that in most cases the online recoverable robust approach improves upon the robust approach. However, this is not always the case and it can be seen that the improvement increases with the number of jobs and the deviation parameter, which indicates the uncertainty level.

Assessing the robustness of the obtained job/group sequences and batching on scenario S

In this part, we fix, for each instance, the best production job sequence, batching and routing found by GSR and TSR (denoted x GSR and x TSR ), and the best production group sequence and batching found by GOSR and TOSR (denoted x GOSR and x TOSR ), found during the optimization on the training set. We evaluate their performance by computing the worst case maximum lateness on the extended scenario set S .

The results presented in Table 2 are also divided into three blocks. The first block concerns the results of the fixed sequence robust solutions (x GSR and x TSR ), the second concerns the online recoverable robust solutions (x GOSR and x TOSR ) and the last one is also meant to compare the two approaches.

Let z(x) = max s∈S {L max (x, s)} for x ∈ {x GSR , x GOSR } the maximum lateness of the fixed sequence robust solutions on scenario set S and let z A (x) = max s∈S {L max (A(x, s), s)} the maximum lateness of the online recoverable sequence robust solutions on the same scenario set, where A(x, s) is the online algorithm used on the realized scenario (i.e. ERD for the production part and NN for the routing part). Columns ∆ give the average gap between the initial solution and the tabu search solution, which is equal to (z(x GSR ) -z(x TSR ))/z(x GSR ) in the first block and (z A (x GOSR ) -z A (x TOSR ))/z A (x GOSR ) in the second.

However considering each scenario in S as a possible outcome outside the training set S, the maximum performance on this set of each method separately has a limited interest and we would rather like to know how each solution issued from a different method behaves on the same unexpected scenario. Columns ∆ give the following average comparisons where

∆ = average s∈S L max (x GSR , s) -L max (x TSR , s) L max (x GSR , s)
for the fixed sequence robust approach and ∆ = average s∈S L max (A(x GOSR , s), s) -L max (A(x TOSR , s), s) L max (A(x GOSR , s), s)

for the online recoverable robust approach. The comparison between the two robust approaches is given in the third block where the first two columns ∆ (1) and ∆ (1) concern the initials solutions and ∆ (2) and ∆ (2) concern the tabu search solutions. ∆ (1) and ∆ (1) give the gap between the initials solutions when a solution is the worst case in S , and the average gap relative to each scenario, respectively. ∆ (2) and ∆ (2) indicate the same information for the tabu search solutions.

∆ (1) = z(x GSR )-z A (x GOSR ) z(x GSR ) ∆ (1) = average s∈S Lmax(x GSR ,s)-Lmax(A(x GOSR ,s),s) Lmax(x GSR ,s) ∆ (2) = z(x TSR )-z(x TOSR ) z(x TSR ) ∆ (2) = average s∈S Lmax(x TSR ,s)-Lmax(A(x TOSR ,s),s) Lmax(x TSR ,s)
Looking first at columns ∆ (1) and ∆ (1) , it appears that the solutions issued from the online recoverable robust approaches (either greedy or tabu) are far more robust than the solutions issued from their standardly robust counterparts, both in the worst case and average setting. Again the ability of the group structure to react to realized scenarios in a controlled way, even with a simple online algorithm, is demonstrated. Looking now at columns TSR and TOSR, if the improvement brought by TSR over GSR is again assessed, the average comparison shows that TOSR may provide solutions that perform worse than those of GOSR in average on scenario S . To exhibit a possible explanation of this behavior, Table 3 gives an aggregation of the results presented in Table 2. These results are divided into two blocks in which all algorithms are compared relatively to TOSR. More precisely, ∆ (1) and ∆ (4) compare GSR and TOSR by computing for each instance, the gap average s∈S L max (x GSR , s) -L max (A(x TOSR , s), s) L max (A(x TOSR , s), s) ∆ (2) and ∆ (5) compare GOSR and TOSR by computing for each instance, the gap average s∈S L max (A(x GOSR , s), s) -L max (A(x TOSR , s), s) L max (A(x TOSR , s), s)

TSR TOSR

n S ω ∆ (%) ∆ (%) ∆ (%) ∆ (%) ∆ (1) (%) ∆ (1) (%) ∆ (2) (%) ∆ (
and ∆ (3) and ∆ (6) compare TSR and TOSR by computing for each instance, the gap average s∈S L max (x TSR , s) -L max (A(x TOSR , s), s) L max (A(x TOSR , s), s)

The first block (columns ∆ (1) , ∆ (2) and ∆ (3) ) present the above-defined gaps averaged over all ω values (i.e. independently of the uncertainty level), while the second block (columns ∆ (4) , ∆ 3, where color blue represents the columns ∆ (1) and ∆ (4) , color red represents the columns ∆ (2) and ∆ (5) , and color green represents the columns ∆ (3) and ∆ (6) . From top to bottom in Figure 7 we have the results when s = 2, s = 5 and s = 10. In the same way, in Figure 8, we have the results when ω = 20%, ω = 40% and ω = 60%. Table 3 andFigures 7,8 show first that the greedy fixed sequence robust approach GSR is always outperformed by the other methods. Also the improvement brought by the online robust recoverable tabu search TOSR on the standard fixed sequence tabu search TSR is stable independently of ω and |S|. However the improvement brought by TOSR on the greedy recoverable robust approach is stable with the number of training scenarios but it decreases with the number of jobs and with the uncertainty level ω. GOSR is even slightly better than TOSR for 100 jobs and |S| = 5 as well for ω = 60. One explanation would be that when the uncertainty level is too high, the effort for optimizing the production group sequence becomes useless.

Conclusion and future research directions

In this paper we investigated a novel use of the concept of groups of permutable jobs, introduced in [Billaut and Roubellat, 1996], for an online recoverable robustness approach in an integrated production scheduling and delivery routing problem when uncertainty is modeled via discrete scenarios. The main novelty of the method is that it builds sequences of groups of production and delivery jobs that are fixed for all scenarios, while the job sequence inside each production or delivery group is constructed by an online greedy algorithm that uses in real time the available information on the realized scenario. With the aim to compare the standard robust scheduling approach with the new approach, we have proposed mixed integer linear programming formulations, greedy algorithms and tabu search methods for both approaches. On randomly generated instances of reasonable size, we compared the robustness of the heuristic approaches on a small size training scenario set and we compared the robustness of the solutions issued from the training phase on an extended scenario set. In both cases, the computational experiments show that for comparable CPU times the online recoverable robust approaches outperform the standard robust approaches. However, there seems to to be no more interest in optimizing the group sequences when the uncertainty level becomes too high.

-

Figure 2 :

 2 Figure 2: Robust solution for two scenarios

Figure 4 :

 4 Figure 4: Delivery group sequence {J 3 }, {J 1 , J 4 } and {J 2 , J 5 } yielding different routes on two possible scenarios

) -Figure 5 :

 )5 Figure 5: Production and transportation phases under the two possible scenarios

Figure 8 :

 8 Figure 8: Comparison for different Ω

  ). A site is assigned to exactly one tour (35) and one group (36), and a group can belong to only one tour (37). Constraints (38) serve to have consecutive tours, stating that if a job is delivered at position k by a tour r, the job delivered at position k+1 must be either in tour r or r+1. Constraints (38) serve to have consecutive delivery groups using the same principle. Constraints (40) ensure consistency between position assignment to groups and tours, stating that if a group h is assigned to a tour r and if group h contains a job delivered at position k, then tour r must also deliver the job at position k. Constraints (41) assign the job at the first position to the first group and to the first tour, and assign the first group to the first tour. Last, Constraints (42) bound L max .

	n j=1	λ s j,k = 1, ∀k ∈ {1, ..., n}, ∀s ∈ S	(33)
	n k=1	λ s j,k = 1, ∀j ∈ {1, ..., n}, ∀s ∈ S	(34)
	n h=1	γ k,h = 1, ∀k ∈ {1, ..., n}	(35)

n r=1 δ k,r = 1, ∀k ∈ {1, ..., n}

Table 1 :

 1 Comparisons of GSR, GOSR, TSR and TOSR methods for |S| ∈ {2, 5, 10}

				TSR		TOSR		TSR Vs TOSR
	n	S	ω	ttb(s)	∆(%)	ttb(s)	∆(%)	∆ 1 (%)	∆ 2 (%)
	10	2	20	0.01	50.42	0.08	23.02	38.31	5.58
			40	0.01	50.09	3.67	23.01	36.08	1.95
			60	0.01	51.39	0.20	26.98	28.75	-4.32
		5	20	7.79	50.27	0.18	22.78	35.77	1.30
			40	0.29	41.71	10.36	18.65	29.67	0.54
			60	0.05	33.46	0.13	17.32	27.46	10.24
		10	20	0.04	48.41	3.02	17.74	38.26	2.40
			40	0.08	41.12	0.64	16.31	31.19	-0.24
			60	0.21	37.97	8.47	13.74	26.85	-0.83
	20	2	20	2.70	70.48	17.76	32.90	56.38	2.44
			40	16.07	61.49	8.74	27.47	47.68	1.32
			60	19.33	61.67	6.30	22.15	51.43	3.78
		5	20	3.68	65.69	2.01	28.20	54.60	3.92
			40	9.52	56.03	10.87	21.66	47.08	4.72
			60	10.45	56.69	9.51	17.87	47.79	2.00
		10	20	28.82	62.76	14.32	30.32	48.98	4.88
			40	30.38	53.34	16.68	20.46	45.03	5.62
			60	23.60	50.09	5.55	14.81	44.12	5.01
	30	2	20	17.36	68.89	9.96	27.44	55.82	0.88
			40	30.93	64.96	11.63	22.54	55.30	4.50
			60	12.51	67.87	6.78	26.11	57.20	1.44
		5	20	25.88	66.36	27.59	28.28	54.56	3.59
			40	22.26	58.97	7.51	16.81	52.19	4.22
			60	42.83	59.24	24.24	22.11	49.20	3.43
		10	20	15.28	63.25	22.79	23.42	54.50	5.01
			40	25.35	56.45	16.90	16.47	48.80	2.49
			60	21.26	54.39	15.54	20.25	45.54	5.46
	50	2	20	24.47	71.33	2.49	21.29	64.22	2.31
			40	25.45	67.42	5.26	21.25	59.17	3.10
			60	14.56	65.27	10.02	18.55	59.16	4.59
		5	20	30.13	68.33	16.19	24.15	59.45	4.27
			40	24.16	61.62	27.54	19.81	54.49	5.26
			60	32.86	56.53	31.75	16.09	51.00	5.95
		10	20	50.60	65.60	32.79	21.16	59.24	6.76
			40	58.49	59.56	56.41	21.99	49.48	3.19
			60	47.13	54.79	55.31	18.01	48.37	6.44
	100	2	20	52.73	70.49	21.22	21.14	62.70	1.30
			40	59.31	64.46	25.72	14.83	60.54	5.41
			60	56.97	62.26	9.65	8.61	61.09	4.90
		5	20	109.51	65.84	64.50	15.02	61.21	4.39
			40	110.73	60.82	66.68	20.45	52.08	3.13
			60	113.04	54.41	64.40	13.80	49.11	3.77
		10	20	119.02	58.36	116.47	16.91	58.37	16.15
			40	120.63	55.47	117.90	16.01	51.64	9.28
			60	122.56	44.55	127.80	13.58	47.90	17.82

Table 2 :

 2 Assessing the robustness of the GSR, GOST, TSR and TOSR solutions on scenario set S

	2) (%)

Table 3 :

 3 5) and ∆ (6) ) give the gaps averaged over over all |S| values (i.e. independently of the training set size level). Aggregated comparisons of the solutions on scenario set S Figures 7 and 8 synthesize table

	n	S	∆ (1) (%)	∆ (2) (%)	∆ (3) (%)	ω	∆ (4) (%)	∆ (5) (%)	∆ (6) (%)
	10	2	63.90	10.13	16.43	20	95.77	21.89	14.34
	20		106.61	3.42	23.40		193.07	27.76	25.68
	30		110.08	2.82	29.07		179.02	18.63	28.28
	50		120.51	1.11	38.31		187.30	13.24	28.61
	100		120.05	2.74	42.43		172.49	3.75	30.63
	10	5	69.23	11.37	18.33	40	68.10	12.83	20.05
	20		128.59	10.81	22.72		98.96	5.25	18.32
	30		129.54	5.52	28.77		110.32	0.06	26.66
	50		133.52	5.82	26.72		123.56	5.35	31.45
	100		120.44	-2.36	24.93		115.30	1.35	31.23
	10	10	72.94	12.95	16.35	60	42.20	-0.27	16.72
	20		129.86	15.52	20.50		73.04	-3.27	22.62
	30		136.92	9.61	20.90		87.21	-0.74	23.79
	50		149.37	11.65	23.63		92.53	-0.01	28.59
	100		132.97	2.67	32.42		85.68	-2.05	37.92

They called this structure an "ordered sequence of subsets", which is equivalent to the "sequence of groups of