Learning GeoTemporal NonStationary Failure and Recovery of Power Distribution
Abstract
Smart energy grid is an emerging area for new applications of machine learning in a nonstationary environment. Such a nonstationary environment emerges when largescale failures occur at power distribution networks due to external disturbances such as hurricanes and severe storms. Power distribution networks lie at the edge of the grid, and are especially vulnerable to external disruptions. Quantifiable approaches are lacking and needed to learn nonstationary behaviors of largescale failure and recovery of power distribution. This work studies such nonstationary behaviors in three aspects. First, a novel formulation is derived for an entire life cycle of largescale failure and recovery of power distribution. Second, spatialtemporal models of failure and recovery of power distribution are developed as geolocation based multivariate nonstationary queues. Third, the nonstationary spatialtemporal models identify a small number of parameters to be learned. Learning is applied to two reallife examples of largescale disruptions. One is from Hurricane Ike, where data from an operational network is exact on failures and recoveries. The other is from Hurricane Sandy, where aggregated data is used for inferring failure and recovery processes at one of the impacted areas. Model parameters are learned using real data. Two findings emerge as results of learning: (a) Failure rates behave similarly at the two different provider networks for two different hurricanes but differently at the geographical regions. (b) Both rapid and slowrecovery are present for Hurricane Ike but only slow recovery is shown for a regional distribution network from Hurricane Sandy.
Nonstationarity, queuing model, mixture model, real data
I Introduction
Nonstationary modeling and learning have been widely applied to many applications [1][2]. This work contributes a new application in an emerging area of smart energy grid. The application is on learning from failure data how distributed power networks respond to external disturbances such as hurricanes. Learned knowledge provides understanding how power networks fail and recover in severe weather. Such understanding is a prerequisite of modernizing our power infrastructure.
Power distribution networks lie at the edge of the energy grid, delivering medium and low voltages to residence and organizations [3]. Distribution networks consist of leaf nodes of the energy infrastructure and are thus susceptible to external disturbances. For example, natural disasters cause widespread destructions and service disruptions to distribution networks [4][5]. There were about 16 major hurricanes and severe storms occurred in north America in the past 5 years [6], each of which disrupted electricity services from 500,000 to several million customers for days [6].
Existing approaches rely primarily on empirical approaches for largescale failures of power distribution. For example, empirical studies have been conducted on assessing damages from largescale power failures (see [7] and references therein). Monitoring systems have been used in power industry to respond to failures (see [8] as examples). As hurricanes and severe storms appear to occur frequently and at a largescale [6], empirical approaches become inadequate for real time failure assessment in a wide geographical area [9]. Furthermore, recovery from largescale power failures is even less understood. This is evidenced by how difficult it was for utilities to provide accurate recovery time to customers [9]. Overall, quantifiable approaches are lacking and needed for characterizing how power distribution networks respond to external disturbances. This is important for discovering and mitigating vulnerabilities for enhancing the power infrastructure [10][11].
Unique challenges emerge for quantifying how power distribution networks respond to largescale external disturbances. The first is randomness. External disturbances such as hurricanes exhibit random behaviors. The resulting power failures occur randomly also. The second is dynamic nature of failures and recoveries due to evolution of external disturbances. For example, a hurricane usually has a landfall with a strong force wind, and then gradually dies down when moving in land. Hence, nonstationarity (randomness and dynamics) is an intrinsic characteristic of largescale failures.
Nonstationary learning is a natural approach for quantifying nonstationary largescale failure and recovery of power distribution induced by external disturbances. However, an additional challenge for learning is lack of data. This may appear to be a paradox: A largescale external disturbance such as a hurricane often results in thousands of power failures, which amounts to a lot of data. However, in the space of external disturbances, a hurricane generates only one sample, i.e., a snapshot of network failures and recoveries from one external disturbance. Hence, data from an individual disturbance is valuable and should be used to enable learning. Note that using real data for studying largescale power failures and recoveries is not yet a common practice for the power infrastructure. Real data on power failures from external disturbance is rare [12][13]. A recent work shows the strength of combining algorithmic approaches with real data on geographically correlated power failures [14]. The focus there is on power transmission rather than distribution.
Incorporating all challenges, a basic research question we intend to answer is how to learn nonstationary behaviors of largescale failure and recovery for distributed power distribution, using real data from one external disturbance? Combining modelbased and datadriven methods is a viable approach for limited samples [15]. A model identifies pertinent quantities that determine nonstationary random processes of failure and recovery. We first derive a problem formulation to obtain a model. What remains unknown are model parameters, which can be learned from data. Such a combination of modelbased and datadriven approaches directs learning to a small number of functions or parameters, and thus makes effective use of data. In addition, a combination of modelbased and datadriven approaches makes learning explanatory: Learned model parameters bear physical meanings on how distributed power distribution responds to an external disruption.
Our formulation focuses on power failures and recoveries induced by exogenous weather. The time scale of such failures and recoveries is considered to be a minute to be consistent to that of a hurricane (see Section V for details). Power failures can also occur in bursts at a small time scale of seconds or less [16]. Such bursty failures are usually due to an internal network structure (see Section V) and not studied in this work. Selfrecoveries often occur at the small time scale of subseconds [16] whereas recovery by field crews occur in minutes or beyond. Hence our model at the time scale of a minute focuses on weather induced failures and recoveries that can not be repaired through selfhealing. Such a model provides understanding how distributed power infrastructure responds to external disturbances.
Our formulation begins with the spatial scale of network nodes and the temporal scale of a minute. As the data from an external disturbance is insufficient to completely specify a detailed temporalspatial model [17], we aggregate spatial variables into groups. A group can be a city that consists of nodes from a small geographical area. The resulting model thus characterizes an entire nonstationary lifecycle of largescale failure and recovery in time and at geolocations. Such a spatialtemporal model is multivariate generalization of queues [18] to include geolocations. ’s and ’s are arrival (failure) processes and departure (recovery) processes for individual geographical area. “” means that it is possible for recovery to occur immediately after a failure, e.g., less than a minute in this work. Hence, multivariate ’s and ’s constitute our model that completely specify nonstationary behaviors of largescale failure and recovery at a power distribution network.
We consider one simplified characterization of queues to the expected values [18]. What to learn then becomes clear: A small number of pertinent parameters of and at different geolocations, i.e., failure rates and recovery time distributions. We first obtain detailed data on largescale power failures from a real life example of a natural disaster, Hurricane Ike. Ike caused power failures in the south states of US and affected more than 2 million users in 2008. We devise learning for two scenarios using the real data. The first learns only temporal processes of nonstationary failure and recovery by aggregating over spatial variables of nodes in an entire network. The second learns geolocation based spatialtemporal processes by aggregating nodes in cities. We show the modeling facilitates learning where model parameters can be easily estimated using the failure data. We then apply the model to another data set from Hurricane Sandy. Hurricane Sandy caused widespread power failures to more than 8 million people in the northeast of US in 2012. The data set consists of aggregated rather than detailed power failures in one of the impacted areas. Our approach is shown to be applicable to the aggregated data for estimating failure and recovery rates. Our approach also shows what can not be learned using aggregated data.
In summary, the contribution of this work consists of the following: (a) a novel model based on nonstationary random processes and dynamic queues for weatherinduced largescale failure and recovery of power distribution, (b) simple learning approaches for estimating parameters of the nonstationary model, (c) applications of the model and nonstationary learning to real data from two hurricanes at different locations.
The rest of the paper is organized as follows. Section II provides background knowledge and an example of largescale failures at a power distribution network. Section III and IV develops a problem formulation of spatialtemporal nonstationary random processes. Section V describes the real data from Hurricane Ike and learns a geotemporal model. Section VI studies nonstationary failure and recovery using parts of real data from Hurricane Sandy. Section VII discusses our findings. Section VIII concludes the paper.
Ii Background and Example
We now provide examples on the temporal scale, and nonstationarity of failure and recovery.
Iia Time Scale of Failure and Recovery
We first discuss the time scale for modeling weather induced failures and recoveries. A power distribution network consists of components such as substations, feeders, transformers, power circuits, circuit breakers, transmission lines, and meters. An example power distribution system is illustrated in Figure 1, with a commonly used radial topology. Three types of components are shown for illustration: A primary substation, three secondary power sources, and loads. Links correspond to power lines. Assume that either a component or a link can fail during a hurricane. Assume that the substation is used as a primary source during normal operation. The secondary sources, that can be distributed renewable sources, are used for backup when the primary source fails [19]. Then the following scenarios can occur for failure and recovery:
(a) If all the sources fail due to an external disturbance, there is no electricity supply to any loads. Hence, the loads experience dependent failures that can occur instantaneously. The scenario of dependent failures also applies to other components upstream in a radial topology that cause loss of electricity at nodes downstream. Dependent failures are often experienced by loads within subseconds.
(b) If a link that connects a load to the network fails due to an external disturbance, there is no electricity supply to the load. Such link failures can occur independently due to fallen trees or power lines. Thus loads experience independent loss of electricity. As such independent failures are caused by exogenous weather, they are assumed to occur at a time scale of a minute or beyond. Such a time scale can be estimated through how rapidly a hurricane force wind passes a city. Consider a small city of acres as an example. Based on the IEEE standard (IEEE/ASTM SI 101997)[20], an approximated “diameter” of the city is about miles. Consider the speed of the force wind at miles per hour. It takes about minutes for the wind to pass the city. This provides a basis of using a minute as a time scale of weatherinduced failures.
(c) Recovery depends on the types of failures and recovery schemes. Certain failures can be repaired through selfrecovery [16]. For example, if the primary substation fails, the electricity supply to all loads can be recovered when the three secondary sources are in operation. In general, selfrecovery and automated reconfiguration built in power distribution usually operate at a time scale of subseconds or seconds [16]. However, failures due to external disturbances, e.g., falling trees and power lines, often require manual repair by field crews. Recovery time depends on not only restoration schemes but also environmental constraints, and is thus considered as random in this work. Such manual recovery time is in either minutes or hours or days from failures.
In summary, failures and selfrecoveries at a small timescale of seconds or subseconds depend on detailed network structure and selfrecovery schemes. Failure and recovery at a larger time scale of a minute and beyond are often due to external disturbances that evolve dynamically and randomly.
IiB Example of NonStationary Failure and Recovery
To gain intuition on an entire life cycle of failure and recovery of a distribution network, we consider a reallife example of largescale power failures occurred during Hurricane Ike in 2008. Figure 2 shows a histogram on failure occurrence time and duration at an operational distribution network before, during and after the hurricane. Each bin has length (failure occurrence time) of hour^{1}^{1}1CDT is used for all plots for Hurricane Ike. and width (duration) of hours. The height of each bin represents the number of failures that occur at time and last for duration . Figure 3 shows geographical distributions of failure occurrences at two different time epochs, where failure occurrence is evidently nonstationary across geographical regions. Hence,
(a) Failure occurrence is nonstationary, i.e., random and timevarying;
(b) Recovery time is nonstationary, i.e., obeys different probability distributions for failures occurred at different time;
(c) Failure occurrence and recovery time are also nonstationary spatially, i.e., exhibit different distributions for different geolocations.
Hence, samples on failure occurrence time and duration are not identically distributed but exhibit geotemporal nonstationarity.
IiC NonStationary Learning
Nonstationary random processes have been studied in the context of drifting concepts (see [21][22][23][24] and references therein). Samples for learning are dynamically drawn from a nonstationary environment. An issue arises on the sample size, i.e., whether data is sufficient for characterizing underlying drifts of distributions.
The problem of learning nonstationary processes in this work exhibits unique challenges in terms of sample size. For simplicity, batch data is assumed to be collected for learning an entire nonstationary life cycle of failure and recovery processes offline. A challenge here is that there is only one snapshot of a distribution network in space and time from one external disturbance. The number of data sets is often small, i.e., from a few severe storms. Therefore, combining modelbased and datadriven approaches becomes important, where data can be used to learn a small number of modelparameters from one external disturbance at a time [15]. In addition, combining modelbased and datadriven approaches for learning is required by the problem: Learned model parameters need to exhibit physical meaning for generic network behaviors upon external disturbances.
Iii Stochastic Model
We now formulate largescale failure and recovery based on nonstationary random processes. We begin with the detailed information on nodal statuses in a distribution system. We then aggregate the spatial variables of nodes to obtain temporal evolution of failure and recovery across geographical areas.
Iiia Failure and Recovery Probability
A geotemporal random process provides a theoretical basis for modeling largescale failures. The temporal variable is time that is assumed to be continuous at the scale of a minute. The spatial variable can be either geo or networklocation of a node. For simplicity, this work considers geolocation as a spatial variable to focus on locationbased failures induced by severe weather. Nodes can be components in a distribution system such as substations, feeders, hubs, transformers, transmission lines, and distributed energy sources. A shorthand notation is used to specify the index of node located at . for a power distribution network with nodes. An underlying network topology is assumed to be radial so that cascading failures occurred in mesh networks are not considered.
Let be the status of the th node at time for . We assume for simplicity that nodes only exhibit two states: if the th node is in a failure mode, i.e., without power supply. if the node is in normal operation. 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, , the probability that node stays in failure mode in is,
(1) 
Equation 1 assumes Markov temporal dependence, and can be applied to nodes in a distribution network. The equations together form a geotemporal model of a network. Note that statistically dependent failures at the small time scale less than a minute are not considered here, as such failures are often caused by an internal network structure rather than exogenous weather. Spatial dependence is embedded in the model but will be studied explicitly in subsequent work.
IiiB Aggregated GeoTemporal Process
When largescale failures are caused by one external disturbance, information available is from one “snapshot” of temporal spatial network statuses, and thus insufficient for specifying a complete temporalspatial model at the node level. Hence, nodes are aggregated over a geographical region (), resulting in
(2) 
Here , where is an indicator function. if event A occurs, and otherwise. We can define a geotemporal process as follows.
Definition: is a geotemporal process where the spatial variables (’s) are aggregated for all nodes in a predefined region . is the number of nodes in failure state at time located in ,
(3) 
Combining Equations 2 and 3, we have,
(4) 
where is an increment of the number of failed nodes in a certain region. is the result of either newlyfailed or newlyrecovered nodes. Hence, we define a failure process and a recovery process respectively.
Definition: 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 failure or recovery occurs at most once to a node during . The increments on a failure process and a recovery process satisfy respectively,
(5) 
where . Similarly, for a sufficiently small , it can be assumed that at most one recovery occurs during . Hence, Equation 2 is simplified as,
(6) 
Furthermore, we assume at time , , , and . Aggregating increments in Equation 6 from to , we have,
(7) 
Hence, the expected number of nodes in the failure state equals to the difference between the expected failures and the expected recoveries. We now group a distribution network of nodes into geographical regions , , based on their geolocations. A city, e.g., a subdivision, is an example of a geographical region widelyused by utilities. Then the failurerecovery process for the entire distribution network is defined as,
(8) 
where characterizes how local power distribution in region responds to an external disturbance.
Iv NonStationary Failure and Recovery
We now derive nonstationary characteristics on failure and recovery. Our derivation reveals pertinent quantities that completely model the behaviors of largescale power failures and recoveries in expected values. This is pertinent to learning a small number of parameters in Section V.
Iva Failure Process
A failure process can be characterized to the first moment by failure rate functions. Let , , , be a vector that consists of the rate function of a failure process, where is the expected number of new failures per unit time at epoch and region , ,
(9) 
The larger is, the faster failures occur in at time . is referred to as the rate function of the failure process . Hence, failure rate quantifies the intensity of failure occurrence. An nonstationary failure process has a timevarying intensity function across geolocations. Assuming a failure process begins at , we have , where
(10) 
for .
IvB Recovery Process
A recovery process can be characterized by recovery rate function , where . is the expected number of new recoveries per unit time at epoch and region ,
(11) 
An nonstationary recovery process has a timevarying rate function. Assuming the temporal failure process begins at , we have for ,
(12) 
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 and where a failure occurs as illustrated in Figure 2. Such nonstationarity of recovery is characterized by which is a conditional probability density function of failure duration given failure time at region . For a given threshold , the conditional probability that a duration is bounded by for failures occurred at time is
(13) 
When is sufficiently small, this probability characterizes rapid recovery that occurs shortly after failures. For a given , the larger is, the more rapid recovery dominates a recovery process. 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: 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 .
IvC Joint FailureRecovery Process
A joint failurerecovery process characterizes an entire life cycle of a failurerecovery process (FRP), and represents the total number of nodes in failure state at time in region (Equation 3). The expected number of nodes in failure can be expressed in rate functions,
(14) 
Failureandrecovery process can be viewed as a birthdeath process. However, commonly used birthdeath processes have a stationary distribution of failure duration and assume independence between failure occurrence and failure duration [27]. Here, these two assumptions do not hold. 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 falling debris and power lines. We shall further elaborate this through the reallife examples in Sections V and VI.
A recovery process is related to a failure process through a probability density function of failure durations.
Theorem Let be an independent increment (failure) process with a rate function , . Let be the duration of a failur occurred at time and region . has a conditional probability density function , where , . Then recovery rate satisfies
(15) 
where , with and being the failure time and recovery time respectively.
The theorem is a corollary of the Transient Little’s Theorem [18]. Intuitively, can be viewed as the probability that a failure occurred at time and region lasts duration. is the average number of failures per unit time recover after duration, i.e., the recovery rate by definition. Aggregating over all failures occurred prior to time results in Equation 15. The detailed proof is given in [28].
IvD What to Learn
What to learn now becomes apparent. Failure rate functions and probability density functions of recovery time completely specify our model to the first moment, i.e.,

