Resilience-Driven Road Network Retrofit Optimization Subject to Tropical Cyclones Induced Roadside Tree Blowdown

Abstract

This article focuses on decision making for retrofit investment of road networks in order to alleviate severe consequences of roadside tree blowdown during tropical cyclones. The consequences include both the physical damage associated with roadside trees and the functional degradation associated with road networks. A trilevel, two-stage, and multiobjective stochastic mathematical model was developed to dispatch limited resources to retrofit the roadside trees of a road network. In the model, a new metric was designed to evaluate the performance of a road network; resilience was considered from robustness and recovery efficiency of a road network. The proposed model is at least a nondeterministic polynomial-time hardness (NP-hard) problem, which is unlikely to be solved by a polynomial time algorithm. Pareto-optimal solutions for this problem can be obtained by a proposed trilevel algorithm. The random forest method was employed to transform the trilevel algorithm into a single-level algorithm in order to decrease the computation burden. Roadside tree retrofit of a provincial highway network on Hainan Island, China was selected as a case area because it suffers severely from tropical cyclones every year, where there is an urgency to upgrade roadside trees against tropical cyclones. We found that roadside tree retrofit investment significantly alleviates the expected economic losses of roadside tree blowdown, at the same time that it promotes robustness and expected recovery efficiency of the road network.

1 Introduction

An increasing concern for the transport sector is that natural disasters can cripple the functionality of a transportation system, especially for road transportation (Mitsakis et al. 2014). Roadside trees are a necessary component in a complete road system.¹ Well-distributed roadside trees have many observed benefits. Such trees can enhance the amenity of driving through the landscape by alleviating the sun and wind and preventing snow drifting, can relieve the intensity of flooding by slowing down and absorbing storm water runoff (Neale 1949), and can improve the air quality on roadways and the visibility of roadsides via cutting down the dust levels on roadways (Zeigler 1986; Wolf 2003). Roadside trees also increase the occurrence of traffic accidents such as tree crashes, however, which has been a major hazard source threating road safety in the United States (Wolf and Bratton 2006) and Poland (Budzyński et al. 2017). In the face of tropical cyclones associated with strong winds, roadside trees will be easily blown down, causing severe injuries or fatalities and disrupting roadways. Prior to the occurrence of tropical cyclones, roadside trees should be planted to improve their survivability and mitigate the disruption of road segments as a kind of pre-disaster investment decision making.

This article first explores optimal retrofit investment in road networks subject to roadside tree blowdown induced by tropical cyclones. We follow a four-step procedure to pursue optimal decision making: (1) tropical cyclone modeling; (2) roadside trees fragility analysis; (3) economic loss and network performance assessment; and (4) optimization of road network retrofit investment in areas facing roadside tree blowdown.

In the optimization step, a trilevel, two-stage, and multiobjective stochastic model was developed to make optimal decisions for retrofit investment. The upper level optimizes retrofit decision making, the middle level determines the optimal repair sequence, and the lower level assigns traffic on the road network. The proposed model at least is a nondeterministic polynomial-time hardness (NP-hard) problem, which can be solved by a trilevel algorithm. A nondominated sorting genetic algorithm (NSGA II) was used to determine a set of nondominated optimal retrofit schemes at the upper level of the algorithm. A genetic algorithm (GA) was implemented to determine the optimal repair sequence of damaged road segments under a specific damage pattern in the middle level. The Frank-Wolfe algorithm (Frank and Wolfe 1956) was adopted to solve the traffic assignment model given a road network at the lower level. To decrease the computational burden, the random forest method was applied to reduce the trilevel algorithm into a single-level algorithm. The highlights of this article include three aspects: (1) a new metric is provided to assess the performance of a road network in which both the damaged status and the nondamaged status are included in one metric; (2) resilience is formulated in two objective functions from two aspects, namely robustness and recovery efficiency; and (3) the combination of the trilevel algorithm and random forest method saves a great deal of running time.

The remainder of the article is organized as follows. Section 2 is a literature review while Sect. 3 introduces the study area and data collection. Section 4 describes the four-step research procedure. The computed results are illustrated and analyzed in Sects. 5 and 6 presents the conclusions of this article.

2 Literature Review

The state of the art in road network retrofit optimization mainly focuses on the optimal allocation of resources for highway bridges subject to seismic hazards (Chang et al. 2012; Rokneddin et al. 2013; Dong et al. 2015; Morbin et al. 2015; Zanini et al. 2016; Fawad et al. 2019; Akiyama et al. 2020; ChienKuo et al. 2020), because highway bridges are vulnerable to earthquakes and crucial to road systems. Some research considered other road infrastructure retrofit s, such as signs, signals, and lights, which are vulnerable to tropical cyclone impacts (Lambert et al. 2013). To the best of our knowledge, however, there is no existing work that considers the impact of roadside tree blowdown on road network performance in the face of tropical cyclones.

The criteria considered for evaluating retrofit schemes in the existing literature include retrofit investment costs (Fan et al. 2010; Huang et al. 2014; Dong et al. 2015; Morbin et al. 2015; Lu et al. 2018; Döyen and Aras 2019), transportation travel costs (Fan et al. 2010; Huang et al. 2014), topological aspects (Rokneddin et al. 2013), deteriorated bridge fragility (Rokneddin et al. 2013), component performance and systematic resilience (Venkittaraman and Banerjee 2013), economic loss (Dong et al. 2014), sustainability (Dong et al. 2015), and travel time reliability (Zhou et al. 2019). In this article, we also consider the limited investment budget, economic risk of roadside tree blowdown, and road network resilience in optimal retrofit decision making.

Previous studies treated road network resilience as an integrated criterion; we consider resilience through two criteria: robustness and recovery efficiency. Low levels of robustness of damaged road networks can cause serious deterioration in road network performance. If the repair process cannot be launched immediately, adverse social and economic consequences will be induced. Recovery efficiency is an integrated product of the post-disaster damage status and the repair scheme. Retrofit schemes are usually implemented before the realization of hazards, which influences the economic loss caused by roadside tree blowdown and the resilience of road networks after the occurrence of hazard impacts, thus further affecting recovery efficiency. Separate consideration of the enhancement of robustness and recovery efficiency can find a better retrofit scheme to make road networks perform better and be more resilient.

