Philippe Lacomme

Aziz Moukrim

Alain Quilliot

Marina Vinot

Supply Chain Optimization with both Production and Transportation Integration: Multiple Vehicles for a Single Perishable Product

Keywords: supply chain coordination, transportation, scheduling, vehicle routing problem, integration

HAL is a

Introduction and literature review

Introduction to integrated problems

The Production and Transportation Scheduling Problem (PTSP) modelled a supply chain problem with production and transportation where product perishability must be addressed (which is common in the food, chemical and pharmaceutical industries). The PTSP is a problem where both scheduling and routing must be jointly solved to have a proper coordination between the production on a single production facility and the transportation, taking the product lifespan into consideration. In this problem, once a lot of products is produced, it must be directly transported to various customer sites within its limited lifespan. A solution to a PTSP is composed of production operations (starting times have to be computed) and transportation operations that define sub-trips (terminology used in the routing community) from the depot (production facility) to a set of customers.

Supply chain management is a cross-functional approach that includes managing the movement of raw materials into an organisation, the processing of materials into finished products (production), and the movement of finished products out of the organisation and toward the end consumers (transport). All these functions (production, storage, transportation) can be seen as independent subproblems or can be integrated under one plan in supply chains. In the literature, a tremendous amount of research considers the production and transportation subproblems successively, and in the vast majority of companies, the scheduling problem is solved first and the routing problem is then addressed, although this kind of approach does not lead to an optimal solution.

Production and transportation stages at the planning level are often linked by intermediate inventory stages whose management strongly depends on the proper coordination between production and transportation as an integrated resolution. Note that integrated production and distribution has a significant impact on customer service, favouring product delivery in due time. The connection between the production (processing stations) and transportation is a highly desirable goal as stressed in numerous publications, including [START_REF] Belo-Filho | An adaptive large neighbourhood search for the operational integrated production and distribution problem of perishable products[END_REF], with a problem with several production facilities and perishable products.

The integration of production and transportation can be relevant when the products are perishable. The product perishability is related to an "expiration date", which commonly appears on product labels to indicate the date from which the quality of the product is no longer guaranteed by the manufacturer (Rivera and Lallmahomed, 2015). [START_REF] Roscoe | Supply chain segmentation in the sporting goods industry[END_REF] underline that perishable goods are quite specific and influence the packing, the storage and the transportation mode. The management of perishable products is critical for integrating production and transportation scheduling in a coordinated manner.

As stressed by [START_REF] Chen | Integrated Production and Outbound Distribution Scheduling: Review and Extensions[END_REF], the integrated production and transportation models at a detailed scheduling level are fairly recent and the majority of models attempt to jointly optimise job-byjob production and transportation by considering the customer service level at the individual job level. Integrated production and routing problems are receiving more and more attention, and several surveys are available, including [START_REF] Moons | Integrated production scheduling and vehicle routing decisions at the operational decision level: A review and discussion[END_REF], [START_REF] Chen | Integrated Production and Outbound Distribution Scheduling: Review and Extensions[END_REF] and Saramiento et al. (1999). Let us note that integrated production and transportation problems fall into two categories of problems according to Kuhn and Liske ( 2017): (a) a Production Distribution Problem; or (b) a Sourcing Production Problem. A large number of publications consider the Production Distribution Problem, i.e., problems where the production is completed first and the products are distributed afterwards. [START_REF] Chen | Integrated Production and Outbound Distribution Scheduling: Review and Extensions[END_REF] provides a survey on integrated production and outbound distribution scheduling problems and classifies these existing problems into several different classes. These problems often have different names: Integrated Production and Outbound Distribution Scheduling (IPODS), Integrated Production and Distribution Problem (IPDP) or, in our case, Production and Transportation Scheduling Problem (PTSP). In these problems, a solution is characterised by a production schedule that plans the date and the location where each demand is processed, and a delivery schedule that satisfies the orders in each shipment, the number of vehicles and the delivery date. Problems can vary depending on the number of production facilities with machine configuration (single-machine, parallel-machine, flow shop, etc.), the number of customers, and customers' characteristics, including delivery time windows and deadlines. The number of available vehicles and their capacities may be limited or not. [START_REF] Chen | Integrated Production and Outbound Distribution Scheduling: Review and Extensions[END_REF] classifies existing IPODS problems into five classes considering individual and immediate delivery, batch delivery to a single customer by direct shipping, batch delivery to multiple customers by direct shipping, batch delivery to multiple customers by routing, and fixed delivery departure dates.

Intregrated problem with one production facility and one perishable product

Perishability naturally occurs in the food, chemical and pharmaceutical industries and its management creates several complications in the scheduling of production and transportation. Special efforts are required to reduce the waste and cost of the storage and transportation of perishable products. [START_REF] Chen | Integrated Production and Outbound Distribution Scheduling: Review and Extensions[END_REF] mentions studies on the perishable products with timesensitive constraints [START_REF] Armstrong | A zero-inventory production and distribution problem with a fixed customer sequence[END_REF][START_REF] Devapriya | Optimal fleet size of an integrated production and distribution scheduling problem for a perishable product[END_REF][START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF]. The problems studied in these articles belong to the same class: batch delivery to multiple customers with routing. The time-sensitive constraints in these problems are due to industrial chemical compounds that must be delivered within a certain time limit once produced. In the survey of [START_REF] Moons | Integrated production scheduling and vehicle routing decisions at the operational decision level: A review and discussion[END_REF] on integrated production scheduling and vehicle routing problems, two studies with one perishable product and one production facility are mentionned [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF][START_REF] Devapriya | Integrated production and distribution scheduling with a perishable product[END_REF]. [START_REF] Armstrong | A zero-inventory production and distribution problem with a fixed customer sequence[END_REF] study the production and transportation problem with a single vehicle and a single production facility. The order in which customers may receive deliveries is fixed. Each customer requests a delivery quantity with a time window for receiving it. The lifespan of the product begins as soon as the production of a customer's order is completed. The problem then turns into a routing problem, minimising the number of non-serviced customers for a subset of customers. [START_REF] Devapriya | Optimal fleet size of an integrated production and distribution scheduling problem for a perishable product[END_REF] focus on a problem with one production facility and a large fleet size. The lifespan of the perishable products immediately begins when the last order of a batch has been completed. The delivery of a batch order therefore begins as soon as the production is achieved. The objective is to minimise the total transportation cost, i.e., the cost of delivery time and the number of vehicles required to satisfy all the demands. Table 1. Integrated problem with one production facility and one perishable product. et al. (2008) address a problem (PTSP) with a network of multiple customers, a single facility with a constant production rate 𝑟 > 0, and a single vehicle with capacity 𝑄. The objective is to minimise the time required to serve all customers, commonly referred to as the makespan. Furthermore, due to a product with a lifespan 𝐵, waiting time must be addressed and the problem can be considered as a two-machine flow shop with maximal time lags. Let us note that in [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF], the problem focuses on a situation where the lifespan is defined by an expiration delay between the end of production and customer delivery. To jointly solve the scheduling and routing problem, the authors introduce a genetic and a memetic algorithm with the algorithm of [START_REF] Gilmore | Sequencing a one-state variable machine: A solvable case of the traveling salesman problem[END_REF] and several lower bounds. The problem has been proved to be NP-hard in the strong sense.

• • • • • • • • Armstrong et al. (2008) • • • • • • • • • • • Karaoğlan and Kesen (2017) • • • • • • • • • Devapriya (2006) • • • • • • • • • Devapriya et al. (2017) • • • • • • • • • •

Geismar

The problem described by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] was recently solved by [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF] using a branch-and-cut (B&C) algorithm. The algorithm uses several valid inequalities taken from the existing literature, and a local search based on a simulated annealing approach is used to improve upper bounds. Numerical experiments are achieved on the same instances as [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] and prove that the B&C is strongly competitive. The objective function consists in the minimisation of the makespan (Table 1).

Recently, [START_REF] Devapriya | Integrated production and distribution scheduling with a perishable product[END_REF] addressed an extension (referred to as IPDSP) of the PTSP introduced by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] where the fleet size is a decision variable. [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] base their approach on an indirect representation of the solutions with giant trips. The difficulty of integration is addressed with the dedicated [START_REF] Gilmore | Sequencing a one-state variable machine: A solvable case of the traveling salesman problem[END_REF] algorithm to solve the no-wait flow shop problem and obtain an optimal order of subtrips. The final solution is then deduced by relaxing the no-wait constraint, keeping the constraints on product perishability. Let us note that for the PTSP addressed by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF], the solution is composed of one trip with several ordered sub-trips with increasing starting times since there is no assignment problem of vehicles to be solved.

Solution representation for scheduling and routing problems

