Modeling the Impact of Modified Inertia Coefficient (H) due to ESS in Power System Frequency Response Analysis

: The advantages of increased penetration of distributed generation are also accompanied by several challenges, low inertia being one of them, which threatens the grid stability. An emerging approach to confront this problem is the introduction of so ‐ called virtual inertia (VI) provided by energy storage systems (ESS). In contrast to the already available literature which considers a conventional load frequency control (LFC) model, this work concentrates on a modified LFC model as the integration of a large portion of ESS changes the inertia constant ( H ) of a power system. A sensitivity function is derived that shows the effects of changes in H on the power system’s frequency response. With the help of the developed mathematical model and simulation results, it is shown that a difference in the actual and calculated values of H can deteriorate the system performance and economy. As one of the reasons for this difference is improper modeling of ESS in the LFC model, therefore, the study signifies the accurate calculation of H in the power systems having enlarged penetration of ESS.


Background and Motivation
A stable profile of the power system frequency response (PSFR) is among the primary requirements for secure and reliable operation of the power system. The PSFR depends on the balance between supply and demand of real power. Despite the fact that power system continuously experiences the fluctuation in power generation and demand, the PSFR is required to be maintained within a specific operating range [1,2]. Conventionally, the objective of frequency regulation (FR) is achieved with the help of spinning and nonspinning reserves during normal conditions. However, to avoid larger blackouts, underfrequency load-shedding or overfrequency protection schemes are deployed during contingency stages [3][4][5]. There are several studies that explain the frequency regulation process in the conventional power systems (e.g., [1,6]). However, the conventional methods of frequency regulation are neither ecofriendly nor cost-effective.
On the other hand, the modern power systems have an increased penetration of distributed energy resources to address the reliability, economy, and environmental concerns. However, the benefits of DERs are accompanied by several challenges, and low inertia is one of them. In the conventional power systems, the inertia inherently provided by the massive machines aids in the system's stability. The DERs, on the other hand, is based on the stationary or light rotating masses. Therefore, they lack in providing enough inertia to the power system when compared to the conventional power plants. To overcome this problem, a concept known as a virtual synchronous machine (VISMA) was introduced [7,8]. A VISMA is essentially a control system that uses an energy storage system (ESS) and a power converter to mimic the inertial and damping properties of a real synchronous machine. The inertia produced by VISMA is known as virtual inertia (VI), and it forms a basis to accommodate a larger share of DERs in the modern power systems without conceding the system's stability. A conceptual block diagram of VISMA is shown in Figure 1.
There are several studies that show the effectiveness of VISMA for FR services, and it is expected that modern power systems will have an increased role of VI for the grid's stability [9]. For example, Ireland and the United Kingdom (UK) aim to increase the introduction rate of DERs up to about 40% by 2020. Their plan also includes to export a huge amount of wind power to the UK from Ireland [10]. However, due to integration of several inverters, it is projected that the inertial constant of the UK power grid could reduce to 70% by 2030 [11]. In this context, VI and associated challenges are a topic of great interest among researchers. A working group supported by the 6 th European Research Framework program presented useful research on this topic by providing the mathematical models, control system designs, and experimental results for VISMA hardware in the loop [12][13][14]. A comprehensive description of existing topologies for VISMA systems is presented in [15]. An improved dynamic response and dispatchable behavior of wind turbines are reported in [16]. This study presents a VISMA control method for the full converter wind turbine and assumes minutelevel energy storage in the dc link as the energy buffer. A recent study [17] models the role of VI and demand response programs for the FR process of a power system. The results of this study show that the VISMA can play a vital role in stabilizing the frequency response of a power system. In [18], a model predictive control based on virtual inertia is proposed for the FR process. The proposed control methodology demonstrates a robust performance even in the presence of a high share of intermittent RES. A control scheme is proposed in [19] to suppress the deviation in a system's voltage and frequency using VISMA. The study in [20] extends the VISMA control for voltage-source converters to modular multilevel converters and shows accurate results with improved efficiency.

