H2 model reduction of linear network systems by moment matching and optimization

TL;DR

This paper develops a novel model reduction method for linear network systems using moment matching and optimization to produce stable, reduced models that preserve structure and minimize H2 norm error.

Contribution

It introduces an optimization-based framework for model reduction that combines moment matching with structural constraints, including two solution approaches for nonconvex problems.

Findings

The first method uses semidefinite programming with a block diagonal Gramian assumption.

The second method employs a gradient projection technique for local optimality.

Application to a power network demonstrates the effectiveness of the proposed methods.

Abstract

In this paper we study the problem of model reduction of linear network systems. We aim at computing a reduced order stable approximation of the network with the same topology and optimal w.r.t. H2 norm error approximation. Our approach is based on time-domain moment matching framework, where we optimize over families of parameterized reduced order models matching a set of moments at arbitrary interpolation points. The parameterization of the low order models is in terms of the free parameters and of the interpolation points. For this family of parameterized models we formulate an optimization-based model reduction problem with the H2 norm of error approximation as objective function while the preservation of some structural and physical properties yields the constraints. This problem is nonconvex and we write it in terms of the Gramians of a minimal realization of the error system. We propose two solutions for this problem. The first solution assumes that the error system admits a block diagonal observability Gramian, allowing for a simple convex reformulation as semidefinite programming, but at the cost of some performance loss. We also derive sufficient conditions to guarantee block diagonalization of the Gramian. The second solution employs a gradient projection method for a smooth reformulation yielding (locally) optimal interpolation points and free parameters. The potential of the methods is illustrated on a power network.

Tables

Table 1. Table I: H 2 subscript 𝐻 2 H_{2} error norm of four reduction methods.

Method	SDP	Gradient	Projection [24]	Positive BT [21]
$H_2$ error	2.813	$5.075\cdot {10}^{-3}$	1.9103	1.7689

Equations

$$\dot{x}$$

$$K(s)=C(sI-A{)}^{-1}B.$$

$$A\Pi +BL=\Pi S,$$

$$\dot{\xi }$$

$$HP=C\Pi ,$$

$${\hat{\Sigma }}_{(S,G,L)}:\dot{\xi }$$

$$\hat{K}(s)=C\Pi (sI-S+GL{)}^{-1}G,$$

$$\Vert K{\Vert }_{H_2}=\sqrt{{\int }_{-\infty }^{\infty }|K(j\omega ){|}^{2}d\omega }.$$

$$\Vert K{\Vert }_{H_2}^2=C^{T}WC=B^{T}MB,$$

$${\dot{x}}_i$$

$${\dot{x}}_4=A_{44}{x}_4+A_{43}{x}_3+B_{4}u\text{and}{A}_{41}=A_{42}=0.$$

$$\left[\begin{array}[]{c|c}A&B\\ \hline\cr C&0\end{array}\right]=\left[\begin{array}[]{cccc|c}A_{11}&A_{12}&\dots&A_{1N}&B_{1}\\ A_{21}&A_{22}&\dots&A_{2N}&B_{2}\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ A_{N1}&A_{N2}&\dots&A_{NN}&B_{N}\\ \hline\cr C_{1}&C_{2}&\dots&C_{N}&0\end{array}\right],$$

$$A_{ij}=0\forall i\in [N],j\notin {\mathcal{N}}_i.$$

$${\dot{\xi }}_i$$

$$\psi =H\xi ,$$

$$\left[\begin{array}[]{c|c}F&G\\ \hline\cr H&0\end{array}\right]=\left[\begin{array}[]{cccc|c}F_{11}&F_{12}&\dots&F_{1N}&G_{1}\\ F_{21}&F_{22}&\dots&F_{2N}&G_{2}\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ F_{N1}&F_{N2}&\dots&F_{NN}&G_{N}\\ \hline\cr H_{1}&H_{2}&\dots&H_{N}&0\end{array}\right],$$

$${\mathcal{K}}_e=K-\hat{K},$$

$${\mathcal{K}}_e(s)={\mathcal{C}}_e(sI-{\mathcal{A}}_e{)}^{-1}{\mathcal{B}}_e,$$

$${\mathcal{A}}_e=\left[\begin{array}{cc}A& 0\\ 0& S-GL\end{array}\right],{\mathcal{B}}_e=\left[\begin{array}{c}B\\ G\end{array}\right],{\mathcal{C}}_e=C\left[\begin{array}{cc}I& -\Pi \end{array}\right].$$

$${\mathcal{A}}_e\mathcal{W}+\mathcal{W}{\mathcal{A}}_e^T+{\mathcal{B}}_e{\mathcal{B}}_e^T$$

$${\mathcal{A}}_e^T\mathcal{M}+\mathcal{M}{\mathcal{A}}_e+{\mathcal{C}}_e^T{\mathcal{C}}_e$$

$$\mathcal{W}=\left[\begin{array}{cc}{W}_{11}& W_{12}\\ W_{12}^T& W_{22}\end{array}\right],\mathcal{M}=\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{12}^T& M_{22}\end{array}\right].$$

$$\mathcal{G}=\{(S,G):(S-GL{)}_{ij}=0\forall i\in [N],j\notin {\mathcal{N}}_i\}.$$

$$\mathcal{R}=\left\{(S,G):\sigma (S-GL)\subset {\mathbb{C}}^{-},(S,G)\in \mathcal{G}\right\}.$$

$$\underset{(S,G)\in \mathcal{R}}{\min }\Vert {\mathcal{K}}_e{\Vert }_2^2$$

$$=\underset{(S,G)\in \mathcal{R},\mathcal{M}\text{s.t.}\text{eqref}{M}_{L}yap}{\min }{\left[\begin{array}{c}B\\ G\end{array}\right]}^T\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{12}^T& M_{22}\end{array}\right]\left[\begin{array}{c}B\\ G\end{array}\right].$$

$$\underset{(S,G,\mathcal{M},\Pi )}{\min }\text{Trace}({\mathcal{B}}_e^T\mathcal{M}{\mathcal{B}}_e)$$

$$\text{s.t.:}A\Pi +BL=\Pi S,\sigma (S-GL)\subset {\mathbb{C}}^{-}$$

$$(S,G)\in \mathcal{G},{\mathcal{A}}_e^T\mathcal{M}+\mathcal{M}{\mathcal{A}}_e+{\mathcal{C}}_e^T{\mathcal{C}}_e=0.$$

$$\underset{(S,G,\mathcal{M})}{\min }\text{Trace}({\mathcal{B}}_e^T\mathcal{M}{\mathcal{B}}_e)$$