The PTSP with the constraint of perishable products can be modelled as a permutation flow shop scheduling problem on two machines with maximal time lags. [START_REF] Fondrevelle | Permutation Flowshop Scheduling Problems with Maximal and Minimal Time Lags[END_REF] prove that this problem is strongly NP-hard even when all the maximal time lags have the same positive value. Nevertheless, note that this problem is a mix between the classical and the nowait flow shop scheduling problem that can be polynomially solved with the algorithm of [START_REF] Johnson | Optimal Two and Three Stage Production Schedules with Set-Up Times, Included[END_REF] and with the algorithm of [START_REF] Gilmore | Sequencing a one-state variable machine: A solvable case of the traveling salesman problem[END_REF], respectively. The first research on the two-machine flow shop problem was proposed by [START_REF] Johnson | Optimal Two and Three Stage Production Schedules with Set-Up Times, Included[END_REF] who defined rules to obtain the optimal schedule. It is commonly accepted that most of the scheduling problem-solving approaches take advantage of a modelling based on the disjunctive graph of [START_REF] Roy | Les problèmes d'ordonnancement avec contraintes disjonctives[END_REF]. Numerous authors have introduced approaches for graph generation and search space exploration. [START_REF] Bierwirth | A generalized permutation approach to job-shop scheduling with genetic algorithms[END_REF] proposal for the job shop remains within the global trend of indirect representation schemes and proves that it is possible to express the machine selections as a vector by repetitions that define a topological order of nodes. Similar remarks hold for the routing problem that has received a considerable amount of attention in recent years. The idea of splitting a giant trip was introduced by [START_REF] Beasley | Route-first cluster-second methods for vehicle routing[END_REF] and was first included in a global framework for the routing problem by [START_REF] Lacomme | Competitive Memetic Algorithms for the Capacitated Arc Routing Problem and its Extensions[END_REF]. The total number of methods that take advantage of such an approach has strongly increased in recent years (as mentioned in the state of the art of [START_REF] Prins | Order-first split-second methods for vehicle routing problems: A review[END_REF]). The second key feature used in routing approaches lies in the local search procedure, which is commonly based on a swap within a trip, a 2-Opt within a trip, a swap between two trips and a 2-Opt between two trips (see [START_REF] Lacomme | Competitive Memetic Algorithms for the Capacitated Arc Routing Problem and its Extensions[END_REF] and [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]).

The remainder of this paper is structured as follows. The following section introduces the problem of interest with multiple vehicles. Section 3 presents a GRASP×ELS framework designed to address the key features for an integrated resolution. Section 4 provides the computational evaluation of the problem by considering a single vehicle. Section 5 illustrates the computational evaluation of the problem with multiple vehicles. The last section provides some conclusions and avenues for further research.

Problem definition for the PTSPm

Definition

The PTSPm (PTSP with multiple vehicles) is an extension of Geismar et al.'s (2008) To solve this problem, three subproblems must be jointly solved: the assignment, the routing and the scheduling. When the assignment and the routing problems are solved, the residual problem can be modelled as a two-stage hybrid flow shop with maximal time lags where the machine in the first stage is the single production facility and where each vehicle defines a machine in the second stage. The [START_REF] Johnson | Optimal Two and Three Stage Production Schedules with Set-Up Times, Included[END_REF] algorithm and the [START_REF] Gilmore | Sequencing a one-state variable machine: A solvable case of the traveling salesman problem[END_REF] algorithm can no longer be used to solve the scheduling part of the PTSPm. A solution could be composed by a set of jobs 𝐽, with 𝑐𝑎𝑟𝑑(𝐽) = 𝑛𝑗 with 𝑛𝑗 ≤ 𝑛. A job 𝑗 includes a set of customers to be served and is composed of two operations 𝑂 𝑗1 and 𝑂 𝑗2 𝑘 .

Operation 𝑂 𝑗1 is a production operation assigned to the production facility, modelling the production of the demands of the customers in the job 𝑗. The operation 𝑂 𝑗2 𝑘 is a transportation operation assigned to the vehicle 𝑘, modelling the sub-trip that makes it possible to deliver the demands of the customers in job 𝑗 from the production facility, also referred to as the depot.

The trip of vehicle 𝑘 is defined by a sequence of all the transport operations (sub-trips) assigned to vehicle 𝑘. The following constraints must hold:  each transportation operation begins and ends at the production facility, also referred to as the depot or depot node, depending on the routing terminology;  the order of the production operations and of the transportation operations should be identical for a given vehicle between two jobs;  a vehicle achieves at most one trip;  all customers of a sub-trip 𝑂 𝑖2

𝑘 must be served within 𝐵 time units after the end of the production operation 𝑂 𝑖1 ;  the total demands in a sub-trip 𝑂 𝑖2 𝑘 cannot exceed the vehicle capacity ∑ 𝑞 𝑗 𝑗∈𝑂 𝑖2 𝑘 ≤ 𝑄;

 deliveries must not be split (each customer must be delivered only once).

The following notations are used to caracterise a solution: 𝐶 𝑗 set of customers included in job 𝑗 𝑂 𝑗1 production operation of job 𝑗 𝑂 𝑗2

𝑘 transportation operation of job 𝑗 assigned to vehicle 𝑘 𝑝 𝑗1 duration of the production operation of job 𝑗 𝑝 𝑗2 duration of the transportation operation of job 𝑗 𝑝 𝑗2 ′ duration of the transportation operation of job 𝑗 without the empty transport to the depot 𝑠 𝑗1 starting time of the production operation of job 𝑗 𝑠 𝑗2 starting time of the transportation operation of job 𝑗 𝑓 𝑗1 finishing time of the production operation of job 𝑗 𝑓 𝑗2 finishing time of the transportation operation of job 𝑗 Figure 1 represents a solution of a PTSPm with two vehicles and five customers. In this solution, the customers are divided into three groups to define three jobs:  the first job (with the operations 𝑂 11 and 𝑂 12 1 ) includes customers 1 and 2;  the second job (with the operations 𝑂 21 and 𝑂 22

2 ) contains customer 3;  the third job (with the operations 𝑂 31 and 𝑂 32 1 ) includes customers 4 and 5; The trip of vehicle 1 is composed of two transportation operations (𝑂 12 1 and 𝑂 32 1 ) and the trip of vehicle 2 is composed of one transportation operation 𝑂 22 2 . The solution presented on the Gantt chart uses the following values for the demands of the customers, 𝑞 1 = 1, 𝑞 2 = 1, 𝑞 3 = 6, 𝑞 4 = 1, 𝑞 5 = 3, and the transportation times are equal to 𝜏 0,1 = 𝜏 2,0 = 𝜏 0,4 = 𝜏 4,5 = 3, 𝜏 1,2 = 𝜏 0,3 = 2 and 𝜏 5,0 = 1. For the first job, the duration of the production operation 𝑝 1 = 2 time units since the production rate 𝑟 = 1 and 𝑞 1 = 𝑞 2 = 1, and the duration of the transportation operation 𝑡 1 = 8 time units since 𝜏 0,1 + 𝜏 1,2 + 𝜏 2,0 = 8. The lifespan of the product is repected: 𝑓 𝑗2 -𝑓 11 -𝜏 2,0 = 10 -2 -3 = 5 ≤ 𝐵 = 7.

Let us note that τ x,y in this problem denotes a transportation time between customer x and customer y that can be precomputed by considering either the minimal distance traveled by the vehicle between customers or the minimal time necessary to travel between two customers that are equal since the speed factor for the vehicles is equal to one. The jobs have to be defined and ordered in order to minimise the makespan 𝐶 𝑚𝑎𝑥 , i.e., the arrival time of the last vehicle at the depot.

p 11 = (q 1 + q 2 )/r p 12 = τ 0,1 + τ 1,2 + τ 2,0 p 21 = q 3 /r p 22 = τ 0,3 + τ 3,0 p 31 = (q 4 + q 5 )/r p 32 = τ 0,4 + τ 4,5 + τ 5,0 Numerical example with: Q = 10, B = 7 and r = 1

1 2 3 4 5 0 q 1 q 2 q 4 q 5 τ 0,1 τ 2,

Linear formulation of the PTSPm

The combined production scheduling and vehicle routing problem is formulated as a mixed integer linear programming model based on five binary variables (𝑦 𝑏𝑖 , ℎ 𝑏 , 𝑥 𝑖𝑗 𝑏 , 𝑧 𝑏𝑐 , 𝑎 𝑏𝑘 ) and seven integer variables (𝑠 𝑏1 , 𝑠 𝑏2 ,𝑝 𝑏1 , 𝑝 𝑏2 , 𝑝 𝑏2 ′ , 𝑞 𝑖 -, 𝐶 𝑚𝑎𝑥 ). The binary decision variables are related to the disjunctions between operations and the integer variables are related to the starting and finishing times or duration of the operations.

