A Bilevel Multiobjective Optimisation Approach for Solving the Evacuation Location Assignment Problem

In this paper, a bilevel multiobjective optimisation model is proposed to solve the evacuation location assignment problem. +e model incorporates the two decision-makers’ spaces, namely, urban planners and evacuees. In order to solve the proposed problem, it is first reformulated into a single-level problem using the Karush–Kuhn–Tucker conditions. Next, the problem is linearised into a mixed-integer linear programming model and solved using an off-the-shelf solver. A case study is examined to showcase the applicability of the proposed model, which is solved using single-objective and multiobjective lexicographic optimisation approaches. +e model provides planners with an ability to determine the best locations for placement of shelters in such a way that the evacuees’ traffic assignment on the existing network is optimised.


Introduction
Failure in infrastructure and buildings during disasters is an important phenomenon to understand when attempting to investigate the response of masses of people to the occurrence of extreme events. Devising evacuation strategies that are effective enough to handle the large masses required to be moved from one point to another during extreme events is presented as an avenue to manage disasters. e complexity of managing an evacuation plan lies in mapping the underling factors that need to be accounted for when forming new strategies. Some of these factors include the behaviour of the evacuees, the response of the underlying structure that needs to accommodate the masses being evacuated, the area that has been impacted by the extreme events, and the availability of shelters that can be used as temporary rescue facilities [1,2].
Studies that optimise the evacuation process through mathematical programming are wide and varied [3,4]. Location-allocation models have been proposed for prescribing evacuation plans during hurricane occurrences [5][6][7]. e evacuation involved in city emergency disasters has also been modelled as a shortest path network flow problem [8,9]. Acts of intention or natural disasters occurring in confined spaces such as stadiums, museums, and shopping centres have also been examined through proposing evacuation optimisation routing problems [10]. Identification of optimal routes and shelter locations during fire disasters is also a topic that has been previously explored [11].
In the literature, studies that investigate disaster management can be classified according to the planning stage being examined. In the predisaster phase, the emphasis is on planning for handling the extreme events, including the focus on building and infrastructure reenforcement [12,13].
is also incorporates the prepositioning of fast relief distribution centres [14] and locating evacuation shelters [15].
In the postdisaster stage, the literature has focused on the evacuation process of the masses to shelters, the distribution of relief, and casualty transportation [13]. In terms of mass transport to shelters, some of the common models proposed include linear programs [16], mixed-integer programs [17,18], cell transmission models [19,20], location routing [17], classical vehicle routing [21], and dynamic network flows [22]. Objectives that are optimised in these models include total travel time [18,19], total evacuation time [21], total travel distance [23], waiting time of evacuees [24], and cost of travel flow [25,26].
is paper proposes a novel approach where a model is formulated to address both the predisaster phase and the postdisaster phases simultaneously, by optimising the positioning of shelters in an urban region and by handling the mass evacuation that occurs after disaster through solving a traffic assignment problem. e novelty of the proposed approach presented in this paper lies in the consideration of the impact of traffic on the movement of people during the evacuation process, through presenting the problem as a bilevel model [39]. In addition, the problem of locating shelters is integrated with the problem of allocating people to designated shelters that have been placed in the region. In order to model a realistic problem, several factors are considered when it comes to planning the pre-and postdisaster stage. is is achieved by modelling multiple objectives at the 2 decision-making levels. Within the upper level decision-making, the objectives that are modelled include the minimisation of the construction costs of shelters, and the maximisation of coverage of these shelters to the surrounding residential zones, through minimising the total system travel time on routes during the evacuation process. At the lower level, the individual travel times of the evacuees on the network are minimised, considering equilibrium conditions in the overall system, where users can no longer reduce their travel efforts.
A multiobjective [40] bilevel model is thus presented in this paper, where the problem of locating shelters and the routing of traffic that results due to the occurrence of an extreme event is solved using a mathematical optimisation approach. e bilevel model attempts to capture the leader and the followers' decision space and response. e leader in the problem examined will be the authorities responsible for the urban planning of the region. In response to the decisions made by the urban planners on locations for shelters, the followers will attempt to optimise their decisions in terms of choosing travel routes with least traffic delays.

