Denis Nikitin

Carlos Canudas-De-Wit

Paolo Frasca

Control of Average and Deviation in Large-Scale Linear Networks

Keywords: Networks Control, Large-Scale Control, Extremum Seeking

This paper deals with the problem of controlling the average state of a large-scale linear network to a constant reference value. We design an output-feedback controller such that no information about state vector or system matrices is needed. For this controller to have arbitrary positive gains, it is sufficient that only a sign condition on system matrices should be satisfied. To assure that the states of the network are close to the average state, the problem of deviation minimization is solved in addition, using a novel extremum seeking algorithm.

I. INTRODUCTION

T HE control of large-scale dynamical networks is a chal- lenging and important question in modern control theory. Its importance originates from its multiple applications to vital systems in technology, nature and society, such as urban traffic networks [START_REF] Bianchin | Gramian-based optimization for the analysis and control of traffic networks[END_REF], social networks [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: Convergence conjectures and theorems[END_REF], fleets of robots [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: Algorithms and theory[END_REF] and power networks [START_REF] Dorfler | Synchronization and transient stability in power networks and nonuniform kuramoto oscillators[END_REF]. Its difficulty originates from the large dimensionality of these real-world networks, where the number of states can reach millions. These large networks challenge the scalability of control methods from several points of view. First, the computation of traditional control algorithms becomes too expensive. Second, the structure and the detailed dynamics of the network may not be fully known. Third, the number of actuators and sensors is often much lower than the number of nodes, so that state feedback is not possible: see for instance biological neural networks, where only an average neuronal activity is measured by electrodes. Another difficulty is that the amount of energy needed to control all nodes of the network can grow unbounded as the number of nodes increases: the growth is actually exponential for some network structures [START_REF] Yan | Controlling complex networks: How much energy is needed?[END_REF], [START_REF] Liu | Controllability of complex networks[END_REF], [START_REF] Cowan | Nodal dynamics, not degree distributions, determine the structural controllability of complex networks[END_REF]. In assessing energyefficient control problems, it has been shown that several aspects such as network centrality [START_REF] Bof | On the role of network centrality in the controllability of complex networks[END_REF], [START_REF] Lindmark | A driver node selection strategy for minimizing the control energy in complex networks[END_REF] and number of control nodes [START_REF] Yan | Spectrum of controlling and observing complex networks[END_REF] play a fundamental role, but still energy analysis of large-scale network controllability is a very hard task to perform [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF].

The control community has often approached the issue of network control by looking for distributed control algorithms, in which the control is applied locally at all nodes and uses only local information. Instead, in this paper we choose to work in a centralized setting, where an external operator has limited information about the network and limited access to few nodes for sensing or actuation purposes. In view of these limitations, the operator shall aim at controlling some aggregate function of the network state, rather than controlling D. Nikitin, C. Canudas-de-Wit and P. Frasca are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. all of its individual nodes. A natural choice for such aggregate function is the average of the node states, which has indeed been defined as a control objective in some prior work that was motivated by opinion dynamics in social networks [START_REF] Vassio | Message passing optimization of harmonic influence centrality[END_REF], [START_REF] Rossi | On the convergence of message passing computation of harmonic influence in social networks[END_REF]. More broadly, the control of a generic output of a largescale network was studied in [START_REF] Klickstein | Energy scaling of targeted optimal control of complex networks[END_REF], [START_REF] Wittmuess | Optimal actuator design for optimal output controllability[END_REF], [START_REF] Klickstein | Locally optimal control of complex networks[END_REF], [START_REF] Commault | Functional target controllability of networks: structural properties and efficient algorithms[END_REF], [START_REF] Casadei | Controllability of large-scale networks: An output controllability approach[END_REF]: the energy required to control aggregated outputs instead of all states is much less.

In this paper, we focus on the problem of controlling the average state of the network, together with the concurrent minimization of its deviation. The average state of the network is defined as an average over all node states of the network, while the squared deviation is defined as an average over all squared differences between node states and the average state. While using a controller for the average state, it is natural to desire the system states to be close to the average: this behavior can be obtained by minimizing their squared deviation. In opposition to previous work, it is assumed that the only values that are measured and regulated are the values of the system outputs, i.e. the average state and the squared deviation. Moreover, we make sure that the system model is not used in the controller. Thus, the controller directly utilizes only system outputs and reference point and the equilibrium of internal states is never computed explicitly. Not using the system model circumvents all issues about computational complexity and uncertainties that affect large networks.

In our problem formulation we assume that the network operator has knowledge of the average and the squared deviation values. There are many physical examples of systems where the average and the squared deviation can be measured without measuring the states of the nodes. Here we briefly mention four examples:

Urban traffic networks. Consider a network of roads in a city, where the state of each node is the number of cars on the corresponding road section. The total number of cars in the city can be estimated either directly or indirectly. A direct estimation of the number of cars can be performed by visionbased methods, by processing images taken from satellites, as in [START_REF] Eslami | Automatic traffic monitoring from satellite images using artificial immune system[END_REF], [START_REF] Palubinskas | Model based traffic congestion detection in optical remote sensing imagery[END_REF]. Although every car is counted independently, the estimation error is defined as a discrepancy in the overall number of cars, therefore these methods effectively reconstruct the total number of cars. An indirect estimation can be based on the vehicle emissions: combustion engines produce CO 2 , which then goes into the atmosphere. The polluted atmosphere changes its reflection properties based on the amount of CO 2 , thus this amount can be measured using infrared sensors mounted on satellites [START_REF] Boynard | First simultaneous space measurements of atmospheric pollutants in the boundary layer from iasi: A case study in the north china plain[END_REF]. Therefore, the number of cars can be reconstructed from the satellite measurements. The total number of cars divided by the number of road sections equals to the average state of the network.

Biological neural networks. A widely known method of monitoring the brain activity is the electroencephalography, with electrodes placed usually along the sculp of the person being monitored. Each electrode measures voltage fluctuations of group of neurons under the surface, therefore it is directly related to the average of individual states of neurons, which obviously cannot be measured independently [START_REF] Van Veen | Localization of brain electrical activity via linearly constrained minimum variance spatial filtering[END_REF].

Dynamics of gas. Every gas consists of a huge number of particles colliding with each other, therefore it can be seen as a dynamical network with neighbouring particles whose interaction depends on their velocities. Thus, we can define the states being the velocities of each individual particle. The gas temperature can be easily measured, but at the same time it corresponds to the internal kinetic energy:

E k = 3 2 k B T = 1 2 mv 2 rms .
Here k B is the Boltzmann constant, T is the temperature and m is the mass of one particle. The variable v 2 rms represents the mean squared deviation of the velocities of particles with respect to the flow velocity. The flow velocity itself is a "wind speed", which represents the average state of the system and can be also directly measured.

Density of a fluid. Fluids also consist of a huge number of particles, and one way to write the dynamical model of a fluid is to consider a space partitioned into individual cells with states being defined as the densities of the fluid inside each cell. In this case, the average state would be the average density in the system: this density can indeed be measured for cryogenic fluids by measuring permittivity by a technique called electrical capacitance tomography [START_REF] Hunt | High-speed density measurement for lng and other cryogenic fluids using electrical capacitance tomography[END_REF].

Another relevant setup is the scale-free control approach to large-scale networks [START_REF] Canudas-De | Erc: Scale-free control for complex physical network systems[END_REF] illustrated in Fig. 1: in this approach, the goal is to control the average state and the deviation of the "hub" regions and the control is applied to the boundaries of the hubs. Output controllability of a scale-free network is studied in [START_REF] Casadei | Model reduction based approximation of the output controllability gramian in large-scale networks[END_REF], the dual problem of reconstruction of an averaged state is solved in [START_REF] Niazi | Average observability of large-scale network systems[END_REF], [START_REF]Average state estimation in large-scale clustered network systems[END_REF], and methods to reduce a large network to a scale-free one, through the suitable aggregation of hub regions, are presented in [START_REF] Martin | Large-scale network reduction towards scale-free structure[END_REF], [START_REF] Martin | Subgraph detection for average detectability of lti systems[END_REF].

The first part of our work is devoted to the problem of controlling the average state of the linear network, where three stability results are presented. Theorem 1 provides conditions on the integral controller gains for the stability of the closedloop system, and Theorem 2 simplifies the conditions under the assumption that the system is positive. Then, Theorem 3 is our main contribution which gives a simple sufficient sign condition on the system matrices which guarantees stability of any positive integral controller for controlling the system output to a constant reference point without knowledge of the system matrices.

Most of our results regarding the output regulation problem of a large linear network system are presented under the assumption that the system is stable and positive (that is, the system matrix has positive elements outside the main diagonal). Network systems with stable dynamics and positive edge weights belong to this class. More generally, positive systems are an important class of systems for which the synthesis of large-scale control algorithms can be greatly simplified. Their impulse response is bounded by their static gain [START_REF] Rantzer | Distributed control of positive systems[END_REF], optimal [START_REF]Scalable control of positive systems[END_REF] and robust [START_REF] Briat | Robust stability and stabilization of uncertain linear positive systems via integral linear constraints: L1-gain and l∞-gain characterization[END_REF] feedback control laws can be easily designed using linear programming, and the state feedback output regulation problem can be explicitly solved [START_REF] Nogueira | Contributions to positive systems control: mathematical methods and applications[END_REF]. From the passivity analysis in the classical control theory it is known that the feedback interconnection between a linear operator with an integral controller is stable irrespective of the gain (has an infinite gain margin) if the linear operator is strictly positive real (SPR) [START_REF] Sepulchre | Constructive nonlinear control[END_REF], [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF]. From this point of view our analysis provides a new simple sufficient condition for the positive system to be SPR, which is summarized in Theorem 4.

The second part of the work focuses on the deviation minimization problem, when the system should be driven to the particular average state while the control inputs should be balanced in such a way that the squared deviation of the states takes the smallest possible value. To solve this problem we use an extremum seeking scheme [START_REF] Ariyur | Real time optimization by extremum seeking control[END_REF], [START_REF] Tan | Extremum seeking from 1922 to 2010[END_REF] which is an adaptive model-free algorithm for the minimization/maximization of a nonlinear steady-state output characteristic. We augment this algorithm with an additional subsystem such that both tasks are accomplished simultaneously: the average is driven to the particular value while the squared deviation is minimized. Theorem 5 proves that the system approaches any small neighbourhood of the optimum state provided the gains of the controller are small enough.

An incomplete account of this research appeared as [START_REF] Nikitin | Boundary control for output regulation in scale-free positive networks[END_REF]: the latter short version is much narrower than the present version as it focuses on the statement of Theorem 3 and on the interpretation of the sufficient condition therein (without providing Theorems 1, 2, 4 or 5, nor the detailed proof of Theorem 3).

Notation: Along this work several types of vector inequalities are used. If x ∈ R n , then x 0 means x i 0 ∀i ∈ {1..n}, x > 0 means x i 0 ∀i ∈ {1..n} and there exists j ∈ {1..n} : x j > 0, and x 0 means x i > 0 ∀i ∈ {1..n}. Matrix I denotes the identity matrix of an appropriate size, while column vector of ones is denoted as 1. The set of complex numbers with non-negative real parts is C + . We use the standard vector norm x =

√

x T x for x ∈ R n throughout the paper.

II. PROBLEM FORMULATION

x 1 : boundary nodes We start the problem formulation with an example: assume the system we need to control is the network given by the graph G = (V , E ), where V is the set of vertices and E is the set of edges. The number of vertices |V | is denoted by n.

On each node v i ∈ V the state x i is defined. Each edge e ∈ E , where e = {v i , v j }, corresponds to the flow between nodes v i and v j . Matrix A ∈ R n×n represents flow ratio. The set of nodes V is split into two parts V 1 and V 2 with state vectors x 1 and x 2 respectively (see Fig. 2). The set V 1 consists of the nodes which are directly controlled by the control action u. We call these nodes "boundary". The set V 2 is a set of uncontrolled nodes, which we call "inner nodes". The average state y = 1 T x/n and the squared deviation V = x 2 /ny 2 are measured. Thus the network depicted on Fig. 2 can be viewed as one particular hub from the scale-free network on Fig. 1.

The evolution of the states x, the average state y and the squared deviation V is given by the following linear timeinvariant system

             ẋ1 = A 11 x 1 + A 12 x 2 + u, ẋ2 = A 21 x 1 + A 22 x 2 , y = 1 n 1 T x, V = x 2 /n -y 2 .
(

) 1
Most of real-world networks are internally stable, so we further assume A being stable. Also in Theorems 2 and 3 we assume A is a Metzler matrix (which means its offdiagonal elements are non-negative), which means all edges have positive weights. Such choice of system matrix together with the fact that B > 0 and C > 0 means that the system (1) belongs to the class of positive systems.

It is useful to analyse more general case than (1), with general stable matrix A ∈ R n×n , C ∈ R 1×n , B ∈ R n×k , and symmetric positive semi-definite matrix P ∈ R n×n defining the quadratic output:

     ẋ = Ax + Bu, y = Cx, V = x T Px. (2)
System (1) can be written in form of (2) using

A = A 11 A 12 A 21 A 22 , B = I 0 , C = 1 n 1 T , P = 1 n I - 1 n 2 11 T
In general, the control goal is to stabilize average state y over the whole network to some desired constant state y d without the explicit knowledge of system matrices. It is assumed that the number of states is too large that it is impossible to use full-state feedback or to use matrix A explicitly.

In the following sections two problems will be solved: 1) Control of average: find the control law u = u(y) for the system (2) such that lim t→∞ y(t) = y d .

2) Control of average and deviation: Find the control law u = u(y,V ) for the system (2) such that lim t→∞ y(t) = y d , and V is minimized. 1Section III is devoted to the first problem of average control, while Section IV shows the solution for the second problem of simultaneous control of average and deviation.

