Kateryna Omelchuk

Alexandre Moncef Stambouli

Alexandre Chagnes

Diphasic acid-basic properties of new organophosphorus extractants

published or not. The documents may come

Diphasic acid-basic properties of new organophosphorus extractants

Acidic organophosphorus compounds have been extensively studied as extractants for numerous metals and have been shown to exhibit a good extraction and separation efficiency [1]. The overall distribution equilibrium of an organophosphorus extractant HL between organic and aqueous phases involves a great number of individual equilibria [2]. However, only three of them are usually taken into account in the literature: the dimerization of HL in organic phase, the partition equilibrium of HL between the two phases and its acid dissociation in aqueous phase [3]-[4]. The metal extraction by organophosphorus extractants requires to increase the pH of acid leaching solution by adding of alkaline reactive such as sodium hydroxide for instance. In some cases, saponification of extractants can occurs and may lead to change in extraction equilibria of metals.

In the present work, saponification of the following extractants diluted in kerosene has been investigated by performing biphasic titration:

The titration curves have been modelled by means of a physicochemical model including the equilibria reported in Figure 2. The formation constant β 1 , β ' , β '' and β n , the dimerization constant (k 2 ) and the dissociation constant (k a ) of the acidic extractant have been used in the model. The values of β 1 , β', β '' and β n have been determined by fitting the experimental data of biphasic titration with the model while k2 and ka values have been determined by VPO measurement and monophasic titration, respectively.

Finally, the physicochemical model is used to investigate the speciation of the extractants as a function of the pH.

[1] Z. Kolarik, Solvent Extr. Ion Exch., vol. 28, no. 6, pp. 707-763, Oct. 2010. [2] C. Sella and D. Bauer, Solvent Extr. Ion Exch., vol. 6, no. 5, pp. 819-833, 1988. [3] X. Fu, Z. Hu, Y. Liu, and J. A. Golding, Solvent Extr. Ion Exch., vol. 8, no. 4-5, pp. 573-595, Jun. 1990.

[4] R. K. Biswas, M. A. Habib, and M. N. Islam, Ind. Eng. Chem. Res., vol. 39, no. 1, pp. 155-160, Jan. 2000.

Key words: ORGANOPHOSPHORUS EXTRACTANTS, MODEL, SAPONIFICATION, EQUILIBRIA

Figure 1 .

 1 Figure 1. Structure of organic extractants synthesized in the laboratory.

Figure 2 .

 2 Figure 2. Equilibria included in the physicochemical model.


Christian Artigues

Emmanuel Hébrard

Yannick Pencolé

Andreas Schutt

Peter J Stuckey

Laas

Data Instance

published or not. The documents may come

Introduction

The overall objective of the GEO-SAFE project [START_REF]Geo-safe -geospatial based environment for optimisation systems addressing fire emergencies[END_REF] is to develop methods and tools enabling to set up an integrated decision support system to assist authorities in optimizing the resources during the response phase to a wildfire (fire suppression, life and goods protection). One critical and crucial part of this integrated decision support is the ability to perform largescale evacuation planning [START_REF] Shahparvari | Enhancing emergency evacuation response of late evacuees: Revisiting the case of australian black saturday bushfire[END_REF]. While in practice most evacuation planning is principally designed by experts using simple heuristic approaches or scenario simulations [START_REF] Veeraswamy | The simulation of urban-scale evacuation scenarios: Swinley forest fire[END_REF], more recently optimization approaches to evacuation planning have been addressed, using a variety of optimization technology, as surveyed recenly in [START_REF] Bayram | Optimization models for large scale network evacuation planning and management: A literature review[END_REF]. This paper presents a challenging variant of the evacuation planning problem in case of wild fire issued from exchanges with practitioners in the context of the GEO-SAFE project and from a specific literature review. A large amount of work has been carried out, notably at NICTA 1 , [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF][START_REF] Even | NICTA evacuation planner: Actionable evacuation plans with contraflows[END_REF][START_REF] Pillac | A pathgeneration matheuristic for large scale evacuation planning[END_REF][START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF][START_REF] Even | Convergent plans for large-scale evacuations[END_REF][START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF][START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF][START_REF] Romanski | Benders decomposition for large-scale prescriptive evacuations[END_REF] mainly in the context of floodings, which can be transposed under some adaptations to evacuation in case of fires. Evacuation planning in case of wildfires is indeed much harder. Wildfire propagations are inherently less predictable than floods. While flood levels mostly rely on the fixed topology of the area and rainfalls, wildfire mainly depends on the wildland fuels [START_REF] Rothermel | A mathematical model for fire spread predictions in wildland fuels[END_REF][START_REF] Anderson | Aids to determining fuel models for estimating fire behavior[END_REF], on the slope of the burning ground and more importantly on the speed and direction of the wind that can suddenly change at any time [START_REF] Weise | A qualitative comparison of fire spread models incorporating wind and slope effects[END_REF][START_REF] Stepanov | Modeling wildfire propagation with delaunay triangulation and shortest path algorithms[END_REF]. Therefore, evacuation planning dedicated to wildfires must be much more robust to difference future scenarios. A good evacuation plan in case of wildfire must not only minimize the evacuation time of the population but also maximize the spatial and temporal safety margin between the evacuees and the actual or potential wildfire front.

This paper present a new optimization model for evacuation planning in the event of wildfire as well as a problem instance generator. On these instances, basic 0-1 integer programming and constraint programming formulations enhanced with a simple compression heuristic are compared. In Sec-tion 2, we provide a literature review and we define the considered problem. Section 3 presents the instance generator. The basic 0-1 integer programming, constraint programming formulations and the heuristic are proposed in Section 4. Computational experiments are given in Section 5.

2 State-of-the-art review and problem definition

Basic evacuation data

We adopt the notation an terminology given in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF]. There is a directed graph G(N = E ∪ T ∪ S, A) representing e.g. a road network in a region that must be evacuated. The graph is made of:

• the set of evacuation nodes E. An evacuation node represents a zone where people to be evacuated are regrouped),

• the set of safe nodes S. A safe node represent a safe geographical zone that people located in the evacuation nodes must reach during the planning horizon

• the set of transit nodes T . A transit node represent an intersection in the road network that can be traversed by the vehicles carrying evacuated people from the evacuation zones to the safe zones.

Suppose the time is discretized in minutes. Each evacuation node k ∈ E is associated with a number of evacuees d k . Each arc q ∈ A has a capacity u e in evacuees/minute, a travel time t e giving the number of minutes that a vehicle takes to traverse the arc and a end time b e , which gives the time at which the arc becomes unavailable due to the fire propagation. Other characteristics may appear in the variants considered in the cited papers.

Previously considered evacuation problems

We limit our review to the the NICTA papers as they explain the genesis of the final model w.r.t. practical considerations. The first considered problem and solution methods were presented by Pillac et al. in 2013 in a research report that was later published in [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF]. One of the practical motivation of the work was that in an urgency situation, the possibility of choices in a road network, such as a fork, generate stress among the population. This is why it is preferable to define a single evacuation path for each evacuation zone. Hence the problem considered in [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF] can be described as a two-level optimization problem, which can be solved in an integrated way or by a twophase approach. At the first level, exactly one evacuation path is determined from each evacuation node to a single safe node. At the second level the flow of evacuees is scheduled trough the path. A MILP (called the restricted flow model RF) is proposed to solve the problem in an integrated way. It is based on time-expanded graph, in which each node is duplicated for each time period. Arc capacities u e in the time expanded network ensure that the traffic does not exceed the road capacity, especially when several paths use the same arc. The model includes continuous flow variable per evacuation node per arc and in the expanded network and also a binary variable per edge in the (nonexpanded) graph G to ensure the uniqueness of the path for each evacuation node. As the MILP is intractable on a set of instances (HN) derived from a real case scenario involving 70000 evacuees in the Hawkesbury-Nepean floodplain, located North-West of Sydney, an iterative column-generationlike two-phase heuristic is proposed. Given a set of potential paths, a master problem solves the path selection and flow scheduling problem with another MILP (called the conflict-based heuristic path generation master problem, CPG-MP) involving a reduced number of binary constraints (one per path). A sub-problem finds new paths based on a subset of critical nodes by solving a multiple-origin, multiple-destination shortest path problem. A lexicographic objective is considered. The main objective is to maximize the number of evacuated people during the time horizon, while a secondary objective aims at maximizing the evacuation start time, based on the practical relevance. Note that the second objective is only indirectly tackled by weighting each arc in the time-expanded graph by a cost inversely proportional to the arc time, hence solving a min-cost flow problem.

Subsequent papers describe variants of the models and new methods. A more sophisticated variant of the heuristic was proposed in Pillac et al. (2014) [START_REF] Pillac | A pathgeneration matheuristic for large scale evacuation planning[END_REF] for the same problem with, however, a different secondary objective for the min-cost flow, aiming at minimizing the evacuation time (called the clearance time), which is indirectly obtained by weighting the arc in the time-expanded graph by a cost proportional to the arc time.

In [START_REF] Even | NICTA evacuation planner: Actionable evacuation plans with contraflows[END_REF] [START_REF] Even | NICTA evacuation planner: Actionable evacuation plans with contraflows[END_REF], an additional degree of freedom is introduced, giving the possibility of selecting contraflow roads, which consists in reversing the direction of some major roads. In practical cases, this possibility can highly increase the network capacity for evacuation. This is done by introducing a subset A c of arcs that can be used in contraflows and by modifying the CPG-MP, interestingly without the need of introducing additional binary variables. A software called the NICTA Evacuation Planner is also presented, with new instances with up to 1 000 000 evacuees, which are solved requiring up to 30 min. of CPU time.

In Even et al. (2015) [START_REF] Even | Convergent plans for large-scale evacuations[END_REF], the case where the selected paths must form a convergent subgraph, i.e. for which each node has at most one outgoing arc in the set of paths, is considered. This is motivated by the fact that convergent evacuation paths can be better controlled. Indeed, even if in the previous model, a single path was assigned to each evacuation node, in the case where two paths are merged and then are separated again, a driver can be confronted to a choice and take a decision that would not correspond to what was planed. Hence the set of paths now form a tree rooted at the safe node. The solution method is still a two-phase method but does not follow anymore the column generation principle. A tree is first built by a MIP working on an aggregated graph (without time discretization). The second phase is the flow scheduling problem given the computed tree, which is a maximum flow problem. The first phase is now able to produce good upper bounds on the maximum number of evacuees and on the minimum clearance time. A simulation study shows that the convergent model outperforms the general model when the presence of a fork in the network generated an hesitation for the driver that is capture by a 0.75 second delay.

