A stability result for network reduced power systems including virtual friction control

We prove that virtual friction control can stabilize a power grid containing several virtual synchronous machines (VSMs), connecting line impedances and loads. Virtual friction is a torque added to the swing equation of each VSM, proportional to the deviation of its frequency from the overall center of inertia (COI) frequency. Our analysis is based on the network reduced power system (NRPS) model. We support our results with simulations for a two-area network of four VSMs, looking at the transients induced by a change of tie-line impedance and an asymmetric load change. We compare the results for the NRPS model with the corresponding results using detailed models of synchronverters and line impedances. We find that virtual friction has a strong stabilizing effect both for the NRPS model and for the full grid model.


I. INTRODUCTION
Virtual synchronous machine (VSM) are switched power converters that are controlled such that they mimic (towards the power grid) the behaviour of synchronous generators, providing frequency and voltage droop as well as (virtual) inertia. They seem to be a very promising way for the development of inverters for the grid integration of distributed power generators, see for instance [3], [6], [16], [18], [19], [20], [25], [32], [39], [40]. VSMs are more flexible than synchronous generators (SGs), as their parameters can be tuned and even changed on-line. Moreover, the frequency and voltage droops can act instantly, while in SGs connected to prime movers, they have large time constants. The main advantage of VSMs, as compared to current source inverters common in renewable energy sources, is that they can form stable grids, just like SGs. Their potential of providing grid frequency support by means of virtual inertia and frequency droop is shown, for instance, in [2], [14], [20], [21] and [30]. Of course, the stability of a (micro)grid comprising VSMs, loads and transmission lines is not guaranteed, a careful analysis and parameter tuning are needed to achieve a stable grid. Thus, the stability analysis of power grids is an important and timely topic, see for instance [1], [4], [8], [9], [11], [13], [15], [17], [25], [28], [31], [36].
Investigating a single inverter connected to an infinite bus, [29] and [30] have proven global asymptotic stability based on hybrid angle control, a combination of DC-matching control and a non-linear angular damping feedback, based on angular differences evolving on a circle with radius 4π. This type of control allows the phase differences of connected VSMs to reach and exceed 180 • . The paper [26] investigates the stability of grids where several VSMs are operated in parallel, by means of µ-analysis. The paper [5] investigates the stability of a power network employing VSMs based on the singular perturbation method, separating the system dynamics into slow, fast and very fast dynamics. It derives sufficient stability conditions for the inverter parameters to guarantee stability of such a network.
A major issue with VSMs employing frequency droop control is that they have to inject large excess power during disturbances where the grid frequency drops below its nominal value. Hardware power ratings of the electronic components and/or the limited availability of excess DC power impose an upper bound for the frequency droop constant and/or they force the droop torque to saturate, which is bad for stability. To avoid reaching saturation frequently, the frequency droop constant (per unit) employed in inverters is typically much smaller than that imposed on the prime mover of an SG [19].
In order to decouple stability requirements from hardware requirements, we propose to introduce virtual friction (VF), as an additional viscous friction torque, acting on the virtual rotors of the connected VSMs. This torque employs the center of inertia (COI) frequency as a reference. (The concept of COI has been introduced in [27] and has been used extensively, see for instance [38]). VF has been introduced in [7] and [35] but only for two VSMs (or two areas), so that the friction torque is proportional to the difference in frequency between these two, without the need for the COI frequency. Simulation experiments in [7] show that VF greatly enhances the stability of a two-area power system with weak links between the two areas. In damping oscillations, the effect of VF is similar to that of the frequency droop, but without the need to inject large excess power if the grid frequency is below the nominal value. The paper [12] has extended this concept and implemented a PID-controller based on VF.
It is important to state that the calculation of the COI frequency requires the transmission of angular velocity information of the VSMs across the power network. The transmission time delays that occur in the communication infrastructure complicate the stability analysis of such systems. In particular, time delays may not be constant over time and may have different values for each communication link, see for example [23] for a stochastic modelling of time delays. In a real grid, time delays may be due to (see [37]): • Transmission delay: Time a packet needs for trans-mission via the employed protocol, dependent on the bandwidth. • Propagation delay: Time for the data to propagate through the medium (copper wire or fiber optic cable). • Processing delay: Time required by the data acquisition unit to prepare and send a frame + time required by the VSM to receive and decode the received data. In order to bypass time delays, a decentralized approach was proposed in [9], where an additional control loop is implemented in the VSM algorithm, based on angular acceleration and active power output. A small-signal approach is used to tune the control parameters. In order to understand the impact of time delay on the VF mechanism, we have conducted exhaustive simulations. Our surprising conclusion is that the system behaviour remains largely unchanged for constant but not necessarily equal time delays up to around 100ms, at least for the simple two area grid with 4 VSMs studied in this paper.
The stability analysis of larger power grids must of course involve simplified models of the generators, transmission lines and loads, because working with full models (high order non-linear systems) makes any analysis intractable, see for instance [24], [33]. An important and successful direction of research in this field has been the study of the network reduced power system (NRPS) model (see [11], [15] and [27]) and its inertia-less limit version, the nonlinear Kuramoto (NK) model. In the NRPS model, each generator is represented by a system of order two (the swing equation) and the transmission lines and loads are static and are represented by complex impedances without state variables. Stability results for the NK and NRPS models have been derived in [15] and improved in [34]. They give -under reasonable assumptions on the impedances of the power network and for small enough inertias -a local exponential stability result on the grounded rotor angles of the synchronous generators in the network.
In this paper we extend the NRPS model by incorporating VF. In [15] and in the stronger version of the main result that is in [34], it was proven that, under suitable assumptions, the NRPS model becomes stable if the inertias of the generators are sufficiently small. However, making the inertias very small comes at a price: The system frequency will vary a lot more and much faster. The advantage of VF is that we no longer have to make the inertia terms small, instead, to achieve stability, it suffices to increase the VF coefficients.
We test VF in simulations of a micro grid, using full models for every component, and we also test the corresponding reduced model which is the NRPS model. For each of the two models (the full and the reduced one) we simulate a change of tie-line impedance between two areas, as well as an asymmetric change of load, and we observe the resulting transients with and without VF. Our results show that indeed, the addition of VF stabilizes the system even if the network is chosen such that it is unstable without the additional damping term. Moreover, with VF, the transients decay much faster. This paper is structured as follows: In Section II we recall the NK and NRPS models and we introduce the new model that has VF. We also state the stability result for the model with VF. Section III contains the proof for the main theorem of this paper, in which we consider separately the dynamics of the COI frequency and then we show that, on a different time scale, the new model with VF reduces to the NRPS model. This allows us to perform a singular perturbation analysis, as in [34]. In Section IV we report our simulations, that show the influence of the VF torques on the second order SG models and compare the evolution of the grounded rotor angles for the same scenarios of the grid. In addition, a more realistic simulation based on the synchronverter detailed in [22] is performed on the same grid, also employing the full model for the transmission lines and loads.