Many existing studies used performance metrics to evaluate the functionality of road networks. For example, Latora and Marchiori (2005) employed the length of the shortest path to assess transportation network performance. Equilibrium travel time was selected as the performance measure of a highway transportation network by Ukkusuri and Yushimito (2009). Günneç and Salman (2011) assessed the functionality of a transportation network facing a disaster by developing multiple probabilistic metrics of connectivity and expectation of travel time/distance. Yang et al. (2018) evaluated the performance of road network using three metrics, namely weighted connectivity, distance efficiency, and total weighted travel time. Jenelius and Mattsson (2012) adopted travel time to evaluate the performance of the original shortest route and the new shortest route during disruption. Almoghathawi et al. (2017) employed several flow-based measures to reflect network functionality. Overall, most of the researchers measured road network performance through computing connectivity, travel distance, traffic flow, and travel time. Travel time is an integrated metric that considers connectivity, travel distance, traffic flow, and user behavior within its computation. However, the existing literature does not consider travel demand and travel time simultaneously in an appropriate way. For severely damaged road networks, a large portion of road segments are dysfunctional, which causes disconnectivity across the network. Only a small part of travel demand can be satisfied. However, the reduced travel demand could release traffic congestion in some road segments that have large traffic in normal situation. In this case, the travel time of the satisfied demand will be reduced. For the less damaged road networks that are still able to serve all original travel demand, travel time could be increased due to more congestion in road network. Only considering travel time or treating travel demand and travel time with the same weight cannot reflect above phenomenon. There is a realistic need to propose a suitable metric for road network performance under disasters.

A useful tool to measure the uncertainty of adverse consequences from hazardous events is risk assessment. Expectation of loss (Peeta et al. 2010; Salmerón and Apte 2010; Günneç and Salman 2011; Du and Peeta 2014) and conditional value at risk (Noyan 2012; Lu et al. 2018) are two widely used risk measurements. Expectation of loss determines the location of average loss distribution, and conditional value at risk (CVaR) can reshape the loss distribution (Sarykalin et al. 2008). The CVaR is a variant of value at risk, which quantifies the scale of expected losses once the cutoff point (value at risk given a significant level α) has been breached (Rockafellar and Uryasev 2000). In our study, expectation of loss is used as a measurement to quantify economic risk of roadside tree blowdown. The CVaR can be used to measure the functional risk of a road network to reduce the chance of having poor network robustness after the occurrence of a hazard. Robustness is the opposite of the functional loss experienced by a road network, which is defined as the surplus functionality of road network after roadside trees are blown down. The lower the robustness, the higher the functional loss. We applied a similar idea from CVaR to measure the expectation of a specific percentile of the worst robustness, and named it as conditional value at risk in terms of robustness. Figure 1 shows the concept of CVaR in terms of robustness, which represents the average of extreme low value of network robustness below cutoff point.

Fig. 1

Demonstration of the conditional value at risk (CVaR) concept in terms of robustness

3 Study Area and Data Collection

Hainan Island, which is in the southern part of China, was selected as our study area. Many highway segments on Hainan Island are planted with different species of roadside trees. The average distance between two adjacent roadside trees is around 20 meters. Hainan Island suffers severely from tropical cyclones with 162 cyclones reaching landfall from 1949 to 2015; 20% of these storms were strong or super-strong typhoons. During tropical cyclones, strong winds often lift or blow down roadside trees, impeding the use of some road segments. For example, on 18 July 2014, Typhoon Rammasun made a landfall on Hainan Island in South China with a central wind speed up to 60 m/s. In Haikou City, more than 100,000 roadside trees were destroyed, disrupting 133 road segments (Haikou Municipal Administration Bureau of Haikou City 2014). There is an urgency on Hainan Island to upgrade the survivability of roadside trees in order to decrease roadside tree blowdown and improve road network performance under tropical cyclone conditions.

The tropical cyclone best-track database for the western North Pacific (1949-2015) was downloaded from the China Meteorological Administration (Ying et al. 2014). The topology of the road network on Hainan Island was digitalized from a map (Star Map Press 2014). An origin-destination traffic matrix was compiled from the 2014 China City Statistical Yearbook (National Bureau of Statistics 2014), which includes passenger traffic and cargo data. Average daily traffic data in 2014 were obtained from the Transport Information Center (TIC) of Hainan Province. The TIC also provided roadside tree data, which include spatial distribution of tree species, and survey data on tree blowdown over provincial highways induced by Typhoon Rammasun.

We focus on roadside tree retrofit at the provincial highway level instead of the street level. The three reasons are: (1) speed limits of provincial highways are much higher than street roads, and drivers may not have enough response time to stop their vehicles when there exists roadside tree blowdown, which threatening the lives of people; (2) When traffic jams happen on provincial highways, drivers cannot easily change origin lanes that are blocked by roadside trees and the capacity of affected highways decreases significantly; (3) Redundancy of a provincial highway network is much lower than a street road network. Drivers in provincial highway networks find it more difficult to determine a short detour than when they are in a street-scale road network. Less redundancy degrades the performance of provincial highway networks.

4 A Four-Step Procedure for Roadside Tree Retrofit

A four-step procedure is used for decision making of roadside tree retrofit, which includes tropical cyclone track modeling, fragility analysis of roadside trees, economic loss and network performance assessment, and optimization of road network retrofit investment. The following subsections introduce each step in details.

4.1 Tropical Cyclone Track Modeling

In this part, we model the cyclone tracks across Hainan Island. The process of cyclone track modeling is presented in Fig. 2, which involves a number of steps:

Fig. 2

Process of cyclone track modeling

