Real-Time Estimation of Microgrid Inertia and Damping Constant

The displacement of rotational generation and the consequent reduction in system inertia is expected to have major stability and reliability impacts on modern power systems. Fast-frequency support strategies using energy storage systems (ESSs) can be deployed to maintain the inertial response of the system, but information regarding the inertial response of the system is critical for the effective implementation of such control strategies. In this paper, a moving horizon estimation (MHE)-based approach for online estimation of inertia constant of low inertia microgrids is presented. Based on the frequency measurements obtained in response to a non-intrusive excitation signal from an ESS, the inertia constant was estimated using local measurements from the ESS’s phase-locked loop. The proposed MHE formulation was first tested in a linearized power system model, followed by tests in a modified microgrid benchmark from Cordova, Alaska. Even under moderate measurement noise, the technique was able to estimate the inertia constant of the system well within ±20% of the true value. Estimates provided by the proposed method could be utilized for applications such as fast-frequency support, adaptive protection schemes, and planning and procurement of spinning reserves.


I. INTRODUCTION
The inertia of modern power systems is in continuous decline due to the increased utilization of renewable energy sources (RESs) displacing traditional rotational synchronous generation. Low inertia power systems are prone to large rateof-change-of-frequency (ROCOF) following a power imbalance, causing large frequency deviations as a result. Proper assessment of the system inertia is of critical importance for system operators to deploy effective strategies to supplement the lost inertia. The inertial response of a power system is the product of instantaneous release/absorption of The associate editor coordinating the review of this manuscript and approving it for publication was Elisabetta Tedeschi . power, typically from synchronous generators, in response to a power imbalance (measured as the deviation in frequency from nominal). The inertia of a power system is thus dependent on the synchronous generators that are online at any given time. With non-synchronous generation sources such as photovoltaics (PVs) and wind continuously displacing the traditional generators, it is challenging to estimate the inertia of a power system at any given time [1]. The variable nature of these RESs adds another degree of uncertainty as the inertia constant of the power system becomes time-varying [2].
To counteract the decline in system inertia, various fast-frequency support mechanisms have been proposed [3], [4]. These strategies inject/absorb power into/from the power system in response to a frequency event, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ much like during the inertial response stage of a synchronous generator. When deploying such strategies, an energy reserve is required to emulate the inertial behavior [5], which could come in the form of energy storage systems (ESSs) and/or even curtailed operation of RESs [6]. For a system operator, accurate knowledge of the inertial and damping behavior of the system is desirable to coordinate appropriate control strategies and protection schemes. Furthermore, it is imperative for the system operator to have enough reserves in place for any plausible contingency in the system [7]. The Electric Reliability Council of Texas (ERCOT) has proposed the concept of ''critical inertia'' as the minimum inertia at or above which the system can be operated reliably warranting the need to actively monitor the system inertia [8]. Situational awareness regarding the inertial response available in the system is thus critical information for the system operator.
Most methods proposed in the literature for inertia estimation are offline (post-event), developed for large interconnected power systems [1], [9]- [12]. In these methods, the data from large disturbances in a power system are logged and, after significant signal-processing, used to estimate the inertia constant. For example in [9], the inertia constant of the system is estimated by solving the swing equation based on the frequency transient measurements. A polynomial approximation with respect to time is applied to the transient frequency measurements to isolate the oscillations and noise in the measurement. The damping of the system has been neglected in this approach. Furthermore, these techniques are known to be susceptible to the identification of the exact time of frequency event onset and the order of the polynomial approximation [10]. System identification approaches using ambient frequency measurements have also been developed [1], but these techniques require computationally expensive offline processing to extract the explicit value of the inertia constant. Although beneficial in traditional synchronous generation based power systems, offline estimation approaches are not suitable for real-time adaptive control systems/protection schemes as the information may arrive late. Furthermore, there may be cases where the system operator does not have enough visibility/information regarding the inertial response of all the sources in the power system. For the stochastic and low-inertia nature of future power systems, online estimation of system inertia and damping constants is key.
Several online approaches have also been presented in the literature. Online estimation of system parameters such as inertia and damping of the power system provides insights into designing adaptive control systems and protection systems, such as under frequency load shedding [13]. Moreover, this allows system operators to optimize the resources to maintain the reliability and resiliency of the power system [14]. An online technique for estimating the inertia constant that uses sliding window-based approach to filter the noisy frequency measurements from phasor measurement units (PMUs) was proposed in [15]. The proposed method again uses the swing equation as the basis to determine the inertia constant. Similarly, in [16], online inertia estimation based on an auto-regressive moving average exogenous (ARMAX) input model is proposed. However, the method requires the data to be collected over a certain time-period (30 s recommended in the paper), which can introduce significant delays in the estimate. In [17], system identification techniques relying on ambient frequency measurements are proposed. However, techniques that depend on ambient frequency measurements may be unable to distinguish between voltage and frequency dynamics and even other electro-mechanical dynamics in the power system [16]. In [18], online identification of inertia response of a system is performed by using a system identification approach on PMU recordings. The method is also shown to require powerful computer resources which is not suitable for embedded application for an ESS in microgrids. The method proposed in this paper is based on frequency estimations from a local phase locked loop (PLL) of the ESS instead of relying on frequency measurements collected over an extensive PMU network as in most other methods. It is expected that in future power systems (including microgrids) ESSs would be widely deployed and could be leveraged for inertia estimation through local PLL based frequency measurements. Similarly, methods based on statistical models have also been explored, but they require significant offline training and calibration using historical data with a number of frequency events [2], [19].
In this paper, we propose a real-time inertia estimation technique based on moving horizon estimation (MHE) for low-inertia microgrids that utilizes local frequency measurements from a PLL of an ESS. The proposed inertia estimator can be implemented within the controller of an ESS and the estimates can be leveraged by the ESS operator to provide fast-frequency services. The main contributions of this paper are: 1) Designed an online inertia estimation technique for low-inertia microgrids to provide the ESS system operator with situational awareness regarding the inertial response of the microgrid. The proposed technique utilizes a MHE approach that relies only on local measurements from a PLL of an ESS. In our previous work [20], the MHE approach was utilized for frequency estimation in microgrids while demonstrating the computational feasibility of implementing MHE algorithms in applications with sampling times in the range of ms. In this paper, the MHE framework is further expanded for inertia and damping constant estimation. 2) Demonstrated the proposed approach in a linearized power system model utilizing a non-intrusive active power excitation signal to induce small changes in the frequency of the microgrid. Measurement noise can have a significant impact on the accuracy of the estimates, especially in the case of microgrids. As most existing work in the literature does not analyze the impact of measurement noise on the inertia estimates, the performance of the proposed technique is analyzed under noisy conditions. Even under noisy PLL measurements, it is demonstrated that the proposed technique can estimate the time-varying system inertia and damping constant. Furthermore, the proposed technique is also shown to perform well when the noise distribution is not Gaussian. 3) Introduced a modified microgrid benchmark from Cordova, Alaska where the proposed approach was tested. The paper is organized as follows -Section II introduces a linearized model representing the frequency dynamics of a microgrid, Section III then introduces the concepts of the MHE approach. The simulation setup is described in Section IV followed by the results and analysis in Section V. The proposed approach is implemented in a microgrid benchmark from Cordova, Alaska, in Section VI. The paper is concluded in Section VII.