Evacuation Location Assignment Problem
In this part of the paper, the multiobjective bilevel programming model for addressing the evacuation location assignment problem (ELAP) is presented. In particular, the adopted notation and the network representation of the model are defined. e ELAP considered lends itself to a class of optimisation problems where its representation can be modelled as a leader-follower or Stackelberg game. Such problems fall into the field of bilevel programming, their main characteristic being that they contain an optimisation problem that is embedded into another one [41], resulting in a hierarchical representation of the decision-making ladder.
Given that the proposed model is split into two levels, namely, an upper level representing the urban planners' decision when it comes to locating the evacuation shelters and a lower-level program related to the formation of the objective and decision space of evacuees on the travel networks linking the underlying region considered, where each of these will be described separately. e lower-level program is modelled under the assumption of user equilibrium [42] and is formulated as a single-commodity flow problem [43]. us, each evacuee is modelled as having the same destination on the network (i.e., all are moving towards the shelters).

Notation and Network Representation.
In this section, the representation of the problem in the form of a transport network is discussed. e network is composed of several types of nodes: the first type of node represents the residential zones in the region; the second type of node represents the potential areas for locating the shelters; the last type of node is a dummy node representation that acts as a sink for all travel towards the shelters in the region. In order to construct the travel network that underlies the region being considered within the evacuation planning problem, geographic information system (GIS) can be utilised. In particular, the use of GIS in evacuation planning has gained extensive popularity in recent years [44,45].
Given that the location of the shelters is a decision variable which is unknown in the proposed evacuation problem, the flow of evacuees towards the shelters cannot be referenced via a specified destination node (as this is unknown at the start). Trips that evacuees need to make to shelters are considered as the single-commodity modelled in the traffic assignment problem of the follower (lower level). As a result, a dummy node acts as a sink for all travel that is heading towards any of the shelters that are yet to be placed; the dummy node is denoted as S in the network model, as shown in Figure 1. e destination for all evacuees in the single-commodity flow problem modelled is thus specified as the dummy node S. Figure 1 also displays the notation description and set adopted in the proposed model.
In order to present the notation adopted in the proposed model, Table 1 is produced.

Upper Level: Leader Model.
e upper level of the bilevel model focuses on capturing the decision by urban planners in selecting appropriate locations for the shelters. A multiobjective model is formulated, where the focus of the urban planners in terms of positioning the shelters is achieved via consideration of 2 main objective functions. e first objective function modelled, Eq. (1), minimises the construction costs associated with building the shelters. e second objective function maximises the coverage of the positioned shelters to surrounding residential zones from which evacuees will be leaving towards the shelters during disastrous events; this is done through minimising the total disruptions of the underlying travel network during the evacuation period to ensure su cient coverage: Within Eq. (1), the cost of construction of each shelter will depend on the region in which the shelter is to be constructed. All shelters are assumed to be of the same type and so the construction cost parameter C i is independent of the shelter type. e binary variable z i is used to indicate the location of shelter i.
Equation (2) minimises the total system travel time (TSTT) of the whole network [46,47]. TSTT re ects the overall travel time spent by network users (evacuees) during the evacuation process, traveling on the links of the network; it is calculated by multiplying the tra c on a given link of the network heading towards the evacuation sink, x ij by the time function, t ij (x ij ). TSTT accounts for tra c congestion, induced by the shelter layout adopted on the network. Eq.
(2) thus acts as a proxy of the level of coverage of the shelters in the region. It is important to note that since this article examines an evacuation problem, the tra c assignment modelled herein is only associated with evacuation ow, with no inclusion of background tra c. e number of shelters to be positioned in the region is not prespeci ed and is left to the model to optimise. An upper bound however is set as the maximum number of nodes available on the network for the placement of the shelters. e nature of the objective functions modelled is such that a con ict between the optimised solutions for each separate function is likely to arise. For instance, the solution that minimises construction costs will be the one associated with the lowest coverage of the shelters to the evacuees requiring to use the shelters. A multiobjective optimisation approach is thus required to handle the two formulated functions.
ere are no constraints formulated for the upper-level decision-maker, apart from the de nition of the domain of the variables controlled by the leader in the Stackelberg game modelled, which is de ned via Eq. (3). e layout of shelters adopted will impact the tra c conditions due to movement of evacuees on the network links. is is because changes to the shelter locations in the urban region examined will create a change in the ow patterns and congestion levels induced by these facilities. Such decisions are re ected in the follower's decision space, Binary variable, which equals 1 if shelter is placed at location p ∈ P and 0 otherwise Travel cost variable, de ned as a positive and continuous function of total ow on link Auxiliary binary variables tij, λ h ij Auxiliary continuous variable Advances in Civil Engineering modelled as the evacuees in this Stackelberg game. In the next section, a lower-level model is formulated for determining the routing of the flow of evacuees through the transport network.

