Nicolò Gionfra
email: nicolo.gionfra@centralesupelec.fr

Guillaume Sandou
email: guillaume.sandou@centralesupelec.fr

Houria Siguerdidjane
email: houria.siguerdidjane@centralesupelec.fr

Damien Faille
email: damien.faille@edf.fr

A Distributed PID-like Consensus Control for Discrete-time Multi-agent Systems

Keywords: Multi-agent Systems, Consensus Control, Discrete-Time LTI Systems

The problem of discrete-time multi-agent systems governed by general MIMO dynamics is addressed. By employing a PID-like distributed protocol, we aim to solve two relevant consensus problems, namely the leaderless consensus under disturbances and leader-follower under timevarying reference state ones. Sufficient conditions for stability as well as two LMI approaches to tune the controller gains are provided. The latter are either based on a H∞ formulation of the problem or on fast response to a reference exogenous signal. Numerical simulations give some insight of which tuning should be considered according to the problem addressed.

INTRODUCTION

In recent years much research effort has been devoted to the area of multi-agent cooperative control because of its wide range of applications and potential benefits. Cooperation of a coordinated multi-agent network is sought via distributed algorithms as they present some interesting advantages over their centralized counterpart, e.g. avoiding single point of failure, reducing communication and computational burden, etc. The main problem in distributed coordination, known as consensus problem, is the one of achieving an agreement on some variables of interest of each agent via local interactions. These variables evolve according to a prescribed dynamics describing the physics of the problem, while interactions among agents are defined by a given communication graph. Finding a distributed protocol to solve the aforementioned problem has been extensively treated for single and double integrator dynamic agents, e.g. [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF]. However, in a more general framework, general dynamics need to be considered in order to describe the agents behavior.

The consensus problem for this latter case has been discussed for both continuous and discrete-time multi-agent systems. In addition, it can be further divided in two main classes of problems, namely leaderless and leader-follower ones. As far as the former is concerned, the most employed distributed protocol is given by a static state feedback law, also called P-like distributed control. One can cite, for instance, [START_REF] Xi | Consensus problems for high-order linear time-invariant swarm systems[END_REF], [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF], [START_REF] Yang-Zhou | Partial stability approach to consensus problem of linear multi-agent systems[END_REF] for the continuous-time framework, and [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF], [START_REF] You | Network topology and communication data rate for consensusability of discrete-time multi-agent systems[END_REF], [START_REF] Su | Two consensus problems for discrete-time multi-agent systems with switching network topology[END_REF], [START_REF] Ge | State consensus analysis and design for highorder discrete-time linear multiagent systems[END_REF] for the discrete one, where the consensus problem is led back to the one of simultaneously stabilizing multiple LTI systems. References [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF], and [START_REF] Su | Two consensus problems for discrete-time multi-agent systems with switching network topology[END_REF]) also solve a leader-follower problem where the leader has an autonomous time-invariant dynamics. Another interesting problem is the one of finding the optimal P-like protocol gain in order to improve consensus under system uncertainties, as in [START_REF] Li | Distributed robust control of linear multi-agent systems with parameter uncertainties[END_REF], and disturbances as in [START_REF] Oh | Disturbance attenuation in a consensus network of identical linear systems: An approach[END_REF], [START_REF] Li | On h∞ and h 2 performance regions of multi-agent systems[END_REF], for continuous time systems, and [START_REF] Wang | H∞ consensus control for discrete-time multi-agent systems with switching topology[END_REF] for discrete-time ones. The proposed approaches usually make use of some H 2 or H ∞ constraints to be respected, and they are in general more involved than the one of simultaneously stabilizing multiple systems. For instance, [START_REF] Li | On h∞ and h 2 performance regions of multi-agent systems[END_REF] provide necessary and sufficient conditions, for the continuous-time case to solve the consensus problem while guaranteeing some properties on the aforementioned norms. On the other hand, for discrete-time systems only sufficient conditions are provided using results from robust control as in [START_REF] Wang | H∞ consensus control for discrete-time multi-agent systems with switching topology[END_REF]. Dynamic distributed controllers are also proposed for consensus achievement based on local output measurements, e.g. [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF]). In the continuoustime framework, [START_REF] Xi | Output consensus analysis and design for high-order linear swarm systems: partial stability method[END_REF] provide a controller with limited energy, while a general full order one is presented in [START_REF] Liu | Dynamic output feedback control for consensus of multi-agent systems: an h∞ approach[END_REF] to achieve some H ∞ performance. Other possible structures have been explored too. Indeed, given the common P-like controller, one can easily think of a more general PID-like structure. In continuoustime, for instance, [START_REF] Carli | A pi consensus controller for networked clocks synchronization[END_REF] propose a PI-like distributed algorithm for single integrator dynamic agents, and [START_REF] Ou | Distributed h∞ pid feedback for improving consensus performance of arbitrary-delayed multi-agent system[END_REF] provide a PID-like controller for general high-order SISO systems. Similar control design is applied to solve a leader-follower consensus under timevarying reference state, as in [START_REF] Ren | Multi-vehicle consensus with a timevarying reference state[END_REF], and in its sampled-data counterpart [START_REF] Cao | Distributed discrete-time coordinated tracking with a timevarying reference state and limited communication[END_REF], where a PD-like protocol is given. Even though the presented literature review is nowhere near exhaustive, one can remark that poorer attention has been devoted to discrete-time dynamic protocols for general LTI MIMO systems, and this is where we wish to place our contribution.

In this paper we propose a PID-like distributed controller for the aforementioned systems, where the agents can communicate on a connected undirected graph, and we provide two possible ways of tuning the controller parameters, based on the solution of LMIs. To the best of our knowledge this distributed control structure has never been fully treated for the mentioned class of dynamic systems. The approach we propose is used to solve two different problems, namely the leaderless consensus under the presence of disturbances, and the leader-follower consensus under a time-varying reference state. Our main results are based on the work of [START_REF] Wu | Lmibased multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF], which we adapted for distributed coordination purposes. The fundamental feature of the aforesaid work is that MIMO PID parameter tuning can be performed via LMIs, avoiding in this way, the need for solving BMIs. Furthermore, in both the analyzed consensus problems the measurement matrix is kept general, allowing a more general problem formulation for the case in which the agents cannot directly measure the variables on which agreement is sought. Eventually, concerning the leaderless consensus, agreement can be focused on particular variables of interest via a proper selection of the controlled output matrix. As for classic control, the PID controller allows good performance despite being rather simple. Concerning the leaderless consensus problem, for instance, it enhances the disturbance rejection, and achieves results that a simple P-like protocol would not permit if the dynamics of the agents are general. Similar conclusions hold for the leader-follower consensus problem with a time-varying reference state, where a P-like control would undoubtedly reach lower performance.

The reminder of this paper is organized as follows. In Section 2 some preliminaries on graph theory are provided and the two main problems are stated. In Section 3 we provide sufficient conditions to solve a leaderless and a leader-follower consensus problem, and we give an LMI approach to tune the distributed PID controller gains. We carry out simulations to test the effectiveness of the proposed controller in Section 4. The paper ends with conclusions and future perspectives in Section 5.

PRELIMINARIES AND PROBLEM STATEMENT

Graph Theory

An undirected graph G is a pair (V, E), where V = {1, . . . , N } is the set of nodes, and E ⊆ V ×V is the set of unordered pairs of nodes, named edges. Two nodes i, j are said to be adjacent if (i, j) ∈ E.

Under the assumption of undirected graph, the latter implies that (j, i) ∈ E too. An undirected graph is connected if there exists a path between every pair of distinct nodes, otherwise is disconnected. The adjacency matrix A = [a ij ] ∈ R N ×N associated with the undirected graph G, considered in this paper, is defined by a ii = 0, i.e. selfloops are not allowed, and

a ij = 1 if (i, j) ∈ E. The Laplacian matrix L ∈ R N ×N is defined as L ii = j =i a ij and L ij = -a ij , i = j.
Considering an undirected graph we make use of the following Lemma 1. [START_REF] Ren | Consensus seeking in multiagent systems under dynamically changing interaction topologies[END_REF] The Laplacian matrix has the following properties: (i) L is symmetric and all its eigenvalues are either strictly positive or equal to 0, and 1 is the corresponding eigenvector to 0; (ii) 0 is a simple eigenvalue of L if and only if the graph is connected.

We will also make use of another Laplacian matrix, according to the following Lemma 2. [START_REF] Lin | Distributed leadless coordination for networks of second-order agents with time-delay on switching topology[END_REF]

Let L = lij ∈ R N ×N be a Laplacian matrix such that lij = N -1 N if i = j,
U LU = I N -1 0 (N -1)×1 0 1×(N -1) 0 Λ, U LU = L 1 0 (N -1)×1 0 1×(N -1) 0
where L 1 ∈ R (N -1)×(N -1) is symmetric and positive definite if the graph is connected.

In addition we employ the Kronecker product ⊗, for which we have Lemma 3. [START_REF] Graham | Kronecker products and matrix calculus with applications[END_REF]

