Two-stage Robust Optimization Approach for Enhanced Community Resilience Under Tornado Hazards

Mehdi Ansari
American Airlines
Fort Worth, TX
[email protected] &Juan S. Borrero
School of Industrial Engineering and Management
Oklahoma State University
[email protected] &Andrés D. González
School of Industrial and Systems Engineering
University of Oklahoma
[email protected]

Abstract

Catastrophic tornadoes cause severe damage and are a threat to human wellbeing, making it critical to determine mitigation strategies to reduce their impact. One such strategy, following recent research, is to retrofit existing structures. To this end, in this article we propose a model that considers a decision-maker (a government agency or a public-private consortium) who seeks to allocate resources to retrofit and recover wood-frame residential structures, to minimize the population dislocation due to an uncertain tornado. In the first stage the decision-maker selects the retrofitting strategies, and in the second stage the recovery decisions are made after observing the tornado. As tornado paths cannot be forecast reliably, we take a worst-case approach to uncertainty where paths are modeled as arbitrary line segments on the plane. Under the assumption that an area is affected if it is sufficiently close to the tornado path, the problem is framed as a two-stage robust optimization problem with a mixed-integer non-linear uncertainty set. We solve this problem by using a decomposition column-and-constraint generation algorithm that solves a two-level integer problem at each iteration. This problem, in turn, is solved by a decomposition branch-and-cut method that exploits the geometry of the uncertainty set. To illustrate the model’s applicability, we present a case study based on Joplin, Missouri. Our results show that there can be up to 20% reductions in worst-case population dislocation by investing $15 million in retrofitting and recovery; that our approach outperforms other retrofitting policies, and that the model is not over-conservative.

Keywords Robust Optimization $\cdot$ Column and constraint generation $\cdot$ Tornado $\cdot$ Community Resilience

1 Introduction

Tornadoes are vertical columns of rotating air that are spawned from supercell thunderstorms caused by rotating updrafts that form because of the shear in the environmental wind field [NOAA, 2022d]. Catastrophic tornadoes are common natural disasters that happen in many populated regions around the world and are particularly concerning in dense high-risk urban areas because of their damage potential. On average, more than 1200 tornadoes happen in the US annually which cause around 60 fatalities [NOAA, 2022a]. Tornadoes have caused damage costs that range between $183 million to $9.493 billion per year in the US within the past two decades [NOAA, 2022e]. The intensity of tornadoes is measured by the Enhanced Fujita (EF) scale, which ranks a tornado from 0 (weakest) through 5 (most violent) based on an estimation of their average wind speed. The average wind speed, in turn, is estimated by comparing the observed tornado damages with a list of Damage Indicators and Degrees of Damage [NOAA, 2022f].

Fatalities and injuries from tornadoes have decreased in the past few decades thanks to warning systems based on weather forecasts [Standohar-Alfano et al., 2017, Koliou and van de Lindt, 2020]. For example, the Storm Prediction Center at the National Weather Service keeps a Day 1-8 Convective Outlook for all the US and issues watches or warnings based on the possibility or observation of tornadoes [NOAA, 2022b]. Tornado warnings are issued as soon as a rotating supercell has been identified by a weather radar; these warnings instruct people in the affected area to take shelter immediately. Even though the warning systems have caused a great reduction in fatalities [FEMA, 2022a], their short lead times (of the order of minutes) do not allow to prevent physical structures from possible destruction. This is further aggravated by the fact that more than 80% of building stocks in the US are wood-frame buildings that are highly vulnerable to wind damage [van de Lindt and Dao, 2009].

Besides weather forecasting, other alternatives can be employed to reduce the impact of tornadoes. In the past decade, several studies have shown that retrofitting strategies with simple and inexpensive actions, e.g., improving the roof cover or enhancing the roof sheathing nailing pattern of a house, can improve building codes to make wood-frame buildings more resistant to tornado damage, particularly from damaging tornadoes of EF2 intensity or less [Simmons et al., 2015, Ripberger et al., 2018, Masoomi et al., 2018, Koliou and van de Lindt, 2020, Wang et al., 2021]. In particular, Ripberger et al. [2018] shows that a 30% or more reduction in lifetime damage can be expected by enhancing existing wood-frame buildings with simple retrofitting actions. Similarly, recent research has studied different recovery strategies and has shown that they have the potential to significantly improve the restoration time and reduce population dislocation [Masoomi and van de Lindt, 2018, Farokhnia et al., 2020, Koliou and van de Lindt, 2020].

At the federal level in the US, the use of retrofitting has been identified as a valuable strategy to reduce the impact of natural disasters [FEMA, 2021, 2022b, The White House, 2022]. At the local level, however, such strategies have not been implemented yet. In fact, there are federal programs that provide general guidelines for local administrators (e.g., the emergency management agencies of cities and towns) to design and fund retrofitting plans [FEMA, 2021]. Motivated by this need, the main goal of this paper is to present a model to guide retrofitting and recovery investments in order to enhance the community resilience of an urban area subject to tornado hazards. Here, community resilience is understood as “the ability to prepare for anticipated hazards, adapt to changing conditions, and withstand and recover rapidly from disruptions” [McAllister et al., 2015].

Specifically, we assume that there is a decision-maker (a government agency or a public-private consortium; for instance the emergency management department of a city and home insurance companies) that can invest a limited budget to retrofit and recover structures across a geographical area of interest, in order to maximize the community resilience due to an uncertain tornado. We model this decision problem as a two-stage robust optimization model. The first-stage decisions are made before the realization of a tornado and determine what structures have to be retrofitted and to what extent. The second-stage decisions determine recovery strategies for the locations that are damaged by the tornado. Without loss of generality, in this article we consider total population dislocation, which is the involuntary movement of people from their residential sites after a severe tornado, as a metric for community resilience to be minimized (therefore, maximizing community resilience would be equivalent to minimizing population dislocation). Our proposed formulation, however, allows for other metrics that assess community resilience.

We formulate the first- and second- stage problems as integer programming problems (IP) whose decision variables select retrofitting strategies and recovery strategies, respectively, for each location. Both stages share the same budget, which models that, at least a priori, the decision-maker does not know what proportion of the budget should be spent on retrofitting and what proportion on recovery. On the other hand, as an accurate prediction of the location/time where a tornado forms, the motion direction, and its magnitude, is out of reach of the current technology [NOAA, 2022c], we approach the uncertainty using a robust lens. That is, we assume that given any retrofitting plan, the corresponding worst-case tornado damage happens. In order to have an adequate balance between accuracy and tractability, we model tornado paths as arbitrary line segments on the plane and assume that an area is affected by the tornado if it is sufficiently close to the line segment. This modeling approach results in a mixed-integer non-linear uncertainty set.

In order to solve the problem, we employ an exact column-and-constraint generation (C&CG) algorithm based on the method discussed by Zeng and Zhao [2013]. Here, at each iteration, a master problem solves a relaxation of the original problem over a subset of possible tornado scenarios. New scenarios are iteratively generated in a subproblem and added to the master until the objective value does not change. In contrast with previous models that employ the C&CG algorithm, for instance Zeng and Zhao [2013], Yuan et al. [2016], Velasquez et al. [2020], our subproblem is a challenging max-min non-linear integer problem that cannot be solved by standard ‘dualize and combine’ techniques. Thus, the subproblem is solved by a decomposition branch-and-cut algorithm that exploits the geometric properties of the uncertainty set. Particularly, we derive ‘double’ and ‘triple’ conflict constraints to define the master relaxation and we check the feasibility of master solutions by using a stabbing line algorithm and a continuous non-convex optimization problem.

Using geographical and population data from IN-CORE [2022], we use the proposed model to determine optimal retrofitting and recovery actions in Joplin, MO. The results of the case study show that retrofitting actions should be performed in most locations in the geographical center of the city, which are also the locations with higher population density. They also show that spending budget in retrofitting is prioritized over recovery, and that only when a ‘critical’ set of locations in the center of the city are retrofitted, there should be expenditures in recovery. By performing simulations we also show that the proposed model outperforms other retrofitting policies, and that the model is not over-conservative: the optimal worst-case population dislocation can be within 10% of the maximum population dislocation observed in the simulations.

To summarize, in this paper we make the following contributions:

•

We propose a novel two-stage optimization model under uncertainty to aid decision-makers in the allocation of resources to retrofit and recover residential wood-frame buildings in tornado-prone regions.
•

We explicitly model tornado paths as arbitrary line segments on the plane and formulate the problem as a two-stage robust optimization problem with integrality requirements in the first and second stages and in the uncertainty set.
•

We develop an exact algorithm that embeds a decomposition branch-and-cut algorithm within a column-and-constraint generation method. By exploiting the geometric properties of line segments, we develop initialization and separation procedures to effectively implement the embedded decomposition branch-and-cut algorithm.
•

Using real data we provide optimal retrofitting and recovery strategies for Joplin, MO. The results show that there can be up to 20% reductions in worst-case population dislocation by investing $15 million; that our approach outperforms other retrofitting policies, and that the model does not suffer from over-conservativeness.

The remainder of this paper is organized as follows. Section 2 reviews the literature. Section 3 describes the two-stage robust optimization problem. Section 4 presents a customized version of the C&CG algorithm. We describe the decomposition branch-and-cut method to solve the subproblem in C&CG algorithm in Section 5. In Section 6, we conduct the numerical experiments on Joplin data as our case study. Lastly, Section 7 concludes the paper.

2 Literature review

2.1 Optimization models to retrofit and recover under tornado hazards

The determination of retrofitting and recovery strategies for tornado hazards, particularly using optimization, is a relatively new topic in the literature. Wen [2021] proposes a multi-objective optimization model to retrofit against tornado hazards. Their model, however, is single-stage, does not considers recovery, does not consider uncertain tornado paths, and aggregates all of the uncertainty of the problem into the parameters of the deterministic model. Simulation, on the other hand, is more commonly used to study the behavior of tornadoes, particularly to evaluate the impact of retrofitting strategies, see for example Strader et al. [2016], Wang et al. [2017], Masoomi and van de Lindt [2018], Fan and Pang [2019], Wang et al. [2021], Stoner and Pang [2021]. These works, however, do not explicitly consider the decision problem of allocating resources to retrofit and recover.

The closest models in the literature to the present work are two-stage robust optimization models for (general) disaster planning, see Yuan et al. [2016], Ma et al. [2018], Matthews et al. [2019], Velasquez et al. [2020], Cheng et al. [2021]. These models, however, are geared towards supply and distribution problems in networks and do not capture the specific patters of the disasters (such as, e.g., paths or ripple effects). Consequently, these models cannot be adapted to deal with tornadoes and the specific retrofitting and recovery setting we study here.

2.2 Robust optimization

Robust optimization methods have been extensively developed in the past two decades to mathematically model and study different systems under uncertainty. Standard robust optimization models assume that the uncertain parameters belong to a convex set [Ben-Tal and Nemirovski, 1999, Bertsimas and Sim, 2004], which have many applications such as scheduling [Lin et al., 2004], supply chain [Ben-Tal et al., 2011], transportation [Yao et al., 2009], power system problems [Xiong et al., 2017], among others. Mixed-integer version of uncertainty sets have also been studied. Exact methods to solve a general robust optimization problem with a mixed-integer uncertainty set are similar to the integer L-shaped method [Laporte and Louveaux, 1993] in which a master problem iteratively solves a restricted version of the original problem under a few uncertain scenarios and for the solution vector. Here, a subproblem generates a new scenario by finding the worst-case of objective within a mixed-integer uncertainty set [Mutapcic and Boyd, 2009, Ben-Tal et al., 2015, Ho-Nguyen and Kılınç-Karzan, 2018, Borrero and Lozano, 2021, Ansari et al., 2022].

Ben-Tal et al. [2004] extend the scope of standard single-stage robust optimization problems by introducing the adjustable robust optimization methodology, where a part of decision variables must be determined before the realization of the uncertainty set. The rest of the decision variables are chosen when the uncertain parameters are revealed. The two-stage approach provides a less conservative framework to deal with the uncertainty and it can model a broad range of applications in different areas such as transportation [Gabrel et al., 2014], networks [Atamtürk and Zhang, 2007], investment [Takeda et al., 2008], power systems problems [Shams et al., 2021, Li et al., 2021], among others.

In general, two-stage robust optimization problems are computationally expensive. To efficiently solve them, decomposition methods are widely employed, under the assumption that the second-stage problem is a linear programming problem [Thiele et al., 2009, Zhao and Zeng, 2012, Jiang et al., 2012, Gabrel et al., 2014]. Zeng and Zhao [2013] presents the column and constraint generation method (C&CG) as an alternative to address this class of robust problems. They showed the C&CG can outperform previous decomposition for many problem settings. The C&CG approach has become a popular tool to solve two-stage robust optimization in the past decade [An et al., 2014, An and Zeng, 2015, Jabr et al., 2015, Neyshabouri and Berg, 2017, Ding et al., 2017, Yuan et al., 2016, Matthews et al., 2019, Velasquez et al., 2020, Cheng et al., 2021]. An important feature of existing applications of the C&CG method is that the second-stage optimization problem is convex. In contrast, in our problem the second-stage problem is an IP problem and the uncertainty set is mixed-integer non-linear. Therefore, the standard ‘dualize and combine’ approach that is used in the literature to solve the subproblem in the C&CG does not apply to our case.

3 Model formulation

Next, we define the two-stage robust optimization model to minimize the total population dislocation under an uncertain tornado and present a mathematical formulation of the problem.

3.1 Two-stage robust optimization formulation

Recall that the first stage decisions seek to determine what locations to retrofit before the realization of a tornado. The second stage decisions, that happen after the uncertainty is revealed, select what recovery strategies should be implemented in the locations that are affected by the tornado.

Formally, suppose $S$ denotes the set of retrofitting strategies and $P$ denotes the set of recovery plans. We consider a set of locations of interest $L$ , hereafter referred to as locations for simplicity, which are placed in a 2-dimensional plane. Let $w_{\ell s}$ be the population dislocation estimation pre-tornado at location $\ell\in L$ under retrofitting strategy $s\in S$ (in most cases $w_{\ell s}=0$ ) and let $g_{\ell sp}$ be the population dislocation post-tornado at location $\ell\in L$ assuming that the retrofitting strategy $s\in S$ is used and that the recovery plan is $p\in P$ .

We assume that the decision-maker has a limited budget of $A\geq 0$ to invest in the retrofitting and recovery plans. This assumption allows the model to allocate the budget endogenously, rather than having the decision-maker setting the proportions before hand. If required, the model is able to handle fixed budgets by adding a linear constraint in the first-stage problem. Let $d_{\ell s}$ be the retrofitting cost of all the buildings in location $\ell\in L$ under strategy $s\in S$ , and let $c_{\ell sp}$ denote the cost of using recovery plan $p\in P$ on the buildings in location $\ell\in L$ if the retrofitting strategy $s\in S$ is applied. All parameters in both stages are computed in Section 6.1 in a real example, by using the publicly available data and engineering models for wood-frame buildings in the literature.

Define the decision variables $f_{\ell s}$ and $r_{\ell sp}$ as

f_{\ell s}=\begin{cases}1,&\text{ if the buildings in location $\ell\in L$ are retrofitted using strategy $s\in S$}\\ 0,&\text{ otherwise},\end{cases}

r_{\ell sp}=\begin{cases}1,&\text{ if the buildings in $\ell\in L$ are retrofitted with strategy $s\in S$}\\ &\text{ and recovered with plan $p\in P$}\\ 0,&\text{ otherwise},\end{cases}

and let $z_{\ell}$ , $\ell\in L$ , be the binary variables that represent the coverage of a tornado across the locations, that is:

z_{\ell}=\begin{cases}1\text{ if location $\ell\in L$ is affected by the tornado}\\ 0\text{ otherwise}.\end{cases}

Then, the two-stage robust optimization model is defined as follows


$\displaystyle v=\min$	$\displaystyle\sum_{\ell\in L}\sum_{s\in S}w_{\ell s}f_{\ell s}+\max_{z\in\mathcal{U}}Q(z,f)$		(1a)
s.t.	$\displaystyle\sum_{s\in S}f_{\ell s}=1$	$\displaystyle\forall\ell\in L$	(1b)
	$\displaystyle f\in\{0,1\}^{\|L\|\|S\|},$		(1c)