This model is further developed by [START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF] [START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF] to incorporate network design aspects in the problem to model possible infrastructure enhancement decisions, as for the west Sydney case concerns about the capacity of the road network were expressed. Two additional decisions are introduced : adding lanes to a road (arc) and elevating a road (arc). The first decision results in an increase of the capacity of the arc while the second decision postpones the flooding time, both coming at a cost. These decisions can be incorporated in the tree design problem, which is also the master problem of a Benders decomposition method. The objective is still to maximize the number of evacuated people under a budget constraint that limits the infrastructure upgrades. The maximum flow subproblem is used to obtain a feasible flow schedule and also to generate optimality cuts that are reinjected in the master problem. To increase the performance of the Benders method, pareto-optimal cuts are generated.

All previous approaches assumed that each individual vehicle of an evacuation (flow unit) moves independently from the other vehicles of the same zone, and exactly as prescribed by the maximum flow model. Pillac et al. [START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF] propose to use the concept of response curves to incorporate behavioral models in the problem. In practice, instead of assigning a start time to each evacuee, the authorities can influence the evacuation start time of a zone and the level of resources mobilized to increase the evacuation rate (e.g. number of agents knocking on people's door), to which people answer according to a behavioral model abstracted by a response curve. It follows that to each evacuation zone k ∈ E, a set of response curves F k is given. For each evacuated zone, a start time δ k and a response curve f f ∈ F k has to be assigned such that the flow of evacuees leaving zone k at time t is givn by

φ t k = 0 if t < δ f (t -δ) if t ≥ δ.
A column generation approach is proposed where the set of all evacuation plans of a zone is considered, where a evacuation plan is a path from the evacuated zone to a safe node, the start time and the response curve. The master problem selects a plan for each evacuation zone in such a way that the network capacity is satisfied and the total cost of the plans is minimized (without ensuring that a convergent is obtained). The plan generation subproblem is solved either via a MIP or via a shortest path approach. The methods experiences difficulties is solving realistic instances due to a longtail effect.

Another drawback of the flow model is that it generated preemptive evacuation plans. As flow units are routed independently there are time periods in which the evacuation of a zone may be stopped and reinitiated later. Even et al. (2015) [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] report that this creates serious implementation issues for the evacuation plans. They propose another model in which the evacuation rate is a decision variable that remains fixed as soon as the evacuation starts and as long as the zone is not fully evacuated or the time horizon is reached. They assume that the evacuation path of each evacuation zone is already determined and concentrate only on the scheduling problem. We present this model in the case that all people must be evacuated and the objective is to minimize the clearance time. Let δ k denote the start time of evacuation of zone k and let λ k denote the evacuation rate of zone k.

Let A k = (e 1 k , . . . , e |A k | k
) the evacuation path of zone k given as its list of edges. An evacuation plan defines a task J i k for each edge e i k with a start time S i k = δ k + i-1 q=1 t q k where t i k is the travel time of edge e i k . The evacuation has a variable duration p k and the total number of evacuees is p k λ k with p k λ k = d k since all people must be evacuated. Each edge task e i k ∈ A k has duration p k . Let J e the set of tasks that use edge e ∈ A. We also denote as u k the maximum evacuation rate of an evacuation node k, which can be seen as a node capacity. The model can be written as follows :

min C max (1) s.t. C max ≥ S |A k | k + p k ∀k ∈ E (2)
S i k + p k ≤ b e i k ∀k ∈ E, ∀i = 1, . . . , |A k | (3)
p k λ k = d k ∀k ∈ E (4)
S i k = δ k + i-1 q=1 t q k ∀k ∈ E; ∀i = 1, . . . , |A k | (5)
J i k ∈Je,S i k ≤t<S i k +p k λ k ≤ u e ∀e ∈ A, ∀t ≥ 0 (6)
δ k ≥ 0 ∀k ∈ E (7) u k ≥ λ k ≥ 0 ∀k ∈ E (8) p k ≥ 0 ∀k ∈ E (9)
This is a no-wait total work-and resource-constrained project scheduling problem where (4) are the total work constraints, (5) are the no-wait constraints, as the start time of the evacuation task on a arc of the evacuation path starts exactly at the decided evacuation time δ k plus the total travel time along the path toward edge e. Constraints ( 6) are the capacity constraints on each edge e. The problem was efficiently solved in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] via constraint programming for both the clearance time minimization version and the maximization of the number of evacuated people variant.

A new evacuation planning problem in case of wild fires

In the work of Even et al. (2015) [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] and previous studies, the time b e at which arc e becomes unavailable comes from a flood propagation model, which is pretty accurate. In the case of fire, even if precise propagation models can be obtained, they depend on multiple parameters. Among them, the wind has a great variability. The subject of promising further studies would be to consider explicitly uncertainty via robust of stochastic approaches, but this would require the definition of various scenarios possibly associated with probability distributions. An alternative to modeling uncertainty of the unavailability dates would be to consider an objective function that seeks to maximize the length of the time interval [C e , b e ], where C e is the completion time of the last task using edge e. It follows that in this paper we consider the following optimization problem, with objective function [START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF] giving the maximum (possibly negative) slack on each evacuation arc, weighted by the population.

