Minimization of Frequency Deviations in Power Network by using majorant functions
Oleg O. Khamisov
math.OC










Optimization of Power Network Frequency Control via Analytic Bounds on the Objective Function
    
O. O. Khamisov

Center for Energy Systems,
 Skolkovo Institute of Science and Technology,
 Moscow, Russia


    
================================================================================================================



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

 



Frequency control in power networks is designed to maintain power balance by adjusting generation that allows to keep frequency at its nominal value (i.e. 50 Hz). If power disturbance occurs, it leads to frequency oscillations that are consequently suppressed by the control. Contemporary control schemes depend on a number of parameters. They can be adjusted to increase control efficiency. The main factor that defines control efficiency is maximal frequency deviation (nadir). However corresponding objective function (value of the nadir) is non-smooth and has multiple local extremums, thus it is difficult to optimize. The aim of this work is to present an analytic conservative estimate of the nadir and then optimize the estimate instead of optimizing nadir, since in practice global minimums of both functions coincide. Nadir and its estimate are non-smooth, therefore standard zero order method is applied. Numerical experiments show that the estimate has unique minimum while nadir has a set of local minimums. Therefore optimization of the estimate, unlike optimization of the nadir, not only computationally easy, but also returns desired global optimum. Such optimization allows to reduce nadir up to two times compared to the nadir with default control parameters. Moreover it allows to limit frequency oscillations in order to reduce wear of the equipment.




§ INTRODUCTION


Power networks are susceptible to power imbalances due to changes in power demand. Additionally generator failure may result in significant disturbance in power balance. As a result, electrical power drained from the generators exceeds mechanical power supply and generators start to slow down. Consequently frequency of the power network drops <cit.>, <cit.>, <cit.>. Large frequency oscillations may result in equipment damage or its emergency shutdown. 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 (50 Hz or 60 Hz).

It is known <cit.>, that maximal frequency oscillation (nadir) happens during first 30 seconds after power disturbance appearance. Currently used frequency control utilises proportional controller called Droop control (or Primary Frequency control) to reduce this initial frequency drop. Droop control 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 <cit.>-<cit.> 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 (nadir) without loss of system's stability.

Classical linearized model of transmission power network is considered. It is assumed that several buses suffer from a step change of power generation or consumption, which results in frequency oscillations. As a result dynamics of the network as defined by a system of linear differential equations

    ẋ=A(r)x+P^D.

Here vector x(t) includes subvector of bus frequency deviations ω(t). It is required to find vector of droop control parameters r=(r_1,…,r_n)^T, that would minimize nadir F(r)=max_l={1,n}max_t≥0|ω_l(t)| for all buses l∈ N, where N is the set of buses. Primary frequency control, considered in this paper stabilizes frequency at equilibrium during the first several tens of seconds <cit.>. Therefore in numerical experiments maximization is done for the first 100 seconds. Vector-function ω(t) is a part of the solution of linear differential equations, therefore it is given by oscillatory functions with infinite number of local extremum. Calculation of the nadir is done by d.c. approach which is computationally expensive. Nadir as a function of control parameters is non-smooth, therefore its optimization is done by Hooke and Jeevse zero order method. Finally nadir function F(r) has multiple local extremum, therefore it is not possible to guarantee that global optimum will be found by such approach.

Alternatively, analytical estimate for nadir is derived in order to improve optimization efficiency and exclude local extremums. Firstly we derive frequency majorants ℳ_l(t)≥ω_l(t) for each bus l that in practice have unique optimum, therefore their maximum G(r)=max_l={1,n}max_t≥0ℳ_l(t)≥ F(r) (estimation of the nadir) can be found via golden section method, which is computationally much faster, then d.c. optimization applied to calculate F(r). Function G(r) is non-smooth and is optimized by Hooke and Jeeves method. However in practice it has unique minimum which coincides with global minimum of F(r). As a result optimization of G(r) is computationally more efficient than optimization of F(r), moreover optimization of G(r) delivers global maximum for both G(r) and F(r), since G(r) does not have local extremums.

Numerical experiments are done in order to compare efficiency of estimate minimization with direct optimization of nadir.

The article is organized in the following way. In section <ref> network model and optimization aim are described.
In section <ref> estimate of the nadir is derived.
In section <ref> optimization of both the nadir and its estimate are described.
In section <ref> results of numerical experiments are presented.
In section <ref> final observations and directions of the future work are discussed.



§ PROBLEM STATEMENT



 §.§ Notations