III. AVERAGE STATE REGULATION ONLY

In this section, we solve Problem 1 of controlling the average. We first describe a general controller structure and then prove stability conditions for low and high gains. Further we discuss these conditions and present examples of the controller performance.

A. Controller Structure

Define transfer function of the system (2) equation is

W s (s) = C(sI -A) -1 B, (3)
y(s) = W s (s)W c (s)e(s), (4)
or, solving for y,

y(s) = W s (s)W c (s) 1 +W s (s)W c (s) y d . (5)
Define closed-loop transfer function

W (s) = W s (s)W c (s) 1 +W s (s)W c (s) . (6)
In the following we investigate what properties do W s and W c have and what properties should W have in order to be stable. Both W s and W c are vector-valued, because controller input e and system output y are scalars, while u which is controller output and system input is a vector, u ∈ R k . Let us look at the i-th component of W s and W c , where i ∈ {1..k}, and define polynomials α(s), β i (s), δ (s), γ i (s) such that

W s (s) i = β i (s) α(s) and W c (s) i = γ i (s) δ (s) . (7)
It is obvious that α(s) is a polynomial of degree n. Moreover, our system is strictly stable, non-casual and matrix D (which directly connects y and u) is zero, thus

deg β i (s) < deg α(s) = n and α(s) = 0 ∀s ∈ C + . ( 8
)
We can choose α(s) and β i (s) such that α(s) ∈ R for s ∈ R and α(s) > 0 for s 0.

Then, the controller W c should also be stable and non-casual, which means

deg γ i (s) deg δ (s) and δ (s) = 0 ∀s ∈ C + . (9)
Again, it is possible to choose δ (s) and γ i (s) such that δ (s) ∈ R for s ∈ R and δ (s) > 0 for s > 0. Now we can rewrite W (s) in terms of polynomials:

W (s) = ∑ i β i (s)γ i (s) α(s)δ (s) + ∑ i β i (s)γ i (s) . (10)
The closed-loop transfer function W (s) should have the following property: for the constant input y d it should give the same output y, thus W (0) = 1. This means that α(0)δ (0) = 0, which is possible only if δ (0) = 0, so δ (s) cannot contain free term. The simplest possible controller that satisfies this necessary condition is the integral controller given by