Suppose that U ∈ R p×p , V ∈ R q×q , X ∈ R p×p , and Y ∈ R q×q . The following hold: (i) (U ⊗ V ) (X ⊗ Y ) = U X ⊗ V Y ; (ii) suppose U , and V invertible, then (U ⊗ V ) -1 = U -1 ⊗ V -1 .

Problems Formulation

Problem 1. We consider N identical agents governed by general discrete-time linear dynamics, according to

     x + i = Ax i + B 2 u i + B 1 ω i , i = 1, • • • , N z i = C 1 x i y i = C 2 x i (1) where A ∈ R n×n , B 2 ∈ R n×l , B 1 ∈ R n×h , C 1 ∈ R r×n , C 2 ∈ R m×n , x i x i (k) ∈ R n and x + i x i (k + 1) ∈ R n are

respectively the agent state at the current step k, and at the next step

k + 1, u i u i (k) ∈ R l is the agent control, ω i ω i (k) ∈ R h its disturbance, z i z i (k) ∈ R r
the variable on which agreement among the agents is sought, and y i y i (k) ∈ R m is the measured output. For the sake of leaderless consensus, a priori we do not require A to be Schur stable. Indeed, as shown by [START_REF] Ge | State consensus analysis and design for highorder discrete-time linear multiagent systems[END_REF], A has a role in determining the consensus function to which the agents converge under proper control. Here it can be thought to be assigned by a previous control design step. The agents can communicate on an undirected connected graph whose Laplacian matrix L has positive minimum nonzero and maximum eigenvalues respectively equal to λ L , and λL . At this point, we can state the problem in a general way as the one of finding a distributed control law for

u i such that z i -z j is minimized for i, j = 1, • • • , N with respect to the disturbance ω [ω 1 , • • • , ω N ] .
In this work though, as previously stated, we focus on local controllers of the form

x + c i = A c x c i + B c s i , i = 1, • • • , N u i = C c x c i + D c s i (2)
where

x c i x c i (k) ∈ R 2l
is the agent controller state, and

A c = I l I l 0 l×l 0 l×l 2l×2l B c = (K i -K d ) K d 2l×m C c = I l 0 l×l l×2l D c = [(K p + K i + K d )] l×m (3)
where K p , K i , K d ∈ R l×m are gain matrices to be tuned, and where

s i s i (k) ∈ R m : s i N j=1 a ij (y i -y j ) (4)
Thus the closed-loop system for agent i has dimension n n + 2l. As shown by [START_REF] Wu | Lmibased multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF], system ( 2) is a state representation of the discrete-time PID MIMO controller, whose ztransform is

u i (z) s i (z) = K p + K i z z -1 + K d z -1 z (5)
The problem can now be restated as the one of finding matrices B c , and D c such that the effect of disturbance ω on the consensus is minimized.

The second problem studied in this paper is the following Problem 2. Consider N + 1 discrete-time linear agents, whose dynamics are described by

     x + 0 = Ax 0 + B 1 u 0 z 0 = C 1 x 0 y 0 = C 2 x 0      x + i = Ax i + B 2 u i , i = 1, • • • , N z i = C 1 x i y i = C 2 x i (6) where A ∈ R n×n , B 1 ∈ R n×h , B 2 ∈ R n×l , C 1 ∈ R r×n C 2 ∈ R m×n , x 0 x 0 (k) ∈ R n is the state of the N + 1 agent, called leader, y 0 y 0 (k) ∈ R m is its measured output, u 0 u 0 (k) ∈ R h is
a timevarying unknown control acting on the leader dynamics, and z 0 z 0 (k) ∈ R r is the variable on which we want the follower controlled outputs z i to converge. Concerning the remaining N follower agents, system description similar to (1) holds. The followers are assumed to communicate on an undirected connected graph whose Laplacian matrix is L. The leader can pass information to a subset of followers. If agent i receives information from the leader, then we set a i0 to 1, and 0 otherwise. Thus we define M L -diag(a 10 , • • • , a N 0 ), which is symmetric and positive definite, and we name λ M , and λM respectively its minimum and maximum eigenvalue. Without loss of generality we consider A to be Schur stable. The aim of the present problem is indeed not the one of stabilizing each single agent, but rather to steer the follower agents state to the leader one despite the presence of u 0 , which makes the leader dynamics time-varying. In order to accomplish such objective we aim to employ the controller of form ( 2), ( 3), where we consider a modified variable s i to take into account the communication with the leader agent, according to

s i = N j=1 a ij (y i -y j ) + a i0 (y i -y 0 ) (7)
Intuitively such a controller is not capable of solving the leader-follower tracking problem, i.e. lim k→∞ z i -z 0 = 0 for i = 1, • • • , N , and for any vector signal u 0 , because the latter acts as an unknown exogenous signal for the overall system including the N + 1 agents. This is why we will focus on tuning the controller matrices B c , and

D c such that z i -z 0 is minimized for i = 1, • • • , N .

MAIN RESULT

H ∞ Output Consensus

In order to state our main result we introduce the following definition, similar to the one given in [START_REF] Wang | Consensus for discretetime multiagent systems[END_REF]. Definition 1. System (1) is said to achieve an H ∞ output consensus with a performance index γ ∈ R + if, for any initial condition, lim k→∞ z iz j = 0 for i, j = 1, • • • , N when ω = 0, and the H ∞ norms of the transfer function matrices, for

i = 1, • • • , N , between ω and z i - 1 N N j=1 z j are inferior to γ.
The following result is based on Theorem 3 in [START_REF] Wu | Lmibased multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF], reported in Theorem 4 in the Appendix.

Theorem 1. Given N agents described by ( 1) on an undirected connected graph; consider the distributed protocol of equations ( 2),( 3),( 4); then the agents achieve H ∞ output consensus with performance index γ if there exist two symmetric positive definite matrices P , P ∈ R n×n such that the LMI conditions of Theorem 4 are simultaneously satisfied for two LTI systems whose matrices are respectively (A, B 2 , λ L C 2 ), and (A, B 2 , λL C 2 ), and they both have controlled output matrix C 1 , and disturbance input matrix B 1 .

Proof. The closed-loop dynamics for the generic agent i, by using (1),(2), and by defining the augmented state

ξ i x i , x c i ∈ R n, and ma- trices C2 [C 2 0 m×2l ], C1 [C 1 0 r×2l ], B B 1 0 h×(2l)
is given by

ξ + i = Âξ i + B N j=1 a ij (ξ i -ξ j ) + Bω i z i = C1 ξ i where  = A B 2 C c 0 A c , B = B 2 D c C2 B c C2 (8)
δ i ξ i - 1 N N j=1 ξ j , thus ζ i = C1 δ i . Note that if ζ i = 0 for i = 1, • • • , N then z i = z j , i.e. out- put consensus is achieved. If now we name ξ ξ 1 , • • • , ξ N , δ δ 1 , • • • , δ N , and ζ ζ 1 , • • • , ζ N , we have that ζ = I N ⊗ C1 δ, and δ = ξ -1 ⊗ 1 N N j=1 ξ j = L ⊗ I n ξ
, where L satisfies the conditions of Lemma 2. Gathering together the equations of the closed-loop agents dynamics, we obtain

ξ + = I N ⊗ Â + L ⊗ B ξ + I N ⊗ B ω ζ = I N ⊗ C1 L ⊗ I n ξ = L ⊗ C1 ξ
We now consider the following change of coordinates δ = L ⊗ I n ξ, which yields

δ + = L ⊗ I n I N ⊗ Â + L ⊗ B ξ+ L ⊗ I n I N ⊗ B ω = = L ⊗ Â + LL ⊗ B   δ + 1 ⊗ 1 N N j=1 ξ j   + L ⊗ B ω = L ⊗ Â + LL ⊗ B δ + L ⊗ B ω
where we used points (i) of Lemma 2, and 3. According to the (ii) point of the former, we employ the orthogonal matrix U ∈ R N ×N to define the change of coordinates: δ

U ⊗ I n δ, ω U ⊗ I h ω, ζ U ⊗ I m ζ
, so that the system equations in the new coordinates are given by

                 δ+ = U ⊗ I n L ⊗ Â + LL ⊗ B (U ⊗ I n) δ + U ⊗ I n L ⊗ B ω = Λ ⊗ Â + ΛU LU ⊗ B δ + Λ ⊗ B ω ζ = U ⊗ I m I N ⊗ C1 (U ⊗ I n) δ = I N ⊗ C1 δ (9)
As shown in Lemma 2, being the last rows of Λ, and U LU zeros, we can split the dynamics (9) in two parts by dividing the system variables as

δ = [ δ 1 , δ 2 ] , ω = [ω 1 , ω 2 ] , and ζ = [ ζ 1 , ζ 2 ]
. The dynamic equation of the second variable is then δ+ 2 = 0, and it does not influences δ1 . It follows that we can study the reduced order system

δ+ 1 = I N -1 ⊗ Â + L 1 ⊗ B δ1 + I N -1 ⊗ B ω1 ζ1 = I N -1 ⊗ C1 δ1
From Lemma 2, it exists an orthogonal matrix