Lower Level: Follower Model.
Each user within the transportation network will attempt to reduce their individual travel time to counter the impacts brought on by the introduction of the demand-inducing sheltering facility to the underlying region. To model this behaviour of evacuees, a user equilibrium (UE) traffic assignment model is formulated as the lower-level model, based on Wardrop's first principle [42]. In particular, equilibrium in the system is achieved when the link flow pattern is such that travel time on all used paths connecting the origin and destination nodes is less than or equal to travel time on the unused paths. e objective function of the lower-level model incorporates a link cost function, t ij (x ij ), dependent on the flow of the network, and it aims at minimising the total of travel time on all the links of the network. In the bilevel model, the link cost function adopted is the travel time function developed by Bureau of Public Roads (BPR) [48], which is given in the following equation: where T 0 ij denotes the free flow travel time, while k 0 ij denotes the existing capacity of link (i, j), respectively. e proposed lower level program of the bilevel model is defined by the following equations: subject to x iS ≤ z i e∈E q eS ∀i ∈ P, e objective function of the lower-level model, Eq. (5), is based on minimising the individual user's travel times, t ij (ω), on all links of the network. Eq. (6) ensures flow conservation at each node of the network, by making use of the flow variable x ij . In Eq. (6), the parameter q iS indicates the demand flow, originating at node i and heading towards the sink node S. Equation (7) ensures that flow from potential shelter locations to the sink node, x iS occurs only in the case that a shelter is placed at that node, i.e. z i � 1. Eq. (8) sets the summation of all flow that moves from the potential shelter locations towards the sink node, p x pS , to be equal to the total demand for shelters by evacuees throughout the region, e q eS . Eq. (9) defines the domain of the follower's decision variables.

Solution Procedure
In this section, an explanation of the methodology adopted to solve the proposed bilevel model is provided. e application of the Karush-Kuhn-Tucker (KKT) conditions to reformulate the model into a single-level representation is highlighted, making use of the convexity of the lower-level program.
e single-level mixed-integer nonlinear programming (MINLP) model is then linearised, through implementing a scheme that is based on piece-wise approximation of the convex BPR function. Finally, a multicriteria analysis that can be applied to deal with the multiple objectives defined for the ELAP is discussed.

Equivalent Lower-Level Model.
e user equilibrium conditions of the lower-level program can be represented by a set of first-order equivalent constraints, namely, the Karush-Kuhn-Tucker (KKT) conditions, as described in [39]. e dual variable μ i is defined for Eq. (6) and can be thought of as the minimum travel time required to evacuate people from node i. e variable t ij captures an approximation to the time function t ij (x ij ).
Equations (5)-(6) and Eq. (9) of the lower-level program can then be replaced by the following equations: e complementary slackness conditions of KKT are enforced by Eqs. (10) and (11). Note that the constraints of the lower-level model, Eqs. (12) and (13), are defined as part of the KKT conditions.
Since the complementary condition, Eq. (11), is nonlinear, the single-level model cannot be solved using a linear solver. To address this, an appropriate linearisation scheme to reformulate the latter equations needs to be applied, as demonstrated in the next section.

Linearising the KKT Conditions.
Let ω ij be an auxiliary binary integer variable, which equals 1 if t ij − μ i + μ j � 0 and 0 otherwise. Eq. (11) is replaced with the following set of 4 Advances in Civil Engineering constraints, Eqs. (15)- (17), resulting in the linearisation of the complementary conditions: where Q and Q ′ are the large positive constants.

Linearising the BPR Function.
Given that the BPR function is nonlinear, a chain of linked special ordered set (SOS) conditions is implemented for the linearisation procedure [49]. e principle idea behind the linearisation scheme adopted is shown in Figure 2. e feasible domain of x ij is partitioned into h ∈ H segments and a continuous nonnegative variable, namely, ψ ijh is associated with each segment. As shown in Figure 2, each segment represents a straight line between two points on the BPR curve. At the start/end of each segment, a predetermined ow is de ned by the variable A ijh . e ow variables can then be represented by the following equation: e BPR function is approximated by the following equation: e conditions imposed on ψ ijh are given by the following equations: Equation (21) is the usual convex combination requirement in piece-wise linear approximation. Eq. (22) states that the variable ψ ijh is of a special ordered set of Type 2 (SOS2), where a maximum of two of the latter variables, which are adjacent, can be nonzero [50]. Most commercial solvers can handle the SOS2 variables; however, in order to improve computations, the SOS2 conditions are linearised.
is is achieved by the introduction of binary variables ζ ijh which are de ned for each segment within a given link, along with a set named H h ∈ H : ord(h) < |H| { }. is set is composed of all segments in the set H, less one.
Equation (21) is replaced with the following equations: It is important to note that as the segmentation of the BPR function increases so does the accuracy of the approximation. However, this will come at the expense of higher computation costs.