where the second stage problem $Q(z,f)$ for given $z=(z_{\ell}:\ell\in L)$ and $f=(f_{\ell s}:\ell\in L,s\in S)$ is


$\displaystyle Q(z,f)=\min$	$\displaystyle\sum_{\ell\in L}z_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}$		(2a)
s.t.	$\displaystyle\sum_{\ell\in L}\sum_{s\in S}\sum_{p\in P}c_{\ell sp}r_{\ell sp}\leq A-\sum_{\ell\in L}\sum_{s\in S}d_{\ell s}f_{\ell s}$		(2b)
	$\displaystyle\sum_{p\in P}r_{\ell sp}=f_{\ell s}$	$\displaystyle\forall\ell\in L,s\in S$	(2c)
	$\displaystyle r\in\{0,1\}^{\|L\|\|S\|\|P\|}.$		(2d)

The first term of the objective function (1a) measures the pre-tornado dislocation across all locations due to implementing the retrofitting strategies. The second term measures the population dislocation after the worst-case tornado with respect to the retrofitting strategies selected at the first stage. Constraint (1b) ensures that each location picks exactly one retrofitting strategy. We make the assumption that there is a “do-nothing” strategy in $S$ with zero cost.

The second stage optimization problem (2) selects the recovery plans that minimize the population dislocation given a retrofitting strategy vector $f$ and the tornado coverage vector $z$ . The budget constraint (2b) ensures that the overall cost for the execution of recovery and retrofitting strategies is $A$ . Constraint (2c) forces each location to pick exactly one recovery strategy corresponding to the selected retrofitting strategy. As before, we assume there exist a “do-nothing” strategy in $P$ with zero cost. Next, we model tornado paths by ensuring that the decision vector $z$ belongs to the uncertainty set $\mathcal{U}\subseteq\{0,1\}^{|L|}$ which contains all feasible tornado damages over locations in $L$ .

3.2 Uncertainty set formulation

We make three assumptions to model a tornado and its associated damage:

A.1.

Tornado paths are represented by line segments of at most a length $E$ in the (real) plane; $E>0$ is a parameter set by the decision maker, and can be easily modified to study the effects of different tornado lengths.
A.2.

Each location $\ell\in L$ represents a cluster of buildings that are geographically close. Such clusters can depict individual buildings, blocks, neighborhoods, counties, or any other arrangement of interest to the decision-maker.
A.3.

A location $\ell\in L$ is damaged by the tornado if an aggregate measure of the coordinates of all buildings in $\ell$ (e.g., the centroid of all the coordinates of the buildings in that location), is within $\Delta$ units of the tornado path, $\Delta>0$ being a parameter set by the decision-maker.

Assumption A.1. follows from considering historical records of tornadoes, such as evidenced in Figure 1, which depicts the tornado paths in Oklahoma County between 1950 and 2020. Observe that, with a few exceptions, the paths can be adequately approximated by line segments. This behavior is also observed across the US, not only Oklahoma, see the damage assessment tool in NOAA [2023]. Whereas A.1. might be seen as restrictive, it is a reasonable and uncluttered mathematical representation of a complex weather phenomenon.

Refer to caption — Figure 1: Map of Oklahoma County tornado paths between 1950-2020 [NWC, 2023].

Assumption A.2. is made to avoid politically inviable policies. Allowing full granularity (by associating a location with each building) might result in policies that prescribe retrofitting and recovery actions for a given building, but no retrofitting nor recovery to its immediate neighbors. Clearly, such a policy might raise fairness concerns, particularly if the decision-maker is a public entity. We note that our model works for any clustering approach; all clustering information required by the model is summarized via the parameters of the model.

Assumption A.3. defines the damage due to the tornado. Considering the variable nature of tornado damage, $\Delta$ can be seen as an upper-bound on the width observed in historical data for a given area to study the worst-case scenario. Due to the uncertain nature of tornado behavior, defining specific rules for modeling the coverage during a tornado remains a challenging task.

We assume that the length of tornado $q$ is at most a given parameter $E$ , and that the two endpoints of $q$ are located within a sufficiently large rectangle $R=[R_{L}^{(1)},R_{U}^{(1)}]\times[R_{L}^{(2)},R_{U}^{(2)}]\subseteq\mathbb{R}^{2}$ that contains all locations in $L$ . Let $e_{0}=(e_{0}^{x},e_{0}^{y})\in R$ and $e_{1}=(e_{1}^{x},e_{1}^{y})\in R$ be the endpoints of line segment $q$ (see Figure 2). The uncertainty set $\mathcal{U}$ can be defined by the following mixed-integer non-linear programming (MINLP) formulation


$\displaystyle\mathcal{U}=\bigg{\{}z\in\{0,1\}$	${}^{\|L\|}:\exists e_{0},e_{1}\in R,t\in\mathbb{R}^{\|L\|},$	(3a)
	$\displaystyle\\|e_{0}+t_{\ell}(e_{1}-e_{0})-(x_{\ell},y_{\ell})\\|\leq\Delta+M_{\ell}(1-z_{\ell}),\ \forall\ell\in L,$	(3b)
	$\displaystyle 0\leq t_{\ell}\leq 1,\ \forall\ell\in L,$	(3c)
	$\displaystyle\\|e_{0}-e_{1}\\|\leq E\bigg{\}},$	(3d)

where $\|.\|$ calculates the Euclidean distance between the vector coordinates and where $(x_{\ell},y_{\ell})$ denote the coordinates of location $\ell\in L$ on the plane.

The value of $\|e_{0}+t_{\ell}(e_{1}-e_{0})-(x_{\ell},y_{\ell})\|$ in Constraint (3b) calculates the distance between a point $e_{0}+t_{\ell}(e_{1}-e_{0})$ on the line segment $q$ and location $\ell$ . This constraint defines $z$ -variables such that if $z_{\ell}=1$ for some $\ell\in L$ , then the location $\ell$ is covered by tornado. The value of $M_{\ell}$ is sufficiently large to ensure the constraint is trivially satisfied if $z_{\ell}=0$ (the value of $M_{\ell}$ is set to be the largest value among the Euclidean distances from location $\ell\in L$ to the corners of rectangle $R$ ). Constraint (3c) enforces the range of continuous variable $t_{\ell}$ to be between 0 and 1, which indicates that $q$ is a finite segment and not an infinite line. Constraint (3d) also ensures that the length of a tornado, which is the Euclidean distance between the two endpoints, does not exceed $E$ . Note that the problem maximizes over the $z$ -variables and that all of them have non-negative coefficients in the objective function. This observation and Constraint (3b) imply that $z_{\ell}=1$ if and only if location $\ell\in L$ is covered by the tornado path.

Remark 1.

The parameters $\Delta$ , $E$ , $R_{L}^{(i)}$ , and $R_{U}^{(i)}$ , $i=1,2$ , are selected by the decision-maker. Also, we note that $E=\infty$ is a valid choice for the length of the tornado path. In this case we refer to the path as a line rather than by a segment. This assumption represents situations where the path is sufficiently long to traverse the entire region of interest. In this case, constraints (3c) and (3d) are not necessary in the formulation of $\mathcal{U}$ .

We close this section by noting that problem (1) belongs to the class of NP-hard problems, see Proposition 3 in Appendix 1. The proof follows by a reduction from the well-known knapsack problem and applies even if there is only one tornado path.

4 A Solution Method Based on the C&CG Framework

We present a method to solve the two-stage robust optimization problem (1), based on the C&CG algorithm in Zeng and Zhao [2013]. First, we provide a linear one-level reformulation of problem (1) by means of an epigraphic reformulation of the worst-case second-stage objective. To this end, note that $\mathcal{U}$ is a finite set and write it as $\mathcal{U}=\{z^{1},\ldots,z^{n}\}$ , where $n\geq 1$ is the cardinality of $\mathcal{U}$ and suppose $r^{i}$ is the vector of second-stage variables corresponding to scenario $z^{i}$ . The one-level linear mixed-integer programming (MIP) reformulation of (1) is:


$\displaystyle v=\min$	$\displaystyle\sum_{\ell\in L}\sum_{s\in S}w_{\ell s}f_{\ell s}+\theta$		(4a)
s.t.	$\displaystyle\sum_{s\in S}f_{\ell s}=1$	$\displaystyle\forall\ell\in L$	(4b)
	$\displaystyle\theta\geq\sum_{\ell\in L}z^{i}_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}^{i}$	$\displaystyle\forall i\in\mathcal{I}$	(4c)
	$\displaystyle\sum_{\ell\in L}\sum_{s\in S}\sum_{p\in P}c_{\ell sp}r^{i}_{\ell sp}+\sum_{\ell\in L}\sum_{s\in S}d_{\ell s}f_{\ell s}\leq A$	$\displaystyle\forall i\in\mathcal{I}$	(4d)
	$\displaystyle\sum_{p\in P}r^{i}_{\ell sp}=f_{\ell s}$	$\displaystyle\forall\ell\in L,s\in S,i\in\mathcal{I}$	(4e)
	$\displaystyle f\in\{0,1\}^{\|L\|\|S\|}$		(4f)
	$\displaystyle r^{i}\in\{0,1\}^{\|L\|\|S\|\|P\|}$	$\displaystyle\forall i\in\mathcal{I},$	(4g)

where $\mathcal{I}=[n]:=\{1,\ldots,n\}$ . Constraint (4c) implies that the minimum value of $\theta$ is equal to the maximum population dislocation among all tornadoes in $\mathcal{U}$ . Constraint (4d) and (4e) ensure that the recourse vectors $r^{i}$ corresponding to each scenario, meet the constraints of the second stage.

Remark 2.

Solving problem (4) directly is not practical in general, given that the uncertainty set $\mathcal{U}$ might have a large number of scenarios in terms of the number of locations. If $E=\infty$ , then there is a polynomial number of scenarios in $\mathcal{U}$ which can be constructed based on the stabbing line algorithm in Section 5.3.1 (in this case, nevertheless, using (4) to directly solve the problem remains impractical because there would be at least a quadratic number of scenarios in terms of the number of locations). If $E<\infty$ , it remains an open question to determine whether there are polynomially many scenarios in $\mathcal{U}$ .

We employ a modified version of the C&CG to solve formulation (4). The algorithm is set to iteratively solve a relaxation of (4) over a subset $\mathcal{I}^{k}\subseteq\mathcal{I}$ of scenarios at iteration $k$ to obtain an optimal solution $(\theta^{k},f^{k})$ . If given $f^{k}$ , there exists a tornado coverage $z^{i_{k}}\in\mathcal{U}$ that produces a dislocation greater than $\theta^{k}$ , then we add scenario $z^{i_{k}}$ to the master problem at the next iteration. Scenario $i_{k}$ is found by solving the following subproblem, evaluated at $f^{k}$ :

\displaystyle\Phi(f)=\max_{z\in\mathcal{U}}\min\left\{\sum_{\ell\in L}z_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}:r\in\mathcal{R}(f)\right\},

(5)

where $\mathcal{R}(f)$ is the set of feasible solutions of problem $Q(z,f)$ ; see constraints (2b)–(2d). Specifically, if the optimal solution is such that $\Phi(f^{k})>\theta^{k}$ , then $z^{i_{k}}$ is the tornado coverage $z$ that attains $\Phi(f^{k})$ . Note that when the variables and constraints associated with scenario $i_{k}$ are added at iteration $k+1$ , the solution $(\theta^{k},f^{k})$ is no longer feasible in the master relaxation.

Let $v^{k}$ be the value of the master problem at iteration $k$ , that is, problem (4) with $\mathcal{I}$ replaced by $\mathcal{I}^{k}$ , with $\mathcal{I}^{k}=\{i_{0},i_{1},\ldots,i_{k-1}\}$ and let $\Phi(f^{k})$ be the subproblem for the retrofitting strategy $f^{k}$ found at iteration $k$ of the algorithm. The C&CG method is presented in Algorithm 1.

Data: Set

\mathcal{U},L,S

, and

P

Result:

f^{k}

1 Set

k=0

\mathcal{I}^{1}=\{{i_{0}}\}\subseteq\mathcal{U}

LB=-\infty

, and

UB=\infty

;

2 while $UB-LB>0$ do

3 Set

k=k+1

;

4 Solve the master MIP

v^{k}

and let

f^{k}

be the optimal retrofitting strategy. Update

LB=v^{k}

;

5 Solve the subproblem

\Phi(f^{k})

and let

{i_{k}}

be the index of the

z

-optimal solution in

\mathcal{I}

. Update

UB=\min\{UB,\sum_{\ell\in L}\sum_{s\in S}w_{\ell s}f^{k}_{\ell s}+\Phi(f^{k})

} and

\mathcal{I}^{k+1}=\mathcal{I}^{k}\cup\left\{{i_{k}}\right\}

;

6 end while

Algorithm 1 C&CG algorithm to solve

v

Algorithm 1 begins with an arbitrary feasible scenario ${i_{0}}\in\mathcal{I}$ . In each iteration $k$ , the relaxed MIP problem $v^{k}$ is solved. The value of $v^{k}$ then provides a lower bound for the value of $v$ . We note that the values of $v^{k}$ , $k\geq 0$ , are non-decreasing in $k$ . Then, the subproblem $\Phi(f^{k})$ generates a new feasible scenario ${i_{k}}$ for the next iteration and updates the upper bound. The algorithm terminates when the subproblem does not find a scenario that violates the lower bound value, which implies the upper bound is equal to the lower bound. Zeng and Zhao [2013] prove that the C&CG algorithm converges in finite time.

Note that Step 1 in Algorithm 1 is expensive computationally because it requires solving the MINLP problem $\Phi(f)$ in (5). The approach to solve the subproblem in standard applications of the C&CG relies on using duality or the Karush Kuhn Tucker (KKT) conditions on the second-stage problem [Zeng and Zhao, 2013]. Notice that $\Phi(f)$ is a bilevel problem with integrality requirements at both levels. Therefore, we cannot use either of these approaches. Alternatively, we propose a decomposition branch-and-cut method next.

5 A Decomposition Branch-and-Cut (DBC) Algorithm to Solve $\Phi(f)$

The DBC method is based on a one-level reformulation of $\Phi(f)$ . We replace the $\min$ in the objective by its epigraphic reformulation and replace the uncertainty set requirements in terms of linear (although exponentially many in the worst-case) constraints. Next, we discuss how to initialize the constraints in the master problem and how to perform the “two types” of separations.

5.1 Overview of the DBC algorithm

Suppose $\mathcal{C}$ is the set of combinations of locations that cannot be covered by a tornado path , i.e., for each $z\notin\mathcal{U}$ , there exists a set of locations $C=\{\ell\in L:z_{\ell}=1\}\in\mathcal{C}$ . A one-level linear reformulation of $\Phi(f)$ is given by:


$\displaystyle\Phi(f)=\max\$	$\displaystyle\eta$		(6a)
s.t.	$\displaystyle\eta\leq\sum_{\ell\in L}z_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}$	$\displaystyle\forall r\in\mathcal{R}(f)$	(6b)
	$\displaystyle\sum_{\ell\in C}z_{\ell}\leq\|C\|-1$	$\displaystyle\forall C\in\mathcal{C}$	(6c)
	$\displaystyle z\in\{0,1\}^{\|L\|}.$		(6d)

The correctness of the reformulation follows from the maximization sense in (6). Note that $r$ in (6b) represents a fixed value (not a variable); and that there is such an $r$ for each element of $\mathcal{R}(f)$ .

Observe that $\mathcal{R}(f)$ and $\mathcal{C}$ are finite sets with potentially exponentially many elements, thus $\Phi(f)$ cannot be solved directly by a commercial solver. The DBC starts with a relaxation of $\Phi(f)$ by replacing $\mathcal{R}(f)$ and $\mathcal{C}$ with subsets $\mathcal{R}^{0}(f)$ and $\mathcal{C}^{0}$ , respectively. The resulting problem is referred as the master relaxation. The standard branch-and-cut (BC) algorithm then starts solving the master and prunes nodes by bound and infeasibility. When at some node $h$ an integral solution $z^{h}$ is found (i.e., a solution that is feasible in the master), its feasibility with respect to the original formulation in (6) must be verified.

