NonStationary Random Process for LargeScale Failure and Recovery of Power Distributions
Abstract
A key objective of the smart grid is to improve reliability of utility services to end users. This requires strengthening resilience of distribution networks that lie at the edge of the grid. However, distribution networks are exposed to external disturbances such as hurricanes and snow storms where electricity service to customers is disrupted repeatedly. External disturbances cause largescale power failures that are neither wellunderstood, nor formulated rigorously, nor studied systematically. This work studies resilience of power distribution networks to largescale disturbances in three aspects. First, a nonstationary random process is derived to characterize an entire life cycle of largescale failure and recovery. Second, resilience is defined based on the nonstationary random process. Close form analytical expressions are derived under specific largescale failure scenarios. Third, the nonstationary model and the resilience metric are applied to a real life example of largescale disruptions due to Hurricane Ike. Real data on largescale failures from an operational network is used to learn timevarying model parameters and resilience metrics.
Nonstationarity, resilience, nonhomogeneous Poisson process, mixture model, real data
I Introduction
Our power grid is a vast interconnected network that delivers electricity from suppliers to consumers. At the edge of the grid lies power distribution networks [1]. Power distribution networks provide medium or low voltages to residence and organizations, and thus serve a unique role of connecting the grid to end users.
Distribution networks consist of ”leaf nodes” of the energy infrastructure and are thus susceptible to external disturbances. For example, natural disasters repeatedly cause devastating destructions and service disruptions to distributions networks [2][3]. There were 6 major hurricanes and more than 10 snow and ice storms occurred in north America in the past 5 years [4]. Each natural disaster disrupted power services to more than 500,000 customers for days [4]. Largescale power outages are more prevalent and damaging in developing countries, whose energy demands are rapidly increasing but power infrastructures are still in development [5].
A fundamental research problem pertaining to the real problem is the resilience of power distribution to largescale external disruptions. Resilience corresponds to the ability of distribution networks to withstand external disturbances and to recover rapidly from failures. Up to date, tremendous attention has been directed to resilience of the core that consists of major power generation and transmission systems of high voltages [6][7][8]. As an external disturbance often affects power distribution networks in a wide geographical span, the service disruptions usually remain local at the edge [9]. The resilience of the edge, however, is understudied [10]. As the demand on energy is growing, the edge of the grid is becoming thicker. For example, a utility provider often serves millions of customers in America. Damages from external disturbances on power distribution can thus result in a profound impact to a large number of users. Hence, the issue of resilience of power distribution networks is much needed to be investigated.
Empirical studies have been conducted on how to access damages from largescale power outages (see [11] and reference therein). Monitoring systems have been developed and used in power industry to respond to failures due to natural disasters (see [12] as an example). These empirical methods predict the degree of damages, e.g., the maximum number of outages upon a natural disaster [2], often through observations and experiences. Quantifiable approaches are lacking and needed to leverage practical experiences to general principles on resilience of power distribution [5][13][10]. In other words, resilience, as a concept about an overall power distribution network, needs to be learned systematically from real data.
Two research challenges emerge. One is what to learn from real data. Real data consists of samples on responses of a distribution network to external disturbances. External disturbances such as hurricanes often occur suddenly and unpredictably. The resulting largescale failure and recovery at distribution networks thus exhibit random and dynamic behavior. For example, recovery depends on multiple random factors such as environmental conditions, available resources and preparation. Prior work models network failures and recoveries using finite state Markov processes [14]. These models belong to the general birthanddeath process [15] which include randomness but assume stationary failure and recovery. Nonstationary failure such as a nonhomogeneous Poisson point process is provided in the context of queues [16][17] with time varying Poisson parameters. General timedependent infiniteserver arrival/service processes are studied in [18][19]. However, nonstationary recovery is rarely included in network resilience. In fact, it is an open issue how to characterize largescale nonstationary failure and nonstationary recovery for power distribution networks.
The second challenge is how to define and characterize resilience for an entire nonstationary lifecycle of failure and recovery. Prior work provides general discussions that resilience should include multiple attributes, such as both failure and recovery [20][21][22]. However, most prior works mainly study resilience in terms of maintaining services, and thus consider failures only [23]. In fact, recovery is rarely included in defining resilience of power distribution networks. It is an open issue how to derive resilience metrics for characterizing nonstationarity in both failure and recovery.
Hence it is necessary to formulate, from ground up, an entire life cycle of largescale failure and recovery. This can prevent us from choosing a model and/or a learning approach subjectively. We first develop a problem formulation of largescale failures at the finest level of network nodes based on temporalspatial stochastic processes. However, as an individual external disturbance results in one ”snapshot” of network responses, information (data) from one snapshot is insufficient to completely specify such temporalspatial model [24]. We thus derive temporal models of an entire distribution network by aggregating spatial variables. The resulting temporal process models an entire nonstationary lifecycle of largescale failures. The model applies to general failure process that can be dependent and with an arbitrary distribution. Two distinct recoverycharacteristics emerge from our model. One is infant recovery that reflects the ability to recover rapidly from failures. The other is aging recovery that corresponds to prolonged failures. We define the resilience as the probability of infant recovery for a power distribution network. We then derive analytical expressions for special cases of failure and recovery.
The model and the resilience metric are studied in a real life example of large scale service disruptions of power distribution. Power failures occurred during a major natural disaster, Hurricane Ike in 2008. Pertinent resilience parameters are learned using real data from an operational distribution network.
The rest of the paper is organized as follows. Section II provides background knowledge and an example of largescale failures at power distribution networks. Section III develops a problem formulation of failureandrecovery processes. Section IV characterizes nonstationary failure and recovery individually and jointly; then derives analytical expressions for special cases of failure and recovery. Section V defines network resilience based on the nonstationary model. Section VI learns pertinent resilience parameters using largescale real data. Section VII discusses our findings and concludes the paper.
Ii Background and Example
A power distribution network is at the edge of the grid from substations to users. A power distribution network consists of components including substations, feeders, transformers, poles, and transmission lines, and meters.
A large number of such devices in a distribution network are often in the open, and thus susceptible to natural disasters such as hurricanes, ice and snow storms. For example, a fallen pole can cause a short circuit, and other devices to fail subsequently. An external disturbance such as a hurricane can cause a large number of failures. The failures interrupt electricity service to end users. A largescale external disruption can affect one or many distribution networks in a wide geographical area. Failure recovery is often done by dispatching crews to the field. Promptness of failure recovery thus depends on environmental constraints, preparedness, and resources.
To gain intuition on an entire life cycle of failure and recovery, we consider an example of largescale power failures occurred during Hurricane Ike. Hurricane Ike had a landfall in 2008 and affected densely populated areas in Texas and Louisiana. Figure 1 shows a histogram on failure occurrence time and duration at an operational distribution network before, during and after the hurricane. The example provides the following observations:
(a) Both failure occurrence and recovery are timevarying, i.e., nonstationary.
(b) Samples on failure occurrence time and duration are not identically distributed. Instead, recovery time characterized by failure duration is different for failures occurred at different time.
This shows that failure occurrence and recovery are statistically dependent, and nonstationary.
Iii Stochastic Model
We formulate largescale failure and recovery based on nonstationary random processes. We begin with the detailed information on nodal statuses in a distribution network. We then aggregate the spatial variables of nodes to obtain temporal evolution of failure and recovery of an entire network.
Iiia Failure and Recovery Probability
A temporalspatial random process provides a theoretical basis for modeling largescale failures at the finest scale of nodes. The temporal variable is time . The spatial variable is the location of a node in the grid. Here, nodes can be any components in a distribution network such as substations, feeders, hubs, transformers, transmission lines, and power circuits. A shorthand notation is used to specify both the index of node and its corresponding location, where for a power distribution network with nodes.
Let be the status of the th node at time for . if the th node is in a failure mode. if the node is in normal operation. For a distribution network, an example of a node failure includes a failed circuit, a fallen pole, a broken link, and an nonoperational substation.
Failures caused by external disturbances exhibit randomness. Whether and when a node fails is random. Whether and when a failed node recovers is also random. Hence, random processes can be used to characterize failure and recovery for all nodes in a network.
Given time , characterizes the probability that node is failed in the near future , where is a small time increment. Assume a node changes state, i.e., from failure to normal and vice versa. Then for the th node, ,
(1) 
Eq.1 provides a model for an individual node in a network. The model includes Markov temporal dependence and spatial dependence among nodal statuses. Such a model can be applied to a heterogeneous grid where nodes experience in general different failure and recovery processes. There are such temporalspatial equations for nodes in a distribution network. The equations together form a temporalspatial model for a network.
IiiB Temporal Process
When largescale outages are caused by individual external disturbances, information available is from ”snapshots” of temporal spatial network statuses. A snapshot corresponds to spatialtemporal nodal statuses with respect to one external disturbance, e.g., one hurricane. As there are usually only a few such snapshots available, information is insufficient for specifying a complete temporal spatial model at the node level. However, spatial variables can be aggregated out from Eq.1, making the information sufficient for temporal characteristics, where
(2) 
Furthermore, the probability can be related to an indicator function,
(3) 
Then we can define a temporal process as follows.
Definition: A temporal processes is a special case of the temporalspatial process where the spatial variables (’s) are aggregated for all nodes in a network. is the number of nodes in failure state at time ,
(4) 
where is an indicator function, i.e., if event A occurs, and otherwise.
Combining Equations 2, 3, and 4, we have,
(5) 
where . Hence, an expected increment of the number of failed nodes in a network equals to the total change of the aggregated probabilities on nodal statuses.
An increment in the total number of nodes in failure state can result from either newly failed or newly recovered nodes. To further characterize the temporal process , we define a failure process and a recovery process respectively.
Definition: Failure and recovery processes: Failure process is the number of failures occurred up to time . Recovery process is the number of recoveries occurred up to time .
Assume is sufficiently small so that at most one failure occurs during , then an increment on the number of failures satisfies
(6) 
where . Similarly, for a sufficiently small , it can be assumed that at most one recovery occurs during . Then,
(7) 
where . Hence, Eq.2 is simplified as,
(8) 
Furthermore, we assume at time , , , and . Aggregating increments in Eq.8 from to , we have,
(9) 
Hence, the expected number of nodes in failure state equals to the difference between the expected failures and the expected recoveries. This characterizes how the status of a distribution network changes from the present to the near future.
For practicality, the empirical processes as the sample means , , and can be used to estimate the true expectations , , and , respectively. Eq.9 can then be represented by the empirical processes,
(10) 
The empirical processes and the failureandrecovery equation allow learning from field data, which shall be elaborated in Section VI.
Iv NonStationary Failure and Recovery
We now focus on the temporal processes to derive nonstationary characteristics on failure and recovery individually and jointly. We derive close form expressions for special cases of failure and recovery processes. Our study reveals pertinent parameters which shall be used to quantify resilience in Section V.
Iva Failure Process
Let be the intensity function of the failure process. is the expected number of new failures per unit time at epoch , i.e.,
(11) 
is also referred to as the rate function of the failure process . The larger is, the more failures occur in a unit time duration. Hence, failure rate quantifies the intensity of failure occurrence. An nonstationary failure process has a timevarying intensity function . Assuming a failure process begins at , we have
(12) 
The probability density function of failure occurrence time at can be obtained from , where
(13) 
As a special case of a general failure process, , as the first moment, can completely determines the failure process . Consider an independent failure process where the number of failures is a counting process with independent increment . Among many such random processes, nonhomogeneous Poisson process (NHPP) [15] captures the nonstationary nature of failures in a parametric form. The parameter is the intensity function, . The nonstationarity refers to a common characteristic of largescale external disruptions where power failures occur at different intensity at different time. The definition of a NHPP, applied to a failure process, is provided below.
Definition [15]: A continuous time counting process is a NonHomogeneous Poisson Process with a timevarying rate function , if satisfies the following conditions:

;

has independent increments;

;

.
In Section VI, we shall show using real data that a failure process from a hurricane follows a NHPP with as the failure rate function.
IvB Recovery Process
We now define intensity function for a recovery process. is the expected number of new recoveries per unit time at epoch ,
(14) 
is also referred to as the rate function of the recovery process .
An nonstationary recovery process has a timevarying intensity function, i.e., . Assuming the temporal failure process begins at , we have
(15) 
The recovery rate characterizes how rapidly recovery occurs, which is measured by failure duration . For an nonstationary recovery process, a failure duration depends on when a failure occurs as illustrated in Figure 1. Such nonstationarity of recovery is characterized by which is a conditional probability density function of failure duration given failure time . For a given threshold , the conditional probability that a duration is bounded by for failures occurred around time is
(16) 
where is a random failure duration.
When is sufficiently small, this probability characterizes rapid recovery that occurs shortly after failures. For a given , the larger is, the more dominating the rapid recovery is. Given desired value of probability , the smaller is, the more dominating the rapid recovery is.
Rapid recovery is referred to as infant recovery. This terminology is borrowed from infant mortality in survivability analysis [25]. Infant recovery is a desirable characteristic of the smart grid. In contrast, slow recovery is referred to as aging recovery in analogous to aging mortality [26]. Infant and aging recovery can be formally defined as follows.
Definition: Infant and aging recovery: Let be a threshold value. If a node remains in failure for a duration less than ; a recovery is an infant recovery. Otherwise, the recovery is aging recovery. Infant recovery is characterized by . Aging recovery is characterized by .
Note here is a function of failure occurrence time. As we shall show through a reallife example in Section VI, failures occurred at different time may experience infant and/or aging recovery of different degrees, showing the nonstationarity of a recovery process.
IvC Joint FailureandRecovery Process
A joint failure and recovery process characterizes an entire life cycle of a failureandrecovery process (FRP), and represents the total number of nodes in failure state at time (Eq.4). The expected number of nodes in failure can be written in terms of rate functions,
(17) 
Failureandrecovery process can be viewed as an alternative form of the birthanddeath process [15]. However, commonlyused birthanddeath processes have a stationary distribution of failure duration and assume independence between failure occurrence and failure duration . Here, these two assumptions are not needed. This implies that failures occurred at different time can last different duration. For example, under strong and sustained hurricane wind, failures that do not happen in daytoday operation can occur due to, e.g., falling debris and power lines. We shall further elaborate this through the reallife example in Section VI.
A failure process and a recovery process are related by failure durations. As a special case when failure follows a nonhomogeneous Poisson process, the following theorem shows that the recovery process is also a nonhomogeneous Poisson process parameterized by recovery rate . In addition can be expressed by failure rate with the help of the distribution of failure duration .
Theorem
Assume
(a) Independent failure occurrence obeys a nonhomogeneous Poisson process with intensity function ;
(b) Failure duration follows a conditional probability density function for , .
Then the recovery time is drawn from a nonhomogeneous Poisson process with recovery intensity function,
(18) 
The proof of the theorem is given in Appendix A.
Hence, if a failure process is nonhomogeneous Poisson, our general formulation of the temporal random process reduces to a queue. indicates timedependent Poisson failures. indicates nonstationary distribution of failure duration. is for the infinite number of servers for repair; thus recovery can occur right after failure. Furthermore, when is stationary, reduces to the convolution of the failure rate and service time distribution . The model reduces to queue developed in [17].
IvD Special Cases
We consider two special cases where failure, recovery and joint processes exhibit simple analytical expressions. These expressions provide insights to an entire nonstationary life cycle of failureandrecovery process under different failure scenarios.
Case 1: Failure and recovery in daytoday operation:
When there are no significant external disruptions, power outages are assumed to occur randomly and sporadically. The failure intensity remains at a constant level, i.e., , where does not vary with time. The number of failures thus follows a homogeneous Poisson process in daytoday operation. The recovery rate in daytoday operation in Eq.18 reduces to,
(19) 
As , . This suggests that the recovery process in daytoday operation is also a homogeneous Poisson process at a long time horizon. Then the expectation of joint failureandrecovery process in daytoday operation is
(20) 
Case 2: Surging failure during a natural disaster:
Now consider a disaster scenario where failures occur suddenly and intensely. increases from a small value to a large value in a short time duration . The failure rate can then be written as,
(21) 
where is the unit step function, . This corresponds to an extreme case when a disaster causes sudden failures at time and then weakens right after.
in Equation 18 becomes,
(22) 
When , . Furthermore, when is sufficiently small, i.e., a surge of failures upon the disaster only lasts a short time,
(23) 
When the peak of surging failures , the recovery rate in Eq.23 is approximately proportional to the surging failure rate and distribution of recovery time immediately after the disaster. In the long run, reduces to .
Substituting Eq.21 and 23 to Eq.17, we have the expectation of the failureandrecovery process in surging failures as,
(24) 
where is the conditional cumulative density function (cdf) of failure duration given time .
Additional characteristics of recovery can further reduce the above expressions. When infant recovery completely dominates a recovery process, when . Hence, for an impulselike surge of failures at , the resulting recovery from the surge lasts duration, i.e.,
(25) 
Eq.24 in dominating infant recovery becomes,
(26) 
Here failure duration is assumed to be larger than the end of failure process to simplify the expression. Eq.25 and Eq.26 show that when a recovery process consists of only infant recovery, all failures due to disasters recover by after the failure eruption. After the distribution network resumes daytoday operation.
When aging recovery dominates the recovery process, for . This implies that recovery begins with delay after a surge of failures. The corresponding recovery rate reduces from Equation 23, where
(27) 
Eq.24 in dominating aging recovery becomes,
(28) 
Here we also assume . Eq.28 shows another extreme case where recovery does not begin until delay from the failure eruption. Then the failures start to recover slowly. Hence, aging recovery shows the difficulty in resuming service to customers, and is thus an undesirable characteristic for the smart grid.
In the general failureandrecovery process, recovery begins as soon as failures occur. During the disaster, the failure process dominates and increases rapidly. Afterwards, the recovery process dominates. At the large time scale and when disaster lasts for a short period of time, the overlap between the failure process and recovery process can be neglected. Thus the failureandrecovery process splits into individual processes, failure and recovery, respectively.
V Resilience
We now derive network resilience using the pertinent parameters for an entire life cycle of nonstationary failure and recovery.
Va Definition
An intuition on how to define resilience results from the previous section. Resilience should be characterized by virulence of failures , speed of recovery and threshold . These three parameters determine infant recovery. A distribution network experiences a combination of infant and aging recovery in general. A network is intuitively more resilient when exhibiting more infant recovery. Hence we define the resilience as the probability of infant recovery.
Definition: Resilience: Given a threshold value on failure duration, the resilience of a power distribution network is defined as
(29) 
At a highlevel, measures the resilience of the grid from largescale disruptions, and reflects the ability for a distribution network as a whole to survive largescale failures. Here we use rather than time dependent conditional probability . This is because resilience, when viewed as a networkwide quantity, should include an entire lifecycle of failure and recovery but not depend on specific time epoch .
VB Resilience Parameters
Resilience can be expressed explicitly by the three pertinent parameters as follows,
(30) 
where , and is the probability density function of failure time given in Equation 13. Eq.30 results in an alternative expression for resilience. There, is the expected value of the probability of infant recovery averaged over the nonstationary failure process. This expression shows explicitly how resilience is determined by the three pertinent parameters.
VC Threshold
is one other pertinent parameter that determines the resilience. For a given value of , if infant recovery dominates the recovery process over the entire failure duration, . If aging recovery dominates the recovery process over the entire failure process, . Thus, a larger represents a more resilient power distribution network. When , the infant and the aging recovery take equal weight.
How to determine the value for ? can result from practical considerations or customer requirements. For example, when recovering within hours is regarded as acceptable for a disaster scenario, defines infant recovery.
A value for can also be determined more objectively. An intuition results from the fact that if a network is dominated by infant recovery, most failure durations are less than . Hence, the slope of is relatively large for and relatively small for . This implies that corresponds to a pertinent change point of the slope of , i.e., a deep valley in the second derivative . This is illustrated in in Fig.2(a). In contrast, when a network is dominated by aging recovery, most of the failure durations are larger than . Thus the slope of is small for and large for . then corresponds to a positive peak in as illustrated in Fig.2(b). Hence, is determined by the largest magnitude of the second derivative of the ,
(31) 
We shall provide an example using real data in Section VI.
Vi LargeScale Outages due to Hurricane Ike
We now apply the above framework on nonstationary failurerecovery processes to a reallife example of largescale utilityservice disruptions caused by a hurricane. Our focus is on using real data to learn the three resilience parameters , and and then to estimate the resilience of an operational power distribution network.
Via Real Data and Processing
Hurricane Ike was one of the strongest hurricanes occurred in 2008. Ike caused large scale power failures, resulting in more than 2 million customers without electricity, and marked as the second costliest Atlantic hurricane of all time [27][28].
Reported by National Hurricane Center [29], the storm started to cause power outages across the onshore areas in Louisiana and Texas on September 12, 2008 prior to the landfall. Ike then made a landfall at Galveston, Texas on 2:10 a.m. Central Daylight Time (CDT), September 13, 2008, causing strong winds, flooding, and heavy rains across Texas. The hurricane weakened to a tropical storm at 1:00 p.m. September 13 and passed Texas by 2:00 a.m. September 14.
Widespread power outages were reported across Louisiana and Texas starting September 12 [28]. A major utility provider collected data on power outages from more than ten counties. The outages include various component failures in a distribution network such as failed circuits, fallen poles, and nonoperational substations. The raw data set consists of 5152 samples. Each sample consists of the failure occurrence time () and duration () of a component () in a distribution network. The accuracy for time is a minute. The data set contains failures occurred from September through , 2008.
The 5152 samples on the failure occurrence time are plotted in Fig.3, where the length of each bin is two hours. There were significantly more failures occurred from 7 a.m. September 12 to 4 a.m. September 14, during which Hurricane Ike made the landfall in Texas. There are 2005 samples in this time period. Those 2005 samples are considered as failures due to Hurricane Ike.
Furthermore, the data set contains groups of failures that occurred within a minute. As a minute is the smallest time scale for each sample, the groups are considered as dependent failures. Dependent failures are grouped as one failed entity (), with a unique failure occurrence time and duration . After such preprocessing, the resulting data set has 465 failed entities. Two outliers with negative failure duration are further removed. The remaining 463 failed entities from 7 am September 12 to 4 am September 14 are referred to as nodes. is the data set we use for the rest of the paper.
ViB Empirical Failure Process
We now focus on studying the empirical failure process using data set.
ViB1 Learning failure rate
First, we use the data set to learn failure rate function . The empirical rate function is estimated using a simple moving average [30]: , where . The resulting rate function is overlaid with the samples on the number of failures in Figure 4, where each bin is of duration 1 hour, hours.
The failure rate function shows a timevarying, i.e., nonstationary, intensity of new failure occurrence, and can be described as follows.

