1 Introduction

Communities in coastal areas are at an ever-increasing risk to prolong power outages caused by hurricanes. Each year during hurricane season, hurricanes cause catastrophic failures in the distribution network and threaten U.S. coastal communities at exponentially higher rates with the increase of populations in these areas. Therefore, communities become more dependent on electric power, and the determinant characteristics of hurricanes (i.e., strength and frequency) change with a changing global climate. As costs resulting from power outages increase, assessing the resilience of the distribution network against hurricanes becomes an emerging task of distribution system operators (DSOs) [1]. In the hours and days leading up to hurricane landfall, the DSO must make logistical and operational decisions crucial to the resilience of the power grid and the safety of its customers. The DSO must quickly anticipate the impacts of hurricane on the grid components and accurately assess the grid resilience with its existing resources given the immediate threat of hurricane-induced outages. These functionalities enable the DSO to reschedule grid operation, and therefore minimize the negative impacts of hurricane on the customers. This analysis additionally provides useful foresight on the effectiveness of its existing resources and hardening measures in preparation for future catastrophic events. Within this context, the aim of this paper is to develop a comprehensive framework to assess the resilience of coastal distribution networks against the threat of outage following hurricane landfall. Specifically, this paper focuses on the benefit of dispatchable energy storage units to store grid-powered energy for post-impact use.

Multiple literature has studied the service restoration of distribution systems in the hours following a catastrophic event. The service restoration methods proposed in the literature inlcude optimal crew assignment problem [2, 3] and deployments of mobile generators [4]. In [5], the service restoration of distribution network with transportable energy storage is also investigated wherein transportable energy storage is proposed to be transported over railway networks to support different islanded areas of a degraded distribution network. The authors assume pre-defined outage scenarios, which neglect the uncertainty and wide variety of outage scenarios that can occur due to hurricane events. While simulation results show that the grid resilience can be improved significantly via rail-transported energy storage, the use of railway network is not always applicable for generic urban distribution networks, e.g., our case study in Houston. In fact, energy storage systems play an important role in enhancing grid resilience, especially for hurricane outage mitigation, and will do so at an increasing value to customers as a result of improving charging efficiency, capacity, and increasing penetration of energy storage in distribution networks [6,7,8,9]. Additionally, as a result of previous hurricane events, many urban areas frequently affected by hurricanes have outlined methods for proper protection, placement, and operation of energy storage units for resilience enhancement. Other researches suggest to take sectionalizing and networked microgrids as a method for enhancing the resilience of power grid following hurricane impact. References [10, 11] present studies in which the distribution network is sectionalized into multiple microgrids for load restoration. These studies consider a known outage event without specifying what causes the outages and how the outages occur. References [12, 13] study the coordination of networked micorgrid as well as power grids and district heating systems for enhancing power grid resilience.

Reference [14] proposes an adaptive robust energy management framework for distribution network against the intermittent interruption from the upstream transmission network while the remaining distribution network topology is not affected. However, the outage scenarios in this work may be too simplified to accurately reflect hurricane impacts. In addition, most of these works employ the linearized DistFlow model [15], which has some inaccuracies in modeling the operation of distribution network in pre-outage and post-outage events.

Preparedness ahead of hurricane-induced outages has also been studied in the literature, which requires modeling the impacts of hurricanes on the distribution network [16, 17]. For example, a negative binomial regression model based on outage data caused by three specific hurricanes occurred in U.S. from 1996 to 1999 is proposed in [16]. The efforts to improve the hurricane model [16] to be applicable for the whole U.S. coastal area are presented in [17]. However, the models proposed in [16, 17] focus on modeling hurricane impacts in long-term planning problems and are not flexible in capturing the spatio-temporal changes of uncertain parameters of hurricanes in short-term notification. In particular, due to strong wind gusts and the possibility of flood, multiple distribution lines may also be damaged and a quick assessment of the grid resilience against the threat of spatio-temporal hurricane-induced outage becomes an emergent task of the DSO. Using long-term hurricane impact models, different resilience enhancement based distribution network expansion strategies such as hardening distribution lines [18, 19] or placements of distributed energy resources (DERs) have been proposed. Moreover, as mentioned above, the linearized DistFlow model employed in [18, 19] might be not accurate to efficiently assess the impacts of hurricane-induced outages on the operation of distribution networks. The resource preparation planning for grid resilience enhancement, such as preparing energy storage units, distributed generators, and crew prepositions, is studied in [20]. However, the proposed assessment methods using DC power flow models are only applicable for modeling transmission networks and are not accurate for modeling distribution networks.