To this end, Algorithm 2 takes the master feasible solution $(\eta^{h},z^{h})$ at node $h$ of the BC tree and verifies its feasibility: first, it checks whether $z^{h}$ satisfies all $\mathcal{C}$ -constraints (6c), i.e., it checks if $z^{h}\in\mathcal{U}$ . If $z^{h}$ does not pass the first check, then let $C^{h}=\{\ell\in L\colon z_{\ell}^{h}=1\}$ be the active locations in $z^{h}$ . Observe that $C^{h}$ is a combination of locations that cannot be covered by a tornado path, thus a cut (6c) with $C=C^{h}$ is added on the fly. If $z^{h}$ passes the first check, then the algorithm checks whether $z^{h}$ satisfies all $\mathcal{R}(f)$ -constraints (6b). If $z^{h}$ also passes the second check, then $z^{h}$ is feasible in $\Phi(f)$ ; if not, then a cut (6b) with $r=r^{h}$ is added on the fly, $r^{h}$ being the optimal solution of the second-stage problem associated with $f$ and $z^{h}$ . The performance of Algorithm 2 depends on the quality of $\mathcal{C}^{0}$ and $\mathcal{R}^{0}(f)$ and on having fast separation routines. We discuss these issues next.

Data: Master feasible solution

(\eta^{h},z^{h})

Result: Generating cuts on the fly, if

(\eta^{h},z^{h})

is infeasible

1 if $z^{h}\notin\mathcal{U}$ then

2 Define

C^{h}=\{\ell\in L:z_{\ell}^{h}=1\}

;

3 Add constraint (6c) with

C=C^{h}

on the fly;

5 else if $\eta^{h}>Q(z^{h},f)$ then

6 Define

r^{h}=\arg\min\{\sum_{\ell\in L}z^{h}_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}:r\in\mathcal{R}(f)\}

;

7 Add constraint (6b) with

r=r^{h}

on the fly;

9 else

(\eta^{h},z^{h})

is feasible.

11 end if

Algorithm 2 Separation Procedure for the DBC

5.2 Definition of $\mathcal{C}^{0}$

The set $\mathcal{C}^{0}$ is constructed by identifying pairs and triples of locations that cannot be covered by a single tornado path. It is worth mentioning that the valid constraints associated with these infeasible cases can also be used to tighten the relaxation of $\mathcal{U}$ by the elimination of potentially many infeasible fractional solutions for $z$ -variables; we evaluate the performance of DBC with $\mathcal{U}$ enhanced by these valid cuts in Appendix 5.

5.2.1 Cuts from infeasible pairs.

Because the length of a tornado is at most $E$ and only locations within a distance of $\Delta$ from a tornado line segment are covered, a tornado cannot cover two locations $\ell_{1},\ell_{2}\in L$ if they are at a distance of more than $E+2\Delta$ . Consequently we have the following result:

Proposition 1.

Let $\Omega=\left\{(\ell_{1},\ell_{2})\in L^{2}:\|(x_{\ell_{1}},y_{\ell_{1}})-(x_{\ell_{2}},y_{\ell_{2}})\|>2\Delta+E\right\}$ be the set of infeasible pairs. Then any $z\in\mathcal{U}$ satisfies that

z_{\ell_{1}}+z_{\ell_{2}}\leq 1\qquad\forall(\ell_{1},\ell_{2})\in\Omega.

(7)

The proof of Proposition 1 is given in Appendix 2.

5.2.2 Infeasible triple cuts.

We also identify infeasible triples of locations that cannot be covered by a tornado path. To introduce the infeasible set of triples, we assume $E=\infty$ which implies a line representation for a tornado path. Observe that if a line does not cover a subset of locations $Z\subseteq L$ , then there is no segment of the line to cover the locations in $Z$ either. For any locations $\ell_{1}$ and $\ell_{2}$ , the analysis below partitions the plane $({\mathbb{R}}^{2})$ in two regions, denoted by $R^{\prime}_{\ell_{1},\ell_{2}}$ and $R^{\prime\prime}_{\ell_{1},\ell_{2}}$ . Region $R^{\prime}_{\ell_{1},\ell_{2}}$ is the set of points $\ell\in\mathbb{R}^{2}$ for which there exist a line that covers $\ell_{1},\ell_{2}$ and $\ell$ , and $R^{\prime\prime}_{\ell_{1},\ell_{2}}=\mathbb{R}^{2}\setminus R^{\prime}_{\ell_{1},\ell_{2}}$ . For simplicity, in the next discussion we interchangeably refer to a point in ${\mathbb{R}}^{2}$ by its name, as $\ell\in{\mathbb{R}}^{2}$ , or by its two-dimensional coordinates, as in $(x_{\ell},y_{\ell})\in{\mathbb{R}}^{2}$ .

Observe that a line $\lambda\subset\mathbb{R}^{2}$ covers both $\ell_{1}$ and $\ell_{2}$ if and only if $\lambda\cap B_{1}\not=\emptyset$ and $\lambda\cap B_{2}\not=\emptyset$ , where $B_{i}$ is the ball with center $\ell_{i}$ and radius $\Delta$ , $i=1,2$ . Let $\mathcal{L}_{\ell_{1},\ell_{2}}$ be the collection of these lines, and let $P_{\ell_{1},\ell_{2}}$ be the set of points that belong to any such line:

P_{\ell_{1},\ell_{2}}=\{\ell\in\mathbb{R}^{2}\colon\exists\lambda\in\mathcal{L}_{\ell_{1},\ell_{2}}\text{ s.t. }\ell\in\lambda\}.

(8)

For any set $T\subseteq\mathbb{R}^{2}$ let $T(\Delta)$ be the set of points that are at a distance at most $\Delta$ from $T$ , that is

T(\Delta)=\{(x_{\ell},y_{\ell})\in\mathbb{R}^{2}\colon\exists(x_{t},y_{t})\in T\text{ s.t. }||(x_{\ell},y_{\ell})-(x_{t},y_{t})||\leq\Delta\}.

(9)

It is readily seen that $R^{\prime}_{\ell_{1},\ell_{2}}=P_{\ell_{1},\ell_{2}}(\Delta)$ .

Now, the set $P_{\ell_{1},\ell_{2}}$ can be characterized by the four lines $\tau_{1}$ , $\tau_{2}$ , $\sigma_{1}$ and $\sigma_{2}$ that are tangent to both $C_{1}$ and $C_{2}$ , where $C_{1}$ and $C_{2}$ are the circles enclosing $B_{1}$ and $B_{2}$ , see Figure 3(a). Specifically, let $\tilde{P}_{\ell_{1},\ell_{2}}$ be the set enclosed by these tangents, see Figure 3(b), we next show that $\tilde{P}_{\ell_{1},\ell_{2}}=P_{\ell_{1},\ell_{2}}$ .

Indeed, let $\lambda_{0}$ be the line that passes through $\ell_{1}$ and $\ell_{2}$ and suppose that a point $\ell$ is in between $\tau_{1}$ and $\tau_{2}$ . Then, the line $\lambda$ that passes through $\ell$ and is parallel to $\lambda_{0}$ intersects both $B_{1}$ and $B_{2}$ . Suppose that $\ell$ belongs to the region left of $\ell_{1}$ between $\tau_{2}$ and $\sigma_{2}$ . Then, the line $\lambda$ that passes through $\ell$ and the point $\sigma_{2}\cap C_{2}$ intersects both $B_{1}$ and $B_{2}$ . Notice that this argument can be repeated for the remaining regions of $\tilde{P}_{\ell_{1},\ell_{2}}$ , and thus we can conclude that $\tilde{P}_{\ell_{1},\ell_{2}}\subseteq P_{\ell_{1},\ell_{2}}$ .

Conversely, let $\ell\in P_{\ell_{1},\ell_{2}}$ and let $\lambda\in\mathcal{L}_{\ell_{1},\ell_{2}}$ be the line intersecting $B_{1}$ and $B_{2}$ such that $\ell\in\lambda$ . If $\lambda$ is one of the tangents, then it is clear that $\lambda\subseteq\tilde{P}_{\ell_{1},\ell_{2}}$ and it follows that $\ell\in\tilde{P}_{\ell_{1},\ell_{2}}$ . Thus, assume that $\lambda$ is not one of the tangents and let $p_{1}$ and $p_{2}$ be two arbitrary elements of $\lambda\cap B_{1}$ and $\lambda\cap B_{2}$ . Observe that the segment $s$ of $\lambda$ between $p_{1}$ and $p_{2}$ must be completely contained in between $\tau_{1}$ and $\tau_{2}$ . Moreover, as $\lambda$ is not one of the tangents, $s$ intersects both $\sigma_{1}$ and $\sigma_{2}$ . Because two straight lines can intersect at most once, the previous observations imply that the $\lambda$ cannot intersect $\sigma_{1}$ nor $\sigma_{2}$ at any point left of $p_{1}$ and right of $p_{2}$ ; in other words, it must be the case that $\lambda$ is completely contained in $\tilde{P}_{\ell_{1},\ell_{2}}$ , which implies that $\ell\in\tilde{P}_{\ell_{1},\ell_{2}}$ , as desired.

Given the previous considerations, in order to characterize $R^{\prime}_{\ell_{1},\ell_{2}}$ it is necessary to characterize $\tilde{P}_{\ell_{1},\ell_{2}}(\Delta)$ . From the shape of $\tilde{P}_{\ell_{1},\ell_{2}}$ , it is clear that $\tilde{P}_{\ell_{1},\ell_{2}}(\Delta)$ looks very similar to $\tilde{P}_{\ell_{1},\ell_{2}}$ , with the difference that its border starts $\Delta$ units up and $\Delta$ units down, see Figure 4(a)

Next, we use geometric facts to derive an explicit equation for lines defining the boundaries of $\tilde{P}_{\ell_{1},\ell_{2}}(\Delta)$ for a given $\ell_{1}$ and $\ell_{2}$ . For ease of exposition hereafter we assume that $x_{\ell_{1}}\leq x_{\ell_{2}}$ and $y_{\ell_{1}}\leq y_{\ell_{2}}$ . Define $\theta=\arctan{\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}}$ to be the angle of slope of the line $\lambda_{0}$ which connects $\ell_{1}$ and $\ell_{2}$ with respect to the horizontal axis and let $\alpha=\arcsin{\frac{2\Delta}{D(\ell_{1},\ell_{2})}}$ , where $D(\ell_{1},\ell_{2})$ is the Euclidean distance between $\ell_{1}$ and $\ell_{2}$ . Then we have that $\widetilde{P}_{\ell_{1},\ell_{2}}(\Delta)=\{(x,y)\in\mathbb{R}^{2}\colon$ Eqs (10a)–(10f) hold $\}$ , where


$\displaystyle y$	$\displaystyle\leq\tan{\theta}\cdot(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos{\theta}}$	$\displaystyle(\lambda_{1})$	(10a)
$\displaystyle y$	$\displaystyle\geq\tan{\theta}\cdot(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos{\theta}}$	$\displaystyle(\lambda_{2})$	(10b)
$\displaystyle y$	$\displaystyle\leq\tan{(\theta+\alpha)}\cdot(x-x_{\ell_{1}})+y_{\ell_{1}}$	$\displaystyle(\lambda_{3})$	(10c)
$\displaystyle y$	$\displaystyle\geq\tan{(\theta+\alpha)}\cdot(x-x_{\ell_{2}})+y_{\ell_{2}}$	$\displaystyle(\lambda_{4})$	(10d)
$\displaystyle y$	$\displaystyle\leq\tan{(\theta-\alpha)}\cdot(x-x_{\ell_{2}})+y_{\ell_{2}}$	$\displaystyle(\lambda_{5})$	(10e)
$\displaystyle y$	$\displaystyle\geq\tan{(\theta-\alpha)}\cdot(x-x_{\ell_{1}})+y_{\ell_{1}}$	$\displaystyle(\lambda_{6}).$	(10f)

The lines $\lambda_{i}$ , $i=0,\ldots,6$ defined in Equations (10a)– (10f) are depicted graphically in Figure 4(b).

Now, let

\Gamma_{\ell_{1},\ell_{2}}=\bigl{\{}\ell\in L\setminus\{\ell_{1},\ell_{2}\}:(x_{\ell},y_{\ell})\in\mathbb{R}^{2}\setminus\tilde{P}_{\ell_{1},\ell_{2}}(\Delta)\bigr{\}}.

(11)

The previous discussion is summarized in the following result, which states that if a location belongs to the infeasible set $\Gamma_{\ell_{1},\ell_{2}}$ then at most two locations among $\ell_{1},\ell_{2}$ , and $\ell_{3}$ can be chosen to be covered by a tornado. A more formal proof of this result is given in Appendix 3.

Proposition 2.

Let $\ell_{1}$ and $\ell_{2}$ be two elements of $L$ and suppose that $D(\ell_{1},\ell_{2})>2\Delta$ . Then any $z\in\mathcal{U}$ satisfies that

z_{\ell_{1}}+z_{\ell_{2}}+z_{\ell}\leq 2\qquad\forall\ell\in\Gamma_{\ell_{1},\ell_{2}}.

(12)

It is worth mentioning that we do not need to consider all possible combinations of locations to construct the infeasible triple cuts because if $\ell_{3}\in\Gamma_{\ell_{1},\ell_{2}}$ , then $\ell_{1}\in\Gamma_{\ell_{2},\ell_{3}}$ and $\ell_{2}\in\Gamma_{\ell_{1},\ell_{3}}$ . Also, we note that the pair-wise and triple valid constraints introduced so far do not completely characterize $\mathcal{U}$ , as shown in the next remark. However, they can be used to define the initial subset $\mathcal{C}^{0}\subseteq\mathcal{C}$ of the master relaxation in the DBC as follows:

\mathcal{C}^{0}=\big{\{}(\ell_{1},\ell_{2})\in\Omega\big{\}}\cup\big{\{}(\ell_{1},\ell_{2},\ell_{3}):\ell_{3}\in\Gamma_{\ell_{1},\ell_{2}},\ell_{1},\ell_{2}\in L,D(\ell_{1},\ell_{2})>2\Delta\big{\}}.

(13)

Remark 3.

Define $\mathcal{V}=\{z\in\{0,1\}^{|L|}:\sum_{\ell\in C}z_{\ell}\leq|C|-1,\forall C\in\mathcal{C}^{0}\}$ . The set $\mathcal{V}$ does not necessarily determine all the feasible solutions in the uncertainty set $\mathcal{U}$ ; specifically $\mathcal{U}\subset\mathcal{V}$ and in general there exist points in $\mathcal{V}$ that do not belong to $\mathcal{U}$ . For example, consider a case of three locations $\ell_{1},\ell_{2}$ , and $\ell_{3}$ with coordinates $(0,0),(4,0)$ , and $(2,1.1)$ on a plane under a tornado with the maximum length of $E=2$ and the width of $\Delta=1$ (See Figure 5). It can be easily observed that there exists no tornado path that covers all three locations because the only tornado path covering $\ell_{1}$ and $\ell_{2}$ is the line segment with endpoints $(1,0)$ and $(3,0)$ . Obviously, the location $(2,1.1)$ is not covered by this path. However, $z_{\ell_{1}}=z_{\ell_{2}}=z_{\ell_{3}}=1$ is a feasible solution in $\mathcal{V}$ because it satisfies both the infeasible pair and triple cuts.

5.3 Feasibility check in $\mathcal{C}$

Next, we present a method to determine whether a master feasible solution $z^{h}$ belongs to $\mathcal{U}$ . We start our analysis for the case when $E=\infty$ and then present a method for the case when $E<\infty$ .

5.3.1 Feasibility check for $E=\infty$ using the stabbing line algorithm (SLA).

The stabbing line problem consists of a collection of (disjoint) circles $\mathcal{O}=\{O_{1},O_{2},\ldots,O_{n}\}$ and the objective is to find the line on the plane that intersects the most circles. The SLA is an algorithm that solves this problem and reports $c(\mathcal{O})$ : the maximum number of circles in $\mathcal{O}$ that can be intersected with a line. Consider the collection of circles $\mathcal{O}^{h}$ with centers at the locations in $C^{h}$ and radius $\Delta$ . Observe that a tornado path with $E=\infty$ covers all locations in $C^{h}$ if and only if the $c(\mathcal{O}^{h})=|C^{h}|$ (see Figure 6). Therefore, if $c(\mathcal{O}^{h})<|C^{h}|$ , then $z^{h}\not\in\mathcal{U}$ whereas if $c(\mathcal{O}^{h})=|C^{h}|$ then $z^{h}\in\mathcal{U}$ .