Prior to 7 p.m. September 12, the intensity was low, i.e., fewer than 5 new failures occurred per hour. Hence where is considered as the failure rate in daytoday operation.

At 7 p.m. September 12, the intensity increased sharply first to 25 new failures per hour. In the next 6 hours, the intensity reached the peak value of nearly 50 new occurrences per hour. This is consistent to the weather report [29] that the strong wind about 145 mph and flooding impacted the onshore areas even prior to the landfall. The time when the peak occurred coincides with the landfall at 2:10 a.m 9/13 CDT.

After staying at the high level for about 12 hours (from 7 p.m. September 12 to 7 a.m. September 13), the intensity decreased rapidly back to a low level of less than 5 new failures per hours.
ViB2 NonHomogeneous Poisson Model
We now consider a hypothesis that these failure occurrences are governed by a nonhomogeneous Poisson process (NHPP). We perform Pearson’s test [31] on hypothesis that the empirical failure process follows a nonhomogeneous Poisson Process (NHPP) with intensity function . The test is on two aspects: (a) independence of the outages occurred at the large time scale of minutes, and (b) the sample mean, i.e., the empirical rate function , is sufficient for characterizing the failure process.
We divide the time duration from 7 a.m. September 12 to 4 p.m. September 14 into 400 intervals. In each interval, the number of new failures is compared with the expectation. The sum of the square errors from all intervals results in a chisquare statistic. The chisquare statistic is , with a degree of freedom of 2. Given a confidence level at , a threshold value is obtained as , where . The chisquare statistic obtained from the data is below the threshold . Hence hypothesis is not rejected. The detailed procedure of Pearson’s test is given in Appendix B.
However, not rejecting is insufficient for accepting the hypothesis. The goodness of fit of NHPP to the data is further validated through QuantileQuantile (QQ) plot given in Figure 5. There, the samples are compared with an nonhomogeneous Poisson process with intensity function . The figure shows that the nonhomogeneous Poisson process with the learned intensity function indeed exhibits a good fit to the data. Based on the result from Pearson’s hypothesis testing and the QQ plot, these 463 power failures occur independently, and obey a nonhomogeneous Poisson distribution.
ViC Empirical Recovery Process
We now learn the empirical recovery process characterized by , the conditional probability density function of failure duration given failure occurrence time. Our objective here is to identify infant and aging recovery.
ViC1 Data
The samples on the failure durations in our data set correspond to the failures occurred from 7 a.m. September 12 to 4 p.m. September 14. These samples result in a joint empirical distribution depicted in the 3D Fig.1. Each bin is of length (initial failure time) of hour and width (duration) of hours. The height of each bin located at initial failure time and failure duration , represents the number of failures that occur at and last for .
ViC2 Mixture Model
Given failure time , we select a mixture model as the probability density function for duration ,
(32) 
where is the number of mixtures at time , () is a weighting factor for the th mixture function , and . All these quantities vary with failure time for a nonstationary recovery process.
A mixture model is chosen since its parameters exhibit interpretable physical meaning. A parametric family of Weibull mixtures is particularly appealing as the parameters correspond to infant and aging recovery directly [15]. Weibull distributions have been widely used in survival analysis [26][25] and reliability theory [15], but not in characterizing recovery from largescale disturbances and resilience parameters. Specifically, a Weibull distribution is
(33) 
where , and are the shape and scale parameters respectively. A shape parameter is especially important for determining the type of recovery signified by a mixture component. The smaller is, the faster the decay rate of , the shorter the failure duration and thus the faster the recovery. Hence, corresponds to infant recovery whereas corresponds to aging recovery. Weighting factor signifies the importance of a component . For a nonstationary recovery process, these parameters are varying with failure time .
ViC3 Learning Mixture Parameters
The parameters of the mixtures of Weibull distribution are learned from the data. For simplicity, we use a piecewise homogeneous function to approximate . We divide the failure time into intervals. Within an interval , is assumed to be a time homogeneous function that does not vary with failure time . For different intervals, ’s have different parameters for nonstationarity, where
(34) 
1  2  3  