$$\text{s.t.:}(S,G)\in \mathcal{R}\left(\iff \sigma (S-GL)\subset {\mathbb{C}}^{-},(S-GL{)}_{ij}=0j\notin {\mathcal{N}}_i\right)$$

$${\mathcal{A}}_e^T\mathcal{M}+\mathcal{M}{\mathcal{A}}_e+{\mathcal{C}}_e^T{\mathcal{C}}_e=0,$$

$$\underset{(S,G)\in \mathcal{G},\mathcal{M}⪰0}{\min }\text{Trace}({\mathcal{B}}_e^T\mathcal{M}{\mathcal{B}}_e)$$

Full text

I Introduction

Complex network systems consist of multiple interacting dynamical subsystems, interconnected through a graph enabling the subsystems to share information, coordinate their activities and have self-control mechanisms [16]. However, the corresponding models of network systems are too complex and difficult to analyze, rendering it is almost impossible to systematically develop operating and/or open/closed-loop control algorithms. Therefore, we need approximation models to do analysis, simulation and control.

State-of-the-art: The problem of model reduction of interconnected systems has been long studied in different frameworks, see e.g. [20] and references therein for a survey. There are two main existing approaches.

A first approach stemming from mathematics considers network systems as static mathematical objects. The reduction is treated with topological objectives, focusing on obtaining a reduced network abstracting a large-scale network by merging groups of nodes into super-nodes (so-called clustering) [4]. For example, [15] aims at preserving stability and synchronisation of the system, [5] preserves an interconnection structure and synchronization by aggregating subsystems with similar frequency responses, [10] provides a reduced system with a dynamical behaviour close to the initial system while preserving several properties for control purpose, and [14] proposes a network reduction method preserving the flow network property and the reduced graph to be scaled-free.

A second approach comes from systems and control theory. The aim is to reduce the network system by preserving consistency/structure in the network. Here, one category of results are in the framework of stability preserving balanced truncation, see e.g. [23, 22], where the balancing yields the so-called structured Hankel singular values (invariants showing the importance of subsystems states with respect to a chosen input-output map for the whole network). The states with the lowest structured Hankel singular values are truncated directly from the full model, resulting in a low order stable network satisfying the given interconnection map. A second category of results are based on interpolatory methods [3], as e.g. in [13, 19, 26], see also [12, 6] for earlier results. Here, structured Krylov projections are applied directly on the entire network to preserve the topology. Note that none of the presented results reduce the number of nodes (subsystems) or alter the interconnection map of the network.

Motivation: In this paper we also consider the second approach of approximating the subsystems of a network. If the number of subsystems is large, we first reduce their number using clustering techniques and then perform subsystem approximation. To the best of our knowledge, in the time-domain moment matching framework [2, 11], finding a reduced model of the network optimally w.r.t. the $H_{2}$-norm of the approximation error in the family of $\nu$ order models that matches a set of $\nu$ moments, while preserving the network topology and stability, is an open question. Some initial progress has been made recently in [17] for general linear systems. However, a direct application of this approach on the full model of the network system does not preserve network topology. Hence, this unsolved problem motivates our work here.

Contributions: In this paper we provide a systematic procedure for approximating the subsystems of a network optimally while preserving the network topology and stability. The proposed procedure is based on time-domain moment matching, where families of parametrized low order models matching a set of moments at arbitrary interpolation points are computed. Here we use the free parameters and the interpolation points defining the parameterization to find the optimal approximation of the network measured in terms of the $H_{2}$-norm of the error system. We formulate an optimization problem with the $H_{2}$-norm of the approximation error as objective function, while the preservation of some structural and physical properties yields the constraints. The problem is nonconvex and we prove that it can be written in terms of the controllability/observability Gramians of a minimal realization of the error system. We propose two solutions for solving this problem. The first solution assumes that the error system admits a block diagonal observability Gramian, allowing for a simple convex reformulation as semidefinite programming, at the cost of some performance loss. Also, sufficient conditions are derived to guarantee block diagonal Gramians. The second solution employs a gradient projection method for a smooth reformulation yielding $H_{2}$ (locally) optimal reduced order models. Both solutions provide (optimal) stable low order network models, parameterized in the interpolation points and in the free parameters, matching a set of moments and preserving the interconnection map of the original network. The efficiency of the methods is illustrated on a positive network and on a multiple area power system.

Content: The paper is organized as follows. In Section II we briefly review the time-domain moment matching model reduction of linear systems. In Section III we formulate the optimization-based $H_{2}$-norm moment matching model reduction problem with a Gramian-type cost function. In Sections IV and V, we propose two numerical optimization methods for solving this model reduction problem and extensions are given in Section VI. In Section VII we illustrate the efficiency of our theory on several network examples.

Notation: $\mathbb{R}$ and $\mathbb{C}$ denotes the set of real and complex numbers, respectively. For a positive integer $N$ we denote by $[N]=\{1,\cdots,N\}$. For a matrix $A\in\mathbb{R}^{n\times n}$, $\sigma(A)$ denotes the set of its eigenvalues and $A_{ij}$ indicates the $(i,j)$ matrix block of $A$ of appropriate dimension.

II Preliminaries

In this section, we briefly present the main results on time-domain moment matching for linear systems [11].

II-A *Moments and moment matching

Consider the linear time-invariant system:

[TABLE]

where $x\in\mathbb{R}^{n}$ is the state of the system, $u\in\mathbb{R}^{m}$ is the input and $y\in\mathbb{R}^{p}$ is the output, respectively. Consequently, system matrices $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{n\times m}$ and $C\in\mathbb{R}^{p\times n}$. Throughout the paper we assume that the system is stable (i.e. $\sigma(A)\subset\mathbb{C}^{-}$) and that (1) is a minimal (i.e. controllable and observable) realization of the transfer function:

[TABLE]

The moments of linear system (1) at a point $s\in\mathbb{C}$ on the complex plane are defined as follows:

Definition 1

[1, 2, 7]** The $k$-moment of the system (1) at the point $s\notin\sigma(A)$, along the direction $\ell\in\mathbb{C}^{m}$ is $\eta_{k}(s;\ell)={(-1)^{k}}/{k!}\ {\operatorname{d\!}^{k}\!K(s)}/{\operatorname{d\!}s^{k}}\ell\in\mathbb{C}^{p},$ with $k\geq 0$ integer.

