Locating line and node disturbances in networks of diffusively coupled dynamical agents

A wide variety of natural and human-made systems consist of a large set of dynamical units coupled into a complex structure. Breakdown of such systems can have a dramatic impact, as in the case of neurons in the brain or lines in an electric grid, to name but a few. Preventing such catastrophic events requires in particular to be able to detect and locate the source of disturbances as fast as possible. We propose a simple method to identify and locate disturbances in networks of coupled dynamical agents, relying only on time series measurements and on the knowledge of the (possibly Kron-reduced) network structure. The strength and the appeal of the present approach lies in its simplicity paired with the ability to precisely locate disturbances and even to differentiate between line and node disturbances. If we have access to measurement at only a subset of nodes, our method is still able to identify the location of the disturbance if the disturbed nodes are measured. If not, we manage to identify the region of the network where the disturbance occurs.


Introduction
Despite recent advances in network science [1,2], the understanding of large, networked dynamical system is still incomplete, even though such systems can play a major role in our daily life. Large scale events therein should be monitored or mitigated, if not completely avoided. Preventing catastrophic events in a networked system requires, in particular, an accurate and reliable assessment of its state in real time. Indeed, the early detection of a source of disturbances and its precise location allow preventive measures to be taken.
In networks, the most fundamental disturbances can affect either of their two basic constituents, namely either nodes or lines. Nodal perturbations, occurring mostly as variations of intrinsic velocities of the disturbed agents, act as an additive perturbations, whereas disturbances on lines, which are modifications of the coupling parameters, are represented as multiplicative disturbances in the model, which are much harder to tackle analytically. This partly explains the scarcity of the literature covering line disturbances. Nevertheless, such perturbations are at least as important as nodal perturbations and, in our opinion, should be investigated in details.
Systems that can be modeled as a set of diffusively coupled dynamical agents where state estimation is needed for security, comfort, or diagnostic reasons are extremely diverse. To name but a few, they range from neurons in the brain to transportation networks and disease spreading models [3][4][5]. The case of high voltage electrical grids is an important example of such systems and has attracted a lot of work along the last years, especially in the context of the ongoing energy transition in many parts of the world [6,7]. We will systematically refer to the power grid and often use its terminology in this manuscript for illustrative purposes, even though our considerations extend to a much larger class of systems.
In the context of power grids, recent applications of results from network science and dynamical systems theory allowed the development of new techniques aiming at a fast assessment of the system's state. Such improvements are based on the measurement of some quantities such as voltage amplitudes, phases, and frequencies, which are nowadays widely accessible on a high time resolution thank, for instance, to phasor measurement units (PMUs). Part of the work on state estimation, such as references [8][9][10], is mostly focused on the estimation of the current operating state of the whole system. This is performed using various techniques, such as estimation of the eigenmodes of the network through probing signals [8] or resonance methods [10], or based on the measurement of the response of the system to ambient noise [9]. Closer to our interest in this manuscript, references [11][12][13][14][15] proposed some approaches aiming at locating the source of a disturbance. For networks composed of areas with weak inter-area connections and strong intra-area connections, reference [11] uses the residues of an estimated transfer function in order to locate a nodal disturbance. Other approaches rely on discrete wavelet transform [12,13] or logistic regression [14] of PMU measurements in order to identify the source of a disturbance. Most of these methods are designed to locate a nodal disturbance only.
Despite their relevance in electrical networks [16][17][18], the case of line disturbances led to fewer results compared to nodal disturbances [19,20], and less focus was put on assessing the impact of a line perturbation. Methods based on the disturbance propagation [15] would be able to determine the area where a disturbed line is located, but to the best of our knowledge, there are no such methods specifically designed to locate line disturbances.
Recent works analyzed the case of line failures in power grids, e.g., following an adversary attack [21][22][23][24][25][26][27]. These approaches assumes full knowledge of the system prior to the line failure, and analyzes the mismatch between the expected and actual voltages and powers to recover the set of failed lines, by optimizing over the potential fault locations. On top of the line failures, the authors of [22][23][24][25] assume that the attacker is able to block data acquisition in the area of the fault, and sometimes even sends some spurious data to the monitoring device. In the approach of [26,27], the authors assume only partial access to measurements, and find the most likely location for the failure. However, in these approaches, the authors assume that the system operator is aware of the attack, and their analysis follows the line failure, which means that dramatic events already happened.
In this manuscript, we propose to detect and locate disturbances whose impact on the system is not (yet) threatening its operation, without any a priori knowledge of the faults' characteristics or approximate location. Our method relies on the measurements of the state of the dynamical agents composing the system and on the knowledge of the interactions between them when the system is unperturbed. This means that, in theory, our method can be applied online, without requiring to know if the system is perturbed or not. Conveniently, we cover the case of disturbed nodes as well as disturbed lines in a unified framework. Assuming that the amplitude and the rate of change of these disturbances are not too large, we propose a way to locate the faulty elements (nodes and/or lines) in the network, based on measurement of the nodes' trajectories. In a similar spirit as references [22][23][24][25], but in a more general framework, we locate the source of the disturbances by inspection of the mismatch between expected and actual agents' velocities. If we have access to measurements at the faulty elements, we are able to precisely locate it, even though such slow disturbance will spread throughout the whole network and have a similar impact on each agent. In case of partial measurements, we are nevertheless able to determine the area of the network where the faulty element is located. Furthermore, our implementation is adapted to a continuous monitoring of the system and can be used to preventively detect and locate deficient elements before the occurrence of a more significant event, such as a cascading failure, and that at a rather low cost as our method does not rely on any optimization but only on a single matrix-vector product at each time step.
The manuscript is organized as follows. We recall some preliminary tools and give a description of the type of models considered in section 2. In sections 3 and 4, we detail the method to locate the faulty element for a single and multiple disturbances, and numerical illustrations are given in section 5.