0.486  0.257  0.257  
0.710  14.400  211.830  
1.000  10.533  10.679  
1  2  3  4  
0.321  0.206  0.019  0.454  
2.680  7.640  21.220  173.580  
0.370  2.910  46.230  3.090  
1  2  3  
0.143  0.472  0.385  
0.530  12.300  135.072  
2.500  15.201  4.424  
1  2  
0.323  0.677  
11.041  112.245  
5.310  12.398  
1  2  3  
0.273  0.159  0.568  
2.479  21.555  134.053  
0.987  1.702  5.070 
Referring to Fig.1, we divide the time into 5 intervals, , the boundaries of each interval are depicted by dashed lines in Fig.1.
Within each interval, the parameters of the mixture Weibull distribution are learned through maximum likelihood estimation [32], and shown in Table I. Each mixture represents one cluster of failure durations. For example, distribution of durations in is depicted in Fig.6, where we observe 3 clusters.
For different intervals, the failure duration exhibits distinct distribution. For example, (7 a.m. September 12 to 7 p.m. September 12) corresponds to the time before hurricane. In this interval, the network was not yet impacted by hurricane Ike and was under daytoday operation. Most failure durations were as short as a few hours. (7 p.m. September 12 to 3 a.m. September 13) corresponds to the time right before the landfall and hurricane began to cause largescale failures. In this interval, more prolonged failures occur and recovery became more difficult than daytoday operation.
ViD Overall FailureandRecovery Process
We now compare the empirical failureandrecovery process with the learned process . is obtained by directly adding up the number of failed nodes from the actual data. is obtained by reconstructing the temporal process with learned and through Eq.17. Fig.7 shows the comparisons. The dotted is obtained with the assumption that is stationary over time. The dashed is obtained with the piecewise stationary given in Table I. All estimated processes are able to capture the trend in the data. However, the stationary distribution of failure durations deviates significantly from the actual sample path . This shows that the piecewise stationary better approximates the actual failureandrecovery process.
ViE Resilience
To evaluate resilience, we first calculate the probability of infant recovery within each interval as shown in Table I. The resilience curve is then obtained by averaging these probabilities over failure time. Fig.8 depicts together with its second derivative . is a concave function. Hence, the threshold is obtained as,
(35) 
we compute the resilience
(36) 
Hence, of the recoveries occurred within 13 hours whereas occurred later than 13 hours. Such resilience close to 0.5 reveals the combined nature of infant and aging recovery in the distribution network. Specifically, indicates a slight dominance of aging recovery over infant recovery.
Vii Discussion and Conclusion
A nonstationary random process has been derived to model largescale failure and recovery of a power distribution network under an external disturbance. The nonstationary network failure and recovery are characterized by timevarying failure rate and probability distribution of recovery time. These two quantities, together with a threshold on recovery time, define network resilience as the probability of rapid recovery. Analytical expressions have been derived to characterize the nonstationary failure and recovery under extreme conditions.
Real data has been obtained from an operational distribution network for a largescale external disturbance, Hurricane Ike. Model parameters of the failurerecovery process have been estimated through nonstationary learning using the real data. Pearson’s statistical test is applied to validate the assumptions of the failure model. The resilience metric has been evaluated, and shown for the operational network that of failures recovered within hours whereas the remaining prolonged as long as 12 days. This provides a quantitative measure of the resilience of the distribution network and a baseline for improvement.
How to improve resilience to reduce failures needs to be further incorporated in the resilience metric. Nonstationary learning algorithms need to be further developed for online learning of timevarying data.
Viii Acknowledgement
The authors would like to thank Chris Kung, Jae Won Choi, Daniel Burnham and Xinyu Dai at Georgia Tech for data processing, Anthony Kuh at University of Hawaii for helpful discussions on distribution networks. Support from National Science Foundation (ECCS 0952785) is gratefully acknowledged.
Appendix A
Proof of Theorem: In order to prove that the recovery process is a Poisson process, we just need to show the increment of recovery process, i.e., recoveries occurs in , is an independent Poisson random variable. Now consider nonoverlapping intervals , …, . We say a failure is type , , if it recovers in the interval . The number of the type events equals to the number of recoveries occur in the interval . By Proposition 5.3 in [15], it follows that the number of recoveries occur in the interval , i.e., the increment of recovery process, is an independent Poisson random variable with mean,
(37) 
where is the probability that a failure is type . Thus, we prove that is a Poisson process.
Then we compute and show the nonhomogeneity. Consider interval and a failure occurs at time , where . The probability that this failure recovers during is . Eq.37 becomes,
(38) 
where . We can rewrite Eq.38 as,
(39) 
The righthandside is just (Eq.14). Hence, the recovery rate function is,
(40) 
Due to the nonhomogeneity of , is a timevarying function, so is also timevarying function. Hence, we show the nonhomogeneity of , which completes the proof.
Appendix B
Pearson’s Hypothesis Test: The hypothesis test is based on a chisquare statistic which compares the failure occurrence times with their sample mean. The details of testing that failure occurrence times are drawn from a NHPP are given here.
1. Compute the estimated failure rate . At time , count the number of failure occurrences in , then .
2. Divide the failure occurrence times into intervals. Count the number of failure occurrences in each interval. Let denote the number of occurrence in interval , where . Let .
3. Count , which is the observed number of intervals with failure occurrences, where .
4. Use the estimated to compute , which is the expected number of intervals with failure occurrences. , where is the th time interval, .
5. Compute the sum . is a chisquare statistic, with degree of freedom (number of independent parameter fitted). Since one parameter is fitted, the degree of freedom is .
6. Given a confidence level, for instance , we obtain a threshold value . The hypothesis is rejected if ; otherwise, cannot be rejected.
References
 [1] S. M. Kaplan, “Smart grid: Electrical power transmission: Background and policy issues.” Congressional Research Service, Tech. Rep., April 2009.
 [2] “Emergency situation report,” Department of Energy, Tech. Rep., 2011.
 [3] R. Albert, I. Albert, and G. L. Nakarado, “Structural vulnerability of the north american power grid,” Phys. Rev. E, vol. 69, p. 025103, Feb 2004.
 [4] List of power outages. [Online]. Available: http://en.wikipedia.org/wiki/List$_$of$_$power$_$outages
 [5] H. Rudnick, “Natural disasters their impact on electricity supply,” IEEE Power and Energy Magazine, vol. 9, no. 2, pp. 22–26, March/April 2011, 2011.
 [6] G. A. Pagani and M. Aiello, “The power grid as a complex network: a survey,” Power, pp. 1–29, 2011.
 [7] P. Hines, J. Apt, and S. Talukdar, “Large blackouts in north america: Historical trends and policy implications,” Energy Policy, vol. 37, no. 12, pp. 5249–5259, December 2009.
 [8] Office of the Manager National Communications System, “Supervisory control and data acquisition (scada) systems,” Communication Technologies Inc., Tech. Rep., October 2004.
 [9] U.S. Department of ENERGY. Distributed energy. [Online]. Available: http://energy.gov/oe/technologydevelopment/smartgrid/distributedenergy
 [10] A. Kuh, personal communication.
 [11] “Comparing the impacts of the 2005 and 2008 hurricanes on u.s. energy infrastructure,” Department of Energy, Tech. Rep., February 2009. [Online]. Available: http://www.oe.netl.doe.gov/docs/HurricaneComp0508r2.pdf
 [12] “Forces of nature come and go. entergy¡¯s preparation never stops,” Entergy, Tech. Rep., November 2011.
 [13] W. H. Hooke, “Engineering for the threat of natural disasters,” BRIDGE WASHINGTON NATIONAL ACADEMY OF ENGINEERING, vol. 37, no. 1, 2007.
 [14] Y. Liu and K. S. Trivedi, “A general framework for network survivability quantification,” in 12th GI/ITG Conf. Measuring, Modeling, and Evaluation of Computer and Communication Systems, 2004.
 [15] S. M. Ross, Introduction to Probability Models, 10th ed. Academic Press, 2010.
 [16] W. A. Massey, “Networks of infiniteserver queues with nonstationary poisson input,” Queueing Systems, vol. 13, no. 1, pp. 183–250, 1993.
 [17] S. G. Eick, W. A. Massey, and W. Whitt, “Mt/g/ queues with sinusoidal arrival rates,” Management Science, vol. 39, no. 2, pp. 241–252, 1993.
 [18] B. L. Nelson and M. R. Taaffe, “The queueing system: Part  the single node,” INFORMS Journal on Computing, vol. 16, no. 3, pp. 266–274, 2004.
 [19] ——, “The queueing system: Part  the multiclass netowrk,” INFORMS Journal on Computing, vol. 16, no. 3, pp. 275–283, 2004.
 [20] M. AlKuwaiti, N. Kyriakopoulos, and S. Hussein, “Network dependability, faulttolerance, reliability, security, survivability: A framework for comparative analysis,” in International Conference on Computer Engineering and Systems, Cairo, Egypt, November 57 2006.
 [21] R. J. Ellison, D. A. Fisher, R. C. Linger, and et al, “Survivable network system: An emerging discipline,” Software Engineering Institute, Carnegie Mellon University, Pittsburgh, PA, Technical Report CMU/SEI97TR013, November 1997.
 [22] X. Lin, R. Xu, H. Xiong, and et al, “A framework of quantitative analysis for information system survivability,” Journal of Electronics and Information Technology, vol. 28, no. 9, pp. 1721–1726, 2006.
 [23] S. C. Liew and K. W. Lu, “A framework for characterizing disasterbased network survivability,” IEEE J. Sel. Areas Commun., vol. 12, pp. 52–58, 1994.
 [24] K. Kant. Surviving large scale internet failures, tutorial. [Online]. Available: http://www.kkant.net/Tutorials/LargeFailures/
 [25] D. W. Hosmer and S. Lemeshow, Applied Survival Analysis: Regression Modeling of Time to Event Data, 2nd ed. WileyInterscience New York, 2008.
 [26] J. D. Kalbfleisch and R. L. Prentice, The statistical analysis of failure time data, 2nd ed. New York: John Wiley and Sons, 2002.
 [27] E. S. Blake and C. W. Landsea, “The deadliest, costliest, and most intense united states tropical cyclones from 1851 to 2010,” National Hurricane Center, Tech. Rep., August 2011.
 [28] J. Colley and S. M. DeBlasio Sr, “Hurricane ike impact report,” Governor¡¯s Office of Homeland Security, Tech. Rep., December 2008.
 [29] “A digital record of the complete best track data,” National Hurricane Center, Tech. Rep., 2008.
 [30] H. L. V. Trees, Detection, Estimation, and Modulation Theory. New York: Wiley, 1971, vol. 1.
 [31] R. L. Plackett, “Karl pearson and the chisquared test,” International Statistical Review / Revue Internationale de Statistique, vol. 51, no. 1, pp. 59–72, April 1983.
 [32] L. Le Cam, “Maximum likelihood ¡ª an introduction,” ISI Review, vol. 58, no. 2, pp. 153–171, 1990.