Vehicle Routing Optimisation in Humanitarian Operations: A Survey on Modelling and Optimisation Approaches

The growing field of humanitarian operations is driven by frequent events of disasters seen in the world today. Within this field, Operations Research (OR) plays a critical role in alleviating the suffering of victims that are impacted by disasters. This paper focuses on the branch of a well-known OR problem, known as the Vehicle Routing Problem (VRP), within the selected scope of humanitarian operations. A total of 123 papers of the last decade are reviewed and classified under the humanitarian operations of supply and delivery, evacuation as well as rescue operations. Besides specific disaster management phases and disaster types, various modelling challenges are highlighted, hinting towards a richer and more complex VRP seen under selected model characteristic classifications. Furthermore, established solution approaches, including hybrid solutions, are highlighted and classified, discussing how they are applied in the context of these humanitarian operations. The inclusion of a machine learning solution approach under the same classification is proposed. Finally, the trend and future outlook of VRP for the suggested humanitarian operations are discussed


Introduction
Following increased events of disasters in the past decades, the urgent need for effective Disaster Management (DM) has been the focus in various academic fields. Being the centre of this topic, humanitarian operations are the key to alleviating suffering when disasters strike. This includes Humanitarian Logistics (HL) as the key in handling postdisaster distribution operations, ensuring that the critical aids are efficiently and effectively delivered to the victims [1]. Aside from HL, humanitarian operations, such as evacuation and rescue operations, are of equal importance within the scope of DM. Humanitarian operation, according to the UN Office for the Coordination of Humanitarian Affairs (OCHA), is operations conducted in order to alleviate human suffering, focusing on civilians, due to extenuating circumstances where service support is often time lacking or non-existence [2].
However, existing theoretical approaches for optimization methods in humanitarian operations rarely account for the real situation of the event, which often includes the disruption of critical facilities, such as the transportation network consisting of the roads and highways. Moreover, chaotic situations could arise due to an uncertain volume of demands and critical conditions that victims may suffer from, such as being trapped in a damaged building or that the road may become congested due to unplanned panic evacuation. This gives rise to various theoretical approaches within the field of Operations Research (OR) in addressing the humanitarian operations in such challenging settings. One of the prominent scopes within the field of OR is the Vehicle Routing Problem (VRP), where the effort to address the humanitarian operations could be traced back as early as 1988 with the work of [3].
Despite almost three decades of research, deterministic problems continue to dominate this field of research, despite the advancement of solution approaches, such as the exact solution, heuristic and metaheuristic solution approachesm as well as various hybrid solution approaches. However, the deterministic problem is not practical for problems that involve uncertain parameters and periodic changes of parameters as one would perceive during the time of disaster. Such a real-world problem demands that the dynamic and stochastic elements of the problem be addressed [4]. Nonetheless, limited literature addressed the rich VRP (RVRP) for the humanitarian operations that, often, would include the stochastic and dynamic elements as well as other attributes, such as multi-objective, multi-trip operations, split deliveries, and time windows, among others.
However, it is challenging to apply exact solution approaches for RVRP that best describes the humanitarian effort operations. This stems from the fact that even the basic VRP is an NP-hard problem, especially for practical problem instances. The same problem is faced by Machine Learning (ML) that is based on the dynamic programming approach that suffers the three curses of dimensionality [5]. At this point, the computation for the exact solution approach or dynamic programming to obtain one true optimal solution would become prohibitive and, eventually, a trade-off needs to be made between the quality of solution and computation ability. For such cases, most researchers and practitioners would turn towards heuristic or metaheuristic solution approaches, where a near-optimal solution could be obtained within a reasonable computation time. Instead of iteratively searching for all potential solutions and evaluating them, as would be done in exact solution approach and dynamic programming, heuristic and metaheuristic work by attempting to avoid evaluating all of the potential solutions in the solution pool while using a controlled randomised solution selection and try to reach a near-optimal value by iteratively improving the current potential solution until no improvement is detected within a reasonable computation time. Through this principle, the exploration and exploitation of potential solutions play an important role in reaching the near-optimal solution. These approaches have been popular ever since, due to the applicability of the solutions in real-world practice in resulting in satisfying improvements.
In the world of DM, the real world VRPs are often more complex and challenging to be solved by any conventional means available. For such a problem, uncertainty and the dynamic of parameters are the main ingredients in addressing a complete and realistic model. ML approach provides such a platform by naturally incorporating dynamic and stochastic elements of the problem through the Markov Decision Processes (MDPs) modelling framework. However, the size of the problem increases exponentially, due to the curse of dimensionality, as more details are incorporated into the model. To this end, neither the ML nor the exact solution approach would provide a satisfactory result within a reasonable computation time. The same is true for heuristic and metaheuristic approaches in converging towards a near-optimal solution, as the size of solution space increases due to the complexity of the problem. Thus, the idea of overcoming the limitation of existing solution approaches through a hybrid solution is very appealing.
This paper aims to extend the preliminary survey in [6] through improved classifications, research methodology, and analysis of the results to classify various published works in the last decade leading to the year 2020 in the attempt of highlighting the trends and focus of VRP in the scope of humanitarian operations in terms of modelling aspects and solution approaches applied. The main contribution of this work is the literature review on recent VRP studies in the scope of selected humanitarian operations focusing on the disaster types, disaster management phases, modelling aspects, and solution approaches. To the best of our knowledge, no work has considered reviewing the applications of VRPs that are based on the suggested humanitarian operations in times of disaster according to the presented classifications. To the best of our knowledge, no such work with the specific humanitarian scope and combination of classifications has been published. From this paper, the challenges in addressing the selected humanitarian operations through specific model characteristics as well as the solution proposed are discussed and the current trend of modelling and solution applied is observed and highlighted. Furthermore, ML is also incorporated as part of the classified solution approach, complementing related past review papers that advocate its application. Finally, the future trends are discussed, providing insights towards research gaps that could be explored.
The remaining of the paper is organised in the following fashion. In Section 2, previous related review papers are discussed. In Section 3, the research methodology is explained. Various aspects of VRP application in the selected humanitarian operations are reviewed according to the proposed classifications presented in Section 4, while Section 5 elaborates the current trends and analysis of the survey findings. The future directions and research gaps are highlighted in Section 6 and, finally, a conclusion is drawn in Section 7.

Previous Review Papers
A number of review or survey papers discussed various VRPs from different scopes and perspectives. However, only a handful of survey papers focused on the scope of the disaster and humanitarian operations. Only three review papers are found addressing VRPs in terms of the application for humanitarian operations excluding our preliminary work [6]. Both [7,8] focused solely on the delivery problem, while [9] looked into the twostage stochastic programming modelling approach of the humanitarian operations. Recent work conducted by [10] focused on casualty management while using OR techniques, which also includes supply and rescue operations, among others, but not limited to only routing problems. Beyond the humanitarian operation's scope, but related to this study, ref. [11] addressed the eco-friendly, unmanned bus routing problem, while [12] proposed and reviewed various RVRP within the field. Ref. [13] discussed the 50 years of progress in VRP since [14], and the same focus is given in [15], presenting a comprehensive study of 277 papers from the year 2009 to 2015, extending the work of [16]. Meanwhile, Ref. [17] discussed the stochastic and dynamic VRPs for various applications. In general, most of the aforementioned review papers (except for [10,15]) only surveyed a moderate number of papers and did not explicitly present the methodology of their review, such as the years in which the reviewed papers are published. Instead, they mainly focused on the findings. Having said that, these review papers all echoed the same suggestions, which include: • the need to address the real-world problems, • the need for the RVRP, and • the importance of stochastic and dynamic elements of the studied problem.
This suggests the persistent gap in RVRP or that the existing work may not meet the criteria of a real-world problem despite the theoretical connotation given. Nonetheless, Ref. [17] highlighted the emerging trend of RVRP, which is strongly related to the stochastic and dynamic elements of routing problem. The advancement of computation power could be one of the key catalysts. Finally, all of the mentioned review papers in this section are presented in Table 1, highlighting the focus of their research as well as the setting of their survey in terms of the year of publication and number of papers reviewed.  [12] Proposed definition for RVRP Real life attributes including stochastic, dynamic, heterogeneity, multi-periodicity Classified based on the constraint classification from [18,19] 1997-2013 50 [13] Decomposition algorithm for exact solution and metaheuristic Attribute learning through naturally inspired algorithm 1959-2009 not specified [17] Addressed Stochastic and Dynamic VRP Comparison with each pure stochastic and dynamic counterpart based on uncertain parameters according to [20] Classified paper based on online, offline computation and hybrid of the two 2001-2015 not specified [15] Surveyed 277 papers based on taxonomy proposed by [16] Shown that metaheuristic is applied the most among the other solutions: exact, heuristic, online solution & simulation based solution Covered multiple attributes of VRPs: split deliveries, multi period VRP, SVRP, DVRP, etc.

Material Collection
This study focuses initially on the review of VRP's application in humanitarian operations with which the search for materials began. From the search, several recurring themes namely the "supply and delivery", "evacuation", and "rescue" operations were noted. Thus, we narrow the scope of VRP's application further into these three specific humanitarian operations. To our knowledge, no study has been conducted on these three specific humanitarian operations put together in the scope of VRP.
Based on the refined scope of the study, reputable review papers are referred to note potential papers that might fit into the scope of this study. Apart from this, an independent search for review materials is also performed through the Google Scholar search engine and the subscribed databases that are guided by the following strategies: 1.
Potential papers are searched with all possible permutations of specific keywords, such as "VRP", "dynamic", "stochastic", "disasters", "search", "rescue", "emergency", "delivery", "evacuation", "vehicle routing", "humanitarian", and "supply". The VRP is set apart from other problems, such as covering tour problem that addresses routing, but not particularly addressing the VRP explicitly. However, known combinations of the OR problems, which include routing problems, such as Allocation Location Routing Problem (ALRP) and Location Routing Problem (LRP), are included within the study. Additionally, this also includes intermodal network problem involving multiple types of vehicle as well as cooperative operations among vehicles, such helicopter and land vehicle.

2.
The search results are filtered by year, in a decreasing order starting with the year 2020 until the year 2010. For each year, the search results are examined page-by-page of minimal 15 search results pages, even though the topic of the searched papers no longer show any relevance to the scope of the study. Otherwise, the search continues in the ongoing pages until no more topic of relevance is obtained in the search. This process is repeated for the year 2020 at three days interval throughout the search duration in order to ensure that any related work is included in the study.
In addition to the papers in [6], the search duration lasted for close to six months (May 2020-October 2020) while reviewing the collected papers in parallel.

3.
All of the selected papers are deemed to be potential papers for review judging by the topic, description in the abstract, few last paragraphs of the introduction section, a glance of the methodology and analysis sections, as well as the conclusions. Furthermore, only English written papers are considered for generating a systematic and consistent study. 4.
From these promising papers, the literature review section or equivalent sections are particularly screened through in order to search for other promising papers. 5.
The resulting papers are then compiled in a Microsoft Excel spreadsheet, where the type and the source of the papers are listed. The source of the papers are then crossed check with the Web of Science (WoS) and Scopus database in order to obtain details of the source, such as the impact factor of a journal and the respective latest quartile performance, with which the quality of the sources is determined. 6.
Materials not belonging to any typical sources, such as Ph.D. theses or technical reports, are not immediately dismissed. Instead, a proper check and study towards these materials are performed in order to judge the relevance, quality, and validity of the work. This ensured that the finding of the review study would cover all good materials from various sources.
Based on the following material collection method, a total of 123 papers are obtained for review. From a considerably huge bulk of potential materials, many are filtered out, as not all materials are explicit in their work in addressing the application. This is especially true for papers addressing VRPs that could be applied to both humanitarian operations and commercial purposes. In this case, unless the humanitarian setting is mentioned in detail, such a paper is dismissed and focus is only paid on those that specifically addressed classified humanitarian operations through the application of VRP in the event of disasters. This is deemed to be appropriate, as the solution and modelling essence of such papers would also incorporate disaster elements that are suitable for the scope of this study.
- : 1 Paper)  0  0  0  2  3  5  6  5  28 15 7  71   TOTAL  0  0  2  3  5  6  12 9  39 28 19 123 From Figure 1, it could be briefly said that the interest in publishing papers addressing the application of VRPs in classified humanitarian operations is gradually increasing from the year 2012 and reaching its peak in the year 2018. The same trend could also be seen in each of the sources of publication type namely the journal and conference proceedings. (Figure 2). When considering all of the sources of publications individually, it could also be observed that each publication contributing roughly to the same amount of papers, ranging from the value of 1 to 4. Thus, no clear source of publication is seen dominating the scope of this review by a large margin. Finally, it could be said that the papers that are published in journals contribute to roughly twice the amount of conference papers, suggesting a healthy amount of comprehensive work published when compared to preliminary studies.

Category Selection
The review is structured and constructed based on the goal of addressing the existing works in a systematically manner, as well as determining the research gaps and trends of works within the scope of the study. Furthermore, in order to integrate the study with existing and future works, existing classifications are noted, compiled, and improved where needed, complimenting previous studies, especially previous review papers that highlight certain gaps within the field. The collected materials are then classified according to the selected classifications that are guided by definitions in the literature.
In this study, the definition of a disaster is as in [8], which is defined as an extraordinary event that could be caused naturally or due to man-made activity. The disaster that could occur with or without warning may also lead to another disaster bringing destruction, despair, loss, and suffering at a scale that could not be addressed by common emergency services.
Based on this definition and guided by the motivation of the study, the materials collected are classified in details, based on the following attributes:

2.
Disaster: types of disasters and phases of DM.

4.
Solution approach: groups of methods to solve the problem.
From these attributes, classifications are then performed in two stages sequentially: 1. Classify papers addressing VRPs based on the applications of three specific humanitarian operations (attribute 1).

2.
Classify the resulting papers under the three umbrellas further in terms of disasters (attribute 2), model characteristics (attribute 3), and solution approaches (attribute 4).
The three specific humanitarian operations that are observed in the literature are: 1. supply and delivery operation, 2. evacuation operation, and 3. rescue operation.
Through the first stage of classification, this study complements the works of [7,8], which reviewed the applications of VRP in disaster, focusing solely on the delivery of the relief supply. The same could be said for [10], which focused on broader OR operations to address the casualty management.
These three applications of VRP in humanitarian operations are further detailed by specifying which disaster phase that each paper addressed in the second stage of classification. These phases are classified into four categories, as proposed by [23] and applied worldwide: (1) Response, (2) Recovery, (3) Preparedness, and (4) Mitigation, under Disaster Phase (DP). Furthermore, disaster type is also included for classification labelled as "Typ" consists of: (i) rapid onset disaster and (ii) slow onset disaster. Through these additional classifications, specific insights in terms of trends could be obtained in relation to the nature of disasters addressed as well as the DM perspectives of the selected humanitarian operations.
In terms of attribute 3, the papers are further classified based on the modelling aspect of complex or rich (R) VRP, deterministic (D) VRP, stochastic (S) VRP, dynamic (Dy) VRP, as well as multi-objective (MO) VRP. Ref. [17] is referred to for the definition of the stochastic and dynamics in the classification. Meanwhile, the definition of RVRP is based on [12]. However, the extended and complex version of typical variants of VRPs, such as CVRP, VRPTW, and VRPPD, are also considered to be rich. The classification of RVRP resonates well with the suggestions from [7][8][9] and ref. [15], which emphasised the importance of RVRP. The DVRP is also included in the classification due to the large number of papers collected addressing such routing problems. Moreover, the classification for DVRP is added in order to investigate the possible trend of research within the context of this study. As is often included in other review papers, the multi-objective attribute of modelling is also of importance and applied in many review papers. As such, the classification for this attribute is also performed.
The final attribute 4 addressed the solution approaches that were adopted by the various papers. Within the scope of this study, a hybrid (Hy) solution is defined as a cross combination in regards to algorithm operators among ML, heuristic (H), metaheuristic (M), and exact solution (E). ML is defined as a solution approach that requires an agent to learn for optimal or sub-optimal decisions or policy, based on supervised, unsupervised, or reinforcement learning. On the other hand, the exact solution approach is defined as a solution that can be solved based on an existing commercial solver, such as CPLEX, Gorubi, and GAMS, which looks into all possible solutions in producing an optimal result as opposed to the sub-optimal results. We note that such a platform also provides heuristic and metaheuristic solution as a choice and took notes on the difference between algorithms, such as branch and bound, and other available solution approaches when making this review. The well-known standard definition of heuristic and metaheuristic is used in the classification ( [24,25], respectively).
The complete abbreviation that is used in the Tables A4-A6 (in Appendix H) is provided in Table 3. Based on these categories, the frequency and focus of the literature available between 2010 and 2020 are presented in Tables A4-A6. Table A4 summarises the characteristics of each paper that researched supply and delivery operations based on the classification that is introduced in Table 3. This is consistent for Tables A5 and A6, for evacuation operations and rescue operations, respectively. Each paper surveyed is then placed under the appropriate categories according to the specified entries in the tables. It should be noted that some fulfill more than one criteria and they are recorded under all relevant categories in this study.

Material Evaluation
The validity of the materials is achieved by continuous checks and critics from the experts in the field studied. This review is extensively and substantially revised and improved based on the previous critical and constructive feedback and review from the experts on the preliminary study [6]. Furthermore, rigorous checks that were assisted by Microsoft Excel allowed for cross-checking, resulting in reduced human error in sorting and placing the relevant papers in their classified groups. Previous review papers are also referred to when constructing the methodology and searching for materials. The materials are collected from reputable journals and conference proceedings, mostly being indexed by WoS and Scopus. However, addressing the small scope as it is, proper and relevant materials are also indiscriminately included. This is to ensure that the resulting trends would cover a wider range of scope instead of trends resulting from a strictly confined scope of sources. Finally, the manuscript is also undergoing various levels of approval processes in terms of its validity, quality, and contribution before its publication.

Review of Literature
This section presents a detailed and comprehensive review of previous published work that is based on the proposed classifications under the umbrella of three specific selected humanitarian operations. In each of the subsections that follows, various works are discussed in terms of modelling, the challenges in modelling the problem, methods in overcoming the problem, as well as the solution proposed along with the challenges that were tackled.