Notations and framework
Throughout this manuscript, we denote by e i the ith vector of the canonical basis, with 1 at index i and zero everywhere else, and e ij = e i − e j . We write 1 n (resp. 0 n ) the vector of ones (resp. zeros) of length n, where the subscript is omitted when the dimension is clear from the context. Finally, we denote by I the identity matrix.

Networked diffusive systems
We model the dynamics of a set of n ∈ N diffusively interacting agents as where the real numbers m i , d i , and ω i are respectively the inertia, the damping, and the natural velocity of agent i, a ij ∈ R is the (i, j)th element of the weighted adjacency matrix of the interaction graph, and f : R → R is an odd, differentiable coupling function, such that f (0) > 0. We assume the interaction graph to be undirected, i.e., a ij = a ji (the directed case is discussed in the appendix A). Note that we will be interested in fixed points of equation (1), i.e., when the left-hand side vanishes, meaning that all our results apply to models with higher order time derivatives. Assume equation (1) has a fixed point, denoted x * ∈ R n . When the system is not too heterogeneous, i.e., the spread of {ω i } is small relatively to the magnitude of the couplings, the fixed point components are small (in absolute value) and are well-approximated by the linearized equation where L ∈ R n×n is the Jacobian matrix defined as, where F i (x) is given by the sum in the right-hand side of equation (1). This Jacobian matrix satisfies the properties (i) j L ij = 0 and (ii) L ii = − j L ij ∀i, it is then a (weighted) Laplacian matrix corresponding to the interaction graph of the system equation (1), with weights related to the strength of the couplings between agents. This approximation is common in the context of power grids, where it is called DC power flow equations [28, section 9.7]. Obviously, for linear coupling functions, equation (2) is not an approximation.
Remark. Throughout this manuscript, we will refer to dynamical agents as nodes and to the connections between them as lines. In the context of graph theory, these two objects are usually called vertices and edges respectively.

Disturbances
We consider disturbances either at nodes, yielding the time varying natural velocity or on lines, described as varying coupling where we take ξ n,l (t) to be any function of time, whose rate of change is typically smaller than any intrinsic time scale of the system equation (1). Formally, we require that the disturbance's rate of change is smaller than the rate at which its effect is damped at the agents and smaller than the speed at which it spreads throughout the network, i.e., for all i, j, where 0 = λ 1 < λ 2 < · · · < λ n are the (real) eigenvalues of the Jacobian matrix L. This first requirement is reasonable in the sense that a disturbance that varies too fast will not spread in the network, and thus either its source is rather easy to locate or it has no influence on the observed nodes.
In order to guarantee that the system remains in the vicinity of the initial fixed point, we also require that the amplitude of ξ n,l is not too large, These last assumptions allow us to consider that the system remains in the linear approximation regime. Again, this assumption is reasonable because a disturbance with too large amplitude would compromise the normal operation of the system, and its location then becomes a secondary problem. We will assume equations (6) and (7) to be satisfied throughout the manuscript.