Consider the matrix $S\in\mathbb{R}^{\nu\times\nu}$, with $\nu\leq n$ and $\sigma(S)=\{s_{i}:i\geq 0\;\text{integer}\}\subset\mathbb{C}$. Let $L=[\ell_{1}\cdots\ell_{\nu}]\in\mathbb{R}^{m\times\nu}$ be such that the pair $(L,S)$ is observable. Since the system is assumed minimal, the Sylvester equation:

[TABLE]

has the unique solution $\Pi\in\mathbb{R}^{n\times\nu}$ with ${\rm rank}\ \Pi=\nu$ provided that $\sigma(A)\cap\sigma(S)=\emptyset$ [1]. Then, the moments of a system can be characterized as follows:

Proposition 1

[11]** Consider the system (1) and let $\Pi$ be the unique solution of equation (3). Then, at the interpolation points $s_{i}\in\sigma(S)$, the moments of the system (1) along directions $\ell_{i}$, $\eta_{k}(s_{i};\ell_{i})$, with $i,k\geq 0$ integers, are characterized by the matrix $C\Pi$.

We now present the moment matching property and the reduced order model satisfying it:

Proposition 2

[2, 11*]**

Consider the $\nu$ order linear system:*

[TABLE]

with the state $\xi\in\mathbb{R}^{\nu}$, input $u\in\mathbb{R}^{m}$ and output $\psi\in\mathbb{R}^{p}$. Here, $\nu\leq n$, $F\in\mathbb{R}^{\nu\times\nu},G\in\mathbb{R}^{\nu\times m}$ and $H\in\mathbb{R}^{p\times\nu}$. Assuming $\sigma(F)\cap\sigma(S)=\emptyset$, then the reduced order system (4) matches the moments of (1) at $\sigma(S)$ if and only if:

[TABLE]

where the invertible matrix $P\in\mathbb{R}^{\nu\times\nu}$ is the unique solution of the Sylvester equation $FP+GL=PS.$

Note that the invertible matrix $P$ in Proposition 2 is merely a coordinate transformation. Hence, taking $P=I_{\nu}$ yields a parameterized $\nu$ order model (with the free parameters $(G,L)$ and the interpolation points matrix $S$) achieving moment matching at $\sigma(S)$, as shown in the next result:

Proposition 3

[2, 11*]**

Assume that $(L,S)$ is observable and $\sigma(A)\cap\sigma(S)=\emptyset$. Consider the $\nu$ order linear system:*

[TABLE]

with $\nu\leq n$ and the transfer function:

[TABLE]

where $\Pi$ is the unique solution of (3). Assuming that $\sigma(S-GL)\cap\sigma(S)=\emptyset$, then the system (6), with the transfer function (7), is a reduced order model of (1) parametrized in $S,G$ and $L$, matching the moments $C\Pi$ of system (1) at $\sigma(S)$.

Remark 1

The system $\widehat{\Sigma}_{(S,G,L)}$ in (7) describes a $\nu$ order approximation of (1) that achieves moment matching at $\sigma(S)$. Since $L$ is only used in $\widehat{\Sigma}_{(S,G,L)}$ to ensure observability of the pair $(L,S)$ and since observability is generic, then, without loss of generality, we fix a priori matrix $L$, i.e., we fix the directions $\ell_{i}$ to compute moments along, see e.g. [2, 11]. Hence, in the rest of the paper we consider $\widehat{\Sigma}_{(S,G,L)}=\widehat{\Sigma}_{(S,G)}$, defining a family of $\nu$ order models matching $\nu$ moments along fixed directions $\ell_{i}$ of system (1) at $\sigma(S)$, for all $G$, such that:

1. i)

$\widehat{\Sigma}_{(S,G)}$* is parametrized in $(S,G)$* 2. ii)

$\sigma(S-GL)\cap\sigma(S)=\emptyset$.

II-B $H_{2}$-norm based on the Gramians of linear systems

We now briefly recall the definition and computation of the $H_{2}$-norm of a linear system. For minimal stable system (1) with transfer function (2), the ${H}_{2}$-norm is defined as [8]:

[TABLE]

This norm can be written explicitly in matrix form as [8]:

[TABLE]

where $W$ is the controlability Gramian and $M$ is the observability Gramian of the linear system (1).

III Optimal $H_{2}$ model reduction formulation of linear network systems

In this section we formulate a model reduction problem, yileding a family of parametrized models for each subsystem of a linear network without altering its structure. To determine the best approximation in terms of the $H_{2}$-norm of the error, we propose an optimization formulation, with the $H_{2}$-norm of the error as objective function, while stability and structure are imposed as constraints.

III-A Linear network systems

We perform model reduction for linear network systems consisting of $N$ interconnected subsystems, with dynamics defined by the linear state space equations:

[TABLE]

where $x_{i}\in\mathbb{R}^{n_{i}}$ represents the state of the $i$th subsystem, $u\in\mathbb{R}^{m}$ is the common input, $A_{ii}\in\mathbb{R}^{n_{i}\times n_{i}}$, $B_{i}\in\mathbb{R}^{n_{i}\times m}$, and $A_{ij}\in\mathbb{R}^{n_{i}\times n_{j}}$. The index set ${\mathcal{N}}_{i}\subseteq[N]$ contains the index $i$ and all the indices of the subsystems which interact with the subsystem $i$. Thus, in (9) we consider that each subsystem is influenced through the states of the neighboring subsystems. For a more general network description see Section VI.

For example, consider the network system in Figure 1, where the arrows indicate the interactions between the subsystems $\Sigma_{1},\Sigma_{2},\Sigma_{3}$ and $\Sigma_{4}$. If we consider the fourth subsystem $\Sigma_{4}$, we have ${\mathcal{N}}_{4}=\{3,4\}$ and hence:

[TABLE]

For model reduction, we also express the dynamics of the entire network system in the compact form (1), $\dot{x}=Ax+Bu$, where $x=[x_{1}^{T}\cdots x_{N}^{T}]^{T}\in\mathbb{R}^{n}$, with $n=\sum_{i=1}^{N}n_{i}$, denotes the states of the entire network and the input $u\in\mathbb{R}^{m}$. As output of the network system we consider a linear combination of the states of each subsystem: $y=Cx$, where $C\in\mathbb{R}^{p\times n}$. Then, the system matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$ and $C=[C_{1}\dots C_{N}]\in\mathbb{R}^{p\times n}$, with $C_{i}\in\mathbb{R}^{p\times n_{i}}$, are given by:

[TABLE]

where we have that the matrix block $(i,j)$ of $A$ satisfies:

[TABLE]

Recall that we assume stable networks, i.e. $\sigma(A)\subset\mathbb{C}^{-}$. Note that the dimension $n$ of entire network system is usually too large, so that it is almost impossible to develop open or closed-loop control algorithms in a systematic way. Therefore, we need to obtain approximation models of the network system (10), useful for analysis, simulation and control. Unfortunately, moment matching-based model reduction techniques, such as in [17], do not preserve the network structure. However, working with an approximation violating basic network constraints it always leaves the question of how conclusive the results on this basis are.

III-B Optimal $H_{2}$ moment matching-based model reduction problem preserving network structure

In this section, we formulate the component model reduction problem of the linear network system (9) with the network structure given in (10). Recall that $A_{ij}=0$ in (10) if $i\in[N],\;j\not\in{\mathcal{N}}_{i}$ and the dimension of the entire network system is $n=\sum_{i=1}^{N}n_{i}$. The goal is to perform model order reduction such that the network structure is preserved, i.e. compute reduced order models for each subsystem of the form:

[TABLE]

where $\xi_{i}\in\mathbb{R}^{\nu_{i}}$, with $\nu_{i}\leq n_{i}$, represents the reduced state of the $i$th subsystem, $u\in\mathbb{R}^{m}$ is the common input, $F_{ij}\in\mathbb{R}^{\nu_{i}\times\nu_{j}}$ and $G_{i}\in\mathbb{R}^{\nu_{i}\times m}$. Moreover, we want to preserve the network structure, that is $F_{ij}=0$ if $i\in[N],\;j\not\in{\mathcal{N}}_{i}$. If the number of subsystems $N$ is large, we first reduce their number using existing clustering techniques (see Section I) and then perform the subsystem approximation procedure described below. Note that the dimension of the whole reduced model $\dot{\xi}=F\xi+Gu$ is $\nu=\sum_{i=1}^{N}\nu_{i}$, where $\xi=[\xi_{1}^{T}\cdots\xi_{N}^{T}]^{T}\in\mathbb{R}^{\nu}$ denotes the full state of the reduced model. We also define the output of the reduced network:

[TABLE]

where $H\in\mathbb{R}^{p\times\nu}$. The matrices of the reduced network system (11) are written in a compact form as:

[TABLE]

where $H_{j}\!\in\!\mathbb{R}^{p\times\nu_{j}}$ and consider also $L=[L_{1}\dots L_{N}]$, with blocks $L_{i}\in\mathbb{R}^{m\times\nu_{i}}$. We want $F_{ij}=0$ for $i\in[N],\;j\not\in{\mathcal{N}}_{i}$. Based on the parametrizations of the reduced model (6) given in terms of the interpolation points matrix $S$ and of the free parameters $G$ (i.e. $F=S-GL,G$ and $H=C\Pi$) and using the $H_{2}$-norm of the approximation error as objective function and physical and structural restrictions as constraints, we derive below an optimization problem to determine the minimal approximation. More precisely, the optimal $H_{2}$ model reduction problem by moment matching is formulated as:

Problem 1

Given a linear network system (1) with the subsystem matrices (10) and the transfer function $K(s)$ as in (2) and the directions of moments $L$, find a reduced order linear network $\widehat{\Sigma}_{(S,G)}$ of the form (6) with subsystem matrices (12) and the transfer function $\widehat{K}(s)$ as in (7), parametrized in the interpolations matrix $S$ and the free parameters $G$, that matches $\nu$ moments of (1) at $\sigma(S)$ and satisfies the constraints:

1. (i)

the $H_{2}$-norm of the error system $\|K-\widehat{K}\|_{2}$ is minimal 2. (ii)

the reduced model $\widehat{K}$ is stable (i.e. $\sigma(S-GL)\subset\mathbb{C}^{-}$) 3. (iii)

the matrix $F=S-GL$ preserves the network topology of $A$ (i.e for all $i\in[N]$, $(S-GL)_{ij}=0$ if $j\notin{\cal N}_{i}$) 4. (iv)

$\sigma(S)\cap\sigma(A)=\emptyset$, $\sigma(S)\cap\sigma(S-GL)=\emptyset$ and the pair $(L,S)$ observable.

However, it is difficult to deal with the restrictions (iv): $(S,L)$ observable, $\sigma(S)\cap\sigma(A)=\emptyset$, and $\sigma(S)\cap\sigma(S-GL)=\emptyset$. One possibility is to fix $S$ and $L$ such that (iv) is automatically satisfied (e.g., without loss of generality, take $S=\operatorname{diag}(s_{1},\dots,s_{\nu})$, with $s_{i}\in\mathbb{C}^{+}\;\forall i$, and $L=[\ell_{1}\dots\ell_{\nu}]$, with $\ell_{i}\not=0\;\forall i$) and search only for the free parameters $G$. All our results hold for this choice. Another possibility, which we also follow in this paper, is to fix $L$. Note that since model reduction procedures usually render $S$ unstable, while $A$ and $S-GL$ are stable, the first two constraints in (iv) are automatically satisfied. Moreover, since observability is generic, by Remark 1, also observability of $(L,S)$ holds. Hence, the constraints (iv) are not imposed in the numerical algorithms, but will be checked after yielding a solution to Problem 1. Therefore, in the sequel we propose an optimization formulation of Problem 1, without constraints (iv) in the unknowns $S$ and $G$, while $L$ is fixed a priori. Under these settings Problem 1 can be recast in terms of the Gramians of the realization of the error system:

[TABLE]

with $\widehat{K}$ from (7), parameterized in $(S,G)$. Let $(\mathcal{A}_{e},\mathcal{B}_{e},\mathcal{C}_{e})$ be a state-space realization of the error transfer function $\mathcal{K}_{e}$:

[TABLE]

where

[TABLE]

Denote the controllability and the observability Gramians of (13) by $\mathcal{W}$ and $\mathcal{M}$, respectively. They are solutions of the Lyapunov equations [1]:

[TABLE]

Since we assume $\sigma(A),\sigma(S-GL)\subset\mathbb{C}^{-}$, then the matrix $\mathcal{A}_{e}$ is also stable, i.e., $\sigma(\mathcal{A}_{e})\subset\mathbb{C}^{-}$. Hence, there exist unique positive semidefinte solutions $\mathcal{W}$ and $\mathcal{M}$ of equations (14), respectively. We partition $\mathcal{W}$ and $\mathcal{M}$ following the two block structure of the error matrix $\mathcal{A}_{e}$:

[TABLE]