Let i be imaginary unit.
Let ℝ be set of real numbers, cardinality of a finite set S is defined as (S).
For an arbitrary matrix (vector) X its transpose is denoted by X^T.
For an arbitrary vector x∈ℝ^n and set I⊆{1,…,n} we define subvector x_I.
For an arbitrary matrix X∈ℝ^n× and set I⊆{1,…,n} we define row submatrix X_I.
Let 0 be zero matrix of the corresponding size, I is identity matrix of the corresponding size, 0 is zero vector of the corresponding size, ρ is vector of ones of the corresponding size.
For vector x=(x_1,…,x_n) we denote by (x)=(x_1,…,x_n) diagonal matrix with elements x_l, l∈{1,…,n}.
Operations X, X, |X| are considered to be elementwise, if X is a vector or a matrix.



 §.§ Model description

Classical generator model <cit.>, <cit.> is used. The power transmission network <cit.> is described by a directed graph (N,E),  where N is the set of n buses defined by their indices: N={1,…,n}, E = {(l,j), l,j∈ N is set of m lines. It is assumed, that the network is connected.

Dynamics of the power transmission network is defined by the system of linear differential equations <cit.>-<cit.>. Kron reduction <cit.>-<cit.> is applied in order to exclude load buses, therefore system has the following form:


    m_lω̇_l=   -d_lω_l+∑_j:(j,l)∈ Ep_ji-∑_j:(l,j)∈ Ep_lj-
       -p^M_l+p^D_l, ω_l(0)=0, l∈ N,
    ṗ_lj=   b_lj(ω_l-ω_j), p_lj(0)=0, (l,j)∈ E,
    
    t^G_lṗ^M_l=   -p_l^M+v_l, p^M_l(0)=0, l∈ N,
    
    t^B_lψ̇_l=   -ψ_l-r_lω_l, ψ_l(0)=0, l∈ N.


Variables of the system have the following meanings:

  
  
  * ω_l, l∈ N are deviations of bus frequencies from nominal value,
  
  * p_lj, (l, j)∈ E are line active power flows,
  
  * p^M_l, l∈ N are mechanic power injections at generators,
  
  * ψ_l, l∈ N are positions of valves,
  

Parameters of the system:

    
  * m_l>0, l∈ N are generators inertia constants,
    
  * d_l>0, l∈ N are steam and mechanical damping of generators and frequency-dependent loads,
    
  * p^d_l, l∈ N are unknown disturbances (assumed constant),
    
  * b_lj>0, (l, j)∈ E are line parameters that depend on line susceptances, voltage magnitudes and reference phase angles,
    
  * t^G_l>0, l∈ N are time constants that characterize time delay in fluid dynamics in the turbine,
    
  * t^B_l>0, l∈ N are time constants that characterize time delay in governor response.
    
  * r_l>0, l∈ N are control parameters.


Equations of the system:

    
  * Equations (<ref>) are generator swing equations,
    
  * Equations (<ref>) are equations of direct current linearized power flows,
    
  * Equations (<ref>) are turbine dynamics,
    
  * Equations (<ref>) are governor dynamics.













Control parameters r are chosen non-negative 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.
















 §.§ Matrix Representation

In order to introduce nadir estimate the system (<ref>) is presented in the following form:

    ẋ=A(r)x+P^D, x(0)=0,



    A(r)=([   -MD   -MC     M     0;  BC^T     0     0     0;     0     0  -T^B   T^B; -T^VR     0     0  -T^V;       ]), x=([   ω(t);   p(t); p^M(t);   v(t) ]), P^D=([ Mp^D;    0;    0;    0 ]).


Here ω(t)=(ω_1(t),…,ω_n(t))^T, p(t)=(p_1(t),…,p_m(t))^T, p^M(t)=(p^M_1(t),…,p^M_n(t))^T, v(t)=(v_1(t),…,v_n(t), D= (d_1,…,d_n)≻0, M=(1/m_1,…,1/m_n)≻0, B= (1/b_1,…,1/b_m)≻0, T^G= (1/t^G_1,…,1/t^G_n)≻0, T^V= (1/t^V_1,…,1/t^V_n)≻0, R= (r_1,…,r_n)≻0, C∈ℝ^n× m is the incidence matrix <cit.> of the system graph G, p^D=(p^D_1,…,p^D_n)^T is the inhomogeneity vector of disturbances.



 §.§ Optimization aim

In the models, describing energy systems, matrix A is diagonalizable <cit.>-<cit.>. We are using this observation in the further results.
Let 

    A(r) = S(r)Λ(r)S^-1(r)

be its eigenvalues decomposition, where Λ(r) diagonal matrix of eigenvalues that has the following form

    Λ(r) =(
                         [ Λ^1(r)      0;      0      0;        ]),

where Λ^1 is diagonal matrix of nonzero eigenvalues.
Let us define nadir as a function of control parameters:

    F(r)=max_l∈ Nmax_t∈[0,t_1]|ω_l(t)|.

Nadir minimization problem  has the following form:


    min_r≥0F(r),

s.t.

    max_i∈ Nλ_i(r)≤0,


    max_l∈ N| λ_l(r)/ λ_l(r)|≥ξ, ξ>0.


where ω(t) is part of the solution of the system (<ref>), that depends on r. Primary frequency control, considered in this paper stabilizes frequency at equilibrium during the first several tens of seconds <cit.>. Therefore in the experiments t_1 is taken equal to 100 seconds. Constraint (<ref>) is used to ensure system's stability, constraint (<ref>) is used to suppress small frequency oscillations. Its effects are discussed in the numerical results section.

Functions ω_l(t) are oscillatory and have infinite amount of extremums as a solution of system of linear differential equations. In addition function max_t∈[0,t_1]|ω_l(t)| is non-smooth, therefore calculation of any value of F(r) is computationally difficult. In addition, as it is shown in numerical results section, function F(r) has local minimums. Thus approximating nadir with  function G(r)≥ F(r) of simpler structure is done. Then problem

    min_r≥0G(r),

is considered instead of (<ref>).




§ ESTIMATE OF THE NADIR

Firstly let us show that there always there always exists set of control parameters R for which frequency deviations ω(t) converge to a constant value.

Let R=0. System (<ref>) is asymptotically stable over ω.


If R=0 then matrix of the system (<ref>) as an upper triangular matrix:

    A(r)=([ 1-21|c-MD    1c|-MC         M         0;   1|cBC^T      1c|0         0         0;      1-30         0  1|c|-T^G       T^G;      3-40         0         0  1|c|-T^V;       4-4 ]).

Variables ψ do not depend on any other variables, and is defined by the equation

    ψ̇=-ψ,

thus is asymptotically stable. Similarly p^m is asymptotically stable as it depends only on ψ.
Diagonal block for the remaining variables ω and p has form

    ([   -MD    MC; -BC^T     0;       ])

and the corresponding homogeneous system is


    ω̇=    -MDω+MCp,
    ṗ =    BC^T.


Let us consider the following Lyapunov function:

    V(ω,p) = 
    ([ ω; p;   ])^T
    ([ M^-1    0;    0 B^-1;      ])([ ω; p;   ]).

Its derivative is given by

    V̇(ω,p) = -2ω^TDω.

That according Barbashin-Krasovskii-LaSalle theorem <cit.> proves asymptotic stability of ω.

Let us introduce the new eigenvalues matrix 
    Λ(r) =(
                         [ Λ^1(r)      0;      0     -I;        ]),

and new system matrix

    A(r) = S(r)Λ(r) S^-1(r).

Then the following lemmas can be derived:

Let

    x(t)=([ ω(t); z(t) ])

be solution of the system (<ref>) and

    x(t)=([ ω(t); z(t) ])

be solution of the system

    ẋ=A(r)x+P^D, x(0)=0.

then

    ω(t)=ω(t)



Solution of (<ref>) is given by

    x(t) = ∫_0^te^A(r)(t-τ)P^Ddτ = S(r)e^Λ(r) t∫_0^te^-Jτdτ S^-1(r)P^D =


    = S(r)e^Λ t(
               [ (Λ^1(r))^-1(I-e^-Λ^1(r)t)                         0;                         0                        It;                           ])S^-1(r)P^D = S(r)(
               [ (Λ^1(r))^-1(e^Λ^1(r)t-I)                        0;                        0                       It;                          ])S^-1(r)P^D.

Similarly

    x(t) = S(r)(
               [ (Λ^1(r))^-1(e^Λ^1(r)t-I)                        0;                        0                   I-e^-t;                          ])S^-1(r)P^D.

Let us consider  their difference

    x(t)-x(t) = S(r)(
               [         0         0;         0 It-I+e^-t;           ])S^-1(r)P^D =


    =(
               [ S^1_N(r)        0; S_Z^1(r) S_Z^0(r);          ])(
               [         0         0;         0 It-I+e^-t;           ])S^-1(r)P^D =


    =(
               [                   0                   0;                   0 S_Z^0(r)(It-I+e^-t);                     ])S^-1(r)P^D.

First n elements of this vector, corresponding to ω(t)-ω(t) are zero.



Vector function ω(t) has the following form:

    ω(t)=Y_N(r)e^Λ(r)ρ t-ω^*,

where

    Y=S(r)Λ(r) (S^-1(r)P^D)



From (<ref>)

    x(t) = ∫_0^te^A(r)(t-τ)P^Ddτ=S(r)Λ^-1(r)e^Λ(r)tS^-1(r)P^D-A^-1(r)p^D =


    = S(r)Λ^-1(r) (S^-1(r)p^D)e^Λ(r)ρ t-A^-1(r)P^D.

Since (r)lim_t→∞(SΛ^-1(r) (S^-1(r)P^D))_N e^Λ(r)ρ t=0 we have (A^-1(r)P^D)_N=ω^*. As is shown in the previous lemma, ω(t) does not depend on the bottom right block of Λ(r), therefore 
    (S(r)Λ^-1(r) (S^-1(r)P^D))_N=(S(r)Λ(r) (S^-1(r)P^D))_N.
 Combining this facts we have

    ω(t)=ω(t)=

    =(SΛ^-1(r) (S^-1(r)P^D))_N e^Λ(r)ρ t-(A^-1(r)P^D)_N=(S(r)Λ(r) (S^-1(r)P^D))_N e^Λ(r)ρ t-(A^-1(r)P^D)_N=

    =Y_N(r) e^Λ(r)ρ t-ω^*.






For the frequency deviation ω(t) the following estimation exists:


    |ω(t)|≤|Y_N(r)|e^ Λ(r) ρ t+|ω^*|



From

    |ω(t)-ω^*|=Y_N(r) e^Λ(r)ρ t+ω^*-ω^*=|Y_N(r) e^Λ(r) ρ t|≤|Y_N(r)|e^Λ(r) ρ t.

We have

    |ω(t)-ω^*|≤ M^1(t,r)=|Y_N(r)|e^Λ(r) ρ t

which gives theorem result.




    |ω(t)|≤|Y_N(r)|e^ Λ(r) ρ t+|ω^*|.




For the frequency deviation ω(t) in (<ref>) we have the following estimation:

    |ω_l(t)|≤ M_l^2(t,r)=∑_j=1^3n+m| y_lj(r)|e^λ_j(r)tmin{|λ_j(r)t|,1}+


    +∑_j=1^3n+m| y_lj(r)|min{|d/dtf_j(t^1_j,r)|t,|f_j(t_j^0,r)-1|},

where

    f_j(t,r)=e^λ_j(r)tcos(λ_j(r)t)-1,


    t_j^0={[ π/|λ_j(r)|+1/λ_j(r)arctanλ_j(r)/λ_j(r),                                λ_j(r)0,;                                       0                                         ].
 

    t^1_j={[                                                                0,                                                λ_j(r)0λ_j(r) = 0,; π/|λ_j(r)|+1/|λ_j(r)|arctan((λ_j(r))^2-(λ_j(r))^2/2λ_j(r)λ_j(r)),                                         (λ_j(r)0­λ_j(r))λ_j(r)≥0,;            1/|λ_j(r)|arctan((λ_j(r))^2-(λ_j(r))^2/2λ_j(r)λ_j(r)),                                         (λ_j(r)0­λ_j(r))λ_j(r)<0, ].




We can use the following notation

    ω(t)=ω^sin(t)+ω^cos(t),

where

    ω^sin(t)=Im(Y_N(r))e^Λ(r)tsin(Λ(r)t)ρ, ω^cos(t)=Re(Y_N(r))e^Λ(r)tcos(Λ(r)t)ρ+ω^*.

We will approximate each of this function separately.

    |ω^sin_l(t)|=|∑_j=1^3n+m y_lj(r)e^λ_j(r)tsin(λ_j(r)t)|≤∑_j=1^3n+m| y_lj(r)|e^λ_j(r)t|sin(λ_j(r)t)|≤


    ≤∑_j=1^3n+m|y_lj(r)|e^λ_j(r)tmin{|λ_j(r)t|,1}.


To approximate ω_c(t) we will use the following expression:

    |ω^cos_l(t)|=|∑_j=1^3n+m y_ljf_j(t,r)|,

Since ω(0)=0, and ω^sin(0)=0, we have ω^cos(0)=0.

Let t_j^0 be a solution of the problem

    max_τ>0|f_j(τ,r)-1|,

then

    t_j^0={[ π/|λ_j(r)|+1/λ_j(r)arctanλ_j(r)/λ_j(r),                                λ_j(r)0,;                                       0                                         ].
 
Let t^1_j be the solution of the problem

    max_τ≥0|f'_j(τ,r)|,

then

    t^1_j={[                                                                0,                                                λ_j(r)0λ_j(r) = 0,; π/|λ_j(r)|+1/|λ_j(r)|arctan((λ_j(r))^2-(λ_j(r))^2/2λ_j(r)λ_j(r)),                                         (λ_j(r)0­λ_j(r))λ_j(r)≥0,;            1/|λ_j(r)|arctan((λ_j(r))^2-(λ_j(r))^2/2λ_j(r)λ_j(r)),                                         (λ_j(r)0­λ_j(r))λ_j(r)<0, ].

Then we have the following estimation

    |ω_c(t)|=|∫_0^td/dt(∑_j=1^3n+m y_lj(r)f_j(η,r)dη)|≤∑_j=1^3n+m| y_lj(r)||∫_0^tf'_j(η,r)dη|=


    =∑_j=1^3n+m| y_lj(r)|min{∫_0^t|f'_j(η,r)|dη,|∫_0^tf'_j(η,r)dη|}=


    =∑_j=1^3n+m| y_lj(r)|min{∫_0^tmax_τ≥0|f'_j(τ,r)|dη,|f_j(t^0_j,r)-1|}≤


    ≤∑_j=1^3n+m| y_lj(r)|min{max_τ≥0|f'_j(τ,r)|t,|f_j(t^0_j,r)-1|}≤


    ≤∑_j=1^3n+m| y_lj(r)|min{|f'_j(t^1_j,r)|t,|f_j(t_j^0,r)-1|}.






    |ω(t)|≤ M(t,r)=min{M^1(t,r),M^2(t,r)},

where

    M^2(t,r)=(M_1^2(t,r),…,M_n^2(t,r))^T, M(t,r)=(M_1(t,r),…,M_n(t,r))^T,



Similarly to nadir F(r), estimate of the nadir is denoted as

    G(r)=max_l∈ Nmax_t∈[0,t_1]M_l(t,r).




§ MINIMIZATION OF THE NADIR



 §.§ Calculations of values F(r) and G(r)

In practice each functions M_l(t,r), l∈ N have unique maximum. Taking this observation into consideration, we assume, that G(r) can be found using golden section algorithm.

For nadir calculation d.c. approximation of frequency deviations is used. Functions ω_l(t) can be represented in the following way:

    ω_l(t)=h_l(t)-q_l(t), l∈ N,

where

    h_l(t) = ω_l(t)+1/2k_l(r)t^2, q_l(t)=1/2k_l(r)t^2.

Here k_l is obtained as follows.

    d^2/dt^2ω_l(t)=∑_j=1^3n+m
    e^ λ_j(r)t(( λ_j(r))^2 y_lj(r)-

    -2 λ_j(r) λ_j(r) y_lj(r)-( λ_j(r))^2 y_lj(r))cos λ_j(r)t+


    .(( λ_j(r))^2 y_lj(r)-2 λ_j(r) λ_j(r) y_lj(r)-( λ_j)^2 y_lj(r))sin λ_j(r)t)≤


    ≤∑_j=1^3n+m((( λ_j(r))^2 y_lj(r)-2 λ_j(r) λ_j(r) y_lj(r)-( λ_j(r))^2 y_lj(r))^2+.


    .+(( λ_j(r))^2 y_lj(r)-2 λ_j(r) λ_j(r) y_lj(r)-( λ_j(r))^2 y_lj(r)))^1/2=k_l(r).

Global maximum of each function |ω_l(t)|, l∈ N is obtained using branch and bound method with concave overestimators and d.c. approximation. Convergence of the method is given in <cit.>. Based on this optimization process values of F(r) are found.



 §.§ Optimization of F(r) and G(r)


Although it is known <cit.>, that eigenvalues and eigenvectors are continuous functions of matrix entries, they cannot be calculated analytically. Therefore Hooke and Jeeves algorithm <cit.> is used optimize both function F(r) and G(r) subject to the following constraints:

  
  * Parameters r must be non-negative.
  
  * Parameters r have to be chosen so, that system of differential equations will remain stable.
  
    max_l∈ Nλ_loverestimators=0.

  
  * System's matrix does not have purely imaginary eigenvalues.
  
  * Oscillation of the frequency oscillations must decrease at reasonable speed. There is no standardized constraints on the frequency oscillations decreases, therefore here the following constraint is used:
  
    max_l∈ N| λ_l(r)/ λ_l(r)|≥ξ, ξ>0.


If those constraints are violated, we take

    F(r)=G(r)=∞.

Points found by Hooke and Jeeves algorithm are denoted by r^F and r^G.



 §.§ Optimization algorithms

Optimization algorithm for F(r) has the following form:
 Minimization of F(r).


  * Currently applied (default) in the optimized power system vector of control parameters r^0 is taken as a starting point.

  * Hooke and Jeeves method is applied to optimize F(r). Each value of F(r) is calculated via d.c. optimization.


Algorithm for optimization G(r) has the same form with only difference in usage of golden section method for calculation of values G(r) instead of d.c. optimization. It is given below:
 Minimization of G(r).


  * Currently applied (default) in the optimized power system vector of control parameters r^0 is taken as a starting point.

  * Hooke and Jeeves method is applied to optimize F(r). Each value of F(r) is calculated via golden section method.




§ NUMERICAL RESULTS

The algorithms were coded in Matlab. Computations were made in PC with Intel Core i7 /2.4GHz / 16GB. Numerical results are presented in the table <ref>. Here column System shows for which power system experiments were held, γ is taken 0.01. Two systems are considered <cit.>: 3 generators system and system of New England, that consists of 10 generators. For each system 100 test with different vectors of disturbances. Averaged results are presented in the table. Column Optimization time represents time in seconds, required for the Hooke and Jeeves method, applied to functions F and G. Next column contains number of function (F or G) calculations required for the method. Last column contains values of function F (maximal absolute value of frequency deviations) at the starting point, after optimization of F and after optimization of G. Although in the last case function G is optimized, aim of the algorithm is to minimize maximal frequency deviations, therefore values of F are provided. Parameters in starting point r are taken from <cit.>.

Time, required for optimization of F is bigger than for G, due to the fact that, during every calculation of F algorithm has to solve n d.c. optimization problems, while during calculation of F golden section computations are required.
As can be seen from the table, optimization of G might give better set of parameters, than optimization of F. Moreover less calculations of the function are required. This effect can be seen in figures <ref> and <ref>. Here New England system is considered, all parameters r_l are frozen with the exception of the first two. Red points on figures <ref> and <ref> represent G(r^G) and F(r^G) respectively, green point on figure <ref> represents G(r^G). In numerical experiments F(r^G) is global optimum of F, so optimization of majorants gives better result, than optimization of maximal values of frequency deviations directly. This happens due to the fact, that function G is smoother. Consider subregion, containing both points G(r^G) and F(r^F) on figures <ref> and <ref>. Function F in this case have unique minimum F(r^F), while G has local minimums, and G(r^G) is one of them.

Dynamics of the frequencies and majorant for New England System are given in figures <ref>, <ref>, <ref> and <ref>. Figure <ref> represents behavior of the system with starting values of the parameters. Figure <ref> represents behavior of the system with starting values of the parameters. Figure <ref> represents behavior of the system with starting values of the parameters. Figure <ref> system dynamics, after optimization of F without suppression of oscillations (<ref>). As can be seen, here maximal frequency deviations are smaller, than in <ref>, however they do not decay, during the observation time.

















§ CONCLUSION


The paper proposes method, that allows to minimize maximum of absolute values of frequency deviations (nadir) in power network under droop control. Control parameters are obtained by this method in a way, that ensures system's stability and good decay rate of the oscillations.

Since frequency oscillations are described by oscillatory functions with infinite number of local extremums, it is necessary to apply d.c. optimization method in order to find nadir for any particular set of control parameters, which is computationally expensive. Moreover nadir as a function of control parameters is non-smooth and have several local extremums. As a result it is not possible to ensure obtaining of global optimum.

Within this paper analytic majorants of frequency oscillations are derived. In practice they have unique maximum, therefore estimate of the nadir based on this majorants can be obtained via golden section, which is faster, than d.c. approach. Additionally, as numerical experiments show, this estimate as a function of control parameters have unique minimum which coincides with global minimum of the nadir.

Efficiency of nadir minimization and estimate minimization is compared. Unlike in nadir minimization, optimization algorithm applied to the estimate cannot finish its work in local minimum and always finds global optimum. Additionally minimization of estimate is less computationally difficult due to exclusion of d.c. optimization.

The obtained approach allows to reduce nadir up to two times compared with the default values of control parameters, used in the power networks. Moreover it keeps frequency oscillations within acceptable limits, which prevents wear of the equipment.

At the future expansion of this technique is planned in order to optimizer secondary frequency control together with the primary frequency control.



§ APPENDIX


System

    A(r)x^*+P^D=0,

where

    x^*=([     ω^*;     p^*; (p^M)^*;     v^* ])

is consistent.


System (<ref>) can be represented as



    -MDω^*+MCp^*+M(p^M)^*+Mp^D=0,
       -BC^Tω^*=0,
       -(p^M)^*+v^*=0
       -Rω^*-v^*=0.


Here ∑_l=1^n(d_l+r_l)ω^*_l=∑_l=1^np^D_l, and ω^*_j=ω^*_l.
Now we need to show, that solution in p exists. It is given by the first set of equations (<ref>). From it we have

    MCp^*=-(D+R)ω^*-p^D.

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

    C^Ty=0

we have y^T((D+R)ω^*-p^D)=0. But C^Ty=0 has only solutions of the form y=kρ, k∈ℝ, therefore y^T(-(D+R)ω^*-p^D)=∑_l=1^n(-(d_l+r_l)ω^*_l-p^D_l)=0.



Every stationary solution

    x^*=([     ω^*;     p^*; (p^G)^*;     v^* ])

of the system (<ref>) satisfies

    ω^*_l=-∑_l=1^np^D_l/∑_l=1^n(d_l+r_l), l∈ N.



From (<ref>) ω^*_l=ω^*_j, sum of rows of (<ref>) gives

    0=∑_l=1^n((d_l+r_l)ω^*_l-p^D_l),

which gives corollary statement.



Kernel of A has the following form:

    (A)=(
            [ 0; p; 0; 0;   ]),

where p∈(C).

 Eigenvectors x from (A) have to satisfy the following system:


    -Dω+Cp+p^M=0,
       -C^Tω=0,
       -p^M+v=0,
       -Rω-v=0.


From (<ref>), (<ref>) p^M=v=-rω. Substituting the latter into (<ref>) we get

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

Sum of the rows of this equation gives

    ∑_l=1^n(d_l+r_l)ω_l=0.

But form (<ref>) ω_l=ω_j for all l∈ N, j∈ N, therefore ω=0 and Cp=0.


Matrix of the system (<ref>) does not have purely imaginary eigenvalues for r∈ U_, where U_ is some neighbourhood of 0.


Firstly let us consider case, when r=0. By contradiction let x be eigenvector and ik be purely imaginary eigenvalue. Then by definition Ax=ikx or in a subvector form


    -MDω+MCp+M(p^G)+Mp^D=ikω,
       -BC^Tω=ikp,
       -p^M+v=ikp^M,
       -v=ikv.


From (<ref>) v=0, from (<ref>) p^M=0 therefore (<ref>) and (<ref>) are equivalent to

    -Dω+i/kCC^Tω-ikω=0.

Separating real and imaginary parts of this equation we get

    (
        [          -D -1/kCC^T+kI;  1/kCC^T-kI          -D;             ])(
               [ ω; ω;   ])=(
                       [ 0; 0;   ]).

Let us consider determinant of the system's matrix by Schur complement formula we have

    (
        [          -D -1/kCC^T+kI;  1/kCC^T-kI          -D;             ])=(D)(D+(1/kCC^T+kI) D^-1(1/kCC^T+kI))0.

First determinant of the product is nonzero, since D is not singular by definition. Second determinant is nonzero as a sum of two positive definite matrices. As a result system (<ref>) does not have any solutions.

As is shown in lemma <ref>, kernel does not depend on the values of r. Therefore there exists some neighbourhood U_ of 0 such that for any r∈ U_ matrix A does not have purely imaginary eigenvalues.


Recall, that  A=S(r)Λ(r)S^-1(r) in the eigenvalues decomposition, where Λ(r) is diagonal matrix of eigenvalues with eigenvalues λ_1(r),…,λ_3n+m(r). Without loss of generality, it is assumed, that they are ordered with respect to real parts increase: λ_l(r)≤λ_l+1(r), l=1,3n+m-1.




    Λ(r) =(
               [ Λ^1(r)      0;      0      0;        ]),

where Λ^1(r) is diagonal nonsingular matrix, corresponding eigenvector matrix is given by

    S(r)=(
               [ S^1(r) S^0(r);        ])=(
               [ S^1_N(r)        0; S_Z^1(r) S_Z^0(r);          ]),

here S^1 is the submatrix of eigenvectors corresponding to nonzero eigenvalues, S^0 is the submatrix of vectors, corresponding to 0 eigenvalue.
Then the following lemma is correct.

Let

    x(t)=([ ω(t); z(t) ])

be solution of the system (<ref>) and

    x(t)=([ ω(t); z(t) ])

be solution of the system

    ẋ=A(r)x+P^D, x(0)=0.

then

    ω(t)=ω(t).



Solution of (<ref>) is given by

    x(t) = ∫_0^te^A(r)(t-τ)P^Ddτ = S(r)e^Λ(r) t∫_0^te^-Jτdτ S^-1(r)P^D =


    = S(r)e^Λ t(
               [ (Λ^1(r))^-1(I-e^-Λ^1(r)t)                         0;                         0                        It;                           ])S^-1(r)P^D = S(r)(
               [ (Λ^1(r))^-1(e^Λ^1(r)t-I)                        0;                        0                       It;                          ])S^-1(r)P^D.

Similarly

    x(t) = S(r)(
               [ (Λ^1(r))^-1(e^Λ^1(r)t-I)                        0;                        0                   I-e^-t;                          ])S^-1(r)P^D.

Let us consider  their difference

    x(t)-x(t) = S(r)(
               [         0         0;         0 It-I+e^-t;           ])S^-1(r)P^D =


    =(
               [ S^1_N(r)        0; S_Z^1(r) S_Z^0(r);          ])(
               [         0         0;         0 It-I+e^-t;           ])S^-1(r)P^D =


    =(
               [                   0                   0;                   0 S_Z^0(r)(It-I+e^-t);                     ])S^-1(r)P^D.

First n elements of this vector, corresponding to ω(t)-ω(t) are zero.


In system (<ref>) ω(t) is asymptotically stable, for r∈ U_st, where U_st is some neighbourhood of 0.



The solution of the (<ref>) can be presented in the following form:

    x(t)=S(r)Λ^-1(r) (S^-1(r)p^D)e^Λ(r)ρ t-A^-1(r)P^D,






    x(t) = ∫_0^te^A(r)(t-τ)P^Ddτ=S(r)Λ^-1(r)e^Λ(r)tS^-1(r)P^D-A^-1(r)p^D =


    = S(r)Λ^-1(r) (S^-1(r)p^D)e^Λ(r)ρ t-A^-1(r)P^D.



Vector function ω(t) has the following form:

    ω(t)=Y_N(r)e^Λ(r)ρ t-ω^*,

where

    Y=S(r)Λ(r) (S^-1(r)P^D)



Since (r)lim_t→∞(SΛ^-1(r) (S^-1(r)P^D))_N e^Λ(r)ρ t=0 we have (A^-1(r)P^D)_N=ω^*. As is shown in the previous lemma, ω(t) does not depend on the bottom right block of Λ(r), therefore (S(r)Λ^-1(r) (S^-1(r)P^D))_N=(S(r)Λ(r) (S^-1(r)P^D))_N. Combining this fact we have

    ω(t)=ω(t)=

    =(SΛ^-1(r) (S^-1(r)P^D))_N e^Λ(r)ρ t-(A^-1(r)P^D)_N=(S(r)Λ(r) (S^-1(r)P^D))_N e^Λ(r)ρ t-(A^-1(r)P^D)_N=

    =Y_N(r) e^Λ(r)ρ t-ω^*.



99
BV
A. R. Bergen, V. Vittal,
Power Systems Analysis, 3nd ed. Prentice
Hall, 2000.
MBB
J. Machowski, J. Bialek, and J. Bumby,
Power System Dynamics: Stability and Control, 3nd ed.
NJ, US: John Wiley & Sons, Inc.,
2008.
WW
A. J. Wood, B. F. Wollenberg,
Power Generation, Operation, and Control, 3nd ed.
NJ, US: John Wiley & Sons, Inc.,
1996.
BH
A. R. Bergen, D. J. Hill,
A Structure Preserving Model for Power System Stability Analysis,
IEEE Transactions on Power Apparatus and Systems,
Vol. PAS-100, no. 1,
1981, pp 25-35.
WWW
A. J. Wood and B. F. Wollenberg,
Power Generation, Operation, and Control, 3nd ed.
NJ, US: John Wiley & Sons, Inc., 1996.
KN
P. Kundur,
Power System Stability And Control, McGraw-Hill, 1994
ZML
C. Zhao, E. Mallada, S. Low
Distributed generator and load-side secondary frequency control in power networks,
2015 49th Annual Conference on Information Sciences and Systems (CISS).
2015, pp. 1-6.
MZL
E. Mallada, C. Zhao, S. Low,
Optimal load-side control for frequency regulation in smart grids,
arXiv preprint arXiv:1410.2931v3,
2015.
ZMLB
C. Zhao, E. Mallada, S. Low; J. Bialek,
A unified framework for frequency control and congestion management,
2016 Power Systems Computation Conference (PSCC),
2016, pp. 1-7.
BRB
R. B. Bapat,
Graphs and Matrices, Springer-Verlag London,
2014.
IEEE  Power Systems Test Case Archive; data available at https://www2.ee.washington.edu/research/pstca/
PST J. Chow and G. Rogers, Power System Toolbox; software available at http://www.eps.ee.kth.se/personal/vanfretti/pst/Power_System_Toolbox_Webpage/PST.html
K
G. Kron,
Tensor Analysis of Networks,
Wiley,
1939.
K2
G. Kron,
A set of principles to interconnect the solutions of physical systems,
Journal of Applied Physics,
Vol. 24,

DB
F. Dörfler, F. Bullo,
Kron Reduction of Graphs with Application to Electrical Networks,
IEEE Transactions on Circuits and Systems I: Regular Papers,
Vol. 60, no. 1,
2013, pp. 150-163.
Voevodin V. V. Voevodin, Yu. A. Kuznetsov,
Matrices and Computations,
Nauka, 1984 (in Russian).
Bellman R. Bellman,
Introduction to Matrix Analysis, 3d ed.,
University of Southern California, 1967.
Schur
R. A. Horn, F. Zhang,
Handbook of Linear Algebra, 3d ed.,
Chapman and Hall, 2013.
HT
R. Horst, H. Tuy
Global Optimization Deterministic Approaches,
3d ed.,
Springer,
1996.
BSS
M. S. Bazaraa, H. D. Sherali, C. M. Shetty,
Nonlinear Programming Theory and Algorithms, 3d ed.,
NJ, US: John Wiley & Sons, Inc.,
2006.
T
E. E. Tyrtyshnikov,
A Brief Introduction to Numerical Analysis,
Springer,
1997.
HC
Haddad  W.M.,  Chellaboina V.
Nonlinear Dynamical Systems and Control, A Lyapunov-Based Approach.
Princeton University Press,
2008, 948 p.

