Zakaria Chekoubi
email: zakaria.chekoubi@univ-lorraine.fr

Wajdi Trabelsi
email: wajdi.trabelsi@icn-artem.com

Nathalie Sauer
email: nathalie.sauer@univ-lorraine.fr

Keywords: Integrated Production-Inventory-Routing Problem, Lot-sizing, Remanufacturing, Inventory control, VRP with Simultaneous Pick-up and Delivery, Reverse logistics, Closed-Loop Supply Chain, Mixed Integer Linear Programming. I

Production-planning, inventory control with returns management, and vehicle routing are critical and highly related problems in the design of closed-loop supply chain networks. In this paper, we address an extension of the Integrated Production-Inventory-Routing Problem (IPIRP) since we consider recovery and remanufacturing of End-Of-Life (EOL) products as a sustainable aspect. A mathematical programming formulation is proposed for the considered problem. The objective is to determine both the optimal amounts of goods (lot-sizes) to be manufactured, remanufactured and stored while meeting customers' requests (deliveries and pick-ups) with a minimum total cost due to the integrated operations considered. The originality of our study lies in the fact of jointly treating the capacited dynamic lot-sizing with remanufacturing problem and the vehicle routing problem with simultaneous pick-up and delivery. Some numerical examples to illustrate the relevance of our proposed model are given.

The Integrated Production-Inventory-Routing problem in the context of Reverse Logistics: the case of collecting and remanufacturing of EOL products INTRODUCTION Managing reverse flow within a close-loop supply chain is considered as a business opportunity for many companies. Today, many firms have initiated efforts to incorporate the concept of reverse logistics in their regular production, inventory and distribution decision systems for several reasons such as the increasing concern for environment, restrictions of applicable government regulations and laws on recycled products and waste disposal, the growing energy consumptions, and the stiff competition between the firms. The integration of these three decisions on the context of reverse supply chains offers tremendous cost savaging opportunities to firms. Moreover, integrating production-inventory-distribution decisions into a single problem is almost indispensable and relevant for some type of goods in particular perishable or timesensitive goods. The Integrated Production-Inventory-Routing Problem (IPIRP) is considered as a new OR field which aims to jointly optimize these three decision levels. It can be described as follows: a plant produces and stocks single or multiple products and distributes them to multiple customers (even by direct shipment or by a special customer sequence) over a multi-period time horizon by jointly considering all the three decisions. Most existing studies on this problem consider the strategic or tactical level [1]. At the strategic level, decisions about facility location and plant capacity are taken [START_REF] Moons | Integrating production scheduling and vehicle routing decisions at the operational decision level: a review and discussion[END_REF]. The tactical level deals with production lot sizes, inventory levels, delivery quantities and routing paths. A recent review of integrated tactical level problems can be found in [START_REF] Díaz-Madroñero | A review of tactical optimization models for integrated production and transport routing planning decisions[END_REF]. The operational level deals with the integration of machine scheduling and distribution decisions. A recent review on the integrated operational level can be found in [START_REF] Moons | Integrating production scheduling and vehicle routing decisions at the operational decision level: a review and discussion[END_REF]. Therefore, the IPIRP can be seen as an integration of three-level lot-sizing, inventory control and routing problems, each of them is quite difficult to solve. Because the IPIRP is NP-hard, there are only few exact algorithms available in the literature ([1], [START_REF] Adulyasak | Formulations and branch-and-cut algorithms for multivehicle production and inventory routing problems[END_REF] and [START_REF] Archetti | Analysis of the maximum level policy in a production-distribution system[END_REF]) and these algorithms can only solve small-scale instances. Moreover, many heuristics have been proposed for several variants of the IPIRP ([6], [START_REF] Bard | Heuristics for a multiperiod inventory routing problem with production decisions[END_REF], [START_REF] Bard | A branch-and-price algorithm for an integrated production and inventory routing problem[END_REF]). For a recent review on the IPIRP, readers can refer to [START_REF] Adulyasak | The production routing problem: A review of formulations and solution algorithms[END_REF], [START_REF] Coelho | Thirty years of inventory routing[END_REF] and [START_REF] Bard | The integrated productioninventory-distribution-routing problem[END_REF]. Furthermore, studying the IPIRP in the context of reverse logistics makes it more challenging issue especially at the tactical level. Our study considers this relevant application area by considering an extension of the IPIRP in which we jointly consider the wellknown Capacited Dynamic Lot-Sizing with Remanufacturing (CDLSR), at production-inventory decision level, and Vehicle Routing Problem with Simultaneous Pickup and Delivery, at routing decision level. We called resulting integrated problem as CDLSR-VRPSPD problem. The first part of the integrated CDLSR-VRPSPD problem deals with the single lot-sizing problem with remanufacturing and inventory control system. The concept of remanufacturing consists on transforming endof-life returned products (reverse flows) into usable products through refurbishment, repair or upgrading. Today many firms voluntary collect end-of-life products to recover residual value and might even design their products to maximize this value [START_REF] Decroix | Optimal policy for a multiechelon inventory system with remanufacturing[END_REF]. The remanufactured products have usually the same quality as the new products (serviceables) and are sold for the same price, but, they are almost loss costly. Also, these products can also be used to fulfil customers' requests. This is why remanufacturing plays, nowadays, a crucial role in the overall production, inventory and distribution systems. Furthermore, one of the significant problems that arises in such kind of production systems is the dynamic lot sizing problem. In [START_REF] Wagner | Dynamic version of the economic lot size model[END_REF],