Supply and Delivery in Routing Optimisation Problems
The VRP in humanitarian operations is commonly addressed as the delivery problem in multiple parts of the supply chain, which includes the delivery to the victims, delivery to the emergency shelters and medical centres, as well as delivery to distribution points that are strategically placed near the disaster area.
This implies the crucial role vehicle routing in delivering critical relief aid materials in alleviating human sufferings, as highlighted by [26]. According to [26], the surgical dressings are among the most crucial materials that are needed during a disaster, as it is used to prevent further infection and contamination, as could be seen in Nepal during the 2015 earthquake. This problem is, in fact, more dire, as there will be degradation of hygiene as well as sanitation facilities and care during such an outbreak [27]. Moreover, the lack of transportation, vehicles, or the proper management of these assets during disasters could prevent aid supplies from being efficiently distributed to the victims [28].
In a response to the global need of efficient vehicle routing in humanitarian operations, various publications can be found addressing the problem. In fact, a large portion of the literature that is related to the topic of routing optimisation during disasters focuses on the supply and delivery problem (as can be seen in Table A4). Table A4 presents a summary of the classification of each work in supply and delivery operations. It can be a guide for reading this section.

Machine Learning Methods in Supply and Delivery VRPs
Among all of the specified VRP solution approaches classified, the adoption of the ML solution approach is very limited, due to the curse of dimensionality that is inherited by ML. This is observed in Tables A4-A6. However, some papers did employ the ML solution approach on a realistic problem instance focusing on a specific case study. For instance, ref. [29] incorporated an assessment on sight within a distribution problem by introducing corporation between three humanitarian agents: Relief Assessment Team (RAT), Emergency Response Team (ERT), and Decision Making Agent (DMA) in order to deliver relief supply more efficiently. The DMA learns to make a routing decision for both RAT and ERT through reinforcement learning, which incorporates both the dynamic and stochastic nature of the demand and route access. The initial assessment is made by RAT at a location that is decided by DMA via route also computed by the DMA. RAT, in return, would update the demand and route information back to the DMA, which could be used in order to decide on demand distribution via ERT, which communicates in the same manner as RAT. The highlight of this methodology is the ability of DMA to learn continuously in both online and offline modes. In the offline mode, simulation is performed to generate a scenario in allowing the DMA to learn. Meanwhile, in the online mode, real information is updated and included in the learning process of DMA. Thus, this particular work highlights the strength of the ML solution approach in addressing complex coordination while addressing both dynamic and stochastic elements of the problem in a flexible manner.

Exact Methods in Supply and Delivery VRPs
In [30], a Generalised VRP with flexible fleet size (GVRP-flex) is addressed. Common assumptions are made, such that the vehicles are homogeneous with one depot and several demand nodes that are clustered into a number of clusters. Only one vehicle is allowed per cluster and the trip to demand nodes should only be made once. The problem is solved while using a Column Generation (CG) and the performance is evaluated and compared with two metaheuristics that are inspired by the iterative local search (ILS). Ref. [30] addressed the problem that is inspired by HL, where the road conditions are often poor, due to the effect of disaster that requires the trip to be done by aeroplane. As opposed to [31], which allows for multiple trips, the trip needs to be done once at each disaster location, thus being ideal for the delivery problem with flexible fleet size, although it is somewhat unrealistic when considering the limited amount of transport during the real world application.
The application of exact method solutions, particularly in the supply and delivery using VRPs, are comparable in number when compared to heuristic and metaheuristic solution approaches (Table A4).
Nonetheless, some work, such as [32] (mark with two stars, ), only addressed a single vehicle for the humanitarian delivery operation involving only 20 to 40 affected locations at the maximum. The problem of multi-trip Cumulative Capacitated Single VRP (mt-CCSVRP) is inspired by the response phase of humanitarian operations operation allowing for a single vehicle with limited capacity to perform successive trips to destinations in order to meet the urgency of delivering medical supplies. For a small instance of 20 sites to visit, ref. [32] proposed two Mixed Integer Linear Programming (MILP) model, while, for a larger instance, a resourced constrained shortest path problem model is adopted to minimise the sum of total arrival time. It is evident that the exact method can actually only be applied up to a certain number of instances based on their complexity in terms of constraints and objectives.
A complex representation of a VRP as inter-modal networks is described by [33]. The problem is inspired from the joint effort of relief distribution coming from both local and international participants in the humanitarian operations. This translates into a routing problem for an inter-modal network that includes maritime and land transport. Applying GAMS solver, ref. [33] is also able to address a multi-objective problem, which further increases the complexity of their model. Moreover, they also stressed the importance and potential of inter-modal transport in delivering relief aid supplies during a disaster event by also considering the vulnerability of the road during the disaster, which is based on a case of the aftermath of an earthquake at the Istanbul seaport. Ref. [34], in his work, also considered multi-objective by adopting fuzzy and crisp formulation in solving the routing problem consists of multiple depots that are assigned with dispatch vehicles, distribution points, and emergency shelters. The problem of split delivering supply towards shelters is addressed while using a heterogeneous fleet during the response phase of disaster. However, due to limited resources of vehicles, the model aims to minimise the numbers of vehicle, while also minimising the total arrival time of vehicle to a depot or to the emergency shelter. Unlike the common VRP, which deals with network involving depots and shelters, this study dealt with three layers of network consisting of depot, distribution centre, as well as emergency shelters, thereby making sure that the distribution centre is able to accommodate deliveries to emergency shelters with vehicles only allowed to go to any one shelters and remain stationary. Thus, the TH method is applied to solve for the objectives introduced in the model.
Similarly, the multi-objective problem is also addressed by [35]. In this routing problem, the efficiency of maximising the delivered supplies through a short time window is considered. Furthermore, the fairness of delivery is also addressed among emergency points. Apart from having multiple depots, the model, known as dynamic continuous time flow network, tackled the complex problem through the dynamic constraint, which depends on time incorporating changes of demand, donation, as well as the condition of the road. In this problem, the decisions involved are the dispatching decision as well as the routing decision. Furthermore, a continuous-time dynamic network flow algorithm that is based on the augmentation of the continuous-time dynamic residual network (CTDRN) involving four steps is proposed in order solve the problem.
Ref. [36] seeked to minimise the cost of travel, in the event of arc failure. The arc failure could be reflected in many scenarios, especially in times of disaster. This study also considered the condition where a route failure is unexpected, where the on board routing software is not available, as well as when communication with a control centre is difficult, due to the damage on communication facilities. It allows for a driver to pick a primary road along with its alternative should an arc failure occur. To solve the problem, ref. [36] proposed a modified label-correcting algorithm where the primary path and alternate secondary path are computed in two parts, respectively. Additionally, through the dominance rule proposed, the huge number of possible instances representing failure in each arc are effectively reduced.
A more detailed model of delivering supplies in a disaster event that is addressed by the exact solution approach is shown in [37]. The Mixed Integer Programming (MIP) model highlights an interesting concept of sharing relief supply among emergency shelters. This concept works by allowing for sharing according to the assigned priorities as well as the number of evacuees. Ref. [37] addressed the problem criteria based on the Japan government's white paper on the Great East Japan Earthquake in 2011. The objective of the model is to minimise the total transportation cost as well as the penalty that is imposed for every failure in satisfying the demand. The goods' priority are represented by each weight that is associated with the good types. Unlike the typical VRP model, shelters that have homogeneous delivery cars with fixed capacity take the role of a depot, while also having their own demand. As such, the sharing concept among shelters could be modelled. Furthermore, each shelter is also associated with priority according to urgency. Here, real data of Chofu City in Japan are applied to generate real data for the parameter of the problem. The CPLEX solver is utilised for solving the problem.
The fairness of relief distribution for event such as earthquake is studied in [38]. The fairness in this work is considered based on the arrival time of vehicle. The satisfaction level of the victims is considered as a utility function with the objective of minimising the range of delivery times from he transportation vehicle. The proposed model is then compared with model for different objectives, namely one that minimises the arrival time as well as the one that minimises the transportation cost. The 2011 Japan earthquake is taken as a case study focusing on the city of Fukushima and Noda village. The results shows that the satisfaction level of victims improved, even with a slight improvement in the time range of delivery. Furthermore, the proposed model improved the fairness of delivery at the price of a higher cost of operation as compared to the two models compared.
Ref. [39] focused on the integration of VRP in the Decision Support System (DSS). However, the embedded VRP considers the extension of conventional VRP. such as the Heterogeneous VRP (HVRP), SDVRP, and the Multi Travel VRP (MTVRP). The prototype model in DSS also considered partitioning the deliveries, where the vehicle has the choice of not delivering all aid within its capacity. Two objectives are considered in this problem: the total time travel per route is minimised for all vehicles as well as to minimise the arrival time to each demand nodes. These two objectives are solved based on priorities, as defined by the user. A small test instance based on hypothetical scenario was applied to validate the model where the problem is solved optimally.
Meanwhile, ref. [40] considered the uncertainty of travel time and demand in CVRP with the robust optimisation approach. The objective is to minimise the total travel time as well as the unmet demand. In this model, two depots are iconsiderednstead of just one, whereas the rest of the formulation is similar to formal CVRP formulation. Finally, based on the robust counterpart obtained, the model is then converted into a deterministic MILP model, which is verified through the commercial solver Gurobi. As expected, the solver struggled for larger instance of the problem, which, according to [40], is largely due to the increasing sub-tour elimination constraints.
An interesting approach of dynamically updating route condition is shown in [41], where a dynamic routing is executed during a disaster. In this problem, the realistic assumption that roads may be affected by disaster with no communication available is addressed. Two problems are specifically addressed in this problem: to obtain an optimal route in delivering relief aid, and to be informed on the route condition while performing the humanitarian operation. Two algorithms are proposed to be solved for the problem, namely the Information Gathering Algorithm (IGA) and Route Construction Algorithm (RCA). This process of information gathering is executed through the IGA. Meanwhile, based on the information that is gathered by IGA, RCA is then dynamically constructed at every query interval while using the Dijkstra shortest path algorithm. The objective is minimise the cost incurred for all vehicles when travelling. A simulation test is performed in order to validate the model and solution based on the rectangular area pivoted from the Nakagyo fire department in Kyoto. The results are compared with the existing method of Geocast, where an improvement can be observed in terms of gathering real time information and the arrival time.
Ref. [42], on the other hand, addressed the stochastic problem of LRP under a fuzzy environment. In this problem a multi-commodity, multi-depot, with heterogeneous vehicle with different capacity and speed are considered. A split delivery is allowed, illustrating a more realistic version of the problem. In this problem, three objectives are addressed. The first objective is to minimise the total expenditure cost, including opening distribution centre, transportation cost, travel expanses, as well as the transportation cost, from suppliers and the distribution centre. The maximum travel time is minimised as the second objective, while the third addresses the unmet demand that should be minimised. Moreover, in this stochastic problem, the available traffic, transportation cost, volume of relief aid, fix cost of providing distribution centre, the distribution points, as well as the normal speed of the vehicle is uncertain. A MIP model is converted to its deterministic counterpart via a crisp model to further address the fuzzy parameters, which are the capacity of the distribution centre as well as the transportation cost. Finally, an exact solution approach is proposed that is based on the proposed crisp model as well as the method presented by [43], which is solved by the GAMS software while using the Baron Solver.
The LRP with uncertain demand and travel time for distributing delivery aid in the event of a disaster is also addressed in [44]. With the objective of minimising the total arrival time of vehicles, the problem is formulated with the associated LRP constraints. This proposed model is applied to small experimental cases, where the commercial solver of LINGO is implemented in order to solve the problem.
Ref. [45] is one of the few that studied the problem of preparing emergency supplies for a disaster instead of responding towards a disaster through the utilisation of a neighbourhood disaster station. The idea is to place the neighbourhood disaster station at a strategic location while restocking the station with emergency supplies that can be distributed to the affected neighbouring areas in the time of a disaster. However, it is unknown when a disaster would strike; thus, to prevent the supplies from expiring, they are to be sold to maximise returns. Once they are sold, the station would need to be replenished. Their objective is to maximise the profit by determining the optimal replenishment policy for each supply items as well as to provide an optimal route for the vehicles that will minimise the cost. A heterogeneous vehicle is considered and the replenishment of the stations is done according to their capacity from a single depot with time constraints. Here, the MIP formulation of the problem is solved by CPLEX solver for a small instance. Because all of the information is known, the problems are denoted as a deterministic problem.
Ref. [46] addressed the scheduling and open vehicle routing of pick up and drop off aid delivery problem through the perspective of a humanitarian coordinator. The coordinator as the decision maker is to decide how many vehicles should be rented from the third party logistics service provider while taking into consideration the limited duration of the rent in the aftermath of a hurricane. This includes tasking each vehicle to a set of clients in an optimised manner and routing the vehicles optimally, such that the total transportation cost is minimised. The problem is considered to be an Open Vehicle Routing Problem with Time Window (OVRPTW). Two models are presented for the problem. The MIP model is formulated by defining an artificial source and sink, and it is solved via a commercial solver. The MIP model is further reformulated into a path based integer model. Here, the routing of vehicle along with the rent period constraint are addressed by the column generation approach. The obtained results indicate that the column generation approach outperformed the commercial solver, especially when the instance size grew.
A DSS is proposed by [47] for aiding the distribution during a disaster event where an embedded VRP module is discussed, resulting in a model driven based decision support system. Here, a deterministic model addressing a heterogeneous fleet that performed split delivery as well as multi-trip is presented. In order to validate the model, a case study was done based on the earthquake and tsunami in the coast of central Chile on 2010, where the smaller instance version is solved by the GAMS solver as well as the prototype solver of the decision support system.
Ref. [48] combined the commodities distribution in the LRP in order to ensure efficient aid distribution in the post disaster phase. However, the location of distribution centres as well as commodities distribution is done during the pre-disaster phase. The problem that focuses on time deprivation in emergency logistics model further reflects the critical scenario, where there is a shortage of supplies. The time deprivation serves as the loss cost incurred during the humanitarian operation. The optimisations are done where the commodities distribution and site selection are considered at pre-disaster phase. The first objective is then to minimise the all logistics cost, which includes the operating cost of the relief centre, the inventory cost, the transportation cost pre-and post-disaster, along with the deprivation cost that is incurred by the unmet demand. The second objective on the other hand is to minimise the total transport time for any edges used. The problem is initially modelled as two stage nonlinear bi-objective integer programming and was further simplified into a single objective integer programming model. Based on a LINGO 10 commercial solver, small randomly generated instances of the model ares applied to validate the model.
In [49], a more realistic problem is considered that is based on the observations done on the Nepal 2015 earthquake. According to the [49], the hilly landscape of Nepal often reduced the road network to a tree network due to a disaster. This often affected the second stage of humanitarian operations especially in delivering medical supplies where mountainous road are often disrupted as the course of the disaster is prolonged. Furthermore, it is more realistic to allow for split delivery representation, as this is, in fact, the case in Nepal earthquake 2015. Based on this observation, a Capacitated Vehicle Routing Problem on Trees with Split Deliveries (CVRPT-SD) is proposed to address the issues of secondary stage in humanitarian operation. A modified MIP model and a heuristic solution based on decomposition scheme are proposed. The primary objective is to minimise the unsatisfied demand, while the secondary objective is to minimise the travel cost. As a solution strategy, the model is decomposed into three models: a delivery model, a capacitated primary VRP model, and a capacitated secondary VRP model. The delivery and capacitated secondary VRP network each have their own designed algorithm, while the capacitated primary VRP network could be solved while using CPLEX efficiently. Both the toy and real dataset that are based on Sindhupalchow district are applied to validate the model and solution. The computation using the real dataset that is based on a general model proved to be prohibitive using CPLEX. The decomposed models are solved instead with their respective proposed exact solution approach.
More realistic depictions of relief distribution during disaster events are shown in [50], where a multi-conflicting objective problem with multi-trip, multi-depot allowing for split delivery among vehicle is addressed. The issue of limited vehicles during a disaster is stressed upon, justifying the multiple trips and split delivery among the limited vehicles, as in [47]. Furthermore, dynamic updates are allowed and solved through the rolling horizon mechanism. Meanwhile, the conflicting objectives are solved through the Pareto optimality method, where the Pareto optimal set is obtained through the augmented constraint method. In this problem, three objectives are considered, namely to minimise the amount of unsatisfied demand based on priority of relief aid, to minimise the total travel time, as well as to minimise the sum of absolute deviation of the unsatisfied demand with the latter aim to ensure fairness in relief distribution. The multi-objective MILP model is solved while using GAMS/CPLEX per each rolling horizon stage. The solution and model is finally validated while using real world data from Tehran for an earthquake case.