II. POWER FLOW MODEL WITH VIRTUAL FRICTION
Notation. We follow the notation from [34] (see also [15]). We represent all angles on the unit circle T by numbers in (−π, π] (modulo 2π). Addition and subtraction are performed modulo 2π. For two angles θ 1 , θ 2 ∈ T the geodesic distance between them is Let n ∈ N be fixed. For any γ ∈ (0, π] we define The above equation means that there exists an open arc of the length γ in T that contains all angles θ j , and θ 0 is the midpoint of this arc. For γ ∈ (0, π] we define ∆(γ) as the closure of ∆(γ) in T n . As in [15], we introduce the mapping grnd : T n → R n−1 as δ = grnd(θ ) where: The above representation of angle differences loses one degree of freedom, since for n angles n − 1 differences are defined. Where required in the following, δ n = 0. For any γ ∈ (0, π] the grounded set ∆ grnd (γ) ⊂ R n−1 is defined as: We denote the closure of this grounded set by ∆ grnd (γ).
Parameters of the models. The matrix [a jk ] ∈ R n×n is such that a jk = a k j > 0 for 1 ≤ j, k ≤ n, j = k and a j j = 0 for j ∈ {1, ... n}. We define an array of phase shifts ϕ jk ∈ T that satisfies The diagonal values ϕ j j are not relevant. We also introduce four sets of n real numbers: M j > 0, D j > 0, F j > 0 and p j ∈ R, j ∈ {1, 2, . . . n}.
The NK model is a collection of n first order differential equations on T n , defined for 1 ≤ j ≤ n. It has been proposed in [15] as an extension to the Kuramoto model, describing 4229 synchronization between a set of coupled oscillators connected through a lossy network: The state of this system is θ = (θ 1 , ...θ n ), which evolves in the state space T n . We denote ω j =θ j . As a simplified model describing a power network, the NRPS model has been derived in [27] (see also Chapter 6 in [10] and [34]). In this model, n SGs are connected via a passive network and θ 1 , . . . θ n are the rotor angles of the generators. The model consists of n second order differential equations on T: for 1 ≤ j ≤ n, We will refer to M j as inertia terms, to D j as damping coefficients and to p j as power terms. The precise connection between these parameters and the parameters of an electric grid is derived in Section 3 of [34]. Let us denote ω j =θ j . The state of the dynamical system (2) is: where, again, ω j =θ j . Let ω c be the center of inertia frequency, defined as in [27]: In this paper we are interested in a modification of the NRPS model, where a new viscous friction term F j (ω j − ω c ) is added to the left-hand side of (2). This leads to: We will refer to F j as a virtual friction coefficient. We assume that the ratios of damping and virtual fricion coefficients over inertia terms are the same for all the generators in the network. This enables us to define the parameters d and f as: Note that the actual values of d and f would be defined beforehand and individual inverters would be configured to fulfil this requirement prior to their integration into the grid.
With the above, we define the following model that we shall call the friction enhanced power system (FEPS) model for convenience: The state of the system (7) is again defined as in (3).
In the following, we show the improved stability properties of this model over the NRPS model. For this, we impose the exact same assumptions on the parameters D j , p j , a jk and ϕ jk as in the main result of [34].