V ∈ R (N -1)×(N -1) : V L 1 V Λ = diag(λ 1 , • • • , λ N -1 ), where 0 < λ L ≤ λ i ≤ λL for i = 1, • • • , N -1.
Thus we can define a further change of coordinates, such that δ1

V ⊗ I n δ1 , ω1 V ⊗ I h ω1 , and ζ1 V ⊗ I m ζ1 . The latter yields

δ+ 1 = I N -1 ⊗ Â + Λ ⊗ B δ1 + I N -1 ⊗ B ω1 ζ1 = I N -1 ⊗ C1 δ1 (10)
It is easy to see that the transfer function matrix of (10) satisfies

T ζ1 ω1 (z) ∞ = T ζ1 ω1 (z) ∞ = T ζ ω (z) ∞ = T ζω (z) ∞ (11)
It follows that we can impose an H ∞ constraint on transfer function matrix T ζω (z) by acting on T ζ1 ω1 (z). We can now separate equation (10) in N -1 subsystems, each of them being governed by 12) can be equivalently seen as the closed-loop form of the two following systems

               δ+ 1 i = (A + B 2 D c (λ i C 2 )) B 2 C c B c (λ i C 2 ) A c δ1 i + B 1 0 ωi ζ1 i = C 1 x1 i (12) where δ1 i [x 1 i x 1,c i ] . System (
     x+ 1 i = Ax 1 i + B 2 ūi + B 1 ωi ȳ1 i (λ i C 2 )x 1 i ζ1 i = C 1 x1 i x+ 1,c i = A c x1,c i + B c ȳ1 i ūi C c x1,c i + D c ȳ1 i (13)
Thus, we can reformulate the problem as the one finding matrices B c , and D c such that for i = 1, • • • , N -1 the closed-loop system of ( 13) is Schur stable when ω i = 0, and to guarantee that T zi ωi (z) ∞ < γ. A sufficient condition to prove the existence of such a solution and a relatively simple way to calculate the controller matrices are obtained by employing Theorem 4. In the latter it is proved that if it exists a symmetric positive definite matrix P i ∈ R n×n such that if a given LMI condition is satisfied, then closed-loop system (12) using controller (2),(3),( 4) is such that

δ 1 i (k + 1)P i δ1 i (k + 1) -δ 1 i (k)P i δ1 i (k) < γ 2 ω i (k)ω i (k) -z i (k)z i (k)
It is important to stress that such LMI condition is affine in the system matrices, variables and matrix P i . We make use of this fact to provide sufficient conditions for which it exists a controller of the considered form such that the mentioned LMI is simultaneously verified for

i = 1, • • • , N -1. Since the generic eigenvalue of L 1 : λ i is such that λ L ≤ λ i ≤ λL , then it always exists α i ∈ R : 0 ≤ α i ≤ 1 so that λ i = α i λ L +(1-α i ) λL .
Notice that the systems to be stabilized, appearing in the first set of equation in (13), can be seen as one single system with an uncertain measurement matrix, whose parameter is λ i . In other words, C 2 i λ i C 2 , and ∃α i :

C 2 i = α i C 2 min + (1 -α i )C 2 max ,
where C 2 min λ L C 2 , and C 2 mix λL C 2 , i.e. it can be written as a convex combination of the extreme matrices C 2 min , and C 2 max . Thus, as in [START_REF] Wang | H∞ consensus control for discrete-time multi-agent systems with switching topology[END_REF], the proof makes use of classic results of robust linear control, and in particular by introducing an affine parameter dependent Lyapunov matrix P (α i ) α i P + (1 -α i ) P , where P , P are Lyapunov matrices solution of simultaneous LMI of Theorem 4 written for respectively C 2 min , and C 2 max . Eventually, it is easy to show that if P , P exist, then the controller solves the problem ∀λ ∈ R : λ L ≤ λ ≤ λL , and in particular for λ = λ i , for i = 1, • • • , N -1. Such a controller is easily found from the solution of the aforementioned LMI condition. Indeed among the LMI variables there are matrices B c , and D c , from which it is easy to calculate the PID gains K p , K i , and K d by employing relations in (3).

Remark 1. Note that the mentioned LMI conditions, if satisfied, guarantees that the consensus error is minimized with respect to the disturbance. However the latter still have a role in determining the consensus function to which the agents converge.

Leader-Follower Consensus under time-varying reference

The result of Subsection 3.1 can be easily adapted for the sake of leader-follower consensus via an H ∞ formulation of the problem. Thus, we give the following Definition 2. System ( 6) is said to achieve an H ∞ output leader-follower consensus with a performance index γ ∈ R + if, for any initial condition, lim k→∞ z i -z 0 = 0 for i = 1, • • • , N when u 0 (k) = 0 ∀k ∈ N, and the H ∞ norms of the transfer function matrices, for i = 1, • • • , N , between u 0 and z i -z 0 are inferior to γ.

Theorem 2. Given the system described by ( 6), where N follower agents can communicate on an undirected connected graph, and one leader can communicate with a non-empty subset of followers; consider the distributed protocol of equations ( 2),( 3),( 7); then the systems achieve H ∞ output leader-follower consensus with performance index γ if there exist two symmetric positive definite matrices P , P ∈ R n×n such that the LMI conditions of Theorem 4 are simultaneously satisfied for two LTI systems whose matrices are respectively (A, B 2 , λ M C 2 ), and (A, B 2 , λM C 2 ), and they both have controlled output matrix C 1 , and disturbance input matrix -B 1 .

Proof. The proof is similar to proof of Subsection 3.1. By defining error e i x i -x 0 , ξ [e i , x c i ] , and ζ i C 1 e i the closed-loop system for the generic follower agent i is given by

ξ + i = Âξ i + B N j=1 a ij (ξ i -ξ j ) + a i0 ξ i + Bu 0 ζ i = C1 ξ i
where Â, B, C1 are defined in (8), and B [-B 1 0 h×2l ] . Defining u 0 1 ⊗ u 0 , we then gather the N agent equations together

ξ + = I N ⊗ Â + M ⊗ B ξ + I N ⊗ B u 0 ζ = I N ⊗ C1 ξ (14)
From the definition of M in Section 2, there exists a orthogonal matrix U :

U MU Λ = diag(λ 1 , • • • , λ N ), where λ i ∈ R : λ i > 0 for i = 1, • • • , N , so that we can define the change of co- ordinates ξ (U ⊗ I n) ξ, u 0 (U ⊗ I l )û 0 , ζ (U ⊗ I r ) ζ.
By applying similar calculation as in the previous subsection, the global system in the new coordinates ξ+ =

I N ⊗ Â + Λ ⊗ B ξ + I N ⊗ B û0 ζ = I N ⊗ C1 ξ (15) As in (11), it results that T ξ û0 (z) ∞ = T ξu 0 (z) ∞ , i.e.
we can minimize the effect of u 0 on the consensus error by acting on system (15). Similar to the passage from equations (10) to (12), splitting (15) in N subsystems yields the following equation for subsystem

i                ξ+ i = (A + B 2 D c (λ i C 2 )) B 2 C c B c (λ i C 2 ) A c ξi + -B 1 0 û0 ζi = C 1 êi (16)
where ξi [ê i x c i ] . Equivalently, it can be described as the connection of the two following systems     

ê+ i = Aê i + B 2 ûi -B 1 û0 ŷi (λ i C 2 )ê i ζi = C 1 êi x+ c i = A c xc i + B c ŷi ûi C c xc i + D c ŷi (17)
The rest of the proof is equivalent to the last part of Subsection 3.1, and it is concluded by invoking Theorem 4, whose LMI conditions have to be simultaneously satisfied for the two systems at the vertices of the polytope having matrices respectively (A, B 2 , λ M C 2 ), and (A, B 2 , λM C 2 ), and same controlled output, and disturbance input matrices C 1 , -B 1 . From the solution of the aforementioned LMIs the controllers gains are easily found as in the proof of Theorem 1. If such a solution exists, then the system is stable.

Having employed a PID structure for the distributed controller suggests that consensus should be reached for any u 0 (k) = ū0 , where ū0 is any constant vector. However this is not automatically guaranteed in the MIMO case by the mentioned LMI conditions, and in this framework it is verified a posteriori. Nonetheless, if such LMI has a solution then, according to the well-known Francis equation, a necessary conditions for the proposed controller to reject constant exogenous signals is that l ≥ r.

In the leader-follower consensus framework a different tuning of the PID controller gains with respect to Theorem 2 could lead to better performance, as shown in section 4. Thus, by proposing the following definition we aim to focus on system fast response rather than imposing some H ∞ constraint. For this last development, we further consider r = m, thus we simply name C C 1 = C 2 . Definition 3. System ( 6) is said to achieve fast leader-follower consensus with performance index τ ∈ R + if for u 0 (k) = 0, and any initial condition, lim k→∞ y i -y 0 = 0 for i = 1, • • • , N , and (1e -1 )% of consensus is achieved in a maximum number of steps equal to τ .

Note that the same kind of definition can be considered for sampled-data systems, by saying that system (6) achieves fast leader-follower consensus with a time constant inferior to τ T s , where T s is the system sampling time. The result we present in the following is based on Theorem 2 in [START_REF] Wu | Lmibased multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF], reported in Theorem 5 in the Appendix.