Notations for the binary variables:

𝑦 𝑏𝑖 = { 1 0 if customer 𝑖 is in job 𝑏 otherwise 𝑖 ∈ 𝐸 𝑏 ∈ 𝐽 ℎ 𝑏 = { 1 0 if job 𝑏 is composed of at least one customer otherwise 𝑏 ∈ 𝐽 𝑥 𝑖𝑗 𝑏 = { 1 0 if customer 𝑖 is serviced immediately before customer 𝑗 in job 𝑏 otherwise (𝑖, 𝑗) ∈ 𝐸 2 𝑏 ∈ 𝐽 𝑧 𝑏𝑐 = { 1 0 if job 𝑏 is scheduled before job c otherwise (𝑏, 𝑐) ∈ 𝐽 2 𝑎 𝑏𝑘 = { 1 0 if job 𝑏 is assigned to vehicle 𝑘 otherwise 𝑏 ∈ 𝐽 𝑘 ∈ 𝑉
Job assignment requirement: this set of constraints ensures that one customer is assigned to one and only one job.

∀𝑖 ∈ 𝐸 ∑ 𝑦 𝑏𝑖 𝑏∈𝐽 = 1 (1)
Vehicle capacity requirements: this constraint ensures that the customers assigned to a job can be serviced by one vehicle considering the total amount of demands of customers.

∀𝑏 ∈ 𝐽 ∑ 𝑦 𝑏𝑖

𝑖∈𝐸

. 𝑞 𝑖 ≤ 𝐶 (2) Job definition: this constraint ensures that if job 𝑏 encompasses no customer (ℎ 𝑏 = 0), no customer is assigned to the job (𝑦 𝑏𝑖 = 0) ∀𝑏 ∈ 𝐽, ∀𝑖 ∈ 𝐸 𝑦 𝑏𝑖 ≤ ℎ 𝑏 (3) Job duration for the production: this constraint defines the job processing time 𝑝 𝑏 on the production considering the set of customers assigned to the job. If customer 𝑖 is assigned to job 𝑏, then 𝑦 𝑏𝑖 = 1 and ∑ 𝑞 𝑖 . 𝑟 𝑖∈𝐸 is the sum of 𝑞 𝑖 moderated with the production rate 𝑟. ∀𝑏 ∈ 𝐽 𝑝 𝑏 = ∑ 𝑦 𝑏𝑖 𝑖∈𝐸 . 𝑞 𝑖 . 𝑟 (4) Vehicle assignment to one vehicle: this constraint ensures that if job 𝑏 is used (there are one or several customers assigned to 𝑏), i..e., ℎ 𝑏 = 1, then one variable 𝑎 𝑏𝑘 is assigned to 1, i.e., one vehicle 𝑘 is assigned to job 𝑏. ∀𝑏 ∈ 𝐽 ℎ 𝑏 = ∑ 𝑎 𝑏𝑘

𝑘∈𝑉

(5) Customer order in the jobs: this set of constraints ensures that each customer has a single predecessor and a single successor

∀𝑗 ∈ 𝐸 ∑ ∑ 𝑥 𝑖𝑗 𝑏 𝑖∈𝐸∪{0} 𝑖≠𝑗 𝑏∈𝐽 = 1 (6) ∀𝑖 ∈ 𝐸 ∑ ∑ 𝑥 𝑖𝑗 𝑏 𝑗∈𝐸∪{0} 𝑗≠𝑖 𝑏∈𝐽 = 1 (7)
Depot definition: this constraint ensures that if job 𝑏 encompasses at least one customer (ℎ 𝑏 = 1), the sub-trip starts and ends at the depot.

∀𝑏 ∈ 𝐽 ℎ 𝑏 = ∑ 𝑥 0𝑗 𝑏 𝑗∈𝐸 (8) ∀𝑏 ∈ 𝐽 ℎ 𝑏 = ∑ 𝑥 𝑖0 𝑏 𝑖∈𝐸 (9)
Customer assignment to a job: this constraint ensures that if a customer 𝑖 is assigned to job 𝑏 (𝑦 𝑏𝑖 = 1), then customer 𝑖 has one predecessor, i. e., ∃𝑗/𝑥 𝑖𝑗 𝑏 = 1. Similar remarks hold for constraint 13, considering the successor.

∀𝑖 ∈ 𝐸, ∀𝑏 ∈ 𝐽 𝑦 𝑏𝑖 = ∑ 𝑥 𝑖𝑗 𝑏 𝑗∈𝐸∪{0} 𝑗≠𝑖 (10) ∀𝑖 ∈ 𝐸, ∀𝑏 ∈ 𝐽 𝑦 𝑏𝑖 = ∑ 𝑥 𝑗𝑖 𝑏 𝑗∈𝐸∪{0} 𝑗≠𝑖 (11)
Sub-trip eliminations: these constraints are the Miller-Tucker-Zemlin (MTZ) constraints. This constraint uses 𝑞 𝑗 -, referred to as the total amount of products remaining in the vehicle after

servicing customer 𝑗. ∀𝑏 ∈ 𝐽, ∀(𝑖, 𝑗) ∈ 𝐸 2 , 𝑖 ≠ 𝑗 𝑞 𝑗 --𝑞 𝑖 -+ 𝑥 𝑖𝑗 𝑏 . 𝐶 ≤ 𝐶 -𝑞 𝑗 (12) ∀𝑏 ∈ 𝐽, ∀(𝑖, 𝑗) ∈ 𝐸 2 , 𝑖 ≠ 𝑗 𝑞 𝑖 -≤ 𝐶 -𝑞 𝑖 (13
) Job duration for transport: this constraint defines the sub-trip duration to service all customers assigned to job

𝑏 ∀𝑏 ∈ 𝐽 𝑝 𝑏2 = ∑ ∑ 𝑥 𝑖𝑗 𝑏 𝑗∈𝐸∪{0} 𝑗≠𝑖 . 𝜏 𝑖𝑗 𝑖∈𝐸∪{0} ( 14
)
Job duration for transport: this constraint defines the trip duration to service all customers assigned to job 𝑏 not considering the empty transport from the last customer to the depot. ∀𝑏 ∈ 𝐽 𝑝 𝑏2 ′ = 𝑝 𝑏2 -∑ 𝑥 𝑖0 𝑏 . 𝜏 𝑖0 𝑖∈𝐸 (15) Disjunctive constraints at the production end: this constraint ensures that the production operation of job 𝑐 and the production operation of job 𝑑 cannot be performed at the same time at the production facility.

∀(𝑏, 𝑐) ∈ 𝐽 2 , 𝑐 ≠ 𝑏 𝑠 𝑏1 + 𝑝 𝑏1 ≤ 𝑠 𝑐1 + (1 -𝑧 𝑏𝑐 ). 𝐻 (16) ∀(𝑏, 𝑐) ∈ 𝐽 2 , 𝑐 ≠ 𝑏 𝑠 𝑐1 + 𝑝 𝑐1 ≤ 𝑠 𝑏1 + 𝑧 𝑏𝑐 . 𝐻 (17)
If 𝑧 𝑏𝑐 = 1 (job 𝑏 is scheduled before job 𝑐), constraint (18) can be rewritten as 𝑠 𝑏1 + 𝑝 𝑏1 -𝑠 𝑐1 ≤ 0, meaning that 𝑠 𝑏1 + 𝑝 𝑏1 ≤ 𝑠 𝑐1 , ensuring that production of 𝑐 on the production facility cannot start before the end of the production operation of job 𝑏. If 𝑧 𝑏𝑐 = 0, constraint (18) holds, and constraint ( 19) can be rewritten as 𝑠 𝑐1 + 𝑝 𝑐1 ≤ 𝑠 𝑏1, meaning that job 𝑐 is processed first and job 𝑑 second for the production operations.

Disjunctive constraints for transport: this constraint ensures that the transport operations of two jobs 𝑏 and 𝑐 are not performed at the same time by the same vehicle 𝑘.

∀(𝑏, 𝑐) ∈ 𝐽 2 , 𝑐 ≠ 𝑏 𝑠 𝑏2 + 𝑝 𝑏2 ≤ 𝑠 𝑐2 + (3 -𝑧 𝑏𝑐 -𝑎 𝑏𝑘 -𝑎 𝑐𝑘 ). 𝐻 (18) ∀(𝑏, 𝑐) ∈ 𝐽 2 , 𝑐 ≠ 𝑏 𝑠 𝑐2 + 𝑝 𝑐2 ≤ 𝑠 𝑏2 + (2 -𝑧 𝑏𝑐 + 𝑎 𝑏𝑘 + 𝑎 𝑐𝑘 ). 𝐻 (19) If 𝑎 𝑏𝑘 = 𝑎 𝑐𝑘 = 1,
this means that the two transport operations are assigned to the same vehicle 𝑘 and the constraints can be rewritten as: 𝑠 𝑏2 + 𝑝 𝑏2 ≤ 𝑠 𝑐2 + (1 -𝑧 𝑏𝑐 ). 𝐻 (20) and 𝑠 𝑐2 + 𝑝 𝑐2 ≤ 𝑠 𝑏2 + 𝑧 𝑏𝑐 . 𝐻 (21). If either 𝑎 𝑏𝑘 or 𝑎 𝑐𝑘 are not assigned to 1, the two constraints hold, regardless of the value of 𝑧 𝑏𝑐 . Precedence constraints per operation: this constraint ensures that the transport operations of one job are performed according to the production-operation sequence first, followed by the transport-operation sequence.

Lifespan products: this constraint ensures that the lifespan of the products are addressed. The product lifespan defines a maximal delay between the delivery of the last customer of a subtrip of a job b, i.e., 𝑠 𝑏2 + 𝑝 𝑏2 ′ , and the finishing time of the production referred to as 𝑠 𝑏1 + 𝑝 𝑏1 : the difference is upper bounded by 𝐵. ∀𝑏 ∈ 𝐽 𝑠 𝑏2 + 𝑝 𝑏2 ′ ≤ 𝑠 𝑏1 + 𝑝 𝑏1 + 𝐵 (21) Makespan. When minimising the makespan, the following constraints should be added to define new integer variables. This constraint ensures that

𝐶 𝑚𝑎𝑥 = 𝑚𝑎𝑥 𝑏∈𝐽 (𝑑 𝑏2 + 𝑝 𝑏2 ) ∀𝑏 ∈ 𝐽 𝑠 𝑏2 + 𝑝 𝑏2 ≤ 𝐶 𝑚𝑎𝑥 (22)

Modelling

In the PTSPm, a sub-trip is fully defined by an ordered sequence of operations: (1) starting at the depot with a pickup operation;

(2) defining an ordered sequence of delivery operations; and

(3) finishing at the depot. The loading of the vehicle can be represented as a decreasing function (see Fig. 2) during each sub-trip. Figure 2 gives some details of the solution presented in Fig. 1. Figure 2 represents both the vehicle load (on the top) and the proper coordination between production and transportation over time, with an explicit modelling of pickup and delivery operations (at the bottom). Time constraints on sub-trips Time dependency only exists between successive sub-trips of the same trip since all subtrips of one trip are assigned to the same vehicle. Concerning the trip of vehicle 1 in Fig. 2, the previous remark implies that 𝑠 32 ≥ 𝑓 12 . Moreover, the earliest starting time of 𝑂 31 depends on the starting time 𝑂 21 , 𝑠 31 ≥ 𝑠 21 + 𝑝 21 . For the trip of vehicle 2, the earliest starting time of the transportation operation 𝑂 22 2 does not depend on the finishing time of 𝑂 12 1 (operations in the diagonal rectangle in Fig. 2). Vehicle 2 can potentially start the transportation operation 𝑂 22 2 before the transportation operation 𝑂 12 1 , as depicted in Fig. 2 Conjunctive and disjunctive arcs to link production and transportation with maximal time lags The delivery nodes where products are unloaded have a time window constraint that defines the maximal arrival time acceptable for the customer, which is correlated to the product perishability. The time windows are related to the earliest starting time of the production on the facility. They are modelled by maximal time lags that are defined using a negative arc between one delivery node and the production node in this special case.

O 22 2 ]-∞,x] ]-∞,z]
As shown in Fig. 2, the time window ] -∞; 𝑦] on node 2, which corresponds to a delivery node, means that the maximal duration between 𝑠 11 and the arrival time at node 2 must not exceed 𝑦 units of time. Similarly, node 3 must be served within a delay of 𝑥 units of time after 𝑠 21 . One positive arc gives the minimal delay between the earliest starting time of two operations. For example, τ 1,2 models a minimal duration required between departure time at node 1 and arrival time at node 2. Normally, these so-called conjunctive arcs are valuated with the shortest path value between two nodes modelling transport, and with the duration of production between two production nodes. Disjunctive arcs are required to define the order of sub-trips in a trip. This type of time windows arises: (1) in routing problems including but not limited to the Dial-A-Ride Problem where both maximal route duration and maximal customer riding time must be taken into account; and (2) in scheduling problems.

