Marc Demange
email: marc.demange@rmit.edu.au

Virginie Gabrel
email: virginie.gabrel@dauphine.fr

Mohamed Haddad

Cecile Murat
email: cecile.murat@dauphine.fr

Marcel A Haddad
email: marcel.haddad@dauphine.fr

Cécile Murat

M A Haddad

C Murat

A robust p-Center Problem under Pressure to locate Shelters in Wildfire Context

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the prevention phase of wildfires, an important issue is to determine the location of fire proof shelters in a given area. The problem is basically to locate p shelters minimizing the maximum distance people have to cover in order to reach one of these shelters in case of fire. From an Operations Research perspective, the problem of locating p centers (warehouses, shelters . . . ) has been extensively studied [?]. In particular, mathematical programming approaches enable to determine optimal solutions for the classical p-Center problem [?].

In our model, the territory, typically with low density habitat, is represented by an adjacency graph G = (V, E). Each node corresponds to a zone. Two nodes i and j are linked by an edge if and only if it is feasible to go directly from one zone to the other without passing through another area. We assume this is a symmetric relation, which makes this graph undirected. The zones are defined by geographical and spatial criteria. In particular, natural barriers like a river or a cliff are natural boundaries. Of course the size of the zones vary depending on the population density. Typical examples of zones are villages and surrounding suburbs or, in the case of very sparse habitat, a homogeneous area with easy interior circulation. We assume that these zones have been defined a priori. Each edge (i, j) is valued with a positive number l i j that can be seen as a distance or a travelling time. Since nodes represent zones, possibly large, the distance (or travelling time) between adjacent nodes can be measured between median points in the two related areas or as a maximum distance; this choice has no incidence on the model. If adjacent zones correspond to a fragmentation of a large homogeneous territory without clear natural boundaries, then the related edges are just a discrete model of a continuous reality. In other cases however, edges represent an accurate representation of the real environment. According to our model, shelters must be placed on nodes. Considering a shelter on node i and a zone represented by a node j, the distance d ji from j to i is the length of a shortest path from j to i in G. Given a solution defined as a set of p nodes, the evacuation distance of zone j is the shortest distance between j and one of the p shelters. The value of a solution corresponds to the longest evacuation distance of nodes. Thus the p-Center problem is to compute a minimum value solution.