II. MODELING OF MICROGRID FREQUENCY DYNAMICS
In this section, an approximate model for the frequency dynamics of a microgrid is introduced. Such a model shows similar behavior to the real frequency dynamics observed in a power system even though the derived model makes several simplifying assumptions and is not a comprehensive representation of microgrid system dynamics. An approximate model is attractive for use with an online MHE to reduce computational cost and ensure optimization convergence [20]. It is shown later that such an approximate model is still able to provide good parameter estimates in a realistic benchmark. Rotational generation is considered to be the major power source and renewable generators supplement the power demand. This allows the frequency dynamics of an isolated microgrid to be modeled using the swing equations based on the concept of an equivalent generator [21], [22]. The frequency dynamics of the entire microgrid (possibly with a number of generators) can thus be represented with a single equivalent generator with appropriate parameter approximation. One can then estimate the inertia and damping constant of this equivalent generator to estimate the overall inertia and damping constant of the microgrid.

A. INERTIA AND DAMPING CONSTANT
When subjected to a disturbance, the generators in a microgrid experience small variations over an average frequency. The average frequency is defined for the center of inertia (COI) of the system as follows: (1) where N represents the total number of generators in the microgrid, H i represents the inertia constant of the i th generator, and ω i is the angular frequency of the i th generator. Similarly, the total inertia constant of the system H can be defined as: where S i is the apparent power rating of the i th generator. Based on these definitions, the linearized swing equation of the microgrid frequency dynamics is given by: where ω 0 is the nominal system frequency, ω represents the ROCOF, P m is the change in the total mechanical power, P is the change in output power of an ESS unit, and P l is the net change in the electrical load of the aggregated microgrid system. Next, if we consider the frequency dependency of the aggregated loads in the microgrid system and use a per-unit (p.u.) system, the swing equation takes the following form: where ω is the change in system frequency, M = 2H is the inertia constant, and D is the damping constant of the system. The damping constant represents the frequency dependence of the power system loads [23]. For instance, a D of 2% means that a 1% change in frequency would result in a 2% change in system load. The net load change P l is assumed to be a disturbance from the point-of-view of the ESS unit. Also, as this is an isolated microgrid, the tie-line power flows are not modeled [21]. This differential equation shows that the inertia and the damping constants define the inertial response of the microgrid system. Estimating these constants in real-time will thus allow the implementation of fast-frequency support mechanisms for ESSs and the design of adaptive control schemes.