Heuristics and Local Search in Supply and Delivery VRPs
As an alternative to the exact method approaches, solving VRPs through a heuristic and metaheuristic approach, could be considered to be a typical choice. This includes the problem of supply and delivery for HL where various heuristic approaches are proposed. Ref. [51], for example, applied heuristic to solve for the large instances of their Last Mile Distribution Problem (LMDP). This problem looks into the effect of considering multiple and conflicting objectives, which are the efficacy, efficiency, and equity, which are appropriate in assessing effective distribution during a chaotic event of disasters. The efficacy is regarded as the amount distributed per node during the time spent, while the efficiency is measured based on the time of delivery and finally the equity is measured based on the disparity among nodes in terms of efficacy. In order to investigate the effect of applying these objectives on the resulting routes, the problem is restricted to one trip operation per vehicle while allowing split delivery among vehicles. Because of this limitation, no vehicle carries capacity enough to fully served one demand point while the total demand is sufficiently supported by the depot. Two decision variables are proposed to govern the routing as well as the allocation of supply, respectively. The effect on route structure is observed based on preferred objectives among three different combinations of three objectives, computed together, which derived the model further into three models. The problem that is formulated with the set partitioning formulation is solved by a commercial solver for small instances and heuristic based on the framework of Greedy Random Adaptive Search Procedure (GRASP) for large instances. The first two models are solved by heuristics that are based on the constructive phase and improvement phase at each iteration for a set number of repetitions. The constructive heuristic is based on an insertion heuristic of randomly restricted selected node candidates. This solution is then improved by performing a local search. This solution approach, however, could not be applied to the model that focuses more on the equity issue of distribution. For this model, a two-phase heuristic approach is utilised. These two phases comprise of setting partial allocation to each node at the initial phase, and the second phase considers adding a node to the route or adding the supply in the initial phase to serve the remainder of the demand. The results concluded with several practical insights, such as that all three objectives could be achieved better by utilising a large number of small vehicles as opposed to commercial practice, which prefers large carriers. This is due to the flexibility and mobility of small vehicles in accessing remote destinations as well as the ability to traverse the affected transportation network. This also improved the equity, as more demands could be covered and served more fairly. However, this comes with non-excessive addition of travel and operation costs that are considered to be reasonable in times of urgency.
Ref. [52] proposed a DSS incorporating both the optimisation platform and simulation, making use of available graphical database. It addressed the coordination problem between private and relief organisation during the response phase of a disaster. In this problem, the lead time between request and delivery is minimised, depending on the road availability. A MIP model is proposed for this problem. The transfer points, demand points, supply points, as well as depots are depicted as nodes or vertices. To solve the model, a Star algorithm that is piloted by a heuristic method is applied, resulting in a near-optimal route. Meanwhile, the facility location problem is addressed by TS. By solving the location problem, the average service time (from request to delivery point) is minimised through relief supply coordination.
A large scale disaster response by coordinating an inter-modal network of both helicopters and vehicles is studied by [53]. Because of the limited number of helicopters, supplying directly to Medical Aid Points (MAPs) is considered impossible leading to the cooperation between helicopters and vehicles in delivering the emergency supplies. Referring to their previous work, refs. [53,54] further stated that the classical Fuzzy C Means (FCM), as applied in their earlier work, tends to ignore the balance of capacity, which could lead to huge wastage of medical supplies. Motivated by these problems of employing FCM, ref. [53] proposed two improved FCM by introducing FCM with Capacity Constraint (FCMwCC) as well as FCM with Number Constraint (FCMwNC). The FCMwCC works with additional proposed heuristics in assigning MAPs to each EDC, since the newly introduced capacity constraint prevents the usual iteration of FCM. FCMwNC, on the other hand, assigns a specific value to the numbers of MAPs assigned to each EDC preventing the unbalance assignment, as in their earlier work. Simulated instances are generated in order to validate the solution and model showing balance distribution, which leads to improved fairness as well as efficient distribution.
The HL problem is presented in [55], where the total travel cost is minimised by synchronising the on site visit with delivery of medical supplies. The real case of the Democratic Republic of Congo (DR Congo) is taken to present the two problems that are faced by the hospital administration in maintaining the pharmacy in a rural hospital as well as the high cost of dropping and picking up a supervisor for an on site visit to rural hospital. The problem is addressed by CVRP with synchronised Pick-ups and Drop-offs (CVRP-syncPD). A heuristic method is developed that is based on a cluster first and route later technique with the ability to adapt the synchronisation between pick up and drop off for an aircraft carrier and on site supervisor. The hospitals are clustered according to the capacity constraint, creating a set of medical routes. The synchronised route is later chosen from this set of medical routes by comparing the synchronisation requirement for on site supervision. These synchronised routes are then evaluated in order to minimise the total cost incurred. The test instance for validating the model and solution is based on the Bandundu province, where 41 hospitals are located.
The Inventory Slack Routing Problem (ISRP) during a disaster where the inventory of relief item is kept under tight observation for fear of shortage of aid is considered in [56]. The problem is solved by decomposing it into two sub-problems to deal with the routing problem and the inventory problem respectively. From the routing side, the Hierarchy TSP concept is applied where the delivery of one vehicle to a cluster of dispensing sites should be done based on priority on the urgency of the demands. By decomposing the problem, the objective of maximising the minimum slack is now converted to a new objective, which is to minimise the delivery time weighted by priority of the dispensing site. The routing is then solved through the hierarchy routing heuristic or sweep algorithm. Consequently, the allocation problem is then solved directly by solving the resulting converted LP model of the allocation problem. From the computational results, the sweep algorithm outperformed the hierarchy routing heuristic, judging from the resulting total travel time as well as the minimal slack.
Ref. [57] focused on the fair allocation of perishable food to ease the problem of hunger in a humanitarian crisis by addressing the problem through the perspective of the ALRP. The problem is based on PDVRP, where waste is avoided while ensuring fair allocation to welfare organisation via optimal route. The envy free allocation could be obtained through the known min-max fair share algorithm. However, it is shown in their work that this algorithm may not be able to satisfy the operational constraint, such that the total cost is minimised. Instead, the fair division problem is assumed for the allocation problem, while the routing component of the problem is formulated as a MILP model. The bi-objective of the model is to minimise the maximum deviation of envy free allocation and the total travel cost respectively. In order to obtain an optimal solution, the problem is decomposed while using Bender's decomposition method, where the relaxed version of the original problem is obtained. However, for larger instances, a hybrid heuristic of local search and a greedy algorithm is proposed. This proposed method consists of pickup node clustering, assignment of delivery nodes, tour construction, and, finally, the improvement of the solution. A numerical experiment is based on the real case, showing the capability of the proposed heuristic in obtaining a balanced solution between quality and computation efforts.
Ref. [58], in his work, argued the necessity of addressing the HL from the perspective of minimising the arrival time, rather than the cost minimising problem to reflect the real urgency of delivering critical supply to the victims. The problem is modelled as the selective cumulative multi-VRP with stochastic travel time, where the priority of the affected area is considered based on the destruction that is caused by the disaster as well as the population of the affected area. Because the fleet of vehicles is limited, only selective areas are served from a single depot. Unlike typical literature addressing the humanitarian problem, homogeneous vehicles are considered to have unlimited capacity. The objective is to maximise the utility that represents the decreasing function of the arrival time. To cater for the urgency of planning, an iterated greedy heuristic is proposed that consists of construction and destruction stages. Benchmark instances are used and the results shows good quality solutions that are computed within a reasonable and comparable time.
The strategy to decompose VRP into multiple TSPs are proposed in [59] through the approach of cluster first and route second. The problem is formulated as a MIP with the objective to minimise the total transportation cost, the pipeline inventory cost, as well as holding cost with back order for the scenario of hurricane Irma in Florida, where the gas supplies for gas stations were urgently needed due to people fleeing from the storm. Prior to solving the routing decision, the optimal replenished amount of gas stations and inventory level or route served by the vehicle are first computed. The solution of the continuous approximation serves as a guidance for the clustering, which would decomposed the VRP into multiple TSPs. Two clustering methods are employed and compared, namely the local observation based clustering and the K-means clustering algorithm to determine the routes. Once the clustering is performed, the TSP is then solved using Christofide's method.
Ref. [60], on the other hand, used the ability of UAVs in order to aid land transport in transporting relief aid during flood. The min cost flow model is applied with a single source and a single destination to routing model, which would also account for the capacity of the arc. Furthermore, the Hortonion Overland flood model is applied to model the problem in a more realistic fashion. Two heuristics are employed to solve the model, depending on the phase of routing. Initially, the quickest flow algorithm is employ to construct the initial routing. However, when flood data affecting transport network are updated, the vehicle is forced to reroute to a new route computed by the earliest arrival flow algorithm. At this point, the model is converted into a multiple source single destination problem. Both the virtual network and real network of Bangalore were applied in order to validate the model.
An Emergency Supply using the Heterogeneous Fleet (ESHF) problem is addressed by [61], where a multiple supply point with limited inventories is used to supply relief aid to multiple delivery points through a multi-trip operation. The problem is similar to the VRPPD with a split delivery, but with an improved representation of the problem addressed. The total operational time that includes the travel time, loading time, as well as the unloading time is reduced as the objective of the model, and the problem is modelled through a MIP formulation. All of the proposed constraints are based on the routing constraints, capacity constraints, timing constraints, and other constraints. The heuristic approach is applied to solve the problem for realistic instance. A case study that is based on the forest fire in the Province of Teruel, Spain was used, through which an efficient plan is obtained when applying the proposed heuristic.
More recently, ref. [62], in their work, addressed a drought scenario. A Multi-Depot VRP (MDVRP) is considered to reflect multiple water sources where water is to be delivered to nodes with demands. The objective is to minimise the total travel time, given the unloading time of the water supply at each demand point that is assumed to be the same. In order to solve the problem, ref. [62] first allocated each vehicle a water source with a set of their unique demand points. Subsequently, the Clarke and Wright algorithm with 2-Opt moves is applied to solve the CVRP with single depot for each vehicle. A case study on the Brazilian scenario showed potential for improvement for practical planning in terms of management and financial.
Ref. [63] planned for helicopter routing, as opposed to ground transport, due to the disrupted road and bridges during a disaster event. Five stages of planning for relief distribution and transportation of personnel in operation are specified where the first stage is to reassess the transport network, followed by the estimation of requirement based on the demands. Meanwhile the air transport is reassessed based on capacity, stability and limitation. Subsequently, the number of helicopters is estimated for the operation. The first operation was a combined distribution for demands which exceed the capacities of the helicopters. The subsequent operation is then optimised based on CVRP with time window where the total time of operation is minimised. The problem is solved using the Clark and Wright algorithm with 2-Opt improvement. The proposed model is applied to the 2011 flood in Serrana Mountain Region, Rio de Janeiro, Brazil for validation.
In [64], a 2-echelon VRP is proposed, where a two-delivery level via drones is addressed during the time of a disaster as a PDVRP. Delivery is made to various first aid locations by the same drone. The model that is based on the mathematical formulation from [65] is solved through a fast constructive heuristic based on a biased randomised algorithm. Interestingly, parallel computing is utilised to fulfil the requirement of fast delivery, which made agile optimisation possible. The node addition process in the constructive phase is improved through biased randomisation, instead of depending on a pure greedy mechanism and local search is applied in order to improve the quality of the solution. Instances from [65] are applied to validate the proposed model, where the commercial problem instances are converted into HL operation instances. Increased performance is observed through parallel computation.