max min k∈E,∀i=1,...,|A k | d k (b e i k -S i k + p k ) (10)
s.t.( 4), ( 5), ( 6), ( 7), ( 8), (

3 Realistic data instance generator Catastrophic wildfire requiring large population evacuation are, thankfully, rare events. However, it means that obtaining useful data is difficult, and indeed this a key problem within the GEO-SAFE project. A significant part of the project revolves around simulation tools such as EXODUS [START_REF] Galea | The EXODUS Model[END_REF], however, even simulated data was hard to come by. Therefore, we opted for taking advantage of the project environment to contribute to this effort by generating our own "realistic" dataset. On the one hand, this approach may introduce biases since we must use models to generate realistic road networks and simulate wildfires. On the other hand, we believe that it will make it much more convenient for benchmarking algorithms in the future. As it turns out, the generated instance are challenging even though relatively modest in size, thus being interesting from an academic viewpoint as well.

Generation of road networks

The first step is to generate a graph standing for the road network. To this end, we used the quadtree model described in [START_REF] Eisenstat | Random road networks: the quadtree model[END_REF]. In a nutshell, this model starts with a single square formed by four nodes and four edges. At each iteration, a square is chosen and five nodes are added, one in the center of the square, and one on each edge connected by a perpendicular edge to the center node. A parameter r controls the sprawl, that is, the preference for splitting larger squares (r < 1) or smaller square (r > 1). The graphs generated in this manner share a many features with real road networks: they are planar, embedded in an Euclidean plane, have similar density distributions, path lengths are within a constant factor of the Euclidean distance, and the number of turns is logarithmic with high probability. An example of random quadtree network is shown in Figure 1a. The colors on the edges correspond to road capacity. To allocate capacities, we first compute a minimum Steiner tree spanning three randomly chosen nodes in high density areas ("cities") and connect these cities to the nearest corner of the outer square. The corresponding set of edges are given the highest capacity and are coloured in blue in Figure 1a. A second set of edges, forming a grid are given an intermediate capacity, they are coloured in green.

Simulating wildfire

The second step consists in determining safety due dates for every edge of the evacuation tree, that is, a time after which the edge become unsafe. To this purpose we use a relatively simple fire propagation model. We chose to use a simple model based on two parameters: a constant intensity γ representing the type of fuel material as well as the temperature, and a wind direction. Indeed, the goal is not so much to accurately predict fire propagation, as it is to generate safety due dates consistent with a wildfire. Of course, should the authorities use this type of planning tools during a real event, then correctly predicting fire propagation would be among the most important factor. The land area is discretized into squares of fixed size (we use another parameter to control this size) which can be in three states: untouched, burning and burned. The fire starts as a single burning square, then at each iteration, any untouched square adjacent of a burning square catches fire with probability γ( π-A π ) 2 , where γ stands for the intensity of the fire, and A is the angle between the wind an a vector going from the center of burning square to the center of the untouched square. Moreover, any burning square that did not propagate stop burning with probability γ 2 . Figure 1b illustrate the state of the simulated wildfire, with burning squares in red and burned squares in black.

Generating evacuation plans

The third step consists in generating the actual evacuation plan, that is, an embedded tree connecting a set of evacuation nodes E to a safe node r. Here again, the goal is not to compute the best evacuation plans, however they must be representative of what would be actual plans.

We first randomly pick a predefined number of evacuation nodes among the nodes of the graph that are in the state burned or burning of the simulated fire. Then we use the convention that the safe zone is the furthest corner from the center of the fire. The evacuation tree is computed simply by using a shortest paths algorithm, however with respect to an arc labeling taking into account first the safety due date of the arc, and only then its length and its capacity.

At this point we have all the information we need to define a fire evacuation problem as defined in Section 2.3.

The tools we developed as well as the benchmarks instances we used in this paper can be accessed here: https://github.com/ehebrard/evacsim.

Formulations and heuristic

In this section, we propose two formulations for the problem: a 0-1 integer linear programing formulation and a constraint programming formulation. Then, we describe a simple compression heuristic able to find quickly an initial solution.

0-1 linear programming formulation

Let H denote an upper bound on the latest evacuation completion time on the evacuation nodes (the time by which the last evacuee leaves the evacuation node). We propose an integer 0-1 linear formulation that makes a discrete approximation of the problem. The set of discrete evacuation possible start times is equal to

H k = {0, . . . , H -p k } where p k = d k u k
is the smallest possible integer evacuation processing time. For an evacuation start time t ∈ H k , the set of possible evacuation processing times is P kt = {p k , . . . , min(d k , H -t)}. Let p mk ∈ P kt denote the m th smallest possible processing time for m = 1, . . . , |P kt |. We introduce a 0-1 variable x ktm equal to 1 if and only if

δ k = t, p k = p mk et λ k = d k p mk .
By analogy to multimode scheduling problems, set {1, . . . , |P kt |} represent the set of processing modes available for scheduling an evacuation task that starts at time t. For a given mode m and a given time t, we denote by I i ktm ⊆ H k the maximal discrete time interval such that

∀τ ∈ I i ktm , t ∈ [τ + i-1 q=1 t q k , . . . , τ + i-1 q=1 t q k + p mk -1]
i.e. the set of evacuation start times τ in mode m that make evacuation on edge e i k in process at time t. This interval is previsely:

I i ktm = {t - i-1 q=1 t q k , . . . , t - i-1 q=1 t q k -p mk + 1} ∩ H k
Given these elements , the problem can be expressed as the following 0-1 integer linear program:

max W min (12) s.t. W min ≤ d k (b e i k - t∈H k |P kt| m=1 (t + p mk )x ktm - i-1 q=1 t q k ) ∀k ∈ E, ∀i = 1, . . . , |A k | ( 13
)
t∈H k |P kt| m=1 x ktm = 1 ∀k ∈ E (14)
J i k ∈Je |P k0 | m=1 τ ∈I i ktm λ km x ktm ≤ u e ∀e ∈ A, ∀t = 0, . . . , H -1 (15)
x ktm ∈ {0, 1} k ∈ E, t ∈ H k , m ∈ P kt (16)

Constraint Programming formulation

In [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF], the NEPP was modeled using standard cumulative constraints. We adapt here this model for our problem. Let x (resp. x) denote the largest (resp. smallest) value in the domain of a variable x. Given a set of tasks J with start time variable

s i ∈ [s i , s i ], processing time variable p i ∈ [p i , p i ],
height variable λ i ∈ [λ i , λ i ] and a resource r of constant capacity u r , recall that cumulative((s i , p i , λ i ) i∈J , u r ) enforces the relations

i∈J|s i ≤t<s i +p i λ i ≤ u r ∀t ∈ H
Consequently, to model the problem, it suffices to associate a task J i k to each arc on the evacuation path of each evacuation node k ∈ E, with height variable

λ k ∈ [1, u k ], start time variable S i k ∈ [0, H -dv uv ], duration variable p k ∈ [ dv uv , d v ].
A resource is defined per arc e ∈ A, with capacity u e . The baseline constraint program for the evacuation planning problem is obtained by replacing constraints [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] in the problem formulation of Section 2.3 by: cumulative((S i k , p k , λ k ) J i k ∈Je , u e ) ∀e ∈ A (17)

Heuristic

We propose a simple compression heuristic to find an initial upper bound.

The heuristic is based on the assumption that scheduling all evacuation tasks at time 0 with the minimum evacuation rate yields a feasible solution, with a high cost. Starting from this solution (∀k ∈ E, we set the start time s k := 0, the end time e k := d k , and the rate λ k := 1). Now an iterative process starts where, at each iteration, the critical evacuation tasks, i.e. the one that minimizes the cost on some edge, is identified. Then, its duration is decreased and its height is consequently increased until (i) either no more height increase/duration decrease can be performed without exceeding an edge capacity or (ii) the task is not critical anymore and another task becomes critical. If case (i) occurs the process stops, otherwise, if case (ii) occurs, the compression process restarts with the new critical task unless the objective increase is smaller than a predefined parameter , in which case the process also stops. Due to the possibility of only left shifting a task k by d k at each iteration, this descent heuristic is of pseudo polynomial computational complexity.

Computational experiments

We generated 240 benchmark instances following the protocol described in Section 3. They are organized into three types of road networks: Dense, Medium and Sparse where the density refers to the number of intersections (respectively 400, 800 and 1200) in the land area. Notice that the graph has always 4 edges per node, so this corresponds to graph size. The impact on the instance is that larger graphs allow for more choices for the shortest paths and therefore longer independent paths. For every type of road network, we generated 4 classes of instances, with respectively 10, 15, 20 and 25 evacuation nodes. Finally, for every class we simulated 20 random wildfires and the subsequent evacuation trees. We used CPLEX 12.7 to solve the MILP formulation with default settings and CPOptimizer 12.7 for solving the CP formulation. The heuristic was used to provide an initial solution to both solvers. We ran every method on every instance of the dataset with a time limit of 45 minutes on 4 cluster nodes, each with 35 Intel Xeon CPU E5-2695 v4 2.10GHz cores running Linux Ubuntu 16.04.4.

We provide here e few implementation details. We only used discrete evacuation rates as it was simpler to implement in the CP solver. It follows that Constraints (4) was implemented as an inequality (≥). The number of edges on which to check the cumulative constraint was reduced thank to the observation that in a path on an evacuation only one edge is a bottleneck. It follows that only one edge per path has to be considered as a limited resource. Last, the opposite of the weighted slack was actually minimized, which amounts to a maximum weighted lateness objective. The results are displayed in table 1, where we give for each solver the average upper bound on each instance family and the optimality ratio, i.e. the percentage of verified optimal solutions found. We can first remark that the generated instances are hard to solve optimally for both solvers: on average only 51% of the instances are solved to optimality with the same behavior for the families : almost all the instances with10 evacuation node and a large part of the instances with 15 evacuation nodes car be solved to optimality, while instances with 20 evacuation nodes becomes much harder and the ones with 25 activities are intractable. As an outcome, our generator is able to produce computationally challenging instances.

In terms of comparison between integer and constraint programming, the integer program is significantly outperformed by CP, both in terms of optimality ratio and of upper bounds on maximum weighted lateness. This is both due to slower convergence time and memory issued due to the huge size of the IP model. As a typical example, the instance medium 10 30 3 2 has 328147 binary variables and 3198 constraints after CPLEX preprocessing.

In terms of the obtained objective function values, on the 240 instances only 45 have negative values, meaning that in a majority of instances the evacuation could not be performed on some edge before the expected deadline. Interestingly, all these 45 instances were solved to optimality, which represents 36% of the 123 instances solved to optimality. 74% of the remaining instances, correspond to pessimistic scenarios where the evacuation road network is unable to ensure the evacuation of the the whole population before the traversal of some route segment would become critical. If such situation occurred in actual road network this could give helpful support to the authorities for increasing the capacity of specific road segments or to build better prescribed evacuation routes.

Concluding remarks

We have proposed a data instance generator and optimization frameworks for a computationally challenging evacuation planning problem, with an objective function tailored to the event of wildfire. This generator could be improved by incorporating more sophisticated fire propagation models and actual road networks. The generation of evacuation routes is also an optimization problem in itself. Feedback from the evacuation planning, in the case, which often occurred in our experiments, where obtained safety margin are not sufficient should be used to modified the evacuation routes accordingly. In terms of the evacuation planning problem we have proposed new integer and constraint programming formulation. To obtain competitive results with IP, one should obviously consider decomposition approach as the problem is huge. Continuous models could also be designed to reduce the number of variables. Even if CP obtains much better results, the vast majority of medium size instances could not be solved to optimality. As future research directions, we will specific global constraint that better capture the structure of the problem as well as dedicated search strategies. Finally, we believe that coupling our approach with simulation and/or stochastic-robust optimization will lead to useful decision support tools in case of response to wild fires.

Figure 1 :

 1 Figure 1: An example of generated instance

  dense 15 (20) 161100 0.70 346404 0.00 dense 20 (20) 311676 0.45 659395 0.00 dense 25 (20) 531016 0.00 1044248 0.00 medium 10 (20) 48591 1.00 103151 0.55 medium 15 (20) 124921 0.70 306552 0.05 medium 20 (20) 276101 0.25 581273 0.00 medium 25 (20) 488282 0.10 955786 0.00 sparse 10 (20) 3747 1.00 14359 0.70 sparse 15 (20) 120173 0.65 288432 0.05 sparse 20 (20) 235771 0.35 554373 0.00 sparse 25 (20) 437618 0.00 868438 0.00 avg (240) 230191 0.51 481763 0.17

Table 1 :

 1 CPO vs MIP: upper bound and optimality ratio

Acknowledgement

This work is partially funded by the H2020-MSCA-RISE-2015 European project GEO-SAFE (id 691161)


Christian Artigues

Emmanuel Hébrard

Yannick Pencolé

Andreas Schutt

Peter J Stuckey

Emmanuel Hebrard

Yannick Pencole

A Study of Evacuation Planning for Wildfires

published or not. The documents may come

L'archive ouverte

Introduction

In EU and Australia, every year thousands of square miles of forests and other lands burn due to wildfires. These fires cause important economic and ecological losses, and often, human casualties as for instance during the Black Saturday bushfires across the Australian state of Victoria in February 2009 [START_REF] Shahparvari | Enhancing Emergency Response in Short-notice Bushfire Evacuation[END_REF]. The overall objective of the GEO-SAFE project [START_REF]Geo-safe -geospatial based environment for optimisation systems addressing fire emergencies, MSCA-RISE-2015 -Marie Sk lodowska-Curie Research and Innovation Staff Exchange (RISE) European[END_REF] is to develop methods and tools enabling to set up an integrated decision support system to assist authorities in optimizing the resources during the response phase to a wildfire (fire suppression, life and goods protection). One critical and crucial part of this integrated decision support is the ability to perform large-scale evacuation planning. As detailed in [START_REF] Shahparvari | Enhancing emergency evacuation response of late evacuees: Revisiting the case of australian black saturday bushfire[END_REF], there are commonly three categories of evacuees: the ones that leave early, the ones that shelter in refuge and the ones who stay at their properties and fight. This work focuses on the evacuation of this third group (late evacuees) which is called the late evacuation planning problem.

While in practice most evacuation planning is principally designed by experts using simple heuristic approaches or scenario simulations [START_REF] Veeraswamy | The simulation of urban-scale evacuation scenarios: Swinley forest fire[END_REF], more recently optimization approaches to evacuation planning have been addressed, using a variety of optimization technology [START_REF] Bayram | Optimization models for large scale network evacuation planning and management: A literature review[END_REF]. The usual test case for most of this work is flood evacuation planning [START_REF] Pillac | A conflict-based path-generation heuristic for evacuation planning[END_REF][START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF][START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF]. Given accurate measurements of rainfall and topology, flood evacuation planning can make use of very accurate predictions of future water levels, and therefore has a very accurate model of what infrastructure will be available at each stage of the evacuation.

Evacuation planning in case of wildfires is much harder. Wildfire propagations are inherently less predictable than floods. While flood levels mostly rely on the fixed topology of the area and rainfalls, wildfire mainly depends on the wildland fuels [START_REF] Rothermel | A mathematical model for fire spread predictions in wildland fuels[END_REF][START_REF] Anderson | Aids to determining fuel models for estimating fire behavior[END_REF], on the slope of the burning ground and more importantly on the speed and direction of the wind that can suddenly change at any time [START_REF] Weise | A qualitative comparison of fire spread models incorporating wind and slope effects[END_REF][START_REF] Stepanov | Modeling wildfire propagation with delaunay triangulation and shortest path algorithms[END_REF]. Therefore, evacuation planning dedicated to wildfires must be much more robust to difference future scenarios. A good evacuation plan in case of wildfire must not only minimize the evacuation time of the population but also maximize the spatial and temporal safety margin between the evacuees and the actual or potential wildfire front.

In this paper, we consider that the authorities already identified the accessible routes and shelters and estimated the population of the late evacuees. Evacuation takes place in individual vehicles and each center has a known population and a single predefined suggested evacuation route. All routes issued from each center converge toward a safe place, so congestion may appear on route segment shared by several evacuation paths. The interest of convergent evacuation plans has been underlined by several studies such as in [START_REF] Even | Convergent plans for large-scale evacuations[END_REF] since it avoids congestion issued by driver slow-downs at forks. Furthermore fire propagation models give a deadline on each route segment beyond which taking this segment comes at high risk. To mitigate these risks, the authorities may delay the evacuation start time and rate for each center. Indeed in practice, besides the possibility of assigning a start time to each evacuated zone, the authorities can also mobilize different levels of resources to increase the evacuation rate (e.g. number of agents knocking on people's door), to which people answer according to a behavioral model abstracted by a response curve [START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF]. Based on these concepts, the evacuation planning model is close to the one proposed in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] called the non-preemptive evacuation planning problem (NEPP) in the context of a flood. In contrast with previously proposed models, the NEPP considers non preemptive evacuation. Once the evacuation of a zone has started, it cannot be interrupted. Indeed, considering preemptive evacuation would make the problem much easier to solve but is hard to implement in practice. In case of wildfire, this would notably cause undesirable stress on evacuees. The major difference of the model considered in this paper with the NEPP model lies in the objective function. In [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF], the main objective was to maximize the number of evacuees, and a secondary objective was to minimize the makespan, i.e. the total evacuation time, while enforcing deadlines on route segments. In case of wildfire, the high variability in the fire propagation makes it necessary to avoid taking a route segment close to the expected deadline. Hence, we consider a single objective by maximizing the minimum gap for all evacuated zones and all route segments between the deadline and the time by which the last evacuee leaves a route segment, weighted by the population of the evacuated zone.

The evacuation planning problem

We are given an a tree G = (E∪T ∪{r}, A) rooted in r standing for the evacuation routes from evacuation nodes E (leaves) to the safe node r through transit nodes T . Every leaf/evacuation node v ∈ E is associated with a population count w v . Every arc u, u has a length l uu and a capacity non-ambiguously denoted q u as G is a tree. Moreover the length of the path from a node u to a node u is also written l uu .

Let H = [0, H] be the time span of the evacuation. We want to associate every leaf/evacuation node v ∈ E to a real s v representing the delay of the evacuation notice and to a "response curve" φ v describing the evacuation flow out of node v over time (starting from s v ). Population flows out of leaf/evacuation node v at rate φ v (t) that is zero before time s v ≥ 0, and such that:

H 0 φ v (t)dt = w v
Let p(u) denote the parent of u in G, p(u) its ascendants, C(u) its children, Ĉ(u) its descendants and L(u) those of its descendants that are leaves of G. We assume that the evacuees never stop, so the flow in the downstream arcs can be computed by summing the flows in the incoming arcs and the flow φ u (t) for any arc ∈ E ∪ T at any time t ∈ H is

φ u (t) = v∈L(u) φ v (t -l uv -s v )
In this paper we consider as in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] a simple response curve. The flow out of an evacuation node v ∈ E is a continuous decision variable and remain constant during the evacuation process, that is, the flow out of node v ∈ E is equal to h v within the time interval [s v , e v ], with e v = s v + wv hv and zero otherwise. Therefore, the flow out of a node u ∈ E ∪ T at time t ∈ H is:

φ u (t) = v∈L(u),sv+luv≤t<sv+luv+ wv hv h v
It follows that the considered evacuation planning problem can be formally defined as the following constraint optimization problem.

Variables: We have a set of |E| non-preemptive tasks, one for every evacuation node. For every task standing for an evacuation node v ∈ E we need two variables, one for the constant rate h v ∈]0, q v ] at which the evacuation will proceed and one for its start time s v ∈ [0, H -wv qv ] at which the evacuation will start.

