Power-grid modelling via gradual improvement of parameters
Bálint Hartmann Géza Ódor Kristóf Benedek István Papp


AIP/123-QED




hartmann.balint@ek.hun-ren.hu
Institute of Energy Security and Environmental Safety, HUN-REN Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary


Institute of Technical Physics and Materials Science, HUN-REN Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary


Institute of Technical Physics and Materials Science, HUN-REN Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary

Department of Theoretical Physics, Budapest University of Technology and Economics, Budafoki út 8, H-1111 Budapest, Hungary

Institute of Technical Physics and Materials Science, HUN-REN Centre for Energy Research, P.O. Box 49, H-1525 Budapest, Hungary




§ ABSTRACT 
The dynamics of electric power systems are widely studied through the phase synchronization of oscillators, typically with the use of the Kuramoto equation. While there are numerous well-known order parameters to characterize these dynamics, shortcoming of these metrics are also recognized. To capture all transitions from phase disordered states over phase locking to fully synchronized systems, new metrics were proposed and demonstrated on homogeneous models. In this paper we aim to address a gap in the literature, namely, to examine how gradual improvement of power grid models affect the goodness of certain metrics. To study how the details of models are perceived by the different metrics, 12 variations of a power grid model were created, introducing varying level of heterogeneity through the coupling strength, the nodal powers and the moment of inertia. The grid models were compared using a second-order Kuramoto equation and adaptive Runge-Kutta solver, measuring the values of the phase, the frequency and the universal order parameters. Finally, frequency results of the models were compared to grid measurements. We found that the universal order parameter was able to capture more details of the grid models, especially in cases of decreasing moment of inertia. The most heterogeneous models showed very low synchronization and thus suggest a limitation of the second-order Kuramoto equation. Finally, we show local frequency results related to the multi-peaks of static models, which implies that spatial heterogeneity can also induce such multi-peak behaviour.


Power-grid modelling via gradual
      improvement of parameters
    I. Papp
    September 4, 2024
================================================================




Modeling power-grid systems has got a major importance in present days as transformation to renewable energy sources requires the complete re-design of energy transmission.
Renewable energy sources can be located quite far from their consumption points because urban and industrial structures do not follow physical constraints and capabilities. Important examples are the sea coast vs
inland divisions in the case of wind power. Ill-constructed high-voltage (HV) power grids can cause catastrophic damages to economies as it was demonstrated in recent history via the emergence of large blackout events<cit.>. The probability distributions of such events was found to be fat-tailed, exhibiting power-law (PL) tails very often
 <cit.>. To understand them, self-organized critical direct current (DC) models have been constructed <cit.> and have been shown to describe well the PL exponents of empirical values. However, many details could not be understood as power-grids work with alternating currents (AC) in which phase differences are the primary causes of the power-flows.




§ INTRODUCTION






AC modelling of power-grids have been proposed since the equivalence of swing equations to the second-order Kuramoto model was shown <cit.>.
Failures leading to blackouts have been studied by composite Kuramoto and threshold models <cit.> and the PL tailed cascade failures could be modelled by them <cit.>. Network topological features, which lead to desynchronization by network fragmentation and Braess paradox phenomena, have been identified <cit.>. We have shown that these are basically consequences of quenched heterogeneity, which can be mitigated by the enhanced fluctuations, that arise naturally in the neighborhood of synchronization transition points, where power grids self-organize themselves by the competition of supply and demand <cit.>.

In power grid systems, the focus is often on the phase synchronization of the individual oscillators <cit.> since their steady state is usually a stable limit cycle. To study the synchronization dynamics, several order parameters are used to characterize the dynamic state of the system. In the literature, there are numerous well-known Kuramoto order parameters, such as the complex order parameter <cit.>, the local order parameter measuring the phase coherence and its global variant <cit.>, a mean-field variant of the complex order parameter <cit.> or the one respecting network topology <cit.>.

Shortcomings of these metrics have been highlighted in a number of papers, most importantly by Ref.  <cit.>, who claim that existing order parameters are not fully suitable to characterize complex oscillator networks as they don't capture all transitions from incoherence over phase locking to full synchrony for arbitrary, finite networks. Hence a universal order parameter was also introduced, which captures partial phase locking, respects the topology of the network, and has been shown to increase monotonically with the coupling strength.

