Minimization of Frequency Deviations in Power Network by using majorant functions

# Minimization of Frequency Deviations in Power Network by using majorant functions

O. O. Khamisov
Center for Energy Systems,
Skolkovo Institute of Science and Technology,
Moscow, Russia
###### Abstract

Frequency control in power networks is designed to maintain power balance by adjusting generation, what allows to keep frequency at its nominal value (i.e. 50 Hz). If power disturbance occurs, it leads to frequency oscillations and deviation form nominal value, that are suppressed by the control. Behavior of the existing control depends on a number of parameters, that are currently chosen form the set, that guarantees control stability. But they can be adjusted within that set to increase control efficiency. The two main factors, that define efficiency are maximal frequency deviations and frequency convergence rate to equilibrium point. However frequencies functions on buses are highly oscillatory and have infinite amount of extremum. The aim of this work is firstly to present analytical conservative estimations of absolute values of frequencies deviations, in order to approximate dependence of frequency behavior on the control parameters and secondly present zero order optimization algorithm that would improve system dynamics by adjusting control parameters.

Keywords: power network; eigenvalues; linear differential equations; frequency control; d.c. optimization; Hooke and Jeeves method.

## 1 Introduction

Power networks are susceptible to power imbalances due to changes in power demand. Additionally generator or line failure may result in significant disturbance in power balance and power flows. As a result, frequency of the power network oscillates [1], [2], [3], and large oscillations may result in equipment damage or emergency shutdown or lines overload. In order to counter this effects frequency control is used. It adjusts power generation in order to restore power balance and deliver frequency to its nominal value (e.g. 50 Hz).

Currently used traditional version of frequency control consists of three different parts. Its first part called Droop control (or Primary Frequency control), which is aimed to reduce initial frequency drop after power disturbance appearance and works at timescale of tens of seconds. During this period the system is most vulnerable due to frequency oscillations. There exist various other versions of frequency control schemes [7]-[9] however they are not implemented in power systems, therefore they are not considered within this paper.

Droop control has a set of control parameters, which are currently chosen to ensure system’s stability. However it is possible to adjust these parameters to reduce frequency oscillations without loosing stability of the system. There is no agreed way to describe system’s behavior with a set of particular parameters without simulations. In particular maximal frequency deviations (nadir) from the nominal value are the key factors, that influence system’s reliability. Aim of the paper is development of approach, that would minimize maximal absolute value of frequency oscillations without loss of system’s stability.

Dynamics of power network are described by system of linear differential equations. Frequency deviations on each bus of the system is highly oscillatory and has infinite number of extremums, therefore calculation of maximum of their absolute values among all buses is computationally difficult.

Additionally, as a function of control parameters, this maximum is defined by eigenvectors and eigenvalues of the system’s matrix, which cannot be calculated analytically. As a result maximum of frequency deviations is a complicated function with many extremum. Applying zero order method for its minimization results in obtaining local minimum.

In order to counter this effects majorants for absolute values of frequency deviations are derived. These majorants still depend on eigenvalues and eigenvectors of the system’s matrix, but have simpler structure. As a result, zero order method often finds not local, but global optimum and does that in a much faster time, since there is no need in maximization of oscillatory function for every bus.

Classical linearized model of transmission power network is considered. It is assumed that several buses (both loads and generators) suffer from a step change of power generation or consumption, what results in frequency oscillations.Aim of this work is derivation of analytical estimations nadir and convergence rate of frequency to the nominal value.

Since control and system’s model are described by linear system of differential equations, eigenvalue analysis is applied to obtain necessary estimations.

The article is organized in the following way. In section 2 network model and optimization aims are described, in section 3 matrix representation of power network dynamics and description of vector variables are introduced, in section 4 axillary lemmas, necessary for derivation of majorants are given, in section 5 approach for search of absolute values of maximal frequency deviations on each bus is described, in section 3 majorants of absolute values of frequency deviations are introduced, in section 7 zero order optimization method is presented, in section 8 results of numerical experiments are presented, in section 9 final observations and directions of the future work are discussed.

## 2 Problem Statement

### 2.1 Notations

Let be set of real numbers, cardinality of a finite set is defined as . Imaginary one is defined as . For an arbitrary matrix (vector) its transpose is denoted by . For an arbitrary vector we define subvector . 0 is zero matrix of the corresponding size, is identity matrix of the corresponding size, is zero vector of the corresponding size, is vector of ones of the corresponding size. For vector we denote by diagonal matrix with elements . Operations , , are considered to be elementwise, if is a vector or a matrix.

### 2.2 Model description