W c (s) i = κ γ i s , (11)
where γ ∈ R k is the vector of gains, defining relative control force applied to different actuated nodes, and κ is the overall gain. The following sections will be devoted to the integral controller and its properties.

B. Stability of the integral controller

Assume we apply the integral controller (11) to the system (2). The closed-loop system may be unstable, and in general in order to prevent this one needs to carefully choose controller gains κ and γ in [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF].

Theorem 1. The system (2) with applied integral controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF] is asymptotically stable if -CA -1 Bγ > 0 and κ ∈ (0, κ * ) for some small κ * ∈ R.

Proof of Theorem 1. Applying the integral controller, the transfer function of the closed-loop system is given by

W (s) = κ ∑ i β i (s)γ i α(s)s + κ ∑ i β i (s)γ i . ( 12
)
For the stability of the closed-loop system W (s) should have no poles on the right-hand side of the complex plane C + . Decompose the denominator:

α(s)s   1 + κ ∑ i β i (s)γ i α(s)s   = 0. Denote Q(s) = ∑ i β i (s)γ i / (α(s)s).
Any point such that α(s) = 0 or s = 0 leads to W (s) = 1, thus the poles of the transfer function can arise only in the roots of 1 + κQ(s). We will prove that there exists κ * such that

∀κ ∈ (0, κ * ) : ∀s ∈ C + \ {0} Re{1 + κQ(s)} > 0. (13) Choose ε, R ∈ R such that ε < |λ i (A)| and R > |λ i (A)| for all i ∈ {1..n}.
Thus all the roots of α(s) lie in a ring between ε and R in a left half-plane. We split complex right half-plane C + into three parts:

H + 0,ε = {s : Re s 0, |s| < ε} H + ε,R = {s : Re s 0, |s| ε, |s| R} H + R,∞ = {s : Re s 0, |s| > R} First we analyse H + 0,ε . Function Q(s)
has a pole at zero, thus it can be written using Laurent series with coefficients Q n :

Q(s) = Q -1 s + ∞ ∑ n=0 Q n s n = Q -1 s + P(s), where P(s) is an analytic function. The residual Q -1 = ∑ i β i (0)γ i /α(0) = -CA -1 Bγ > 0, thus Re Q -1 s 0 ∀s ∈ H + 0,ε \ {0}, while P(s) is analytic in C and thus has a minimum in H + 0,ε . Next we analyse H + R,∞ . This set is contained into a set H R,∞ = {s : |s| > R}. If R is big enough, Q(s) is analytic in H R,∞
, but it vanishes at infinity, therefore by a maximum modulus principle Q(s) is bounded from below in H R,∞ by values on its boundary, and consequently it is bounded in H + R,∞ . Finally, set H + ε,R is compact and does not contain zeros or roots of α(s). Therefore Q(s) is analytic in it and thus bounded. We obtained that Re Q(s) is bounded from below in C + \ {0}. Denoting this bound as Q in f , we see that choosing [START_REF] Rossi | On the convergence of message passing computation of harmonic influence in social networks[END_REF] and therefore proves the theorem.

κ * = -1/Q in f in case Q in f < 0 or κ * = +∞ in case Q in f 0 as- sures satisfaction of
Theorem 2. The system (2) with applied integral controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF] is asymptotically stable if the system (2) is positive, γ > 0, -CA -1 Bγ = 0 and κ ∈ (0, κ * ) for some small κ * ∈ R.

Proof of Theorem 2. Positivity of the system (2) means that all elements of matrices B and C are greater or equal than zero, and matrix A is a Metzler matrix. Now we introduce the notion of M-matrix:

Definition 1 (M-matrix, [START_REF] Plemmons | M-matrix characterizations. i-nonsingular mmatrices[END_REF]). An n × n matrix M that can be expressed in the form M = αI -L, where L = (l i j ) with l i j 0, 1 i, j n, and α λ (L) where λ (L) is the maximum of the moduli of the eigenvalues of L, is called an M-matrix.

From this definition it follows immediately that a negative of a Metzler stable matrix is an M-Matrix. The main property of any M-matrix M is that its inverse M -1 is a positive matrix, thus (M -1 ) i j 0 for all i, j. [START_REF] Fan | Topological proofs for certain theorems on matrices with nonnegative elements[END_REF] Matrix -A is an M-Matrix which means that -A -1 has its all elements nonnegative, therefore -CA -1 B is a positive vector. By the theorem statement γ > 0, and having -CA -1 Bγ = 0 means -CA -1 Bγ > 0 and, applying Theorem 1, this leads to the asymptotic stability of the closed-loop system for small enough κ.

By Theorem 2, if the system is positive it is enough to choose γ = 1 (or any γ 0, just to satisfy -CA -1 Bγ = 0), and then pick up small enough overall gain κ.

C. Controller with arbitrary large gains

It appears although that there exists simple criteria on the system matrices which says whether the closed-loop system will converge irrespectively of controller gain values, provided that they are positive. This result is one of the main contributions of this paper, and it is formulated as follows:

Theorem 3. The system (2) with applied integral controller (11) is asymptotically stable for arbitrary large positive controller gains κ and γ if the system (2) is positive, CA 2 > 0 and CA 2 Bγ > 0.

This result means that, irrespective of the gains, the integral controller will preserve stability for a very large class of systems. One of the important types of large-scale networks for which Theorem 3 is satisfied is a general consensus network (e.g. for social interactions), see the example below.

Example (Damped consensus). Assume the system (2) is given by matrices A = -L -αI, where L is a Laplacian matrix of some network with n nodes, α > 0 means additional damping to the system to preserve stability, and C = 1 T /n represents average state of the network. Then A is a Metzler stable matrix, and C is the eigenvector of A with corresponding eigenvalue -α, thus CA 2 = α 2 C > 0. Then any controller with positive gains κ and γ will lead to the convergence, provided Bγ > 0.

One should notice that the condition CA 2 Bγ > 0 on the control matrix B is very non-restrictive, because by choosing appropriate vector gain γ it is always possible to make Bγ > 0, and hence, provided CA 2 > 0 and CA 2 B = 0, we will have CA 2 Bγ > 0. The reason for this is the fact that the regulation variable is a single scalar output. Proposition 1. Condition CA 2 Bγ > 0 is a sufficient condition for the output controllability of the system (2).

Proof. Indeed, Kalman rank test for the output controllability of (2) can be written as

rank{C B AB A 2 B ... A n-1 B } = 1,
and by CA 2 Bγ > 0 we have CA 2 B = 0, which means that the rank test is satisfied.

Note that the analogue of this Proposition can be proven for Theorems 1 and 2.

Corollary 1. The positive system (2) with CA 2 0 is asymptotically stable for any integral controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF] with positive gains applied to any single boundary node. Therefore, it is enough to control only one node.

Before proving the third theorem we need to state three technical lemmas.

Lemma 1. Suppose we have a matrix M = M + ibI, which is a complex matrix with real part M and imaginary part bI, with b ∈ R and I an identity matrix. Assume M being invertible and having no eigenvalues on the imaginary axis. Denote L = M -1 = L + i L. Then the real part of L is given by

Re L = L = (M + b 2 M -1 ) -1 . ( 14
)
Proof. See Appendix A.

Lemma 2. Let M be an M-matrix. Let C be a row-vector such that CM 2 > 0. Then

C(M + tM -1 ) -1 > 0 ( 15
)
for any t 0.

Proof. See Appendix B.

Lemma 3. Let M be an M-matrix. Let C be a row-vector such that CM 2 > 0 and CM 2 Bγ > 0. Then

C(M + tM -1 ) -1 Bγ > 0 ( 16
)
for any t 0.

Proof. See Appendix C.

Proof of Theorem 3. Applying the integral controller and multiplying nominator and denominator by s, the transfer function of the closed-loop system is given by