B. MODELING OF THE FREQUENCY CONTROL LOOPS
The typical frequency control loops in an isolated power system are shown in Fig. 1. The primary frequency control is implemented through governors in the system to stabilize the frequency change. The secondary frequency controller then removes the steady-state error as described in [23]. The type of governor and the dynamics of the turbine itself affect the dynamics of the frequency response of the system. It is assumed that in general the turbine-governor dynamics and the secondary control loop can be represented by the following set of differential equations [24]: where T g is the turbine-governor time constant of the aggregated generators, R p is the aggregated droop constant, K i represents the integral gain of the secondary control loop, P s is the secondary power, and δ is the change in the phase angle of the equivalent generator. The turbine-governor dynamics model used may be more complex if the system is dominated by other generation types. The proposed estimation technique remains equally applicable to other modeling assumptions. Based on (4), (5), and (6), the following differential equation describing the overall frequency dynamics of the isolated power system can then be derived: The derivative term Ṗ is neglected to simplify the predictive model to be used by the MHE. Such an approximation ensures optimization convergence and efficient computation. Moreover, because the MHE only requires an approximate model, such an assumption will not affect performance [20], [24]. The state-space representation of the differential equation is then given by:

III. MOVING HORIZON ESTIMATION AND PROBLEM FORMULATION
In this paper, the predictive model described in (8) is used in combination with online measurements for frequency and ROCOF from a PLL of an ESS to identify the equivalent inertia and damping constants of the system. There are various approaches to perform state and parameter estimation of dynamic systems, such as Kalman filters and extended/unscented Kalman filters [25], [26]. Most of these techniques make an assumption on the noise distribution of the measurement signals. Particle filters have also been widely discussed for online state and parameter estimation of dynamic systems, but they are known to be sensitive to errors in the initial guess, and thus may have larger convergence times [27]. Furthermore, if there are physical constraints on the estimated parameters, these techniques lose accuracy [28]. MHE, on the other hand, does not require assumptions on noise distribution and can explicitly handle constraints. In this section, the fundamental concept of MHE is first introduced. This is followed by the problem formulation when implementing MHE for inertia and damping constant estimation.

A. CONCEPTS OF MOVING HORIZON ESTIMATION
MHE is an online optimization-based estimation technique.
In MHE, past measurements are collected over a finite horizon and then the state and parameters at the current sampling time are estimated based on minimizing a cost-function while satisfying the constraints on the states and the parameters. A number of related work [27], [29], [30] show the wide applicability of MHE in linear and non-linear systems. Recent advances in efficient algorithms and optimization solvers have made the real-time implementation of these dynamic optimization problems feasible [28], making MHE suitable for embedded applications. The concept of MHE is illustrated in Fig. 2. The crosses (X) represent the measurements taken from the system, while the dashed line is the system output based on a predictive model. A measurement window of fixed length L is shifted one-step ahead in each sample time.
At the initialization stage, the window is initially empty so the window starts to move once the appropriate amount of data is collected. As the window moves forward, the newest measurement is incorporated into the estimation window while the oldest measurements are removed. The moving window allows the MHE to approximate a full-information estimator while maintaining computational tractability [28].