The communication graph between subsystems imposes a constraint on the admissible parameterizations. If for subsystem $i\in[N]$, $j\notin{\cal N}_{i}$, then there is no communication link from subsystem $j$ to subsystem $i$. Thus, in the reduced model, the subsystem $i$ cannot be also influenced by the subsystem $j$. This leads to the structured constraint $\mathcal{G}$, defined by:

[TABLE]

Let us define the feasible set for the reduced model:

[TABLE]

By (8) and the simplification stated above, Problem 1 becomes:

[TABLE]

Using the matrix notations above, optimization problem (16) can be written in matrix form explicitly as:

[TABLE]

Note that, the Sylvester equation $A\Pi+BL=\Pi S$ need not be solved since, by [11, Lemma 1], we can take $\Pi=VT$, with $V$ a certain Krylov projection and $T$ some non-singular matrix. Therefore, we get the following simplified nonconvex optimization formulation for Problem 1:

[TABLE]

In the rest of the paper we derive several numerical procedures for solving the nonconvex problem (18), whose optimal solution yields a stable reduced order model of dimension $\nu$ of the linear network system, which preserves the topology of the network and minimizes the $H_{2}$-norm of the error system.

IV Convex model reduction using block diagonal Gramians

The nonconvex problem (18) can be written equivalently in terms of matrix inequalities (semidefinite programming):

[TABLE]

where $L$ is fixed a priori. Clearly, semidefinte program (SDP) (19) is not convex since it contains bilinear matrix inequalities (BMIs). However, our next result shows that we can obtain a suboptimal solution through convex SDP using a simple assumption that the error system admits a block diagonal observability Gramian. While diagonal Gramians have recently been exploited in the balanced truncation model reduction of positive systems [9], the application of block diagonal Gramians on structured moment matching model reduction of general network systems is discussed in our paper.

Theorem 1

If the convex SDP relaxation:

[TABLE]

has a solution, then we can recover a suboptimal solution of the model reduction Problem 1 expressed in terms of the SDP problem (19) through the relations:

[TABLE]

Proof:

Using the block form of the Gramian $\mathcal{M}$, (19) yields the equivalent SDP problem (21).

Note that problem (21) is not convex since if we assume $M_{12}\not=0$, then we cannot convexify the previous BMIs since we need to define $M_{12}G=Z_{12}$ and $M_{22}G=Z_{22}$ and require $\mathcal{M}\succeq 0$. However, if we assume for the Gramian $\mathcal{M}$ the block $M_{12}=0$ and $M_{22}$ block diagonal, then problem (21) can be recast as a convex SDP. More precisely, if we introduce additional variables, then we get the following SDP:

[TABLE]

Letting $Z_{22}=M_{22}G,\Theta_{22}=M_{22}S$ and using the Schur complement, problem (22) becomes the convex SDP (1). However, this change of variables $M_{22}(S-GL)=\Theta_{22}-Z_{22}L$ is in general not suitable when imposing the structured constraints $\mathcal{G}$ on $S-GL$. Although the constraint on the parameterization $(S,G)\in\mathcal{G}$ is linear and thus convex, the corresponding constraint on $\Theta_{22},Z_{22}$ and $M_{22}$ (i.e. $M_{22}^{-1}(\Theta_{22}-Z_{22}L)\in\mathcal{G}$) is nonlinear and consequently nonconvex. If we restrict the structure of $M_{22}$, assuming it is block diagonal with the block sizes compatible to those of the reduced subsystems, i.e., $M_{22}=\text{diag}(M_{22}^{1},\dots,M_{22}^{N})$, with $M_{22}^{i}\in\mathbb{R}^{\nu_{i}\times\nu_{i}}$, the structured constraints are naturally guaranteed:

[TABLE]

Note that the block diagonal assumption on $M_{22}$ is a sufficient condition for (23) given an arbitrary network structure $\mathcal{G}$. Moreover, we can recover a suboptimal solution of the original problem through the relations: $G=M_{22}^{-1}Z_{22},S=M_{22}^{-1}\Theta_{22}$ and $\mathcal{M}=\text{diag}(M_{11},M_{22})$. ∎

Clearly, (1) is a suboptimal solution of the original SDP problem (19) since we restrict the Gramian matrix $\mathcal{M}$ to have the blocks $M_{12}=0$ and $M_{22}$ to be block diagonal. Hence, (1) is a convex SDP relaxation of the original problem (19).

IV-A Sufficient conditions on block diagonal Gramians

As we can see from the proof of Theorem 1, the block diagonal assumption on the observability Gramian is crucial. In this section we derive sufficient conditions for the feasibility of the SDP relaxation (1) in the case of general linear systems, i.e. sufficient conditions to guarantee that the error system admits a block diagonal observability Gramian. For this we need the following result that holds for any two vectors $u$ and $v$:

[TABLE]

This inequality follows from the relation $(u-Pv)^{T}P^{-1}(u-Pv)\succeq 0$ for $P\succ 0$. We are interested in deriving sufficient conditions to guarantee that the SDP (19) admits a feasible triplet $(S,G,\mathcal{M})$ with $\mathcal{M}$ of block diagonal form and consequently the SDP relaxation (1) is well-defined.

Theorem 2

Given the stable minimal system (1) there exists a stable reduced order model (6) such that the error system admits a block diagonal observability Gramian $\mathcal{M}=\operatorname{diag}(M_{11},M_{22})$, with $M_{22}=\text{diag}(M_{22}^{1},\dots,M_{22}^{N})$, if the following conditions hold:

[TABLE]

for some matrices $M_{11},M_{22},P\succeq 0$ and $(S,G)\in{\cal G}$.

Proof:

Note that the feasible set of SDP (19) is nonempty if $\mathcal{M}\succeq 0$, $(S,G)\in{\cal G}$ and the following inequality holds:

[TABLE]

which, using that $\mathcal{M}=\operatorname{diag}(M_{11},M_{22})$, is equivalent to

[TABLE]

Using now (24) in the last term we get that:

[TABLE]

for all $x,\xi$. Consequently, if the inequality:

[TABLE]

holds for all $x,\xi$, then (26) also holds. This proves the sufficient conditions (25). ∎

This theorem provides sufficient conditions and a procedure for constructing a reduced order network model for which the corresponding error system admits a bock diagonal observability Gramian. Indeed, let us, for example, fix $L,M_{22}=I_{\nu}$, and some $P\succ 0$. Then, the existence of a solution $(S,G,\Pi)$ satisfying $(S,G)\in{\cal G}$ of the system of equations:

[TABLE]

together with the existence of an $M_{11}\succ 0$ satisfying $A^{T}M_{11}+M_{11}A+C^{T}C+C^{T}(C\Pi)P^{-1}(C\Pi)^{T}C\preceq 0$ guarantee that we have a reduced order model for which the corresponding observability Gramian of the error system is block diagonal (i.e. $M_{12}=0$ and $M_{22}=I_{\nu}$ is also block diagonal). For example, we can fix matrices $S$ , $L$ and $G$ such that $S-GL$ is stable and has the same network structure as $A$ (just take stable diagonal matrix $S$ and $G=0$). Let $\Pi$ be the solution of $A\Pi+BL=\Pi S$ and define $P=-(S-GL)^{T}-(S-GL)-(C\Pi)^{T}(C\Pi)$. If the resulting $P\succ 0$ and if there exists $M_{11}\succ 0$ satisfying $A^{T}M_{11}+M_{11}A+C^{T}C+C^{T}(C\Pi)P^{-1}(C\Pi)^{T}C\preceq 0$, then we obtain a stable reduced model preserving the network structure and for which the corresponding error system admits a block diagonal Gramian $\mathcal{M}=\text{diag}(M_{22},I_{\nu})$.

From our best knowledge, the most common dynamical systems that admit block diagonal Gramians are the positive systems. The system matrices for these systems satisfy: all off-diagonal elements of the matrix $A$ and all the entries of the matrices $B$ and $C$ are non-negative. Positive systems occur in modelling of applications with special structures from, e.g., biomedicine, economics, networks (see [24, 21]).

V Nonconvex model reduction using projected gradient

The SDP approximation (1) offers a convex way to solve the nonconvex problem (18) of Problem 1 at the cost of some performance loss in general. We can also derive a numerical approach based on projected gradient method to solve (18). Here, our idea is to use a partial minimization approach (see Appendix -A) to (18) leading to a smooth reformulation, and apply the gradient projection method to get a (locally) optimal solution for (18). More precisely, consider the nonconvex optimization problem (18), where $L$ fixed a priori and $\mathcal{A}_{e}=\mathcal{A}_{e}(S,G)$. Then, the following partial minimization holds for problem (18):

[TABLE]

However, if $S-GL$ and $A$ are stable, then there exists unique $\mathcal{M}=\mathcal{M}(S,G)$ positive semidefinte solution of the Lyapunov equation:

[TABLE]

Hence, for any pair $(S,G)$ stable, the partial minimization in $\mathcal{M}$ leads to an optimal value $f(S,G)=\min_{\mathcal{M}:\mathcal{A}_{e}^{T}\mathcal{M}+\mathcal{M}\mathcal{A}_{e}+\mathcal{C}_{e}^{T}\mathcal{C}_{e}=0}\text{Trace}(\mathcal{B}_{e}^{T}\mathcal{M}\mathcal{B}_{e}))$, which can be written explicitly as:

[TABLE]

where $\mathcal{M}(S,G)$ is the unique positive semidefinite solution of the Lyapunov equation:

[TABLE]

Thus, we get the following equivalent reformulation for (18):

[TABLE]

For solving the equivalent problem (28) we can apply any first- or second-order optimization method. Hence, we need to compute the gradient and eventually the Hessian of the objective function $f$. In the sequel, we show that we can compute the gradient of the objective function of (28) solving two Lyapunov equations. Indeed, since $\text{Trace}(MN)=\text{Trace}(NM)$ for any matrices $M,N$ of compatible sizes, the objective function of (28) becomes:

[TABLE]

Theorem 3

The objective function $f$ of optimization problem (28) is differentiable on the set of stable matrices ${\cal D}=\left\{(S,G):\;\sigma(S-GL)\subset\mathbb{C}^{-}\right\}$ and the gradient of $f$ at any pair of matrices $(S,G)\in{\cal D}$ is given by $\nabla f(S,G)=[\nabla_{S}f(S,G)\;\nabla_{G}f(S,G)]\in{\mathbb{R}}^{\nu\times(\nu+m)}$ with:

[TABLE]

where $\mathcal{M}(S,G)$ solves the Lyapunov equation (V) and $\mathcal{W}(S,G)$ solves the Lyapunov equation from below:

[TABLE]

Proof:

To compute the gradient $\nabla f(S,G)$, we write the derivative $f^{\prime}(S,G)\operatorname{d\!}{(S,G)}$ for some $\operatorname{d\!}{(S,G)}\in\mathbb{R}^{\nu\times(\nu+m)}$ in gradient form using the trace. We introduce the gradient as:

[TABLE]

Then, we have:

[TABLE]

We compute separately the two terms in the above expression. Let us define:

[TABLE]

Since $(S,G)\in{\cal D}$ and ${\cal D}$ is an open set, then $\Phi_{\mathcal{M}}(S,G,\mathcal{M})\operatorname{d\!}\mathcal{M}$ given by:

[TABLE]

is surjective and also we have:

[TABLE]

Since $\Phi(S,G,\mathcal{M})+{\cal C}_{e}^{T}{\cal C}_{e}=0$, the Implicit Function Theorem yields the differentiability of $\mathcal{M}(S,G)$ and the relation:

[TABLE]

Moreover, by (14a) the Gramian $\mathcal{W}(S,G)$ is the unique solution of the Lyapunov equation (30). Subtracting (31) multiplied by $\mathcal{W}(S,G)$ to the left from (30) multiplied by $\mathcal{M}^{\prime}(S,G)$ to the right, taking the trace, and reducing the appropriate terms, we get the relation:

[TABLE]

Similarly, for the second term we get:

[TABLE]

Hence, combining (32) and (33), using the block structure of $\mathcal{W}$ and $\mathcal{M}$, and the definition of trace, we obtain easily the closed form expression for the gradient from (3). ∎

The previous theorem also yields the necessary optimality condition for the model reduction optimization problem (28):

Lemma 1

Let the block presentations $S=(S_{ij})_{i,j=1:N}$, $G=(G_{i})_{i=1:N}$ and $L=(L_{i})_{i=1:N}$, with $L$ fixed. If $(S^{*},G^{*})\in{\cal D}$ solves the optimization problem (28), then

[TABLE]

where the expressions of $\nabla_{S}f$ and $\nabla_{G}f$ are given in (3).

We can replace the open set ${\cal D}$ with any sublevel set:

[TABLE]