W (s) = κC(sI -A) -1 Bγ s + κC(sI -A) -1 Bγ . ( 17
)
It is sufficient to show that the real part of denominator is strictly greater than zero in the right half-plane. Since Re s > 0 in the right half-plane, it is enough to show that Re κC(sI -A) -1 Bγ > 0.

Denote Re s = α and Im s = β , so matrix (sI A) is an M-matrix and matrix M is an M-matrix too. Moreover, condition CA 2 > 0 implies CM 2 > 0 and CA 2 Bγ > 0 implies CM 2 Bγ > 0. Applying Lemma 1 we conclude that Re κC(M

-A) -1 = ((αI -A) + iβ I) -1 . Denote M = αI -A. Matrix A is a Metzler stable matrix, thus (-
+ iβ I) -1 Bγ = κC(M + β 2 M -1 ) -1 Bγ. ( 18
)
By Lemma 3 C(M + tM -1 ) -1 Bγ > 0 for any t 0, and assuming κ > 0 we trivially obtain a sufficient condition on positivity of the real part of the denominator. Proof. The SPR transfer function should by definition satisfy Re{C(sI -A) -1 B} > 0 for Re S > 0, see [START_REF] Sepulchre | Constructive nonlinear control[END_REF], [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF]. Therefore the proof of this theorem follows the same steps as the second part of the proof of Theorem 3.

From the point of view of linear systems theory, Theorem 4 is the main result of our paper. Indeed, such a simple condition for passivity for positive systems appears in literature for the first time.

Moreover, it appears that Theorem 3 is a direct consequence of Theorem 4, as the following reasoning shows. Assume we fix an input gain vector γ and define a new controller output α such that u = γα. Then the system (2) becomes SISO with respect to input variable α. Define H 1 (s) = κ/s and H 2 (s) = C(sI -A) -1 Bγ. It is possible to construct control loop with feedback interconnection as depicted on Fig. 4. The closed-loop system input is defined as y d and the system output is α. It is known that for L 2 stability of the system with feedback interconnection it is sufficient that the transfer function of one of the blocks is positive real (PR, which is equivalent to passivity) and another is strictly positive real (SPR) [START_REF] Sepulchre | Constructive nonlinear control[END_REF]. Passivity of an integral controller H 1 (s) is obvious, and Theorem 4 is used to prove that H 2 (s) is SPR. Therefore the closed-loop system is L 2 stable. Now, it remains to prove that y → y d , which is obvious if one recalls that an output of a stable system with constant input converges to a constant value, thus for any constant y d there exists α * such that α → α * . But convergence of an output of an integral controller means that its input converges to zero, which reads as e → 0, which is exactly y → y d .

D. Practical use of the theorems

Theorems 1-3 presented in the previous sections provide the same result, the stability of the closed-loop system (2) with controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF]. However, they differ in their assumptions. During the derivation of the controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF] we assumed that the system matrices are not known. However, usually one has a knowledge about some general properties of the system, such as positivity. These properties can sometimes be induced from the nature of the problem itself and do not rely on the particular topology. Therefore the results presented in our work can be used to analyse the stability based on these properties.

Theorem 1 requires -CA -1 Bγ > 0 for the integral controller to be stable for κ ∈ (0, κ * ). This scalar condition essentially means that the direction of adaptation of integral controller forms an acute angle with the zero frequency gain of the system. In practice one usually knows the direction of the zero frequency gain. At worst, it is enough to change the sign of γ once.

Theorem 2 exploits positivity of the system: the zero frequency gain of the positive system is positive. Therefore it is enough to use positive gains for the integral controller, and the condition -CA -1 Bγ > 0 can be loosened just to -CA -1 B = 0. However the gain κ still should satisfy κ ∈ (0, κ * ).

In Theorem 3 the small-gain condition κ ∈ (0, κ * ) is removed at the cost of adding the vector inequality CA 2 > 0. This inequality can be used to determine stability without knowing the particular matrices for some classes of systems, such as damped consensus, see the example after Theorem 3. For other systems, this condition should be interpreted as a constraint on the system parameters.

In the remainder of this section we will analyse the condition CA 2 > 0 more closely. Namely, first of all we will prove that this condition cannot be relaxed, since weaker conditions would not assure the stability for all gains. Then, we will provide some graph-theoretical intuition and rewrite the condition in terms of quadratic constraint on the node selfdampings.

If A is a Metzler stable matrix, all elements of A -1 are nonpositive. Multiplication of a positive vector by a matrix with nonpositive elements renders negative vector, therefore right multiplying the condition CA 2 > 0 by A -1 one obtains CA < 0, and the same argument provides C > 0. The condition CA 2 > 0 is new and it is used in the Lemmas 2 and 3 (substituting M = αI -A as in the proof of the theorem). When one looks at the statement of the Lemma 2, one might think that it would be enough to require the less restrictive condition CA < 0 (This condition is obtained from the statement of the Lemma 2 by letting t → +∞) and has been proposed for the full state static feedback output control of positive systems by [START_REF] Nogueira | Contributions to positive systems control: mathematical methods and applications[END_REF]).

However, let us show that condition CA 2 > 0 is significant and CA < 0 is not sufficient. An example of the positive system with CA < 0 but CA 2 > 0 would be

A =   -1 0 0 1 -1 0 0 1 -1   , B =   1 0 0   , C = 1 T . ( 19
)
For this system CA = (0, 0, -1), but CA 2 = (0, -1, 1). We can then show that this system is not SPR. To check this, by definition we take a pole in the complex right half-plane s = 0.01 + 2i, which results in

ReC(sI -A) -1 B = -0.0047. ( 20
)
Since the transfer function value is negative, the system is not positive real and thus it is not passive. Moreover, there exists an integral controller which makes this system unstable, for example one with the control vector γ = 1 (since only one node can be controlled) and the gain κ = 3 (although with κ = 2 the system is still stable). This confirms our understanding that the novel CA 2 > 0 condition is meant to ensure stability using any arbitrary boundary node and arbitrary positive gain κ.

Going deeper to understand topological properties of the condition CA 2 > 0, we first start with more intuitive one, CA < 0, which is implied by CA 2 > 0.

Define matrices D and E such that A = E -D, with D being diagonal and E having all diagonal elements zero. Thus both D and E have all their entries positive. Matrix E can be viewed as adjacency matrix of the network, with element E i j meaning influence of node v j on node v i . Matrix D consists of the selfdamping powers on the diagonal. Therefore condition CA < 0 reads as CD > CE. This condition states some kind of diagonal dominance in the network.

Assume some C i = 0. Then (CD) i = 0 because D is diagonal. Thus (CE) i should be also zero, which means that for every index j either C j = 0 or E ji = 0.

Corollary 2. If node v i is not included in the aggregated output (C i = 0), then its reachable set should not be included either.

• For a strongly connected graph this means that all nodes should be included in the aggregated output. • If the network is divided into "boundary" nodes and "inner" nodes, and the goal is to control the average of the inner nodes, then at least one of the boundary nodes should also be included into the average.

In the same manner it is possible to see this condition as a lower bound on the damping of each node:

D ii ∑ j C j E ji /C i .
Thus the bigger is the influence of the node's neighbours in the output, the bigger should be the node's damping.

We can use the same decomposition A = E -D in order to understand the condition CA 2 > 0 and conclude that

CE 2 +CD 2 > C(ED + DE). (21)
Being a quadratic inequality, this condition bounds damping of each node from above and below with respect to dampings of other nodes. We will see examples in the following section.

E. Examples for average control

Here we present three examples of networks, namely a star, a line and an Erdős-Rényi graph, and analyse the condition CA 2 > 0 for them.

1) Network with star topology: To begin with we choose network with star topology with one central node and n leafs, average state of which we want to control. Let nodes 1...n be the leafs and node n + 1 be the center. Assume the center and the first leaf belong to the boundary node set and thus can be controlled (see Fig. 5). Dynamics of this network can be written as system (2) with matrices