In the Gantt chart in Fig. 1, the first job provides products for customers 1 and 2 for a production duration of 𝑝 11 = 2 time units, the earliest starting time of 𝑂 11 is equal to 0 and vehicle 1 is available at time 0. Therefore, on node 1 (Fig. 2), the time window is defined by ] -∞; 2 + 𝛿 1 ] and ] -∞; 2 + 𝛿 2 ] for node 2, where 𝛿 𝑖 is the lifespan of products for customer 𝑖. The earliest starting time of a transportation operation is the maximal value between the finishing time of the previous transportation operation assigned to the same vehicle and the finishing time of the production operation. Similar considerations make it possible to define the earliest starting time 𝑂 22 2 at time 8 with a time window ] -∞; 8 + 𝛿 3 ] for nodes 3. Because the problem being considered is a single product PTSPm, the lifespan is equal for all nodes, i.e., 𝛿 𝑖 = 𝐵, ∀𝑖, implying that all time windows of all delivery nodes of the same sub-trip are identical and equal to ] -∞; 𝑥 + 𝐵], where 𝑥 is the duration of the production linked to the transport.

Figure 3. Sub-trip with maximal time lag simplification.

The earliest arrival time of a vehicle increases within a sub-trip. Therefore, if an earliest arrival time exists for the last delivery operation of the sub-trip that meets the time window constraints, all previously scheduled operations with time window constraints will hold. In Fig. 3, sub-trip 𝐺 with three delivery nodes takes the three maximal time lags modelling the product lifespan explicitly into account. For PTSPm with a single product, which consequently defines the same lifespan for all customers, a new graph 𝐺′ can be used to model only the time window of the last delivery operation of each sub-trip.

The feature that differentiates this problem from previous ones is the combination of the scheduling decisions with the limited product lifespan and the vehicle routing decisions. These interdependent decisions lead to the possibility that the product may expire before it reaches a customer if an unprofitable scheduling solution is chosen. The objective consists in solving the scheduling and routing problem by minimising the makespan to comply with the classical objective function introduced by Geismar et al. ( 2008) and providing a semi-active solution, i.e., a left-shifted solution.

A PTSPm resolution based on a GRASP ELS

This proposal is based on a GRASP×ELS, which introduces a new splitting algorithm for the assignment and the routing problem, and a new resolution of the scheduling problem using an approach based on a disjunctive graph. The disjunctive graph is specially designed to efficiently take the perishability constraint into account using maximal time lags.

Key features for a PTSPm resolution

The key point for the PTSPm resolution is to alternate between solutions encoded by giant trips (TSP -Traveling Salesman Problem -permutation list on the 𝑛 customers), set of trips (VRP -Vehicle Routing Problem -on the 𝑛 customers) that comply with the ordered set of customers defined by the giant trip, and a flow shop resolution by considering jobs linked to the trips (PTSPm on the 𝑛 customers). The approach is a combination of three search space representations that favour partial enumeration of the whole search space (Fig. 4).

The iterative search space exploration takes advantage, first, of the indirect representation of solutions by using a giant trip. Second, a set of feasible sub-trips minimising the total transportation time is computed using a split-based approach on the sequence of customers that has been defined by the giant trip. Note that the specific local search for node routing can be applied to obtain a high-quality local routing solution.

The resolution framwork is expressed by the following key features:  a local search for routing with 2-Opt or swap within a sub-trip, and 2-Opt or swap between two sub-trips, based on classical VRP neighbourhoods;  a splitting procedure for proper coordination between production and transportation;  a scheduling procedure based on a disjunctive graph with the sub-trips made by the splitting procedure;  a concatenation procedure to obtain a new giant trip. The GRASP×ELS approach is very different from Geismar et al.'s proposal and differs in both scope and the resolution scheme that are summarised in Table 3.

The GRASP×ELS Principle