where $(S_{0},G_{0})\in{\cal D}$ is any pair of initial stable reduced order system matrices. By arguments as in [25] we can show that ${\cal L}(S_{0},G_{0})$ is a compact set. Then, the theorem of Weierstrass implies that for any given matrix $(S_{0},G_{0})\in{\cal D}$, the model reduction Problem 1 given by optimization formulation (28) has a global minimum in the sublevel set ${\cal L}(S_{0},G_{0})$. We can also show that the gradient $\nabla f(S,G)$ is Lipschitz continuous on the compact sublevel set ${\cal L}(S_{0},G_{0})$. We briefly sketch the proof of this statement. First we observe that $\mathcal{M}(S,G)$ and $\mathcal{W}(S,G)$ are continuous functions, since they are solutions of some algebraic linear systems. Moreover, there exists finite $\ell_{M}>0$ such that for all $(S,G),(S^{\prime},G^{\prime})\in{\cal L}(S_{0},G_{0})$:

[TABLE]

Then, using the expression of $\nabla f(S,G)$, compactness of ${\cal L}(S_{0},G_{0})$, continuity of $\mathcal{M}(S,G)$ and $\mathcal{W}(S,G)$, and the previous relation we conclude that there exists $\ell_{f}>0$ such that for all $(S,G),(S^{\prime},G^{\prime})\in{\cal L}(S_{0},G_{0})$:

[TABLE]

This property of the gradient is useful when analyzing the convergence behavior of the projected gradient algorithm we propose below for solving optimization problem (28).

V-A Projected gradient method

We have proved that the nonconvex optimization problem (28) has differentiable objective function and its gradient is given by (3). Moreover, the gradient is smooth (i.e. Lipschitz continuous) on any compact set. Then, we can apply the projected gradient method for obtaining a (local) optimal solution of (28). Starting from an initial stable matrix pair satisfying the structured constraints, $(S_{0},G_{0})\in{\cal R}$, we consider the following update rule:

[TABLE]

where the stepsize $\alpha_{k}$ can be chosen by a backtracking procedure or constant in the interval $(0,\ 2/\ell_{f})$ (where $\ell_{f}$ denotes the Lpschitz constant of the gradient) [18]. Here $\Pi_{\cal G}\left(\nabla f\right)$ denotes the projection of the gradient of the objective function $\nabla f$ onto the convex set ${\cal G}$ describing the network structure. Note that the projection of $\nabla f$ onto the convex set ${\cal G}$ is straightforward and computationally cheap: using the expressions of the gradient from (3), we only need to set the blocks of the gradient for all $i\in[N]$ and $j\notin{\cal N}_{i}$ as

[TABLE]

Based on this update law and with these choices for the stepsize, using the Lipschitz gradient property for the objective function we can easily prove that the sequence of value functions $f(S_{k},G_{k})$ is nonincreasing [18]:

[TABLE]

for some constant $\Delta>0$ for all $k\geq 0$. Therefore all the iterates remain in the compact sublevel set ${\cal L}(S_{0},G_{0})$. Moreover, we define gradient mapping of $f$ at $(S_{k},G_{k})$ [18]:

[TABLE]

which, according to Lemma 1, represents the natural measure of optimality for the constrained problem (28). Since $f$ is bounded from below by zero, then for any positive integer $k>0$ it is straightforward to prove from the previous descent inequality the following global convergence rate for the gradient mapping:

[TABLE]

where $f^{*}$ is the optimal value of problem (28). Moreover, under some mild assumptions, such as the Hessian of $f$ at a local minimum is positive definite and bounded, then starting sufficiently close to this local optimum the gradient iteration converges linearly to this solution [18]. Therefore, the speed of convergence of this iterative process depends on the starting point. Hence, to obtain a good initial stabilizing pair of matrices $(S_{0},G_{0})$ satisfying the network conditions $\mathcal{G}$, we can solve the structured SDP problem (1) and initialize with its solution. Note that for computing the gradient we first need to find the Gramians $\mathcal{M}_{k}$ and $\mathcal{W}_{k}$, solutions of the Lyapunov equations (V) and (30) in $(S_{k},G_{k})$. The solvability of $\mathcal{M}_{k}$ and $\mathcal{W}_{k}$, unique positive semidefinite solutions of (V) and (30), is implied by the stability of the error matrix $\mathcal{A}_{e}(S_{k},G_{k})$.

Hence, we consider the following algorithmic procedure:

Algorithm 1** ($H_{2}$ optimal network reduction algorithm)**

Let $(S_{0},G_{0})\in\mathcal{G}$ (e.g., $(S_{0},G_{0})$ - solution of SDP (1)). 2. 2.

Perform update (34) until $\|\Gamma_{\cal G}(S_{k},G_{k})\|\leq\epsilon$.

Note that the Gramians $\mathcal{M}_{k}$ and $\mathcal{W}_{k}$ are in general dense, which means the block diagonal assumption from Section IV is relaxed during each iteration of our projected gradient method. In fact, the block diagonal Gramian can be viewed as an intermediate step between a diagonal Gramian in [9] on positive linear systems and a full one as provided by our gradient method on general linear network systems.

VI Extensions

Note that our approach is general and flexible, allowing to tackle network systems with even more structure. For example, we can consider that coupling among subsystems is through both, states and inputs, i.e. we modify the dynamics (9) as:

[TABLE]

where $x_{i}$ and $u_{i}$ represent the states and inputs of the $i$th subsystem. Our optimization-based model reduction framework allows to easily incorporate the additional structured constraints, inherited from $B$ (i.e. $B_{ij}=0\;\forall i\in[N],j\notin{\cal N}_{i}^{u}$), on the matrix $G$. For example in the convex SDP (1) we just need to add the additional convex constraint $(Z_{22})_{ij}=0\;\forall i\in[N],j\notin{\cal N}_{i}^{u}$. Similarly, in the projected gradient method we just need to set the corresponding block components of the gradient $\nabla_{G}f$ to zero, i.e. $\nabla_{G_{ij}}f=0\;\forall i\in[N],j\notin{\cal N}_{i}^{u}$. Similarly, for positive networks we can easily incorporate positivity constraints on the system matrices of the reduced model.

VII Illustrative examples

In this section, we illustrate the efficiency of the proposed methods numerically. In particular, we compute and compare reduced order models for two network examples achieving (possibly) the minimum $H_{2}$ norm approximation.

VII-A Random positive network system