Constraints:

We have a single type of constraints to avoid "jams", i.e., flow exceeding the capacity of an arc. For each node u ∈ T , we have a cumulative resource constraint of capacity q u ensuring that φ u (t) ≤ q u , which is written:

v∈L(u),sv+luv≤t<sv+luv+ wv hv h v ≤ q u , ∀u ∈ T , ∀t ∈ H
Objective: To every transit node u is associated a due date d u , standing for the time at which the next road portion becomes unsafe. The objective is to minimize the maximum lateness of any task, i.e., the difference between the time at which it leaves a node u and d u , weighted by the population. Hence the objective is: min max u∈T ,v∈L(u)

s v + h v /w v + l uv -d u
Dominance rules and problem formulation A first observation is that we can simplify the objective function by retaining the transit node u minimizing d u -l uv . Let d v = min u∈T {d u -l uv } denote this value.

The objective can therefore be rewritten as the maximum of |E| expressions:

min max v∈E s v + h v /w v -d v
The second following observation, also made in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF], allows to reduce the number of cumulative constraints to consider. Observations 1. For any two transit nodes u, u ∈ T , if q u ≤ q u and u ∈ Ĉ(u), then a jam in u entails a jam in u.

We thus only need to check jams in nodes u such that ∀v ∈ p(u), q v > q u . In practice, it means that given an evacuation node v ∈ E we can explore the nodes in the route from v to r in reverse order, and for every stretch of road without branch, keep only the arc of minimal capacity, called the critical arc. Let T denote the reduced set of critical transit nodes to consider. The problem can thus be formulated as follows. minimize max

v∈E w v {s v + h v /w v -d v } (1)
subject to:

v∈L(u),sv+luv≤t<sv+luv+ wv hv h v ≤ q u , ∀u ∈ T , ∀t ∈ H (2)
h v ∈]0, q v ], s v ∈ [0, H - w v q v ], ∀v ∈ E (3)
3 Baseline approach using cumulative constraints

In [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF], the NEPP was modeled using standard cumulative constraints. We consider this model for our problem as the baseline approach. Let x (resp. x) denote the largest (resp. smallest) value in the domain of a variable x. Given a set of tasks J with start time variable s i ∈ [s i , s i ], processing time variable

p i ∈ [p i , p i ], height variable h i ∈ [h i , h i ]
and a resource r of constant capacity q r , recall that cumulative((s i , p i , h i ) i∈J , q r ) enforces the relations

i∈J|si≤t<si+pi h i ≤ q r ∀t ∈ H
Consequently, to model the problem, it suffices to associate a task v to each evacuation node v ∈ E, with height variable h v ∈]0, q v ], start time variable s v ∈ [0, H -wv qv ], completion time variable e v ∈ [ wv qv , H], and to duplicate and translate this task for each critical transit node u on its path towards the safe node. For each critical arc u ∈ T , let i uv denote the duplicate for evacuation node v ∈ L(u). A resource is defined per critical arc u ∈ T , with capacity q u .

The baseline constraint program for the evacuation planning problem is obtained by replacing constraints (2) in the problem formulation by: cumulative((s iuv , e iuv -s iuv , h iuv ) v∈L(u) , q u ) ∀u ∈ T ( 4)

w v = h v (e v -s v ) ∀v ∈ E (5)
s iuv = s v + l uv ∀u ∈ T , ∀v ∈ L(u) (6)
e iuv = e v + l uv ∀u ∈ T , ∀v ∈ L(u) (7)
h iuv = h v ∀u ∈ T , ∀v ∈ L(u) (8)
In the model above, the information that tasks have a fixed energy is lost to the cumulative constraint propagator. Given a task v with energy w v start time s v , completion time e v and rate h v , the algorithm will consider a task of duration wv h v and height wv e v -s v . When the upper bounds on consumption and duration are large, these values tend to 0, thus greatly hindering constraint propagation.

Example 1. Consider a resource of capacity 4 and four tasks as shown in Table 1. Given the total energy (second column) possible ranges for the rates, minimum start times and maximum completion times (3rd to 5th columns), bounds consistency on the duration and demand variables yields the ranges shown in the 6th and 7th columns, respectively.

Therefore, a standard cumulative constraint will use the lower bounds and consider four tasks of durations and consumptions 2, 3 2 , 1 and 1, respectively. The classic cumulative constraint will not adjust the domains further, as shown in Figure 1. The set of 3 solutions illustrated in this figure shows that every

w i q i s i e i duration height v 1 6 3 0 3 [2, 3] [2, 3] v 2 6 4 1 5 [ 3 2 , 4] [ 3 2 , 4] v 3 4 4 1 5 [1, 4] [1, 4] v 4 4 4 1 5 [1, 4] [1, 4] 0 1 2 3 4 5 v 1 v 3 v 2 v 4
(a)

0 1 2 3 4 5 v 1 v 3 v 2 v 4 (b) 0 1 2 3 4 5 v 1 v 3 v 2 v 4 (c)

The Energetic Cumulative Constraint

The standard cumulative constraint is weak for this class of problems since it is unable to reason about total energy of each task. We can take advantage of the knowledge that the product (duration × rate) is a constant by using a more global constraint.