Metaheuristics in Supply and Delivery VRPs
Looking into the multi-objective VRP in this problem domain, ref. [66] proposed a delivery model that takes into account the congestion level that is based on the destruction of the road network and a chaotic environment involving the road network. The road capacity is assumed to be affected by the disaster and it can be associated with the congestion level that is related to the flow capacity of the road. Both the total delivery time as well as the total cost incurred using the weight coefficient transformation method are optimised. In addition, more priority is given according to the weight assigned. Genetic algorithm (GA) is applied in order to solve the capacitative vehicle routing problem with time window.
Ref. [67], on the other hand, proposed a model that minimised the total cost of delivery and transportation, while, at the same time, maximising the satisfaction degree of the recipient based on the anticipation, endurance, and failure level. The problem is stochastic as the actual road capacity is uncertain, modelled through an uniform distribution as compared to its designed fixed capacity. The satisfaction level is expressed in a trapezoidal fuzzy function that is based on the anticipation degree, endurance degree, as well as the failure degree. The travel time is computed based on the Bureau of Public Road impedance function. A Simulated Annealing (SA) is employed to solve the problem.
Ref. [68] is of the opinion that selective VRP is more appropriate in representing the routing problem with uncertainties as opposed to conventional VRPs. They argued that hard choices are sometimes unavoidable and need to be made, regardless of the extenuating circumstances. This translates into the fact that not all of the nodes could be visited, thus the Selective VRPs. They proposed three selective VRP formulation based on fuzzy Stochastic Vehicle Routing Problem (SVRP), reliable SVRP, as well as robust SVRP. These proposed models fill the gap, where, previously, only deterministic Selective VRP are addressed. In order to solve these models, ref. [68] proposed three improved Parallel GA. The three algorithms differ from each other based on their communication and diversification strategies. Their performances, together with conventional GA, are then compared by testing them on both small and large network. Among the proposed model formulation, Fuzzy SVRP shows better potential in addressing large problems, due to the fuzzy nature of various uncertain elements in the operation.
A common approach of clustering first and routing later can be seen in [69], which addressed the problem of routing vehicle to large affected destination due to disaster while having a large amount of concentrated medical relief supply at distribution centre. In this problem, a single depot and a fixed number of homogeneous vehicles are assumed with multiple affected areas for delivery distribution. Furthermore, in order to address the large quantities of relief aid, several parameters are proposed based on the quality, handling, and packaging of the distributed materials, which are incorporated into the constraint of the model. The objective is to minimise the total travel time of the delivery operation. In order to solve the problem efficiently, a K-Means clustering algorithm is applied to cluster various emergency location into groups before computing the route for each cluster while using Particle Swarm Optimisation (PSO).
Ref. [70] viewed the critical delivery of emergency supplies from the supply and delivery perspective by considering the supply that is needed for each centre in order to distribute to selected emergency points. However, in disaster time, such as an earthquake, the demand is usually uncertain. This leads to penalty added should the supplies be insufficient for delivery. Strategically, the problem is divided into two, where clustering is needed to assign emergency points to each distribution centre. The objective is to minimise the total cost, which includes the penalty for delivering insufficient demand, penalty for servicing more than demand needed, as well as transportation and operation costs. The first problem is tackled by converting the multi-depot VRP into multiple single VRP through the nearest assignment and average distance method. Meanwhile, the second problem is solved through their proposed improved version of GA. The tests are then conducted based on two numerical examples where the first consist of a simulated instance and the latter is based on the water-logging disaster in Changsha City, Hunan.
Meanwhile, an interesting approach of the ACO can be found in [71], where the problem of vehicle allocation and dynamic demand are considered. The problem placed the government as the decision maker that needs to decide on routing solution, depending on loading capacity of the vehicle, the time taken for vehicles to gather at a distribution centre, as well as the time that is taken for a vehicle to reach disaster areas. These three variables changes over time turning this problem into a dynamic VRP with changing demands and numbers of vehicles based on information updates. The problem is formulated as MDPs, where the objective is to minimise the total longest transportation of all vehicle per each dynamic phase. To allow for the continuous update of information, dynamic solution approach is needed in order to compute the route allocation for each vehicle online via the dynamic programming approach. This new update is embedded in the route allocation algorithm, where the assigned route is computed while using the ACO algorithm. A case study based on the earthquake event in Lushan County in 2013 is adopted in order to validate the model and solution.
A SA is applied in [72] in order to address a multi-objective blood donation problem. The problem is in two-fold: blood mobile location problem and the shuttle routing problem. These problems addressed the issue of perishable platelets, uncertain blood donor's arrival, as well as supply chain planning, which involves both facility location problem and VRP, whereby maximising the blood collection while minimising the cost is difficult. The first problem of locating blood-collection mobile is solved by the crisp multi-objective linear model, obtained through the conversion of the respective multi-objective fuzzy programming model. The output of the blood donation mobile location problem is then applied to the VRP (second model) that is solved by the SA in an attempt to maximise blood donation while minimising the cost that is incurred by both blood mobile and blood donation shuttle. Several test instances are utilised to validate the model and solution method. The first model is solved by CPLEX, while the second model is solved by the proposed SA. The results shows applicability of the proposed solution.
The problem of relief distribution under disruption risk through the LRP is addressed by [73]. Particular to this problem, the risk of route disruption, the number of vehicles, as well as temporary distribution centre's capacity is considered in the first 72 h of emergency response while trying to maintain fairness in delivering relief aid. Multi-commodity relief items are considered to be delivered from Temporary Distribution Centers (TDCs) to Victims Temporary Residence (VTRs) via heterogeneous vehicles on the ground. In order to account for the disruption, each respective disruption from the route, vehicle disruption, as well as the TDC capacities disruption would incur an additional cost, which would need to be minimised as part of the multi-objective problem. The operational cost which includes the setup of TDC cost, the transportation cost, the violated time window penalty, as well as unmet demand penalty are also minimised. Furthermore, the risk of disruptions is measured via the measured of Conditional Value at Risk (CVaR), which is easy to incorporate in a LP model for optimisation purposes as compared to the measure of Value at Risk (VaR). Adapting the CVaR into the objective function, the model is then linearised to address the nonlinearity of the initial model. The problem is then solved by GA based on the new proposed representation of the chromosome specific to the problem. Better results are obtained when compared to the optimal solution that was obtained by CPLEX in terms of computation time.
An integrated LRP is presented in [74] in order to govern the UAVs application in relief distribution. The optimal location of a warehouse is determined, followed by the routing of UAVs with the objective of minimising the humanitarian operations' cost as well as the number of hubs that are utilised by the UAVs. Interestingly, in the model is the introduction of a penalty for any loss of life due to the delay in delivering relief aid. Besides that, the operational costs for setting up a depot as well as the operational cost of UAVs are accounted for in the objective. The small instance of the model is solved by CPLEX using the branch and bound algorithm. However, to solve for the real instance, the bi-layered Large Neighbourhood Search (LNS) is employed with the inclusion of Informed Rapidly explored Random Tree (IRRT), such that the path is re-planned and improved. The existing benchmark instances are based on the damages and relief demand in the central region of Nanto.
Meanwhile, a multi-conflicting optimisation criteria regarding the last mile distribution of relief aid is addressed in [75], where the issues of safety, efficiency, equity, and effectiveness are considered to allow for more a realistic representation of the problem. Specifically, the problem search for an optimal route to deliver relief aid from a distribution centre to a set of demand nodes. Here, the multi-modal transport, multi-depot, multi-loading, as well as split delivery are considered. The problem is formulated as a compromised programming for the last mile distribution, where non-dominating criteria are achieved as a solution with the optimal trade-off between all the criteria. The utilisation of metaheuristic is warranted due to the complexity of the problem. Evolutionary algorithm is avoided, since slight changes in potential solution might render the solution infeasible. Ref. [76], on the other hand, focused on the multi-objective, multi-period, multicommodity open LRP where special care is given to the condition and the short term repair of the affected road while distributing supplies and repair crew during the post disaster phase. In this problem, all the cost associated with setting up distribution centre, travel cost as well as repair cost is minimised. Furthermore, the maximum travel time of the route is also minimised while maximising the quality of the route by ensuring that the service of repair equipment is done along the way. What is unique in this problem is the fact that the safety of crewman is considered within the problem model, which is rarely the case in literature. Because of the multiple objectives needed to be solved by the problem, a multi-objective metaheuristic that could achieve Pareto optimal solution is desired. A Non-dominated Sorting Genetic Algorithm-II (NSGA-II) is utilised, along with the multi-objective PSO (MOPSO), to solve the problem.
The LRP with Simultaneous Pickup and Delivery (LRPSPD) under uncertainty of victims to be evacuated during disaster is solved by [77]. The objective is to reduce the total cost of operations, which includes the transportation cost, the setup of an emergency logistics centre cost, and fixed cost for vehicles. The problem is also treated as a multi-depot problem, as is usually the case with LRP. A numerical experiment is conducted by varying the Decision maker Preference Index (DPI) and the resulting objective value is observed based on the risk that the decision maker is willing to make. The problem is formulated while using the chance constraint programming approach and finally solved by GA.
The sustainability issue concerning the social, economic, and environment concerns when relief aid delivery operation is performed in times of disaster is addressed by [78]. The sustainability issue is regarded in terms of the access, equity, and needs of fulfillment. What links these factors in the proposed VRP is the amount of relief, as well as the time of arrival, which is reflected by the Victim's Perception of Satisfaction (VPS) with regards to sustainability. As such, the objective is to maximise the minimum VPS as well as to minimise the largest difference of VPS among victims at an interval. The problem is formulated as a Mixed Integer Non Linear Programming (MINLP) model, where split delivery is allowed in delivering multi-commodities via heterogeneous vehicles from multi-Relief Distribution Center to multiple Affected Specific Areas. For this problem, a GA is adopted as the solution approach, due to its balanced trade-off between solution quality and computation time, which is shown in the test case of the Wenchuan Earthquake.
Meanwhile, ref. [79] addressed the single layer LRP for relief commodity distribution, where priority is given to commodity and demand while considering the road damage. The priority to relief commodity is given through a set of real numbers that differentiate them from urgent relief material and general life relief material. The demand priority is reflected by the relationship between the urgency of demand and their time tolerance. A satisfaction function is introduced through a soft time window of each node and the time sensitivity coefficient where the link to urgency of demand is highlighted. Additionally, weight is assigned to each road signifying the level of damage sustained by the road. The two objective functions are to maximise the satisfaction level and the penalty incurred by violating time window is minimised. The formulated model is solved through the NSGA-II where the first objective function is simplified to ease the analysis. Two layers of coding are introduced with regards to the location decision as well as the routing decision. Finally, Solomon's instance is applied in order to validate and evaluate the model and solution.
Ref. [80] focused on the location inventory routing problem (LIRP) to prepare for both pre and post disaster by looking into three stages of the supply chain. This problem focuses on the earthquake in Tehran as a case study. The model is divided into two addressing the pre-and post-DM. In the first stage, the ideal location of warehouses is considered noting the hospital as a relief points as well as potential affected areas due to disaster. The cost of setting up warehouses, the cost of transporting the supplies from supplier to warehouses, as well as the inventory cost are addressed. Meanwhile, in the second stage, the distribution of the relief aid is considered based on the Multi-Echelon, Multi-Depot VRP (MEMDVRP). The first echelon reflects transporting the supply between warehouses and hospitals. The second echelon focuses on the transport and routing between hospitals and the affected areas. Three objective functions have been considered for the problem, namely to minimise: (i) the operational cost caption in the pre-disaster operation, (ii) the cost that is incurred by the post-disaster humanitarian operation, and (iii) the total relief time during the relief operation in the post-disaster phase. The problem is formulated as a MILP model, where the small instance is solved by the constraint method. The larger scale of the problem is solved by the NSGA-II as well as the Reference Point Based NSGA-II (RPBNSGA-II), where they outperformed each other for the case of smaller instance and larger instance, respectively.
The LIRP to supply multi-commodity to disaster sites is addressed in [81] while using two proposed models. The former model is based on the the assumption where one affected area only needs to be served by one distribution centre, while the latter deals the problem with the flexibility of split delivery allowing for more than one distribution centre to serve one affected site. In this problem, a two level planning approach is considered. First, the location of distribution centre is sought for with different capacity to be considered. The second level focused on distribution of the commodity through an optimised route. A multi-objective is defined for the problem that is based on this level of operations. The first level sought to minimise the cost of setting up warehouses and distribution centre, as well as the cost of storing commodity. While the second phase aims to minimise the travel cost incurred by routing the vehicle for relief distribution, the travel time as well as maximising the route reliability between crucial nodes. As researched in [76], the model is solved by the NSGA-II and MOPSO. Both of the solutions were compared through 35 instances, whereby the NSGA-II proved to be better when addressing model with uncertain parameter. Meanwhile, MOPSO is more efficient when solving the deterministic model.
The problem of scheduling the Disaster Medical Assistance Team (DMAT) and the optimal relief distribution simultaneously is considered in [82]. The intention is to coordinate both operations so as to allocate limited resources in the most efficient way in times of disaster. The objective is to minimise the total service completion time. The problem is further structurally analyses to decompose it into two types of problem formulation: P1 and P2. Looking at the DMAT routing alone without the relief aid allocation and distribution, the initial problem would be reduced to Time Dependant VRP (TDVRP) (P1) with the objective to reduce the total completion time. P1 could be solved to near optimality through the Rolling Horizon (RH) algorithm. The problem is then reduced to a special case of the original problem (P2), where the objective is to minimise the total tardiness of delivering relief aid. Based on this analysis, a solution approach is proposed, where P1 is solved by the modified Artificial Bee Colony (ABC) algorithm. Next, by applying the DMAT route to the initial solution, P2 is then solved by the RH algorithm. This work concludes with a numerical experiment that is based on the 2016 Kumamoto Earthquake in Japan, where the solution was obtained in practical time.
Ref. [83], on the other hand, addressed both the CVRP as well as the Split Delivery VRP (SDVRP), while trying to minimise with five separate models each: the total number of vehicles, the total travel time, the latest arrival time, the demand weighted arrival time, and the summation of arrival time. The problem is first modelled as a deterministic model serving as a basic model for the robust counterpart. Depending on the objective functions, an adjustment is made to the deterministic model by changing the objective function and adding corresponding constraints if necessary. The same process is applied to both the deterministic CVRP and SDVRP model. Similarly, respective robust counterparts are also developed to adopt the robust optimisation approach. A two stage heuristic-metaheuristic is proposed based on a known construction heuristic and tabu search. An extended insertion heuristic is applied to the complicating constraint, namely the capacity constraint in robust CVRP formulation. Meanwhile, two additional algorithms are added in order to account for robust SDVRP constraint revolving the increase demand and travel time. A good initial solution is then obtain for tabu search to improved upon. Through their work, it is shown that the robust SDVRP model was able to reach the goal better when compared to CVRP when dealing with uncertain travel time and demands.
Ref. [84] focused on the cooperative delivery among vehiclew, where product transshipment is allowed in a dynamic setting for disaster relief aid. In the dynamic environment, a mew demand might arise in addition to a prior demand, ahen all of the vehicles are at the depot. In order to satisfy the new demand, a main vehicle with sufficient capacity to supply for any shortage of other heterogeneous vehicles is proposed. A penalty is incurred for any unsatisfied demands in order to ensure that most demand is satisfied. A mix integer nonlinear model is presented with the objective to minimise the total travel time, the penalty costs, as well as the cost that is incurred, including the loading and unloading cost. In order to solve for a realistic instance, a GA is employed, where the chromosome is adapted to the dynamic demand. The problem is decomposed by iterating the GA at every period to account for any new updates.
Scheduling and routing problems are studied in [85] in order to address both the evacuation and delivery operations post disaster. Two sets of vehicles are considered each, with their own respective task of evacuating or distributing relief aid. Both of the operation allows for split delivery when considering the limited capacity of both type of vehicles. The objective is to minimise the total arrival time, where the return time to depot is ignored neither that such actions are specified in the model. A small test instance is performed to compare the proposed integrated MIP model's performance with the sequential models, such that the evacuation phase serves as an input for the distribution phase. A Memetic Algorithm (MA) is adopted in order to solve for real problem instance. The chromosome is adapted to the two operations simultaneously resulting an integrated solution candidates.
Meanwhile, ref. [86] proposed a model that addressed the case where insufficient supply is presented instead. Because of the insufficient supply, the fairness in distributing supply among emergency locations is maximised instead. Specifically, the maximum difference between two demand points is minimised. Furthermore, the VRPTW also minimised the total travel time. In order to establish the insufficiency of supply, each vehicle is made to only serve one location, where the demand is higher than the vehicle's capacity. Given the scenario, the optimal routing of vehicles is computed while maximising fairness in the supply distribution among emergency locations. This problem is then solved by GA, which is then evaluated and verified by a small test instance.
Ref. [87] focused on distributing aids while using a vehicle that could originate from any parking area. Here, the parking area is not treated as depot, whereby the vehicle upon arrival will not immediately return. Instead, the vehicle upon the completion of task will resume a new task once ordered. Thus, a multi-trip operation is allowed. The goal is to minimise the total operation time, including the travel time, loading time, and unloading time. An Enhanced Monarch Butterfly Optimisation (EMBO) is adopted, where the modification is made on the crossover operator as well as the self adaptive strategy in order to address the problem of the solution being stuck at the local optima. The results are compared with the original algorithm (MBO) as well as seven other intelligent algorithms. A set of random generated test instance is developed in order to validate the model and solution showing better performance as compared to the original MBO and other algorithms compared.
A Fuzzy Low Carbon Open Location Routing Problem (FLCOLRP) is considered in [88], which focuses on the carbon emission while performing emergency relief distribution during a disaster, such as earthquake. In this problem, multiple distribution centres are located to serve multiple demand nodes through different types of vehicles with different capacities. The objectives are to minimise the total travel cost, total carbon emissios, as well as delivery time. The formulations of the model are done by first considering the multi-objectives from which adequate constraints are derived. Because of the uncertain demand during chaos, triangular fuzzy numbers are applied in order to estimate the demand centring around the most optimistic, pessimistic, and, most likely, demand. A penalty cost is incurred in the model should the delivery be insufficient, due to the limited capacity of vehicles per delivery. Moreover, the carbon emissions are derived from the fuel consumption, which is dependent on the vehicle speed, vehicle load, the inclination of the roads, as well as the travelled distance. A combined PSO and TS metaheuristic mechanism is applied in order to solve the problem. The initial solution is first computed by the PSO, which is later improved by the TS. Although the mechanism is regarded as hybrid in their work, the method is considered as two separate metaheuristics, since no internal operator of respective metaheuristic is changed in any manner.
The problem of LRP in post disaster phase is addressed by [89], where the medical facilities need to be setup or located at an optimal location within the vicinity of disaster affected area along with their fixed emergency medical vehicle. It then needs to route each medical fleet respectively in an optimal manner such that the delivery aid is delivered to the demand nodes minimising the time spent on providing humanitarian aid. Each fleet belonging to each k erected site is heterogeneous allowing for a more complex representation of the real case. The Emergency k-Location Routing Problem (EkLRP) next deals with routing the heterogeneous vehicle with varying capacity in such a way that the wait time is minimised. In order to execute the proposed SA algorithm, the initial solution is first obtained through the Greedy Randomised Algorithm (GRA) that decomposed the problem into two sub-problem namely the location problem and the routing problem. The GRA is also applied to adjust or improve the solution that is computed by SA in the consequent steps.
Ref. [90] focused on the shortage of supplies scenario, where fairness in distributing relief aid is demanded. The problem is modelled as LRP where the distribution centres are built just after the occurrence of a disaster. Partial road damage is addressed in the model as well as the inclusion of urgent window constraint. In this problem, the maximum demand loss is minimised, along with the total loss of all the nodes. In order to account for efficiency in delivery, another objective is to consider where the maximum duration of demand points in receiving relief aid is minimised. The problem is formulated through the lexicographic order multi-objective method. A Hybrid Heuristic Algorithm (HHA) is employed that consists of a greedy algorithm as well as ACO. Here, the fair allocation distribution is done by a greedy algorithm where the loss of demand nodes are computed iteratively. The demand nodes are also allocated to the each distribution centre via a greedy algorithm, which is also employed in order to determine the mode of delivery. Finally both the helicopter and land vehicle routing are solved by the ACO. The case of Great Wenchuan Earthquake was chosen as the study case.
In [91], the issue of relief distribution and redistribution is highlighted during a disaster to ensure that there are no severe shortages of supplies in the needed area. An IRP is used to address both distribution and redistribution issues. The first period of the operation is distribution by routing vehicles utilising the VRPPD to assign locations where the redistribution is done next according to the Network Flow Problem (NFP). Dynamic changes of updated data are considered at each multi-period planning increasing distribution and redistribution efficiency and ensure fairness in the case of increase casualties of destruction, due to the ongoing disaster. A Specialised Simulated Annealing (SSA) is proposed in order to solve the problem. The partial solution is obtained by partially solving the problem (turning it into a sub-problem due to some feasible variable fixing). This partial solution, now reduced in complexity, is solved by CPLEX for an optimum sub-solution. This subsolution is then iterated further by SA with another crossover mechanism, besides local search to reach a near optimal solution for the whole problem.
We also consider [92] as a VRP, although the proposed multi-period relief distribution model also consists of scheduling problem. Here, the problem of ideal depot location as well as the availability of transport network post disaster are being addressed for the humanitarian operation of relief distribution. Apart from a set of multi-periods and introduction of penalty update function, as well as a specified objective function, a standard modelling approach for VRP is adopted and referred to in [93] with shelters, depots and road network along with capacitated vehicles with the aim to minimise the costs and penalties. Because of the complexity of the problem, a GA is applied with a two tiered chromosome representing two decision variable: a sequence of shelter and numbers of relief distribution. The solutions could be at different lengths, where the mutation and crossover are sufficiently applied over a thousand generations that stem from 250 population size. The model and solution is then validated while using a case study based on the Nankai earthquake affecting the Minato ward of Osaka.
Ref. [94], on the other hand, raised the issue of a possible assault when transporting relief aid due to the scarcity of supply in a post disaster area where people are desperate for emergency supplies. The problem is modelled as the last mile distribution problem with known numbers for relief aid prior to delivery, which involves a network of unavail-able roads due to damages. Six criteria are considered as the objective of the model, namely: time, cost, equity, priority, security, and reliability. In order to balance all six objectives, compromise programming is applied instead of the typical approach of assigning weights on each of the objectives. An ACO is adopted to solve the problem that is based on the novel representation of ants, representing both vehicle and aid kits and effective pheromone mechanism, among four other types of pheromones. Here, the solution is constructed based on both the pheromones and a greedy approach to diversify the solutions. Using the effective pheromone and greedy function in an additive equation leads to more accuracy when constructing a solution, as opposed to the conventional multiplicative equation between the greedy and pheromone function. The model is tested on real cases of the Haiti earthquake 2010 and the Niger famine 2005.
Most recently, ref. [95] complemented the critical delivery for temporary shelter construction in times of disaster by including drones in the Multi-Objective VRP with Time Window and Drones (MO-VRPTWD) with the objectives of minimising the total energy consumption of both the trucks and drones, as well as the total number of trucks. The energy consumption is reflected by different capacities that each truck carries, the distance travelled, as well as the different speed of each drones and the distance travelled while operating. As commonly applied in various literature, a penalty is also imposed should a time window be violated. An Improved Artificial Bee Colony (IABC) is proposed to solve the problem. An improvement is made in the bee scouting strategy to improve the global search. Solomon's benchmark is applied to test the model and superior performance is observed by comparing the proposed solution with TS, GA, and VNS.