• 1. Obtain the landfall locations of historical cyclones. In previous studies, scientists employed a buffer zone to define the impact area of a cyclone (Lin et al. 2010; Malmstadt et al. 2010). For a specific study area, its buffer zone was designed according to its shape and size. For Hainan Island, 113 km is set as the average radius of historical wind speed whose Beaufort level is equal to or larger than 10. For each historical cyclone track, its origin is the coordination of movement when the cyclone began to affect the island.

• 2. Model relevant factors of origin. Based on statistical analysis of historical cyclone tracks on Hainan Island, the annual frequency of historical cyclones, cyclone occurrence times, and central pressure differences follow general extreme value distributions (GEV). There is an obvious positive correlation between the latitude and longitude of origins (Kendall correlation coefficient is 0.4996), which is fitted by a joint empirical copula probability distribution. Central maximum wind speeds of historical cyclones were modeled by their own empirical distribution, because no existing distribution fitted them. Inverse transform sampling was employed to generate random numbers for these cyclone factors.

• 3. Model coordination of movement. The statistical model used to calculate the coordination of cyclone movement was identified from the literature (Rumpf et al. 2007).

• 4. Model central pressure difference and central maximum wind speed. The central pressure difference (CPD) and central maximum wind speed (CMWS) are the two main factors of cyclone intensity. Affected by the surface roughness, both CPD and CMWS will exponentially diminish after the landfall of a cyclone. Based on statistical analysis, CPD and CMWS respectively follow the attenuation model proposed by Vickery (2005) and Kaplan and Demaria (1995).

• 5. Model ending locations of cyclones. Most of the cyclones affecting Hainan Island passed through the island and only a few of them ceased on the island. Here we define the ending location as the coordination of movement when the cyclone began to leave Hainan Island. If the central maximum wind speed of a cyclone satisfies one of the following conditions, the movement will stop: the central maximum wind speed is zero; there is no historical record within a preset range (330 km); and the time step is larger than the simulated cyclone duration.

• 6. Validating cyclone track simulation. The simulated cyclone tracks can be verified following two aspects: First, by examining cyclone tracks. Five longitude lines (18.5°, 19°, 19.5°, 20°, and 20.5°) and four latitude lines (109°, 109.5°, 110°, and 110.5°) were selected as basic lines, which exactly cover the whole of Hainan Island. Comparison of the number of the historical and simulated cyclone tracks that encounter each basic line: for longitude basic lines, the absolute error is 17.4 and the relative error is 13.2%; for the latitude basic lines, the absolute error is 31 and relative error 20.5%. Overall, the absolute error is 23.4 and the relative error is 16.7%. Second, we examined the intensity of cyclone tracks. An accumulated cyclone energy index (ACE) was adopted to evaluate cyclone intensity. The ratio of simulated ACE to historical ACE is between 0.689 and 1.32. Comparing this with existing studies (Hall and Jewson 2007; Rumpf et al. 2007), the accuracy of cyclone tracks is acceptable.

Following the above procedure of cyclone tracks generation with Monte Carlo simulation, we produced a set of 1,020 tropical cyclone tracks for 500 years. The China Meteorological Administration released typhoon meteorological service regulations in 2001, which extended the Beaufort scale up to Level 17. Table 1 lists wind speeds at Beaufort Levels 10 to 17. Given a cyclone track and according to the Miller wind field model (Miller 1967), the impact area at each Beaufort level during a tropical cyclone can be estimated (Yang et al. 2016). For example, Fig. 3 displays Typhoon Rammasun’s impact areas (radius of wind, RW) at Beaufort levels from 10 to 17 in 2014.

Table 1 Wind speed at different Beaufort levels

Fig. 3

Impact area of Typhoon Rammasun on Hainan Island at Beaufort wind levels from 10 to 17. RW radius of wind

4.2 Fragility Analysis of Roadside Trees

Winds associated with tropical cyclones are strong enough to uproot or blow down roadside trees, which may impede traffic on road segments. According to the data collected by the Transport Information Center (TIC) of Hainan Province, there are 15 most commonly planted roadside tree species on the island. Table 2 shows the tree species and their properties. These tree species need to be further categorized to determine their fragility curves. Based on tree height, leaf area index, and leaf type, the Federal Emergency Management Agency (FEMA) classified the daily roadside tree species into 30 categories in total and each category corresponds to a fragility curve (FEMA 2003). The curves illustrate that hurricane wind has a probability to blow down trees when its velocity is higher than or equal to 50 knots (at Beaufort level 10). In this study, there are 14 categories corresponding to all tree species on Hainan Island. The fragility curves of the fourteen 14 on Hainan Island are sorted based on the value of probability of tree blowdown for Beaufort levels from 10 to 17. Table 2 also lists the rank of fragility curve of each tree species on Hainan Island. The less vulnerable a tree species is to blowdown, the higher rank the fragility curve has.

Table 2 Roadside tree species and their properties on Hainan Island

To verify the applicability of the FEMA fragility curves, we used a survey result on roadside tree losses of provincial highways induced by Typhoon Rammasun, which was provided by the TIC of Hainan Province. That survey summarized the number of blowdown trees of different species on different road segments that encountered different intensities of wind induced by Typhoon Rammasun. The rate of blowdown trees over different species encountering different intensities of wind can be calculated, as listed in Table 3, in comparison with probabilities from the FEMA fragility curves. From this table, the rate of blowdown trees by Typhoon Rammasun is basically consistent with the FEMA fragility curves.

Table 3 Comparison of historical rate of blowdown trees with probabilities from fragility curves

Figure 4a displays the fragility curves of tree species on Hainan Island, labelled by fragility curve rank. Figure 4b shows the spatial distribution of the fragility curve rank of roadside tree species. The most vulnerable roadside tree species are in the northeastern and southern parts of the road network, and the most robust tree species are concentrated in the center of the road network.

Fig. 4

Roadside tree species on Hainan Island. a Fragility curves of different species; b Spatial distribution of roadside tree species based on fragility curve rank