Given a set of tasks J with start time

s i ∈ [s i , s i ], completion time e i ∈ [e i , e i ], height h i ∈ [h i , h i ],
constant energy w i and a resource r of constant capacity q r , energetic cumulative((s i , e i , h i , w v ) i∈J , q r ) enforces the relations:

i∈J|si≤t≤ei h i ≤ q r ∀t ∈ H w i = h i (e i -s i ) ∀i ∈ J
This constraint is NP-hard since the particular case where the variables h i are fixed is a cumulative constraint. However, because of the malleable nature of the tasks, the problem becomes easier under some assumptions. In particular, we consider here that we do not have minimum values on the height variable h i of any task i, and we denote this particular case energetic cumulative\{h}. We consider three further relaxations:

Let "r" denote the release dates, whose relaxation means that s i = 0 for every task i ∈ J; "d" denote the due dates whose relaxation implies e i = H for every task i ∈ J and some constant H; "h" denote the maximum height whose relaxation implies h i = q r for every task i ∈ J. We write energetic cumulative\{S} for the case where the subset of constraints S ⊆ {h, h, r, d} are relaxed.

Theorem 1. energetic cumulative\{h, x, y} is in P for any x = y ∈ {r, d, h} sketch. Because of the lack of space, we give only the algorithm for energetic cumulative\{h, r, d} as the other two are trivial: without bound on the height, we can vertically slice the tasks as thin as possible while satisfying the capacity q r and stack them w.r.t. to their release or due dates (depending on which was relaxed). Algorithm 1 schedules the tasks so that every tasks starts at time 0 and the latest completion time is minimum while satisfying h i ≤ h i ∀i ∈ J Thus, the constraint is satisfiable iff this latest completion time is less than H.

Algorithm 1: UsageSorted W.L.O.G, let J = {1, . . . , n} be such that i < j =⇒ w i q i ≥ w j q j ; K ← n i=1 w i ; Q ← qr; foreach i ∈ [1, n] do s i ← 0; h i ← min(q i , Qw i K ); K ← K -w i ; Q ← Q -h i ;
When Algorithm 1 ends the completion time of a task i is wi hi . Moreover, it can be shown that i < j =⇒ wi hi ≥ wj hj (we omit that part of the proof). Therefore 1 has the latest completion time. There are two cases, either this completion time is w1 q1 (in which case it cannot be reduced) or it is n i=1 wi qr (in which case the resource is fully used).

Theorem 2. energetic cumulative\{h, x} is NP-complete for any x ∈ {r, d, h} Proof. Membership: checking the energy equation is trivial and checking the resource equation is not more difficult than for the cumulative constraint.

We first prove hardness of energetic cumulative\{hd} by reduction π from an instance x = {x 1 , . . . , x n } of the SubsetSum problem of deciding the following proposition:

∃I | i∈I x i = s. Let K = n i=1 x i . For every i ∈ [1, n]
we create a task i with w i = x i , q i = 1 and r i = 0. Then, we create three tasks L, M and R with q L = w L = 2Ks, q R = w R = s + K, q M = w M = 2K and r L = 0, r M = 1, r R = 2. Finally, we set q = 2K and we ask if all the tasks can be scheduled with H ≤ 3.

We show first that a solution of x entails a solution of π(x). We set s L = 0, s M = 1, s R = 2, h L = w L , h M = w M and h R = w R . Now, let I be the subset of {1, . . . , n} such that i∈I x i = s. For all i ∈ I, we set s i = 0 and for all i ∈ {1, . . . , n} \ I we set s i = 2 and for all v ∈ {1, . . . , n} we set h i = 1. Now suppose that there is a solution of π(x). Observe first that in any solution, we have s R = 2 and h R = w R because s R ≥ 2 and h R ≤ w R and any solution such that s R > 2 or h R < w R implies e R > 3. Now, suppose that h M < w M , since e R ≤ 3 then M and R must overlap, and therefore

h M ≤ q -w R = 2K -(s + K) = K -s. So we have p M = 2K K-s > 1.
By contradiction, we have h M = w M and therefore s M = 1. By a similar reasoning we have s R = 0 and h R = w R . The remaining tasks are therefore distributed so that they exactly fill two disjoint areas of size s and K -s. The sum of the energies of the tasks in the leftmost area is therefore equal to s.

The proof is the same for energetic cumulative\{hr} except that the constraints r L = 0, r M = 1 and r R = 2 are replaced by d L = 1, d M = 2 and d R = 3, and for every other task i, we have

d i = 3.
The proof is the same for energetic cumulative\{hh} except that the constraints ∀i, q i = 1 are relaxed, the constraints q L = w L , q R = w R , q M = w M are replaced by d L = 1, d M = 2 and d R = 3, and for 1 ≤ i ≤ n, we have

d i = 3.

Flow Global Constraints

Note that a global constraint similar to energetic cumulative, called the continuous energy constraint, has been proposed in [START_REF] Nattaf | Cumulative scheduling with variable task profiles and concave piecewise linear processing rate functions[END_REF]. In this variant the rate of the activity is no longer constant and may vary over time. A flow-based propagator proposed for this constraint can be adapted for the energetic cumulative. The propagator works by building a flow network relaxation f (D) of the problem, before it propagates, using the current domain D.

Given current domain D, let T J = {s i , e i , s i , e i | i ∈ J} be the set of O(|J|) minimum or maximum start and completion times of the tasks of J and let t q be the q-th largest element of T J . We partition the time line into the set of consecutive intervals

I(J) = {[t q , t q+1 ] | 1 ≤ q < |T J |}
Next, we create a flow f (D) network as follows. We have source node S and a first layer of task nodes J, the flow from S to i is w i for each i ∈ J. We create a second layer of time interval nodes I q = [t q , t q+1 ], 1 ≤ q ≤ |J| -1. There is an edge from i to

I q if [t q , t q+1 ] ⊆ [s i , e i ].
The edge is bounded by 0..(t q+1 -t q )×h i . If, in addition, [t q , t q+1 ] ⊆ [s i , e i ], the lower bound of the edge can be increased to h i . We create a final layer to a sink node E. There is an edge from each interval node I q to E with capacity bounded by (t q+1 -t q ) × q r . Theorem 3. Any solution to energetic cumulative((s i , e i , h i , w i ) i∈J , q r ) given current domain D, is extendible to a solution of the flow network f (D).

Proof. Given a solution of the constraint we extend it to a solution of the flow network f (D) as follows. The flow from S to each node i is set to w i . The flow from each node i to I q is set to (min(e i , t q+1 ) -max(s i , t q )) × h i The flow from each node I q to E is given by the sum of the incoming flows to I q . We show that these flows obey the bounds and flow balance equations. Examining each node i, the flow in is w i and the flow out is h i ×(e i -s i ). These must be equal by equation [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF]. Examining each node I q the flow in is i∈E ((min(e i , t q+1 )-max(s i , t q ))×h i but by equation ( 4) at no time is there more than q r resource being used, hence this is no more than (t q+1 -t q ) × q r , the capacity of the outgoing arc. The flow balance at I q holds by construction.

Example 2. If we consider the circumstances explained in Example 1, the constraint energetic cumulative generates the flow network shown in Figure 2.

While the standard propagator can determine nothing, the flow network is infeasible, so the energetic cumulative propagator immediately fails. 1.

S v 2 v 1 v 3 v 4 [0, 1] [1, 3] [3, 6]
The natural way to implement this propagation for the evacuation planning problem at hand is to define one energetic cumulative constraint per critical transit node u ∈ T . The special structure of the evacuation planning problem allows us to consider simultaneously all energetic cumulative constraints and to integrate the flow propagators of each critical transit nodes into a single global flow propagator. The basic remark is that all vehicles passing in an interval [t q , t q+1 ] on a transit node u of an evacuation path will appear on a ascendant transit node u ∈ p(u) exactly l uu time units later, i.e. in interval [t q + l uu , t q+1 + l uu ]. Considering the evacuation tree reduced to the evacuation nodes and the critical transit node, the global network flow gf (D) is built as follows. We create the source S and one layer of nodes for the evacuation nodes E, the flow from S to v ∈ E being w v . Let T v = {s v , e v , s v , e v } the time events associated to evacuation node v ∈ E and let

I(v) = {[t v q , t v q+1 ] | 1 ≤ q < |T v |}
the corresponding consecutive intervals. We connect node v to each interval I v q ∈ I(v) by an edge of capacity bounded by 0..h v × (t v q+1 -t v q ) (or with a lower bound of h v under the conditions described above). We assume that the critical nodes are sorted in the topological order from direct successors of the evacuation nodes to the direct predecessors of the safe nodes. Then, for each critical node u ∈ T taken in this order, we build recursively the set of events T u from the events of its parent nodes in the evacuation tree and the corresponding consecutive interval set I(u):

T u = ∪ u ∈p(u) {t u q + l uu | 1 ≤ q < |T u |} and I(u) = {[t u q , t u q+1 ] | 1 ≤ q < |T u |}.
We add an edge between each interval of the parent node I u q = [t u q , t u q +1 ] ∈ I(u ) and each interval of the child node I u q = [t u q , t u q+1 ] ∈ I(u) such that t u q + l uu ≤ t u q and t u q +1 + l uu ≥ t u q+1 with a capacity bounded by 0..(t u q+1 -t u q ) × q u . Finally, we build an edge of unlimited capacity from each interval of each parent node u ∈ p(r) of the safe node to the target E. The global flow has O(|E| + u∈ T L(u)) nodes. The following theorem shows that is there is no feasible global flow, then energetic cumulative fails.

Theorem 4. Any solution to the conjunction of energetic cumulative((s iuv , e iuv , h iuv , w v ) v∈L(u) , q u ) for all u ∈ T given current domain D, is extendible to a solution of the flow network gf (D).

Proof. The proof uses similar arguments as in the proof of Theorem 3.

Heuristic Upper Bound

We propose a simple compression heuristic to find an initial upper bound. Assume w.l.o.g. that the evacuation rate initial domain is D(h v ) = [1, q v ] for all evacuation node The heuristic is based on the assumption that scheduling all evacuation tasks at time 0 with the minimum evacuation rate yields a feasible solution, with a high cost. Starting from this solution (∀v ∈ E, s v := 0, e v = w v , h v := 1), the set of minimum and maximum start end end time events T u is built for each transit node, as well as the corresponding consecutive intervals I(u) as for the global flow constraint. In addition, the resource consumption profile Γ u where ρ qu ∈ Γ u denotes the resource consumption level right after time t u q is also computed for each transit node u. At each iteration, the heuristic finds the critical task v * , i.e. the one that sets the objective, as well as the second most critical task v . The principle is to left shift the end of the critical task enough so that v becomes the new critical task with an absolute gap of below the new cost of v * . This yields a target completion time e * := c v - wv * -d v * for v * . However a left shift of e v * means an increase of its height. On a transit node u, the evacuation task v ends at time e v * + l uv * . Let I u q + denote the interval starting with the evacuation completion time on u: t u q + = e v * + l uv * . Let I u q - denote the interval starting with the evacuation start time on u: t u q -= s v * +l uv * . The maximum height ∆ h increase is equal to ∆ h = min q=q -,••• ,q + -1 (q q -ρ qu ). Left shifting the completion time of the critical task v * on u to the beginning of the preceding interval t u q + -1 is possible if

ww t u q + -1 -s v * -l uv * ≤ h v * +∆ h . Otherwise,
the critical task v * can only be shifted to s v * + l uv * + ww h v * +∆ h ∈]t u q + -1 , t u q + ] and the maximum left shifted end time on transit node u has been found. In the first case, the process can be iterated by tentatively ending the task in interval t u q + -2 till he maximum left shifted end time on u is found. Repeating this check on each resource yields a minimum possible left shifted completion time e. If e ≤ e * , then the left shift to e * is possible and v becomes the new critical task. The compression process restarts with v * := v . Otherwise, the left shit is only performed to e. v * is still critical and the heuristic stops. Due to the possibility of only left shifting a task by at each iteration, the heuristic is pseudo polynomial.