B. FORMULATION OF MOVING HORIZON ESTIMATION
The MHE, which will be implemented in an ESS, is formulated in this section. To ensure the convexity of the optimization problem, the MHE is formulated as a quadratic program (QP). This allows for the use of efficient solvers to be implemented to solve the online optimization problem [31]. Let us define L as the backward time-horizon and x k = [ δ k ω k ω k ] as the states of the system at the discrete-time instant k. Then the MHE problem will take the following form: where J L is the cost-function to be optimized and y k is a vector representing measured states. The measured power output (the excitation signal injected by the ESS) is denoted by p k , while p k represents the predicted power output from the ESS unit. The discretized state-space matrices A d , B d , and C d are obtained using the Zero-Order Hold (ZOH) method. The output matrix is given as C d = diag(0, 1, 1) as only the frequency deviation and the ROCOF are measured using a PLL. Equation (9b) defines the dynamics constraints used by the MHE while constraint (9c) and (9d) limit the search range of the inertia estimates to realistic values. The cost function consists of two terms weighted by matrices V and W . The matrix V is defined as V = diag(0, V 22 , V 33 ). The element V 11 is set to zero as the δ k is not utilized when implementing the MHE. The matrix W has only one term as there is only a single control input. The first term in the cost function J L penalizes the difference between the measured outputs and the predicted outputs. Similarly, the second term accounts for actuation errors [32]. Solving this optimization problem at each sampling time yields the parameter estimatesM andD. The MHE estimator also yields p k which is the estimated power output from the ESS.

C. IMPLEMENTATION OF MOVING HORIZON INERTIA ESTIMATION
The MHE is implemented using the ACADO Toolkit for MATLAB [33]. ACADO for MATLAB is an open-source toolbox that provides a general framework to implement dynamic optimization problems. The algorithm for the online estimation procedure is presented in Algorithm 1. Every time a new estimate is desired, the microgrid frequency dynamics are excited using an excitation signal p k (in this case a pulse train signal). The frequency of acquiring these new estimates will depend on the application. Once the time horizon L and the sampling time are defined, the frequency and ROCOF measurements from the PLL of the ESS are sampled at each Algorithm 1 Online Moving Horizon Estimation Input: L, τ , y k , p k , V , and W Output:M andD initialisation: 1: i ← 0 2: repeat 3: append frequency, ROCOF, and power measurements to the measurement matrix obtain initial guess for the solver 4: if (i = 0) then 5: get initial guess fromM andD of previous iteration 6: else 7: get initial guess from offline estimation 8: i ← 1 9: end if call the MHE solver 10: call MHE solver generated by ACADO Toolkit to estimateM andD 11: until excitation signal is stopped sampling instant during excitation. The MHE does not produce an estimate until the first L data points are collected. After the first L datapoints, any new subsequent data point is fed into the estimator while the oldest measurement is neglected. The MHE solver is then implemented to estimate the true states (the frequency and the ROCOF) and the parameters (the inertia and the damping constant). In the subsequent iterations, the latest estimate from the previous estimation is used as the initial guess of the state and the parameters for that particular time instant. This is done as ACADO is an iterative solver that converges faster with an accurate initial guess.

IV. CASE STUDY IN A LINEARIZED POWER SYSTEM MODEL A. SIMULATION SETUP
The proposed technique is first tested in the linearized power system model derived in Section II. The setup used for the simulation case study is shown in Fig. 3 and was implemented in MATLAB/Simulink. Table 1 specifies the parameters that were used for the simulation. An excitation signal p k is used to perturb the frequency of the system. This excitation signal is a pulse train with an amplitude of 0.1 p.u. and frequency of 1 Hz as shown in Fig. 4. This signal can be generated from the ESS where the proposed estimator will be implemented. If desired, an energy neutral square wave signal that charges and discharges the ESS by the same amount could also be used. However, the suitability of such an excitation signal will depend upon the topology of the ESS and will not impact the parameter estimates. The amplitude of the excitation signals also has an impact on the performance of the estimator. In this case, for the level of noise considered in the measurements, a pulse train of amplitude 0.1 p.u. (i.e., 10% of the total size of the system) was found to be sufficient. Reducing the amplitude below 10% of the microgrid size deteriorates It is assumed that a PLL is used to measure the signals ω k and ω k . For the size of the power systems considered, the measurements at any one location are sufficient as the observed inertia constant would remain relatively similar regardless of the measurement point. The measurements from the PLL are typically quite noisy, especially in low-voltage microgrid systems, due to factors such as unbalanced operations, harmonics, etc. Hence, to check the performance of the MHE under noisy measurement conditions, different levels of measurement noise are added to these signals. There has been no consensus on the measurement noise distribution of PLLs and PMUs in the literature. Different studies provide inconclusive evidence on the suitable distribution to be used for estimation approaches [34], [35]. Thus, the proposed technique is tested under different noise characteristics. First, the performance of the algorithm is validated when there is no measurement noise. Then, the performance of the technique is also tested for noise with different amplitudes and distribution characteristics. The signals ω k and ω k are the inputs to the MHE algorithm along with the excitation signal p k . Based on these local measurements, the MHE algorithm produces the frequency estimate ω k , the ROCOF estimates ω k , and estimates for the inertia constantM and damping constantD of the power system.

