Topology identification of heterogeneous networks

This paper addresses the problem of identifying the graph structure of a dynamical network using measured input/output data. This problem is known as topology identification and has received considerable attention in recent literature. Most existing literature focuses on topology identification for networks with node dynamics modeled by single integrators or single-input single-output (SISO) systems. The goal of the current paper is to identify the topology of a more general class of heterogeneous networks, in which the dynamics of the nodes are modeled by general (possibly distinct) linear systems. Our two main contributions are the following. First, we establish conditions for topological identifiability, i.e., conditions under which the network topology can be uniquely reconstructed from measured data. We also specialize our results to homogeneous networks of SISO systems and we will see that such networks have quite particular identifiability properties. Secondly, we develop a topology identification method that reconstructs the network topology from input/output data. The solution of a generalized Sylvester equation will play an important role in our identification scheme. © 2020TheAuthors.PublishedbyElsevierLtd.ThisisanopenaccessarticleundertheCCBYlicense


Introduction
Graph structure plays an important role in the overall behavior of dynamical networks.Indeed, it is well-known that the convergence rate of consensus algorithms depends on the connectivity of the network topology.In addition, many properties of dynamical networks, like controllability, can be assessed on the basis of the network graph (Chapman & Mesbahi, 2013;Jia et al., 2019;Liu et al., 2011).Unfortunately, the graph structure of dynamical networks is often unknown.This problem is particularly apparent in biology, for example in neural networks and genetic networks (Julius et al., 2009), but also emerges in other areas such as power grids (Cavraro & Kekatos, 2018).
To deal with this problem, several topology identification methods have been developed.Such methods aim at reconstructing the topology (and weights) of a dynamical network on the basis of measured data obtained from the network.et al. ( 2019), Ramaswamy et al. (2018), Van den Hof et al. (2013), van Waarde et al. (2018), van Waarde et al. (2018), along with the joint topology and dynamics recovery problem (Ioannidis et al., 2019;Wai et al., 2019).
The goal of this paper is to provide a comprehensive treatment of topology identification for linear MIMO heterogeneous networks, with no assumptions on the network structure such as sparsity or regularity.Most existing work on topology identification emphasizes the role of the network topology by considering relatively simple node dynamics.For example, networks of single integrators have been studied in Hassan-Moghaddam et al. (2016), Morbidi and Kibangou (2014), Nabi-Abdolyousefi and Mesbahi (2010), van Waarde et al. (2019b).In addition, the papers (Suzuki et al., 2013) and Shahrampour and Preciado (2015) consider homogeneous networks comprised of identical single-input single-output systems.Nonetheless, there are many examples of networks in which the subsystems are not necessarily the same, for example, mass-spring-damper networks (Koerts et al., 2017), where the masses at the nodes can be distinct.Heterogeneity in the node dynamics has also been studied in the detail in synchronization problems, see e.g.Wieland et al. (2011), Yang et al. (2014).
We study topology identification for the general class of heterogeneous networks, where the node dynamics are modeled by general, possibly distinct, MIMO linear systems.We divide our analysis in two parts, namely the study of identifiability and the development of identification algorithms.The study of identifiability of the network topology deals with the question whether there exists a data set from which the topology can be uniquely identified.Identifiability of the topology is hence a property of the node systems and the network graph, and is independent of any data.Topological identifiability is an important property.Indeed, if it is not satisfied, then it is impossible to uniquely identify the network topology, regardless of the amount and richness of the data.After studying topological identifiability, we will turn our attention towards identification algorithms.Our two main contributions are hence the following: (1) We provide conditions for topological identifiability of general heterogeneous networks.Our results recover an identifiability result for the special case of networks of single integrators (Paré et al., 2013;van Waarde et al., 2019b).We will also see that homogeneous networks of singleinput single-output systems have quite special identifiability properties that do not extend to the general case of heterogeneous networks.(2) We establish a topology identification scheme for heterogeneous networks.The idea of the method is to reconstruct the interconnection matrix of the network by solving a generalized Sylvester equation involving the Markov parameters of the network.We prove that the network topology can be uniquely reconstructed in this way, under the assumptions of topological identifiability and persistency of excitation of the input data.
A preliminary version of our work was presented in van Waarde et al. (2019a).The contributions of the current paper are significant in comparison to van Waarde et al. (2019a) for two reasons.First, the identifiability results presented here are more general as they are applicable in situations when not all network nodes are excited.Also, the necessary conditions for identifiability of single-integrator networks are shown to carry over to the more general class of homogeneous networks of single-input singleoutput systems.Secondly, the topology identification approach is new, and attractive in comparison to van Waarde et al. (2019a) since the network interconnection matrix is computed directly and without the use of auxiliary variables.Our approach is also suitable for ''parallelization'' in the sense that each row block of the interconnection matrix can be computed independently.
The paper is organized as follows.In Section 2 we formulate the problem.Section 3 contains our results on topological identifiability.Subsequently, we describe our topology identification method in Section 4. Finally, we state our conclusions in Section 5.