Classical generator model [4], [5] is used. The power transmission network [1] is described by a directed graph , where is the set of buses, is set of lines. It is assumed, that the network is connected.

Dynamics of the power transmission network is defined by the system of linear differential equations [7]-[9]. Kron reduction [13]-[15] is applied, therefore system has the following form:

 ˙θl= ωl,l=¯¯¯¯¯¯¯¯1,n (1a) Ml˙ωl= −dlωl+∑j:(j,l)∈Epjl−∑j:(l,j)∈Eplj− −pMl+pdl,ωl(0)=0,l=¯¯¯¯¯¯¯¯1,n, (1b) plj= blj(θl−θj),plj(0)=0,(l,j)∈E, (1c) TGl˙pMl= −pMl+pCl,pMl(0)=0,l=¯¯¯¯¯¯¯¯1,n, (1d) TBl˙ψl= −ψl+pCl,ψl(0)=0,l=¯¯¯¯¯¯¯¯1,n. (1e)

Variables of the system have the following meanings:

• are deviations of bus voltage angles from their nominal values,

• are deviations of bus frequencies from nominal value,

• are deviations of line active power flows from their reference values,

• are mechanic power injections at generators,

• are positions of valves,

• are control values.

Parameters of the system:

• are generators inertia constants,

• are steam and mechanical damping of generators and frequency-dependent loads,

• are unknown disturbances (assumed constant),

• are line parameters that depend on line susceptances, voltage magnitudes and reference phase angles,

• are time constants that characterize time delay in fluid dynamics in the turbine,

• are time constants that characterize time delay in governor response.

Equations of the system:

• Equations (1b) are generator swing equations,

• Equations (1c) are equations of direct current linearized power flows,

• Equations (1d) are turbine dynamics,

• Equations (1e) are governor dynamics.

Equation (1c) can be differentiated. Variables can be substituted from (1a). As a result, system can be described purely by differential equations:

 Ml˙ωl= −dlωl+∑j:(j,l)∈Epjl−∑j:(l,j)∈Eplj− −pMl+pdl,ωl(0)=0,l=¯¯¯¯¯¯¯¯1,n, (2a) ˙plj= blj(ωl−ωj),plj(0)=0,(l,j)∈E, (2b) Tl˙pMl= −pMl+pCl,pMl(0)=0,l=¯¯¯¯¯¯¯¯1,n, (2c) TBl˙ψl= −ψl+pCl,ψl(0)=0,l=¯¯¯¯¯¯¯¯1,n. (2d)

### 2.3 Control description

Droop control is applied to counter frequency drop, that occurs in a case of power deficit. In practice it reduces maximal frequency deviations from the nominal value (e.g. 50Hz) of the system. Droop control is given buy the following formula

 pCl(t)=−rlωl(t),l=¯¯¯¯¯¯¯¯1,n, (3)

here

 ¯¯¯rl≥rl>0,l=¯¯¯¯¯¯¯¯1,n (4)

is control parameter. These parameters are chosen with respect to (4) in order to ensure system’s stability. It is possible to adjust this parameters in order to reduce maximal frequency deviations even further and still keep system stable.

### 2.4 Optimization aims

Problem statement has the following form:

 maxl=¯¯¯¯¯¯1,nmaxt∈[0,t1]|ωl(t)|→minr,¯¯¯r≥r>0,

where is part of the solution of the system (2), that depends on . Primary frequency control, considered in this paper stabilizes frequency at equilibrium during the first several tens of seconds [6]. Therefore in the experiments is taken equal to seconds.

Main purpose of the work is to minimize maximal frequency deviations among all buses by adjusting control parameters , . Firstly majorants for the absolute values of frequency deviations are derived based on the eigenvectors and eigenvalues of the system’s matrix. In practice they have unique maxima therefore golden section method can be applied to find them. Then Hooke and Jeeves method is used to minimize majorants, thus reducing maximal absolute value of frequency oscillations among all buses of the system.

## 3 Matrix Representation

The system (2) with control , substituted from (3) can be presented in the following form:

 ˙x=Ax+PD,x(0)=0, (5)
 A=ZQ (6)
 Z=⎛⎜ ⎜ ⎜ ⎜⎝M0000B0000TG0000TB⎞⎟ ⎟ ⎟ ⎟⎠,Q=⎛⎜ ⎜ ⎜ ⎜⎝−D−CI0CT00000−II−R00−I⎞⎟ ⎟ ⎟ ⎟⎠,
 x(t)=⎛⎜ ⎜ ⎜ ⎜⎝ω(t)p(t)pM(t)ψ(t)⎞⎟ ⎟ ⎟ ⎟⎠,PD(t)=⎛⎜ ⎜ ⎜⎝pD000⎞⎟ ⎟ ⎟⎠,