We explain next the SLA following Saba [2013]. Assume first a set of disjoint circles $\mathcal{O}=\{O_{1},O_{2},\ldots,O_{m}\}$ , $m>0$ . Given any two disjoint circles $O_{a}$ and $O_{b}$ , find the two external and two internal tangent lines to each pair of circles $(O_{a},O_{b})$ and record the intersection points on both circles; see Figure 7. Let $p_{ab}^{E},q_{ab}^{E}$ and $p_{ab}^{I},q_{ab}^{I}$ be the intersection points of these external and internal tangent lines, respectively, with $O_{a}$ . Note that a tangent line to $O_{a}$ that touches any point $p$ in the arc $[p_{ab}^{E},p_{ab}^{I}]$ (or $[q_{ab}^{I},q_{ab}^{E}]$ ) intersects $O_{b}$ , see Figure 7.

Given $O_{a}$ , the SLA computes the arc-intervals $[p_{ab}^{E},p_{ab}^{I}]$ and $[q_{ab}^{I},q_{ab}^{E}]$ for all $b\not=a$ ; note there exist exactly $2(m-1)$ such intervals for $O_{a}$ . Let $\pi_{a}$ be an endpoint of any such arc-intervals that belong to the maximum amount of arc-intervals. Then, the tangent line that touches $O_{a}$ at $\pi_{a}$ is a line intersecting $O_{a}$ that intersects the most circles in $\mathcal{O}$ . Iterating this process starting from each $a=1,\ldots,m$ , will determine a line that intersects the most circles. The SLA can be implemented $O(m^{2}\log m)$ time; see Saba [2013] for further details on the algorithm.

In our case the circles in $\mathcal{O}^{h}$ might be non-intersecting, as the locations in $C^{h}$ might be at a distance less than $2\Delta$ . We thus adapt the SLA, noting that if $O_{a}\cap O_{b}\not=\emptyset$ then the internal tangents are not necessary and that doing the analysis with the arc-interval $[p_{ab}^{E},q_{ab}^{E}]$ is sufficient for this case.

5.3.2 Feasibility check if $E<\infty$ .

In this case, if $c(\mathcal{O}^{h})<|C^{h}|$ , then we can again conclude that there is no finite line segment that covers $C^{h}$ and thus $z^{h}\not\in\mathcal{U}$ . Otherwise, if $c(\mathcal{O}^{h})=|C^{h}|$ , then there is an infinite line that covers all locations in $C^{h}$ , but there is no guarantee that a line segment of length at most $E$ covers all locations of $C^{h}$ . Thus, we complete the feasibility check in two steps: first, we use the line resulting from the SLA and projection ideas to determine whether the line can be trimmed to be of length $E$ without compromising its coverage, the details are in Appendix 4. If even after this step the feasibility of $z^{h}$ cannot be guaranteed, then we check the feasibility of $z^{h}$ algebraically using a simplified version of $\mathcal{U}$ .

Specifically, consider the set $\tilde{\mathcal{U}}^{h}$ defined as


$\displaystyle\tilde{\mathcal{U}}^{h}=\bigg{\{}(e_{0},e_{1})\in$	$\displaystyle R^{2}:\exists t\in[0,1]^{\|C^{h}\|},\\|e_{0}-e_{1}\\|\leq E$	(14a)
	$\displaystyle\\|e_{0}+t_{\ell}(e_{1}-e_{0})-(x_{\ell},y_{\ell})\\|\leq\Delta,\ \forall\ell\in C^{h}\bigg{\}}.$	(14b)

Observe that $z^{h}\in\mathcal{U}$ if and only if $\tilde{\mathcal{U}}^{h}\not=\emptyset$ and that (14) is far simpler than $\mathcal{U}$ because it has fewer variables and does not involve binary variables nor ‘big-M’ values. Checking whether $\tilde{\mathcal{U}}^{h}$ is non-empty can be done using a quadratic (non-convex) continuous optimization solver. The such check can be significantly speed-up by adding linear cuts in $\tilde{\mathcal{U}}^{h}$ , see the details in Appendix 4.

5.4 Definition of $\mathcal{R}^{0}(f)$ and feasibility check for $\mathcal{R}(f)$

To initialize the master problem we assume that $\mathcal{R}^{0}(f)=\{r^{0}\}$ where $r^{0}$ is an optimal solution for the second-stage problem $Q(z,f)$ assuming no location is hit, i.e., $r^{0}$ is a solution of $Q(\mathbf{0},f)$ . On the other hand, once it has been checked that $z^{h}\in\mathcal{U}$ , we check that $z^{h}$ satisfies all $\mathcal{R}(f)$ -constraints in (6b) by solving the second-stage problem evaluated at $z^{h}$ , i.e., by solving $Q(z^{h},f)$ . Let $r^{h}$ be an optimal solution of $Q(z^{h},f)$ . If $\eta^{h}>Q(z^{h},f)$ then $(\eta^{h},z^{h})$ is not feasible in (6b) and we add the constraint $\eta\leq\sum_{\ell\in L}z_{\ell}\sum_{s\in S}\sum_{p\in P}g_{\ell sp}r_{\ell sp}^{h}$ at all active nodes of the BC tree. Conversely, if $\eta^{h}\leq Q(z^{f},f)$ , then $(\eta^{h},z^{h})$ is a feasible solution of (6) and node $h$ is pruned by feasibility after updating the lower bound of the BC. Observe that the initialization and separation procedure involve solving the $NP$ -hard problem $Q(z,f)$ ; however, this problem is relatively simple and can be solved quickly using commercial MIP solvers, see our computational results in Section 6.

6 Case study: retrofitting/recovery residential buildings in Joplin, MO

Next, we perform several numerical experiments with the proposed model. Our objectives are to study the performance of the two-stage optimization problem, to analyze how the optimal solutions change with the length of the path, to gather high-level insights into the optimal retrofitting and recovery activities, and to compare the optimal policies with other benchmark policies. The experiments are based on location and cost data from Joplin, MO. Joplin was the location of one of the most severe tornadoes in US history that occurred on May 22, 2011. The tornado was rated EF-5 with an average speed of more than 200 mph and crossed 22.1 miles on the ground with a width of one mile. This catastrophe caused the loss of 161 lives and more than a thousand injuries. The total damage was approximated at 3 billion dollars as around 7500 residential structures were destroyed [Kuligowski et al., 2014].

6.1 Definition of Parameters

The geographic and demographic data are based on information from Joplin and extracted from the online platform IN-CORE [2022]. Particularly, this dataset divides buildings in about 1500 “block IDs” (which are contiguous blocks of the city) and includes over 20,000 buildings. For our experiments, we use only the single or multi-story residential wood-frame building archetypes.

We assume an EF-2 tornado with a wind speed of 134 mph, a width of 0.75 miles (thus $\Delta$ =0.75/2), and a length of 5 and $\infty$ miles. Following Wang et al. [2021], we consider the retrofitting strategies presented in Table 1; these strategies are analyzed in Masoomi et al. [2018] with respect to the component fragility of wood-frame residential building archetypes. The retrofitting activities in Table 1, have minimal impact in the livability of the buildings; thus we assume that the first stage population dislocation parameter is zero for all locations and strategies, that is, $w_{\ell s}=0$ for all $\ell\in L,s\in S$ . This assumption, however, does not change the complexity of solving the problem, and if necessary could be relaxed to better depict a particular context or instance of interest.

Table 1: Retrofitting strategies for the experiments.

Strategy $s$	Description
$R_{1}$	Roof covering with asphalt shingles, using 8d nails to attach roof deck sheathing panel to rafters spaced at 6/12 inch, and the selection of two 16d toenails for roof-to-wall
$R_{2}$	Roof covering with asphalt shingles, using 8d nails to attach roof deck sheathing panel to rafters spaced at 6/6 inch, and the selection of two H2.5 clips for roof-to-wall
$R_{3}$	Roof covering with clay tiles, using 8d nails to attach roof deck sheathing panel to rafters spaced at 6/6 inch, and the selection of two H2.5 clips for roof-to-wall.

The cost of implementing these strategies are computed from the data of IN-CORE [2022]. As mentioned in Section 3, a ‘do-nothing’ strategy with zero cost is also included in the retrofitting strategy set $S$ , that is, $S$ ={do-nothing, R1, R2, R3}. On the other hand, only two recovery strategies are assumed in our experiment, namely, recover ( $p=1$ ) and do-nothing ( $p=0$ ); thus $P=\{0,1\}$ .

The second stage population dislocation parameters, $g_{\ell sp}$ , are defined as the expected population dislocation $60$ days after the tornado (Appendix 6 has additional experiments for 30 days of recovery). These values are evaluated based on the state of damage $d\in D$ immediately after the tornado hits, where $D$ is the set of damage states. Let $X_{\ell}$ be the random variable denoting the functionality of $\ell$ after 60 days; where $X_{\ell}=0$ , only if $\ell$ reaches 100% functionality after $60$ days (no dislocation) and $X_{\ell}=1$ otherwise (dislocation). Also, let $Y_{\ell}\in D$ be the random variable representing the damage at $\ell$ immediately after the tornado, and let $\sigma_{\ell}$ denote the retrofitting strategy. Then, the expected population dislocation after recovery is

g_{\ell s1}=N_{\ell}\times E[X_{\ell}|\sigma_{\ell}=s],

(15)

where $N_{\ell}$ is the population of location $\ell\in L$ . Since the value of $X_{\ell}$ depends on $Y_{\ell}$ , using conditional probability and the (reasonable) assumption that $X_{\ell}$ is independent of $\sigma_{\ell}$ given the knowledge of $Y_{\ell}$ , we have that

E[X_{\ell}|\sigma_{\ell}=s]=\sum_{d\in D}P[X_{\ell}=1|Y_{\ell}=d]P[Y_{\ell}=d|\sigma_{\ell}=s].

(16)

Following Koliou and van de Lindt [2020], the damage states are given by $D=\{$ Minor, Moderate, Extensive, Complete $\}$ . We compute $P[X_{\ell}=u|Y_{\ell}=d]$ using the repair time distributions described in FEMA [2003]. These distributions are lognormal and their parameters can be found in Koliou and van de Lindt [2020]. The probabilities $P[Y_{\ell}=d|\sigma_{\ell}=s]$ are computed from the building-level tornado fragility-curves for wood-frame residential buildings in Masoomi et al. [2018], Koliou and van de Lindt [2020].

The cost of recovery $c_{\ell s1}$ is the expected percentage of the total location replacement cost: suppose $R_{\ell}$ is the area of location $\ell\in L$ and $\alpha$ is the cost of replacement per $m^{2}$ . For the wood residential archetype, we set $\alpha=\$862.0/m^{2}$ and, the percentage of the location replacement cost to be $r_{minor}=0.5\%$ , $r_{moderate}=2.3\%$ , $r_{extensive}=11.7\%$ , and $r_{complete}=23.4\%$ , see Koliou and van de Lindt [2020]. The following formula is used to compute the cost of recovery for location $\ell$ if it is retrofitted by strategy $s$ :

c_{\ell s1}=\alpha\times R_{\ell}\times\sum_{d\in D}r_{d}P[Y=d|\sigma=s].

(17)

We assume the decision-maker takes charge of all recovery expenses.

On the other hand, if a location does not receive resources from the decision-maker, then the population dislocation will be higher and depends on the resident’s actions. We assume the population dislocation for the do-nothing strategy is a factor of the population dislocation resulting from the recovery actions of the decision-maker, i.e.,

g_{\ell s0}=\mu_{\ell s}g_{\ell s1}\quad\forall\ell\in L,s\in S.

(18)

The factor $\mu_{\ell s}$ can vary such that $g_{\ell s0}\in[g_{\ell s1},N_{\ell}]$ : $g_{\ell s0}=g_{\ell s1}$ implies that the residents pay a full amount of expenses by themselves to recover their buildings as quickly as the decision-maker, and $g_{\ell s0}=N_{\ell}$ means that the residents take no recovery action on their own and therefore the whole population is dislocated from location $\ell$ . For our experiments, we fix $g_{\ell s0}=(g_{\ell s1}+N_{\ell})/2$ which is the midpoint in the range. Note that the cost of the do-nothing strategy is $c_{\ell s0}=0,\forall\ell\in L,s\in S$ from the decision-maker’s point of view.

We remark that various parameters of our case study we selected rather arbitrarily; for instance, the assumption that recovery means that the building recovers 100% functionality; that the last day of recovery is 60 days; or that ‘half’ the residents recover their property without assistance. These are reasonable values that are selected for illustration purposes and can be changed by the decision-maker as desired. Generally, modifications of these parameters should not impact the model’s performance, and as long as there is data available that account for the modifications, our model should be able to run; see for example Appendix 6 where experiments with a recovery time of 30 days are reported.

6.2 Robust retrofitting and recovery strategies in Joplin, MO

We implement the proposed two-stage optimization model to find robust retrofitting and recovery strategies for locations in Joplin. All algorithms are coded in Python 3.9 and solved using Gurobi Optimizer 9.1.1 on a 64-bit Windows operating system with Intel(R) Core(TM) i7-8550U CPU and 8.00 GB of RAM under a time limit of 3600 seconds. The source code is available in online repository Ansari [2023]. We use Algorithm 1, along with the DBC method described in Section 5 with the setup that yields the best performance; see our comparative experiments in Table 6 in Appendix 5.

In order to define the locations in $L$ , we group Joplin “block IDs” in 100 center locations using the $k$ -means method [MacQueen, 1967] and aggregate the data accordingly. Particularly, we assume that when the tornado hits a location then all of its buildings are hit, and that the tornado hits a location only if the location’s centroid is within $\Delta$ units of the path. Table 2 compares the population dislocation for when the tornado length is $E=5$ miles or $\infty$ and the budgets are $A=0,15,30$ million USD. Columns 1 and 2 in Table 2 show $E$ and $A$ . Column 3 has the optimal population dislocation under the worst-case tornado for each set of parameters. Columns 4 and 5 report the CPU run time of C&CG Algorithm 1 in seconds and its number of iterations, respectively. Column 6 shows the CPU runtime of the DBC method in seconds to solve subproblem $\Phi(f)$ . Column 7 compares the CPU time that the DBC algorithm spends in a callback function to verify the feasibility of integral solutions. Column 8 shows the total CPU time to find a feasible solution in the set $\tilde{\mathcal{U}}^{h}$ . Columns 9 and 10 present the proportion of the total budget that is spent for retrofitting and recovery, respectively.

Table 2: Solving the two-stage robust optimization model with different parameters of tornado length and available budget, after 60 days of recovery for 100 locations group in Joplin, MO.

Length Budget Population CCG CCG DBC Callback $\tilde{\mathcal{U}}^{h}$ feasibility Retrofitting cost Recovery cost ($) dislocation run time (sec.) iteration run time (sec.) run time (sec.) run time (sec.) /budget /budget 5 0 M 16318 127 2 126 3 1 - - 15 M 13408 307 3 306 11 6 0.6 0.4 30 M 13091 300 3 300 11 5 0.3 0.7 $\infty$ 0 M 17293 568 2 568 3 0 - - 15 M 14236 716 2 716 5 0 0.6 0.4 30 M 13839 986 3 985 5 0 0.4 0.6

Table 2 shows that the solver is able to find optimal solutions for every experiment in at most 986 seconds. It also shows that solving the subproblems is the most time-consuming step (which is virtually 99% of the total runtime). In general, the algorithm solves the subproblem for 5 miles line segment faster than the (infinite) line; this fact can be justified by the effect of the conflict constraints (7) which are not present when $E=\infty$ . In addition, for all experiments, the callback function (i.e., Algorithm 2) is very fast, which implies that the bottleneck of the entire algorithm is the branch-and-bound that solves the relaxation of the subproblem. As expected, for the infinite line, the results show zero seconds to check the feasibility using $\tilde{\mathcal{U}}^{h}$ , because in this case, the SLA alone can verify the feasibility.