the authors addressed the dynamic lot sizing problem for the first time. The so-called Wagner-Whitin problem assumes a finite planning horizon; demand is given for each interval and may vary over time; capacity constraints are not considered and a single item of a single-level product is produced. The problem is to determine the number of items to be produced in each time period in order to satisfy the demands while minimizing the total production, setup and inventory costs. Moreover, in the recent literature, there has been an increasing interest in the study of lot-sizing problems with remanufacturing due to their large applicability in reverse logistics. Two recent reviews on the modelling of lot-sizing and reverse logistics inventory systems are given by [START_REF] Aloulou | A bibliography of non-deterministic lot-sizing models[END_REF] and [START_REF] Bazan | A review of mathematical inventory models for reverse logistics and the future of its modeling: An environmental perspective[END_REF]. Mainly, the single-item dynamic lotsizing problem with remanufacturing (SDLSR) has attracted the attention of several authors during the last years. One of the most important extensions of the SDLSR problem, is the Capacited Dynamic Lot-Sizing with Remanufacturing (CDLSR) problem in which limited production capacity is considered. Several extensions of the CDLSR problem can be found in [START_REF] Masmoudi | Multiitem capacitated lot-sizing problem in a flow-shop system with energy consideration[END_REF] and [START_REF] Zhang | A Lagrangian relaxation based approach for the capacitated lot sizing problem in closed-loop supply chain[END_REF]. Recent, works on heuristics approaches can be found on [START_REF] Piñeyro | The economic lot-sizing problem with remanufacturing: analysis and an improved algorithm[END_REF], [START_REF] Schulz | A new Silver-Meal based heuristic for the single-item dynamic lot sizing problem with returns and remanufacturing[END_REF] and [START_REF] Sifaleras | Variable neighborhood search for the economic lot sizing problem with product returns and recovery[END_REF]. For recent review of mixed integer programming techniques, readers can refer to [START_REF] Cunha | A computational comparison of formulations for the economic lot-sizing with remanufacturing[END_REF] and [START_REF] Retel Helmrich | Economic lot-sizing with remanufacturing: complexity and efficient formulations[END_REF].