4.3 Economic Loss and Network Performance Assessment

Roadside tree blowdown induced by hurricane winds can cause economic loss. The procedure for evaluating the physical damage associated with economic loss of roadside tree blowdown includes: (1) For each cyclone event, determine the percentage of each road segment within the impacted area of each Beaufort wind level by GIS overlap analysis; (2) Based on fragility curves and the spatial distribution of Beaufort wind levels on road segments, estimate the expected total amount of trees that will collapse during each cyclone event. The spatial distribution of the expected number of damaged trees on the road network subject to a specific cyclone event is called the damage pattern of the cyclone; (3) The economic loss of roadside tree blowdown on each road segment is the sum of clean up costs and replanting fees for the damaged trees as shown in Eq. 5, assuming that replanted tree species are the same as those that collapsed.

Roadside tree blowdown causes traffic disruption on road segments, further diminishing the performance of the road network. Scientifically evaluating road network performance is the fundamental issue for calculating robustness and recovery efficiency. Travel demand and travel time are two important aspects that need to be considered when designing a performance metric for a road network. A damaged road network is still connected if there is at least one path connecting each origin and destination (OD) pair. Travel demand is assumed to be stable if the damaged road network is connected. Travel demand from disrupted OD pairs will be removed if the damaged road network is not connected. Intuitively, the damaged road network will perform better if it can serve more travel demand and have less travel time costs. Unfortunately, there are some numerical barriers to making an easy judgement. On the one hand, Bureau of Public Roads (BPR)’s link performance function is usually employed to estimate the travel time of each road segment, as shown in Eq. 15. This function means the travel time of each road segment increases exponentially as traffic flow increases. On the other hand, generally traffic flow assignment follows a user equilibrium (Patriksson 2015) in which all users pursue an individual optimum by choosing a path with the least travel time, so that each individual user cannot find a better path via unilateral change. Travel times based on the individual optimum will not be less than system optimum. The above two reasons imply that travel time increases much faster than travel demand. It is difficult to judge which road network performs better when comparing one road network with less satisfied travel demand and less travel time to another road network with more satisfied travel demand and more travel time. To better consider the impact of travel time and travel demand on the performance of a damaged road network, we propose a metric, expressed in Eq. 7. When a damaged road network cannot satisfy all the original travel demand, the value of the performance metric equals the percentage of the satisfied travel demand; otherwise, the performance metric also considers the ratio of the reciprocal of original travel time for the intact road network to the reciprocal of travel time for the damaged road network.

4.4 Optimization of Road Network Retrofit Investment

In this part, a mathematical model is first developed for describing the operation optimization for retrofitting roadside trees. Then a solution algorithm is designed to solve the problem.

4.4.1 Mathematical Model

The retrofit strategy we adopted is to replant the roadside tree species on each road segment with a species that have a higher fragility curve rank. An ideal retrofit scheme should largely alleviate the economic risk of roadside tree blowdown and enhance resilience of the road network after disasters. The proposed model is organized into three levels: upper level, middle level, and lower level.

At the upper level, both pre-disaster retrofit stage and the post-disaster recovery stage are considered. The prior stage determines the retrofit scheme for roadside trees of road segments; The later stage estimates economic loss of roadside tree blowdowns, road network robustness, and recovery efficiency under each possible tropical cyclone scenario. The economic risk of blowdown roadside trees is defined as the expected economic loss of roadside tree blowdown. The risk of low robustness can be measured by CVaR, the expectation of a specific percentile of worst robustness. Three objectives optimized at the upper level are: (1) minimize the economic risk of roadside tree blowdown; (2) maximize CVaR of robustness of the road network; and (3) maximize the expected recovery efficiency of the road network. Usually, the retrofit cost is constrained within a specific budget. The pareto-optimal solutions will be provided at the end of the estimation process, and the decision makers can choose the most suitable scheme according to their preferences.

The middle level determines the optimal repair sequence and calculates the recovery efficiency of the damaged road network under a damage pattern. Recovery efficiency is an integrated product of any post-disaster damage status and repair scheme. The optimal recovery efficiency is represented by the area under the curve of the recovery trajectory of the damaged road network, computed in Eq. 10.

The lower level model is a traffic assignment model, which follows the user equilibrium principle in which users minimize their own individual travel time during route choices. When disasters happen, some travel demands will be cancelled because of disconnectivity; some travel demands will be delayed due to detours replacing the initial paths. Therefore, total travel demand will decrease. Users will settle into a new user equilibrium.

Necessary assumptions are made in the model that: (1) drivers follow user equilibrium with full information of road status; (2) unsatisfied travel demand will be removed due to disconnectivity, and the travel demand of an OD pair is inelastic if it is still connected; (3) roadside trees on the same road segment will be repaired at the same time; and (4) only one road segment will be repaired at a time after a disaster.

The parameters and decision variables are defined as follows:

• (1) Parameters for the road network:

• N: set of nodes in the road network, i and j are its elements.

• A: set of arcs in the road network, arcs (i, j) are its element.

• L_(i,j): length of arc (i, j).

• (2) Parameters for roadside trees:

• Q: average distance between two spatially adjacent roadside trees.

• K: set of roadside tree species, k ∊ K.

• M(k): price of each roadside tree species k.

• Z(k): removal fee for roadside tree species k.

• Ψ(k): rank of fragility curve of roadside tree species k.

• T_(i,j): roadside tree species on arc (i,j).

• (3) Parameters for the retrofit process associated with different cyclone scenarios:

• S: set of cyclone events, s ∊ S; Υ number of cyclone events.

• TT_s: total repair time spent on road network under damage pattern s.

• IT_ɛ: repair time spent on damaged arc ɛ under damage pattern s.

• V: set of all strategies $V_{{\left( {i,j} \right)}} = \{ v|v \in \left\{ {0, 1, \ldots , q} \right\} \cap \varPsi \left( {\varXi \left( v \right)} \right) \ge \varPsi \left( {T_{{\left( {i,j} \right)}} } \right)\}$, q is the total number of strategies, Ξ(v) is the replanted roadside tree species when adopting strategy v.

