In this paper, a bilevel multiobjective optimisation model is proposed to solve the evacuation location assignment problem. The 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. 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]. Location-allocation models have been proposed for prescribing evacuation plans during hurricane occurrences [57]. 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], 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].

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 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. 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 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. The lower-level program is modelled under the assumption of user equilibrium [42] and is formulated as a single-commodity 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 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 in the network model, as shown in Figure 1. The destination for all evacuees in the single-commodity 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 upper-level decision-maker, 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 lower-level 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 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. The objective function of the lower-level 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 lower-level 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 single-level representation is highlighted, making use of the convexity of the lower-level program. The 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.

3.1. Model Reformulation
3.1.1. Equivalent Lower-Level Model

The 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]. 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 lower-level 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 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.

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 piece-wise 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. Single-Level MILP

The final linearised single-level 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 single-objective optimisation problem will conflict with one another. The 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 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.

end for
Lexicographic minimiser:

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 1st stage. The 2nd 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 free-flow 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 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.

4.1. 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. 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%.

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 decision-maker, 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 1st 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 single-level one. The resulting MILP was then solved using an off-the-shelf 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 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. 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.