The second part of the integrated CDLSR-VRPSPD problem considers the VRP with simultaneous pickup and delivery which is one of the most known variants of the classical VRP. In VRPSPD, the vehicles (homogeneous or heterogenous fleet) are not only required to deliver goods to customer, but also pickup goods at customer locations [START_REF] Karaoglan | Vehicle Routing Problem with Simultaneous Pickup and Delivery: Mixed Integer Programming Formulations and Comparative Analyses[END_REF]. A general assumption in VRPSPD is that all delivered goods must be originated from the depot, all pickup goods must be transported back to the depot. Delivery and pickup goods must be met simultaneously when each customer is visited only once by a vehicle and unloading is carried out before loading at the customers without violating vehicle capacity constraints [START_REF] Chen | Vehicle routing problem with simultaneous deliveries and pickups[END_REF]. Also, the VRPSPD is very used and more relevant in reverse logistics application area. Thus, after VRPSPD has been introduced by [START_REF] Min | The multiple vehicle routing problem with simultaneous delivery and pick-up points[END_REF], the research on this problem has been mainly focused on heuristic and metaheuristic approaches ( [START_REF] Polat | A perturbation based variable neighborhood search heuristic for solving the vehicle routing problem with simultaneous pickup and delivery with time limit[END_REF], [START_REF] Tchoupo | Ant Colony Optimization Algorithm for Pickup and Delivery Problem with Time Windows[END_REF], [START_REF] Kalayci | An ant colony system empowered variable neighborhood search algorithm for the vehicle routing problem with simultaneous pickup and delivery[END_REF], and [START_REF] Euchi | Genetic scatter search algorithm to solve the one-commodity pickup and delivery vehicle routing problem[END_REF]). An exact method, based on branch and price approach, has been developed by [START_REF] Dell'amico | A branch-and-price approach to the vehicle routing problem with simultaneous distribution and collection[END_REF]. Interested readers can refer to [START_REF] Berbeglia | Static pickup and delivery problems: a classification scheme and survey[END_REF] and [START_REF] Parragh | A survey on pickup and delivery problems[END_REF] for further investigated VRPSPD extensions and reviews.

Despite the increasing of studies describing a variety of variants and models for the IPIRP, few studies consider this problem in the context of reverse logistics. Moreover, to our best knowledge, there is no study considering jointly the CDLSR and VRPSPD problems. In this paper, a particular closed-loop supply chain network integrating production planning, inventory control and distribution decisions is considered. The remainder of the paper is organized as follows. In Section 2, we provide a brief description of the considered problem as well as the assumptions adopted to develop our mathematical model. In section 3, we formally define the problem and present a MILP formulation. Computational experiments and the obtained results are described in section 4. Finally, we conclude our paper with some perspectives for future research work.

II. PROBLEM DESCRIPTION

This section provides a formal description of the CDLSR-VRPSPD problem and then presents the assumptions we considered to develop the MILP formulation model of the considered problem.

Let G = (N, A) be a complete directed graph network with = {0} ∪ {1, 2, … , } = {0} ∪ the set of nodes and = {( , ): , ∈ , ≠ } the set of arcs. The node {0} represents the central depot and the others represent the costumers. Each arc ( , ) ∈ has a non-negative symmetric routing cost that represents the cost for reaching node j from node i and satisfying triangular inequality (i.e., + ≥ ). Each customer ∈ has a periodic demand including delivery and pickup (we assume = = 0 ∀ ∈ ). We consider a single vehicle with capacity (we assume that > ∈ {∑ ∈ ; ∑ ∈ } ) and a fixed cost . At the central depot are located two production systems, respectively, for manufacturing and remanufacturing of new products (serviceables) and return-products. The production strategy of both systems is based on the economic lot-sizing with capacity of a single product for which the demand requirements of each period over the finite planning horizon must be satisfied either by manufacturing new products or by remanufacturing end-oflife products returned to the returns inventory. The studied closed-loop supply chain can be modelled as shown in Fig. 1.

The objective of the CDLSR-VRSPD problem is to determine the optimal integrated planning that minimizes the sum of setup, manufacturing, remanufacturing, holding and transportation costs by determining the optimal amounts of goods to be manufactured, remanufactured, stored and to be distributed through optimal routes. The following assumptions are considered:

FIG. 1: CONSIDERED CLOSED-LOOP SUPPLY CHAIN SYSTEM

