Reduced Order Modeling of Diffusively Coupled Network Systems: An Optimal Edge Weighting Approach
Xiaodong Cheng, Lanlin Yu, Dingchao Ren, Jacquelien M.A. Scherpen
11 Reduced Order Modeling of Diffusively Coupled Network Systems:An Optimal Edge Weighting Approach
Xiaodong Cheng, Lanlin Yu, Dingchao Ren and Jacquelien M.A. Scherpen
Abstract —This paper studies reduced-order modeling of dy-namic networks with strongly connected topology. Given a graphclustering of an original complex network, we construct aquotient graph with less number of vertices, where the edgeweights are parameters to be determined. The model of thereduced network is thereby obtained with parameterized systemmatrices, and then an edge weighting procedure is devised,aiming to select an optimal set of edge weights that minimizes theapproximation error between the original and the reduced-ordernetwork models in terms of H -norm. The effectiveness of theproposed method is illustrated by a numerical example. I. I
NTRODUCTION
With the growing complexity of spatially interconnected dy-namic systems, the importance of understanding and managingdynamic networks has been widely recognized. An importantclass of dynamic networks is given by the so-called diffusivelycoupled networks, which are commonly used for describingdiffusion processes, e.g., information or energy spreading innetworks. The examples can be found in vehicle formations,electrical networks, synchronization in sensor networks, andopinion dynamics in social networks, see e.g. [1]–[4]. Thespatial structures of such systems are usually complex andresult in high-dimensional models that cause challenges foranalysis, control, and optimization. To effectively capture thecollective behaviors of dynamics over networks, it is desirableto simplify the structure of a complex network without asignificant loss of accuracy.Different from model reduction problems for other typesof dynamic systems, the one considered in this paper putsemphasis on the preservation of the network structure, whichis necessary for applications e.g., distributed controller designand sensor allocation [5]–[7]. Conventional model reductionmethods, e.g., balanced truncation and moment matching,merely focus on approximating the input-output behavior ofa given dynamic system [8], while the preservation of thenetwork structure is barely guaranteed. Although a generalized
This work of X. Cheng is supported by the European Research Council(ERC), Advanced Research Grant SYSDYNET, under the European UnionsHorizon 2020 research and innovation programme (Grant Agreement No.694504). This work of L. Yu is supported by the National Natural ScienceFoundation of China Under Project 61761136005.Xiaodong Cheng is with Control Systems Group, Department of ElectricalEngineering, Eindhoven University of Technology, 5600 MB Eindhoven, TheNetherlands [email protected]
Lanlin Yu is with the Institute of Advanced Technology, Westlake Institutefor Advanced Study, Westlake University, Hangzhou 310024, China andDingchao Ren is with the Department of Automation, University of Scienceand Technology of China, Hefei 230026, China. { yulanlin1992,dcrenustc } @gmail.com Jacquelien M.A. Scherpen is with Jan C. Willems Center for Systems andControl, Engineering and Technology Institute Groningen, Faculty of Scienceand Engineering, University of Groningen, Nijenborgh 4, 9747 AG Groningen,the Netherlands. [email protected] balanced truncation approach in [9] is able to construct anaccurate reduced-order model with a network interpretation,the relation between the original and obtained new typolo-gies is not yet clear. Singular perturbation approximation, asalternative network structure-preserving approach, has beenapplied to the reduction of electric circuits [10] and chemicalreaction networks [11]. This class of methods mainly relies ontime-scale separation of the states in an autonomous networksystem, while the external inputs are not considered explicitly.Besides, the resulting reduced topology hardly retains sparsity.Recently, clustering-based methods have been intensivelystudied and become the mainstream methodology for reducingnetwork systems, see e.g., [12]–[20]. With graph clustering,the vertices in a large-scale network are partitioned into severaldisjoint clusters. This class of methods has a clear advantage inretaining the consensus property [12], [19], system positivity[15], and scale-free property [20] in reduced-order models.The model reduction procedure can be implemented via thePetrov-Galerkin projection framework, where the projectionmatrix is formed based on the vertex clusters. However, all thecurrent clustering-based methods put their main focus on find-ing suitable clusters. After clusters are found, reduced-ordernetwork models are then directly determined by the projectionframework, while the freedom to construct a reduced-ordernetwork model with higher accuracy is overlooked.In this paper, we will explore the latter freedom and providea novel method for reduced-order modeling of directed net-works. We do not aim to find an optimal clustering. Instead, weassume that the clustering of a network is given, which leadsto a quotient graph . A parameterized reduced-order model isestablished based on this quotient graph, in which the edgeweights are free variables to be optimized. Then, the majorproblem in this paper follows: How to tune the edge weightsin the parameterized reduced-order model to minimize theapproximation error?This problem can be formulated as an optimization problemwith the objective to minimize the H -norm of the reductionerror between original and reduced network systems, in whichthe edge weights of the reduced network are variables to beoptimized. This edge weighting problem is subject to a bilinearmatrix inequality (BMI) constraint, which is computationallyexpensive. Therefore, we devise a novel edge weighting al-gorithm based on the convex-concave decomposition , whichlinearizes the nonconvex constraint as a convex one in theform of a linear matrix inequality (LMI). An iterative schemeis implemented to search for a set of optimal weights. Theconvergence of this algorithm is theoretically ensured, andthus at least a local optimum can be reached. Moreover,we initialize the edge weights as the outcome of clustering-based projection, such that the obtained reduced-order network a r X i v : . [ m a t h . O C ] M a r model is guaranteed a better approximation accuracy than theclustering-based projection methods.The rest of this paper is organized as follows. In Section II,we recap some preliminaries on graph theory and introduce theproblem setup. In Section III, the parameterized reduced-ordermodel is formulated, and an edge weighting algorithm is pro-posed to minimize the approximation error. In Section IV, theproposed method is illustrated by an example, and Section Vfinally makes some concluding remarks. Notation:
The symbol R and R + denote the set of realnumbers and positive real numbers, respectively. Let S n bethe set of real symmetric matrices of size n × n . I n is theidentity matrix of size n , and n represents the vector in R n of all ones. The cardinality of a set S is denoted by |S| . Fora real matrix A , the columns of A ⊥ form a basis of the nullspace of A , that is, AA ⊥ = 0 .II. P RELIMINARIES AND P ROBLEM S ETTING
This section provides necessary definitions and concepts ingraph theory used in this paper, and we refer to [21] for moredetails. The model of a dynamical network is then introducedand the model reduction problem is formulated.
A. Graph Theory
A directed graph G := ( V , E ) consists of a finite andnonempty node set V := { , , · · · , n } and an edge set E ⊆ V × V . Each element in E is an ordered pair of V , andif ( i, j ) ∈ E , we say that the edge is directed from vertex i tovertex j . A directed graph G is called simple , if G does notcontain self-loops (i.e., E does not contain any edge of theform ( i, i ) , ∀ i ∈ V ), and there exists only one edge directedfrom i to j , if ( i, j ) ∈ E .Next, we introduce several important matrices for charac-terizing a directed simple graph. Let m := |E| , the incidencematrix B ∈ R n × m is defined by B ij = +1 if edge j is directed from vertex i, − if edge j is directed to vertex i, otherwise . If each edge is assigned a positive value (weight), then the weighted adjacency matrix of G , denoted by A , is definedsuch that A ij ∈ R + denotes the weight of edge ( j, i ) ∈ E ,and A ij = 0 if ( j, i ) / ∈ E . In the case of a simple graph, A isa binary matrix with zeros on its diagonal. Then, the Laplacianmatrix L ∈ R n × n of the graph G is defined as L ij = (cid:26) (cid:80) nj =1 ,j (cid:54) = i A ij if i = j, −A ij otherwise. (1)Clearly, L = 0 . The diagonal entries of L are strictly positive,and the off-diagonal entries are non-positive. Alternatively,we can characterize the Laplacian matrix using the incidencematrix of G as L = B W B (cid:62) , (2)where B is a binary matrix obtained by replacing all “ − ”entries in the incidence matrix B with zeros, and W := Diag( w ) , with w = (cid:2) w w · · · w |E| (cid:3) (cid:62) , and w k the positive weight associated to the k -th edge, for all k ∈ { , , · · · , |E|} .For a vertex in a weighted graph, the indegree and outdegree of the vertex are computed as (cid:80) j ∈V A ij and (cid:80) i ∈V A ij ,respectively. A strongly connected graph G is called balanced if the indegree and outdegree of each vertex in G is equal.From (1), the following lemma is immediate. Lemma 1.
A weighted strongly connected graph G is balancedif and only if one of the following conditions hold. The edge weights of G satisfies B w = 0 . The Laplacian matrix of G satisfies ker( L ) =ker( L (cid:62) ) = span( ) . The strong connectivity implies that there is only one zeroeigenvalue of L [21], and the balance of G then indicates thatboth the row and column sums of L are zero. Remark 1.
Undirected graphs can be viewed as special bal-anced directed graphs with bidirectional edges. The Laplacianmatrix of an undirected graph is L = B W B (cid:62) , where B is an incidence matrix obtained by assigning an arbitraryorientation to each edge of the undirected graph, and W is apositive diagonal matrix representing edge weights. Next, we recapitulate the notion of graph clustering, whoseconcept can be found in e.g., [12]–[17].
Definition 1.
Let G := ( V , E ) be a directed graph. Then, agraph clustering is a partition of V into r nonempty disjointsubsets C , C , · · · , C r covering all the elements in V , where C i is called a cluster of G . Let {C , C , · · · , C r } be a clustering of G with n vertices.This graph clustering can be characterized by a binary char-acteristic matrix Π ∈ R n × r , whose rows and columns arecorresponding to the vertices and clusters, respectively: Π ij := (cid:40) if vertex i ∈ C j , otherwise . Remark 2.
Note that all the clusters are nonoverlapping, i.e.,each vertex can be not assigned to distinct clusters. Therefore,each row of the characteristic matrix Π only has one nonzeroelement. Specifically, we have Π r = n and (cid:62) n Π = [ |C | , |C | , · · · , |C r | ] . (3) B. Problem Setup
In this paper, we consider a network system evolving overa directed graph G , which is simple, weighted and stronglyconnected. The dynamics of each vertex is governed by ˙ x i ( t ) = − n (cid:88) j =1 A ij [ x i ( t ) − x j ( t )] + p (cid:88) j =1 f ik u k ( t ) , (4)where x i ( t ) ∈ R is the state of vertex i , and A ij is the ( i, j ) -thentry of the adjacency matrix of G , representing the strengthof the coupling between vertices i and j . u k ( t ) ∈ R is theexternal input, and f ij ∈ R is the gain of the j -th input actingon vertex i , which is zero if and only if u j has no effect onvertex i . Let F ∈ R n × p be the matrix such that F ij = f ij . We then present the dynamics of the overall network in a compactform as Σ : (cid:40) ˙ x ( t ) = − Lx ( t ) + F u ( t ) ,y ( t ) = Hx ( t ) , (5)with x ( t ) := [ x , x , · · · , x n ] (cid:62) ∈ R n and u :=[ u , u , · · · , u p ] (cid:62) ∈ R p . The vector y ∈ R q collects theoutputs of the network, and H is the output matrix.This paper aims for structure-preserving model reduction ofdiffusively coupled networks in form of (5), and the reduced-order model not only approximates the input-output mappingof the original network system with a certain accuracy but alsoinherits an interconnection structure with diffusive couplings.To this end, we adopt graph clustering to build up a reduced-order network model. Specifically, the problem addressed inthis paper is as follows. Problem 1.
Given a network system Σ as in (5) and a graphclustering {C , C , · · · , C r } , find a reduced-order model ˆΣ : (cid:40) ˙ˆ x = − ˆ Lx + ˆ F u ˆ y = ˆ H ˆ x (6) with ˆ x ∈ R r , r (cid:28) n , such that ˆ L is the Laplacian matrix of areduced directed graph, and the reduction error (cid:107) Σ − ˆΣ (cid:107) H is minimized. ˆ L ∈ R r × r , ˆ F ∈ R r × p , ˆ H ∈ R q × r are matricesdepending on the graph clustering. It is worth emphasizing that Problem 1 does not aim to findan appropriate graph clustering of the network G . Instead, wefocus on how to establish a “good” reduced-order model withgiven clusters. Thus, it is an essentially different problem frome.g., [14], [15], [17], [19], and we do not apply the Petrov-Galerkin projection framework.III. M AIN R ESULT
In this section, a novel model reduction approach fornetwork systems is presented with two steps. In the fist step,a parameterized model of a reduced network is constructed onthe basis of graph clustering. Then, the second step computes aset of parameters in an optimal fashion such that the H -normof approximation error is minimized. A. Parameterized Reduced-Order Network Model
Given a graph clustering of the original network, we presenta parameterized model for the reduced network, whose in-terconnection topology is determined by the clustering. Animportant property of this parameterized model is that itguarantees the boundedness of the reduction error (cid:107) Σ − ˆΣ (cid:107) H for all positive edge weights.To derive a parameterized reduced-order network modelwith such a property, we first convert the system (5) to its balanced graph representation as follows. Lemma 2.
If the underlying graph of Σ in (5) is stronglyconnected, then there exists a diagonal M ∈ R n × n withpositive diagonal entries such that Σ is equivalent to (cid:26) M ˙ x ( t ) = − L b x ( t ) + F b u ( t ) ,y ( t ) = Hx ( t ) , (7) where F b = M F and L b = M L is the Laplacian of abalanced graph.
The proof follows directly from [19]. Next, we establish areduced-order model using the representation (7) to guaranteea bounded reduction error (cid:107) Σ − ˆΣ (cid:107) H .Let G b be the balanced graph of G . Note that G and G b have the same incidence matrix B . Given a graph clustering {C , C , · · · , C r } , the quotient graph ˆ G b is r -vertex directedgraph obtained by aggregating all the vertices in each clusteras a single vertex, while retaining connections between clustersand ignoring the edges within clusters. Specifically, if there isan edge ( i, j ) in G b with vertices i, j in the same cluster, thenit will not be presented as an edge in ˆ G b . If there exists anedge ( i, j ) with i ∈ C k and j ∈ C l , then there is an edge ( k, l ) in ˆ G b .Let ˆ B be the incidence matrix of the quotient graph ˆ G b .Algebraically, it can be verified that ˆ B is obtained by removingall the zero columns of Π (cid:62) B , where B is the incidence matrixof G b (or G ). Furthermore, we denote ˆ W = Diag( ˆ w ) , with ˆ w = (cid:2) ˆ w ˆ w · · · ˆ w m (cid:3) (cid:62) , (8)as the edge weight matrix of G b , where ˆ w k ∈ R + , and m isnumber of edges in ˆ G b . In order to maintain ˆ G b as a balancedgraph, we impose the constraint on its edge weights as ˆ B ˆ w = 0 , (9)according to Lemma 3. Thereby, the dynamics on the balancedquotient graph ˆ G b is then obtained as (cid:26) ˆ M ˙ x ( t ) = − ˆ L b ( ˆ W ) x ( t ) + ˆ F b u ( t ) ,y ( t ) = ˆ Hx ( t ) , (10)with the reduced matrices ˆ M = Π (cid:62) M Π , ˆ L b ( ˆ W ) = ˆ B ˆ W ˆ B (cid:62) , ˆ F b = Π (cid:62) F b , and ˆ H = H Π , (11)where ˆ B is the binary matrix obtained by replacing all the“ − ” entries with zeros in ˆ B , and ˆ M b ∈ R r × r , ˆ F b ∈ R r × p ,and ˆ H ∈ R q × r are reduced matrices determined by the givenclustering of G b . Since the graph clustering is given, i.e., Π isknown, the only parameters to be decided are the weights in ˆ W , which satisfy the constraint (9).From the reduced graph balanced representation (10), weimmediately construct a parameterized reduced-order modelin the form of (6) with the reduced matrices ˆ L ( ˆ W ) = ˆ M − ˆ B ˆ W ˆ B (cid:62) , ˆ F = ˆ M − Π (cid:62) F b , ˆ H = H Π , (12)where ˆ L represents a reduced weighted graph ˆ G . In (12), onlythe weight matrix ˆ W is to be determined, which is selectedfrom the following set M := { W = Diag( ˆ w ) | ˆ w ∈ R m + , ˆ B ˆ w = 0 } . (13)In the following example, we demonstrate the parameterizedmodeling of a simplified dynamic network. Example 1.
Consider an network example in vehicle forma-tion [1], where the formation topology G is depicted in Fig. 1a.Clearly, G is balanced, i.e., G = G b , with the incidence matrix B = − − − − − −
10 0 0 0 − − − − . Suppose that each vehicle is modeled as a first-order integratorwhich has the identical mass, i.e., M = I . An external control u is applied on vertex 4, and the vertex 1 is measured as theoutput signal y . Then, the network model is obtained in theform of (7) with L b = − − − − − − − − − − , F b = , and H = [ ] . Consider a clustering of G b as C = { , } , C = { , , } , C = { } , which leads to the characterization matrixas Π = (cid:62) The topology of the quotient graph ˆ G b is shown in Fig. 1b withthe incidence matrix ˆ B = − − − − . All the edge weights of ˆ G b are positive parameters to be deter-mined, as labeled in Fig. 1b, which leads to the parameterizedLaplacian matrix as ˆ L b ( ˆ W ) = ˆ w − ˆ w − ˆ w ˆ w − ˆ w − ˆ w ˆ w + ˆ w . The weights satisfy the constraint ˆ B ˆ w = 0 , namely, ˆ w = ˆ w − ˆ w , ˆ w = ˆ w , such that ˆ G b is balanced. The other matrices inthe reduced-order model (10) are computed as ˆ M = Π (cid:62) M Π = , ˆ F b = Π (cid:62) F b = , and ˆ H = H Π = (cid:2) (cid:3) . Then, in the parameterized reduced-order model 6, we have ˆ L ( ˆ W ) = ˆ M − ˆ L b ( ˆ W ) = ˆ w − ˆ w − ˆ w ˆ w − ˆ w − ˆ w ˆ w + ˆ w , with ˆ w = ˆ w − ˆ w , and ˆ w = ˆ w . The corresponding reducedgraph is depicted in Fig. 1c, which is no longer balanced. Remark 3.
The physical interpretation of the reduced matri-ces in (11) are explained. ˆ M is constructed such that the massof vertex k in ˆ G b is equal to the mass sum of all the verticesin C k in G b . The expression of ˆ F b means that if a vertex in a
54 3 6 12 u y (a)
2' 3'1' u ^ y ^ w ^ w ^ w ^ w (b)
2' 3'1' u ^ y ^ w ^ w ^ w ^ w (c) Fig. 1: (a) A directed balanced network G b consisting of 6 ver-tices, in which vertex 3 is controlled and vertex 4 is measured.The vertices in different clusters are indicated by distinctcolors. (b) The quotient graph ˆ G b consisting of 3 vertices,where the edge weights are parameters to be determined. Thequotient graph is balanced, when the constraints ˆ w = ˆ w − ˆ w and ˆ w = ˆ w are imposed. (c) The resulting reduced graph ˆ G . cluster C k of G b is controlled by an external input, then vertex k in ˆ G b is also controlled. Analogously, ˆ H indicates that avertex k in ˆ G b is measured if there is a measurement takenfrom a vertex in C k . With the reduced matrices in (12) and the constraint in (9),an important property of the reduced-order network model ˆΣ is that it guarantees the H reduction error between the originalsystem Σ in (5) and ˆΣ is always bounded.The computation of the reduction error amounts to find the H norm of the following error system: G e ( s ) = C e ( sI − A e ) − B e , (14)where A e = − (cid:20) L
00 ˆ L (cid:21) , B e = (cid:20) F ˆ F (cid:21) , C e = (cid:2) H − ˆ H (cid:3) . Note that A e is not Hurwitz, since L and ˆ L are both Laplacianmatrices containing zero eigenvalues. Thus, (cid:107) G e ( s ) (cid:107) H cannotbe calculated directly using the state space representation (14).Here, we employ the following matrices S n = (cid:20) − I n − (cid:62) n − (cid:21) ∈ R n × ( n − , S r = (cid:20) − I r − (cid:62) r − (cid:21) ∈ R r × ( r − (15)which are independent of system dynamics and satisfy S (cid:62) n n = 0 , and S (cid:62) r r = 0 . Let their left pseudo inverses be S + n : = ( S (cid:62) n M − S n ) − S (cid:62) n M − ∈ R ( n − × n ,S + r : = ( S (cid:62) r ˆ M − S r ) − S (cid:62) r ˆ M − ∈ R ( r − × r . Then, using the matrices in (15), we show the following result.
Lemma 3.
Consider the network system Σ in (5) and thereduced-order network model ˆΣ in (6) with matrices in (12) .Then, Σ − ˆ Σ ∈ H holds for all ˆ W ∈ M . Proof.
With S n and S r in (15), we construct a nonsingular ( n + r ) × ( n + r ) matrix as U e = (cid:20) σ − M n M − S n σ − M r M − S r (cid:21) , (16)where σ M := (cid:62) n M n = (cid:62) r Π (cid:62) M Π r = (cid:62) r ˆ M r . Theinverse of U e is given as U − e = (cid:62) n M (cid:62) r ˆ MS † n M S † r ˆ M . Note that L = M − L b and ˆ L = ˆ M − ˆ L b , where both L b and ˆ L b are the Laplacian matrices of balanced graphs, satisfying (cid:62) n L b = 0 , L b n = 0 , and (cid:62) r ˆ L b = 0 , ˆ L b r = 0 . Using theseproperties, we obtain G e ( s ) = C e U e ( sI − U − e A e U e ) − U − e B e = (cid:2) ¯ C e C e (cid:3) (cid:20) sI − A e ) − (cid:21) (cid:20) ¯ B e B e (cid:21) , = ¯ C e ¯ B e + C e ( sI − A e ) − B e , (17)where A e = − (cid:20) S + n L b M − S n S + r ˆ L b ˆ M − S r (cid:21) ,B e = (cid:20) S + n F b S + r ˆ F b (cid:21) , C e = (cid:2) HM − S n − ˆ H ˆ M − S r (cid:3) , (18)and ¯ B e := (cid:20) (cid:62) n F b (cid:62) r ˆ F b (cid:21) , and ¯ C e := σ − M (cid:2) H n − ˆ H r (cid:3) . (19)It follows from (3) that (cid:62) r ˆ F b = (cid:62) n F b , and ˆ H r = H n ,which yield ¯ C e ¯ B e = σ − M ( H n (cid:62) n F b − ˆ H r (cid:62) r ˆ F b ) = 0 . Thus,(17) becomes G e ( s ) = C e ( sI − A e ) − B e . (20)It is not hard to verify that both the matrices − S + n L b M − S n and − S + r ˆ L b ˆ M − S r are Hurwitz. Consequently, G e ( s ) in (20)is asymptotically stable, i.e., Σ − ˆ Σ ∈ H , for all ˆ W ∈ M .Next, we discuss the consensus property of the reduced-order network (6) with the matrices in (12). Consensus is atypical property of diffusively coupled networks, and it impliesthe nodal states converge to a common value in the absenceof the external input. More precisely, the network system in(5) reaches consensus if lim t →∞ [ x i ( t ) − x j ( t )] = 0 holds for all i, j ∈ V and all initial conditions. Proposition 1.
Consider the network system Σ in (5) whichreaches consensus. Then, the reduced-order model ˆΣ in (6) also reaches consensus, for any clustering Π and ˆ W ∈ M .Proof. It can verified that the parameterized Laplacian matrix ˆ L defined in (12) characterizes a strongly connected graph.Thus, ˆ L has only one zero eigenvalue. Then, the proof imme-diately follows from [19], [21]. The parameterized modeling of the reduced dynamic net-work using the graph balanced representation in (7) guaranteesthe stability of the error system (14), whose H norm canbe evaluated via the transfer function (20) with the Hurwitzmatrix A e . Note that in (20), the matrices S n and S r in (15) isonly dependent on the sizes of the networks, and Π is knownfor a given graph clustering, then the weights in ˆ W becomethe only unknown parameters to be determined in the follow-up procedure. In the following section, we aim for an optimalselection of the edges weights in the reduced network. B. Optimal Edge Weighting
In this section, we aim for an optimization scheme fordetermining ˆ W ∈ M that minimizes the approximation error (cid:107) G e ( s ) (cid:107) H . Thereby, the following problem is addressed. Problem 2.
Consider the original network system Σ in (5) .Given a graph clustering Π , find a ˆ W ∈ M such that (cid:107) G e ( s ) (cid:107) H is minimized, where ˆΣ is the reduced networkmodel defined in (6) with the matrices (12) . To solve this problem, we apply an optimization techniquebased on the convex-concave decomposition , which can beimplemented to search for a set of optimal weights itera-tively. A fundamental step toward the implementation is todevelop a necessary and sufficient condition for characteriz-ing (cid:107) G e ( s ) (cid:107) H , which leads to suitable constraints for theoptimization problem. Theorem 1.
Given the network system Σ in (5) . There exists areduced-order network model ˆΣ in (6) such that (cid:107) G e ( s ) (cid:107) H < ˆ γ if and only if there exist matrices ˆ Q = ˆ Q (cid:62) > withdimension ˆ Q ∈ R ( n + r − × ( n + r − , ˆ R = ˆ R (cid:62) > withdimension R ∈ R q × q , ˆ W ∈ M , and ˆ δ ∈ R + , such that thefollowing inequalities are satisfied, ˆ Q ¯ A + ¯ A (cid:62) ˆ Q ˆ QB e ˆ QEB (cid:62) e ˆ Q − ˆ δI E (cid:62) ˆ Q + − ¯ A (cid:62) r ¯ A r A (cid:62) r A r − I < , (21) (cid:20) ˆ Q ˆ δC (cid:62) e ˆ δC e ˆ R (cid:21) > , (22) tr( ˆ R ) < ˆ γ, (23) where B e , C e are defined in (18) , and ¯ A = (cid:20) − S + n L b M − S n
00 0 (cid:21) , E = (cid:20) I (cid:21) ¯ A r = (cid:20) − S + r ˆ B ˆ W ˆ B (cid:62) ˆ M − S r (cid:21) . (24) Proof.
Consider the error system G e ( s ) in (20), whichis asymptotically stable. Following e.g., [22], we have (cid:107) G e ( s ) (cid:107) H < γ , with γ ∈ R + , if and only if there existmatrices Q = Q (cid:62) > and R = R (cid:62) > such that (cid:20) QA e + A (cid:62) e Q QB e B (cid:62) e Q − I p (cid:21) < , (25) (cid:20) Q C (cid:62) e C e R (cid:21) > , (26) tr( R ) < γ, (27) where A e , B e , C e are defined in (18).In the following, we prove that the three inequalities areequivalent to (21), (22), and (23), respectively. First, it is nothard to verify that (25) is equivalent to QA e + A (cid:62) e Q QB e QEB (cid:62) e Q − I E (cid:62) Q − δI < (28)for a sufficiently large scalar δ ∈ R + , where E is defined in(24). Consider a nonsingular matrix T = I I − ¯ A r I . Pre- and post-multiplying by T (cid:62) and T , respectively, (28) thenbecomes (21), where the equation A e = ¯ A + E ¯ A r , and thesubstitutions ˆ δ = δ > , ˆ Q = δ Q > are used.Next, we observe that the following implications hold. (cid:20) Q C (cid:62) e C e R (cid:21) > ⇔ (cid:20) δ ˆ Q C (cid:62) e C e R (cid:21) > ⇔ (cid:20) ˆ Q ˆ δC (cid:62) e ˆ δC e ˆ R (cid:21) > R ) < γ ⇔ δ tr( ˆ R ) < γ ⇔ tr( ˆ R ) < ˆ γ, with ˆ R = ˆ δR and ˆ γ = ˆ δγ . As a result, (22) and (23) areequivalent to (26) and (27), respectively.Based on Theorem 1, we reformulate Problem 2 moreexplicitly as the following minimization problem min ˆ Q> , ˆ W ∈M tr( R ) (29)s.t. (21) and (22) hold , where ˆ R = ˆ δR with a given ˆ δ ∈ R + . Note that the constraint(22) can be solved efficiently using standard LMI solvers,while (21), due to the nonlinearity term ¯ A (cid:62) r ¯ A r , is a bilinearmatrix inequality, which causes the major challenge in solvingthe problem (29).To handle the bilinear constraint (21), we adopt the tech-nique called psd-convex-concave decomposition [23]. Definition 2.
A matrix-valued mapping
Φ : R n → S (cid:96) is calledpositive semidefinite convex concave (psd-convex-concave) if Φ can be expressed as Φ = Φ − Φ , where Φ k , with k = 1 , ,are positive semidefinite convex (psd-convex), i.e., Φ k ( λw + (1 − λ ) w ) ≤ λ Φ k ( w ) + (1 − λ )Φ k ( w ) , (30) holds for all λ ∈ [0 , and w , w ∈ R n . The pair (Φ , Φ ) is called a psd-convex-concave decomposition of Φ . Consider the bilinear inequality (21), and define the follow-ing matrix-valued mapping:
Φ( ˆ Q, ˆ δ, ˆ W ) = ψ ( ˆ Q, ˆ δ ) + ϕ ( ˆ W ) , (31)where ψ ( ˆ Q, ˆ δ ) = ˆ Q ¯ A + ¯ A (cid:62) ˆ Q ˆ QB e ˆ QEB (cid:62) e ˆ Q − ˆ δI E (cid:62) ˆ Q , (32) ϕ ( ˆ W ) = − ¯ A (cid:62) r ¯ A r A (cid:62) r A r − I . (33) Then, the following lemma shows that the pair ( ψ, − ϕ ) is apsd-convex-concave decomposition of Φ . Lemma 4.
The matrix-valued mapping
Φ( ˆ Q, ˆ δ, ˆ W ) in (31) ispsd-convex-concave.Proof. Note that the matrix ψ ( ˆ Q, ˆ δ ) in (31) is linear withrespect to ˆ Q and ˆ δ . Thus, it is immediate that ψ ( ˆ Q, ˆ δ ) is psd-convex. Then, the claim holds if − ϕ ( ˆ W ) in (31) is psd-convex.With the structure of ¯ A r in (24), the only nonlinear subma-trix in − ϕ ( ˆ W ) can be expressed as ¯ A (cid:62) r ¯ A r = (cid:20) ϕ a ( ˆ W ) (cid:21) , (34)with ϕ a ( ˆ W ) := S (cid:62) r ˆ M − ˆ B ˆ W ˆ B (cid:62) ( S + r ) (cid:62) S + r ˆ B ˆ W ˆ B (cid:62) ˆ M − S r . Then, showing the psd-convexity of − ϕ ( ˆ W ) in (33) is equiv-alent to prove that ϕ a ( ˆ W ) is psd-convex.Let W , W ∈ M , and denote W λ = λ ˆ W + (1 − λ ) ˆ W .For any λ ∈ [0 , , we have ϕ a ( W λ ) − λϕ a ( ˆ W ) − (1 − λ ) ϕ a ( ˆ W )= S (cid:62) r ˆ M − ˆ B W λ ˆ E (cid:62) ( S + r ) (cid:62) S + r ˆ E W λ ˆ B (cid:62) ˆ M − S r − tS (cid:62) r ˆ M − ˆ B ˆ W ˆ E (cid:62) ( S + r ) (cid:62) S + r ˆ E ˆ W ˆ B (cid:62) ˆ M − S r − (1 − λ ) S (cid:62) r ˆ M − ˆ B ˆ W ˆ E (cid:62) ( S + r ) (cid:62) S + r ˆ E ˆ W ˆ B (cid:62) ˆ M − S r = − λ (1 − λ )( V ( ˆ W − ˆ W ) V ( ˆ W − ˆ W ) V (cid:62) ) ≤ , (35)where V = S (cid:62) r ˆ M − ˆ B , and V = ˆ B (cid:62) ( S + r ) (cid:62) S + r ˆ B . Since − t (1 − t ) ≤ and V ( ˆ W − ˆ W ) V ( ˆ W − ˆ W ) V (cid:62) ≥ , itholds that ϕ a ( λ ˆ W + (1 − λ ) ˆ W ) ≤ λϕ a ( ˆ W ) + (1 − λ ) ϕ a ( ˆ W ) , which implies that the mapping ϕ a ( ˆ W ) is psd-convex from(30), i.e., − ϕ ( ˆ W ) is psd-convex. As a result, it follows fromDefinition 2 that the matrix-valued mapping Φ( ˆ Q, ˆ δ, ˆ W ) in(31) is psd-convex-concave.The psd-convex-concave decomposition in (31) allows usto linearize the optimization problem (29) at a stationarypoint ˆ W ∈ M . To simplify the optimization procedure, weintroduce a new optimization variable µ to eliminate theequality constraint ˆ B ˆ w = 0 in (13), where ˆ W = Diag( ˆ w ) ,and m is the number of edges in the reduced network.Let ¯ B ∈ R ¯ r × m be a full row rank matrix obtained byremoving linearly independent rows of the ˆ B ∈ R r × m , andit still holds that ¯ B ˆ w = 0 . Then, there exists a columnpermutation matrix P ∈ R m × m such that ¯ B ˆ w = (cid:2) ¯ B a ¯ B b (cid:3) P ˆ w = ¯ B a µ a + ¯ B b µ = 0 , with P ˆ w = (cid:20) µ a µ (cid:21) , where ¯ B a ∈ R ¯ r × ¯ r is full rank. µ ∈ R ¯ m + , and ¯ m = m − ¯ r , isdefined as the new optimization variable. Note that ˆ w = P (cid:62) (cid:20) − ¯ B − a ¯ B b I ¯ m (cid:21) µ, (36)which projects the weights ˆ w into ker( ˆ B ) . Thereby, we rewritethe constraint ˆ W ∈ M as µ ∈ R ¯ m + (36). Now, we redefine thematrix-valued mapping ϕ ( ˆ W ) in (31) as φ ( µ ) = ϕ ( ˆ W ) , (37) Algorithm 1
Iterative Edge Weighting
Input: L , F , H , Π , and a small scalar ˆ δ ∈ R + Output: ˆ W ∗ . Compute the incidence matrix ˆ B of the quotient graph ˆ G b . Choose an initial vector µ (0) ∈ R ¯ m + . Set iteration step: k ← . repeat Solve (39) to obtain the optimal solution µ ∗ . k ← k + 1 , and µ ( k ) ← µ ∗ . until | f ( µ ( k +1) ) − f ( µ ( k ) ) | ≤ ε . Compute ˆ w ∗ using (36), and return ˆ W ∗ ← Diag( ˆ w ∗ ) .which remains psd-convex due to the linear relation in (36).The derivative of the matrix-valued mapping φ ( µ ) at µ is alinear mapping Dφ : R ¯ m + → S (cid:96) , with (cid:96) = n + 2 r + p + pq − ,which is defined as Dφ ( µ )[ h ] = ¯ m (cid:88) i =1 h i ∂φ∂µ i ( µ ) , ∀ h ∈ R ¯ m . (38)Given a point µ ( k ) ∈ R ¯ m + , the linearized formulation of theproblem (29) at µ ( k ) is formulated as min ˆ Q> ,µ ∈ R ¯ m + f ( µ ) = tr( R ) (39)s.t. (cid:20) ˆ Q ˆ δC (cid:62) e ˆ δC e ˆ R (cid:21) > , ˆ δ ∈ R + , ˆ R = ˆ δR > ψ ( ˆ Q, ˆ δ ) + ϕ ( ˆ W ( k ) ) + Dφ ( µ ( k ) )[ µ − µ ( k ) ] < , where the derivative of φ ( µ ( k ) ) is given as Dφ ( µ ( k ) )[ µ − µ ( k ) ] := m (cid:88) i =1 ( µ i − µ ( k ) i ) ∂φ∂µ ( k ) i ( µ ( k ) ) , with j = 1 , · · · , ¯ m. Notice that the optimization problem(39) is convex , of which the global optimum can be solvedefficiently using standard SDP solvers e.g., SeDuMi [24].Based on Lemma 4 and (39), we are now ready to present analgorithmic approach for solving the minimization problem in(29) in an iterative fashion, see Algorithm 1, in which ε ∈ R + is a prefixed error tolerance determining whether to terminatethe iteration loop.The initial condition µ (0) can be chosen as an arbitraryvector with all strictly positive entries. With (36), it willguarantee ˆ W ∈ M , i.e., the reduced graph is balanced.Furthermore, we can also initialize µ using the outcome ofgraph clustering projection in [18], [19]. Specifically, from agiven clustering Π , we construct an initial reduced Laplacianmatrix in (11) as ˆ L (0) b := Π (cid:62) L b Π , with L b the Laplacianmatrix of the balanced graph G b . Then, the initial weight ofthe edge ( i, j ) in the quotient graph ˆ G b is the ( i, j ) -th entry of ˆ L (0) b . By doing so, µ (0) can be formed such that ˆ W ∈ M .The convergence analysis of Algorithm 1 follows naturallyfrom [23], and it means that a local optimum can be obtained.More importantly, if we select the initial condition from theclustering-based projection, it is guaranteed that, through iter-ation, the approximation accuracy of reduced-order networkmodel with the weights obtained by Algorithm 1 will be improved. In this sense, the approximations obtained by theproposed method is at least better than the ones obtained byclustering-based projection methods in e.g., [18], [19]. We willshow this merit from a numerical example in the next section.IV. I LLUSTRATIVE E XAMPLE
To illustrate the effectiveness of the proposed edge weight-ing approach, we implement it to a sensor network examplefrom [3], [18]. The topology of the this network is shownin Fig. 2, which consists of 14 strongly connected vertices,and all the edge weights are 1. In this example, two externalinput signals are injected into the network via vertices 2 and 7,respectively, and the states of vertices 9 and 10 are measured.Suppose that 5 clusters are given for this directed networkas C = { , , , } , C = { } , C = { , , } , C = { } , and C = { , , , , } , which leads to the quotientnetwork in Fig. 3, with incidence matrix ˆ B := − − − − − − − − , There are 8 edges in the quotient graph, and each edge isassigned with a symbolic weight as labeled in Fig. 3. Thesevariables, determining the reduction error, are to be determinedby our optimization approach.First, the parameterized reduced model in (10) of thequotient graph is generated with matrices ˆ L b = ˆ w + ˆ w − ˆ w − ˆ w − ˆ w ˆ w − ˆ w w + ˆ w + ˆ w − ˆ w − ˆ w − ˆ w ˆ w − ˆ w w , ˆ M = [ . . . . . ] , ˆ F b = (cid:104) . . (cid:105) (cid:62) , ˆ H = (cid:104) (cid:105) , and the weight vector ˆ w in (8) satisfy the following constraintsfor a balanced graph: ˆ w = ˆ w , ˆ w = ˆ w + ˆ w , ˆ w =ˆ w , ˆ w = ˆ w . Next, we implement Algorithm 1 to solve the optimizationproblem (29) with µ = [ ˆ w , ˆ w , ˆ w , ˆ w ] (cid:62) ∈ R as theoptimization variable. Particularly, the SeDuMi solver[24] is adopted to solve the convex problem (39).We choose the initial edge weights obtained by theclustering-based projection [18], [19], which gives ˆ w (0) =[0 . , . , . , . , . , . , . , . and the approximation error (cid:107) G e ( s ) (cid:107) H = 0 . .With ˆ δ = ε = 10 − , Algorithm 1 stops after 72iterations. The convergence trajectory of the resulting H reduction error is shown in Fig. 4. The finalsolution of the edge weights are given as ˆ w ∗ =[0 . , . , . , . , . , . , . , . ,which provides the approximation error (cid:107) G e ( s ) (cid:107) H = 0 . .Through iteration, the edge weighting method further reducesthe error by 41.93%, compared to the clustering-basedprojection. Therefore, our method can provide a reducednetwork systems with a better H approximation error. u u Fig. 2: A connected directed sensor network containing 14vertices, in which the red vertices are controlled, and theshadowed ones are measured. u u ^ w ^ w ^ w ^ w ^ w ^ w ^ w ^ w Fig. 3: The quotient graph obtained by clustering, where thecontrolled vertices are labeled with red color, and measuredvertices are indicated by shadow. The weights of the edges areparameters to be determined.
Interation
Edge WeightingClustering-based Projection
Fig. 4: Approximation errors of clustering-based projectionand the proposed edge weighting method.V. C
ONCLUSIONS
In this paper, the H model reduction problem for dy-namical networks consisting of diffusively coupled agentshas been formulated as a minimization problem, in whichthe edge weights in the reduced network are parameters tobe chosen. Necessary and sufficient conditions have beenproposed for constructing a set of optimal edge weights. Aniterative algorithm has been provided to search for the desirededge weights such that the H norm of the approximationerror is small. Finally, compared with the projection-basedmethod in [12], the feasibility of this method is illustrated byan example. The advantage of this proposed model reductionmethod is that not only the structure of the original networkhas been preserved but also the approximation error has beenoptimized. For future works, we will improve the effectivenessof the iterative algorithm such that the obtained solution isnot restricted to a local optimum. Moreover, an extension tonetworked high-order linear subsystems are also of interest. R EFERENCES[1] J. A. Fax and R. M. Murray, “Graph Laplacians and stabilization ofvehicle formations,” 2001.[2] F. D¨orfler, J. W. Simpson-Porco, and F. Bullo, “Electrical networks andalgebraic graph theory: Models, properties, and applications,”
Proceed-ings of the IEEE , vol. 106, no. 5, pp. 977–1005, 2018.[3] G. Scutari, S. Barbarossa, and L. Pescosolido, “Distributed decisionthrough self-synchronizing sensor networks in the presence of propa-gation delays and asymmetric channels,”
IEEE Transactions on SignalProcessing , vol. 56, no. 4, pp. 1667–1684, 2008.[4] A. V. Proskurnikov, A. S. Matveev, and M. Cao, “Opinion dynamics insocial networks with hostile camps: Consensus vs. polarization,”
IEEETransactions on Automatic Control , vol. 61, no. 6, pp. 1524–1536, 2015.[5] T. H. Summers and J. Lygeros, “Optimal sensor and actuator placementin complex dynamical networks,”
IFAC Proceedings Volumes , vol. 47,no. 3, pp. 3784–3789, 2014.[6] A. J. Gates and L. M. Rocha, “Control of complex networks requiresboth structure and dynamics,”
Scientific Reports , vol. 6, no. 1, pp. 1–11,2016.[7] T. Ishizaki, A. Chakrabortty, and J.-I. Imura, “Graph-theoretic analysisof power systems,”
Proceedings of the IEEE , vol. 106, no. 5, pp. 931–952, 2018.[8] A. C. Antoulas,
Approximation of Large-Scale Dynamical Systems .Philadelphia, USA: SIAM, 2005.[9] X. Cheng, J. M. A. Scherpen, and B. Besselink, “Balanced truncationof networked linear passive systems,”
Automatica , vol. 104, pp. 17–25,2019.[10] E. Bıyık and M. Arcak, “Area aggregation and time-scale modeling forsparse nonlinear networks,”
Systems & Control Letters , vol. 57, no. 2,pp. 142–149, 2008.[11] S. Rao, A. J. van der Schaft, and B. Jayawardhana, “A graph-theoretical approach for the analysis and model reduction of complex-balanced chemical reaction networks,”
Journal of Mathematical Chem-istry , vol. 51, no. 9, pp. 2401–2422, 2013.[12] N. Monshizadeh, H. L. Trentelman, and M. K. Camlibel, “Projection-based model reduction of multi-agent systems using graph partitions,”
IEEE Transactions on Control of Network Systems , vol. 1, pp. 145–154,Jun. 2014.[13] B. Besselink, H. Sandberg, and K. H. Johansson, “Clustering-basedmodel reduction of networked passive systems,”
IEEE Transactions onAutomatic Control , vol. 61, no. 10, pp. 2958–2973, Oct 2016.[14] T. Ishizaki, K. Kashima, J. I. Imura, and K. Aihara, “Model reductionand clusterization of large-scale bidirectional networks,”
IEEE Transac-tions on Automatic Control , vol. 59, pp. 48–63, 2014.[15] T. Ishizaki, K. Kashima, A. Girard, J.-i. Imura, L. Chen, and K. Aihara,“Clustered model reduction of positive directed networks,”
Automatica ,vol. 59, pp. 238–247, 2015.[16] H.-J. Jongsma, P. Mlinari´c, S. Grundel, P. Benner, and H. L. Trentelman,“Model reduction of linear multi-agent systems by clustering with h and H ∞ error bounds,” Mathematics of Control, Signals, and Systems ,vol. 30, no. 1, p. 6, 2018.[17] X. Cheng, Y. Kawano, and J. M. A. Scherpen, “Model reduction ofmulti-agent systems using dissimilarity-based clustering,”
IEEE Trans-actions on Automatic Control , vol. 64, no. 4, pp. 1663–1670, April 2019.[18] X. Cheng and J. M. A. Scherpen, “A new controllability Gramian forsemistable systems and its application to approximation of directednetworks,” in
Proceedings of IEEE 56th Conference on Decision andControl , Melbourne, Australia, Dec 2017, pp. 3823–3828.[19] X. Cheng and J. M. A. Scherpen, “Clustering-based model reduction oflaplacian dynamics with weakly connected topology,”
IEEE Transactionson Automatic Control , 2019.[20] N. Martin, P. Frasca, and C. Canudas-de Wit, “Large-scale networkreduction towards scale-free structure,”
IEEE Transactions on NetworkScience and Engineering , vol. 6, no. 4, pp. 711–723, 2018.[21] M. Mesbahi and M. Egerstedt,
Graph Theoretic Methods in MultiagentNetworks . Princeton University Press, 2010.[22] G. Pipeleers, B. Demeulenaere, J. Swevers, and L. Vandenberghe,“Extended LMI characterizations for stability and performance of linearsystems,”
Systems & Control Letters , vol. 58, no. 7, pp. 510–518, 2009.[23] Q. T. Dinh, S. Gumussoy, W. Michiels, and M. Diehl, “Combiningconvex-concave decompositions and linearization approaches for solvingBMIs, with application to static output feedback,”
IEEE Transactions onAutomatic Control , vol. 57, no. 6, pp. 1377–1390, 2011.[24] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimizationover symmetric cones,”