Distribution of blackouts in the power grid and the Motter and Lai model
Yosef Kornbluth, Gabriel Cwilich, Sergey V. Buldyrev, Saleh Soltan, Gil Zussman
aa r X i v : . [ phy s i c s . s o c - ph ] A ug Distribution of blackouts in the power grid and the Motter and Lai model.
Yosef Kornbluth , Gabriel Cwilich , and Sergey V. Buldyrev , Saleh Soltan , , and Gil Zusman Department of Physics, Yeshiva University, 500 West 185th Street, New York, New York 10033 Princeton University, Princeton, NJ, USA Amazon.com, New York, NY, USA Department of Electrical Engineering, Columbia University, New York, NY, USA (Dated: August 5, 2020)Carreras, Dobson and colleagues have studied empirical data on the sizes of the blackouts in realgrids and modeled them by computer simulations using the direct current approximation. Theyhave found that the resulting blackout sizes are distributed as a power law and suggested that thisis because the grids are driven to the self-organized critical state. In contrast, more recent studiesfound that the distribution of cascades is bimodal as in a first order phase transition, resulting ineither a very small blackout or a very large blackout, engulfing a finite fraction of the system. Herewe reconcile the two approaches and investigate how the distribution of the blackouts changes withmodel parameters, including the tolerance criteria and the dynamic rules of failure of the overloadedlines during the cascade. In addition, we study the same problem for the Motter and Lai model andfind similar results, suggesting that the physical laws of flow on the network are not as importantas network topology, overload conditions, and dynamic rules of failure.
I. INTRODUCTION
Cascading failures in the power grids continue to hap-pen in spite of efforts to make power grids more re-silient [1–3]. The standard criterion of resiliency is the N − N − et al. [4]suggested that the power grid is driven to the SOC stateby recursive upgrading of the power grid with constantlygrowing power demand, applying the N − N − II. POWER-LAW VERSUS BIMODALDISTRIBUTIONS
In this section we will review the conditions underwhich the distribution of cascade sizes in overload mod-els are power law or bimodal, and compare these overloadmodels with the topological models [22–29], in which theoutcome of the cascade depends only on the initial topol-ogy of the network and location of the initially damagedelements. Topological or overload models can be investi-gated for two different classes of initial attacks.In the first well studied class [26, 27, 30, 32, 33], thecascade is started by a massive attack, after which onlya fraction p of elements survive. For this class, an in-teresting problem is to study the behavior of the orderparameter S ( p ), the fraction of the surviving nodes at theend of the cascade, as a function of the initially surviv-ing elements p . Many topological and overload modelsexhibit a first order phase transition at p = p t , at whichthe order parameter S ( p ) exhibits a step discontinuity,or a second order transition at p = p c , at which the orderparameter is continuous but its derivative with respect to p exhibits a step discontinuity. For finite systems, near p = p t the distribution of the order parameter becomesbimodal, and the transition point p c is defined as thepoint at which the population of the two peaks is equal.For the second order transition the distribution of the or-der parameter is always unimodal. For some topologicaland overload models, the line of the first order transi-tions, p t , in a plane of two parameters may transforminto the line of second order phase transitions, p c at acertain value of the second parameter [27, 30, 31, 33].In the second class of initial attacks, such as in theabove mentioned discussion of power grids, the cascadeis started by the failure of a single element. Here, it isinformative to study the distribution of the “blackout”sizes, i.e., the number of nodes that failed during thecascade and how such a probability distribution dependson the model parameters. This type of initial attacksmay exist in two flavors.In the first flavor, the system first undergoes a cas-cade of failures caused by a massive initial attack with p > p t and the order parameter stabilizes at S ( p ). Thisis analogous to the first class of initial attacks. After-wards, the avalanche is started by the removal of one ad-ditional element. In the well-studied topological models,it is known [32, 33], that S ( p ) − S ( p t + 0) ∼ ( p − p t ) / , (1)where S ( p t + 0) is the limit of S ( p ) for p → p t abovethe step discontinuity. Because in the topological modelsthe order of removal elements does not play any role, theremoval of one node after the system is stabilized at p is identical to the initial removal of pN + 1 nodes, where N is the total number of elements in the system. In eithersituation, the average avalanche size under this initialconditions is s = N [ S ( p + 1 /N ) − S ( p )] ≈ ∂S ( p ) /∂p ∼ ( p − p t ) − / . Hence, the avalanche size diverges above thetransition point p t , with the divergence characterized bythe critical exponent γ = − /
2. In network theory thistransition is called a hybrid transition, but in fact this be-havior is completely equivalent to the mean field behaviornear the spinodal of the first order phase transition, suchas the behavior of the isothermal compressibility nearthe spinodal of the gas-liquid phase transition. This fol-lows from the fact that the isothermal compressibility isequivalent to the susceptibility in the Ising model, whichin turn is equivalent to the average cluster size at thepercolation transition. The divergence takes place onlyfor p > p t , while for p < p t the avalanches remain of finitesize or do not exist at all, since S ( p ) = 0 for p < p t . For p → p t + 0, the distribution of the avalanche sizes devel-ops a power law behavior with an exponential cutoff S ∗ : P ( s ) ∼ s − τ exp( − S/S ∗ ) with the mean field exponent τ = 3 / S ∗ ∼ ( p − p t ) − /σ with σ = − γ = − /
2, satisfy a usual percolationscaling relation γ = (2 − τ ) /σ , but with different γ and σ than in the percolation theory [36, 37].For the second flavor of the one-element-failure initialattack, which is typical for most power grid simulations,the system exists in a state characterized by some set ofparameters, but the parameter p describing the size ofthe initial attack is not applicable at all. In this case,both in topological and overload models, the cascadingfailures evolve as a branching process in which the failureof one element leads to the failure of b elements, whichdepend on the failed element in the topological models,or get overloaded due to the removal of the failed ele-ment in the overload models. The critical point of thesimplest branching process with a fixed distribution ofthe branching factor P ( b ) is defined by h b i = 1 [34]. Forthe mean-field variant of the fixed distribution model,the distribution of the avalanche sizes is a power lawwith an exponential cutoff P ( s ) ∼ s − τ exp( − s/s ∗ ), where s ∗ ∼ ( h b i − − /σ . Here σ = 1 / τ = 3 / γ = (2 − τ ) /σ = 1 and β = ( τ − /σ = 1, where γ governs the divergence of theaverage finite avalanche size for h b i → β governsthe probability of an infinite avalanche (i.e., failed net-work) for h b i >
1. These infinite avalanches, which areapproximated in the finite system by finite avalanchesthat distintegrate the network, and whose size fluctuatesaround µN , form the second peak of the bimodal dis-tribution. The finite avalanches, which do not scale asa finite fraction of the network, form a first peak of thedistributions, which is separated from the second peakby a huge gap because there are no avalanches of size s , such that s ∗ << s << µ (Fig. 1). Thus, this simplebranching process model predicts the existence of regimeswith bimodal ( h b i >
1) and power-law distributions with s=x -5 -4 -3 -2 -1 P ( s > x ) =1.10, N=100000=1.05, N=100000=1.00, N=100000=0.95, N=100000=0.90, N=100000=1.00, N=200000=1.05, N=200000=1.00, N=200000=0.95, N=200000=0.90, N=200000slope =-1/2 FIG. 1. Cumulative distributions of the avalanche sizesin the failure branching process with the Poisson branch-ing degree distributions, with the average branching factors h b i = 0 . , . , , . , . N = 100000 and N = 200000 nodes. One can see that the bimodal distributionof avalanches, characterized by a wide horizontal plateau ofthe cumulative distribution, emerges for h b i >
1. For h b i < h b i = 1, the distribution is a power law with afinite-system cutoff. The power law region with the exponent τ − . exponential cutoff ( h b i ≤ N − h b i small: Electriccompanies want to eliminate large catastrophic blackoutswhich exist for h b i >
1, but trying to save money on theinfrastructure, are willing to accept finite blackouts. Ob-viously, this compromise is achieved for h b i = 1 − ǫ . Inthe following sections we will check this failure-branchingprocess model by computer simulations. III. SIMULATIONS OF THE POWER GRIDS
In order to verify that the N − n p generator nodes which produce power, n c loads whichconsume power and n t transmitter nodes, which do notproduce or consume power but redistribute the flow be-tween several power lines. We assume that all genera-tors or loads produce or consume equal power. We usethe model, setting n p = n c = 1000 and n t = 8000. All N = n p + n c + n t nodes are randomly placed on a square ofedge L and periodic boundaries. In the RR model, eachnode is connected to its k nearest neighbors. In order tomake sure that each line connected to exactly k nodes,we make a list of all nearest neighbor pairs and start toselect pairs from this list in ascending order of length, cre-ating a link between them if both of the nodes in the pairhave less than k links. At the end of the process, the ma-jority of nodes have k links with a few exceptions whichhave only k − r = L p h k i / ( πN ) andis connected to all the nodes within this circle. Sincethe nodes are randomly distributed geographically, thedegree distribution of nodes becomes a Poisson distribu-tion with average degree h k i . We used k = 4 for RR and h k i = 5 . N − n L links, we obtain the current in a given link ij after theremoval of each one of the other ( n L −
1) links, and we findthe maximal current I ∗ ij through that link ij obtainedover the set of those ( n L −
1) currents. We use this value I ∗ ij as the maximal possible load that the link ij may sus-tain. We start each simulation with the removal of a ran-dom pair of links, i a j a and i b j b , and recompute the cur-rents in all remining links, I ij . If in some links I ij > I ∗ ij ,we find the link with maximal overload I ij /I ∗ ij , remove itand find the new currents I ij (one-by-one update rule).We continue this process until, after removal of n links,for all remaining links, I nij ≤ I ∗ ij . If at a certain stage,the grid splits into several disconnected clusters, we ap-ply the power production-consumption equalization foreach cluster using the minimal production-consumptionrule described in [13] We compute the blackout size as afraction of the consumed power lost in the cascade.In order to carefully assess the importance of the N − I ∗ ij = ( α + 1) I i,j , where I i,j is the current in line i, j in the unperturbed network. Thus, for each networkmodel, we study four cases: N − N − N − α . Thisis not surprising because for all four models of the powergrid, the effective tolerance of line ij under the N − α ij ≡ I ∗ ij /I ij − I ij is the initial current inthe line ij ) is a wide distribution with almost 5% of lineshaving α i,j >
1. Thus, the fraction of runs that resultin any additional line overload is about 0.03; for mostpairs of initially attacked lines, no further lines are over-loaded. For the tolerance model with a uniform α , thislevel of protection is achieved only for α >
1. Thus, im-plementation of the N − h k i = 5 . h k i = 4 . N − h k i = 5 .
0, for which the N − N − N − blackout fraction, x P [I -I f ) / I > x ] USWIDADARR, k=4RML,
FIG. 2. Cumulative distributions of blackout sizes, measuredas a function of the fraction of the consumed power lost inthe blackout for different grid topologies described in the text,with the DC current approximation and the N − h k i = 5 . τ − . -4 -3 -2 -1 fraction of deleted lines, f d -5 -4 -3 -2 -1 P (f d > x ) USWI, N_1, all-at-onceUSWI, N_1, one-by-oneLBSPG, N-1, all-at-onceLBSPG, N-1, one-by-one (a) -5 -4 -3 -2 -1 fraction of deleted lines, f d -4 -3 -2 -1 P (f d > x ) α =0.1, al-at-once α =0.1, one-by-one α =0.8, all-at-once α =0.8, one-by-one α=1.0, all-at-once α=1.2, all-at-once α=1.4, all-at-once α=1.6, all-at-once α=1.8, all-at-once α=2.0, all-at-once (b) -5 -4 -3 -2 -1 fraction of deleted lines, f d -4 -3 -2 -1 P (f d > x ) α=0.0, all-at-once α =0.4, all-at-once α =0.6, all-at-once α =0.8, all-at-once α =0.9, all-at-once α=1.0, all-at-once α=1.1, all-at-once α=1.2, all-at-once α=2.0, all-at-once α=2.8, all-at-once α=0.6, one-by-one α=0.8, one-by-one α=0.9, one-by-one α=1.0, one-by-one α=1.1, one-by-oneslope -0.54 (c) FIG. 3. Cumulative distributions of blackout sizes, measuredas a function of the fraction of failed lines for the USWI andLBSPG models. (a) Both models for the N − α and bothupdate rules. One can see that for sufficiently large valuesof the tolerance, α > .
2, the distribution becomes approxi-mately power-law with large values of τ − ≥
1. The imple-mentation of the one-by-one update rule with the uniform tol-erance model shows that the power law distribution emergesfor smaller value of the tolerance, α = 0 . depends only on the tolerance level α [Fig. 3(b),(c)]. Asin [13], we start a cascade with an initial removal of asingle line. In agreement with the original results of [13],the bimodality disappears for α > α = 0 .
8, and a significant reduction of the sizes of the largest blackouts(Fig. 3(b)).All in all, our simulations suggest that neither the N − α .The N − α , but is not strictly nec-essary. Another important aspect is the network topol-ogy. For sparse homogeneous networks whose failuresare close to the percolation threshold such as our RMLmodel, the power-law distribution of blackouts emergesfor smaller levels of protection, while for the other modelssimilar failure rules still yield a bimodal distribution. IV. BETWEENNESS CENTRALITYOVERLOADS WITH THE N − CONDITION
As one can see, the models of the power grids we studyare relatively complex and have many different parame-ters, the role of which is difficult to clarify. Therefore,it is desirable to study a simpler model of overloadsin which the role of each parameter becomes clear. Itis also interesting to see if the emergence of the powerlaw distribution is a universal phenomenon for differentoverload models. A good candidate for such a modelis the betweenness centrality model, which displays aclear first order transition for the all-at-once removal ruleand uniform tolerance condition employed by Motter andLai [15]. However, we can expect that, similarly to thepower grid models, the implementation of the N − P ( k ) and compute the shortestpaths between each pair of nodes. The length of eachedge is assumed to be equal to 1 + ǫ , where ǫ is a nor-mally distributed random variable with a small standarddeviation σ ǫ ≪
1. This precaution is taken in order tomake sure that each pair of nodes ij has a unique shortestpath connecting them. The betweenness centrality b ( k )of each node k is computed as the number of all short-est paths ij passing through node k , such that k = i , k = j . The N − i is deleted, the rest of the nodes remain connectedand the betweenness centrality of each remaining node j , b i ( j ) does not exceed its maximal possible load b ∗ ( j ).Thus b ∗ ( j ) = max i = j b i ( j ) (2)and the original graph must be biconnected. Note thatin the original Motter and Lai model with a uniform tol-erance α > b ∗ ( j ) = (1 + α ) b ( j ) , (3)where b ( j ) is the initial betweenness of node j in a com-pletely intact network.To implement this model, we first construct a randomlyconnected graph of N nodes with a given degree distri-bution P ( k ) by the Molloy-Read algorithm [41] and findits largest biconnected component with N b nodes usingthe Hopcroft-Tarjan algorithm [42]. After that, we com-pute the betweenness centrality of each node and repeat N b simulations with each node removed in turn to find b ∗ ( i ) for each node. To find the distribution of cascadesof failures we remove a pair of random nodes and find thenodes i with overloads z = b ( i ) /b ∗ ( i ) >
1. In the firstmodel of all-at-once removal, at each stage of the cas-cade we remove all the nodes with overloads, recomputethe betweenness centrality of the remaining nodes andrepeat this process until no overloaded nodes are left. Inthe second model of one-by-one removal, at each stage ofthe cascades we remove only one of the overloaded nodes,the one with the largest overload z , we then recomputeall the centralities, and repeat this process until no over-loaded nodes are left. We study several cases of Erd¨osR´enyi (ER) networks with different average degrees h k i and different number of nodes in the biconnected com-ponent N b . The size of the blackout is computed as thefraction of removed nodes at the end of the cascade. Wealso study the distribution of the number of nodes dis-connected from the giant component.Again, the N − h k i = 1 . . < α < .
2, we seethe transition from a bimodal distribution of the cascadesto a power-law distribution with an exponential cut-off[Fig. 4(a)] as α increases. This is somewhat similar to theeffect observed for large α and small fraction p of nodesthat survive the initial attack in Ref. [30], where the bi-modal distribution of the cascades ceases to exist. Infact, small p effectively brings the network to the perco-lation transition equivalent to small h k i →
1. Replacingthe all-at-once update rule by the one-by-one rule shiftsthe transition from bimodal to a power-law distributionat a smaller α [Fig. 4(b)], but both behaviors can still beseen.Simultaneous application of the N −
10 100 1000 x=S b -3 -2 -1 P ( S b > x ) α=0.01α=0.1α=0.2α=0.3 x=S b -4 -3 -2 -1 P ( S b > x ) one-by-one α=0.1 one-by-one α=0.01 all-at-once α=0.1 all-at-once α=0.01 FIG. 4. (a) Cumulative distributions of blackout sizes, mea-sured as the fraction of failed nodes in the Motter andLai model with different values of tolerance α and all-at-once update rule. Black lines indicate graphs for values of α = 0 . , . , ... .
19. (b) Comparison of all-at-once and one-by-one update rules for selected values of the tolerance. Onecan see the emergence of the power law distribution in theone-by-one removal case for much smaller values of tolerance. values of h k i (Fig. 5). For small h k i , close to the perco-lation threshold and small system sizes, the bimodalityof the blackout distribution disappears and the distri-bution of blackouts becomes an approximate power lawwith an exponential cutoff. As the size of the networkincreases, the exponential cutoff disappears and the bi-modality emerges again, with a peak for large blackoutsemerging. The beginning of the cumulative distributionfor all values of h k i and all sizes of the system exhibits thesame slope of − /
2, which is the mean-field value for theSOC models [11]. For large h k i , we observe only a smallfraction of blackouts distributed as a power law, whilethe entire distribution remains bimodal; many networkscompletely disintegrate. However, the increase of the sys-tem size leads to the relative increase of the range of thepower law behavior. As the system size increases, we ob-serve similar changes in the distribution of the blackoutsfor all values of h k i studied (Fig. 5). For small systemsizes, we observe a smaller initial negative slope ( > − / − /
2, and the cutoffshifts to larger values. For even larger sizes, the slope ofthe initial part of the cumulative distribution increaseseven further ( < − /
2) but an inflection point appears inthe distribution, indicating the start of bimodality. Asthe system size continues to increase the negative initialslope continues to grow and the plateau region of thecumulative distribution separating large blackouts fromsmall becomes longer and longer.Interestingly, the behavior of the cumulative distribu-tion of the blackouts for fixed h k i and increasing systemsize is analogous to the behavior of the cumulative distri-bution of the sizes of avalanches in the failure-branchingprocess on a system of fixed size when the branching fac-tor h b i increases (Fig. 1). The cumulative distribution forthe avalanche distribution for h b i < τ −
1, and has a strong exponential cutoff.As h b i increases, the slope increases and the exponentialcutoff shifts to larger values. When h b i approaches thecritical point, h b i = 1, the distribution develops an in-flection point and a plateau that separates finite clustersand the giant component. For h b i >
1, the initial slopebecomes larger than τ . In contrast to this, the behaviorof the blackout size distribution, for fixed size, as h k i ap-proaches the critical value of the percolation transitionis analogous to increase of the system size in the distri-bution of avalanches in the failure-branching process. As h k i →
1, the range over which we see the power behaviorwith the mean field critical exponent τ = − . N − N − N real numbers,“fitnesses”, uniformly distributed between0 and 1 is created. At each time step the smallest ele-ment in the array and m randomly selected elements arereplaced by new values selected from the same uniformdistribution. As a result, after certain transition period,the majority of elements have fitness above the criticalvalue f c = 1 / ( m + 1). The number of active elements number of deleted nodes -3 -2 -1 P ( N d > x )
2. We see[Fig. 6(a)] that in our model the behavior of the numberof overloaded nodes is similar to the behavior of the ran-dom walk, which is consistent with the observed valueof τ ≈ .
5. However, the DFA analysis [43, 44] of thenumber of overloaded nodes during the cascade gives thevalue of the Hurst exponent ≈ .
1, which is much smallerthan the exponent for the uncorrelated random walk ofthe mean-field SOC model. The number of overloadednodes scales as ln N b , or a very small power law. In ourmodel, the fitness of a node is its relative overload z i , andthe level of z i is not arbitrary as in the SOC models, butis defined initially as z i = 1. When, after implementingthe N − N − k = 2 over nodes with k = 3.At the beginning of the cascade the removals take placein the biconnected component, reducing it by the lengthof a chain of nodes with k = 2 in between a pair of nodeswith k ≥
3. Each such removal only reduces the sizeof the giant component by one. But as more and morenodes are removed, further removals may happen in thesingly connected part of the network. Figure 6 (b) showsthe size of the biconnected component to which the re-moved node belongs as function of the cascade stage. If itbelongs to a singly connected part of the remaining net-work, we plot zero instead of the size of the biconnectedcomponent. When a node that is part of the singly con-nected part fails, the network separates into disconnectedparts and the size of the giant component drops signif-icantly. At the beginning, the smaller part is usually adangling end and the drop in size is relatively small, buteventually, two approximately equal parts get separated.At this point the betweenness of all nodes decreases dra-matically (approximately by factor of 4) and the cascadeends. Thus, although the cascade dynamics has somesimilarities with SOC, (at each time step the node withthe largest overload is removed and the overloads of therest of the nodes are changed), our model does not stayat the critical state but moves away from it.The comparison between a mean-field variant of SOCand our network overload dynamics is also instructive tohelp us understand when a bimodal distribution occursand when a power law distribution does. A power lawdistribution, such as in a branching process, occurs whenthere is a series of failures, each of which has a chance ofprecipitating b = 1 additional failures. Since there is thepossibility, at each stage of the cascade, of no new nodesfailing and the process stopping, the chance that exactlyone more node will fail decreases exponentially. At thecritical point, the cascade begins with b = 1, but as thenetwork disintegrates b changes. The all-at-once updaterule causes the network to disintegrate in fewer steps,since multiple nodes are removed at each step, and thus b increases during the cascade of failures without manychances for the cascade to end prematurely. Thus, aftera few steps b >
1, and the process will not stop itself withan intermediate level of destruction. However, there aretwo different conditions that will prevent b from increas-ing drastically as the network begins to disintegrate. Oneis a high level of overall protection, whether through theimplementation of the N − α . High overall protection means that minornetwork destruction will not cause widespread overloadand drive b away from one; only severe network destruc-tion can do that. Thus, the overload process has a chanceto stop spontaneously before it reaches that stage. Theother conditions is if the network is near the percolation nu m b e r o f ov e r l o a d s time step s i ze o f b i c onn ec t e d c o m p . size of biconnected comp.number of nodes in GC FIG. 6. (a) Dependence of the number of overloaded nodesas function of the cascade step for the Motter and Lai modelfor h k i = 1 . N = 500000, with the N − point, where the failure of a single node has a high chanceof separating a small, but sizable component from thenetwork. This will drastically decrease the overall loadand thus stabilize the network, arresting the cascade offailures [30]. In cnclusion, a one-by-one update rule and ahigh overall protection both contribute to the formationof a power law. The effect of the system size discussedabove can also be understood in light of the compari-son to a branching process. The larger the system, themore nodes (as an absolute number, not a fraction) failfollowing the initial attack. Thus, the chance that noneof the overloaded nodes will cause another node to over-load (i.e. a spontaneous end to the cascade) becomesincreasingly unlikely as the number of overloaded nodesincreases, even if b remains close to one for many stepsof the cascade. Thus, the chance that the cascade offailures will end with an intermediate amount of dam-age becomes more unlikely when the size of the originalsystem increases. V. DISCUSSION AND CONCLUSION
We have shown that the DC model of the power gridhas remarkable similarities to the Motter and Lai modelof betweenness centrality, suggesting that the exact phys-ical laws governing the flux do not play as big a role asthe network topology and the dynamics of the cascadepropagation. We have shown that both overload mod-els have similarities and differences with the topologicalmodels with dependency links [27]. In the topologicalmodels the outcome of the cascade depends only on thetopology of the network and the location of the initialdamage, while in the overload models the outcome sig-nificantly depends on the dynamics of the cascade, whichitself depends on the order of removal of overloaded el-ements and the relative speed of failure of overloadedelements and redistribution of flux after elements fail-ure. In general, if the failure is slow compared to theredistribution of flux (one-by-one removal rule) the cas-cades become smaller and the distribution of blackoutsbecomes unimodal, while for fast failures of overloadedelements (the all-at-once removal rule) the distributionof blackouts is bimodal. The shape of the distribution ofthe blackout sizes becomes approximately a power law asthe topology of the network approaches the percolationtransition for large level of protection. In summary, thegeneral picture of the behavior of the overload models isconsistent with the simple model of the failure branchingprocess. For large protection levels (small h b i ) the distri-bution of blackout sizes is a power law with exponentialcutoff, while for low protection levels (large h b i ) the dis-tribution of blackout sizes is bimodal, with a fraction oflarge avalanches which decreases as the protection levelincreases.The empirical observation of the power law distribu-tion of the blackout sizes in real power grids indicatesthat these grids are near the critical point of the cor-respondent failure branching process. However, this isunlikely to be caused by a SOC mechanism in which themodel is driven to the critical point by a certain dynamicrule. One hypothesis would be that it is the N − h b i < h b i to be as close as possible to the critical pointfrom below, so that limited but occasionally very largeblackouts are still possible. It is still desirable to developa model of long time evolution of power grids similarto [4], which reproduces the power-law distribution ofavalanches. The observation of power law distributionof the blackouts in the sparse networks close to the per-colation point also emphasizes the role of islanding forpreventing large blackouts. VI. ACKNOWLEDGEMENTS
Our research was supported by HDTRA1-14-1-0017,HDTRA1-19-1-0016, and HDTRA1-13-1-0021. SVB andGC acknowledge the partial support of this researchthrough the Dr. Bernard W. Gamson ComputationalScience Center at Yeshiva College. [1] “U.S-Canada Power System Outage Task Force report onthe August 14, 2003 blackout in the United States andCanada: Causes and recommendations.” 2004.[2] “Report of the enquiry committee on grid disturbance inthe Northern region on 30th July 2012 and in Northern,Eastern and North-Eastern region on 31st July 2012.”Aug. 2012.[3] “Report on the Grid Disturbance on 30 th July 2012 andGrid Disturbance on 31 st July 2012.” 2012.[4] H. Ren, I. Dobson, and B. A. Carreras, IEEE transac-tions on power systems 23, 1217 (2008).[5] B. A. Carreras, D. E. Newman, I. Dobson, IEEE Trans-actions on Power Systems 31, 4406 ( 2016 )[6] B. A. Carreras, V. E. Lynch, I. Dobson, and D. E.Newman. “Critical points and transitions in an electricpower transmission model for cascading failure black-outs.” Chaos , 985-994 (2002).[7] B. A. Carreras, V. E. Lynch, I. Dobson, and D. E. New-man. “Complex dynamics of blackouts in power trans-mission systems.” Chaos
3, 643-652 (2004).[8] I. Dobson and L. Lu. “Voltage Collapse Precipitated bythe Immediate Change in Stability When Generator Re-active Power Limits are Encountered.”
IEEE Trans. CAS vol. 39, No. 9, pp. 762-766. Sept. 1992.[9] P. Bak and K. Sneppen, Phys. Rev. Lett. , 4083 (1993).[10] M. Paczuski, S. Maslov, and P. Bak, Phys. Rev. E ,414 (1996). [11] H. Flyvbjerg, K. Sneppen,and P. Bak, Phys. Rev. Lett. , 4087 (1993).[12] R. Albert, I. Albert, and G. L. Nakarado. Phys. Rev. E, , 3694.[15] A. E. Motter and Y.C. Lai. Phys. Rev. E, , 065102(2002).[16] A. E. Motter. Phys. Rev. Lett., , 098701 (2004).[17] S. Soltan, D. Mazauric, and G. Zussman, Cascading fail-ures in power grids Analysis and algorithms, in Proc.ACM e-Energy14, 2014.[18] H. Cetinay, S. Soltan, F.A Kuipers, G. Zussman, and P.Van Mieghem, IEEE Transactions on Network Scienceand Engineering, (4), 301-312 (2017).[19] S. Soltan, A. Loh, and G. Zussman, IEEE Transactionson Control of Network Systems, (3), 1424-1433 (2017).[20] J.-W. Wang, L.-L. Rong, Safety Sci., , 5766 (2002). [23] J. P. Gleeson, and D. J. Cahalane, Proc. SPIE ,66010W (2007).[24] G.J. Baxter, S.N. Dorogovtsev, A.V. Goltsev, and J.F.F.Mendes, Phys. Rev. E , 011103 (2010).[25] G.J. Baxter, S.N. Dorogovtsev, A.V. Goltsev, and J.F.F.Mendes, Phys. Rev. E , 051134 (2011).[26] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, andS. Havlin, Nature , 1025 (2010).[27] R. Parshani, S. V. Buldyrev, and S. Havlin, Proc. Natl.Acad. Sci. U. S. A. 108, 1007 (2011).[28] M. A. Di Muro, S. V. Buldyrev, H. E. Stanley, and L. A.Braunstein, Phys. Rev. E , 042304 (2016).[29] M. A. Di Muro, L. D. Valdez, H. H. Aragao Rego, S. V.Buldyrev, H. E. Stanley and L. A. Braunstein, Sci. Rep. , 15059 (2017).[30] Y. Kornbluth, G. Barach, Y. Tuchman, B. Kadish, G.Cwilich, and S. V. Buldyrev, Phys.Rev. E 97, 052309(2018).[31] R. Parshani, S. V. Buldyrev, and S. Havlin, Phys. Rev.Lett. 105, 048701 (2010)[32] S.N. Dorogovtsev, A.V. Goltsev, J.F.F. Mendes, “k-core architecture and k-core percolation on complex net-works”, Physica D s by counting the number of clusters of differentsizes N ( s ) in a finite network, which is characterized bythe critical exponent τ = 5 /
2. An alternative definition,consistent with the notation in the SOC literature, is theprobability p ( s ) of a randomly selected site to belong toa cluster of size s . This probability distribution is charac-terized by a power law with an exponent τ ′ = τ − / , 2014.[39] S. Soltan, A. Loh, G. Zussman, “A learning-basedmethod for generating synthetic power grids,” IEEE Sys-tems Journal, , no. 1, pp. 625634, Mar. 2019.[40] S. Soltan, A. Loh, and G. Zussman, “Columbia Uni-versity Synthetic Power Grid with Geographical Co-ordinates,” https://doi.org/10.17041/drp/1471682, Jan.2018[41] Molloy, M. and Reed, B. The size of the giant componentof a random graph with a given degree sequence. Combin.Probab. Comput. 7, 295-305 (1998).[42] J. Hopcroft and R. Tarjan, Commun. ACM 16, 372(1973).[43] C. K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley and A. L. Goldberger, Mosaic Organizationof DNA Sequences, Phys. Rev. E , 1685-1689 (1994).[44] S. V. Buldyrev, A. L. Goldberger, S. Havlin, R. N. Man-tegna, M. E. Matsa, C.-K. Peng, M. Simons, and H.E. Stanley, Long-Range Correlation Properties of Cod-ing and Noncoding DNA Sequences: GenBank Analysis,Phys. Rev. E51