Although there has been multiple recent researches in the literature, the assessment of distribution network resilience against the threat of a hurricane is an on-going and challenging research subject. In particular, the spatio-temporal relationship of hurricanes over the course of their progression are vital for accurate scenario modeling and subsequent resilience studies [1]. In addition, it is important to incorporate the distribution network operation constraints such as the conic network power flows [21, 22] that accurately model the operation of radial distribution networks.

This paper aims to fill the gap in literature by providing an effective resilience assessment framework for distribution networks equipped with energy storage in the impending threat of an approaching hurricane. More specifically, this paper develops a resilience assessment framework that utilizes the high-quality predictions produced by the proposed spatio-temporal hurricane impact analysis (STHIA) model, and a second-order conic program (SOCP) optimization method to assess the impacts of energy storage systems on the resilience of distribution network in coastal areas. The proposed model utilizes rich sources of hurricane data provided by the National Hurricane Center (NHC) of the U.S. National Weather Service. The major contributions of this paper are summarized as follows:

  1. 1)

    This paper develops a comprehensive spatio-temporal model of hurricane simulation and grid outage analysis. The proposed model can be used to anticipate the hurricane-induced outages in distribution networks after a hurricane makes the landfall in a coastal area.

  2. 2)

    This paper presents an optimization based assessment framework to analyze the impact of energy storage systems on the resilience of distribution systems given the spatio-temporal outages induced by a hurricane. The proposed model incorporates the grid operation in both pre-outage and post-outage events, which is casted as a SOCP and can be solved by available commercial solvers.

The rest of the paper is organized as follows. Section 2 presents and explains the spatio-temporal hurricane model and the damage analysis on the grid components after the hurricane landfalls. Section 3 provides the mathematical formulation for the resilience assessment of the distribution network considering the hurricane-induced outages and the use of energy storage units. Numerical results conducted on the modified IEEE 33-bus distribution network uisng real hurricane data analysis in the U.S. coastal area is presented in Section 4 to highlight the effectiveness of the proposed method. Section 5 concludes the paper.

2 Spatio-temporal hurricane-induced outage prediction

Figure 1 illustrates the key aspects of STHIA for power grid resilience studies. The proposed model is composed of two major sequential steps, i.e., spatio-temporal hurricane simulation and grid component impact analysis, the output of which is a simulated grid outage scenario including both the time and locations of the faults.

Fig. 1
figure 1

Stochastic STHIA model

Full size image

2.1 Spatio-temporal hurricane simulation

The spatio-temporal path and strength of the hurricane are the most important factors that affect the distribution network. The pressure difference of the eye of hurricane \(\Delta P\) is widely used by atmospheric scientists to capture the size and strength of the hurricane [23, 24], and can be determined from a number of common observation methods including aircraft fly-over, satellite imagery, etc. The pressure of the eye of hurricane \(P_{eye}\) can then be calculated by \(\Delta P\) and the average sea level pressure \(P_{atm}\) (which is approximately 1013 mbar) as follows:

$$\begin{aligned} P_{eye} = P_{atm} - \Delta P \end{aligned}$$
(1)

Consequently, the corresponding maximum sustained wind speed \(\overline{V}\) is then approximated by interpolating in the Saffir-Simpson table [25] given the hurricane central pressure. The parameters for this computation are summarized in Table 1.

Table 1 Parameters for different hurricane categories
Full size table

The wind speed of the hurricane at a particular location is defined as a function of its distance from the eye and the maximum sustained wind speed \(\overline{V}\) as follows [26, 27]:

