Abstract
In this paper, a bilevel multiobjective optimisation model is proposed to solve the evacuation location assignment problem. The model incorporates the two decisionmakers’ spaces, namely, urban planners and evacuees. In order to solve the proposed problem, it is first reformulated into a singlelevel problem using the Karush–Kuhn–Tucker conditions. Next, the problem is linearised into a mixedinteger linear programming model and solved using an offtheshelf solver. A case study is examined to showcase the applicability of the proposed model, which is solved using singleobjective and multiobjective lexicographic optimisation approaches. The 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.
1. 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. The 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]. Locationallocation models have been proposed for prescribing evacuation plans during hurricane occurrences [5–7]. The 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]. This 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], mixedinteger 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].
In terms of relief distribution, some of the model types proposed include dynamic network flows [27, 28], uncapacitated network flow [29, 30], and classical vehicle routing [31, 32]. The objectives optimised include the vehicle arrival time [33], cumulative unserved causalities [34], expected travel time [35], weighted unmet demand [27], equity of satisfied demand [36], maximal covering demand [37], and operation costs [38].
This 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.
The 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. This is achieved by modelling multiple objectives at the 2 decisionmaking levels. Within the upper level decisionmaking, 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. The bilevel model attempts to capture the leader and the followers’ decision space and response. The 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.
2. 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.
The ELAP considered lends itself to a class of optimisation problems where its representation can be modelled as a leaderfollower 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 decisionmaking 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 lowerlevel 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. The lowerlevel program is modelled under the assumption of user equilibrium [42] and is formulated as a singlecommodity flow problem [43]. Thus, each evacuee is modelled as having the same destination on the network (i.e., all are moving towards the shelters).
2.1. Notation and Network Representation
In this section, the representation of the problem in the form of a transport network is discussed. The 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 singlecommodity 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 in the network model, as shown in Figure 1. The destination for all evacuees in the singlecommodity flow problem modelled is thus specified as the dummy node . 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.
2.2. Upper Level: Leader Model
The 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. The first objective function modelled, Eq. (1), minimises the construction costs associated with building the shelters. The 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 sufficient 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 is independent of the shelter type. The binary variable is used to indicate the location of shelter .
Equation (2) minimises the total system travel time (TSTT) of the whole network [46, 47]. TSTT reflects 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 traffic on a given link of the network heading towards the evacuation sink, by the time function, . TSTT accounts for traffic 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 traffic assignment modelled herein is only associated with evacuation flow, with no inclusion of background traffic.
The number of shelters to be positioned in the region is not prespecified 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.
The nature of the objective functions modelled is such that a conflict 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.
There are no constraints formulated for the upperlevel decisionmaker, apart from the definition of the domain of the variables controlled by the leader in the Stackelberg game modelled, which is defined via Eq. (3).
The layout of shelters adopted will impact the traffic conditions due to movement of evacuees on the network links. This is because changes to the shelter locations in the urban region examined will create a change in the flow patterns and congestion levels induced by these facilities. Such decisions are reflected in the follower’s decision space, modelled as the evacuees in this Stackelberg game. In the next section, a lowerlevel model is formulated for determining the routing of the flow of evacuees through the transport network.
2.3. 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 demandinducing sheltering facility to the underlying region. To model this behaviour of evacuees, a user equilibrium (UE) traffic assignment model is formulated as the lowerlevel 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. The objective function of the lowerlevel model incorporates a link cost function, , 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 denotes the free flow travel time, while denotes the existing capacity of link , respectively.
The proposed lower level program of the bilevel model is defined by the following equations:
The objective function of the lowerlevel model, Eq. (5), is based on minimising the individual user’s travel times, , on all links of the network. Eq. (6) ensures flow conservation at each node of the network, by making use of the flow variable . In Eq. (6), the parameter indicates the demand flow, originating at node and heading towards the sink node .
Equation (7) ensures that flow from potential shelter locations to the sink node, occurs only in the case that a shelter is placed at that node, i.e. . Eq. (8) sets the summation of all flow that moves from the potential shelter locations towards the sink node, , to be equal to the total demand for shelters by evacuees throughout the region, . Eq. (9) defines the domain of the follower’s decision variables.
3. Solution Procedure
In this section, an explanation of the methodology adopted to solve the proposed bilevel model is provided. The application of the Karush–Kuhn–Tucker (KKT) conditions to reformulate the model into a singlelevel representation is highlighted, making use of the convexity of the lowerlevel program. The singlelevel mixedinteger nonlinear programming (MINLP) model is then linearised, through implementing a scheme that is based on piecewise 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.
3.1. Model Reformulation
3.1.1. Equivalent LowerLevel Model
The user equilibrium conditions of the lowerlevel program can be represented by a set of firstorder equivalent constraints, namely, the Karush–Kuhn–Tucker (KKT) conditions, as described in [39]. The dual variable is defined for Eq. (6) and can be thought of as the minimum travel time required to evacuate people from node . The variable captures an approximation to the time function .
Equations (5)–(6) and Eq. (9) of the lowerlevel program can then be replaced by the following equations:
The complementary slackness conditions of KKT are enforced by Eqs. (10) and (11). Note that the constraints of the lowerlevel model, Eqs. (12) and (13), are defined as part of the KKT conditions.
Since the complementary condition, Eq. (11), is nonlinear, the singlelevel 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.
3.2. Linearising the KKT Conditions
Let be an auxiliary binary integer variable, which equals 1 if and 0 otherwise. Eq. (11) is replaced with the following set of constraints, Eqs. (15)–(17), resulting in the linearisation of the complementary conditions:where and are the large positive constants.
3.3. 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]. The principle idea behind the linearisation scheme adopted is shown in Figure 2. The feasible domain of is partitioned into segments and a continuous nonnegative variable, namely, 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 flow is defined by the variable . The flow variables can then be represented by the following equation:
The BPR function is approximated by the following equation:
The conditions imposed on are given by the following equations:
Equation (21) is the usual convex combination requirement in piecewise linear approximation. Eq. (22) states that the variable 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. This is achieved by the introduction of binary variables which are defined for each segment within a given link, along with a set named . This set is composed of all segments in the set , 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.
3.4. Linearising the TSTT Objection Function
To linearise the TSTT objective function, Eq. (2) can be replaced by the following equivalent equation:where represents the dual variable defined for Eq. (6) above and is the total number of evacuees leaving node . Both Eqs. (2) and (25) are equivalent, as discussed in [46].
3.5. SingleLevel MILP
The final linearised singlelevel ELAP is given as the following equations:
3.6. Lexicographic Optimisation
Given that the ELAP examined is multiobjective in nature, it is expected that the solutions to each singleobjective optimisation problem will conflict with one another. The concept of single optimality adopted in singleobjective optimisation problems is therefore replaced with that of Pareto optimality. In particular, a solution to a multiobjective optimisation problem is said to be Pareto optimal if there does not exist another feasible solution such that and for at least one index , where is the set of objective functions solved in the multiobjective optimisation problem [51].
There 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.