B. SETUP FOR THE MHE USING ACADO
The following setting was used when implementing the MHE solver based on the formulation described in (9). The matrices defining the system dynamics A d and B d in Constraint (9b) are obtained based on the system parameters   Table 1. In this particular study, it is assumed that the inertia constant M lies between 1 s and 10 s while the damping constant D can vary between 0.5% and 2%. The weighting matrix V is set based on the co-variance of the measurement noise. If g represents the co-variance of the measurement noise on the states, the weighting matrix is set to V = diag(0, g −0.5 , g −0.5 ). The weighting matrix W is set to 100 × g −0.5 as the measurement of the excitation signal is not expected to have significant actuation errors. This ensures that there is a relatively higher weight on the noise that is expected to be small in the cost function hence ensuring all the terms in the cost function contribute relatively equally and badly scaled problems are avoided. The final parameter that needs to be set is the backward time-horizon L. The selection of this parameter depends on the time-scale of the dynamics of interest (in this case the inertial response of the system) and on the amount of noise in the measured signals, which will be discussed in the subsequent sections.

V. RESULTS AND ANALYSIS
The performance of the MHE is first evaluated assuming that there is no noise in the measurement signals. Fig. 5 shows the performance of the MHE when estimating the inertia and the damping constant of the system. Initially, the inertia constant of the system is set to 4 s and the damping constant to 1.5%. It is assumed that these are unknown quantities. The inertia constant of the system is changed from 4 s to 2 s at a simulation time of 25 s while the damping constant remains constant throughout the simulation. The change in inertia could be a result of a generation loss in the microgrid or a rotational generator shutting down and being supplemented by a PV system without any inertial response. In this first case study, the inertia constant is changed by manually modifying the matrices A d and B d . However, in the second case study with the microgrid benchmark, the inertia constant is changed by connection/disconnection of generators in the system as discussed in Section VI. The initial guesses supplied to the online estimator were 2 s forM and 2% forD. The backward time-horizon L was set to 25 samples (0.25 s) in this case. Once initialized, the MHE waits for 25 samples before providing the first estimate. The MHE estimates converge to the true value within 2-3 s. As the inertia of the system does not typically change in such short timescales, this delay is acceptable for real-time estimation of inertia. The MHE was able to estimate both the inertia and damping constant accurately when there was no measurement noise involved. When the inertia constant of the system changes to 2 s at the simulation time of 25 s, the MHE estimates the new value fairly quickly within 1-2 s. During this time, there is a slight error in the estimate of the damping constant, but the MHE converges back to the true value within 2 s. Thus the MHE was able to estimate a drastic change in the inertia constant within 2 s.

A. ESTIMATION WITH GAUSSIAN NOISY MEASUREMENTS
In this case, the performance of the MHE is evaluated in the presence of measurement noise. Initially, it is assumed that the measurement noise has a Gaussian distribution. The performance was evaluated for different noise levels. Measurement noise that had a mean of 0 and different levels of co-variance was added to the measurements. The performance of the estimator was tested for co-variance of 1e-8, 1e-7, and 1e-6 corresponding to signalto-noise (SNR) ratios of 75, 65, and 55 dB, respectively. These are typical SNR metrics found for PLL measurements in the literature [36]. Fig. 5 shows the performance of the estimator when the measurement noises correspond to SNR of 75 dB and 55 dB. The estimates are compared against the case when there was no measurement noise in the system.  In both cases, the MHE was able to track the true inertia constant of the system although with small variations around the true value. However, the variations are small and within ±20% of the true value, providing sufficient accuracy to be used for adaptive control and protection schemes. A closer inspection shows that increasing the noise in the measurements results in constant small errors in the estimate especially in the case of the damping constant estimate.