, for ,

, for .
In general, the forms and the parameters of these two functions are unknown, and need to be learned from real data. The learned functions and the parameters can then be used to estimate the empirical processes. The empirical processes are the sample means , , and that estimate the true expectations , , and , respectively.
V Hurricane Ike
We first apply learning to a reallife example of largescale utilityservice disruptions caused by a hurricane.
Va Data From Hurricane Ike
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 [29][30].
Reported by National Hurricane Center [31], the storm started to cause power failures 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. (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.
A major utility provider collected data on power failures from more than ten cities. The failures include failed circuits, fallen poles and power lines, and nonoperational substations. The raw data set has of 5152 samples. Each sample consists of the failure occurrence time () and duration () of a component () in a distribution network from September through , 2008. The accuracy for time is a minute.
VB Data Processing
The data set contains bursts of failures that occurred within a minute. As a minute is the smallest time scale for each sample, the bursts 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 learning.
Spatial variables ’s can be either chosen a priori or through learning from data. In this work, we choose ’s to be small cities to include a natural living environment of customers and this method is widelyused by utility providers. There are cities in the data set as illustrated in Figure 7.
VC Temporal Failure Process
We first study the temporal nonstationarity of the failureandrecovery process. Spatial variables are aggregated across the entire network. This is equivalent to reducing multiple geographical areas to one entire impactregion from the hurricane. Then the geotemporal failurerecovery process reduces to a temporal process. For notational simplicity, spatial variables are omitted for temporal processes.
The empirical rate function is estimated using a simple algorithm based on moving average [32]: , where is chosen to be hours. 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.
The learned failure rate function shows a timevarying rate of new failure occurrence:
(a) Prior to 7 p.m. September 12, the rate was low, i.e., fewer than 5 new failures occurred per hour. Hence per hour is considered as the failure rate in daytoday operation.
(b) At 7 p.m. September 12, the rate increased sharply first to 25 new failures per hour. In the next 6 hours, the rate reached the peak value of nearly 50 new occurrences per hour. This is consistent to the weather report [31] that the strong wind about 145 mph and flooding impacted the onshore areas prior to the landfall. The time of the peak coincides with the landfall at 2:10 a.m 9/13 CDT.
(c) After staying at the high level for about 12 hours (from 7 p.m. September 12 to 7 a.m. September 13), the rate decreased rapidly back to a low level of less than 5 new failures per hours.
VD Temporal Recovery Process
We now learn the empirical recovery process characterized by , the conditional probability density function of failure duration given failure occurrence time . As the spatial aggregation removes the geolocation variables, is the conditional density function of failure duration of an entire network.
We use the samples on the failure durations and occurrences in our data set. These samples result in a joint empirical distribution in Figure 2. The height of each bin located at represents the number of failures that occur at time and last for duration . Figure 2 shows nonstationarity of failure durations. For example, a large number (217) of failures occurred between 7 p.m. September 12 and 8 a.m. September 13 lasted for more than a day. This indicates that many failures occurred during the surge of the hurricane were difficult to recover. Hence, a nonstationary distribution for is an appropriate assumption.
Given failure occurrence time , we observe that the distribution of duration is a combination of two components: Infant and aging recoveries. We thus select a mixture model for the probability density function where ,
(16) 
where is the number of mixtures at time , () is a weighting factor for the th mixture function , and . Weighting factor signifies the importance of the th component . For a nonstationary recovery process, these parameters vary with failure time .
A mixture model is chosen since its parameters exhibit interpretable physical meaning [33][34][17]. A parametric family of Weibull mixtures is particularly appealing as the parameters correspond to infant and aging recovery directly. Weibull distributions have been widely used in survival analysis [26][25] and reliability theory [27], but not in characterizing recovery from largescale external disturbances. Specifically, a Weibull distribution is
(17) 
where , and are the shape and scale parameters respectively. Hence, th component in Equation 16 is .
Shape and scale parameters, and , are pertinent for characterizing the type of recovery. The smaller and are, the faster the decay of , the shorter the failure duration and thus the faster the recovery. Hence, and moderate (e.g., or smaller) correspond to infant recovery. and large (e.g., ) correspond to aging recovery.
For simplicity, we use a piecewise homogeneous function to approximate . The failure time is divided into intervals shown in Figure 2. Within interval for , is assumed to be stationary that does not vary with failure time . For different intervals, ’s have different parameters for nonstationarity,
(18) 
The parameters of the Weibull mixtures within each interval are learned through maximum likelihood estimation [17] from the data. Failure durations obey different distributions for failures occurred at different intervals, showing the nonstationarity. For example, the first duration (7 a.m. September 12 to 7 p.m. September 12) is when the network was not yet impacted widely by Hurricane Ike. Three Weibull mixtures are learned from the data, with the shape, the scale and weighting parameters as , and . The first two components result in dominating infant recovery, where of failures recovered within a day. In contrast, the third duration (3 a.m. September 13 to 3 p.m. September 13) is when the largescale failures continued to occur after the landfall. Two Weibull mixtures are learned from the data. The shape, the scale and weighting parameters are and , showing dominating aging recovery. As the result, only of failures recovered within a day. The second duration (7 p.m. September 12 and 8 a.m. September 13) is around the hurricane landfall, where about a half of the failures occurred experienced infant recovery within a day (see Figure 5 for the three Weibull mixtures). For durations overall, the probability of infant recovery within a day changes over time, showing the nonstationary of failurerecovery processes.
We then reconstruct the empirical temporal failurerecovery process with learned and through Equation 14. Figure 6 shows the comparisons between and , the reconstructed and the actual sample paths of the failurerecovery process repectively. The closeness between the two sample pathes shows that the piecewise stationary approximates well the actual failureandrecovery process.
VE GeoTemporal Failure Process
We now incorporate geolocation variables to learn the geotemporal nonstationarity. Failure process is a geotemporal process with multiple attributes from geographical regions, . The empirical failure rate functions for are estimated using the same algorithm of moving average. The resulting rate vector is multivariate, consisting of timevarying functions. Due to the small sample size, there are 6 out of 13 cities shown in Figure 7, each of which has sufficient samples ranging from to . Figure 8 shows the failure rates of the 6 cities. The multivariate failure rates exhibit the following characteristics:
(a) Temporal nonstationarity: At a given geographical region , is a timevarying function similar to the bellshaped curve obtained for the entire network. Consider as an example. The failure rate was low (few than 5 failures) prior to 7 p.m. September 12. Then, the rate increased sharply and reached the maximum value of 25 new failures per hour, at about 1 a.m. September 13. After that, the rate decreased rapidly to few than 5 failures.
(b) Spatial nonstationarity: At a given time , is a spatiallyvarying function. The peak values of failure rates vary from 1.5 to 27 per hour across the 9 cities. The time when the rate reached the peak value varies between 8 p.m. September 12 to 7 a.m. September 13, and is depicted as a dashed line at the bottom in Figure 8.
(c) Spatial temporal nonstationarity: The regions are then labeled with respect to the order of failure rates that reached the maximum value in Figure 8. For example, the failure rate at City reached the peak value first, followed by the failure rates at City through City . The figure shows the geotemporal characteristic that failure rates at different city reached their peak values approximately from the coast to inland. This appears to be consistent to the movement of the hurricane track (Figure 7).
VF GeoTemporal Recovery Process
To learn the geotemporal nonstationary recovery, we extend the mixture model (Equation 16) to a geotemporal bivariate mixture, where for ,
(19) 
Again our learning focuses on the 6 cities with sufficient samples. Dependencies of failure durations among cities are not studied in this work because of the small sample size.
We apply the piecewise homogeneous distribution function in Equation 18 to each region ,
(20) 
Here, each component is a Weibull distribution . Mixture ’s and their coefficients vary with respect to not only failure occurrence time (temporal nonstationarity) but also geolocations ’s (spatial nonstationarity).
Applying the maximum likelihood estimation [17], we obtain the estimated parameters of Weibull distributions in the 6 cities. Note that due to the small sample size in some of the regions, the parameters of distributions of failure duration have to be assumed, in our implementation, not varying with failure occurrence time within a region. The probability of infant recoveries is also computed accordingly. Three cities (1, 4, 6) show a similar percentage of infant recovery from to whereas the remaining cities (3, 5, 8) have infant recovery from to . Table I shows the learned model parameters for two example cities. Figure 9 shows the geographical distribution of infant and aging recoveries for the 6 cities.
The probability of infant recovery as well as model parameters vary across different geographical regions, showing the spatial nonstationarity of the recovery process. Examining more details, adjacent cities (e.g., 1 and 3) that are close to the coast can exhibit different percentages of infant recovery. Faraway cities (e.g., city 8 which is far in land and city 5 which is close to the coast) can also exhibit a similar percentage of infant recovery. Hence, recovery processes seem to be complex and require further study.
Vi Hurricane Sandy
We now learn using real data from another reallife example of largescale disruptions caused by Hurricane Sandy. This provides an understanding how our model and learning approach can be generalized to other hurricanes.
Via Data
Hurricane Sandy had a landfall at Northeastern United States on October 28, 2012. Hurricane Sandy resulted more than 6 million customers without electricity for days. The state with the most customers without power was New Jersey, where about 1.98 million customers lost power supplies [9].
A utility company, reported the number of failures (outages) in more than 10 counties in New Jersey from October 28, 2012 to November 22, 2012. The aggregated number of reported outages is a sample in our data set. Each sample consists of a given geolocation and time at the scale of minutes (the reporting interval). The geolocation variable corresponds to a county in New Jersey for . The data set consists of such samples, i.e., for time from October 28 to November 22, 2012. Figure 11(a) plots the data. Note that such aggregated data does not provide accurate occurrence time nor duration of each power failure.
1  2  3  

0.3478  0.3188  0.3333  
0.0045  12.1893  197.0316  
0.2490  2.7891  3.7629  
1  2  3  
0.3000  0.1500  0.5500  
0.0650  12.2138  129.7408  
0.2897  3.9992  2.8037 
ViB Empirical Failure Process
Learning now begins with the aggregated number of failures for , from which failure and recovery rates are estimated accordingly. This is a reverse process to learning from detailed failure data in Hurricane Ike.
To learn the failure rate, we recall that from Equation 10, and from Equation 14. This suggests that a lower bound on the failure rate can be estimated from the aggregate number of failures at time as
(21) 
where is a time epoch when increases.
To determine how to obtain such an estimate, we examine characteristics of raw (time series) data at the county level. Figure 10 shows two examples of the number of aggregated failures at two different counties in New Jersey. shows sharp increases and sharp decreases. A sharp increase occurs when the failure rate exceeds the recovery rate whereas a sharp decrease happens when recovery rate exceeds the failure rate. Hence, a change point in can be used to identify a lower bound for either a failure rate or a recovery rate. In addition, a sharp increase/decrease indicates a salient rather than noisy change point, where a lower bound can be obtained accurately.
We first obtain the positive increments from for each region using Equation 21. We then aggregate the increments over the regions to obtain a lower bound for the failure rate of the utility network. , the estimated lower bound on the number of failures up to time , can then be obtained by integrating , which is shown in Figure 11(b) ^{2}^{2}2 EST is used for plots in regard to Hurricane Sandy..
ViC Empirical Recovery Process
To learn the empirical recovery rate, we apply Equation 21 except that corresponds to the time epoch of a decrease in the number of failures. Figure 11(c) shows an estimated lower bound for recovery rate and the cumulative number of recoveries respectively.
Since the aggregated data from Hurricane Sandy does not contain detailed recovery time for each failure, it is impossible to learn the timevarying distribution of failure duration . Nevertheless, the aggregated data can be used to estimate a stationary distribution of recovery time, i.e., . As the detailed information on failure duration is not available from the data, we consider a simple distribution with one Weibull mixture . Applying discrete samples to Theorem IVC, reconstructed recovery rate can be related with and as
(22) 
where minuets is the step size, and is the discrete time. Weibull parameters and are then estimated to minimize the estimation error . Figure 12 shows the estimated Weibull distribution, where the shape parameter and the scale parameter . The resulting stationary distribution of failure durations is then used to reconstruct a lower bound for the recovery rate. Figure 12 shows the estimated from the data set and the reconstructed . Reconstructed thus provides a profile on how the recovery varies with time.
Vii Findings and Discussions
Viia Findings
Learning from Hurricane Ike and Hurricane Sandy results in the following findings.
ViiA1 Failure process
Failure rates are timevarying for both Hurricane Ike and Hurricane Sandy. The corresponding failure processes are nonstationary in time and geographical regions. However, the failure rates exhibit different characteristics at the county level for Hurricane Ike and Hurricane Sandy: The failure rates for Hurricane Ike appear to vary gradually. However, the failure rates for Hurricane Sandy exhibit sharp changes, showing that failures occurred in groups ^{3}^{3}3The cause shall be sought for when more detailed data becomes available.. When aggregated over geographical regions, failure rates for both hurricanes exhibit similar characteristics, i.e., first rapidly increasing and then decreasing. More quantitative study is needed to further compare the failure processes for different hurricanes at different spatial scales.
ViiA2 Recovery process
Learned recovery rates from Hurricane Ike and Hurricane Sandy are both timevarying. For Hurricane Ike, the learned probability distributions of failure durations exhibit nonstationarity in time and geolocations, i.e., depend on when failures occur. Such distributions constitute both infant and aging recovery, as shown in Table I and Figure 9. The degree of infant recovery, however, is different at different cities. Three out of the six chosen cities recovered more rapidly then the rest. Failures with infant and aging recoveries are also interleaving in geolocations.
The recovery for the provider network from Hurricane Sandy shows a nearly steady rate of 7000 recoveries per hour. In addition, the estimated Weibull distribution of the failure duration exhibits stronger aging recovery than infant recovery. A lack of infant recovery for this utility provider may indicate that power distribution networks suffered virulent disruptions during Hurricane Sandy. The recovery can thus be difficult. Yet, detailed rather than aggregated failure data is needed for accurately estimating distributions of failure durations.
Note that failures and recoveries can occur simultaneously within a minute interval. That is why the amount of increase in is a lower bound of the actual failure rate . When the number of failures increased rapidly, e.g., from October 28 to October 31, recovery appeared to be minor. When the hurricane passed the area after October 31, recovery dominated. This is shown by the lower bounds of the failure and the recoveryrate in Figure 11 and 12.
ViiB Discussions
The type of available data is important for learning nonstationary behaviors of power distribution in response to external disruptions. The accurate failure data from Hurricane Ike characterizes an entire life cycle of failure and recovery processes. Data from Hurricane Sandy is aggregated and thus lack of exact information on individual failure occurrence and duration. Hence, learning is to infer failure and recoveryprocesses, which is a reverse process to that for Hurricane Ike. The 15minute sampling interval seems to be sufficient for estimating the lower bounds of failure and recoveryrates from Hurricane Sandy. The aggregated data is insufficient for characterizing a nonstationary distribution of failure duration but can be used to learn a stationary distribution as an approximation.
To deal with the small sample size, a rule of thumb is used where training samples should be several times more than parameters [17]. For Hurricane Ike, 20 or more samples seem to be sufficient for estimating temporal characteristics of failure and recoveryrates but insufficient when the spatial nonstationarity is studied. This suggests that the algorithm need to be enhanced, e.g., to identify spatial scales appropriate for aggregation.
Our model assumes an underlying radial topology, where failures can be considered as independent increments at large temporal spatial scales (minutes, cities). Detailed network configuration is yet to be included in our model. For example, topology and power flows [35][36] are two possible characteristics to be included for failures and recoveries. Failure and recoveryprocess at a small time scale of subseconds then need to be considered accordingly. A challenge is much increased complexity and innetwork measurements at temporal spatial scales.
Viii Conclusion
This work shows that nonstationary geotemporal random processes naturally model largescale failure and recovery of power distribution induced by hurricanes. In particular, multivariate geolocation based queues provide such nonstationary failure and recovery processes. The nonstationary failure and recovery can be completely characterized to the expected values by timevarying failure rate and probability distribution of recovery time across geographical regions.
Real data from two hurricanes have been used to learn failure and recovery processes. Learning detailed failure data from Hurricane Ike reveals that the failure process across different geographical regions follows a similar trend to that of the hurricane. However, the failure and recoveryprocesses exhibit different infant and aging recovery across geographical regions. Learning aggregated data from an impact area by Hurricane Sandy shows that our model can infer failure and recovery rates using aggregated data. The failure rates have more significant discrete components for Hurricane Sandy than for Hurricane Ike at geographical regions. The recovery process is dominated by aging recovery for one utility network from Hurricane Sandy but consists of a significant component of infant recovery for another utility from Hurricane Ike. This shows that model is indeed needed for general failure and recoveryprocesses in dynamic queues. Note that these findings are for power distribution through open rather than underground networks.
These findings call for subsequent research on how distributed power distribution are impacted by external disturbances. For example, power failures and recoveries are yet to be studied at all impact areas for Hurricane Sandy. Spatial temporal dependencies among power distribution networks at different geographical regions need to be studied explicitly. This requires combining detailed configurations of power distribution with the dynamic model. These studies shall provide further understanding on how to enhance the distributed power infrastructure.
Ix Acknowledgement
The authors would like to thank Chris Kung, Jae Won Choi, Daniel Burnham, Xinyu Dai and Michael Perez for data processing, Amanda Cox for providing parts of the data and helpful discussions, Anthony Kuh for helpful discussions on distribution networks, anonymous reviewers for valuable comments, and Associate Editors for helpful suggestions. Support from National Science Foundation (ECCS 0952785) is gratefully acknowledged.
References
 [1] C. Alippi, G. Boracchi, and M. Roveri, “A JustinTime Adaptive Classification System Based on the Intersection of Confidence Intervals Rule,” Neural Networks, vol. 24, no. 8, pp. 791–800, 2011.
 [2] O. Schabenberger and C. A. Gotway, Statistical Methods for Spatial Data Analysis. Chapman and Hall/CRC, December 2004.
 [3] S. M. Kaplan, “Smart Grid: Electrical Power Transmission: Background and Policy Issues,” Congressional Research Service, Tech. Rep., April 2009.
 [4] Department of Energy, “Emergency Situation Report,” Tech. Rep., 2011. [Online]. Available: http://www.oe.netl.doe.gov/emergency_sit_rpt.aspx
 [5] 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.
 [6] List of power outages. [Online]. Available: http://en.wikipedia.org/wiki/List_of_power_outages
 [7] “Comparing the impacts of the 2005 and 2008 hurricanes on u.s. energy infrastructure,” Department of Energy, Tech. Rep., February 2009.
 [8] “Forces of nature come and go. entergy¡¯s preparation never stops,” Entergy, Tech. Rep., November 2011.
 [9] U.S. Department of Energy, “Emergency Situation Reports: Hurricane Sandy,” Tech. Rep., October 2012.
 [10] 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.
 [11] W. H. Hooke, “Engineering for the Threat of Natural Disasters,” Bridge Washington National Academy of Engineering, vol. 37, no. 1, 2007.
 [12] D. Zhu, “Electric Distribution Reliability Analysis Considering TimeVarying Load, Weather Conditions and Reconfiguration with Distributed Generation,” Ph.D. dissertation, Virginia Polytechnic Institute and State University, 2007.
 [13] Y. Liu and C. Singh, “A Methodology for Evaluation of Hurricane Impact on Composite Power System Reliability,” IEEE Transactions on Power Systems, vol. 26, no. 1, February 2011.
 [14] A. Bernstein, D. Bienstock, D. Hay, M. Uzunoglu, and G. Zussman, “Power Grid Vulnerability to Geographically Correlated Failures Analysis and Control Implications,” arXiv:1209.1099v1, June 2012.
 [15] S. Geman, E. Bienenstock, and R. Doursat, “Neural Networks and the Bias/Variance Dilemma,” Neural Comput., vol. 4, no. 1, pp. 1–58, Jan. 1992.
 [16] M. Amin and J. Stringer, “The Electric Power Grid: Today and Tomorrow,” MRS Bulletin, vol. 33, pp. 399–407, Apr. 2008.
 [17] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, 2nd ed. Wiley, November 2000.
 [18] D. Bertsimas and G. Mourtzinou, “Transient Laws of NonStationary Queueing Systems and Their Applications,” Queueing Syst. Theory Appl., vol. 25, no. 1/4, pp. 115–155, Jan. 1997.
 [19] T. Perry, “Solar Sandy Project Brings Panels to the People,” IEEE Spectrum, Tech. Rep., November 2012.
 [20] “IEEE Standard for Use of the International System of Units (SI): the Modern Metric System,” IEEE/ASTM SI 101997, pp. i–, 1997.
 [21] A. Kuh, T. Petsche, and R. L. Rivest, “Learning TimeVarying Concepts,” in Proceedings of the 1990 conference on Advances in neural information processing systems 3, ser. NIPS3. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1990, pp. 183–189.
 [22] G. Widmer and M. Kubat, “Learning in the Presence of Concept Drift and Hidden Contexts,” Machine Learning, vol. 23, pp. 69–101, 1996, 10.1007/BF00116900.
 [23] J. Gama, P. Medas, G. Castillo, and P. P. Rodrigues, “Learning with Drift Detection,” in Advances in Artificial Intelligence  SBIA 2004, 17th Brazilian Symposium on Artificial Intelligence, Sao Luis, Maranhao, Brazil, September 29  October 1, 2004, Proceedings, ser. Lecture Notes in Computer Science, vol. 3171. Springer, 2004, pp. 286–295.
 [24] R. Elwell and R. Polikar, “Incremental Learning of Concept Drift in Nonstationary Environments,” IEEE Transactions on Neural Networks, vol. 22, no. 10, pp. 1517–1531, 2011.
 [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] S. M. Ross, Introduction to Probability Models, 10th ed. Academic Press, 2010.
 [28] Y. Wei, C. Ji, F. Galvan, S. Couvillon, G. Orellana, and J. Momoh, “Dynamic Resilience for Power Distribution: Modeling and Learning from LargeScale Failure and Recovery,” IEEE Journal on Selected Areas in Communications: Smart Grid Communications Series, October 2012, submitted.
 [29] 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.
 [30] J. Colley and S. M. DeBlasio Sr, “Hurricane ike impact report,” Governor¡¯s Office of Homeland Security, Tech. Rep., December 2008.
 [31] “A Digital Record of the Complete Best Track Data,” National Hurricane Center, Tech. Rep., 2008. [Online]. Available: ftp://ftp.nhc.noaa.gov/atcf/archive/2008/bal092008.dat.gz
 [32] H. L. V. Trees, Detection, Estimation, and Modulation Theory. New York: Wiley, 1971, vol. 1.
 [33] S. P. Chatzis and Y. Demiris, “Nonparametric Mixtures of Gaussian Processes With PowerLaw Behavior.” IEEE Trans. Neural Netw. Learning Syst., vol. 23, no. 12, pp. 1862–1871, 2012.
 [34] W. Fan, N. Bouguila, and D. Ziou, “Variational Learning for Finite Dirichlet Mixture Models and Applications.” IEEE Trans. Neural Netw. Learning Syst., vol. 23, no. 5, pp. 762–774, 2012.
 [35] M. E. Baran and F. F. Wu, “Network Reconfiguration in Distribution Systems for Loss Reduction and Load Balancing,” IEEE Transactions on Power Delivery, vol. 4, no. 2, pp. 1401–1407, 1989.
 [36] Y. Zhao, R. Sevlian, R. Rajagopal, A. Goldsmith, and H. V. Poor, “Outage Detection in Power Distribution Networks with OptimallyDeployed Power Flow Sensors,” to appear in IEEE Power and Energy Society General Meeting, 2013.