Inhibiting failure spreading in complex networks
IInhibiting failure spreading in complex networks
Franz Kaiser,
1, 2
Vito Latora,
3, 4, 5 and Dirk Witthaut
1, 21
Forschungszentrum J¨ulich, Institute for Energy andClimate Research (IEK-STE), 52428 J¨ulich, Germany Institute for Theoretical Physics, University of Cologne, K¨oln, 50937, Germany School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, UK Dipartimento di Fisica ed Astronomia,Universit`a di Catania and INFN, 95123 Catania, Italy The Alan Turing Institute, The British Library, London NW1 2DB, UK (Dated: September 8, 2020)
Abstract
In our daily lives, we rely on the proper functioning of supply networks, from power grids towater transmission systems. A single failure in these critical infrastructures can lead to a completecollapse through a cascading failure mechanism. Counteracting strategies are thus heavily soughtafter. In this article, we introduce a general framework to analyse the spreading of failures incomplex networks and demonstrate that both weak and strong connections can be used to containdamages. We rigorously prove the existence of certain subgraphs, called network isolators, thatcan completely inhibit any failure spreading, and we show how to create such isolators in syntheticand real-world networks. The addition of selected links can thus prevent large scale outages asdemonstrated for power transmission grids. a r X i v : . [ phy s i c s . s o c - ph ] S e p omplex networked systems are subject to external perturbations, damages or attackswith potentially catastrophic consequences [1, 2]. The loss of even a single edge can causea blackout in a power grid [3, 4], the dieback of a biological network [5], or the collapse ofan entire ecological network [6]. It is thus essential to understand how the structure of anetwork determines its response to perturbations and its global resilience [7–11]. Here, wepropose a general framework to model the redistribution of flows in a complex network thatfollows a small and local failure, and we suggest novel and more efficient strategies to improvenetwork resilience. Our findings reveal that propagation of damages can counterintuitivelybe better limited by adding selected links instead of removing links and can turn very usefulto construct more robust networks or to improve existing ones.The division of a network into weakly coupled parts provides the most intuitive methodto inhibit the spreading of failures, thus improving system resilience [12–15]. An example isshown in Fig. 1(a) for an elementary supply network with two weakly connected modules.The response to an edge failure is strong locally, but it is reduced in the other module of thenetwork which has only few links to the part where the failure happened. A similar effectis observed in a real-world case, the Scandinavian power grid in Fig. 1(d). The study ofcommunity structures in both natural and man-made systems is an integral part of networkscience: a variety of methods has been developed to define and identify the weakly connectedmodules of a network [16–18], and the important role of community structures in networkdynamics is today well recognised.Limiting connectivity for the sake of additional security is, however, not always desirable.For instance, microgrid concepts and intentional islanding are heavily discussed in energysystems research [21, 22], but the overall demand for electric power transmission actuallyincreases [23, 24]. Other methods to contain perturbations or damages in complex networksare thus needed. Indeed, an exceptionally strong inter-connectivity between two modules canalso suppress failure spreading as shown in Fig. 1(b,e). Notably, a strong inter-connectivitycan be realised in different ways. In the random network example in Fig. 1(b), a high numberof links connects a subset of nodes of the two modules. In real vascular networks of leafs,the suppression of failure spreading occurs naturally because the central vein between theleft and right parts has an exceptionally large capacity (Fig. 1(e), cf. also [25]).Remarkably, failure spreading can be completely stopped by certain subgraphs which werefer to as network isolators in the following, an example being shown in Fig. 1(c). The2 eak connectivity | ∆ F | (a) strong connectivity | ∆ F | (b) isolator | ∆ F | (c) | ∆ F | (d) | ∆ F | (e) | ∆ F | (f) − − − − − − − − − − − − − − − − − − Figure 1.
Different network structures inhibit the spreading of failures in complexnetworks.
We simulate the impact of a single failing link (red) for different network structures;resulting flow changes are colour coded. (a,b) Both a weak and a strong inter-connectivity cansuppress the spreading of failures between two modules of a complex network. (c) Failure spreadingis prevented completely by a network isolator (blue shading); flow changes on the grey links areexactly zero. (d) The Scandinavian powergrid consists of three weakly connected modules whichsuppresses failure spreading between the modules [19]. (e) The vascular network of a
Burserahollickii leaf contains a strong central vein [20], which suppresses failure spreading between thetwo sides of the leaf. (f) Same as in (d) but with the addition of two links (blue shading) to createa network isolator. failure of an edge in the right part of the network does not affect the flows in the left part atall. Real world networks can be made perfectly resistant to failure propagation by the ad-hoc addition of few links to create network isolators, as demonstrated for the Scandinavianpower grid in Fig. 1(f) consisting of three weakly coupled modules. The suppression offailure spreading is readily generalised to networks with more than two modules [26].Our results are based on a general framework that allows a theoretical analysis of theinterplay of network connectivity and robustness in the context of supply or transportationnetworks. Many such systems can in fact be modelled by linear flow networks where the3ow over an edge ( i, j ) depends linearly on the gradient of a potential function across theedge, F i → j = K ij · ( ϑ i − ϑ j ). In particular, this description applies to power transmissiongrids [27], where F is the real power flow, ϑ i denotes the nodal voltage phase angle and K ij is given by the line susceptance, and to vascular networks [28, 29], where F is the flow ofwater or nutrients, ϑ i is the local pressure and K ij the edge’s capacity (see SupplementalMaterial [26]). The generalisation to nonlinear systems will be discussed below. Further-more, equivalent problems arise in the linearisation of general diffusively coupled networksof dynamical systems around an equilibrium or limit cycle [30].The impact of a damage in linear flow networks can be calculated analytically. Assumethat an edge (cid:96) = ( r, s ) fails, and summarise the response at all nodes i = 1 , . . . , N by thevector of changes ∆ (cid:126)ϑ = (∆ ϑ , . . . , ∆ ϑ N ) (cid:62) . One then finds [26] L ∆ (cid:126)ϑ = q (cid:96) (cid:126)ν (cid:96) , (1)where L is the Laplacian matrix of the weighted network, (cid:126)ν (cid:96) is a vector with +1 at position r and − s , and q (cid:96) = 1 − K rs (cid:126)ν (cid:62) (cid:96) L − (cid:126)ν (cid:96) is a source strength [31]. We thus findthat the response of a network to failures is essentially determined by the Laplacian L . Thismatrix incorporates the network structure and is defined as follows; L ij = − K ij if i is connected to j, (cid:80) ( i,k ) ∈ E ( G ) K ik if i = j, . (2)Notably, the Laplacian matrix is also useful to infer the large scale connectivity and thecommunity structure of a given network [32].To quantify the effect of connectivity on failure spreading, we have studied the impact ofdifferent failures in a variety of synthetic networks as well as in several real-world networks.For a given initial failure of an edge (cid:96) , we compute the flow changes ∆ F i → j = K ij · (∆ ϑ i − ∆ ϑ j )for all edges ( i, j ) in a given subgraph G (cid:48) of the network. Furthermore, we must take intoaccount that the impact of a failure generally decreases with distance [31, 33, 34]. As anoverall measure of the impact of a failure we thus consider the expression (cid:104)| ∆ F i → j |(cid:105) ( i,j ) ∈ G (cid:48) d ,which gives the magnitude of flow changes averaged over all edges ( i, j ) ∈ G (cid:48) at a givendistance d to the edge (cid:96) . The prime question is now whether the impact differs substantially4 solator Figure 2.
Effectiveness and robustness of shielding network structures. (a,b) Adjacencymatrices for the graphs shown in Fig. 1(a,b). Two random graphs G (30 , .
4) are inter-connected viaa fraction c of their nodes chosen at random, and links are added with probability µ , interpolatingbetween weak (a) or strong (b) inter-connectivity. (c) Adjacency matrix for the 6-regular graphshown in Fig. 1(c) containing a network isolator. (d) The average ratio of flow changes R ( (cid:96) ) in thetwo components (Eq. 3) is strongly suppressed for both high and low inter-connectivity µ . The blueline represents the median value over all distances and the shaded region indicates the 0 .
25- and0 . R vanishes for a perfect network isolator described bythe condition ξ ( A ) = 0 and increases algebraically with coherence parameter ξ when perturbed. between the communities or moduli of a network. We thus plot the ratio R ( (cid:96), d ) = (cid:104)| ∆ F i → j |(cid:105) ( i,j ) ∈ O d (cid:104)| ∆ F i → j |(cid:105) ( i,j ) ∈ S d . (3)between the remote part of the network G (cid:48) = O and the part G (cid:48) = S containing the failingedge (cid:96) . If this ratio approaches or reaches zero, this is indicative of a very strong suppressionof failure spreading into the other part of the network.To study how the impact of failure spreading depends on the network structure, weconsidered synthetic graphs obtained by connecting two Erd˝os-R´enyi (ER) random graphs toeach other at preselected, randomly chosen vertices with a tunable probability µ ∈ [0 ,
1] [26,35]. The resulting graphs have a connectivity structure ranging from two weakly connectedcommunities for low values of µ shown in Figure 1(a) to strongly connected parts shownin panel (b). In the limit µ = 1, the two modules are connected via a complete bipartite5 solator Grid modi fi cation Figure 3.
Network isolators can contain cascading failures in power grids. (a) Lineloading (colour code) on the Scandinavian grid in units relative to maximal loading before theinitial failure of a single line (coloured red). (b) The initial failure results in a cascade of overloads(red coloured lines) until the grid disconnects into several parts. (c) Magnification of the gridstructure in Eastern Norway. A small modification of the grid enables a network isolator. (d)Introducing the network isolator completely suppresses the spreading of failures from WesternNorway to the rest of the grid thus inhibiting the cascade observed in (b). graph as shown in Fig. 1(c). This is a possible realization of a network isolator , since itcompletely suppresses flow changes as we will explain in the following. The correspondingadjacency matrices clearly indicate the connectivity structure, revealing the strong or weakcoupling between the two modules of the networks (Fig. 2(a,b,c)). Remarkably, evaluatingthe quantity R ( (cid:96) ), obtained by averaging the ratio over flow changes R ( (cid:96), d ) over all distances d for a specific trigger link (cid:96) , for a varying connectivity structure tuned by µ , we find thatthe spreading of failures is largely suppressed for both weak and strong connectivity betweenthe two modules as shown in Fig. 2(d). Distance plays a minor role for the ratio of flowchanges R ( d ) as illustrated in a separate Figure in the Supplemental Material [26].Network symmetries are known to play an important role for the dynamics and syn-6hronisability of a network [36–38]. Network isolators as a specific connectivity structurecompletely inhibit the spreading of failures from one network module to another. They man-ifest also as particular, symmetric patterns in the region of the adjacency matrix describingthe connectivity between the two parts of the network as we have seen in Fig. 2(c). They arecharacterised by the following theorem, which we prove in the Supplemental Material [26] Theorem 1.
Consider a linear flow network composed of two modules 1,2 and let A denote the weighted adjacency matrix of the mutual connections. An edge failure in onemodule does not affect the flows in the other module if rank( A ) = 1 . For unweightednetworks this criterion is fulfilled if A describes a complete bipartite graph. Note that, while network isolators prevent failure spreading, we found that they do notinfluence network controllability [26].Since most real world examples of networks do not contain perfect network isolators,we have studied the robustness of a network isolator against modifications of the topology.Starting from a unit rank matrix, we perturb the adjacency matrix A iteratively. In eachstep we choose one of the matrix columns (cid:126)a i , i = 1 , . . . , m at random and perturb it accordingto (cid:126)a (cid:48) i = (cid:126)a i + (cid:126)e (cid:107) (cid:126)a i (cid:107) . The elements of the perturbation vector (cid:126)e are chosen uniformly at randomfrom the interval [ − β, β ], where β is a small parameter, here β = 0 .
05. The deviation of theperturbed matrix A from a unit rank matrix is quantified using its coherence statistics [39], ξ ( A ) = 1 − min i,j (cid:104) (cid:126)a i , (cid:126)a j (cid:105)(cid:107) (cid:126)a i (cid:107)(cid:107) (cid:126)a j (cid:107) . (4)The performance of the isolator is then measured by calculating the ratio of flow changes R , which is obtained from R ( (cid:96), d ) by averaging over all possible trigger links and distances.A perfect isolator is characterised by ξ ( A ) = 0 and enables a complete containment offailure spreading such that R = 0. For perturbed isolators, we find that R increases approx-imately algebraically with ξ ( A ), see Fig. 2(e). Hence, the isolation effect persists for smallperturbations, albeit with reduced efficiency.Perfect network isolators can be easily constructed to improve the resilience of complexnetworked systems. As a practical example we show an application to electric power grids,where large scale blackouts are typically triggered by the outage of a single transmissionelement which leads to a cascade of failures [3, 40]. We demonstrate the impact of networkisolators against cascading failures in the case of the Scandinavian grid. In the original grid7 ertex cut = 2 edge cut = 2 edge cut = 3 edge cut = 3initial failure Figure 4.
Different ways of constructing isolators and their effects in full non-linear ACload flow. (a) to (d) Alternative methods of creating an isolator in a given network. We show thenetwork structure before (top left) and after (top right) the addition of a network isolator, as wellas corresponding adjacency matrices (bottom) with the different shades of blue representing thecapacity K ij of the respective edge. A lower prior connectivity simplifies the creation of isolators asmeasured by the vertex cut (a) or edge cut (b,c,d) which is visible in the adjacency matrix (entriescolored red). The creation of network isolators results in characteristic patterns in the adjacencymatrix in terms of the capacities of the isolator edges (shades of blue). (e) An initially failing linkwith unit flow (red) in the British grid results in changes of real power flow (colour code) throughoutthe whole network, as obtained by computing a non-linear full AC power flow [19, 26]. (f,g) Afterintroducing a network isolator based on the strategy presented in panel (a), failure spreading isperfectly inhibited in the linear power flow approximation, and still significantly reduced in thenon-linear full AC load flow. (h) Median absolute flow changes are significantly lower in the uppermodule of the grid shielded by the isolator (light blue, straight line) compared to the situationbefore introducing the isolator (light blue, dotted line), whereas flow changes in the lower moduleare almost the same (dark blue lines). layout, failure spreading between different areas of the grid is reduced due to the presenceof only few connections between different areas of the network – but it is possible. A failurein one area can thus spread to other areas and cause a global cascade, as demonstrated in8igure 3(a,b) for a cascade emerging in Western Norway. This spreading may in principlebe prevented by decoupling different areas of the grid, but this is highly undesired. Infact, future energy systems will require more connectivity, not less, to transmit renewableelectric power [23, 24]. In contrast, building a network isolator can completely inhibit failurespreading at increased connectivity. A perfect isolator can be realised with moderate effortby reconstructing two substations in Norway, such that they effectively form two nodeseach. The new nodes must be linked by internal connections and one additional two-circuitoverhead line, whose parameters are optimised such that the condition rank( A ) = 1 issatisfied (Fig. 3(c)). A simulation for such an optimised grid layout shows that the spreadingof the cascade is completely suppressed (Fig. 3(d)). The grid remains connected and loadshedding is no longer necessary as a containment strategy [2, 3].Network isolators are not limited to the particular situation shown in Figure 1, or to linearflow models only. In Figure 4(a,d) we identify several subgraphs that allow to easily introducenetwork isolators into existing topologies. For subgraphs with a prior low connectivity, asmeasured by a small vertex cut (Fig. 4(a)) or a small edge cut (Fig. 4(b,c,d)), networkisolators may be introduced with small grid modifications - by adding (a,b,d) or removingand adding (c) selected links. The concept of network isolators has been established forlinear flow networks, but can be generalised in two ways. (1) We can rigorously proof theexistence of isolators for an important class of nonlinear network dynamical systems [26]. (2)For many nonlinear systems of practical importance, the impact of failures or perturbationsis well described by a linearisation around an equilibrium or limit cycle [30]. We demonstratethat the strong isolation of line outages in power grids persists beyond the linearised flowequations via direct numerical simulations in Fig. 4(e to h).In conclusion, connectivity determines the resilience of complex networks in manifoldways. As expected, a division of a network into weakly coupled modules suppresses thespreading of failures from one module to the others. Remarkably, we have found that a stronginter-connection can equally well suppresses the spreading in both flow networks and innetworks of nonlinear dynamical systems. We have demonstrated that an even stronger effectcan be created by certain subgraphs called isolators, which inhibit the spreading completely.Isolators can be applied to mitigate cascading failures, for instance in electric power grids,while enabling an arbitrary degree of connectivity between the grid parts. These resultswiden our perspective on the large scale organisation of complex networks in general, showing9hat very diverse structural patterns can exist that discriminate functional units and improvenetwork resilience.We thank Tom Brown and Jonas H¨orsch for help with the processing of power grid data,Eleni Katifori and Henrik Ronellenfitsch for providing the leaf data and Raissa D’Souzaand J¨urgen Kurths for valuable discussions. We gratefully acknowledge support from theFederal Ministry of Education and Research (BMBF grant no. 03EK3055B D.W.), theHelmholtz Association (via the joint initiative ”Energy System 2050 – a Contribution of theResearch Field Energy” and grant no. VH-NG-1025 to D.W.). V.L.’s work was funded bythe Leverhulme Trust Research Fellowship CREATE. [1] R. Albert, H. Jeong, and A.-L. Barab´asi, Error and attack tolerance of complex networks,Nature , 378 (2000).[2] Y. Yang, T. Nishikawa, and A. E. Motter, Small vulnerable sets determine large networkcascades in power grids, Science , eaan3184 (2017).[3] P. Pourbeik, P. Kundur, and C. Taylor, The anatomy of a power grid blackout - root causesand dynamics of recent major blackouts, IEEE Power and Energy Magazine , 22 (2006).[4] D. Witthaut, M. Rohden, X. Zhang, S. Hallerberg, and M. Timme, Critical links and nonlocalrerouting in complex supply networks, Physical Review Letters , 138701 (2016).[5] L. Sack, E. M. Dietrich, C. M. Streeter, D. Sanchez-Gomez, and N. M. Holbrook, Leaf palmatevenation and vascular redundancy confer tolerance of hydraulic disruption, Proceedings of theNational Academy of Sciences , 1567 (2008).[6] J. A. Dunne and R. J. Williams, Cascading extinctions and community collapse in model foodwebs, Philosophical Transactions of the Royal Society B: Biological Sciences , 1711 (2009).[7] S. H. Strogatz, Exploring complex networks, Nature , 268 (2001).[8] S. Sahasrabudhe and A. E. Motter, Rescuing ecosystems from extinction cascades throughcompensatory perturbations, Nature Communications , 170 (2011).[9] A. E. Motter, Cascade Control and Defense in Complex Networks, Physical Review Letters , 098701 (2004).[10] R. M. D’Souza, Curtailing cascading failures, Science , 860 (2017).[11] Y. Lin, K. Burghardt, M. Rohden, P.-A. No¨el, and R. M. D’Souza, Self-organization of dragon ing failures, Physical Review E , 022127 (2018).[12] D. B. Stouffer and J. Bascompte, Compartmentalization increases food-web persistence, Pro-ceedings of the National Academy of Sciences , 3648 (2011).[13] L. J. Gilarranz, B. Rayfield, G. Li˜n´an-Cembrano, J. Bascompte, and A. Gonzalez, Effects ofnetwork modularity on the spread of perturbation impact in experimental metapopulations,Science , 199 (2017).[14] C. Hens, U. Harush, S. Haber, R. Cohen, and B. Barzel, Spatiotemporal signal propagationin complex networks, Nature Physics , 403 (2019).[15] C. D. Brummitt, R. M. D’Souza, and E. A. Leicht, Suppressing cascades of load in interde-pendent networks, Proceedings of the National Academy of Sciences , E680 (2012).[16] F. Radicchi, C. Castellano, F. Cecconi, V. Loreto, and D. Parisi, Defining and identifyingcommunities in networks, Proceedings of the National Academy of Sciences , 2658 (2004).[17] M. Girvan and M. E. J. Newman, Community structure in social and biological networks,Proceedings of the National Academy of Sciences , 7821 (2002).[18] M. E. J. Newman, Communities, modules and large-scale structure in networks, Nature Physics , 25 (2012).[19] J. H¨orsch, F. Hofmann, D. Schlachtberger, and T. Brown, PyPSA-eur: An open optimisationmodel of the european transmission system, Energy Strategy Reviews , 207 (2018).[20] H. Ronellenfitsch, J. Lasser, D. C. Daly, and E. Katifori, Topological Phenotypes Constitutea New Dimension in the Phenotypic Space of Leaf Venation Networks, PLOS ComputationalBiology , e1004680 (2015).[21] R. H. Lasseter, Microgrids, in , Vol. 1 (IEEE, 2002) pp. 305–308.[22] M. Mureddu, G. Caldarelli, A. Damiano, A. Scala, and H. Meyer-Ortmanns, Islanding thepower grid on the transmission level: less connections for more security, Scientific Reports ,34797 (2016).[23] D. P. Schlachtberger, T. Brown, S. Schramm, and M. Greiner, The benefits of cooperation ina highly renewable european electricity network, Energy , 469 (2017).[24] T. Tr¨ondle, J. Lilliestam, S. Marelli, and S. Pfenninger, Trade-offs between geographic scale,cost, and infrastructure requirements for fully renewable electricity in europe, Joule (2020).[25] T. Gavrilchenko and E. Katifori, Resilience in hierarchical fluid flow networks, Physical Review , 012321 (2019).[26] See Supplemental Material for a proof of the main theorem on network isolators, extendedFigures on their applicability, a discussion of their relationship to network controllability anddetailed methods, which include Refs. [41–57].[27] K. Purchala, L. Meeus, D. V. Dommelen, and R. Belmans, Usefulness of dc power flow foractive power flow analysis, in IEEE Power Engineering Society General Meeting (2005) pp.454–459 Vol. 1.[28] E. Katifori, G. J. Sz¨oll˝osi, and M. O. Magnasco, Damage and fluctuations induce loops inoptimal transport networks, Physical Review Letters , 048704 (2010).[29] D. A. Coomes, S. Heathcote, E. R. Godfrey, J. J. Shepherd, and L. Sack, Scaling of xylemvessels and veins within the leaves of oak species, Biology Letters , 302 (2008).[30] D. Manik, M. Rohden, H. Ronellenfitsch, X. Zhang, S. Hallerberg, D. Witthaut, andM. Timme, Network susceptibilities: Theory and applications, Physical Review E , 012319(2017).[31] J. Strake, F. Kaiser, F. Basiri, H. Ronellenfitsch, and D. Witthaut, Non-local impact of linkfailures in linear flow networks, New Journal of Physics , 053009 (2019).[32] S. Fortunato, Community detection in graphs, Physics Reports , 75 (2010).[33] F. Kaiser, J. Strake, and D. Witthaut, Collective effects of link failures in linear flow networks,New Journal of Physics , 013053 (2020).[34] S. Kettemann, Delocalization of disturbances and the stability of ac electricity grids, PhysicalReview E , 062311 (2015).[35] P. Erd˝os, On the evolution of random graphs, Publications of the Mathematical Institute ofthe Hungarian Academy of Sciences , 17 (1960).[36] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, Cluster syn-chronization and isolated desynchronization in complex networks with symmetries, NatureCommunications , 1 (2014).[37] F. Sorrentino, L. M. Pecora, A. M. Hagerstrom, T. E. Murphy, and R. Roy, Complete charac-terization of the stability of cluster synchronization in complex dynamical networks, ScienceAdvances , e1501737 (2016).[38] V. Nicosia, M. Valencia, M. Chavez, A. D´ıaz-Guilera, and V. Latora, Remote SynchronizationReveals Network Symmetries and Functional Modules, Physical Review Letters , 174102 , 948 (2010).[40] B. Sch¨afer, D. Witthaut, M. Timme, and V. Latora, Dynamically induced cascading failuresin power grids, Nature Communications , 1975 (2018).[41] M. E. J. Newman, Networks: An introduction (Oxford University Press, 2010).[42] B. Bollob´as,
Modern graph theory , Graduate texts in mathematics No. 184 (Springer, 1998).[43] F. D¨orfler, J. W. Simpson-Porco, and F. Bullo, Electrical networks and algebraic graph theory:Models, properties, and applications, Proceedings of the IEEE , 977 (2018).[44] A. J. Wood, B. F. Wollenberg, and G. B. Shebl´e,
Power Generation, Operation and Control (John Wiley & Sons, New York, 2014).[45] A. E. Motter and Y.-C. Lai, Cascade-based attacks on complex networks, Physical Review E , 065102 (2002).[46] P. Crucitti, V. Latora, and M. Marchiori, Model for cascading failures in complex networks,Physical Review E , 045104 (2004).[47] Z. Yuan, C. Zhao, Z. Di, W.-X. Wang, and Y.-C. Lai, Exact controllability of complex net-works, Nature Communications , 2447 (2013).[48] P. Van Mieghem, K. Devriendt, and H. Cetinay, Pseudoinverse of the laplacian and bestspreader node in a network, Physical Review E , 032311 (2017).[49] J. Gao, Y.-Y. Liu, R. M. D’Souza, and A.-L. Barab´asi, Target control of complex networks,Nature Communications , 5415 (2014).[50] Y.-Y. Liu, J.-J. Slotine, and A.-L. Barab´asi, Controllability of complex networks, Nature ,167 (2011).[51] N. Biggs, Algebraic potential theory on graphs, Bulletin of the London Mathematical Society , 641 (1997).[52] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Self-organized synchronization in decen-tralized power grids, Phys. Rev. Lett. , 064101 (2012).[53] T. Nishikawa and A. E. Motter, Comparative analysis of existing models for power-grid syn-chronization, New Journal of Physics , 015012 (2015).[54] D. Manik, D. Witthaut, B. Sch¨afer, M. Matthiae, A. Sorge, M. Rohden, E. Katifori, andM. Timme, Supply networks: Instabilities without overload, The European Physical Journal pecial Topics , 2527 (2014).[55] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, The kuramoto model in complexnetworks, Physics Reports , 1 (2016).[56] R. A. Wurbs and W. P. James, Water Resources Engineering (Pearson, 2001).[57] J. Reichold, M. Stampanoni, A. L. Keller, A. Buck, P. Jenny, and B. Weber, Vascular graphmodel to simulate the cerebral blood flow in realistic vascular networks, Journal of CerebralBlood Flow & Metabolism , 1429 (2009). upplemental Material FLOW NETWORKS
In this section, we briefly review the theory and applications of linear flow networks.
Mathematical description
In this work, the main model of interest is a linear flow network model which we introducemore formally in this section. Consider a connected graph G = ( E, V ) consisting of N = | V | nodes and L = | E | edges. Assign to each node in the network a potential ϑ n ∈ R , n ∈ V ( G )and to each edge a capacity K ij ∈ R + , (cid:96) = ( i, j ) ∈ E ( G ). Now we assign a flow F i → j ∈ R to each link (cid:96) = ( i, j ) ∈ E ( G ) in the network that is assumed to be linear in the potentialdrop F i → j = K ij · ( ϑ i − ϑ j ) = − F j → i . (5)Suppose that there are sources and sinks attached to the nodes of the networks P i ∈ R , i ∈ V ( G ). In this case, the in- and outflows at each node have to balance with the sources andsinks P i = N (cid:88) k =1 F i → k . (6)This equation is known as continuity equation or Kirchhoff ’s current law . If the sourcesand sinks P i are given, Eqs. (6) and (5) completely determine the potentials in the network(up to a constant shift to all potentials). In a power grid, the sources and sinks are thepower injections or withdrawals as a result of power production or consumption. Whenlooking at the stable, operational fixed point of a power grid they are balanced such that (cid:80) i P i = 0 – we therefore assume this to hold in the following sections. The theory of linearflow networks applies resistor networks, as well as AC power grids in the DC approximation,hydraulic networks and networks of limit cycle oscillators, which will be discussed in detailin section .Now we introduce a compact, vectorial notation which facilitates the analysis of pertur-bations or damages to the network. Note that the flow is a signed quantity that dependson the orientation of the edges that we arbitrarily fix for this purpose and say that the flow15s directed from node i to node j in this case. We can write the flows in vectorial notation (cid:126)F = ( F , ..., F L ) (cid:62) ∈ R L as follows; (cid:126)F = KI (cid:62) (cid:126)ϑ. (7)Here, K = diag( K , ..., K L ) ∈ R L × L is the graph’s weight matrix that collects the edgeweights and I (cid:62) is the transpose of the the graph’s edge-node incidence matrix I ∈ Z N × L that determines the orientation of the graph’s edges by the following relationship I j(cid:96) = +1 if edge (cid:96) ˆ= ( j, k ) starts at node j, − (cid:96) ˆ= ( j, k ) ends at node j , . (8)Furthermore, (cid:126)ϑ = ( ϑ , ..., ϑ N ) (cid:62) ∈ R N is a vector of potentials or voltage phase angles. Wecan also define a vector of power injections (cid:126)P = ( P , ..., P N ) (cid:62) ∈ R N such that the continuityequation reads as (cid:126)P = I (cid:126)F . (9)In this expression, the correspondence between the power balance and Kirchhoff ’s currentlaw becomes manifest: it states that the in- and outflows at each node have to balancethe injections and withdrawals of power. Combining Equations (7) and (9), we may find arelationship between angles (cid:126)ϑ and power injections (cid:126)P , thus defining the graph’s weighted
Laplacian matrix L = IKI (cid:62) ∈ R N × N , by (cid:126)P = IKI (cid:62) (cid:126)ϑ = L (cid:126)ϑ. (10)The weighted Laplacian matrix used here has the following entries L ij = − K (cid:96) if i is connected to j via (cid:96) = ( i, j ) , (cid:80) (cid:96) =( i,k ) ∈ E ( G ) K (cid:96) if i = j, . (11)The Laplacian matrix plays an important role in graph theory [41]. If the underlying graph G is connected, it has one zero eigenvalue λ = 0 with corresponding eigenvector (cid:126)v = (cid:126) / √ N .Therefore, the matrix is not invertible. In many cases, it would nevertheless be desirableto invert the matrix, e.g. in order to find the phase variables given the power injections inEquation (10). This problem is typically overcome by making use of the matrix’s Moore-Penrose pseudoinverse L † . It may be used to invert Equation (10) in the same way as for the16rdinary matrix inverse in the case of balanced power injections [48]. The Moore-Penrosepseudoinverse of the graph Laplacian L allows for the following representation: using L ’seigenvalues sorted by magnitude λ = 0 , λ ≤ ... ≤ λ N with corresponding eigenvectors (cid:126)v = (cid:126) / √ N , (cid:126)v , ..., (cid:126)v N , we can express its pseudoinverse L † as [43] L † = ( (cid:126)v , (cid:126)v , ..., (cid:126)v N ) ... λ − ... ... ... ... ...... ... ... λ − N ( (cid:126)v , (cid:126)v , ..., (cid:126)v N ) (cid:62) . The second eigenvalue λ is usually referred to as Fiedler eigenvalue or algebraic connectivity and is an indicator of the graph’s overall connectivity. If we assume the overall graph to beconnected, this eigenvalue is strictly greater than zero λ >
0. Importantly, a large differencebetween second and third eigenvalue λ − λ implies a strong modularity in the graph andthus indicates the existence of a community structure [16–18].Before we proceed, let us briefly fix the notation for the following sections: we will referto an edge (cid:96) = ( (cid:96) , (cid:96) ) ∈ E ( G ) and its index (cid:96) in the ordered set of all edges interchangeablyor refer to it by its terminal nodes (cid:96) and (cid:96) . If we assume the edge space to be spannedby vectors in the two element field GF (2), we may express the edge by a unit vector (cid:126)l (cid:96) =(0 , ..., (cid:124)(cid:123)(cid:122)(cid:125) l , .. (cid:62) ∈ GF (2) L which we refer to as the edge’s indicator vector. The edge-nodeincidence matrix I then maps this unit vector to the corresponding unit vectors in the fieldof vertices GF (2) N . We thus get the following result for the edge expressed in terms of itsstarting vertex (cid:96) and terminal vertex (cid:96) : (cid:126)ν (cid:96) = I · (cid:126)l (cid:96) = (cid:126)e (cid:96) − (cid:126)e (cid:96) = ... ... − ... } (cid:96) } (cid:96) , (cid:126)e (cid:96) and (cid:126)e (cid:96) are basis vectors in GF (2) N (cid:126)e (cid:96) = ... ...... } (cid:96) , (cid:126)e (cid:96) = ...... ... } (cid:96) . This formulation allows us to easily switch between the edges expressed in edge space andthe nodes corresponding to its terminal ends.
Applicability of linear flow models
The theoretical framework in the last section has many different applications. We willdemonstrate its applicability to the following systems in this section:1. Power grids, [31, 44]2. Resistor networks [42],3. Hydraulic networks [5, 28],4. Limit cycle oscillators [30].
Application to power grids
The power flow equations describing the steady state of a power system at an arbitrarynode i are given by [44] P i = N (cid:88) k =1 | V i || V k | ( G ik cos( ϑ i − ϑ k ) + B ik sin( ϑ i − ϑ k )) Q i = N (cid:88) k =1 | V i || V k | ( G ik sin( ϑ i − ϑ k ) − B ik cos( ϑ i − ϑ k )) . (12)Here, P i and Q i are the real and reactive power generated or consumed at node or bus i , ϑ i is the voltage angle at the same bus and | V i | is the voltage magnitude. The matrices18 ∈ R N × N and B ∈ R N × N with elements G ij and B ij , respectively, are the real part andthe complex part of the complex nodal admittance matrix Y = G + i B ∈ C N × N . Notethat the matrices B and G are not actually matrices of susceptances and conductances,respectively. Instead, their entries read as follows B ij = − b ij if ( i, j ) ∈ E ( G ) , i (cid:54) = j,b shunt i + (cid:80) ( i,k ) ∈ E ( G ) b ik if i = j, , where b shunt i denotes the shunt susceptance of node i and b ij is the susceptance of the circuitconnecting node i to node j . G has an analogous structure with elements G ij = − g ij if ( i, j ) ∈ E ( G ) , i (cid:54) = j,g shunt i + (cid:80) ( i,k ) ∈ E ( G ) g ik if i = j, , where g ij are the conductances of the circuit between nodes i and node j . The matrices B and G thus have the structure of a Laplacian matrix except for the diagonal entrieswhich contain additional terms given by the shunt susceptances and conductances. Theoff-diagonal elements of the nodal admittance matrix thus read as Y jk = − y jk , ∀ j (cid:54) = k ; y jk = g jk + ib jk = 1 r jk + ix jk , with the circuit’s reactance x jk and resistance r jk . Note that line susceptances b (cid:96) = − xr + x are thus negative. The Equations (12) reduce to the lossless power flow equations in thecase where the real part of the nodal admittance matrix is negligible G ≈ , i.e. lines arepurely inductive.We will focus on the so called DC approximation of this full AC power flow equations.This approximation is based on three assumptions [44]:1. Voltages vary little, i.e., | V i | ≈ const , ∀ i with respect to their base values,2. Angular differences are small, i.e., sin( ϑ i − ϑ j ) ≈ ϑ i − ϑ j , ∀ ( i, j ) ∈ E ( G ),3. Transmission lines are purely inductive, i.e., B ij (cid:29) G ij , ∀ ( i, j ) ∈ E ( G ).19ypically, these assumptions are fulfilled for high voltage transmission grids if the line loadingis not too large [27]. Using these approximation, Equation (12) reduces to P i = N (cid:88) k =1 | V i || V k | B ik (cid:124) (cid:123)(cid:122) (cid:125) K ik ( ϑ i − ϑ k ) , thus revealing the analogy to Equation (6). Application to resistor networks
Resistor networks are another example which may be described using linear flow net-works. They have been studied for a long time leading to many fundamental results ofgraph theory [42]. We will briefly introduce the theory of resistor networks and use the sym-bol ˆ= to refer to the corresponding quantity in the mathematical framework of linear flownetworks as introduced in section . For resistor networks, the flow along the graph’s edgesis a current flow (cid:126)i ∈ R L ˆ= (cid:126)F between nodes of different voltage (cid:126)V ∈ R N ˆ= (cid:126)ϑ . The line weightsare given by the inverse resistances, i.e. the conductances, of the lines G ∈ R L × L ˆ= K suchthat Equation (7) reads in this case (cid:126)i = GI (cid:62) (cid:126)V , where I is again the node-edge incidence matrix. Along the same lines, Equation (9) trans-lates to (cid:126)i in = I (cid:126)i. Here, (cid:126)i in ∈ R N ˆ= (cid:126)P is a vector of currents injected at the graphs’ nodes and the Equation isagain a manifestation of Kirchhoff’s current law. We may thus apply the same theoreticalframework to resistor networks. Applications to hydraulic networks
The same formalism can also be shown to apply water transport networks that we referto as hydraulic networks or pipe networks. Consider a hydraulic network consisting of pipesthat connect to each other at junctions. Then we form the underlying graph by assigning avertex to each of the junctions and put an edge between two vertices if they are connected20ia a pipe. The nodal quantity of interest in this case is the pressure (cid:126)p ∈ R N ˆ= (cid:126)ϑ . If weassume the pipes to be much longer than their radius r (cid:28) L and the flow across all pipesin the network to be laminar with a Newtonian, incompressible fluid flowing through it, wecan approximate the fluid flow (cid:126)Q ∈ R L ˆ= (cid:126)F across a pipe (cid:96) = ( i, j ) by the Hagen-Poiseuilleequation Q (cid:96) = K (cid:96) · ( p i − p j ) . Here, we collected different parameters describing the pipe and the fluid in the line parameter K (cid:96) = πr (cid:96) µL (cid:96) with the pipe radius r (cid:96) , the pipe length L (cid:96) and the fluid’s dynamic viscosity µ . Conservationof mass then requires that inflows and outflows balance as in Equation (6). Importantapplications of this framework are blood vessels in humans and animals [57], the vascularsystem of plants [28] or hydraulic networks [56]. For vascular networks, the system does notconsist of pipes but rather of smaller vascular bundles such that the scaling of line parameter K with the radius r does not necessarily exactly hold [29]. Applications to limit cycle oscillators
The linear flow model may be regarded as a linearisation of the
Kuramoto model whichnaturally appears in many cases, in particular when approximating weakly coupled oscillatorsystems near a stable limit cycle [30].Consider a connected, simple graph G = ( E, V ). The Kuramoto model describes a set ofweakly coupled oscillators with phase angles (cid:126)ϑ ∈ R N attached to the graph’s vertices that arecoupled via the graph’s edges through coupling constants K ij , ( i, j ) ∈ E ( G ), see e.g.[55]. Theoscillators’ tendency to synchronise through the coupling is counteracted by each oscillator’snatural frequency ω j that is written compactly as a vector (cid:126)ω = ( ω , ..., ω N ) (cid:62) ∈ R N . Thenthe dynamics of the phase angle ϑ i attached to node i, i ∈ { , ..., N } , reads˙ ϑ i = ω i − (cid:88) k K ik sin( ϑ i − ϑ k ) . As before, we fix an orientation of the graph’s edges and summarise the coupling coefficientsfor all edges ( i, j ) ∈ E ( G ) in the diagonal coupling matrix K ∈ R L × L , such that the21ectorised dynamics reads ˙ (cid:126)ϑ = (cid:126)ω − IK sin( I (cid:62) (cid:126)ϑ ) . (13)Here, I is again the graph’s node-edge incidence matrix (8) and the sine function is under-stood to be taken element-wise, i.esin( I (cid:62) (cid:126)ϑ ) = (sin([ I (cid:62) (cid:126)ϑ ] ) , ..., sin([ I (cid:62) (cid:126)ϑ ] L )) (cid:62) . Fixed points of the dynamics are defined by a vanishing time derivative ˙ (cid:126)ϑ = (cid:126)
0. Therefore,the equation characterising the phase angles at the fixed point (cid:126)ϑ ∗ reads (cid:126)ω = IK sin( I (cid:62) (cid:126)ϑ ∗ ) . If the angular differences on all edges are small, we may reduce this to the linear equationsin( I (cid:62) (cid:126)ϑ ) ≈ I (cid:62) (cid:126)ϑ , again retrieving an expression analogous to the discrete Poisson equa-tion (10). The second-order Kuramoto model
An extension of the Kuramoto model presented in Equation (13) is given by the second-order Kuramoto model that is also used frequently in power systems analysis to describesynchronising generators [52–54], where it is also referred to as Kuramoto model with inertia.The model contains an additional second-order time derivative of phase angles representingthe generators’ inertia and reads as¨ (cid:126)ϑ = − α ˙ (cid:126)ϑ + (cid:126)ω − IK sin( I (cid:62) (cid:126)ϑ ) . (14)Here, α = diag( α , ..., α N ) ∈ R N × N is a diagonal matrix incorporating the generators’ inertiaand friction coefficients [52] and the other quantities are defined the same way as for thefirst order Kuramoto model (13). The vector of frequencies in this model corresponds tothe power injections (cid:126)ω ∈ R N ˆ= (cid:126)P . Fixed points of the second order model with phase angles (cid:126)ϑ ∗ are characterized by both, first and second order time derivative vanishing ¨ (cid:126)ϑ = ˙ (cid:126)ϑ = (cid:126) (cid:126)ω = IK sin( I (cid:62) (cid:126)ϑ ∗ ) . Again, this model reduces to the linear flow model if phase differences at the fixed point aresmall sin( I (cid:62) (cid:126)ϑ ∗ ) ≈ I (cid:62) (cid:126)ϑ ∗ . 22 escription of link failures In this section, we will briefly review the analysis of link failures within the linear flowtheory setting. We will first demonstrate how the effects of a link failure may be approachedon the nodal level [31]. Assume that a link k = ( r, s ) with preoutage flow ˆ F k fails, whichdoes not disconnect the graph. This induces a change in the potentials (cid:126)ϑ (cid:48) = (cid:126)ϑ + ∆ (cid:126)ϑ by virtue of the discrete Poisson equation (10). Here, we introduced the vector of potentialchanges ∆ (cid:126)ϑ ∈ R N and a vector of potentials after the failure (cid:126)ϑ (cid:48) ∈ R N . The correspondingequation for the new grid reads as (cid:126)P = ( L + ∆ L )( (cid:126)ϑ + ∆ (cid:126)ϑ ) . Here, ∆ L is the change in the Laplacian matrix due to the removal of link k and takes theform ∆ L = K k I (cid:126)l k ( I (cid:126)l k ) (cid:62) . If we subtract the discrete Poisson equation for the old grid beforethe failure of link k from this equation, we arrive at the expression∆ (cid:126)ϑ = − ( L + ∆ L ) † ∆ L (cid:126)ϑ. Finally, we can use the Woodbury Matrix identity to rewrite the expression into the followingform [31] L ∆ (cid:126)ϑ = q k (cid:126)ν k , (15)where q k = (1 − K k ( I · (cid:126)l k ) (cid:62) L † I · (cid:126)l k ) − ˆ F k is a source term. Similar expressions appear naturally when analysing resistor networks andhave been studied, for example, in Refs. [48, 51]. After calculating the potential changesbased on this equation, the flow changes on a link (cid:96) = ( (cid:96) , (cid:96) ) are given by the followingequation ∆ F (cid:96) → (cid:96) = K (cid:96) · (∆ ϑ (cid:96) − ∆ ϑ (cid:96) ) . NETWORK ISOLATORS INHIBIT FAILURE SPREADING COMPLETELY
In this section we formally establish the existence of network isolators. To this end wefirst fix some notation. 23 undamentals and notation
We consider a linear flow network consisting of two parts, i.e. its vertex set V is writtenas V = V + V . We now label the nodes in V as follows without loss of generality1 , . . . , m : nodes in V that are connected to V m + 1 , . . . , n : nodes in V that are not connected to V n + 1 , . . . , n + m : nodes in V that are connected to V n + m + 1 , . . . , n + n : nodes in V that are not connected to V . Then the weighted adjacency matrix of the network can be written as A = A A A (cid:62) A A = a
00 0 with A ∈ R n × n , A ∈ R n × n , A ∈ R n × n and a ∈ R m × m . Furthermore, we definethe degree matrices D , D and d associated with the adjacency matrices A , A and a ,that is d kl = (cid:80) p a kp for k = l k (cid:54) = l , and the Laplacian matrices L = D − A of subnetwork 1, L = D − A of subnetwork2 and L of the whole system. Main theorem on network isolators
We can then formulate the main theorem on network isolators.
Theorem 2 (Network isolators completely suppress failure spreading between modules) . Consider a linear flow network consisting of two parts with vertex sets V and V and assumethat a single link in the induced subgraph G ( V ) fails, i.e. a link ( r, s ) with r, s ∈ V . If theadjacency matrix of the mutual connections has unit rank rank( A ) = 1 , then the flows onall links in the induced subgraph G ( V ) are not affected by the failure, that is ∆ F (cid:96) ,(cid:96) ≡ ∀ (cid:96) , (cid:96) ∈ V . The subgraph corresponding to the mutual interactions is referred to as network isolator . roof. Assume that the adjacency matrix of the mutual connections has unit rank rank( A ) =rank( a ) = 1. We first proof that for any vector (cid:126)y ∈ R n the following statement holds (cid:126)x = d − a
00 0 (cid:126)y = c , (16)where c ∈ R is some real number. This result can be obtained by writing (cid:126)x ∈ R n incomponents. For all j ∈ { , . . . , m } we have x j = (cid:80) k a jk y k (cid:80) k a jk . Since a has unit rank all its rows are linearly dependent such that we can write a jk /a k = a j /a for all k ∈ { , . . . , n } , such that a jk = a k a j /a . Hence, x j = a j /a × (cid:80) k a k y k a j /a × (cid:80) k a k = (cid:80) k a k y k (cid:80) k a k = x =: c, and all elements of the vector are equal. The remaining n − m elements of the vectorvanish, x j = 0 , ∀ j ∈ { m + 1 , ..., n } , because the corresponding adjacency matrix A hasonly zero entries at the respective positions.We now compute the impact of a failure of link k in G ( V ) via the discrete Poissonequation (15) L ∆ (cid:126)ϑ = q k (cid:126)ν k . We decompose this equation as well as the vectors ∆ (cid:126)ϑ and (cid:126)ν into two parts correspondingto the two parts of the network∆ (cid:126)ϑ = ∆ (cid:126)ϑ ∆ (cid:126)ϑ , (cid:126)ν = (cid:126)ν (cid:126) , where ∆ (cid:126)ϑ , (cid:126)ν ∈ R n and ∆ (cid:126)ϑ , (cid:126)ν ∈ R n . Then the lower part of Equation (15) correspondingto the vertices n + 1 , . . . , n + n reads L + d
00 0 ∆ (cid:126)ϑ = a
00 0 ∆ (cid:126)ϑ (17)25sing the notation established above. Using the prior result (16) and multiplying by thematrix d −
00 1 this equation can be rewritten as d −
00 1 L + ∆ (cid:126)ϑ = d − a
00 0 ∆ (cid:126)ϑ = c Now one can easily check via a direct calculation that∆ (cid:126)ϑ = c is a solution to this equation. Furthermore, this solution is unique as the linear system ofequation has full rank. This is most easily seen for Equation (17), as the matrix on the lefthand side is normal and positive definite.We have thus shown that the nodal potentials in V are shifted by the same constant c when a link in G ( V ) fails. Hence the flow changes are given by∆ F (cid:96) → (cid:96) = K (cid:96) (∆ ϑ (cid:96) − ∆ ϑ (cid:96) ) = 0 ∀ (cid:96) , (cid:96) ∈ V . Corollary 1 (Complete bipartite graphs are network isolators) . Consider a linear flownetwork consisting of two modules with vertex sets V and V and assume that a single linkin the induced subgraph G ( V ) fails, i.e. a link ( r, s ) with r, s ∈ V . If the subgraph G (cid:48) ofmutual connections between the two modules is a complete bipartite graph with uniform edgeweights K = K (cid:96) = K m , ∀ (cid:96), m ∈ E ( G (cid:48) ) , then the subgraph is a network isolator. If the wholegraph is unweighted, G (cid:48) always has uniform edge weights, thus a complete bipartite graph ofmutual connections always is a network isolator for any unweighted network. roof. If the subgraph G (cid:48) is complete and bipartite (ignoring all connections within bothinduced subgraphs G ( V ) and G ( V )), its adjacency matrix takes the form A (cid:48) = K · m × m (cid:62) m × m . We can immediately see that this matrix has unit rank, such that by theorem 2, G (cid:48) is anetwork isolator. Network isolators in non-linear systems
We will now demonstrate how to extend the concepts of network isolators from linearsystems to a certain class of non-linear networked systems with diffusive coupling. Let (cid:126)f ( L (cid:126)x ) = ( f ([ L (cid:126)x ] ) , ..., f N ([ L (cid:126)x ] N )) (cid:62) : (cid:126)x ∈ R N → (cid:126)f ( L (cid:126)x ) ∈ R N be a continuous function on the real numbers that depends on the product of Laplacianmatrix L and vector (cid:126)x . Here, [ L (cid:126)x ] j denotes the j -th row of the standard matrix-vectorproduct L (cid:126)x . We assume that the underlying network topology is again separated into twosubgraphs G ( V ) and G ( V ), see the beginning of this section. We further assume that f j (0) = 0 , ∀ j ∈ { , ...N } , i.e. each of the functions vanishes at the origin. Note that the functions f j ([ L (cid:126)x ] j ) can bedifferent and non-linear, as long as they vanish at the origin. Consider a dynamical systemof the form ˙ (cid:126)x = (cid:126)f ( L (cid:126)x ) (18)that admits a fixed point solution (cid:126)x ∗ with vanishing time derivative ˙ (cid:126)x = (cid:126) (cid:126) (cid:126)f ( L (cid:126)x ∗ ) . (19)Now add a perturbation vector ∆ (cid:126)P = ∆ (cid:126)P (cid:126) (20)to the system that has non-zero entries only at the nodes of the first induced subgraph G ( V )and assume that the dynamical system 18 relaxes to a new fixed point (cid:126)x (cid:48) with∆ (cid:126)P = (cid:126)f ( L (cid:126)x (cid:48) ) . (21)27hen the following corollary holds Corollary 2 (Isolation in non-linear systems) . Consider a non-linear dynamical networkedsystem of the form (18) that consists of two modules with vertex sets V and V which areconnected by a network isolator as of Theorem 2. Assume that the system admits a fixedpoint solution as given in Eq. 19. Assume that a small perturbation as in Eq. 20 is appliedto the nodes in the first induced subgraph G ( V ) and that the system relaxes to a new fixedpoint as in Eq. 21. Then the new fixed point has the following form (cid:126)x (cid:48) = (cid:126)x (cid:48) c(cid:126) , where c ∈ R is a constant. The second module is thus isolated against perturbations in the first module and viceversa in the sense that a perturbation in one module results in a constant shift in the othermodule.
Proof.
The proof is analogous to the proof of theorem 2. Applying the function (cid:126)f to Equa-tion 17 describing the fixed point in the non-perturbed subgraph G ( V ), we see that thesystem is still solved by (cid:126)x (cid:48) = ∆ (cid:126)x (cid:48) c(cid:126) . Approximate isolation for weakly non-linear systems
Even if not rigorously valid, we find that strong network isolation persists for an evenlarger class of non-linear systems that we will discuss in this section. Note that our analysishere closely follows a linear response theory analysis of Kuramoto oscillators that can befound in Ref. [30].Consider a networked non-linear dynamical system of the form˙ (cid:126)x i = (cid:126)f ( (cid:126)x i ) i + N (cid:88) k =1 A ik g ( x i − x k ) . (22)28ere, (cid:126)x ∈ R N is a vector of nodal dynamical variables, (cid:126)f is a differentiable function of self-interactions of these variables and (cid:126)g ( (cid:126)x ) is a differentiable function that depends only on thedifferences of nodal variables at neighbouring nodes. The strength of interactions is encodedin the graph’s adjacency matrix A . Assume that the system relaxes to a fixed point with˙ (cid:126)x i = 0 where (cid:126)x ( t ) = (cid:126)x ∗ . If we perturb the network locally at a node or an edge, we cancompute the change in this fixed point using linear response theory [30]: to leading order,we obtain a linear system as above.Assume that we perturb a single edge ( n, m ) by modifying its edge weight by a small number∆ A ij such that A ij → A (cid:48) ij = A ij + ∆ A ij ∆ A ij = i, j ) (cid:54) = ( n, m )∆ A if ( i, j ) = ( n, m ) . Assume that this modification causes a change of the fixed point by x ∗ j → x (cid:48) j = x ∗ j + ∆ x j , ∀ j ∈ { , ..., N } . We can expand the dynamics to leading order in terms of the new fixed point ∂f ( x ∗ j ) ∂x j ∆ x j + N (cid:88) k =1 A jk ∂g ( x ∗ j − x ∗ k ) ∂x j (∆ x j − ∆ x k ) + s j = 0 . Here, s j is a source term that vanishes if node j is not part of the edge ( n, m ), j (cid:54) = n, m . Thesum in this expression may be compactly written in terms of an effective Laplacian matrix˜ L N (cid:88) k =1 A jk ∂g ( x ∗ j − x ∗ k ) ∂x j (∆ x j − ∆ x k ) = [ ˜ L ∆ (cid:126)x ] j , where the Laplacian matrix has the off-diagonal entries˜ L jk = − A jk ∂g ( x ∗ j − x ∗ k ) ∂x j . Thus, if the underlying graph contains a network isolator, we can apply Theorem 2 to thesystem and see immediately that each component is (approximately) isolated against smallperturbations in the other one. In particular, this description applies to Kuramoto oscillators(see section ) perturbed at a few nodes or edges and powergrids described by AC load flowequations 12 subject to a link failure. We can thus get approximate isolation in both modelsas shown in Figure 4(e to h) for the AC load flow model.29
INEAR CONTROLLABILITY OF COMPLEX NETWORKS
We now turn to a different theoretical concept in complex networks research: the con-trollability of a network. In this section, we briefly analyse the influence of network isolatorson the controllabilty of complex systems with a linear dynamics. In general, we find thatintroducing a network isolator to a complex network has no generic influence on its control-lability.Consider a linear dynamical system on a network with N nodes with a state vector (cid:126)x ∈ R N whose dynamics is given by [47] ˙ (cid:126)x = A (cid:126)x + B (cid:126)u. (23)Here, A ∈ R N × N denotes the graph’s adjacency matrix, (cid:126)u ∈ R m is a (potentially time-varying) input vector that is supposed to achieve control of the network and B ∈ R N × m isthe control matrix. Then one definition of controllability is the following: Can we find a setof m driver nodes identified by the controllability matrix B such that the system may bedriven from any initial state (cid:126)x to any final state (cid:126)x f in finite time? If yes, the system is saidto be controllable and a measure of its controllability is given by the minimum number ofdriving nodes N d ≤ N necessary to achieve full controllability [47, 49, 50].We identify this set of driver nodes necessary for exact controllability for a small samplenetwork using a method due to Yuan et al. [47] who demonstrated that the minimum numberof driver nodes N d can be found by determining the multiplicity of the eigenvalues of thegraph’s adjacency matrix A [47]. Assume that the underlying network is undirected suchthat its adjacency matrix is symmetric as for the networks studied in this manuscript. Inthis case, we can calculate the algebraic multiplicity δ ( λ i ) for all eigenvalues λ i of this matrixto calculate the minimum number of driver nodes, N D , necessary to control the network (cf.Eq.4, Ref. [47]) N D = max i [ δ ( λ i )] . (24)This approach has the advantage that the driver nodes necessary to control the network, i.e.the controllability of a network, may immediately be identified, which is more complicatedwhen using the classical Kalman rank condition [47].In Figure S8, we illustrate a potential application of this formalism to network isolators.30he adjacency matrix of the graph reported in panel A has the eigenvalue λ M = − δ ( λ M ) = 2, while all other eigenvalues have multiplicity one. An eigenvalue λ M = − δ ( λ i ) = 1 , ∀ i , which implies that the graph can becontrolled by a single node (colored red). Therefore, in this case, the controllability of thenetwork is increased after constructing the isolator. We emphasize that the network isolatorprevents only flow changes, but not flows from passing as demonstrated in panels B-C andE-F.For the remaining network isolators constructed in throughout this manuscript, we didnot find any influence of the introduction of network isolators on the controllability of theunderlying network and thus conclude that isolators do not generically influence networkcontrollability. COMPUTATIONAL METHODSCreating graphs with strong or weak inter-module connectivity
We introduce a model to create ensembles of graphs consisting of two subgraphs withweak or strong interconnectvity, see Figures 1 and 2. We start with two disconnected Erd˝os-R´enyi random graphs G ( N , p ) and G ( N , p ), where N denotes the number of nodes inthe grid and p the probability that two randomly chosen nodes are connected by an edge[35]. Then we randomly choose n = [ c · N ] nodes v = { v , ..., v n } in G and n = [ c · N ]nodes w = { w , ..., w n } in G . Here, c ∈ [0 , ⊂ R is a constant representing the shareof nodes connecting to the other subgraph and [ · ] denotes the nearest integer. Out of allpossible edges e = { ( v , w ) , ..., ( v n , w ) , ..., ( v n , w n ) } between the two sets of nodes v and w , we randomly add a share of µ ∈ [0 , µ controls the connectivity of thetwo subgraphs G and G : They remain disconnected for µ = 0 and they are connectedvia a complete bipartite graph for for µ = 1. For c = 1 and µ = p = p we recover a31ingle Erd˝os-R´enyi random graph with N = N + N nodes. Note that this procedure is inprincipal not limited to ER random graphs. We apply it to study other types of graphs asshown in Supplemental Figure S1. Perturbing network isolators
The robustness of network isolators to structural perturbations is analysed as follows.Let G = ( E, V ) be a graph whose nodes are split into two subsets V and V . Furthermore,let A be the part of the graph’s weighted adjacency matrix that encodes the mutualconnections between the two parts as described in theorem 2. Without loss of generality wecan order the nodes of the network in such a way that the matrix has the structure A = (cid:126)a · · · (cid:126)a m (cid:126) · · · (cid:126) (cid:126) · · · (cid:126) (cid:126) · · · (cid:126) . (25)According to theorem 2, a perfect network isolator is found if rank( A ) = 1, i.e. if all vectors (cid:126)a , . . . , (cid:126)a m are linearly dependent.To investigate the robustness of network isolators, we start from a unit rank matrixrank( A ) = 1 and perturb it iteratively. In each step we choose one of the vectors (cid:126)a i , i =1 , . . . , m at random and perturb it according to (cid:126)a (cid:48) i = (cid:126)a i + (cid:126)e (cid:107) (cid:126)a i (cid:107) . The elements of theperturbation vector (cid:126)e are chosen uniformly at random from the interval [ − β, β ], where β isa small parameter, here β = 0 . A from a unit rank matrix is quantified usingits coherence statistics [39], ξ ( A ) = 1 − min i,j (cid:104) (cid:126)a i , (cid:126)a j (cid:105)(cid:107) (cid:126)a i (cid:107)(cid:107) (cid:126)a j (cid:107) , (26)where (cid:104)· , ·(cid:105) denotes the standard scalar product on R n and (cid:107)·(cid:107) denotes the (cid:96) -norm. For amatrix A of unit rank we have ξ ( A ) = 0 as all vectors are linearly dependent. For vectorsdeviating from linear dependence, the measure increases until it reaches its maximum valueif two vectors are linearly independent with ξ ( A ) = 1.To create Figure 2(e), we repeated this process 1000 times starting from the perfectisolator shown in panel c. Edge weights were randomly chosen from a normal distribution N (10 ,
1) with mean µ = 10 and variance σ = 1 except for the isolator, where we choose32our groups of four edges having the same weight such that initially rank( A ) = 1. For eachperturbed network, we evaluate ξ ( A ) and the ratio of flow changes R according to Eq. (3)averaged over all possible trigger links (cid:96) and distances d . For a perfect isolator, this ratiovanishes due to a vanishing numerator. Power grid data and cascade model
Power grid data has been extracted from the open European energy system model PyPSA-Eur, which is fully available online [19]. The model includes the topology as well as thesusceptance b (cid:96) and the line rating F max i → j for each high voltage transmission line in Europe.We consider the Scandinavian synchronous grid spanning Norway, Sweden, Finland andparts of Denmark. This grid is coupled to other synchronous grids (central European grid,British grid and Baltic grid) only via high voltage DC transmission lines. Power flow on theselines are actively controlled and can thus be considered constant, thus leading to constantreal power injections at the coupling nodes. The Scandinavian grid has 269 nodes and 370edges, counting multiple-circuit lines only once.Cascading failures are simulated for fixed power injections P i for each node correspondingto an economic dispatch for the entire PyPSA-Eur model that includes a security margingiven by the constraint | F i → j | ≤ . · F max i → j . The cascade is triggered by the failure of asingle line ( r, s ) which is effectively removed from the grid. The simulation then proceedsstep-wise; In each step, we first calculate the nodal phase angles ϑ i and real power flows F i → j for all nodes and lines, respectively, by solving the continuity equation P i = (cid:80) j F i → j with F i → j = K ij ( ϑ i − ϑ j ). Then we check for overloads: Any line ( i, j ) with | F i → j | > F max i → j undergoes an emergency shutdown and is removed from the grid. The simulations arestopped when no further overload occurs or when the grid is disconnected.Note that this mechanism for cascading failures is different from the cascading failuremechanism typically analysed in node capacity load models (see e.g. Refs. [45, 46]). Theredistribution of nodal loads or flows after failures in such models is typically based onneighborhood of nodes, on shortest path betweenness measures or on other ‘intelligent’ re-distribution schemes whereas the redistribution of flows after failures in linear flow networksor power grids studied using AC load flow analysis are given by the physical laws governingelectrical networks. Furthermore, usually nodes - not edges - are assumed to fail, which is33ot the typical case in real power grids. Processing leaf data
The leaf venation network is based on a microscopic recording of a leaf of the species
Bursera hollickii provided by the authors of Ref. [20]. Edge weights K ij are assumed toscale with the radius r ij of the corresponding vein ( i, j ) as K ij ∝ r ij according to the Hagen-Poisseuille law, see Ref. [29] for a detailed discussion. We used the radius in pixel scannedat a resolution of 6400 dpi. 34 . . . . . . µ − − R (a) ER graphs G (50 , .
3) 0 . . . . . . µ − − R (b) ER graphs G (30 , . . . . . . . µ − − R (c) BA graphs N = 40 , k = 4 0 . . . . . . µ − − R (d) Random 4-regular graphs µ = 0 . µ = 0 . µ = 0 . µ = 0 . Figure S1.
Averaged ratio of flow changes decays with high and low connectivityfor different random graphs.
All panels show ratio of flow changes R averaged over all linksand distances against connectivity parameter µ (see methods) along with corresponding graph forlow values of the connectivity parameter. (a) Two ER graphs with parameters N i = 50 , p i = 0 . µ = 0 .
02 at a randomly chosen share of c = 0 . N i = 30 , p i = 0 . , µ = 0 . , c = 0 .
2. (c) A similar scaling isobserved if two BA random graphs with parameters N i = 40 , k i = 4 are connected with probability µ = 0 .
016 at a randomly chosen share of c = 0 . N = 50 , µ = 0 . , c = 0 .
2. Blue linerepresents median value over all distances and shaded region indicates 0 .
25- and 0 . − − − − hh | ∆ F p → q | i ( p , q ) ∈ O / S d i ‘ µ = 0 . (a) ‘ ∈ S‘ ∈ O − − − − µ = 0 . (b) ‘ ∈ S‘ ∈ O − − − − µ = 0 . (c) ‘ ∈ S‘ ∈ O d − − − R ( d ) (d) d − − − (e) d − − − (f) Figure S2.
Ratio of flow changes depends weakly on distance.
We examine the scalingof link flow changes with distance for two ER random graphs G (120 , .
02) that are connected at c = 0 . µ = 0 .
02 (left), µ = 0 . µ = 0 . R ( d ) revealsa weak dependence of the ratio on distance. Blue line represents median value over all distancesand shaded region indicates 0 .
25- and 0 . igure S3. Increasing or decreasing connectivity between more than two modulesreduces failure spreading equally well.
Here, we demonstrate a possible extension of thesynthetic network model described in the Methods section to more than two modules. For eachpanel, we simulate a single link failure (red) that results in flow changes (colour coded). (a to c)Three ER random graphs G (30 , .
3) (right), G (50 , .
2) (bottom) and G (40 , .
4) (top left) that aremutually interconnected with probability µ = 0 .
05 at 20 percent, i.e. c = 0 .
2, thus resulting inthree mutually weakly connected modules. (d to f) Connecting the same modules as shown in (Ato C) with probability µ = 0 .
85, thus resulting in strong inter-module connectivity, reduces failurespreading equally well. igure S4. Networks isolators can be generalised to network consisting of more thantwo modules. (a) Topology of a network consisting of three ER random graphs G (40 , .
4) (left), G (20 , .
3) (top) and G (30 , .
2) (bottom right) that are mutually connected through network iso-lators. (b to d) Link failures in each of the individual subgraphs (red lines) do not change flows(colour code) in any of the other subgraphs. a) − − − ξ ( A )10 − − − R (b) Figure S5.
Robustness of network isolators shows the same scaling with perturbationsfor different graphs.
Robustness of network isolators measured by ratio of flow changes R averaged over all links against measure of perturbations to network isolators ξ ( A , ). (a) Graphcreated from the graph ensemble and shown in Fig. S1C was modified in such a way that it containsa network isolator connecting five nodes from one part to five nodes of the other part througha bipartite connectivity structure. Edge weights are drawn randomly from a normal distribution N (10 ,
1) except for the network isolator where the randomly chosen weights of five edges starting inthe same node and connecting to all connecting nodes in the other part were chosen as basis weightsfor all other connections between the two parts. (b) The isolator robustness shows qualtiativelythe same scaling as for the 6-regular graph shown in Fig. 1(c). Perturbations were applied in 1000repetitions choosing a perturbation strength of α = 0 .
05. Dotted line takes into account the factthat the curve goes through the point ξ = R = 0 for a perfect isolator. ∆ F | (a) | ∆ F | (b) − − − − − − − − Figure S6.
Network isolators do not increase grid vulnerability. (a) Failure of a linkwith unit flow in the Scandinavian grid before the construction of the network isolator yields astrong response in terms of absolute flow changes | ∆ F | . (b) After adding two links to create anetwork isolator (blue shaded region, see Figure 3(c)), we simulate a failure of one of the links in the isolator. We observe that both, the failure within the isolator (b) as well as a failure in theinitial grid in close proximity to the location where the isolator is constructed (a) yield a similareffect. We thus conclude that introducing the network isolator will not make the network morevulnerable compared to the network without the isolator. However, a failure in the isolator maypotentially affect the whole network. ∆ F | (a) | ∆ F | (b) | ∆ F | (c) | ∆ F | (d) | ∆ F | (e) | ∆ F | (f) − − − − − − − − − − − − − − − − − − − − Figure S7.
Network isolators may be realised in various real-world power grids.
Allgrid topologies and line susceptances were extracted from the open European energy system modelPyPSA-Eur, which is fully available online[19]. (a,c,e) Initial failure of a link (red) with unit flowresults in flow changes in the whole network for Scandinavia (a,c) as well as the central Europeangrid (e). (b,d,f) After introducing network isolators to the grids, failure spreading to other partsof the network is completely stopped. The construction of isolators follows the “recipes” illustratedin Figure 4, namely with (b) following the construction in Fig. 4(b), (d) following the constructionin Fig. 4(a), and (f) following the construction in Fig. 4(d). a) | F | (b) | ∆ F | (c)(d) | F | (e) | ∆ F | (f) − − − − − − − − Figure S8.
Isolators do not generally prevent the controllability of a network. (a) Anexample of an undirected network with two weakly connected components that requires N D = 2driving nodes (in orange) to be controlled. This can be calculated from the graph adjacency matrix,which has, by construction, an eigenvalue λ M = − δ ( λ M ) = 2 (SeeEq. 24 and Ref. [47]). (d) After adding a few links to create a network isolator, we have N D = 1 andonly one node (colored orange) is necessary to control the entire network, i.e. the network isolatorhas in this case increased the controllability of the network. (b) We show the flows obtained byour linear flow model for a single source of power P = 1 at the node colored in red and a singlesink with P = −1 at the node colored in blue. The resulting (absolute) flows are color-coded:The flow can easily reach from the red node to the blue node. (e) Adding the isolator, flow canstill propagate freely from the source node (red) to the target node (blue) in the same way as inpanel B. Hence, the isolator does not prevent the propagation of flows. (c) Simulating the failureof a single link (red), we observe that flows do also change in the other part of the network. (f)Conversely, the isolator does prevent propagation of flow changes caused by a link failure in theright part of the network to its left part.