For both tornado lengths, there are decreases in the dislocation by increasing the budget because decision-makers are able to retrofit and recover more locations. For example, increasing the budget from $0 to $15M results in around 3000 fewer dislocations in damaged neighborhoods. Note that this translates in reductions of around 20% in population dislocation for both tornado lengths by investing $15M. However, adding another $15 million will reduce dislocations only by around 500 people because of the great cost of recovery. In fact, the relationship between budget and dislocation exhibits a clear diminishing returns behavior, see Figure 8, where it can be seen that after around $10M there is a small decrease in population dislocation. From a decision-making point of view, these results implies that it is not reasonable to spend more than $10M-$15M in retrofitting and recovery, as the return of investment is not worth it.

More details on the optimal strategies can be found in Figure 9, which maps the locations, their optimal strategies, and the worst-case tornado scenarios for various budgets. It shows that in all experiments the decision-maker uses most of the extra amount of budget from $15M to $30M to recover a few locations rather than in retrofit locations that are not in the path of the worst-case tornado. Notice that recovery costs are generally more expensive than retrofitting costs. For example, the recovered location in the $\infty$ -$15M map in Figure 9 has 470 residents; 382 of which are estimated to be dislocated by spending $194K for retrofitting. However, the decision-maker has to spend $5.5M in recovery to further reduce the estimated dislocation to 295 people for this location. Columns 9-10 in Table 2 show that the proportion of the budget that is assigned to recovery increases by having more budget. In addition, we can observe in Figure 9 that the path of worst-case tornadoes does not change by much in different experiments, suggesting that the location of the tornado path is largely driven by the population density in the city center area.

These results suggest, from a worst-case perspective, that it is optimal to invest only in locations in or close to where the worst-case tornado hits and that it is better to invest both in retrofitting and recovery, rather than just in retrofitting. The results are somewhat surprising in the sense that, at least a priori, it would be reasonable to invest all budget in retrofitting across a wider area (particularly because recovery is much more expensive than retrofitting). As such, the results suggest that there is a critical set of clustered locations, in this case the densely populated center of the city, where most of the investments must be concentrated. We suspect, however, that there is a different structure for the optimal actions in data sets with a more homogeneous population distribution.

6.3 Advantages of robust retrofitting

In this section we compare the optimal two-stage robust retrofitting and recovery strategies against a policy that makes retrofitting and recovery decisions at random. We assume a $E=5$ miles tornado path and a total budget of $A=$ $15 million. The decision-maker is allowed to spend $0, $3, $9, and $15 million of the budget in each experiment to randomly retrofit locations with random selections in set $S$ ; the remaining budget is reserved for the recovery actions. Once a retrofitting strategy $f$ is determined, then the corresponding worst-case second-stage happens (i.e., the value of $\Phi(f)$ is determined, see Equation (5)). Table 3 compares the average, maximum, minimum, and standard deviation of population dislocation for 10 random replications. Note that for the $0M case, we do not need multiple replications.

Table 3: Comparison of population dislocation statistics for 10 random retrofitting plans with spending a certain amount of budget out of a total $15M.

Budget to randomly retrofit Population Dislocation Average Maximum Minimum Standard deviation $0M (0%) 16318 - $3M (20%) 16057 16263 15717 163 $9M (60%) 15596 15905 15306 196 $15M (100%) 15143 15451 14764 189 Robust optimal retrofitting 13408 -

The results in Table 3 show that with random policies the average, maximum, and minimum dislocation after the corresponding worst-case tornado decrease with the percentage of the budget spent in retrofitting. These results are depicted graphically in the Box plots in Figure 10. Note that the proposed robust optimization model results in 13408 dislocations, which is noticeable better than the random policies. To further emphasize the value of optimization, recall that the optimal solution with our proposed method in Table 2 spends $9M out of $15M (60% of budget) in retrofitting. By contrast, even the minimum population dislocation with random retrofitting using the same budget proportion is 15306 (a 14% increase).

We conduct additional simulations to randomly retrofit locations with the full amount of $15M and $30M budget as before, but this time the worst-case tornado found in Section 6.2 is used. The results across 10 replications are shown in Table 4, which shows that for both cases the robust optimal population dislocation is significantly less than the minimum of population dislocation with random retrofitting. The plot in Figure 11 provides a visualization of the significant gap between robust optimal value (stars) and the distribution of values for random retrofitting plans. We can also observe that population dislocation reduces by increasing the amount of budget.

Table 4: Comparison of population dislocation statistics using ten different random retrofitting plans and budgets of $15M and $30M.

Budget Population dislocation Robust optimal retrofitting Average Maximum Minimum Standard deviation $15M 13408 15314 15941 14545 331 $30M 13091 14795 15248 14145 294

The results of this section show that deciding using the proposed two-stage robust method yields significant reductions in population dislocation against random policies when evaluating in terms of worst-case outcomes. In the next section we complement these experiments by analyzing how the optimal policy behaves when random tornadoes, instead of the worst-case scenario, happen.

6.4 Robust strategies in random scenarios

In this section we fix the optimal retrofitting-recovery plan and compare the robust population dislocation values resulting from having 10 simulated 5-miles tornadoes, see Table 5 and Figure 12.

Table 5: Comparison of population dislocation statistics for the worst-case tornado and simulated tornadoes by using the optimal retrofitting plan

Budget Population dislocation Worst-case tornado Average Maximum Minimum Standard deviation $0M 16318 10324 13201 5233 2835 $15M 13408 9634 12720 5031 2679 $30M 13091 9416 12480 4900 2607

The results for $0M in Table 5 and Figure 12 show that there is a large gap between the worst-case tornado 16318 (blue star) and the maximum value among simulation results which is 13201. As expected, by having more budget, the average population dislocations in the simulations reduce from 10324 to 9643 for $15M and to 9416 for $30M. In addition, the gap between the robust and maximum simulation values for $15 and $30M are noticeable smaller than the gap for the no-budget case.

These observations show the following about the optimal retrofitting-recovery policy: first, that there can be a large variability for population dislocation assuming any tornado is possible. Second, that the variability seems to decrease as more budget is available, and third, that the worst-case scenario is not that far way from simulated scenarios. This last point is further emphasized by the fact that only 10 tornadoes are simulated; one should expect that with more simulated tornadoes the gap between the worst-case and the maximum dislocation observed reduces even more. Further, this last point also shows that the worst-case is not as rare as one might anticipate, and thus that the proposed robust optimization model does not suffer from ‘over-conservatism,’ which is a common drawback of deciding in a robust manner.

7 Conclusions

Motivated by the need of data-driven methods to allocate resources in retrofitting for tornado-prone areas, we consider a two-stage robust optimization model to determine retrofitting and recovery strategies before and after the realization of an uncertain tornado disaster. We assume that the decision-maker (a combination of public-private agencies) has to allocate a budget in retrofitting and recovery actions, to minimize population dislocation. We propose a two-stage optimization model where the first-stage variables determine the retrofitting strategies and the second-stage variables determine the recovery strategy, assuming that for any retrofitting action, the corresponding worst-case tornado happens. Given that the uncertainty corresponds to the possible tornado paths, the resulting problem is a two-stage mixed-integer robust optimization problem with a mixed-integer non-linear uncertainty set.

We show that the proposed problem is NP-hard. A standard approach to solve the problem is a column-and-constraint generation algorithm that needs to solve a max-min subproblem at each iteration. We embed a decomposition branch-and-cut algorithm in the column-and-constraint solution method to address the non-convex structure of the bilevel subproblem. We propose two sets of conflict constraints and a custom separation procedure to implement the decomposition branch-and-cut algorithm. Using data from INCORE, we use our model to provide the optimal retrofitting and recovery strategies for the city of Joplin, Missouri against the worst-case scenario. The results show that the optimal policy exhibits a diminishing returns behavior that sets-in quickly; prioritizes the high-density city center; and prefers to retrofit and recover than to only retrofit. The proposed policy also outperforms other benchmark and does not suffer from ‘over-conservatism.’

Future work will focus on modeling tornado damages assuming that there might be several levels of disruptions. Additionally, future research would analyze other social vulnerability indices to measure community resilience under the risk of a tornado. One line of research might study the population dislocation in residential areas with different income levels after certain days of a tornado occurrence, in order to recognize more vulnerable locations that require immediate agencies’ attention. Another line of research might take our proposed framework to study other natural disasters like hurricanes with minor modifications in the uncertainty model and the set of retrofitting strategies.

Acknowledgment