Here , , , , , , , , , , is the incidence matrix [10] of the system graph , is the inhomogeneity vector of disturbances.

## 4 Auxiliary lemmas

In the models, describing energy systems, matrix is diagonalizable [11]-[12]. Let be eigenvalues decomposition, where is diagonal matrix of eigenvalues with eigenvalues ordered with respect to real parts increase: . Further for each column (eigenvector or generalized eigenvector) of transition matrix we will use the following notation:

 s=⎛⎜ ⎜ ⎜ ⎜⎝ωspspMsψs⎞⎟ ⎟ ⎟ ⎟⎠.
###### Lemma 4.1.

Kernel of has the following form:

 ker(A)=⎛⎜ ⎜ ⎜⎝0ps00⎞⎟ ⎟ ⎟⎠,

where .

###### Proof.

Eigenvectors from have to satisfy the following system:

 −Dωs+Cps+pMs=0, (7a) −CTωs=0, (7b) −pMs+ψs=0, (7c) −Rωs−ψs=0. (7d)

From (1d), (1e) . Substituting the latter into (7a) we get

 −(D+R)ω+Cp=0.

Sum of the rows of this equation gives

 n∑l=1(dl+rl)ωl=0.

But form (7b) for all , , therefore and . ∎

###### Lemma 4.2.

There exists diagonal matrix such, that for all for which and for any corresponding matrix the following properties are satisfied:

1. Matrix is negative semi-definite.

2. Matrix does not have purely imaginary eigenvalues.

###### Proof.

(1) Consider the following matrices:

 G=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝−I−D−12CI0CTD−1200000−II−RD−1200−I⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,P=⎛⎜ ⎜ ⎜ ⎜⎝D−120000I0000I0000I⎞⎟ ⎟ ⎟ ⎟⎠.

Then . Let us consider quadratic for with :

 xTGx=−ωTω+ωTpM−(pM)TpM+ψTpM−ψTψ−ωTRD−1ω=
 =−12(ω−pM)T(ω−pM)−12(pM−ψ)T(pM−ψ)−12ωTω−12ψTψ−ωTRD−1ψ.

Therefore if , and, as a result [16], matrices and are negative semi-definite. We will take

 ¯¯¯¯R1=D.

(2) From the contradiction. Let us assume, that is a purely imaginary eigenvalue. Than we have the following system:

 −MDωs+MCps+MpMs=ikωs, (8a) −BCTωs=ikps, (8b) −TGpMs+TGψs=ikpMs, (8c) −TBRωs−TBψs=ikψs. (8d)

Excluding all variables except we have

 Hωs=0, (9)

where

 H=−MD+ikMCBCT−
 −MR(I−k2(TG)−1(TB)−1−ik((TG)−1+(TB)−1))((I+k2(TG)−2)(I+k2(TB)−2))−
 −ikI.

Here

 ReH=−MD−MR(I−k2(TG)−1(TB)−1)((I+k2(TG)−2)(I+k2(TB)−2)).

Matrix

 V(k)=(I−k2(TG)−1(TB)−1)((I+k2(TG)−2)(I+k2(TB)−2))

is diagonal with diagonal elements equal

 Vl(k)=1−k2TGlTBl(1+k2(TGl)2)(1+k2(TBl)2),l=¯¯¯¯¯¯¯¯1,n.

It can be seen, that

 limk→−∞Vl(k)=limk→∞Vl(k)=0,l=¯¯¯¯¯¯¯¯1,n.

Functions are continuous, therefore they have infima . Let us choose such that

 −MD−M¯¯¯¯R2V––≻0.

Than for all element wise we have

 ReH≻0.

Complex system of linear equations (9) is equivalent to

 (ReH−ImHImHReH)(ReωsImωs)=(00). (10)

Real and imaginary parts of are symmetrical matrices, and have real eigenvalues [17]. is diagonal nonsingular matrix. Determinant of the matrix from (10) can be found using Schur Complement [18].

 det(ReH−ImHImHReH)=detReHdet(ReH+ImH(ReH)−1ImH)=
 =detReHdet(I+ImH(ReH)−1ImH(ReH)−1)detReH=
 =detReHdet(I+(ImH(ReH)−1)2)detReH.

Here , therefore , , hence and . Matrix of the system (10) is not singular and . As a result .

Taking

 ¯¯¯¯R=min{¯¯¯¯R1,¯¯¯¯R2}

we ensure that statements of the lemma will be fulfilled. ∎

###### Lemma 4.3.

System

 Ax+pD=0 (11)