Hybrid Optimisation in Supply and Delivery VRPs
The number of literature employing hybrid methods in solving supply delivery VRPs seems to be almost comparable to those employing heuristic methods. This is probably related to the fact that, in this specified problem, the hybrids applied are often those that combine different types of heuristics and metaheuristics or merging both heuristic and metaheuristic methods into a new version of a hybrid method. However, some problems require the merging of exact solution with either the heuristic or metaheuristic method.
Ref. [96] addressed the problem of HL in terms of routing and allocation in the response phase of disaster by taking the various capacity of heterogeneous vehicle, storage capacity. as well as vehicle's working time into account. The storage capacity of the distribution points is considered while being served periodically, making it a VRP problem with dynamic demand. Due to critical nature of the delivery, a penalty is incurred for any delay based on certain penalty factors. The main objective is thus to minimise the penalty incurred thereby ensuring satisfactory delivery. The problem also seeks to minimise the unbalance distribution, making it a multi-objective problem. A hybrid solution is proposed in order to solve the problem. This algorithm executes an exact solution approach while performing perturbation in SA, ensuring the simultaneous optimisation of routing and allocating demands. The sub-problem is formulated as integer programming and it is solved using the Gorubi software.
Additionally, Cumulative Capacitated VRP (CCVRP) versions are optimised through hybrid methods in [31] that looks into multi-trip CCVRP (mt-CCVRP). Ref. [31] addressed the problem with large scale disaster, where the road could be compromised and affected for which helicopters would be an ideal choice to make a first evaluation on the affected area. Because of critical arrival time imposed on each helicopter, the sum of arrival time is reduced when compared to conventional CCVRP problem that minimises cost and distance. The problem is solved by a Multi Start Iterated Local Search (MS-ILS) with a Variable Neighbourhood Descent (VND). Each initial solution obtained by greedy randomised heuristic is optimized by VND within the working framework of MS-ILS. Through this mechanism, the exploration of MS-ILS is effectively utilised, while the VND acts as an exploitation tool for each of the initial solutions.
A similar structure is applied in [97] that studies the relief operation during the response phase of disaster. The authors went further to consider in detail the different priority of supplies, various injury level of evacuees which corresponds to their priority as well as the different roles of relief workers along with their different priorities. The objectives of the problem are to minimise the total unsatisfied demand, non-evacuated evacuees, as well as non-transported relief workers. Because of the thorough detail taken by the author, such a complex problem could only be solved by a specific heuristic tailored to the requirements detailed for larger instances. A greedy heuristic and a local search are both applied in an algorithm consisting of four phases. Furthermore, the second phase of the algorithm is executed through an exact solution approach. Based on the humanitarian operation addressed, ref. [97] could be considered as both applied to supply and delivery as well as evacuation processes.
In [98], the problem of cumulative Multi Depot VRP (cum-MDVRP) is solved by incorporating within the conventional framework of an ACO and a 2-Opt heuristic. The authors addressed the problem of the post disaster relief operation during the response phase where multi-depot facilities allows for multi-trip among vehicles. The objective function reflects the urgency of time where the delay time of total relief delivery per trip is minimised excluding the trip to depot. The problem is formulated as a MIP model, for which a hybrid solution approach, known as multiple round ACO (mr-ACO), is proposed. A 2-Opt heuristic is applied to each the best selected and constructed path created by the ants to increase the speed of the convergence of the ACO. A Tabu list is created specifically to prevent the usage of non-permitted move encountered by the ants while constructing the routes within the ACO iterations. A more stable solution is achieved with the proposed hybrid when compared with existing heuristics.
Meanwhile, through the formation of hybrid Cellular Automata (CA) and a GA in [99], a complex stochastic VRP is solved, addressing the stochastic road travel time as well as the demands. In this study, the demands and travel time during the chaotic event of a disaster are assumed to be unknown. One depot is assumed with homogeneous vehicles with the goal of minimising the total arrival time. In the proposed CGA, an additional step is added after the selected chromosome undergoes mutation, where it is further mutated into a cellular cell to expand the exploitation capability of the CGA. As a result, the convergence efficiency of the CGA is increased. This also ensures that the constraints are satisfied before the CGA resumes for the next generation.
Ref. [100] argued that the probabilistic theory should not be applied to represent the degree of belief that some uncertain events might occur. Instead, the uncertainty theory would be a better choice. Based on this argument, an uncertainty theory is utilised in order to represent two models that addressed the routing problem in distributing urgent medical provision to critical destinations. In this problem, the demands and delivery time are uncertain and their expected values are minimised or maximised in order to obtain an optimal solution of the models. A hybrid solution that is based on a proposed uncertain simulation and a GA is presented. The proposed uncertain simulation is utilised to compute the expected values, the inverse uncertainty distribution as well as the uncertainty measure from the chromosomes that are generated by GA. A numerical example is presented, showing the applicability of the proposed model, which leads to improved fairness in delivery as opposed to the conventional model that relies heavily on operational time.
The problem of emergency grain delivery in post disaster is considered by [101] where three objectives are optimised through the proposed model and hybrid solution. The satisfaction level is maximised based on the arrival time of the vehicles. The total cost consisting of demand cost, penalty cost, and travel cost are to be minimised. Additionally, the sum of transportation time is also minimised. A time window is assigned to each demand point, resulting in a penalty cost should the time window be violated. This problem is solved by an improved hybrid solution Immune ACO (IACO), where the Immune Algorithm (IA) is applied in order to distribute the pheromones based on its fast global convergence and random mechanism, while the solution is optimised based on the exploration ability and positive feedback by the ACO. The non-dominant sorting algorithm is applied to design the immunity of the ACO in order to address the convergence and diversification of Pareto optimisation. A numerical experiment that is based on Solomon's benchmark data is conducted in order to compare the hybrid performance, which outperformed the results obtained from ACO, IA, and GA.
The problem of distributing large relief material via Ship Routing Problem (SRP) in the event of disaster is addressed in [102]. The problem is regarded as an open Split Delivery VRP with the addition of varying ship's speed that decreases over time, depending on the geographical condition. The problem is formulated as a mixed non-linear programming model and the objective is to minimise the total travel time. A Multi-criteria Fuzzy Clustering algorithm is then proposed in order to obtain suitable clustering of the affected area based on the computation of dissimilarity of the demand points, the resulting fuzzy correlation matrix and finally clustering. A hybrid GA is presented as the solution of the model, where the Dijkstra shortest path algorithm is adopted in the decoding phase and solution construction. The solution scheme and model are then validated and evaluated through a numerical experiment that is based on 2008 Wenchuan Earthquake in China.
A dynamic and stochastic last mile distribution for disaster relief response is addressed in [103]. In the final stage of distribution, the coordinator faces the challenge of high demand, due to shortage of supply at the end of distribution. Furthermore, the damage that is sustained by various infrastructure, including the road, worsens at the last stage of delivery, leading to the uncertainty of demand and travel time. This leads to dynamic routing, where a pre-computed route needs to be rerouted based on the updated info. The formulated problem aims to minimise the total travel time. The problem is solved by the proposed hybrid SA and VNS. The initial solution is obtained by first allocating the nodes to the route, which is later constructed by performing nearest neighbour insertion heuristic. In order to deal with the dynamic nature of the problem, the problem is treated as a multi-stage static problem. This means that any changes, i.e adding node to the pre-computed node will be done via the Clark and Savings algorithm. This complete initial solution is then improved by the hybrid SA and VNS, where the local search in SA is performed by stochastic VNS. A real case study based on the data of Yogyakarta Indonesia is applied and the performance of the proposed solution is compared with SA, VNS, and ACO, showing comparable results.
Ref. [104], on the other hand, addressed both covering tour and the dynamic routing problem, where temporary relief shelters are erected to cover the affected disaster area while the dynamic routing aerial transport vehicle is determined in order to distribute relief aid in the aftermath of the disaster. Through the proposed model, the affected region is covered by at least one temporary relief shelter and uncertain demand, and the affected region is addressed. A virtual end node is added to the model in order to model the objective of minimising the arrival time of last assigned shelter. In solving the problem, two algorithms are applied, including the proposed hybrid algorithm. First, the known GENI based heuristic that is composed of two combined insertion heuristics is used to create route for the vehicles. Subsequently, the partition of disaster areas served by the temporary relief shelter and the assigned vehicles are determined by the hybrid Scatter Search (SS) and VNS. The local search operator is substituted with the VNS instead within the SS mechanism. Randomly generated instances are applied for a small scale problem. For a larger instance, an existing benchmark instance from [105] is modified and adopted.
The hybrid between metaheuristic and metaheuristic could be seen in [106], where the problem of fairness in relief distribution in HL is addressed in order to prevent discontent outbreak, as well as to save as many victim as possible. Four measurement of fairness are considered that are based on previous studies: absolute fairness, the min-max method, fair sequence, and multi-objectives. A Peer-induced Fairness Capacitated (PFCVRP) model is proposed in order to solve the problem. The fairness coefficient is integrated in the objective function. An improved ACO based on tabu list is proposed and considered the least used vehicle rather than applying all vehicle at disposal. The initial solution from the improved ACO is used by VNS to improve the solution further and, afterwards, local search is applied, leading to a near optimal solution. The model is validated based on randomly generated small instances which are inspired by [107]. The larger instances are also randomly generated and applied to compare the performance with ACO, GA, and ACO with 2-opt search, showing comparable results.

Evacuation in Routing Optimisation Problems
In Table A5, it could be seen that quite a number of literature addressing VRP in regards to evacuation. This section covers various evacuation types, from evacuation from bush fire through a shuttle bus via PDVRP, planning for the flow of vehicle evacuation when considering congestion traffic to in-building evacuation through path planning optimisation. This problem rarely can be seen while incorporating a combination of OR problems, such as LRP or IRP, as opposed to VRP for supply and delivery operations. However, the focus on the time window and numbers of victim saved becomes more prevalent when compared to supply and delivery operation as the time urgency is much more critical. Because there is no application of ML solution approach under this classification, the discussion of the survey is directed to the exact solution method below instead.

Exact Methods in Evacuation VRPs
The evacuation through bus routing optimisation is addressed in [108], which uniquely minimised the exposed risk of evacuees to danger by introducing a VRPTW. The objective is to minimise the time taken by evacuees to reach a safe place free from the risk of danger exposure. Similar to [109], a constructive heuristic approach is adapted in order to achieve an upper bound for the Lagrangian Relaxation (LR) method in solving the problem. The LR problem is converted into a min-cost flow problem by creating virtual arcs between pick up and delivery nodes representing potential visits. In conclusion, this work stressed the theoretical perspective of the work, rather than a practical solution, since the algorithm becomes prohibitive as the problem becomes larger, to which the application of the aforementioned heuristic is advocated.
The dynamics of real time transit vehicle routing for evacuation operations is considered by [109], where the model is updated with real information that is obtained based on the conditions of the disaster in affecting the routing network. The PDVRP is extended in order to address the transit operation to evacuate victims without transportation to emergency shelters. The resulting model is embedded in the real time transit evaluation system, SmartEvac, for humanitarian operations minimising the total travel time. The model consists of both offline and online planning. One is applied before the operation begins as a static Capacitated VRP with Pickup and Delivery (CVRPPD), while the other is applied during the real time operation as Capacitated Dynamic VRP with Pickup and Delivery (CDVRPPD), where the information updates are reflected within. In order to solve the problem, a few measures are taken in [109] to reduce the computation burden, such as introducing Column Generation to decompose the problem into master and sub-problem to address the exponential increase of variables, due to the constraint reduction. Furthermore, the dynamic multi-depot CDVRPPD was also reduced to a single depot CDVRPPD by merely introducing virtual pickup points and, eventually, a good initial solution was applied, computed by the Clarke and Wrights Saving Algorithm [110], which reduced the computation burden further. A real case study based on the Gustavo Hurricane evacuation in Gulfport was adopted and simulated while using the CORSIM network development to validate the model.
An in-building evacuation problem is studied by [111]. The Building Information Model (BIM) is converted to a graph network in order to approach the evacuation problem from the perspective of Time Dependent VRP (TDVRP). In order to determine the closest exits for each evacuees, a single shortest path algorithm is applied, such as Bellman-Ford's or the Dijkstra algorithm. The route from the group of evacuees to each of their exits is then computed from the TDVRP model. If the numbers of evacuees within the group is more than allowed, clustering method k-medoids is applied to partition evacuees in that particular group to each assigned exits. In this problem, the exit is considered as depots, while the demand points are where the evacuees are located. The problem is again solved by solving the TDVRP model formulated as a MIP with the objective of minimising the arrival of the last returning rescue team.
A more complex evacuation problem is addressed in [112], where the stochasticity of a storm scenario and loss of life, as well as uncertain travel time, is considered. Focus is given in the two-stage stochastic programming apart from the proposed deterministic model for allocation budgets for various operations that are needed during the post wildfire debris flow. The objective function carries the objective from the deterministic model, which is to minimise the expected total damage for all clustered area, where probability distribution determines the uncertain storm. In addition, the numbers of casualties and cost for picking up evacuees are added in the former objectives function in the two stage stochastic programming model. The model is then solved with Baron solvers from GAMS in order to obtain an exact solution, where a case study based on Santa Barbara Jesusita fire in 2009 is adopted.
Another location and routing problem could be found in [113], which addresses the problem of rescuing victims in post-disaster scenario. In the multi-agent system, the agents represent the centre of coordination, centre of emergency, fire station including the fire engines, ambulances, and the wounded victims. The allocations of emergency centres to the wounded are performed based on the proposed Multiplicatively Weighted Network Voronoi Diagram (MWNVD) where the number of ambulances of the emergency centre are considered as the weight of the emergency centre. A PSO-MWNVD is proposed in order to further optimise the allocation of ambulance as well as minimising the gap between the estimated victim and the expected population in the Voronoi Diagram. Based on this assignment, the ambulance is dispatched from each of emergency centre to cover the service, which are denoted by the PSO-MWNVD via the optimise path computed by the Bi Level Dijkstra algorithm. District 5 in Tehran is taken as an example of rescue implementation in order to validate the solution and model.
Ref. [114] focused on providing more than one optimal path while combining GIS information as well as the DSS framework. In this DSS, the spatial danger rating computed based on three different models (McArthur Fire Danger Index, Carvill Hurricane Index, and Earthquake Damage Estimation) is applied to the GIS for each three different disasters: fire, hurricane, and earthquake. From the danger rating, a buffer zone can be generated and overlapped onto the road network in order to determine which roads are compromised and blocked as well as possible safe destination to route to. Only the chosen road based on these considerations are taken for computation while using shortest path Dijkstra algorithm. The weight or cost reflecting the changing travel time is updated, depending on the parameters affecting the road conditions. The application of Autonomous Vehicles (AVs) for evacuation is considered in [115]. The deterministic model describes a set of AVs, persons to be evacuated, as well as a set of safe shelters. The problem is formulated as an integer linear programming aiming to minimise the evacuation time where a Dynamic Programming (DP) approach is applied to solve the problem. A divide and conquer strategy is utilised, making use of the fact that a small problem can be solved to optimality exactly. The division is achieved through clustering via K-means algorithm, where the evacuees are divided into groups, depending on the capacity of the AV and each group is served akin to TSP problem with one AV. This TSP problem is then solved by DP.
Ref. [116] extended the work of [117] in addressing both the vehicle routing as well as the flow of arrival of evacuees pattern similar to [118], which, according to the author, would improve the effectiveness of the evacuation planning. A network consisting of two sub-network is proposed where the former addressed the maximum flow as well as the earliest arrival pattern flow of evacuees, while the latter is attributed to the vehicle routing of buses in transporting evacuees to a safe shelters with limited capacity. The objective is to minimise the total time that is taken to evacuate victims from various known pickup locations. The problem is solved through the branch and bound algorithm making use of various versions of greedy heuristic in determining the upper bound in order to improve the computation time as opposed to the conventional branch and bound algorithm.
An agent based modelling is adopted by [119] in modelling the evacuation of both pedestrian and vehicle putting more focus on the evacuation by vehicle. Making full use of agent based modelling, the complex behaviour of both pedestrian and vehicles are modelled based on their respective parameters guided by realistic yet optimistic assumption. In addition, the moving speed, congestion and casualty estimation were considered in the model adding to the problem complexity. As is often with agent based modelling literature, much emphasise is put on the modelling where the model is validated based on the 2011 Tohoku Earthquake and Tsunami within Tagajyo City, Japan.
Most recently, a branch of LP that addressed the multi-objective with clear priority problem also known as the lexicographic goal programming is studied by [120]. Various conflicting goals or objectives are addressed when evacuating a group of evacuees with respective classifications. In particular, the paper addressed the objectives of evacuating as many victims as possible, reducing operating time as well as minimising cost where no trade-off should be considered. In their model, the flow of evacuees and the vehicle assignment to a set of nodes are made compatible guided by the constraints proposed. To reduce the trade-off between the conflicting objectives, the goal lexicographic goal programming is adopted. The evacuation of those with specific needs are considered to have highest priority followed by maximising the evacuees with normal conditions, minimising evacuation time and finally reducing the operation cost. In this goal programming which is solved by GAMS 23.7, the optimal value of highest priority leads to a new constant that is considered for the next goal and will continue until all objectives are achieved according to priority level. Finally, the model is tested with the case of the 2018 earthquake and tsunami in Palu, Indonesia.