This research is partially funded by the National Science Foundation (NSF) (Awards ENG/CMMI # 2145553 and GEO/RISE # 2052930) and by the Air Force Office of Scientific Research (AFORS) (Award # FA9550-22-1-0236). This research was also partially funded by the National Institute of Standards and Technology (NIST) Center of Excellence for Risk-Based Community Resilience Planning through a cooperative agreement with Colorado State University [Awards # 70NANB20H008 and 70NANB15H044]. The contents expressed in this paper are the views of the authors and do not necessarily represent the opinions or views of the National Institute of Standards and Technology (NIST), the National Science Foundation (NSF), or the Air Force Office of Scientific Research (AFORS).

References

An and Zeng [2015] Y. An and B. Zeng. Exploring the modeling capacity of two-stage robust optimization: Variants of robust unit commitment model. IEEE Transactions on Power Systems, 30(1):109–122, 2015. doi: 10.1109/TPWRS.2014.2320880.
An et al. [2014] Y. An, B. Zeng, Y. Zhang, and L. Zhao. Reliable p-median facility location problem: two-stage robust models and algorithms. Transportation Research Part B: Methodological, 64:54–72, 2014.
Ansari [2023] M. Ansari. Github repository. Online at https://github.com/mehdi-ansari/Two_Stage_Robust_Tornado_Problem, Accessed September 2023, September 2023.
Ansari et al. [2022] M. Ansari, J. S. Borrero, and L. Lozano. Robust minimum-cost flow problems under multiple ripple effect disruptions. INFORMS Journal on Computing, 2022.
Atamtürk and Zhang [2007] A. Atamtürk and M. Zhang. Two-stage robust network flow and design under demand uncertainty. Operations Research, 55(4):662–673, 2007.
Ben-Tal and Nemirovski [1999] A. Ben-Tal and A. Nemirovski. Robust solutions to uncertain linear programs. Operations Research Letters, 25:1–13, 1999.
Ben-Tal et al. [2004] A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski. Adjustable robust solutions of uncertain linear programs. Mathematical programming, 99(2):351–376, 2004.
Ben-Tal et al. [2011] A. Ben-Tal, B. Do Chung, S. R. Mandala, and T. Yao. Robust optimization for emergency logistics planning: Risk mitigation in humanitarian relief supply chains. Transportation research part B: methodological, 45(8):1177–1189, 2011.
Ben-Tal et al. [2015] A. Ben-Tal, E. Hazan, T. Koren, and S. Mannor. Oracle-based robust optimization via online learning. Operations Research, 63(3):628–638, 2015.
Bertsimas and Sim [2004] D. Bertsimas and M. Sim. The price of robustness. Operations Research, 52(1):35–53, 2004.
Borrero and Lozano [2021] J. S. Borrero and L. Lozano. Modeling defender-attacker problems as robust linear programs with mixed-integer uncertainty sets. INFORMS Journal on Computing, 33(4):1570–1589, 2021.
Cheng et al. [2021] C. Cheng, Y. Adulyasak, and L.-M. Rousseau. Robust facility location under demand uncertainty and facility disruptions. Omega, 103:102429, 2021.
Ding et al. [2017] T. Ding, C. Li, Y. Yang, J. Jiang, Z. Bie, and F. Blaabjerg. A two-stage robust optimization for centralized-optimal dispatch of photovoltaic inverters in active distribution networks. IEEE Transactions on Sustainable Energy, 8(2):744–754, 2017. doi: 10.1109/TSTE.2016.2605926.
Fan and Pang [2019] F. Fan and W. Pang. Stochastic track model for tornado risk assessment in the us. Frontiers in built environment, 5:37, 2019.
Farokhnia et al. [2020] K. Farokhnia, J. W. van de Lindt, and M. Koliou. Selection of residential building design requirements to achieve community functionality goals under tornado loading. Practice Periodical on Structural Design and Construction, 25(1):04019035, 2020.
FEMA [2003] FEMA. Multi-hazard loss estimation methodology: Earthquake model. Hazus MH 2.1. 2003.
FEMA [2021] FEMA. Natural hazard retrofit program toolkit. https://www.fema.gov/sites/default/files/documents/fema_natural-hazards-retrofit-program-tookit.pdf, Accessed January 2023, February 2021.
FEMA [2022a] FEMA. Tornado–alerts and warnings. Online at https://community.fema.gov/ProtectiveActions, Accessed November 2022, November 2022a.
FEMA [2022b] FEMA. Resources for repairing, retrofitting and rebuilding after a tornado. https://www.fema.gov/fact-sheet/resources-repairing-retrofitting-and-rebuilding-after-tornado-0, Accessed January 2023, March 2022b.
Gabrel et al. [2014] V. Gabrel, M. Lacroix, C. Murat, and N. Remli. Robust location transportation problems under uncertain demands. Discrete Applied Mathematics, 164:100–111, 2014.
Ho-Nguyen and Kılınç-Karzan [2018] N. Ho-Nguyen and F. Kılınç-Karzan. Online first-order framework for robust convex optimization. Operations Research, 66(6):1670–1692, 2018.
IN-CORE [2022] IN-CORE. Interdependent networked community resilience modeling environment. https://incore.ncsa.illinois.edu/, Accessed October 2022, October 2022.
Jabr et al. [2015] R. A. Jabr, I. Džafić, and B. C. Pal. Robust optimization of storage investment on transmission networks. IEEE Transactions on Power Systems, 30(1):531–539, 2015. doi: 10.1109/TPWRS.2014.2326557.
Jiang et al. [2012] R. Jiang, M. Zhang, G. Li, and Y. Guan. Benders’ decomposition for the two-stage security constrained robust unit commitment problem. In IIE Annual Conference. Proceedings, page 1. Institute of Industrial and Systems Engineers (IISE), 2012.
Koliou and van de Lindt [2020] M. Koliou and J. W. van de Lindt. Development of building restoration functions for use in community recovery planning to tornadoes. Natural Hazards Review, 21(2):04020004, 2020.
Kuligowski et al. [2014] E. D. Kuligowski, F. T. Lombardo, L. Phan, M. L. Levitan, and D. P. Jorgensen. Final report, national institute of standards and technology (nist) technical investigation of the may 22, 2011, tornado in joplin, missouri. 2014.
Laporte and Louveaux [1993] G. Laporte and F. V. Louveaux. The integer L-shaped method for stochastic integer programs with complete recourse. Operations Research Letters, 13:133–142, 1993.
Li et al. [2021] Y. Li, F. Zhang, Y. Li, and Y. Wang. An improved two-stage robust optimization model for cchp-p2g microgrid system considering multi-energy operation under wind power outputs uncertainties. Energy, 223:120048, 2021.
Lin et al. [2004] X. Lin, S. L. Janak, and C. A. Floudas. A new robust optimization approach for scheduling under uncertainty: I. bounded uncertainty. Computers & chemical engineering, 28(6-7):1069–1085, 2004.
Ma et al. [2018] S. Ma, B. Chen, and Z. Wang. Resilience enhancement strategy for distribution systems under extreme weather events. IEEE Transactions on Smart Grid, 9(2):1442–1451, 2018.
MacQueen [1967] J. MacQueen. Classification and analysis of multivariate observations. In 5th Berkeley Symp. Math. Statist. Probability, pages 281–297, 1967.
Masoomi and van de Lindt [2018] H. Masoomi and J. W. van de Lindt. Restoration and functionality assessment of a community subjected to tornado hazard. Structure and Infrastructure Engineering, 14(3):275–291, 2018.
Masoomi et al. [2018] H. Masoomi, M. R. Ameri, and J. W. van de Lindt. Wind performance enhancement strategies for residential wood-frame buildings. Journal of Performance of Constructed Facilities, 32(3):04018024, 2018.
Matthews et al. [2019] L. R. Matthews, C. E. Gounaris, and I. G. Kevrekidis. Designing networks with resiliency to edge failures using two-stage robust optimization. European Journal of Operational Research, 279(3):704–720, 2019.
McAllister et al. [2015] T. P. McAllister et al. Community resilience planning guide for buildings and infrastructure systems, volume i. 2015.
Mutapcic and Boyd [2009] A. Mutapcic and S. Boyd. Cutting-set methods for robust convex optimization with pessimizing oracles. Optimization Methods & Software, 24(3):381–406, 2009.
Neyshabouri and Berg [2017] S. Neyshabouri and B. P. Berg. Two-stage robust optimization approach to elective surgery and downstream capacity planning. European Journal of Operational Research, 260(1):21–40, 2017.
NOAA [2022a] NOAA. National centers for environmental information, tornadoes statistics. Online at https://www.ncei.noaa.gov, Accessed September 2022, September 2022a.
NOAA [2022b] NOAA. Storm predication center - national weather service. Online at https://www.spc.noaa.gov/, Accessed October 2022, October 2022b.
NOAA [2022c] NOAA. Tornadoes frequently asked questions. Online at https://www.weather.gov/lmk/tornadoesfaq, Accessed November 2022, November 2022c.
NOAA [2022d] NOAA. Supercells. Online at https://www.spc.noaa.gov/misc/AbtDerechos/supercells.htm, Accessed November 2022, November 2022d.
NOAA [2022e] NOAA. Weather related fatality and injury statistics. Online at https://www.weather.gov/hazstat/, Accessed November 2022, November 2022e.
NOAA [2022f] NOAA. National weather service, the enhanced fujita scale (ef scale). https://www.weather.gov, Accessed September 2022, September 2022f.
NOAA [2023] NOAA. Damage assessment toolkit. Online at https://apps.dat.noaa.gov/stormdamage/damageviewer/, Accessed March 2023, March 2023.
NWC [2023] NWC. Historical data of oklahoma county tornadoes. Online at https://www.weather.gov/oun/tornadodata-county-ok-oklahoma, Accessed January 2023, January 2023.
Ripberger et al. [2018] J. T. Ripberger, H. C. Jenkins-Smith, C. L. Silva, J. Czajkowski, H. Kunreuther, and K. M. Simmons. Tornado damage mitigation: Homeowner support for enhanced building codes in oklahoma. Risk analysis, 38(11):2300–2317, 2018.
Saba [2013] S. Saba. Line intersecting maximal number of circles (circle “stabbing” problem). https://sahandsaba.com/line-intersecting-maximal-number-of-circles.html, Accessed October 2022, October 2013.
Shams et al. [2021] M. H. Shams, M. Shahabi, M. MansourLakouraj, M. Shafie-khah, and J. P. Catalão. Adjustable robust optimization approach for two-stage operation of energy hub-based microgrids. Energy, 222:119894, 2021.
Simmons et al. [2015] K. M. Simmons, P. Kovacs, and G. A. Kopp. Tornado damage mitigation: Benefit–cost analysis of enhanced building codes in oklahoma. Weather, climate, and society, 7(2):169–178, 2015.
Standohar-Alfano et al. [2017] C. D. Standohar-Alfano, J. W. van de Lindt, and B. R. Ellingwood. Vertical load path failure risk analysis of residential wood-frame construction in tornadoes. Journal of Structural Engineering, 143(7):04017045, 2017.
Stoner and Pang [2021] M. Stoner and W. Pang. Tornado hazard assessment of residential structures built using cross-laminated timber and light-frame wood construction in the us. Natural Hazards Review, 22(4):04021032, 2021.
Strader et al. [2016] S. M. Strader, T. J. Pingel, and W. S. Ashley. A monte carlo model for estimating tornado impacts. Meteorological Applications, 23(2):269–281, 2016.
Takeda et al. [2008] A. Takeda, S. Taguchi, and R. Tütüncü. Adjustable robust optimization models for a nonlinear two-period system. Journal of Optimization Theory and Applications, 136(2):275–295, 2008.
The White House [2022] The White House. Fact sheet: Biden–harris administration launches initiative to modernize building codes, improve climate resilience, and reduce energy costs. https://www.whitehouse.gov/briefing-room/statements-releases/2022/06/01/fact-sheet-biden-harris-administration-launches-initiative-to-modernize-building-codes-improve-climate-resilience-and-reduce-energy-costs/, June 2022.
Thiele et al. [2009] A. Thiele, T. Terry, and M. Epelman. Robust linear optimization with recourse. Rapport technique, pages 4–37, 2009.
van de Lindt and Dao [2009] J. W. van de Lindt and T. N. Dao. Performance-based wind engineering for wood-frame buildings. Journal of Structural Engineering, 135(2):169–177, 2009.
Velasquez et al. [2020] G. A. Velasquez, M. E. Mayorga, and O. Y. Özaltın. Prepositioning disaster relief supplies using robust optimization. IISE Transactions, 52(10):1122–1140, 2020.
Wang et al. [2017] J. Wang, S. Cao, W. Pang, and J. Cao. Experimental study on effects of ground roughness on flow characteristics of tornado-like vortices. Boundary-Layer Meteorology, 162:319–339, 2017.
Wang et al. [2021] W. Wang, J. W. Van De Lindt, N. Rosenheim, H. Cutler, B. Hartman, J. Sung Lee, and D. Calderon. Effect of residential building wind retrofits on social and economic community-level resilience metrics. Journal of Infrastructure Systems, 27(4):04021034, 2021.
Wen [2021] Y. Wen. Development of Multi-Objective Optimization Model of Community Resilience on Mitigation Planning. PhD thesis, University of Oklahoma, 2021.
Xiong et al. [2017] P. Xiong, P. Jirutitijaroen, and C. Singh. A distributionally robust optimization model for unit commitment considering uncertain wind power generation. IEEE Transactions on Power Systems, 32(1):39–49, 2017.
Yao et al. [2009] T. Yao, S. R. Mandala, and B. Do Chung. Evacuation transportation planning under uncertainty: a robust optimization approach. Networks and Spatial Economics, 9(2):171, 2009.
Yuan et al. [2016] W. Yuan, J. Wang, F. Qiu, C. Chen, C. Kang, and B. Zeng. Robust optimization-based resilient distribution network planning against natural disasters. IEEE Transactions on Smart Grid, 7(6):2817–2826, 2016.
Zeng and Zhao [2013] B. Zeng and L. Zhao. Solving two-stage robust optimization problems using a column-and-constraint generation method. Operations Research Letters, 41(5):457–461, 2013.
Zhao and Zeng [2012] L. Zhao and B. Zeng. Robust unit commitment problem with demand response and wind energy. In 2012 IEEE power and energy society general meeting, pages 1–8. IEEE, 2012.

APPENDIX

1 Proof of $NP$ -hardness

Proposition 3.

The optimization problem (1) is NP-hard.

Proof: Consider an arbitrary instance of the knapsack problem with $N$ items, weights $b_{i},\forall i\in N$ , values $q_{i},\forall i\in N$ , capacity $B$ , and value $Q$ . We build an instance of problem (1) that is equivalent to the knapsack problem. Consider an instance of problem (1) with $|S|=1$ and $|P|=2$ . Observe that because $|S|=1$ , $f$ has to be equal to a vector of ones. Suppose that the coordinates of the locations and the values of $E$ and $\Delta$ are such that $z_{\ell}=1$ for all $\ell\in L$ is feasible in $\mathcal{U}$ . For simplicity, we drop index $s$ from the notation. Also, note that $r_{\ell 2}$ can be replaced by $1-r_{\ell 1}$ for all $\ell\in L$ because only two recovery strategies are assumed. The specific case of the problem can be rewritten as


$\displaystyle v=\min$	$\displaystyle\sum_{\ell\in L}(g_{\ell 1}-g_{\ell 2})r_{\ell 1}+\sum_{\ell\in L}(w_{\ell}+g_{\ell 2})$	(19a)
s.t.	$\displaystyle\sum_{\ell\in L}(c_{\ell 1}-c_{\ell 2})r_{\ell 1}\leq A-\sum_{\ell\in L}(d_{\ell}+c_{\ell 2})$	(19b)
	$\displaystyle r_{\ell 1}\in\{0,1\}^{\|L\|}.$	(19c)

Consider the example of problem (19) where $L=N$ , $c_{\ell 1}-c_{\ell 2}=b_{\ell}$ and $g_{\ell 2}-g_{\ell 1}=q_{\ell}$ for each item $\ell\in L$ , and $A-\sum_{\ell\in L}(d_{\ell}+c_{\ell 2})=B$ . The answer of knapsack problem for some instance $L^{*}\subseteq N$ is YES (i.e., $\sum_{i\in L^{*}}q_{i}\geq Q$ ), if and only if ${v}\leq\sum_{\ell\in L}(w_{\ell}+g_{\ell 2})-Q$ .

2 Proof of Proposition 1

Proof: Suppose $(\ell_{1},\ell_{2})\in\Omega$ , and for the sake of contradiction that $z_{\ell_{1}}=z_{\ell_{2}}=1$ . For the constraint in (3b), we have


	$\displaystyle\\|(1-t_{\ell_{1}})e_{0}+t_{\ell_{1}}e_{1}-(x_{\ell_{1}},y_{\ell_{1}})\\|\leq\Delta$		(20a)
	$\displaystyle\\|(1-t_{\ell_{2}})e_{0}+t_{\ell_{2}}e_{1}-(x_{\ell_{2}},y_{\ell_{2}})\\|\leq\Delta.$		(20b)

The sum of two above inequalities results

\|(1-t_{\ell_{1}})e_{0}+t_{\ell_{1}}e_{1}-(x_{\ell_{1}},y_{\ell_{1}})\|+\|(x_{\ell_{2}},y_{\ell_{2}})-(1-t_{\ell_{2}})e_{0}-t_{\ell_{2}}e_{1}\|\leq 2\Delta,

(21)

and by applying the triangle inequality $\|X+Y\|\leq\|X\|+\|Y\|$ , we conclude that

\|(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})+(t_{\ell_{2}}-t_{\ell_{1}})(e_{0}-e_{1})\|\leq 2\Delta.

(22)

Finally, we employ the fact that $\|X\|-\|Y\|\leq\|X-Y\|$ to reach

\|(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})\|-\|(t_{\ell_{1}}-t_{\ell_{2}})(e_{0}-e_{1})\|\leq 2\Delta.

(23)

Now, we show (23) cannot be true because $(\ell_{1},\ell_{2})\in\Omega$ , i.e.,

\|(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})\|>2\Delta+E.

(24)

By constraint (3d) and using $\|aY\|=|a|\|Y\|$ , we have

\|(t_{\ell_{1}}-t_{\ell_{2}})(e_{0}-e_{1})\|=|t_{\ell_{1}}-t_{\ell_{2}}|\|e_{0}-e_{1}\|\leq|t_{1}-t_{2}|E\leq E.

(25)

Note that $|t_{1}-t_{2}|\leq 1$ . So, the two above inequalities imply that

		$\displaystyle\begin{cases}\\|(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})\\|>2\Delta+E\\ -\\|(t_{1}-t_{2})(e_{0}-e_{1})\\|_{1}\geq-E\end{cases}$		(26)
		$\displaystyle\quad\implies\\|(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})\\|-\\|(t_{\ell_{1}}-t_{\ell_{2}})(e_{0}-e_{1})\\|>2\Delta.$		(26)

Observe that (23) violates the true inequality (LABEL:contradiction_res). Therefore, the cut $z_{\ell_{1}}+z_{\ell_{2}}\leq 1$ is valid because $z_{\ell_{1}}=z_{\ell_{2}}=1$ is an infeasible solution under the assumption $(\ell_{1},\ell_{2})\in\Omega$ .

3 Proof of Proposition 2

This section provides a more detailed discussion (and a more detailed proof) for Proposition 2. Proof: We are interested to identify the division of space $R$ into feasible ( $R_{\ell_{1},\ell_{2}}^{\prime}$ ) and infeasible ( $R_{\ell_{1},\ell_{2}}^{\prime\prime}$ ) subspaces for two locations $\ell_{1}$ and $\ell_{2}$ with coordinates $(x_{\ell_{1}},y_{\ell_{1}})$ and $(x_{\ell_{2}},y_{\ell_{2}})$ , respectively. For the sake of simplicity, we assume that $x_{\ell_{1}}\leq x_{\ell_{2}}$ throughout the proof; the results can be easily extended to the case $x_{\ell_{1}}>x_{\ell_{2}}$ .

Lemma 1.

Lines $\lambda_{1},\ldots,\lambda_{6}$ in (10a)-(10f) determine boundaries of subspaces.

Proof: We find feasible slope $m$ and y-intercept $q$ of lines in the form of $y=mx+q$ with a distance $\Delta$ from both locations $\ell_{1}$ and $\ell_{2}$ which means the following equalities must simultaneously hold,


	$\displaystyle\frac{\|y_{\ell_{1}}-mx_{\ell_{1}}-q\|}{\sqrt{1+m^{2}}}=\Delta,$		(27a)
	$\displaystyle\frac{\|y_{\ell_{2}}-mx_{\ell_{2}}-q\|}{\sqrt{1+m^{2}}}=\Delta.$		(27b)

The above equalities imply two possible cases:

	$\displaystyle\frac{\|y_{\ell_{1}}-mx_{\ell_{1}}-q\|}{\sqrt{1+m^{2}}}$	$\displaystyle=\frac{\|y_{\ell_{2}}-mx_{\ell_{2}}-q\|}{\sqrt{1+m^{2}}}$		(28)
	$\displaystyle\implies$	$\displaystyle\begin{cases}y_{\ell_{1}}-mx_{\ell_{1}}=y_{\ell_{2}}-mx_{\ell_{2}}\ (\textbf{Case 1})\\ y_{\ell_{1}}-mx_{\ell_{1}}-q=mx_{\ell_{2}}+q-y_{\ell_{2}}\ (\textbf{Case 2})\end{cases}.$		(28)

Case 1.

m=\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}=\tan(\theta)\ (\text{by the definition}).

(29)

By replacing $m=\tan(\theta)$ in (27a), we will have


	$\displaystyle\frac{\|y_{\ell_{1}}-\tan(\theta)x_{\ell_{1}}-q\|}{\sqrt{1+\tan^{2}(\theta)}}=\Delta$		(30a)
$\displaystyle\implies$	$\displaystyle\|y_{\ell_{1}}-\tan(\theta)x_{\ell_{1}}-q\|\|\cos(\theta)\|=\Delta$	$\displaystyle(\sqrt{1+\tan^{2}(\theta)}=\frac{1}{\|\cos(\theta)\|})$	(30b)
$\displaystyle\implies$	$\displaystyle\|y_{\ell_{1}}-\tan(\theta)x_{\ell_{1}}-q\|=\frac{\Delta}{\cos(\theta)}$	$\displaystyle(x_{\ell_{1}}\leq x_{\ell_{2}}\implies\cos(\theta)\geq 0)$	(30c)
$\displaystyle\implies$	$\displaystyle\begin{cases}q^{\prime}=y_{\ell_{1}}-\tan(\theta)x_{\ell_{1}}-\frac{\Delta}{\cos(\theta)}\\ q^{\prime\prime}=y_{\ell_{1}}-\tan(\theta)x_{\ell_{1}}+\frac{\Delta}{\cos(\theta)}\end{cases}$		(30d)