B. EFFECT OF ESTIMATION TIME HORIZON ON PARAMETER ESTIMATES
Next, the effect of increasing the estimation time horizon of the MHE on the accuracy of parameter estimates is analyzed. For this, the measurement noise was assumed to have a co-variance of 1e-7 (SNR = 65 dB). Again, it was assumed the measurement noise had a Gaussian distribution. Using L < 25 did not provide good parameter estimates resulting in large deviations from the true values. Hence, Fig. 6 shows the inertia and damping constant estimates starting from L = 25. It can be observed that increasing the time horizon greatly improves the accuracy of the estimates. For instance, with L = 25 there are large variations in the inertia constant estimate around the true value and the damping constant estimate shows significant deviations compared to the true value. However, increasing the time horizon to 50 and 75 samples significantly improves both the inertia and the damping constant estimates. Increasing the time horizon, however, increases the computation cost of the estimator. Thus a balance needs to be found on the length of the selected prediction horizon based on the level of measurement noise and computation resource available.

C. ESTIMATION WITH NON-GAUSSIAN NOISY MEASUREMENTS
Contradicting conclusions have been made in the literature in regards to the distribution of the measurement noise from PLLs and PMUs. Furthermore, measurement error distribution can vary through time with changes/deterioration of potential transformers, current transformers, and change in communication channels [35]. This is where traditional estimation techniques such as Kalman filters and/or extended Kalman filters may fail. The deviation of the measurement noise distribution is characterized using the skewness and kurtosis metrics.
In this case, the performance of the proposed estimator is evaluated for different measurement noise distributions. The estimation is performed for a horizon length of L = 25. Fig. 7 shows the inertia and the damping constant estimates with skewness of +0.5 and kurtosis of 7.0. This provides skewed, non-Gaussian measurement noise and the proposed technique can be tested for cases when the measurement noise is non-Gaussian. The measurement noise distribution compared against a Gaussian distribution of the same mean and co-variance is shown in Fig. 7 (a). Figs. 7(b) and 7(c) show the inertia and damping constant estimation, respectively. The MHE can still estimate the inertia and the damping constant, although there is a slight increase in the variation in both estimates. A similar observation can be made when the measurement noise has negative skewness in Fig. 8. These results illustrate that the MHE performs well even though there are deviations from the usual assumption of a Gaussian noise distribution.

D. CHECK FOR PERSISTENCY OF EXCITATION
In [37] and [38], the authors have highlighted that absence of persistency of excitation can lead to incorrect estimation of parameters and in some cases exhibit periods of unstable (possibly oscillatory) behavior. To illustrate that the lack of persistency is not an issue in the proposed scheme, a long simulation was performed where both the inertia and the damping constant of the system do not change after the initial change at 25 s as described in previous tests. The simulation was performed for SNR = 65 dB and L = 25. The results shown in Fig. 9 show that the estimates do not diverge leading to the conclusion that persistency of excitation is not an issue in the proposed scheme. Thus, the excitation signal beyond the perturbation signal used is not required for the proposed scheme to work.

VI. VERIFICATION IN MICROGRID BENCHMARK A. BENCHMARK DESCRIPTION
To test the proposed inertia estimation technique, the remote microgrid test system from Cordova, Alaska, was modified as shown in Fig. 10. Only two substations are consideredthe ORCA 12.47 kV substation with four diesel generators (ORCA 3, ORCA 4, ORCA5, and ORCA 6) and the Humpback Creek (HBC) 12.47 kV substation, with three hydro generators (HBC 1, HBC 2, and HBC 3). Detailed parameters used to model each generator are listed in Appendix A1.
A 1 MW ESS unit was assumed to be installed at the ORCA substation. This ESS was used to inject an excitation signal similar to the one described in Fig. 4. The PLL of the ESS unit is then used to measure the change in frequency and ROCOF of the system based on the voltage measurements v abc from point of connection of the ESS. In reality, these measurements are noisy, so Gaussian noise of various noise levels was added with SNR of 75, 65, and 55 dB. These measurements were then fed into the MHE unit which then estimates the inertia constant of the system along with the filtered estimates of the change in frequency and ROCOF. The weights were set in a similar fashion to the case studies defined in Section IV. The backward time horizon L was set to 25 samples.