Heuristics and Local Search in Evacuation VRPs
An interesting problem is addressed by [121] where a livestock evacuation plan is optimised based on the heuristic version of a Lagrange Relaxation method. Inspired by the earthquake and tsunami events in Fukushima in 2011, the livestock is evacuated via vehicles to prevent contamination from radioactive exposure known for causing harmful effects on bone marrow and thyroid. The problem is formulated into two models each to maximise evacuated livestock as well as minimising the vehicle numbers respectively. The first model addressed both dynamic setting and uncertain road accessibility which is addressed through probability obtained by simulations. Meanwhile, the second model which is said to be more realistic is based on a deterministic VRP model employing graph theory as opposed to the network flow model of the first model. One decision variable is proposed namely the number of evacuated livestock traversing the arc at time t in the first model. Meanwhile, the second model addressed two decision variables where the decision of vehicle travelling an arc and decision to allocate capacity to each demand node are proposed respectively. In this research, ref. [121] further argued that in times of emergency, a sub-optimal evacuation plan would be more realistic as opposed to optimal solution computation which would consume valuable time. Thus, the complex integer programming models are both solved by utilising a modified Lagrange Multiplier heuristic which is considered simple yet efficient in solving the problem. This method solved the problem via Lagrange Relaxation partially until a satisfying evacuated livestock number is reached. The computation is then halted to allow for the rounding up phase where the fractional computed flow is rounded based on the worst-case scenario. A similar concept is also applied to solve the second model. The models and solutions are then validated through generated instances with comparisons made based on the solutions computed by Gorubi. A case study based on the earthquake and tsunami events in Fukushima in 2011 is applied, where the transportation network and the number of livestock are referred to. A safety distance is determined as well as the time limit for evacuation in the first formulation. The second formulation, on the other hand, is based on a fixed fleet of vehicles that is determined based on the feasibility of the solution for the second model. Near-optimal solution quality is observed, which is comparable to the solution that is provided by Gorubi.
The evacuation planning focusing on both the routing and flow of evacuation on the generated optimised path is considered by [122]. An optimised flow of evacuation eliminates the occurrence of congestion that could hinder an effective evacuation plan as a whole, even if the routing is optimised. The problem is modelled through an evacuation graph with nodes representing evacuated node, safe nodes and transit nodes as well as edges representing roads. The objective is to maximise the number of evacuees reaching the safe nodes via vehicles. This problem is solved by a Conflict based Path Generation (CPG) heuristic that was developed to decouple the problem into sub-problem and master problem. Here, the sub-problem generates an optimised path and, as a result, adds a new column and row to the master problem. The master problem, on the other hand, solved the optimised flow in response to the generated path via the scheduling of the evacuation. The cycle of sub-problem and master problem iteration is repeated until convergence is observed based on a determined characteristic. The Hawkesbury Napeon case study is adopted, where a flood evacuation is necessary once every 200 years involving 70,000 victims. The CPG with and without contraflow is compared, showing greater improvement in the former.
A CVRP with simulated demand as part of the three stages in developing a DSS integrating the GIS is studied in [123]. In this problem, the simulated demand is obtained through the population estimation algorithm that estimates both the criticality of the demand as well as the demand itself. However, the demand or the population is deemed to be uncertain, as in the case of the population in Phuket, which mostly consists of vacationers and locals, depending on the season. As such, the GIS enables the population estimation to be applied in order to execute the Monte Carlo sampling to estimate the demand of a node. In this model, the total distance travelled is minimised. The tsunami prone Phuket is identified as potential emergency location, where its population is estimated as a form of evacuation demand. The problem is solved by constructive heuristic while using the Clark and Wright Savings algorithm.
Meanwhile the aerial vehicle routing for an evacuation plan is addressed in [124]. The characteristic of the helicopters as rescue vehicle and the urgency level of the victim is captured through the concept of agent-based modelling. Based on this concept, the rules of interaction between the rescue agent-also named as the vehicle agent-and the demand agent-or the victim-is proposed. A realistic representation of the problem is reflected based on these characteristics. A nonlinear optimisation problem is formulated, where the objective is to maximise the number of saved evacuees. This problem is solved by the proposed heuristic namely the Cooperative Multi Agent (CMA) based algorithm where the rules of interaction is considered in solving the problem. The initial solution computed by the vehicle agents based on the information obtained from the demand agents and the safe location is evaluated by the demand agents in order to improve the quality of the solution by means of replying what is to be removed or inserted in the task list. Vehicle agents then update its task list by means of removing or inserting demands into the task list and this new task list is further evaluated. This process continues for as long as there is a request of updating task list from the demand agents. A numerical experiment is performed in order to validate the model that is based on the 2011 earthquake and tsunami in Japan.
The wild fire evacuation optimisation process is considered in [125], by taking the dynamic changes of the route due to spreading fire into consideration. Based on the flow model depiction, a single destination is considered from multiple evacuation locations. All of the evacuees are assumed to start evacuation at the same time, which is regarded as group evacuees, where they share the same computed route for evacuation. The Capacity Aware Shortest Path Evacuation Router (CASPER) is used interchangeably with the Capacity Aware Reverse Map Analyser (CARMA) in order to compute the route for all evacuees. The objective is to minimise the final departure time of evacuees with CASPER being assigned to the load dynamic and static input of the problem, while CARMA transforms the network GIS data into a graphical network and seeks for the solution based on the single Dijkstra algorithm of each shortest path tree. While the computation is improved in searching for the routes, the observed egress time is not minimised, thus an iterative heuristic is proposed where a repair mechanism is applied to a pre-computed solution with disadvantages with a new paths until an applicable solution is obtained. In this heuristic, the conflicting and competing two potential paths over edges is defined. A worst case scenario of wild fire in California is adapted in order to test the model and algorithm.
In [126], the problem of evacuating evacuees with disabilities is addressed based on their categories by a fleet of heterogeneous vehicles. Here, the vehicle ranged from ambulances, vehicles that are capable of carrying only normal evacuees, as well as vehicles for carrying evacuees with wheelchairs. A DVRP is considered for this problem, which is formulated as an MILP with the aim to minimise the total evacuation time, where, for the case of multiple optimal route, the route with the least travel cost is opted. The problem is solved by a proposed heuristic that is based on two concepts. The first concept develops route sequentially for vehicles according to their capacities in descending order until every vehicle exhausts its capacity. On the other hand, the second concept assigns vehicles to pickup nodes in parallel with descending capacity to corresponding descending demand at pick up nodes. The two concept heuristics are compared in both experimental instances as well as a real case of a forest fire in Teruel, Spain. The results shows that the first concept is better in terms of the total distance covered, while the second concept is more efficient in minimising the total evacuation time.
A complex representation in addressing the evacuation during a bush fire is shown in [127] through the Multi Destination Capacitated Dynamic VRPTW (MDCDVRPTW). In this study, the problem of evacuating late evacuees to a safe shelter with the aim of maximising the numbers of late evacuees through a low risk route is addressed. The late evacuees reflect a realistic view on the event, where not all evacuees could board the vehicle to the shelter within the specific time window, due to the capacity of the vehicle. A GA and two heuristics are presented to solve for this complex problem. A comparison made shows that the constructive based heuristic exceeds the GA in terms of solution quality. Additionally, the exact algorithms are applied to be served as a comparative test, which again highlights the advantages of heuristic over exact solutions when dealing with larger instances. This work could also be classified under the Metaheuristic in Section 4.2.4 because this work applied both GA and heuristics as the solution methods.
The evacuation planning between the warning time of flood and the time when the flood occurred is addressed in [128]. In this problem, the evacuees are evacuated via public transit bus to a safer locations with the aim of maximising the number of evacuees saved with minimum cost. Two models are presented based on the planning steps that are divided as pre-and during the flood event addressing the problem of choosing a pickup location as well as the problem of computing an optimal route from the pickup location to a safer place. In the FLP, a few candidates from a pickup location are considered, such as the bus stop, which everyone is familiar with, the location where huge population is detected, and a minimum number of existing bus stops. The FLP is formulated, such that the distance between pickup locations and their corresponding shelter is minimised. Meanwhile, the routing problems are formulated as the Split Delivery VRP, where one single shelter is denoted as the depot. Pickup locations corresponding to the shelter are the demand points and the bus is only allowed to return when its capacity is full. The problem is solved via sweep heuristic algorithm, where the simulations are executed and the performance is evaluated.
In [129], a two-stage problem is formulated in order to address the issue of locating, scheduling, and routing evacuees during evacuation. The first stage dealt with determining the optimal pickup location based on the distance of the depot and safe shelter. The second stage then solved the problem of scheduling and routing of the vehicle in transporting evacuees from the computed pickup locations to the shelters. Three objectives are achieved through this solution approach, where, in the first stage, the walking time of the evacuees to the pickup locations is minimised. While, in the second stage, the total transportation time and the number of vehicles used are minimised. The first stage solution stems from clustering method that clusters the evacuees into pickup points. This problem is based on the uncapacitated multi-facility location problem, which is also known as the multi-Weber problem. These computed pickup points from the IP model is then used as a parameter, together with their corresponding time window to the second stage problem, which is a MILP model. A Hybrid GA (HGA) that is based on the clustering algorithm and a GA is applied to improve the clustering quality of some other methods that are discussed in their work. This proposed algorithm addressed the problem of noise sensitivity of Fuzzy C Means (FCM) in performing clustering. Furthermore, three aspects are considered when solving the second stage problem where the route path should be determined, the evacuation flow should be selected as well as scheduling the vehicles. An Interval Round Trip based Routing and Scheduling heuristic is then proposed based on these considerations.
A Noncombatant Evacuation Operation (NEO) is considered by [118] to facilitate the decision makers in evacuating civilians as well as other noncombatant in South Korea to a safe location. A time staged network model is proposed as a MILP model with the objective of minimising the required evacuation time. The network structure consists of features from the perspective of a second infantry division that is usually tasked with overseeing the evacuation process. From the MILP model, the computed optimal evacuation flow of every arc is then applied in order to solve the heterogeneous CVRP using the Clark and Wright heuristic involving 2-Opt improvement. In the subsequent VRP, constraints are introduced to be applied to each transportation mode (bus, train, and helicopter)m with the objectives of minimising the total evacuation cost where the cost is equivalent to the travel time. Furthermore, the determination of the transportation mode at every interval of fix time leads to a dynamic VRP.
Recently, ref. [130] solved the problem of evacuees evacuating with the assumption that each could evacuate via their own vehicle. The multiple disaster locations are addressed with multiple evacuation pointsm leading to multiple sources and multiple sinks, respectively, in the resulting pre-min cost flow model formulation. The min cost flow problem is then obtained by converting the multi-source and multi-sink to a supersource and a supersink. Because of the uncertain travel time and road capacities, which depend on different scenarios, a robust optimal routing model that is based on the two-stage stochastic programming model is proposed. However, the model consists of a hard coupling constraint that could not be solved in a reasonable time. This leads to decomposing the problem into two single stage sub-problem via a heuristic based on the Lagrangian relaxation. In order to validate the model and the Lagrangian relaxation based heuristic algorithm, the Chicago sketch and grid network are applied for experimental test.

Metaheuristics in Evacuation VRPs
The problem of evacuation of victims with disabilities is addressed in [131]. This problem is classified as an Overburdened VRP (OBVRP), where not all victims can be evacuated, due to the overburdening of resources. This is reflected in the objective function that tried to minimise the total number of unevacuated victims. Consequently, typical swapping heuristic operations are no longer feasible in addressing the problem. This is due to the fact that each possible route must now contain unevacuated victims. To deal with the problem, different versions of ACO extensions that consist of Min-Max Ant System (MMAS), Elitis Ant System (EAS), as well as neighbourhood ant heuristic and greedy ant heuristic were employed. Each of them corresponds to creating upper and lower bounds for the pheromones, providing an emphasis on a strong solution from the start of the search serving as a good initial solution. Because there is no available and existing standard benchmark for validating the model and solution, newly generated instances are derived from [132]. The obtained results highlighted the importance of proper parameter selection for the algorithm.
Ref. [133], on the other hand, addressed the dynamic VRP, where the evacuation routing from a local shelter to a long term safe settlement is updated at each 15-minute interval. Here the travel time taken along the route, the capacity of vehicles, and the numbers of evacuees left at the local shelter are dynamically updated in order to capture the essence of dynamic changes during disaster. The route planning need to be executed in such a way that minimised the total travel time which results in reduced evacuation time. The two-stage solution procedure assists the evacuation routing, where the clustering of shelters to each depot is performed in the first stage. The second stage executes a routing optimisation via a SA with an initial solution obtained by the nearest neighbour heuristic. Stage 1 and 2 are both updated at each 15 min. interval, eliminating empty local shelters and updating vehicle capacities. The transportation network of north west of Tehran was utilised as a case study to test the model where the static and dynamic model's solutions are compared.
A Multi-Objective Multi-Trip CVRP (MO-MTCVRP) is discussed in [134] when addressing the evacuation operation with limited resources when considering multiple, yet conflicting, objectives. The assumption is made that these limited resources should be used as many times as possible during a critical time. Furthermore, a fast response in the last mile problem is ensured by the objective of minimising the maximum latency, which reflects the waiting time of the last affected victim in the operation. Additionally, the total cost and number of vehicles are minimised in order to ensure an efficient relief operation. To solve the complex problem, a Multi Start algorithm with Intelligent Neighbourhood Selection (MSINS) is proposed. Here, a number of neighbourhood structures undergo the two phases of a Multi-Start Algorithm, where the feasible solution is first created and later improved. Based on the principle of optimality that was introduced in the previous study, a multi objective local search is then executed in the two phases of a Multi-Start Algorithm. This solution approach is compared with the NSGA-II algorithm and the results indicate that the proposed solution outperformed the NSGA-II in both computation time and solution quality.
Ref. [135], on the other hand, raised an interesting notion of using the autonomous heterogeneous bus fleet to evacuate victims from an affected area to a safer location. The objective is to minimise the total travel cost that consists of the cost incurred by the travel time as well as the cost related to the risk while evacuating victims. In this problem, a conventional formulation of VRP is presented with the exception of a heterogeneous fleet. The two associated costs are then combined into one cost that is based on different weights assigned to each cost. An evolutionary algorithm is proposed to solve the problem. A case study is used to evaluate both the modal and solution based on Hillsborough County, Florida, where hurricane Irma occurred in 2017. From the results, an encouraging number of cost reduction is obtained as compared to the manual evacuation effort.
The evacuation problem due to typhoon is considered in [136] by looking into the secondary re-occurrence and the chain effect of the typhoon. The evacuation process is divided into two phases. In the first phase, the evacuees are transported to the emergency shelter, while, in the second phase, the injured victim will be transported to the hospital. Heterogeneous vehicles are made in order to scatter and their priority are based on their tasks given. Furthermore, the secondary disaster, such as the resulting mud slide and flood, are reflected through the damage of the road affecting the travel time of the vehicle. This phenomenon is denoted by the road impedance coefficient which is directly correlated to the damage of the respective road. The aim of the model is to minimise the time that is spent by vehicle throughout the operation. A matrix encoding based GA is proposed in order to solve the problem. The problem model and solution are then evaluated by a simulation that is based on the real Meranti Typhoon that occurred in 2016, affecting as far as Xianmen, China as a super typhoon.
Meanwhile, ref. [137] addressed the issue of solving hte minimax problem for evacuating evacuees due to a disaster in two stages, resulting in two objective functions. The first stage involves transporting evacuees from an affected area to a temporary shelter resulting in an objective of minimising the maximum evacuation time of evacuees. The second stage involves transporting the victims from a temporary shelter to a safe area where from the perspective of decision maker, the total evacuation time for all evacuees is minimised. No direct connection is given between the affected areas, the temporary shelters and the safe areas. Solving the nonlinear minimax optimisation problem is challenging, where a unique solution is sought for in this specific problem as compared to other problems in the literature. The reason for this being that a multiple solution namely the routing guide would confuse the evacuees who are obtaining the guide from the decision maker. As such, a lexicographic minimax MIP model version of the problem is considered, where the tabu search is employed. The model was then validated by the real network of Honolulu, Hawaii while using two test scenarios where the demand of evacuees is different.
The efficiency of evacuation by simultaneously optimising the location of evacuation hub as well as the routing of public transport via two layers target model is studied by [138]. A public transit network is optimised for evacuation by introducing the reversible lane for the bus transit as well as the private vehicles for evacuation. Assumptions were made, such that the capacity of evacuation hub as well as the demand is based on refugee numbers. It is also assumed that there is no delay in the evacuation process. In this problem, any suitable trunk road can be considered as the bus reversible lane increasing the capacity of the flow in order to improve the evacuation. Moreover, the far right of the lane is reserved as the reversible bus lane, so the activity of evacuation can be done back and forth. The problem is modelled as the min cost flow model, where the total travel distance and the total evacuation time are minimised. A multi-objective GA is then performed, where the objectives are to improve simultaneously. A case study of the Beilun district in Ningbo is adopted for the case of typhoon evacuation in order to validate the model.

Hybrid Optimisation in Evacuation VRPs
In [139], the evacuation planning is addressed through the LRP with Time Windows (LRPTW), where the problem is solved through a hybrid Memetic Algorithm (MA). The setting presented is the evacuation process in an urban area, where a bus pickup is arranged for evacuees in the event of an earthquake. Each emergency shelter has homogeneous capacity that is represented by nodes that are analogous to the multi-depot in VRP, while the evacuation points each have a certain number of evacuees represented as demand. In this problem, the decision on location of evacuation point along with their respective evacuees number and assignment to which emergency shelter is considered. An optimal route is computed in order to minimise the total evacuation time. The Hybrid MA is employed where the early convergence is avoided through the introduction of a heuristic and randomness. The known neighbourhood operator CustMove is adopted within the local search procedure to improve the offspring. The model is validated through a case study that is based on a seismic disaster in Bucaramanga in Columbia.

Rescue in Routing Optimisation Problems
Besides routing for supply and delivery as well as routing for evacuation operations, VRP for rescue operation also dominates the literature in terms of addressing humanitarian operations during disasters. This can be seen in Table A6. However, their numbers pale in comparison when compared to the two aforementioned humanitarian operations. Despite that, in the 2015 Nepal earthquake disaster, for example, the calls for improving Search and Rescue (SAR) mission could be found in [26] where road blockage and landslide proved to be the largest obstacle in SAR missions. Ref. [27], on the other hand, pointed out the need for better coordination among local authorities in SAR missions, despite the fact that the aid operations are largely dependent on international aid. Following the call for more efficient SAR mission during disasters, various publications could be found addressing the routing problem involving SAR missions, which are further elaborated in the sections below.

Machine Learning Methods in Rescue VRPs
Similar to Section 4.1.1 under the supply and delivery VRP, this review only encountered one work that solves the rescue VRP through the ML approach. As the case of [29] in Section 4.1.1, the application of ML here also addressed the stochastic and dynamic of the problem. In the case of SAR, such an application can be found in [140], where a Cooperative Receding Horizon (CRH) based VRP is formulated in order to allow for cooperation among agents. This application of a solution concept is applied to an emergency case, where the ambulance is routed to specific emergency cases under the condition of uncertain road network or cases due to a disaster or other hazard. In this ML approach, the agents learned to optimise the route by accumulating rewards. The sum of rewards are maximised in a dynamic and stochastic environment by providing stationary trajectories for agents that represent the ambulances to their assigned targets under certain conditions, following the work of [141]. The authors expanded this application by applying the concept on graphbased routing as well as developing a behaviour switching algorithm for the agents to adapt to the changing targets over time, allowing for cooperative rescue or service among agents. Furthermore, the switching algorithm attempts to address the occasional oscillation that might occur when applying CRH, which causes a non-stationary trajectory for the agents. It is worth noting that the proposed algorithm tried to make use of both versions of CRH, namely the CRH and the target oriented CRH (tCRH), in achieving the optimal effect by determining the trigger point to switch behaviour, based on the distance between agents and the targets. Meanwhile to exploit both aforementioned versions, switching CRH (sCRH) is proposed. Several real scenario's data obtained from the local emergency agency are applied in order to validate the model where the results are compared between CRH, tCRH, and sCRH, showing best performance in sCRH solutions.