Generating a Realistic Data Set

Catastrophic wildfire requiring large population evacuation are, thankfully, rare events. However, it means that obtaining useful data is difficult, and indeed this a key problem within the GEO-SAFE project. A significant part of the project revolves around simulation tools such as EXODUS [START_REF] Galea | The EXODUS Model[END_REF], however, even simulated data was hard to come by.

Therefore, we opted for taking advantage of the project environment to contribute to this effort by generating our own "realistic" dataset. On the one hand, this approach may introduce biases since we must use models to generate realistic road networks and simulate wildfires. On the other hand, we believe that it will make it much more convenient for benchmarking algorithms in the future. As it turns out, the generated instance are challenging even though relatively modest in size, thus being interesting from an academic viewpoint as well.

Generation of road networks

The first step is to generate a graph standing for the road network. To this end, we used the quadtree model described in [START_REF] Eisenstat | Random road networks: the quadtree model[END_REF]. In a nutshell, this model starts with a single square formed by four nodes and four edges. At each iteration, a square is chosen and five nodes are added, one in the center of the square, and one on each edge connected by a perpendicular edge to the center node. A parameter r controls the sprawl, that is, the preference for splitting larger squares (r < 1) or smaller squares (r > 1). The graphs generated in this manner share many features with real road networks: they are planar, embedded in an Euclidean plane, have similar density distributions, path lengths are within a constant factor of the Euclidean distance, and the number of turns is logarithmic with high probability. An example of random quadtree network is shown in Figure 3a. The colors on the edges correspond to road capacity. To allocate capacities, we first compute a minimum Steiner tree spanning three randomly chosen nodes in high density areas ("cities") and connect these cities to the nearest corner of the outer square. The corresponding set of edges are given the highest capacity and are coloured in blue in Figure 3a. A second set of edges, forming a grid are given an intermediate capacity, they are coloured in green.

Simulating wildfire

The second step consists in determining safety due dates for every edge of the evacuation tree, that is, a time after which the edge become unsafe. To this purpose we use a relatively simple fire propagation model. We chose to use a simple model based on two parameters: a constant intensity γ representing the type of fuel material as well as the temperature, and a wind direction. Indeed, the goal is not so much to accurately predict fire propagation, as it is to generate safety due dates consistent with a wildfire. Of course, should the authorities use this type of planning tools during a real event, then correctly predicting fire propagation would be among the most important factors. The land area is discretized into squares of fixed size (we use another parameter to control this size) which can be in three states: untouched, burning and burned. The fire starts as a single burning square, then at each iteration, any untouched square adjacent of a burning square catches fire with probability γ( π-A π ) 2 , where γ stands for the intensity of the fire, and A is the angle between the wind an a vector going from the center of burning square to the center of the untouched square. Moreover, any burning square that did not propagate stop burning with probability γ 2 . Figure 3b illustrate the state of the simulated wildfire, with burning squares in red and burned squares in black.

Generating evacuation plans

The third step consists in generating the actual evacuation plan, that is, an embedded tree connecting a set of evacuation nodes E to a safe node r. Here again, the goal is not to compute the best evacuation plans, however they must be representative of what would be actual plans.

We first randomly pick a predefined number of evacuation nodes among the nodes of the graph that are in the state burned or burning of the simulated fire. Then we use the convention that the safe zone is the furthest corner from the center of the fire. The evacuation tree is computed simply by using a shortest paths algorithm, however with respect to an arc labeling taking into account first the safety due date of the arc, and only then its length and its capacity.

At this point we have all the information we need to define a fire evacuation problem as defined in Section 2. However, we use Observation 1 to remove redundant arcs. For every evacuation path, we explore the arcs from the safe node r to an evacuation node v. For any section of the path on which nodes have single child, all the nodes of that section have the same set of descendant leaves in the tree. Therefore the same set of evacuees will go through these nodes with the same relative delay. Therefore, we conserve only the arc of minimum capacity in that section, and only if this capacity is strictly smaller than that of the previous section. It follows that an instance with n evacuation tasks has at most O(n) "jam" constraints. Moreover, the objective can be stated as the maximum of O(n) expressions as shown in Section 2.

The tools we developed as well as the benchmarks instances we used in this paper can be accessed at anonymizedfortheblindreviewprocess.

Experimental Results

We generated 240 benchmark instances following the protocol described in Section 7. They are organized into three types of road networks: Dense, Medium and Sparse where the density refers to the number of intersections (respectively 400, 800 and 1200) in the land area. Notice that the graph has always 4 edges per node, so this corresponds to graph size. The impact on the instance is that larger graphs allow for more choices for the shortest paths and therefore longer independent paths. For every type of road network, we generated 4 classes of instances, with respectively 10, 15, 20 and 25 evacuation nodes. Finally, for every class we simulated 20 random wildfires and the subsequent evacuation trees. We used CPLEX to solve the flow formulation, which turned out to be extremely costly. Therefore we did not systematically call the propagator at each node. Instead we implemented a heuristic method to decide whether we should solve the flow using two criteria:

The first criterion that we use is the total load Ω(r) of a resource r or capacity q r on the tasks J:

Ω(r) = i∈J w i q r (max i∈J e i -min i∈J s i )
We only call the flow when max u∈ T Ω(u) ≥ 0.95, this value was empirically chosen. Moreover, observe that when Ω(r) > 1 the constraint is trivially infeasible. Therefore we can return an early failure in this case.

The second criterion is the size of the remaining search space. Indeed solving the flow might not be worthwhile even it is infeasible, if brute-force search would have been faster. Here we empirically chose 2 50 as minimal search space, defined as the product of the ranges of all start time and height variables. That is, we solve the flow only when:

i∈E log |D(s i )| + log |D(h i )| ≥ 50.
However, even with this approach, using the global flow constraint made CP Optimizer about one order of magnitude slower (measured by number of fails per second). Therefore we also tested a version that never solves the flow, and only fails when the overall load is too large for the resource's capacity (Ω(r) > 1).

We ran every method on every instance of the dataset with a time limit of 45 minutes on 4 cluster nodes, each with 35 Intel Xeon CPU E5-2695 v4 2.10GHz cores running Linux Ubuntu 16.04.4. Table 2 show results of running CP Optimizer with the strategy Depth first search (DFS). We use ' h ' to denote the use of the heuristic upper bound described in Section 6, ' + ' to denote the use of the overall load check described above, and ' f ' for the flow propagator. The columns "#s", gives the ratio of instances for which the method found a feasible solution, the column "ub" the mean upper bound, and we highlight the best outcomes using colors. We can see that DFS does not find a feasible solution in every case, whereas our upper bound heuristic always finds one. Moreover, starting from a good quality solution pays off and the best solutions are significantly better than that of the baseline method. However, the flow propagator (DFS f ) and even the simple load checking rule (DFS + ) do not help and in fact decrease the performance of CP Optimizer.

In fact, the results are even more negative when using CP Optimizer's default strategy. Table 3 show the results in this setting, using the same conventions. The number of feasible solution is here replaces by the number of optimality proofs, since a feasible solution was found in every case. Here we can see that again the extra propagation does not improve the results. Moreover, the impact of the heuristic upper bound is milder: it helps find better solution in some cases, but overall it is actually slightly worse than not using it.

We find the empirical results puzzling. In the case of the flow propagator, we expect a tradeoff, not necessarily positive, between the extra pruning and the computational overhead. However, it is much more surprising that the load 4 shows the average number of nodes required to solve an instance, the average number of fails due to the flow, the ratio of fails per node and the "speed" in number of branches explored per second. On the one hand, we observe that the flow propagator decreases the speed by a factor 6, while failing in only about 1% of the calls. Recall, however, that the flow is solved only when the remaining search space is large enough, so these fails are likely to cut a relatively large subtree. On the other hand, we can see that the impact of the load checker on the search speed is rather limited, although visible, and that is does causes early failures (though admittedly rarely).

Conclusion

We have introduced a variant of the flood evacuation planning problem, where the objective function is adapted to the context of wildfires. We have designed a generator of challenging and realistic instances for this problem. We have introduced a new type of cumulative constraint where tasks have flexible duration and height, but a constant energy. We have analysed the complexity of this constraint, and proposed several solution approaches, including an effective greedy upper bound and a global propagation method.