III. THE PROOF OF THEOREM 2.1
Step 1. The COI frequency is a function of the frequencies ω j . In order to analyse the system stability, we have to investigate the stability of ω c . By differentiating (4) and using (7), we obtain that where P e, j (θ ) is defined by 4230 and its interpretation is the electrical power output of generator j. By the definition of ω c , n ∑ j=1 M j (ω c − ω j ) = 0. Thus (11) can be written as: Dividing by n ∑ k=1 M k gives a first order linear differential which shows that ω c converges to some limit ω * c if the functions P e, j (θ ) converge to some limits.
Step 2. The FEPS model (7) corresponds to the NRPS model with the additional disturbance term f ω c and an enhanced damping coefficient M j ( f + d) (in place of M j d).
Lemma 3.1: Let (θ , ω) be a solution of (7). Let (θ 0 , ω 0 ) be the solution of the related equation where P e, j (θ 0 ) is defined as in (12) and ω 0 j =θ 0 j , starting from the same initial conditions: Then the function ζ is independent of j and it is the unique solution oḟ starting from ζ (0) = 0, where ω c is defined in (4).
Notice that the only real difference between the equations (7) and (13) is thatp j has been replaced with p j .
Step 3. We introduce new variables x j ( j ∈ {1, 2, ...n}) and a new timescale τ as follows: Substituting for t and ω 0 j in (13) gives: We choose β = d + f so that ε = 1 β 2 by (8). This allows rewriting equation (13) in the form of a standard perturbation problem with ε as the singular perturbation parameter: For ε → 0, the reduced system dynamics are given by: Writing the above using the grounded angles from equation (1) gives where δ 0 j = θ 0 j − θ 0 n and δ 0 n = 0. The equations (18) are now of the same form as the equations (19), (20) from [34], with x j in place of ω j , and with M j in place of m j and also in place of D j from [34]. Thus we can apply Theorem 4.1 from [34] to conclude that there exist suitable conditions (as stated in Theorem 2.1) for which the following holds: It follows thus with (16) and the fact that θ 1 = θ , that P e, j (δ ) =P e, j (δ 0 ), whencẽ P e, j (δ ) →P e, j (δ 0 * ).

IV. SIMULATION RESULTS
This section presents simulation results for a two area network where each area consists of two generators and one load (See Fig. 1). An additional load L 3 in parallel with load L 2 serves to introduce a load step change in area 2. The line impedance between the two areas can be changed by operating switch S 1 , (dis-)connecting transmission line Y 5 . We present the simulations based on the full VSM algorithm from [18], [22] as well as results for the NRPS and FEPS model. All simulations were conducted in the MATLAB/Simulink (R) environment.