In this paper, we aim to address a research gap in the literature, namely to examine how different modeling assumptions regarding the heterogeneities of a power grid are captured by the different order parameters. In the different scenarios, to analyze and compare the dynamic behavior of the various models, we will use the frequency spread, the global order parameter, and the newly proposed universal order parameter by ref. <cit.> as the main measures. We also present a frequency analysis of the simulation results and confirm q-Gaussian distributions, matching real data distributions presented in one of our earlier work <cit.>.

The remainder of the paper is structured as follors. Section <ref> introduces the synchronization model and the twelve grid models. Section <ref> presents the results organized around five aspects. Finally, these results are discussed in Section <ref>, and conclusions are drawn.



§ DATA AND MODELS




 §.§ The synchronization model


Modeling power-grid systems come in different flavors, but at the heart of most approaches describing the time evolution lies the so-called swing equations <cit.>, set up for mechanical elements (e.g. rotors in generators and motors) with inertia. Mathematically it is formally equivalent to the second-order Kuramoto equation <cit.>, for a network of N oscillators with phases θ_i(t).

To investigate the effects of different parametrizations and to facilitate benchmarking with previous results, we used a more specific form <cit.>, which includes dimensionless electrical parametrization and approximations for unknown ones:

    θ̈_i+α θ̇_i=P_i/I_iω_S
    +P_i^max/I_i ω_S ∑
    _j=1^NW_𝑖𝑗 sin(θ_j-θ_i) + Ω_i  .

In this equation θ_i is the phase angle, ω_i=θ̇_̇i̇, is the frequency of node i, α is the dissipation or damping factor, W_ij is the coupling strength and P_i is the source/load power. Furthermore, I_i denotes the rotation inertia, ω_S the system frequency, and P_i^max the maximal transmitted power in the system. Note that we can also have an intrinsic frequency of nodes Ω_i = 50 Hz (in Europe), but it can be transformed out in a rotating frame and we have omitted it in the calculations. Our frequency results show the deviations from this value. 

If we know more details of the electrical parameters we can cast this into the form with real physical dimensions:

    ω̇_i = -D_i ω_i/M_i ω_S + L_i/M_i ω_S + ∑_j=1^NY_ij V_i V_j/M_i ω_Ssin(θ_j - θ_i),

where D_i has dimension of [·^2/^2] and describes the damping effect of element i in the system, L_i [·^2/^3] is the power capacity of node i, Y_ij=1/X_ij [1/] is the susceptance of lines, the inverse of reactance, V_i [] is the nodal voltage level and M_i [·^2] is the moment of inertia.

Topological heterogeneity of power grids is the result of two factors, (i) the structure and connectivity of the grid itself, and (ii) the heterogeneity of power line capacities and nodal behaviors, as it was presented in our recent work <cit.>. In the second-order Kuramoto equation, these varying properties are represented by the parameters L_i, Y_ij and M_i.

The time step resolution of the calculations was set to be Δ t= 0.25 and α=0.1 was used, similarly as in Refs <cit.>. Ref. <cit.> also used  α=0.4. In this way, the results of the Kuramoto equation become dimensionless.

To model station fluctuations, we have added a multiplicative, quenched noise to the source/sink terms of <ref>

    η_i,j = 0.05 ξ_j P_i/I_i ω_s ,

where ξ_j ∈ N(0.1) is drawn from a zero centered Gaussian distribution. To solve the equations of motion we used an adaptive, Runge-Kutta-4, from the package Numerical Recipies.

We investigated the standard synchronization measures of the phases R(t) and the frequency spread Ω(t), called the frequency order parameter. We measured the Kuramoto phase order parameter:

    z(t_k) = r(t_k) exp[i θ(t_k)] = 1 / N ∑_j exp[i θ_j(t_k)]  .

Sample averages over different initial fluctuations for the phases

    R(t_k) = ⟨ r(t_k)⟩

and for the variance of the frequencies

    Ω(t_k) = 1/N⟨∑_j=1^N (ω(t_k)-ω_jt_k))^2 ⟩

were determined, where ω(t_k) denotes the mean frequency within each respective sample at time step
t_k = 1 + 1.08^k, k=1,2,3.... 
Sample averages were computed from the solutions with hundreds of independent self-frequency realizations (i.e. η_i,j) for each control parameter.

