Distributed Integral Action: Stability Analysis and Frequency Control of ...
Distributed Integral Action: Stability Analysis and Frequency Control of Power Systems Martin Andreasson∗† , Henrik Sandberg∗ , Dimos V. Dimarogonas∗ and Karl H. Johansson∗ ∗ ACCESS
Linnaeus Center, KTH Royal Institute of Technology, Stockholm, Sweden.
Abstract—This paper analyzes distributed proportionalintegral controllers. We prove that integral action can be successfully applied to consensus algorithms, where attenuation of static disturbances is achieved. These control algorithms are applied to decentralized frequency control of electrical power systems. We show that the proposed algorithm can attenuate step disturbances of power loads. We provide simulations of the proposed control algorithm on the IEEE 30 bus test system that demonstrate its efficiency.
I. I NTRODUCTION Distributed control is in many large-scale systems the only feasible control strategy. To attenuate unknown disturbances, proportional-integral (PI) controllers are often employed in such plants. However, it is still an open problem under which conditions distributed PI-control can stabilize a certain plant in general [1]. One important and increasingly interesting class of systems that require integral action to attenuate disturbances are electric power systems. Due to disturbances in form of load and generation changes, a proportional frequency controller is not able to reach the desired reference frequency asymptotically. To eliminate static errors, integrators are employed. Due to the inherent difficulties with distributed PI control, automatic frequency control of power systems is typically carried out at two levels. In the inner control loop the frequency is controlled with a proportional controller against a dynamic reference frequency. At the outer loop, the reference frequency is controlled with a PI controller to eliminate static errors [2]. While this control architecture often works satisfactorily today, future power system developments might render it unsuitable. For instance, large scale penetration of renewable power generation increases generation fluctuations, creating a need for fast and local disturbance attenuation. Decentralized control of power systems might also provide efficient anti-islanding control and self-healing features, even when communication is unavailable [3], [4]. In this paper we propose a fully decentralized frequency control algorithm for electrical power systems. We model the power system by the swing equation, used in transient stability [5] and fault detection [6]. We apply and further develop the theory of distributed consensus to prove the stability of This work was supported in part by the European Commission by the Hycon2 project, the Swedish Research Council (VR) and the Knut and Alice Wallenberg Foundation. The 3rd author is also affiliated with the Centre for Autonomous Systems at KTH and is supported by the VR 2009-3948 grant. † Corresponding author. E-mail: [email protected]
the proposed algorithms. These results could be of interest not only for power systems, but also for other distributed or decentralized dynamical systems with first- or second-order dynamics. We give a short review of the consensus problem to position our theoretical results with respect to the literature. The coordination of autonomous agents based solely on local interactions and decentralized control algorithms [7], has applications in multi-vehicle control [8], formation-control [9], flocking [10], [11], rendezvous [12] and distributed estimation [13], [14], [15]. Consensus with integral action has been studied in [16] for agents with single integrator dynamics. It was shown that the proposed consensus protocol can attenuate constant and some time-varying disturbances to a certain degree. The main theoretical contributions of this paper is the extension of consensus for agents with single integrator dynamics under integral action, and consensus for agents with damped double integrator dynamics under integral action. We study consensus algorithms both with and without absolute position measurements. The rest of this paper is organized as follows. In section II we introduce the mathematical notation. In section III we analyze consensus protocols with PI controllers for agents with single integrator dynamics. In section IV we propose two consensus protocols for agents with double integrator dynamics, one without and one with absolute position measurements. We also provide av example to support our theoretical results. In section V, based on the developed theory, we propose a decentralized control algorithm for frequency control of electrical power systems, and compare the performance to traditional control algorithms. The paper ends by some concluding remarks in section VI. II. P RELIMINARIES Let G be a graph. Denote by V = {1, . . . , n} the vertex set of G, and by E = {1, . . . , m} the edge set of G. Let Ni be the set of neighboring nodes to i. Denote by B = B(G) the vertex-edge adjacency matrix of G, and let L be the Laplacian matrix of G. In this paper we will only consider static, undirected and connected graphs. For the application of frequency control of power systems, this is a reasonable assumption as long as there are no power line failures. For undirected graphs it holds that L = BB T . Let C− denote ¯ − its closure. We the open left half complex plane, and C will denote the position of agent i as xi , and its velocity as vi , and collect them into column vectors x = (x1 , . . . , xn )T ,
v = (v1 , . . . , vn )T . We denote by cn×m a vector or matrix of dimension n × m whose elements are all c. In denotes the identity matrix of dimension n. III. D ISTRIBUTED I NTEGRAL ACTION FOR S INGLE I NTEGRATOR DYNAMICS A. Without absolute position measurements The objective is to reach asymptotic consensus on the states xi , i.e., limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V. We consider the following consensus protocol for agents with single integrator dynamics: (1)
z˙i = xi x˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj ) + di
(2)
where α ∈ R+ , β ∈ R+ and di ∈ R is an unknown scalar disturbance. zi represents the local integral state of agent i’s position. This consensus protocol is different from the one proposed in [16], as it does not require absolute position measurements. Furthermore, the above consensus protocol does not require communication of the integral state between the agents, as it suffices for each agent to measure its neighbors states. Theorem 1. Under the dynamics (1)–(2) we have that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V for any constant disturbances di . If di = 0P ∀i ∈ V, the agents converge to the common value x∗ = n1 i∈V xi (0). Proof: We first consider the case where di = 0 ∀i ∈ V. By introducing the state vector z = [z1 , . . . , zn ]T we may rewrite (1) – (2) in vector form. Since we are only interested in analyzing the x-dynamics, we may write z˙ 0n×n In z = . x˙ −αL −βL x {z } | ,A
By elementary column operations we note that the characteristic polynomial of A is given by 0 = det (sIn ) · det (βs + α)L + s2 In The first factor yields the multiple eigenvalue s = 0. By comparing the second factor with the characteristic polynomial of L: 0 = det (L − κIn ) with solutions κ = λi ≥ 0, we obtain the equation 2
0 = s + λi βs + λi α. The above equation has two solutions s = 0 if λi = 0, and solutions s ∈ C− if λi > 0. Since the above equation has exactly two solutions for every λi , it follows that the algebraic multiplicity of the eigenvalue 0 must be two. It is well-known that for connected graphs G, λ1 is the only zero-eigenvalue of the Laplacian L. By straightforward calculations we obT tain the that e11 = 11×n 01×n is an eigenvector and
T e21 = 01×n 11×n is a generalized eigenvector of A corresponding to the eigenvalue 0. It can also be verified that v1 = n1 11×n 01×n and v2 = n1 01×n 11×n are a generalized left eigenvector and an eigenvector of A corresponding to the eigenvalue 0, and that v1 e11 = 1, v2 e21 = 1 and v2 e11 = 0, v1 e21 = 0. If we let P be an orthonormal matrix consisting of the normalized eigenvectors of A, we can chose the first columns of P to be e11 and e21 , and the first rows of P −1 to be v1 and v2 respectively. Since all other eigenvalues of A have strictly negative real part we obtain lim eAt = lim P eJt P −1 t→∞ 1 t 01×(2n−2) −1 0 1 01×(2n−2) = P lim P 0 t→∞ 0(2n−2)×1 0(2n−2)×1 eJ t 1 t 01×(2n−2) 0 1 01×(2n−2) P −1 = lim P t→∞ 0(2n−2)×1 0(2n−2)×1 0(2n−2)×(2n−2) 1 1n×n t1n×n = lim 1n×n t→∞ n 0n×n t→∞
Given an initial position x(0) = x0 , we obtain 1X lim xi (t) = x0,i ∀ i ∈ V t→∞ n i∈V
i.e. the agents converge to the average of their initial positions. We now turn our attention to the case where d = [d1 , . . . , dn ]T is a constant, nonzero disturbance. Let us define the output of the system as # " yz z BT 0(n−1)×n · = x yx 0(n−1)×n BT {z } | ,C
Let us consider the linear coordinate change # " h i √1 11×n 1 n×1 n S u x = √n 1 x u= ST " # i h √1 11×n 1 n×1 n S w z = √n 1 w= z ST h i where S is a matrix such that √1n 1n×1 S is an orthonormal matrix. In the new coordinates, the system dynamics (1)–(2) become w˙ = u " u˙ =
0(n−1)×1 "
+
0
0 0(n−1)×1
# 01×(n−1) w −αS T LS # " # 1 01×(n−1) 11×(n) n u+ d −βS T LS ST
We note that the states u1 and w1 are both unobservable and uncontrollable. We thus omit these states to obtain a minimal realization by defining the new coordinates u0 =
T T u2 , . . . , u n and w0 = w2 , . . . , wn , we obtain the system dynamics # " # 0 " (n−1)×(n−1) w˙ w0 ST d 0 I(n−1) = + u˙ 0 0(n−1)×1 −αS T LS −βS T LS u0 By Lemma 10 in [16], S T LS is invertible and we may define # 00 0 " 1 T w w (S LS)−1 S T d α = 0 − u00 u 0(n−1)×1 It is easily verified that 0(2n−2)×1 is the only equilibrium of the system dynamics, which in the new coordinates are given by # 00 " 0(n−1)×(n−1) I(n−1) w00 w˙ 00 = T T u˙ −αS LS −βS LS u00 | {z } ,A00
By [17], the eigenvalues of A00 are given by det κ2 I(n−1) + (βκ + α)S T LS
Theorem 2. Under the dynamics (3)–(4), the agents converge to a common value x∗ for any constant disturbance di . If P di = 0 ∀i ∈ V, the agents converge to x∗ = n1 i∈V xi (0). Proof: If di = xi (0) = 0 ∀i ∈ V, the system dynamics (1) – (2) may be written in vector form as z z˙ 0 In = n×n x˙ −αL −βL − δI x {z } | ,A
By elementary column operations, the characteristic polynomial of A may be written as 0 = det (sIn ) det (βs + α)L + (s2 + δs)In By similar arguments used in the proof of Theorem 1, A has a simple eigenvalue 0, T with the corresponding eigen1 0 vector e = and the left eigenvector v1 = 1×n 1×n 1 1 1 0 , whereas all other eigenvalues have negative 1×n 1×n n real part. We see that v1 e1 = 1, and hence it follows that lim eAt = lim P eJt P −1 = # " t→∞ 1 1n×n 1 01×(2n−1) −1 0 P = P lim t→∞ 0(2n−1)×1 eJ t n 0n×n
t→∞
By comparing this with the characteristic equation det sI + B T LB , which by lemma 10 in [16], has solutions −si < 0, we know that the eigenvalues of A00 must satisfy
lim x(t) = 0
t→∞
with solutions κ ∈ C− . Thus A00 is Hurwitz. B. With absolute position measurements Note that Theorem 1 only guarantees relative consensus, i.e. limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ G under disturbances. It does not guarantee any properties of the agents absolute positions. As a matter of fact limt→∞ |xi (t)| = ∞ ∀i ∈ G under the consensus protocol (1)–(2) in general. This can be very undesirable in many applications, such as sensor fusion and rendezvous of mobile robots. The goal might be to have both limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V and xi (t) → x∗ ∀i ∈ V for any disturbance d, for some x∗ ∈ R. In order to satisfy these requirements, we introduce absolute position measurements to the consensus protocol. We thus consider the following consensus protocol for agents with single integrator dynamics:
x˙ i = −
(3) X j∈Ni
Given any initial position x(0) = x0 , it immediately follows that
κ2 + si βκ + si α = 0
z˙i = xi
0n×n 0n×n
α(zi − zj ) + β(xi − xj ) − δ(xi − xi (0)) + di (4)
where α ∈ R+ , β ∈ R+ , δ ∈ R+ and di ∈ R is an unknown scalar disturbance. Remark 1. Note that we require each agent to know its initial position xi (0). This will allow the agents to converge to their initial average positions under certain cond
Now let di 6= 0 and x(0) 6= 0 for at least one i ∈ V. Defining the output as # " yz BT 0(n−1)×n z = x yx 0(n−1)×n In {z } | ,C
the proof follows analogously to the proof of the second part of theorem 1. Finally, if di = 0 ∀i ∈ V, the stationarity of x(t) implies: lim 11×n −αLz(t) − βLx(t) − δx(t) + δx(0) = 0 t→∞ X ⇒ nx∗ = xi (0) i∈V
which concludes the proof. IV. D ISTRIBUTED I NTEGRAL ACTION FOR DOUBLE INTEGRATOR DYNAMICS WITH DAMPING
A. Without absolute position measurements We consider the following consensus protocol for agents with double integrator dynamics: (5)
z˙i = xi
(6)
x˙ i = vi v˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj ) − γvi + di
(7)
where α ∈ R+ , β ∈ R+ , γ ∈ R+ and di ∈ R is an unknown scalar disturbance.
Theorem 3. Under the dynamics (5)–(7) we have that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V for any constant disturbance di if and only if α < βγ. If di = P0 ∀i 1∈ V, the ∗ agents converge to the common value x = i∈V n xi (0) + P 1 v (0). i i∈V γn Proof: We first consider the case where δi = 0 ∀i ∈ V. Let also di = 0 ∀i ∈ V. By introducing the state vector z = [z1 , . . . , zn ]T we may rewrite (5) – (7) in vector form as: z˙ z 0n×n In 0n×n x˙ = 0n×n 0n×n In x v˙ v −αL −βL −γIn {z } | ,A
By elementary column operations it is easily shown that the characteristic polynomial of A can be written as 0 = det (sI) · det s2 I · det (α + βs)L + s2 (s + γ)I where I is the identity matrix of appropriate dimensions. The first two factors obviously give the solutions s = 0. Comparing the third factor with the characteristic polynomial of L, we get
We define the output of the system as BT 0(n−1)×n 0(n−1)×n z yz yx = BT 0(n−1)×n x 0(n−1)×n v yv 0(n−1)×n 0(n−1)×n BT | {z } ,C
Now let us consider the linear coordinate change " # i h √1 11×n 1 n×1 0 0 n √ 1 S z z = z= z n ST " # i h √1 11×n 1 n×1 0 0 n S x x = x = √n 1 x ST " # i h √1 11×n 1 n×1 0 0 n S v v = v = √n 1 v ST i h where S is a matrix such that √1n 1n×1 S is an orthonormal matrix. In the new coordinates the system dynamics are z˙ 0 = x x˙ 0 = v " 0
v˙ =
0 = s3 + γs2 + λi βs + λi α where λi is an eigenvalue of L. If λi > 0, the above equation has all its solutions s ∈ C− if and only if α < βγ, and α, β, γ > 0 by the Routh-Hurwitz stability criterion. Since G by assumption is connected, λ1 = 0 and λi > 0 ∀i = 2, . . . , n. For λ1 = 0, the above equation has the solutions s = 0, s = −γ. Bystraightforward calculations T it can be shown that e11 = 11×n 01×n 01×n and T 2 e1 = 01×n 11×n 01×n are an eigenvector and a generalized eigenvector corresponding to the eigenvalue 0. Furthermore v1 = γ 21n γ 2 11×n 01×n 11×n and v2 = 1 01×n γ11×n 11×n are a generalized left eigenvector γn and a left eigenvector of A corresponding to he eigenvalue 0. Furthermore v1 e11 = 1, v2 e21 = 1 and v2 e11 = 0, v1 e21 = 0. Hence the first columns of P can be chosen as e11 and e21 , and the first rows of P −1 can be chosen to be v1 and v2 . Since all other eigenvalues of A have strictly negative real part we obtain lim eAt = lim P eJt P −1 t→∞ t→∞ 1 t 01×(3n−2) −1 0 1 01×(3n−2) = P lim P t→∞ J 0t 0(3n−2)×1 0(3n−2)×1 e 1 t1n×n 1+tγ γ 2 1n×n 1 n×n 1 = lim 0n×n 1n×n γ 1n×n t→∞ n 0n×n 0n×n 0n×n Given any initial position x(0) = x0 , v(0) = v0 , we obtain 1X 1 X lim xi (t) = x0,i + v0,i ∀ i ∈ V t→∞ n γn i∈V
i∈V
0
0(n−1)×1 "
+
0 (n−1)×1
0
# 01×(n−1) 0 z −αS T LS # " # 1 01×(n−1) 11×(n) 0 0 n x − γv + d ST −βS T LS
We note that the states z10 , x01 and v10 are unobservable and uncontrollable. We thus omit these states to obtain a minimal realization by new coordinates z 00 = 0 0defining0 the 0 T T 0 T 00 00 z2 , . . . , zn , x = x2 , . . . , xn and v = v2 , . . . , vn0 we obtain the system dynamics 00 0(n−1)2 I(n−1)2 0(n−1)2 z˙ z 00 00 x˙ 00 = 0(n−1)2 I(n−1)2 x 0(n−1)2 v˙ 00 v 00 −αS T LS −βS T LS −γI(n−1)2 | {z } ,A00
ST d + 0(n−1)×1 0(n−1)×1 We now shift the state space by defining: 1 z (3) z 00 (S T LS)−1 S T d α (3) 00 0(n−1)×1 x = x − v 00 0(n−1)×1 v (3) It is easily verified that 0(2n−2)×1 is the only equilibrium of the system dynamics, and that the dynamics in the new coordinates are also characterized by the matrix A00 . By a similar argument used when showing that A is negative semidefinite, we may show that A00 is negative semi-definite. But since S T LS is full-rank, A00 must also be full-rank, and hence A00 is Hurwitz. This implies that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V even in the presence of constant disturbances di . It is also clear that whenever α ≥ βγ, at
least one eigenvalue will have non-negative real part, and that its (generalized) eigenvector will be distinct from e11 and e21 . This implies that consensus cannot be reached. Remark 2. If additionally di = vi (0) = 0 ∀i ∈ V, the agentsPconverge to the average of their initial positions, i.e. x∗ = i∈V n1 xi (0). B. With absolute position measurements As for the consensus algorithm for single integrator dynamics proposed in section III-A, in general we have limt→∞ |xi (t)| = ∞ ∀i, j ∈ V under the consensus protocol (5)–(7) We now propose a consensus algorithm for double integrator dynamics that guarantees bounded consensus even under disturbances, by introducing absolute position measurements. Moreover, with the proposed consensus algorithm the agents reach the average of their initial positions for arbitrary initial velocities, in the absence of disturbances. Consider the following consensus protocol for agents with double integrator dynamics: (8)
z˙i = xi
(9)
x˙ i = vi v˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj )
− γvi − δ(xi − xi (0)) + di
(10)
where α ∈ R+ , β ∈ R+ , γ ∈ R+ , δ ∈ R+ and di ∈ R is an unknown scalar disturbance. Theorem 4. Under the dynamics (8)–(10), the agents converge to a common value x∗ for any constant disturbance di . P If di = 0 ∀i ∈ V, the agents converge to x∗ = n1 i∈V xi (0) for arbitrary vi (0). Proof: If di = xi (0) = 0 ∀i ∈ V, the dynamics (8)–(10) can be written in vector form as z˙ 0n×n In 0n×n z x˙ = 0n×n 0n×n In · x v˙ −αL −βL − δI −γIn v | {z }
with one solution s = 0, and two solutions s ∈ C− . As all other eigenvalues are strictly positive, the solutions s corresponding to strictly positive λi satisfy s ∈ C− iff λi α + δ < λi βγ. This is satisfied if α < T βγ. It can be verified that e1 = 11×n 01×n 01×n and v1 = T 1 δ11×n γ11×n 11×n are a right and left eigenvecn tor of A, corresponding to the eigenvalue 0. Furthermore v1 e1 = 1. Since all other eigenvalues have strictly negative real part, we have lim eAt = lim P eJt P −1 t→∞ " # 1 01×(3n−1) 0 = P lim P −1 t→∞ 0(3n−1)×1 eJ t δ1 γ1n×n 1n×n 1 n×n 0n×n 0n×n 0n×n = n 0n×n 0n×n 0n×n
t→∞
Given any initial position x(0) = x0 , v(0) = v0 , we obtain lim x(t) = 0
t→∞
Defining the output of the system as z yz BT 0(n−1)×n 0(n−1)×n yx = BT 0(n−1)×n · x 0(n−1)×n v yv 0n×n 0n×n In {z } | ,C
the proof follows analogously to the proof of the second part of theorem 3. Finally, assume di = 0 ∀i ∈ V. The stationarity of v(t) implies: lim 11×n −αLz − βLx − δx + δx(0) − γv = 0 t→∞ X ⇒ nx∗ = xi (0) i∈V
which concludes the proof.
Example 1. Let the dynamics of the agents be given by (5)– (7) with β = 5 and γ = 3. The setup we will consider consists By performing elementary column operations the character- of a string of 5 mobile robots, whose communication topology istic polynomial of A can be written as is a string graph. We consider the system with a constant disturbance d = [1, 0, 0, 0, 0], for α = 1 and for α = 0, i.e. 0= with and without integral action. The initial conditions are det (sI) det (s + γ)I det (α + βs)L + (s3 + γs2 + δs)I given by x(0) = [5, −6, 8, 4, 5], v(0) = [0, 0, 0, 0, 0]T . By Theorem 3 stability is guaranteed when α > 0. In The first two factors give the solutions s = 0 and s = −γ Figure 1 the state trajectories are shown for both α = 0 and respectively. By comparing the third factor with the characterα = 1. We observe that asymptotic consensus is only reached istic polynomial of L, it can be seen that the other eigenvalues when α = 1. We conclude that introducing integral action is s of A satisfy important for reaching consensus under disturbances. 3 2 0 = s + γs + (δ + λi β)s + λi α V. F REQUENCY CONTROL OF POWER SYSTEMS where λi ∈ spec(L) Since G by assumption is connected, In this section we show that a similar protocol to the one L has a single simple eigenvalue λ1 = 0, which gives the proposed in section IV-B, can be employed for automatic characteristic polynomial frequency control in power systems. Let us consider a power 0 = s(s2 + γs + δ) system modeled by a graph G = (V, E). Each node, here ,A
frequency is often measured at a specific bus. This will typically lead to longer delays, since disturbances need to propagate through the system before control action can be taken. The centralized controller architecture is illustrated in figure 2.
x(t)
5 0 −5
0
5
10 t
15
20
+ ω
ref
Σ
x(t)
5
ω ˆ
Cc
0
C1
−
Σ 0
5
10 t
15
20
Figure 1. The upper figure shows the state trajectories of (5)–(7) when α = 0, and the lower figure shows the state trajectories when α = 1.
referred to as a bus, is assumed to obey the linearized swing equation [2] X mi δ¨i + di δ˙i = − kij (δi − δj ) + pm (11) i + ui
Figure 2.
ω1 Bus 1
.. . +
−5
u1
Cn
−
un
Bus n
ωn + +
Σ
The figure illustrates the centralized control architecture.
ω ref
u1
C10
Bus 1
ω1 , δ1
j∈Ni
where δi is the phase angle of bus i, mi and di are the inertia and damping coefficient respectively, pm i is the power load at bus i and ui is the mechanical input. kij = |Vi ||Vj |bij , where Vi is the voltage of bus i, and bij is the susceptance of the edge (i, j), here referred to as line. By defining the state vectors δ = δi , . . . , δn and ω = δ˙ = ω1 . . . , ωn , we may rewrite (11) in state-space form as " # 0n×n In δ 0 0n×1 δ˙ = + n×1 + (12) −M Lk −M D ω M pm Mu ω˙ where M = diag( m11 , . . . , m1n ), D = diag(d1 , . . . , dn ), Lk is the weighted Laplacian with edge weights kij , pm = m T T p1 , . . . , p m , u = ui , . . . , u n . n A. Centralized PI control Traditionally, automatic frequency control of power systems is carried out at two levels [2]. In the inner control loop the frequency is controlled with a proportional controller against a reference frequency. At the outer loop, the reference frequency is controlled with a PI controller to eliminate static errors. We model the inner controller of bus i as ui = α(ˆ ω − ωi )
and the outer PI controller is assumed to be given by 1X ωi zˆ˙ = ω ref − n i∈V 1X ω ˆ˙ = ω ref − ωi + β zˆ n
(13)
(14) (15)
i∈V
where we have assumed that the average frequency of the buses is measured by the central controller. In reality the
.. . ω ref
Figure 3.
un
Cn0
Bus n
ωn , δn
The figure illustrates the decentralized control architecture.
Theorem 5. The power system described by (12) where ui is given by (13)–(15), is stable for α > 0, β > 0, γ > 0 satisfying βγ < α. Furthermore limt→∞ ωi (t) = ω ref . Proof: As we are only interested in studying ω, we may rewrite (12) with u given by (13)–(15), as ω ˆ ω ˆ˙ − nγ 11×n 0 − nβ 11×n z · = 0n×n In z˙ 0n×1 ω ω˙ α1n×1 −M L −M D − αIn | {z } ,A
ref
γω + −ω ref 1n×1 m p
(16)
We now consider the the matrix A0 defined as 0 − nβ 11×n − nγ 11×n A0 , 0n×1 0n×n In α1n×1 −M L −αIn By performing elementary column operations on A0 , we may
write the characteristic polynomial of A0 as s3 γ + s2 α + sα + αβ · det M L + (s2 + sα) = 0 The first factor has solutions s ∈ C− if βγ < α. Comparing the second factor with the characteristic polynomial of M L, we have that s must satisfy s2 + sα + ti = 0 where ti ≥ 0 is an eigenvalue of M L. Also this equation has solutions s ∈ C− and one solution s = 0. However as A0 is full-rank, 0 cannot be an eigenvalue of A0 , so A0 is Hurwitz. But then clearly also A is Hurwitz. Finally, stationarity implies that the equilibrium solution must satisfy ω = ω ref 1n×1 . B. Decentralized PI control In this section we analyze a completely decentralized PI controller, where each bus controls its own frequency based only on local information. Thus, no frequency measurements need to be sent to a central controller, and there is no need to send control signals. This architecture might be favorable due to security concerns when sending unencrypted frequency measurements and control signals over large areas. Another benefit could be better performance when the tripping of one or several power lines causes the network to be split up into two or more sub-networks, so called islanding. The controller of node i is assumed to be given by z˙i = (ω ref − ωi )
ui = α(ω ref − ωi ) + βzi
(17) (18)
The controller architecture is illustrated in Figure 3. This decentralized control architecture is not practically feasible with only frequency measurements available at the generation buses. Even the slightest frequency measurement error will be integrated and prevent the frequencies from reaching consensus [2]. With recent advances in phasor measurement unit (PMU) technology however, phase measurements are likely to be available to all generator buses [18]. By employing optimal PMU placement, the number of PMU:s needed for complete observability can be drastically reduced [19]. By integrating (17) we obtain zi = ω ref t − δi
This implies that in order to accurately estimate the integral state zi , each generator bus needs access only to accurate time and phase measurements, both provided by PMU:s. The decentralized controller architecture is illustrated in figure 3. Theorem 6. The power system described by (12) where ui is given by (17)–(18), is stable for any choice of α > 0, β > 0. Furthermore limt→∞ ωi (t) = ω ref . Proof: As we are only interested in studying ω, we may rewrite (12) with ui is given by (17)–(18), as " ref # z˙ 0n×n In z −ω = · + (19) ω˙ −M L − βIn×n −M D ω M pm | {z } ,A
where the i’th element of β is βi . Let m = mini mi and d = mini di . We can now write: M D = m dIn + D0 where D0 is a diagonal matrix with non-negative entries. We now consider the the matrix A0 defined as 0n×n In 0 A , −M L − βIn×n −m · dIn By [17] the eigenvalues of A0 are given by det (s2 + sm d)In + M L + βIn×n
Comparing this with the characteristic polynomial of M L + diag(β), we conclude that s must satisfy s2 + sm d + ti = 0 where ti ≥ 0 is an eigenvalue of M L + diag(β). By [20] M L + diag(β) is positive definite, and hence the above equation has all its solutions in C− . It follows that also A is Hurwitz. Now consider the coordinate shift 0 z z z = − o ω0 ω ωo where z0 = (βIn + M L)−1 (M Dω ref + pm ), ω0 = ω ref 1n×1 . In the translated coordinates, 02n×1 is the only equilibrium of the system. Hence limt→∞ ωi (t) = ω ref ∀i ∈ V. C. Simulations The centralized and decentralized frequency control algorithms were tested on the IEEE 30 bus test system [21].The line admittances were extracted from [21] and the voltages were assumed to be 132 kV for all buses. M and D were set to reasonable numerical values. The power system is initially in an operational equilibrium, until the power load is increased by a step of 200 kW in the buses 2, 3 and 7. This will immediately result in decreased frequencies at the extra load buses. The frequency controllers at the buses will then control the frequencies towards the desired frequency of ω ref = 50 Hz. In the centralized controller the parameters were set to α = 0.8, β = 0, γ = 0.04, while in the decentralized control architecture the parameters were α = 0.8, β = 0.04. The step responses of the frequencies is plotted in figures 4–5. We note that if there is a centralized PI-controller for the reference frequency, the generation is increased uniformly among the generators. If however the integral action is distributed amongst the generators, some generators will increase their generation more than others. VI. D ISCUSSION AND C ONCLUSIONS In this paper we have studied consensus protocols with integral action for agents with double and single integrator dynamics. We have proved that the proposed consensus protocols reach asymptotic consensus even in the presence of constant disturbances. If we allow for absolute position measurements, the proposed consensus protocol converges asymptotically to a common state. In the absence of disturbances, the proposed consensus protocols asymptotically
φ(t) [Hz]
systems, with reasonable performance and possible benefits over centralized frequency control.
50
R EFERENCES
49.95 49.9
0
20
40
60
80
100
P (t) [kW]
t [s] 60 40 20 0
0
20
40
60
80
100
t [s]
φ(t) [Hz]
Figure 4. The upper figure shows the bus frequencies with centralized frequency control, while the lower figure shows the control signals at all buses. We observe that the frequencies converge quickly to a common frequency. The frequencies are then slowly regulated towards the reference frequency by the central PI controller.
50 49.95 49.9
0
20
40
60
80
100
60
80
100
P (t) [kW]
t [s] 60 40 20 0
0
20
40 t [s]
Figure 5. The upper figure shows the bus frequencies with decentralized frequency control, while the lower figure shows the control signals at all buses. We observe that the frequencies converge quickly to a common frequency. The frequencies are then slowly regulated towards the reference frequency by the PI controllers present at each bus.
solve the initial average consensus problem. Furthermore we have demonstrated that a similar consensus protocol can be employed for decentralized frequency control of power
[1] M. Morari and E. Zafiriou. Robust Process Control. Prentice Hall, Englewood Cliffs, 1989. [2] J. Machowski, J.W. Bialek, and J.R. Bumby. Power System Dynamics: Stability and Control. Wiley, 2008. [3] N. Senroy, G.T. Heydt, and V. Vittal. Decision tree assisted controlled islanding. Power Systems, IEEE Transactions on, 21(4):1790 –1797, nov. 2006. [4] Bo Yang, V. Vittal, and G.T Heydt. Slow-coherency-based controlled islanding: A demonstration of the approach on the august 14, 2003 blackout scenario. Power Systems, IEEE Transactions on, 21(4):1840– 1847, nov. 2006. [5] Florian Doerfler and Francesco Bullo. on the critical coupling for kuramoto oscillators. SIAM Journal on Applied Dynamical Systems, 10(3):1070–1099, 2011. [6] Iman Shames, Andr M.H. Teixeira, Henrik Sandberg, and Karl H. Johansson. Distributed fault detection for interconnected second-order systems. Automatica, 47(12):2757 – 2764, 2011. [7] R. Olfati-Saber and R.M. Murray. Consensus problems in networks of agents with switching topology and time-delays. IEEE Transactions on Automatic Control, 49(9):1520 – 1533, sept. 2004. [8] W. Ren and R.W. Beard. Distributed Consensus in Multi-vehicle Cooperative Control. Springer, 2008. [9] J.A. Fax and R.M. Murray. Information flow and cooperative control of vehicle formations. IEEE Transactions Automatic Control, 49(9):1465 – 1476, sept. 2004. [10] Liu, Y. and Passino, K.M. and Polycarpou, M.M. Stability analysis of w-dimensional asynchronous swarms with a fixed communication topology. IEEE Transactions Automatic Control, 48:76–95, 2003. [11] H.G. Tanner, A. Jadbabaie, and G.J. Pappas. Stable flocking of mobile agents, part i: Fixed topology. In IEEE Conference on Decision and Control, volume 2, pages 2010 – 2015 Vol.2, dec. 2003. [12] F. Bullo, J. Cort´es, and S. Mart´ınez. Distributed Control of Robotic Networks. Applied Mathematics Series. Princeton University Press, 2009. [13] L. Xiao, S. Boyd, and S. Lall. A scheme for robust distributed sensor fusion based on average consensus. In International Symposium on Information Processing in Sensor Networks, pages 63 – 70, april 2005. [14] R. Olfati-Saber. Distributed Kalman filter with embedded consensus filters. In IEEE Conference on Decision and Control, European Control Conference, pages 8179 – 8184, dec. 2005. [15] R. Carli, A. Chiuso, L. Schenato, and S. Zampieri. Distributed Kalman filtering based on consensus strategies. IEEE Journal on Selected Areas in Communications, 26(4):622 –633, may 2008. [16] R.A. Freeman, Peng Yang, and K.M. Lynch. Stability and convergence properties of dynamic average consensus estimators. In IEEE Conference on Decision and Control, pages 338 –343, dec. 2006. [17] John R. Silvester. Determinants of block matrices. The Mathematical Gazette, 84(501):pp. 460–467, 2000. [18] A.G. Phadke. Synchronized phasor measurements in power systems. Computer Applications in Power, IEEE, 6(2):10 –15, april 1993. [19] R.F. Nuqui and A.G. Phadke. Phasor measurement unit placement techniques for complete and incomplete observability. Power Delivery, IEEE Transactions on, 20(4):2381 – 2388, oct. 2005. [20] Andr´e Teixeira. Multi- agent systems with fault and security constraints. Master’s thesis, Royal Institute of Technology (KTH), June 2009. [21] Power systems test case archive-30 bus power flow test case. available in: http://www.ee.washington.edu/research/pstca/pf30/pgtca30bus.htm.
Linnaeus Center, KTH Royal Institute of Technology, Stockholm, Sweden.
Abstract—This paper analyzes distributed proportionalintegral controllers. We prove that integral action can be successfully applied to consensus algorithms, where attenuation of static disturbances is achieved. These control algorithms are applied to decentralized frequency control of electrical power systems. We show that the proposed algorithm can attenuate step disturbances of power loads. We provide simulations of the proposed control algorithm on the IEEE 30 bus test system that demonstrate its efficiency.
I. I NTRODUCTION Distributed control is in many large-scale systems the only feasible control strategy. To attenuate unknown disturbances, proportional-integral (PI) controllers are often employed in such plants. However, it is still an open problem under which conditions distributed PI-control can stabilize a certain plant in general [1]. One important and increasingly interesting class of systems that require integral action to attenuate disturbances are electric power systems. Due to disturbances in form of load and generation changes, a proportional frequency controller is not able to reach the desired reference frequency asymptotically. To eliminate static errors, integrators are employed. Due to the inherent difficulties with distributed PI control, automatic frequency control of power systems is typically carried out at two levels. In the inner control loop the frequency is controlled with a proportional controller against a dynamic reference frequency. At the outer loop, the reference frequency is controlled with a PI controller to eliminate static errors [2]. While this control architecture often works satisfactorily today, future power system developments might render it unsuitable. For instance, large scale penetration of renewable power generation increases generation fluctuations, creating a need for fast and local disturbance attenuation. Decentralized control of power systems might also provide efficient anti-islanding control and self-healing features, even when communication is unavailable [3], [4]. In this paper we propose a fully decentralized frequency control algorithm for electrical power systems. We model the power system by the swing equation, used in transient stability [5] and fault detection [6]. We apply and further develop the theory of distributed consensus to prove the stability of This work was supported in part by the European Commission by the Hycon2 project, the Swedish Research Council (VR) and the Knut and Alice Wallenberg Foundation. The 3rd author is also affiliated with the Centre for Autonomous Systems at KTH and is supported by the VR 2009-3948 grant. † Corresponding author. E-mail: [email protected]
the proposed algorithms. These results could be of interest not only for power systems, but also for other distributed or decentralized dynamical systems with first- or second-order dynamics. We give a short review of the consensus problem to position our theoretical results with respect to the literature. The coordination of autonomous agents based solely on local interactions and decentralized control algorithms [7], has applications in multi-vehicle control [8], formation-control [9], flocking [10], [11], rendezvous [12] and distributed estimation [13], [14], [15]. Consensus with integral action has been studied in [16] for agents with single integrator dynamics. It was shown that the proposed consensus protocol can attenuate constant and some time-varying disturbances to a certain degree. The main theoretical contributions of this paper is the extension of consensus for agents with single integrator dynamics under integral action, and consensus for agents with damped double integrator dynamics under integral action. We study consensus algorithms both with and without absolute position measurements. The rest of this paper is organized as follows. In section II we introduce the mathematical notation. In section III we analyze consensus protocols with PI controllers for agents with single integrator dynamics. In section IV we propose two consensus protocols for agents with double integrator dynamics, one without and one with absolute position measurements. We also provide av example to support our theoretical results. In section V, based on the developed theory, we propose a decentralized control algorithm for frequency control of electrical power systems, and compare the performance to traditional control algorithms. The paper ends by some concluding remarks in section VI. II. P RELIMINARIES Let G be a graph. Denote by V = {1, . . . , n} the vertex set of G, and by E = {1, . . . , m} the edge set of G. Let Ni be the set of neighboring nodes to i. Denote by B = B(G) the vertex-edge adjacency matrix of G, and let L be the Laplacian matrix of G. In this paper we will only consider static, undirected and connected graphs. For the application of frequency control of power systems, this is a reasonable assumption as long as there are no power line failures. For undirected graphs it holds that L = BB T . Let C− denote ¯ − its closure. We the open left half complex plane, and C will denote the position of agent i as xi , and its velocity as vi , and collect them into column vectors x = (x1 , . . . , xn )T ,
v = (v1 , . . . , vn )T . We denote by cn×m a vector or matrix of dimension n × m whose elements are all c. In denotes the identity matrix of dimension n. III. D ISTRIBUTED I NTEGRAL ACTION FOR S INGLE I NTEGRATOR DYNAMICS A. Without absolute position measurements The objective is to reach asymptotic consensus on the states xi , i.e., limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V. We consider the following consensus protocol for agents with single integrator dynamics: (1)
z˙i = xi x˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj ) + di
(2)
where α ∈ R+ , β ∈ R+ and di ∈ R is an unknown scalar disturbance. zi represents the local integral state of agent i’s position. This consensus protocol is different from the one proposed in [16], as it does not require absolute position measurements. Furthermore, the above consensus protocol does not require communication of the integral state between the agents, as it suffices for each agent to measure its neighbors states. Theorem 1. Under the dynamics (1)–(2) we have that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V for any constant disturbances di . If di = 0P ∀i ∈ V, the agents converge to the common value x∗ = n1 i∈V xi (0). Proof: We first consider the case where di = 0 ∀i ∈ V. By introducing the state vector z = [z1 , . . . , zn ]T we may rewrite (1) – (2) in vector form. Since we are only interested in analyzing the x-dynamics, we may write z˙ 0n×n In z = . x˙ −αL −βL x {z } | ,A
By elementary column operations we note that the characteristic polynomial of A is given by 0 = det (sIn ) · det (βs + α)L + s2 In The first factor yields the multiple eigenvalue s = 0. By comparing the second factor with the characteristic polynomial of L: 0 = det (L − κIn ) with solutions κ = λi ≥ 0, we obtain the equation 2
0 = s + λi βs + λi α. The above equation has two solutions s = 0 if λi = 0, and solutions s ∈ C− if λi > 0. Since the above equation has exactly two solutions for every λi , it follows that the algebraic multiplicity of the eigenvalue 0 must be two. It is well-known that for connected graphs G, λ1 is the only zero-eigenvalue of the Laplacian L. By straightforward calculations we obT tain the that e11 = 11×n 01×n is an eigenvector and
T e21 = 01×n 11×n is a generalized eigenvector of A corresponding to the eigenvalue 0. It can also be verified that v1 = n1 11×n 01×n and v2 = n1 01×n 11×n are a generalized left eigenvector and an eigenvector of A corresponding to the eigenvalue 0, and that v1 e11 = 1, v2 e21 = 1 and v2 e11 = 0, v1 e21 = 0. If we let P be an orthonormal matrix consisting of the normalized eigenvectors of A, we can chose the first columns of P to be e11 and e21 , and the first rows of P −1 to be v1 and v2 respectively. Since all other eigenvalues of A have strictly negative real part we obtain lim eAt = lim P eJt P −1 t→∞ 1 t 01×(2n−2) −1 0 1 01×(2n−2) = P lim P 0 t→∞ 0(2n−2)×1 0(2n−2)×1 eJ t 1 t 01×(2n−2) 0 1 01×(2n−2) P −1 = lim P t→∞ 0(2n−2)×1 0(2n−2)×1 0(2n−2)×(2n−2) 1 1n×n t1n×n = lim 1n×n t→∞ n 0n×n t→∞
Given an initial position x(0) = x0 , we obtain 1X lim xi (t) = x0,i ∀ i ∈ V t→∞ n i∈V
i.e. the agents converge to the average of their initial positions. We now turn our attention to the case where d = [d1 , . . . , dn ]T is a constant, nonzero disturbance. Let us define the output of the system as # " yz z BT 0(n−1)×n · = x yx 0(n−1)×n BT {z } | ,C
Let us consider the linear coordinate change # " h i √1 11×n 1 n×1 n S u x = √n 1 x u= ST " # i h √1 11×n 1 n×1 n S w z = √n 1 w= z ST h i where S is a matrix such that √1n 1n×1 S is an orthonormal matrix. In the new coordinates, the system dynamics (1)–(2) become w˙ = u " u˙ =
0(n−1)×1 "
+
0
0 0(n−1)×1
# 01×(n−1) w −αS T LS # " # 1 01×(n−1) 11×(n) n u+ d −βS T LS ST
We note that the states u1 and w1 are both unobservable and uncontrollable. We thus omit these states to obtain a minimal realization by defining the new coordinates u0 =
T T u2 , . . . , u n and w0 = w2 , . . . , wn , we obtain the system dynamics # " # 0 " (n−1)×(n−1) w˙ w0 ST d 0 I(n−1) = + u˙ 0 0(n−1)×1 −αS T LS −βS T LS u0 By Lemma 10 in [16], S T LS is invertible and we may define # 00 0 " 1 T w w (S LS)−1 S T d α = 0 − u00 u 0(n−1)×1 It is easily verified that 0(2n−2)×1 is the only equilibrium of the system dynamics, which in the new coordinates are given by # 00 " 0(n−1)×(n−1) I(n−1) w00 w˙ 00 = T T u˙ −αS LS −βS LS u00 | {z } ,A00
By [17], the eigenvalues of A00 are given by det κ2 I(n−1) + (βκ + α)S T LS
Theorem 2. Under the dynamics (3)–(4), the agents converge to a common value x∗ for any constant disturbance di . If P di = 0 ∀i ∈ V, the agents converge to x∗ = n1 i∈V xi (0). Proof: If di = xi (0) = 0 ∀i ∈ V, the system dynamics (1) – (2) may be written in vector form as z z˙ 0 In = n×n x˙ −αL −βL − δI x {z } | ,A
By elementary column operations, the characteristic polynomial of A may be written as 0 = det (sIn ) det (βs + α)L + (s2 + δs)In By similar arguments used in the proof of Theorem 1, A has a simple eigenvalue 0, T with the corresponding eigen1 0 vector e = and the left eigenvector v1 = 1×n 1×n 1 1 1 0 , whereas all other eigenvalues have negative 1×n 1×n n real part. We see that v1 e1 = 1, and hence it follows that lim eAt = lim P eJt P −1 = # " t→∞ 1 1n×n 1 01×(2n−1) −1 0 P = P lim t→∞ 0(2n−1)×1 eJ t n 0n×n
t→∞
By comparing this with the characteristic equation det sI + B T LB , which by lemma 10 in [16], has solutions −si < 0, we know that the eigenvalues of A00 must satisfy
lim x(t) = 0
t→∞
with solutions κ ∈ C− . Thus A00 is Hurwitz. B. With absolute position measurements Note that Theorem 1 only guarantees relative consensus, i.e. limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ G under disturbances. It does not guarantee any properties of the agents absolute positions. As a matter of fact limt→∞ |xi (t)| = ∞ ∀i ∈ G under the consensus protocol (1)–(2) in general. This can be very undesirable in many applications, such as sensor fusion and rendezvous of mobile robots. The goal might be to have both limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V and xi (t) → x∗ ∀i ∈ V for any disturbance d, for some x∗ ∈ R. In order to satisfy these requirements, we introduce absolute position measurements to the consensus protocol. We thus consider the following consensus protocol for agents with single integrator dynamics:
x˙ i = −
(3) X j∈Ni
Given any initial position x(0) = x0 , it immediately follows that
κ2 + si βκ + si α = 0
z˙i = xi
0n×n 0n×n
α(zi − zj ) + β(xi − xj ) − δ(xi − xi (0)) + di (4)
where α ∈ R+ , β ∈ R+ , δ ∈ R+ and di ∈ R is an unknown scalar disturbance. Remark 1. Note that we require each agent to know its initial position xi (0). This will allow the agents to converge to their initial average positions under certain cond
Now let di 6= 0 and x(0) 6= 0 for at least one i ∈ V. Defining the output as # " yz BT 0(n−1)×n z = x yx 0(n−1)×n In {z } | ,C
the proof follows analogously to the proof of the second part of theorem 1. Finally, if di = 0 ∀i ∈ V, the stationarity of x(t) implies: lim 11×n −αLz(t) − βLx(t) − δx(t) + δx(0) = 0 t→∞ X ⇒ nx∗ = xi (0) i∈V
which concludes the proof. IV. D ISTRIBUTED I NTEGRAL ACTION FOR DOUBLE INTEGRATOR DYNAMICS WITH DAMPING
A. Without absolute position measurements We consider the following consensus protocol for agents with double integrator dynamics: (5)
z˙i = xi
(6)
x˙ i = vi v˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj ) − γvi + di
(7)
where α ∈ R+ , β ∈ R+ , γ ∈ R+ and di ∈ R is an unknown scalar disturbance.
Theorem 3. Under the dynamics (5)–(7) we have that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V for any constant disturbance di if and only if α < βγ. If di = P0 ∀i 1∈ V, the ∗ agents converge to the common value x = i∈V n xi (0) + P 1 v (0). i i∈V γn Proof: We first consider the case where δi = 0 ∀i ∈ V. Let also di = 0 ∀i ∈ V. By introducing the state vector z = [z1 , . . . , zn ]T we may rewrite (5) – (7) in vector form as: z˙ z 0n×n In 0n×n x˙ = 0n×n 0n×n In x v˙ v −αL −βL −γIn {z } | ,A
By elementary column operations it is easily shown that the characteristic polynomial of A can be written as 0 = det (sI) · det s2 I · det (α + βs)L + s2 (s + γ)I where I is the identity matrix of appropriate dimensions. The first two factors obviously give the solutions s = 0. Comparing the third factor with the characteristic polynomial of L, we get
We define the output of the system as BT 0(n−1)×n 0(n−1)×n z yz yx = BT 0(n−1)×n x 0(n−1)×n v yv 0(n−1)×n 0(n−1)×n BT | {z } ,C
Now let us consider the linear coordinate change " # i h √1 11×n 1 n×1 0 0 n √ 1 S z z = z= z n ST " # i h √1 11×n 1 n×1 0 0 n S x x = x = √n 1 x ST " # i h √1 11×n 1 n×1 0 0 n S v v = v = √n 1 v ST i h where S is a matrix such that √1n 1n×1 S is an orthonormal matrix. In the new coordinates the system dynamics are z˙ 0 = x x˙ 0 = v " 0
v˙ =
0 = s3 + γs2 + λi βs + λi α where λi is an eigenvalue of L. If λi > 0, the above equation has all its solutions s ∈ C− if and only if α < βγ, and α, β, γ > 0 by the Routh-Hurwitz stability criterion. Since G by assumption is connected, λ1 = 0 and λi > 0 ∀i = 2, . . . , n. For λ1 = 0, the above equation has the solutions s = 0, s = −γ. Bystraightforward calculations T it can be shown that e11 = 11×n 01×n 01×n and T 2 e1 = 01×n 11×n 01×n are an eigenvector and a generalized eigenvector corresponding to the eigenvalue 0. Furthermore v1 = γ 21n γ 2 11×n 01×n 11×n and v2 = 1 01×n γ11×n 11×n are a generalized left eigenvector γn and a left eigenvector of A corresponding to he eigenvalue 0. Furthermore v1 e11 = 1, v2 e21 = 1 and v2 e11 = 0, v1 e21 = 0. Hence the first columns of P can be chosen as e11 and e21 , and the first rows of P −1 can be chosen to be v1 and v2 . Since all other eigenvalues of A have strictly negative real part we obtain lim eAt = lim P eJt P −1 t→∞ t→∞ 1 t 01×(3n−2) −1 0 1 01×(3n−2) = P lim P t→∞ J 0t 0(3n−2)×1 0(3n−2)×1 e 1 t1n×n 1+tγ γ 2 1n×n 1 n×n 1 = lim 0n×n 1n×n γ 1n×n t→∞ n 0n×n 0n×n 0n×n Given any initial position x(0) = x0 , v(0) = v0 , we obtain 1X 1 X lim xi (t) = x0,i + v0,i ∀ i ∈ V t→∞ n γn i∈V
i∈V
0
0(n−1)×1 "
+
0 (n−1)×1
0
# 01×(n−1) 0 z −αS T LS # " # 1 01×(n−1) 11×(n) 0 0 n x − γv + d ST −βS T LS
We note that the states z10 , x01 and v10 are unobservable and uncontrollable. We thus omit these states to obtain a minimal realization by new coordinates z 00 = 0 0defining0 the 0 T T 0 T 00 00 z2 , . . . , zn , x = x2 , . . . , xn and v = v2 , . . . , vn0 we obtain the system dynamics 00 0(n−1)2 I(n−1)2 0(n−1)2 z˙ z 00 00 x˙ 00 = 0(n−1)2 I(n−1)2 x 0(n−1)2 v˙ 00 v 00 −αS T LS −βS T LS −γI(n−1)2 | {z } ,A00
ST d + 0(n−1)×1 0(n−1)×1 We now shift the state space by defining: 1 z (3) z 00 (S T LS)−1 S T d α (3) 00 0(n−1)×1 x = x − v 00 0(n−1)×1 v (3) It is easily verified that 0(2n−2)×1 is the only equilibrium of the system dynamics, and that the dynamics in the new coordinates are also characterized by the matrix A00 . By a similar argument used when showing that A is negative semidefinite, we may show that A00 is negative semi-definite. But since S T LS is full-rank, A00 must also be full-rank, and hence A00 is Hurwitz. This implies that limt→∞ |xi (t) − xj (t)| = 0 ∀i, j ∈ V even in the presence of constant disturbances di . It is also clear that whenever α ≥ βγ, at
least one eigenvalue will have non-negative real part, and that its (generalized) eigenvector will be distinct from e11 and e21 . This implies that consensus cannot be reached. Remark 2. If additionally di = vi (0) = 0 ∀i ∈ V, the agentsPconverge to the average of their initial positions, i.e. x∗ = i∈V n1 xi (0). B. With absolute position measurements As for the consensus algorithm for single integrator dynamics proposed in section III-A, in general we have limt→∞ |xi (t)| = ∞ ∀i, j ∈ V under the consensus protocol (5)–(7) We now propose a consensus algorithm for double integrator dynamics that guarantees bounded consensus even under disturbances, by introducing absolute position measurements. Moreover, with the proposed consensus algorithm the agents reach the average of their initial positions for arbitrary initial velocities, in the absence of disturbances. Consider the following consensus protocol for agents with double integrator dynamics: (8)
z˙i = xi
(9)
x˙ i = vi v˙ i = −
X j∈Ni
α(zi − zj ) + β(xi − xj )
− γvi − δ(xi − xi (0)) + di
(10)
where α ∈ R+ , β ∈ R+ , γ ∈ R+ , δ ∈ R+ and di ∈ R is an unknown scalar disturbance. Theorem 4. Under the dynamics (8)–(10), the agents converge to a common value x∗ for any constant disturbance di . P If di = 0 ∀i ∈ V, the agents converge to x∗ = n1 i∈V xi (0) for arbitrary vi (0). Proof: If di = xi (0) = 0 ∀i ∈ V, the dynamics (8)–(10) can be written in vector form as z˙ 0n×n In 0n×n z x˙ = 0n×n 0n×n In · x v˙ −αL −βL − δI −γIn v | {z }
with one solution s = 0, and two solutions s ∈ C− . As all other eigenvalues are strictly positive, the solutions s corresponding to strictly positive λi satisfy s ∈ C− iff λi α + δ < λi βγ. This is satisfied if α < T βγ. It can be verified that e1 = 11×n 01×n 01×n and v1 = T 1 δ11×n γ11×n 11×n are a right and left eigenvecn tor of A, corresponding to the eigenvalue 0. Furthermore v1 e1 = 1. Since all other eigenvalues have strictly negative real part, we have lim eAt = lim P eJt P −1 t→∞ " # 1 01×(3n−1) 0 = P lim P −1 t→∞ 0(3n−1)×1 eJ t δ1 γ1n×n 1n×n 1 n×n 0n×n 0n×n 0n×n = n 0n×n 0n×n 0n×n
t→∞
Given any initial position x(0) = x0 , v(0) = v0 , we obtain lim x(t) = 0
t→∞
Defining the output of the system as z yz BT 0(n−1)×n 0(n−1)×n yx = BT 0(n−1)×n · x 0(n−1)×n v yv 0n×n 0n×n In {z } | ,C
the proof follows analogously to the proof of the second part of theorem 3. Finally, assume di = 0 ∀i ∈ V. The stationarity of v(t) implies: lim 11×n −αLz − βLx − δx + δx(0) − γv = 0 t→∞ X ⇒ nx∗ = xi (0) i∈V
which concludes the proof.
Example 1. Let the dynamics of the agents be given by (5)– (7) with β = 5 and γ = 3. The setup we will consider consists By performing elementary column operations the character- of a string of 5 mobile robots, whose communication topology istic polynomial of A can be written as is a string graph. We consider the system with a constant disturbance d = [1, 0, 0, 0, 0], for α = 1 and for α = 0, i.e. 0= with and without integral action. The initial conditions are det (sI) det (s + γ)I det (α + βs)L + (s3 + γs2 + δs)I given by x(0) = [5, −6, 8, 4, 5], v(0) = [0, 0, 0, 0, 0]T . By Theorem 3 stability is guaranteed when α > 0. In The first two factors give the solutions s = 0 and s = −γ Figure 1 the state trajectories are shown for both α = 0 and respectively. By comparing the third factor with the characterα = 1. We observe that asymptotic consensus is only reached istic polynomial of L, it can be seen that the other eigenvalues when α = 1. We conclude that introducing integral action is s of A satisfy important for reaching consensus under disturbances. 3 2 0 = s + γs + (δ + λi β)s + λi α V. F REQUENCY CONTROL OF POWER SYSTEMS where λi ∈ spec(L) Since G by assumption is connected, In this section we show that a similar protocol to the one L has a single simple eigenvalue λ1 = 0, which gives the proposed in section IV-B, can be employed for automatic characteristic polynomial frequency control in power systems. Let us consider a power 0 = s(s2 + γs + δ) system modeled by a graph G = (V, E). Each node, here ,A
frequency is often measured at a specific bus. This will typically lead to longer delays, since disturbances need to propagate through the system before control action can be taken. The centralized controller architecture is illustrated in figure 2.
x(t)
5 0 −5
0
5
10 t
15
20
+ ω
ref
Σ
x(t)
5
ω ˆ
Cc
0
C1
−
Σ 0
5
10 t
15
20
Figure 1. The upper figure shows the state trajectories of (5)–(7) when α = 0, and the lower figure shows the state trajectories when α = 1.
referred to as a bus, is assumed to obey the linearized swing equation [2] X mi δ¨i + di δ˙i = − kij (δi − δj ) + pm (11) i + ui
Figure 2.
ω1 Bus 1
.. . +
−5
u1
Cn
−
un
Bus n
ωn + +
Σ
The figure illustrates the centralized control architecture.
ω ref
u1
C10
Bus 1
ω1 , δ1
j∈Ni
where δi is the phase angle of bus i, mi and di are the inertia and damping coefficient respectively, pm i is the power load at bus i and ui is the mechanical input. kij = |Vi ||Vj |bij , where Vi is the voltage of bus i, and bij is the susceptance of the edge (i, j), here referred to as line. By defining the state vectors δ = δi , . . . , δn and ω = δ˙ = ω1 . . . , ωn , we may rewrite (11) in state-space form as " # 0n×n In δ 0 0n×1 δ˙ = + n×1 + (12) −M Lk −M D ω M pm Mu ω˙ where M = diag( m11 , . . . , m1n ), D = diag(d1 , . . . , dn ), Lk is the weighted Laplacian with edge weights kij , pm = m T T p1 , . . . , p m , u = ui , . . . , u n . n A. Centralized PI control Traditionally, automatic frequency control of power systems is carried out at two levels [2]. In the inner control loop the frequency is controlled with a proportional controller against a reference frequency. At the outer loop, the reference frequency is controlled with a PI controller to eliminate static errors. We model the inner controller of bus i as ui = α(ˆ ω − ωi )
and the outer PI controller is assumed to be given by 1X ωi zˆ˙ = ω ref − n i∈V 1X ω ˆ˙ = ω ref − ωi + β zˆ n
(13)
(14) (15)
i∈V
where we have assumed that the average frequency of the buses is measured by the central controller. In reality the
.. . ω ref
Figure 3.
un
Cn0
Bus n
ωn , δn
The figure illustrates the decentralized control architecture.
Theorem 5. The power system described by (12) where ui is given by (13)–(15), is stable for α > 0, β > 0, γ > 0 satisfying βγ < α. Furthermore limt→∞ ωi (t) = ω ref . Proof: As we are only interested in studying ω, we may rewrite (12) with u given by (13)–(15), as ω ˆ ω ˆ˙ − nγ 11×n 0 − nβ 11×n z · = 0n×n In z˙ 0n×1 ω ω˙ α1n×1 −M L −M D − αIn | {z } ,A
ref
γω + −ω ref 1n×1 m p
(16)
We now consider the the matrix A0 defined as 0 − nβ 11×n − nγ 11×n A0 , 0n×1 0n×n In α1n×1 −M L −αIn By performing elementary column operations on A0 , we may
write the characteristic polynomial of A0 as s3 γ + s2 α + sα + αβ · det M L + (s2 + sα) = 0 The first factor has solutions s ∈ C− if βγ < α. Comparing the second factor with the characteristic polynomial of M L, we have that s must satisfy s2 + sα + ti = 0 where ti ≥ 0 is an eigenvalue of M L. Also this equation has solutions s ∈ C− and one solution s = 0. However as A0 is full-rank, 0 cannot be an eigenvalue of A0 , so A0 is Hurwitz. But then clearly also A is Hurwitz. Finally, stationarity implies that the equilibrium solution must satisfy ω = ω ref 1n×1 . B. Decentralized PI control In this section we analyze a completely decentralized PI controller, where each bus controls its own frequency based only on local information. Thus, no frequency measurements need to be sent to a central controller, and there is no need to send control signals. This architecture might be favorable due to security concerns when sending unencrypted frequency measurements and control signals over large areas. Another benefit could be better performance when the tripping of one or several power lines causes the network to be split up into two or more sub-networks, so called islanding. The controller of node i is assumed to be given by z˙i = (ω ref − ωi )
ui = α(ω ref − ωi ) + βzi
(17) (18)
The controller architecture is illustrated in Figure 3. This decentralized control architecture is not practically feasible with only frequency measurements available at the generation buses. Even the slightest frequency measurement error will be integrated and prevent the frequencies from reaching consensus [2]. With recent advances in phasor measurement unit (PMU) technology however, phase measurements are likely to be available to all generator buses [18]. By employing optimal PMU placement, the number of PMU:s needed for complete observability can be drastically reduced [19]. By integrating (17) we obtain zi = ω ref t − δi
This implies that in order to accurately estimate the integral state zi , each generator bus needs access only to accurate time and phase measurements, both provided by PMU:s. The decentralized controller architecture is illustrated in figure 3. Theorem 6. The power system described by (12) where ui is given by (17)–(18), is stable for any choice of α > 0, β > 0. Furthermore limt→∞ ωi (t) = ω ref . Proof: As we are only interested in studying ω, we may rewrite (12) with ui is given by (17)–(18), as " ref # z˙ 0n×n In z −ω = · + (19) ω˙ −M L − βIn×n −M D ω M pm | {z } ,A
where the i’th element of β is βi . Let m = mini mi and d = mini di . We can now write: M D = m dIn + D0 where D0 is a diagonal matrix with non-negative entries. We now consider the the matrix A0 defined as 0n×n In 0 A , −M L − βIn×n −m · dIn By [17] the eigenvalues of A0 are given by det (s2 + sm d)In + M L + βIn×n
Comparing this with the characteristic polynomial of M L + diag(β), we conclude that s must satisfy s2 + sm d + ti = 0 where ti ≥ 0 is an eigenvalue of M L + diag(β). By [20] M L + diag(β) is positive definite, and hence the above equation has all its solutions in C− . It follows that also A is Hurwitz. Now consider the coordinate shift 0 z z z = − o ω0 ω ωo where z0 = (βIn + M L)−1 (M Dω ref + pm ), ω0 = ω ref 1n×1 . In the translated coordinates, 02n×1 is the only equilibrium of the system. Hence limt→∞ ωi (t) = ω ref ∀i ∈ V. C. Simulations The centralized and decentralized frequency control algorithms were tested on the IEEE 30 bus test system [21].The line admittances were extracted from [21] and the voltages were assumed to be 132 kV for all buses. M and D were set to reasonable numerical values. The power system is initially in an operational equilibrium, until the power load is increased by a step of 200 kW in the buses 2, 3 and 7. This will immediately result in decreased frequencies at the extra load buses. The frequency controllers at the buses will then control the frequencies towards the desired frequency of ω ref = 50 Hz. In the centralized controller the parameters were set to α = 0.8, β = 0, γ = 0.04, while in the decentralized control architecture the parameters were α = 0.8, β = 0.04. The step responses of the frequencies is plotted in figures 4–5. We note that if there is a centralized PI-controller for the reference frequency, the generation is increased uniformly among the generators. If however the integral action is distributed amongst the generators, some generators will increase their generation more than others. VI. D ISCUSSION AND C ONCLUSIONS In this paper we have studied consensus protocols with integral action for agents with double and single integrator dynamics. We have proved that the proposed consensus protocols reach asymptotic consensus even in the presence of constant disturbances. If we allow for absolute position measurements, the proposed consensus protocol converges asymptotically to a common state. In the absence of disturbances, the proposed consensus protocols asymptotically
φ(t) [Hz]
systems, with reasonable performance and possible benefits over centralized frequency control.
50
R EFERENCES
49.95 49.9
0
20
40
60
80
100
P (t) [kW]
t [s] 60 40 20 0
0
20
40
60
80
100
t [s]
φ(t) [Hz]
Figure 4. The upper figure shows the bus frequencies with centralized frequency control, while the lower figure shows the control signals at all buses. We observe that the frequencies converge quickly to a common frequency. The frequencies are then slowly regulated towards the reference frequency by the central PI controller.
50 49.95 49.9
0
20
40
60
80
100
60
80
100
P (t) [kW]
t [s] 60 40 20 0
0
20
40 t [s]
Figure 5. The upper figure shows the bus frequencies with decentralized frequency control, while the lower figure shows the control signals at all buses. We observe that the frequencies converge quickly to a common frequency. The frequencies are then slowly regulated towards the reference frequency by the PI controllers present at each bus.
solve the initial average consensus problem. Furthermore we have demonstrated that a similar consensus protocol can be employed for decentralized frequency control of power
[1] M. Morari and E. Zafiriou. Robust Process Control. Prentice Hall, Englewood Cliffs, 1989. [2] J. Machowski, J.W. Bialek, and J.R. Bumby. Power System Dynamics: Stability and Control. Wiley, 2008. [3] N. Senroy, G.T. Heydt, and V. Vittal. Decision tree assisted controlled islanding. Power Systems, IEEE Transactions on, 21(4):1790 –1797, nov. 2006. [4] Bo Yang, V. Vittal, and G.T Heydt. Slow-coherency-based controlled islanding: A demonstration of the approach on the august 14, 2003 blackout scenario. Power Systems, IEEE Transactions on, 21(4):1840– 1847, nov. 2006. [5] Florian Doerfler and Francesco Bullo. on the critical coupling for kuramoto oscillators. SIAM Journal on Applied Dynamical Systems, 10(3):1070–1099, 2011. [6] Iman Shames, Andr M.H. Teixeira, Henrik Sandberg, and Karl H. Johansson. Distributed fault detection for interconnected second-order systems. Automatica, 47(12):2757 – 2764, 2011. [7] R. Olfati-Saber and R.M. Murray. Consensus problems in networks of agents with switching topology and time-delays. IEEE Transactions on Automatic Control, 49(9):1520 – 1533, sept. 2004. [8] W. Ren and R.W. Beard. Distributed Consensus in Multi-vehicle Cooperative Control. Springer, 2008. [9] J.A. Fax and R.M. Murray. Information flow and cooperative control of vehicle formations. IEEE Transactions Automatic Control, 49(9):1465 – 1476, sept. 2004. [10] Liu, Y. and Passino, K.M. and Polycarpou, M.M. Stability analysis of w-dimensional asynchronous swarms with a fixed communication topology. IEEE Transactions Automatic Control, 48:76–95, 2003. [11] H.G. Tanner, A. Jadbabaie, and G.J. Pappas. Stable flocking of mobile agents, part i: Fixed topology. In IEEE Conference on Decision and Control, volume 2, pages 2010 – 2015 Vol.2, dec. 2003. [12] F. Bullo, J. Cort´es, and S. Mart´ınez. Distributed Control of Robotic Networks. Applied Mathematics Series. Princeton University Press, 2009. [13] L. Xiao, S. Boyd, and S. Lall. A scheme for robust distributed sensor fusion based on average consensus. In International Symposium on Information Processing in Sensor Networks, pages 63 – 70, april 2005. [14] R. Olfati-Saber. Distributed Kalman filter with embedded consensus filters. In IEEE Conference on Decision and Control, European Control Conference, pages 8179 – 8184, dec. 2005. [15] R. Carli, A. Chiuso, L. Schenato, and S. Zampieri. Distributed Kalman filtering based on consensus strategies. IEEE Journal on Selected Areas in Communications, 26(4):622 –633, may 2008. [16] R.A. Freeman, Peng Yang, and K.M. Lynch. Stability and convergence properties of dynamic average consensus estimators. In IEEE Conference on Decision and Control, pages 338 –343, dec. 2006. [17] John R. Silvester. Determinants of block matrices. The Mathematical Gazette, 84(501):pp. 460–467, 2000. [18] A.G. Phadke. Synchronized phasor measurements in power systems. Computer Applications in Power, IEEE, 6(2):10 –15, april 1993. [19] R.F. Nuqui and A.G. Phadke. Phasor measurement unit placement techniques for complete and incomplete observability. Power Delivery, IEEE Transactions on, 20(4):2381 – 2388, oct. 2005. [20] Andr´e Teixeira. Multi- agent systems with fault and security constraints. Master’s thesis, Royal Institute of Technology (KTH), June 2009. [21] Power systems test case archive-30 bus power flow test case. available in: http://www.ee.washington.edu/research/pstca/pf30/pgtca30bus.htm.