A =     -1 -β 0 • • • 1 0 -1 -β • • • 1 . . . . . . . . . . . . 1 1 • • • -n -α     , B =     0 1 0 0 . . . . . . 1 0     , (22)
C = 1 T /(n + 1). ( 23) Such choice of system matrices corresponds to the undirected network with star topology and damping α > 0 for central node and β > 0 for all other nodes. The choice of B explores both cases of controlling leaf and center. It allows for maximum generality, moreover the controllability is guaranteed by Corollary 1. Integral controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF] with γ = (1, 0) T would correspond to the control applied only to the center, and controller with γ = (0, 1) T would correspond to the control of the first leaf.

x n+1 x 1 x 2 x n u 1 u 2 1 1 1 -α -β -β -β
Calculating CA and CA 2 gives

CA = -β -β • • • -α /(n + 1) < 0, (CA 2 ) 1...n = β 2 + (β -α) /(n + 1), (CA 2 ) n+1 = α 2 + n(α -β ) /(n + 1). (24)
CA 2 > 0 means then α 2 +n(α -β ) 0 and β 2 +(β -α) 0 with at least one of these inequalities being strict. Solving this for damping of leaf nodes we obtain

α + 1 4 - 1 2 β α + α 2 n , (25)
thus β is bounded from both sides with respect to α. Moreover, as n → ∞, we obtain limit inequality β α, which means that damping for leafs should be lower than damping for the center. Simulation results for both cases, γ = (1, 0) T and γ = (0, 1) T , and for n = 20 leafs are given on Fig. 6, with dampings α = 2, β = 1.1, desired output value y d = 5 and integral controller gain κ = 12. On Fig. 6a it is clearly seen that controlling the central node and controlling the leaf has almost the same effect on the output y.

2) Line network: Now we explore an example of a directed line network with n nodes. This network is depicted on Fig. 7. As usual, we are interested in controlling average state of the network, and it is assumed that we can control only the input node x 1 of the system. System matrices for n nodes are given as follows: x 1

A =      -1 -β 0 • • • 0 1 -1 -β • • • 0 . . . . . . . . . . . . 0 0 • • • -α      , B =      1 0 . . . 0      , (26)
C = 1 T /n, γ = 1. ( 27
) (a) (b)
x 2 x 3 x n u 1 1 -β -β -β -α
CA = -β -β • • • -α /n < 0, (CA 2 ) 1..n-2 = β 2 /n > 0, (CA 2 ) n-1 = (β 2 + β -α)/n, (CA 2 ) n = (α 2 + α -β )/n, (28)
CA 2 > 0 means then α 2 +(α -β ) 0 and β 2 +(β -α) 0. Solving this for damping of leaf nodes we obtain

α + 1 4 - 1 2 β α + α 2 , ( 29
)
thus β is bounded from both sides with respect to α.

In order to validate the conclusions about this example, we take directed line networks with 4 and 100 nodes and check whether they are stable or unstable for different κ.

Fix α = 0.2, therefore for condition CA 2 > 0 to hold one needs √ 0.45 -0.5 β 0.24. On Fig. 8 simulation results are shown for κ = 12, y d = 5 and for two values of β , the first, β = 0.2, satisfies the condition, and the second β = 0.002 does not. In the case β = 0.2 and n = 100 it is very interesting to see what are limit values of the state variables x. It appears that they decrease exponentially starting from the controlled node x 1 , while preserving their average equal to y d . This is due to the fact that in the steady state all node states except the first one and the last one should satisfy relation x i-1 -(1+β )x i = 0.

3) Random Erdős-Rényi graph: Here we present a simulation results for an integral controller for random Erdős-Rényi graph with n = 4000 nodes and probability of creating an edge p = 0.01. Vector C = 1 T /n represents the average, and the system matrix A is the negative of a Laplacian of this ER graph with an additional random damping on every node taken from the uniform distribution U(6, 7) such that CA 2 > 0. Matrix B is chosen to be a random vector of zeros and ones with equal probability. With such setup for any κ > 0 the system converges to the desired output reference y d = 5, see Fig.

IV. DEVIATION MINIMIZATION

A. Problem motivation

The previous sections showed that the control of average can be performed rather easily and without any knowledge of the system. But in some cases controlling the average with arbitrary control direction γ can lead to a poor performance: although the average state y → y d , states themselves can be very far from y d . As an example one can look at Fig. 8c-d, where the spread of states x is shown. Although the average value is 5, most of the states in the steady state are almost zero, while some states are much larger, around 80. This dispersion between states is captured by the squared deviation V . The smaller it is, the closer are the states to the average value.

Therefore, it makes sense to find the control law u = u(y,V ) for the system (2), which solves simultaneously two problems: assures lim t→∞ y(t) = y d and minimizes the squared deviation V .

Preliminary, let us make the following observation. Controlling the average to the desired value in the steady state means that the system should satisfy one-dimensional constraint -CA -1 Bu * = y d , thus if the dimension of the steadystate control vector u * is k, there are still k -1 degrees of freedom left for optimizing the control direction in sense of minimization of the squared deviation.

B. Explicit solution

Let us assume that the desired steady state is reached, and try to find it. Denote x * and u * as the state vector and the control vector respectively in the steady state. Also denote the steady-state squared deviation as V * . Then the equations for the steady state, obtained from the system (2), assuming

y → y d , are      0 = Ax * + Bu * , y d = Cx * , V * = x * T Px * . (30)
Our problem can be seen as a linear constrained quadratic minimization problem:

minimize V * = x * T Px * , subject to Ax * + Bu * = 0, Cx * = y d . (31)
In comparison to the standard linear-quadratic regulator, note that problem [START_REF]Scalable control of positive systems[END_REF] is formulated for the steady state, thus there is no more dynamics in it, as well as no optimization of the transient process.

Assume for a moment that all the system matrices are known. Using the fact that the matrix A is stable, we can take the inverse and thus obtain the steady state vector x * = -A -1 Bu * . Denoting S = B T A -T PA -1 B and η = -B T A -T C T , we can write the minimization problem [START_REF]Scalable control of positive systems[END_REF] in terms of u * :

minimize V * = u * T Su * , subject to η T u * = y d . ( 32
)
Solution for the constrained problem is found using the Lagrangian:

L(u) = u T Su + λ (η T u -y d ). ( 33
)
Minimizing it over the control variable and solving for the Lagrange multiplier, we find that the explicit solution to the minimization problem is given by

             λ * = - 2y d η T S -1 η , u * = y d η T S -1 η S -1 η, x * = - y d η T S -1 η A -1 BS -1 η. ( 34
)
Without loss of generality we will assume that S is positive definite for the future analysis. The solution (34) cannot be used explicitly due to the fact that the system matrices are assumed to be unknown. But the next section introduces an algorithm which is able to stabilize the system in the arbitrary small neighbourhood of this solution.

C. Extremum seeking for squared deviation minimization

Extremum seeking is a form of adaptive control where the steady-state input-output characteristic is optimized, without requiring any explicit knowledge about this input-output characteristic other than that it exists and that it has an extremum [START_REF] Ariyur | Real time optimization by extremum seeking control[END_REF], [START_REF] Tan | Extremum seeking from 1922 to 2010[END_REF]. This algorithm, developed in the first part of XX century, explores the control space with small oscillations and provides an approximation of the gradient, which then can be integrated in order to find the optimum.

In standard realisations of extremum seeking, one adds to a current control input an oscillating signal, which should be small and slow in comparison with the system dynamics. Further, multiplying the output by the same oscillating signal, it is possible to recover an estimate of the gradient of the output with respect to the input. This standard algorithm is unfortunately not usable for us, since we want to perform a constrained optimization [START_REF] Briat | Robust stability and stabilization of uncertain linear positive systems via integral linear constraints: L1-gain and l∞-gain characterization[END_REF] with a constraint that the average steady state should be equal to the desired one. However, if we modify the algorithm so to minimize the Lagrangian (33) instead of the squared deviation itself, we will optimize the original squared deviation while preserving the average state constraint. This modification leads to the introduction of Lagrange multiplier λ , which can be reconstructed by an additional integrator. Assume the control law for the system (2) is given by

     u = -κωr(ωt)(V + λ (y -y d )), λ = κaωκ λ (y -y d ), u = ū + ar(ωt), (35)