Contribution and Organization
From the previous discussion, it is evident that the magnified capacity of ESS modifies the inertia constant ( H ) of a power system. This observation calls for an in-depth investigation of the impact of the inertia constant of the power system on the FR process, as H is an important parameter of the load frequency control (LFC) model and the power balance equation. In solving the LFC problem, H is typically considered as of constant value. However, this assumption is not justified for the power systems having an enlarged capacity of ESS and can result in undesired results, as explained in the subsequent sections. Hence, in comparison to the existing literature, the contribution, assumptions, and highlights of the presented work are summarized as follows: 1. The major part of existing literature that addresses the FR process in the presence of ESS considers a conventional load frequency control (LFC) model [8,15,18,21]. Some other studies, such as [17,22], adopt modified LFC models to attain the objective of FR while considering the connection or disconnection of a large load or generation. In these studies, however, there is no to very little contemplation over the intrinsic changes to the power system due to the addition of energy storage devices. In contrast, the proposed work investigates the FR process by introducing a novel stability analysis which considers frequency deviation not only as a function of load disturbance but also a function of H . 2. A sensitivity function is defined and derived to show the impact of variations in H on the PSFR in the developed LFC model. The results are compared with the conventional model. It is notable that the sensitivity function developed in this work emphasizes on H in comparison to studies such as [23][24][25][26][27], where sensitivity function is developed for the load-damping coefficient or other system parameters. 3. The results show that errors in the calculation of the inertia constant affect the system performance and stability. Although, for essentially stable systems, the effects of miscalculated inertia constants are overshadowed by external disturbances, and the existing controls system of the power system is robust enough to overcome these effects. However, the consequences are adverse when the system is unstable or has a low stability margin. Alternatively, it means that the frequency protection devices may malfunction due to miscalculated H . 4. Please note that the presented work mainly focuses on the modified LFC model considering changes in H and analytical results obtained thereof in a MATLAB/Simulink environment. To find out how wind and solar systems affect the value of the inertia constant and similar topics, please refer to other resources, such as [28][29][30][31]. Further, implementation of the proposed model in a hardware-in-the-loop setup or using a detailed machine model is not reported in this work.
The paper is organized as, in Section 2, the derivation of the sensitivity function with respect to H for the frequency response analysis is presented for two typical types of power systems. A detailed analysis of PSFR using the total differential equation is presented in Section 3 considering two cases: when the system is essentially stable and when the system is unstable. The simulation results are presented and discussed in Section 4. Section 5 delivers a summary and conclusion of this study. A description of the symbols used in the paper is provided in the Abbreviations and Nomenclature section at the end of the paper.

Derivation of Sensitivity Function
In this section, a sensitivity function is defined and derived to show the impact of variations in the inertia constant on the frequency deviation of a power system. A change in the inertia constant is denoted by H . For the sake of generality, two typical categories of power systems (i.e., with primary control only and with primary and supplementary control) are considered for the derivation.

Power System with Primary Control
For the frequency control analysis, the block diagram of a power system containing a non-reheat turbine and a primary or droop control mechanism is shown in Figure 2 [15,23]. This model of the power system does not have a supplementary control, which will be explained in the subsequent sections. A mathematical relation for the calculation of the frequency deviation due to a step load disturbance is provided in Equation (1), as follows [1]: Figure 2. Block diagram of a power system with primary control only.
The term t T and g T represent the time constants of the non-reheat turbine and governor, respectively. The mathematical calculations and derivations carried out in this work are equally valid for other types of turbines. Considering other types, such as hydro or reheated turbines, will require a modification in the turbine transfer function (TF), as described in Equation (2) [32].
H is calculated to demonstrate a potential variation in frequency due to a change in H . For the old conventional power systems where the penetration level of ESS was small, the effect of a change in H was not significant in comparison to the external disturbances, and the available control of a power system was robust enough to handle the minor change in H due to small-sized ESS. However, this difference cannot be overlooked in the present and future power systems where the inertia constant changes significantly during a time period due to an intensified involvement of the ESS.
A detailed derivation of Equation (3) is presented in Appendix A.
Next, a unit less sensitivity function ( ) f H  is defined in Equation (4) for the relative change in the frequency deviation with respect to a relative change in the inertia constant.
For the power system without supplementary control,  f H is calculated according to Equation (5), as follows:

Power System with Primary and Supplementary Controls
The LFC model of a power system containing supplementary control in addition to primary control is shown in Figure 3, [15,23].
Following the steps of the previous subsection, the frequency deviation and the sensitivity function are calculated according to Equation (7) and Equation (8), respectively. Appendix B provides a detailed derivation of the sensitivity function given in Equation (8) for the power system shown in Figure 3.
Again, the unit less sensitivity function ( ) f H  is calculated for the power system with supplementary control, as given in Equation (9).
The formulation of sensitivity functions for the two models of the power system is the same as shown in Equations (3) and (8). Similarly, Equations (5) and (9)

Stability Analysis
For the sake of simplicity, the conventional models consider the frequency deviation to be associated with an external load disturbance while the fixed values are of the power system coefficients [33,34]. In terms of a differential equation, this approach is represented by Equation (10), as follows: However, the effect of changes in H on the PSFR is not considered in Equation (10) which can lead to undesired results, as shown in the following section. Based on the discussion presented in the previous sections, this study recommends that  ( ) f s must be a function of H and  ( ) Considering a mutual independence between  ( ) d P s and H , formulation of the frequency deviation is modified to include the role of, H as presented in Equation (11): Using the formulation of   ( ) / f s H as already calculated in Equations (3) and (8), the following relationship is obtained for the modified model: control only using Equations (6) and (7), as follows. It is notable that the given derivation is also applicable to the system with primary and supplementary controls and results in the same formulation, as shown in Equations (13)-(15): Combining Equations (12) and (13), Equation (14) is obtained as follows: The time-domain representation of the frequency deviation is presented in Equation (14) by integrating the inverse Laplace transformation of Equation (13), as follows: Having a time-domain representation of the frequency deviation, the stability analysis is performed for the PSFR considering two situations of the power system, as explained below.