is consistent.

###### Proof.

System (11) can be represented as

 −MDω+MCp+M(pG)+MpD=0, (12a) −BCTω=0, (12b) −(pG)+ψ=0 (12c) −Rω−ψ=0. (12d)

Similarly to lemma 4.1 we have , and . Now we need to show, that solution in exists. It is given by the first set of equations (12a). From it we have

 MCp=−(D+R)ω−pD. (13)

Here sum of the elements of the left hand side vector is equal . By the Fredholm theorem of alternative (13) has solution if for any such that

 CTy=0 (14)

we have . But has only solutions of the form , , therefore . ∎

###### Corollary 4.4.

Every stationary solution

 x∗=⎛⎜ ⎜ ⎜ ⎜⎝ω∗p∗(pG)∗ψ∗⎞⎟ ⎟ ⎟ ⎟⎠

of the system (1) satisfies

 ω∗l=∑nl=1pDl∑nl=1(dl+rl),l=¯¯¯¯¯¯¯¯1,n. (15)

If system’s graph is not a tree, then is singular and eigenvalues matrix can be presented in the following form

 Λ=(Λ1000), (16)

where is diagonal nonsingular matrix, corresponding eigenvector matrix is given by

 S=(S1S2)=⎛⎜ ⎜ ⎜ ⎜⎝ccΩs10Ps1Ps2PGs10Ψs10⎞⎟ ⎟ ⎟ ⎟⎠,

here is the submatrix of eigenvectors corresponding to nonzero eigenvalues, is the submatrix of vectors, corresponding to eigenvalue. Let us introduce the following matrix

 ¯¯¯¯Λ=(Λ100Λ2),

where is an arbitrary matrix and matrix

 ¯¯¯¯A=S¯¯¯¯ΛS−1.

Then the following lemma is correct.

###### Lemma 4.5.

Let

 x(t)=⎛⎜ ⎜ ⎜ ⎜⎝ω(t)p(t)pG(t)α(t)⎞⎟ ⎟ ⎟ ⎟⎠

be solution of the system (5) and

 ¯¯¯x(t)=⎛⎜ ⎜ ⎜ ⎜⎝¯¯¯ω(t)¯¯¯p(t)¯¯¯pG(t)¯¯¯¯α(t)⎞⎟ ⎟ ⎟ ⎟⎠

be solution of the system

 ˙¯¯¯x=¯¯¯¯Ax+pD,¯¯¯x(0)=0. (17)

then

 ω(t)=¯¯¯ω(t).
###### Proof.

Solution of (5) is given by

 x(t)=∫t0eA(t−τ)pDdτ=SeΛt∫t0e−ΛτdτS−1pD=
 =SeΛt(Λ−11(I−e−Λ1t)00It)S−1pD=S(Λ−11(eΛ1t−I)00It)S−1pD.

Similarly

 ¯¯¯x(t)=S(Λ−11(eΛ1t−I)00∫t0eΛ2(t−τ)dτ)S−1pD.

Let us consider their difference

 x(t)−¯¯¯x(t)=S(000∫t0eΛ2(t−τ)dτ)S−1pD=
 =(S1S2)(000∫t0eΛ2(t−τ)dτ)S−1pD=
 =⎛⎜ ⎜ ⎜ ⎜⎝Ωs10Ps1Ps2PGs10Ψs0⎞⎟ ⎟ ⎟ ⎟⎠(000∫t0eΛ2(t−τ)dτ)S−1pD=
 =⎛⎜ ⎜ ⎜ ⎜⎝000Ps2∫t0eΛ2(t−τ)dτ0000⎞⎟ ⎟ ⎟ ⎟⎠S−1pD.

First elements of this vector, corresponding to are zero. ∎

Further it is assumed, that is negative definite diagonal matrix.

###### Corollary 4.6.

Solution of (1) always converges to some stationary pointy.

###### Proof.

It is a consequence of the fact, that solution of (2) is the same as solution of asymptotically stable system (17). ∎

###### Lemma 4.7.

The solution of the (1) can be presented in the following form:

 ¯¯¯x(t)=S¯¯¯¯Λ−1diag(S−1pD)e¯¯¯Λρt+x∗,

where

 x∗=⎛⎜ ⎜ ⎜ ⎜⎝ω∗p∗(pM)∗ψ∗⎞⎟ ⎟ ⎟ ⎟⎠=−¯¯¯¯A−1pD.
###### Proof.
 ¯¯¯x(t)=∫t0e¯¯¯¯A(t−τ)pDdτ=S¯¯¯¯Λ−1e¯¯¯ΛtS−1pD−¯¯¯¯A−1pD=
 =S¯¯¯¯Λ−1diag(S−1pD)e¯¯¯Λρt−¯¯¯¯A−1pD.