A. The full model
The VSM algorithm was run at 10kHz while the power part is simulated at a step size corresponding to 100kHz. We use the ode4 (Runge-Kutta) solver for these simulations. Loads are modelled as constant impedance and transmission lines as RL circuits. Switches in the model have negligible breaker resistance and snubber resistance of 1kΩ and snubber capacitance of 200nF. All synchronverters have nominal power 9kW and have identical parameters with M j = 188.5Ws 2 and D j = 94.2Ws and F j = 3141.6Ws. The set powers of all inverters were P set = 5kW and Q set = 0kVar.
A 1ms round trip time delay for the COI frequency to the synchronous machines is used. The transmission network admittances Y 1 to Y 5 consist of R = 0.1Ω and L =1mH and Y 6 is equivalent to R = 10Ω and L =100mH. The inverter output are connected to symmetric LCL-filters Y F with each coil L = 1.1mH, R = 0.1Ω and a capacitor of C = 20µF. These filters Y F are parts of our network. Loads 1 and 2 consume 12kW active power, while load L 3 consumes 5kW. No reactive power is consumed by the loads. We use a reactive power control loop as in [18], [22] and we inject random current and voltage measurement errors of a realistic size and power spectral density, see [18] for the details. As the transmission line impedance between the two generators within the same area is relatively low, their coupling is strong. Simulation results show that the effect of VF has 28 Fig. 2. The plots of active power output (in W) following the disconnection of load L 3 starting at t = 30s. Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF.
a beneficial impact when the load imbalance between the two areas is large and a strong tie line current is required to balance the two areas. In the case where the two areas are relatively balanced and therefore almost zero power flow is required between the two areas, a change in tie line impedance has little effect, as expected. We thus describe our results for the following two scenarios: • The tie line impedance is high and a sudden change in load occurs in one of the two areas. • The tie line impedance sees a step change at a high load imbalance.

1)
Step down change in load at high tie line impedance: Figures 2 and 3 show simulation results for the four generators after a sudden disconnection of load L 3 leading to a step change of the active power consumption in area 2 at t = 30s. This change occurs when switch S 1 is open, thus at a weak tie line between the two areas. When the VF is inactive, the active power output of the two generators in area 2 suddenly decreases by about 900W. Strong oscillations between the two areas are induced (See Fig. 2a). These oscillations do not occur when the VF is active (see Fig. 2b). Instead, within 3 seconds, the power output of the 4 generators converges without oscillations to equal values. In both cases, the overall active power output of all four generators has a transient lasting for about 7s, an effect that is due to the reactive power control loop [22]. Fig. 3 shows the generator frequencies for the same experiment with and without VF. In both cases, the frequency of the four generators increases at a slower rate dictated again by the rate of change of the excitation current control. The inter area oscillations show as modulations of the generator frequencies with an initial amplitude of 1 rad. These are fully damped when VF is active. Fig. 4 shows the grounded angles of generators G 2 to G 4 with respect to generator G 1 . The oscillations initially reach an amplitude of ±0.35rad. 2) Step down change in tie line impedance at high load imbalance: Data is shown for a sudden change of the tie line impedance between area 1 and area 2 (closing of S 1 at t = 40s). Fig. 5 shows the active power output of the generators. The change in tie line impedance from a high value to a low value triggers strong power oscillations between the two 4232   28   29  30  31  32  33  34  35  36  37  38   292   294   296 298 300 (a) VF inactive 28 Fig. 3. The plots of frequencies (in rad/s) following the disconnection of load L 3 starting at t = 30s. Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF.  Fig. 4. The plots of grounded rotor angles (in rad) with respect to generator G 1 , following the disconnection of load L 3 starting at t = 30s. Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF.
areas which slowly increase as long as the load imbalance persists. When the VF is activated, the active power output of the generators converge within 1s with no oscillations. In Fig. 6 we see that the amplitude of the frequency oscillations reaches ±1rad/s. A slow decay of the system frequency is observed over the course of 10s. Again VF efficiently damps the oscillations.
The grounded angles of the four generators can be seen in Fig. 7 to have an angular difference of 0.35rad before the reconnection of tie line Y 5 , which is due to the high phase shift induced by the weak inter area link between the two areas. In the absence of VF control, oscillations vary by   The plots of grounded rotor angles (in rad) with respect to generator G 1 , following a tie-line impedance step down starting at t = 40s. Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF.
±0.35rad. With VF, the grounded angles converge to zero within 2 to 3 seconds.
Note that the steady state power output of the VSMs only depends on P set and ω. With P set being equal for all 4 generators, the active power outputs of the synchronized generators reach equal values even if the load is not balanced between the areas (see [22]). Note also, that the value for D is chosen relatively low and that stronger damping of course can also be achieved by increasing D. This however comes at the cost of requiring the VSMs to supply large excess power in the case of ω deviating from its nominal value. An adequate choice of d and f should be made such that for a given (micro)grid d can stabilize the grid in normal operation, while f helps improving speed of convergence and decreases oscillations.
The above simulations show frequency oscillations of 0.55 − 0.65Hz in the case of a high tie-line impedance and twice faster oscillations for lower tie line impedance.
We mention that we did simulations with different communication time delays and we found that the behaviour remains good for non-uniform time delays up to 100ms. The details of this will be reported in the journal version of this paper.