• U: set of Beaufort wind levels from 10 to 17, u is its element.

• π ^u_(i,j)s : percentage of length of arc (i, j) within the impact area of Beaufort wind level u under cyclone event s.

• τ_(i,j)s(π ^u_(i,j)s ,$\varXi \left( {V_{{\left( {i,j} \right)}} } \right)$): probability of roadside tree blowdown on arc (i, j) after adopting strategy V_(i,j) under damage pattern s.

• (4) Decision variables for the retrofit process associated with different cyclone scenarios:

• RE_s: recovery efficiency of the road network under damage pattern of cyclone event s (we call it damage pattern s below).

• PD_s: the economic loss of roadside tree blowdown under damage pattern s.

• Ru_s: robustness of road network under damage pattern s; x is the variable of robustness of the road network under all damage patterns.

• $\Delta_{s}$: optimal repair sequence of damaged arcs under damage pattern s, ɛ is its index; E_s is the number of damaged arcs under damage pattern s; $\Delta_{s}$ is the decision variable at the middle level of the model.

• Ω: set of intact arcs on the current road network. If the road network is not damaged, Ω = Ω₀ = A; if the road network is damaged and not repaired yet, $\varOmega = \varOmega_{s} = A/\left\{ {\Delta_{s} } \right\}$; if the damaged road network has been repaired r damaged arcs, $\varOmega = \varOmega_{s}^{r} = A/\left\{ {\Delta_{s} /\left\{ {\Delta_{s}^{1} , \Delta_{s}^{2} , \ldots ,\Delta_{s}^{r} } \right\}} \right\}$. We use Ω to represent the status of the road network.

• P_Ω: the performance of road network $\varOmega$.

• (5) Parameters for traffic flow in the road network:

• W_Ω: set of origin-destination (OD) pairs that can be connected by at least one path in the road network Ω, w is its element.

• O ^Ω_w : origin of OD pair W_Ω.

• D ^Ω_w : destination of OD pair W_Ω.

• Y ^Ω_w : travel demand of OD pair W_Ω.

• t ⁰_(i,j) : free flow travel time of arc (i, j).

• c_(i,j): capacity of arc (i, j).

• (6) Decision variables for traffic flow in the road network:

• ϱ(Ω): total travel demand of OD pairs.

• f ^Ω_(i,j) : traffic flow of arc (i, j).

• t ^Ω_(i,j) : travel time that users spend passing through arc (i, j).

• X ^Ω_(i,j)w : traffic flow of arc (i, j) for OD pair w in the road network Ω; a decision variable at the lower level of the model.

• λ(Ω): total travel time of the road network Ω, including three types of scenarios Ω₀, Ω_s, and Ω ^r_s .

• Z_Ω: objective value of user-equilibrium based on traffic assignment in the road network Ω.

The proposed model is constituted by three levels: upper level, middle level, and lower level. The upper level deals with optimal retrofit decision making; the middle level solves the optimal repair sequence problem; and the lower level reflects traffic demand on a road network.

(1) The upper level:

$${\text{Min}}\frac{1}{\Upsilon} \times \mathop \sum \limits_{s \in S} PD_{s}$$

(1)

$${\text{Max}}\quad CVaR_{\alpha } \left( x \right) = E[x|x < \hbox{max} \{ z|F_{x} \left( z \right) \le \alpha \} ],\quad x \in \left\{ {Ru_{s} ,\forall s \in S} \right\}$$

(2)

$${\text{Max}}\quad \frac{1}{\Upsilon} \times \mathop \sum \limits_{s \in S} RE_{s}$$

(3)

Subject to:

$$\mathop \sum \limits_{{\left( {i,j} \right) \in A}} \text{sgn} \left( {V_{{\left( {i,j} \right)}} } \right) \times M\left( {\varXi \left( {V_{{\left( {i,j} \right)}} } \right)} \right) \times \frac{{L_{{\left( {i,j} \right)}} }}{Q} \le B$$

(4)

$$PD_{s} = \mathop \sum \limits_{{\left( {i,j} \right) \in A}} \mathop \sum \limits_{u = 10}^{17} \tau_{{\left( {i,j} \right)s}} \left( {\pi_{{\left( {i,j} \right)s}}^{u} , \varXi \left( {V_{{\left( {i,j} \right)}} } \right)} \right) \times \left( {M\left( {\varXi \left( {V_{{\left( {i,j} \right)}} } \right)} \right) + Z\left( {\varXi \left( {V_{{\left( {i,j} \right)}} } \right)} \right)} \right) \times \frac{{L_{{\left( {i,j} \right)}} }}{Q}\quad \forall s \in S$$

(5)

$$Ru_{s} = P\left( {\varOmega_{s} } \right),\quad \forall s \in S$$

(6)