GRASP×ELS is a hybridisation [START_REF] Prins | A GRASP × Evolutionary Local Search Hybrid for the Vehicle Routing Problem[END_REF] of a GRASP (Greedy Randomised Adaptive Search Procedure) [START_REF] Feo | Greedy Randomized Adaptive Search Procedures[END_REF], with an ELS (Evolutionary Local Search) [START_REF] Wolf | Evolutionary local search for the super-peer selection problem and the p-hub median problem[END_REF][START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]) that takes advantage of both methods. The multi-start approach of the GRASP, which provides 𝑛𝑝 initial solutions, is based on a greedy randomised heuristic, and solutions are then improved by a local search procedure. The second metaheuristic can be expressed as an extension of the ILS (Iterated Local Search) [START_REF] Lourenço | Iterated local search[END_REF], referred to as ELS and proposed by [START_REF] Wolf | Evolutionary local search for the super-peer selection problem and the p-hub median problem[END_REF]. The solution space investigation is achieved by the GRASP that favours diversity, and the intensification phase is devoted to the ELS via a proper local search investigation into the local search space. In addition, to combine GRASP with ELS, another important feature is the alternation between solution spaces as stressed by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] and [START_REF] Prins | Order-first split-second methods for vehicle routing problems: A review[END_REF]. Converting a PTSPm solution into a giant trip is achieved by a concatenation procedure, and the operation converting a giant trip into a VRP solution is achieved by a dedicated splitting procedure (Split). Algorithm 1 is composed of a loop from line 13 to 36, which iterates on a new starting solution for the GRASP.

For each initial solution, the procedure Generation_of_initial_solution(), using a greedy heuristic creates a new initial solution submitted to the Split procedure (line 16) to obtain a routing solution that is improved (line 17) by the local search procedure Local_Search_on_Routing(S, nl). Secondly, the Scheduling() procedure computes a solution (line 18). If the cost of the solution 𝑆, 𝑓(𝑆) is better than the cost of the best PTSPm solution found, 𝑆*, 𝑓* (line 20), 𝑆 becomes the new best solution. The while loop from lines 23 to 35 is the ELS and encompasses the loop (lines 25 to 32) for neighbourhood generations. According to the common Split approach, the mutation operator (line 26) is defined on the giant trip T and not on the solution.

The random heuristics, referred to as Generation_of_initial_solution(), generates initial solutions. This method uses two heuristics to build the initial solution.

The first one, with a probability of 0.9 is based on a greedy randomised heuristic based on a path-scanning-like approach. The first heuristic builds sub-trips one-by-one, starting from the depot, and extends each sub-trip customer-by-customer using two criteria:  the extension step at node 𝑖 considers that the sub-trip moves to the nearest customer not yet served. If vehicle load ≥ 𝑄, the next customer to be served is selected to minimise the distance to the depot;  the next customer to be served is the customer that minimises the transportation time.

The second heuristic with a probability of 0.1 classifies the customers by decreasing transportation time. f* := ; O := Ø 13.

for p := 1 to np do // GRASP loop 14.

S := call Generation_of_initial_solution () 15.

T := call Concat (S) 16.

S := call Split (T) 17.

S := call Local_Search_on_Routing (S, nl) 18.

S := call Scheduling (S) 29.

T := call Concat (S) 20.

if (f(S) < f*) then f* := f(S); S* := S; // f: the cost of a solution 21. endif 22.

i, r := 0 23. while (i < ne) and (r < nr) do // ELS loop 24.

i := i + 1; f" :=  25.

for j := 1 to nd do // mutation loop 26.

T' := call Mutation (T) 27.

S' := call Split (T') 28.

S' := call Local_Search_on_Routing (S', nl) 29.

S' := call Scheduling (S') 30.

T' := call Concat (S') 31.

