Research Article  Open Access
Dynamic Phasors Estimation Based on TaylorFourier Expansion and Gram Matrix Representation
Abstract
The paper presents a new approach to estimation of the dynamic power phasors parameters. The observed system is modelled in algebra of matrices related to its TaylorFouriertrigonometric series representation. The proposed algorithm for determination of the unknown phasors parameters is based on the analytical expressions for elements of the Gram’s matrix associated with this system. The numerical complexity and algorithm time are determined and it is shown that new strategy for calculation of Gram’s matrix increases the accuracy of estimation, as well as the speed of the algorithm with respect to the classical way of introducing the Gram’s matrix. Several simulation examples of power system signals with a timevarying amplitude and phase parameters are given by which the robustness and accuracy of the new algorithm are confirmed.
1. Introductions
The very intense and rapid development of electronic technologies, including the production of renewable energy sources (commercial solar, wind power plants, and biomass power plants) and power devices, caused a change in the structure of the traditional power system. Electrical variables, such as the basic harmonics amplitude and its frequency, can be significantly altered in micronetworks, weak networks, and networks of the island type due to reduced capacity in terms of shortcircuit currents. In addition, the harmonic is injected into the system using electrical electronic equipment [1]. Under these circumstances, an ideal algorithm for assessing and estimation of the pharos must enable quick, accurate, and steady monitoring of the changing parameters of the electrical signals that are contaminated by the present harmonic components.
In the last few years, many algorithms have been proposed in the literature for evaluating the phases in different dynamic states and conditions that can be controlled in the network. In general, they can be classified into methods based on the application of discrete Fourier transform (DFT) [2–4] and nonDFTbased methods. Each of the algorithms requires a harmonic/pharos model and uses some of the specific techniques that model the parameters of the observed system. In accordance with the modelling method, the measurement method, the algorithms can be divided into two main classes: algorithms that rely on the pure sinusoidal signal model (static model) and algorithms based on a nonsinusoidal model (dynamic model) [5].
A simultaneous estimation of phasor and frequency based on the application of the fast recursive GaussNewton algorithm was proposed in [6], while the method based on the modified Fourier transformation in order to eliminate the DC offset was presented in [7]. The paper [8] proposes a new method based on an adaptive bandpass filter in order to evaluate the phasor. On the other hand, the authors of [9] introduce a special angleshifted energy operator to separate the value of the instantaneous amplitude of the phasor. The least square curve fitting approach was proposed in [10] in order to eliminate the unwanted effect of saturation of the current transformer. Measurement of the phase and frequency during the transient processes was considered in [11], while the dynamic estimation of the phase based on the application of maximally flat differentiators was proposed in [12] and the phasorlet in [13]. The paper [14] introduces an approach to estimate the parameters of the phasor based on recursive wavelet transformation. TaylorFourier algorithm [15] approximates the dynamic phase with the secondorder Taylor expansion and the least square observer. The TaylorKalman method based on the Kalman observer in order to achieve nondelayed dynamic phase estimation is described in [16]. This method is further enhanced by the development of the state space for harmonic infiltration in [17] and is called Tailor Fourier Kalman. The Prony algorithm [18] was used to estimate the phasor described by the exponential amplitude and the linear phase, with no Taylor expansion. There are other methods that are modifications of earlier methods [19, 20]. The Shank method is another method for estimating a dynamic phasor, using the least square and consecutive delays of the unit response system [21]. The algorithm [22] for TFT (TaylorFourier transform) harmonic analysis uses recursive multipleresonatorbased computational techniques, which enable the reduction of both the computational cost and the memory requirements of the algorithm.
In [23] a more precise dynamic model describing the complex trajectory of the dynamic phasor is derived owing to a revised state transition equation of the rotating phasor and its derivatives. Based on the improved dynamic model [23], a modified TaylorKalman Filter for instantaneous dynamic phasor estimation is developed which is coincident with the selfadaptive nature of the Kalman filter principle. An extension of the Taylor Weighted Least Squares (TWLS) algorithm for the estimation of the phasor, frequency, and ROCOF, RateOfChangeOfFrequency, parameters of an electric waveform has been analyzed in [24]. Paper [25] introduces a dynamic phasor estimation method for PMUs/IEDs based on a modified time domain hybrid method. Dynamic phasor is estimated using the Taylor expansion, where frequency deviation is derived directly from fitting parameters to avoid magnification of fitting errors. A double suboptimalscaling factoradaptive strongtracking Kalman filter (DSTKF)based phasor measurement unit algorithm which can meet the accuracy requirement of the IEEE standard C37.118.1 under the dynamic condition was proposed in [26]. This method uses a kth Taylor polynomial to linearize the complex exponential of the signal model and estimates the dynamic phasor using DSTKF. In [27], a combination of the least squarebased Prony analysis and Taylor expansion called Taylor–Prony is proposed to estimate the dynamic phasor.
In this paper, a special transformation (representation) of TaylorFourier expansion and corresponding Gram matrix is used to estimate the dynamic phasors from input phasorlets in real domain, based on synchronously sampled values of the input signal. The dynamic system’s behaviour is described in algebra of matrices related to their TaylorFouriertrigonometric series representation. The explicit formulas for calculation of elements of the associated Gram’s matrix are given. In this way, we are in a position to improve the computational performance of the proposed algorithm, as regards precision and speed. Namely, the proposed algorithm is based on matrix multiplication. It is a wellknown fact that mathematical definition of multiplication of matrices gives an algorithm that takes time to multiply a matrix of type with a matrix of type (in the case of square matrices, it takes time). If we use the explicit formulas for elements of Gram’s matrix in our algorithm, then one matrix multiplication is avoided, and the task is reduced to consecutive calculations of the values of specific polynomials. Moreover, each element of Gram’s matrix depends on the number of samples (it is expressed in terms of a polynomial of variable ), but the time necessary for this calculation depends on the degree of that polynomial. Consequently, the time needed for calculation of Gram’s matrix remains the same even if the number of samples is increased. Another important reason for using analytical expressions for elements of Gram’s matrix instead of their interpretation through the usual matrix multiplication can be found in the fact that rounding of exact numbers is almost unavoidable when reporting many computations and these rounding errors generally accumulate (explicit formulas for elements of Gram’s matrix induce smaller roundoff errors in the case of large number of samples).
The accuracy of the proposed method is compared with some wellknown methods from literature, illustrating the capability of tracking dynamic phasors, especially in comparison with other procedures based on the least squares method. Sinusoidal and step changes of the amplitude and phase, harmonic condition, frequency tracking test, and computation time are different tests which are used to validate the proposed method, according to the definition and test cases in the standard [28]. The potential of the proposed approach is demonstrated by simulating various numerical signals in MATLAB. The proposed algorithm is particularly suitable for the integration of distributed generating sources with microgrids, when fast detection of faults and the islanding condition is required.
2. Dynamic Signal Model and Algorithm Description
The behaviour of a power system under oscillation is usually modelled by trigonometric series where and are the amplitude and the phase functions which represent variations of the mth dynamic harmonic over time, and ( is frequency of fundamental component). In an equivalent form, the behaviour of such power system can be represented aswhere
Clearly, the following holdsand
Let and be approximated with two Taylor series within a short period of time near the reference time , i.e.,and
Therefore, we have
In a discrete version, this can be written as the system of equationswhere .
For each , the mth harmonic representation is
In order to perform the necessary and proposed calculations for the estimation of unknown phasors parameters from input phasorlet, samples of the processed signal are first transformed through the Fourier nonrecursive algorithm [29, 30]. After that they are introduced in the form of input parameters to the algorithm proposed in this paper. The mathematical model of the algorithm for digital filtering is defined aswhere is nth filtered sample of the mth harmonic voltage or current signal, is the processed signal–sample in the time moment , is the sampling period, and is the number of samples, . After the extraction we obtain the vector of the filtered samples, the length of , on which to apply the proposed procedure. For the realization of this filter transformation, the FIR structure described in [31] can be used, but it was chosen (11) for the reason of obtaining a simpler and faster procedure for estimating the processed phasors.
Let us observe that for , we obtain the following system of equations
The traditional algorithm for estimating the value of the phasor is based on the subsystem of (12), in which only two column vectors are taken into account; in this way the dynamic phasor is approximated by Tailor's zeroorder polynomial over the interval in which observations are made. This generates a staircase function, with a variable step from one interval to another. Such a model is accurate only when the input signal is in the stationary state. This is certainly not enough in a situation where oscillations in the power system occur, in which the first and second derivatives are as relevant as the constant term. A matrix representation of this discrete system (12) has the formwhere And
Our goal is to minimize the sum of squares of the residuals i.e., to find the solution to the matrix equation
The matrix is called the Gram’s matrix. The best solution (in the least squares sense) exists provided that the Gram's matrix is invertible, which is fulfilled when the column vectors H are linearly independent [32]. The Grammian inversion depends on the size of the interval N and the order describing the model of the signal itself that is the subject of processing. In our approach we assume that , for some (in the case that , the number of samples N covers the whole period).
It is clear that phasors are given by the inverse transform of the phasorlets and that the best solution is given by
The pseudoinverse matrix depends only on the parameters of the adopted signal model. The estimation of the phasors in the center of the evaluation interval where Tailor's error is zero is correct if the input signal is well described by the adopted model for which the least mean squared error (LMS) is also zero. An LMS error would affect the estimate if the signal was outside the projection subspace of the LMS algorithm.
3. Gram’s Matrix of Dynamic Signal Model
In this section, we will investigate the Gram's matrix of the dynamic signal model (12), and for each element of that matrix, an explicit formula will be given.
The Gram’s matrix has the form where , and , and for all and , the elements , , and are finite trigonometric series
With the goal of calculating the elements , , and , we will consider finite series
Clearly, for the following holds
The series , , and can be represented equivalently as
It is a wellknown fact that an analytic solution for a sum of powers of positive integers is where is Riemann zeta function [33, 34], is Hurvitz zeta function [35], and is generalized harmonic number [36]. The Swiss mathematician Jacob Bernoulli (16541705) derived the formula for the finite sum of powers of consecutive positive integers [37] representing as polynomial in of degree . In the case that the power takes values in the set , the corresponding polynomials are presented in Table 1.
Now, let us consider the finite sums involving trigonometric functions
The following holds for their exponential formwhere
For arbitrary , let functions be defined in the following way:
Clearly, for hold and
Therefore, we have and implicitly we obtain
The series , , and can be represented, in the terms of and , in the following way:
Also, these series can be represented in the terms of , , and in the following simple notation:
In the case that takes values in the set , the explicit forms of functions obtained by consecutive derivations according to formula (35) are presented in Table 2.
In the main section of this article, the synchronous model is observed, i.e., the number is chosen with respect to the relation . For and , this leads to the values given in Table 3.

If we represent each in in its expanded form as a polynomial in and of degree (the degree of each term in this polynomial in two variables is the sum of the exponents in each term), then can be transformed into the form given in Table 4.

In a similar way, by representing each in in its expanded form as a polynomial in and of degree , can be transformed into the form given in Table 5.

Now, by formula (42), we obtain values , , given in Table 6.
Using formula (43) we obtain values , for , given in Table 7.
Finally, by (44) we obtain values , for , given in Table 8.

Now, by (27) we obtain the elements of Gram's matrix . Namely, for all and , the elements are given by
Further, for all and , the elements are given by
Finally, for all and , the elements are given by