The Sherman-Morrison-Woodbury formula
Our results rely heavily on the Sherman-Morrison-Woodbury formula [29, section 2.1.4], giving an explicit formulation for the inverse of the rank-k perturbation of an invertible matrix A, where U ∈ R n×k and V ∈ R k×n characterize the rank-k perturbation. We emphasize that the Sherman-Morrison-Woodbury formula applies to the pseudoinverse of Laplacian matrices (even though such matrices are singular), provided that the columns (resp. rows) of U (resp. V) are orthogonal to the kernel of A.
In the case where k = 1, equation (8) reduces to where u and v are vectors characterizing the rank-1 perturbation. Equation (9) is usually referred to as the Sherman-Morrison formula.

The Kron reduction
In the context of electrical networks, it is possible to rewrite the equation governing the currents (the power flow equations [30]) with respect to a restricted number of voltage variables through Kron reduction [31]. This is done by taking the Shur complement [29, section 3.2.11] of the coupling matrix with respect to a subset of nodes and adapting the velocities accordingly. Indeed, this reduction is not restricted to power grids and can be applied to the coupling matrix of any network. Partitioning the nodes in two sets I g = {1, . . . , n g } (the nodes that are not reduced) and I c = {n g + 1, . . . , n g + n c } (the ones that are reduced), we write the elements of equation (2) in block form (reordering indices if necessary) The Kron-reduced coupling matrix is then the Schur complement of L cc in L, and applying a similar reduction to the vector of natural velocities gives This yields the Kron-reduced version of equation (2), which is restricted to the non-reduced nodes which allows to solve the equations for the subset of nodes I g .
Note that if the coupling matrix is Laplacian, then its Kron-reduced counterpart is also Laplacian.

Locating a single disturbance
We assume the disturbance to vary sufficiently slowly, meaning that its effect will propagate throughout the network. The detection of the disturbance can then be performed at any point of the system (e.g., through Fourier transform of the signal), but locating it requires a more elaborate procedure. We use the same method to identify disturbances at nodes or on lines. In order to perform this we need to know the Jacobian matrix L of the system, possibly Kron-reduced if some nodes are not measured, i.e., we need to be informed about the unperturbed state of the system. Under our assumption of relatively small amplitude in ω, a reasonable approximation of a stable fixed point of equation (1) is where † denotes the Moore-Penrose pseudoinverse of the matrix [29, section 5.5.2]. Note that it then satisfies LL † = I − n −1 11 because the Kron reduction L r of L is a Laplacian matrix. The effect of the disturbances we consider translates into time-varying (reduced) Jacobian matrixL r (t) and/or velocity vectorω r (t). As we assumed the amplitude of the disturbance to be relatively small, the evolution of the state of the system is reasonably approximated as x g (t) ≈ [L r (t)] †ωr (t). If we have access to the time series of x g (t) and the knowledge of L r , we can compute which we refer to as the frequency mismatch. We show now that the time series of ψ(t) are able to locate the disturbed element with as much accuracy as can be expected. We distinguish three cases: a disturbance at a node (which can be reduced or not); a line disturbance between two nodes that are not Kron-reduced; and a line disturbance with at least one end-node that is Kron-reduced.

Nodal disturbance
This case is the easiest to treat. The Jacobian matrix is constant in timeL r (t) = L r and only one component of the vector of velocities is time-varying. There are now two possible cases. If the perturbed node is not reduced, we can locate it exactly. Indeed, in this case, the reduced velocity vector has the form which, once plugged into equation (15), yields Provided that n is not too small, the amplitude of ψ(t) will be significantly larger at the perturbed node, allowing then to identify it. If the perturbed node is Kron-reduced, the reduced velocity vector is expressed as Again, plugging this into equation (15), we see that all non-reduced nodes that are neighbors of the reduced component are significantly impacted by the disturbance. We can then identify the area of the network in which the faulty node is located, but we cannot identify the node exactly. This is not surprising as we do not have a good enough accuracy in the measurements.