Notation
We denote the Kronecker product by ⊗.The direct sum of matrices A 1 , A 2 , . . ., A k is the block diagonal matrix defined by Moreover, the concatenation of matrices Finally, let A(z) be an n × m rational matrix.Then the constant

Problem formulation
We consider a network model similar to the one studied by Fuhrmann and Helmke (Fuhrmann & Helmke, 2015, Ch. 9).Specifically, we consider networks composed of N discrete-time systems of the form where x i (t) ∈ R n i is the state of the ith node system, v i (t) ∈ R m i is its input and w i (t) ∈ R p i is its output for i = 1, 2, . . ., N.
The real matrices A i , B i and C i are of appropriate dimensions.We occasionally use the shorthand notation (A i , B i , C i ) to denote (1).
The coupling between nodes is realized by the inputs v i (t), which are specified as where u(t) ∈ R m is the external network input and Q ij and R i are real matrices of appropriate dimensions.In addition, let S i be a real p × p i matrix and consider the external network output y(t) ∈ R p , defined by Then, by introducing the block diagonal matrices and the matrices we can represent the network dynamics compactly as (3) Here x(t) = col(x 1 (t), x 2 (t), . . ., x N (t)) ∈ R n where n is defined as n := ∑ N i=1 n i .We emphasize that the coupling of the node dynamics is induced by the matrix Q , which we will hence call the interconnection matrix.
There are a few important special cases of node dynamics (1) and resulting network dynamics (3).If A i = A 0 , B i = B 0 and C i = C 0 for all i = 1, 2, . . ., N, the dynamics of all nodes in the network are the same and the resulting dynamical network is called homogeneous.The more general setting in which the node dynamics are not necessarily the same is referred to as a heterogeneous network.Another special case of node dynamics occurs when m i = p i = 1 for all i = 1, 2, . . ., N. In this case, the node systems are single-input single-output (SISO) systems, and the resulting dynamical network is referred to as a SISO network. 1 Topology identification of homogeneous SISO networks has been studied in Suzuki et al. (2013) and Shahrampour and Preciado (2015).In addition, topology identification has been well-studied (see e.g (Gonçalves & Warnick, 2008;Hassan-Moghaddam et al., 2016;Nabi-Abdolyousefi & Mesbahi, 2010;van Waarde et al., 2019b)) for networks of so-called single-integrators, in which the node dynamics are described by ẋi (t) = v i (t).This type of node dynamics can be seen continuous-time counterpart of (1) where The purpose of this paper is to study topology identification for general, heterogeneous dynamical networks of the form (3).Although we focus on discrete-time systems, our results can be stated for continuous-time systems as well.In order to make the problem more precise, we first explain what we mean by the topology of (3).Let G = (V, E) be a weighted directed graph with We refer to G as the topology of the dynamical network (3).With this in mind, the problem of topology identification concerns finding G (equivalently, finding Q ) using measurements of the input u(t) and output y(t) of (3).We assume knowledge of the local node dynamics (i.e., the matrices A, B and C ) as well as the external input/output matrices R and S. 2  At this point, we may ask the following natural question: is it possible to uniquely reconstruct the topology of (3) from input/output data?To formalize and answer this question, we define the notion of topological identifiability.Let y u,x 0 ,Q (t) denote the output of (3) at time t, where the subscript emphasizes the dependence on the input u(•), the initial condition x 0 = x(0) and interconnection matrix Q .The following definition is inspired by Grewal and Glover (1976) and defines the notion of distinguishability of interconnection matrices.