The notation adopted for the lexicographic optimisation process is established via the term , where given 2 objective functions, and , the priority is given to which is minimised in the 1^{st} stage. The 2^{nd} stage involves minimising subject to the constraint , where is the optimal solution of obtained at the initial stage. This process is generalised for the case where more than 2 objective functions are involved.
4. 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. The 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.
The 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 freeflow travel time and capacities of each link of the network. The cost of locating each of the shelters considered is given in Table 4, as derived from local builders in Australia.
The 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. The analysis is divided into 2 main parts. The 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. This analysis is conducted based on optimising each individual objective function.
The second part of the analysis is related to the examination of the tradeoffs 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.
4.1. Examination of Evacuees on Links and Allocation to Shelters
In order to display the tradeoff that exists between the multiple objectives solved for the ELAP, each objective is singly optimised first. The 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 traffic assignment; the cost drops by almost 57% moving from planned to optimised, while TSTT reduces by almost 68%.
(a)
(b)
The flow 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 flow.
Some of the insight that can be gained from analysing the traffic 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. Third, 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.
4.2. Pareto Optimality
The 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:
The 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 decisionmaker, through minimising the total system travel time on the network to avoid travel delays. This 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. The ELAP is hence solved over 2 stages.
The 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. The Pareto optimum solution corresponds to opening shelters at 8, 9, 10, and 11. The 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.
5. 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 singlelevel one. The resulting MILP was then solved using an offtheshelf solver for a practical case study. The 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. The proposed model can therefore be put into effective use for enhancing the decisionmaking 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. Third, 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
The author declares no conflicts of interest.