Line disturbance between non-reduced end-nodes
This is the most interesting case because we are able to locate the faulty line exactly. It also covers the case where no Kron reduction is applied. Here the faulty line has no influence on the reduced velocity vector and is a rank-1 perturbation of the reduced Jacobian matrix, where we adapted the size of e ij and the indices i and j in order to comply with the Kron reduction. Introducing equation (19) into the computation of the frequency mismatch yields where we used Sherman-Morrison equation (9) at the second line and defined where we used that LL † e ij = e ij which follows from the property of the pseudoinverse discussed below equation (14). One sees that the only time varying components of ψ(t) are at the two end points of the perturbed line, which allows to locate it exactly. This case is illustrated in figure 1(d), which shows the time series of ψ i (t) for nodes 1 to 9 of the network shown in figure 1(a), when the coupling of the green line is varying with time. One sees that the perturbed line is unambiguously identified.
Remark. We notice that the larger the weight of the disturbed line (in the Jacobian matrix L), the more accurate the identification. Indeed, if the value L ij is small, its variation will hardly be noticed on the  (15). The correspondence between indices and colors is given by (1,2,3,4,5,6,7,8,9) ↔ (blue, orange, green, red, purple, brown, pink, gray, yellow). (b) Disturbance at the non-reduced node 6. (c) Disturbance at the reduced node 11. (d) Disturbance on line (7,8). None of the end-nodes are reduced. (e) Disturbance on line (2, 10) where node 10 is reduced. (f) Disturbance on line (11,12), where both nodes are reduced. dynamics of the agents. This means that our method is the most reliable for disturbances occurring on the most important lines of the system.

Line disturbance with at least one reduced end-node
In this situation, the exact location of the perturbed line cannot be determined based on the measurements because at least one of the end points is hidden to the observer. We detail the case where the two end points of the faulty line are in the Kron-reduced set. The mixed case where only one of the end points is reduced is very similar. All of our computations are inspired from [19, sections V.1 and V.2] and rely on an appropriate use of the Sherman-Morrison formula equation (9) and of the formulae for the Kron reduction equations (11) and (12).
Assuming that i and j are in the Kron-reduced set, one gets time-varying reduced velocity vector and time-varying reduced Jacobian matrix where we defined and e ij is adapted according to the Kron reduction, and we used that L is symmetric, i.e., L gc = (L cg ) . Now in the same spirit as before, the frequency mismatch is where again, we used the Sherman-Morrison formula, equation (9), and where we gathered all time-varying quantities in γ(t).
For the case where node i is not reduced and node j is, a similar calculation gives with e i and e j being adapted to the Kron-reduced system, and where, again, we assembled all time-varying quantities inγ(t).
In these two cases, the effect of the perturbed line will be measured [through ψ(t)] on the non-reduced nodes connected to the reduced component containing the perturbed line. This can be verified by inspection of the submatrix L gc and is seen in figures 1(e) and (f), showing the time series of ψ i (t), when the disturbed line is at the blue and orange lines respectively. One cannot identify precisely its location, but it is possible to identify the reduced component to which it belongs, which is the most we could expect from the available measurements. In figure 1(e), where the perturbed line connects a reduced node to a non-reduced one, we observe that the amplitude of the frequency mismatch ψ(t) is much larger at the non-reduced end of the disturbed line than at any other non-reduced node. This makes sense if line couplings are all similar, the ith component ofṽ (equation (28)) is likely to be significantly larger than its other components, which will translate as a larger amplitude of variation in ψ i (t). However, in full generality, we cannot guarantee that the amplitude of ψ i (t) will be larger than all other components of the frequency mismatch.
In summary, by simply computing the frequency mismatch, we are able to locate node and line disturbances, at locations in the system that are not necessarily measured.