B. The simplified model
In order to see if we get comparable results for the simplified model (the FEPS model), we have computed the Y -matrices of the grid by network reduction from the parameters of the full model. In order to calculate a i j , ϕ i j 4233 28  The stability criterion in Theorem 2.1 predicts stable operation of the model for a low tie line impedance. For L 3 disconnected, Γ min = 564s −1 and Γ crit = 394s −1 , while for L 3 connected, Γ min = 549s −1 and Γ crit = 368s −1 . This gives γ min = 42.7 • and γ max = 137.3 • and γ min = 40.8 • and γ max = 139.2 • respectively. At high tie line impedance, Γ min < Γ crit . This shows that the stability and synchronization criterion shown in the theorem is conservative, because the NRPS model shows stable behaviour even at high tie-line impedance. Note however, that the full model was unstable at low tie-line impedance and high load imbalance. 1) Step down change in load at high tie line impedance: Fig. 8 shows the frequency of the 4 generators for a load step change of −5000W in area 2. Startingt at t = 30s, the generator frequency increases from below 300rad/s to 310rad/s. Inter area oscillations appear with a period of 1.5s. VF efficiently damps these oscillations. (Note that the convergence of the frequencies to the new higher value is faster than in in the full model.) 39 40 41  42  43  44  45  46  47  48  49   295   300   305   39  40  41  42  43  44  45  46  47  48  49   295   300   305   G1 G2 G3 G4 Fig. 10. The plots of frequencies (in rad/s) following a tie-line impedance step down (starting from steady state). Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF.  Fig. 11. The plots of grounded rotor angles (in deg) with respect to generator G 1 following a tie-line impedance step down (starting from steady state). Subfigure (a) shows the plots without VF, while subfigure (b) shows the much smoother transition with VF Fig. 9 shows grounded rotor angles relative to G 1 . Here the NRPS model exhibits slightly higher amplitude than the full model. Note that the grounded angles of G 3 and G 4 converge to approximately +0.14rad/s, whereas in the full model this difference is close to zero. 2) Step down change in tie line impedance at high load imbalance: We show results for the simplified model switching at t = 40s from high to low tie line impedance while S 2 is closed. As in the full model, the change entails strong oscillations with frequencies around 3Hz. VF damps the oscillations of the network efficiently as seen in Fig. 10b and Fig. 11b. A small overshoot can be observed shortly after t = 40s.

V. CONCLUSION
This paper provides a proof of the stability of the NRPS model when employing VF, leading to the FEPS model. A singular perturbation analysis of the FEPS model shows that the same stability boundaries apply as in [34], however, the new result does not require to work with very small inertia. Our theoretical results are supported by simulations of the NRPS and FEPS models for a small groid comprising four VSMs of 9kW each. Comparison with the simulations for the full model of the grid shows good agreement between the behaviours, both with and without VF. In both cases VF was highly effective in damping oscillations. For a VF coefficient of F j = 3141.6Ws, the system showed transients decaying within 1-2 seconds.