Besides, we measured a more complex order parameter suggested for the second-order Kuramoto model, which claimed to accurately track the degree of partial phase locking and synchronization<cit.>

    r_iui(t_k) = 1/(∑_i,j^N w_ij) ∑_i,j^N w_ijcos(θ_i-θ_j)

and it's sample and temporal average in the steady state:

    R_uni=⟨ r_iui(t_k) ⟩


The fluctuations of the order parameters are measured by the standard deviations of the sample and temporal averages in the steady state, typically after 100 s transient time.



 §.§ The grid model

We chose the Hungarian high-voltage (132, 220, and 400kV) network to create our grid models. 
The network consists of 387 nodes and 640 edges, and its most important features are presented in Table <ref>. 
Cross-border transmission lines were reduced to their domestic terminals as sources or consumers, 
thus resulting in a standalone synchronous system. In the modeled loading state, 351 nodes behave as consumers and 36 as sources.

 

In order to study how modeling depth is perceived by the different order parameters, we created 12 different variations of the Hungarian network. These variations introduce heterogeneity to the parameters W_ij, L_i,  and M_i. The resulting representations thus range from completely homogeneous networks, which are the most widely covered in related literature, to completely heterogeneous ones, where electric parameters and nodal behaviors are defined using the actual data and measurements of the Hungarian system.
The following assumptions are used for the three parameters.


    
  * W_ij, coupling strength:
        
            
  * Identical value for each edge, the value corresponding to the largest thermal capacity (ampacity) limit in the system (approx. 1400 MW). This option represents the benchmark used by e.g. Refs. <cit.>.
            
  * Unique value for each edge, depending on their actual thermal capacity limits (range between 40 and 1400 MW).
            
  * Unique value for each edge, depending on their actual admittance Y_ij and voltage level.
        
    
  * L_i, nodal power:
        
            
  * The sources (L_i<0) are distributed equally among the nodes representing power plants and the consumers (L_i>0) are distributed equally among nodes truly representing consumption. This is a slightly modified assumption of Ref. <cit.>, where half the nodes correspond to consumers (L_i>0), while the other half to power sources (L_i<0).
            
  * Every L_i value is uniquely assigned, based on measured data (SCADA).
        
    
  * M_i, moment of inertia:
        
            
  * Aligning with the literature, we set a constant value for M_i, corresponding to a 400 gas turbine power plant as in <cit.>. 
            
  * We evenly distribute the moment of inertia among 400 and 220 nodes, which host the majority of synchronous generators (conventional power plants).
            
  * We evenly distribute the moment of inertia along all nodes of the model.
            
  * We set unique values based on measured data and whether the node actually hosts a synchronous machine or not.
        


In the following the 12 scenarios will refer to the different models as shown in Table <ref>.





§ RESULTS


In the following, results of the synchronization studies are presented in a sequential, interdependent way. First, the coupling strength was varied to compare how different order parameters display criticality (Section <ref>). Then the transient behaviour of the three order parameters, R, Ω and R_uni are analysed in sections <ref>, <ref> and <ref>, respectively. Section <ref> compares the frequency data of the simulations to grid frequency data presented in Ref. <cit.>.

Note that to assist the interpretation of the results, not all scenarios are displayed in all figures.


 §.§ Dependence of the critical point on the models


We assume that power-grid-like systems operate near the state of self-organized criticality (SOC) <cit.>,. This means that the system is not operated on 100% load capacity, but usually at a lower level. To mimic the not fully loaded behavior we can cast equation (<ref>) in the following form:


    θ̈_i+α θ̇_i=P_i/I_iω_S
    +P_i^max/I_i ω_S λ∑
    _j=1^NW_𝑖𝑗 sin(θ_j-θ_i)  .


The multiplicative factor λ in front of the interaction term is the chosen constant for the system and its value corresponds to different load levels. Alternatively it can be also said, that the initial P_i^max/I_i ω _S term transforms to  P_i/I_i ω_S, with arbitrary P_i power. In a mathematical sense this maps to a changing coupling strength, which allows us to identify the SOC behavior by analyzing the standard deviation of the Kuramoto order parameter.