if (f(S') < f") then f" := f(S'); T" := T'; S" := S'; endif 32. endfor 33.

if The local search on the routing (Local_Search_on_Routing()) is achieved using several classical VRP neighbourhood moves to improve the initial solution, namely 2-Opt within a subtrip and insertion within a sub-trip. At each iteration, the first improved move is executed but requires verification of the sub-trip feasibility by considering the capacity of the vehicles and the lifespan. All solutions are converted into a giant trip by random concatenation of their subtrips and then evaluated by the Split procedure. The key point for the efficiency of GRASP×ELS is to alternate between solutions encoded as giant trips and PTSPm solutions.

Split procedure

Giant trips are evaluated via the Split procedure that minimises the total trip duration subject to the vehicle capacity, as reported by [START_REF] Beasley | Route-first cluster-second methods for vehicle routing[END_REF] and [START_REF] Prins | Order-first split-second methods for vehicle routing problems: A review[END_REF]. Split is a key procedure used to convert a giant trip into a VRP solution (with respect to the sequence) and it is based on the classical Split procedure adjusted to address the specific PTSPm constraints. The algorithm provides an efficient solution by building an auxiliary graph 𝐺 = (𝑋, 𝐴), where X represents n + 1 nodes numbered from 0 to n. Node 0 is a dummy node, while the nodes 1 … n correspond to the sequence of the giant trip 𝐺𝑇 = (𝜎 1 , … , 𝜎 𝑛 ). An arc (𝑖, 𝑗) belongs to A if a sub-trip serving customers 𝜎 𝑖+1 to σ j (inclusive) is both weight-feasible and lifespanfeasible. An initial label is set at node 0 and the labels are propagated from node to node in G using arcs. The best label at node 𝑛 is kept as the best solution.

The routing problem has a resource constraint since there is a limited number of available vehicles, each vehicle defining one specific machine to be scheduled. The Split procedure for the PTSPm is an extension of the split introduced for the HVRP (Heterogeneous Vehicle Routing Problem) by [START_REF] Duhamel | A hybrid evolutionary local search with depth first search split procedure for the heterogeneous vehicle routing problems[END_REF] because there is the earliest finishing time of the production to take into account as well as the earliest finishing time of vehicles for the previous sub-trip. Computation of the resource-constrained shortest path in the graph is typically achieved by a label-correcting algorithm. According to Desrocher (1988), several labels per node have to be handled and the key point consists in defining the label structure with the cost and the resource availability. The generic Split procedure of [START_REF] Duhamel | A hybrid evolutionary local search with depth first search split procedure for the heterogeneous vehicle routing problems[END_REF] is based on three key points:  a label description based on the resource availability;  a dominance rule between labels to improve the running time by keeping only promising labels on the node;  a propagation rule to define a label from node 𝑖 to node 𝑗 depending on the customer sequence (σ i+1 , … , σ j ).

Let 𝐿 𝑖 𝑝 = (𝑑 𝑃𝑖 , 𝑑 𝑉1𝑖 , … , 𝑑 𝑉𝑁𝑖 ) be the 𝑝 𝑡ℎ label assigned to node 𝑖 when the number of vehicles 𝑉 = 𝑁, where L i p (𝑗) is the earliest finishing time of vehicle 𝑗 for 𝑗 = 2. . 𝑁 + 1 and + 𝜏 𝜎 𝑗 ,0 . The max operator, which appears in the propagation rule, is explained in Fig. 5. In the first case, 𝑂 𝑗1 has a duration (𝑝 𝑗1 ) greater than the duration between the end of the production operation 𝑂 𝑖1 at time 𝑑 𝑃𝑖 and the end of the previous sub-trip assigned to vehicle 𝑉1 (𝑑 𝑉1𝑖 ). . In the other case, if 𝑂 𝑗1 has a duration (𝑝 𝑗1 ) lower than the duration between the end of the production operation 𝑂 𝑖1 at time 𝑑 𝑃𝑖 and the end of the previous sub-trip assigned to vehicle 𝑉1 (𝑑 𝑉1𝑖 ), then to be compliant with the no-wait constraint, the constraint 𝑑 𝑃𝑗 = 𝑑 𝑉1𝑖 must hold and the production planning will encompass a period of inactivity. Since a large number of labels may be generated, dominance criterion is used to improve the running time, keeping only non-dominated labels on the node thanks to the following dominance propriety.

L i p (

Dominance property

Label L = (𝑑 𝑃 , 𝑑 𝑉1 , … , 𝑑 𝑉𝑁 ) dominates label P = (𝑑 𝑃 , 𝑑 𝑉1 , … , 𝑑 𝑉𝑁 ) when the two labels are sorted and is therefore equal to (𝑑 𝜎 1 𝑝 , … , 𝑑 𝜎 𝑁+1

𝑝

) and (𝑑 𝜎 1 𝑞 , … , 𝑑 𝜎 𝑁+1 𝑞 ), respectively, and all of the following conditions hold:

 ∃𝑖 * ∈ {1, … , 𝑁 + 1}, d σ 𝑖 * 𝑝 < d σ 𝑖 * 𝑞 ;  ∀𝑗 ∈ {1, … , 𝑁 + 1}, d σ 𝑗 𝑝 ≤ d σ 𝑗 * 𝑞 .
The dominance rule limits the number of labels stored at each node but several authors, including [START_REF] Duhamel | A hybrid evolutionary local search with depth first search split procedure for the heterogeneous vehicle routing problems[END_REF], have reported that a large number of labels could still be generated. Another time-saving approach consists in limiting the maximal number 𝑁 𝑚𝑎𝑥 of labels generated during the split process and the maximal number 𝑁𝐵 𝑚𝑎𝑥 of labels stored at each node. This can reduce the CPU time, but discarding some labels may cause the algorithm to miss the optimal split. Such an approach offers a compromise between split quality and the number of labels kept on nodes and provides time-efficient sub-optimal split solutions.

The Split procedure is detailed in Algorithm 2 and uses low-level procedures to handle the label on the graph:  Propagation_label_to_vehicle(L,v,𝑖) propagates the label 𝐿 using the vehicle v at node 𝑖 and returns a new label;  CheckDomination(L, i, Ni) performs the dominance check between the label 𝐿 and all labels stored at node 𝑖, 𝑁 𝑖 . It returns 0 if 𝐿 is dominated by at least one label from 𝑖. Thus, 𝐿 must not be saved. It returns 1 if 𝐿 dominates at least one label from 𝑖. It returns 2 otherwise;  InsertLabel(L, i, CD, Ni, NBmax)attempts to add the label 𝐿 at node 𝑖 if 𝐶𝐷 ∈ {1,2} and by considering that the number of labels 𝑁 𝑖 must not exceed the maximal number 𝑁𝐵 𝑚𝑎𝑥 . If 𝐶𝐷=1, all labels from i dominated by 𝐿 are deleted. The list of labels is ordered by decreasing cost;  Extract_trip() checks the shortest path into the graph for the last node to the first node, and returns the set of trips.

The Split algorithm (Algorithm 2) is composed of two parts: the initialisation part where local variables are initialised (lines 14-15) and a for loop from lines 16 to 52, which iterates for the ordered set of customers defined by the giant trip T. The while loop from lines 19 to 51 iterates and makes it possible to evaluate the partial sequence (𝜎 𝑖+1 , … , 𝜎 𝑗 ). If 𝑖 = 𝑗, a new sub-trip is created and the initial sub-trip costs are assigned (lines 22-23) and are updated in lines 25-26.

The condition line 28 makes it possible to determine if the constraints hold and the loop is stopped early thanks to the condition on the Boolean condition (stop=true). Because a set of labels is saved at each node, the for loop of line 29 iterates over all labels of node i and addresses the label 𝐿 𝑖 𝑃 .

Algorithm 2. Split 1. procedure Split 2. input parameters 3. T: giant tour 4. output parameters 5.

S: VRP-PTSPm solution 6. global parameter 7.

Q : maximal vehicle weight capacity 1.

B : lifespan of the product 8.

qi : total items ordered by customer i 9. 𝜏ij : cost from customer i to j 10.

n : number of customers 11.

N : number of vehicles 12.

NBmax : maximal number of labels in each node 13. begin 14.

L0 1 :=(0,0,…,0), S :=  15.

for i := 1 to n do NBi := 0 endfor 16.

for j := 0 to n do 17. j := i 18.

stop := false 19.

while (j < n and stop=false) 20.

customer := Tj 21.

if (j = i) then 22.

production_cost := dcustomer/r 23.

transport_cost := cdepot,customer + ccustomer,depot 24. else 25.

production_cost += dcustomer/r 26.

transport_cost += ccustomer-1,customer +ccustomer,depot -ccustomer-1,depot 27. endif 28.

if ((production_cost*r < Q) and (transport_cost -ccustomer,depot ≤ B)) then 29.

for p := 1 to NBi do 30.

if (dPi-1 = 0) then 31.

insertion

Scheduling procedure

After the splitting procedure, each job is fully defined by an ordered sequence of customer demands that are assigned to a specific vehicle. The previous sections refer to a problem that can be defined as a two-stage hybrid flow shop with one machine at the first stage and 𝑁 machines at the second stage. At the first stage, the term 'machine' is a general term that refers to the production (production machine) and to the transportation for the second stage (also referred to as transportation machine).

A job simultaneously models the production (the first operation on the production machine) and a sub-trip (the transportation operation on one transportation machine among the 𝑁 machines available).

The problem is modelled as a disjunctive graph model first defined by [START_REF] Roy | Les problèmes d'ordonnancement avec contraintes disjonctives[END_REF] using a directed graph G = (V, A, E P , E T ), where V represents the set of nodes that contains one element for each operation 𝑂 𝑖 , a source node 0 connected to the first operation of each job, and a sink node * linked to the last operation of each job (Fig. 6). The set A represents the set of conjunctive arcs, E P the set of pairs of disjunctive arcs between the production nodes, and E T the set of pairs of disjunctive arcs between the transportation nodes. With the specific graph characteristic, the classical disjunctive graph can be adapted using maximal time lags between the starting time of the production operation 𝑂 𝑖1 and the starting time of the transportation operation 𝑂 𝑖2 , which is consistent with the perishability constraints. The value of the maximal time lag defines the larger gap between the starting time of 𝑂 𝑖1 and the starting time of 𝑂 𝑖2 , and is defined by -(𝐵 + 𝑝 𝑖1 -𝑝′ 𝑖2 ).

Conjunctive arcs are used to connect each pair of consecutive operations of the same job. Each pair of disjunctive arcs on the production connects operation 𝑂 𝑖1 to 𝑂 𝑗1 (belonging to different jobs), and has a duration 𝑝 𝑖1 . Each pair of disjunctive arcs on the transportation connects two operations 𝑂 𝑖2 , 𝑂 𝑗2 in this order belonging to different jobs that are to be processed on the same transportation machine (vehicle), and has a duration 𝑝 𝑖2 (Fig. 6). A feasible solution corresponds to an acyclic graph, and an evaluation procedure can be defined using a Bellmanlike longest path algorithm to obtain the earliest starting time of each operation, including the makespan 𝐶 𝑚𝑎𝑥 . In Fig. 6, there are two disjunctive arcs on the transportation due to jobs using the same vehicle. For vehicle 1 assigned to jobs 1 and 4, the value on the disjunctive arc is equal to p 21 . Figure 7 is a graphical representation of relations between both the earliest and finishing times of operations with the lifespan constraint. The Bierwirth vector used in the scheduling procedure gives an order on the job using the sub-trips created by the splitting procedure. An efficient local search algorithm can be defined considering the critical path including the neighbourhood of [START_REF] Van Laarhoven | Jobshop scheduling by simulated annealing[END_REF] and [START_REF] Grabowski | A block approach for single machine scheduling with release dates and due dates[END_REF] The neighbourhood of Laarhoven consists of swapping two consecutive operations assigned to the same machine along the critical path, leading to a modification of the machine disjunction. Due to the indirect representation of solutions, the operation swap is achieved in the Bierwirth vector by switching the two job numbers of the two operations in the disjunctive graph. The permutation of operations based on the block definition and located at the critical path is included in a depth first search local search investigates the move from the end of the critical path to the dummy node. If one permutation leads to a lower cost solution, the critical path exploration is restarted to the end of the graph. Note that due to computation time considerations, a maximal number of iterations ns is reached. A full example is given in detail on the Web page.

Hash function for PTSPm

To avoid premature convergence of the global algorithm and to stimulate the investigation of the part of the search space not previously investigated, it is necessary to use an efficient clone detection system. Previous research into scheduling/routing problems based on indirect representation approaches has proved that efficient clone detection must be defined considering the indirect representation or the solution. For the giant trip 𝐺𝑇, the related hash function 𝐻(𝑡) is defined as follows:

𝐻(𝑡) = ∑ 𝐺𝑇[𝑖] × 𝑖 𝑚𝑜𝑑 𝐾 𝑛 𝑖=1
. 𝐺𝑇[𝑖] is the 𝑖 𝑡ℎ customer in the giant trip and 𝐾 is a constant. The value of 𝐾 impacts the size of the map and the probability of collision (consequently, 𝐾 must be as large as possible considering the memory available).

Lower bound on the makespan

Three PTSP lower bounds were introduced by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF], leading to a poor evaluation with multiple vehicles. Considering sub-trips with only one customer associated with immediately available vehicles, 𝑁 ≥ 𝑛, a lower bound can be defined using Jackson's algorithm (1955) that considers the one-machine scheduling problem with tails. To obtain the lower bound, the proposal is to extend Jackson's earliest due date rule. The lower bound based on the Jackson rule creates an optimal schedule that minimises the maximal lateness by ordering the jobs by non-decreasing due dates because the release dates are equal to 0.

Line 18 in Algorithm 3 makes it possible to define the sum of the processing times of previously scheduled jobs. At line 19, the release date of the job is updated as the maximum between the current release date and the sum of the processing time plus the due date of the job. These two lines do not appear in Jackson's original algorithm; they have been added to the algorithm in order to match the lower bound of the PTSPm. J: set of job/sub-trip 4. global parameter 5. rj : release date of the job/sub-trip j 6. pj : processing time of the job/sub-trip j 7.

q : sum of the processing time already scheduled 8. dj : due date of the job/sub-trip j 9.

n : number of customers/sub-trip 10. output parameters 11.

Cmax : lower bound 12. begin 13.

L := {1, …, n}, q := 0 14.

if (J ≠ ) then 15.

while (L ≠ ) do 16.

u := minj{rj}, j in L 17.

choose i in L such as ri = u and di maximal 18.

q := q + pi 19.

for all k in L do rk:= max{rk, q + di}, Cmax := rk 20.

L := L\{i} 21. endwhile 22. endif 23. end

To apply Algorithm 3 on the PTSPm, a parallel should be drawn between the processing time and the production time of the PTSPm, as well as between the due date and the PTSPm transportation time. The lower bound is then given by Cmax at the end of the algorithm.

Computational evaluation for the PTSP

The benchmark introduced by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] is composed of three small-scale instances with 40 customers and three medium-scale instances with 50 customers, and customers' demands are uniformly distributed between 100 and 300. The locations of the customers in the three instances with 40 or 50 customers are randomly generated in the square of 200, 300 and 400 side lengths, which denotes 12 sets of parameters (𝑆𝑃 is detailed in Table 5a). As a consequence, there are six datasets (𝐷𝑆 is is detailed in Table 5b). For each instance, the value of three parameters is set: 𝑟 ∈ {1,2,3} is the production rate; 𝑄 ∈ {300,600} is the truck capacity; and 𝐵 ∈ {300,600} is the product lifespan. We therefore have 72 resulting instances in the benchmark. Note that the distance is the transportation time between customers by considering the Euclidean distance and a vehicle speed of 1 unit of distance per unit of time.