Theorem 3. Given the system described by ( 6), where N follower agents can communicate on an undirected connected graph, and one leader can communicate with a non-empty subset of followers; consider the distributed protocol of equations (2),( 3),( 7); then the systems achieve fast leader-follower consensus with performance in-

dex τ = - 1 log(R)
, where R ∈ R : 0 ≤ R < 1, if there exist two symmetric positive definite matrices P , P ∈ R n×n such that the LMI conditions of Theorem 5 are simultaneously satisfied for two LTI systems whose matrices are respectively (A, B 2 , λ M C), and (A, B 2 , λM C), and where the real constants (a, b) to be set in Theorem 5 are chosen to be (a, b) = (0, R).

Proof. The proof employs the same change of coordinates as in the previous one, so that we can restate the problem as the one of stabilizing the top system of equation ( 17), for i = 1, • • • , N , with the bottom system in (17), i.e. a PID controller whose matrices are defined in (3). Unlike Theorem 2, as previously mentioned, we invoke Theorem 5, where it is stated that given two real constants (a, b), if there exists a symmetric positive definite matrix P i such that a given LMI condition is satisfied, then system ( 16) is stable with all its eigenvalues λ laying in the complex plane region defined by

F D ( [λ], [λ]) : ( [λ] + a) 2 + [λ] 2 < b 2 .
As for the two previous proofs, we employ classic results of linear robust control to impose that this condition is simultaneously satisfied for two systems at the vertices of the polytope whose matrices are respectively (A, B 2 , λ M C), and (A, B 2 , λM C). If such a solution exists then the eigenvalues of system ( 14) are guaranteed to lie in F D . In this framework we are interested in speeding up the system response to u 0 . For this reason we set a = 0, and b = R, where R : 0 ≤ R < 1. Thus, all system eigenvalues are guaranteed to have a module inferior to R. As a result, the system has the slowest time-constant inferior to -T s log(R)

. In terms of number of iterations it is easy to see that such performance is equal to a maximum value -1 log(R) of iterations. Eventually, from the LMI solution, the PID gains are found as in the two previous proofs.

Remark 2. In this latter problem too, having imposed a PID structure does not directly guarantee achievement of consensus for any constant u 0 (k) in the general MIMO case. According to Francis equations, if the mentioned LMI has a solution, a necessary condition though is given by l ≥ m.

Simulation Examples

First of all we carry out a numerical simulation to test the H ∞ output consensus control. We consider a network of 5 agents as show in Fig. 1.a. Note that ( 18) is not Schur stable because two of its eigenvalues lay on the unit circle. Each agent is perturbed by a disturbance of the form ω i (k) = 0.8ν i (k) + c i , where ν i is an aleatory variable with uniform distribution of probability in [0, 1], and c i is some constant value. The PID gains found via LMIs allow an H ∞ performance index of γ = 0.18. Fig. 2 shows the 5 agents trajectories (colored dashed lines) as well as their average (blue continuous line). Then we compare the two proposed PID gain tuning for a leaderfollower consensus problem. For this example we consider the graph of Fig. 1.b, where agent 0 is the leader. The system dynamics is governed by equation ( 6 18). The controller tuned following Theorem 2 allows an H ∞ performance index of γ = 2, while the one tuned according to Theorem 3 guarantees a performance index of τ = 6.1531. In Fig. 3 we simulate the system step response for a value of u 0 = 3. For ease of comparison, we plot here the only output associated to matrix C 1 . As mentioned in section 3, fast consensus (green dashed lines) outperforms the H ∞ one (red dashed-dotted lines). Indeed, even if the latter respects Theorem 2, its consensus error goes slowly to zero with respect to the former one. Eventually, in Fig. 4 it is shown the system behavior for fast consensus tuning of PID gains when u 0 is a time-varying vector signal, and where we set C = C 1 = C 2 = 1 0 0 0 1 0 . The

Conclusion

We presented a PID-like distributed protocol for general LTI MIMO discrete-time agents communicating on an undirected connected graph. By employing LMIs we showed how the controller gains can be tuned to solve two different, yet similar, problems, namely a leaderless under system disturbances and a leader-follower under timevarying reference state consensus problem. Treating the system disturbances in the H ∞ framework revealed good performance, whereas a gain tuning based on fast response seems to be preferable when dealing with a leader-follower problem.

Our results are based on robust control to deal with the problem of simultaneous stabilization of a given number of systems. The given conditions are sufficient and therefore conservative. In near future work we are interested in studying less restrictive conditions when treating the discretetime consensus problem in the H ∞ framework.

APPENDIX

In the following we report the two cited theorems in [START_REF] Wu | Lmibased multivariable pid controller design and its application to the control of the surface shape of magnetic fluid deformable mirrors[END_REF]. Consider the system of equations

x + = Ax + B 1 ω + B 2 u z = C 1 x, y = C 2 x ( 19
)
where R r×(n-3l) , H 1 ∈ R n×r , H 2 ∈ R n×r , H 3 ∈ R h×r , H 4 ∈ R r×r , Y ∈ R q×m , and

A ∈ R n×n , B 2 ∈ R n×l , B 1 ∈ R n×h , C 1 ∈ R r×n , C 2 ∈ R m×n ,
F 22 ∈ R q×3l , 1 ≤ q ≤ 3l, G 1 [G 11 0] ∈ R n×n , G 11 ∈ R n×(n-3l) , G 2 [G 21 0] ∈ R h×n , G 21 ∈ R h×(n-3l) , G 3 [G 31 0] ∈ R r×n , G 31 ∈
N 1 = 0 (n-q)×n 0 (n-q)×2l Y C 2 0 q×2l
and we further name and if exists K such that F 22 K = Y , then the H ∞ norm of the closed-loop system given by ( 19) and

Ψ 11 P -F T -(F T ) , Ψ 21 N 1 + (F T Ã) -G 1 T + (H 1 C cl ) , Ψ 22 -P + G 1 T Ã + (G 1 T Ã) + H 2 C cl + (H 2 C cl ) , Ψ 31 (F T B) -G 2 T , Ψ 32 G 2 T Ã + H 3 C cl + (G 1 T B) , Ψ 33 -γ 2 I + G 2 T B + (G 2 T B) , Ψ 41 -G 3 T -H 1 , Ψ 42 G 3 T Ã + H 4 C cl -H 2 , Ψ 43 G 3 T B -H 3 ,
x + c = A c x c + B c y u = C c x c + D c y (21)
satisfies T zω ∞ < γ.

Theorem 5. Consider system (19). If there exists a positive definite matrix P ∈ R n×n , and a matrix J = J 11 0 (n-q)×3l J 21 J 22

J 22 ∈ R 3l×3l , and X ∈ R 3l×m , and we further name (23) and if J is nonsingular, then by choosing K = J -1 22 X, the eigenvalues of the following matrix