Stable Power System
The classical methods of the root locus plot or Routh-Hurwitz stability criterion provide sufficient information about the stability of a linear system, and they are frequently used in control system applications [35]. The characteristic polynomial of a stable power system is a Hurwitz polynomial, which has all the roots located in the left half of the complex plane. It means that a finite or bounded power disturbance  ( ) bounded. It is mathematically represented by Equation (16), which means that  ( ) f t is within a small band, as the power system is essentially stable.
Similar to the conventional model, the norm of TF is bounded in the developed model (Equation (11) (17) represents a small number whose values are less than its unity for a stable power system and ( ) u t represents a step function.
The maximum frequency deviation for the conventional model of Equation (10) in the timedomain is bounded by the term shown in Equation (17) for a stable power system.
It is evident from Equations (17) and (18) that  ( ) f t is bounded for the two models of a stable power system; however, their boundaries are dissimilar.

Unstable Power System
If the TF of a power system has the right-half plane poles (i.e., the characteristic polynomial of the power system has zeros in the right-half plane in the s -domain), the power system is considered to be unstable. It implies that a finite power disturbance  ( )  (19), where  represents a very large number whose value is greater than its unity.
For the conventional LFC model, which considers the frequency deviation as a function of load disturbance (  ( ) d P s ) only, the lower bound of the frequency deviation in the case of an unstable system is calculated using Inequality (20) [23]. It is notable that the frequency deviation is directly proportional to  and  ( ) On the other hand, the lower limit of the frequency deviation for an unstable power system can be calculated using Inequality (21). A complete derivation of Inequality (21) is provided in Appendix C and follows the steps explained in the preceding subsection.
Comparing Inequality (20) and (21), it is clear that the frequency deviation in the time-domain is very large even if the values of H and  ( ) d P s are small per unit. This observation is true for the two LFC models (i.e., the conventional model and the model proposed in this work). However, the boundaries of the frequency deviation are not the same for the two models. In the conventional model, the lower bound of the frequency deviation in the case of an unstable system is much smaller than the lower bound given in Inequality (21) that has a  2 term. In other words, the effect of changes in H on the frequency deviation cannot be ignored in this situation.

Calculation of Maximum Frequency Deviation
Though the frequency deviation is bounded for an essentially stable power system, it is always desirable to calculate the maximum frequency deviation, which is denoted by  max ( ) f t in this work.
The value of  max ( ) f t provides vital information about the power system. From a design and operational point of view, the criterion of the power system stability must fulfill the requirements associated with  max ( ) f t .
For a stable power system, the partial time-derivate of ( ) f t must be zero at the point of the maximum deviation [5,21]. This is mathematically shown in Equation (22), which assumes no initial frequency deviation (i.e.,  For the two types of power systems with and without supplementary control, the maximum frequency deviation is calculated according to Equations (23) and (24), respectively, as follows: For an oscillatory response of  ( ) f t , there are several values of t that satisfy Equations (23) and (24). In this case, the smallest value of t must correspond to  max ( ) f t as the first swing in a stable power system corresponds to the maximum frequency deviation.

Case Studies
The simulation setup is divided into four case studies, as described below: 1. Case I: A stable power system with primary control. 2. Case II: An unstable power system with primary control. 3. Case III: A stable power system with primary and supplementary controls. 4. Case IV: An unstable power system with primary and supplementary controls.
The parameters of the test systems are shown in Table 1 [36]. Appendix D shows the MATLAB implementation of the test system. A step load disturbance of 0.05 p.u. is applied in all the test cases. The changes in H are made with a step of 15%, and the possibility of positive and negative changes is considered in this work. As already discussed, several systems are facing the reduction of H as high as more than 50%. Therefore, the simulation cases consider up to ±60% variations in the inertia constant.  Figure 4 shows the effects of changes in H on the PSFR for a stable power system having primary control only. Since the power system does not have an integral control action, therefore, the frequency deviation does not go to zero, as shown in the results. The key points of this simulation study are presented as follows:

Case I: A Stable Power System with primary control
1. In a stable power system with only primary control, the level of frequency deviation changes with H , as shown in Figure 4a. For the same load disturbance, the maximum frequency deviation is different from the one calculated using a fixed value of H (i.e., the case when ). 2. The observation presented in the preceding point has significant implications over the conventional model, which does not consider the effect of changes in H over the frequency deviation. From Figure 4a, it is evident that the error in calculation of the amount of frequency deviation could be as high as 43% for the extreme cases (i.e., ∆ 60% and ∆ 60%) and 30% for ∆ 0% and ∆ 60%. It is notable that, from a practical perspective, only a major portion of ESS devices can introduce such a significant difference. However, the results presented in Figure 4a indicate the importance of considering variation in the inertia constant. 3. Additionally, the observation presented in the previous points has practical applications as well.
For example, the capacity determination of spinning/nonspinning reserves is also dependent on the calculated frequency deviation; therefore, an error in the H estimation can lead to over/undercapacity reserves. Similarly, the reliability of frequency relays is also jeopardized if H is miscalculated or there are changes in its value due to the virtual inertia provided by ESS. It is once again iterated that, in conventional LFC models, the value of H is considered to be constant and does not change during an operational time-period. 4. The rate of change of the frequency deviation in Figure 4b is zero when the frequency deviation in Figure 4a touches a maximum or minimum point. As discussed earlier, the first zero-crossing of   ( ) / f t t corresponds to  max ( ) f t in a stable power system. This is demonstrated for a simulation run with 60%. The maximum frequency deviation occurs at 0.692 s , as verified by the curve shown in Figure 4b. For the sake of clarity, this observation is highlighted for only one case, although it is valid for all simulation cases.

Case II: An Unstable Power System with Primary Control
The effects of changes in H on the PSFR for an unstable power system primary control only are shown in Figure 5. Since the power system is not stable, therefore, the load disturbance takes the frequency deviation to higher levels. The instability is modeled by converting the left-half plane pole of the turbine TF to a pole in the right-half plane (i.e., Further, the generation rate constraint (GRC) of 10% does not allow extreme changes to the power, though the frequency response remains unstable. The key points of this simulation study are presented as follows: 1. Although the frequency response is unstable in all cases, however, the level of the frequency deviation changes with H , as shown in Figure 5a. In some cases, the response of the frequency deviation is quicker than the other cases. 2. The observation presented in the preceding point also has practical considerations. For example, the frequency relays can mal-operate if H is miscalculated or there are changes in its value due to the virtual inertia provided by ESS. 3. As the frequency deviation goes much beyond the operational limits, the shown simulation time of 10 s in this study provides enough information. Since the power system is unstable, the frequency deviation should (theoretically) keep on increasing, and its time derivative will not be zero. In any case, the practical power plants have control, and protective devices (e.g., frequency relays or generation rate limiters) that do not allow their operation after certain limits on  ( ) f t and   ( ) / f t t are met.  Figure 6 shows the effects of variation in H on the PSFR for a stable power with primary and supplementary controls. Again, the changes in H are made with a step of 15%. In contrast to previous cases, the power system has an integral control action; therefore, the frequency deviation goes to zero, as shown in Figure 6a. The summary of observations for this simulation study is presented as follows:

3. Case III: A stable power system with primary and supplementary controls
1. For the same load disturbance, the level of frequency deviation changes with H . As shown in Figure 6a,  max ( ) f t and max t are different for all cases. This observation calls for the need of accurate information about the power system's inertia constant. It is notable that, as long as the power system reserves have enough capacity, the stability of a power system is not jeopardized.
In any case, a significant difference in the actual and calculated value of H may result in insensitive or oversensitive behavior of the protective relays. It is once again explained that a huge amount of ESS can introduce such a significant difference. 2. Whenever  ( ) f t reaches a maximum or minimum point,    ( ) / 0 f t t , as shown in Figure 6b.
As discussed earlier, the first zero crossing of   ( ) / f t t corresponds to  max ( ) f t in a stable power system. This is demonstrated for a simulation run with 60%. The maximum frequency deviation occurs at 0.66 s, as verified by the curve shown in Figure 6b. For the sake of clarity, this observation is highlighted for only one case, although it is valid for all simulation cases.

Case IV: An Unstable Power System with Primary and Supplementary Controls
The effects of changes in H on the PSFR for an unstable power system with primary and supplementary controls are shown in Figure 7. The instability is modeled by converting the left-half plane pole of the turbine TF to a pole in the right-half plane, as explained earlier. Since the power system is not stable, therefore, the load disturbance keeps on increasing the frequency deviation. The following points summarize the important observations of this simulation study: 1. Although the frequency response is unstable in all cases, however, the level of the frequency deviation changes with H , as shown in Figure 7a. In some cases, the response of the frequency deviation is quicker than the other cases. In other words, the frequency relays can trip earlier or later than expected if there is a difference between the actual and calculated H . 2. As shown in Figure 7b, the rate of change in the frequency deviation will not go to zero for an unstable power system. This is due to the existence of a right-half plane pole in the system TF, which makes the error grow. It is once again explained that, due to the deployment of the GRC, the response does not shoot instantaneously. The results shown in Figure 7 are for a demonstrative purpose to show the significance of an accurate calculation of H for an improved PSFR. Further, the practical systems do have watchful protection that does not allow the generator to operate beyond set limits [6].

Conclusion
The modern power systems have increased the role of ESS, which can provide VI with the help of power electronics-based power converters. These sources of VI are seen as an alternative of the synchronous generator in the context of the power grid's stability enhancements. There exist several simulation studies that show the promising role of VI in frequency profile improvements of power systems using the well-established LFC model. These studies, however, do not consider the changes in the inertia coefficient due to significant penetration of ESS. This work, on the other hand, provides a modified LFC model which includes the role of increased VI. The sensitivity function derived in this work shows the importance of the inertia coefficient in the LFC model. The stability analysis and simulation results show that the effects of differences in the calculated and actual value of H could be detrimental, as it can result in the mal-operation of frequency relays or improper calculations of power reserves. The findings of this study call for the proper modeling of ESS for VI and an accurate and dynamic calculation method of the inertia coefficient of the power system.
As a future work, experimentation of the proposed model in a hardware-in-the-loop configuration could lead to interesting results. Moreover, the presented work can be extended by taking into account the nonlinearities of the power system and considering complex multi-machine, multi-area configurations. Additionally, the presented model can be augmented with consideration of changes in the damping coefficient of the power system.

Appendix D
A screenshot of a simulation setup showing the deadband and the generation constraints of the governor and the turbine, respectively. The values of these parameters are provided in Table 1.