Table 4 proposes speed factors that are used in the following comparative tables to ensure a fair comparative study and that were established by previous research articles including [START_REF] Dongarra | Performance of various computers using standard linear equations software[END_REF], and on http://www.roylongbottom.org.uk/linpackresults.htm. Since the MIPS performance is not the only influence on the CPU time, Table 4 also provides the information available about the computer, the operating system and the language. 5a gives the definition of the 12 sets of parameters 𝑆𝑃 = (𝑟, 𝑄, 𝐵), and Table 5b gives the six dataset definition 𝐷𝑆 = (𝐶, 𝑆). In Table 6, each value is an average value over the 12 sets of parameters (𝑆𝑃 = ⟦1,12⟧), for the six datasets (𝐷𝑆 = ⟦1.6⟧) and over the five replications for our proposal. Note that in Table 6, Table 7 andTable 8, the values reported by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF], [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF], and for our proposal are rounded off to the nearest integer value to ensure a fair comparative study. The genetic algorithm of [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] has an average time of 169 seconds, while the average computing time of our method is 141 seconds (average time scale of Table 6). The average computing time of the branch-and-cut of [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF] is significantly larger, with approximately 3249 seconds. ). The analysis of the full solution set (available on the Web page) makes it possible to confirm that the GRASP×ELS provides better solutions in 31 instances, equal solutions in 35 instances and worse solutions in only six instances. In general, when the solutions are improved, they are improved by approximately 4.7%, whereas when they are worsened, they are worsened by approximately 0.04%. The efficiency of the GRASP×ELS algorithm is more significant for the instances in which 𝑄 = 600 (𝑆𝑃 from 7 to 12) particularly for instances with the sets of parameters 𝑆𝑃 = 9 and 𝑆𝑃 = 12 that have a high production rate (Table 7). For instances from 𝑆𝑃 = 7 to 𝑆𝑃 = 12, the transportation time dominates the production time, as underlined by [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF]. Therefore, improvements on routing phase lead to better results for these instances. In Table 7, with 𝑆𝑃 = 12, the gap has been significantly reduced from 45.2% [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] to 39.3% for our proposal. With the method used by [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF], the gap has also been reduced to 39.2%. For the PTSP, our method provides 18 optimal solutions vs. 16 optimal solutions for the dedicated method of [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF], as illustrated in Table 7. The method also outperforms the method of [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF], providing 18 solutions that reached the lower bound vs. 13 solutions for [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF].

Computational evaluation for the PTSPm

Benchmarks for the PTSPm

The instances of [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] are extended to tackle up to six vehicles and to evaluate the impact of routing with regard to the scheduling solutions (Table 8). It can be observed that by increasing the number of vehicles from one to two, the average makespan for all the instances is reduced from 25% (decreasing from 7,884 to 5,820) and then decreased by 6% (from 5,820 to 5,478) when increasing the number of vehicles from two to three. Table 8. Average results with the GRASP×ELS for the PTSPm.

After intensive numerical experiments, we found that a fleet of six vehicles is sufficient to reach the lower bound of each instance with the strongest constraints, 𝑟 ∈ {1,2,3}, 𝑄 = 300 and 𝐵 = 300. The instances of [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] and the details of our solutions are available on the Web page. A careful analysis of the ratio between the average customer demand and the vehicle capacity remains constant (approximately 0.3) for all instances. Consequently, the average number of customers per sub-trip remains quite low (approximately two on average).

Proposal of new large-scale instances with multiple vehicles

A new set of instances with a broad range of parameters is defined considering:  customer distribution to model both urban areas and rural surroundings;  random depot node location (centered, peripheral location, etc.);  large vehicle capacities and, consequently, larger sub-trips. Each instance is characterised by:  the number of customers that follows a discrete uniform distribution 𝑈(50,100) for the large instances or 𝑈(100,200) for the very large instances;  an interest area that contains all the customers and the plant. The length X and the width Y of this area can vary between 100 and 500, which represents 25 possibilities (100 × 100, 100 × 200, …, 500 × 400 and 500 × 500). This area is centered in (0,0) and divided into four sectors [-𝑋, 0] × [0, 𝑌], [0, 𝑋] × [0, 𝑌], [0, 𝑋] × [-𝑌, 0] and [-𝑋, 0] × [-𝑌, 0];  the number of customer centers 𝑐 = {1,2,3}.

Moreover, the following parameters are necessary to define an instance:  the capacity of the trucks is equal to 1000;  the demand of each customer 𝑖, 𝑑 𝑖 = {50, 100, 200, 300, 500} with the associated probabilities {0.1, 0.2, 0.4, 0.2, 0.1};  the position of each customer is uniformly distributed in one customer center;  the position of each customer center is randomly selected among the sector centers with a variation of ±10%;  the size of each customer centre is equal to the minimal distance between the customer centre and the edge of the interest area with a variation of ±20%;  the position of the depot has a probability of 0.7 to be in the centre of one sector centres with a variation of ±10%, and 0.3 to be centred in the interest area with a variation of ±10%.

With all these parameters, we have 150 instances. An instance is caracterised by a triplet (𝑥, 𝑧, 𝑐) with: 𝑥 ∈ {𝐿, 𝑉𝐿} to indicate if the size of the instance is large 𝐿 or very lage 𝑉𝐿, 𝑧 ∈ ⟦1,25⟧ to define the size of the interest area, and 𝑐 = {1,2,3} to give the number of customer centres.

Computational results on the new instances

Five replications are achieved for each instance and Tables 9 gives the average of the best solutions found over the five replications and on the three customers centers (column ℎ * , * ) with the average total CPU time (column 𝑡𝑡 * , * ) and the average time required to find the best solution (column 𝑡 * , * ). Note that the lower bound (𝐿𝐵 * column) refers to the lower bound introduced in Section 3.6. The parameter setting remains identical for all instances and was obtained after preliminary experiments:  np/ne Number of GRASP/ELS/ neighbourhood iterations 50/10/5  lr/ns Number of local searches on the routing/scheduling 40/500  NBmax Maximal number of labels per node 5

The GRASP×ELS on the (𝐿, * , * ) instances has an average time ranging from 1.13 to 0.1 seconds, depending on the number of vehicles, whereas the average computing time on the (𝑉𝐿, * , * ) instances is equal to 6.63 seconds for one vehicle and 0.86 seconds for two vehicles (Tables 9a and9b). The computing time has thus been multiplied by six between the large and the very large instances. Moreover, on average, the gap between our results and the lower bound is very small, with 0.09% on large instances and 0.03% on very large instances for one vehicle. With two vehicles, the gap between our results and the lower bound is even smaller, with 0.005% on large instances and 0.003% on very large instances

The analysis of trips makes it possible to prove that the trip design is linked to the sub-graph defined by the densities of nodes, referred to as customer centres. All the trips and a graphical representation for instance 1 are available on the Web site.

Concluding remarks

Proper integration of production planning and routing is the key feature in a supply chain since the coordination of these two functions has a significant impact on the customer service level. This paper addresses the PTSP with multiple vehicles (PTSPm) in order to extend the PTSP with a single vehicle in [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] and can be efficiently used to compute coordinated solution in supply chain. This approach proves it is possible to solve the two problems in a coordinate way and permits to obtain better solution than classical approach where the two problems are solved sequentially.