Production level: we consider the mono-item dynamic lot sizing policy for both manufacturing and remanufacturing production systems. Each production line has a production capacity should not to be exceeded in each period. Furthermore, each production system has its own set-up cost. In addition, a unitary production cost occurs for manufacturing (resp., remanufacturing) a new product (resp., returned product) in each period. We also assume that manufacturing and remanufacturing lead times are zero. Moreover, we consider that all demands (delivery and pickup) are known for all periods of planning horizon. Also, the remanufactured products will be of the same quality as new products and backlogging is not allowed.

Inventory level:

We assume that the initial stocks of serviceables and returns are both predefined (zero initial stocks conditions) and that there us a positive demand in the first period. Moreover, we assume that new manufactured (resp., remanufactured) products are stored in the serviceables (resp., returns inventory) without exceeding its periodic capacity. We consider also that demand is meet either from newly manufactured products or from the remanufactured returns or both. We also assume that each inventory has its own periodic holding cost.

Vehicle routing level: we assume that the vehicle route starts and ends at the central depot, and every customer is visited at most once by the vehicle in each time period and each customer with zero requests should not be visited ( > 0 and ≤ ). Also, the vehicle should perform pickup and delivery simultaneously. Furthermore, the vehicle total load should not exceed its capacity.

III. MATHEMATICAL FORMULATION

In this section, we provide a mixed integer linear program (MILP) formulation for the CDLSR-VRPSPD problem which we used to solve some small instances.

A. Parameters

To formulate the problem, following parameters are used:

: set of all nodes (node 0 represents production facilities and central depot) = \{0}: set of nodes (customers) without node 0 : number of time periods (horizon planning), ∈ {1, … , } : delivery demand of customer C ∈ i N at period : pick-up demand of customer C ∈ i N at period ℎ : unitary holding cost of one serviceable product at period ℎ : unitary holding cost of one return-product at period : set-up cost for manufacturing at period : set-up cost for remanufacturing at period : unitary production cost for manufacturing at period : unitary production cost for remanufacturing at period : maximal production capacity of manufacturing system at period : maximal production capacity of remanufacturing system at : maximum storage capacity of serviceable products at : maximum storage capacity of return-products at : vehicle capacity; : fixed vehicle cost; : transportation cost between node i and node j (assume = and = 0 ∀( , ) ∈ )

B. Variables

Used variables in this model are as follows:

: number of products manufactured in period

:
= + + -∑ ∈ ∀t ∈ T (2) = -+ ∑ ∈ ∀t ∈ T (3) = ; = (4) 0 ≤ ≤ ∀t ∈ T (5) 0 ≤ ≤ ∀t ∈ T (6) ≤ min ∑ ∑ ∈ ; . ∀t ∈ T (7) ≤ min ∑ ∑ ∈ ; ∑ ∑ ∈ ; . ∀t ∈ T (8) ∑ y ∈ ≤ 1 ∀j ∈ N ; ∀t ∈ T (9) ∑ ∈ -∑ ∈ = 0 ∀j ∈N ; ∀t ∈ T (10) ∑ ∈ ≤ 1 ∀t ∈ T (11) + ≤ 1 (∀i, j ∈ N: i ≠ j); ∀t ∈ T (12) ∑ ∈ -∑ ∈ = ∀j ∈ N ; ∀t ∈ T (13) ∑ ∈ -∑ ∈ = ∀j ∈ N ; ∀t ∈ T (14)
≤ ∀t ∈ T ( 15) 16) 17) 18)