$$\begin{aligned}&V(r) = \Bigg (\frac{AB \Delta P{\mathrm {e}}^{-A/r^B}}{\rho r^B }+ \frac{r^2f^2}{4}\Bigg )^{\frac{1}{2}}\! - \!\frac{rf}{2} \end{aligned}$$
(2)
$$\begin{aligned}&{\left\{ \begin{array}{ll} A = R^{B}_{\max }\\ B = \frac{{\overline{V}}^2 \rho f}{\Delta P}\\ f = 2 \varOmega \phi \end{array}\right. } \end{aligned}$$
(3)
$$\begin{aligned}&R_{\max } = \text{exp}(2.556 - m_1 \Delta P^2 + m_2 \phi ) \end{aligned}$$
(4)

where A and B are the shape parameters of the wind speed profile of the hurricane; r is the radius from hurricane center; \(\rho\) is the density of air; \(R_{\max }\) is the radius of the maximum wind; \(\varOmega\) is the angular velocity of the Earth; \(\phi\) is the storm latitude; \(m_1 =0.000050255\); and \(m_2=0.042243032\).

The temporal hurricane path is defined by the following time series based parameters: ① approach angle \(\varPhi _t\), ② landfall location \(\Delta x_t\), and ③ translation velocity \(c_t\) that represents the translation movement speed of the eye of hurricane and determines the temporal dimension of the hurricane path. For simplicity, hurricane paths are assumed to be linear with a constant-speed translation c.

Given the values of \(\varPhi _t\), \(\Delta x_t\) and \(c_t\), the temporal location of the hurricane can be calculated as follows:

$$\begin{aligned}&Y_{t} = ct - \frac{Y_{tot}}{2} \end{aligned}$$
(5)
$$\begin{aligned}&Y_{\mathrm {tot}} = c N_T \end{aligned}$$
(6)

where \(Y_{t}\) is the temporal location of hurricane at time t; \(Y_{tot}\) is the total distance which the hurricane traveled in the scheduling horizon; and \(N_T\) is the total number of time slots in the scheduling horizon. This paper considers the spatio-temporal hurricane model for emergency preparation of the grid. Hence, the initial location of the hurricane is considered to be half of its total distance traveled such that the hurricane makes landfall approximately half way through the horizon.

2.2 Impact analysis of grid component

This step can be summarized as follows. Firstly, the locations and structural properties of the grid components are used for quantifying the impacts of a hurricane on a proximate power distribution network. Then, a vector between the location of the grid component and the changing location of the eye of hurricane is used to determine the component-to-eye distance. We then map the obtained distance to the radial wind speed curve in [26] to calculate the wind speed experienced by each component over the time horizon considered. The details of the impact analysis are presented in the following.

The spatio-temporal distance between the grid component and the hurricane, which is called the grid component-to-eye vector, is used to map the wind speed analysis in Section 2.1 to the location of the grid component. In this paper, we neglect the Earth’s curvature (the earth’s curvature will be considered in future works for capturing large-scale distribution networks spanning in a considerably large area) and consider that all grid components can be mapped on a 2-dimensional plane. We consider the landfall location as the origin of the coordinate system. Figure 2 illustrates a Cartesian plane aligned on an assumed straight coastline that determines the locations of grid components relative to the eyes of hurricanes. For simplicity, we do not reference the cardinal directions (north, south, etc.) in this paper. All component locations are initially specified in terms of their proximity to the expected landfall location as \(x_i\) and \(y_i\), where i is the index of grid component. However, by shifting the coordinate system by the approach angle \(\phi\), the locations of grid components can be expressed in terms of rotated variables \(x'_{i }\) and \(y'_{i }\), which are aligned with the hurricane path. Hence, the magnitude of the component-to-eye variable between the eye of hurricane and the component \(D_{i,t}\) can be calculated by applying Pythagorean theorem on \(x'_{i}\) and \(y'_{i}\) as follows:

$$\begin{aligned}&x'_{i} = (x_i + \Delta x_i) \cos(\phi) + y_i \sin(\phi) \end{aligned}$$
(7)
$$\begin{aligned}&y'_{i} = -(x_i + \Delta x_i)\sin(\phi) + y_i \cos(\phi) \end{aligned}$$
(8)
$$\begin{aligned}&D_{i,t} = \sqrt{(Y_{t} + y'_{i})^2 + (x'_{i})^2} \end{aligned}$$
(9)
Fig. 2
figure 2

Cartesian coordinate system and variables representing locations of grid components relative to eye of hurricane

Full size image

The impact of landfall event on the wind speed of the hurricane and consequently the distribution network outage prediction should be considered carefully. In general, the hurricane is assumed to have a constant central pressure and thus constant radial wind profile. Hence, we can simply generate a time-based wind speed curve by plugging in the hurricane-to-component variable \(D_{i,t}\) into the radial wind speed profile formula in (2) while the hurricane remains over water. However, the strength of the hurricane will decay after the landfall occurs, which is captured by the following equation [28].

$$\begin{aligned}&V^{\mathrm {s}}_{i,t} = V_b +(V (D_{i,t})R - V_b) \text{exp}(-\alpha (t-T_{lf})) \end{aligned}$$
(10)

where \(\alpha\) is the decay constant; \(V_b\) is the background sustained wind speed in the hours following a hurricane event; R is a wind speed reduction factor; and \(T_{lf}\) denotes the time period in which the hurricane makes landfall. The values of parameters \(V_b\), \(\alpha\), R vary by region [28]. This paper considers the Gulf Coastal area where \(\alpha =0.095\) per hour, \(V_b\) = 13.735 m/s, and R = 0.9. Note that neglecting wind speed decay when the hurricane is landing could result in overestimating the wind speed of the hurricane and the prediction results of grid outages. For each time period t, the wind speed value from the wind speed curve of the component is mapped to the fragility model of component. For simplicity, this paper only considers the line outage induced by the hurricane, which can be simplified as follows:

$$\begin{aligned} \alpha ^{ij} = {\left\{ \begin{array}{ll} 1 &{} V^{\mathrm {gust}}_{ij,t} \ge v_{ij}^{\mathrm {limit}}\\ 0 &{} V^{\mathrm {gust}}_{ij,t} < v_{ij}^{\mathrm {limit}} \end{array}\right. } \end{aligned}$$
(11)

where \(v_{ij}^{\mathrm {limit}}\) is the limit value of wind speed that the line can sustain; and \(V^{\mathrm {gust}}_{ij,t} \!=\! 1.287 V^{\mathrm {s}}_{ij,t}\) [27] is three-second gust wind speed at the middle of line ij.

3 Grid resilience assessment

3.1 Mathematical formulation for grid operation under hurricane-induced outage

Based on the information of hurricane-induced outage that can be obtained from the proposed model, the DSO can assess the grid resilience by simulating optimal grid operation in both pre-outage and post-outage events. The optimal operation of distribution network under hurricane-induced outage can be modeled by the following optimization model:

$$\begin{aligned} \min \sum \limits _{t \in \mathcal {T}}\left( \lambda ^{\mathrm {grid}}_t P_t^{\mathrm {grid}} + \sum \limits _{i \in \mathcal {I}} \sum \limits _{c \in \mathcal {C}} \lambda ^{\mathrm {curtail}}_{c,i,t} \left( 1-\gamma _{c,i,t} \right) P_{c,i,t}\right) \tau \end{aligned}$$
(12)

s.t.

$$\begin{aligned}&\underline{E}_{i} \le E_{i,t} \le \overline{E}_{i} \quad \forall i \in \mathcal {I}, \forall t \end{aligned}$$
(13)
$$\begin{aligned}&E_{i,t+1 } = E_{i,t} + \left( \eta ^c_i C_{i,t} - \frac{1}{\eta ^d_i} D_{i,t}\right) \tau \nonumber \\&\quad\quad\quad\quad\quad \forall i \in \mathcal {I}, ~ 1 \le t \le T_{\text{end}} -1, \forall s \end{aligned}$$
(14)
$$\begin{aligned}&{\left\{ \begin{array}{ll} 0 \le C_{i,t}\quad \forall i \in \mathcal {I}, \forall t\\ D_{i,t} \le \overline{B}_{i} \quad \forall i \in \mathcal {I}, \forall t \end{array}\right. } \end{aligned}$$
(15)
$$\begin{aligned}&- \overline{Q}^{\mathrm {e}}_{i} \le Q^{\mathrm {e}}_{i,t} \le \overline{Q}^{\mathrm {e}}_{i} \quad ~\forall i \in \mathcal {I}, \forall t \end{aligned}$$
(16)
$$\begin{aligned}&E_{i,1 } = \overline{E}_{i} \quad \forall i \in \mathcal {I} \end{aligned}$$
(17)
$$\begin{aligned}&P_{t}^{\mathrm {grid}} + D_{1,t} - C_{1,t} - \sum \limits _{c \in \mathcal {C}_1} \gamma _{c,1,t} P_{c,1,t} = \sum \limits _{j \in \mathcal {N}(1)} P^{1j}_{t} \quad \forall t \end{aligned}$$
(18)
$$\begin{aligned}&{Q_{t}^{\text {grid}}} {+} {Q^{\text{e}}_{1,t}} {-} {\sum _{c \in \mathcal {C}_1} \gamma _{c,1,t} Q_{c,1,t}} {=} {\sum \limits _{j \in \mathcal {N}(1)}} {Q^{1j}_{t}} \quad \forall t \end{aligned}$$
(19)
$$\begin{aligned}&D_{i,t} - C_{i,t} - \sum _{c \in \mathcal {C}_i} \gamma _{c,i,t} P_{c,i,t} = \sum \limits _{j \in \mathcal {N}(i)} P^{ij}_{t} \quad i \in \mathcal {I}_{t}, \forall t \end{aligned}$$
(20)
$$\begin{aligned}&Q^{\mathrm {e}}_{1,t}- \sum _{c \in \mathcal {C}_i} \gamma _{c,i,t} Q_{c,i,t} = \sum \limits _{j \in \mathcal {N}(i)} Q^{ij}_{t} \quad i \in \mathcal {I}_{t}, \forall t \end{aligned}$$
(21)
$$\begin{aligned}&{\left\{ \begin{array}{ll} 0\le u^{ij}_{i,t} \le \frac{\overline{V}^2_{i}}{\sqrt{2}} \alpha ^{ij}_{t} \quad \forall t, \forall ij \in \mathcal {L} \\ 0 \le u^{ij}_{j,t} \le \frac{\overline{V}^2_{j}}{\sqrt{2}}\alpha ^{ij}_{t} \quad \forall t, \forall ij \in \mathcal {L} \end{array}\right. } \end{aligned}$$
(22)
$$\begin{aligned}&0\le \nu _{i,t} - u^{ij}_{i,t} \le \frac{\overline{V}^2_{i}}{\sqrt{2}}\left( 1-\alpha ^{ij}_{t}\right) \quad \forall t, \forall ij \in \mathcal {L} \end{aligned}$$
(23)
$$\begin{aligned}&0 \le \nu _{j,t} - u^{ij}_{j,t} \le \frac{\overline{V}^2_{j}}{\sqrt{2}}\left( 1-\alpha ^{ij}_{t}\right) \quad \forall t, \forall ij \in \mathcal {L} \end{aligned}$$
(24)
$$\begin{aligned}&P^{ij}_{t} = \sqrt{2}g^{ij} u^{ij}_{i,t} - g^{ij} R^{ij}_{t} - b^{ij} T^{ij}_{t} \quad \forall t, \forall ij \in \mathcal {L} \end{aligned}$$
(25)
$$\begin{aligned}&Q^{ij}_{t} = -\sqrt{2} \left( b^{ij}+\frac{b_\mathrm {sh}^{ij}}{2}\right) u^{ij}_{i,t} + b^{ij} R^{ij}_{t}- g^{ij} T^{ij}_{t} \quad \forall t, \forall ij \in \mathcal {L} \end{aligned}$$
(26)
$$\begin{aligned}&{\left\{ \begin{array}{ll} 2 u^{ij}_{i,t} u^{ij}_{j,t} \ge \left( R^{ij}_{t}\right) ^2 + \left( T^{ij}_{t}\right) ^2 \quad \forall t, \forall ij \in \mathcal {L} \\ &R^{ij}_{t} \ge 0 \quad \forall t, \forall ij \in \mathcal {L} \end{array}\right. } \end{aligned}$$
(27)
$$\begin{aligned}&\left( P^{ij}_{t}\right) ^2 + \left( Q^{ij}_{t}\right) ^2 \le \left( S^{ij}_{\max } \right) ^2 \quad \forall t, \forall ij \in \mathcal {L} \end{aligned}$$
(28)
$$\begin{aligned}&\frac{\underline{V}^2_{i}}{\sqrt{2}}\le \nu _{i,t} \le \frac{\overline{V}^2_{i}}{\sqrt{2}} \quad \forall t, \forall i \in \mathcal {I} \end{aligned}$$
(29)

where \(\mathcal {T}\) is the set of time slots; \(\mathcal {I}\) is the set of system buses; \(\mathcal {L}\) is the set of distribution lines; \(\mathcal {C}\) is the set of customers; \(\mathcal {N}(i)\) is the set of buses connecting to bus i by a line in the original network; \(\mathcal {C}_{i}\) is the set of customers c at bus i; \(\tau\) is the length of time interval (15 min); \(\alpha ^{ij}_{t}\) is the load connection status, equal 1 if line ij is connected in time t in scenario s, 0 if otherwise; \(\lambda ^{\mathrm {grid}}_{t}\) is the local marginal price forecast at time t; \(\lambda ^{\mathrm {curtail}}_{c,i,t}\) is the value of lost load type c of bus i in time t; \(P_{c,i,t}\) and \(Q_{c,i,t}\) are the active and reactive power demands of customer c at bus i in time t, respectively; \(\underline{E}_{i}\) and \(\overline{E}_{i}\) are the minimum and maximum energy levels of the battery at bus i; \(B_{i}\) is the rating power of the battery at bus i; \(\eta ^c_i\) and \(\eta ^d_i\) are the charging and discharging efficiencies of energy storage at bus i, respectively; \(b_{\mathrm {sh}}^{ij}\) is the shunt susceptance of line ij; \(g^{ij}\) and \(b^{ij}\) are the resistance and susceptance of line ij, respectively; \(\underline{V}_{i}\) and \(\overline{V}_{i}\) are the minimum and maximum voltages of bus i, respectively; \(C_{i,t }\) and \(D_{i,t }\) are the charging and discharging power of energy storages at bus i in time t, respectively; \(Q^{\mathrm {e}}_{i,t }\) is the reactive power dispatched from energy storage at bus i at time t; \(Q^{\mathrm {e}}_i\) is the maximum reactive power dispatched from energy storage at bus i; \(E_{i,t }\) is the energy level of battery at bus i in time t in scenario s; \(P_{t}^{\mathrm {grid}}\) and \(Q_{t}^{\mathrm {grid}}\) are the active and reactive power supplied from upstream network at time t, respectively; \(P^{ij}_{t}\) and \(Q^{ij}_{t}\) are the active and reactive power flows from bus i to bus j over line ij at time t, respectively; \(R^{ij}_{t}\) and \(T^{ij}_{t }\) are the auxiliary variables used for calculating power flows in SOCP formulation; \(S^{ij}_{\max }\) is the maximum line capacity; \(\nu _{i,t}\) and \(u^{ij}_{i,t}\) are the auxiliary variables used for calculating bus voltages in SOCP formulation; and \(\gamma _{c,i,t}\) is the percentage of the load served of customer c at bus i in time t.

The objective function (12) represents the cost of energy imported from the upstream grid and the cost of energy not served due to outages.

Constraints (13)-(17) express the operation constraints of the battery including battery minimum and maximum energy level constraints (13), temporal dynamics of energy level (14), charging and discharging power limit constraints (15), reactive power dispatch constraint (16), and initial battery level condition (17). In this paper, we assume that the battery has been fully charged before the hurricane makes the landfall and affects the grid area. The power injection constraints are represented in (18)-(21) where active and reactive power injections at the substation bus (bus 1 in this paper) are captured in (18) and (19), and active and reactive power injections at other buses are captured in (20) and (21).

The set of constraints (22)-(28) represent the distribution network power flows in conic form considering the connection status of the line \(\alpha ^{ij}_t\) [21] using auxiliary variables \(R^{ij}_{t}, T^{ij}_{t}, \nu _{i,t}, u^{ij}_{i,t}, u^{ij}_{j,t}\). In particular, \(R^{ij}_{t},T^{ij}_{t}\) and \(\nu _{i,t}\) are defined as follows:

$$\begin{aligned}&R^{ij}_{t} = V_{i,t} V_{j,t} \cos {(\theta _{ij})} \end{aligned}$$
(30)
$$\begin{aligned}&T^{ij}_{t} = V_{i,t} V_{j,t} \sin {(\theta _{ij})} \end{aligned}$$
(31)
$$\begin{aligned}&\nu _{i,t} = \frac{V_{i,t}^2}{\sqrt{2}} \end{aligned}$$
(32)

where \(V_{i,t}\) is the bus voltage magnitude; and \(\theta _{ij}\) is the phase angle. In addition, two auxiliary variables \(u^{ij}_{i,t}, u^{ij}_{j,t}\) are defined for two buses i, j respectively to capture the spatio-temporal hurricane-induced outage in line ij as follows:

$$u_{{i,t}}^{{ij}} = \left\{ {\begin{array}{*{20}c} {\frac{{V_{i}^{2} }}{{\sqrt 2 }}} & {\alpha _{t}^{{ij}} = 1} \\ 0 & {\alpha _{t}^{{ij}} = 0} \\ \end{array} } \right.$$
(33)
$$u_{{j,t}}^{{ij}} = \left\{ {\begin{array}{*{20}l} {\frac{{V_{j}^{2} }}{{\sqrt 2 }}} \hfill & {\alpha _{t}^{{ij}} = 1} \hfill \\ 0 \hfill & {\alpha _{t}^{{ij}} = 0} \hfill \\ \end{array} } \right.$$
(34)

It is easy to observe that, if the line ij is connected, the following equality holds:

$$\begin{aligned} 2 u^{ij}_{i,t} u^{ij}_{j,t} = \left( R^{ij}_{t}\right) ^2 + \left( T^{ij}_{t}\right) ^2 \end{aligned}$$
(35)

Remark: the values of phase angles and bus voltage magnitudes can be easily recovered from auxiliary variables \(R^{ij}_{t},T^{ij}_{t},\nu _{i,t}, u^{ij}_{i,t},u^{ij}_{j,t}\) based on their definitions (30)-(34).

The distribution network power flow constraints (22)-(28) can be explained as follows. Constraints (22)-(24) represent the definition (33) and (34) of \(u^{ij}_{i,t}, u^{ij}_{j,t}\) based on the line connection status \(\alpha ^{ij}_{t}\). The active and reactive power flows \(P^{ij}_{t}, Q^{ij}_{t}\) are calculated from auxiliary variables \(R^{ij}_{t},T^{ij}_{t},u^{ij}_{i,t}, u^{ij}_{j,t}\) as in (25) and (26). Constraint (27) which defines a rotated quadratic cone is the relaxation of equality condition (35). Constraint (28) represents the line thermal constraints. The bus voltage limit requirement \(\underline{V}_{i} \le V_{i,t} \le \overline{V}_{i}\) is rewritten in (29) using the new auxiliary variable \(\nu _{i,t}=V_{i,t}^2/\sqrt{2}\).

The impact of line connection status on the network power flow constraint can be explained as follows. When \(\alpha ^{ij}_{t}=1\) (i.e. line ij is connected), we can obtain \({u^{ij}_{i,t}}\,=\,{ \nu _{i,t}}{=}{V^2_{i,t}/\sqrt{2}}\) and \({u^{ij}_{j,t}}\,=\,{\nu _{i,t}}{=}{V^2_{i,t}/\sqrt{2}}\) based on constraints (22)-(24). Hence, constraints (25)-(28) express the regular branch flow equations for the connected line:

$$\begin{aligned} P^{ij}_{t }\,=\, & {} g^{ij} {\left( V_{i,t }\right) ^2} - g^{ij} V_{i,t } V_{j,t }\cos {(\theta _{ij})} - b^{ij} V_{i,t } V_{j,t } \sin {(\theta _{ij})} \end{aligned}$$
(36)
$$\begin{aligned} Q^{ij}_{t}= & {} - \left( b^{ij}+\frac{b_{\mathrm {sh}}^{ij}}{2}\right) {V_{i,t}^2} + b^{ij} V_{i,t } V_{j,t } \cos {(\theta _{ij})} \nonumber \\&- g^{ij} V_{i,t } V_{j,t } \sin {(\theta _{ij})} \end{aligned}$$
(37)

The conic power flow model (22)-(24) is exact for radial distribution network [21]. Hence, the set of constraints (22)-(24) accurately model the operations of radial distribution networks in pre-outage events. When \(\alpha ^{ij}_{t}=0\) (i.e. line ij is disconnected), \({u^{ij}_{i,t}}\,=\,{u^{ij}_{j,t}}\,{=}\,{0}\), constraints (22)-(24) enforce \({u^{ij}_{i,t}}\,=\,{u^{ij}_{j,t}}\) to be zero. Hence, constraint (27) results in \({R^{ij}_{t}}\,{=}\,{0}\) and \({T^{ij}_{t}}\,{=}\,{0}\), which means active and reactive power flows \(P^{ij}_{t},Q^{ij}_{t}\) calculated in (25) and (26) are zero. In other words, there are no active and reactive power flows in the disconnected line ij. That means the constraints (22)-(24) accurately model the operation of radial distribution networks in post-outage events.

The proposed grid assessment problem (12)-(29) is a SOCP, which can be solved efficiently by available convex optimization solvers such as CONOPT with reduced gradient methods, by MOSEK with interior-point methods, or by CPLEX with barrier methods [21, 22]. Specifically, constraint (27) is a rotated cone constraint [21], constraint (28) is a quadratic cone constraint, the objective (12) and remaining constraints are all linear. Hence, the problem (12)-(29) can be solved directly by CPLEX with quadratic constraint programming mode (QCP) [29].

4 Case study

In this paper, the modified IEEE 33-bus distribution system is used to represent the distribution network as shown in Fig. 3. Network data is given in [15]. We map the locations of 33 buses of the network to the Houston Metro system to analyze the impact of the approaching hurricane on the distribution network. Specifically, we consider a typical hurricane scenario in the Gulf-Coast region. The longitude and latitude of (29.7604 N, 95.3698 W) in Texas were initialized as the expected landfall location. We consider a 24-hour time horizon starting from the time when the hurricane makes landfall. An interval of 15 min is used in this paper to ensure the gust wind speed can be captured in the hurricane simulation and outage impact analysis. Table 2 summarizes key inputs used in the spatio-temporal hurricane model which is based on the study in [27]. Detailed explanation of these parameters can be found in [27]. The impact analysis of hurricanes on distribution system components is based on the STHIA model proposed in [30]. Spatio-temporal visualization of few major recent hurricanes and their impacts on power grids are presented in [31]. The impacts of the hurricane Harvey on power grid in August 26, 2017 18:00 Coordinated Universal Time (UTC) is illustrated in Fig. 4 [31]. The distribution network is assumed to connect to the upstream grid at bus 1. The price of energy imported from the upstream grid \(\lambda ^{\mathrm {grid}}\) is based on the CAISO LMP on February 2, 2016, which is shown in Figure 5. The price exhibits a typical “duck” curve with two peaks.

Table 2 Hurricane characteristic inputs
Full size table

In this paper, we investigate two case studies:

  1. 1)

    Case 1: medium hurricane scenario with moderate wind speed decay factor after making landfall, \(\alpha =0.095\) per hour.

  2. 2)

    Case 2: strong hurricane scenario with small or ignored decay factor after making landfall, \(\alpha =0\).

The analysis of hurricane-induced outages is presented in Table 3. Since Houston is a coastal city, Case 1 can be considered as the case that the hurricane does not pass through the city directly. Although the grid is still affected, the strength of wind speed decays since the hurricane does not make landfall in the city. Case 2 can be considered as the extreme scenario, i.e., the eye of the hurricane directly passes through the center of the city and the grid components are subjected to the full strength of the hurricane. It can be observed that more lines are damaged by the hurricane and the outages occur sooner in Case 2 than in Case 1. Clearly, distribution network components laying directly in, or near to, the path of the hurricane will be more severely affected as they are subjected to higher wind speeds, and thus a higher probability of failure. This result highlights the correctness of the STHIA model developed in this paper.

Fig. 3
figure 3

Distribution network diagram

Full size image
Fig. 4
figure 4

Spatio-temporal visualization of Hurricane Harvey

Full size image
Fig. 5
figure 5

CAISO LMP on February 2, 2016

Full size image
Table 3 Hurricane-induced outages
Full size table

Figure 6 captures the operation of the distribution network in both pre- and post-outage events where the network connection is affected by the wind speed of the hurricane. In particular, the outage scenario affects the energy resource scheduling of the grid. A new peak is observed before the wind speed of the hurricane starts affecting the grid area, i.e., the DSO needs to import more energy from the upstream grid to ensure all energy storage units are fully charged. Note that the negative power flow in Fig. 7 denotes energy storage discharge whereas positive power flow denotes energy storage charging. Once the outage occurs, the energy storage units in the islanded network begin to discharge energy. Note that in Case 1, energy storage units E1 and E3 remain connected to the upstream grid, which allows them to be charged at the end of the scheduling horizon. It is worth mentioning that the power imported from the upstream grid is smaller than the total demand. This can be explained by the fact that parts of the network are isolated from the main grid and can only be served by energy storage units E2 and E4. In addition, the energy reserved in energy storage units E1 and E3 in the pre-outage event also helps the DSO to flexibly reduce the imported energy. These results distinguish the operation of distribution network in pre-outage and post-outage events.

Fig. 6
figure 6

Distribution grid operation results

Full size image

The impact of energy storage capacity in the grid is illustrated in Fig. 7. In particular, we increase the size \(E_{i,\max }\) by multiplying its value with a scale factor 0, 0.5, 1, 1.5, and 2, respectively. It can be observed that when there is no energy storage (i.e., scale factor is 0) in the system, the total cost, cost of energy not served, and energy not served increase since the network has no backup resources to support the islanded parts of the network. As the scale factor increases, the distribution network has more backup resources, which results in the reduction of the cost of energy not served and the energy not served. Consequently, it leads to the reduction of total cost. However, it is worth mentioning these values become saturated at some points as the scale factor increases, expressing the limits of using energy storage units in supporting the load after outages. This illustrates the need of active network management such as network reconfiguration to support the customers’ demands after outage events and utilize the resources more effectively.

Fig. 7
figure 7

Impact of energy storage capacity on grid resilience

Full size image

5 Conclusion

This paper presents a comprehensive framework to analyze and assess the impacts of hurricane-induced outages and energy resources on distribution system operation during impending hurricane landfall. Specifically, we provide a framework to analyze the impact of energy storage in improving the resiliency of the distribution network with the spatio-temporal outages induced by hurricanes. Our proposed framework combines a spatio-temporal model of hurricane progression based on empirical weather parameters, as well as a second-order conic optimization model for power grid operation during hurricane events. Our framework allows the accurate anticipation of outage scenarios caused by the gust-wind speed of hurricane and captures both pre-outage and post-outage optimal operation of distribution network. Numerical simulations are performed on the modified IEEE 33-bus distribution system with real hurricane data in the Houston shows that the proposed framework can quickly and accurately assess the effectiveness of using energy storage in boosting the grid resilience.

Future works include extending the proposed resilience assessment framework to distribution network planning problems for placement of energy storage units as well as decisions on adopting alternative resilience enhancement strategies, e.g. distribution automation.