To identify the cross-over point to synchronization (cf. fig. <ref>), we varied the λ parameter from 0.1 to 1.0 with steps of 0.1. For the sake of completeness, we performed analysis starting the system from a phase-ordered state, i.e. all oscillators have the same initial phase and some noise, or from a disordered phase, i.e. all the oscillators have random initial phase assignment. In the figures of this section, dash-dotted lines with unfilled markers will represent the phase-ordered states, and continuous lines with filled markers are the disordered states.  

We have chosen model scenarios 1, 5, and 9 for finding the optimal λ value, as these scenarios represent homogeneous nodal behavior with three different options for defining the coupling strength. We performed the statistical calculation on a minimum 2500 sample and a maximum 7500 for each λ value. The results are shown in Figs. <ref>-<ref> for R, Ω and R_uni, respectively.

Fig. <ref> shows that the behavior of the phase-ordered curves  (dash-dotted lines with unfilled markers) is non-trivial. As the curves of scenarios 1, 5, and 9 overlap, this implies that R does not capture the difference in modeling the coupling strength. Very similar results were obtained for Ω-s, as shown in Fig. <ref>.

Fig. <ref> displays that the R_uni order parameter indeed increases monotonically with the exception of the case of scenario 9, phase ordered initial condition, large λ values. Also, the standard deviation of phase-disordered results shows a peak at λ=0.4 and λ=0.3 for scenarios 5 and 9 and respectively 1. Phase-ordered results show now trivial behavior, with decreasing tendency with the exception of scenario 9. Finally, it is important to notice that R_uni captures the difference in modeling the coupling strength. Scenario 5, representing the unique thermal capacity of the lines shows a steep increase as the function of λ, showing that a heterogeneous grid structure tends to desynchronize under less severely weakened states.










 §.§ Kuramoto order parameter transient behavior


Fig. <ref> shows the Kuramoto order parameter, R in the transient from different initial conditions. It can be seen that the scenarios align in three curves. The highest R values are reached by scenarios 1, 5, and 9 (R≈0.7 for phase disordered initial conditions), then 2, 6, and 10 (R≈0.5), and finally by 3, 7, and 11 (R≈0.4). These results indicate that changing the coupling strength, W_ij does not cause visible changes in steady-state values of R if all other parameters are fixed. This implies that R is unable to describe the synchronization behavior by changing the characteristics of the transmission lines.
It can also be seen that increasing the heterogeneity of the moment of inertia, M_i, decreases the value of the order parameter, as expected.





 §.§ The Ω order parameter transient


Fig. <ref> shows the evolution of the order parameter Ω. Similarly to the results of R, the scenarios form three distinct curves. However, the order of these curves is the opposite of the Kuramoto order parameter. Scenarios 1, 5, and 9 exhibit the lowest results (ΔΩ≈6 for phase disordered initial conditions), followed by 2, 6, and 10 (ΔΩ≈40), and by 3, 7, and 11 (ΔΩ≈5000). These results again imply that changing the coupling strengths, W_ij does not affect the steady-state values of Ω too much if all other parameters are fixed, thus heterogeneity in the edge behavior seems to be irrelevant.





 §.§ The new order parameter R_uni


Fig. <ref> 
shows the evolution of R_uni from disordered (left) and from ordered (right) initial conditions. This time all nine curves are separated, thus it is validated that R_uni is able to capture different heterogeneity of the model. Similarly to the other order parameters, the highest synchronization is shown by the completely homogeneous model, Scenario 1. However, the results of the other scenarios are significantly different than before. In general, unique nodal behavior modeled through L_i and M_i decreases synchronization, which is an important advantage compared to the use of R, shown in Section<ref>.