Multiple disturbances
Interestingly, we observe that our approach allows to locate multiple disturbances acting simultaneously, provided they have sufficiently different characteristics so one is able to distinguish them. We treat here the case of two simultaneous disturbances, the case with more disturbances extends naturally from it. We show the cases where the disturbances occur in non-reduced areas. The reduced case works similarly as in the single disturbance case.
For two nodal disturbances ξ 1 (t) and ξ 2 (t), at non-reduced nodes i and j respectively, the calculation of the frequency mismatch yields, Provided the disturbances ξ 1 and ξ 2 have sufficiently different characteristics, one can identify their location if n is not too small. If one of the nodes is reduced, only the area of the disturbance can be located. Two simultaneous line disturbances can be modeled as the following rank-2 perturbation of the Jacobian matrix,L where lines (i, j) and (k, ) are perturbed. As in the single disturbance case, one can then calculate ψ(t), which yields, using the Sherman-Morrison-Woodbury formula equation (8), Recall that ω, the columns of U, and the row of V t are orthogonal to the kernel of L.
Parsing equation (31), we see that Furthermore, we see numerically that the matrix X is close to be diagonal. Namely, the off-diagonal elements of X are smaller (in absolute value) than its diagonal elements by at least one order of magnitude. This yields where X pp is the pth diagonal element of X. The disturbances can then be located unequivocally, provided they have sufficiently different characteristics. Even though we have not been able to certify analytically that X is almost diagonal, we explain it heuristically as follows. Analytically, everything boils down to showing that V t L † U is (close to) diagonal. Let us forget the time dependence of V for now (ξ 1 = ξ 2 ≡ 1) and compute the (α, β)th term of this matrix, where we denote by (i α , j α ) the disturbed lines, where λ γ is the γth eigenvalue of L and u γ,k denotes the kth component of the associated eigenvector. At first sight, one could think that the main contribution to equation (34) comes from eigenmodes with small eigenvalues, i.e., γ close to 2. However it appears that, for the corresponding eigenvectors, the components of neighboring nodes are rather similar (see figure 2, showing some eigenmodes on a large-scale network modeling the European power transmission grid). It is particularly true in large-scale networks where these low laying eigenmodes spread over the whole network, at odds with high laying ones that are rather localized on a few nodes. One could think then, that the contribution comes from these high laying eigenmodes as the components of neighboring nodes might differ significantly [32]. However, these contributions also remain small due to the denominator in equation (34) where large values for λ γ bring down the contribution. The diagonal terms in X however are not vanishing due to the identity matrix appearing in equation (31). We performed a numerical investigation of the matrix VL † U for the European grid. On average over all possible pairs of lines, the absolute value of its off-diagonal term is 5.2 × 10 −6 and it reaches 5.9 × 10 −2 for the worst configuration. This seems to validate that the assumption we made on X holds, at least for interaction graphs based on infrastructures.
Even though we have not been able to ground our observations in analytical results, it is remarkable that such a naive approach is able to isolate various disturbances in a complex system.

Implementation and numerical validation
In order to validate our approach, we simulated three instantiations of the dynamics equation (1) with three different types of line disturbances, on three different networks.

Algorithm
We propose here and summarize in figure 3 the course of an algorithm that would implement our method.
First, with the Jacobian matrix L and the times series {x(t)} T t=1 as input, one can straightforwardly compute the frequency mismatch {ψ(t)} T t=1 . Then, the algorithm needs to identify groups of outlier nodes, whose frequency mismatch trajectory are significantly different from the bulk of all measured nodes. On top of that, these outlier nodes need to be grouped in sets (I 1 , . . . , I p ) of qualitatively similar trajectories, e.g., blue and green (resp. orange and red) curves are grouped together in figures 5 and 6. Defining what 'significantly different' means is not always trivial and might be problem dependent, but it is necessary in order to identify disturbed nodes. Furthermore, identifying 'significantly different' curves and classifying them by similarity might by a challenging problem for data analysts, but these issues are beyond the scope of our expertise in general and of this manuscript in particular.
Once the nodes have been classified, each subset of nodes should correspond to a disturbance in the system. In figure 3, we propose a way to determine the type of disturbance associated to each subset I j .