B. RESULTS AND ANALYSIS
For this benchmark, only the inertia constant of the system was estimated to simplify the MHE implementation. For typical power systems, the damping constant is usually in the range of 1-2%. First, the estimates are performed based on average frequency measurement calculated based on (1). The results with various noise levels in the measurements are shown in the form of violin plots in Fig. 11. The first three violin plots correspond to the case when only the generators at the ORCA substation are online, while the last three violin plots are for the case when all the generators at both the ORCA and the HBC substation are online. It can be observed that the inertia estimate is close to the inertia estimated based on (2), which is 2.34 when the noise is low at SNR of 75 dB. With increased noise (for SNR of 65 dB and 55 dB), the mean estimate deviated from this value and the spread of the estimates in general also increased. It is worth noting that, in the last three violin plots, when the generators at the HBC substation were added, the MHE was able to detect the change  in the equivalent inertia of the system. The inertia computed based on (2) in this case was 2.57.
Next, the same simulations were repeated using PLL-based frequency and ROCOF measurements instead of the average measurements. The details of the PLL are provided in Appendix A2. Similar observations as in the previous case can be made. The spread of the distribution of estimates increases with an increase in the noise level, and the increase in equivalent inertia with the addition of generators from the HBC substation was also detected by the MHE. However, it is interesting to note that the inertia constant is overestimated in general in all cases. This can be attributed to the dynamics of the low-pass filter (LPF) in the PLL. The LPF results in delays in frequency and ROCOF signals which the MHE interprets as the microgrid system having higher inertia thus leading to overestimation of the inertia constant.

VII. DISCUSSIONS AND CONCLUSION
A technique to estimate the inertia and damping constant of a microgrid system using MHE was introduced. The proposed technique was able to estimate the constants in real-time, which would allow an ESS operator to deploy fast-frequency support strategies. The estimates were obtained based on local measurements. A non-intrusive excitation signal was used to perturb the frequency of the system to perform the estimates. The MHE showed good performance under typical PLL measurement noise amplitudes and distributions. Since there has not been a clear consensus on the measurement noise distribution for frequency measurements, the ability of the estimator to work on different noise distributions is a major contribution. The proposed technique was also tested in a realistic benchmark from Cordova, Alaska. Leveraging an approximate model of the microgrid system, the MHE was able to closely estimate the inertia constant of the system. From a practical standpoint, it is important to note that in some cases introducing active power pulse trains as proposed in this paper may lead to inter-area oscillations especially in multi-area power systems. However, inter-area oscillations are extremely rare especially in microgrids. These aspects need to be carefully considered when implementing the proposed technique. Instead of pulse trains, Pseudo-Random Binary Signals (PRBS), multi-sine wave signals, etc. could be utilized to prevent such issues. However, from a theoretical perspective, the active power pulse train gives the best performance [39]. The technique also assumes a certain amount of the ESS energy is reserved for inertia estimation, however, as the estimation is performed within a short-time interval, the energy consumption is expected to be low. Furthermore, the energy spent in one estimation phase (with discharging pulse trains) may be recovered in the consequent estimation phase (with discharging pulse trains) through proper coordination. Another key assumption that the paper makes is that the microgrid is still dominated by rotational generators. The applicability of the underlying predictive model used by the MHE needs to be further evaluated in microgrids with large amounts of converter-based generation. We hypothesize that adjustments to the predictive model and re-tuning of the weighting parameters might be required in case the algorithm does not perform well in this specific situation. Our future work will investigate how this framework can be adapted for larger interconnected (multi-area) power systems and power system with large amounts of converter-based generation.

APPENDIX A1
The parameters used to model the generators in the Cordova, Alaska microgrid benchmark are provided in Table 2.

APPENDIX A2
The standard three phase PLL available in the MATLAB/Simulink library was used to measure the frequency and ROCOF in the Cordova, Alaska microgrid benchmark. The parameters of the PLL are provided below: 1) Regulator gains (K p , K i , K d ) = (180, 3200, 1) 2) Time constant for derivative action (T d ) = 0.0001 s 3) Filter cut-off frequency (f c ) = 10 Hz