Definition 1.
Let y u,x 0 ,Q (•) and y u,x 0 , Q (•) denote the output trajectories of two systems of the form (3) with interconnection matrices Q and Q and initial conditions x 0 and x0 , respectively.We say that Q and Q are indistinguishable if there exist initial conditions x 0 , x0 ∈ R n such that for all input functions u.Moreover, Q and Q are said to be distinguishable if they are not indistinguishable.
With this in mind, the topology of (3) is said to be identifiable if Q is distinguishable from all other interconnection matrices.More formally, we have the following definition.
1 Here we emphasize that 'SISO' refers to the node systems of the network.
The overall network dynamics (3) can still have multiple external inputs and outputs.
2 This assumption is standard in the literature on topology identification, see, e.g., Shahrampour and Preciado (2015) and Suzuki et al. (2013).Without knowledge of the node dynamics, topology identification becomes a full system identification problem.
Definition 2. Consider system (3) with interconnection matrix Q .The topology of system (3) is said to be identifiable if Q and The importance of topological identifiability lies in the fact that unique reconstruction of Q from input/output data is only possible if the topology of (3) is identifiable.Indeed, if this is not the case, there exists some Q ̸ = Q that is indistinguishable from Q , meaning that both Q and Q explain any input/output trajectory of (3).Topological identifiability is hence a structural property of the system (3) that is independent of a particular data sequence and that is necessary for the unique reconstruction of Q from data.
Following Grewal and Glover (1976), it is straightforward to characterize topological identifiability in terms of the transfer matrix from u to y.This transfer function will be denoted by The topology of the networked system (3) is identifiable if and only if the following implication holds: Although Proposition 1 provides a necessary and sufficient condition for topological identifiability, the condition involves the arbitrary matrix Q .Hence, it is not clear how to verify the condition of Proposition 1. Instead, in this paper we want to establish conditions for topological identifiability in terms of the local system matrices A, B and C and the matrices Q , R and S.This is formalized in the following problem.
Problem 1. Find necessary and sufficient conditions on the node dynamics A, B, C , the external input/output matrices R, S and the interconnection matrix Q under which the topology of (3) is identifiable.
Our second goal is to identify Q from input/output data.
Problem 2. Develop a methodology to identify the interconnection matrix Q from measurements of the input u(t) and output y(t) of system (3).