A cl (A + B 2 D c C 2 ) B 2 C c B c C 2 A c lie in the region F D ( [λ], [λ]) : ( [λ] + a) 2 + [λ] 2 < b 2 .

  Similar to[START_REF] Liu | Dynamic output feedback control for consensus of multi-agent systems: an h∞ approach[END_REF]), and ((Wang and Gao, 2011)), we defineζ i z i -1 N N j=1 z j , and

  C 2 as in (

Figure 1 :

 1 Figure 1: (a): leaderless communication graph (left); (b): leader-follower communication graph (right).

Figure 2 :

 2 Figure 2: H∞ output consensus.

Figure 3 :

 3 Figure 3: Step response comparison for H∞ output and fast leader-follower consensus.

Figure 4 :

 4 Figure 4: Leader-follower consensus under timevarying reference.

  x x(k) ∈ R n and x + x(k + 1) ∈ R n are respectively the system state at the current step k, and at the nextstep k + 1, u u(k) ∈ R l is the control input, ω ω(k) ∈ R h is an exogenous input signal, z z(k) ∈ R r the controlled output, and y y(k) ∈ R m is the measured one. Define the matrices C cl C 1 0 r×(2l) , B B 1 0 h×(2l) , K D c B c ,andà A B 2 C c 0 2l×n A cwhere A c , B c , C c , and D c are defined in (3).Assuming B 2 to be of full column rank without loss of generality, there exists an invertibleT b ∈ R n×n T b B 2 = 0 l×(n-l) I l×l . Finally define T T b 0 n×2l 0 2l×n I 2l×2lThus, we have the following theorems Theorem 4.Consider system (19). If there exists a positive definite matrix P ∈ R n×n , where n n + 2l, matrices F = F 11 0 (n-q)×3l F 21 F 22

  following LMI has a solution bP * Ω + JT Ã + aJT b(JT + (JT ) -P ) > 0

  and Ψ 44 I -H 4 -H 4 , such that the following LMI has a solution

		Ψ 11	*	*	*	
	 	Ψ 21 Ψ 22 Ψ 31 Ψ 32 Ψ 33 *	* *	  < 0	(20)
		Ψ 41 Ψ 42 Ψ 43 Ψ 44

ACKNOWLEDGEMENTS

This study has been carried out in the RISEGrid Institute (www.supelec.fr/342p38091/risegriden.html), joint program between CentraleSupélec and EDF ('Electricité de France') on smarter electric grids.


N Gionfra

G Sandou

H Siguerdidjane
email: houria.siguerdidjane@centralesupelec.fr

P Loevenbruck
email: philippe.loevenbruck@edf.fr

D Faille

A Novel Distributed Particle Swarm Optimization Algorithm for the Optimal Power Flow Problem

The distributed optimal power flow problem is addressed. No assumptions on the problem cost function, and network topology are needed to solve the optimization problem. A distributed particle swarm optimization algorithm is proposed, based on Deb's rule to handle hard constraints. Moreover, the approach enables to treat a class of distributed optimization problems in which the agents share a common optimization variable. Under mild communication assumptions, agents are only required to know local variables, cost function, and constraints to solve a common optimization problem. A simulation example is provided, based on a 5-bus electric grid.

I. INTRODUCTION

Optimal power flow (OPF) is a well-known engineering problem that has a central role in the power dispatching and planning of the electric grid. Since its first formulation [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF], the problem has been widely studied and a great number of centralized algorithms have been proposed to solve it (see for instance [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF] and references therein). Such research effort is motivated by the complexity of the problem itself, which happens to be highly nonconvex, mainly due to the power flow equations. Nonetheless, nowadays OPF is still object of great attention. Distributed generation, free electricity market as well as the increasing penetration of renewable energy sources in the grid require the OPF solutions to adapt to a definitely complexified electric network. In particular, the requirements for scalability and efficiency make centralized solutions no longer tractable and justify the research for distributed ones [START_REF] Liu | Consensus-based distributed optimal power flow algorithm[END_REF].

In this regard, numerous solutions are available in the literature. Generally speaking, the latter usually split the problem in two steps. Firstly, it is reduced to a convex one, either by approximation or exact relaxation. Secondly, the obtained optimization problem is distributed among the agents, i.e. the buses of the grid. This step usually capitalizes on alternating direction method of multipliers because it proved good convergence and scalability properties (e.g. [START_REF] Dall'anese | Distributed optimal power flow for smart microgrids[END_REF], [START_REF] Peng | Distributed optimal power flow algorithm for balanced radial distribution networks[END_REF], [START_REF] Magnússon | A distributed approach for the optimal power-flow problem based on admm and sequential convex approximations[END_REF]). Moreover, typically, the optimization problem is formulated in a distributed way via the introduction of an additional constraint called consistency or consensus constraint, which is responsible for the consistency, or agreement, of the variables which are shared among the agents of the network (e.g. [START_REF] Peng | Distributed optimal power flow algorithm for balanced radial distribution networks[END_REF], [START_REF] Magnússon | A distributed approach for the optimal power-flow problem based on admm and sequential convex approximations[END_REF], [START_REF] Liu | Consensus-based distributed optimal power flow algorithm[END_REF]). Considerable attention has been devoted to exact convex relaxation methods (e.g. [START_REF] Peng | Distributed optimal power flow algorithm for balanced radial distribution networks[END_REF], [START_REF] Dall'anese | Distributed optimal power flow for smart microgrids[END_REF]) because they enable to find the OPF global optimum in polynomial time. Unfortunately they are restricted to special network classes, such as radial ones. Other approaches (e.g. [START_REF] Sun | Fully decentralized ac optimal power flow algorithms[END_REF], [START_REF] Magnússon | A distributed approach for the optimal power-flow problem based on admm and sequential convex approximations[END_REF]), even if removing the assumptions on the network topology, are based on some convex approximations of the original OPF. In [START_REF] Magnússon | A distributed approach for the optimal power-flow problem based on admm and sequential convex approximations[END_REF], for instance, the authors rely on convexification of nonconvex constraints. Thus, neither optimality of the solution nor feasibility of the original problem are guaranteed. In [START_REF] Liu | Consensus-based distributed optimal power flow algorithm[END_REF] a consensus based distributed OPF is proposed by employing available convex optimization tools. In addition the consensus constraint is taken into account via the method of the penalty function. Thus, it inherits the problem of distortion of the original objective function [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF].

In this paper, as in [START_REF] Magnússon | A distributed approach for the optimal power-flow problem based on admm and sequential convex approximations[END_REF], we aim at scalability of the solution and refer to an OPF problem formulation which is as general as possible. To accomplish such goals we capitalize on particle swarm optimization (PSO) algorithm as it showed good performance for the centralized OPF problem (e.g. [START_REF] Abido | Optimal power flow using particle swarm optimization[END_REF], [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF], [START_REF] Syai'in | A distribution power flow using particle swarm optimization[END_REF]). In particular we suggest a new distributed version of the aforementioned algorithm, leading mainly to a two-fold contribution. On the one hand, to the authors' knowledge, PSO was never applied to solve OPF in a distributed manner. On the other hand, the proposed algorithm has a general formulation that makes it applicable to a class of distributed optimization problems for which such solution is not available yet. For instance, related works to distributed PSO (DPSO) can be found in [START_REF] Wakasa | Distributed particle swarm optimization based on primal-dual decomposition architectures[END_REF], and [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF]. Unfortunately these approaches cannot be applied to OPF because they solve a problem in which the agents do not share common variables. Approaches in [START_REF] Hereford | A distributed particle swarm optimization algorithm for swarm robotic applications[END_REF] and, [START_REF] Akat | Decentralized asynchronous particle swarm optimization[END_REF] are based on considering each agent as a PSO particle of a common optimization problem. On the one hand this in turns implies that each agent has knowledge of the common cost function, which is typically not the case in optimization problems such as the OPF one. On the other hand, each agent decides on the whole optimization variable, which could be unnecessary, and computationally demanding for the problem addressed in our work. Eventually, while [START_REF] Wakasa | Distributed particle swarm optimization based on primal-dual decomposition architectures[END_REF] makes use of a perturbed primal-dual method, and [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF] capitalizes on projection functions, we choose to handle constraints by employing Deb's rule [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF], by following the centralized PSO proposed by [START_REF] Aguirre | Copso: Constrained optimization via pso algorithm[END_REF].

The remainder of this paper is organized as follows. In Section II the optimization problem is stated in a general way. The main result is presented in Section III, where the main features of DPSO algorithm are detailed. Section IV is concerned with the self-configuration of a finite-time average consensus, which is needed in the proposed DPSO algorithm. OPF and its formulation allowing a problem solution via DPSO is addressed in Section V. We carry out a simulation example to test the effectiveness of the approach in Section VI. The paper ends with conclusions, and future perspectives in Section VII.

II. PROBLEM STATEMENT

Let us provide a general formulation of the optimization problem. The specific case of OPF is derived in Section V. Consider a group of N agents, each of which disposes of a private control variable x i ∈ R n i , i = 1, • • • , N, private cost function, inequality and equality constraints. By private we mean that the latter are only known by the associated agent i. We provide the following useful definition.

Definition 1: Agent i is physically coupled to agent j if at least one among its private cost function, inequality and equality constraints depend on agent j private control variable. Agent j is thus said to be a physical neighbor of agent i. We then consider a graph G p = (V p , E p ) that keeps track of the physical relations among the agents. In particular V p = {1, . . . , N} is the set of the agents, i.e. the nodes of the graph, and E p ⊆ V p × V p is the set of edges among them, where the edge (i, j) ∈ E p if and only if agent j is a physical neighbor of agent i. Note that according to the given definition of physical neighbor, G p generally defines a digraph. The set of physical neighbors of agent i is defined as N p i j ∈ V p : (i, j) ∈ E p . By defining x x x i j x j ∈ R n j : j ∈ N p i as the set of physical neighbors variables of agent i, we are able to represent the generic agent i cost function, inequality, and equality constraints respectively as f i (x i , x x x i j ), g i (x i , x x x i j ) ≤ 0 0 0, and h i (x i , x x x i j ) = 0 0 0, where f i :

R n i + ∑ j∈N p i n j → R, g i : R n i + ∑ j∈N p i n j → R m i , h i : R n i + ∑ j∈N p i n j → R q i ,
and where 0 0 0 is a vector of proper dimension with all zero entries. Note that no assumptions were made concerning these functions. We are now able to state the optimization problem as min

x x x [x 1 •••x N ] F(x x x) min {x i ,i=1,••• ,N} N ∑ i=1 f i (x i , x x x i j ) subject to g i (x i , x x x i j ) ≤ 0 0 0, i = 1, • • • , N h i (x i , x x x i j ) = 0 0 0, i = 1, • • • , N (1)
Thus, agents have to cooperatively minimize a common cost function F : R ∑ N i n i → R while sharing the common optimization variable x x x. In the next section we provide a distributed way to solve [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF].

III. DISTRIBUTED PSO

Before providing the DPSO algorithm we need to clarify how we intend to handle constraints in [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF], and how we need to set the distributed communication among the agents.

A. Constraints Handling

As known, one of the major difficulties when employing PSO to solve a constrained optimization problem is the lack of an explicit method to direct the optimum search towards the feasible region [START_REF] Aguirre | Copso: Constrained optimization via pso algorithm[END_REF]. Having made no constraints assumption we choose a generic method to handle them. In particular, as previously stated, we make use of Deb's rule, which belongs to the penalty function approaches but it does not require any penalty parameter. Moreover, the aforementioned rule allows avoiding any cost function distortion that may occur when incorporating the constraints in the problem via penalty functions. Deb's rule consists of a tournament selection in which, when comparing two solutions of (1), the following criteria is adopted [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF]: (i) any feasible solution is preferred to any infeasible solution; (ii) among two feasible solutions, the one having better objective function value is preferred; (iii) among two infeasible solutions, the one having smaller constraint violation is preferred. In addition, equality constraints are handled via a transformation into inequality constraints with the introduction of a positive

threshold ε ∈ R + such that g i+N (x i , x x x i j ) |h i (x i , x x x i j )|-ε ε ε ≤ 0 0 0, i = 1, • • • , N,
where ε ε ε is a vector of proper dimension with all its entries equal to ε. PSO is a population based algorithm in which a given number of particles in the search space can evaluate the cost function. In order to implement Deb's rule, though, we consider the following modified cost function, called fitness function, which particles need to evaluate to decide on their search direction [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF].

F(x x x)    F(x x x) if g i (x x x) ≤ 0 0 0 i = 1, • • • , 2N f max + N ∑ i=1 m i ∑ k=1 χ g i,k (x x x) + q i ∑ k=1 χ g i+N,k (x x x) otherwise (2)
being χ : R → R the function

χ(y) = y if y > 0 0 otherwise
and where g i,k (x x x), g i+N,k (x x x) are respectively the k-th component of g i (x x x), and g i+N (x x x). Constraint violation is evaluated via the sum of violated constraints appearing in the second equation of (2). f max is the maximum value of F among the available feasible solutions. If some information about the common cost function is known a priori, for ease of implementation, f max can be set as a sufficiently high positive value.

B. Communication Settings

Assumption 1: There exists a communication graph

G 1 c = (V 1 c , E 1 c ) such that G 1 c = G p .
This basically means that each agent being a physical neighbor of another one can directly communicate to it. We now consider the undirected graph associated to it:

Ḡ 1 c = ( V 1 c , Ē 1 c ), where V 1 c = V 1 c , and Ē 1 c is such that if (i, j) ∈ E 1 c then (i, j), ( j, i) ∈ Ē 1 c . Ḡ 1 c is not necessarily connected. Thus, it can be divided in M connected undirected subgraphs Ḡ 1 c,m , m = 1, • • • , M, such that ∪ M m=1 Ḡ 1 c,m = Ḡ 1 c , and Ḡ 1 c,m ∩ Ḡ 1 c,r = ( / 0, / 0), ∀m, r = 1, • • • , M, m = r. Assumption 2: There exist M undirected connected graphs Ḡ 2 c,m , m = 1, • • • , M such that Ḡ 2 c,m ∩ Ḡ 2 c,r = ( / 0, / 0), ∀m, r = 1, • • • , M, m = r, and Ḡ 2 c ∪ M m=1 Ḡ 2 c,m = ( V 2 c , Ē 2 c ) where V 2 c = V 1 c .
The latter assumption simply states that the communication among the M groups of agents defined via Ḡ 1 c,m , m = 1, • • • , M, are not required to communicate via Ḡ 1 c itself, as long as the considered communication Ḡ 2 c keeps the same M groups of agents connected within each group. A graph example is reported in Fig. 1.

Remark 1: For ease of implementation one can simply consider only one communication graph

G c = (V c , E c ) such that V c = V 1 c , E 1 c ⊆ E c , and there exist M subgraphs G c,m , m = 1, • • • , M such that ∪ M m=1 G c,m = G c , G c,m ∩ G c,r = ( / 0, / 0), ∀m, r = 1, • • • , M, m = r.
Remark 2: The considered graphs of communication keep two different groups of agents disconnected whenever there is no physical coupling among them. If the disconnected groups are M, it basically implies that (1) can be split in M separate optimization subproblems, whose solutions can be found separately with no inter-group communication.

C. Variables Settings

Each agent has N p particles associated to its private variable, x i,p ∈ R n i , p = 1, • • • , N p . Moreover, to each of them we associate a variable that keeps memory of the personal best position visited in the search space, named b i,p ∈ R n i , p = 1, • • • , N p . Each x i,p (resp. b i,p ) has a spy copy in each agent for which it is a physical neighbor,

x j i,p (resp. b j i,p ) ∈ R n i ∀ j : i ∈ N p j , p = 1, • • • , N p .
The role of the spy variables is to allow evaluation of the private cost functions, and constraints. In addition, each agent particle (resp. personal best) has its personal estimate of the average common cost function, and average sum of violated constraints evaluated in the optimization variable

x x x p = [x 1,p • • • x N,p ] (resp. b b b p [b 1,p • • • b N,p ]
), i.e composed by all agents particles (resp. personal bests) having the same p index. Respectively we name the former h 1 i,p

1 /N ∑ N l=1 f l (x l,p , x x x l j,p ) resp. h b,1 i,p 1 /N ∑ N l=1 f l (b l,p , b b b l j,p
) , the latter . We keep index i of h 1 i,p , h 2 i,p , (resp. h b,1 i,p , h b,2 i,p ) to stress that agent i estimate of the average value can have a small difference with respect to the other agents one. Furthermore, each particle is endowed with a speed value v i,p , and it has access to the global best position, g i,p , among a subset of particles S i p , i.e. the subset of particles from which x i,p can retrieve information about their personal best. We additionally require such subset to be the same for each particle among the agents having same p index. Thus we can drop its index i: S i p = S p , i = 1, • • • , N. Note that S p can be simply chosen as the set of all agent particles, so that g i = g i,p would be agent i best visited position. Nonetheless it might be convenient to choose S p as a smaller set of particles to prevent the algorithm from premature convergence.

h 2 i,p 1 N N ∑ l=1 1 m l m l ∑ k=1 χ g l,k (x l,p , x x x l j,p ) + 1 q l q l ∑ k=1 χ g l+N,k (x l,p , x x x l j,p resp. h b,2 i,p 1 N N ∑ l=1 1 m l m l ∑ k=1 χ g l,k (b l,p , b b b l j,p ) + 1 q l q l ∑ k=1 χ g l+N,k (b l,p , b b b l j,p

D. PSO Update Law

PSO equations determine the particle motion in the search space in order to find the global optimum. The latter is a wellknown system of stochastic difference equations. For ease of notation, at iteration k of the algorithm we note x i,p x i,p (k), and x + i,p x i,p (k + 1). Same notations hold for v i,p , b i,p , and g i,p . PSO equations, for each particle, are given by

v + i,p = ωv i,p + φ 1 i,p (g i,p -x i,p ) + φ 2 i,p (b i,p -x i,p ) x + i,p = x i,p + v + i,p (3)
(b + i,p , h b i,p ) = (x + i,p , h + i,p ) if h + i,p < h b i,p (b i,p , h b i,p ) otherwise (4)
g + i,p = arg min b + i,p :p∈S p {h b i,p } (5)
where

h i,p h i,p (k) = h 1 i,p if h 2 i,p = 0 f max + h 2
i,p otherwise and h + i,p = h i,p (k + 1). A similar definition holds for h b i,p , h b,1 i,p , h b,2 i,p . Equation ( 4) is the decision step of the algorithm in which the new particles are compared to the according personal best in terms of average fitness function evaluation.

φ 1 i,p φ 1 i,p (k) ∼ U(0, c 1 ), φ 2 i,p φ 2 i,p (k) ∼ U(0, c 2
), i.e. they are two aleatory variables with uniform distribution of probability in the respective intervals [0, c 1 ], and [0, c 2 ], where c 1 , c 2 ∈ R + . Exhaustive works in the literature show how to set parameters c 1 , c 2 , and ω, called inertial factor, in order to ensure convergence of the algorithm, (e.g. [START_REF] Sandou | Using particle swarm optimization for reduced order h synthesis[END_REF]). Here, convergence is meant as ∀i

= 1, • • • , N: lim k→∞ x i,p (k) - x i,t (k) = 0, ∀p,t = 1, • • • , N p ,
as convergence to the global optimum is not guaranteed.

Remark 3: In this paper we additionally consider lower and upper bounds of the variables of the form:

x i ≤ x i ≤ xi , x i , xi ∈ R n i , i = 1, • • • , N,
by limiting the maximal speed of each particle, and forcing x i,p to stay within the aforementioned bounds. This is simply obtained by adding the following equations to (3)

v + i,p = max{min{v + i,p , vi }, -vi } (6) x + i,p = max{min{x + i,p , xi }, x i } (7)
where we place respectively (6) right after the first equation in [START_REF] Liu | Consensus-based distributed optimal power flow algorithm[END_REF], and ( 7) after the second equation in [START_REF] Liu | Consensus-based distributed optimal power flow algorithm[END_REF], and where vi 1 /2( xix i ).

E. Diversity and Dynamic Tolerance

For the case of OPF we consider two additional features of PSO algorithm that enhance the optimality of the solution. These are descripted in [START_REF] Aguirre | Copso: Constrained optimization via pso algorithm[END_REF]. The reader may refer to it for further details. The first one consists in keeping diversity of the flock of particles, which prevents from premature convergence. This is obtained by applying two perturbation operators to the particle personal bests b i,p , namely the C-Perturbation, and the M-Perturbation, which will not be described in this paper. However, in order to reduce the communication burden we consider a unique perturbation, obtained by sequential composition of the aforementioned ones, and we name it CM-Perturbation. Eventually, we employ an additional constraint handling technique, named dynamic tolerance. It consists in letting the threshold ε, defined in Subsection III-A, be iteration-varying, typically linearly decreasing from an initial value ε to a final target value ε.

F. Distributed Algorithm

Before providing the overall algorithm we need to introduce a distributed subroutine that allows the evaluation of the average common cost function and sum of constraints, and that has to be run each time that a decision step, as in (4), has to be performed. Whenever the latter is required, we refer to Algorithm 1 to evaluate the problem functions in the generic point x

x x p [ x1,p • • • xN,p ] ,
introduced to be consistent with the notation. Via the aforementioned subroutine, each of its component xi,p is associated with the estimation of the average fitness function ĥi,p . It is important to notice that its value is dependent on the chosen value of equality constraints threshold ε. In general, when the chosen point to be evaluated is one of the particles we refer to the notation h i,p to indicate the corresponding function estimation. In all the other cases we introduce a superscript to the associated function evaluation with the same name of the considered point, for instance b i,p ↔ h b i,p . The overall distributed algorithm to be run by each agent on each of its particles x i,p is shown in Algorithm 2, where iterations are performed until a prescribed number max iter, and where we introduced r, and s, respectively an aleatory variable r ∼ U(0, 1), and a linearly decreasing value from 1 to 0. In Algorithm 2, lines 10 -16 are used to keep diversity of the flock. In particular t i,p is a temporary value obtained by perturbation of the new computed personal best b + i,p . Lines 17 -21 are only necessary if ε is iterationvarying. Indeed if ε decreases at each algorithm step, then the personal bests need to be reevaluated to see how they fit the new, ε-dependent, fitness function. If not performed, a b i,p computed at the beginning of the iterations has higher probability to have better fitness function value with respect to one computed afterwards. Because of the structure of the proposed algorithm we are able to conclude the following Proposition 1: Consider the centralized PSO algorithm obtained by substitution of Algorithm 1 with direct computation of ĥi,p . If it converges, so Algorithm 2 does. In addition, they converge with same performance.

Remark 4: It is important to note that thanks to Assumption 1 we are able to update the spy variables with one communication step (line 1 in Algorithm 1). Moreover we are able to avoid adding any consensus or consistency constraint to [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF].

Remark 5: Proposition 1 is essentially verified thanks to the average consensus step in Algorithm 1, and to assumption 1. The latter can be relaxed by requiring that each agent that is a physical neighbor to another one has at least a communication path to it. Each of the agents belonging to this path is required to have spy variables of the mentioned physical neighbor. Then, for instance, line 1 of Algorithm 1 can be modified by considering a leader-follower consensus algorithm to allow the update of the spy variables.

IV. FINITE-TIME AVERAGE CONSENSUS

The main communication burden is due to the need for performing the average consensus in Algorithm 1 to evaluate the problem functions. Unfortunately this step cannot be approximated as a one-step communication as shown in [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF], for the accuracy of consensus would not be sufficient to perform decision steps on particles, as (4). Nonetheless we employ a finite-time average consensus algorithm to reduce the communication burden as much as possible. In addition, being interested in a whole-distributed solution, we capitalize on a distributed self-configuration of the aforementioned algorithm. Thus, in this section we aim at providing the configuration approach to tune the average consensus appearing in step 5 of Algorithm 1. The following has to be performed once, and offline.

A. Self-configuration

We consider the self-configuration problem for the generic undirected graph Ḡ 2 c,m = ( V 2 c,m , Ē 2 c,m ), associated to the msubproblem of (1), defined in Subsection III-B. We rename it G = (V , E ) for ease of notation, where V = 1, • • • , N m . Let us state the following graph features.

• The distance between two nodes (i, j) ∈ V : dist(i, j) is the length of the shortest path between i, and j. • The eccentricity of a node i ∈ V : ecc(i) max j∈V dist(i, j). • The radius of G : r(G ) min i∈V ecc(i).

• The diameter of G : d(G ) max i∈V ecc(i).

• The neighborhood of agent i:

N i = { j ∈ V : (i, j) ∈ E }. The state at time t associated to each agent in V is h i (t) ∈ R, i = 1, • • • , N m . The overall network state is h h h [h 1 • • • h N m ] ∈ R N m .

Each agent updates its own state according to

h i (t) = w t ii h i (t -1) + ∑ j∈N i w t i j h j (t -1) (8)
where w t ii , w t i j : j

∈ N i ∈ R, i = 1, • • • , N m , ∀t > 0,
are the weights at time t to be chosen. Equation [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF] in matrix form can be written as h h h(t) = W t h h h(t -1). It is well-known that finite-time average consensus in D steps can be solved if we are able to find D matrices W t , t = 1, . . . , D such that the following is verified 1

∏ t=D W t = 1 N m 1 1 11 1 1 (9)
where 1 1 1 ∈ R N m has all its entries equal to 1. In general D is known to be a value in [d(G ), 2r(G )]. In this paper we refer to the work of [START_REF] Dung | Distributed design of finite-time average consensus protocols[END_REF], where the authors suggest a gradient back-propagation method to solve the factorization problem (9) in a distributed way. Please refer to the mentioned reference for details. This is done by solving the following optimization problem min

{W 1 ,••• ,W D } N m ∑ i=1 Q ∑ q=1 (h i,q (D) -y q ) 2 (10)
where h i,q (D) is obtained via iteration of (8), and h i,q (0), y q , i = 1, • • • , N m , q = 1, • • • , Q are the chosen learning sequences. Nonetheless, being the approach gradientbased, and the optimization problem [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF] nonconvex, the solution is strongly dependent on the algorithm initialization of matrices W t . We notice that (10) can be written as [START_REF] Carpentier | Contribution to the economic dispatch problem[END_REF], with no constraints, by taking the agent i private cost function as ∑ Q q=1 (h i,q (D)y q ) 2 , and common cost function as the objective function of [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF]. Thus, we propose to use Algorithm 2 to provide an initial point to the algorithm shown in [START_REF] Dung | Distributed design of finite-time average consensus protocols[END_REF], as it could help finding a better solution. In this self-configuration step, clearly, average consensus in Algorithm 1 cannot be performed via the finite-time algorithm, which is the aim of this stage. For this reason we consider the following update law for ĥi,p variables at step 5 in Algorithm 1.

ĥ+ i,p = ϕ ii ĥi,p + ∑ j∈N i ϕ i j ĥ j,p (11)
ϕ i j =        1 max{d i , d j } + 1 if (i, j) ∈ E 1 -∑ j∈N i 1 max{d i , d j } + 1 if i = j (12)
where [START_REF] Wakasa | Distributed particle swarm optimization based on primal-dual decomposition architectures[END_REF] are Metropolis-Hastings weights, and d i |N i | is the degree of node i in G . Equation [START_REF] Syai'in | A distribution power flow using particle swarm optimization[END_REF] only requires each agent to know little information about its neighbors, and it has to be run for a sufficient number of steps in order to ensure proper convergence of Algorithm 2. In addition, for the particular problem of [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF], there is no need to use spy variables neither. Agents do share common variables, but this link is kept implicit in the computation of h i,q (D), appearing in [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF], via [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF].

B. Example

We consider MATLAB ® bucky graph in Fig. 2, to show the effectiveness of the approach. This graph has radius, and diameter r = d = 9, and a total number of agents equal to 60. Good results were found for a value of D = 11, i.e. 11 steps are needed to find the average consensus on the network with an acceptable small error. Fig. 3 shows the results of the combined DPSO and gradient back-propagation algorithm versus the solution provided by the only gradient backpropagation, where the dashed-dotted black line represents the average of the initial network state 1 /N m ∑ N m i=1 h i (0). From this example it is clear that DPSO not only can enhance the optimality of the solution of [START_REF] Wannakarn | Optimal power flow problem solved by using distributed sobol particle swarm optimization[END_REF], but it also helps finding one when gradient back-propagation fails. Indeed, in this example, weights found via gradient back-propagation do not let convergence to the average consensus value (black dashed-dotted line) as it shows a steady state constant error. On the other hand, this is attained by employing the algorithm initialization provided by DPSO.

V. OPF MODEL AND PROBLEM FORMULATION

OPF is concerned with minimizing a given cost function while satisfying the electric grid constraints. These are mainly given by the limits on the voltage value, typically allowed in the interval [0.95, 1.05]p.u., and by the power flow equations

P i = U i ∑ j∈N p i U j (G i j cos(θ i -θ j ) + B i j sin(θ i -θ j )) Q i = U i ∑ j∈N p i U j (G i j sin(θ i -θ j ) -B i j cos(θ i -θ j )) (13)
which show how active power P, and reactive power Q at every bus i of the considered electric grid are function of its own voltage values of amplitude U i , and phase θ i , and of its physical neighbors voltages (U j , θ j ) : j ∈ N p i . G i j , B i j are parameters depending on the impedance of the line between buses i, and j. If we consider (U i , θ i ) to be the private variable of bus i, then it is easy to see that ( 13) are already written in the form of the constraint functions of (1). Among the possible cost functions usually employed for OPF, in this paper we choose to minimize the sum of power losses in the transmission lines min

{(U i ,θ i ),i=1,••• ,N} 1 2 N ∑ i=1 ∑ j∈N p i P i j loss + Q i j loss , where (14)
P i j loss = g i j U 2 i +U 2 j -2U i U j cos(θ i -θ j ) Q i j loss = -b sh i j U 2 i +U 2 j -b i j U 2 i +U 2 j -2U i U j cos(θ i -θ j ) ( 15
)
where N is the number of buses, g i j , b i j , b sh i j are other parameters depending on the impedance of the line between i, and j, and

g ii = b ii = b sh ii = 0, i = • • • , N .
The reader may refer to [START_REF] Andersson | Modelling and analysis of electric power systems[END_REF] for a detailed explanation of ( 13), [START_REF] Akat | Decentralized asynchronous particle swarm optimization[END_REF]. In the sequel we neglect the shunt susceptance b sh i j , for it is usually verified that |b sh i j | |b i j |. Note that considering a weighted sum of the kind of P i j loss + αQ i j loss , α ∈ R + , would not change the solution of the optimization problem, as P i j loss , and Q i j loss are basically the same function (being g i j , -b i j > 0). For the same reason one could only consider P i j loss , or Q i j loss and obtain the same solution. We can restate ( 14), [START_REF] Akat | Decentralized asynchronous particle swarm optimization[END_REF], allowing the formulation of the problem in the form of (1), by considering the following min

{(U i ,θ i ),i=1,••• ,N} N ∑ i=1 P i loss + Q i loss , where (16)
P i loss ∑ j∈N p i g i j U 2 i -g i j U i U j cos(θ i -θ j ) Q i loss ∑ j∈N p i -b i j U 2 i + b i j U i U j cos(θ i -θ j ) (17)
Clearly ( 16), [START_REF] Sandou | Using particle swarm optimization for reduced order h synthesis[END_REF] are equivalent to [START_REF] Hereford | A distributed particle swarm optimization algorithm for swarm robotic applications[END_REF], [START_REF] Akat | Decentralized asynchronous particle swarm optimization[END_REF]. We are now able to formulate the chosen OPF problem. Consider an electric grid of N buses, N G of which are generators, here considered PU buses, N L are loads, here considered PQ buses, and one is a Uθ bus, used to balance generation, load and losses. Each bus has the private cost function: P i loss + Q i loss . Thus, the common cost function is the objective function of [START_REF] Aguirre | Copso: Constrained optimization via pso algorithm[END_REF]. Moreover they all have their voltage amplitude constrained in the previously mentioned allowed interval. This kind of constraint is treated via ( 6), [START_REF] Sun | Fully decentralized ac optimal power flow algorithms[END_REF]. The Uθ bus also gives its θ as a reference to the grid, and it is simply set as 0 • . Every generator bus i presents an equality constraint on P i , set as Pi , which is a problem data, via the first equation of ( 13), and an inequality constraint on Q

i , Q i min ≤ Q i ≤ Q i max ,
employing the second equation of [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF]. Eventually, every load bus j presents two equality constraints given by Pj , and Q j , data of the problem, via [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF]. Overall, each bus of the network presents two generalized nonconvex private inequality constraints. Thus the optimization problem is [START_REF] Aguirre | Copso: Constrained optimization via pso algorithm[END_REF], [START_REF] Sandou | Using particle swarm optimization for reduced order h synthesis[END_REF] subject to for i = 1, • • • , N 0.95 ≤ U i ≤ 1.05 via [START_REF] Wakasa | Distributed particle swarm optimization using an average consensus algorithm[END_REF] :

P i = Pi , Q i = Qi if i is a PQ bus P i = Pi , Q i min ≤ Q i ≤ Q i max if i is a PU bus VI. SIMULATIONS
We illustrate the performance of the proposed algorithm on a 5-bus grid shown in Fig. 4, which is a modified scheme of a grid example in [START_REF] Andersson | Modelling and analysis of electric power systems[END_REF]. The chosen communication graph (red dashed-dotted lines) simply traces the electric lines among the buses. The considered grid has nominal voltage equal to 63 kV, and the grid cables impedance are set to be R = 0.15 Ω/km, and X = 0.21 Ω/km. Bus 5 is the Uθ bus, and it can be considered as a source substation. The generators are prescribed to produce P1 = 50 MW, P3 = 70 MW, and load buses require P2 = 60 MW, P4 = 85 MW. As far as the reactive power is concerned, load buses are operated at constant power factor, cos ϕ, equal to 0.97. Generators have to keep the reactive power within given limits with respect to their active power production. Indeed the reactive power supplied above a Q /P ratio of (tan ϕ) max or absorbed under (tan ϕ) min is subject to the application of a charge. In France, for instance, it amounts to 16.3 e/MVARh for electric grids whose nominal voltage is in [50, 130] kV. In particular (tan ϕ) max , (tan ϕ) min are set respectively as 0.4, and -0.35. According to the TURPE 41 , the cost for the active power losses is set to 53 e/MWh. As far as the DPSO set-up is concerned, we chose N p = 70 particles for each bus, and max iter = 600. Having normalized the power and voltage values, we considered a threshold linearly decreasing from ε = 1 to ε = 10 -5 until the 90% of the iterations. For the considered graph, finitetime average consensus algorithm to be run in Algorithm 1, only requires 2 steps to be solved. Fig. 5 shows the particles behavior during the iterations. A particle is assigned with a sum of constraints value whenever it is not feasible, and with a common cost function value otherwise. At the beginning, being randomly initialized, they are very likely to be infeasible. Nonetheless they reach feasibility quite quickly since ε is still enough large. As iterations grow, ε becomes smaller, thus, the common cost function respecting constraints is updated, and this explains its increasing pace from about iteration 100. At the very end all particles move according to their sum of constraints value since ε reaches a very small value. However, this fact does not have to be interpreted as nonconvergence. Indeed, at this point, the particles would only adjust their position in a small neighborhood of their last feasible visited position. By naming g best the best among all g i,p found at the end of the iterations, this fact is confirmed by a maximum convergence quadratic error: max i,p x i,p (max iter)g best I. Equality constraints are presented in Table II. This shows that good performance is achieved since the ideal ratio P/P d pso , and Q/Q d pso should be 1. The total active power losses are 5.95 MW, and the total reactive power losses are 8.34 MVAR. Eventually, in order to evaluate the optimality of the given solution, we compare it with the solution to the OPF problem having another common cost function. For instance, let us consider min U i ∑ N i=1 (U i -1) 2 , i.e. we want to keep voltages at every bus as close as possible to 1 p.u.. In this case, the active power losses amount to 6.57 MW, and reactive ones to 9.19 MVAR, which shows the effectiveness of the optimal solution of Table I. VII. CONCLUSIONS A novel distributed approach to solve OPF was presented. Being based on a population algorithm, no assumptions were made concerning either the cost function of the optimization problem, and the network topology. Moreover the presented algorithm enables to solve general distributed optimization problems that can be written in the form of (1). Potentiality of the algorithm was also shown when solving a distributed factorization problem to tune the finite-time average consensus algorithm.

The main drawback of the approach is due to the communication burden, in turns due to the need for performing a finite-time average consensus algorithm for each step of DPSO. This basically raises the time for convergence especially for communication graph having high diameter, and radius. We are currently investigating an algorithm modification to avoid performing the average consensus step.

Fig. 1 .

 1 Fig. 1. Communication graphs example.

Fig. 2 .

 2 Fig. 2. Bucky communication graph.

Fig. 3 .

 3 Fig. 3. Finite-time average consensus. (up): combined DPSO-gradient back-propagation solution; (down): gradient back-propagation solution.

Fig. 4 . 5 -

 45 Fig. 4. 5-bus electric grid, and its communication graph G c .

2 3 •

 3 10 -7 , which shows convergence of the algorithm. Optimal variables are shown in Table

Tarifs d'Utilisation du Réseau Public de Distribution d' Électricité

ACKNOWLEDGMENT

This study has been carried out in the RISEGrid Institute (www.supelec.fr/342p38091/risegrid-en.html), joint program between CentraleSupélec and EDF ('Electricité de France') on smarter electric grids.