Exact Methods in Rescue VRPs
The problem of overwhelming demand for emergency medical help that is snowballing during the time of disaster is studied in [142]. A bi-objective location allocation and ambulance routing is addressed in order to determine the location, the number of temporary emergency shelters, and the optimal route for ambulance in rescuing the victim that need medical assistance. In this problem, the demand is divided into two types of patients that are coded with red and green colour. The red type patient is considered to be critically injured and should be transported to the temporary emergency station. The green code patient, on the other hand, receives treatment on site. The objective is to minimise the complete service time of treating both red and green code patients, as well as to minimise the cost of setting up the temporary emergency stations. By first linearising the non linear constraint, the constrain method is then applied to solve the problem, where the test instance is based on the case from Tehran city.
Ref. [143], in his work, made use of the public monitoring system service provider, known as the Golden Eye, to provide proof of real time monitoring of the transport network condition to guide the emergency vehicle into the heart of disaster affected area during the post disaster event. Based on the real time data pre-and post-disaster, a custom AI makes a prediction on the traffic and road condition in guiding the decision makers. Moreover, the internal Golden Eye program analysis is utilised to fill in the blind spot due to damaged monitoring equipment on site. Based on the existing expert input on damage and possibility of repair, an optimal route for the emergency vehicle is then dynamically computed while using an improved Bidirectional Dijkstra optimal path selection algorithm. For a path that failed to be monitored due to faulty equipment, a prediction is then made that is based on data pre-and post-disaster at a regular interval. As a result, these predicted paths could also be included when computing an optimal path. The rich topological demographic of Fujian province in China was considered when verifying the model. The problem of emergency vehicles dealing with road conditions, such as the connectivity reliability of road, the computed rescue time travel, as well as the index road security is considered in [144]. For the travel time, the model is built considering the reduced capacity of a link as well as the increased flow of the link due to a sudden disaster or emergency event. Two correction coefficients are introduced to incorporate into the existing Bureau of Public Road (BPR) function formulation to address the sudden nature of a disaster onto the traffic travel time. The connectivity reliability is defined as the measure of the possibility of the link being used under specified conditions and times. The third variable addressed the safety of the emergency vehicle while in transit during disaster event, especially those occurring in intersection, which is frequently the case. The Pareto optimal solution is applied in order to address the multi-objective problem. The depth search is adopted to search for feasible solutions from which the sub-problem (based on the three consideration's value function) is updated. Furthermore, the alternate solution is obtained through the non-dominated sorting method that is based on the predefined conditions. Eventually the problem and solution method is validated based on the traffic accident that occurred at the intersection of Yunnan Road and Beijing Westroad, Nanjing China under heavy downpour.
Ref. [145] focused on the efficiency of rescue routing by integrating two separate models namely the Data Envelopment Analysis (DEA) model with the Incomplete Linguistic Group decision Constraint Cone (ILGCC) and an IP model to compute for optimal rescue routing. The former model incorporated five decision variables to determine the rescue efficiency. These variables consisted of input variables representing the distance of the arc, the risk of travelling along the arc, as well as the output variable that represents the affected population, the level of damages sustained by buildings at the site, and the possibility of a second disaster occurrence. Taking the missing of information and data from the experts into account, the ILGCC model was developed and adopted in developing the DEA model. The ILGCC reflected the expert opinions and decisions that are often incomplete and vague as a first response when a disaster occurred. A rescue efficiency could be computed and applied to the IP model with the five decision variables and by applying the DEA model with ILGCC. Furthermore, in the dynamic setting, the changes of data are expected and applied for every period. In the IP model, the total rescue efficiency is maximised with an ongoing decrease of efficiency, as the operation time is prolonged. The case study of Wenchuan earthquake that hit Sinchuan, China is taken to generate instances in order to verify the model. The efficiency of SAR is highlighted as compared to the time-based solution approach.
The problem of saving victims that are trapped in various critical locations due to a disaster that occurred by allocating and routing the SAR team in the post disaster phase is addressed in [146]. In reflecting a realistic scenario, both dynamic and stochastic elements are presented in the proposed model. The SAR team is heterogeneous with different operation capacities. Additionally, two versions of the online problems are proposed with two different objectives. The first problem tried to minimise the operation time (known as makespan) as well as minimising the total weighted latency before the victim is finally rescued. The weights are computed based on the percentage of numbers of victims at a critical nodes. The second proposed online model tried to minimise the total waiting time of the victim until they are finally rescued. These two online models are then compared to the corresponding offline model, where all of the information is known. Finally, a deterministic online strategy is proposed, making use of the solution obtained through the offline MIP version of the problem. Ten generated networks are applied to validate the model; they are obtained from Gorubi, each with instances up to 500 nodes.
More recently, the problem of both allocating resources and routing the resources for SAR operation via a two-stage decomposition approach was studied by [147]. The model also addressed the risk of a secondary disaster occurrence. Furthermore, the uncertainties are also considered through the environment setting, which is addressed by the interval based robust optimisation approach. In the first stage, there are several affected districts that consist of several demand locations with varying demands, depending on the damage afflicted. The model addresses the fairness in distributing allocation by computing the minimum total demand coverage through estimating the SAR demand in each afflicted district. Each demand is weighted by its importance weight and the objective is to maximise the weighted sum of the demand minimum's coverage for all types of demand, which, in turn, maximises the fairness allocation. The output of the first stage, namely the optimal SAR resource allocation is then forwarded to the second stage in the routing problem model. The total weighted demand completion time is minimised. Solving both of the MIP problems, an interval data robust optimisation is adopted as a solution approach, where each robust counterpart is defined for each problem model to address the uncertain parameters. A few scenarios are considered from the case study of earthquake in Iran in order to verify the models where the robust models are compared with a deterministic counterpart in terms of the solutions obtained.

Heuristics and Local Search in Rescue VRPs
Since 2010, there have only been two studies found in the literature that used heuristic and local search in the domain of rescue VRPs. In [148], the operation stage of the five stages of SAR are addressed, where the maritime victim survivors in the ocean are searched for and rescued by nearby ships and helicopters operating simultaneously and dispatched by the Maritime Rescue Coordination Center (MRCC). In this problem, the longer the arrival time of the rescue team, the more risk of the survivors being dispersed and drifted apart although; in this problem, the location of the survivor is known. A Generalised VRP (GVRP) is adopted in order to represent the problem where multiple depots are considered. The objective is to maximise the profits that are associated to the lives saved during the rescue operation. Because the exact solution performs poorly for a larger instance of the problem, two constructive heuristics are presented to solve the problem, with one being greedier than the other. The first heuristic, which is the sequential construction heuristic, adds a possible next destination that is based on the constraints given. However, the second heuristic constructs the route by considering all feasible next destinations and finally choosing one that has the highest return; hence, the term greedy. These heuristics were applied and they resulted in a pilot heuristic, where, instead of adding and constructing the route one-by-one, a partial solution computed by the pilot method was extended through a sequence of possible clusters. A historical incident that is handled by the local navy was taken as the case study to validate the proposed models.
The dynamic or multi stage static CVRP with deadlines to emphasise the importance of rescuing the victims in floods was considered by [149]. The deadline is presented as a hard constraint, as opposed to a time window that only incurred a penalty in the objective function. A time flow between demand request and treatment time is added to the objective function weighted by priorities. Additionally, four priorities are defined, namely: the victim that can remain on the spot, victim that need to be rescued within 12 h, victims that need to be rescued within 6 h, and victims that need to be rescued immediately. From these priorities, the dateline and weight factors were determined. The objective of the problem is to minimise the total time flow, which is weighted by the priority of respective demands. The problem is formulated in LP and tested with various scenarios. In this dynamic setup, a multi-static problem is solved offline based on the new updated demand at each stage. Because of the complexity of the real scenario reflected, a MILP formulation is presented and solved by the proposed heuristic. The problem is solved by the proposed Best Time Flow Insertion heuristic, where it is adapted from the well known Best Fit algorithm in the bin packing problem.

Metaheuristics in Rescue VRPs
A Large Neighbourhood Search (LNS) is used in [150] in order to address the critical time in achieving optimal ambulances routing during critical event of rescuing injured patients. In this problem, the injured patients consist of those heavily (red code) and slightly (green code) injured. These two groups of patients will be treated in hospital or on field, demanding different courses of actions. With the aim to minimise the weighted sum of the latest service times of both red code and green code patients, considerations, such as hospital capacities and priorities of the two groups of patients, are taken into consideration. The service is considered to be complete if the red code patient arrives at any hospital or the green code patient is treated with first aid. The service time addresses the preparation time for a red code patient for transportation to the hospital. Meanwhile, the service time for a green code patient is denoted as the treatment time at the scene. The problem can be considered to be a multi-depot with multiple hospitals equipped with homogeneous ambulances. The solution is addressed in two mathematical models with two index decision variables and three index decision variables, respectively. The former aims to reduce the burden of computation. A LNS is applied in order to solve for the complex problem in a fast manner and robustly. Two initials solutions were obtained from the insertion and constructive heuristics, leading to the excellent quality of solutions. No existing benchmark is available for the problem addressed. Hence, a large set of instances is generated to investigate the performance of the models with better performance in favor of the two indexed formulation model.
Additionally, ref. [151], in their work, also focused on the ambulance routing problem through a minor extension of basic VRP. The deterministic problem aims to minimise the total travel cost, which is limited by each ambulance's budget for operations. Based on the emergency medical services to victim rescued or Points of Care (PoC) and the number of ambulances, a feasible solution is obtained respecting the capacity and limitation constraint imposed by the model. A GA was implemented to solve for near optimal routing. The test case based on the information that was obtained from a public hospital in Jendouba, Tunisia was adopted where the occurrence of a disaster, such as floods or winter storms are common.
The SAR mission utilising UAVs as working robots with battery limitations is studied in [152]. The limitation is handled through locating the charging station, where the UAVs is required to pause their activity upon reaching critical battery level. Hence, both the location of the charging station and the trajectories of the UAVs are optimised in this problem. In order to reduce the computation burden, a two-layer level of grid map is utilised to address the trajectories as well as placement of the charging station. A GA is applied, where genes represent mission points in the sequence of the order of the genes that need to be visited, while the chromosome is evenly divided to accommodate multiple UAVs. Based on the computation, a single battery limitation can indicate the final visit that a robot is able to make. In this problem, the completion time and the operation cost, as well as the distance between the charging points, are minimised, while a penalty is incurred for violating the energy limitation via the weighted sum method. The solution method is evaluated via Monte Carlo simulation in Matlab. Furthermore, a Gazebo robot simulation is applied to evaluate the performance further.
The problem of road damage and congestion during a disaster is addressed in [153] by incorporating the short term traffic flow prediction. This unique solution approach consists of multiple solutions and modelling stages. The real traffic data are managed through the internet based data gathering program, as well as other existing real data gathering systems, such as the traffic control, the vehicle tracking program, and the road toll system. These rich data are converted into a fuzzy model, where the attributes or indicators are represented by the fuzzy set, whose membership degree is denoted by the value of the data attributes. Through this fuzzy information system, the traffic congestion can be predicted and applied when computing the optimal route for the emergency rescue vehicle via the guided local search with a variable penalty function. In order to validate the proposed model, a test data based on 13 of the main roads of Lanzhau is utilised. In addition, the test instance from previous literature was applied and the results were compared with different solution algorithms, including the traditional local guided search. The stability and quality of the proposed algorithm was proven as well as showing superior results when compared to the traditional guided local search.
Ref. [154] focused on improving the saving of lives and reducing human suffering in the allocation and routing problem. In this problem, three humanitarian objectives are formulated and captured with regards to lifesaving utility, fairness, and delay cost. The rescue operation was considered within the first 72 h of post disaster event. Two decision variables are introduced, each addressing the allocation and routing component, respectively. In this problem, the allocation decision variable addressed the number of rescue victims that should be loaded onto the vehicle at the demand point. The objective was given to minimise the total travel time, which is also the total rescue time, in order to minimise the delay cost that is associated with the deprivation cost and maximise the lifesaving utility reflecting priority of wounded victims for efficient lifesaving. This problem was solved by two metaheuristic algorithms, namely the NSGA-II as well as the proposed NSGA-II with ACO. Both complex and simple examples are applied in order to validate both methods. In the analysis, it is shown that the proposed mutation and crossover operators lead to a more quality solution when compared to their original counterpart.
The ambulance allocation routing problem was considered by [155], where the rescue operation of the patients was based on three categories, depending on the severity and the required services. In order to emphasise the urgency in saving the patients, a soft time window constraint is introduced, whereby the thresholds are based on the survival function of the patient category. A penalty was added for the late arrival and the objective is to minimise the weighted sum of latest service completion time over all patients treated, the penalty time imposed while treating the red coded patient, and the penalty time that is imposed to treat the yellow coded patient. The three different categories of patients are colour coded, symbolising non-critical patients that can be treated on the spot, severe patients that should be transferred to hospital, as well as severe patients that need to be transferred to hospital while sustaining his or her life with equipment within a certain type of ambulance, respectively. Four types of ambulance are defined, depending on the severity of patient. GA and TS are both applied to solve for the realistic problem size, due to the vast and complex problem that is being addressed. Both of the solutions are compared with a modified greedy heuristic and CPLEX solver to evaluate the results obtained. Better performance and efficiency are observed in the proposed TS method.
Ref. [156] addressed the emergency relief routing for injured victims by rescuing and dispatching them to the medical centre. Two humanitarian points are considered and compared. The first point is the varying degree of injury as opposed to the assumption of all injuries are identical, which is depicted as one of the objectives. The second point is the varying degree of suffering tolerance while en-route to the medical centre reflected by the varying upper bound of time window constraints. Human sufferings are reflected in the model through the absolute deprivation cost as well as the relative deprivation cost. Based on the problem setting, three models are proposed with regards to identical injury level, diverse injury level, and diverse injury level with a hybrid transportation strategy for mitigating the issue of the diverse injury level of victims. In this model, the transportation cost and the absolute and the relative deprivation cost are minimised. The second model considers a different priority of a severely injured person when compared to a slightly injured person, while the hybrid model incorporates strategy that allocates some space for a slightly injured person to be rescued along with the severely injured victim if the capacity of the vehicle allows it. All of the proposed models are validated by comparing the model with and without the proposed additional objectives, as well as the models with and without proposed varying time window constraint. The problem is solved by ACO in order to verify and evaluate the model showing the improved level of fairness in rescuing victims as well as the necessity to consider various degrees of injury when performing a rescue operation.
Recently, Ref. [157] addressed the problem of distributing injured victims over three transportation network consisting of road, rail, and air transport while considering the allocation and location problem of the operation. Different levels of effects from a disaster are depicted through the circled based approach. Through this estimation approach, the travel time, percentage of injured victims as well as the damage to the link is estimated depending on the radius of the circle from the disaster origin point and their relation formulation, respectively. Some assumptions are made, auch as capacitated and limited vehicle depending on region, non-regional residence that is also considered as candidate for additional aid stations, and the homogeneous type of injury. The objective of the LRP, formulated in the integer nonlinear programming model, is to minimise the operation time, consisting of establishing total aid station time as well as the relief time. A GA and a discrete imperialist competitive algorithm are proposed and compared when solving the problem. A case study of supposed earthquake in Tabriz, Iran is referred to with a discrete imperialist competitive algorithm showing better performance when compared to GA.

Hybrid Methods in Rescue VRPs
The Urban SAR Team Deployment Problem (USAR-TDP) is addressed by [158]. For cases, such as rescuing victims trapped in a damaged building due to the disaster, it is appropriate to model the problem with stochastic demands, travel time, as well as service time. This reflects the nature of such key characteristics that cannot usually be known a priory. As the USAR operations prolong, the situational awareness of the post disaster environment increases while the survival of the victims on the other hand reduces to reflect a more realistic scenario. The approximation technique can usually be applied to reduce the Multi-stage Stochastic Programming (MSP) into a single stage. The MSP is decomposed into a sequence of inter related two-stage stochastic programming with recourse where the previous decisions and information of the former two-stage programming with recourse is used for the subsequent two-stage problem. A column generation method is devised in order to reduce the large two-stage stochastic programming at every decision epoch where the inclusion of each possible tour for a team is regarded as one column. The decomposed sub-problems are then solved through VNS and, if needed, through an exact solution from the branch and price method. As such, we regard the solution as a hybrid solution strategy. The earthquake that struck Port-au-Prince in Haitian Capital in 2010 is taken as a case study to validate the model applied.
The problem of rescuing injured victims in the time of a disaster is studied by [159] in the form of non fix Multi Depot Pickup and Delivery VRP (MDPDVRP), where the survival rate of the victim is maximised. In this problem, simple assumptions are taken such as the homogeneity of the vehicle with fixed speed and capacity. An Emergency Reply Center (ERC) plays the role of assigning vehicles to patients to be rescued, depending on the distance. Furthermore, ERC also provides critical information regarding the road condition in the aftermath of the disaster. In contrast to the typical approach of maximising the service, the focus here is to maximise the number of saved life. The model is solved by a proposed hybrid GA, where the chromosome signifies the nodes and the order of the genes represents different hospitals. An improvement is made on the feasible solution from the permutation through the local search. A case study that is based on Tehran city is applied in order to validate the solution and model. GA and the proposed hybrid GA are both compared where better performance could be seen in the proposed solution.
The problem of rescuing victim during storm flood disaster in urban cities is considered in [160], where the chance constraint programming specifying a confidence level under which the stochastic travel time constraint is guaranteed to hold. In this problem, a multi-objective with an assigned weight is addressed to minimise the total travel time budget of all vehicles as well as the maximum travel time budget of a single vehicle. A simple model is assumed, where there is only one depot, and that the vehicle could only visit a site once with no return trip is allowed. A lognormal distribution is applied to model the uncertain travel time based on the real data obtained. In this problem, the vehicle unloads a rescue operative as well as facilities before resuming its route onto other disaster sites. As such, this problem could also be considered as a supply and delivery problem.
To solve a realistic problem size, a hybrid of ACO-TS is proposed, where the addition of TS is applied in order to improve the solution that was obtained by the ants before the pheromone is updated. By doing so, the immature local optima obtained by the ACO is avoided. The model is then validated through a real storm flood that occurred in the urban area of Guangzhou, China.

Current Trends
Throughout the survey, various disasters are considered, such as earthquakes, floods, tsunami, drought, mudslide, hurricane, typhoon, war, and bush fires, with some work covering the possibility of secondary subsequent disaster, such as in [136,147].
By observing Tables A4-A6, it could be observed that the application of VRP in humanitarian operations largely revolves for supply and delivery problems with 73 papers, followed by evacuation problems with 32 papers and, finally, rescue operation problems with 22 papers. Figure 4 also clearly depicts this. Out of the 123 papers reviewed, only three papers solve the routing problem involving recovery phase. Two other papers addressed the disaster preparedness phase operations, while one paper looked into the mitigation phase operation. These trends in DM phases are observed in Figures 5 and 6. Meanwhile, the breakdown of DM phase according to supply and delivery, evacuation, and rescue operation is depicted in Figures 7-9, respectively.      Furthermore, only 16 papers address the slow onset disaster type, such as flood (Figures 10 and 11). Six papers addressed both rapid and slow onset disaster types, four of which applied for supply and delivery operation, while two addressed the rescue operation, as illustrated in Figures 12 and 13, respectively. Meanwhile, Figure 14 shows the break down of disaster type for evacuation operation. The rest of the papers addressed the sudden onset or rapid onset disaster, as can be seen in Figure 10. Only 28 papers could be considered as RVRP, according to the definition that is presented in Section 3.3, with a breakdown of 13 papers (46.43%) in supply and delivery operation, 12 (42.86%) in evacuation operation, as well as three (10.71%) in rescue operation addressing real world problem in a complex setting ( Figure 15). Meanwhile, 79 papers addressed deterministic VRP, with one work focusing only for one vehicle. Moreover, only 26 works addressed SVRP, followed by 30 addressing DVRP, of which only seven addressed the SDVRP. The overall distribution of these model characteristics according to selected humanitarian operations could be seen in Figure 16, while Figure 17 shows the distribution of all reviewed papers against model characteristics. Surprisingly, the breakdown of the these papers shows that more rescue operation addressed both stochastic and dynamic problem (four papers), followed by supply and delivery operations with just two and evacuation operation with the remaining one paper (Figure 18). This suggests that the setting for rescue mission might be more complicated and complex, due to the life that is involved in the operation, which resonates well with the study presented in [10]. Additionally, only 49 literatures addressed multiobjectives VRP, with more than half (69.39%) applied in supply and delivery operation, followed by the rescue operation (10.20%), and finally, the evacuation operation (20.41%) ( Figure 19). In all three humanitarian operations, deterministic models are addressed slightly more than one-third of the papers each. This is depicted in Figures 20 and 21, as well as in Figure 22, where the summary of all three is seen in Figure 17. Additionally, the classification of papers reviewed based on deterministic VRP, SVRP, and DVRP is shown in Figures 23-25, respectively.                Meanwhile, from the solution approach perspective, only two papers solved their VRP through the ML solution approach (Figure 26), each for supply and delivery operation and rescue operation, followed by 16 papers that adopted the hybrid solution approach with breakdown of 11 papers for supply and delivery operation and three each for evacuation and rescue operations ( Figure 27). In addition, 35 papers employed the exact solution approach with slightly more than half addressing supply and delivery operation, followed by seven papers addressing rescue operation, and nine for the evacuation operation ( Figure 28). Meanwhile, 25 papers adopted the heuristic solution, which saw an almost equal share of both supply and delivery operation and evacuation operation with 12 and 11 papers, respectively, followed by only two papers for rescue operation. Finally, 50 papers applied a metaheuristic to solve the VRP with both evacuation and rescue operation having an almost equal share of papers, with 10 and nine each, respectively, leaving supply and delivery operation with 31 papers. When looking at the heuristic and metaheuristic individually, the breakdown percentage could be observed in Figures 29 and 30, respectively. Interestingly, the results show that the number of works adopting hybrid solution is only slightly lower when compared to heuristic solution in supply and delivery operation (Figure 31), while metaheuristic continues to dominate the works that address VRP in disaster applications, followed by the exact solution approach (Figure 32). Furthermore, these numbers echoed the trend in [15] in terms of ranking by mostly applied solution approaches if hybrid solution and ML are excluded. The only difference is the reduced dominance of metaheuristic usage among the reviewed papers. That being said, this review focuses on the specific application of VRP in humanitarian operations, as opposed to the broad field of applications covered by [15]. Figure 33 provides the trend of solution approaches based on the selected humanitarian operation.        When looking into each of the operations individually, the supply and delivery operations show a clear dominance by the metaheuristic in terms of solution approach with roughly 42 percent of the papers belonging to this class. This is followed by the exact solution approach with approximately 26 percent of the papers, leaving both heuristic and hybrid solutions with almost an equal share of percentage. As can be expected, the ML solution approach only accounts for about one percent when compared to all solutions applied. In details, the finding is illustrated in Figure 31. Likewise, a similar trend could be observed for the rescue operation in Figure 34, leaving the evacuation operation showing a different trend. Interestingly, here, the heuristic is applied more, albeit only three percent more, than the metaheuristic, followed by the exact solution that roughly triples the value of hybrid adoption in percentage. This is shown in Figure 35. Finally, the current trends of modelling and solution approaches applied in the scope of this survey based on year are presented below in Figures 36-39 for clearer depiction.      In general, the modelling of classified VRPs in this survey can be divided by those adopting graph theory (such as the minimum cost flow model and spanning tree model), agent based model, such as in [113,119,124], among others, the set partitioning model, as can be seen in [51], and the MDP model. There is also a trend of the combined OR problem that also addressed the VRP to fully capture the important mechanism of humanitarian operations, such as the LRP, which can be seen in [42,44,48,73,74,77,79,89,138,139], followed by the LIRP in [80,81], scheduling and routing in [85,92], ARP in [96], as well as ALRP in [57].
From these models, various interesting and complex problems are addressed. As a general trend, it could be seen that most work addressed split delivery, such as in [34,42,46,47,[49][50][51]61,71,73,75,76,78,81,83,85,87,94,96,102,116,123,128,158] and multi-trip VRP that can be seen in [31,47,50,61,76,88,90,94,97,98,124,126,129], ref. [134,149,159] normally stressed on the problem of vehicle limitations such as the numbers of vehicle or the fuel limitation. However some works addressed the VRP in this scope by assuming flexible fleet size such as [30,87]. Additionally, heterogeneous fleet as can be seen in [75] is usually assumed for works that deal with complex requirements for the selected routing operations. It is also not uncommon to see multi-depots, such as in [70,75,80,87,94,98], and ref. [82] being applied to address a more realistic scenario rather than depicting only one depot. Understandably most literature involving ambulance routing are often addressing VRP for rescue operation, which could be seen in [150,151,155]. These compliment the works that address urgent medical needs of victim that can be seen in [113,150,151,[155][156][157]. In particular, road condition is one of the most important key factors when addressing routing in humanitarian operation. Various conditions of road can best be modelled through a flow model that details the flow of the road, capacity, as well as impedance. Such detail can be seen in [121,125,136,153], to cite a few. Some works make use of the dispersion of disaster effects, as can be seen in [157] that link the road's condition based on circled areas of the affected region. Some works assume binary capacity of road, such as [114,146]. On the other hand, uncertain road capacity could be observed in [67,75]. Meanwhile, ref. [76] addressed the problem of repairing some roads while executing delivery operation.
Furthermore, various applications of the fuzzy method are observed, such as in [42,44,68,77] and [88], which address uncertainty on the parameters into the fuzzy parameter. The fuzzy method is also applied to address multi-objectives such as in [34]. Meanwhile, a multi-criteria fuzzy clustering is applied to estimate certain parameter as can be seen in [102].
What are also worth noting are works that address inter-modal network, such as in [33,148], involving maritime and land transport as well as air transport, respectively. Furthermore, ref. [53] addressed the cooperation between air and land transport. More depictions of inter-modal networks are seen in [33,52,90,104,118]. Additionally, the routing problem is not unique to just land vehicles. Some works, such as [60,64,74,152], only address drone routing, while [63,124] addresses the helicopter routing. To that note, the trend of cooperation between drones and vehicles in enabling more sophisticated emergency operation could also be observed, as can be seen in [64,95]. Interestingly, autonomous vehicles are also addressed, such as in [115] as well as ship routing in [102].
Meanwhile, several works addressed their work as an explicit DSS. For example, ref. [39,47,52] in supply and delivery operation as well as [114,123] in evacuation operations. These types of applications usually also utilise the aid of external information, such as [52,123], which applied the GIS, as well as [111], which aid evacuation through the BIM. Other research, such as in [102], make use of the geographical topography information to compute decreasing speed during the humanitarian operation.
Some other works worth mentioning are [37], which addresses the fairness in terms of sharing relief aid among shelters. The importance of fairness is also highlighted in [38,51,90,106]. The importance of including livestocks in an evacuation plan is highlighted in [121], while the mitigation phase of a disaster is addressed by [45] in restocking earthquake shelters. Ref. [55] synced the supervisory visit at a medical centre with the resupply of the pharmacy through PDVRP. Ref. [60] tackles the problem of dynamic route conditions via the incorporation of drones and land vehicles with the drone acting as the eye. Ref. [62] addressed the delivery process during drought, which is uncommon in the literature, with [72] addressing the blood donation problem. Meanwhile, a unique consideration is taken in [94], which addressed the possibility of an attack from desperate victims by constraining the delivery operation only via convoy. Ref. [112] addresses the flow of debris, which causes harm during evacuation, while [126,131] focuses on evacuating victims with disability. On the other hand, ref. [116] first addressed the flow of evacuees during evacuation before considering the routing problem and, finally, refs. [143,144] addressed both the connectivity network and the respective emergency vehicle.
The application of ML as solution approach is truly lacking in this scope of research. A complete approach of ML could only be seen in [29], which incorporates MDP as part of the reinforcement learning approach. Ref. [140], on the other hand, introduced learning by maximising rewards and applied this to agent based modelling instead. The exact solution approach on the other hand shows a frequent application of the Dijkstra algorithm, CG method, A star algorithm, as well as many other unspecified commercial solvers from CPLEX, Gorubi, and LINGO, among others. It is interesting to note that all of the papers applied commercial solver instead of using the open source platform, such as GLPK, COIN-OR, and SCIP. Furthermore, it could be observed that the classified VRP that is solved by exact solution approach ranges from a simple problem, such as [32] with only one vehicle, to complex problem solving for multi-objectives, inter-modal, multi-trip, fuzzy environment, stochastic, dynamic problem, as well as a combination of the OR problem, such as LRP and ARP.
The same could be said for the heuristic solution approaches in solving problems for various ranges of complexity. Among the well known heuristics applied in this scope of survey include Clarke and Wright heuristic, greedy heuristic, sweep heuristic, local search, and sequential insertion heuristic. Other heuristics that are not commonly applied in this scope include Best Time Flow Insertion heuristic and Lagrangian Relaxation based heuristic. Furthermore, a tailored heuristic could be seen in [124], with the cooperative multi-agent based algorithm heuristic.
GA is commonly applied as a metaheuristic solution approach. Here, various versions of GA could be seen, such as NSGA-II, matrix encoding based GA, and parallel GA. Some metaheuristic are especially applied for solving multi-objectives such as NSGA-II and MOPSO. Furthermore, ACO is popularly adapted, especially for the case where swapping and mutation could not yield feasible solutions. Such is also the case for objectives with strong associations as can be seen in [94]. Additionally, it could be noted that heuristics are commonly applied for an initial solution prior to improvement by the metaheuristic seen in [83,133]. However, it is not uncommon to see metaheuristic being used in a subsequent manner, such as in [88,89,154]. Meanwhile, TS is commonly used as part of the improvement phase of other metaheuristics, rather than used in isolation. Other metaheuristics that were applied within the scope of survey include SA, ABC, MA, GRASP, LNS, GRA, PSO, and MBO.
From the perspective of hybrid approaches, ACO remains popular as the main component for a hybrid solution between metaheuristic approach. Especially, the application of TS, commonly as an improvement operator for ACO, can be seen in [106,160]. Meanwhile, the hybrid of ACO and other Immune Algorithms are seen in [101]. Furthermore the hybrid of GA and other heuristic can be seen in [159]. Additionally, a hybrid between metaheuristic and heuristic is also observed, such as in [98,139]. Furthermore, the hybrid of exact solution and metaheuristic or heuristic are seen in [96,97,102].

Future Direction and Research Gaps
The survey findings indicate that the metaheuristic will continue to dominate the solution approaches for the VRP, in particular for the application of humanitarian operations. This is due to the tendency of the resulting model becoming more complex, as more important criteria should be addressed to describe the chaotic situation of the humanitarian operations. Naturally, metaheuristic would be the obvious choice to solve these complex models in a practical sense. To echo the work of [15], a considerable number of RVRPs have started to emerge, especially in this specific field of application. It is expected that more RVRPs would surface within the field ensuring prolonged relevance of metaheuristic solutions.
The advancement of computing power, on the other hand, will ensure that the exact solution remains relevant, as the strong appeal for obtaining an optimal solution rather than settling for a near-optimal solution will always be tempting. For this very reason, the growing number of adoption of an exact solution will continue to increase proportional to the advancement of computing power.
In this paper, the reviewed literature are differentiated between cooperative mechanism among solution approaches and hybrid solution approaches, where one solution approach might be incorporated within another solution approach. Nonetheless, the hybrid solution approach is slowly gaining traction due to its relevance in overcoming the limitation of existing solution approaches. The dilemma between optimality and practicality remain prevalent within the context of optimisation. The idea of obtaining an optimal solution for practical practices will continue to be desired. Rather than waiting for the advancement of computation in the near future, making use of existing solution approaches towards reaching this goal is also an obvious step. Such efforts would surely give rise to more hybrid versions of solution approaches.
Meanwhile, the natural modelling of both dynamic and stochastic by MDP, which often leads to the ML solution approach, could no longer be ignored. In representing the chaotic situation in pre-and post-disaster settings, the uncertainty of various parameters will remain an important issue. Furthermore, the changing variable within an interval is commonly observed and utilised in humanitarian operations. Looking at the trend, it will soon be a must that DSVRP be addressed for all humanitarian operation purposes, should effective delivery, evacuation, and rescue operations be executed. As such, it is anticipated that more hybrid versions of ML would emerge in addressing the curse of dimensionality haunting the ML solution approach.
In terms of DSS, various online tools, such as GPS, BIM, and GIS, will soon be widely incorporated, as can be seen in [52,111,123]. This would ensure the continued relevance of DVRP. Perhaps more online tools covering various aspects of humanitarian operations will soon emerge, giving rise to more potential research to cater for this new development. Moreover, the advancement of computation also allows for ML to advance further in assisting humanitarian operations in terms of vehicle routing. Such a decision support system would perhaps be able to imitate the human decision-making process, as observed in the game of Go.
Finally, we curiously note the abundance of work covering the routing problem in terms of delivering critical aids to victims. When compared to such problems, evacuation routing problems and rescue routing problems are rather scarce. In the RVRP setting, such problems are even rarer. Furthermore, most of the reviewed papers addressed the response phase of DM, leaving the other three phases with, at most, three papers reviewed out of 123 papers. A similar observation could be made for papers addressing different disaster types. Research addressing slow onset disaster in this scope is still lacking. From the modelling perspective, the deterministic problem is addressed the most. This could point to the missing link between the theoretical and practical side of the routing in the selected humanitarian operations, as it could be argued that the deterministic model rarely accounts for a real-life problem that is dynamic and stochastic in nature and mostly applied for theoretical purposes. Finally, as it is often mentioned in this study, the adoption of the ML approach is very much lacking as far as the research gap is concerned. This leads to some research questions for future research: 1.
Could the ML modelling and solution approaches be utilised to its full potential (overcoming the curses of dimensionality) to address a more complex RVRP in detailing accurate depiction of humanitarian operations for practical instances? 2.
How to overcome the limitations of the existing solution approaches utilising hybrid solution in addressing RVRP for practical instances? 3.
How to incorporate the emerging online observation and communication tools in the optimisation computation of user-friendly DSS for the applications of humanitarian operations? 4.
How to incorporate as many elements of the chaotic environment as possible into the evacuation and rescue routing problems for the humanitarian operations in the setting of RVRP? 5.
Could a general model be developed that applied to all disaster phases as well as disaster types?
Guided by these research questions, more focus should be paid in developing a much richer VRP that incorporates as many elements as possible in the event of disasters. To this end, more work should be done in making use of the MDP modelling approach that naturally addresses the respective dynamic and stochastic elements, as would be desired for the humanitarian operations. As pointed out by [15], we have yet to see more multiple stochastic and dynamic parameters being addressed simultaneously. As such, more effort should be made to look into such problems. A richer RVRP should entail an intermodal network, as the humanitarian operations normally utilised all available modes of transportation in achieving effective outcomes. All of these qualities demanded by such RVRP would no doubt lead to a more complex model where solving it would be challenging. A hybrid solution approach could be the key in solving such a problem. Concurrently, the ML solution approach could also play the same role. Despite the potential of Deep Neural Network in handling a complex task, we have yet to see one, such as in [161], addressing the RVRP in the settings of the humanitarian operations. Hence, we should also expand the optimisation world of OR beyond the exact, heuristic, metaheuristic, and hybrid solution approaches alone. The emerging ML approaches should also receive an equal amount of focus given to the rest of the existing solution approaches. Furthermore, more online and communication tools should be incorporated and translated into the dynamic and stochastic elements of RVRP. These online tools serve as a source for dynamic and stochastic elements that are obtained from the location of disasters in depicting the real situation of the affected locations. Finally, evacuation and rescue routing problems both deserve an equal amount of attention among the academicians, if not more. The huge gap between these problems and the delivery routing problems being addressed shows a disconnection between the theoretical and practical worlds.

Concluding Remarks
In this review, a collection of 123 papers is reviewed and classified accordingly based on the proposed classifications. Out of the 123 papers, four papers are classified twice under the classified scope of VRP routing operations. Furthermore, current trends, future research, and research gaps are also discussed after reviews of all the collected papers. The findings show that metaheuristic is applied the most for the selected routing problems. On the other hand, the exact solution approach remains relevant to the field, despite the computation challenges. This is attributed to the advancement of computation power. More hybrid solutions are seen applied in order to improve and overcome the limitations of the existing solutions. In the same manner, more technological tools are incorporated in order to holistically address the routing problem by providing more inputs and realtime observation. This leads to a more challenging model addressing stochastic and dynamic VRP, giving rise to the notion of applying a machine learning solution approach. Few research gaps are also identified from the perspective of modelling and solution approaches, such as the incorporation of existing tools in providing a more effective decision support system, modelling richer problem, especially for evacuation and rescue operations that are still lacking, as well as the inclusion of both dynamic and stochastic elements in the model. Finally, slow onset disasters, such as war, famine, and drought, deserve equal attention, if not more, to further contribute to humanitarian operations.