Application of Kalman filter based estimation techniques to electric power systems

º Ph.D. Dissertation Program: Sistemas de Energía Eléctrica Department of Electrical Engineering Escuela Técnica Superior de Ingeniería Universidad de Sevilla Seville, 2020 Application of Kalman filter – based estimation techniques to electric power systems Author: Miguel Ángel González Cagigal Advisors: José Antonio Rosendo Macías Antonio Gómez Expósito iii iv Ph. D. Dissertation Program: Sistemas de Energía Eléctrica Application of Kalman filter – based estimation techniques to electric power systems Author: Miguel Ángel González Cagigal Advisors: José Antonio Rosendo Macías Antonio Gómez Expósito Department of Electrical Engineering Escuela Técnica Superior de Ingeniería Universidad de Sevilla Seville, 2020 v vi Ph. D. dissertation: Application of Kalman filter – based estimation techniques to electric power systems Author: Miguel Ángel González Cagigal Advisors: José Antonio Rosendo Macías y Antonio Gómez Expósito The thesis committee, composed by the following members: Agree to the calification: Seville, 2020 The commitee Secretary. vii viii A mi familia A mis maestros ix xvi LIST OF FIGURES Figure 1-1. Variation of the resistivity with the temperature for different materials . 2 Figure 2-1. Diagram of the recursive estimation process of the KF ............................. 10 Figure 3-1. Diagram of the two-stage implementation for the UKF estimator .......... 20 Figure 3-2. Voltage magnitude and angle at the generator bus .................................... 23 Figure 3-3. Considered measurements for the estimation process............................... 24 Figure 3-4. Estimation of the modified parameters in the first stage ........................... 26 Figure 3-5. Estimation of the modified parameters in the second stage ..................... 29 Figure 4-1. Wind speed variation in the base case ........................................................... 40 Figure 4-2. Modified inertia estimation. Base case ........................................................... 41 Figure 4-3. Estimation of the modified parameters from the pitch angle controller and the VSC. Base case ........................................................................................................... 41 Figure 4-4. Representation of the wind gust modelled as a Mexican hat wavelet ....... 44 Figure 4-5. Estimation of the modified shaft inertia. Disturbance 1 ............................. 44 Figure 4-6. Estimation of the modified parameters. Disturbance 1 .............................. 45 Figure 4-7. Voltage dip considerd in the second disturbance ....................................... 47 Figure 4-8. Example of the estimation results for the modified parameters. Disturbance 2 ............................................................................................................................ 47 Figure 4-9. Example of the estimation results for the modified parameters. Disturbance 3 ............................................................................................................................ 48 Figure 4-10. Performance comparison in the estimation of shaft inertia by CKF, UKF and EKF ............................................................................................................................ 58 xvii Figure 5-1. Block diagram of the proposed methodology for the CPI problem ....... 63 xviii LIST OF TABLES Table 3-1. Distribution of inputs and measurements for the synchronous generator ...................................................................................................................................................... 21 Table 3-2. Relative errors for the parameters of the synchronous generator ............. 27 Table 3-3. Relative error in the estimation of the parameters of the machine regulators ................................................................................................................................... 29 Table 3-4. Estimated values of the parameters and relative errors for 𝑅𝑤=10−6 ... 32 Table 3-5. Estimated values of the parameters and relative errors for 𝑅𝑤=10−4 ... 33 Table 3-6. Estimated values of the parameters and relative errors for 𝑅𝑤=10−3 ... 34 Table 3-7. Mean and máximum relative errors for different values of 𝑅𝑤 ................. 35 Table 4-1. Distribution of inputs and measurements for the wind turbine ................ 38 Table 4-2. Relative errors in the parameter estimation. Base case ................................ 43 Table 4-3. Estimation errors for the parameters with a 3% measurement error ....... 50 Table 4-4. Estimation errors for the parameters with a 5% measurement error ....... 51 Table 4-5. Estimation errors for the parameters with a 7% measurement error ....... 52 Table 4-6. Maximum relative error for increasing levels of measurement error....... 52 Table 4-7. Estimation errors for the parameters with a 2% model error ..................... 54 Table 4-8. Estimation errors for the parameters with a 5% model error ..................... 55 Table 4-9. Estimation errors for the parameters with a 10% model error ................... 56 Table 4-10. Estimation errors for the parameters with a 15% model error ................. 57 Table 4-11. Maximum relative error for increasing levels of model error .................. 57 Table 5-1. Estimation results for the UKF. Scenario 1 ...................................................... 64 Table 5-2. Estimation results for the CKF. Scenario 1 ...................................................... 65 xix Table 5-3.Estimation results for the EnKF. Scenario 1 ..................................................... 65 Table 5-4. Estimation results for the UKF. Scenario 2 ..................................................... 66 Table 5-5. Estimation results for the CKF. Scenario 2 ...................................................... 67 Table 5-6. Estimation results for the EnKF. Scenario 2 .................................................... 67 Table 5-7. Estimation results for the UKF. Scenario 3 ..................................................... 68 Table 5-8. Estimation results for the CKF. Scenario 3 ...................................................... 69 Table 5-9. Estimation results for the EnKF. Scenario 3 .................................................... 69 Table 5-10. Estimation results for the UKF. Scenario 4 ................................................... 70 Table 5-11. Estimation results for the CKF. Scenario 4 ................................................... 71 Table 5-12. Estimation results for the EnKF. Scenario 4 ................................................. 71 Table 5-13. Rate of correct assignments for the PCA, the LASSO-based and the proposed KF-based methods ................................................................................................ 73 1 1 INTRODUCTION his introductory chapter is aimed to establish the context embracing this doctoral thesis, which is split into three parts subsequently described. Along with the technical framework, a thorough revision of the state of the art is presented regarding the considered dynamic state estimator applications, allowing to assess the contributions of this study. 1.1. Parameter estimation in electric power plants. Context and motivation Electric power systems (EPS) are currently experiencing a remarkable period of change, given the inclusion of new technologies based on power electronics, and the consolidation of the renewable sources of energy, which have modified the existing paradigm in the conventional energy production system. The latter aspects make more difficult the operation and control in electrical networks, increasing the vulnerability of the system towards the possible disturbances that might appear. Moreover, the mentioned tasks of operation and control require a deep knowledge of the dynamics associated to each of the elements comprising the grid, so that the response of this components can be predicted, and precautionary measures can be taken. The transient behavior is determined by a set of differential-algebraic equations, establishing the dynamic models for the different units. These equations involve a collection of characteristic parameters, for which an adequate knowledge is essential in order to correctly identify the state of the system and take the T 2 Application of Kalman filter – based estimation techniques to electric power systems corresponding actions. The values of these parameters are sometimes provided by the manufacturer of the equipment, while in several cases this information is not available. Regardless the situation, a usual practice is assuming the parameters to be invariable with respect to the operating conditions. A particular example might be the widely studied dependence of some parameters, as the electrical resistivity, with the temperature (Figure 1-1, [1]). Figure 1-1. Variation of the resistivity with the temperature for different materials Besides the temperature, other factors might affect significantly the values of the parameters involved in the dynamic models of the components of an EPS, such as: - The operating point of the elements, determined by the voltage and the exchanged power in a certain moment. Introduction 3 - The degradation caused by unusual working conditions or aging. The variation of the parameters due to the latter factors is not studied as deeply as that caused by the temperature, given its potential complexity, and resulting in the previously mentioned assumption of neglecting those changes. Finally, the absence of information provided by the manufacturer may also lead to a lack of knowledge on their equipment parameters. In light of the above mentioned aspects, it is clear the usefulness of a method that allows to establish a reliable estimation of the parameters involved in the dynamics of the components of an EPS, given the external factors previously highlighted. This is the main goal of parts 1 and 2 of this thesis, which will approach the described issues with a parameter estimation technique for two particular elements of the system: the synchronous generator and the variable speed-wind turbine. 1.2. Electrical phase assignment in distribution grids Focusing on distribution networks, operation and control tasks are essential in order to guarantee the good quality of the service provided to the consumers. In this context, grid operators require an unambiguous knowledge of the load distribution in each of the three phases of the system, so that this information allows to state a correct load equilibrium in those phases. The importance of this problem is even more remarkable when considering renewable sources of energy connected to the distribution grids, since the phase identification helps to establish a better generation-consumption balance in each phase. In this regard, despite the efforts undertaken by distribution companies, they frequently lack enough information about the phase connection of their singlephase customers, owing, for instance, to network reconfiguration after faults, phase switching derived from improper maintenance, or inaccurate recording of the true load-to-phase connectivity. In these circumstances, a method must be developed to estimate as accurately as possible the actual phase to which a customer is connected in LV feeders, which is known as the customer-phase identification (CPI) problem. The CPI problem has been approached in several ways by previous works. In [2], a signal processing perspective is applied to voltage observations, which are also 4 Application of Kalman filter – based estimation techniques to electric power systems used in [3], for a correlation-based methodology. Two techniques are highlighted in this section: a method proposed in [4], using Least Absolute Shrinkage and Selection Operator (LASSO) and a novel approach in [5] for phase idenfication using graph theory and principal component analysis (PCA). The latter methodologies will be included in a comparative study with the estimation technique proposed in this thesis for the CPI problem. 1.3. Dynamic state estimators. Kalman filter A dynamic state estimator (DSE) is aimed to identify the state of the system which is unobservable through direct measurement. This state, defining the system under study, is considered to be time-variant, unlike the traditional static estimators. An example of DSE might be the Luenberger observer, [6], based on the definition of an equivalent system for which it is possible the direct observation of the state variables whose knowledge is desired. However, this work is focused on other DSE, the so-called Kalman filter (KF). Chapter 2 includes a brief historical context of this estimator and the filter equations, while the particular applications considered will be described in chapters 3-5. The KF has been widely used in several works related to EPS. A major división can be made regarding those studies related to large-scale estimation (the reader may refer to [7–9]), and those works focused on the characterization of single elements of the system, being the latter field of investigation closer to the objectives of this document. The original formulation of the KF, [10], assumes a linear dynamic of the system, being this condition quite uncommon in the equations of the elements in EPS. To approach this issue, different formulations have arised to deal with non-linear dynamics. Among them, the extended KF (EKF), is based on the calculation of the jacobian in order to linearize the state equations. This scheme has been widely used in several works related to EPS, as that in [11], where this estimator is used for state estimation with unkown inputs. Regarding synchronous machines, [12] makes use of the EKF to estimate the initial rotor position of permanent magnets synchronous machines (PMSMs), while the flux losses are included in the estimation process of [13], using vector control, and in [14], where demagnetization problems are Introduction 5 considered. A study in [15] concludes that the EKF and another non-linear formulation, the unscented KF (UKF), are the most computationally efficient when the synchronous generator dynamic is considered. The latter formulation has been used in [16] for state estimation in single machines and in [17–19] for multimachine systems. In the particular case of synchronous generators, the joint estimation of state and parameters has been approached in several works using the UKF (see [20] and [21]). However, the application of those estimation techniques in real systems is hindered by the use of internal measurements. To avoid these difficulties, different strategies have been developed, as that in [22], using a set of measurements obtained from phasor measurement units (PMUs). The method proposed in [23] broadens the problem by including in the model the parameters of the speed governor (SG) and the automatic voltage regulator (AVR), studying the effect of a variable error (up to ±30%) in the initialization of those parameters. The power system stabilizer (PSS) is introduced in the KF model in [24], although the parameters of this regulator are taken as known inputs of the system. This assumption might not be accurate enough given the variability of these parameters and the possible lack of information. With the above-mentioned considerations, and to the author’s knowledge, the publication included in annex II.1 is the first attempt to consider a fully regulated generation set, including the synchronous machine itself along with the whole set of regulators (SG, AVR and PSS). Both state and parameters from these elements are jointly estimated using UKF and measurements obtained at the genartor terminal bus, which could be provided by a PMU in real applications. Furthermore, a recent formulation of the KF, the so-called cubature KF (CKF) has proven a correct performance when it is used for the dynamic state estimation in synchronous generators, [25]. To see the theoretical aspects of this scheme, the reader may refer to [26], where some limitations of other KF formulations, such as UKF, which are not suffered by the CKF, are highlighted. Joint state and parameter estimation using CKF is studied in [27], applied to a vehicle model, and in [28], to permanent magnet synchronous motors. In the second part of this thesis, the CKF is considered for the state and parameter 12 Application of Kalman filter – based estimation techniques to electric power systems The continuous-time counterpart of the original KF formulation, the so-called Kalman-Bucy filter [30], is used for continuous systems represented by the following state and measurement equations: 𝑥󰇗(𝑡)= 𝐴·𝑥(𝑡)+𝐵·𝑢(𝑡)+𝑤(𝑡) (7) 𝑧(𝑡)= 𝐶·𝑥(𝑡)+𝑣(𝑡) (8) being 𝑢(𝑡) the system input. However, since the measurements are collected at discrete time instants in real applications of the estimator, a continuous-time, discrete-measurement framework is defined with the following equations: 𝑥󰇗(𝑡)= 𝐴·𝑥(𝑡)+𝐵·𝑢(𝑡)+𝑤(𝑡) (9) 𝑧𝑘= 𝐶𝑘·𝑥𝑘+𝑣𝑘 (10) Where a discretization of the state equation (9) is used to calculate the a priori estimation of the state vector at instant 𝑘. Finally, in most applications related to EPS, and specifically those presented in this thesis, the KF is applied to discrete-time, discrete-measurement nonlinear systems. In the discrete-time framework, the associated equations may be expressed as, 𝑥𝑘=𝑓(𝑥𝑘−1,𝑢𝑘−1)+𝑤𝑘 (11) 𝑧𝑘= 𝑔(𝑥𝑘,𝑢𝑘)+𝑣𝑘 (12) where 𝑓(·) and 𝑔(·) are nonlinear functions. Kalman filter context 13 In order to approach the nonlinearities in the system equations, several formulations of the KF have been developed, from which three of them are presented in the next sections, with their corresponding equations. 2.3. Unscented Kalman filter At instant 𝑘, a cloud of 2𝐿+1 vectors, the so-called σ-points, is obtained from the previously estimated state vector, 𝑥𝑘−1 (dimension 𝐿), and the covariance matrix of the state estimation error, 𝑃𝑘−1, as follows [31]: 𝑋𝑘−1 0=𝑥𝑘−1 (13) 𝑋𝑘−1 𝑖=𝑥𝑘−1 +[√(𝐿+𝜆)𝑃𝑘−1]𝑖 𝑖 =1,…,𝐿 𝑋𝑘−1 𝑖=𝑥𝑘−1 +[√(𝐿+𝜆)𝑃𝑘−1]𝑖 𝑖 =1,…,𝐿 [√(𝐿+𝜆)𝑃𝑘−1]𝑖 being the 𝑖𝑡ℎ column of the corresponding matrix, and λ a scaling factor calculated as 𝜆 = 𝛼2(𝐿+𝜅)−𝐿, where α and κ are two filter parameters to be tuned. These σ-points are evaluated through equation (11), yielding 2𝐿+1 vectors, 𝑋𝑘 𝑖−, from which the a priori estimations 𝑥𝑘 − and 𝑃𝑘 − are obtained: 𝑥𝑘 −= ∑𝑊𝑚𝑖𝑋𝑘 𝑖− 2𝐿 𝑖=0 (14) 𝑃𝑘 −=∑𝑊𝑐𝑖(𝑋𝑘 𝑖− 2𝐿 𝑖=0 −𝑋 𝑘 −)(𝑋𝑘 𝑖− −𝑋 𝑘 −)𝑇+𝑄𝑘 (15) where the weighting vectors 𝑊𝑚 and 𝑊𝑐 are calculated from 14 Application of Kalman filter – based estimation techniques to electric power systems 𝑊𝑚0 =𝜆 𝐿+𝜆 (16) 𝑊𝑐0 =𝜆 𝐿+𝜆+1−𝛼2+𝛽 𝑊𝑚𝑖 =1 2(𝐿+𝜆) 𝑖 =1,…,2𝐿 β being another tunable parameter. On the basis of the a priori information, the correction stage starts with the calculation of a new cloud of vectors, 𝑋𝑘 𝑖, which are evaluated with the measurement function 𝑔(·) in equation (12), and weighted with the vectors 𝑊𝑚, yielding 𝛾𝑘 𝑖− = 𝑔(𝑋𝑘 𝑖, 𝑢𝑘) i= 0,…,2L (17) 𝑧𝑘 −=∑𝑊𝑚𝑖𝛾𝑘 𝑖− 2𝐿 𝑖=0 (18) Then, the covariance matrix of the measurement estimation error, 𝑃𝑧𝑘 −, and the cross-covariance matrix of state and measurements, 𝑃𝑥𝑧𝑘 −, are obtained using the vector 𝑊𝑐 as follows: 𝑃𝑧𝑘 −= ∑𝑊𝑐𝑖(𝛾𝑘 𝑖− 2𝐿 𝑖=0 −𝑧𝑘 −)(𝛾𝑘 𝑖− −𝑧𝑘 −)𝑇+𝑅𝑘 (19) 𝑃𝑥𝑧𝑘 −=∑𝑊𝑐𝑖(𝑋𝑘 𝑖− 2𝐿 𝑖=0 −𝑥𝑘 −)(𝛾𝑘 𝑖− −𝑧𝑘 −)𝑇 (20) Kalman filter context 15 The correction stage concludes with the a posteriori predictions, 𝑥𝑘=𝑥𝑘 −+𝐾𝑘(𝑧𝑘−𝑧𝑘 −)𝑇 (21) 𝑃𝑘=𝑃𝑘 −−𝐾𝑘𝑃𝑧𝑘 −𝐾𝑘 𝑇 (22) which are based on the 𝑎 𝑝𝑟𝑖𝑜𝑟𝑖 predictions at instant 𝑘 and the so-called Kalman gain, 𝐾𝑘, calculated from 𝐾𝑘= 𝑃𝑥𝑧𝑘 −(𝑃𝑧𝑘 −)−1 (23) 2.4. Cubature Kalman filter The CKF formulation, more recent than the UKF, was developed in 2010 as a technique to estimate the positioning in the navigation field. This estimator uses 2𝐿 cubature points, which are calculated from the previous information, 𝑥𝑘−1 and 𝑃𝑘−1, through the following equations, [32]: 𝑆𝑘 −𝑆𝑘 −𝑇 =𝑃𝑘 − (24) 𝑥𝑘 𝑖− = 𝑆𝑘 −ξ𝑖 √𝐿+𝑥𝑘 − 𝑖 =1,...,2𝐿 (25) where 𝑆 is a positive-definite square root matrix of 𝑃, and ξ𝑖 is the 𝑖𝑡ℎ cubature node, obtained as the intersection of the unit sphere and the ℝL axis. The state function 𝑓(·) in equation (11) is evaluated for the set of cubature points, yielding 2𝐿 vectors 𝑋𝑘 𝑖−, from which the a priori estimation is computed, 16 Application of Kalman filter – based estimation techniques to electric power systems 𝑥𝑘 −=1 2𝐿∑𝑋𝑘 𝑖− 2𝐿 𝑖=1 (26) 𝑃𝑘 −=1 2𝐿∑𝑥𝑘 𝑖−𝑥𝑘 𝑖−𝑇 2𝐿 𝑖=1 −𝑥𝑘 −𝑥𝑘 −𝑇 +𝑄𝑘 (27) For the correction stage, the covariance matrix is factorized in order to calculate both the matrix 𝑆𝑘 − 𝑆𝑘 −𝑆𝑘 −𝑇 =𝑃𝑘 − (28) and a new set of 2𝐿 cubature points, 𝑋𝑘 𝑖−, at which function 𝑔(·) is evaluated to obtain 𝛾𝑘 𝑖−. Then, the measurement estimation, 𝑧𝑘 −, its covariance matrix, 𝑃𝑧𝑘 −, and the crosscovariance matrix of state and measurement, 𝑃𝑥𝑧𝑘 −, are calculated as follows: 𝑧𝑘 −=1 2𝐿∑𝛾𝑘 𝑖− 2𝐿 𝑖=1 (29) 𝑃𝑧𝑘 −=1 2𝐿∑𝛾𝑘 𝑖−𝛾𝑘 𝑖−𝑇 2𝐿 𝑖=1 −𝑧𝑘 −𝑧𝑘 −𝑇 +𝑅𝑘 (30) 𝑃𝑥𝑧𝑘 −=1 2𝐿∑𝑥𝑘 𝑖−𝛾𝑘 𝑖−𝑇 2𝐿 𝑖=1 −𝑥𝑘 −𝑧𝑘 −𝑇 (31) Kalman filter context 17 The a posteriori predictions of the state vector, 𝑥𝑘, and the covariance 𝑃𝑘 are calculated with the Kalman gain using the same equations (21) and (22) as in the UKF. 2.5. Ensemble Kalman filter Finally, the EnKF formulation, [33], is a Monte Carlo approximation of the original KF, which has proven accurate enough in high-dimensional state-space problems. The ensemble is represented by a 𝐿x𝑁 matrix, 𝑁 being the number of samples considered. This ensemble is first propagated through the state and measurement functions [34], 𝑥𝑘 𝑖− = 𝑓(𝑥𝑘−1 𝑖,𝑢𝑘−1) (32) 𝑧𝑘 𝑖− = 𝑔(𝑥𝑘 𝑖,𝑢𝑘) 𝑖 =1,…,𝑁 and then the mean values are calculated: 𝑥𝑘 −     =1 𝑁∑𝑥𝑘 𝑖− 𝑁 𝑖=1 (33) 𝑧𝑘 −    =1 𝑁∑𝑧𝑘 𝑖− 𝑁 𝑖=1 (34) 18 Application of Kalman filter – based estimation techniques to electric power systems The EnKF correction stage is based on the calculation of the intermediate matrices 𝑃𝑘    𝐻𝑘 𝑇=1 𝑁∑(𝑥𝑘 𝑖− −𝑥𝑘 −     )(𝑧𝑘 𝑖− −𝑧𝑘 −    )𝑇 𝑁 𝑖=1 (35) 𝐻𝑘𝑃𝑘    𝐻𝑘 𝑇=1 𝑁∑(𝑧𝑘 𝑖− −𝑧𝑘 −    )(𝑧𝑘 𝑖− −𝑧𝑘 −    )𝑇 𝑁 𝑖=1 (36) which allow the Kalman gain and the updated values of each sample in the ensemble to be obtained: 𝐾𝑘= 𝑃𝑘    𝐻𝑘 𝑇(𝐻𝑘𝑃𝑘    𝐻𝑘 𝑇+𝑅)−1 (37) 𝑥𝑘 𝑖=𝑥𝑘 𝑖− +𝐾𝑘(𝑧𝑘−𝑧𝑘 𝑖−) 𝑖 =1,…,𝑁 (38) Finally, the corrected covariance matrix, 𝑃𝑘, is calculated as follows: 𝑥𝑘    =1 𝑁∑𝑥𝑘 𝑖 𝑁 𝑖=1 (39) 𝑃𝑘=1 𝑁∑(𝑥𝑘 𝑖− −𝑥𝑘    )(𝑥𝑘 𝑖− −𝑥𝑘    )𝑇 𝑁 𝑖=1 (40) 19 3. APPLICATION 1: FULLY REGULATED SYNCHRONOUS GENERATOR his section presents the application of a KF-based technique to the estimation of the parameters involved in the dynamic equations of a synchronous generator and its regulators (SG, AVR and PSS). First, the proposed implementation of the UKF formulation is described, as well as the tuning of the estimator. The considered case studies are presented subsequently, showing the corresponding estimation results. 3.1. Kalman filter implementation The UKF scheme is used in this part of the thesis for the joint estimation of the state and parameters of a fully regulated synchronous generator. The algorithm of this formulation has been coded using MATLAB. First attempts with a vanilla implementation of the joint estimator, directly using the whole set of parameters in the system equations, led to unacceptable convergence behavior. Therefore, as proposed in [35] for a simpler case with singleaxis generator and only frequency regulation, a new modified parameter vector has been considered using different combinations with their original values. Finally, all the convergence problems were solved using a novel two-stage formulation, as illustrated in Figure 3-1, so that the modified parameters from the synchronous machine are identified in the first stage. These results are incorporated T 20 Application of Kalman filter – based estimation techniques to electric power systems to the full model in the second stage for the estimation of the parameters from the machine regulators. Figure 3-1. Diagram of the two-stage implementation for the UKF estimator The UKF tuning is determined by the following points: - Values of the UKF tuning parameters. - Definition of the system inputs and measurements. - Initial estimation. Application 1: fully regulated synchronous generator 21 - Process and measurement noise covariance matrices. First, the UKF algorithm involves a set of tuning parameters, α,β and κ. As proposed in [36], where the influence of these values in the performance of the estimator is analyzed and accordingly the values α=10−4 y β = 2 have been taken in this work. Regarding the optimal tuning for κ, it is obtained using the following expression: κ+𝐿 = 3 being L the size of the state vector. With respect to the measurements taken from the system, five magnitudes can be easily measured at the generator terminal bus, namely: - Voltage magnitude, 𝑉 - Voltage angle, θv - Angular speed, ω - Active power exchanged, 𝑝𝑒 - Current magnitude, 𝐼 Those magnitudes are divided, as proposed in [37], into inputs, 𝑢(𝑡), and measurements, 𝑧(𝑡), as it is summarized in Table 3-1. Table 3-1. Distribution of inputs and measurements for the synchronous generator Input vector, 𝑢(𝑡) Measurement vector, 𝑧(𝑡) 𝑉 𝐼 θ𝑉 𝑝𝑒 ω 28 Application of Kalman filter – based estimation techniques to electric power systems 3.2.2. Second estimation stage In the second stage of the estimation process, the modified parameters of the three machine regulators considered (SG, AVR and PSS) are included in the augmented state vector of the KF, to be estimated. In this case, the information obtained from the first stage is used, given that the initial values of the parameters of the synchronous machine are taken as their previous estimation (Table 3-2), with a high confidence in the components of 𝑥0, meaning low values for the corresponding elements of the matrix 𝑃0. The system inputs and the measurements considered for the UKF implementation remain the same as those used in the first stage, including the considered variations in the terminal voltage magnitude. For the representation of the modified parameters, the same considerations as in the first stage have been taken, showing the estimated value jointly with a ±3𝜎𝑖 deviation. Figure 3-5 illustrates the evolution of the modified parameters of the considered regulators. The graphics corresponding to the estimation of the machine parameters are not included, since they suffer no significant variation with respect to the initial estimation and, therefore, no additional information is provided. In this figure it can be observed, as with the generator modified parameters, that those of the three regulators reach stable values close to the simulated ones, with a reduced covariance, giving evidence of the good performance of the proposed estimation methodology. Finally, Table 3-3 summarizes the estimated values of the original parameters of the regulators, with their corresponding relative errors. The fact that the maximum error is only slightly higher than 3% also shows that the accuracy of the estimation technique is more than acceptable. Application 1: fully regulated synchronous generator 29 Table 3-3. Relative error in the estimation of the parameters of the machine regulators Parámetro Valor estimado Error relativo (%) 𝑅 0.1032 3.20 𝑇𝑟 0.1006 0.64 𝐾𝑖 49.75 0.498 𝑇𝑎𝑣𝑟 0.988 1.186 𝑇𝑒 0.2002 0.101 𝐾0 2.479 0.84 𝐾𝑤 30.78 2.596 𝑇𝑤 10.135 1.35 Figure 3-5. Estimation of the modified parameters in the second stage 30 Application of Kalman filter – based estimation techniques to electric power systems Application 1: fully regulated synchronous generator 31 3.2.3. Variations of the base case To asses the robustness of the UKF-based estimator proposed in this thesis, three additional scenarios have been considered with different values in the standard deviation of the random walk, 𝑅𝑤, defining the magnitude and angle of the external voltage. The comparison is based on the estimated values of the parameters and their corresponding relative errors. Tables 3-4 to 3-6 include these estimated values for 𝑅𝑤=10−6, 10−4 and 10−3 respectively. To make the analysis easier, Table 3-7 summarizes the mean and maximum relative errors in each scenario, including the base case. 32 Application of Kalman filter – based estimation techniques to electric power systems Table 3-4. Estimated values of the parameters and relative errors for 𝑅𝑤=10−6 Parameter Estimated value Relative error (%) 𝑥𝑑 0.985 1.50 𝑥𝑑 ′ 0.298 0.67 𝑥𝑞 0.656 0.92 𝑥𝑞 ′ 0.562 2.18 𝐻 6.345 2.38 𝐷 1.994 0.03 𝑇′𝑑0 0.506 1.20 𝑇′𝑞0 5.116 2.32 𝑅 0.1031 3.05 𝑇𝑟 0.1006 0.60 𝐾𝑖 49.75 0.50 𝑇𝑎𝑣𝑟 0.987 1.31 𝑇𝑒 0.198 1.02 𝐾0 2.479 0.84 𝐾𝑤 30.65 2.17 𝑇𝑤 10.115 1.15 Application 1: fully regulated synchronous generator 33 Table 3-5. Estimated values of the parameters and relative errors for 𝑅𝑤=10−4 Parameter Estimated value Relative error (%) 𝑥𝑑 0.982 1.78 𝑥𝑑 ′ 0.296 1.12 𝑥𝑞 0.661 1.64 𝑥𝑞 ′ 0.563 2.35 𝐻 6.340 2.51 𝐷 1.983 0.84 𝑇′𝑑0 0.510 2.03 𝑇′𝑞0 5.137 2.74 𝑅 0.1041 4.15 𝑇𝑟 0.1011 1.11 𝐾𝑖 49.519 0.97 𝑇𝑎𝑣𝑟 0.984 1.59 𝑇𝑒 0.203 1.33 𝐾0 2.473 1.09 𝐾𝑤 31.203 4.01 𝑇𝑤 10.309 3.10 34 Application of Kalman filter – based estimation techniques to electric power systems Table 3-6. Estimated values of the parameters and relative errors for 𝑅𝑤=10−3 Parameter Estimated value Relative error (%) 𝑥𝑑 0.980 2.04 𝑥𝑑 ′ 0.296 1.33 𝑥𝑞 0.661 1.74 𝑥𝑞 ′ 0.567 2.89 𝐻 6.308 3.04 𝐷 2.023 1.14 𝑇′𝑑0 0.511 2.20 𝑇′𝑞0 5.136 2.72 𝑅 0.1048 4.84 𝑇𝑟 0.1013 1.27 𝐾𝑖 49.363 1.29 𝑇𝑎𝑣𝑟 1.024 2.35 𝑇𝑒 0.204 2.06 𝐾0 2.453 1.88 𝐾𝑤 31.503 5.01 𝑇𝑤 10.407 4.07 Application 1: fully regulated synchronous generator 35 Table 3-7. Mean and máximum relative errors for different values of R_w 𝑅𝑤 10−6 10−5 10−4 10−3 |𝐸𝑟 𝑚𝑎𝑥| (%) 3.05 3.2 4.15 5.01 |𝐸𝑟 𝑚𝑒𝑑| (%) 1.36 1.27 2.02 2.49 As it is observed, the results are similar in the first two cases, with 𝑅𝑤=10−6 and 10−5, becoming the errors higher as the value of 𝑅𝑤 increases. Nevertheless, in all the considered scenarios, the mean and maximum errors remain within an acceptable range, showing the accuracy and robustness of the two-stage estimation technique presented in this thesis using UKF. 36 Application of Kalman filter – based estimation techniques to electric power systems 37 4. APPLICATION 2: VARIABLE-SPEED WIND TURBINE his section presents the application of a KF-based technique to the estimation of the parameters involved in the dynamic equations of a variable-speed PMSM wind turbine connected to the grid through a back-to-back pair of VSCs. First, the proposed implementation of the CKF formulation is described, as well as the tuning of the estimator. The considered case studies are presented subsequently, showing the corresponding estimation results. 4.1. Kalman filter implementation In this part of the thesis, the CKF scheme is used for a joint estimation of the state and parameters involved in the dynamic equations of the wind turbine, the PMSM, the VSC and the pitch angle controller, [37]. The algorithm of the CKF, presented in chapter 2, has also been coded using MATLAB. As described in chapter 3 for the fully regulated synchronous generator, the proposed methodology is based on the modification of the original parameters in the equations. Regarding the CKF tuning, it consists of the same points presented for the UKF, except for the use of tuning parameters, which are not included in the CKF algorithm. The following signals have been taken from a simulated model, representing the measurements obtained from a real system: T 44 Application of Kalman filter – based estimation techniques to electric power systems Figure 4-6 shows the evolution experienced by the rest of the modified parameters in the model. The robustness of the proposed methodology can be concluded given the minimum deviations presented by all the parameters from the VSC and the pitch angle controller. Figure 4-4. Representation of the wind gust modelled as a Mexican hat wavelet Figure 4-5. Estimation of the modified shaft inertia. Disturbance 1 Application 2: variable-speed wind turbine 45 Figure 4-6. Estimation of the modified parameters. Disturbance 1 46 Application of Kalman filter – based estimation techniques to electric power systems 4.2.3. Disturbance 2: Voltage dip The second disturbance considered is an abrupt descent in the voltage magnitude at the connection point of the generation set with the external network, which emulates a possible voltage dip caused by a short-circuit close to this bus of the electrical grid. Specifically, a 70% drop in the voltage is simulated, with a duration of 1 s. The standard deviation of the random walk of this voltage remains at the same value as in the base case. The evolution of the magnitude and angle in this case is represented in Figure 4-7. The wind speed profile in this scenario is that used in the base case. With the simulated signals used as inputs as measurements in the CKF estimator, the process was carried out using the same initial estimation as in the first disturbance. As an example, Figure 4-8 includes the evolution of a modified time constant from the VSC. Since the rest of the parameters show no significant variations they are omitted in this document. Application 2: variable-speed wind turbine 47 Figure 4-7. Voltage dip considerd in the second disturbance Figure 4-8. Example of the estimation results for the modified parameters. Disturbance 2 48 Application of Kalman filter – based estimation techniques to electric power systems Once again, the transient behavior due to the disturbance is soon vanished and the parameter stabilizes in its initial estimation. The nature of this deviation remarks the robustness of the proposed technique. 4.2.4. Disturbance 3: Topological change In the last disturbance considered in this part of the thesis, an abrupt and permanent change is simulated in the impedance modeling the tie-line connecting the generation set and the external network. This might represent a sudden change in the external grid configuration. The wind profile and the grid voltage remains as in the base case. For this scenario, the evolution of the modified time constant of the pitch angle controller, 1/𝑇𝑝, is represented in Figure 4-9 as an example. Figure 4-9. Example of the estimation results for the modified parameters. Disturbance 3 Application 2: variable-speed wind turbine 49 An important point to be highlighted is the fact that the operating point is permanently modified due to this disturbance, meaning that the estimated values of the modified parameters might be modified with respect to their initial estimation. This point leads in some cases to a slight improvement in the accuracy of the estimation. In light of the presented estimation results provided by the CKF-based method in the base case, and the proven robustness in the disturbances studied, the estimation technique proposed in this thesis is confirmed. 4.2.5. Measurement error impact Once the accuracy of the estimator has been proven, the influence of the measurement error is tested in this section. The same external conditions as in the base case are considered in this analysis. With a 1% error in the base case, successive increments in this value have been studied (3%, 5% and 7%). Tables 4-3 to 4-5 collect the estimated values of the original parameters in the model along with their relative errors. To establish a comparison of the different noise levels, Table 4-6 summarizes the maximum relative error in each case. Note that, while low measurement errors lead to sufficiently accurate estimation results, when the errors in the measured signals are large enough, the CKF-based estimator performance deteriorates, in proportion to those errors. 50 Application of Kalman filter – based estimation techniques to electric power systems Table 4-3. Estimation errors for the parameters with a 3% measurement error Parameter Estimated value Relative error (%) 𝐻𝑡𝑚 4.015 0.369 𝐾𝑞𝑐 34.374 1.821 𝑇𝑑𝑐 0.503 0.699 𝐾𝑑𝑐 1.517 1.123 𝑇𝑑𝑠 0.505 0.885 𝑇𝑞𝑠 0.515 2.854 𝐾𝑑𝑠 1.514 0.900 𝑇𝑝 2.990 0.321 𝐾𝑝 2.025 1.237 Application 2: variable-speed wind turbine 51 Table 4-4. Estimation errors for the parameters with a 5% measurement error Parameter Estimated value Relative error (%) 𝐻𝑡𝑚 4.022 0.551 𝐾𝑞𝑐 34.13 2.547 𝑇𝑑𝑐 0.510 2.011 𝐾𝑑𝑐 1.589 5.926 𝑇𝑑𝑠 0.515 3.001 𝑇𝑞𝑠 0.526 5.224 𝐾𝑑𝑠 1.518 1.233 𝑇𝑝 2.941 1.987 𝐾𝑝 2.047 2.355 52 Application of Kalman filter – based estimation techniques to electric power systems Table 4-5. Estimation errors for the parameters with a 7% measurement error Parameter Estimated value Relative error (%) 𝐻𝑡𝑚 4.059 1.466 𝐾𝑞𝑐 33.228 5.333 𝑇𝑑𝑐 0.520 4.021 𝐾𝑑𝑐 1.662 10.800 𝑇𝑑𝑠 0.538 7.698 𝑇𝑞𝑠 0.548 9.577 𝐾𝑑𝑠 1.560 4.000 𝑇𝑝 2.832 5.912 𝐾𝑝 2.121 6.041 Table 4-6. Maximum relative error for increasing levels of measurement error Measurement error |𝐸𝑟 𝑚𝑎𝑥| (%) 1% 2.100 3% 2.854 5% 5.926 7% 10.800 Application 2: variable-speed wind turbine 53 4.2.6. Model error impact In the implemented model for the PMSM coupled to the wind turbine, the dq axis inductances of the machine, 𝐿𝑑 and 𝐿𝑞, has been taken as known. The reason is that these values can be obtained from different tests performed on the considered machine. Nevertheless, the accuracy of this information might not be adequate, leading to a source of errors in the estimation technique. In order to test the performance of the CKF-based methodology with model error, a total of four scenarios have been studied, with increasing levels of error (2%, 5%, 10% and 15%) in the values of 𝐿𝑑 and 𝐿𝑞 included in the KF model with respect to those in the simulation model. Tables 4-7 to 4-10 collect the estimated values of the original parameters in the model with their relative errors. To establish a comparison of the different model error levels, Table 4-11 summarizes the maximum relative error in each case. Note that, while low model errors lead to sufficiently accurate estimation results, when the errors in the known system parameters are large enough, the CKF-based estimator performance deteriorates, in proportion to those errors. 60 Application of Kalman filter – based estimation techniques to electric power systems 61 5. APPLICATION 3: CUSTOMER-PHASE IDENTIFICATION his section presents the application of KF-based estimation techniques applied to the identification of the electrical phase of single-phase clients in distribution grids, denominated as CPI problem. First, the proposed KF methodology is described, as well as the tuning of each of the three formulations. The considered case studies are presented subsequently, showing the corresponding estimation results. 5.1. Kalman filter implementation For the identification of the electrical phase in distribution grids, a comparative analysis is made for three formulations of the KF, the UKF, the CKF and the EnKF. The algorithms of these schemes have been coded using MATLAB framework. To approach the CPI problem with the three mentioned estimators, a novel methodology is proposed, consisting of a sequential selection and assignment of the single-phase clients, using a statistical criteria based on the calculation of a set of coefficients from the estimated values of the state variables and the covariance of the error in those estimations. This methodology is illustrated in the block diagram presented in Figure 5-1 and described in detail in section IV of the publication in annex II.3. The use of this selective method is one of the major contributions of this work, T 62 Application of Kalman filter – based estimation techniques to electric power systems avoiding the inclusion of non-convex equality constraints of the form 𝑥(𝑥−1) =0 in the estimator formulation, which can lead to convergence problems. The proposed technique is novel to the author’s knowledge and can be applied to other binary problems. Regarding the tuning of the KF-based estimators, the covariance matrices 𝑄 and 𝑅 have been taken as diagonal matrices with 𝑄𝑖𝑖 =10−4 and 𝑅𝑖𝑖 =9·10−4 in the three formulations. For the initial estimation, zero values have been used for the state vector, 𝑥0, in the UKF and CKF schemes, whereas in the case of EnKF, a binary random initialization has been considered for the samples in the ensemble. In the three cases, the covariance matrix 𝑃0 has been taken as diagonal with 𝑃𝑖𝑖 =10. The measurements taken from the simulated distribution grid are: - The hourly active energy consumption for each client in the network considered as the system input vector, 𝑢(𝑡). - The hourly active energy consumption for each of the three phases at the substation level, taken as the measurement vector, 𝑧(𝑡), for the correction stage of the KF algorithm. Additionally, while for the simulated distribution grid a load flow has been solved in order to calculate the network losses accurately, a simplified loss model has been incorporated in the filter equations in order to improve the performance of the estimation. This model assumes unitary power factor and nominal voltage for each consumer. Application 3: customer-phase identification 63 Figure 5-1. Block diagram of the proposed methodology for the CPI problem 64 Application of Kalman filter – based estimation techniques to electric power systems 5.2. Results A total of four scenarios have been considered in order to test and compare the accuracy and robustness of the different KF formulation applied to the CPI problem with the methodology developed in this work. Finally, the performance of the KF is compared with that of other existing techniques approaching the electrical phase identification. 5.2.1. Scenario 1: Original measurements In the first scenario, the original measurements from the simulated grid are taken, being the objective to compare the accuracy of the different formulation with increasing number of clients. The hourly energy measurements are obtained from [39], where real hourly data from a European distribution company, comprising smart meter readings for 20 days, are provided. Tables 5-1, 5-2 and 5-3 show the estimation results for the UKF, the CKF and the EnKF, respectively. Table 5-1. Estimation results for the UKF. Scenario 1 Consumption curves Correct assignments Wrong assignments Percentage 50 50 0 100% 100 95 5 95% 200 181 19 90.5% 300 216 84 72% 400 233 167 58.25% Application 3: customer-phase identification 65 Table 5-2. Estimation results for the CKF. Scenario 1 Consumption curves Correct assignments Wrong assignments Percentage 50 50 0 100% 100 100 0 100% 200 195 5 97.5% 300 262 38 87.33% 400 313 87 78.25% Table 5-3.Estimation results for the EnKF. Scenario 1 Consumption curves Correct assignments Wrong assignments Percentage 50 50 0 100% 100 100 0 100% 200 200 0 100% 300 283 17 94.33% 400 346 54 86.5% 66 Application of Kalman filter – based estimation techniques to electric power systems In light of the performance of UKF and CKF with increasing number of loads shown in Tables 5-1 and 5-2, respectively, it can be noticed that, while both estimators correctly assign 100% of customers for a reduced number of loads (50 clients), the behavior of the UKF deteriorates faster as the number of loads increases. As shown in Table 5-3, the EnKF response to increasing system sizes is better than that of the other formulations, which confirms the expected behavior of this KF scheme for high-dimensional problems. 5.2.2. Scenario 2: Analysis of the required amount of data For the results shown above, hourly smart meter readings for 20 days are used, leading to 480 total energy measurements for each consumer. This section analyzes the quality of the estimation as the number of available measurements decreases. Tables 5-4, 5-5 and 5-6 show the estimation results for the UKF, the CKF and the EnKF, respectively, using 200 consumption curves in each case. Table 5-4. Estimation results for the UKF. Scenario 2 Available hourly measurements Hit rate 480 90.5% 400 90.5% 300 90% 200 85% 100 67% Application 3: customer-phase identification 67 Table 5-5. Estimation results for the CKF. Scenario 2 Available hourly measurements Hit rate 480 97.5% 400 97.5% 300 96% 200 90% 100 70% Table 5-6. Estimation results for the EnKF. Scenario 2 Available hourly measurements Hit rate 480 100% 400 100% 300 100% 200 95% 100 72% 68 Application of Kalman filter – based estimation techniques to electric power systems Those results suggest a deterioration of the performance of the KF-based estimation techniques when the number of available measurements is lower than 200 (around 8 days) in this particular scenario. Further tests with feeders comprising different numbers of customers show that the required number of measurement snapshots increases with the number of loads, as intuitively expected. For instance, for 100 clients, the success rate of the three estimators does not deteriorate substantially, even when only 100 energy measurements (around 4 days) are considered for each load. 5.2.3. Scenario 3: Noisy measurements The performance of the different KF formulations is evaluated in a realistic scenario where errors in the measurements are considered. As the objective of these case studies is to determine the robustness of the KF schemes against measurement errors, a relatively low number of loads is considered, namely 100 consumption curves. Gaussian noise is artificially added to each measurement after the load flow is computed. Tables 5-7, 5-8 and 5-9 show the estimation results for the UKF, the CKF and the EnKF, respectively. Table 5-7. Estimation results for the UKF. Scenario 3 Noise level (%) Hit rate 1 95% 2 92% 3 89% 5 83% Application 3: customer-phase identification 69 Table 5-8. Estimation results for the CKF. Scenario 3 Noise level (%) Hit rate 1 100% 2 96% 3 93% 5 89% Table 5-9. Estimation results for the EnKF. Scenario 3 Noise level (%) Hit rate 1 100% 2 100% 3 98% 5 92% 76 Application of Kalman filter – based estimation techniques to electric power systems involved in the model equations. - The difference between the simulation and the estimator models, jointly with the use of external measurements for the UKF, enable the application of the proposed technique to real systems, with a set of measurements obtained from PMUs or a local data logger. - The inclusion of the PSS in the estimator model constitutes an improvement with respect to other related studies. - The estimation technique presented has not only shown good accuracy, but also robustness when modification in the Random Walk of the voltage in the external grid is considered. 6.1.2. Wind turbine In the second part of this thesis, the CKF formulation is considered for the joint estimation of the state and parameters of a wind turbine coupled to a PMSM with full converter control and pitch angle regulation. The following conclusions are obtained from the results presented in previous sections: - As it was remarked for the synchronous generator, the parameter modification methodology has successfully solved the convergence problems observed when the original parameters were included in the estimator model. - The estimation process based on CKF has provided a high degree of accuracy in the identification of the model parameters considering smooth variations in the system for a typical operating point, with the corresponding wind variability. - Three additional case studies were included to prove the robustness of the algorithm, which presented a good performance with the considered disturbances, recovering the steady state estimation for the parameters after a light transient behavior. - It was also observed that the CKF based estimator shows acceptable results when the model and measurement noises are increased. Final Conclusions 77 - Finally, the CKF formulation has shown to be superior to other KF schemes, such as the UKF and the EKF. 6.1.3. Customer-phase identification In the third part of this thesis, a technique based on Kalman filtering has been presented, applied to the identification of the electrical phase which individual loads are connected to in distribution grids, using exclusively hourly energy measurements obtained from smart meters. A total of 4 scenarios have been considered as case studies to prove the good performance of the proposed methodology in different situation. Additionally, the obtained results for the different KF formulations have been compared to those from two previously published methods. From the results presented in the chapter 5 of this document, the following conclusions are drawn: - The UKF scheme is less accurate as the number of loads in the grid increases, showing unacceptable hit rates (<60%) when 400 customers are considered. - The performance of the CKF is better than that of the UKF, with a lower sensitivity to the size of the state vector and measurement and model errors. - The EnKF formulation has proven to be the best of the considered KF schemes for the consumer-phase identification problem, both in terms of success rates and sensitivity to noise. - As it was expected, the proposed methodology is affected by the amount of energy measurements available, relative to the number of clients in the distribution grid under study. As the number of loads in a certain feeder increases, the number of required measurements for a correct phase identification also raises. - Finally, from the comparative study considered it can be concluded that the CKF and the EnKF schemes present better results than those of the 78 Application of Kalman filter – based estimation techniques to electric power systems methods based on PCA and LASSO, particularly when the number of loads is high, being the success rates similar when the grid size is small. 6.2. Future lines of investigation 6.2.1. Regulated synchronous generator Regarding the possible future investigations related to the parameter estimation in fully regulated synchronous generators, several lines are proposed: - A more complex model included in the UKF for the generator and its regulators would reflect in more detail the real behavior of the system under study, contributing to the correct operation and control of the grid. - Other promising future work would be the use of the proposed estimation technique with measurements obtained from real generation sets, given the good performance observed in the simulations. - The proposed estimation technique might be applied to more complex dynamic phenomena with unknown parameters, such as the natural frequency in high dimensional electric power systems. 6.2.2. Wind turbine Related to the parameter estimation in variable-speed wind turbines, the possible future lines of investigation are similar to those presented for the synchronous generator, i.e. the use of more complex models for the elements of the studied system, and the application of the proposed technique to real wind turbines, given the good performance in the simulations and the direct application of the proposed method, which uses exclusively external measurements. Additionally, the proposed methodology is appropriate for the parameter identification in other types of systems using power electronic converters, as might be photovoltaic generation plants or battery energy storage systems. Final Conclusions 79 6.2.3. Customer-phase identification Finally, the following point are remarked regarding the future lines of investigation related to the phase identification problem: - The proposed methodology based on Kalman filtering could be modified in order to include three-phase consumers for which only aggregated hourly energy measurements are available. - Given the potential of the existing measuremet equipments, a possible upgrade of the presented algorithm might be based on the addition of information related to the voltage or the reactive energy consumption, in order to improve the success rates of the implemented estimator. 6.3. Publications First, the publications related to the three parts of this thesis are highlighted: • [40] M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito, “Parameter estimation of fully regulated synchronous generators using Unscented Kalman Filters”, in: Electric Power Systems Research, Volume 168, 2019, Pages 210-217, November, 2018. • [41] M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito, “Parameter Estimation of Wind Turbines with PMSM using Cubature Kalman Filters,” in: IEEE Transactions on Power Systems, Volume .35, 2020, pp. 1796-1804, May 2020. • [42] M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito, “Application of nonlinear Kalman filters to the identification of customer phase connection in distribution grids,” in: International Journal of Electrical Power and Energy Systems, Volume 125, 2021. 80 Application of Kalman filter – based estimation techniques to electric power systems Related to publications [40] and [41], the following presentation at a conference special session was made: • [43] J.A. Rosendo-Macías, M.A. González-Cagigal, A. Gómez-Expósito, “PMU-Based Estimation of Renewable Power Plants Parameters.” In: IEEE Powertech Conference. Milán (Italia). 2019 Finally, as mentioned in the introduction of this document, an additional contribution was made outside the field of power systems, by applying estimation techniques based on KF for the monitorization and possible prediction of the evolution of the Covid-19. We believe some of the novel ideas proposed in this work to model and estimate slowly changing parameters, could be as well of application in the future to related power system problems. A daily report of the situation of the pandemic in Spain was published in the repository idUS, property of the University of Seville: • [44] A. Gómez-Expósito, J.A. Rosendo-Macías and M.A. González- Cagigal, “Modelo y análisis de la evolucion de una pandemia vírica mediante filtros de Kalman: el caso del Covid-19 en España”, In: https://idus.us.es/handle/11441/94508 Additionally, a preliminary version of the paper was submitted to the repository Medrxiv: Final Conclusions 81 • [45] A. Gómez-Expósito, J.A. Rosendo-Macías and M.A. González- Cagigal, “Monitoring and tracking the evolution of a viral epidemic through nonlinear Kalman filtering: Application to the Covid-19 case”, In: https://www.medrxiv.org/content/10.1101/2020.05.11.20098087v1 While the final version, included in annex II.4, and which has been submitted to “IEEE Journal of Biomedical and Health Informatics” is currently under minor revision. 82 Application of Kalman filter – based estimation techniques to electric power systems 83 ANNEX I. KALMAN FILTER APPLIED TO THE MONITORIZATION OF A VIRUS SPREAD n this annex, the motivation and the most relevant results are presented for the estimation technique proposed by the Ph. D. student and his thesis advisors for the monitorization and tracking of the evolution of a viral epidemic. Finally, the conclusions derived from the obtained results are mentioned. I.1 Motivation Despite the spectacular medical advances of the 20th century, and the practical eradication of viral diseases that in the past caused great mortality (e.g., smallpox), modern societies are still very vulnerable to the sudden appearance of new viruses, such as the SARS-CoV-2, cause of the Covid-19, for which, in the summer of 2020, there is still no vaccine. In addition, once a viral outbreak originates in a region of a country (in the case of the Covid-19, the Chinese region of Hubei, where the first reported case was dated on December 2019), the globalization of the economy and mass tourism spread it almost inevitably and quickly to the rest of the world. In the absence of effective treatments, once a certain threshold has been passed, the main and almost sole remedy against the spread of the disease to the entire population is social distancing, the objective of which is to minimize the contact between people, and therefore morbidity [46]. In extreme cases, when the speed of propagation of the outbreak is very high, massive lockdowns of entire countries may be needed, which cannot last indefinitely owing to their drastic impact on the economic activity. For this reason, all the agents involved (governments, international organizations, institutions, companies and individuals) have the greatest interest in knowing how the number of affected and deceased people will evolve over time, with a view, on the one hand, on verifying the beneficial effects of social distancing, and on the other to scheduling the already saturated health resources and taking the economic measures intended to mitigate as far as possible the devastating effects of an epidemic like that of Covid-19. Scientists, engineers, economists, etc. are acquainted with several mathematical and statistical toolkits (recently renamed collectively as "data analytics") for the treatment and filtering of time series, with a view on extracting useful information from the available data, uncertain by definition, such as trends, patterns, average values, expected variances, etc. In the specific case of a viral epidemic, such as that of Covid-19, there are basically two categories of models for processing the information: - Models that try to characterize the "physical" reality explaining the observed data., as those in [47] and [48]. I 84 Application of Kalman filter – based estimation techniques to electric - Models that try to determine explanatory parameters or variables from a purely mathematical point of view, the so-called “black box” approach, used in [49] and [50]. The methodology proposed, which belongs to the second category, uses KF to process both the assumed dynamic model and the information available throughout the outbreak. Reference [51] applies the KF for the estimation of the evolution of AIDS, while several recent studies related to the Covid-19 have arisen. In [52] the KF is used to deal with the estimation of the reproductive number of the virus. A short-term prediction model is proposed in [53], where the time update equations of the estimator are used for future forecasts of the pandemic spread. ARIMA models are combined with a KF in [54] to track the evolution of the Covid-19 in Pakistan. Unlike in those references, where the parameters involved in the state estimation process are supposed to be known, in this work such assumptions are not required. This is the major distinguishing feature of the proposed methodology, compared to the state of the art, and the main contribution of the paper included in annex II.4. I.2 Main results In this section, the most relevant results are presented. The reader may refer to the last version of the paper in annex II.4 for the detailed description of the proposed KF-based estimation technique. II.2.1 Simulation results First, a set of simulated scenarios, where the SIRD model [55-57] is considered, are used to test the performance of the proposed implementation of the KF. Figure I-1 represents the time evolution of the estimated geometric ratio, 𝑟(𝑛), along with the raw noisy measurements provided by the simulation and the actual value of 𝑟(𝑛). Figure I-1. Estimation of r(n) in the simulated scenario It can be noticed that the simulated value of 𝑟(𝑛) remains always within the 5-95% CI, giving evidence of the good performance of the proposed technique. From the estimated 𝑟(𝑛) and the initial testing 85 Annex I. Kalman filter applied to the monitorization of a virus spread ratio considered, 𝑡(0), an estimation is obtained for the evolution of the infectious people, 𝐼(𝑛), compared in Figure I-2 with the simulated value. Figure I-2. Estimation of I(n) in the simulated scenario It is observed how the maximum estimation error (around 4.5%) takes place, as expected, at the peak of the epidemic. Additionally, the effect of two error sources is studied: - An abrupt change in the testing ratio, t(n). - An error in the initial guess of the testing ratio, t(0). Both analysis confirmed the robustness of the estimation technique, which is subsequently applied to real data from different countries. II.2.2 Case studies The case studies are divided in two subsections: 1) the time period when massive lockdowns occurred in most countries, denoted in the media as the “first wave” of the pandemic [58], and 2) the subsequent transient period, once the lockdowns are relaxed, usually through several de-escalation phases, towards the so-called “new normality”. 1. Lockdown period (first wave) A total of four countries have been considered in this period: China, South Korea, Spain and the UK. At the early stage of the pandemic, the information provided by these countries was sufficient to allow the application of the proposed methodology. Contents lists available at ScienceDirect Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr Parameter estimation of fully regulated synchronous generators using Unscented Kalman Filters M.A. González-Cagigal, J.A. Rosendo-Macías ⁎ , A. Gómez-Expósito Department of Electrical Engineering, University of Seville, Spain ARTICLE INFO Keywords: Unscented Kalman Filter Synchronous generator Speed governor Automatic voltage regulator Power system stabilizer Parameter estimation ABSTRACT This paper presents a parameter estimation technique for generation sets including the synchronous machineturbine pair, along with the customary controllers: speed governor, automatic voltage regulator and power system stabilizer. The proposed technique is based on the Unscented Kalman Filter for the joint estimation of the system dynamic state and a modified set of parameters from which the actual model parameters and constants can be computed. To the authors’knowledge, this work is the first attempt to estimate such a full set of dynamic state variables and parameters, using just external measurements taken at the generator terminal bus. 1. Introduction The operation and control of electric power systems require a complete understanding of their associated dynamics, described by a set of equations. Gaining accurate knowledge of the parameters involved in those equations is essential in order to identify the system state and make adequate control decisions for the power system. Usually, those parameters are assumed to be constant, i.e., independent of factors such as exchanged power, terminal voltage, etc., characterizing the operating point. However, this assumption may not be accurate enough, owing to nonlinearities and time dependence (aging), and a method must be established to determine the updated values of system parameters taking the mentioned factors into consideration. A feasibility study is made in [1], which proves the accuracy of dynamic state estimators (DSEs) based on Kalman Filters (KFs) in power systems with nonlinear dynamics. Studies such as [2] use KFs to estimate the state of a synchronous machine, showing inaccurate estimation results in the presence of model uncertainties. Also in this context, [3] presents a comparison of different DSEs based on KFs, concluding that Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF) are computationally the most efficient for the synchronous machine dynamics. In [4], the EKF is used for state estimation in power systems with unknown inputs, but no parameter estimation is included. Particularly, among the different KF formulations, UKF has been used in [5] for state estimation of single synchronous machines, and in [6–8] for multi-machine power systems. Parameter estimation is addressed in [9,10], using a set of internal measurements which are difficult to obtain in practice. The dependence of the available measurements and the generator parameters to be estimated using UKF is studied in [11]. To avoid these difficulties, [12] proposes a method, using phasor measurement units (PMUs), for the estimation of the parameters of a synchronous machine. The most significant internal parameters of these machines are also estimated in [13], whereas [14] considers the machine under saturation conditions. The method proposed in [15] broadens the problem by including in the model the parameters of the Speed Governor (SG) and Automatic Voltage Regulator (AVR), studying the effect of a variable error (up to ± 30%) in the initialization of the parameters. The Power System Stabilizer (PSS) model is considered in [16]. Based on measurements from PMUs, the parameters associated with the generator, SG and AVR are estimated, while PSS parameters are supposed to be known in the proposed UKF implementation. To our knowledge, this work is the first attempt to consider a fully regulated generation set, including the synchronous machine itself along with the whole set of regulators (SG, AVR and PSS). Both state and parameters from these elements are jointly estimated using UKF and measurements obtained at the generator terminal bus. This work is organized as follows: Section 2formulates the equations used in the proposed UKF algorithm. Section 3presents the modeling of the overall generation set. The implementation of the UKF is described in Section 4. Section 5presents a case study to test the accuracy of the proposed estimation technique. Finally, the conclusions obtained from the estimation results are presented in Section 6. https://doi.org/10.1016/j.epsr.2018.11.018 Received 17 July 2018; Received in revised form 23 October 2018; Accepted 28 November 2018 ⁎ Corresponding author. E-mail address: [email protected] (J.A. Rosendo-Macías). Electric Power Systems Research 168 (2019) 210–217 0378-7796/ © 2018 Elsevier B.V. All rights reserved. T 2. Unscented Kalman Filter Kalman Filter implementations involve a set of state equations, including the dynamic and the measurement equations. In the case of continuous-time, discrete-measurement, non-linear systems, these equations can be expressed as =+xt fxt ut wt ˙() ( (), ()) () (1) =+zt gxt ut vt() ((), ()) ( ) kkkk (2) where x(t) is the state vector, u(t) the system input, and z(t k ) the available measurements at instant t k . The model and measurement noises, w t( ) and vt() k , are assumed Gaussian processes with covariance matrices Qand R, respectively. Considering a time step Δt, the above equations have the following discrete counterparts: =+ + −−− xx tfxu wΔ· ( , ) kk k k k 111 (3) =+zgxu v(, ) kkk k (4) which are more appropriate for non-linear Kalman filtering techniques, like the EKF and UKF. Linearizing in (3) the state function, f, has yielded poor results when combined with the approximations made by the EKF, owing to the nature of the equations that describe the behavior of the synchronous machine and associated regulators. Instead, this work makes use of the UKF, which is based on an iterative process with two different stages, as summarized in the sequel [17]. 2.1. Prediction stage For each iteration at instant k, a cloud of 2L+ 1 vectors, called σpoints, is calculated from the estimated expected value of the state vector, − x ˆ k 1 (dimension L), and the covariance matrix of the state estimation error, P k−1 , using the following expression: ⎧ ⎨ ⎪ ⎩ ⎪ = =++ =−+ =… −− −−− − +−−+ xx xx LλP xx LλP iL ˆ ˆ[( ) ] ˆ[( ) ] 1, , kk k ikki k iL kkiL 1 01 111 111 (5) where + − LλP [ ()] k i 1is the ith column of the matrix + − LλP() k1 , and λis a scaling factor calculated as =+− λ αL κ L () 2(6) with αand κbeing two filter parameters to be tuned. Those σ-points are evaluated in (3) obtaining 2L+ 1 resultant vectors, − xk i, from which to obtain the a priori estimations − x ˆkand − P k : ∑ = − = − xWx ˆ k i L k i 0 2 mi (7) ∑ =−−+ − = −−−− P Wx x x x Q(ˆ)( ˆ) k i L k ikk ikT k 0 2 ci (8) where the weighting vectors W m and W c are calculated as follows: ⎧ ⎨ ⎪ ⎩ ⎪ = =+−+ == =… + + + W Wαβ WW i L 1 1, ,2 mλ Lλ cλ Lλ Lλ 0 02 mi ci 1 2( ) (9) and βis another tunable parameter. 2.2. Correction stage On the basis of the a priori estimations, a new cloud of vectors is calculated, ⎧ ⎨ ⎪ ⎩ ⎪ = =+ + =− + =… −− −− − +− −− + xx xx LλP xxLλP iL ˆ ˆ[( ) ] ˆ[( ) ] 1, , kk k ikk i kiL kk iL 0 () (10) and evaluated with the measurement function gin (4), yielding ==… −− γgxui L (,) 0,,2 k ik ik(11) Those values are weighted using the vector W m defined in (9), ∑ = − = − zWγ ˆk i L k i 0 2 mi (12) Then, the covariance matrix of the measurement estimation error, − P z k , and the cross-covariance matrix of state and measurements, − P xzk , are obtained using the vector W c , as follows: ∑ =−−+ − =−−− − P Wγ z γ z R(ˆ)( ˆ) i L ikk ikT k zk 0 2 ci (13) ∑ =−− − = −−−− P Wx x γ z(ˆ)( ˆ) i L k ikk ik T xzk 0 2 ci (14) By using the a priori predictions at instant k, from (7),(8) and (13), (14), both the Kalman gain, =−− − K PP() kxzk zk 1 (15) and the respective a posteriori predictions can be obtained, =+ − −− xxKzz ˆˆ (ˆ) kkkk k T (16) =− −− P PKPK kkkk T zk (17) which are needed for the next algorithm iteration. 3. System dynamic modeling The joint estimation of parameters and state variables requires the knowledge of the dynamic equations modeling the behavior of the synchronous machine and associated regulators. 3.1. Dynamic model of the rotating machine A number of models can be used to represent the dynamics of a synchronous generator. Some of them are excessively complex for the purpose of this work while others do not allow the implementation of voltage control. The two-axis model in [18] is widely used in transient stability studies and has proven to be suitable for this application. After neglecting the rotor resistance, the model equations yield: =−δωω ˙Ω( ) bs (18) =−−− ω HppDωω ˙1 2(() ) me s (19) ′=′−′−− ′+eTexxiv ˙1(( ) ) q d qdddf 0(20) ′=′−′−− ′ eTexxi ˙1(( ) ) d q dqqq 0(21) where Ω b is the base frequency of the system and the state variables related to the generator are: δ: rotor angle. ω: rotor angular speed. ′ e ˙ d : transient electromotive force (emf) in d-axis. ′ e ˙q : transient emf in q-axis. The electric power, p e ,in(19) is calculated as follows: M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 211 =+ p vi vi edd qq (22) =−vδθVsin( ) d (23) =−vδθVcos( ) q(24) Vand θbeing the magnitude and phase angle of the voltage at the generator bus. Note that Eqs. (18)–(21) include a number of parameters, listed in Table 1, to be estimated. The goal of this work is to estimate as many of them as possible, using UKF and resorting to the existing knowledge about their values. 3.2. Speed governor This paper considers primary and secondary controls for the mechanical power produced by the generator turbine, characterized by [19], =−−−− p TpppRωω ˙1(1() ) mrcm s ref (25) =− p Kω ω ˙() cis (26) where the turbine actuation delay has been neglected. The parameters involved in these equations, listed in Table 2, will be also included in the UKF model. After integration, updated values of p m are used in (19). 3.3. Automatic voltage regulator The type-III model, providing a reasonable tradeoffbetween Table 1 Synchronous machine parameters. Symbol Parameter D Damping coefficient H Inertia constant x q q-axis synchronous reactance x d d-axis synchronous reactance ′ xqq-axis transient reactance ′ x d d-axis transient reactance ′ T q 0 q-axis transient time constant ′ T d 0 d-axis transient time constant Table 2 Speed Governor parameters. Symbol Parameter RGovernor droop T r Primary control time constant K i Secondary control gain Table 3 Automatic voltage regulator parameters. Symbol Parameter T avr Controller time constant T e Exciter time constant K 0 AVR gain Table 4 Power system stabilizer parameters. Symbol Parameter T w PSS time constant K w PSS gain Fig. 1. Representation of the simplified system under study. Fig. 2. Evolution of the grid voltage in the base case. Fig. 3. Evolution of the angular speed and terminal voltage. M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 212 accuracy and complexity [20], has been used for the AVR, as follows: =+−−vTKV v V v ˙1(( ) ) rsr avr 0ref (27) =+−vTvv v ˙1( ) f e rf f0(28) being vfthe AVR excitement signal, to be used in the emf equation (20), and vf 0 its initial value. Table 3 collects the AVR parameters, also to be estimated using the UKF. 3.4. Power system stabilizer For this regulator, which modifies the AVR reference signal, V ref ,a simplified type-II model has been implemented, [20]: =− − +vTKω ω v ˙1(( ) ) w ws11 (29) =−+vKωω v() sw s 1 (30) Table 4 collects the parameters of this regulator, also to be estimated using UKF. 4. Implementation of UKF State estimation requires the previous knowledge of the parameters involved in the above equations. However when these parameters are not known, dynamic estimation techniques such as UKF can be used for a joint estimation of state variables and parameters using and augmented state vector [21]. First attempts with a vanilla implementation of the joint estimator, directly using the whole set of parameters listed above, led to unacceptable convergence behavior. Therefore, as proposed in [22] for a simpler case with single-axis generator and only frequency regulation, a new modified parameter vector ψ, to be embedded in the augmented state vector, x a was adopted. This way, =xxψ[, ] aTTT , where x contains the state variables, =′′ xδωeeppvvvv[,,,, ,,,,,] Tqdmcfr s1 and =′′ ′′ ψxxDHT T [dx , dx , , , , , , , ,,,,,,,] Tdq dq dT RT K T KTT K w 0 1 1 1 100 1 1 100 q ri e w 0 avr 0avr where =− ′ xxdx dd d and =− ′ xxdxqqq. So, the size of the augmented state vector is L= 26, and the dynamic system model, replacing (3) and (4),is ⎡ ⎣⎤ ⎦=⎡ ⎣ ⎢+⎤ ⎦ ⎥+ −−− − x ψ xtfxu ψw Δ· ( , ) k k kkk k k 111 1 (31) =+zgxu v(,) kk k ak (32) where wk is now the augmented-model noise vector, including not only the state variable components, but also the parameter components. Five magnitudes can be easily measured at the generator terminal bus, namely: V,θ,ω,Iand p e , where θis always measured with respect to the synchronous reference frame (for instance, with the help of a PMU) an Iis the stator current leaving the generator. Please, note that balanced operating conditions are assumed and, hence, those magnitudes refer to the positive sequence (single-phase model). In the proposed formulation, the five magnitudes are divided into inputs, u=[V, θ], and measurements, z=[ω,I,p e ], [23]. The vector zneeds to be formulated as a function of the augmented state and input vectors, i.e., the function g(x a ,u)in(32). The first z component, the rotor speed ω, is already a state variable. For the rms value of the stator current, I, and the electric power, p e , the following expressions are used: =+Iii dq 22 (33) =+ p vi vi edd qq (34) where =−vδθVsin( ) d (35) =−vδθVcos( ) q(36) =′−′ i ev x d qq d(37) =−′+ ′ i ev x qdd q(38) As explained above, the first estimator model, based on the regular set of state variables and parameters, was prone to convergence problems during the initial tests. Trying to circumvent this problem, this work proposes a two-stage estimation technique. In the first stage, the generator parameters are estimated using a simpler model, that considers only the synchronous machine, with no regulators. In other words, the initial estimator comprises the following reduced model, =′′ xδωee[, , , ] Tqd and =′′ ′ ′ ψ xxDHT T [dx , dx , , , , , , 1] Tdq dq d q 0 0 Then, in the second stage, the generador parameters provided by the simpler model are reintroduced in the fully regulated model above, as an initial estimation. Table 5 Parameter simulation values. Parameter Unit Simulated value x d pu 1 x q pu 0.65 ′ x d pu 0.3 ′ xqpu 0.55 Dpu2 H s 6.5 ′ T q 0 s 0.5 ′ T d 0 s5 R pu 0.1 T r s 0.1 K i pu 50 T avr s1 T e s 0.2 K 0 pu 2.5 T w s10 K w pu 30 Table 6 Initial values of the state variables. Variable Initial value δ0.427 ω1 ′ e q 1.2 ′ e d 0 p m 0.8 p c 0 v f 1.68 v r 0 v 1 0 vs0 M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 213 Fig. 4. Estimated parameters for the synchronous generator at the first stage. M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 214 Fig. 5. Estimated parameters for SG, AVR and PSS. M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 215 5. Case study In this section, the performance of the proposed method is tested with the reduced power system shown in Fig. 1. The simulation has been carried out using Matlab Simulink. The model considered for the generation of measurements, including the rotor resistance, the turbine actuation delay and a set of control parameters in the PSS model, is more accurate than that adopted by the UKF algorithm [20]. In this case study, the generator is connected to an infinite busbar through a tie-line with impedance z L = 0.01 + j0.1 pu. The voltage at the infinite busbar is modeled with magnitude and angle evolving smoothly, as Gaussian random walks with standard deviations = − R 10 w5 . Two abrupt voltage magnitude steps (0.03 pu upwards and 0.05 downwards) are also applied, as shown in Fig. 2. Moreover, initial values of active and reactive power contributions to the infinite busbar are assumed, with complex power s= 0.8 + j0.4 pu, so that p ref and V ref can be accordingly obtained. As an example, the evolution of ωand Vunder these abrupt voltage magnitude variations is represented in Fig. 3. The system base frequency is 50 Hz, so Ω b = 100πrad/s and ω s = 1. The parameters adopted in the simulation models are listed in Table 5, and the time step for simulation and UKF is Δt= 0.01 s. The UKF requires an initial state for x a . The initial values for the state variables, x, are obtained by solving the steady-state case, which is given in Table 6. For the modified parameter vector adopted, ψ, random initial values are considered in the range ± 20% to ± 40% of their real values. In the second stage of the proposed two-stage estimation technique, the generator parameters are initialized with their estimated value from the first stage. The covariance of the initial estimation error is defined by a diagonal matrix = P PPdiag([ , ]) x Tψ T 000 , where the vector P x0 corresponds to the state and P ψ0 to the parameters. For the elements of P x0 , increasing values are adopted, from the generator variables to those of the PSS, resulting the vector: = −−−− −−−− P[10,10,10,10, 10 , 10 , 10 , 10 , 1, 1] x T 04444 3322 The elements of vector P ψ0 for the parameters of the machine and its regulators are tuned in order to achieve a better estimation, obtaining the following result: =−−−−−−− − P[10,10,10,10,10,10,10, 10 ,5,5,5,1,1,1,10,10] ψ T06666666 6 In this work, the UKF has been implemented with α=10 −4 ,κ=0, β= 2, following the values commonly found in works such as [24], where an analysis of the influence of these scaling parameters in the estimation process is made. The covariance matrix Rhas been taken as a diagonal matrix with R ii =10 −4 , corresponding to an error with 1% s.d. (lower accuracy than that typically provided by PMUs). The covariance matrix Qis assumed to be of the form = Q QQdiag([ , ] ) x Tψ T, where the vector Q x corresponds to the state and Q ψ to the parameters. For the elements of Q x , a common value of 10 −6 is considered, while Q ψ is proposed to have increasing values, from the generator parameters to those of the PSS, resulting the vector: = −−−−−− −−−−−− −−−− Q[10,10,10,10,10,10, 10 , 10 , 10 , 10 , 10 , 10 , 10 , 10 , 10 , 10 ] ψ T666666 665555 5544 The proposed algorithm presents a consistent performance in its ability to properly estimate the model modified parameters. Fig. 4 shows the estimation results for the generator modified parameters, as obtained from the first stage of the proposed technique. The graphics include the evolution of the estimation error covariance. For each parameter, the estimated value, x ˆ i , is represented jointly with a deviation ±xP ˆ3 iii . Note that the covariances tend to Q ii , showing the accuracy of the converged estimation. In the second stage, the synchronous generator parameters show no significant evolution from their initial values, obtained at the first stage, but the model modified parameters corresponding to the SG, the AVR and the PSS evolve successfully to the real values, as shown in Fig. 5. Hence, the method achieves a correct estimation of the modified parameters, allowing the final computation of the original parameters in (18)–(30). The relative estimation error of the original parameters is summarized in Table 7. The robustness of the proposed technique for parameter estimation has been tested with different standard deviations of the Gaussian random walks on the grid voltage, R w . For each case, the maximum relative error in the estimated parameters, E | | r max , is shown in Table 8. 6. Conclusions In this paper, an UKF implementation is proposed to perform the joint estimation of the state variables and parameters of a fully regulated two-axis synchronous machine, including SG for the frequency control and AVR and PSS for voltage control. Although many techniques, including UKF, have already been used for the estimation of generator parameters, the contribution of this work lies in the use of a modified set of parameters, comprising up to 16 parameters, which are estimated using only 5 external measurements taken at the generator terminal bus. Also, the use of a two-stage cascaded estimation process has made it possible to deal with bad convergence behavior affecting some generator parameters. A case study has shown that the proposed estimator yields accurate enough results when using different initialization strategies. Four performance tests with different standard deviations of the random walks in the grid voltage show similar relative errors on the estimated values. Acknowledgement The first author thanks the Ministry of Education and Professional Training of Spain for the financial support (grant FPU17/06380). Table 7 Relative error in parameter estimation. Parameter Estimated value Relative error (%) x d 0.987 1.30 x q 0.6583 1.28 ′ x d 0.298 0.67 ′ xq0.559 1.43 D1.999 0.05 H6.341 2.45 ′ T q 0 0.502 0.40 ′ T d 0 5.117 2.34 R0.1032 3.20 T r 0.1006 0.64 K i 49.75 0.498 T avr 0.988 1.186 T e 0.2002 0.101 K 0 2.479 0.84 T w 10.135 1.35 K w 30.78 2.596 Table 8 Maximum relative error for each standard deviation R w. R w 10 −6 10 −5 10 −4 10 −3 E | |(% ) r max 3.05 3.20 4.15 5.01 M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 216 References [1] Z. Huang, K. Schneider, J. Nieplocha, Feasibility studies of applying Kalman Filter techniques to power system dynamic state estimation, International Power Engineering Conference (IPEC), December, 2007. [2] P. Ren, H. Lev-Ari, A. Abur, Robust continuous-discrete extended Kalman filter for estimating machine states with model uncertainties, Power Systems Computation Conference (PSCC), June, 2016. [3] Z. Ning, M. Da, H. Zhenyu, G. Welch, Dynamic state estimation of a synchronous machine using PMU data: a comparative study, IEEE Power & Energy Society General Meeting, January, 2015. [4] E. Ghahremani, I. Kamwa, Dynamic state estimation in power system by applying the extended Kalman filter with unknown inputs to phasor measurements, IEEE Trans. Power Syst. 26 (November (4)) (2011). [5] S. Wang, W. Gao, A.P. Sakis Meliopoulos, An alternative method for power system dynamic state estimation based on unscented transform, IEEE Trans. Power Syst. 27 (May (2)) (2012). [6] J. Qi, K. Sun, J. Wang, H. Liu, Dynamic state estimation for multi-machine power system by unscented Kalman filter with enhanced numerical stability, IEEE Trans. Smart Grid (March) (2016). [7] F.M. Taimah, A.N. Merzah, Power system dynamic state estimation based on Kalman filter, Int. J. Comput. Appl. (0975-8887) 154 (November (11)) (2016). [8] J. Zhao, L. Mili, Robust unscented Kalman filter for power system dynamic state estimation with unknown noise statistics, IEEE Trans. Smart Grid (October) (2017). [9] A. Rouhani, A. Abur, Constrained iterated unscented Kalman filter for dynamic state and parameter estimation, IEEE Trans. Power Syst. 33 (May (99)) (2018) 1, https:// doi.org/10.1109/TPWRS.2017.2764005. [10] M. Huang, W. Li, W. Yan, Estimating parameters of synchronous generators using square-root unscented Kalman filter, Electr. Power Syst. Res. 80 (September (9)) (2010) 1137–1144. [11] G. Valverde, E. Kyriakides, V. Terzija, A non-linear approach for on-line parameter estimation of synchronous machines, 17th Power Systems Computation Conference, August, 2011. [12] Y. Wehbe, L. Fan, UKF based estimation of synchronous generator electromechanical parameters from phasor measurements, North American Power Symposium (NAPS), September, 2012. [13] E.L. Geraldi, T.C.C. Fernandes, R.A. Ramos, An approach to estimate parameters of a three-phase synchronous generator model, 2016 IEEE Power and Energy Society General Meeting (PESGM), November, 2016. [14] G. Valverde, E. Kyriakides, G.T. Heydt, V. Terzija, On-line parameter estimation of saturated synchronous machines, 2011 IEEE Power and Energy Society General Meeting, October, 2011. [15] J.C.N. Pantoja, A. Olarte, H. Daz, Simultaneous estimation of exciter governor and synchronous generator parameters using phasor measurements, Electric Power Quality and Supply Reliability Conference (PQ), June, 2014, pp. 43–49. [16] S. Dutta Chowdhury, N. Senroy, PMU data based online parameter estimation of synchronous generator, International Conference on Power Systems (ICPS), March, 2016. [17] S Simon D., Optimal State Estimation: Kalman H Infinity and Nonlinear Approaches. ISBN: 13 978-0-471-70858-2. [18] Glover J., Sarma M., Overbye T., Power Systems Analysis and Design. ISBN-13: 978- 1-111-42579-1. [19] Sauer P.W., Pai M.A., Chow J.H., Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox 2e. ISBN: 9781119355779. [20] Milano F., Power System modelling and scripting. ISBN: 978-3-642-13668-9. [21] E.A. Wan, R. Van der Merwe, The unscented Kalman filter for nonlinear estimation, Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000. AS-SPCC. The IEEE 2000, October, 2000, pp. 153–158. [22] H. Ghassempour, Z. Miao, L. Fan, W. Jiang, D. Manjure, Identification of synchronous generator model with frequency control using unscented Kalman filter, Electr. Power Syst. Res. 126 (2015) 45–55. [23] S. Julier, The scaled unscented transformation, American Control Conference, 2002. Proceedings of the 2002, IEEE Press, 2002, pp. 4555–4559. [24] J.J. Hyun, L. Hyung-Chul, Analysis of scaling parameters of the batch unscented transformation for precision orbit determination using satellite laser rangingdata, J. Astron. Space Sci. 28 (2011) 183–192, https://doi.org/10.5140/JASS.2011.28.3. 183. M.A. González-Cagigal et al. Electric Power Systems Research 168 (2019) 210–217 217 Annex II.2 M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito, “Parameter Estimation of Wind Turbines with PMSM using Cubature Kalman Filters,” in: IEEE Transactions on Power Systems, Volume .35, 2020, pp. 1796-1804, May 2020. 1802 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 35, NO. 3, MAY 2020 TABLE V RELATIVE ERROR IN PARAMETER ESTIMATION WITH SMOOTH VARIATIONS Fig. 6. Performance comparison in the estimation of shaft inertia by CKF, UKF and EKF. representing the turbine-motor inertia constant Htm, while the VSC and the pitch angle control parameters are represented in Fig.5. Foreach parameter, theestimated value, ˆxi,isrepresented along with the ˆxi±3Sii bounds. Note that the covariances tend to Qii, showing the accuracy of the converged estimation. As can be seen, the method eventually converges to the neighborhood of the correct modified parameters, from which the originalparameterscanberecovered,asperequations(18)–(32). The estimated value of the original parameters and the relative estimation error are summarized in Table V. Note that the largest relative error is 2.1%. The performance of the proposed estimation technique, based on CKF, has been compared with that shown by other KF filter schemes, such as UKF and EKF, as illustrated in Fig. 6, where the estimation of the shaft inertia under soft variations in the system is shown. Note that the EKF provides an unacceptable estimated value for the modified parameter and that the UKF presents convergence problems in the long term. Once the accuracy of the parameter estimation process has been shown under mild operating conditions, the robustness of the proposed estimation technique is tested with three large disturbances separately arising in the system under study. B. Disturbance 1: Wind Gust This test considers an abrupt change in the wind speed, which is modeled as a Mexican Hat Wavelet [31], as shown in Fig. 7. The evolution of the voltage at the point of connection is the same as in the base case. Fig. 7. Representation of the Mexican Hat Wavelet. Fig. 8. Estimation result for 10/Htm with wind gust. In this case, the parameters are initialized with their estimated values, as provided by the base case. As the initial estimation is more accurate, the elements of Pψ0are reduced to 10−4.The value of Qii has also been modified to 10−6. The estimated parameters showed no significant variations during the wind gust, reaching the same values as in the base case after a small transient, once the disturbance vanishes. As an example, the evolution of 10/Htm is provided in Fig. 8, where a zoomed view is added to better visualize the transient behavior. C. Disturbance 2: Voltage Dip The second disturbance considered is a 70% voltage dip at the point of connection that elapses for 1 s (Fig. 9). The tie-line impedance and the standard deviation of the Gaussian random walks are the same as in the previous case. The wind speed is assumed to remain as in the base case. The initial estimation for the parameters and the elements of Pψ0and Qremain the same as in the previous case (Disturbance 1). Except for a small transient arising during the disturbance, the modified parameters showed no significant variations in the presence of the voltage dip, reaching the same estimated values as in the base case. As an example, the evolution of 1/Tdc in Authorized licensed use limited to: Universidad de Sevilla. Downloaded on September 30,2020 at 06:48:51 UTC from IEEE Xplore. Restrictions apply. GONZÁLEZ-CAGIGAL ET AL.: PARAMETER ESTIMATION OF WIND TURBINES WITH PMSM USING CUBATURE KALMAN FILTERS 1803 Fig. 9. Voltage dip at the point of connection. Fig. 10. Estimation result for 1/Tdc with voltage dip. shown in Fig. 10, where the transient behavior can be noticed in the zoomed view. D. Disturbance 3: Topological Change In the last case study considered in this work, an abrupt variation in the tie-line impedance is assumed, representing a sudden topological change in the external system. For this purpose, at instant t=20s, the value of zLis modified to zL=0.01 + j0.5pu permanently. The value of the voltage at the connection point and the wind speed remain the same as in the base case. The initial estimation for the parameters and the elements of the matrices Pψ0and Qremain the same as in the disturbances 1 and 2. As the disturbance considered in this case is permanent, the estimated values are slightly different to those prior to the disturbance. See the evolution of 1/Tpin Fig. 11 as an example, where the scale of the transient period is augmented. Note that theaccuracyoftheestimatedparameterinthiscasehasimproved minimally. Fig. 11. Estimation result for Kpwith topological change. TABLE VI MAXIMUM RELATIVE PARAMETER ERROR FOR INCREASING MEASUREMENT ERRORS TABLE VII MAXIMUM RELATIVE PARAMETER ERROR FOR INCREASING ERRORS IN LdAND LqVALUES VI. MEASUREMENT AND MODEL ERROR IMPACT Once the accuracy of the proposed estimation technique has been proved, the influence of the measurement error on the CKF performance is tested. For this purpose, the smooth system perturbations considered in the base case are repeated by gradually increasing the s.d. of the measurement errors. Table VI summarizes the results obtained in terms of maximum relative error in the parameter estimation. As expected, the lower the measurement accuracy the poorer the estimation results. Note that,ifPMUsareusedasthesourceofmeasurements,witherrors typically lower than 1%, then the maximum parameter error will be less than 2%. Ontheotherhand,thesynchronousgeneratormodelpresented in Section III, as given by (25) and (26), includes the generator inductances in dq axis, Ldand Lq, which are considered to be known for the CKF implementation. The performance of the estimator with increasing error in those parameters is compared in Table VII, where the variation of the maximum relative error is shown. Note that, while low model errors lead to sufficiently accurate estimation results, when the errors in the known system parameters are large enough, the CKF-based estimator performance deteriorates, in proportion to those errors. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on September 30,2020 at 06:48:51 UTC from IEEE Xplore. Restrictions apply. 1804 IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. 35, NO. 3, MAY 2020 VII. CONCLUSION Inthis paper,aCKF isused to performa jointestimation of the state variables and parameters of a variable-speed wind turbine with a direct-drive synchronous generator and a back-to-back VSC. The proposed method includes a set of modified parameters, providing accurate estimation results (2.1% maximum relative error) under normal small variations in the operating point. A comparison of three KF schemes (CKF, UKF and EKF) has shown that the CKF is the most suitable for the parameter estimation of the system under study. The robustness of the technique has been successfully proven with three different performance tests, simulating typical large disturbances that can occur in real life. Finally, the results provided by the CKF, for increasing values of measurement errors, show that the estimator performance is acceptable as long as the quality of the measurements remains within a range typically achieved by existing synchrophasors. REFERENCES [1] Offshore Wind in Europe: Key trends and statistics 2017, European Wind Energy Association (EWEA). https://windeurope.org/wpcontent/uploads/files/about-wind/statistics/WindEurope-Annual- Offshore-Statistics-2017.pdf [2] NERC, “Reliability Guideline: Power Plant Dynamic Model Verification using PMUs,” 2016. Downloadable from: https://www.nerc.com/ comm/PC_Reliability_Guidelines_DL/Reliability%20Guideline%20- %20Power%20Plant%20Model%20Verification%20 using%20PMUs% 20-%20Resp.pdf [3] Z. Huang, K. Schneider, and J. Nieplocha, “Feasibility studies of applying Kalman Filter techniques to power system dynamic state estimation,” in Proc. Int. Power Eng. Conf., Dec. 2007, pp. 376–382. [4] A. Kumar and B. C. Pal, “Dynamic estimation and control of power systems,” ISBN: 978-0-12-814005-5. San Diego, CA, USA: Academic edition, Oct. 2018. [5] J. Zhao and L. Mili, “A robust generalized-maximum likelihood unscented Kalman filter for power system dynamic state estimation,” IEEE J. Sel. Topics Signal Process., vol. 12, no. 4, pp. 578–592, Aug. 2018. [6] A. Rouhani and A. Abur, “Constrained iterated unscented Kalman filter for dynamic state and parameter estimation,” IEEE Trans. Power Syst., vol. 33, no. 3, pp. 2404–2414, May 2018. [7] S. Wang, W. Gao, and A. P. Sakis Meliopoulos, “An alternative method for power system dynamic state estimation based on unscented transform,” IEEE Trans. Power Syst., vol. 27, no. 2, pp. 942–950, May 2012. [8] M. A. González-Cagigal, J. A. Rosendo-Macías, and A. Gómez-Expósito, “Parameter estimation of fully regulated synchronous generators using unscented Kalman Filters,” Electric Power Syst. Res., vol. 168, 2019, pp. 210–217, Nov. 2018. [9] C.Mengnan,Q.Yingning,F.Yanhui,W.Hao,andD.Infield,“Windturbine fault diagnosis based on unscented Kalman Filter,” in Proc. Int. Conf. Renewable Power Gener., Oct. 2015, pp. 1–5. [10] M. Abdelrahem, C. Hackl, and R. Kennel, “Application of extended Kalman filter to parameter estimation of doubly-fed induction generators in variable-speed wind turbine systems,” in Proc. Int. Conf. Clean Electric. Power, Jun. 2015, pp. 226–233. [11] S. P. Azad and J. E. Tate, “Parameter estimation of doubly fed induction generator driven by wind turbine,” in Proc. IEEE/PES Power Syst. Conf. Expo., Mar. 2011, pp. 1–8. [12] Q. Jilong, T. Yantao, G. Yimin, and Zhucheng, “A sensorless initial rotor position estimation scheme and an extended Kalman Filter observer for the direct torque controlled permanent magnet synchronous motor drive,” in Proc. Int. Conf. Elect. Mach. Syst., Oct. 2008, pp. 3945–3950. [13] X. Jiang, P. Sun, and Z. Q. Zhu, “Modeling and simulation of parameter identification for PMSM based on EKF,” in Proc. Int. Conf. Comput., Mechatronics, Control Electron. Eng., Aug. 2010, pp. 345–348. [14] X. Xiao and C. Chen, “Reduction of torque ripple due to demagnetization in PMSM using current compensation,” IEEE Trans. Appl. Supercond., vol. 20, no. 3, pp. 1068–1071, Jun. 2010. [15] L. Otava and L. Buchta, “Permanent magnet synchronous motor stator winding fault detection,” in Proc. IECON - 42nd Ann. Conf. IEEE Ind. Electron. Soc., Oct. 2016, pp. 1536–1541. [16] L. Otava and L. Buchta, “Implementation and verification of the PMSM stator interturn short fault detection algorithm,” in Proc. 19th Eur. Conf. Power Electron. Appl., Sep. 2017, pp. P.1–P.10. [17] G. Feng, C. Lai, J. Tjong, and N. Kar, “Non-invasive Kalman filter based permanent magnet temperature estimation for permanent magnet synchronous machines,” IEEE Trans. Power Electron., vol. 33, no. 12, pp. 10673–10682, Dec. 2018. [18] A. Echchaachouai, S. Elhani, A. Hammouch, and S. Guedira, “Extended Kalman filter used to estimate speed rotation for sensorless MPPT of wind conversion chain based on a PMSG,” in Proc. 1st Int. Conf. Elect. Inf. Technologies Proc., Jul. 2015, pp. 172–177. [19] Z. Zhang, F. Wang, G. Si, and R. Kennel, “Predictive encoderless control of back-to-back converter PMSG wind turbine systems with extended Kalman filter,” in Proc. IEEE 2nd Annu. Southern Power Electron. Conf., Dec. 2016, pp. 1–6. [20] S. A. A. Shahriari, M. Raoofat, M. Mohammadi, M. Dehghani, and M. Saad, “Dynamic state estimation of a doubly fed induction generator based on a comprehensive nonlinear model,” Simul. Model. Pract. Theory, vol. 69, p. 92–112, Dec. 2016. [21] D. G Pillai, A. Vivek, and V. Srikanth, “Non-linear state estimation of PMSM using derivative-free and square-root Cubature Kalman Filter,” in Proc. Int. Conf. Intell. Comput., Instrum. Control Technologies, Jul. 2017, pp. 126–131. [22] I. Arasaratnam, “Cubature Kalman filtering: Theory and applications,” Ph.D. Thesis, School of Graduate Studies of McMaster University, Hamilton, Ontario 2009. [23] Y. Sun and Q. Chen, “Joint estimation of states and parameters of vehicle model using cubature Kalman filter,” in Proc. IEEE Int. Conf. Syst., Man, Cybernet., Budapest, Oct. 2016, pp. 000977–000982. [24] G. R Gopinath and S. P. Das, “Speed and position sensorless control of interior permanent magnet synchronous motor using square-root cubature Kalman filter with joint parameter estimation,” in Proc. IEEE Int. Conf. Power Electron., Drives Energy Syst., Trivandrum, Dec. 2016, pp. 1–5. [25] H. Pesonen and R. Piché, “Cubature-based Kalman filters for positioning,” in 7th Workshop Positioning Navigation Commun., Mar. 2010, vol. 45, no. 49, pp. 11–12. [26] F. Milano, “Power system modelling and scripting,” ISBN: 978-3-642- 13668-9, Springer Publishing Company, Incorporated, 2010. [27] H. Ghassempour, Z. Miao, L. Fan, W. Jiang, and D. Manjure, “Identification of synchronous generator model with frequency control using unscented Kalman filter,” in Electric Power Syst. Res., vol. 126, pp. 45–55, Sep. 2015. [28] R. Gagnon and J. Brochu, “Wind Farm - Synchronous Generator and Full Scale Converter (Type 4) Detailed Model,” Mathworks. https://es.mathworks.com/help/physmod/sps/examples/wind-farm- synchronous-generator-and-full-scale-converter-type-4-detailed- model.html. Last visited: July 2019. [29] S. Julier, “The scaled unscented transformation,” in American Control Conference, 2002. Proceedings of the 2002, Piscataway, NJ, USA: IEEE Press, 2002, pp. 4555–4559. [30] J.G. Slootweg,H.Polinder,and W. L.Kling,“Initializationof windturbine models in power system dynamics simulations,” in IEEE Porto Power Tech Proc. (Cat. No.01EX502), vol. 4, pp. 6, Aug. 2001. [31] IEC 61400-1, Wind Turbines - Part 1: Design Requirements, International Electrotechnical Commission, Geneva, Switzerland, Tech. Rep., August 2005. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on September 30,2020 at 06:48:51 UTC from IEEE Xplore. Restrictions apply. Annex II.3 M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito, “Application of nonlinear Kalman filters to the identification of customer phase connection in distribution grids,” in: International Journal of Electrical Power and Energy Systems, Volume 125, 2021. Contents lists available at ScienceDirect Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes Application of nonlinear Kalman filters to the identification of customer phase connection in distribution grids M.A. González-Cagigal, J.A. Rosendo-Macías, A. Gómez-Expósito Department of Electrical Engineering, University of Seville, Spain ARTICLE INFO Keywords: Ensemble Kalman filter Cubature Kalman Filter Unscented Kalman Filter Phase identification State estimation ABSTRACT This paper presents a state estimation approach to address the problem of identifying the phase to which singlephase customers are connected in three-phase distribution grids. The proposed method performs Kalman filtering on the information provided simultaneously by the smart meter of every customer and the aggregated energy consumption measured at each phase of the secondary substation feeding the set of customers. Different nonlinear formulations of the Kalman filter are tested and their performance compared, showing that the ensemble Kalman filter provides better estimation results when the system size increases. The accuracy, robustness and limitations of the estimator are also tested when measurement errors are considered. 1. Introduction Correct operation monitoring and control in distribution systems are essential in order to assure a good quality of the service provided to the customers. In this context, it is most important for grid operators to unambiguously know which loads are connected to each of the three phases of the system (the European feeder topology is assumed in this work). An accurate connectivity information is a prerequisite to promote the correct phase balance of LV consumers, alleviating in this way the problems derived from feeder unbalances, such as sharper voltage drops, which can even violate the grid codes, and increasing power losses which also affect the lifespan of the equipment due to temperature rise. Moreover, the penetration of renewable energy resources at the distribution level also benefits from the phase identification, since it helps to establish a better production-consumption balance for each phase of the grid. In this regard, despite the efforts undertaken by distribution companies, they frequently lack enough information about the phase connection of their single-phase customers, owing for instance to network reconfiguration after faults, phase switching derived from improper maintenance, or inaccurate recording of the true load-to-phase connectivity. In these circumstances, a method must be developed to estimate as accurately as possible the actual phase to which a customer is connected in LV feeders, which is known as the customer-phase identification (CPI) problem. The CPI problem has been approached in several ways by previous works. In [1], a signal processing perspective is applied to voltage observations, which are also used both in [2], for a correlation-based methodology, and in [3], where a spectral clustering technique is proposed. The connecting phase of underground distribution transformers is determined in [4] through phase voltage measurements. Smart meters have improved the communication between the loads and the substations and can be also used for the CPI problem [5,6]. A method based on Least Absolute Shrinkage and Selection Operator (LASSO) is proposed in [7], also using smart meter data from a LV distribution network. In [8], a novel approach for phase identification using graph theory and principal component analysis (PCA) is tested. The possible missing information in smart meter data is dealt with in [9] through a correlation analysis. In this work, the CPI is addressed by applying a Kalman filtering (KF) state estimation technique to a set of hourly consumption curves obtained from real loads, along with hourly energy measurements taken at each phase of the secondary substation. The proposed method conservatively assumes that other electrical magnitudes potentially provided by smart meters, such as voltage readings or reactive power consumed by each load, are not available. Moreover, a simplified loss model is adopted, allowing the impact of each load on the energy delivered by each phase of the transformer to be estimated for a given topology. The KF is a dynamic state estimator (DSE) widely used in electric power systems in any of its diverse forms [10–12]. Particularly, the application of KF to parameter estimation in power systems has been successfully tested, e.g. in [13–15], providing evidence of the DSE potential for the CPI problem. In this work, three nonlinear KF schemes are tested and compared, namely the so-called unscented Kalman filter (UKF), the cubature Kalman filter (CKF), and the ensemble Kalman filter (EnKF). The proposed technique takes advantage of the hourly information provided by smart meters to sequentially assign customers to the most likely phase. The proposed KF implementation does not explicitly enforce, beforehand, binary constraints for the state variables, but rather adopts a novel, statistically-based inference logic successively rounding state https://doi.org/10.1016/j.ijepes.2020.106410 Received 26 December 2019; Received in revised form 6 July 2020; Accepted 30 July 2020 Electrical Power and Energy Systems 125 (2021) 106410 0142-0615/ © 2020 Elsevier Ltd. All rights reserved. T variables to their nearest binary value. As the input information considered in this work is similar to that assumed in [7,8], the KF-based estimation is also compared to those competing techniques for different scenarios. The paper is organized as follows: Section 2provides a brief background on the different KF schemes considered. Next, the feeder used for testing and the simplified loss model adopted are presented in Section 3. The implementation details of the KF based estimation techniques, as applied to the CPI problem, are given in Section 4.In Section 5, the results obtained in two different scenarios are presented and discussed, while the proposed KF technique is compared in Section 6with other published works dealing with CPI. The conclusions derived from the case studies considered are drawn in Section 7. 2. Kalman filter background In this section, three different Kalman filter formulations are succinctly reviewed, as applied to continuous-time, discrete-measurement nonlinear systems. In the discrete-time framework, the associated equations may be expressed as follows, =+ −− xfxu w(, ) kkk k 11 (1) =+zgxu v(, ) kkk k (2) where x k is the state vector at instant ku, k the system input, and z k the vector of available measurements. Gaussian processes are considered for the model and measurement noises, wk and v k , with covariance matrices Qand R, respectively. The iterative processes of the Kalman filter schemes considered in this work, all of them involving prediction and correction stages, are described below. 2.1. Unscented Kalman Filter At instant k, a cloud of +L2 1 vectors, the so-called σ -points, is obtained from the previously estimated state vector,  − x k 1 (dimension L), and the covariance matrix of the state estimation error, − P k 1 , as follows [16]:    ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ = =++ =−+ =… −− −−− − +−− + xx xx LλP xx LλP iL [( ) ] [( ) ] 1, , kk k ikk i k iL kk iL 1 01 111 111 (3) + − LλP [ ()] k i 1 being the i th column of the corresponding matrix, and λ a scaling factor calculated as follows: =+− λ αL κ L () 2 (4) where α and κ are two filter parameters to be tuned. These σ -points are evaluated through Eq. (1), yielding +L2 1 vectors, − xk i, from which the a priori estimations  − xkand − P k are obtained: ∑ = − = − xWx k i L mi k i 0 2 (5)   ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟−+ − = −−−− P Wx x x x Q() k i L ci k ikk ik T k 0 2 (6) where the weighting vectors W m and W c are calculated from: ⎧ ⎨ ⎪ ⎩ ⎪ = =+−+ == =… + + + W Wαβ WW i L 1 1, ,2 mλ Lλ cλ Lλ mi ci Lλ 0 02 1 2( ) (7) β being another tunable parameter. On the basis of the a priori information, the correction stage starts with the calculation of a new cloud of vectors, − x k , which are evaluated with the measurement function g(· ) in Eq. (2), and weighted with the vectors W m , yielding ==… −− γgxui L (,) 0,,2 k ik ik(8) ∑ = − = − zWγ k i L mi k i 0 2 (9) Then, the covariance matrix of the measurement estimation error, − P z k , and the cross-covariance matrix of state and measurements, − P xz k , are obtained using the vector W c as follows:   ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟−+ − = −−−− P Wγ z γ z R() zk i L ci k ikk ik T k 0 2 (10)  ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟− − = −−−− P Wx x γ z() xzk i L ci k ikk ik T 0 2 (11) The correction stage concludes with the a posteriori predictions,    =+ − −− xxKzz() kkkk k T (12) =− −− P PKPK kkkzk k T (13) which are based on the a priori predictions at instant kand the so-called Kalman gain, K k , calculated from =−− − K PP() kxzk zk 1 (14) 2.2. Cubature Kalman Filter This KF formulation uses a set of L 2cubature points calculated from  − x k 1 and − P k 1 through the following expressions [17]: = −−− SS P kk Tk11 1 (15)  =+=… −−− xSξLxi L1, ,2 k ikik 111 (16) where Sis a positive-definite square root of matrix P(the Cholesky factorization of matrix Pis customarily used), and ξiis the i th cubature node, obtained as the intersection of the unit sphere and the L  axis. The state function f (·) in (1) is evaluated for the set of cubature points, yielding a set of L 2vectors − xk i, from which the a priori estimation is computed, ∑ = − = − xLx 1 2 k i L k i 1 2 (17)   ∑ =−+ − = −− −− P Lxx xx Q 1 2 k i L k ik iT kk T k 1 2 (18) For the correction stage, the covariance matrix − P k is factorized in order to calculate both the matrix − S k , = −− − SS P kk Tk (19) and a new set of L 2cubature points, − x k , at which function g(· ) in (2)is evaluated to obtain − γ k . Then the measurement estimation,  − z k , its covariance matrix, − P z k , and the cross-covariance matrix of state and measurements, − P xz k , are calculated as follows: ∑ = − = − zLγ 1 2 k i L k i 1 2 (20)   ∑ =−+ − = −− −− P Lγγ zz R 1 2 zk i L k i k iT kk T k 1 2 (21) M.A. González-Cagigal, et al. Electrical Power and Energy Systems 125 (2021) 106410 2  ∑ =− − = −− −− P Lxγ xz 1 2 xzk i L k i k iT kk T 1 2 (22) The a posteriori predictions of the state vector,  x k , and the covariance P k are calculated with the Kalman gain using the same Eqs. (14) and (13) as in the UKF algorithm. 2.3. Ensemble Kalman Filter The EnKF [18], is a Monte Carlo approximation of the original KF which has proven accurate enough in high-dimensional state-space problems. The ensemble is represented by an LxNmatrix, Nbeing the number of samples considered. The ensemble is first propagated through the state and measurement functions, ⎧ ⎨ ⎩ = = =… −−− − xfxu zgxu iN (, ) (, ) 1, , k ik ik k ik ik 11 (23) and then the mean values are calculated: ∑ = − = − xNx 1 k i N k i 1(24) ∑ = = − zNz 1 k i N k i 1 (25) The EnKF correction stage is based on the calculation of the intermediate matrices ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟− = −−− P HNxxzz 1() kk T i N k ikk ikT 1 (26) ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟− = −− HPH Nzzzz 1() kk k T i N k ikk ikT 1 (27) which allow the Kalman gain and the updated values of each sample in the ensemble to be obtained: =+ − K PH HPH R() kk k Tkk k T 1 (28) =+ − =… −− xx Kzzi N ()1,, k ik ikk k i(29) Finally, the corrected covariance matrix, P k , is calculated as follows: ∑ = = xNx 1 k i N k i 1 (30) ∑ =⎛ ⎝ ⎜−⎞ ⎠ ⎟− = P Nxxxx 1() k i N k ikk ik T 1(31) 3. Problem statement and modeling In order to assess the ability of the KF-based methodology to address the CPI problem, a typical distribution grid is considered comprising N s single-phase and N tthree-phase customers, not necessarily balanced, resulting in =+ N NN3 cs t total consumption curves in the network. In this work, it is assumed that the number of three-phase customers is 20 % of the total. The proposed methodology, as presented in Section 4,is easily extensible to other load connection arrangements, such as twophase loads, still found in some areas. Fig. 1 shows an example of a distribution network with 100 loads, used in the sequel to test the estimation techniques. The energy consumption of each load iat a certain hour k, denoted by E i k ,, is obtained from [19], where real hourly data from a European distribution company, comprising smart meters readings for 20 days, are provided, leading to a total of 480 energy measurements for each customer. In case of three-phase loads, a single-phase consumption is assigned to each phase. As customers with null consumption cannot be identified (they provide no information), the corresponding curves are removed from the raw data. To fully characterize the energy consumption of the customers, the reactive energy is also obtained from the raw data in [19]. Finally, the resulting hourly curves are randomly associated to a certain phase (a, b or c). Along with the previous customer information, a typical distribution feeder topology is considered, each customer being associated to one of the grid nodes. Then, a load flow can be solved at each hour kin order to obtain the energy delivered by each phase of the MV/LV secondary substation, E E, SS k aSS k b ,, and E SS k c,, which is used by the proposed KF-based estimation technique as a measurement. This fully defines the distribution grid model involved in the estimation process. 4. Kalman Filter Implementation The application of the different KF schemes described in Section 2to the CPI problem is illustrated in the flowchart represented in Fig. 2,as summarized in the sequel. 4.1. Parallel filtering For every iteration of each tested technique, three independent KF- based estimators run in parallel, one for each phase pof the distribution grid, being the state vector composed of N c variables, x i p , with = p abc,,, so that =LN c . Any three-phase customer is characterized by a set of three consecutive state variables, ++ xx x{,,} ipipip 12 , for each phase. Those variables indicate if the consumption iis associated to the corresponding phase p. For these static variables, the state function, f (·), in Eq. (1), is taken as a random walk, independent of the system input, u k , yielding the following expression: =+ − xx w k p k p k p 1(32) where the model noise, w k p , is considered to have the same covariance matrix, Q, for each phase. The system input, u k , used in the measurement function g(· ) in Eq. (2), is determined by the hourly energy consumption of each customer, =… u EE E,,, k TN k 12 c . Regarding the network losses, as no measurements about the voltage and power factor of each load are available, a Fig. 1. Single-line diagram of one of the test networks. M.A. González-Cagigal, et al. Electrical Power and Energy Systems 125 (2021) 106410 3 simplified loss model is adopted by the KF estimators, for which a given rated voltage and power factor =φcos 1 ik, are assumed for all loads points. With these assumptions, an average hourly current, Ii k ,is calculated for each consumer ifrom the corresponding hourly energy consumption, E i k ,: =IE TU· ik ik , , (33) where Uis the voltage magnitude, and Tis the energy integration period. In this work, as the energy measurements are supposed to be obtained hourly, =T h 1. Accordingly, the energy loss attributable to customer iis given by: = E ITrl··· ik loss ik i ,, 2(34) where ris the conductor resistance per unit length and li is the electrical distance from each consumer to the secondary substation, both assumed to be known ( =r0.223 Ω/k m , has been adopted in this paper). The vector u k also includes the variables x i p corresponding to customers which have been previously assigned to a phase, as described below. The measurement of each filter, z k p , is taken as the energy delivered at each hour by the corresponding phase of the secondary substation, =zE k p SS k p ,, so the measurement function g(· ) in Eq. (2) reduces to: ∑ =⎛ ⎝ ⎜+⎞ ⎠ ⎟+ = zxEEv· k p i N ik pik ik loss k p 1 ,,, c (35) where the measurement noise, v k p , is considered to have the same covariance matrix, R, for each phase. 4.2. Initialization and tuning As explained formerly, before the KFs are applied, a set of input data are gathered, including the hourly consumption of each customer and the energy delivered by the MV/LV transformer for the time window available. A separation is made between single-phase and three-phase loads. For the last ones, it is assumed that three energy measurements are available, one per phase, but the phase labels (a, b or c) are unknown. At the end of the estimation processes, the actual distribution of customers among the three phases is used to evaluate the performance of each estimation procedure. The KF formulations require that an initial estimation is adopted, which is determined by the state vector  x 0 and the covariance matrix of the initial estimation error, P 0. For the UKF and CKF implementations, all state variables in the three phases are initialized as 0, assuming a complete lack of knowledge about their real values. In the case of EnKF, the samples are given random binary values as initial estimation, which enhances the convergence of the estimation algorithm. The covariance matrix, P 0is considered as a diagonal matrix with = P 10 ii in the three KF schemes. The covariance matrices Qand Rare defined as diagonal matrices, considering typical values of = − Q 10 ii 4 and = − R 9·10 ii 4 respectively, the last one being equivalent to assuming a s.d. of 3% for the measurement errors. For the estimation based on UKF, it is necessary to define the tunable parameters of the filters introduced in Eqs. (4) and (7). A study is made in [20] over the influence of these parameters in the estimation process, concluding that == − α β10 , 2 4and =− κ N3 c , are reasonable values for good estimation results. Regarding the EnKF, the number of samples in the ensemble is taken as = N N10· c so that it can suit the different sizes of the state vector. 4.3. Candidate selection and assignment Once the three parallel estimators have provided a customer distribution for their corresponding phase, a single consumption is assigned and removed from the subsequent estimation processes. The candidate selection at each iteration of the CPI procedure is based on the estimated values of the state variables,  xi p , and the covariance of their estimation error, P ii. Forcing the state variables to be only 0 or 1, as required by the nature of the problem in hand, would involve equality constraints of the form   −=xx·( 1) 0. However, the application of the KF estimation algorithms to such non-convex model is prone to convergence problems. For this reason, the approach proposed in this work to enforce the binary character of the state variables relies on the confidence level of a given consumption being associated to a certain phase and not to the others. This is quantified through the so-called phase confidence factors, C i p , calculated with the probability density function (PDF) of the state variables, as follows. First, the cumulative density function for >x0. 5 ip, denoted as M i p , is calculated as follows:  ∫ ⎧ ⎨ ⎪ ⎩ ⎪ = = =… ∞− − Medx pabc iN ,, 1, , ip πP c 1 2· · 0.5 ii xx ip Pii () 2 2· Since M i p provides information on the discrete value associated to the corresponding state variable, this coefficient is called phase discretization factor. Fig. 3 includes a graphic representation of this factor, corresponding with the shaded area in the Gaussian density function, considering 0.5 as lower limit. Then, for each energy consumption curve i, the phase confidence factor combine the information derived from the three phases of the corresponding load through the following expression: Fig. 2. Flowchart of the CPI methodology. M.A. González-Cagigal, et al. Electrical Power and Energy Systems 125 (2021) 106410 4 JBHI-00622-2020.R2 1 Abstract— This work presents a novel methodology for systematically processing the time series that report the number of positive, recovered and deceased cases from a viral epidemic, such as Covid-19. The main objective is to unveil the evolution of the number of real infected people, and consequently to predict the peak of the epidemic and subsequent evolution. For this purpose, an original nonlinear model relating the raw data with the timevarying geometric ratio of infected people is elaborated, and a Kalman Filter is used to estimate the involved state variables. A hypothetical simulated case is used to show the adequacy and limitations of the proposed method. Then, several countries, including China, South Korea, Italy, Spain, UK and the USA, are tested to illustrate its behavior when real-life data are processed. The results obtained clearly show the beneficial effect of the severe lockdowns imposed by many countries worldwide, but also that the softer social distancing measures adopted afterwards have been almost always insufficient to prevent the second virus wave. Index Terms - Nonlinear Kalman filtering, parameter estimation, Covid-19, geometric series. I. INTRODUCTION espite the spectacular medical advances of the 20th century, and the practical eradication of viral diseases that in the past caused great mortality (e.g., smallpox), modern societies are still very vulnerable to the sudden appearance of new viruses, such as the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), cause of the coronavirus disease 2019 (Covid-19), for which there is still no vaccine. In addition, once a viral outbreak originates in a region of a country (in the case of the Covid-19, the Chinese region of Hubei, where the first reported case was dated on December 2019), the globalization of the economy and mass tourism spread it almost inevitably and quickly to the rest of the world. In the absence of effective treatments, once a certain threshold has been passed, the main and almost sole remedy against the spread of the disease to the entire population is social distancing, the objective of which is to minimize the contact Manuscript submitted for review on May 11, 2020. The third author thanks the financial support of the Spanish Ministry of Education and Professional Training (grant FPU17/06380). between people, and therefore morbidity [1]. In extreme cases, when the speed of propagation of the outbreak is very high, massive lockdowns of entire countries may be needed, which cannot last indefinitely owing to their drastic impact on the economic activity. For this reason, all the agents involved (governments, international organizations, institutions, companies and individuals) have the greatest interest in knowing how the number of affected and deceased people will evolve over time, with a view, on the one hand, to verifying the beneficial effects of social distancing, and on the other to scheduling the already saturated health resources and taking the economic measures intended to mitigate as far as possible the devastating effects of an epidemic like that of Covid-19. Scientists, engineers, economists, etc. are acquainted with several mathematical and statistical toolkits (recently renamed collectively as "data analytics") for the treatment and filtering of time series, with a view to extracting useful information from the available data, uncertain by definition, such as trends, patterns, average values, expected variances, etc. In the specific case of a viral epidemic, such as that of Covid-19, there are basically two categories of models for processing the information: 1. Models that try to characterize the "physical" reality explaining the observed data. In the case of a viral epidemic, these models [2], [3] consider, for example, what fraction of people are at work, in teaching or travelling, how long it takes for an infected person to manifest symptoms, what is the mortality rate according to age groups, etc. This type of modeling is widely used in engineering, because the dynamics of the underlying systems or devices are generally well characterized, through mathematical relationships obtained from the physical laws that govern them (such is the case, for example, of electrical networks or an artificial satellite). 2. Models that try to determine explanatory parameters or All authors are with the Department of Electrical Engineering, Laboratory of Engineering for Energy and Environmental Sustainability, University of Seville, Spain (e-mails: ag[email protected], [email protected], [email protected]). . Monitoring and Tracking the Evolution of a Viral Epidemic Through Nonlinear Kalman Filtering: Application to the Covid-19 Case A. Gómez-Expósito, Fellow, IEEE, J. A. Rosendo-Macías, Senior Member, IEEE, and M. A. González-Cagigal D JBHI-00622-2020.R2 2 variables from a purely mathematical point of view ("black box" approach), without going into the causes or interactions between components that explain the resulting data. Given uncertain data, which enter the system regularly (in our case, every day), the aim is to characterize its temporal evolution by adjusting the parameters of a mathematical model, so that the differences between what is observed and what is estimated are minimized. Two variants can be considered in this category: 2.a) Mathematical models that do not assume a priori what the shape of the temporal evolution of the involved magnitudes will be, but rather use a state transition equation, which tries to capture the dynamics of the problem in question by relating the variables in an instant of time to the variables in the previous instant. In this case, it is a matter of determining how the coefficients that define this equation evolve over time. In the case of epidemics, among the most popular models are those derived from the SIR (Susceptible-Infected-Recovered) model, [4], such as the one used for example in [5] to analyze the evolution of Covid-19 in Italy. This model is also considered in [6], where the evolution of the epidemic is forecasted using a novel state filtering algorithm. 2.b) Mathematical models based on the assumption that the evolution of infected people, deceased, etc. obeys a predetermined curve (based on the experience of previous epidemics), whose coefficients are estimated based on the time series of reported data. For example, the evolution of the accumulated number of infected people can be satisfactorily approximated by means of a sigmoid curve, as assumed in [7], where the curve proposed by Gompertz [8] is used. The methodology proposed in this work belongs to the second category. As explained in the next section, we depart from the basic SIR model, by considering that the number of susceptible people, being large enough and changing relatively slowly, does not have to be explicitly considered in the model, but can be rather embedded in other equally significant parameters, such as the time-varying ratio of the geometric series characterizing the progression of affected people. Moreover, the proposed model explicitly distinguishes between people who have proved positive in a test, and people actually infectious, who are many more and for whom there is no reliable information available. In this work, a Kalman filter (KF) is used to process both the assumed dynamic model and the information available throughout the outbreak. The KF, proposed for linear dynamic systems in the early 1960s, is considered one of the fundamental tools that allowed men to walk on the moon, as it was successfully used in guiding the Apollo program space missions [9]. This filter, which constitutes a generalization of the technique known as "recursive least squares", estimates the maximum likelihood evolution (that is, the most statistically probable, according to the assumed uncertainties and the observed samples) of the state of a dynamic system, and can be generalized to the non-linear case, including situations where the model parameters are also to be estimated. Reference [10] applies the KF for the estimation of the evolution of AIDS, while several recent studies related to the Covid-19 have arisen. In [11] the KF is used to deal with the estimation of the reproductive number of the virus. A shortterm prediction model is proposed in [12], where the time update equations of the estimator are used for future forecasts of the pandemic spread. ARIMA models are combined with a KF in [13] to track the evolution of the Covid-19 in Pakistan. Unlike in those references, where the parameters involved in the state estimation process are supposed to be known, in this work such assumptions are not required. This is the major distinguishing feature of the proposed methodology, compared to the state of the art, and the main contribution of the paper. II. PROPOSED MODEL We start from the well-known and simple SIR model [4], mathematically described by: 𝑆󰇗𝑐(𝑡) = −𝛽 ∙𝑆𝑐(𝑡)𝐼(𝑡)/𝑁 𝐼󰇗(𝑡) = 𝛽 ∙ 𝑆𝑐(𝑡)𝐼(𝑡)/𝑁 −𝛾 ∙ 𝐼(𝑡) 𝑅󰇗𝑐(𝑡) = 𝛾 ∙𝐼(𝑡) where 𝑆𝑐(𝑡) and 𝑅𝑐(𝑡) are, respectively, the cumulative or total susceptible and recovered people, 𝐼(𝑡) represents the active infectious (not to be confused with cumulative infectious), β and γ are the transmission and recovery rates, and 𝑁 is the total population of the studied region, satisfying 𝑁 = 𝑆𝑐(𝑡)+ 𝐼(𝑡)+𝑅𝑐(𝑡). Note that, in this compact model, the deceased cases are paradoxically included in 𝑅𝑐(𝑡) (alternatively, they could be subtracted from 𝑁). For practical purposes, the discrete counterparts obtained by numerical integration (forward Euler) are rather of interest. Moreover, as dead people are separately reported, they can be explicitly modeled, leading to a discrete-time SIRD (Susceptible-Infected-Recovered-Deceased) model, as used in [14-16]: 𝑆𝑐(𝑛 +1)= 𝑆𝑐(𝑛)−𝛽 ∙ 𝑆𝑐(𝑛)𝐼(𝑛)/𝑁 (1) 𝐼(𝑛 +1)= 𝐼(𝑛)+𝛽 ∙ 𝑆𝑐(𝑛)𝐼(𝑛)/𝑁−(𝛾 +𝜇)𝐼(𝑛) (2) 𝑅𝑐(𝑛 +1)= 𝑅𝑐(𝑛)+𝛾 ∙ 𝐼(𝑛) (3) 𝐷𝑐(𝑛 + 1)= 𝐷𝑐(𝑛)+𝜇 ∙ 𝐼(𝑛) (4) where 𝑛 is the elapsed time (in days) from a given origin, 𝐷𝑐(𝑛) is the cumulative dead and 𝜇 is a mortality ratio. The data publicly reported (available in references such as [17],[18]), typically comprise the following three items: • Fraction of infectious people who, subject to a test, yield a positive outcome. This considers the fact that there may be many more infected than those reported positives, as happens with a large number of asymptomatic people. The cumulative positives will be denoted 𝑃𝑐(𝑛). JBHI-00622-2020.R2 3 • Fraction of recovered people who have been previously identified as positive. For simplicity of notation, the same symbol as in the basic SIR model, 𝑅𝑐(𝑛), will be used in the sequel, even though we are referring here to a subset of 𝑅𝑐. • Cumulative number of deceased, 𝐷𝑐(𝑛), which is assumed to be the same as in the SIRD model, even though the actual number of dead by the virus may differ from the reported figures. Some sources directly provide the active positive cases, 𝑃(𝑛), defined as: 𝑃(𝑛)= 𝑃𝑐(𝑛)− 𝑅𝑐(𝑛) − 𝐷𝑐(𝑛). Note that both 𝐼(𝑛) and 𝑃(𝑛) tend to zero as 𝑛 increases sufficiently (end of the viral outbreak), while the remaining cumulative magnitudes asymptotically reach a maximum or steady-state value. Epidemiologists use the so-called basic reproductive number, 𝑅0 (average number of people infected by a single infectious person during the infective period at the onset of the outbreak) to characterize whether and how fast an epidemic spreads at the very beginning. If 𝑅0>1, then the epidemic will progress exponentially. As time elapses, though, the number of susceptible people decreases, either by the virus evolution itself or as a consequence of social distancing measures, and 𝑅0 is replaced by the effective reproductive number, 𝑅𝑡. In terms of SIRD coefficients, 𝑅0 and 𝑅𝑡 are given by: 𝑅0= 𝛽/(𝛾 + 𝜇) ; 𝑅𝑡(𝑛) = 𝑅0𝑆𝑐(𝑛)/𝑁 However, as thoroughly discussed in [19], the basic reproductive number 𝑅0 is not free from ambiguity and controversy. For instance, it is stated in [19] that "using 𝑅0 as a threshold parameter for a population-level model could produce misleading estimates of the infectiousness of the pathogen, the severity of an outbreak, and the strength of the medical and/or behavioral interventions necessary for control". Moreover, if 𝑅0 is estimated from time series of reported data, as in [11], then there is no way, at least for a new virus such as Covid-19, to subsequently check or contrast the accuracy of the estimates. This probably explains the wide confidence intervals so far reported for 𝑅0 values [20]. Similar arguments apply to 𝑅𝑡. For this reason, instead of or in addition to 𝑅0, we postulate in this work the use of a more intuitive and measurable index, related with the growth rate of the infected class, to duly and unambiguously characterize a viral epidemic. Let the daily evolution of the active infectious be expressed as a geometric time series: 𝐼(𝑛 + 1) = 𝒓(𝒏)∙𝐼(𝑛) (5) where 𝑟(𝑛) is the time-varying ratio of the series. Then the daily growth rate is obtained from: Growth rate (p.u.) = ∆𝐼(𝑛)/𝐼(𝑛) = 𝑟(𝑛) −1 Clearly, as long as 𝑟(𝑛)>1, the viral outbreak will continue its expansion, whereas the disease extinguishes when 𝑟(𝑛)<1. There is no ambiguity in using 𝑟(𝑛) as a threshold, when referred to a whole population. Note however that, if (5) were expressed in terms of cumulative magnitudes, rather than daily or active cases, then 𝑟(𝑛) would tend asymptotically to 1. By direct comparison of (5) with (2), the following relation is obtained, 𝑟(𝑛) = 1 +𝛽 ∙𝑆𝑐(𝑛)/𝑁− (𝛾 + 𝜇) or, in terms of 𝑅𝑡: 𝑟(𝑛)= 1+ (𝛾 +𝜇)[𝑅𝑡(𝑛)−1] Given 𝑟(𝑛), one still would have to guess the values of the parameters involved in the SIRD model (1)-(4), to obtain 𝑅0. We contend that there is no need to worry in the short term about 𝑅0, as 𝑟(𝑛) suffices to duly track the epidemic evolution on a daily basis. This work is aimed at estimating, from the daily reported data, the evolution of 𝑟(𝑛) and, as a consequence, the growth rate of the infectious people. Note that, if 𝑟(𝑛) can be somehow estimated, then equation (1) becomes unnecessary. In our approach, the impact of susceptible people, a factor which varies smoothly, is also embedded into 𝑟(𝑛). In order to take advantage of the reported numbers of positive, deceased, and recovered cases, the following relationships are considered, taking into account (3)-(4): 𝑃(𝑛)= 𝑡(𝑛) ∙𝐼(𝑛) (6) 𝐷(𝑛)= 𝜇(𝑛)∙𝐼(𝑛) (7) 𝑅(𝑛)= 𝑡(𝑛)∙𝛾(𝑛)∙𝐼(𝑛)=𝛾𝑡(𝑛) ∙𝐼(𝑛) (8) where 𝑡(𝑛) is a testing or reporting ratio that models the fraction of those infectious who are subject to tests and yield positive, 𝐷(𝑛)= 𝐷𝑐(𝑛)− 𝐷𝑐(𝑛 − 1) is the daily increase in the number of deaths, and 𝑅(𝑛)= 𝑅𝑐(𝑛)−𝑅𝑐(𝑛 − 1) is the daily variation in the number of recovered cases. From (7) at two consecutive instants, keeping (5) in mind: 𝐷(𝑛 +1)= 𝜇(𝑛 +1)∙𝑟(𝑛)∙𝐼(𝑛) 𝐷(𝑛)= 𝜇(𝑛)∙𝐼(𝑛) and dividing: 𝑟𝐷(𝑛)= 𝑟(𝑛)∙𝐾𝜇(𝑛) (9) where 𝑟𝐷(𝑛)= 𝐷(𝑛 + 1)/𝐷(𝑛) is the ratio of consecutive daily deaths and 𝐾𝜇(𝑛)= 𝜇(𝑛 + 1)/𝜇(𝑛) is in turn a ratio of consecutive mortality ratios. JBHI-00622-2020.R2 4 Similarly, from (6) and (5): 𝑃(𝑛 +1)= 𝑡(𝑛 + 1)∙𝑟(𝑛)∙ 𝐼(𝑛) 𝑃(𝑛)= 𝑡(𝑛)∙𝐼(𝑛) and dividing: 𝑟𝑃(𝑛)= 𝑟(𝑛)∙𝐾𝑡(𝑛) (10) where 𝑟𝑃(𝑛)= 𝑃(𝑛 + 1)/𝑃(𝑛) is the ratio of consecutive daily positives and 𝐾𝑡(𝑛)= 𝑡(𝑛 + 1)/𝑡(𝑛) is a ratio of consecutive 𝑡(𝑛) coefficients. Finally, from (8) and (5): 𝑅(𝑛 +1)=𝛾𝑡(𝑛 + 1)∙𝑟(𝑛)∙𝐼(𝑛) 𝑅(𝑛)=𝛾𝑡(𝑛)∙𝐼(𝑛) and dividing: 𝑟𝑅(𝑛)= 𝑟(𝑛)∙𝐾𝛾(𝑛) (11) where 𝑟𝑅(𝑛)=𝑅(𝑛 + 1)/𝑅(𝑛) is the ratio of daily recovered cases and 𝐾𝛾(𝑛)=𝛾𝑡(𝑛 + 1)/𝛾𝑡(𝑛) is a ratio of consecutive 𝛾𝑡(𝑛) coefficients. III. KALMAN FILTER APPLICATION AND IMPLEMENTATION From the above simple model, given by (9)-(11), a system of state equations can be formulated allowing the sequence 𝑟(𝑛) to be estimated by means of a non-linear KF, such as the Extended KF (EKF), Unscented KF (UKF) or Ensemble KF (EnKF) [9], [21]. This type of filter is capable of estimating the dynamic evolution of both parameters and state variables. Even though its maximum likelihood has only been proven for linear problems, it is applied successfully in nonlinear problems, such as the one in hand. In our proposal, the state vector is composed of the infectious geometric ratio, 𝑟(𝑛), not directly measured, and the gains 𝐾𝜇(𝑛), 𝐾𝑡(𝑛) and 𝐾𝛾(𝑛), all of them initially assumed to evolve with a random walk. Therefore, the state vector is 𝑥 = [𝐾𝜇𝐾𝑡𝐾𝛾𝑟]𝑡 and the state transition equation: 𝑥(𝑛)= 𝑥(𝑛 −1)+𝑤(𝑛) where 𝑤 = [𝑤𝜇𝑤𝑡𝑤𝛾𝑤𝑟]𝑡 is a Gaussian noise vector with a covariance matrix Q(n), which accounts for model errors such as possible time lags not duly considered in the model. The state is to be estimated with the help of a measurement vector, composed of 𝑟𝐷(𝑛), 𝑟𝑃(𝑛) and 𝑟𝑅(𝑛), along with the pseudo-measurements of 𝐾𝜇(𝑛), 𝐾𝑡(𝑛) and 𝐾𝛾(𝑛). Thus, assuming that the coefficients 𝜇(𝑛), 𝑡(𝑛) and 𝛾𝑡(𝑛) change slowly, their variation ratios can be considered to lie around the unity. Therefore, the measurement vector is, 𝑧(𝑛)=[𝑟𝐷(𝑛)𝑟𝑃(𝑛)𝑟𝑅(𝑛)1 1 1]𝑡 and the measurement equation, 𝑧(𝑛)= ℎ[𝑥(𝑛)]+𝑣(𝑛) where 𝑣 = [𝑣𝐷𝑣𝑃𝑣𝑅𝑣𝜇𝑣𝑡𝑣𝛾]𝑡 is a Gaussian noise vector that models the measurement error with a covariance matrix R(n), and the nonlinear measurement function is ℎ(𝑥)=[𝑟 ∙𝐾𝜇𝑟 ∙ 𝐾𝑡𝑟 ∙𝐾𝛾𝐾𝜇𝐾𝑡𝐾𝛾]𝑡 With this formulation, the EKF is able to deal properly with the non-linearity arising in the measurement equation, performing at each iteration a linear prediction step, followed by the non-linear correction step. So, the KF provides the sequence of states estimates 𝑥(𝑛), whose last component, the estimate of 𝑟(𝑛), will be denoted as 𝑟𝐾(𝑛). The sequence 𝑟𝐾(𝑛) incorporates, in a statistically optimal fashion, the noisy information provided by raw data ratios, such as 𝑟𝐷(𝑛), 𝑟𝑃(𝑛) or 𝑟𝑅(𝑛), which results in a more reliable estimation than the raw data ratios themselves. After several initial trials, an enhanced model has been finally implemented, considering two improvements: 1) The random walk in 𝑟(𝑛) can be advantageously replaced, while keeping the linearity of the state equation, by a linear prediction based on the recent past history, e.g., based on 𝑟(𝑛 − 1), 𝑟(𝑛 −2) and 𝑟(𝑛 −3). In this case, 𝑥 = [𝐾𝜇𝐾𝑡𝐾𝛾𝑟 𝑟𝑛−1 𝑟𝑛−2 𝑟𝑛−3]𝑡 and 𝑥(𝑛)= 𝐴 ∙ 𝑥(𝑛 − 1)+ 𝑤(𝑛) (12) where: 𝐴 = [ 1 1 1 4 3 −1 3 1 1 1 ] , 𝑤 = [ 𝑤𝜇 𝑤𝑡 𝑤𝛾 𝑤𝑟 0 0 0 ] 2) The information available at instant 𝑛 can be used to improve previous estimations. For this purpose, the Rauch- Tung-Striebel (RTS) smoother [9], is implemented in two steps: • Forward pass performed with EKF. • Backward recursion smoother based on the linear state JBHI-00622-2020.R2 5 transition equation. Regarding the tuning of the filter, the diagonal elements of 𝑅(𝑛) related to 𝑣𝐷, 𝑣𝑃 and 𝑣𝑅 are self-tuned according to the respective sample variances of the last available days of 𝑟𝐷(𝑛), 𝑟𝑃(𝑛) and 𝑟𝑅(𝑛). The terms related to the pseudomeasurements are set to 10−4 to increase their weights. The diagonal elements of 𝑄(𝑛) are set to constant values 𝑄𝑖,𝑖 = 10−3, according to the observed variance of 𝑟𝐾(𝑛)− 𝑟𝐾(𝑛 −1), except in the case of history variables, where 𝑄𝑖,𝑖 = 0. The initial values of 𝑥(0) and P(0), required by the filter, are not relevant as their effect quickly vanishes. The proposed values are: 𝑟𝐾(0) and history variables initialized to the mean of the raw ratios, 𝐾𝑖 gains initialized to 1, and P(0) initialized to an identity matrix. Additionally, the results of the KF allow the estimation or prediction of two magnitudes of great interest: 1) The day in the future when the infection peak will be reached, 𝑛𝑝, which will occur when 𝑟𝑘(𝑛𝑝) = 1. For this purpose, the sequence 𝑟𝐾(𝑛) can be fitted to a predetermined evolution in order to predict its future behavior from the past history. According to what can be observed empirically, we have fitted the sequence 𝑟𝐾(𝑛) to a decreasing exponential, characteristic of first-order systems. It is worth stressing that, at this point, we are talking about the peak of infected, the peak of positives being a proxy for it. 2) The estimated number of active infectious people, 𝐼𝐾(𝑛). This magnitude can be computed from those infected at a given reference day, 𝑛0, as follows: 𝐼𝐾(𝑛)=∏𝑟𝑘(𝑖) 𝑛−1 𝑖=𝑛0∙𝐼(𝑛0) (13) This means that, besides the errors of the sequence 𝑟𝐾(𝑛), the uncertainty of 𝐼𝐾(𝑛) inherits that of the initial guess 𝐼(𝑛0). Therefore, as happens also with the SIR model, and in fact with any other model based on differential equations, precisely estimating 𝐼𝐾(𝑛) based on the above expression requires that an accurate number of infectious people be known on a given day, which can be very challenging. For our results, we will obtain the initial guess from 𝐼𝐾(𝑛0)= 𝑃(𝑛0)/𝑡(𝑛0), for an assumed value of 𝑡(𝑛0). For instance, 𝑡(𝑛0) = 0.2 if we believe the first day there were 5 infectious people per reported positive. IV. SIMULATION RESULTS In this section, the performance of the proposed implementation of the KF is tested on a set of simulated scenarios, where the SIR model described in Section II is considered for the propagation of a virtual virus. While the total simulation time is 90 days, a lockdown is assumed to take place at day 15. This restrictive policy is aimed to completely stop non-essential public mobility, resulting in a quick reduction of the transmission rate, β (which is assumed to evolve exponentially, from an initial value 0.5 to 0.09) and, therefore, also of the ratio 𝑟(𝑛). The remaining simulation parameters are given in Table I. In all simulated cases presented in the sequel, artificial Gaussian noise has been added to the measurements used by the KF algorithm in order to represent a more realistic scenario where the reported information presents inaccuracies. TABLE I VALUE OF THE PARAMETERS IN THE SIMULATION Parameter Definition Simulation value μ Mortality ratio 0.004 γ Recovery rate 0.1 𝐼(0) Initial infectious 1000 𝑁 Total population 47 Million 𝑡(𝑛) Testing ratio 0.2 A. Base case In the base case, the testing ratio, 𝑡(𝑛), is assumed to remain constant (except for the noise). The time evolution of the estimated geometric ratio, 𝑟(𝑛), is represented in Fig. 1 along with the raw noisy measurements provided by the simulation and the actual value of 𝑟(𝑛). Note that the estimation provided by the KF is very close to the simulated value, giving evidence of the good performance of the proposed method. Fig. 2 shows the benefit attained from the incorporation of the smoothing filter mentioned above to the basic EKF algorithm, in terms of more damped oscillations. Fig. 1. Estimation of r(n) in the base case To compare the proposed KF implementation with other methods customarily employed as filters in these cases, Fig. 3 shows the estimated value of 𝑟(𝑛) along with the results provided by three moving-average filters respectively applied to each of the noisy measurements, 𝑟𝑃, 𝑟𝐷 and 𝑟𝑅. As can be seen, the KF more closely tracks the evolution of 𝑟(𝑛). From the estimated 𝑟(𝑛) sequence, and the initial testing ratio, 𝑡(0), an estimation is obtained for the evolution of infectious people, which is compared in Fig. 4 with the simulated value of 𝐼(𝑛) and the reported positives. The maximum estimation error JBHI-00622-2020.R2 6 (around 4.5%) takes place, as expected, at the peak of the epidemic. Fig. 2. Estimation of r(n) with and without the smoother Fig. 3. Comparison of the proposed method with moving average filters Fig. 4. Estimation of the infectious people in the base case B. Error assessment Once the proposed estimation technique is validated in the base case, where only Gaussian noise is considered, the effect of different error sources is studied in the following scenarios. • Step in 𝑡(𝑛) An abrupt change is simulated in the testing ratio from 𝑡(𝑛)=0.2 to 0.3 at day 25, representing an increase in the availability of the tests (this has been observed in practice in several countries). Fig. 5 shows the estimation of 𝑟(𝑛), along with the simulated value and the measurements of the geometric ratios 𝑟𝑃, 𝑟𝐷 and 𝑟𝑅. Note that the step in 𝑡(𝑛) is observed as an impulse in the ratio 𝑟𝑃, which is quite effectively filtered out by the proposed KF implementation. Fig. 5. Estimation of r(n) with a step on the testing ratio Fig. 6. Estimation of the infectious people with a step on the testing ratio Step on 𝑡(𝑛) seen as an impulse on 𝑟𝑃(𝑛) Step on 𝒕(𝒏) day 25 JBHI-00622-2020.R2 7 The estimation of 𝐼(𝑛) is represented in Fig. 6, where the actual simulated value is again very close to the estimated one, giving evidence of the good performance of the method in the presence of a step in the testing ratio. • Deviations in 𝑡(0) Finally, the last scenario considered in this section shows how errors in the initial guess of the testing ratio, 𝑡(0), with respect to the assumed value 𝑡(0)= 0.2, affect the results. Given that the errors in this factor only affect the estimation of the infectious people, 𝐼(𝑛), the representation of the estimated 𝑟(𝑛) is not repeated. Fig. 7 shows the estimation of 𝐼(𝑛) for 𝑡(0)= 0.15 and 𝑡(0)= 0.25 (±25% error). Fig. 7. Estimation of the infectious people for a range of t(0) values. The results in Fig. 7 clearly show the importance of having an accurate guess of the number of infectious at the onset of the outbreak, as this initial error propagates proportionally up to the peak. Note, however, that any epidemiological model, such as SIR, faces the same challenge. V. CASE STUDIES In this section, the proposed KF-based estimation technique is applied to the real data reported by different countries. For convenience, the results have been divided in two subsections: 1) the time period when massive lockdowns occurred in most countries, denoted in the media as the “first wave” of the pandemic [22], and 2) the subsequent transient period, once the lockdowns are relaxed, usually through several de-escalation phases, towards the so-called “new normality”. 1. Lockdown period (first wave) A total of four countries have been considered in this period: China, South Korea, Spain and the UK. At the early stage of the pandemic, the information provided by these countries was sufficient to allow the application of the proposed methodology. Figs. 8-15 represent the estimated sequence of the geometric ratio, 𝑟(𝑛), and the number of infected people, 𝐼(𝑛), for the four countries. The KF implementation is tuned as described in Section III for the covariance matrices 𝑄 and 𝑅, and the initial values of the vector 𝑥(0) and the covariance matrix 𝑃(0). In order to estimate the number of infectious people, 𝐼(𝑛), according to (13), an initial value for the parameter 𝑡(0) is needed. In absence of a better clue, 𝑡(0) = 0.2 is considered in all cases, expect for Spain (see the discussion of this particular case below). The points for which 𝑟(𝑛)= 1 (peak of the epidemic) are highlighted with a dot. The following remarks can be made from those results: o A different evolution of 𝑟(𝑛) can be observed for the Asian countries (China and South Korea), where the effects of Covid-19 started earlier. Once the geometric ratio 𝑟(𝑛) < 1, the trend for South Korea is to remain roughly constant throughout the considered period, whereas for China a certain rebound can be noticed after March 10. o The estimation results obtained for Spain show an asymptotic trend towards 𝑟(𝑛)=0.95. A slight increase is observed in 𝑟(𝑛) between April 10 and 15, probably influenced by a sudden increase in the number of tests. Regarding the parameter 𝑡(0), in the Spanish case we have taken into account the results of a massive seroprevalence test performed by the government in the first half of May [23], from which it was concluded that the total number of infected people was around 5.2% of the population (approximately 2.3 million people). In view of this valuable information, the initial value 𝑡(0) has been adjusted so that the cumulative number of infected people matches the result of the survey on the date it was released (May 13), leading to 𝑡(0)=0.12. This provides the estimation of 𝐼(𝑛) shown in Fig. 13, where a maximum value of the active infectious people of around 1.3 million can be noticed by mid-April. Fig. 16 represents the estimation of the cumulative infectious people for the Spanish territory, where the total number of infected people matches the results of the survey. o With the available information in mid-May, some countries had already left behind the peak of the epidemic (i.e., 𝑟(𝑛) < 1). For those cases (China, South-Korea and Spain), a rather accurate early forecasting of the epidemic evolution can be made, around 10 to 14 days before the peak, by fitting a decreasing exponential to a window of past estimated data. This prediction is shown with green dotted lines in Figs. 8- 13. o Regarding the UK, where the peak of the number of infectious people had not been reached in the period considered (i.e., 𝑟(𝑛) > 1), an exponential fitting (made between around mid-April and early May) and the corresponding extrapolation is considered for this case. According to such fitting, the peak should take place in the second half of May. JBHI-00622-2020.R2 8 Fig. 8. Estimation of r(n) in China in the first period considered Fig. 10. Estimation of r(n) in South-Korea in the first period considered Fig. 12. Estimation of r(n) in Spain in the first period considered Fig. 9. Estimation of infectious people in China in the first period considered Fig. 11. Estimation of infectious people in South-Korea in the first period Fig. 13. Estimation of infectious people in Spain in the first period considered JBHI-00622-2020.R2 9 Fig. 14. Estimation of r(n) in UK in the first period considered Fig 16. Estimation of cumulative infectious people in Spain Fig 17. Fitted geometric ratios from common threshold Fig. 15. Estimation of infectious people in UK in the first period considered Finally, Fig. 17 shows the evolution of the above-mentioned fitted exponential curves for different countries (Italy has been included in the representation in order to establish a more complete comparison), all of them represented from a common threshold 𝑟(𝑛)= 1.2, so that the corresponding time constants can be easily compared. In light of this representation, it can be noticed that the reduction of the geometric ratio is faster in China (just 13 days from 𝑟(𝑛)= 1.2 to 𝑟(𝑛)= 1), possibly as a consequence of a more severe lockdown, followed by Spain and South Korea (between 25 and 27 days to reach 𝑟(𝑛)= 1), showing similar trends, and finally Italy (40 days to reach 𝑟(𝑛)= 1). 2. Post-lockdown period (second wave) As the pandemic evolves, it becomes more difficult to properly report on a regular basis all the information involved in the estimation of active positives. Many countries (notably Spain) stopped reporting the number of recovered people, probably owing to the remarkable increase in the number of asymptomatic positive cases, which never entered a hospital and hence never counted as recovered or dead. For this reason, it is not possible to accurately update the estimations of the geometric ratios of active positives, 𝑟(𝑛), for some of the countries considered in the early stages. Instead, Figs. 18-25 represent the estimated geometric ratio, 𝑟(𝑛), and the number of active infectious people, 𝐼(𝑛), for four countries (USA, Italy, India and Brazil), all of them specially affected by the pandemic and still reporting the information required by the proposed estimation technique. JBHI-00622-2020.R2 10 Fig. 18. Estimation of r(n) in USA in the second period considered Fig. 20. Estimation of r(n) in Italy in the second period considered Fig. 22. Estimation of r(n) in India in the second period considered Fig. 19. Estimation of infectious people in USA in the second period considered Fig. 21. Estimation of infectious people in Italy in the second period considered Fig. 23. Estimation of infectious people in India in the second period considered