A major difficulty in wildfire management is due to the uncertainty associated with fire outbreaks and spreading. The circumstances which interest us are when a wildfire is expected to spread rapidly and to be of extreme intensity. This phenomenon happens for instance, on hot and windy days with dense and very dry vegetation. Furthermore the way how the fire spreads is hard to predict, since it depends on meteorological conditions (drought, heat, winds..) and on the nature of the vegetation. The black Saturday disaster in Victoria, Australia in February 2009 and in particular the tragic events in Marysville, a small town devastated by flames, is a perfect illustration. With climate change, such fires are an increasing phenomenon and are the subject of the project GEO-SAFE (http://geosafe. lessonsonfire.eu/) under the European Union's H2020 research and innovation program. Our study is part of this project. Under these extreme circumstances, in case of a fire outbreak, the whole population of the area has to be evacuated to shelters and stay there until being rescued. This assumes an efficient early warning system and clear messages to evacuate. In this process, we need to take into account that some routes to shelters may no longer be practicable.

To address this challenge, we define a new version of the p-Center problem considering a set of scenarios. A scenario corresponds to a zone on fire meaning that the zone can no longer be reached. Moreover in the emergency context of wildfires, the evacuation decision is made under pressure and implies a specific evacuation strategy. The evacuation context involves new distances to reach shelters and the value of a solution is computed as the worst evacuation distance. This leads to a specific robust problem called the Robust p-Center problem under Pressure (RpCP).

This paper is organised as follows: in Section 2, we define RpCP and we compare it to the state of the art. In Section 3, we study the complexity of the problem in specific classes of graphs. In particular we establish the NP-hardness of this problem on subgraphs of grids. In Section 4, we propose a new formulation based on 0-1 Linear Programming. Finally in section 5, a new exact algorithm is defined and extensive experimental results show the tractability of our approach.

2 Robust p-Center problem under pressure In Subsection 2.1, we define and motivate the problem, while in Subsection 2.2, we underline the differences between solutions of RpCP and the classical p-Center problem. Finally in Subsection 2.3, we compare RpCP to relevant variants in the literature.

Definition

To model fire hazards, a scenario is associated with a specific fire outbreak in the territory that is represented by a graph G = (V, E) with n nodes. We restrict ourselves to single fire outbreak and consequently, each scenario s corresponds to a single node s on fire. This restriction is motivated by our primary focus on a relatively short time period after outbreak which assumes an efficient early warning system. In this case everybody can escape to a shelter before the fire spreads to adjacent zones. The operational graph associated with the scenario s, denoted by G s , is a directed graph obtained from G as follows. All edges (i, j) are replaced by two arcs (i, j) and ( j, i) except the edges (s, v) incident to s: these edges are replaced by a unique arc (s, v). Consequently, in G s , node s is no longer accessible from another node. In G s = (V, E s ), paths are directed and directed shortest paths from a node x = s in G s cannot go through s. For all i and j, d s i j denotes the length of the directed shortest path from i to j in the graph G s . By convention, for all j ∈ V \ {s}, we have d s js = +∞ which means there is no path from j to s.

Given a set of p nodes, solution to the p-Center problem, the evacuation distance of a node j is the shortest distance between j and its nearest shelter. This evacuation distance is not relevant in our context since such shortest path may no longer be practicable because of fire. For this reason, we consider a more realistic model to evaluate the evacuation distances to shelters. This model has been discussed with final users as part of the Geosafe project 1 and in particular in dedicated open sessions of GEO-SAFE Wildfire conference, (https: //geosafe.lessonsonfire.eu/fireconference/).

In case of a fire outbreak on node s, we have:

1. for people in zone s, two cases have to be considered.

(a) If a shelter is located in zone s, we assume the following: all the persons present in that zone can reach this shelter safely. So the evacuation distance is set to zero. To support this hypothesis, we assume that the shelter location in the zone and the local layout guarantee easy accessibility [?] during the phase after outbreak. It is even reasonable to assume that on and around the shelter, clear signs direct people accurately. However, for people outside the zone, to attempt to reach the shelter could be dangerous. For this reason s is not accessible from other nodes in the operational graph G s . (b) If there are no shelters in zone s, we assume that people in this zone will first flee in any other direction to reach an adjacent zone j, they will evacuate to the nearest shelter from j in G s . Thus, the evacuation distance associated with this zone is the worst out of all the possible adjacent nodes to reach first. Without strong evidence (like a shelter in the zone), people, especially under pressure, may react in very diverse ways: choosing for instance to run in the opposite direction of the fire, or in the opposite direction of the wind or in the most accessible direction. On reaching a safer zone, they are more likely to observe the information provided, as their own decisions have become more rational.

To ensure an acceptable level of risk, we must consider the worst case scenario. 2. for people who are not in zone s, for example in zone j = s, the evacuation distance from j to shelter k corresponds to d s jk in graph G s , i.e., avoiding node s. As already mentioned, we assume that all people in this territory need to be sheltered. Going to the closest accessible shelter is safe. It is important to remember that fires can be very unpredictable, very fast spreading, and can flare up from embers and at a distance. This examination of evacuation distances renders the problem specific compared to the literature and introduces some additional complexity.

For a given p, a solution C of RpCP corresponds to a set of p nodes where to locate the p shelters. For a given solution C and a given scenario s, the evacuation distance of zone j is denoted by r s j (C). If a shelter is located on j, r s j (C) = 0 otherwise we have:

r s j (C) =      min k∈C d s jk if j = s max v∈N + G s (s) {l sv + min k∈C d s vk } otherwise (1)
where N + G s (s) is the set of all nodes v such that (s, v) ∈ E s . Notice that r s j (C) is equal to +∞ if j can't reach any shelter in G s . In this case, the solution C is not feasible and consequently instances may occur without feasible solutions (see Subsection 2.2 for further details).

The radius associated with scenario s is defined as r s (C) = max j∈V r s j (C). A given solution C corresponds to n radius r 1 (C), . . . , r n (C) for n different scenarios. In the context of emergency evacuation, the value r(C) of a solution C is obtained by considering the worst (max) radius: r(C) = max{r 1 (C), . . . , r n (C)}. The problem RpCP is then to determine a solution C * with the minimum worst value:

r(C * ) = min C r(C) = min C max s∈V r s (C) = min C max s∈V max j∈V r s j (C)
Example 1 Let us consider the operational graph associated with scenario 2 in the Figure 1 where the two arcs (i, j) and ( j, i) are represented by a unique continuous line. The solution C = {3, 10} corresponding to shelters located on nodes 3 and 10 is represented by pentagon nodes. In case of fire on node 2 (scenario 2), the modification of the graph and the evacuation strategy induce: -The shortest path value from 1 to 3 is no longer 3 but 23, using the shortest path 1, 6, 7, 8, 3. Consequently, the nearest shelter from node 1 is 10 at a distance of 8. Thus the evacuation distance of 1 in scenario 2, is equal to 8 and node 1 is evacuated to node 10. -To compute the evacuation distance of node 2 in scenario 2, we have to consider three neighbours:

for neighbour 1, the distance to the nearest shelter 10 is 1 + 8 = 9; for neighbour 7, the distance to the nearest shelter 10 is 3 + 9 = 12; for neighbour 3 with a shelter, the distance is 2. Consequently, r 2 2 (C) = 12. -The radius of the scenario 2 is given by r 2 (C) = max j=1,...,14 r 2 j (C) = r 2 13 (C) = 15. Finally, the radius of this solution is defined by r(C) = max s=1,...,14 r s (C) = 35 induced by node 4 in scenario 3.

Comparison with the p-Center problem and feasibility condition

A solution C ⊂ V for RpCP is feasible if it has a finite evacuation distance for all nodes in all scenarios. In some cases, most of the feasible solutions of the p-Center problem are no longer feasible for RpCP. For example, consider the path G given in Fig 2 . When p = 2, any set of 2 nodes is a feasible solution for the 2-Center problem while the only feasible solution for the robust 2-Center problem is the set {1, 5} of the extremity nodes of the path. For all C without node 1 (or equivalently 5), we necessarily have r 2 1 (C) = +∞ (since node 1 cannot reach any shelter) and the solution C is not feasible. We observe that any feasible solution for RpCP must include all nodes of degree 1, called leaf nodes. Consequently, when p is smaller than the number of leaf nodes, RpCP admits no feasible solutions. Moreover, a solution C ⊂ V for RpCP has a finite evacuation distance for all nodes in all scenarios when, in all subgraphs G s , there is a path between any node v ∈ G s and a shelter. This implies the need for a shelter to be located in any connected components of G s for all scenarios s. Let us denote cc s the number of connected components of G s . When p is smaller than max s∈V cc s , RpCP admits no feasible solutions.

Even when any solution of the p-Center problem is feasible for RpCP, the relative error of using an optimal solution of the p-Center problem for RpCP can be arbitrarily large, as shown in the following. Consider the instance represented in Figure 3, for M ≥ 1. The optimal choice for the 2-Center problem is to locate shelters at nodes 2 and 3, for a radius of 1. The value of this solution in R2CP is M + 3, induced by the evacuation distance of node 1 to shelter 3 in scenario 2. However, if we locate the shelters at node 5 and 6, the worst radius value is 3 induced by scenarios 2, 3, 5 and 6. The relative error of using an optimal solution of the 2-Center problem for R2CP is (M+3)-3 3 = M 3 . To the best of our knowledge, RpCP including our specific evacuation strategy has never been studied so far. Some variants of the p-Center problem close to our model are proposed in the literature. In the next section, we present some of them and underline the main differences.

Comparison with the state of the art

The robust version of the p-Center problem is usually defined by introducing uncertainty of edge weights: each weight may vary in an uncertainty set (usually an interval) and the problem is to determine a solution minimizing the worst case or the maximum regret [?,?,?,?]. In [?], authors apply the Bertsimas and Sim approach [?] to model and to solve a p-Center problem with interval associated to edge weights. In these works, weights vary independently one from each other. In our context, this independence hypothesis is not relevant since, if a fire ignites on a node, all weights of the edges incident to this node are modified in the same way. In our context however, uncertainty has to be represented by a set of discrete scenarios. Such approach has already been studied for facility location problems [?,?]. These two papers review the literature on stochastic and robust facility location models. Most deal with scenario-based approaches on generalizations of the p-median facility location problem. For example, in [?], a p-median problem under scenario-based demand uncertainty is considered. In [?], a new robustness measure is introduced and used on a discrete scenario-based approach for the p-median problem. This measure consists in minimizing the expected cost under the constraint that the relative regret is bounded in each scenario. More recently, in [?], authors propose a stochastic evacuation planning model that optimally locates shelters and assigns evacuees to the nearest shelter so as to minimize the expected total evacuation time. The two reviews [?,?] underline that robust p-Center problems are often more difficult than the related p-median problems. This explains why scenario-based approaches were not considered for the p-Center problem until very recently. To our knowledge, the only paper studying the robust p-Center problem with a scenario approach is [?]. They propose a robust model for a reliable p-Center problem. Each scenario corresponds to a set of disrupted facilities and updated demands and costs. Clients are reallocated to the nearest surviving facility. In all these studies, different travel times on the edges are realised for different scenarios. Given a solution C, the radius of node j in scenario s is: min k∈C d s jk . This distance only corresponds to the case j = s in Equation 1 defining the evacuation distance r s j (C). In the other case ( j = s), the evacuation distance is no longer computed with a min operator over C but with a max min since max v∈N + G s (s) {l sv + min k∈C d s vk } = max v∈N + G s (s) min k∈C {l sv + d s vk }. This case is due to our specific evacuation strategy that differentiates our problem from the classical scenario-based approach.

In [?], Chaudhuri, Garg and Ravi consider the k-neighbour p-Center problem. It is a generalization of the p-Center problem where, given a number k, we have to place p centers so as to minimize the maximum distance of any non-center node to its k th closest center. They give a 2-approximation algorithm for this problem, and show it is the best possible. However, the evacuation strategy in RpCP does not correspond to an automatic reassignment to a k th center. For a given solution and for all scenarios s, we do not assign node v to some predetermined k th closest center in G, because, among other things, we have no guarantee that the k th closest center for v in G is accessible in G s . If we consider graph G s , while the nodes v ∈ V \ {s} are evacuated to their closest center in G s , the evacuation strategy of node s depends on its neighbourhood. Thus, there is no simple way to reduce RpCP to the k-neighbour p-Center problem.

In a variant of the p-Center problem for large-scale emergencies considered in [?], the disaster affects a single node s, including any facility (e.g. a shelter) on this node. The difference with our model is twofold. Any facility on an affected node is no longer available but only the population on this node requires evacuation. This model is motivated by different kind of disasters that affect a single node but also by the fact that each node corresponds to a large zone, like an entire city. Our context is really different since all zones must be evacuated in each scenario s and a shelter always secures at least the people from the corresponding area.

So, to the best of our knowledge, the RpCP that we introduce with specific evacuation strategy under pressure has never been studied. In the following section, we analyze the complexity of RpCP on different classes of graphs.

Complexity

RpCP is clearly polynomial (even thought not tractable in practice) for any constant p; so, in this section, we will consider that p is not fixed and thus, it is part of the input.

For this section, we consider the decision version of RpCP with constant threshold. For any constant k, RpCP k takes as input an integer k and an instance of RpCP; it is to decide whether there is a solution C with radius r(C) ≤ k. This problem is clearly in NP: we consider a polynomial number of scenarios and consequently, for any solution C, checking whether r(C) ≤ k can be done in polynomial time. For each scenario, it requires evaluating the evacuation distances of each node using a minimum path algorithm.

In this paper, we will address k = 1 or 2 and all edges of length 1. For k = 1, we outline a close relation between RpCP 1 and the problem of deciding whether a graph has a node cover of size p, i.e., a set of p nodes such that each edge of the graph is incident to at least one node of the set (Minimum Node Cover problem). Then, we consider a subclass of bipartite planar graphs, the class of induced subgraphs of a grid called subgrids, that is a realistic class of real instances. RpCP 1 is polynomially solvable in this class while RpCP 2 reveals to be NP-complete. We conclude with some comments about the complexity of RpCP.

RpCP 1 and the Minimum Node Cover problem

We consider an undirected graph G = (V, E), where all edges are of length 1.

The following elementary proposition leads immediately to a first hardness result:

Proposition 1 Let G be a graph. The radius of a set of nodes C ⊆ V is at most 1 if and only if C is a node cover that includes all leaf nodes.

Proof Suppose C a set of nodes with r(C) ≤ 1. Since its radius is finite, it should include all leaf nodes. Let u ∈ V \C, then for scenario u (u is on fire) and the set C, we have r u u (C) = 1, which means that all neighbours of u are in C. Thus, C is a node cover.

Conversely, if C is a node cover including all leaf nodes, we consider a scenario u and

a node v such that v / ∈ C (in particular v is not a leaf node). If v = u, since v / ∈ C, every neighbour of v is in C. If v = u
and since v is of degree at least 2, v has at least one neighbour in C \{u}. In both cases, the evacuation distance is 1 and consequently r u (C) = 1. Therefore, r(C) ≤ 1, which completes the proof.

In particular, if G has no leaf node, then r(C) is at most 1 if and only if C is a node cover. We deduce immediately:

Corollary 1 1. RpCP 1 is NP-complete in all classes of graphs of minimum degree 2 for which the decision version of Minimum Node Cover problem is NP-complete. The problem RpCP is NP-hard in these classes of graphs. 2. RpCP 1 is polynomial-time solvable in all classes of graphs of minimum degree 2 for which the decision version of Minimum Node Cover problem is polynomial-time solvable.

For hereditary classes of graphs H, we can easily use a pre-processing allowing to reduce the Minimum Node Cover problem in this class to the same problem in the subclass of graphs in H without leaf nodes. This leads to the following corollary.

Corollary 2

1. RpCP 1 is NP-complete in all hereditary classes of graphs for which the decision version of Minimum Node Cover problem is NP-complete. The problem RpCP is NP-hard in these classes of graphs. 2. RpCP 1 is polynomial-time solvable in all hereditary classes of graphs for which the decision version of Minimum Node Cover problem is polynomial-time solvable.

The proof is not essential for this paper and consequently, to make the text easier to read, it is detailed in Appendix A.1.

In particular, deciding whether the minimum radius is 1 (RpCP 1 ) is NP-complete on planar graphs of maximum degree 3 since the decision version of Minimum Node Cover is NP-complete on this hereditary class [?]. The case of planar graphs of low degree is particularly relevant for real case applications. Since the graph G represents the adjacency graph of zones in the territory, it is planar and in most cases, each zone has a small number of adjacent zones. In some cases however, the underlying graph has even a simpler structure. A common case is a rectangular grid or a subgraph of a grid (called subgrid) when the territory has some "holes" corresponding to large spaces not likely to burn like a lake, a clearing or even a stone field. Based on these cases, a natural question is the complexity of RpCP in bipartite planar graphs and even in subgrids. In what follows, we answer this question.

RpCP 2 in subgrids

The main result of this subsection is Theorem 1 that states the hardness of RpCP 2 in subgrids. To make the presentation easier to read and the proof more intuitive, we first propose a weaker result (Proposition 2) stating the hardness of RpCP 2 in bipartite planar graphs of maximum degree 3. Then, we show how to extend the reduction with a more complex construction to prove Theorem 1. This step uses ideas and techniques proposed in [?,?] to prove NP-hardness results in subgrids for a large range of problems known to be hard in planar bipartite graphs with nodes of degree 2 or 3.

Proposition 2 RpCP 2 is NP-complete in planar bipartite graph with nodes of degree 2 or 3, even if all edges have length 1. RpCP is NP-hard in this case.

Proof We revisit a reduction mentioned in [?] for another problem. The reduction is from the decision version of Minimum Node Cover problem in planar graphs with nodes of degree 2 or 3, known to be NP-complete [?]. Since this Proposition is weaker than Theorem 1, we propose in the main text a sketch of proof and report in Appendix A.2 the proof of the main argument (Claim 1). It is not essential for the paper but helps understanding the proof of Theorem 1.

Given a planar graph G = (V, E) with degrees 2 or 3, one builds a bipartite graph G = (V , E ) by replacing each edge (u, v) with a gadget L uv as presented in Figure 4. G has all its nodes of degree 2 or 3. We then have the following claim:

Claim 1 For any t ≤ |V |, G has a node cover of size t if and only if G has a set C of (t + |E|) nodes with r(C) ≤ 2. Moreover, for each edge (u, v) ∈ E, C has exactly one node in {y u uv , y v uv } and none in {s uv ,t uv }.

Proof (of Claim 1) We prove only the necessary condition and prove the sufficient condition in Appendix A.2. Suppose first that G has a node cover U ⊂ V of size t, we will add to it a set U E of nodes to make it the required set of centers. To make the following construction non-ambiguous, we consider an orientation of the graph G. Consider any edge (u, v) oriented from u to v. If u ∈ U, then we add y v uv to U E and if v ∈ U but u / ∈ U, then we add y u uv to U E . Then, U ∪U E is a set of (t + |E|) nodes of radius 2 (nodes in V are of degree at least 2).

Claim 1 states that the decision version of Minimum Node Cover in planar graphs with nodes of degree 2 or 3 polynomially reduces to RpCP 2 in this class of graphs. Since the former problem is NP-complete in this class, so does the latter, which concludes the proof.

In what follows we show how we can adapt the proof of Proposition 2 to prove a stronger result. It seemed to us easier to devise directly a reduction from Minimum Node Cover than reducing RpCP 2 in planar graphs with nodes of degree 2 or 3 to the same problem in a more restrictive class. The proof of Proposition 2 is given only to make this reduction clearer and more intuitive.

Theorem 1 RpCP 2 is NP-complete in subgrids with nodes of degree 2 or 3, even if all edges have length 1. RpCP is NP-hard in this case.

Proof We already have noted that the problem RpCP 2 is in NP. For clarity, the reduction is divided in two steps.

Step 1: This step follows general ideas proposed in [?,?] for proving NP-hardness results in subgrids. Given a planar graph G = (V, E) with nodes of degree 2 or 3, we first embed it in a grid of polynomial size using a result of [?]: nodes are mapped to nodes of the grid and edges are mapped to non-crossing paths in the grid. Embedding can be done in polynomial time. For seek of simplicity, for every node u of G, we will denote as well by u the node of the grid it maps to. Using this convention, we consider that the node set of G is contained in the node set of the embedding of G in the grid. Thus, any edge (u, v) of G is replaced by a path of length uv between u and v in the embedding for some positive integer uv . The resulting graph is not a subgrid but only a partial subgraph of the grid. The next step will make it a subgrid with, in addition, the required properties to ensure the validity of the reduction.

Step 2: The main idea is inspired from the reduction seen in Proposition 2, where a 4cycle was inserted on each edge of the original graph (gadget of Figure 4). After embedding the graph in a grid, an edge (u, v) of the original graph is a path of length uv and the idea is to insert on this path a 4-cycle like in the previous reduction. The only technical difficulty is to manage the length of the paths between u, v and this 4-cycle to ensure the reduction will work. To this aim, we use the gadget H uv represented in Figure 5 with nodes x 1 uv , x 2 uv , y u+ uv , y u- uv , s uv ,t uv , y v+ uv , y v- uv and z i uv , i = 1, . . . , 9. It can be seen as a gadget similar to Figure 4 (nodes y u+ uv , y u- uv , s uv ,t uv , y v+ uv , y v- uv ) with two paths of 3 and 9 edges attached to y u+ uv and y v+ uv , respectively. Note that this gadget can replace a section of 12 consecutive horizontal edges and similarly, a sequence of 12 vertical edges can be replaced with the same gadget rotated by π 2 . To this purpose, the next step after step 1 is to subdivide every edge in 12 edges by inserting 11 new nodes. This has few advantages: it produces another embedding of the original graph in a grid that ensures that every path P uv in the grid associated with an edge (u, v) of the original graph has now 12 uv edges and with its 12 first edges from u (or from v) either all horizontal or all vertical. In addition, such an expansion gives enough space to ensure we can insert gadgets H uv while guaranteeing the resulting graph to be a subgrid. The strategy is to insert H uv by replacing the 12 first edges on one side of P uv . Since it is non symmetric, we will use an orientation of the original graph G. Starting from the graph obtained at step 1, the second step is summarized below and will conclude the reduction:

-Subdivide every edge in 12 edges by inserting 11 new nodes. The resulting graph, G is a subgrid obtained from G by replacing any edge (u, v) of E with a path P uv of 12 uv edges.

-Select an orientation of each edge of G.

-For an edge (u, v) oriented from u to v, replace the 12 first edges of P uv , starting from u, with the gadget H uv represented in Figure 5 Note that for any edge (u, v) of G oriented from u to v, there is a path of length 12( uv -1) in G from the node z 9 uv and v. In particular, if uv = 1, then z 9 uv = v. Else, the nodes of this path are denoted by z 10 uv , . . . , v = z 9+12( uv -1) uv

. We then denote by H uv the subgraph obtained by adding to H uv the path from z 9 uv to v. If uv = 1, then H uv = H uv . The graph H uv has 6 + 12 uv nodes including u and v. The graph G is obtained from G by replacing each edge (u, v) oriented from u to v with H uv . It is an induced subgrid since for every node x of H uv \ {u, v} and for every node y of H u v \ {u , v } with (u, v) = (u , v ), (x, y) is not an edge in the grid. In addition, nodes in V have the same degree in G than in G and all other nodes in G have degree 2 or 3. So, G has |V | + 4|E| + 12L nodes, all of degree 2 or 3, where

L = ∑ (u,v)∈E ( uv ).
This concludes the construction that can be performed in polynomial time. To conclude the proof of Theorem 1, we need to show that the Minimum Node Cover in the graph G reduces to RpCP 2 in the graph G. For this, we establish two claims. Claim 2 is a technical result proved in Appendix A.3. It is used to prove Claim 3 that immediately concludes the proof.

Claim 2 Let C be a set of nodes of radius at most 2, then for every edge (u, v) of G oriented from u to v, the following holds:

1. C includes at least one node from {s uv ,t uv , y u- uv , y v- uv } 2. C includes at least two nodes from {s uv ,t uv , y u- uv , y v- uv , y u+ uv , y v+ uv }

The proof is given in Appendix A.3.

We are now ready to establish Claim 3. The proof requires to define, for an edge (u, v) ∈ E of G oriented from u to v, two disjoint sets of nodes of H uv :

Denote C v+ uv = {z 3i uv , i = 1, . . . , 2 + 4( uv -1)} (last node is z 6+12( uv -1) uv ); Denote C v- uv = {z 3i-1 uv , i = 1, . . . , 3 + 4( uv -1)} (last node is z 8+12( uv -1) uv
, which is linked to v).

C v+

uv and C v- uv include, after a first node z i 0 uv , every third node along the path z i uv , i 0 ≤ i ≤ 8 + 12( uv -1). Roughly speaking, C v+ uv (resp. C v- uv ) corresponds to the optimal position of centers along the path from y v- uv to v in H uv after deciding to implement a center on node y v+ uv (resp. y v- uv ). Note that on a path with 3k + 1 nodes, at least k -1 centers are required in addition to the two extremities to ensure the radius to be at most 2. The only solution using this number of centers is to place centers at the extremities of the path and on every third node in between. We have

|C v+ uv | = 2 + 4( uv -1) and |C v- uv | = |C v+ uv | + 1.
In Figure 5, nodes in C v+ uv and C v- uv are filled with lines and dots respectively.

Claim 3 G = (V, E) has a node cover K of size k = |K| if and only if G has a set C K of (k + 4|L| + |E|) nodes with r(C K ) = 2.
Proof (of Claim 3) Let us first consider a node cover K of cardinality k in G = (V, E). We complete K in C K in G by adding, for every edge (u, v) ∈ E of G oriented from u to v, 4 uv + 1 nodes as follows:

If u ∈ K, then we add to C K nodes {y u+ uv , y v- uv } ∪C v- uv . If u /
∈ K, then we add to C K nodes {x 1 uv , y u- uv , y v+ uv } ∪C v+ uv .

Since we add 1 + 4 uv nodes for each edge, we have:

|C K | = |K| + 4|L| + |E|.
We can check that the radius of C K is at most 2 and it is at least 2 since, in all cases, we have two consecutive nodes in H uv that are not in C K .

To prove the converse, note first that we need at least 1 + 4 uv nodes of H uv \ {u, v} in C K to ensure r(C K ) ≤ 2. Actually, if there are three consecutive nodes of degree 2 that are not in C K , then the radius is at least 3. Consequently, even if u and v are in C K , we need at least 4 uv centers of H uv \ {u, v, s uv ,t uv , y u- uv , y v- uv } to ensure that every three consecutive nodes of degree 2 include at least one center. The best way to do it is to add {y u+ uv , y v+ uv } ∪ C v+ uv , which makes 4 uv centers. Using Claim 2 (first item), at least one additional node from {s uv ,t uv , y u- uv , y v- uv } should be added in any set of nodes C satisfying r(C) = 2. Let us now assume that G includes a set of nodes C of radius 2. As we just noted, for every edge (u, v) of G, C includes at least 1 + 4 uv nodes of H uv \ {u, v} and thus, C includes in all k + |E| + 4L nodes for some non-negative k. Suppose now that, for an edge Proposition 3 There are graph classes for which RpCP 1 is NP-complete while all graphs in the class have a solution of radius 2 using only two centers (thus, RpCP 2 is trivial).

(u, v) of G, oriented from u to v, neither u nor v is in C. Since r(C) ≤ 2, we have {x 1 uv , x 2 uv } ∩C = / 0 and {z 8+12( uv -1) uv , z 7+12( uv -1) uv } ∩C = / 0.
Proof Consider any class C of graphs with all nodes of degree at least 2 and including at least three independent nodes (i.e., not linked by an edge) such that the decision version of Minimum Node Cover problem is NP-complete on C . The condition that at least three independent nodes exist in any graph G = (V, E) of this class is not restrictive since Node Cover can be trivially solved in polynomial time on graphs that do not satisfy this condition.

The condition ensures that the size of a minimum node cover is at most |V | -3 for any graph G = (V, E) ∈ C . We build the class C of all graphs obtained from a graph G ∈ C by adding two nodes u 0 , v 0 completely linked to all nodes of G. For any graph in C , a minimum node cover includes u 0 , v 0 and a minimum node cover of the graph obtained by removing u 0 , v 0 . If u 0 or v 0 is not included in the node cover, then all other nodes should be included which, by hypothesis, is larger than the proposed solution. So, the decision version of Minimum Node Cover problem and RpCP 1 using Corollary 1 are both NP-complete on C . Note however that for any graph in C obtained from G = (V, E) ∈ C by adding u 0 , v 0 , the set {u 0 , v 0 } is a solution of R2CP of radius 2, which completes the proof.

In the following section, we propose an integer linear programming (ILP) formulation for RpCP. This model is inspired by the models proposed in [?,?] for the p-Center problem. We generalize this model to take into account different fire scenarios and our evacu-ation strategy. When handling both robust p-Center problem under pressure and the usual p-Center problem, the latter will be qualified as deterministic, unless no confusion occurs.

4 Integer Linear Programming model for the robust p-Center problem under Pressure

We propose an ILP model with 0-1 variables representing the maximal radius over all the scenarios to be minimized under linear constraints. The starting point of our model is the deterministic p-Center model presented in [?]. We recall this model in the next sub-section.

Integer Linear Programming model for the deterministic p-Center problem

The deterministic p-Center model presented in [?] is similar to the model proposed by Elloumi, Labbé and Pochet in [?]. In [?], the formulation is based on the observation that the optimal value of the p-Center problem corresponds to one of the distances between two nodes in G. Denote D the finite list of distinct distances between nodes, using the distance of shortest path lengths. Starting with the matrix SP = (d i j ) of the shortest path lengths between every couple of nodes (with d ii = 0 and d i j = d ji = +∞ if there is no path between i and j), D is obtained by sorting in increasing order the T different finite values of the matrix SP:

D min = D 1 < D 2 < D 3 < . . . < D T = D max .
In [?] and [?], two kinds of binary variables are introduced. More specifically, in [?], the following variables are used:

for all j = 1, . . . , n, y j is a binary variable with y j = 1 if a shelter is located on j and 0 otherwise, for all t = 1, . . . , T , u t is a binary variable with u t = 1 if the value of the solution is equal to D t and 0 otherwise.

In [?], Calik and Tansel introduce the following formulation P det :

P det                                      min T ∑ t=1 D t u t (1 det ) s.t. n ∑ j=1 y j = p (2 det ) T ∑ t=1 u t = 1 (3 det ) ∑ j:d i j ≤D t y j ≥ t ∑ q=1 u q i = 1, . . . n, t = 1, . . . , T, (4 det ) y j ∈ {0, 1} j = 1, . . . , n u t ∈ {0, 1} t = 1, . . . , T
Constraint (2 det ) fixes the number of shelters to be located. Constraint (3 det ) ensures that exactly one variable u t is equal to 1 and the corresponding D t value is selected as the objective value according to the objective function (1 det ). If u t = 1, then ∑ t q=1 u q = 1 and Constraints (4 det ) ensure for each node i that at least one shelter is located at a distance less or equal than D t .

The number of binary variables is equal to n +T and the number of constraints is nT +2. The size of this model can be huge since it depends on the number T of distinct shortest path lengths. As explained in [?,?], it is possible to reduce this size: knowing a lower bound LB and an upper bound UB for the optimal value of P det , we can delete some variables since:

u t = 0 , ∀t : D t < LB u t = 0 , ∀t : D t > UB
In the following, for an integer linear program P, its relaxed version, where constraints x ∈ {0, 1} are relaxed as 0 ≤ x ≤ 1, is denoted by LP. The optimal value of the programs P and LP are denoted by v(P) and v(LP), respectively.

In order to determine a lower bound for v(P det ), Daskin [?] proposes an algorithm based on the set covering problem with the following formulation denoted by SC r :

SC r            min n ∑ j=1 y j s.t.
= D s 1 < D s 2 < D s 3 < . . . < D s T s = D s max .
For each s, D s is computed in two stages:

step 1: compute SP s = (d s i j ) the matrix of shortest path lengths from i to j in G s ; extract from SP s the different finite shortest path lengths for all i = s and initialize D s with these values; step 2: ∀k = s and ∀v ∈ N + G s (s) compute the distance l sv + d s vk from s to v to k and add it to D s . Finally, we merge all the lists D s in one ordered set

D rob = {D rob 1 , . . . , D rob T rob }, with D rob min = D rob 1 < D rob 2 < D rob 3 < . . . < D rob T rob = D rob max .
In our formulation P rob , the decision variables are the y variables (similar to P det ) and the u variables with the following interpretation: u t = 1 if and only if, t is the minimum index in {1, . . . , T rob } such that, for any given scenario, all the nodes are at a distance to a shelter less than or equal to D rob t . We then introduce the following formulation for RpCP:

P rob                                                      min T rob ∑ t=1 D rob t u t (1) s.c. n ∑ j=1 y j = p (2) T rob ∑ t=1 u t = 1 (3) ∑ j:d s i j ≤D rob t y j ≥ t ∑ q=1 u q s = 1, . . . , n, i = 1, . . . , s -1, s + 1, . . . , n, t = 1, . . . , T rob (4) ∑ j:l sv +d s v j ≤D rob t y j ≥ t ∑ q=1 u q -y s s = 1, . . . , n, ∀v ∈ N + G s (s),t = 1, . . . , T rob (5) y j ∈ {0, 1} j = 1, . . . , n u t ∈ {0, 1} t = 1, . . . , T rob
Constraints (2), and (3) are similar to constraints (2 det ) and (3 det ) with the only difference that T is replaced with T rob . Constraints (4) ensure that for each node i = s, at least one shelter is located at a distance less than or equal to D rob t in every scenario s. Constraints (5) are specific to RpCP and allow to model the chosen evacuation strategy:

if y s = 1, then a shelter is located on s and constraints (5) are relaxed; if y s = 0, then no shelter is located on s and the set of constraints (5) on all neighbours of s ensure that the evacuation distance (worst case value) in scenario s is less than or equal to D rob t

The number of binary variables is equal to n + T rob and the number of constraints is n 2 T rob + 2mT rob with m the number of edges of G. The size of the model P rob depends on the size of the list D rob leading to huge integer linear programs. Consequently, in order to obtain optimal solutions we have to reduce its size (fixing some variables) and to define specific exact algorithm based on a generalization of a binary search. These methods are presented in the following section and experimental results are given.

Computational study

The size of model P rob depends on the size of the list D rob . Similarly to P det , the size of P rob can be reduced knowing a lower bound LB and an upper bound UB for v(P rob ) since some variables can be fixed as follows:

u t = 0 , ∀t : D rob t < LB u t = 0 , ∀t : D rob t > UB
A first challenge is to determine tight upper and lower bounds. We propose in the following section several methods to compute such bounds.

Upper and lower bounds

We propose four different methods to compute such bounds:

-The first method uses an optimal solution of P det with the algorithm proposed in [?].

Obviously, the value of an optimal solution of RpCP can not be less than the optimal value of P det . We denote LB 1 = v(P det ). When this solution is feasible for RpCP, its value gives an upper bound UB 1 for RpCP. -The second method is an extension of the method presented in section 4.1 to compute a lower bound of P det . The formulation SC r is adapted to RpCP with SC rob r given by:

SC rob r                    min n ∑ j=1 y j s.c. ∑ j:d s i j ≤r y j ≥ 1 ∀s, i = 1, . . . , n (6)
∑ j:l sv +d s v j ≤r y j ≥ 1y s ∀s, ∀v ∈ N + G s (s) ( 7)

y j ∈ {0, 1} j = 1, . . . , n
A binary search can be performed on D rob to find the minimum radius r * for which v(LSC rob r * ) ≤ p. A lower bound for P rob is then r * , denoted by LB 2 . -In a third method, we randomly construct solutions and compute their value for RpCP.

The lowest obtained value represents a second upper bound UB 2 . -The fourth method consists in considering P rob without the Constraints ( 4). The obtained model denoted by RP rob corresponds to the problem where only the evacuation distance of the node s is taken into account for scenario s. It reduces the number of constraints by n 2 T rob . The value of the obtained solution is a lower bound LB 3 for RpCP. Similarly to the first method, if RP rob has an optimal solution, which is feasible for RpCP, it gives a third upper bound UB 3 for RpCP.

In our preliminary experiments, despite using these bounds for P rob , the number of constraints and variables were still too high in order to solve exactly the problem, more precisely even to write the LP instance. For example, for an instance with 100 nodes, 200 edges and T rob = 400, the number of constraints exceeds 4 millions. It took us more than 9 hours and 120 gigabytes of memory usage to obtain the optimal solution.

Thus we propose a general scheme using a generalization of binary search algorithm. As we will see in the experimental results section, the same instance using our algorithms can be solved in less than 40 seconds.

Exact solution method

Consider P(D) a linear programming formulation whose objective value (to be minimized) takes a value from an ordered set D = {D 1 , D 2 , . . . , D T }. Denote LB and UB two initial lower and upper bounds for v(P(D)). A σ -quantile search, presented in Algorithm 2, can be used to solve P(D) by solving at most log σ (T + 1) instances of P(D ) with D ⊆ D. While UB = LB, we repeat the following steps:

-First we compute a restricted set D = {D k 1 , . . . , D k σ } ⊆ D using function Fkernel (Algorithm 1): we delete all values in D that are less than LB or greater than UB, then

Algorithm 1 Fkernel

Require:

D = {D 1 , D 2 , . . . , D T }, LB,UB, σ ∈ N, σ ≥ 3 Ensure: Returns a subset of D 1: Find k 1 ∈ {1, . . . , T } such that D k 1 = LB 2: Find k σ ∈ {1, . . . , T } such that D kσ = UB 3: step ← (k 1 + k σ )/(σ -1) 4: for i ← 2 to σ -1 do 5: k i ← k i-1 + step 6: end for 7: Return {D k 1 , . . . , D k i , . . . , D kσ }
Algorithm 2 σ -quantile search Require: P(D), LB,UB, σ ∈ N, σ ≥ 3 Ensure: Returns the optimal value and an optimal solution to P(D)

1: while UB = LB do 2:

D = {D k 1 , . . . , D kσ } ← Fkernel(D, LB,UB, σ ) 3:
Solve P(D ) 4:

Set q ∈ N + such that D kq ← Optimal value of P(D ) 5:

sol ← Optimal solution of P(D ) 6:

if q = 1 then 7:

UB ← LB 8: else 9:

UB ← D kq 10:

LB ← D k q-1 11:

end if 12: end while 13: Return LB and sol D k 1 = LB, D k σ = UB and the intermediate values correspond to a (σ -1)-quantile. So, between every two consecutive values of D there are roughly the same number of values of D.

-Next we solve P(D ). Let v(P(D )) = D k q , then D k q is an upper bound for v(P(D))

as D k q ∈ D. In addition, there is no feasible solution in P(D ) with value D k q-1 and equivalently in P(D), so D k q-1 is a lower bound for v(P(D)). Note that if D k q = D k 1 , then D k 1 is the optimal solution for P(D). -Finally, LB = D k q-1 and UB = D k q .

Note that for σ = 3, the σ -quantile search is actually a binary search.

Our initial exact algorithm to solve P rob is presented in Algorithm 3. In a first step, all lower and upper bounds are computed and in a second step, a σ -quantile search is performed.

In the next section, we evaluate the computational efficiency of the proposed algorithm.

Experimental results

We implement the Exact Algorithm in Python 3.7 for two sets of instances: p-median instances from OR-Library and subgrids randomly generated. We generate the distance matrices SP = (d i j ) and SP s = (d s i j ) for all scenario s ∈ V using networkx library 2.3. We execute our experiments on a server with 254Gb of RAM and 14

Algorithm 3 Exact Algorithm

Require: G = (V, E), x, σ 1 , σ 2 Ensure: Returns the optimal value and an optimal solution to P rob 1: Generate D rob 2: Solve P det and generate LB 1 and UB 1 3: Compute LB 2 4: Compute UB 2 by generating x random solutions 5: LB 3 , sol ← σ -quantile search(RP rob (D rob ), max{LB 1 , LB 2 }, min{UB 1 ,UB 2 }, σ 1 ) 6: UB 3 ← Value of sol for P rob 7: UB ← min{UB 1 ,UB 2 ,UB 3 } and LB ← max{LB 1 , LB 2 , LB 3 } 8: optValue, optSolution ← σ -quantile search(P rob (D rob ), LB,UB, σ 2 ) 9: Return optValue Intel Core (Haswell; no TSX) Processor at 2.3 Ghz. Mathematical programs are solved with CPLEX 12.9 (with MIPEmphasis option set to 0).

Experimental results on OR-Library Instances

The input data used for the computations are the 40 instances of the p-median problem from the OR-Library ([?]) which are used also for solving the p-center problem ([?,?]). n varies between 100 and 900 nodes and p varies between 5 and n/3 . In the following we focus on the instances which could be solved within 5 hours.

Table 1 contains the value of the upper and lower bounds computed for the instances ordered by the number of their nodes and the value of p. We mark with an * the values of the bounds that are equal to optimal value. It stands out that on all the instances considered, LB 3 = UB 3 = v(P rob ). The equality between UB 3 and LB 3 is not mandatory since only evacuation paths of a subset of nodes are considered in LB 3 . In fact, we record some instances in which LB 3 < v(P rob ) (see in Figure 7) and, in these cases, the computation time is more important. The equality between LB 3 and UB 3 means that the evacuation distance of node s induces the radius of an optimal solution in scenario s. It may be due to the fact that OR-Library instances are considerably sparse.

Concerning lower bounds, LB 2 is also a tight lower bound very close and often equal to v(P rob ), while LB 1 is the worst lower bound. Concerning upper bounds, UB 1 is globally better than UB 2 .

The Exact Algorithm computes all pairs of lower and upper bounds. In order to better understand the trade-off between bounds quality and computational time, we compare the processing time of three variants of the Exact Algorithm:

-EA1 is a version of the Exact Algorithm in which only LB 1 and UB 1 are computed, its processing time is TEA1. This variant is mainly based on the resolution of the deterministic p-Center problem. -EA2 is a version of the Exact Algorithm in which only LB 2 and UB 2 are computed, its processing time is TEA2. This variant is adapted from the Daskin's algorithm for the RpCP. -EA3 is a version of the Exact Algorithm in which only LB 3 and UB 3 are computed, its processing time is TEA3. This variant is a new one specific to RpCP and its evacuation strategy. For pmed35 to pmed40, TEA1 and TEA2 exceed five hours: only EA3 can exactly solve all instances in less than 5 hours.

In Figure 6, we represent the processing times for instances from pmed1 to pmed34 for three different values of p: p = n/3 , p = n/10 and p = 10. For each p, three curves represent the processing times TEA1, TEA2 and TEA3, in function of n. We observe that our dedicated algorithm EA3 performs faster than EA1 and EA2. The reason is twofold: in EA1, the poor quality of the lower bound LB 1 (see Table 1) increases the number of iterations in step 8 of the Exact Algorithm. Conversely, in EA2 the quality of the lower bound LB 2 highly decreases the number of iterations in step 8, however computing LB 2 is very time consuming. More precisely, the generation of the input CPLEX instance for SC rob r requires p=n/3 p=n/10 p=10 Fig. 6: Processing time of each variant of the Exact Algorithm on some OR-Library instances a large amount of time (at least 90% of the total computing time for generating LB 2 ). On the other hand, LB 3 and UB 3 can be computed much faster while providing the best quality of bounds. Thus, EA3 is clearly the most efficient of the three algorithms.

Once we have identified EA3 as the best of the three variants, we must verify whether improvements can be made by adjusting parameters σ 1 and σ 2 . Given the quality of LB 3 , increasing σ 2 to values greater than 3 is counterproductive: one iteration with σ -quantile search is enough to prove that the lower bound is an upper bound. In this case, increasing σ 2 would only increase the size of the ILP model constructed in the step 8 of Exact Algorithm.

However, we can potentially improve LB 3 processing time by using other values of σ 1 . Therefore, we compare the performance of algorithm EA3 on the 34 first instances from the OR-Library with different values of σ 1 in a range of values between 3 and 11. For larger instances, computational times are not reported since they exceed 5 hours for some values of σ 1 . The results are given in Table 2, where, for each instance, the processing time of EA3 is standardised with respect to the processing time of EA3 for σ 1 = 3. Then we use the geometric mean to compare the average processing time for the different values of σ 1 . The experiment reveals that, with σ 1 = 4, EA3 is at least 14 percent faster than the other tested values for σ 1 .

Then we perform a qualitative study to better understand the impact of p on the processing time. In Figure 7, the processing time of EA3 is given for three instances with 100 nodes and 200 edges, named pmed1, pmed4 and pmed5. For each instance, the curve describes the evolution of TEA3 in function of p, with p ranging from 2 to 33. For pmed2, the processing time is relatively stable for all values of p. On the contrary, the processing time to solve pmed4 and pmed5 is much more impacted by the variation of p. Precisely, processing times directly depend on the number of iterations to solve P rob with the σ -quantile search (step 8).

When the number of iterations is equal to 1 (which corresponds to the case LB 3 = v(P rob )) the computation time is stable. But when LB 3 < v(P rob ), the number of iterations increases (up to 5) and the computation time significantly increases. These results underline that for a given size of instance, there is no obvious relationship between the value of p and the complexity of solving P rob .

Our experimental results allow to conclude that the original algorithm EA3 is the best one. In particular we are able to precisely tune the parameters defining Algorithm 4 called EA3*. The efficiency of EA3* comes from the quality of the lower bound obtained with a very efficient 4-quantile search algorithm. With such a bound, in most cases, only one iteration of the 3-quantile search is enough to determine an optimal solution of P rob in step 5.

Instance pmed2

Instance pmed4

Instance pmed5 To measure the advantage of using RpCP over a classical deterministic p-Center problem, we compare the value on the objective function of RpCP of an optimal solution of p-Center problem (referred to as deterministic solution) to the optimal solution value of RpCP. Let us recall that UB 1 exactly corresponds to the value, on RpCP, of a deterministic solution. Thus, we compute in Table 3 the ratio UB 1 /v(P rob ) (none of the instances correspond to n=700 and p = n/3 in OR-Library). It appears that the gap can be quite significant, up to 25%. Consequently, the value of a deterministic solution can be far away from the optimal solution value of RpCP and better solutions can be found solving RpCP.

In the following, we report our experimental results on randomly generated subgrid instances.

Experimental results on random subgrid instances

As explained in Section 3.1, subgrids are relevant graphs for real case applications. Thus, we chose to test our algorithms on random subgrid instances. We generated a set of 18 sub- 3: Ratio between UB 1 and the optimal value of RpCP on OR-Library instances grids along the following steps. From an original undirected unit grid of size (l × w), we generate three subgrids by randomly removing a node with probability 0.05, 0.1 and 0.2.

For each subgrid, we remove isolated and leaf nodes being mandatory shelter locations. We apply this process to three grids of size 10 × 10, 20 × 10 and 20 × 20, into nine unit subgrids SG1, . . . , SG9. We then generate nine weighted subgrids wSG1, . . . , wSG9 by randomly assigning length edge values, from 1 to 10, applied to the edges of SG1, . . . , SG9 respectively. The data set is available at [?].

In the following, we present our synthetic analysis of the results. First, we solve P rob using Exact Algorithm with σ 1 = σ 2 = 3 and x = 10. Table 4 gives the values of the upper and lower bounds computed for the instances SG1, . . . , SG9 ordered by the number of their nodes and the value of p. Likewise, the values of the bounds that are equal to optimal values are marked with an * . For bounds corresponding to non-feasible solutions for RpCP the value is 10000. In Table 4, NA denotes that the values cannot be computed within 5 hours. Please note that the best values achieved for the unsolved instances can be found in Table 10 in Appendix A.4. LB 1 is still the worst lower bound while LB 2 performs slightly better than LB 3 and both of them quite often reach optimal values. Concerning the upper bounds, we note that UB 1 and UB 3 are not so good as those of obtained for OR-Library instances.

Table 5 gives the value of the upper and lower bounds computed for the instances wSG1, . . . , wSG9 ordered by the number of their nodes and the value of p.

In the weighted case, LB 2 and LB 3 outperform LB 1 . Both lower bounds are not as good when the value of p is very low. UB 3 is the best upper bound sometimes reaching optimal value. In this case, the solution of RP rob obtained with EA3 is also the optimal solution for P rob .

In Tables 8 and9 in the Appendix A.4, we record the processing time of our algorithms on the unit subgrids and the weighted subgrids, for the fixed values of p (5,10 and 20), and relatives values of p (n/5 and n/3). For processing time exceeding 5 hours, we record NA. For each instance, the underlined value is the best. Overall, for both unit and weighted subgrids, EA3 is the most efficient algorithm.

In Tables 6 and7, we record the ratio UB 1 /v(P rob ) on the unit subgrids and the weighted subgrids. When the deterministic solution is not feasible, the corresponding ratio value is marked ∞. We observe that the gap, up to 167%, is larger than the one recorded for OR-Library instances. In Table 6, when p = n/3, the ratio UB 1 /v(P rob ) equals 1 in all instances. In this case, we also have LB 1 = 1 which means that the optimal solution of the p-Center problem is a dominating set. We established the following result in [?]: if G is a planar triangle-free graph with no leaf nodes and C a dominating set, r(C) ≤ 2. Moreover, when Table 7: Ratio between UB 1 and optimal solution value of RpCP on weighted subgrids sion version of Minimum Node Cover problem is NP-complete. In addition, we prove that RpCP is NP-hard in subgrids which are relevant for forest fire management problems. In order to propose an exact algorithm, we present a 0-1 Linear Program to model the problem and different methods to compute lower and upper bounds. On the basis of first experimental results, we propose a refined exact algorithm, which is the best current exact algorithm to solve RpCP. The results presented in this paper should be considered as preliminary results in view of solving the problem given. Further research directions will include the design of more efficient heuristics and the study of relevant classes of topologies for real case applications. These classes of instances motivate two main challenges: analyzing the problems and their difficulty in such particular cases and how to take into account specific structures of real case instances in the linear programming approaches so as to make the related models much more tractable. Moreover, this linear programming approach can be adapted to solve variants that are interesting from a practical perspective. For example, under specific circumstances, a shelter can be located only in specific places, or a scenario may induce a fire on different nodes. This is left for future research.

Fig. 1 :

 1 Fig. 1: The operational graph associated with scenario 2 with the solution C = {3, 10}

Fig. 2 :

 2 Fig. 2: A path with 5 nodes and shelters on the extremities

Fig. 3 :

 3 Fig. 3: Comparing the optimal solutions of the 2-Center problem and R2CP

Fig. 4 :

 4 Fig. 4: The gadget L uv for an edge (u, v)

Fig. 5 :

 5 Fig.5: The Gadget H uv for the edge (u, v) oriented from u to v (continuous lines are used for H uv 's edges while dashed lines correspond to edges outside H uv )

  ∑ j:d i j ≤r y j ≥ 1 i = 1, . . . , n y j ∈ {0, 1} j = 1, . . . , n Considering the linear relaxation LSC r , we have v(P det ) > r if v(LSC r ) > p. When r increases, v(LSC r ) decreases and the best lower bound is the smallest value of r ensuring v(LSC r ) ≤ p. The algorithm proposed by Daskin performs a dichotomic search on D to determine such lower bound solving several LSC r . Combining Daskin's algorithm with formulation P det gives the best known results as outlined by Elloumi et al. [?] and Calik et al. [?]. In the following section, we present an extension of model P det to RpCP. 4.2 New model for the robust p-Center problem under pressure We extend the formulation of P det to RpCP. First, we have to replace D by D rob , the list of distinct finite distance values in all G s considering the evacuation strategy. The list D rob is obtained by merging and ordering all the ordered sets D s of distinct finite distances between nodes in G s , for s ∈ V . The elements of D s are denoted by D s min

Fig. 7 :

 7 Fig. 7: Processing time of EA3 on three different instances for p ranging from 2 to 32

  Then, since we cannot have three consecutive nodes of degree 2 outside C, |C ∩ {z i uv , i = 1, . . . , 8 + 12( uv -1)}| ≥ 4 uv -1. Using Claim 2 (second item), at least two nodes from {s uv ,t uv , y u- uv , y v- uv , y u+ uv , y v+ uv } should be in C. In all, C has at least 4 uv + 2 nodes in H uv . These nodes can be replaced with {u, y u+ uv , y v- uv } ∪ C v- uv without augmenting the cardinality of the set. By repeating this transformation, we obtain a set of nodes C such that |C | ≤ k + 4|L| + |E|, where C ∩V is a node cover of size at most k. This completes the proof of Claim 3. Claim 3 states that the decision version of Minimum Node Cover in planar graphs with nodes of degree 2 or 3 polynomially reduces to RpCP 2 in this class of graphs. Since the former problem is NP-complete in this class [?], so does the latter. This concludes the proof of Theorem 1. 3.3 Does increasing the radius make the decision problem harder? Subsection 3.2 gives an example of graph classes for which RpCP 2 is hard while RpCP 1 is polynomially solvable. Meanwhile, Subsection 3.1 outlines classes of graphs for which RpCP 1 is hard without any evidence about hardness of RpCP 2 on these classes. These results make natural the question of the hardness of RpCP k when k varies. Is there a reduction allowing to state the hardness of RpCP k+1 on a class if RpCP k is known to be hard on this class? In particular, can we conclude results for larger k values for the classes studied here? The following remark gives evidence of graph classes on which RpCP 1 is hard but RpCP 2 is trivial. This justifies that such a reduction from a k to k + 1 does not exist in the general case and consequently, the hardness of RpCP k on a given graph class requires to be studied for any value of k and cannot be deduced, in general, from hardness results dealing with different values of k. This leaves open avenues for future researches with, in particular, the challenge to devise reductions working for different values of k.

Table 1 :

 1 Optimal solution values and bound values for OR-Library instances

	Instance	n	|E|	p	OPT	LB1	UB1	LB2	UB2	LB3	UB3
	pmed1	100	200	5	222	127	251	221	263	222*	222*
	pmed2	100	200	10	194	98	229	192	251	194*	194*
	pmed3	100	200	10	191	93	226	186	240	191*	191*
	pmed4	100	200	20	157	74	184	156	218	157*	157*
	pmed5	100	200	33	115	48	144	115*	180	115*	115*
	pmed6	200	800	5	180	84	208	180*	205	180*	180*
	pmed7	200	800	10	156	64	163	155	180	156*	156*
	pmed8	200	800	20	143	55	153	143*	188	143*	143*
	pmed9	200	800	40	124	37	136	124*	164	124*	124*
	pmed10	200	800	67	100	20	118	100*	136	100*	100*
	pmed11	300	1800	5	153	59	157	153*	169	153*	153*
	pmed12	300	1800	10	145	51	150	145*	171	145*	145*
	pmed13	300	1800	30	129	36	136	128	154	129*	129*
	pmed14	300	1800	60	116	26	125	115	145	116*	116*
	pmed15	300	1800	100	105	18	118	105*	133	105*	105*
	pmed16	400	3200	5	143	47	147	143*	153	143*	143*
	pmed17	400	3200	10	136	39	139	136*	150	136*	136*
	pmed18	400	3200	40	122	28	127	122*	144	122*	122*
	pmed19	400	3200	80	112	18	119	111	132	112*	112*
	pmed20	400	3200	133	103	13	113	103*	130	103*	103*
	pmed21	500	5000	5	137	40	139	137*	147	137*	137*
	pmed22	500	5000	10	133	38	137	133*	147	133*	133*
	pmed23	500	5000	50	118	22	122	118*	135	118*	118*
	pmed24	500	5000	100	110	15	115	110*	132	110*	110*
	pmed25	500	5000	167	103	11	113	103*	131	103*	103*
	pmed26	600	7200	5	134	38	137	134*	147	134*	134*
	pmed27	600	7200	10	128	32	132	128*	139	128*	128*
	pmed28	600	7200	60	114	18	118	114*	134	114*	114*
	pmed29	600	7200	120	108	13	113	108*	130	108*	108*
	pmed30	600	7200	200	103	9	109	103*	124	103*	103*
	pmed31	700	9800	5	128	30	136	124	135	128*	128*
	pmed32	700	9800	10	127	29	128	123	164	127*	127*
	pmed33	700	9800	70	113	15	119	110	126	113*	113*
	pmed34	700	9800	140	107	11	111	107*	125	107*	107*
	pmed35	800	12800	5	128	30	130	-	136	128*	128*
	pmed36	800	12800	10	125	27	127	-	134	125*	125*
	pmed37	800	12800	80	112	15	121	-	130	112*	112*
	pmed38	900	16200	5	127	29	129	-	138	127*	127*
	pmed39	900	16200	10	122	23	127	-	174	122*	122*
	pmed40	900	16200	90	111	13	113	-	123	111*	111*

Table 2 :

 2 Comparison of the standardized execution time of EA3 for different values of σ 1

Table

			p
	n	5	10	n/10	n/3
	100	1.13	1.18	1.18	1.25
	200	1.16	1.04	1.07	1.18
	300	1.03	1.03	1.05	1.12
	400	1.03	1.02	1.04	1.10
	500	1.01	1.03	1.03	1.10
	600	1.02	1.03	1.04	1.06
	700	1.06	1.01	1.05	-

Table 4 :

 4 Optimal solution values and bound values for unit subgrids p = n/3, there is no solution with radius equals to 1. So, the deterministic solution is an optimal solution for RpCP.

	Instance
	Graph	p	n	|E|	OPT	LB 1	LB 2	LB 3	UB 1	UB 2	UB 3
	SG1	5	73	102	7	4	7*	6	10000	9	10000
	SG1	10	73	102	4	3	4*	4*	10000	7	6
	SG1	20	73	102	2	1	2*	2*	2*	4	2*
	SG1	14	73	102	3	2	3*	3*	8	6	7
	SG1	24	73	102	2	1	2*	2*	2*	4	2*
	SG2	5	92	153	7	4	6	5	8	8	9
	SG2	10	92	153	4	3	4*	4*	6	6	6
	SG2	20	92	153	3	2	3*	3*	4	5	4
	SG2	18	92	153	3	2	3*	3*	4	5	5
	SG2	30	92	153	2	1	2*	2*	2*	4	2*
	SG3	5	96	167	7	4	5	5	8	8	9
	SG3	10	96	167	4	3	4*	4*	6	5	6
	SG3	20	96	167	3	2	3*	3*	4	5	4
	SG3	19	96	167	3	2	3*	3*	4	4	4
	SG3	32	96	167	2	1	2*	2*	2*	3	2*
	SG4	5	143	218	11	6	10	9	20	13	13
	SG4	10	143	218	7	4	6	6	10	10	10
	SG4	20	143	218	4	2	4*	4*	6	7	7
	SG4	28	143	218	3	2	3*	3*	5	6	5
	SG4	47	143	218	2	1	2*	2*	2*	4	2*
	SG5	5	179	304	9	6	8	8	14	11	12
	SG5	10	179	304	6	4	6*	5	8	8	9
	SG5	20	179	304	4	2	4*	4*	7	7	6
	SG5	35	179	304	3	2	3*	3*	4	5	4
	SG5	59	179	304	2	1	2*	2*	2*	4	2*
	SG6	5	190	332	9	6	8	8	11	12	11
	SG6	10	190	332	6	4	5	5	8	9	7
	SG6	20	190	332	4	3	4*	4*	6	7	6
	SG6	38	190	332	3	2	3*	3*	5	5	6
	SG6	63	190	332	2	1	2*	2*	2*	4	2*
	SG7	5	306	483	15	9	13	12	10000	10000	10000
	SG7	10	306	483	9	6	9*	8	10000	17	10000
	SG7	20	306	483	6	4	6*	6*	9	11	8
	SG7	61	306	483	3	2	3*	3*	10000	7	10000
	SG7	102	306	483	2	1	2*	2*	2*	5	2*
	SG8	5	359	620	NA	8	11	10	16	16	16
	SG8	10	359	620	NA	6	7	7	11	12	12
	SG8	20	359	620	NA	4	5	5	11	9	9
	SG8	71	359	620	3	2	3*	3*	5	5	5
	SG8	119	359	620	2	1	2*	2*	2*	4	2*
	SG9	5	383	704	NA	9	11	11	16	16	18
	SG9	10	383	704	NA	5	7	7	11	12	13
	SG9	20	383	704	NA	4	5	3	8	8	9
	SG9	76	383	704	3	2	3*	3*	4	5	4
	SG9	127	383	704	2	1	2*	2*	2*	4	2*

Table 5 :

 5 Optimal solution values and bound values for weighted subgrids 6 ConclusionWe introduce a new version of the p-Center problem motivated by the context of evacuation in case of wildfires. We call it the Robust p-Center problem under Pressure and emphasize its differences with the existing models in the literature. We study the complexity of this problem: we show that RpCP is NP-hard in all hereditary classes of graphs where the deci-

			p
	Instance	5	10	n/5	n/3
	SG1	∞	∞	2.67	1
	SG2	1.14	1.5	1.33	1
	SG3	1.14	1.5	1.33	1
	SG4	1.82	1.43	1.67	1
	SG5	1.56	1.33	1.33	1
	SG6	1.22	1.33	1.67	1
	SG7	∞	∞	∞	1
	SG8	NA	NA	∞	1
	SG9	NA	NA	1.33	1

Table 6 :

 6 Ratio between UB 1 and optimal solution value of RpCP on unit subgrids

			p
	Instance	5	10	n/5	n/3
	wSG1	∞	1.16	1.45	1.54
	wSG2	1.21	1.25	1.19	1.25
	wSG3	1.3	1.32	1.4	1.36
	wSG4	1.3	1.26	1.47	1.92
	wSG5	1.31	1.14	1.41	1.85
	wSG6	1.2	1.47	1.24	1.46
	wSG7	∞	∞	∞	∞
	wSG8	1.3	1.34	1.29	1.31
	wSG9	1.36	1.34	1.29	1.31

Table 8 :

 8 Processing time of each variant of the Exact Algorithm on the unit subgrids

	Instance
	Graph	p	n	|E|	T EA1	T EA2	T EA3
	SG1	5	73	102	8.62	6.42	4.62
	SG1	10	73	102	4.48	6.14	2.6
	SG1	14	73	102	3.37	6.25	2.43
	SG1	20	73	102	0.84	5.42	0.31
	SG1	24	73	102	0.81	5.58	0.32
	SG2	5	92	153	15.02	19.87	13.47
	SG2	10	92	153	5.98	11.83	5.46
	SG2	18	92	153	2.71	10.19	3.8
	SG2	20	92	153	2.86	10.01	2.3
	SG2	30	92	153	1.43	10.23	0.59
	SG3	5	96	167	29.38	38.74	31
	SG3	10	96	167	7.54	9.74	6.08
	SG3	19	96	167	3.21	9.82	2.66
	SG3	20	96	167	3.28	11.87	2.61
	SG3	32	96	167	1.66	10.56	0.67
	SG4	5	143	218	61.04	72.64	32.84
	SG4	10	143	218	29.54	52.4	27.98
	SG4	20	143	218	14.57	39.34	13.16
	SG4	28	143	218	13.34	36.56	10.61
	SG4	47	143	218	4.75	37.96	1.42
	SG5	5	179	304	188.98	281.31	317.06
	SG5	10	179	304	70.04	80	57.87
	SG5	20	179	304	37.37	72.1	23.35
	SG5	35	179	304	15.16	69.52	12.48
	SG5	59	179	304	8.8	73.25	2.61
	SG6	5	190	332	399.16	855.92	486.3
	SG6	10	190	332	143.33	236.57	101.38
	SG6	20	190	332	44.3	88.26	30.09
	SG6	38	190	332	31.71	83.41	25.69
	SG6	63	190	332	10.74	87.64	3.13
	SG7	5	306	483	1298.46	1364.17	1683.74
	SG7	10	306	483	585.89	394.51	924.87
	SG7	20	306	483	205.19	354.8	110.61
	SG7	61	306	483	309.12	304.88	136.95
	SG7	102	306	483	37.03	324.04	8.53
	SG8	5	359	620	NA	NA	NA
	SG8	10	359	620	NA	NA	NA
	SG8	20	359	620	NA	NA	NA
	SG8	71	359	620	355.2	675	258.2
	SG8	119	359	620	88.4	791.4	40.3
	SG9	5	383	704	NA	NA	NA
	SG9	10	383	704	NA	NA	NA
	SG9	20	383	704	NA	NA	NA
	SG9	76	383	704	196.3	776.4	124.4
	SG9	127	383	704	95	837	21.9

Table 9 :

 9 Processing time of each variant of the Exact Algorithm on the weighted subgrids

	Instance
	Graph	p	n	|E|	Best Value = UB	LB
	SG8	5	359	620	14	11
	SG8	10	359	620	9	7
	SG8	20	359	620	7	5
	SG9	5	383	704	13	11
	SG9	10	383	704	11	7
	SG9	20	383	704	7	5

Table 10 :

 10 Best results obtained on the unsolved instances of unit subgrids

We would like to thank particularly the Pau Costa Foundation (www.paucostafoundation.org/), the Fire Organisation of Andalucia (INFOCA, www.juntadeandalucia.es/medioambiente/site/ portalweb/), the Fire Organisation of Corsica (SDIS2B, www.sdis2b.fr) and the Country Fire Authority of Victoria, Australia (CFA, www.cfa.vic.gov.au).

Acknowledgements

We would like to thank two anonymous referees, whose comments have significantly improved the paper.

Work supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 691161

Appendix

A.1 Proof of Corollary 2

Corollary 2 1 RpCP 1 is NP-complete in all hereditary classes of graphs for which the decision version of Minimum Node Cover problem is NP-complete. The problem RpCP is NP-hard in these classes of graphs. 2 RpCP 1 is polynomial-time solvable in all hereditary classes of graphs for which the decision version of Minimum Node Cover problem is polynomial-time solvable.

Proof (1): Suppose H is an hereditary class of graphs for which the decision version of Node Cover is NPcomplete. We show it is reducible to RpCP 1 in this class. Consider (G, k), with G ∈ H and k ∈ N an instance of the decision version of Node Cover. We transform G using the following pre-processing: while the graph has a leaf node x and k ≥ 0, add to the cover the unique neighbour y of x, remove x, y and all edges incident to y and subtract 1 to k. If k becomes negative, then the original instance was a no-instance. If all nodes of G are eliminated and k ≥ 0, then it was a yes-instance. In other cases, if (G , k ) are the resulting graph and constant, the original instance is a yes-instance if and only if (G , k ) is a yes-instance. Since G is in H and has no leaf node, this is equivalent to deciding whether there is in G a set of k nodes and radius at most 1.

(2): Suppose H is an hereditary class of graphs for which the decision version of Node Cover can be solved in polynomial-time and consider an instance (G, k), with G ∈ H and k ∈ N of RpCP 1 . Using Proposition 1, it is a yes-instance if and only if there is a node cover including all leaf nodes and of cardinality at most k. If G has k leaf nodes, with k > k, then it is a no-instance. If k ≤ k now, this is equivalent to say that the subgraph obtained by removing the leaf nodes and the related incident edges has a node cover of size kk . This question can be answered in polynomial-time since the resulting graph is in H.

A.2 Proof of Claim 1 (sufficient condition)

Claim 1 For any t ≤ |V |, G has a node cover of size t if and only if G has a set C of (t + |E|) nodes with r(C) ≤ 2. Moreover, for each edge (u, v) ∈ E, C has exactly one node in {y u uv , y v uv } and none in {s uv ,t uv }.

Proof We have already shown the necessary condition and prove only the sufficient condition in this Appendix. Suppose that a set U satisfies r(U ) ≤ 2. We first note that, for any edge (u, v) oriented from u to v, U ∩ {y u uv , s uv , y v uv ,t uv } = / 0. Else, in the scenario where y u uv is on fire, any evacuation path from y u uv passing through s uv or t uv would be of length at least 3. We perform such a transformation for every edge (u, v) ∈ E one by one in order to compute a set U of radius at most 2 and of cardinality at most t + |E| satisfying, for each edge (u, v) ∈ E, |U ∩ {y u uv , s uv , y v uv ,t uv }| = 1. If it has less than t + |E| nodes, then we can add nodes in V to make it of cardinality t + |E| without increasing its radius (recall t ≤ |V |). Now, we note that U ∩V is a node cover of G of cardinality t, which will complete the proof. Consider an edge (u, v), oriented from u to v and such that u / ∈ U. In the scenario where u is on fire, y u uv needs to be in U, else an evacuation path in the direction y u uv would be of length greater than 2. By construction, in this case we have U ∩ {s uv , y v uv ,t uv } = / 0 and consequently v needs to be in U since in the opposite case, an evacuation path from v through y v uv would be of length 3. This concludes the proof of Claim 1.

A.3 Proof of Claim 2

Claim 2 Let C be a set of nodes of radius at most 2, then for every edge (u, v) of G oriented from u to v, the following holds:

1. C includes at least one node from {s uv ,t uv , y u- uv , y v- uv } 2. C includes at least two nodes from {s uv ,t uv , y u- uv , y v- uv , y u+ uv , y v+ uv }

Proof (1): If C ∩ {s uv ,t uv , y u- uv , y v- uv } = / 0, then if y v- uv is in fire, its evacuation distance is at least 3 (induced by the evacuation path through s uv and y u- uv ). The same occurs for y u- uv .

(2): If s uv ∈ C and y u- uv / ∈ C, then t uv or y v- uv should be in C to ensure an evacuation distance at most 2 for y u- uv when it is on fire. So, |C ∩ {s uv ,t uv , y u- uv , y v- uv , y u+ uv , y v+ uv }| ≥ 2. Similarly, if t uv ∈ C, we have |C ∩ {s uv ,t uv , y u- uv , y v- uv , y u+ uv , y v+ uv }| ≥ 2. Finally, if s uv ,t uv / ∈ C, then we need C ∩{y u- uv , y u+ uv } = / 0 and C ∩{y v- uv , y v+ uv } = / 0 to ensure an evacuation distance at most 2 when t uv is on fire.

A.4 Tables