Linearising the TSTT Objection Function.
To linearise the TSTT objective function, Eq. (2) can be replaced by the following equivalent equation: where μ i represents the dual variable de ned for Eq. (6) above and q iS is the total number of evacuees leaving node i. Both Eqs. (2) and (25) are equivalent, as discussed in [46].

3.5.
Single-Level MILP. e nal linearised single-level ELAP is given as the following equations: x iS ≤ z i e∈E q eS ∀i ∈ P, Advances in Civil Engineering 5 3.6. Lexicographic Optimisation. Given that the ELAP examined is multiobjective in nature, it is expected that the solutions to each single-objective optimisation problem will conflict with one another. e concept of single optimality adopted in single-objective optimisation problems is therefore replaced with that of Pareto optimality. In particular, a solution to a multiobjective optimisation problem z * is said to be Pareto optimal if there does not exist another feasible solution z such that f e (z) ≤ f e (z * ) ∀e ∈ O and f m (z) ≤ f m (z * ) for at least one index m ∈ O, where O is the set of objective functions solved in the multiobjective optimisation problem [51]. ere are multiple approaches that can be adopted to solve the multiobjective ELAP proposed in this study, including the ε-constraint approach [52], lexicographic optimisation [53], and scalarisation [54]. In this study, the lexicographic approach is adopted because it is common for urban planners to have a preference defined over the importance of the objective functions being optimised when deciding on evacuation planning approaches. When objective functions are prioritised, a unique solution of the Pareto hyper surface exists [55]. In addition, lexicographic optimisation has been implemented extensively in the literature [56,57].
Algorithm 1 displays the procedure adopted in lexicographic optimisation. A solution, which is a lexicographic minimiser of the ELAP, is the one where the objective function being minimised can only be reduced further at the expense of at least one of the higher ranked objective functions. e notation adopted for the lexicographic optimisation process is established via the term lex min[B v , B w ], where given 2 objective functions, B v and B w , the priority is given to B v which is minimised in the 1 st stage. e 2 nd stage involves minimising B w subject to the constraint B v ≤ B * v , where B * v is the optimal solution of B v obtained at the initial stage. is process is generalised for the case where more than 2 objective functions are involved.

Case Study
In order to examine the applicability of the proposed model, a realistic evacuation problem is examined. Figure 3 displays the network being considered in this study. A total of 7 residential zones are modelled in the network with 5 other potential locations for shelters identified. e total number of potential locations for positioning the shelters will therefore act as an upper bound to the total number of shelters that can be located in the region examined. Total number of nodes and links in the case considered is 12 and 30 links, respectively. e demand for shelters from each zone is given in Table 2. Because an uncapacitated model is proposed in this paper, each of the shelters positioned is assumed to be able to accommodate an unrestricted number of evacuees. Table 3 displays the free-flow travel time and capacities of each link of the network. e cost of locating each of the shelters considered is given in Table 4, as derived from local builders in Australia. e proposed model was implemented in AMPL [58] and run on a personal computer with Windows 10 as the operating system, Intel Core i7 processor and 16 GB of RAM. CPLEX is adopted as the MILP solver in this study [59]. Solving time for the case study was reported as 89 seconds. e analysis is divided into 2 main parts. e first part is associated with the examination of traffic flow on the links of the network in Figure 3 and the examination of the number of shelters placed along with the number of people served by each placed shelter.
is analysis is conducted based on optimising each individual objective function. e second part of the analysis is related to the examination of the trade-offs that exit between the 2 objective functions considered, namely, total system travel time and construction cost. In this part of the analysis, the problem is solved using the lexicographic approach described above.

