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 Taylor-Fourier-trigonometric 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 time-varying 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 short-circuit 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 non-DFT-based 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 Gauss-Newton 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 band-pass filter in order to evaluate the phasor. On the other hand, the authors of [9] introduce a special angle-shifted 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. Taylor-Fourier algorithm [15] approximates the dynamic phase with the second-order Taylor expansion and the least square observer. The Taylor-Kalman 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 (Taylor-Fourier transform) harmonic analysis uses recursive multiple-resonator-based 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 Taylor-Kalman Filter for instantaneous dynamic phasor estimation is developed which is coincident with the self-adaptive 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, Rate-Of-Change-Of-Frequency, 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 suboptimal-scaling factor-adaptive strong-tracking 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 square-based Prony analysis and Taylor expansion called Taylor–Prony is proposed to estimate the dynamic phasor.
In this paper, a special transformation (representation) of Taylor-Fourier 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 Taylor-Fourier-trigonometric 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 well-known 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 round-off errors in the case of large number of samples).
The accuracy of the proposed method is compared with some well-known 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 zero-order 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 well-known 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 (1654-1705) 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
The presented representation of Gram’s matrix shows that the computation of pseudoinverse can be performed more rapidly and with more precision.
4. Simulation Results
The proposed phasor measurement scheme is employed to estimate coefficients of the power signal under nonstationary scenarios, during which we compared the results obtained in the proposed algorithm with some other well-known procedures. Practically, through the simulation test we estimated the ability of the algorithms to track harmonic changes. Tests have been performed to show the performance of the proposed technique with reference to the conditions of the standard [28] for synchronized phasor measurement systems in power systems. The first is a proposed estimation procedure verified on the basis of the following test signal, which is characteristic for oscillating conditions:where , , , , and .
As can be seen from the above relations, the test signal is sampled with a frequency of 1200 Hz, thus forming 24 samples within a 20 ms wide window, which is the period of the signal with a fundamental frequency of 50 Hz. Figures 1(a) and 1(b) show the size of the estimated amplitude and phase of the fundamental dynamic phasor, using the proposed estimation method. It is clear that the main difference between the proposed method and some other methods based on the application of the least square methods [15, 21, 27], and the methods based on the application of Kalman filters [16, 17, 23, 31] is in the occurrence of the delay in the estimation due to the use of the data window.
(a)
(b)
One of the established criteria for assessing the quality of the estimation is based on the Total Vector Error (TVE), where we can estimate the error size in the estimation of the phasor magnitude and angle, defining it aswhere and are real and estimated values. Figure 2 shows the total vector error of the proposed method in the first ten cycles; one cycle time delay of estimated phasor is compensated as in [27]. It can be concluded that the error occurring is above all the results of the introduced delay, similar to that of [15, 21]. Owing to the modification made in calculating the parameters of the phasors, the magnitudes of this calculated error are smaller than in other methods based on the least squares method, with their size becoming comparable to [16, 17, 23, 38], therefore being much smaller than [27, 39–41]. The proposed approach can calculate the derivatives of the phasor-phasor speed and acceleration, which reduces the estimation error. In addition, it is possible to calculate the frequency of the processed signal and detect faults and power swings.
4.1. Amplitude Oscillation Case
The ability of the proposed approach to track an arbitrary oscillating amplitude is given in the following example whereis the input signal to the estimator, with Hz and nonzero Fourier coefficients (, , and ), and where represents an amplitude oscillation, given by the following second-order polynomialwith coefficients , , and . The corresponding digital signal is sampled at 24 samples/cycle (Hz). The oscillation of main-original signal and its corresponding estimated signal components are shown in Figure 3. Based on this representation, it can be clearly seen that the proposed method for estimation is fully conducive to following the dynamics of all components of the processed input signal. Also, we can see that the quality of estimation decreases with increasing the number of harmonic, which are the results of relatively small and fixed-constant number of samples. For higher harmonics components this number of samples cannot provide precise reconstruction. If the number of samples adaptively changed as the function of the order of harmonics, the accuracy of the estimation would be the same as on the first-fundamental harmonic component. An additional problem is certainly the nonideality of the filter function.
(a)
(b)
(c)
(d)
4.2. Magnitude-Phase Step Test
In this section the amplitude and the phase step tests are performed using the proposed estimation procedure. The test signal is defined as [28] (in the standard, changes in this step are happening separately)where is a unit step function, is the amplitude step size, is the phase step size, and is the step instant. In this test, a signal with the 10° phase angle step and the amplitude step changes from 0.9 to 1 p.u. in the fundamental component which is fed to the algorithm. The parameters are chosen as , , , while the measurement noise is set to be negligible. The results for the amplitude and phase step tests are shown in Figure 4, in which the upper subplot shows the transient response of the amplitude or the phase, while in the lower plot the TVE are illustrated.
(a)
(b)
(c)
It is observed that the proposed approach achieves relatively long transient periods but with small overshoots. Based on the results presented in Figure 4, the proposed method is capable of performing a precise estimation of the amplitude and phase after the expiry of the transition period (the model is more appropriate for smooth amplitude and phase changes), since it represents the modification of Tailor’s-Fourier expansion. The duration of this transition in the proposed approach is shorter than in some other methods that have the same starting point [17, 27], which is the result of the position of the poles of the transfer function with respect to the unit circle in the plane (as the poles are more spaced). In addition, there is no large overshoot in the response, as opposed to [16, 17, 23] due to the amplitude and phase discontinuities at the step change, since no previous assessment is performed concerning a sample window.
In Table 9 the performance indices defined in the standard [28] including response time, delay time, and overshoot are listed for both amplitude and phase step change tests. The performances achieved by the proposed algorithm fulfill the M-Class measurement requirements in the standard.
4.3. Frequency Step Test
During this simulation test, a signal having a 5 Hz frequency step is used to evaluate the possibility of a transient response in a frequency step condition (although 5 Hz frequency step is not likely to happen in a power grid, this test is designed based on IEEE standard [28]).
The frequency was calculated as a ratio of roots of characteristic equation based on extracted unknown coefficients , , , , , from (12) and [15, 27]. The proposed method shows characteristics, as well as the evaluation process, similar to the previous type of signal (magnitude-phase step condition). Figure 5 shows that the proposed processing method can be successfully tracked as +5 Hz frequency step, versus the traditional approach to processing the phasors.
4.4. Harmonic Infiltration Test
Input signal defined with (66) was used to validate the proposed methods in harmonic condition:
Figure 6 represents the output of the proposed method in the harmonic condition. According to Figure 6, fundamental phasor estimation based on proposed modification is free from harmonic, which shows the superiority of harmonic modification of the proposed method in this paper. Since in this test, it is possible to estimate the dynamic phasor of the first harmonics. To increase the range of harmonics, the sampling number per cycle should be increased.
4.5. Frequency Response
In order to evaluate the frequency characteristic of the proposed estimator, its transfer function is defined in the formwhere and are transforms of the input (main signal) and output (estimated dynamic phasor). The frequency response of the proposed method is given in Figure 7. The magnitude responses are unsymmetrical, which is completely understandable since it belongs to complex filters. In the surroundings of the fundamental frequency, complete elimination of the negative one and complete pass of the positive one occur. The flatter transfer characteristic causes a lower distortion in the estimation of the phasor and provides zero-gain in the nonfundamental component. In this way, the harmonics at the output are removed, as opposed to [16, 21, 23].
(a)
(b)
4.6. Error Bounds
In order to determine the boundary of the possible error in the estimation of the phasor, signals with variable envelopes were used:where the time constants (τi) were generated by a uniform random process in the interval of cycles. In a similar way, the three frequencies were randomly generated in the intervals of , , and Hz. The error is calculated bywhere and are real and estimated amplitudes, respectively. Figure 8 shows the histograms of the errors attained by proposed method.
The RMS error belongs to the range of 3×10−3 and 8×10−3 for proposed estimation method. This error is related to elimination of the third derivative term of the Taylor expansion in phasor estimation process, and it is much better than error estimated in [16, 17, 22–24, 27].
4.7. Noise, Subharmonics and Interharmonics Infiltration
It is useful to analyze the estimation error of the proposed method when the noise level changes. A criterion named transient monitor (TM) is used for this evaluation [21]. TM is an index to monitor a sudden change of the measured signal and is calculated bywhere is the real (measured) value of sample and is the recomputed sample obtained by the estimated dynamic phasor. Estimation error , is the difference between these two quantities. TM is calculated by estimation error in every sample and can be used as a quality measure of phasor estimation. Figure 9 shows the TM as a function of noise variance, for the proposed and some well-known methods for phasors estimation. As the noise variance increases, higher errors are obtained, so there is an upward trend when noise level increases. According to Figure 9, there is a stable trend for least squares-based methods and our proposed method, before critical point with variance of about 10−2. After this point, transient monitor increases steeply by increasing error variance. Generally Kalman based methods [16, 17, 23] show lower values compared to least square-based methods; however, our proposed algorithm has the best performance relative to other Taylor expansion based methods, when noise variance is low, according to Figure 9. Also we can conclude that the proposed method depicts a downward trend by increasing the sampling frequency because the estimation accuracy increases with a high sampling rate, and contrary to [21], does not increase its noise sensitivity at a high sampling rate.
The standard [28] does not explicitly define the behaviour of estimation methods in the presence of interharmonic and subharmonic components, which can be found in the installations in which these measurements are applied. In general, interharmonics can be very difficult to detect and isolate if they are located in the band of interest of the phasor dynamic model. In order to obtain an even more complete picture of the performance of the phasor estimation method described here, a subharmonic frequency of 20 Hz, as well as two interharmonic components frequency of 57.5 Hz and 85 Hz [5], with a 10% amplitude with respect to the fundamental component amplitude, has been added to a sinusoidal signal with the aforementioned contents. An out-of-band interference test is defined to test the capability to filter interfering signals that could be aliased into the measurement, and it is considered only for M class and it is eliminated for P class [39].
The precision in estimation of the proposed algorithms is directly related to their filtering capabilities. When the interharmonic is in the passband of the algorithm (designed to include dynamic phenomena of interest), as for interharmonic components on 57.5 Hz, the TVE is about 3%, whereas for interharmonics in the stopband frequencies, TVE depends on the stopband attenuation. In this case, for an interharmonic frequency of 85 Hz, TVE is about 0.5 %. The same situation is for subharmonics components on 20 Hz. The obtained results are better than the results obtained in [5, 39].
The proposed approach in estimation of phasor values was also tested in a situation when processing the input signal contains the third harmonic with amplitude settled to 20% of the fundamental component, and an interharmonic frequency 75 Hz. This level of the harmonics is very commonly present in current waveforms in the control applications such as a current tracking in shunt active power filters (APFs) in marine networks or smart distributed grids with renewable energy sources [39]. TVE in estimation of third harmonics in this test was about 1%. In all above-mentioned tests the number of samples was N=24.
4.8. Computational Complexity and Simulation Time
The inputs of the algorithm are matrix of the type and matrix of the type given by (16) and (13), respectively. The output of the algorithm is column matrix of the type obtained by (19). The proposed algorithm can be decomposed in the following steps:(A1)Calculation of the elements of Gram’s matrix by (45)-(59).(A2)Inversion of Gram’s matrix .(A3)Multiplication of with .(A4)Multiplication of with .
(A1) If the formulas (45)-(59) are used for calculation of the elements of Gram's matrix , then the task is reduced to consecutive calculations of specific values based on polynomials. If denotes the maximum of all computational times needed for executing (45)-(59), then the computational time for does not exceed .
(A2) The solution to (19) includes the inversion of the Gram’s matrix. To find the inverse of -matrix one can use Lower-Upper (LU) decomposition method which takes FLOPs.
(A3) In this step, the matrix of type is to be multiplied by the matrix of the type which takes the FLOPs.
(A4) In this step, the matrix of the type is to be multiplied by the matrix of the type which takes the FLOPs.
Based on the all above-mentioned procedures-steps, the proposed algorithm does not exceed operations-FLOPs, which practically defined its computational costs. If the step (A1) is replaced with
(A’1) Multiplication of with ,
then the computational time of this step is increased with respect to the time needed for performing (A1). Since is -matrix, it follows that is -matrix and multiplication of matrices gives an algorithm that takes time to multiply the matrix of type with the matrix of type . Clearly, the time needed for such an algorithm depends on the number of samples . On the other hand, the time in (A1) is independent of , so the number of samples does not influence the speed of (A1). This enables us to deal with a large number of samples, without decreasing the performance of algorithm. Also, if the time is the maximum of all computational times needed for executing (45)-(59), clearly, it is the time needed for realization of (54). If is regarded as , then its computational time takes 3 FLOPs for basic arithmetic operations. Consequently, the computational time for (54) can be approximated with finite number which is the number of basic arithmetic operations (addition, subtraction, multiplication, and division) in (54).
In order to practically (in real conditions) determine the time needed for the computation of unknown phasors parameters according to the proposed method, the frequency of the processed signal is increased, which can clearly define whether the method can be used for offline or online processing. A comparison was made with well-known methods in the literature and the resulting comparisons were given in Figure 10. The hardware features of the following characteristics were used: Intel(R) Core(TM) i5-4460, 3.2GHz, 16 GB RAM. The proposed modification in the algorithm described here showed that the time required for the calculation was significantly less than in other methods based on the least squares method [24, 27], and approaching the speed by methods [22, 23, 39–41].
5. Conclusions
In this paper we used the dynamic phasor concept to estimate the variable amplitude and phase in processing of oscillating signals in modern power grids. The specific modification of Taylor-Fourier expansion and corresponding Gram matrix give us the possibility to develop more precise and computational attractive algorithm for estimation of all unknown parameters of processing signal, in all conditions defined with standard. Simulation results demonstrate the accuracy and capability of the proposed method in view of TVE, harmonic condition, frequency tracking, and computation time. The proposed method has been investigated under different conditions and found to be a valuable and efficient tool for detection of signal components under dynamic conditions. The simulations results have shown that the proposed technique provides accurate estimates and offers the possibility to track the phase, frequency, and amplitude changes of nonstationary signals. Due to its specific form, which enables the reduction of both the computational cost and the memory requirements of the algorithm the proposed algorithm can be very useful for real-time digital systems.
Data Availability
The numerical data, obtained analytical relations, and presented graphics data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research is supported by the Ministry of Education, Science and Technological Development, Republic of Serbia, Grant Nos. 42009, 172057, and 174013.