###### Corollary 4.8.

For diagonal matrix with diagonal entries

 ^λl={1λlλl≠00λl=0

and

 ^x(t)=YeΛρt−¯¯¯¯A−1pD,

where

 Y=S^Λdiag(S−1pD),

we have

 ^ω(t)=¯¯¯ω(t)=ω(t).
###### Corollary 4.9.

For the following equality is correct

 ω(t)=YωeΛt+ω∗=ReYωeReΛtcosImΛt−ImYωeReΛtsinImΛt+ω∗ (18)

or in element-wise representation

 ωl(t)=3n+m∑j=1eReλjt(ReyljcosImλlt−ImyljsinImλlt)+ω∗l. (19)

## 5 D.c. approximation of frequency deviations

Functions can be represented in the following way:

 ωl(t)=hl(t)−ql(t),l=¯¯¯¯¯¯¯¯1,n,

where

 hl(t)=ωl(t)+12klt2,ql(t)=12klt2.

Here is obtained as follows.

 d2dt2ωl(t)=3n+m∑j=1eReλjt(((Reλj)2Reylj−2ReλjImλjImylj−(Imλj)2Reylj)cosImλjt+
 ((Imλj)2Imylj−2ReλjImλjReylj−(Reλj)2Imylj)sinReλjt)≤
 ≤3n+m∑j=1(((Reλj)2Reylj−2ReλjImλjImylj−(Imλj)2Reylj)2+
 +((Imλj)2Imylj−2ReλjImλjReylj−(Reλj)2Imylj))12=kl.

Global maximum of each function is obtained using branch and bound method with concave overestimators and d.c. approximation. Convergence of the method is given in [19].

We denote obtained maximum of absolute value of frequency deviation as

 ωmaxl=maxt∈[0,t1]|ωl(t)|,l=¯¯¯¯¯¯¯¯1,n.

## 6 Majorants of frequency deviations

###### Theorem 6.1.

For the frequency deviation in (5) the following estimation exists:

 |ω(t)−ω∗|≤M1(t)=|Yω|eReΛρt.
###### Proof.
 |ω(t)−ω∗|=|YωeΛρt+ω∗−ω∗|1:n=|YeΛρt|≤|Yω|eReΛρt.

###### Corollary 6.2.
 |ω(t)|≤|Yω|eReΛρt+|ω∗|
###### Theorem 6.3.

For the frequency deviation in (5) we have the following estimation:

 |ωl(t)|≤M2l(t)=3n+m∑j=1|Imylj|eReλjtmin{|Imλjt|,1}+
 +3n+m∑j=1|Reylj|min{∣∣∣ddtfj(t1j)∣∣∣t,|fj(t0j)−1|},

where

 fj(t)=eReλjtcos(Imλjt),
 t0j=⎧⎪⎨⎪⎩(π+arctanReλjImλj)/Imλj,λj≠0,0otherwise.
 t1j=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩arctan((Reλj)2−(Imλj)22ReλjImλj)Reλj≠0 and Imλj≠0,ReλjReλj≠0 and Imλ=0,0otherwise.
###### Proof.

Since , we can use form

 ω(t)=ωs(t)+ωc(t),

where

 ωs(t)=Im(Yω)eReΛtsin(ImΛt)ρ,ωc(t)=Re(Yω)eReΛtcos(ImΛt)ρ+ω∗.

We will approximate each of this function separately.

 |(ωs(t))l|=∣∣ ∣∣3n+m∑j=1ImyljeReλjtsin(Imλjt)∣∣ ∣∣≤3n+m∑j=1|Imylj|eReλjt|sin(Imλjt)|≤
 ≤3n+m∑j=1|ylj|eReλjtmin{|Imλjt|,1}.

To approximate we will use the following expression:

 |(ωc(t))l|=∣∣ ∣∣3n+m∑j=1Reyljfj(t)∣∣ ∣∣,

Since , and , we have .

Let be a solution of the problem

 maxτ>0|fj(τ)−1|,

then

 t0j=⎧⎪⎨⎪⎩(π+arctanReλjImλj)/Imλj,λj≠0.0otherwise.

Let be the solution of the problem

 maxτ≥0∣∣f′j(τ)∣∣,
 t1j=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩arctan((Reλj)2−(Imλj)22ReλjImλj)Reλj≠0 and Imλj≠0,ReλjReλj≠0 and Imλ=0,0otherwise,

Then we have the following estimation

 |ωc(t)|=∣∣ ∣∣∫t<