So, the tornado paths with formulations


	$\displaystyle y=\tan(\theta)(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{\Delta}{\cos(\theta)},$		(31a)
	$\displaystyle y=\tan(\theta)(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{\Delta}{\cos(\theta)},$		(31b)

are two boundaries that cover both locations (observe that any other $q$ between $q^{\prime}$ and $q^{\prime\prime}$ ( $q^{\prime}\leq q\leq q^{\prime\prime}$ ) corresponding to the slope $m=\tan\theta$ forms a tornado path that has a distance at most $\Delta$ to the locations $\ell_{1}$ and $\ell_{2}$ ). The tornado in (31a) covers all locations in distance $\Delta$ between the lines $y=\tan\theta(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos\theta}$ (see (10b) for line $\lambda_{2}$ ) and $y=\tan\theta(x-x_{\ell_{1}})+y_{\ell_{1}}$ . Similarly, the tornado (31b) covers all points between $y=\tan\theta(x-x_{\ell_{1}})+y_{\ell_{1}}$ and $y=\tan\theta(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos\theta}$ (see (10a) for line $\lambda_{1}$ ).

Replacing $m=\tan\theta$ in (27b) will give the same boundaries as (27a).
Case 2.

q=\frac{y_{\ell_{1}}+y_{\ell_{2}}-m(x_{\ell_{1}}+x_{\ell_{2}})}{2}.

(32)

By replacing $q=\frac{y_{\ell_{1}}+y_{\ell_{2}}-m(x_{\ell_{1}}+x_{\ell_{2}})}{2}$ in (27a), we have

		$\displaystyle\frac{\|y_{\ell_{1}}-mx_{\ell_{1}}-\frac{y_{\ell_{1}}+y_{\ell_{2}}-m(x_{\ell_{1}}+x_{\ell_{2}})}{2}\|}{\sqrt{1+m^{2}}}=\Delta$		(33)
	$\displaystyle\implies$	$\displaystyle\frac{\|y_{\ell_{2}}-y_{\ell_{1}}-m(x_{\ell_{2}}-x_{\ell_{1}})\|}{\sqrt{1+m^{2}}}=2\Delta$		(33)

Instead of solving (33) for $m$ which leads to a complicated quadratic equation, we exploit the form of equality to compute possible values of $m$ . Observe that (33) can be interpreted as a line passing through the origin and has a distance $2\Delta$ from the point $(x_{\ell_{2}}-x_{\ell_{1}},y_{\ell_{2}}-y_{\ell_{1}})$ (see Figure 13).

Now, we discuss that the angles corresponding to the possible slopes of $m$ must be

	$\displaystyle m^{\prime}=\tan{\left(\arctan{\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}}+\arcsin{\frac{2\Delta}{d}}\right)},$		(34)
	$\displaystyle m^{\prime\prime}=\tan{\left(\arctan{\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}}-\arcsin{\frac{2\Delta}{d}}\right)},$		(35)

where $\arctan{\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}}=\theta$ and $\arcsin{\frac{2\Delta}{D(\ell_{1},\ell_{2})}}=\alpha$ by the definition. Therefore,

	$\displaystyle m^{\prime}=\tan{(\theta+\alpha)},$		(36)
	$\displaystyle m^{\prime\prime}=\tan{(\theta-\alpha)}.$		(37)

So, the resulting tornado paths are

	$\displaystyle y$	$\displaystyle=\tan(\theta+\alpha)x+\frac{y_{\ell_{1}}+y_{\ell_{2}}-\tan(\theta+\alpha)(x_{\ell_{1}}+x_{\ell_{2}})}{2}$		(38)
		$\displaystyle=\tan(\theta+\alpha)(x-\frac{x_{\ell_{1}}+x_{\ell_{2}}}{2})+\frac{y_{\ell_{1}}+y_{\ell_{2}}}{2}$		(38)

	$\displaystyle y$	$\displaystyle=\tan(\theta-\alpha)x+\frac{y_{\ell_{1}}+y_{\ell_{2}}-\tan(\theta-\alpha)(x_{\ell_{1}}+x_{\ell_{2}})}{2}$		(39)
		$\displaystyle=\tan(\theta-\alpha)(x-\frac{x_{\ell_{1}}+x_{\ell_{2}}}{2})+\frac{y_{\ell_{1}}+y_{\ell_{2}}}{2}$		(39)

Observe that these tornado paths pass through the midpoint $(\frac{x_{\ell_{1}}+x_{\ell_{2}}}{2},\frac{y_{\ell_{1}}+y_{\ell_{2}}}{2})$ . So, two edges of the tornado path in (38) with distance $\Delta$ will be $y=\tan(\theta+\alpha)(x-x_{\ell_{1}})+y_{\ell_{1}}$ and $y=\tan(\theta+\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}$ (see lines $\lambda_{3}$ and $\lambda_{4}$ in (10c) and (10d)). Likewise, the $\Delta$ edges of (39) are $y=\tan(\theta-\alpha)(x-x_{\ell_{1}})+y_{\ell_{1}}$ and $y=\tan(\theta-\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}$ (see lines $\lambda_{5}$ and $\lambda_{6}$ in (10e) and (10f)). Replacing $q$ in (27b) leads to the same equality in (33), and so the same conclusions.

In the following, we determine the feasible pair $(\beta,q(\beta))$ for the line $y=\tan(\beta)x+q(\beta)$ to cover both locations $\ell_{1}$ and $\ell_{2}$ within the distance $\Delta$ , where $\beta=\arctan(m)$ is the angle of the slope of the line and $q(\beta)$ is the $y$ -intercept of the line given angle $\beta$ . It can be shown that any feasible slope angle $\beta$ must be inside the interval $[\theta-\alpha,\theta+\alpha]$ as follows.


	$\displaystyle\frac{\|y_{\ell_{1}}-\tan(\beta)x_{\ell_{1}}-q(\beta)\|}{\sqrt{1+\tan^{2}\beta}}\leq\Delta,$		(40a)
	$\displaystyle\frac{\|y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}-q(\beta)\|}{\sqrt{1+\tan^{2}\beta}}\leq\Delta.$		(40b)

The summation of (40a) and (40b) results in

\frac{|-y_{\ell_{1}}+\tan(\beta)x_{\ell_{1}}+q(\beta)|+|y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}-q(\beta)|}{\sqrt{1+\tan^{2}\beta}}\leq 2\Delta

(41)

By the triangle inequality $\|X+Y\|\leq\|X\|+\|Y\|$ , we conclude

\frac{|y_{\ell_{2}}-y_{\ell_{1}}-\tan(\beta)(x_{\ell_{2}}-x_{\ell_{1}})|}{\sqrt{1+\tan^{2}\beta}}\leq 2\Delta.

(42)

Based on Case 2 in the proof of Lemma 1, (42) implies that $\beta\in[\theta-\alpha,\theta+\alpha]$ .

Now, we compute the feasible $q(\beta)$ for given $\beta\in[\theta-\alpha,\theta+\alpha]$ . From (40a) and (40b), we have


	$\displaystyle y_{\ell_{1}}-\tan(\beta)x_{\ell_{1}}-\frac{\Delta}{\|\cos(\beta)\|}\leq q(\beta)\leq y_{\ell_{1}}-\tan(\beta)x_{\ell_{1}}+\frac{\Delta}{\|\cos(\beta)\|}$		(43a)
	$\displaystyle y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}-\frac{\Delta}{\|\cos(\beta)\|}\leq q(\beta)\leq y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}+\frac{\Delta}{\|\cos(\beta)\|}.$		(43b)

Note that $\sqrt{1+\tan^{2}\beta}=1/{|\cos(\beta)|}$ . For the sake of simplicity, we assume hereafter that $\theta+\alpha<\pi/2$ and $\theta-\alpha>-\pi/2$ (which implies $\cos(\beta)>0$ ) and then discuss two cases for a feasible $q(\beta)$ followed by the inequalities (43a) and (43b):

$\theta\leq\beta\leq\theta+\alpha\implies\tan(\theta)\leq\tan(\beta)\implies\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}\leq\tan(\beta)$ :

y_{\ell_{1}}-\tan(\beta)x_{\ell_{1}}-\frac{\Delta}{\cos(\beta)}\leq q(\beta)\leq y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}+\frac{\Delta}{\cos(\beta)}.

(44)

$\theta-\alpha\leq\beta\leq\theta\implies\tan(\theta)\geq\tan(\beta)\implies\frac{y_{\ell_{2}}-y_{\ell_{1}}}{x_{\ell_{2}}-x_{\ell_{1}}}\geq\tan(\beta)$ :

y_{\ell_{2}}-\tan(\beta)x_{\ell_{2}}-\frac{\Delta}{\cos(\beta)}\leq q(\beta)\leq y_{\ell_{1}}-\tan(\beta)x_{\ell_{1}}+\frac{\Delta}{\cos(\beta)}.

(45)

Figure 14 visualizes an example of the feasible ranges (44) and (45) for $q(\beta)$ given angle $\beta\in[\theta-\alpha,\theta+\alpha]$ .

Inequalities (44) and (45) imply line $y=\tan(\beta)x+q(\beta)$ for given $\beta\in[\theta-\alpha,\theta+\alpha]$ is feasible (covers $\ell_{1}$ and $\ell_{2}$ within a distance of $\Delta$ ), if and only if a point $(\bar{x},\bar{y})$ on the line meets

\begin{cases}\beta\in[\theta,\theta+\alpha]:\bar{y}\in[\tan(\beta)(\bar{x}-x_{\ell_{1}})+y_{\ell_{1}}-\frac{\Delta}{\cos(\beta)},\tan(\beta)(\bar{x}-x_{\ell_{2}})+y_{\ell_{2}}+\frac{\Delta}{\cos(\beta)}],\\ \beta\in[\theta-\alpha,\theta]:\bar{y}\in[\tan(\beta)(\bar{x}-x_{\ell_{2}})+y_{\ell_{2}}-\frac{\Delta}{\cos(\beta)},\tan(\beta)(\bar{x}-x_{\ell_{1}})+y_{\ell_{1}}+\frac{\Delta}{\cos(\beta)}].\end{cases}

(46)

Remark 4.

Given the angle $\beta\in[\theta-\alpha,\theta+\alpha]$ , a feasible line $y=\tan(\beta)x+q(\beta)$ covers a point $(x_{0},y_{0})$ within the distance $\Delta$ , if and only if

\begin{cases}\beta\in[\theta,\theta+\alpha]:y_{0}\in[\tan(\beta)(x_{0}-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\beta)},\tan(\beta)(x_{0}-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\beta)}],\\ \beta\in[\theta-\alpha,\theta]:y_{0}\in[\tan(\beta)(x_{0}-x_{\ell_{2}})+y_{\ell_{2}}-\frac{2\Delta}{\cos(\beta)},\tan(\beta)(x_{0}-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\beta)}].\end{cases}

(47)

Regarding Remark 47, to obtain an inclusive set of feasible coverage for any line $y=\tan(\beta)x+q(\beta)$ where $\beta\in[\theta-\alpha,\theta+\alpha]$ , define

I(x)=\bigcup_{\theta\leq\beta\leq\theta+\alpha}[\tan(\beta)(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\beta)},\tan(\beta)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\beta)}],

and

J(x)=\bigcup_{\theta-\alpha\leq\beta\leq\theta}[\tan(\beta)(x-x_{\ell_{2}})+y_{\ell_{2}}-\frac{2\Delta}{\cos(\beta)},\tan(\beta)(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\beta)}].

A point $(x,y)$ is covered within distance $\Delta$ of some feasible line if and only if $y\in I(x)\cup J(x)$ .

In the following, we find the upper and lower bounds of the set $I(x)\cup J(x)$ to create the set of infeasible $\Gamma_{\ell_{1},\ell_{2}}$ for some location $b_{3}\in L$ which $y_{b_{3}}\notin I(x_{b_{3}})\cup J(x_{b_{3}}),\forall\beta\in[\theta-\alpha,\theta+\alpha]$ . First, the upper bound of set $I(x)$ is discussed as follows.

Define

x^{I}_{u}=x_{\ell_{2}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta+\alpha)}-\sqrt{1+\tan^{2}(\theta)}}{\tan(\theta+\alpha)-\tan(\theta)}\right),

which is the intersection of two extreme lines $\lambda_{1}$ and $\lambda_{3}$ (see Lemma 1 and note that $y=\tan\theta(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos\theta}$ and $y=\tan(\theta+\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta+\alpha)}$ are other representations of $\lambda_{1}$ and $\lambda_{3}$ ). We claim that line $\lambda_{1}$ for $x\leq x^{I}_{u}$ and line $\lambda_{3}$ for $x\geq x^{I}_{u}$ are the upper bounds of $I(x)$ .

Lemma 2.

Function $f(z)=\frac{\sqrt{1+\tan^{2}z}-\sqrt{1+\tan^{2}a}}{\tan z-\tan a}$ is strictly increasing.

Proof: We show the derivative of $f(z)$ is strictly positive:

$\displaystyle f^{\prime}(z)=$	$\displaystyle\frac{\frac{\tan z\sec^{2}z}{\sqrt{1+\tan^{2}z}}(\tan z-\tan a)-\sec^{2}z(\sqrt{1+\tan^{2}z}-\sqrt{1+\tan^{2}a})}{(\tan z-\tan a)^{2}}$	(48)
$\displaystyle=$	$\displaystyle\frac{\sec^{2}z}{(\sqrt{1+\tan^{2}z})(\tan z-\tan a)^{2}}$
	$\displaystyle\ \times(\tan^{2}z-\tan z\tan a-1-\tan^{2}z+\sqrt{1+\tan^{2}z+\tan^{2}a+\tan^{2}z\tan^{2}a})$
$\displaystyle=$	$\displaystyle\frac{\sec^{2}z}{(\sqrt{1+\tan^{2}z})(\tan z-\tan a)^{2}}$
	$\displaystyle\ \times(\sqrt{(1+\tan z\tan a)^{2}+\tan^{2}z+\tan^{2}a-2\tan z\tan a}-\tan z\tan a-1)$
$\displaystyle=$	$\displaystyle\frac{\sec^{2}z}{(\sqrt{1+\tan^{2}z})(\tan z-\tan a)^{2}}$
	$\displaystyle\ \times(\sqrt{(1+\tan z\tan a)^{2}+(\tan z-\tan a)^{2}}-(\tan z\tan a+1))>0.$

Using Lemma 2, we conclude:

$x\leq x^{I}_{u}$ :

$\displaystyle x$	$\displaystyle\leq x_{\ell_{2}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta+\alpha)}-\sqrt{1+\tan^{2}(\theta)}}{\tan(\theta+\alpha)-\tan(\theta)}\right)$		(49)
	$\displaystyle\leq x_{\ell_{2}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\beta)}-\sqrt{1+\tan^{2}(\theta)}}{\tan(\beta)-\tan(\theta)}\right)$	$\displaystyle\forall\beta\in[\theta,\theta+\alpha]$
$\displaystyle\implies$	$\displaystyle\tan(\theta)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta)}\geq\tan(\beta)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\beta)}\quad$	$\displaystyle\forall\beta\in[\theta,\theta+\alpha]$

$x\geq x^{I}_{u}$ :

$\displaystyle x$	$\displaystyle\geq x_{\ell_{2}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta+\alpha)}-\sqrt{1+\tan^{2}(\theta)}}{\tan(\theta+\alpha)-\tan(\theta)}\right)$		(50)
	$\displaystyle\geq x_{\ell_{2}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta+\alpha)}-\sqrt{1+\tan^{2}(\beta)}}{\tan(\theta+\alpha)-\tan(\beta)}\right)$	$\displaystyle\forall\beta\in[\theta,\theta+\alpha]$
$\displaystyle\implies$	$\displaystyle\tan(\theta+\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta+\alpha)}\geq\tan(\beta)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\beta)}\quad$	$\displaystyle\forall\beta\in[\theta,\theta+\alpha]$

Hence, the upper bound of $I(x)$ will be (see Figure 15)

\begin{cases}\tan(\theta)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta)}\text{ (line }\lambda_{1})&x\leq x^{I}_{u}\\ \tan(\theta+\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta+\alpha)}\text{ (line }\lambda_{3})&x\geq x^{I}_{u}\end{cases}

(51)

Likewise, upper bound of $J(x)$ can be determined by defining:

x^{J}_{u}=x_{\ell_{1}}-2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta)}-\sqrt{1+\tan^{2}(\theta-\alpha)}}{\tan(\theta)-\tan(\theta-\alpha)}\right),

as the intersection of lines $\lambda_{1}$ and $\lambda_{5}$ . Employing Lemma 2, similarly it is concluded that the upper bound of $J(x)$ will be

\begin{cases}\tan(\theta)(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\theta)}\text{ (line }\lambda_{1})&x\geq x^{J}_{u}\\ \tan(\theta-\alpha)(x-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\theta-\alpha)}\text{ (line }\lambda_{5})&x\leq x^{J}_{u}\end{cases}

(52)

The lower bounds of $I(x)$ and $J(x)$ can be determined with a similar approach. Define

x^{I}_{l}=x_{\ell_{1}}+2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta+\alpha)}-\sqrt{1+\tan^{2}(\theta)}}{\tan(\theta+\alpha)-\tan(\theta)}\right),