where a and κ are small gains, ω is a small frequency, κ λ is a relative Lagrange multiplier adaptation gain, and the oscillating signal r(ωt) is defined as

r(ωt) = √ 2 • sin(2πωt) cos(2πωt) sin(4πωt) • • • T .
Inputs to this control scheme (Fig. 10) are the output average y and output squared deviation V . Therefore, the control law does not use any state feedback or the system matrices. This is a multi-variable [START_REF] Moase | Non-local stability of a multi-variable extremum-seeking scheme[END_REF] extremum seeking control scheme, augmented with an additional integrator for the adaptation of the Lagrange multiplier as in primal-dual method [START_REF] Nedić | Subgradient methods for saddle-point problems[END_REF], [START_REF] Simpson-Porco | Input-output performance of linear-quadratic saddle-point algorithms with application to distributed resource allocation problems[END_REF] with κ λ being the relative speed of adaptation.

To begin with the proof of the stability of the control scheme [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF], we first present a preliminary analysis, treating ω, κ and a as "small" parameters. Assuming ω is small, we define a new time-scale τ := ωt, which should be slow enough such that the dynamics of the system (2) is much faster than the dynamics of the adaptation. Under this time-scale the closedloop system equations are

             ω dx dτ = Ax + B( ū + ar(τ)), d ū dτ = -κr(τ) x T Px + λ (Cx -y d ) , dλ dτ = κaκ λ (Cx -y d ). (36)
With small ω the singular perturbation analysis [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF] can be performed, thus system dynamics is substituted with its steadystate input-state mapping, i.e. x * = x * (u) = -A -1 Bu. The reduced dynamics is then approximated by

           d ū dτ = -κr(τ) ( ū + ar(τ)) T S( ū + ar(τ))+ + λ η T ( ū + ar(τ)) -y d , dλ dτ = κaκ λ (η T ( ū + ar(τ)) -y d ). (37)
As a next step we introduce an additional time-scale θ := κτ, using which the system becomes

           d ū dθ = -r(τ) ( ū + ar(τ)) T S( ū + ar(τ))+ + λ η T (( ū + ar(τ)) -y d ) , dλ dθ = aκ λ (η T ( ū + ar(τ)) -y d ). (38)
This system is periodic in τ with unit period. When κ is small, the reduced dynamics can be approximated well by the dynamics averaged over the unit period:

                     d ūav dθ = - 1 0 r(σ ) ( ūav + ar(σ )) T S( ūav + ar(σ ))+ +λ av η T ( ūav + ar(σ )) -y d dσ , dλ av dθ = 1 0 aκ λ (η T ( ūav + ar(σ )) -y d ) dσ . (39)
Recall that by the definition of r(•) the oscillating signal has the following properties: 1 0 r(σ )dσ = 0 and 1 0 r(σ )r(σ ) T dσ = I. Then we can rewrite the system:

     d ūav dθ = -2aS ūav -aηλ av -a 2 R, dλ av dθ = aκ λ η T ūav -aκ λ y d , (40)
where R = 1 0 r(σ )r(σ ) T Sr(σ )dσ . It can be seen that 2Su + ηλ is the gradient of the Lagrangian, thus

     d ūav dτ = -a∇ ūav L + O(a 2 ), dλ av dτ = aκ λ ∇ λ av L, (41)
which converges to O(a) of the explicit solution (u * , λ * ).

Concretely, analysing steady state we obtain

y * av = η T ū * av ≡ y d , λ * av = λ * -a η T S -1 R η T S -1 η , ū * av = u * + a 2 η T S -1 R η T S -1 η S -1 η -S -1 R . (42)
The rigorous stability proof is based on the notion of semiglobal practical asymptotic stability: Definition 2 (SPA stability, [START_REF] Tan | On non-local stability properties of extremum seeking control[END_REF]). Consider the parametrized family of systems:

ẋ = f (t, x, ε 1 , ε 2 , ..., ε l ), (43)
where x ∈ R n and parameters of the system ε i > 0 ∀i = 1, 2, ..., l. The system ( 43) is said to be semi-globally practically asymptotically (SPA) stable in [ε 1 , ε 2 , ..., ε l ] at x * , if there exists β ∈ K L [START_REF] Khalil | Nonlinear systems[END_REF] such that the following holds: for each pair of strictly positive numbers (∆, ν), there exists ε * 1 > 0 and for any

ε 1 ∈ (0, ε * 1 ) there exists ε * 2 = ε * 2 (ε 1 ) > 0 and for any ε 2 ∈ (0, ε * 2 ) there exists ε * 3 = ε * 3 (ε 1 , ε 2 ) > 0, ..., there exists ε * l = ε * l (ε 1 , ε 2 , ..., ε l-1
) > 0 such that for any ε l ∈ (0, ε * l ) the solutions of [START_REF] Simpson-Porco | Input-output performance of linear-quadratic saddle-point algorithms with application to distributed resource allocation problems[END_REF] with the parameters [ε 1 , ε 2 , ..., ε l ] satisfy:

|x -x * | β (|x 0 -x * |, (ε 1 • ε 2 • • • ε l )(t -t 0 )) + ν (44)
for all t t 0 0, x(t 0 ) = x 0 with |x 0 -x * | ∆.

Remark 1. Note that the order of the parameters [ε 1 , ε 2 , ..., ε l ] is very important, because the bound for every parameter depends on the choice of all previous parameters, i.e. ε * 3 depends on the chosen ε 1 and ε 2 .

Theorem 5. System (2) with applied control law (35) is SPA stable in [a, κ, ω] at (x * , u * , λ * ).
Proof. First we see that the system ( 40) is a linear system with the system matrix M = -S -η κ λ η T 0 , multiplied by a.

The matrix M is stable, which can be shown by analysing its eigenvalues. Assuming µ 0 being eigenvalue of M, we show that the characteristic polynomial can have no roots. By the Schur complement:

det(M -µI) = -κ λ η T (S + µI) -1 η -µ det(-S -µI). ( 45
) Matrix S is positive definite, thus S + µI is also positive definite. Therefore det(-S -µI) = 0, η T (S + µI) -1 η > 0 and µ 0, which means that (45) cannot equal zero.

We see that the matrix M have no non-negative eigenvalues, which means that the system [START_REF] Fan | Topological proofs for certain theorems on matrices with nonnegative elements[END_REF] is SPA stable in a.

Then, using Lemma 1 from [START_REF] Tan | On non-local stability properties of extremum seeking control[END_REF], we see that the system (37) is SPA stable in [a, κ], and finally using Lemma 2 from [START_REF] Tan | On non-local stability properties of extremum seeking control[END_REF], we conclude that the original closed-loop system is SPA stable in [a, κ, ω], which concludes the proof.

Note that the stability of the closed-loop system (2)- [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF] heavily depends on the chosen parameters. We proved SPA stability in [a, κ, ω], which by definition means that the bound for κ depends on the chosen α, and the bound for ω depends on the chosen α and κ. Therefore, it is difficult to find any rigorous bounds for how small these parameters should be. We can make only heuristic assumptions, such as requiring that all these parameters should be an order of magnitude smaller than the system impulse response. Remark 2 (Additional integral controller). The extremum seeking scheme [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF] can be enhanced by an array of possible modifications. For instance, an additional integral controller can be added:

     u = -κaωγ(y d -y) -κωr(ωt)(V + λ (y -y d )), λ = κaωκ λ (y -y d ), u = ū + ar(ωt), (46)
where the vector γ such that η T γ > 0 plays the same role as in the first sections of this paper regarding average control. I.e. for positive system it is enough to take γ > 0, as in Theorem 2.

Stability proof for the controller [START_REF] Khalil | Nonlinear systems[END_REF] follows exactly the same steps as in Theorem 5, replacing S > 0 by S + γη T > 0.