$$P_{\varOmega } = \left\{ {\begin{array}{*{20}l} {\frac{1}{2} \times \frac{\varrho \left( \varOmega \right)}{{\varrho \left( {\varOmega_{0} } \right)}}, } \hfill & if \ \varrho \left( \varOmega \right) < \varrho \left( {\varOmega_{0} } \right) \hfill \\ {\frac{1}{2} + \frac{1}{2} \times \frac{{\lambda \left( {\varOmega_{0} } \right)}}{\lambda \left( \varOmega \right)},} \hfill & if \ \varrho \left( \varOmega \right) = \varrho \left( {\varOmega_{0} } \right) \hfill \\ \end{array} } \right.,\quad \varOmega \in \left\{ {\varOmega_{0} , \varOmega_{s} , \varOmega_{s}^{r} } \right\},\quad \forall s \in S$$

(7)

$$\lambda \left( \varOmega \right) = \mathop \sum \limits_{{\left( {i,j} \right) \in A}} f_{{\left( {i,j} \right)}}^{\varOmega } \times t_{{\left( {i,j} \right)}}^{\varOmega } \left( {f_{{\left( {i,j} \right)}}^{\varOmega } } \right),\quad \varOmega \in \left\{ {\varOmega_{0} , \varOmega_{s} , \varOmega_{s}^{r} } \right\},\quad \forall s \in S$$

(8)

$$\varrho \left( \varOmega \right) = \mathop \sum \limits_{w \in W} Y_{w}^{\varOmega } ,\quad \varOmega \in \left\{ {\varOmega_{0} , \varOmega_{s} , \varOmega_{s}^{r} } \right\},\quad \forall s \in S$$

(9)

where α is the significant level of the conditional value at risk, which is set to be 0.05. Equation 1 is the objective function for minimizing the economic risk of roadside tree blowdown. Equation 2 is the objective function for maximizing the conditional value at risk of robustness, F_x is the cumulative distribution function of robustness, and z is the percentile of F_x that is not less than α. Equation 3 is the objective function for maximizing the expected recovery efficiency. Equation 4 is the retrofit cost based on the chosen retrofit scheme. The retrofit cost should be less than the preset budget. Equation 5 is the formulation of the economic loss of roadside tree blowdown under damage pattern s. Equation 6 is the formulation of robustness of the road network under damage pattern s. Equations 7 and 8 are respectively the formulations of the performance and total travel time of road network $\varOmega$. Equation 9 shows the formulation of satisfied travel demand of OD pairs under damage pattern s.

(2) The middle level:

$${\text{Max}}\quad RE_{s} = \frac{1}{{TT_{s} }} \times \mathop \sum \limits_{\varepsilon = 0}^{{E_{s}}} P\left( {\varOmega_{s}^{\varepsilon } } \right) \times IT_{\varepsilon } ,\quad \forall \varepsilon \in \Delta_{s} ,\quad \forall s \in S$$

(10)

Subject to

$$TT_{s} = \mathop \sum \limits_{\varepsilon = 1}^{{E_{s} }} IT_{\varepsilon } ,\quad \forall \varepsilon \in \Delta_{s} ,\quad \forall s \in S$$

(11)

Equations (7)–(9)

where Eqs. 10 and 11 are respectively the formulations of recovery efficiency and total repair time under damage pattern s.

(3) The lower level:

$$\hbox{min} Z_{\varOmega } = \mathop \sum \limits_{{\left( {i,j} \right) \in A}}^{{}} \mathop \int \limits_{0}^{{f_{{\left( {i,j} \right)}}^{\varOmega } }} t_{{\left( {i,j} \right)}}^{\varOmega } \left( v \right)dv$$

(12)

Subject to

$$\mathop \sum \limits_{j \in N}^{{}} X_{{\left( {i,j} \right)w}}^{\varOmega } - \mathop \sum \limits_{j \in N}^{{}} X_{{\left( {j,i} \right)w}}^{\varOmega } = \left\{ {\begin{array}{*{20}l} {Y_{w}^{\varOmega } , \quad i = O_{w}^{\varOmega } } \hfill \\ { 0, \quad i \in N/\left\{ {O_{w}^{\varOmega } , D_{w}^{\varOmega } } \right\}} \hfill \\ { - Y_{w}^{\varOmega } , \quad i = D_{w}^{\varOmega } } \hfill \\ \end{array} } \right.,\quad \forall w \in W_{\varOmega } ,\quad \forall i \in N$$

(13)

$$f_{{\left( {i,j} \right)}}^{\varOmega } = \mathop \sum \limits_{{w \in W_{\varOmega } }}^{{}} X_{{\left( {i,j} \right)w}}^{\varOmega } ,\quad \forall \left( {i,j} \right) \in A$$

(14)

$$t_{{\left( {i,j} \right)}}^{\varOmega } \left( v \right) = t_{{\left( {i,j} \right)}}^{0} \times \left( {1 + \varpi \left( {\frac{{f_{{\left( {i,j} \right)}}^{\varOmega } }}{{C_{{\left( {i,j} \right)}} }}} \right)^{\beta } } \right),\quad \forall \left( {i,j} \right) \in A$$

(15)

$$X_{{\left( {i,j} \right)w}}^{\varOmega } \ge 0,\quad \forall \left( {i,j} \right) \in A,\quad \forall w \in W_{\varOmega }$$

(16)

where the constant parameter ϖ is 0.15 and β is 4. The lower level of the model was used in three scenarios: (1) intact road network; (2) damaged road network; and (3) partially repaired road network. Only the travel demand between the connected OD pairs is considered in the model. With reference to the traffic assignment models provided by Patriksson (2015), and with the premise of satisfying symmetry of the cost function, the objective function of user equilibrium, minimizing the individual user’s travel time, can be converted into Eq. 12. Equation 13 calculates the traffic flow balance. Equation 14 means that the traffic flow of arc (i,j) is the sum of any OD pairs’ traffic flow that pass through arc (i,j). Equation 15 is the BPR’s link performance function. Equation 16 states that the traffic flow of any arc should be nonnegative.

4.4.2 Solution Approach

The upper level problem is a multiobjective, nonlinear optimization problem. The middle level problem is an optimal sequential problem that has been proved to be a NP-hard problem. The lower level problem is equivalent to a user-equilibrium-based traffic assignment model, which is a convex optimization problem. To solve this trilevel optimization problem, we developed a trilevel algorithm as shown in Fig. 5. At the upper level of the problem, the NSGA II algorithmic framework was used to identify a set of nondominated optimal retrofit schemes. Predisaster retrofit schemes can influence the post-disaster damage patterns of roadside trees on road segments. At the middle level, under a specific damage pattern, a GA was implemented to find the optimal repair sequence of damaged road segments. At the lower level, the Frank-Wolfe algorithm was adopted to solve the traffic assignment model for a given road network (intact/damaged/being repaired).

Fig. 5

Structure of the trilevel algorithm and the single-level algorithm

The computation time of the trilevel algorithm is huge, particularly when the road network or the number of damaged patterns is large. As shown in Fig. 5, the output of the GA at the middle level of the algorithm is the recovery efficiency, which is an input factor of fitness functions of NSGA II at the upper level of the algorithm. Solving the middle level problem under each damage pattern costs considerable running time. Outputs of the lower level of the algorithm are traffic flow and travel time of each road segment that can be used to further calculate road network performance. Road network performance is an input factor of the recovery efficiency fitness function of the GA at the middle level and a conditional value at the risk of robustness fitness function of NSGA II at the upper level. Frequently using the lower level of the algorithm also consumes considerable running time. To lighten the computational burden, the random forest method was used to directly find the relationship between the damage patterns and the recovery efficiency of the road network. Therefore, the trilevel algorithm can be transformed into a single-level algorithm, saving a great deal of running time.

The implementation process of the random forest method is described in the literature (Hastie et al 2009). Essential user design elements of random forest include: (1) Independent variables are the total numbers of damaged roadside trees on road segments. One dependent variable is recovery efficiency. Table 4 shows the compiled disaster events dataset including damage patterns and recovery efficiency; (2) The random forest method for regression is used instead of classification due to continuous independent variables; (3) To avoid overfitting, purity improvement was used when pruning unnecessary nodes.

Table 4 Disaster events dataset including damage patterns, robustness, and recovery efficiency

5 Result Analysis

Five main implications are found in this study. The effectiveness of road network performance metric is validated via simulation. Random forest method is successfully leveraged for estimating recovery efficiency. The relationship between robustness and recovery efficiency is examined by simulation results. Optimal solutions obtained by the proposed solution algorithm are analyzed. The marginal effect of a pre-set budget on retrofit schemes is also investigated.

5.1 Road Network Performance Metric Analysis

To verify the practicability of the proposed road network performance metric, we conducted a simulation to explore the relationship between satisfied travel demand and travel time. The number of realizations is Y, and was set at 30,000. Road segments are assumed to have the same failure probability. The failure probability at the yth realization is f = y/Y. For each road segment, a random number rnd (a fraction between 0 and 1) was generated to compare with the failure probability. If rnd was larger than f, the road segment would be dysfunctional. Ultimately, different damaged road networks are generated. This simulation makes sure that all possible damaged road network segments can be evenly sampled. Based on the traffic assignment model, we calculated the satisfied travel demand and travel time for each damaged road network segment, as shown in Fig. 6a. It is implied that the damaged road network with less satisfied travel demand also probably has less travel time. This is because when a large portion of road segments are damaged, disconnectivity will take place within a network, which reduces the satisfied travel demand, as well as the traffic time in some road segments due to less congestion. Figure 6b shows the performance value of all damaged road networks in descending order. The discontinuous curve separates the damaged networks into two groups: one is the networks in which satisfied travel demand equals original travel demand and another is the networks with satisfied travel demand that is less than original travel demand.

Fig. 6

Simulation result analysis. a Relationship between travel time (unit: hours) and satisfied travel demand. b Performance values of all damaged road network in descending order. c Correlation of status of road segments. d Relationship between robustness and recovery efficiency

5.2 Recovery Efficiency Estimation Using Random Forest Method

To prepare the dataset used by the random forest method for predicting recovery efficiency, we designed a simulation. The overall retrofitting percentage of road segments (equivalent to preset budget) was set from 0.1 to 0.9 with intervals of 0.1. The retrofit strategy for each road segment was randomly determined. The above settings make sure that all possible retrofitting road networks can be evenly sampled. For each retrofitted road network, based on the impact area of each Beaufort wind level during each cyclone event and the roadside tree fragility curve of each road segment, the damage pattern of the road network for each cyclone event can be generated. If there exists at least one damaged roadside tree on a road segment, the road segment will be dysfunctional and will need to be repaired. Based on different damage patterns, we can search the optimal repair sequence and calculate the corresponding recovery efficiency of the damaged road networks. The number of realizations in the simulation is 30,000. The damage patterns and calculated recovery efficiency constitute the dataset for the random forest procedure.

We used 70% of the dataset as a training dataset and 30% of the dataset as a test dataset. We split the training dataset into 10 pieces, using 1 piece as a validation dataset and another 9 pieces as a training dataset and implemented cross-validation for 10 times and calculated the average validation error, using the obtained training model to predict the test dataset and estimate the test error. Random forest (RF) was coded using the Scikit-learn package (Pedregosa et al. 2011) in PyCharm (2019.2 community edition).² To validate the effectiveness of random forest, we compared random forest with several common machine learning methods including Decision Tree, Bagging, AdaBoosting, Gradient Boosting, Linear Regression, Partial Least Square, Ridge, Lasso, and Support Vector Machine (SVM). Table 5 presents the errors of the cross-validation and test errors of different methods measured by the mean square error, which shows that random forest has the least mean square error (MSE) in the prediction of recovery efficiency.

Table 5 Comparison among different machine learning models

Figure 6c shows the correlation between the damaged status of any two road segments. In general, the correlations among the status of road segments are low, implying there is no obvious multicollinearity. That is why linear regression and robust regression models (Partial Least Square, Ridge, Lasso, and SVM) do not perform well and have a lower MSE. Random forest has the best learning effect in this case because it is an ensemble learner that consists of many decision trees. Compared with other ensemble learners (Bagging, AdaBoosting, and Gradient Boosting), RF generates each decision tree by choosing a different random subset of features instead of using all features, which decorrelates the decision trees. Therefore, RF has less variance and less bias.

5.3 The Relationship Between Robustness and Recovery Efficiency

Robustness of the damaged road network under each damage pattern can be calculated within the previous simulation. Figure 6d illustrates the linear relationship between robustness and recovery efficiency at a statistical level. While at an operational level, the damaged road network with a larger value of recovery efficiency may have a lower value of robustness, such as Point B compared with Point A, because the topology of the damaged road network has an influence on the optimal repair sequence and further affects recovery efficiency. Low robustness of damaged road networks can cause high deterioration of road network performance. If the repair process cannot be launched immediately, adverse social and economic consequence will be induced. Therefore, considering both robustness and recovery efficiency has practical meaning.

5.4 The Optimal Retrofit Scheme

A set of retrofitting strategies for replanting tree species along damaged portions of the road network on Hainan Island was determined based on Table 2 according to two criteria: lowest cost (the sum of plants’ prices and cleaning up fees) and lowest fragility curve ranking. Each of these replanting tree species is the best choice by both criteria. There are four potential retrofit strategies displayed in Table 6 and one strategy involves no action.

Table 6 Corresponding tree species of the set of retrofit strategies for the road network on Hainan Island

NSGA II was also coded and implemented using PyCharm (2019.2 community edition). The experimental environment was an Intel Core i7 8565U Processor PC with 16 GB memory. Figure 7a illustrates the Pareto-optimal solutions found by NSGA II. We can see that the solution having the least expected economic loss of roadside tree blowdown also has the least expected recovery efficiency. The 1st solution also has a lower CVaR in terms of robustness than the 2nd solution and the 2nd, 3rd, and 4th solutions have equal CVaR in terms of robustness. Four different optimal nondominated solutions correspond to four different retrofit schemes, which are the result of the spatial heterogeneity in road network. Spatial heterogeneity in a road network is reflected by spatial heterogeneity of hazard intensity (wind speed), spatial heterogeneity of fragility curves of roadside trees, and spatial heterogeneity of topological importance of different road segments. For example, Fig. 7b, c shows the spatially distributed retrofit schemes of the 1st and 4th solutions. The 1st solution retrofits road segments that have higher probability of suffering from higher hazard intensity (for example, Road Segment A in Fig. 7b) or whose roadside trees are more fragile (for example, Road Segment B in Fig. 7b). The 4th solution retrofits road segments whose topological positions are more important (for example, Road Segment C in Fig. 7c).

Fig. 7

Pareto-optimal solutions. a Four optimal nondominated solutions. b, c respectively are the 1st and 4th solutions of optimal nondominated solutions

5.5 Comparison Analysis

To study the influence of the marginal effect of a pre-set budget on optimal retrofit schemes, we set budget options to RMB 25,000,000, 50,000,000, and 75,000,000 yuan, respectively. Figure 8a illustrates the optimal nondominated solutions for retrofit schemes subject to three pre-set budgets. Figure 8b–d show the mean values of expected economic loss of roadside tree blowdown, CVaR in terms of robustness and the expected recovery efficiency of optimal nondominated solutions respectively with three pre-set budgets and the option without retrofitting investment. Retrofitting investment largely alleviates the expected economic loss of roadside tree blowdown and enhances robustness and expected recovery efficiency. The expected economic loss of roadside tree blowdown with the RMB 25,000,000 yuan budgeted investment option is only 38.37% of that of the option without retrofitting investment, while its CVaR in terms of robustness and expected recovery efficiency are 1.74 and 14.5 times of those of the option without retrofitting investment, respectively.

Fig. 8

Retrofit schemes comparison. ) The optimal nondominated solutions for retrofit schemes subject to three pre-set budgets. b–d are the mean values of expected economic loss of roadside tree blowdown, CVaR in terms of robustness, the expected recovery efficiency of optimal nondominated solutions respectively with three pre-set budgets

Figure 8b–d show that by increasing budget, the expected economic loss of roadside tree blowdown will decrease while both CVaR in terms of robustness and expected recovery efficiency will increase. It is implied that the resilience of optimally retrofitted road network will be enhanced with a higher budget. The expected economic loss of roadside tree blowdown decreases by 79.18% when the budget increases from RMB 25,000,000 to 50,000,000 yuan, and it decreases 39.81% when the budget increases from RMB 50,000,000 to 75,000,000 yuan. CVaR in terms of robustness and expected recovery efficiency increases by 124.31% and 8.06% when the budget increases from RMB 25,000,000 to 50,000,000 yuan, respectively, and it increases by 19.08% and 0.81% when the budget increases from RMB 50,000,000 to 75,000,000 yuan, respectively. Overall, the marginal effect between RMB 25,000,000 and 50,000,000 yuan is much larger than that between RMB 50,000,000 and 75,000,000 yuan.

6 Conclusion

This article is the first study on optimal retrofit investment in road networks subject to roadside tree blowdown induced by strong winds from tropical cyclones. A trilevel, two-stage, and multiobjective stochastic mathematical model was developed. The upper level deals with optimal retrofit decision making, the middle level solves the optimal repair sequence problem, and the lower level reflects traffic demand on the road network. The proposed model was solved by a trilevel algorithm that can be reduced to a single-level algorithm by random forest, decreasing a great amount of running time.

The road network on Hainan Island was chosen as a study area. Our analyses show that: (1) The proposed performance metric of a damaged road network is better than other existing metrics, since it separates damaged networks into two different groups: one with original travel demand being satisfied and another with original travel demand not being all satisfied; (2) Compared with other machine learning models, the Random Forest method presents the best performance with least mean square error in both validation and test processes when predicting the robustness and recovery efficiency of damaged road networks; (3) Consideration of both robustness and recovery efficiency has realistic meaning because the damaged road network, with a larger value of recovery efficiency, may have a lower value of robustness at operational level; (4) Four optimal nondominated solutions were found as a result of the spatial heterogeneity of road networks; (5) Retrofitting investment significantly diminishes the expected economic loss of roadside tree blowdown and enhances robustness and expected recovery efficiency. Increased investment budgets produce marginal effects for the three objective function values.

The proposed mathematic model can also be adapted to solve the optimal retrofit problem of other infrastructure network systems. The ideas of coupling evolutionary algorithms and machine learning models provides insights into how to reduce the computational burden in heuristic optimization, which can be adopted by many other optimization problems.