The framework we propose takes advantage of an indirect representation of the solutions using a split-based approach with search space alternation between TSP solutions, VRP solutions and PTSPm solutions. This indirect approach is one of the key features of the proposal. The framework efficiency leads to a special disjunctive graph for the trips and two kinds of disjunctive arcs due to a single production facility and the vehicles. The method has proven to be efficient in the one-vehicle instances, providing better solutions than [START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] and [START_REF] Karaoğlan | The coordinated production and transportation scheduling problem with a time-sensitive product: a branch-and-cut algorithm[END_REF], with shorter computation times. The method also generates new solutions for the multi-vehicle extension. In order to ensure fair comparative studies, a new set of instances has been introduced.

Our research is now directed towards a bi-objective resolution where a second criterion could be introduced for the quality of service, which could be defined as the delay between the arrival time of one vehicle at a customer node and the upper bound of the time window. This second criterion should be relevant for the quality of service, provided that there is a large enough difference between the delivery date and the expiration date for customer products modelled by the perishability constraint.

Figure 2 .

 2 Figure 2. Example of a two-vehicle workload.

Figure 4 .

 4 Figure 4. The three search spaces to build a PTSPm solution.

  The finishing time of 𝑂 𝑖1 therefore corresponds to the starting time of 𝑂 𝑗1 and, consequently, 𝑑 𝑃𝑗 = 𝑑 𝑃𝑖 +

Figure 5 .

 5 Figure 5. The two cases of label propagation in the Split procedure.

Figure 6 .

 6 Figure 6. Example of the scheduling procedure with five sub-trips and V=3.

Figure 7 .

 7 Figure 7. Representation of the constraint due to the disjunctive arcs of the transportation.

Algorithm 3 :

 3 A Jackson rule-based lower bound 1. procedure A lower bound based on the Jackson Rule 2. input parameters 3.

Table 2 .

 2 Example of solution for the two-vehicle PTSPm.

		𝑂 11	0 1	𝑂 12 1	2	0	𝑂 21	0	𝑂 22 2	3	0	𝑂 31	0	4	𝑂 32 1	5	0
	Starting time	0					2					8
	Time Windows			]-∞;9]				]-∞;15]					]-∞;19]
	Departure time		2 5		7			8	10			12 15		18
	Arrival time		5		7 10			10 12			15		18 19

  and Table 2. Assignment of

	vehicles to sub-trips is a challenging problem that should be solved avoiding extra waiting time
	on the Gantt chart.
	Vehicle							O 12 1					O 22 2						O 32 1
	loading						Node 1		Node 2		Node 3		Node 4	Node 5
		0				s 12				]-∞,y]	s 22	f 12				s 32 f 22						f 32	Time
		s 11			s 21						s 31
	Prod.	O 11	p 11	O 21			p 21			O 31				p 31
	Trip 1	0	S	0 p 11	P	0	τ 0,1 p 21	D	τ 1,2	-y	D	τ 2,0	0	p 31		0	P	τ 0,4	D	τ 4,5 -z	D	τ 5,0	0	0	*
	∀𝑏 ∈ 𝐽 Trip 2			𝑠 𝑏1 + 𝑝 𝑏1 ≤ 𝑠 𝑏2 P τ 0,3 No time dependency between 0	-x	D	τ 3,0	0			0		(20) End of trip 1
					and the starting time of the finishing time of O 12 1							End of trip 2

Table 3 .

 3 GRASP×ELS approach vs. Geismar et al.'s approach

			(Geismar et al., 2008)	GRASP×ELS
	Scope		PTSP	PTSP + PTSPm
		Scheduling	Gilmore-Gomory	Disjunctive graph
		Routing	Split of the classical VRP	Split with production and transport consideration
	Resolution scheme	Local search for scheduling	/	Based on the proposal of Laarhoven et al. (1992) and Grabowski et al. (1986)
		Local search for routing	2-Opt moves	Inter/extra trips 2-Opt and swap moves
		Solution representation	Chromosome	Giant trip

  with the introduction of blocks.

			B+p 11
			B
	P	O 11	O 14		lifespan constraint
	V1		O 21		O 24
	V2
	V3
	s 11	f 11	s 21	s 21 +p' 21	fs 21

Table 4 .

 4 Comparative performance of processors.Five replications (𝑟𝑒𝑝) are processed for each instance and Tables 6, 7 and 8 give the average of the best solutions found over the five replications with the average CPU time required for each run. The parameters used in the GRASP×ELS are given below. They remain unchanged for all instances and were obtained after preliminary experiments:

		(Geismar et al., 2008)	(Karaoğlan et al., 2017)	Our framework
	Computer	Intel Pentium 4 3.2GHz	Intel Xeon 3.16GHz	Intel Core i7 3.4GHz
	OS	Windows XP Pro.	Windows 7	Windows 7
	Language	BASIC	C++ IBM ILOG CPLEX 12.6	Visual C++
	MFlops	1573	1892	2671
	Speed factor	1.0	1.2	1.7

Table

Table 5 .

 5 Instances characteristics a. Definition of the 12 sets of parameters (SP).

					b. Definition of the 6 datasets (DS).
	SP	r	Q	B	DS	C	S
	1	1	300	300	1	40 200
	2	2	300	300	2	40 300
	3	3	300	300	3	40 400
	4	1	300	600	4	50 200
	5	2	300	600	5	50 300
	6	3	300	600	6	50 400
	7	1	600	300
	8	2	600	300
	9	3	600	300
	10	1	600	600
	11	2	600	600
	12	3	600	600

Table 6 .

 6 Average results on the instances of[START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF] with five replications.

		(Geismar et	(Karaoğlan et		Our proposal
	* 𝐿𝐵 *	al., 2008) ℎ * , * * 𝑡𝑡 * , * *	𝑡 * , * *	al., 2017) ℎ * , * * 𝑡𝑡 * , * *	𝑡 * , * *	ℎ * , * *	𝑡𝑡 * , * *	𝑡 * , * *	ℎ * *	𝑡𝑡 * *	𝑡 * *
	Avg. 5383	8045		7891		7884			7876
	Avg.time	169 N/A	2708 N/A		82.7 40.7		81.1 42.7
	Avg. time scale	169 N/A	3249 N/A		141	69		138	73
	Nb. Opt	16/72		13/72		16/72			18/72
	Our method provides an average value ℎ * , * * of about 7884, which is better than Geismar et al.
	(2008) (average value of approximately 8045) and Karaoğlan et al. (2017) (average value of
	approximately 7891

Table 7 .

 7 Average results on the instances of[START_REF] Geismar | The Integrated Production and Transportation Scheduling problem for a Product with a Short Lifespan[END_REF].

	k ∈ SP	𝐿𝐵 * 𝑘

(Geismar et al., 2008) (Karaoğlan et al., 2017) Our proposal

		ℎ * , * 𝑘		%Gap	ℎ * , * 𝑘	%Gap	ℎ * , * 𝑘	%Gap
	1	8781	10040	12.5	10049	12.6	10039	12.5
	2	4412	9157	51.8	9179	51.9	9148	51.7
	3	2956	9125	67.6	9153	67.7	9118	67.5
	4	8781	10041	12.5	10048	12.6	10026	12.4
	5	4412	9171	51.9	9180	51.9	9151	51.7
	6	2956	9153	67.7	9152	67.7	9122	67.6
	7	8781	8781	0.0	8782	0.0	8782	0.0
	8	4412	5744	23.1	5346	17.5	5347	17.5
	9	2956	5421	45.4	4887	39.5	4919	39.9
	10	8781	8781	0.0	8782	0.0	8781	0.0
	11	4412	5724	22.9	5276	16.4	5300	16.7
	12	2956	5403	45.2	4861	39.2	4875	39.3
	Avg.	5383	8045	33.4	7891	31.4	7884	31.4
	Nb. Opt		16/72		13/72		16/72

Table 9 .

 9 Results on the large and very large-scale instances with one or two vehicle. a. Large instances (𝑳, * , * )

				1 vehicle			2 vehicles
	Inst.	LB * ̅̅̅̅̅	ℎ * , *	𝑡𝑡 * , *	𝑡 * , *	ℎ * , *	𝑡𝑡 * , *	𝑡 * , *

Corresponding author. e-mail addresses: placomme@isima.fr (P. Lacomme); aziz.moukrim@utc.fr (A. Moukrim); alain.quilliot@isima.fr (A. Quilliot); marina.vinot@isima.fr (M. Vinot)

Acknowledgements: This work was carried out and funded within the framework of the ATHENA project (Ref.: ANR-13-BS02-0006) and Labex MS2T. It was supported by the French government through the "Investments for the Future" program managed by the French National Research Agency (Ref.: ANR-11-IDEX-0004-02).