We consider a stable $12$th order linear positive network system as in (9) with matrices $A,B$ and $C$ described in Appendix -B, generated randomly in the interval $(-5,1)$ such that $A$ is stable and positive and satisfying the interconnection map from Figure 1. Here, we select $n_{i}=3$, $m=1$ and $p=1$ for all $N=4$ subsystems. Hence, $A\in\mathbb{R}^{12\times 12}$. We compute an $H_{2}$ optimal reduced order network, with $\nu_{i}=1$ for all $i=1:4$, with the same network structure. For the initialization of the algorithm we consider the solution $(S_{0},G_{0})$ of the structured SDP problem (1), with $L=[0\;0\;0\;1]$ fixed a priori. The output responses of $H_{2}$ SDP $\widehat{K}_{\text{sdp}}(s)$ and $H_{2}$ optimal obtained with projected gradient $\widehat{K}_{\text{grad}}(s)$ are displayed in Figure 2. We observe that the $H_{2}$ optimal output response, corresponding to projected gradient $\widehat{K}_{\text{grad}}(s)$, is almost identical with the response of $K(s)$ in the frequency range we considered $10^{-6}:10^{2}$.

In Table I, we compare the $H_{2}$-norm of the errors yielded by the proposed SDP and projected gradient methods versus the positive preserving projections method in [24] and the positivity preserving balanced truncation (BT) in [21]. Note that our projected gradient method performs best, that is the $H_{2}$ error of the $4$th order network decreased from $2.813$ (for SDP) to $5.075\cdot 10^{-3}$ (for projected gradient). Moreover, the constraints (iv) in Problem 1 are satisfied by the optimal model.

VII-B Power network system

Consider a power system split into $N$ control areas consisting of a generator and a load, with tie-lines providing interconnections between them, as described in [27]. For the area $i$ the simplified model is given by the differential equations:

[TABLE]

The interconnection to the control area $j\neq i$, is made through the tie-line with the power flow modeled by the equation:

[TABLE]

with $\Delta P_{\rm tie}^{ij}=-\Delta P_{\rm tie}^{ji}$. The notation $\Delta$ indicates the deviation from the steady-state of the variable, e.g., $\Delta\omega$ is the deviation of the angular speed from the nominal operating value. The variable are defined as follows:

• •

$\omega$ is the angular speed of the rotor,

• •

$M^{a}$ is the angular momentum,

• •

$D$ is the ratio between the percentage change in load and the percentage change in frequency,

• •

$P_{m}$ is the mechanical power,

• •

$T_{\rm CH}$ is the charging time constant,

• •

$P_{v}$ is the steam valve position,

• •

$P_{\rm ref}$ is the load reference setpoint,

• •

$R^{f}$ is the ratio between the percentage change in frequency and the percentage change in unit output,

• •

$T_{G}$ is the governor time constant,

• •

$P_{\rm tie}^{ij}$ is the tie-line power flow between areas $i$ and $j$,

• •

$T_{ij}$ is stiffness coefficient of the tie-line (between area $i$ and area $j$).

Considering a chain of $N$ areas, equations (VII-B) and (37) yield a $4N-1$th order model of the power network of the form (35) with the states:

[TABLE]

The input of each area is $u_{i}=\Delta P_{{\rm ref}_{i}}$, while as measured outputs we consider the angular speed deviations $\Delta\omega_{i}$ of each area $i=1:N$. Note that the matrix $A$ of the power network model has almost block bidiagonal form as seen in Figure 3.

The values of the parameters are chosen for each of the $N$ areas, randomly in the interval between the lowest and the highest possible physical values extracted from [27]:

• •

$D_{i}\in[0.255,80]$,

• •

$M_{a}^{i},T_{{\rm CH}_{i}}\in[1,5]$,

• •

$R^{f}_{i}\in[0.03,0.07]$,

• •

$T_{G_{i}}\in[4,10]$,

• •

$T_{ij}\in[1.5,2.5]$,

with $i,j=1:N$. With these values the resulting model of the power network is stable. For $N=4$, the matrices of the linear network model are given by (50) and by (51). We aim at selecting $N$ interpolation points, that is the order of each reduced subsystem is $\nu_{i}=1$ for all $i=1:N$, and accordingly the matrices $S,L,G\in\mathbb{R}^{N\times N}$. For solving the model reduction Problem 1 for the power network model from above with a chain of $N$ areas we use the projected gradient algorithm from Section V-A initialized with the SDP solution (1), for $L=I_{N}$ fixed a priori. This approach yields a $N$th order $H_{2}$ optimal approximation of the power network of the form (12) redering an optimal $H_{2}$ norm of the approximation error. Note that the projected gradient method preserves the (block) bidiagonal network topology of $S-GL$ as Figure 4 shows, as well as the stability of the network.

We plot in Figure 5 the optimal $H_{2}$ norm of the error approximation yielded by the projected gradient method for a number of $N$ areas ranging from $4:30$. The variation of the optimal $H_{2}$ norm is consistent with the increase in $N$, that is the $N$th order approximation of the $4N-1$ model is not necessarily more accurate as $N$ increases.

VIII Conclusions

In this paper we have studied the model order reduction for linear network systems. Using moment matching techniques, we have developed an optimization framework to compute parameterized reduced order stable models achieving moment matching, minimizing the $H_{2}$ norm of the error system, and preserving the structure of the to-be-reduced model of the network. For this, we have proposed two numerical procedures, based on SDP and projected gradient, for finding the (optimal) reduced order model of the network. Preliminary numerical simulations have confirmed the efficiency of our approach.

-A Note on partial minimization

The partial minimization is valid for any nonconvex program, i.e. given any function $F(x,y)$ (not necessarily convex) we always have:

[TABLE]

Let us briefly prove this statement. For this, take $x$ fixed and determine, as a function of $x$, that $y$ which minimizes $F(x,y)$ in the second variable. Let us denote this solution by $y^{*}(x)\in\arg\min_{y}F(x,y)$ for any fixed $x$. Further, let us define the partial function $f(x)=F(x,y^{*}(x))$ defined only in the variable $x$. Let $x^{*}\in\arg\min_{x}f(x)$. Then, the pair $(x^{*},y^{*}(x^{*}))$ is an optimal solution of the original problem since for all $(x,y)$ we have:

[TABLE]

-B Parameters of model from Section VII-A

For the positive network example from Section VII-A the numerical values of the $A$, $B$, $C$ matrices generated randomly in the interval $(-5,1)$ such that $A$ is a stable matrix and satisfying the interconnection map from Figure 1 are given in (50) and (51). Furthermore, in our numerical experiments we fix $L=[0\ \ 0\ \ 0\ \ 1]$.

[TABLE]