∑ ∈ = ∑ ∈ ∀t ∈ T (
∑ ∈ = 0 ∀t ∈ T (
∑ ∈ = ∑ ∈ ∀t ∈ T (
∑ ∈ = 0 ∀t ∈ T (19) + ≤ . ∀(i,j : i≠j) ∈ N ; ∀t ∈ T (20)
≥ . ∀(i,j : i≠j) ∈ N ; ∀t ∈ T ( 21)

≤ ( -). ∀(i,j : i≠j) ∈ N ; ∀t ∈ T (22)
≥ . ∀(i,j : i≠j) ∈ N ; ∀t ∈ T ( 23)

≤ ( -). ∀(i,j : i≠j) ∈ N ; ∀t ∈ T (24) ∈ {0, 1} ∀(i,j : i≠j) ∈ N ; ∀t ∈ T (25) ∈ {0, 1} ; ∈ {0, 1} ∀t ∈ T (26) , , , , , ∈ ℤ ∀(i,j : i≠j) ∈ N ; ∀t ∈ T (27)
The objective function (1) minimizes the sum of setup, manufacturing, remanufacturing, holding and transportation costs. Constraints (2)-( 4) are flow balance equations and initial inventories conditions for serviceable and return products respectively. Constraints ( 5)-( 6) define, respectively, the maximum storage capacity of serviceable and return products not to be exceeded during each period. Constraints ( 7)- [START_REF] Bard | A branch-and-price algorithm for an integrated production and inventory routing problem[END_REF] ensure that the quantities of manufactured (resp. remanufactured) products do not exceed manufacturing definite limits. Constraints ( 9)-( 12) are vehicle routing constraints. Constraint [START_REF] Adulyasak | The production routing problem: A review of formulations and solution algorithms[END_REF] ensures that the vehicle visits at most a customer per period. Constraint [START_REF] Coelho | Thirty years of inventory routing[END_REF] represents the flow conservation, while the constraint [START_REF] Bard | The integrated productioninventory-distribution-routing problem[END_REF] guarantees that one vehicle must be used. Constraints [START_REF] Decroix | Optimal policy for a multiechelon inventory system with remanufacturing[END_REF] ensures that any feasible route does not contain sub-tours with only two customers. Constraint ( 13)-( 14) guarantee that the delivery and pickup demands at customer j are satisfied. Constraint [START_REF] Bazan | A review of mathematical inventory models for reverse logistics and the future of its modeling: An environmental perspective[END_REF] ensures that the quantity of remanufactured products at the actual period was held in the returns inventory in the previous period. This constraint guarantees the policy 'remanufacturing after collecting'. Constraint [START_REF] Masmoudi | Multiitem capacitated lot-sizing problem in a flow-shop system with energy consideration[END_REF] ensures that the total pickup load on the last arcs is equal to the total pickup demand, while constraint [START_REF] Zhang | A Lagrangian relaxation based approach for the capacitated lot sizing problem in closed-loop supply chain[END_REF] guaranties that the load on the first arcs is zero. Similarly, constraints [START_REF] Piñeyro | The economic lot-sizing problem with remanufacturing: analysis and an improved algorithm[END_REF] and (19) determine the delivery loads on the first and the last arcs. Constraint [START_REF] Sifaleras | Variable neighborhood search for the economic lot sizing problem with product returns and recovery[END_REF] guaranties that the total load on the vehicle does not exceed its capacity. Constraints ( 21)-( 24) are called as bounding constraints and restrict the upper and lower bounds of the quantities to be delivered and to be collected on all arcs. Finally, constraints ( 25)-( 27) define the nature of the decision variables considered.

IV. COMPUTATIONAL RESULTS

A. Instance generation

We implemented our model and tested all the generated instances using the commercial solver IBM ILOG CPLEX 12.7 on a PC with CPU Core i5 (2.40 GHz) and 4 GB RAM. Moreover, we used the default configuration of the considered solver for all performed tests. The program execution is stopped even when the optimum is found or after a time limit of 3600 s.

The data sets were generated based on some adapted data sets from the literature of lot-sizing with remanufacturing problems and those from the integrated production routing problems in particular research works of [START_REF] Sifaleras | Variable neighborhood search for the economic lot sizing problem with product returns and recovery[END_REF] and [33]. All numerical instances were generated according to the parameters values shown in Table 1. The numerical tests were done on a set of problems that were generated by varying the number of time periods (|T|), the delivery demands ( ) and the pickup demands ( ). We obtain a total of 10*3=30 problems (combinations), and for each problem five instances are randomly generated by drawing the delivery ( ) and the pickup demands ( ) from uniform distribution. Also, we change the travelling cost ( ) from one problem to another one by redrawing new coordinates and calculating the new cost matrix for the new problem. A total of 10*3*5=150 instances were generated and tested by the MILP formulation presented in the previous section.

B. Numerical results and analysis:

To test the efficiency and the limits of our MILP formulation, we analysed the obtained results from CPLEX solver in term of the average CPU time in regards to problems sizes. The summary of computational results is given in Table 2. In this table columns Ave. Tot. Cost, Ave. CPU (s) and #OPT show the average CPU time in seconds, the average of the total cost and the number of instances solved to optimality, respectively. Based on the computational results given in Table 2 and the graphic of Fig. 2, it is easy to observe that the proposed MILP formulation is very sensitive to the number of periods (| |) and number of customers (| |). Indeed, as shown in Fig. 2, when the number of customers rises from 30 to 50, the average CPU time shows a sharp increase for all possible time period values, especially, when | | = 9 and | | ≥ 40. We can also observe that Cplex has solved only 4 instances in the case of problem 26, and just 1 instance in the case of problem 29 out of 5 instances with an average CPU time that exceeded the one-hour time limit. The same fact can be observed for problems 27, 29 and 30 where Cplex was unable to solve no instance within the limit time of one hour.

These observations can be explained by the fact that the developed model induces a large number of decision variables and also by the fact that it considers the joint optimization of the capacited dynamic lot-sizing problem with remanufacturing and the vehicle routing problem with simultaneous pick-up and delivery problem each of them is NP-Hard. In this paper, an extension of the integrated productioninventory-routing problem was studied by considering the case of collecting and remanufacturing end-of-life products within a closed-loop supply chain. We presented a MILP formulation for the considered problem. The proposed model aims to jointly optimize the capacited dynamic lot-sizing with remanufacturing problem and the vehicle routing problem with simultaneous pickup and delivery. This joint optimization in the context of reverse logistics networks represents both an economic interest, as long as we minimize the total cost of the considered operations, and a challenging issue since each of the considered problems is NP-hard. The numerical results show the need to develop approximate resolution methods for resolving large instances. This fact is due to the considerable number of variables induced by the proposed model. Moreover, the next step would be to further develop the proposed model and expand it to cover more realistic real-life situations by considering the multi-item lot-sizing with remanufacturing problem, at production level, and considering a heterogenous fleet of vehicles with divisible deliveries and pickups at distribution level.

  The main scientific contributions of our study are summarized as follows. (a) It provides a modelling approach for the integrated production-inventory-routing problem with returns management. (b) Our study offers a novel MILP formulation in which we conjointly optimize the Capacited Dynamic Lot-Sizing with Remanufacturing and the VRP with Simultaneous Pickup and Delivery. Finally, (c) our study provides managerial insights in the trade-offs among the integrated operational decisions.

FIG. 2 :

 2 FIG. 2: EFFECTS OF NUMBER OF TIME PERIODS (| |) AND THE NUMBER OF CUSTOMERS | | ON THE AVERAGE CPU TIME

[ 1 ]

 1 Ruokokoski, M., Solyali, O. G. U. Z., Cordeau, J. F., Jans, R., & Süral, H. (2010). Efficient formulations and a branch-and-cut algorithm for a production-routing problem. GERAD Technical Report G-2010-66.

TABLE 1 :

 1 GENERATION OF INSTANCES

	Time horizon ( )	Possible values {3, 6, 9}
	Number of customers ( )	Possible values {5, 10, 15, 20}
	Manufactirung production line capacity ( )	Calculated as: 2 • ∑	∈	∑ ∈
	Remanufactirung production line capacity ( )	Calculated as: 2 • ∑	∈	∑ ∈
	Vehicle capacity ( )	>	Calculated as: ∈ {∑ ∈ ; ∑ ∈	}
	Customers delivery
	demand in each period	Uniform distribution U(20;60).
	( )
	Customers pickup demand in each period ( )	Uniform distribution U(5;20).
	Maximum storage capacity of serviceable inventory in each period ( )	Calculated as: 2 • ∑	∈	∑ ∈	• | |
	Maximum strorage
	capacity of retruns inventory in each perdion	Calculated as: 2 • ∑	∈	∑ ∈	• | |
	( )
	Manufacturing setup cost at each period ( )	Possible value 500
	Remanufacturing setup cost at each period ( )	Possible value 500
	Unitary production cost for
	manufacturing in each	Possible value 10
	period ( )
	Unitary production cost for
	remanufacturing in each		Possible value 2
	period ( )
	Unitary holding cost for manufactured product (ℎ )		Possible value 1
	Unitary holding cost of return-product (ℎ )	Possible value 0.8
	Vehicle fixed cost ( )		500
			Calculated as:
	Transportation cost ( )	-	+	-	+ 0.5
	Initial level of serviceable inventory ( )		0
	Initial level of returns inventory ( )		0

TABLE 2 :

 2 SUMMARY OF COMPUTATIONAL RESULTS ((*) SIGNIFICATES THAT THE AVE. CPU (S) WAS CALCULATED BASED ONLY ON THE SOLVED INSTANCES)

	Prob.#	| |	| | #Ins.	Ave. Tot. Cost	Ave. CPU (s)	#OPT
	1		3	5	9644	0,148	5
	2	5	6	5	16681,8	0,432	5
	3		9	5	25725,84	0,48	5
	4		3	5	14418,08	0,964	5
	5	10	6	5	28002,96	3,688	5
	6		9	5	42298,52	1,2	5
	7		3	5	19789,6	2,274	5
	8	15	6	5	39248,36	4,62	5
	9		9	5	57242,44	5,476	5
	10		3	5	26300,84	10,3	5
	11	20	6	5	48405,92	21,16	5
	12		9	5	70064,08	55,318	5
	13		3	5	30230,76	17,272	5
	14	25	6	5	59657,52	111,63	5
	15		9	5	87066,76	102,094	5
	16		3	5	35803,6	60	5
	17	30	6	5	67476,2	160,242	5
	18		9	5	100826,56	184,152	5
	19		3	5	41216,56	136,256	5
	20	35	6	5	77410,4	297,349	5
	21		9	5	114921,44	428,772	5
	22		3	5	46145,4	376,342	5
	23 24	40	6 9	5 5	88316,2 129821,55 (*)	379,292 >3600	5 4
	25 26	45	3 6	5 5	52040,6 96612,15 (*)	299,566 >3600	5 4
	27		9	5	-	>3600	0
	28 29	50	3 6	5 5	56447,92 106717.2 (*)	551,968 >3600	5 1
	30		9	5	-	>3600	0


Emmanuel Kieffer
email: emmanuel.kieffer@uni.lu

Martin Rosalie
email: martin.rosalie@uni.lu

Grégoire Danoy
email: gregoire.danoy@uni.lu

Pascal Bouvry
email: pascal.bouvry@uni.lu

Bayesian optimization to enhance coverage performance of a swarm of UAV with chaotic dynamics

Introduction

We recently propose a mobility model of a swarm of UAVs whose purpose is to cover an unknown area: CACOC (Chaotic Ant Colony Optimization to Coverage) [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. Our algorithm is based on Ant Colony Optimization where chaotic dynamics are used to enhance the exploration part of the algorithm. CACOC mobility model uses repulsive pheromones to guide the UAVs over the area they have to cover. The UAVs share a map of virtual pheromones that indicates recently visited areas when high pheromone concentrations are present. The UAVs then have a higher probability to move to the least recently visited areas. When there is no pheromone to guide the UAVs, the introduction of chaotic dynamics permits to obtain an efficient exploration of the unknown area. Since the chaotic dynamics are obtained from a three differential equations system with parameters, we can tune one parameter to obtain another chaotic dynamic. The aim of this work consists in discovering a chaotic dynamic leading to the best coverage which can be only computed after a CACOC simulation. For this purpose, we consider Bayesian Optimization which has been originally designed for time-consuming black-box optimization. Indeed in order to evaluate the chaotic dynamic for a specific parameter, a full simulation should be realized. S Global optimization techniques (e.g. population-based heuristics) could be very time-consuming contrary to a surrogate model which would minimize the number of simulations to perform to determine the best parameter of the chaotic system. The next section will first introduce the Rossler system and its chaotic system. Section 3 describes how the chaotic system of CACOC has been optimized using Bayesian Optimization. Finally, we will conclude and discuss future works.

Chaotic dynamics of the Rössler system

We used the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] with a parametrization of its parameters [START_REF] Sprott | Asymmetric bistability in the Rössler system[END_REF]

     ẋ = -y -z ẏ = x + ay ż = b + z(x -c) , where      a = 0.2 + 0.09α b = 0.2 -0.06α c = 5.7 -1.18α . (1)
We already provide a detailed analysis of the chaotic dynamics that can be obtained by varying the parameter α [START_REF] Rosalie | Templates and subtemplates of Rössler attractors from a bifurcation diagram[END_REF]. As a chaotic dynamic is a solution of a deterministic process, the periodic orbits of attractors are used to describe the topological structure of chaotic attractors. We used these periodic orbits to obtain recurrent pattern in our exploratory movement. In our original method [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF], we found three patterns: straight lines (period 1 orbit), large turns (period 2 orbit) and serpentine patterns (period 4 orbit) when α = -0.25.

For α ∈ [-0.8; 0.4] we know that the global dynamic is only composed by an increasing branch followed by a decreasing branch. However, when α varies, the orbits are not the same. For instance, this can lead to the apparition of a new orbit of period four, itself leading to a pattern with very large turn. This can be instead of the serpentine pattern or in addition. As it is fastidious to explore all the periodic orbits changes when α varies for this range of values, we choose to use Bayesian optimization to find out the best α for CACOC. We used a metric developed to evaluate CACOC: the slope of coverage. It permits to evaluate how the 10 UAVs spread on the area and then, how they manage to visit unvisited area. For a given simulation, the percentage of covered area is ploted against the number of steps, the slope of coverage is the linear regression of this percentage over the number of step: the higher the value, the best it is.

Bayesian optimization

Bayesian optimization belongs to surrogate-based optimization algorithms. It has been originally designed to tackle black-box optimization problems where the evaluation of a solution could be very time-consuming. Indeed, many optimization problems do not have a formal mathematical expression because it is either unknown or difficult to obtain. To cope with such issues and minimize the number of evaluations, a surrogate model based on Gaussian processes [START_REF] Mockus | Bayesian approach to global optimization: theory and applications[END_REF] is created to approximate the unknown objective function. Bayesian optimization samples promising zones in the feasible region by computing a distribution of the objective function. This distribution give us a prior knowledge on location of the optimal solution. Bayesian optimization is thus characterized by a probability measure on F describing our prior beliefs on F as well as an acquisition function which allows to gain information on the location of the minimum value of the objective function. It has been succesfulluy used to tackle bi-level optimization problems (cite Manu) to reduce dramatically the number of evaluation of the objective function while providing accurate solution values. In this work, the objective function, i.e. coverage cannot be formally described and a simulation has to be run to measure the influence of the α parameter on the chaotic dynamic. Therefore, we attempt to create a surrogate model using Gaussian process which will be refined until no improving alpha can be obtained. Due to the non-linearity of the problem, we consider the Matern32 kernel to define the covariance matrix of the Gaussian processes. Promising area are sample using the optimization of the Lower Confidence Bounds acquisition function.

Conclusion

In this work, we introduced the optimization of CACOC through Bayesian Optimization. CA-COC is based on a chaotic system, i.e. Rossler system whose behavior can be modified by tuning the α parameter. In order to evaluate the performance of CACOC for different value of α, the coverage metric has to be evaluated after simulation. The latter is time-consuming. Therefore, a surrogate-based optimization, i.e. Bayesian Optimization has been privilegied to tackle this issue. Experiments have been conducted on the High Performance Computing (HPC) platform of the university of Luxembourg [START_REF] Varrette | Management of an academic HPC cluster: The UL experience[END_REF]. An analysis of the chaotic system with the obtained α value has been performed to compare the periodic orbits and their associated patterns. Numerical results show that the best α value avoid a waste of time in periodic region of the bifurcation diagram. Future works will focus on more complex chaotic system as well as new application domain of the optimized CACOC approach.