Emergency Response in Complex Buildings: Automated Selection of Safest and Balanced Routes

The extreme importance of emergency response in complex buildings during natural and human‐induced disasters has been widely acknowledged. In particular, there is a need for efficient algorithms for finding safest evacuation routes, which would take into account the 3‐D structure of buildings, their relevant semantics, and the nature and shape of hazards. In this article, we propose algorithms for safest routes and balanced routes in buildings, where an extreme event with many epicenters is occurring. In a balanced route, a trade‐off between route length and hazard proximity is made. The algorithms are based on a novel approach that integrates a multiattribute decision‐making technique, Dijkstra's classical algorithm and the introduced hazard proximity numbers, hazard propagation coefficient and proximity index for a route.


INTRODUCTION
Emergency response in the built environment is being widely studied, with a significant surge of interest in this area after 9/11. The focus of such studies is on rescue and evacuation, which are based on route finding and indoor navigation (Choi and Lee, 2009;Kwan and Lee, 2005;Lay, 2007;Lee, 2007;Zlatanova, 2011, 2012;Vanclooster et al., 2014). Despite many publications in this field in recent times, there is still a lack of appropriate evacuation algorithms and their implementations (Lee, 2007;Lee and Zlatanova, 2008;Meijers et al., 2005;Vanclooster et al., 2010). Kwan and Lee (2005) investigated possible improvements of navigable networks for the particular pur-pose of facilitating quick emergency response to terrorist attacks in the integrated system of the ground transportation system and multistory office buildings. They concluded that extending the standard 2-D GIS (Two-Dimensional Geographic Information System) to a real-time 3-D GIS has "considerable potential for improving the speed of emergency response after terrorist attacks on multilevel structures in urban areas." Meijers et al. (2005) developed a semantic model representing interior spaces in a building. Their model can be used for an intelligent computation of evacuation routes. Lee (2007) reviewed 3-D models and building evacuation models, and developed a pedestrian-based indoor navigation model using 3-D GIS. The human behavior was studied by Choi and Lee (2009) using a social force model. Recently, an advanced configurable crowd model for different behaviors and scenarios was developed by Sun and Wu (2014), in particular, the simulation of evacuation in a building was implemented. Liu and Zlatanova (2011) proposed a new door-to-door approach for finding routes between rooms and also a detailed route in a single room. Vanclooster et al. (2010) developed a capacity constrained flow algorithm on a 3-D geometric network model. In 2014, Vanclooster et al. applied Grum's least risk path algorithm to an indoor space for minimizing risks of getting lost, and proposed several improvements to Grum's algorithm to make it more compatible with indoor networks. The role of elevators and stairs in efficient evacuation was investigated by Lay (2007).
A spatial model always underlies any evacuation algorithms, including those mentioned above. A good example of such a model was given by Kwan and Lee (2005, Figure 5). The structure of a building is represented as a logical network, where the nodes represent spatial objects such as rooms, corridors, and other navigable areas. The edges represent navigable connections between adjacent objects. The network can be further extended to a geometric network to model precise geometric properties (e.g., distance between nodes and their locations) and provide real navigation routes. This representation can be used for graph algorithms such as Dijkstra's or A* algorithms for finding shortest routes in evacuation planning (Dijkstra, 1959;Hart et al., 1968). Also, there is a considerable potential for using Ant Colony Optimization algorithms, which have been successfully applied for finding evacuation routes during a tsunami (Forcael et al., 2014), combining both safety and distance to determine the best route. An advanced and efficient strategy for the shortest path problem with uncertain travel cost can be found in Shahabi et al. (2015). Their approach might be particularly relevant for evacuation algorithms in the built environment because the distribution of people in a building during an extreme event is typically unknown, so that the link travel time function is uncertain.
As a rule, 2-D floor plans are used for reconstruction of horizontal navigable networks, and a 3-D building is obtained by linking contiguous floors at some connection points, for example, staircases. Spatial relationships between the rooms in the vertical direction are not reflected in the model. This solution is sufficient for a simple analysis of indoor human movements, but other important phenomena related to emergency response, for example, fire spread or heat propagation, cannot be simulated in such models. Also, many 3-D navigation models represent buildings with a simple, often regular, structure, which is not sufficient for complex interiors. Liu and Zlatanova (2012) proposed an interesting concept of automatic navigable network generation based on the geometry and semantics of a building. Their method requires a valid spatial model with preserved consistency between the geometry and topology, which is not always readily available.
In their pioneering work Kisko and Francis (1985) proposed the EVACNET+ software application for evacuation scenarios planning. By formulating the task at hand as an optimization problem, they developed a user-friendly interface allowing a user interaction, for the computation of egress times and the determination of problematic building locations. For the sake of efficient emergency evacuation, workload relaxation techniques from querying theory were used by Deng et al. (2008) to deal with the curse of complexity in large buildings. For this purpose, the authors developed tools for modelling and optimizing the occupancy evolution during evacuation scenarios. They derived complexity lower bounds on the evacuation time and consolidated their findings by providing realistic building simulations. Pursals and Garzon (2009) presented a new formulation for the problem of building evacuation, through the development of a model for occupants' movement, allowing a good choice of evacuation paths during emergency scenarios. In contrast to the approach for the proper selection of evacuation routes leading trapped occupants to main building exit points, Park et al. (2009) focused on computing optimal routes leading search and rescue personal to disaster locations. For this purpose, the authors developed a time-dependent optimal routing solution based on a network representing the building configuration, which has been enriched by relevant information about the facility. Cellular automata have been widely used for evacuation simulations. For example, a crowd simulation model for large facilities was given by Abdelghany et al. (2010) for replicating the selection of exit gates, based on a trade-off between travel distance and the level of congestion. In the work of Kirchner and Schadschneider (2002) such models have been applied to analyze evacuation scenarios and reduce egress time for situations with a small number of doors. In the same spirit, Daoliang et al. (2006) developed a 2-D cellular automata prototype. It was used in evacuation scenarios for the simulation of exit dynamics of the occupants of a building, that is, for better choices of egress doors. They also provided some hints on the choice of the dimensions of the building exit points; this can be helpful for good building designs. Inspired by this work, Varas et al. (2007) also applied cellular automata but focused instead on studying the effect of fixed obstacles. The problem of the presence of obstacles in evacuation simulations has also been studied by Huang and Guo (2008), who proposed the model of a floor field and used a rectangular lattice site for the computation of navigation routes within it. One of their main conclusions was that familiarizing the occupants with rooms' inner configurations and exit door locations plays an important role in reducing the evacuation time. Based on available video records of student evacuation from a classroom, typical evacuation characteristics were investigated by Zhang et al. (2008) using an improved multigrid model. The simulation process was modelled and its comparison with real evacuation experiments showed its closeness to reality, thus opening the way for understanding evacuation behavioral aspects. The discussed articles are summarized in Table 1.
In this article, we present an algorithm for finding the safest route in a building, where an extreme event with many epicenters is occurring. Another algorithm produces a balanced route, which is achieved by a trade-off between route length and hazard proximity. The algorithms are based on a novel approach that integrates the recently developed 3-D building model described in the next section, a multiattribute decision- Table 1 Summary of the literature 3-D modelling , Boguslawski (2011, Choi (2009), Kwan (2005), Lee (2007Lee ( , 2008, Lienhardt (1991), Liu (2012), Meijers (2005), Weiler (1988) Route/evacuation optimization Dijkstra (1959), Hart (1968), Kisko (1985), Liu (2011), Park (2009 making technique, Dijkstra's classical algorithm and the introduced hazard proximity numbers, hazard propagation coefficient and proximity index for a route. This study proposes the enhanced concept of hazard proximity with two criteria: distance and the number of obstructions (i.e., walls and floors) between the epicenter of an extreme event and points in the building. The algorithms are validated by testing them on different buildings and discussing the results in Section 5. Note that the underlying 3-D model is constructed automatically and it includes all the semantics (3-D or others) necessary for the aforementioned algorithms.

THE BUILDING INFORMATION MODELLING-GEOGRAPHY INFORMATION SCIENCE (BIM-GIS) MODEL
The research in this article is based on the 3-D building model recently developed by Boguslawski et al. (2015) and Barki et al. (2015). This model is an integration of BIM technology and GIS analysis. Ideally, the Industry Foundation Classes (IFC) format could be used to exchange model between the design environment and GIS analysis tool. However, the information content in IFC is very high for GIS applications such as emergency management. The original model needs to be simplified in the process of generalization to reduce the storage cost and extract building elements and geometry essential for reconstruction of indoor spatial relations. This can be achieved through development of generalization procedures or using simplified models, which can be automatically generated in commercial software packages, for example, Green Building XML (gbXML).
We consider the gbXML format as a simplified BIM model, which is the main data format used as an input in our research: a detailed geometry is simplified to a level sufficient to preserve the adjacency relationship between rooms, which is essential. The information about room volumes and wall surfaces, including openings (i.e., doors and windows), is accompanied by attributes such as room names, etc. The openings are used for conventional navigation and egress routes computation. However, other alternative routes can also be considered in case of direct hazard, for example, walls or partitions can be drilled to get access to adjacent rooms. This requires additional information about construction materials, which can be obtained from the original model. The external links in gbXML allow one to look up the required parameters if the original model is available, for example, the IFC model. For example, suppose that there are two adjacent rooms (without a door between them) modelled by nodes u and v, and we know that the wall between them is thin, which is an attribute of the link uv. Then the link uv is made available for navigation in the network and the user is informed if this link happens to be in the egress route.
Our models are Boundary Representations (B-Rep) of the concerned buildings, where volumes (cells) are enclosed by faces, faces by edges, and edges by points. A B-Rep may be modelled and implemented by various data structures, for example, the radial-edge (Weiler, 1988), G-Maps (Lienhardt, 1991), or recently developed dual half-edges (DHE) (Boguslawski, 2011). The DHE structure has been adopted to encode the model geometry as the underlying navigable network can be simultaneously and automatically constructed together with the 3-D model.
A simple model of a building stored in the Autodesk Revit format was exported to the gbXML model and reconstructed using the primal/dual DHE data structure (see Figure 1). Surfaces stored in the gbXML model represent boundaries of spaces, that is, rooms. Information about adjacent spaces and incorporated openings, for example, doors and windows, is attached to each surface. Original surfaces are represented with DHE as double-sided faces. Adjacency information is used to group surfaces into sets representing boundaries of separate spaces. They are merged together (see Figure 2) along adjacent edges into closed cells representing rooms using the Cardboard and Tape method (Boguslawski, 2011). The resulting model is a cell complex, where rooms are represented as cells in the primal structure with an associated dual node unambiguously representing this cell. Adjacent cells are connected by dual edges bounded by dual nodes (see Figure 3).
A door between two rooms is represented as a zerovolume cell with an associated unique node. Thus, there  (c) room with incorporated openings. are two dual edges connecting the first room node to the door node, and the door node to the second room node. The same idea is applied in case of a door between an internal room and an external space. The latter is represented as a cell or a set of connected cells if it was partitioned. The complete graph of indoor connections is shown in Figure 1b. Some openings, which are not directly connected to the boundary of the enclosing surface, for example, windows, are connected to the surface boundary by bridge edges (dotted lines in Figure  1a). A bridge edge is not a part of the original model, but it is introduced to preserve a valid topology of the B-Rep model.
The structure of the model presented in Figure 1 is simple and was reconstructed without additional improvement or validation. However, models with more complex structure exported to the gbXML format must be processed first to reconstruct a valid naviga-ble network. Some common issues are unclosed and overlapping cells. Such models are valid for most of engineering analyses, but not for GIS, which requires a complete topology with a proper representation of spatial relations among objects. For validated models, the graph of connections reflecting spatial relationships among cells in the complex is created automatically using the DHE construction operators. For further details, see the article by Boguslawski et al. (2015).
Summarizing the BIM-GIS model, a building in general consists of several connected rooms that have volumes (corridors, offices, storage spaces, etc. are considered as rooms too), so they are represented by primal cells. The geometry of a room is modelled by the links and nodes of a cell, and relations between adjacent rooms are represented with dual links connecting the corresponding cells. Those relations are described in terms of access level from one room to the adjacent one: access to the next room is by a door; the next room is not accessible because of a wall, but if the wall is thin a hole can be made. It may not be possible to get directly to the next room if the wall is made of concrete. This is an example of a basic set of attributes that are assigned to connections between rooms and then used as weights in graph traversal algorithms, for example, Dijkstra's algorithm. Rooms are not the only objects in a building that are important. Walls, doors, and windows are essential and included in the model. They are represented as cells with geometry, and some attributes are assigned to them. Figure 4a shows a hypothetical building created in Revit. The building is then exported to gbXML and reconstructed with the DHE structure ( Figure 4b). The latter is the primal model, which includes the geometry of the building. The graph of connections, that is, the dual model, is then automatically created (Figure 4c).  The building has three exits and three stairwells, and all floors have a similar structure, except for the exits at the ground floor. The plan of the ground floor is shown in Figure 5.

TEST BUILDING AND TESSELLATION OF CORRIDORS
For navigation purposes, all corridors and large open spaces in a building should be partitioned automatically to generate a navigable network reflecting real navigation routes. Without a proper partition, an incorrect distance between nodes (cells) might be calculated and hence wrong connections (links) might be selected for evacuation. This is illustrated in Figure 5, where a person has to go from Room 11 to the nearest exit. According to the current model with the corridor modelled by just one node, the person would follow a dashed line, while the actual walking pattern should be different; one possible path is shown in Figure 6 (the bold line). To achieve an appropriate tessellation, we use an approach based on the Voronoi Diagrams (VD). VD partition a space into a set of adjacent cells represented as a graph. The dual graph to VD consists of links connecting adjacent dual nodes, which represent primal cells. These links form a network we use for navigation. It should be noted that for some concave shapes a Voronoi tessellation may produce cells, which are split into several unconnected parts, for example, when a boundary of a corridor overlaps with the cell and some parts of the divided cell are not connected to the cell enclosing the dual node. However, in our implementation this situation does not exist. All dual nodes are enclosed by exactly one Voronoi cell.
For the corridors of the building in Figure 4a, this is illustrated in Figure 6, where one possible tessellation of the corridor is shown. Note that some cells in this tessellation have a node, for example, the node u in Figure 6, which is either a door or a concave corner. Hence, such a cell is already associated with one primal node in the 3-D model. However, for other cells, new tessellation nodes are introduced with the corresponding links, an example is the node v in Figure 6. Note that different tessellations are possible and the density of tessellation may be higher for larger areas.
The choice of the tessellation density depends on the precision, which is required in the navigation route. Let us denote by G the original 3-D dual model of the building without tessellation, that is, the graph of connections between cells, and by G + the 3-D model of the building with all the necessary tessellations. Thus, G + has some additional nodes such as corner points and new nodes added during the tessellation process of corridors and large open spaces. Examples of the graphs G and G + are given in Figure 7, where only navigable links are shown. Note that non-navigable links, which are not shown in Figure 7, represent physical obstructions, for example, the link between the nodes in the left and right top corners. Now, the shortest route from a selected room to the nearest exit can be easily calculated using Dijkstra's algorithm.

ALGORITHMS FOR SAFEST AND BALANCED ROUTES
In this section, two algorithms for finding safest and balanced routes in buildings are presented. The algorithms are designed for complex buildings, for example, high-rise buildings, shopping malls, where an extreme event is occurring. Typically, for people it is safer to stay further away from the epicenters of an extreme event in terms of the distance and the number of obstructions (e.g., walls). In this article, we consider extreme events having this natural property. Examples of such events are fire, and terrorist activities such as bomb attacks or hostage taking situations. For the latter, the location of a bomb or a terrorist is an "epicenter," whereas for fire the "epicenters" can be defined as points (nodes) in the building, where the temperature exceeds a certain threshold. Walls, floors, ceilings are called "obstructions." For a particular point p in the building and an extreme event with the epicenter z, the hazard for that point is a function of the direct distance from p to the epicenter z and the minimum number of obstructions between p and z. Note that an extreme event may have many epicenters. For convenience, we summarize the notation, which will be used in this section: AS(P ρ ) Aggregate score for P ρ AS* Maximal aggregate score for all routes in R D(e) Length of link e in meters D(P ρ ) Length of P ρ in meters G Graph of connections G + Graph G extended by tessellation nodes GM Geometric mean H(v) Hazard proximity number for node v H i (v) Hazard proximity number for node v w.r.t. z i HD(e) Hazard proximity number for link e P ρ (p,q)-route for the propagation coefficient ρ PI(P ρ ) Proximity index for P ρ R Set of (p,q)-routes S D (P ρ ) Distance score for P ρ S PI (P ρ ) Proximity index score for P ρ WGM Weighted geometric mean b(v,z i ) Minimum number of obstructions between v and z i d(v,z i ) Direct distance from v to z i in meters d + (d -) Maximal (minimal) lengths of routes in R l(P ρ ) Number of links in P ρ p + (p -) Maximal (minimal) proximity indices of routes in R r(e) Proximity ratio for link e r i (e) Proximity ratio for link e w.r.t. z i t Hazard tolerance, t = 0, 0.5, or 1 z i ith epicenter ρ Propagation coefficient ρ max Maximal propagation coefficient τ Adapted propagation coefficient, τ = 1 + ρ/100

Safest routes
The first objective is to find the safest route in a building, that is, the total hazard proximity from a route to all the epicenters of the extreme event should be minimized. As explained above, the hazard has two criteria: distance and the number of obstructions. These criteria are very natural and should be considered together. For example, the distance of 40 m between a person and an epicenter of one of the aforementioned extreme events in a direct visibility might be considered not as safe as the distance of 30 m with two concrete walls between the person and the epicenter. It may be pointed out that we consider a generic extreme event.
For a particular event, for example, fire, our algorithms can be further developed by taking into account more precise models for heat propagation, smoke spread, and structural collapse. Thus, we consider an extreme event with k epicenters and two criteria for safety: distance and the number of obstructions. Based on these criteria, Algorithm 1 finds the safest available (p,q)-route in a building. Note that the application of this algorithm is threefold: it can be used by a rescue team to get from one of the entrances to a particular place in the building; as an evacuation algorithm from room p to one of the exits; or for navigation from point p to point q in the building.
Let d(v,z i ) denote the direct distance from node v to the epicenter z i in meters, and b(v,z i ) the minimum number of obstructions between v and z i . At the first stage of Algorithm 1, the direct distances and the number of obstructions are calculated in the DB Procedure. More precisely, using the standard geometric technique, the direct distance d(v,z i ) from node v to the epicenter z i in meters is calculated for all nodes vϵV(G + ) and all epicenters. Then, the Breadth First Search Algorithm is run in the graph G from the node z i . It finds the number of links in shortest paths from z i to all other nodes in G. Since each link in those paths represents a physical obstruction, the number of links in such (z i ,v)-path is the minimum number of obstructions (i.e., walls and floors) between z i and v. This number is denoted by b(v,z i ), and it should be further adjusted for doors and exit doors. For example, if v represents a door with two adjacent rooms r 1 and r 2 , then we put b (v, z i ) = min{b (r 1 , z i ) , b (r 2 , z i )} because the links vr 1 and vr 2 do not represent a physical obstruction. Note that for some nodes, such as corner nodes and the "new" nodes in the tessellation, the numbers b(v,z i ) have not been calculated because those nodes are not in the dual graph G. Hence, for each such node vϵV(G + ) -V(G) we put b(v,z i ) = b(w,z i ), where the node wϵV(G) represents the cell whose tessellation contains v.
It may be pointed out that the calculation of the parameters b(v,z i ) is only possible because the dual graph incorporates all the necessary 3-D information about the building. The second step of Algorithm 1 is to compute the hazard proximity numbers for all links in G + for the given propagation coefficient ρ = ρ max , where ρ represents the degree of hazard propagation, it is explained in more detail below. This is done in the HP Procedure. Initially, for each epicenter z i , the hazard proximity numbers H i (v) are calculated for all nodes in G + . These numbers go from small positive values, which mean "very far from the epicenter," to 100, which stands for "inside the epicenter z i ." The hazard proximity numbers are based on the parameters d(v,z i ) and b(v,z i ), which should be first replaced by one variable representing their average. For two variables with different numerical ranges it is appropriate to use a (weighted) geometric mean. Since d(v,z i ) and b(v,z i ) have different ranges, the weighted geometric mean is applied: where w 1 and w 2 are the relative weights for the direct distance and the number of obstructions. Because the direct distance is as important as the number of obstructions, we can assume that the corresponding weights for the two variables are in proportion 50:50, that is, w 1 = 0.5 and w 2 = 0.5. However, these weights can be adjusted if necessary, for example, for buildings with many large open spaces and few obstructions. Thus, the above formula is simplified to the standard geometric mean: Next, the values of geometric means should be transformed to the scale going from 100 to 0 taking into account the propagation coefficient ρ ࣙ 0. This is achieved by using the following formula: Now, if we denote τ = 1 + ρ 100 , then a well-justified formula for the hazard proximity numbers is obtained: For instance, if ρ = 100, then hazard proximity numbers for nodes propagate quickly from 100 (in the epicenters) to small positive numbers (far from the epicenter). This puts a strong emphasis on the epicenter and the rooms in its close proximity. In contrast, if ρ is a small positive number, then the propagation is slow, thus putting less emphasis on the epicenter and the nearby rooms. In the extreme case ρ = 0 there is no propagation, that is, all hazard proximity numbers for nodes are equal to 100.
Having calculated H i (v) for all the nodes in G + and all values of i = 1,2, . . . ,k, the following formula is used to compute the final hazard values for nodes: Here we assume that the hazard at a particular node is equal to the maximal hazard proximity number at this node for all the epicenters, this approach is justified for many cases. A different formula can be easily incorporated in the algorithm if it is necessary, for example, to take into account the cumulative effect of all the hazard proximity numbers H i (v). Further, for each link e = uv in the graph G + , the hazard proximity number HD(e) for e is determined by calculating the arithmetic average of the hazard proximity numbers of its end-nodes, and then by multiplying the resulting number by the length of e in meters:

H D(e) = 0.5[H (u) + H (v)]D(e)
The last operation is important because G + is not a homogeneous network. For example, let us suppose that one link is 2 m long and another is 10 m long, and they both have the same hazard proximity number, say 10. If they both are used in a navigation route, then it is natural to assume that a travel time for a longer link would be approximately five times longer, so the hazard proximity number for the longer link should be 50. In other words, if we subdivided the longer link into five 2-m-long links to make the network more homogeneous, then those five links would approximately contribute 50 to the total hazard proximity of the route.
The final stage of Algorithm 1 is to run Dijkstra's algorithm in the graph G + from node p with link weights HD(e). It produces the safest available (p,q)route P in the building with two criteria: distance and the number of obstructions. Note that the formula for the hazard proximity numbers is based on the weights w 1 = 0.5 and w 2 = 0.5 for d(v,z i ) and b(v,z i ). However, different weights can be easily incorporated in the formula.

Balanced routes
The second objective is to produce a balanced (p,q)route in a building, where an extreme event is occurring. This is achieved in Algorithm 2, where one of the input parameters is the hazard tolerance coefficient t.
The hazard tolerance is a trade-off between distance and safety, and it can be equal to 0, 0.5, or 1. For example, if t = 1, then the shortest route will be generated. If t = 0, then the algorithm finds the safest available route. If t is not specified and there are enough routes in the set R, then a route with the 50/50 balance of distance/safety will be reported as the balanced route.
Typically, there are many (p,q)-routes. The shortest route might go through the epicenter and be dangerous, whereas the safest route might be the longest one. A member of the rescue team, who is fully protected from the hazard, may wish to use the shortest route even if it is the most dangerous one, that is, their hazard tolerance t is equal to 1. In contrast, an unprotected person with a respiratory disease may want to use the safest evacuation route, whatever is its length, in which case the hazard tolerance t is 0. The hazard tolerance is an optional parameter, and at the moment there are only two values. If it is not given, then by default t = 0.5. The default value of 0.5 simply means that the hazard tolerance has not been specified, and this number will be used as a relative weight for the distance attribute, thus the relative weight for the hazard proximity will be 0.5 too. We use one variable t in this context because the single formula for the aggregate score AS(P ρ ) will be applied for all values of t. Thus, if t is not specified, then a route with a right balance of distance/proximity will be chosen.
The first part of Algorithm 2 runs the DB Procedure to determine all parameters d(v,z i ) and b(v,z i ). Then, the binary search is carried out with respect to ρ. The first run is for ρ = 0, producing the shortest (p,q)route P 0 because all hazard proximity numbers for links are 100D(e). The route P 0 is included in the set R. The next run is for ρ = ρ max . If the resulting route coincides with P 0 , then there is no interval for the binary search and it is terminated. Otherwise, the route is different from P 0 and it is included in R. The next run is for ρ = 0.5ρ max . There are three possibilities here. If the resulting route is a new one, then it is included in R and the binary search continues for two intervals (0; 0.5ρ max ) and (0.5ρ max ; ρ max ). If the resulting route coincides with one of the routes in the set R, then one of the intervals is removed from the search and the other interval is used in the binary search. For example, if the route coincides with P 0 , then the binary search continues for the interval (0.5ρ max ; ρ max ), whereas the interval (0; 0.5ρ max ) is removed. This procedure is terminated if at least one of stopping criteria is satisfied: a specified size of R, a specified length of the widest interval, and a running time. For our test buildings, the procedure goes on until seven routes are found or the length of the widest interval is less than 0.1.
In the next block of Algorithm 2, the total distance of each route in R is calculated. Then, the proximity ratios are computed for each link in a route. They are based on the geometric averages of parameters d(v,z i ) and b(v,z i ) for end-nodes of the link and its length. In contrast to hazard proximity numbers, the proximity ratios do not depend on the propagation coefficient ρ, and a small proximity ratio means a close proximity to one of the epicenters. For a given link, the final proximity ratio r(e) is the smallest proximity ratio for that link: r (e) = min 1<i<k r i (e). The proximity index for a route P is the harmonic mean of r(e)'s for all links in the route: where l(P) is the number of links in P. Note that the proximity index is an average of rates. Also, the proximity index should not be dominated by sections of a route with large proximity ratios, and actually the impact of small proximity ratios is important. Therefore, the harmonic mean is an appropriate measure for the proximity index. Since the proximity index is independent on the propagation coefficient, it can be used for comparison of the routes from the set R.
In the final part of the algorithm, a multiattribute decision-making technique is used to rank the routes in R and choose a balanced (p,q)-route. First of all, the maximal and minimal values of the lengths and proximity indices are calculated for all routes in R: Then, quadratic value functions are applied for rating the routes with respect to two attributes, the distance and the proximity index. Different value functions were tested, and it turned out that the most appropriate one is quadratic. For each route in R, the scores for these Finally, the aggregate score is calculated as a weighted average of the routes' rates, where the weights depend on the tolerance coefficient t. More precisely, the weights are t and 1 -t. The balanced (p,q)-route is one with the highest aggregate score. For example, if t = 1, then the shortest route is chosen. If t = 0, then the algorithm returns the route with the highest proximity index, that is, the safest available route. If t is not specified and there are enough routes in the set R, then a route with the 50/50 balance of distance/proximity will be reported as the balanced route.

TESTING THE ALGORITHMS
We start testing Algorithms 1 and 2 with a virtual building shown in Figure 8. This 10-floor building has five stairwells and five exits, and at each level there are five rooms connected by a long corridor as can be seen in   Figure 9. The relatively large number of stairwells is needed to illustrate the behavior of the algorithms.
In what follows, we put ρ max = 100 for Algorithms 1 and 2. However, further testing is needed to decide which values of the maximal propagation coefficient are appropriate for different buildings. The epicenter of an extreme event (labelled by a star) is on the fourth floor, as can be seen in Figure 10. The starting point (the node p) is located on the top floor, above the epicenter, and the node q is not specified. Thus, we are looking for a route from p to one of the exits. It is not difficult to see that Algorithm 1 is a particular case of Algorithm 2 if we put t = 0 in the latter, that is, the former finds a route with the highest proximity index. Thus, it is enough to test Algorithm 2 for different values of the hazard tolerance t.
For the above scenario, the binary search of Algorithm 2 produces five different routes with the following propagation coefficients: ρ = 0, 7, 10, 11, 12. The lengths and the proximity indices for those routes are summarized in Table 2. As can be seen in the table, Algorithm 2 returns the shortest route P 0 if t = 1, and the safest route P 12 if t = 0. If t is not specified, then the balanced route P 10 is returned by the algorithm.
The route P 0 is shown in Figure 10a. It goes through the first stairwell, which is very close to the epicenter, so it is the most dangerous route, but the shortest one. The second route P 7 goes through the second stairwell; Fig. 11. The scatter plot for routes P 0 -P 12.
it is safer but longer. The routes P 10 , P 11 , and P 12 go through the third, fourth, and fifth stairwells, respectively. The routes P 10 and P 12 are shown in Figures 10b and c. The scatter plot for the two parameters of the five routes is given in Figure 11. It is not surprising that there is a very strong positive correlation (at 1% significance level) between distance and proximity index. Also, the routes form a so-called "efficient frontier" in the sense that no route is "dominated" by another one, that is, for any two routes one of them is safer but longer.
Different floors for the starting room have been tested. In general, Algorithm 2 produces good results; however, in some cases there is an unexpected behavior. For example, if the starting point is located on the sixth floor, then the binary search finds five routes, see Table 3. The first four routes are similar to the results for the top floor; they go through the first four stairwells, respectively. However, the fifth route P 26 does not go directly to the fifth stairwell. According to the algorithm, it is safer to first go upstairs and then use the fifth stairwell as can be seen in Figure 12. This is reflected in the corresponding proximity indices: 3.87 for P 26 and 3.59 for the route that goes directly to the fifth stairwell.  Note that P 26 is much longer compared to other routes, so the balanced route P 25 might be considered as a more reasonable one.
The aforementioned behavior is not necessarily unreasonable; it depends on the type of hazard. For some extreme events it might be deemed as safe, for others, for example, fire, as unsafe. In the latter case, this problem can be rectified differently. The first approach is to include another criterion, route complexity, which will be investigated in a separate article. The complexity of the route P 26 would be rather high because it goes upstairs and uses two staircases, thus decreasing the like-lihood that it will be eventually chosen. Another approach is to use a better model for hazard propagation in a building if the nature of hazard is known, however, such models are out of the scope of this article where we consider a generic extreme event.
For the time being, we can use a simple approach based on a binary input variable x. It is equal to 1 if going upstairs is undesirable; and 0 otherwise. By default, x = 1, in which case we add an additional "penalty" for going upstairs in terms of distance. This penalty increases the hazard proximity numbers for links representing sections of a staircase that go up, thus making them undesirable in the routes. Note that such links are not forbidden completely, because in some cases going upstairs is unavoidable. This adjustment of Algorithm 2 produces the route P 29 instead of P 26 . The new route has a slightly worse proximity index (3.59 versus 3.87), but it is much shorter (64.2 m versus 93.1 m) and does not go upstairs. The adjusted algorithm returns P 22 as a balanced route because the "outlier" P 26 was replaced by a much shorter route P 29 .
Let us now consider the building of Figure 4, which has three stairwells. Instead of looking at simple situations with one or two epicenters, we simulate an extreme event with four epicenters as illustrated in Figure 13. This is an extremely tight situation, because the west stairwell is blocked by two epicenters at the third and seventh floors, whereas another two epicenters are located in the east and north stairwells at the fifth floor. Thus, to avoid the epicenters on the fifth floor, one has to use the west stairwell, which is not safe either. This means that a very safe route from the top floor to one of the exits does not exist, that is, any route is very close to the hazard in such an extreme configuration of epicenters. Note that the hazard proximity numbers for rooms, which are shown in red/orange/yellow hues, correspond to the hazard proximity numbers for nodes in the dual model.   The shortest, balanced, and safest routes produced by Algorithm 2 are shown in Figure 14, where stars represent the epicenters, and the information about the routes is summarized in Table 4. Note that ρ max is now 200 because there are many epicenters. As can be seen in Figure 14, the shortest route goes through two epicenters in the west stairwell, the balanced route uses the north stairwell with one epicenter, and the safest available route tries to stay further away from the epicenters by first using the east stairwell down to the sixth floor, then the west stairwell down to the fourth floor, and finally the north stairwell down to the ground floor.
The above extreme example makes it obvious that not only the global hazard proximity of a route is important, but also a local proximity of route's nodes/links to the hazard should be taken into account. Indeed, the proximity index is a global parameter, which can be used to compare different routes. However, it does not tell us that the route is going through or very close to one of the epicenters. Some threshold values for the proximity ratios r(e) could be used for this purpose if an epicenter is not located in a corridor or a large open space. Otherwise, the direct distance d(v,z i ) from v to the epicenter could be a better measure of local hazard proximity of the node v. The assessment of route's local hazard proximity is out of the scope of this article. However, in the context of the above example, the proximity ratio less than 0.5 means an extreme proximity to an epicenter, whereas the proximity ratio between 0.5 and 1 represents a close proximity. Thus, if a user opts for    a balanced route in our example, then (s)he should be given the information that it is extremely close to the epicenter together with the option to choose the safest available route.

Courtesy of MZ&Partners Architectural & Engineering Consultancy
Finally, let us consider a more realistic example, which is based on the Doha World Trade Centre (DWTC). A typical actual floor of this building is shown in Figure 15, and its DHE reconstruction and tessellation are illustrated in Figure 16. The left and right pictures of Figure 17 demonstrate the shortest and safest routes between two rooms, respectively. Further, the 37-floor building of Figure 18 was generated using the aforementioned floor plan. To add even more complexity, the locations and the number of staircases were modified: there are three staircases between the ground floor and the third floor; three staircases at different locations between Floors 5 and 36; and six staircases on Floor 4 (three going up and three going down). The resulting 37-floor building is not exactly the DWTC, but it is based on the actual floor plan of the DWTC. Figure 18 illustrates the shortest and safest routes in this building, where stars represent three hazard epicenters.
As mentioned above, the algorithms can be used for navigation from point p to point q in the building; this is illustrated in Figure 17. Also, they can be used for evacuation from room p to one of the exits or by a rescue team to get from one of the entrances to a particular place in the building. This is illustrated in Figure 18. From the practical viewpoint, the 3-D model of the building can be created in advance and kept in the cloud. On arrival, the rescue team should detect the event epicenters. For example, in case of fire, the team can scan the building with appropriate equipment or use temperature sensors in the building if they are available. Then the model in the cloud is updated with epicenters and the safest available route is found. Notice that people inside the building will also have access to the updated model in the cloud, and there is a potential for using an indoor navigator.

CONCLUSIONS
In this article, we presented an algorithm for finding the safest route in a building, where an extreme event with many epicenters is occurring. Another algorithm produces a balanced route, in which a trade-off between route's length and hazard proximity is made. Note that the hazard proximity has two criteria: distance and the number of obstructions between the epicenter and a point in the building. The proposed algorithms are based on the BIM-GIS model developed by Boguslawski et al. (2015), and they essentially use the underlying 3-D structure of the model. For example, the calculation of the parameters b(v,z i ) is only possible because the dual graph incorporates all the necessary 3-D information about the building.
As illustrated above, the application of the algorithms is threefold: they can be used by a rescue team to get from one of the entrances to a particular place in the building; for evacuation from room p to one of the exits; or for navigation from point p to point q in the building. In addition to the 3-D BIM-GIS model, the algorithms are based on multiattribute decision-making technique and the introduced hazard proximity numbers, hazard propagation coefficient, and proximity index for a route. The formulae for these parameters are well justified and they have been validated by testing. Also, they can be easily adjusted to incorporate different weights for d(v,z i ) and b(v,z i ).
The results of testing are promising; in many instances the algorithms produce very reasonable balanced and safest routes. However, in some cases the safest available route is rather long and it goes upstairs. One possible extension to overcome this issue is to include another criterion, route complexity, which will be investigated in a separate article. Another approach is to use a better model for hazard propagation in a building if the nature of the hazard is known. For example, in the case of fire, more precise models for heat propagation, smoke spread, and structural collapse can be used. A further limitation is that the algorithms do not take into account multiple agents that create route conflicts and congestion in a building. However, the capacity constraints can be easily included in the 3-D model, so there is a potential for extending the algorithms in this direction.
In the future work, we will consider three criteria: distance, hazard proximity, and route complexity. The multiattribute rating technique will be applied for finding the "best" route, that is, a route which is reasonably short, safe, and simple. It may be pointed out that the information about building material may be taken into account in the corresponding algorithm, as well as available sensor information (if any), for calculating the best egress route. Also, if distribution of people in a building is known, or can be predicted, then the distance criterion may be replaced by the time criterion thus taking into account multiple agents in the building. Another interesting extension would be to take into account the dynamics of the situation, which is particularly important if a hazardous event develops rapidly.