This scheme can provide much faster convergence of the average (see the examples), which means that the bigger adaptation gains can be used without the possibility for λ to diverge. The equilibrium point for the averaged reduced model of the closed-loop system with this scheme is [START_REF] Nedić | Subgradient methods for saddle-point problems[END_REF], exactly the same as in the previous case.

There are a lot of other possible modifications of the extremum seeking scheme that can be usefully included, for example the high-and low-pass filters that are added to pick up the adaptation signal [START_REF] Tan | Extremum seeking from 1922 to 2010[END_REF].

D. Examples for squared deviation minimization

All the algorithms presented before were tested on a graph constructed as a random Erdős-Rényi graph with n = 40, probability of creating an edge (with weight 1) p = 0.1 and self-loops with weight -5. The dimension of the control vector was chosen k = 3, with matrix B ∈ R n×k being filled randomly: each element was set either to 0 or 1 with equal probability. Matrices C and P were chosen such that y corresponds to the average and V to the squared deviation of the states of the system. Desired value for the average was set y d = 5.

In order to compare the speed of different algorithms we calculated characteristic times for the dynamics of average and squared deviation, defined as a negative inverse of the largest eigenvalue of the closed-loop system. Results of the simulations of different algorithms are presented in Figures 1112131415.

To begin with, we apply the integral controller (11) to the system, and the dynamics of the average state y and the squared deviation V are shown on Now we aim to minimize the deviation of the system states together with controlling the average. If the extremum seeking scheme ( 35) is used, the goal is achieved, and the performance is shown in Fig. 12. The gain values are ω = 0.1, κ = 2, a = 1 and κ λ = 0.01. This scheme is rather difficult to tune and also very slow, the characteristic time is T = 2350, three orders of magnitude higher than in Fig. 11. Also, significant oscillations in y, V and u can be seen even after the convergence of the system. To minimize the oscillations in the extremum seeking, one needs to minimize the gains, but this would lead to an increased convergence time. Therefore we may try to improve the extremum seeking controller by adding the timedependence to the gains, making them large at the beginning and decreasing them over time. In adaptation algorithms the gains should decrease slower than 1/t [START_REF] Borkar | Stochastic approximation: a dynamical systems viewpoint[END_REF], otherwise the algorithm does not converge. In the extremum seeking [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF] the gains a and κ are multiplied together, therefore their product should decrease slower than 1/t. In this example we set a = κ = Performance of the extremum seeking scheme ( 46) is shown on Fig. 14. The gain values are ω = 0.15, κ = 2, a = 1 and κ λ = 0.1. With respect to the scheme [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF], the gain κ λ can be chosen larger. This leads to the faster adaptation of the Lagrange multiplier, thus this scheme works faster than [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF]. The dynamics for average and deviation now behave differently due to the additional controller for the average, therefore it makes sense to find separate characteristic times for them. The characteristic time for average is just T y = 3.9, while for squared deviation it is T V = 732. The parameter of the integral controller is the same as in the case of the integral controller for average, κaωγ = [0, 0, 3] T .

Finally, the implementation of the scheme (46) with timedecreasing gains is presented in Fig. 15 characteristic time for squared deviation being T V = 41 and for average T y = 18. We see that the scheme [START_REF] Khalil | Nonlinear systems[END_REF] in general works faster that the scheme [START_REF] Kottenstette | On relationships among passivity, positive realness, and dissipativity in linear systems[END_REF], but both of them are too slow to compare with the simple average controller [START_REF] Pasqualetti | Controllability metrics, limitations and algorithms for complex networks[END_REF]. Their performance can be significantly increased using time-varying gains, although this leads to a large overshoot at the beginning.

V. CONCLUSION

In this paper, we considered the problem of the control of aggregates of a large-scale network system (namely, its average and variance). First we studied a linear output control problem and examined the general properties of the transfer functions of the system and the controller. We then studied the integral controller for the linear output regulation and formulated sufficient condition CA 2 > 0 for the convergence of any positive integral controller, giving in addition an example showing the conservatism of this condition. If the system satisfies this condition, the parameters of the controller can be chosen arbitrarily, and there is no need to have knowledge of the state vector or of the values of the elements of the A matrix. Control of the average state does not mean that the individual system states will be close to the average state. Therefore, in addition to controlling the average it is worth to minimize the deviation of the system states. To solve this problem we used the extremum seeking algorithm augmented with the primal-dual method for the constrained minimization. The stability of this scheme was proven and its performance, as well as that of several modified versions, was tested in the numerical simulations.

To conclude, we would like to further discuss the scope of application of this work. In the introduction, we argued that average and deviation can be directly measured in several practical examples. In an even broader range of cases, however, average and deviation can be estimated through sampling some nodes and constructing suitable observers, as recently illustrated in [START_REF]Average state estimation in large-scale clustered network systems[END_REF], [START_REF] Niazi | State variance estimation in large-scale network systems[END_REF]: the inclusion of such observers in our control scheme should be a topic of future work.

APPENDIX A PROOF OF THE LEMMA 1

By the definition of inverse M L = (M + ibI)(L + i L) = I, which decomposes into real and imaginary parts:

ML -b L = I, M L + bL = 0. ( 47
)
From the second equation L = -bM -1 L, and substitution of L into the first equation gives

ML + b 2 M -1 L = I, (48)
which means L = M + b 2 M -1 -1 .

APPENDIX B PROOF OF THE LEMMA 2

Denote L(t) = (M + tM -1 ) -1 . We need to prove that CL(t) > 0 for all t 0. The idea of the proof is to provide series expansion for CL(t) such that each term in the expansion is positive. First the coefficients of Taylor series will be computed and then summation by parts will be used twice to obtain series with positive terms. The proof of the Lemma is separated into subsections A-F.

A. Series expansion

Matrix M is an M-matrix, which by definition means that there exists some matrix P with P i, j 0, ρ(P) < 1 and scalar s > 0 such that M = s(I -P). Now make the following transformations: L(t) = (M + tM -1 ) -1 = M(M 2 + tI) -1 = s(I -P) s 2 (I -P) 2 + tI .

(

) 49
Multiplier s s 2 +t is always positive, thus it doesn't affect the sign of the result, so in future we will omit it. Now denote α = s 2 s 2 +t . By definition of t and s this variable satisfies 0 < α 1. Case α = 1 is trivial (it corresponds to the case t = 0), thus often in the following we will use 0 < α < 1. Then L = (I -P) I -2αP + αP 2 -1 .

(50)

We aim to find a coefficients in formal series expansion of L in the powers of P:

L = +∞ ∑ k=0 L k P k . (51)

B. Coefficients of the series expansion

We can introduce scalar function F(x) which has the same expansion as (51) and for which recursive computation of series coefficients is possible. Concretely, define

f (x) = 1 -x 1 -2αx + αx 2 , (52)
where x ∈ [0, 1). Writing the same expansion as L:

f (x) = (1 -x) +∞ ∑ k=0 (2αx -αx 2 ) k . ( 53
)
At the same time, f (x) can be expanded as Taylor series centered at 0:

f (x) = +∞ ∑ k=0 f (k) (0) k! x k . ( 54
)
Power series expansion is unique, thus coefficients L k = f (k) (0) k! . The next step is to determine derivatives of f (x) evaluated at x = 0. Let us introduce function g(x) = 1 -2αx + αx 2 .

(55)

It is obvious that f (x)g(x) = 1x. Now take n-th derivative of this multiplication: (56) Function g(x) is a polynomial of the degree 2, thus its derivatives can be explicitly written: g (0) (0) = 1, g (1) (0) = -2α, g (2) (0) = 2α, (57) and all higher derivatives are zero. Moreover, (1x) (0) (0) = 1 and (1x) (1) (0) = -1 with all higher derivatives also zero.

Recall that L n = f (n) (0) n! . Using (56) we have the following recurrent relation for L n :

L n -2αL n-1 + αL n-2 = 0, ∀n 2, ( 58
)
with initial conditions L 0 = 1 and L 1 = 2α -1.

C. Solving the linear recurrent equation