Conditions for topological identifiability
In this section we state our solution to Problem 1 by providing necessary and sufficient conditions for topological identifiability.We start by providing an overview of the results that are proven in this section.In the following table, ''N'' denotes necessary and ''S'' denotes sufficient.
Theorem 2 General N-S conditions Theorem 3 N condition; also S if R has full rank Theorem 5 N condition for homogeneous SISO networks Theorem 7 N-S conditions for homog.SISO networks For analysis purposes, we first rewrite the network transfer matrix F Q (z).Note that Premultiplication by (zI − A) −1 and postmultiplication by the matrix (zI This means that where the transfer matrices of all node systems.Finally, by rearranging terms we obtain (5) Note that the inverse of I − G(z)Q exists as a rational matrix.Indeed, since (zI − A) −1 is strictly proper we see that lim z→∞ (I − G(z)Q ) = I.Therefore, we conclude by ( 5) that the transfer matrix We remark that ( 6) is an attractive representation of the network transfer matrix, since the matrices A, B and C describing the local system dynamics are grouped and contained in the transfer matrix G(z).
Remark 1.By ( 6), we see that the networked system (3) can be represented by the block diagram in Fig. 1.Hence, the problem of topology identification can be viewed as the identification of the static output feedback gain Q , assuming knowledge of the system G(z) and the external input/output matrices R and S.
The following theorem gives necessary and sufficient conditions for topological identifiability.We will use the notation Theorem 2. Consider the networked system (3) and assume that the matrix S has full column rank.The topology of (3) is identifiable if and only if cker By hypothesis, S has full column rank and hence We define ∆ := Q − Q .Then, ( 8) is equivalent to each of the following statements: Equivalently, Next, let vec(M) denote the vectorization of a matrix M. Then ( 9) is equivalent to (10) By ( 10) it is clear that the topology of ( 3) is identifiable if and only Finally, by the block diagonal structure of G(z), this is equivalent to (7) which proves the theorem.□ By Theorem 2, topological identifiability is equivalent to the matrices G i (z) ⊗ H ⊤ Q (z) having zero constant kernel.Note that this condition generally depends on the -a priori unknown -matrix Q .Notably, identifiability is independent of the particular matrix Q whenever all node inputs are excited and all node outputs are measured, as stated in the following theorem.
Theorem 3. Consider the networked system (3).If the topology of for all i, j ∈ V.In addition, suppose that S has full column rank and R has full row rank.Then the topology of (3) is identifiable if and only if (11) holds.
The importance of Theorem 3 lies in the fact that the identifiability condition (11) can be verified without knowledge of Q .This means that, whenever the rank conditions on S and R hold, one can check for topological identifiability before collecting data from the system.Remark 2. A proper transfer matrix T (z) has constant kernel {0} if and only if the matrix col(M 0 , M 1 , . . ., M r ) has full column rank.Here M 0 , M 1 , . . ., M r are the Markov parameters of T (z) and r is greater than or equal to the order of T (z).As such, the conditions of Theorems 2 and 3 can be verified by computing the rank of the Markov parameter matrices associated to the transfer matrices in (7) and (11).
Proof.We first prove the second statement.Suppose that S has full column rank and R has full row rank.Then We define Exploiting the block diagonal structure of G(z), we conclude that the topology of ( 3) is identifiable if and only if (11) holds.□ A consequence of Theorem 3 is that identifiability of the topology of (3) implies that the constant kernel of both G ⊤ i (z) and Based on this fact, we relate topological identifiability and output controllability of the node systems.
Definition 3. Consider the system where x ∈ R n , u ∈ R m and y ∈ R p , and let y u,x 0 (•) denote the output trajectory of (12) for a given initial condition x 0 and input u(•).System ( 12) is called output controllable if for every x 0 ∈ R n and y 1 ∈ R p there exists an input u(•) and time instant T ∈ N such that y x 0 ,u (T ) = y 1 .
Corollary 4. If the topology of (3) is identifiable then the systems Proof.By Theorem 3, identifiability of the topology of (3) implies that the constant kernel of The latter implication holds if and only if the output controllabil- is output controllable (Trentelman et al., 2001, Ex. 3.22).The proof for the necessity of output controllability of (A as an 'excitability' condition.Indeed, it guarantees that we have enough freedom in steering the output w i (t) of each node i ∈ V.
Example 1.We will now illustrate Theorems 2 and 3. Consider a network of N = 10 oscillators of the form Here mod denotes the modulo operation and ≡ denotes congruence.The network nodes are diffusively coupled, and an external input is applied to node 1, that is, otherwise, where This means that the interconnection matrix Q is defined element-wise as Since we only externally influence the first node system, the corresponding matrix R is given by the first column of I. We assume that we externally measure all node outputs, meaning that S = I.
Using Theorem 2, we want to show that the topology of (3) is identifiable.First, note that the transfer function G i (z) of node system i is given by implies that the topology of (3) is identifiable if and only if cker . This is equivalent to the output controllability of the system (A + BQC , BR, C ).It can be easily verified that the output controllability matrix has full row rank.We therefore conclude by Theorem 2 that the topology of (3) is identifiable.Note that the rank of the output controllability matrix (and hence, identifiability) depends on the interconnection matrix Q .
Next, we discuss the scenario in which R = I.In this case, we can externally influence all nodes.Now, identifiability can be checked without knowledge of Q .In fact, by Theorem 3, the topology of ( 3) is identifiable if and only if cker . This condition is satisfied, since all local transfer functions are nonzero scalars.
So far, we have provided a general condition for identifiability in Theorem 2, and we have discussed some of the implications of this result in Theorem 3 and Corollary 4.However, possible criticism of the results may arise from the full rank condition on S in Theorem 2, which, until now, has been left rather unjustified.
It turns out that full column rank of S (or the dual, full row rank of R) is necessary for topological identifiability in case the networked system is homogeneous and SISO.For this important class of networked systems, the rank condition on S in Theorem 2 is hence not restrictive.

Theorem 5.
Consider a homogeneous SISO network, that is, a system of the form (3) Before proving Theorem 5, we state the following lemma.
Lemma 6. Suppose that m i = p i = 1 and A is controllable and (S, Q ) is observable.
Proof.Suppose on the contrary that (S, Q ) is unobservable.Let v ∈ R N be a nonzero vector in the unobservable subspace of (S, Q ), i.e., This implies that SQ k = S(Q + vv ⊤ ) k for all k ∈ N. By (6), the network transfer matrix is given by where G 0 (z) := C 0 (zI − A 0 ) −1 B 0 is a scalar transfer function.Next, by expanding F Q (z) as a formal series Hence, the topology of (3) is not identifiable.The proof for necessity of controllability of (Q , R) is analogous and therefore omitted.□ Proof of Theorem 5. Suppose on the contrary that rank R < N and rank S < N. Then there exist nonzero vectors v 1 , v 2 ∈ R N such that Sv 1 = 0 and v ⊤ 2 R = 0. We assume without loss of generality that v 2 is such that v ⊤ 2 v 1 ̸ = −1.Next, we define By our assumption on v 2 , the matrix T is hence invertible, and .
We define the matrix Now, we distinguish two cases: we obtain , where T := T ⊗ I.Here we have used the fact that p i = m i = 1 for all i ∈ V, as well as the property (X Multiply from right by v 2 and rearrange terms to obtain Trentelman et al. (2001, Ch. 3)).By the previous lemma, this implies that the topology of (3) is not identifiable.□ Theorem 5 is interesting because it shows that the ability to measure all node outputs or to excite all node inputs is necessary for identifiability in the case of homogeneous SISO networks.This result allows us to sharpen Theorem 2 for this particular class of networks.

Theorem 7.
Consider a homogeneous SISO network, that is, a system of the form (3) with m i = p i = 1 and A i = A 0 , B i = B 0 and C i = C 0 for all i ∈ V.The topology of (3) is identifiable if and only if G 0 (z) := C 0 (zI −A 0 ) −1 B 0 ̸ = 0 and at least one of the following two conditions holds: Proof.To prove the 'if'-statement, we first assume that G 0 (z) is nonzero, rank S = N and (Q , R) is controllable.By Theorem 2, the topology of ( 3) is identifiable if and only if cker We expand the latter matrix as a formal series as We claim that by strict properness of G 0 (z), the powers G k 0 (z) (k = 0, 1, 2, . . . ) are linearly independent over the reals.Indeed, suppose α q 0 (z) where p 0 and q 0 are polynomials.
By strict properness of G 0 (z), this is a contradiction since every numerator on the right hand side of (15) has degree less than p k 1 0 (z)q kr −k 1 0 (z).Thus α 1 = 0.In fact, we can repeat the same argument to show α 1 = • • • = α r = 0, proving the claim of independence.It follows from ( 14 where we leveraged the hypothesis that G 0 (z) is nonzero.Now, using the fact that G k 0 (z) (k = 0, 1, 2, . . . ) are linearly independent, we obtain v ⊤ Q k R = 0 for all k ∈ N. We conclude by controllability of the pair (Q , R) that v = 0, hence cker H ⊤ Q (z) = {0}.In other words, the topology of ( 3) is identifiable.The sufficiency of the three conditions G 0 (z) ̸ = 0, rank R = N and (S, Q ) is observable is proven in a similar fashion and thus omitted.
To prove the 'only if'-statement, suppose that the topology of Furthermore, by Theorem 5, either S or R has full rank.□ It is noteworthy that full rank of either R or S is not necessary for topological identifiability of heterogeneous networks, as demonstrated next.
In addition, assume that R = ( ) ⊤ and S = ( ) . It can be easily verified that where Q 11 , Q 12 , Q 21 and Q 22 are the entries of the interconnection matrix ) .
We assume that for some interconnection matrix Q .By comparing the numerators of F Q and F Q we see that Q 21 = Q21 .Moreover, by comparing the coefficients corresponding to z 2 and z in the denominator, we obtain Q 11 = Q11 and Q 22 = Q22 .Finally, by comparing constant terms in the denominator, we see that . Hence, Q = Q and we conclude that the topology of ( 3) is identifiable.However, S does not have full column rank and R does not have full row rank.

Topology identification approach
In this section, we focus on the problem of topology identification, as formulated in Problem 2. The proposed solution consists of two steps: first identify the Markov parameters of the networked system (3), and then extract the matrix Q .There are several ways of computing the Markov parameters on the basis of input/output data, we will summarize some of them in the next section.

Identification of Markov parameters
Consider a general linear system of the form where x ∈ R n is the state, u ∈ R m is the input and y ∈ R p the output.In this section we recap how one can identify the Markov parameters D, CB, CAB, . . ., CA r B for r ∈ N, using measurements of the input and output of ( 16)-( 17).For a given signal f (t) with t = 0, . . ., T − 1, we define the Hankel matrix of depth k as The signal f (0), f (1), . . ., f (T − 1) is said to be persistently exciting of order k if H k (f ) has full row rank.Now suppose that we measure T samples of the input u(t) and output y(t) of ( 16)-( 17) for t = 0, 1, . . ., T − 1.We rearrange these measurements in Hankel matrices of depth n + r + 1.Moreover, we partition ) , where U p and Y p contain the first n row blocks of H n+r+1 (u) and H n+r+1 (y), respectively.The following result from Markovsky and Rapisarda (2008, Prop. 4) shows how the Markov parameters can be obtained from data.
Moreover, the Markov parameters can be obtained as Y f G = col(D, CB, CAB, . . ., CA r B).
Theorem 8 shows how the Markov parameters of the system can be obtained from measured input/output data.The input should be designed in such a way that it is persistently exciting, special cases of such inputs have been discussed in Verhaegen and Dewilde (1992).For u(0), . . ., u(T − 1) to be persistently exciting of order 2n + r + 1 a number of samples T ≥ (m + 1)(2n + r + 1) − 1 is necessary.In fact, there are input functions that achieve persistency of excitation of this order exactly for T = (m + 1)(2n + r + 1) − 1.A refinement of Theorem 8 is possible using the notion of weaving trajectories (Markovsky et al., 2005), which reduces the order of excitation to 2n+1.More generally, one can extend the notion of persistency of excitation to an arbitrary concatenation of multiple trajectories (van Waarde et al., 2020).This is useful in situations where single experiments are individually not sufficiently informative.
Remark 5.In addition to the deterministic setting of Theorem 8, there are approaches to identify the Markov parameters of systems with disturbances, i.e., systems of the form where v and w are zero mean, white vector sequences.In particular, the paper (Oymak & Ozay, 2018) studies the identification of the system's Markov parameters from finite data, and provides statistical guarantees for the quality of estimation.

Topology identification
Subsequently, we will turn to the problem of identifying the topology of (3) from the network's Markov parameters.As in Theorem 2, we will assume that S has full column rank.In fact, to lighten the notation, we will simply assume S = I, even though all results can be stated for general matrices S having full column rank.Under the latter assumption, the Markov parameters of (3) are given by Whenever the dependence of M ℓ (Q ) on Q is clear, we simply write M ℓ .It is not immediately clear how to obtain Q from the Markov parameters since M ℓ depends on the ℓth power of A + BQC .The following lemma will be helpful since it implies that M ℓ can essentially be viewed as an affine function in Q and lower order Markov parameters.
Lemma 9. We have that Proof.First, we claim that for square matrices D 1 and D 2 of the same dimensions, we have for all ℓ = 1, 2, . . . .It is straightforward to prove this claim by induction.Indeed, for ℓ = 1, (18) holds.If (18) holds for ℓ ≥ 1 then proving the claim.Subsequently, by substitution of D 1 = A and D 2 = BQC into (18), we obtain Finally, the lemma follows by pre-and postmultiplication by C and BR, respectively.□ Using Lemma 9, we can come up with a system of linear equations in the unknown interconnection matrix Q .To see this, let us denote K ℓ := M ℓ − CA ℓ BR.Moreover, define the Toeplitz matrix L by where r ≥ 2n − 1.We apply Lemma 9 for ℓ = 1, . . ., r to obtain Next, let L i denote the (i + 1)th column block of L and define the matrix K := col(K 1 , K 2 , . . ., K r ).We can then write (19) in a more compact form as which reveals that Q is a solution to a generalized Sylvester equation.Topology identification thus boils down to (i) identifying the network's Markov parameters, (ii) constructing the matrices K , L i and M i for i = 0, . . ., r −1 and (iii) solving the Sylvester equation.
We summarize this procedure in the following theorem.
Let the matrices K and L i be as before.If the topology of (3) is identifiable then the interconnection matrix Q is the unique solution to the generalized Sylvester equation in the unknown Q.
Proof.Note that the interconnection matrix Q is a solution to (21) by construction.Suppose that Q is also a solution to (21).
We want to prove that Q = Q .Since Q and Q are both solutions to (21), we have for ℓ = 1, 2, . . ., r.Here we have written the dependence of M ℓ−i−1 on Q explicitly, to distinguish between Q and Q .By Lemma 9 we have Clearly, M 0 (Q ) = CBR = M 0 ( Q ).In fact, we claim that M k (Q ) = M k ( Q ) for all k = 0, 1, . . ., r. Suppose on the contrary that there exists an integer s such that 0 < s ≤ r and M s (Q ) ̸ = M s ( Q ).We assume without loss of generality that s is the smallest integer for which this is the case.Then M k (Q ) = M k ( Q ) for all k = 0, 1, . . ., s − 1.By combining ( 22) and ( 23) we obtain using (24).This is a contradiction and we conclude that Finally, as the topology of ( 3) is identifiable, we conclude that Q = Q .This completes the proof.□

Solving the generalized Sylvester equation
In the previous section, we saw that the generalized Sylvester equation ( 21) plays a central role in our topology identification approach.In this section, we discuss methods to solve this equation.One simple approach to the problem is to vectorize Q and write (21) as the system of linear equations in the unknown vec(Q) of dimension ) .
However, a drawback of this approach is that the dimension of vec(Q) is quadratic in the number of nodes N.This means that for large networks, solving ( 26) is costly from a computational point of view.
For the 'ordinary' Sylvester equation of the form there are well-known solution methods that avoid vectorization. 3The general idea is to transform the matrices L 0 and M 1 to a suitable form so that the Sylvester equation is easier to solve.A classic approach is the Bartels-Stewart method (Bartels & Stewart, 1972) that transforms L 0 and M 1 to real Schur form by means of two orthogonal similarity transformations.The resulting equivalent Sylvester equation is then simply solved by backward substitution.A Hessenberg-Schur variant of this algorithm was proposed in Golub et al. (1979).The approach was also extended to be able to deal with the more general equation using QZ-decompositions (Golub et al., 1979, Sec. 7).The problem with all of these transformation methods is that they rely on the fact that the Sylvester equation consists of exactly two Qdependent terms, i.e., r = 1.Therefore, it does not seem possible to extend such methods to solve generalized Sylvester equations of the form (21) for r > 1, see also the discussion in Van Loan (2000, Sec. 2).Nonetheless, we can improve upon the basic approach of vectorization ( 26) by noting that the matrices A, B and C have a special structure.Indeed, recall from (2) that these matrices are block diagonal.This allows us to write down a Sylvester equation for each row block of Q.Let Q (j) denote the jth block row of Q for j ∈ V. Then it is straightforward to show that (21) is equivalent to for all j ∈ V, where L (j) i is the (i + 1)th column block of the matrix L (j) , given by ℓ the jth row block of K ℓ .The importance of ( 27) lies in the fact that each row block of Q can be obtained independently, which significantly reduces the dimensions of the involved matrices.In fact, ( 27) is equivalent to the linear system of equations in the unknown vec ) . Note that the unknown is linear in the number of nodes, assuming that m j and p i are small in comparison to N.