However, the experimental results we carried out suggests that these techniques tend to degrade CP Optimizer's performances. One possible insight that we can draw from these experiments is that the amount of tuning, sophisticated techniques, and the overall black-box nature of CP Optimizer, although making it extremely efficient in practice is also a problem when trying to test preliminary ideas on how to best tackle a problem. In order to improve over the results of CP Optimizer, it is therefore probable that we will need to use more refined algorithms and in particular dedicated flow algorithms which would be incremental and from which we could propagate bounds.

Figure 1 :

 1 Figure 1: Some feasible schedules

Figure 2 :

 2 Figure2: Flow network for the tasks given in Table1.

  (a) road network (b) simulated wildfire (c) evacuation plan

Figure 3 :

 3 Figure 3: An example of generated instance

Table 1 :

 1 Task parameters

Table 2 :

 2 Depth first search: upper bound and optimality ratio

			DFS	DFS h		DFS +	DFS h+		DFS f	DFS hf
		#s	ub	#s	ub	#s	ub	#s	ub	#s	ub	#s	ub
	dense 10 (20) 1.00	26469 1.00 26469 1.00	25485 1.00 25573 1.00	25485 1.00 25573
	dense 15 (20) 1.00 179099 1.00 176391 1.00 189864 1.00 188238 1.00 194737 1.00 195631
	dense 20 (20) 1.00 365596 1.00 367030 1.00 406828 1.00 407116 1.00 486971 1.00 447924
	dense 25 (20) 0.95 697661 1.00 663314 0.85 1093784 1.00 759437 0.85 1201755 1.00 881494
	medium 10 (20) 1.00	49294 1.00 49294 1.00	49567 1.00 49567 1.00	49567 1.00 49567
	medium 15 (20) 1.00 148513 1.00 143999 1.00 153372 1.00 149243 1.00 169840 1.00 170785
	medium 20 (20) 0.95 339103 1.00 339370 1.00 376395 1.00 384381 1.00 507420 1.00 427572
	medium 25 (20) 0.90 1291811 1.00 659269 0.80 1461162 1.00 722167 0.75 1286011 1.00 790511
	sparse 10 (20) 1.00	4474 1.00	4759 1.00	4464 1.00	5951 1.00	4464 1.00	5951
	sparse 15 (20) 1.00 131229 1.00 129276 1.00 135192 1.00 134420 1.00 146487 1.00 146917
	sparse 20 (20) 0.90 293930 1.00 312834 1.00 320080 1.00 323449 1.00 388194 1.00 364291
	sparse 25 (20) 0.95 598501 1.00 574534 0.80 593962 1.00 732975 0.80 684041 1.00 767851
	avg	(240) 0.97 333506 1.00 287212 0.95 369874 1.00 323543 0.95 395298 1.00 356172

Table 3 :

 3 Default search: upper bound and optimality ratio (20) 311676 0.45 312501 0.45 312697 0.15 313114 0.15 317035 0.10 316596 0.00 dense 25 (20) 531016 0.00 529889 0.00 530749 0.00 530975 0.00 539192 0.00 564722 0.00 medium 10 (20) 48591 1.00 48591 0.80 48591 0.95 48591 0.75 48776 0.90 48979 0.70 medium 15 (20) 124921 0.70 125195 0.50 125559 0.60 125110 0.60 126011 0.40 127421 0.40 medium 20 (20) 276101 0.25 276094 0.25 277545 0.10 277485 0.10 279654 0.05 278723 0.10 medium 25 (20) 488282 0.10 488555 0.10 491794 0.00 491423 0.00 493682 0.00 507275 0.00 (20) 120173 0.65 120151 0.65 120162 0.60 119803 0.55 120895 0.55 120576 0.55 sparse 20 (20) 235771 0.35 253807 0.35 237503 0.15 236920 0.20 238300 0.20 240112 0.15 sparse 25 (20) 437618 0.00 437219 0.05 439699 0.00 438160 0.00 441707 0.00 448979 0.00 avg (240) 230191 0.51 231762 0.47 231188 0.42 230991 0.39 233102 0.37 237175 0.33

		CPO		CPO h		CPO +		CPO h+	CPO f		CPO hf
		ub	opt	ub	opt	ub	opt	ub	opt	ub	opt	ub	opt
	dense 10 (20)	23294 0.95 23567 0.80 23354 0.95 23681 0.85 23503 0.90 23807 0.80
	dense 15 (20) 161100 0.70 161700 0.75 162856 0.60 162751 0.55 164524 0.45 165025 0.45
	dense 20 sparse 10 (20)	3747 1.00	3880 0.90	3747 1.00	3880 0.90	3946 0.85	3880 0.85
	sparse 15

Table 4 :

 4 Ratio of failures due to the propagator heuristic, which are both almost computationally free do not help at all, let alone hinders CP Optimizer. Table

	method	#sol	#nodes (total)	#fails (flow)	fails/node	speed
			avg	avg	avg	avg
	F/LN S	1.00	525793	2371	0.01161	5324
	O/LN S	1.00	3042714	2376	0.00137	30157
	LN S	1.00	3084572	0	0.00000	34508
	checker and the

Tasks v 4 and v 3 are symmetric so we omit the bounds for v 4

Acknowledgement

This work is partially funded by the H2020-MSCA-RISE-2015 European project GEO-SAFE (id 691161)


Christian Artigues

Emmanuel Hébrard

Yannick Pencolé

Andreas Schutt

Peter J Stuckey

Geosafe

GeoSafe -Evacuation planning problems

published or not. The documents may come

Introduction

This paper aims at identifying new challenging variants of the evacuation planning problems, especially in case of bush fires, based on exchanges with practitioners in the context of the GeoSafe project [START_REF]Geo-safe -geospatial based environment for optimisation systems addressing fire emergencies[END_REF]. A large amount of work has been carried out at NICTA [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF][START_REF] Even | NICTA evacuation planner: Actionable evacuation plans with contraflows[END_REF][START_REF] Pillac | A pathgeneration matheuristic for large scale evacuation planning[END_REF][START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF][START_REF] Even | Convergent plans for large-scale evacuations[END_REF][START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF][START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF][START_REF] Romanski | Benders decomposition for large-scale prescriptive evacuations[END_REF] mainly in the context of floodings, which can be easily transposed to evacuation in case of fires. The recent survey by Bayram [START_REF] Bayram | Optimization models for large scale network evacuation planning and management: A literature review[END_REF] presents the state-of-the-art on large scale network evacuation planning and management.

Basic evacuation data

We adopt the notation an terminology given in [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF]. There is a directed graph G(N = E ∪ T ∪ S, A) representing e.g. a road network in a region that must be evacuated. The graph is made of:

• the set of evacuation nodes E. An evacuation node represents a zone where people to be evacuated are regrouped),

• the set of safe nodes S. A safe node represent a safe geographical zone that people located in the evacuation nodes must reach during the planning horizon

• the set of transit nodes T . A transit node represent an intersection in the road network that can be traversed by the vehicles carrying evacuated people from the evacuation zones to the safe zones.

Suppose the time is discretized in minutes. Each evacuation node k ∈ E is associated with a number of evacuees d k . Each arc q ∈ A has a capacity u e in evacuees/minute, a travel time t e giving the number of minutes that a vehicle takes to traverse the arc and a end time b e , which gives the time at which the arc becomes unavailable due to the fire propagation. Other characteristics may appear in the variants considered in the cited papers.

Previously considered evacuation problems

We first reviews the NICTA papers. The first considered problem and solution methods were presented by Pillac et al. in 2013 in a research report that was later published in [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF]. One of the practical motivation of the work was that in an urgency situation, the possibility of choices in a road network, such as a fork, generate stress among the population. This is why it is preferable to define a single evacuation path for each evacuation zone. Hence the problem considered in [START_REF] Pillac | A conflictbased path-generation heuristic for evacuation planning[END_REF] can be described as a two-level optimization problem, which can be solved in an integrated way or by a two-phase approach. At the first level, exactly one evacuation path is determined from each evacuation node to a single safe node. At the second level the flow of evacuees is scheduled trough the path. A MILP (called the restricted flow model RF) is proposed to solve the problem in an integrated way. It is based on time-expanded graph, in which each node is duplicated for each time period. Arc capacities u e in the time expanded network ensure that the traffic does not exceed the road capacity, especially when several paths use the same arc. The model includes continuous flow variable per evacuation node per arc and in the expanded network and also a binary variable per edge in the (non-expanded) graph G to ensure the uniqueness of the path for each evacuation node. As the MILP is intractable on a set of instances (HN) derived from a real case scenario involving 70000 evacuees in the Hawkesbury-Nepean floodplain, located North-West of Sydney, an iterative column-generation-like two-phase heuristic is proposed. Given a set of potential paths, a master problem solves the path selection and flow scheduling problem with another MILP (called the conflict-based heuristic path generation master problem, CPG-MP) involving a reduced number of binary constraints (one per path). And a sub-problem finds new path based for a subset of critical nodes by solving a multiple-origin, multiple-destination shortest path problem. A lexicographic objective is considered. The main objective is to maximize the number of evacuated people during the time horizon, while a secondary objective aims at maximizing the evacuation start time, based on the practical relevance. Note that the second objective is only indirectly tackled by weighting each arc in the time-expanded graph by a cost inversely proportional to the arc time, hence solving a min-cost flow problem.

Subsequent papers describe variants of the models and new methods. A more sophisticated variant of the heuristic was proposed in Pillac et al. (2014) [START_REF] Pillac | A pathgeneration matheuristic for large scale evacuation planning[END_REF] for the same problem with, however, a different secondary objective for the min-cost flow, aiming at minimizing the evacuation time (called the clearance time), which is indirectly obtained by weighting the arc in the time-expanded graph by a cost proportional to the arc time.

In Even etl a. (2014) [START_REF] Even | NICTA evacuation planner: Actionable evacuation plans with contraflows[END_REF], an additional degree of freedom is introduced, giving the possibility of selecting contraflow roads, which consists in reversing the direction of some major roads. In practical cases, this possibility can highly increase the network capacity for evacuation. This is done by introducing a subset A c of arcs that can be used in contraflows and by modifying the CPG-MP, interestingly without the need of introducing additional binary variables. A software called the NICTA Evacuation Planner is also presented, with new instances with up to 1 000 000 evacuees, which are solved requiring up to 30 min. of CPU time.