Equation ( 58) is a linear recurrent equation, which solution is found by solving the characteristic polynomial

λ 2 -2αλ + α = 0. ( 59
)
For 0 < α < 1 roots are complex conjugate pair (λ , λ * ) with

λ = α + i α(1 -α), |λ | = √ α. ( 60
)
The general solution to the equation ( 58) is given by L n = Re [zλ n ] , where z is a complex value that should be determined from the initial conditions. From L 0 = 1 we simply recover Re z = 1, and from L 1 = 2α -1 it is found that Im z = 1-α α . Thus the solution to the equation ( 58) is given by

L n = Re 1 + i 1 -α α α + i α(1 -α) n . (61)

D. Back to the matrix equation

It is established that matrix L can be expressed by the series

L = +∞ ∑ k=0 L k P k , (62)
where L k are given by (61). Now it is evident that (62) is a convergent series due to ρ(P) < 1 and |λ | = √ α < 1. Coefficients L k can be both positive and negative, thus in general matrix L should not be positive. But we want to prove positivity of the vector CL:

CL = +∞ ∑ k=0 L k CP k > 0.
(63)

E. Properties of {CP k } sequence Now it is time to use the condition CM 2 > 0. First of all, M is an M-matrix, thus for any vector x inequality xM > 0 implies x > 0. Therefore CM > 0 (and actually C > 0 automatically).

From CM > 0 we obtain C(I -P) > 0, which means C > CP. Moreover, matrix P is positive, thus multiplying both sides of this inequality on P preserves it. Thus the order relation holds:

C > CP > CP 2 > CP 3 > ... > 0. ( 64
)
Therefore sequence {CP k } is monotonically decreasing with a limit zero (because ρ(P) < 1). Now let us use the next condition, CM 2 > 0. Essentially it means C(I -P)(I -P) > 0, or (C -CP) > (CP -CP 2 ). Again, multiplication by P preserves order, so we have

C -CP > CP -CP 2 > CP 2 -CP 3 > ... > 0, (65)
or

CP k -2CP k+1 +CP k+2 > 0. ( 66
)
This implies that sequence {CP k -CP k+1 } is also monotonically decreasing to zero. In some sense this is equivalent to the "convexity" of {CP k } sequence.

F. Summations by parts

For any series (69)

N ∑ k=0 x k y k = x N Y N - N-1 ∑ k=0 (x k+1 -x k )Y k , (67)
Let us calculate H n :

H n = n ∑ k=0 L k = Re z n ∑ k=0 λ k = Re z 1 -λ n+1 1 -λ , (70)
where λ = α + i α(1α) and z = 1 + i 1-α α . Multiply nominator and denominator by (1λ * ):

H n = 1 1 -α Re z(1 -λ * )(1 -λ n+1 ) . (71)

Figure 1 .

 1 Figure 1. Scale-Free network with hubs (shaded in yellow) controlled from the boundary nodes (double circles)

x 2 :

 2 inner nodes

Figure 2 .

 2 Figure 2. Network with boundary and inner nodes separation

Figure 3 .

 3 Figure3. Control loop given by closed-loop transfer function[START_REF] Liu | Controllability of complex networks[END_REF]

Figure 4 .

 4 Figure 4. Feedback interconnection of passive systems

Figure 5 .

 5 Figure 5. Network with star topology with n leafs. Boundary nodes in green.

Figure 6 .

 6 Figure 6. Output control of the star network with n = 20 leafs. α = 2, β = 1.1, κ = 12, y d = 5. (a). Output y for different γ vectors. (b). Spread of states x for γ = (1, 0) T corresponding to the control of the central node. All the leaf states x 1 ..x 20 have the same asymptotic value 4.751 (which is obvious from the symmetry), while the central state converges to 9.978.

Figure 7 .

 7 Figure 7. Directed line network with n nodes. Boundary node in green

Figure 8 .

 8 Figure 8. Output control of the directed line network, α = 0.2, κ = 12, y d = 5. Inset (b). n = 4, β = 0.002. Output y of the network is unstable, CA 2 > 0. Insets (a), (c), (d). n = 100, β = 0.2. Network is stable. (a). Output y. (c). Spread of states x. (d). Values of lim t→∞ x i depending on the number i ∈ {1..100}, which is the distance from the controlled node.

9

 9

  .

Figure 9 .

 9 Figure 9. Output control of Erdős-Rényi graph for n = 4000 nodes, κ = 250. (a). Dynamics of output y. (b). Dynamics of states x.

Figure 10 .

 10 Figure 10. Extremum seeking scheme for constrained minimization

Fig. 11 .Figure 11 .

 1111 Figure 11. Integral controller (11) for average control, T = 3.84. (a). Average state y, with black dashed line denoting y d . (b). Squared deviation V , with black dashed line denoting V * .

Figure 12 .

 12 Figure 12. Extremum seeking control (35), T = 2350. (a). Average state y, with black dashed line denoting y d . (b). Squared deviation V , with black dashed line denoting V * . (c). Lagrange multiplier λ , with black dashed line denoting λ * . (d). Control vector u, with dashed lines denoting u * .

12 (Figure 13 .

 1213 Figure 13. Extremum seeking control (35) with the gains decreasing over time, T = 434. (a). Average state y, with black dashed line denoting y d . (b). Squared deviation V , with black dashed line denoting V * . (c) and (d): Shortterm plots for y and V . (e). Lagrange multiplier λ , with black dashed line denoting λ * . (f). Control vector u, with dashed lines denoting u * .

Figure 14 .

 14 Figure 14. Extremum seeking control (46), T V = 732 for squared deviation and T y = 3.9 for average. (a). Average state y, with black dashed line denoting y d . (b). Squared deviation V , with black dashed line denoting V * . (c) and (d): Short-term plots for y and V . (e). Lagrange multiplier λ , with black dashed line denoting λ * . (f). Control vector u, with dashed lines denoting u * .

Figure 15 .

 15 Figure 15. Extremum seeking control (46) with the gains decreasing over time, T V = 41 and T y = 18. (a). Average state y, with black dashed line denoting y d . (b). Squared deviation V , with black dashed line denoting V * . (c) and (d): Short-term plots for y and V . (e). Lagrange multiplier λ , with black dashed line denoting λ * . (f). Control vector u, with dashed lines denoting u * .

- 1 =

 1 s(I -P) (s 2 + t)I -2s 2 P + s 2 P 2 -1

  k) (x)g (k) (x).

  This transformation is called Abel transformation or summation by parts. We will apply this procedure twice to obtain series with each term positive.Denote H n = n ∑ k=0 L k . Then H n is bounded because L k consists of powers of λ with |λ | = √ α < 1.By ρ(P) < 1 follows lim k→+∞ CP k = 0. Thus lim k→+∞ H k CP k = 0 and we can writeCL = +∞ ∑ k=0 L k CP k = -+∞ ∑ k=0 H k (CP k+1 -CP k ) = = +∞ ∑ k=0 H k (CP k -CP k+1 ).

  CP k -2CP k+1 +CP k+2 ).

Note that the problems are formulated in the steady-state, therefore we will not pursue any optimization of the transient process.

ACKNOWLEDGMENT

The Scale-FreeBack project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement N 694209).

Product z(1λ * ) = i 1-α α , which is purely imaginary, so

and calculate G n :

Multiply nominator and denominator by (1λ * ):

Product w(1λ * ) = 1, thus this function reads as

This means that G n > 0 for any n 0. Furthermore, by convexity of the sequence {CP k } for any k 0 : CP k -2CP k+1 +CP k+2 > 0.

Thus every term in (69) is greater than zero, which concludes the proof.

APPENDIX C PROOF OF THE LEMMA 3

As in previous proof, define L = (M + tM -1 ) -1 and M = s(I -P). By Lemma 2 CL > 0. Using the series expansion (69), we can write

where all G k > 0. Condition CM 2 Bγ > 0 reads as

Then

where the first term is strictly greater than zero and all others a greater or equal than zero. Thus CLBγ > 0, which concludes the proof.