When comparing the scenarios, an interesting behavior can be observed. If only the heterogeneity of the coupling strength, W_ij is increased, better synchronization is found in the case of the scenarios, where the values of W_ij are calculated from admittances (Scenarios 9, 10, 11) as compared to the ones, where actual thermal capacity was considered (5, 6, 7). If we compare curves with the same coupling strength, it is seen that the benchmark of literature, namely assuming the moment of inertia of a large power plant at each node, is a significant factor in reaching high synchronization. (E.g. (R_uni≈0.9 for Scenario 1 and ≈0.8 for Scenario 2.) However such models are demonstrably overly optimistic, especially with the increase of non-synchronous generation in the mix <cit.>. As R_uni is able to capture the different nodal behaviors and also shows their dependency on the coupling strength, its use over R should be favored, especially when using the Kuramoto model for grid case studies.





 §.§ Frequency distributions


Testing the predictive power of the Kuramoto-based modeling is a great challenge and has not been done on the quantitative level on realistic power-grids, according to our knowledge. Here we show node frequency results obtained by different levels of parameter approximations. We calculated the PDF-s at nodes obtained in the steady state form samples at the last 10 time steps and from thousands of independent realizations, which differ in the input/output power by the addition of the small quenched fluctuation values as shown in Eq.(<ref>).




These calculations were done for each scenario and for each λ.
We tried to fit the PDF-s with the 8 most popular distributions: Gaussian, exponential, Student's t, log-normal, Pareto, double Weibull, generalized extreme value, and beta, from the Python distfit package as well as by the q-Gaussian functions as this distribution was commonly fitted well other HV studies of AC electrical data <cit.>. Agreement with the q-Gaussian is remarkably good for lower λ values, see Fig <ref> in case of
Scenario 9.

For λ > λ_c we found multi-peak behavior as displayed on the inset of Fig. <ref> even though the width of the frequency spread decreases by increasing λ. Multi-peak frequency behavior is very common in the European power grid, especially in islands, like Great Britain, Ireland, and Mallorca. <cit.>. A numerical analysis based on the extension of the swing equations with a time-dependent damping factor could reproduce such global frequency fluctuations, suggesting that the system is wandering around the nominal 50Hz peak <cit.>.

Our swing equation solutions on the full Hungarian power-grid network suggest that even with constant parameters multi-peak frequency behavior can emerge when we overload the system with global power transmission above the synchronization point λ_c.
Thus, beyond temporally different behavior we can also find sub-peaks on static power grids, due to the network heterogeneity. Note, that in our previous large-scale simulations, we showed different synchronization behaviors at fixed control parameters in different communities of Europe <cit.>.
This may hint at the dangers in power grids with multi-peaks being out of optimal operation control even if the frequency spread is narrow.

For scenarios other than 9 fat-tailed distributions were also obtained. As they are very numerous publication of them will be published elsewhere.



§ DISCUSSION AND CONCLUSIONS


As it was shown in the comparative analysis of Section <ref>, the use of R_uni, proposed by Ref.<cit.> is encouraged to display the differences of heterogeneous power grid models. We found that in contrast to R and Ω, this order parameter is able to capture the difference of the coupling strength, showing higher synchronization values when W_ij is calculated based on the admittances of power lines. We also found that decreasing inertia of the system is more distinctly presented across the different scenarios. This feature of R_uni is especially advantageous for case studies, as the benchmark models of the literature tend to overestimate the amount of inertia in the system.

Considering the modeling depth, we found that completely heterogeneous models with unique nodal behavior based on SCADA measurements show very low synchronization and thus seem to be inappropriate for Kuramoto models. The underlying reasons could be (i) the ignorance of power losses and (ii) the ignorance of reactive power. To bridge these gaps, a promising approach was presented in Refs <cit.>, where the voltage magnitudes are incorporated in the Kuramoto model.

We provided local frequency results for the whole Hungarian power-grids, agreeing quite well with empirical measurements <cit.>. The calculated PDF-s of ΔΩ, with respect to the nominal 50 Hz exhibit similar width and shapes as those recorded and
published in <cit.>. For λ-s, which drive the system above the synchronization point earlier we observed community dependent synchronization <cit.>.These are related to the frequency multi-peaks of static models we report now. This means, that
they occur not only in systems with time dependent parameters, but spatial heterogeneity can also induce them and warn for over-driven power grids, in the sense that they are away from the optimal SOC behavior.


Bálint Hartmann acknowledges the support of the Bolyai János Research Scholarship of the Hungarian Academy of Sciences. Support from the Hungarian National Research, Development and Innovation Office NKFIH (K146736) is also acknowledged.




§ DATA AVAILABILITY STATEMENT

The data that support the findings of
this study are available from the
corresponding author upon reasonable
request.





§ REFERENCES