In Even et al. (2015) [START_REF] Even | Convergent plans for large-scale evacuations[END_REF], the case where the selected paths must form a convergent subgraph, i.e. for which each node has at most one outgoing arc in the set of paths, is considered. This is motivated by the fact that convergent evacuation paths can be better controlled. Indeed, even if in the previous model, a single path was assigned to each evacuation node, in the case where two paths are merged and then are separated again, a driver can be confronted to a choice and take a decision that would not correspond to what was planed. Hence the set of paths now form a tree rooted at the safe node. The solution method is still a two-phase method but does not follow anymore the column generation principle. A tree is first built by a MIP working on an aggregated graph (without time discretization). The second phase is the flow scheduling problem given the computed tree, which is a maximum flow problem. The first phase is now able to produce good upper bounds on the maximum number of evacuees and on the minimum clearance time. A simulation study shows that the convergent model outperforms the general model when the presence of a fork in the network generated an hesitation for the driver that is capture by a 0.75 second delay.

This model is further developed by Kumar et al. (2016) [START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF] to incorporate network design aspects in the problem to model possible infrastructure enhancement decisions, as for the west Sydney case concerns about the capacity of the road network were expressed. Two additional decisions are introduced : adding lanes to a road (arc) and elevating a road (arc). The first decision results in an increase of the capacity of the arc while the second decision postpones the flooding time, both coming at a cost. These decisions can be incorporated in the tree design problem, which is also the master problem of a Benders decomposition method. The objective is still to maximize the number of evacuated people under a budget constraint that limits the infrastructure upgrades. The maximum flow subproblem is used to obtain a feasible flow schedule and also to generate optimality cuts that are reinjected in the master problem. To increase the performance of the Benders method, pareto-optimal cuts are generated.

All previous approaches assumed that each individual vehicle of an evacuation (flow unit) moves independently from the other vehicles of the same zone, and exactly as prescribed by the maximum flow model. Pillac et al. [START_REF] Pillac | A columngeneration approach for joint mobilization and evacuation planning[END_REF] propose to use the concept of response curves to incorporate behavioral models in the problem. In practice, instead of assigning a start time to each evacuee, the authorities can influence the evacuation start time of a zone and the level of resources mobilized to increase the evacuation rate (e.g. number of agents knocking on people's door), to which people answer according to a behavioral model abstracted by a response curve. It follows that to each evacuation zone k ∈ E, a set of response curves F k is given. For each evacuated zone, a start time δ k and a response curve f f ∈ F k has to be assigned such that the flow of evacuees leaving zone k at time t is givn by

φ t k = 0 if t < δ f (t -δ) if t ≥ δ.
A column generation approach is proposed where the set of all evacuation plans of a zone is considered, where a evacuation plan is a path from the evacuated zone to a safe node, the start time and the response curve. The master problem selects a plan for each evacuation zone in such a way that the network capacity is satisfied and the total cost of the plans is minimized (without ensuring that a convergent is obtained). The plan generation subproblem is solved either via a MIP or via a shortest path approach. The methods experiences difficulties is solving realistic instances due to a longtail effect.

Another drawback of the flow model is that it generated preemptive evacuation plans. As flow units are routed independently there are time periods in which the evacuation of a zone may be stopped and reinitiated later. Even et al. (2015) [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] report that this creates serious implementation issues for the evacuation plans. They propose another model in which the evacuation rate is a decision variable that remains fixed as soon as the evacuation starts and as long as the zone is not fully evacuated or the time horizon is reached. They assume that the evacuation path of each evacuation zone is already determined and concentrate only on the scheduling problem. We present this model in the case that all people must be evacuated and the objective is to minimize the clearance time. Let δ k denote the start time of evacuation of zone k and let λ k denote the evacuation rate of zone k. Let A k = (e 1 k , . . . , e

|A k | k
) the evacuation path of zone k given as its list of edges. An evacuation plan defines a task J i k for each edge e i k with a start time S i k = δ k + i-1 q=1 t i k where t i k is the travel time of edge e i k . The evacuation has a variable duration p k and the total number of evacuees is p k λ k with p k λ k = d k since all people must be evacuated. Each edge task e i k ∈ A k has duration p k . Let J e the set of tasks that use edge e ∈ A. The model can be written as follows :

min C max (1)
C max ≥ S |A k | k + p k ∀k ∈ E (2)
S i k + p k ≤ b e i k ∀k ∈ E, ∀i = 1, . . . , |A k | (3)
p k λ k = d k ∀k ∈ E (4)
S i k = δ k + i-1 q=1 t i k ∀k ∈ E; ∀i = 1, . . . , |A k | ( 5
)
J i k ∈Je,S i k ≤t<S i k +p k λ k ≤ u e ∀e ∈ A, ∀t ≥ 0 (6)
This is a no-wait total work-and resource-constrained project scheduling problem where (4) are the total work constraints, ( 5) are the no-wait constraints and ( 6) are the cumulative constraints. The problem was efficiently solved via constraint programming for both the clearance time minimization version and the maximization of the number of evacuated people variant.

New problems

In this section we identify possibly relevant problems that were not addressed in the literature.

Integrated evacuation routing and non preemptive scheduling

The CP approach [START_REF] Even | A constraint programming approach for non-preemptive evacuation scheduling[END_REF] was limited to the flow scheduling part. Consequently, it could be integrated as the subproblem solver in the benders method proposed in [START_REF] Kumar | Optimizing infrastructure enhancements for evacuation planning[END_REF]. This would allow the integrated generation of convergent evacuation routes and non preemptive schedules. However the subproblem would become a combinatorial optimization problem, which raises the question of generating combinatorial benders cuts.

Objective function

The time b e at which arc e becomes unavailable comes from a fire/flood propagation model. In the case of fire, even if precise propagation models can be obtained, they depend on multiple parameters. Among them, the wind has a great variability. An alternative to modeling uncertainty of the unavailability dates would be to consider an objective function that seeks to maximize the length if the time interval [C e , b e ], where C e is the completion time of the last task using edge e. we can also considered the weighted case, which could allow to model the risk associated to a particular edge. Note that the objective would be to maximize e∈A w e (b e -C e ), which amounts to minimize the weighted resource disposal time e∈A w e C e . A question arise about the complexity of this objective for preemptive evacuation. The preemptive scheduling problem (probably) becomes NP-hard but that has to be checked, for particular networks. For non preemptive evacuation, it is not clear whether the CP method would remain efficient because of the objective change from a max to a sum. It could be interesting to design hybrid MILP/CP method.

Another question raised would we how to incorporate the smoke problem in the model constraints and objective. May we have two values on the arcs one for the burning time and another one for the presence of smoke that can be earlier ?

Uncertainty

We could introduce an explicit modeling of uncertainty either in the arc unavailability dates, and/or in the evacuation rate and/or in the numbers of evacuees. Discrete scenarios could be considered with the objective to find a robust evacuation plan. A complexity question, still for preemptive evacuation scheduling, is it easy to find the best schedule that maximizes the worst number of evacuated people (or minimizes the worst clearance time) on the set of discrete scenarios ? Can we propose an efficient robust optimization algorithm in the general non preemptive case with integrated routing ? In a two-phase method the uncertainty would be tackled in the sub-problem (to find a schedule feasible in the worst case). A certain number of approaches tacking uncertainty in evacuation planning are reported in [START_REF] Bayram | Optimization models for large scale network evacuation planning and management: A literature review[END_REF], namely robust optimization, stochastic optimization and chance constrained programming. Note that in case of uncertainty in the unavailability time of an arc, the chance constraint

P[C e ≤ be ] ≥ α ( 7
)
where be is the stochastic unavailability time variable for arc e and α is the desired guarantee has a deterministic expression provided that the inverse of cumulative distribution function F has an analytical expression. P[C e ≤ be ] = 1 -P [ be ≤ C e ] = 1 -F (C e ). It follows that ( 7) is equivalent to 1 -F (C e ) ≥ α and so

C e ≤ F -1 (1 -α) (8)
This means that the deterministic constraints also offers a probabilistic guarantee in case of a stochastic unavailability date.

Evacuation rates

In non preemptive evacuation schedule, a variable evacuation rate could be introduced. We can imagine a situation where the authorities may regulate the rate of the vehicles that exit the evacuation area. This would allow to find a better compromise between congestion and evacuation speed. The energy (or total work) constrained scheduling model, that allows to define for each task a variable rate from a minimum value λ min k > 0 to avoid preemption to a maximal value λ min k . The rate function is a decision variable. Note that considering such variable rates can be used when the decision maker has actually the control on the rate but also in the case that the rate is an uncertain parameter, whose lower and upper bounds are known. In that case constraint propagation on the total work can be used in a search tree to establish that there no feasible scenario according to the current search status. The cumulative constraint with variable resource usage model proposed in [START_REF] Nattaf | A hybrid exact method for a scheduling problem with a continuous resource and energy constraints[END_REF] and the associated constraint propagation algorithm could be used. Already established results that state that is is dominant to consider rate changes only at the start and completion time of tasks, will probably be transposable here.

Evacuation vehicle routing, location of evacuation and safe nodes

The models presented above implicitly considered that each evacuee uses its own vehicle but we could also consider that some evacuees need to be evacuated in shuttles available in a limited number and with a limited capacity that have to travel back and forth from the evacuation node to the safe nodes, (or even using more complex routes that have to be determined but this will conform no more to the hypothesis of a single path from an evacuation node to a safe node). Another decision level is to decide the location of evacuation and safe nodes. This is mentioned in the survey [START_REF] Bayram | Optimization models for large scale network evacuation planning and management: A literature review[END_REF] as the Shelter location problem also studied by several authors.

Data instances

A crucial issue is the generation/obtention of realistic data instances. The ones used in the NICTA papers are not available. On possibility is to use the SPARK toolkit developed at CSIRO to generate instances [START_REF] Spark | a wildfire simulation toolkit for researchers and experts in resilience field[END_REF].

Acknowledgement

This work is partially funded by the H2020-MSCA-RISE-2015 European project GEO-SAFE (id 691161)