PanTaGruEl
The PanTaGruEl model of the interconnected European grid [32,33] is built on publicly available data of geolocalization of power system elements, and parameters (admittances, productions, loads, . . . ) are Figure 3. Block diagram of a suggested algorithm implementing our method. Note that this manuscript covers the first block only, whose task is to ease the work of downstream blocks. This does not mean that these subsequent blocks are trivial. Reasonable noise was added at each node in order to challenge the robustness of our approach. The disturbed lines are exactly identified as the ones whose end nodes have a significant amplitude for ψ i . reconstructed based on standard assumptions. It consists of 3809 nodes and 7343 lines. For our simulations, we consider an oscillating disturbance, ξ l (t) = ξ 0 sin(ω m t), satisfying the conditions detailed in section 2.2. From a practical point of view, this allows us to tune the time scale of the disturbance, which is the oscillation frequency ω m in order to guarantee that it is sufficiently slow with respect to all time scales of the network. In the following, we compare the disturbance detection using our method relying on ψ(t) or by simply measuring the agents' frequenciesẋ(t) and value x(t). The left panels of figure 4 shows the maximal amplitude of ψ i (t) for the oscillating disturbance of five different lines of the PanTaGruEl model, shown in figure 2. Lines are perturbed one at a time, in different simulations. The right panels show the maximal amplitude of the frequenciesẋ i (t) for the same time series. The red circle indicate the indices of the two extremities of the disturbed line. Our method identifies unequivocally the endpoints of the disturbed lines whereas the time series of frequencies are not able to do so.
In the case of sparse (incomplete) measurements, we have tested, but not shown in this works, our method for faults where both ends of the perturbed line belong to the reduced (non-measured) nodes. In this case, our methods correctly detect the nodes in the vicinity of the faulty elements. However, at this task, our method is not significantly better than the detection proceeding by inspection of the trajectories of x i (t),ẋ i (t).

Euroroad and US airports
In order to illustrate the variety of couplings that our method covers, we applied it to: Figure 5. Phases, velocities, and frequency mismatch for the European road network [34,35] with linear coupling and two noisy lines. The pairs of trajectories blue-green and orange-red correspond to the end nodes of each disturbed lines. The gray area covers the same quantities for all other nodes. The disturbances cannot be seen in the phases and velocities trajectories, whereas they are clearly identified and located with the frequency mismatch. The panels on the right side are the snapshot at the end of the time series of the left side. Figure 6. Phases, velocities, and frequency mismatch for the US airports network [37,38] with third order Kuramoto coupling [36]. One line is disturbed by a step signal and the other by a ramp signal. The pairs of trajectories blue-green and orange-red correspond to the end nodes of each disturbed lines. The gray area covers the same quantities for all other nodes. The disturbances cannot be seen in the phases and velocities trajectories, whereas they are clearly identified and located with the frequency mismatch. The panels on the right side are the snapshot at the end of the time series of the left side.
• A set of n = 1039 agents linearly coupled according to the largest component of the European road network from [34,35]. Two arbitrary lines are perturbed by a noisy signal. Results are shown in figure 5; • A set of n = 1572 agents interacting through a higher-order Kuramoto coupling [36] (couplings of order up to q = 3 are used here), with interaction structure given by the largest connected component of the network of US airports from [37,38]. One line is perturbed by a ramp signal and another one is subject to a step signal. Results are shown in figure 6. In each case, we added a reasonable amount of noise at each node, for sake of generality.
In both figures 5 and 6, the left panels show the time series of x i (t),ẋ i (t), and ψ i (t) from top to bottom. The right panels show the snapshot of the same quantities at the end of the time series on the left, at t = 3.0.
We acknowledge that, even though it is straightforward to identify the end-nodes of the disturbed lines by visual inspection of figures 5 and 6, it might be an algorithmic challenge to automatize this process.

Conclusion
We proposed an elegant method to identify and locate disturbances in networked dynamical systems. Our method relies on time series of the agents degree of freedom and is able to differentiate between line disturbances and nodal disturbances. In the case of partial measurements, we are able to locate precisely the perturbation if the nodes where the it occurs are measured. Otherwise, we can determine the area of the network where the faulty, unmeasured nodes/lines are located, using the Kron-reduced network. The main condition required for our method to work out well is that the disturbance changes much slower than the intrinsic time scales of the system. However it is not restricted to a particular shape for the perturbation or a specific dynamical model, making it rather widely applicable. Interestingly, our method is also able to locate multiple disturbance occurring at the same time, given that they have different characteristics. This represents a remarkable feature as in large scale networks made of thousands or millions of nodes, many disturbances are likely to overlap in the time series.
We believe that the idea raised in this manuscript will help the development of more efficient tools to detect and locate disturbances more accurately. Future work will aim at proposing a scheme for online detection and localization of disturbances.