and

x^{J}_{l}=x_{\ell_{2}}+2\Delta\left(\frac{\sqrt{1+\tan^{2}(\theta)}-\sqrt{1+\tan^{2}(\theta-\alpha)}}{\tan(\theta)-\tan(\theta-\alpha)}\right).

Then, the lower bound of $I(x)$ will be

\begin{cases}\tan(\theta)(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\theta)}\text{ (line }\lambda_{2})&x\geq x^{I}_{l}\\ \tan(\theta+\alpha)(x-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\theta+\alpha)}\text{ (line }\lambda_{4})&x\leq x^{I}_{l},\end{cases}

(53)

and for the lower bound of $J(x)$ , we have

\begin{cases}\tan(\theta)(x-x_{\ell_{2}})+y_{\ell_{2}}-\frac{2\Delta}{\cos(\theta)}\text{ (line }\lambda_{2})&x\leq x^{J}_{l}\\ \tan(\theta-\alpha)(x-x_{\ell_{2}})+y_{\ell_{2}}-\frac{2\Delta}{\cos(\theta-\alpha)}\text{ (line }\lambda_{6})&x\geq x^{J}_{l}.\end{cases}

(54)

In conclusion, for given location $(x_{b_{3}},y_{b_{3}})$ if it is located above the upper bounds in (51) and (52) which means

\begin{cases}y_{b_{3}}>\tan(\theta)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\theta)}(=\tan(\theta)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta)}),\\ y_{b_{3}}>\tan(\theta+\alpha)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta+\alpha)}(=\tan(\theta+\alpha)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}),\\ y_{b_{3}}>\tan(\theta-\alpha)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}+\frac{2\Delta}{\cos(\theta-\alpha)}(=\tan(\theta-\alpha)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}),\end{cases}

(55)

OR below the lower bounds in (53) and (54) which implies

\begin{cases}y_{b_{3}}<\tan(\theta)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\theta)}(=\tan(\theta)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}-\frac{2\Delta}{\cos(\theta)}),\\ y_{b_{3}}<\tan(\theta+\alpha)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}-\frac{2\Delta}{\cos(\theta+\alpha)}(=\tan(\theta+\alpha)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}),\\ y_{b_{3}}<\tan(\theta-\alpha)(x_{b_{3}}-x_{\ell_{2}})+y_{\ell_{2}}+\frac{2\Delta}{\cos(\theta-\alpha)}(=\tan(\theta-\alpha)(x_{b_{3}}-x_{\ell_{1}})+y_{\ell_{1}}),\end{cases}

(56)

then $y_{b_{3}}\notin I(x)\cup J(x)$ . The set of infeasible locations for $\ell_{1}$ and $\ell_{2}$ , $\Gamma_{\ell_{1},\ell_{2}}$ , includes locations in $L$ that meet (55) or (56), i.e.,

	$\displaystyle\Gamma_{\ell_{1},\ell_{2}}$	$\displaystyle=\{b_{i}\in L\setminus\{\ell_{1},\ell_{2}\}:y_{b_{i}}\notin I(x_{b_{i}})\cup J(x_{b_{j}})\}$		(57)
		$\displaystyle=\{b_{i}\in L\setminus\{\ell_{1},\ell_{2}\}:(x_{b_{i}},y_{b_{i}})\in R^{\prime\prime}_{\ell_{1},\ell_{2}}\}.$		(57)

Therefore, locations in $\Gamma_{\ell_{1},\ell_{2}}$ cannot be covered by a line that has a distance at most $\Delta$ from two locations $\ell_{1}$ and $\ell_{2}$ , and so $z_{\ell_{1}}+z_{\ell_{2}}+z_{b_{3}}\leq 1$ for $b_{3}\in\Gamma_{\ell_{1},\ell_{2}}$ is a valid cut.

Remark 5.

The boundaries of the cases $\theta+\alpha\geq\pi/2$ or $\theta-\alpha\leq-\pi/2$ can be easily obtained by rotating coordinates to get $\theta^{\prime}$ such that $\theta^{\prime}+\alpha<\pi/2$ and $\theta^{\prime}-\alpha>-\pi/2$ the same as the conditions in the proof.

4 Details on the Feasibility Check when $E<\infty$

Next we detail the feasibility check of the DBC for a given integer solution $z^{h}$ for the case where $E<\infty$ and that the solution of the SLA generates a line of length greater than $E$ . Let $C\subseteq L$ be a subset of locations in $L$ and let $\lambda$ be a line in $\mathbb{R}^{2}$ . We say that $\lambda$ is a covering line for $C$ if the distance from $\lambda$ to any $\ell\in C$ is at most $\Delta$ . Similarly, we will refer to a finite segment $\xi$ that covers all locations in $C$ a covering segment for $C$ ; the length of segment $\xi$ will be denoted by $u_{\xi}$ . For any $a,b\in L$ , let $\lambda_{a}$ and $\lambda_{b}$ the projections of $a$ and $b$ in $\lambda$ , respectively.

Lemma 3.

Let $C\subseteq L$ be given and let $\lambda$ be a covering line for $C$ . If $D(\lambda_{a},\lambda_{b})\leq E$ for all $a,b\in C$ , then there exists a covering segment $\xi$ of $C$ such that $u_{\xi}\leq E$ .

Proof: Let $w,w^{\prime}\in\operatorname{arg\,max}\{D(v,v^{\prime})\colon v=\lambda_{a},v^{\prime}=\lambda_{b},a,b\in C\}$ . That is, $w$ and $w^{\prime}$ are the projections of the points in $C$ on $\lambda$ that attain the maximum distance. If $\xi=[w,w^{\prime}]$ , then clearly $\xi$ covers $C$ and $u_{\xi}\leq E$ .

Let $C\subseteq L$ be given and let $\lambda$ be a covering line of $C$ . We denote by $C_{\lambda,E}\subseteq C$ the set of building pairs in $C$ whose projections in $\lambda$ are at a distance of more than $E$ , i.e.,

C_{\lambda,E}=\bigl{\{}(a,b)\in C\times C\colon D(\lambda_{a},\lambda_{b})>E\bigr{\}}.

(58)

For any $\ell\in L$ let $\delta_{\ell}=D(\ell,\lambda_{\ell})$ and let $\xi_{\ell}=\sqrt{\Delta^{2}-\delta_{\ell}}$ . For any $(a,b)\in C_{\lambda,E}$ let $\lambda_{a}^{\prime}=\lambda_{a}+\xi_{a}(\lambda_{b}-\lambda_{a})$ and $\lambda_{b}^{\prime}=\lambda_{b}+\xi_{b}(\lambda_{a}-\lambda_{b})$ . Note that both $\lambda_{a}^{\prime}$ and $\lambda_{b}^{\prime}$ are points in $\lambda$ .

Lemma 4.

Let $C\subseteq B$ be given and let $\lambda$ be a covering line for $C$ . If $D(\lambda_{a}^{\prime},\lambda_{b}^{\prime})\leq E$ for all $a,b\in C_{\lambda,E}$ then there exist a covering segment $\xi$ of $C$ such that $u_{\xi}\leq E$ .

Proof: By construction, for any $(a,b)\in C_{\lambda,E}$ , both $\lambda_{a}^{\prime}$ and $\lambda_{b}^{\prime}$ are elements of $\lambda$ and $D(\lambda_{a}^{\prime},a)=D(\lambda_{b}^{\prime},b)=\Delta$ , thus the segment $[\lambda_{a}^{\prime},\lambda_{b}^{\prime}]$ covers both $a$ and $b$ . Let $w,w^{\prime}\in\operatorname{arg\,max}\{D(\lambda_{a}^{\prime},\lambda_{b}^{\prime})\colon a,b\in C_{\lambda,E}\}$ , then by definition $D(\lambda_{a}^{\prime},\lambda_{b}^{\prime})\leq D(w,w^{\prime})$ for all $a,b\in C_{\lambda,E}$ and, moreover, $[w,w^{\prime}]$ covers $C$ . Because $D(w,w^{\prime})\leq E$ then $\xi=[w,w^{\prime}]$ is a covering segment of $C$ with $u_{\xi}\leq E$ .

Step 2 in Algorithm 2 which checks the feasibility of solution $z^{h}$ in $\mathcal{U}$ , can be reframed as finding two endpoints $e_{0}$ and $e_{1}$ of some covering segment $\xi$ for $C^{h}$ with $u_{\xi}\leq E$ . To accelerate the search for two endpoints in step 2, we implement SLA for $C^{h}$ along with two lemmas 3 and 4 as shown in Algorithm 3.

Data: Set

\tilde{\mathcal{U}}^{h},C^{h},E

Result:

feasible

1 set

feasible=FALSE

;

2 run SLA for locations in

C^{h}

. Let

\lambda_{MAX}

be the line intersecting maximal number of circles and

\xi_{MAX}

be the shortest covering segment;

3 if maximum coverage with $\lambda_{MAX}=|C^{h}|$ then

4 if $u_{\xi_{MAX}}\leq E$ then

feasible=TRUE

6 end if

7 else

8 if $\tilde{\mathcal{U}}^{h}\neq\emptyset$ then

feasible=TRUE

10 end if

12 end if

14 end if

Algorithm 3 Feasibility check for solution

z^{h}

E<\infty

Algorithm 3 executes SLA for the set of circles with centers in $C^{h}$ and radius $\Delta$ (step 3). If the maximal stabbing line $\lambda_{MAX}$ does not cover locations in $C^{h}$ , then no segment produces the exiting solution which means $z^{h}$ is infeasible. Else, using the discussion in Lemmas 3 and 4, we check whether the shortest covering segment of $\lambda_{MAX}$ , denoted as $\xi_{MAX}$ , has a length at most $E$ (step 3). If so, we recognize $z^{h}$ feasible by $\xi_{MAX}$ . If the length of the segment is greater than $E$ , we find a feasible solution in $\tilde{\mathcal{U}}^{h}$ . If set $\tilde{\mathcal{U}}^{h}$ is not empty, then there are feasible endpoints for a segment covering $C^{h}$ (step 3).

In step 3 we further enhance the tractability of $\tilde{\mathcal{U}}^{h}$ by restricting the search in ranges for feasible endpoints as follows. Let $\ell_{i}$ and $\ell_{j}$ be two extreme locations that have the maximum distance among all pairs of locations in $C^{h}$ . Then, the endpoints of the line segment $e_{0}$ and $e_{1}$ must be within distance $\Delta$ from each $\ell_{i}$ and $\ell_{j}$ . Particularly, points $e_{0}$ and $e_{1}$ should satisfy that (see Figure 16):


	$\displaystyle x_{\ell_{i}}-\Delta\leq e^{x}_{0}\leq x_{\ell_{i}}+\Delta$		(59a)
	$\displaystyle y_{\ell_{i}}-\Delta\leq e^{y}_{0}\leq y_{\ell_{i}}+\Delta$		(59b)
	$\displaystyle x_{\ell_{j}}-\Delta\leq e^{x}_{1}\leq x_{\ell_{j}}+\Delta$		(59c)
	$\displaystyle y_{\ell_{j}}-\Delta\leq e^{y}_{1}\leq y_{\ell_{j}}+\Delta.$		(59d)

The enhancement, therefore, adds the inequalities (59) in the definition of $\tilde{\mathcal{U}}^{h}$ . The enhanced definition is solved significantly faster.

5 Performance of Subproblem Methods

We compare three different methods to address the one-level subproblem $\Phi(f)$ in Algorithm 1:

1.

ORG: solves the subproblem with original set $\mathcal{U}$ (instead of using (6c)) by DBC and only updates $\mathcal{R}(f)$ on the fly,
2.

AVC: adds valid cuts for infeasible combinations in $\mathcal{C}^{0}$ to $\mathcal{U}$ to enhance ORG performance while updates $\mathcal{R}(f)$ on the fly,
3.

DEC: replace $\mathcal{U}$ by conflict constraints (6c) and update $\mathcal{C}$ and $\mathcal{R}(f)$ on the fly as described in Section 5.1.

Table 6 compares the three approaches in terms of optimality gap, solution time, number of iterations, subproblem run time, and its callback function run time. We generate a testbed of 5 sample problems in which 10 blocks in Joplin are randomly selected as locations. Column 2 shows the method that is implemented to solve the subproblem. Columns 3, 4, and 5 present the best objective value (expected population dislocation), the best lower bound, and the optimality gap found within the time limit of one hour, respectively. Columns 6 and 7 show the CPU run time of C&CG Algorithm 1 in seconds and its number of iterations, respectively. Column 8 compares the CPU run time of the DBC method in seconds for each selection of subproblems. Column 9 also presents the total amount of time that the DBC spends in the callback function to generate cuts on the fly with Algorithm 2; note that the callback run time for ORG and AVC methods does not include the time of checking feasibility in $\mathcal{U}$ as we directly use the original feasible solution.

Sample # Subproblem method Best objective Best bound Gap CCG run time (sec.) CCG iteration DBC run time (sec.) Callbacks run time (sec.) 1 ORG 365 - - 3600 1 3600 42 AVC 185 185 0% 6 3 6 0 DEC 185 185 0% 1 3 1 1 2 ORG 161 153 5% 4012 2 4012 1 AVC 153 153 0% 26 3 26 1 DEC 153 153 0% 1 3 1 0 3 ORG 178 178 0% 4484 3 4484 1 AVC 178 178 0% 9 3 9 0 DEC 178 178 0% 2 3 1 1 4 ORG 56 - - 3600 1 3600 0 AVC 75 75 0% 3549 3 3549 1 DEC 75 75 0% 1 3 1 1 5 ORG 103 - - 3600 1 3600 0 AVC 65 65 0% 36 3 36 1 DEC 65 65 0% 1 3 1 0

Table 6: Comparing C&CG performance for 3 different methods to solve the subproblem

The results in Table 6 show that the ORG method is not able to solve any sample with only 10 locations within the one-hour time limit. Also, based on the results, solving the master problem and the second stage recovery problem in the callback function is very fast and almost all the run time in ORG is spent on solving the MINLP subproblem (5). For samples 1, 4, and 5 the subproblem was not solved to optimality in one hour. For some samples, the ORG method exceeds the time limit while solving the MINLP subproblem. In these cases, we wait until the subproblem finishes its execution before terminating the algorithm, which results in some computational times reported being larger than the one-hour time limit (italic values in the table).

Adding infeasible pairs and triples cuts in AVC significantly reduces the computation times by as much as two orders of magnitude, resulting in all the samples being solved in less than one hour. However, AVC can be potentially time-consuming for problems such as sample 4 in the experiment. The DEC method greatly outperforms the other two methods as it solved all samples within a few seconds. The subproblem runtime for DEC shows how fast the subproblem can lead to finding optimal tornado paths. The callbacks’ runtime also emphasizes that the feasibility check with Algorithm 3 in Appendix 4 which is equipped by the polynomial stabbing line algorithm is a quick task.

6 Additional Experiments for 30 Days of Recovery

We conduct more experiments for locations in Joplin to study the population dislocation after 30 days of recovery. The results are shown in Table 7.

Table 7: Solving the two-stage robust optimization model with different parameters of tornado length and available budget, after 30 days of recovery for 100 locations group in Joplin, MO.

Length Budget Population CCG CCG DBC Callback $\tilde{\mathcal{U}}^{h}$ feasibility Retrofitting cost Recovery cost ($) dislocation run time (sec.) iteration run time (sec.) run time (sec.) run time (sec.) /budget /budget 5 0 M 16465 153 2 153 5 2 - - 15 M 14847 194 2 193 5 2 0.6 0.4 30 M 14688 162 2 161 5 1 0.3 0.7 $\infty$ 0 M 17449 361 2 360 3 0 - - 15 M 15742 527 2 527 5 0 0.6 0.4 30 M 15582 500 2 500 3 0 0.3 0.7

Compared with the results in Table 2, we observe more population dislocation for each set of parameters because of the shorter period to recover the damaged locations. Despite this difference in the magnitude of population dislocation, it is important to note that the information, clarifications, and examinations between both Table 2 and 7 remain consistent. Therefore, Table 7 provides numerical experiments with a focus on the impact of shorter recovery times, complementing the information presented in Table 2 in the main body of the paper.