Examination of Evacuees on Links and Allocation to
Shelters. In order to display the trade-off that exists between the multiple objectives solved for the ELAP, each objective is singly optimised first. e results of the optimised and evaluated functions are displayed in Table 5 and Figure 4. In particular, Table 5 reports on the optimised evaluated results of each objective function, while Figures 4(a) and 4(b) display a contrast between the optimised objective functions, and the solution yielded by an experienced planner. When the construction cost is minimised, i.e., Eq. (1), the TSTT is assessed at 5,916,700 mins. In the case that the TSTT is solely optimised, the resulting TSTT measure falls by 56% compared to the case when construction cost is optimised; additionally, the construction cost increases by 1300%. For both the construction cost and TSTT, the optimised results are always lower that the results yielded by the planner's allocation of shelters and tra c assignment; the cost drops by almost 57% moving from planned to optimised, while TSTT reduces by almost 68%. e ow and distribution of evacuees on the links of the network for each objective function that is optimised is shown in Tables 6 and 7. In particular, Table 6 shows that when the construction cost is minimised, 1 shelter is placed at Node 8 of the network. All 47,000 evacuees of the region are directed towards the 1 open shelter. Flow on links connected to Node 8 thus occupies the greatest number of evacuees. In Table 7, the optimised solution when TSTT is minimised is such that 4 shelters are activated. As a result, the 47000 evacuees are   spread amongst all 4 shelters, with the shelter at Node 8 accommodating the most evacuees, assessed at around 21,000 evacuees. Links surrounding all four shelter nodes are also carrying a large number of ow. Some of the insight that can be gained from analysing the tra c assignment of evacuees on the links of the network are as follows. First, it becomes apparent to designers, the capacity required for existing infrastructure to handle evacuees during extreme events. Second, any upgrade required for the current infrastructure can be determined based on links that occupy a large capacity of evacuees during the solving of the evacuation allocation problem. ird, the distribution of the shelters in such a manner that equity is ensured in terms of access can be determined based on proximity and availability of shelters in the region.

Pareto Optimality.
e ELAP for the case study shown in Figure 3 is now solved using the lexicographic optimisation approach displayed in Algorithm 1.
Assume that the following assignment is implemented: e ranking adopted in the solution procedure is based on As a result, emphasis is placed on increasing coverage of the shelters, as the primary goal of the decision-maker, through minimising the total system travel time on the network to avoid travel delays. is is justifiable as planners will need to ensure that all evacuees can get easy access to a nearby shelter. Second preference is then given to minimising the cost of construction.
e ELAP is hence solved over 2 stages. e results of the lexicographic optimisation runs at each stage are displayed in Table 8. As can be noticed, no change is seen between the two stages involved in the solution process.
e Pareto optimum solution corresponds to opening shelters at 8, 9, 10, and 11. e first stage of the lexicographic optimisation enforces the minimisation of the TSTT of evacuee flow on the network; this yields a shelter configuration where the total construction cost is assessed at $140,000. At the second stage, even though the cost objective is now the one minimised, due to the constraint imposed on the total TSTT of the network, no further reduction in cost is possible while maintaining the TSTT at the same value as that of the 1 st stage. As a result, a unique solution exists in this case, which corresponds to the opening of 4 shelters in the network examined.

Conclusion
A bilevel problem for evacuation planning was proposed in this paper, based on modelling planners and evacuees' decision spaces. Two objective functions were incorporated in the model, namely, a monetary cost objective for shelter construction costs and a travel system objective function that captures the overall coverage of shelters positioned in the network. Impact of induced traffic towards positioned shelters on the network, in terms of congestion levels created, was assessed via the bilevel structure of the model. A solution approach was then proposed based on a linearisation scheme integrated with KKT equivalent conditions to convert the bilevel model into a single-level one. e resulting MILP was then solved using an off-the-shelf solver for a practical case study.
e examined case revealed that the 2 objective functions resulted in conflicting solutions, with differences that were evaluated to be up to 83% in the level of the respective solutions produced. In addition, when contrasted with solutions produced by an experienced planner, the optimised model yielded results that were 57% and 68% less expensive both in terms of construction costs and TSTT. A lexicographic approach for the multiobjective optimisation of the problem revealed no differences in solutions yielded among the 2 solving stages involved, highlighting the uniqueness of the Pareto solution.
e proposed model can therefore be put into effective use for enhancing the decision-making process involved during the precrisis evacuation planning phases.
A number of limitations can be identified in this study. First, background traffic on the network is not considered. Second, enticing evacuees to adopt certain routes via reward and punishment has not been examined. ird, the nature of disaster is assumed as generic; depending on the disaster event type, evacuee's behaviour can change. Future works will involve tackling these limitations by the author.

Data Availability
Data have all been included in the manuscript.

Conflicts of Interest
e author declares no conflicts of interest. Advances in Civil Engineering 9