Robustness analysis
In the case that the Markov parameters M 0 , M 1 , . . ., M r are identified exactly, we can reconstruct the topology by solving the generalized Sylvester equation ( 21), or equivalently, the system of linear Eqs. ( 26).Now suppose that our estimates of the Markov parameters are inexact, and we have access to where the real matrices ∆ ℓ represent the perturbations.Accordingly, we define Kℓ := Mℓ − CA ℓ BR = K ℓ + ∆ ℓ .Let ∆ := col(∆ 1 , ∆ 2 , . . ., ∆ r ).In this case it is natural to look for an approximate (least squares) solution vec( Q ) that solves An obvious question is how the solution Q is related to the true interconnection matrix Q .The following lemma provides a bound on the infinity norm of vec( Q ) − vec(Q ).In what follows, we will make use of the constant where X † denotes the Moore-Penrose inverse of X .
Lemma 11.Consider the network (3) with S = I and suppose that its topology be identifiable.Assume that the solution Q to (30) is unique.Then we have that ) . (31) Note that the bound (31) tends to zero as ∆ 0 , ∆ 1 , . . ., ∆ r tend to zero, so Q is a good approximation of Q for small perturbations.
An overestimate of (31) can be obtained if some prior knowledge is available.In particular, note that α is readily computable from the estimated Markov parameters (29).The first two norms in (31) can be upper bounded if a bound on ∥∆ i ∥ ∞ is given.Iden- tification error bounds on the Markov parameters are derived, e.g., in Oymak and Ozay (2018).Finally, to estimate ∥vec(Q )∥ ∞ one requires a bound on the largest network weight, i.e., an upper bound on the largest (in magnitude) entry of Q .The upper bound (31) is useful in the case that the nonzero weights of the network are lower bounded in magnitude by some known positive scalar γ , an assumption that is common in the literature on consensus networks, cf.LeBlanc et al. (2013, Sec. 3).Indeed, in this case we can exactly identify the graph structure G from noisy Markov parameters if since identified entries smaller than 1 2 γ are necessarily zero.We will further illustrate this point in Example 3.
Proof.We make use of the shorthand notation The hypothesis that Q is unique is equivalent to A E having full column rank.By using ( 26) and the relation Mi = M i + ∆ i , we get Finally, taking infinity norms yields the upper bound (31).This completes the proof.□ Example 3. Consider the networked system in Example 1.We consider the situation in which only the first node of the network is externally excited.We already know by the discussion in Example 1 that the topology of the system is identifiable.Here, our aim is to reconstruct the topology on the basis of the noisy Markov parameters (29), where r = 40.The perturbations are drawn randomly from a normal distribution using the Matlab command randn, and scaled such that   ∆ ⊤ i   ∞ ≤ 10 −5 for all i.Since ∆ i is a vector, this also implies that ∥∆ i ∥ ∞ ≤ 10 −5 .In this example, we assume that the weights of the network (i.e., the entries of Q ) have magnitudes between 1 2 and 1.
We identify the matrix Q by solving (30).To get an idea of the quality of estimation, we want to find a bound on (31).First, we compute α = 464.7040.By the assumptions on the perturbations and network weights, we obtain the bounds ∥∆∥ ∞ ≤ 10 −5 and ∥vec(Q where we have used (Lancaster & Farahat, 1972, Thm. 8 & p. 413) to bound the Kronecker product.Combining the previous bounds, we conclude that (31) is less than or equal to 0.1883.
Since 0.1883 ≤ 0.25 we can round all entries of Q that are less than 0.25 to zero, since the corresponding entries in Q are necessarily zero.The resulting zero/nonzero structure of Q can be captured by a graph Ĝ that we display in Fig. 2. Clearly, the structure of Ĝ is identical to the graph defined in Example 1, and the weights of Ĝ are close to the weights of G. Next, we repeat the experiment for larger perturbations, i.e., for ∥∆ i ∥ ∞ and   ∆ ⊤ i   ∞ bounded by 0.01.We identify Q and use the same rounding strategy as before to obtain a graph Ĝ in Fig. 3.Note that Ĝ resembles the original network structure G.In fact, all links are identified correctly, except for (7, 8) and the spurious link (4, 8).In this case, the bound (31) equals 49.9997, illustrating the fact that (31) can be conservative.

Conclusions
In this paper we have studied the problem of topology identification of heterogeneous networks of linear systems.First, we have provided necessary and sufficient conditions for topological identifiability.These conditions were stated in terms of the constant kernel of certain network-related transfer matrices.We have also seen that homogeneous SISO networks enjoy quite special identifiability properties that do not extend to the heterogeneous case.Subsequently, we have turned our attention to the topology identification problem.The idea of the identification approach was to solve a generalized Sylvester equation involving the network's Markov parameters to obtain the network topology.One of the attractive features of the approach is that the structure of the networked system can be exploited so that each row block of the interconnection matrix can be obtained individually.
The generalized Sylvester equation ( 21) plays an important role in our identification approach.Numerical solution methods are less well-developed for this equation than they are for the standard Sylvester equation (Bartels & Stewart, 1972;Golub et al., 1979).Hence, it would be of interest to further develop numerical methods for Sylvester equations of the form (21).We note that a Krylov subspace method has already been developed in Bouhamidi and Jbilou (2008).Another direction for future work is to study topological identifiability with prior information on the interconnection matrix.For example, from physical principles it may be known that Q is Laplacian.Such prior knowledge could be exploited to weaken the conditions for identifiability in Theorems 2, 3 and 7.

Remark 4 .
Theorem 5 generalizes several known results (seeParé et al. (2013), van Waarde et al. (2019a, 2019b)) for networks of single-integrators.Indeed, in the special case that A 0 = 0, B 0 = C 0 = 1, the node output w i (t) equals the node state x i (t) for all i ∈ V, and Theorem 5 asserts that either full state measurement or full state excitation is necessary for identifiability.This fact has been observed in different setups inParé et al. (2013, Thm.1), van  Waarde et al. (2019b, Rem. 2), and van Waarde et al. (2019a, Thm.  5).