# Dynamic Phasors Estimation Based on Taylor-Fourier Expansion and Gram Matrix Representation

Authors: Predrag Petrović; Nada Damljanović
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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, andROCOF, 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 akth 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 takespqr time to multiply a matrix of type p×q with a matrix of type q×r (in the case of square p×p-matrices, it takes p3 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 pr consecutive calculations of the values of specific polynomials. Moreover, each element of Gram’s matrix depends on the number of samples N (it is expressed in terms of a polynomial of variable N), 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(1)xt=∑m=0Mamtsin⁡mωt+ψmt,where am(t) and ψm(t) are the amplitude and the phase functions which represent variations of themth dynamic harmonic over time, and ω=2πf1 (f1 is frequency of fundamental component). In an equivalent form, the behaviour of such power system can be represented as(2)xt=∑m=0MAmtsin⁡mωt+Bmtcos⁡mωt,where (3)Amt=amtcos⁡ψmt,Bmt=amtsin⁡ψmt.Clearly, the following holds(4)ψmt=arctgBmtAmtand(5)amt=Am2t+Bm2t.LetAm(t) and Bm(t) be approximated with two Taylor series within a short period of time near the reference time tref=0, i.e.,(6)Amt=∑k=0KAmk0k!tk=∑k=0KAmktkand(7)Bmt=∑k=0KBmk0k!tk=∑k=0KBmktk.Therefore, we have(8)xt=∑m=0M∑k=0KAmktksin⁡mωt+∑k=0KBmktkcos⁡mωt.In a discrete version, this can be written as the system ofN equations(9)xnT=∑m=0M∑k=0KAmknkTksin⁡mωnT+∑k=0KBmknkTkcos⁡mωnT,where n=1,2,…,N.For eachm=1,2,…,M, themth harmonic representation is(10)xmnT=∑k=0KAmknkTksin⁡mωnT+∑k=0KBmknkTkcos⁡mωnT.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 as(11)xmfilnT=1N∑k=0N-1xk+nTcos⁡kmωT,n=1,2,…,N,where xmfil(n) isnth filtered sample of themth harmonic voltage or current signal, x((k+n)T) is the processed signal–sample in the time moment (k+n)T, T is the sampling period, and N is the number of samples, N=f/f1  (f=1/T). After the extraction we obtain the vector of the filtered samples, the length of N, 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 forK=2, we obtain the following system of equations(12)xmfilT=Am0+Am1T+Am2T2sin⁡mωT+Bm0+Bm1T+Bm2T2cos⁡mωT,xmfil2T=Am0+Am12T+Am222T2sin⁡mω2T+Bm0+Bm12T+Bm222T2cos⁡mω2T,xmfil3T=Am0+Am13T+Am232T2sin⁡mω3T+Bm0+Bm13T+Bm232T2cos⁡mω3T,⋮xmfilNT=Am0+Am1NT+Am2N2T2sin⁡mωNT+Bm0+Bm1NT+Bm2N2T2cos⁡mωNT.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 form(13)xmfil=H·ym,where (14)xmfil=xT,x2T,x3T…,xNTT,(15)ym=Am0,Am1,Am2,Bm0,Bm1,Bm2T,And(16)H=sin⁡mωTTsin⁡mωTT2sin⁡mωTcos⁡mωTTcos⁡mωTT2cos⁡mωTsin⁡mω2T2Tsin⁡mω2T22T2sin⁡mω2Tcos⁡mω2T2Tcos⁡mω2T22T2cos⁡mω2Tsin⁡mω3T3Tsin⁡mω3T32T2sin⁡mω3Tcos⁡mω3T3Tcos⁡mω3T32T2cos⁡mω3T⋮⋮⋮⋮⋮⋮sin⁡mωNTNTsin⁡mωNTN2T2sin⁡mωNTcos⁡mωNTNTcos⁡mωNTN2T2sin⁡mωNT.Our goal is to minimize the sum of squares of the residuals(17)Ry~m=xmfil-H·y~m,i.e., to find the solution to the matrix equation(18)HT·H·y~m=HT·xmfilThe matrixHT·H 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 vectorsH are linearly independent [32]. The Grammian inversion depends on the size of the intervalN and the order describing the model of the signal itself that is the subject of processing. In our approach we assume that 2mωNT=2sπ, for some s∈N (in the case that s=1, the number of samplesN 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(19)y~m=HT·H-1·HT·xmfil.The pseudoinverse matrix(HT·H)-1·HT 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 matrixHT·H has the form (20)HT·H=ACCB,where A=aiji,j=13, B=biji,j=13 and C=ciji,j=13, and for all i,j=1,2,3 and p=i+j-2, the elements aij, bij, and cij are finite trigonometric series(21)aij=∑n=1Nnp·sin2mωnT,(22)bij=∑n=1Nnp·cos2mωnT,(23)cij=∑n=1Nnp·sin⁡mωnT·cos⁡mωnT.With the goal of calculating the elementsaij, bij, and cij, we will consider finite series (24)ApN,x=∑n=1Nnp·sin2nx,(25)BpN,x=∑n=1Nnp·cos2nx,(26)CpN,x=∑n=1Nnp·sin⁡nx·cos⁡nx.Clearly, forx=mωT the following holds(27)aij=Tp·ApN,mωT,bij=Tp·BpN,mωT,cij=Tp·CpN,mωT.The seriesAp(N,x), Bp(N,x), and Cp(N,x) can be represented equivalently as (28)ApN,x=12·∑n=1Nnp-12·∑n=1Nnp·cos⁡2nx,(29)BpN,x=12·∑n=1Nnp+12·∑n=1Nnp·cos⁡2nx,(30)CpN,x=12·∑n=1Nnp·sin⁡2nx.It is a well-known fact that an analytic solution for a sum of powers of positive integers is(31)SpN=∑n=1Nnp=ζ-p-ζ-p,N+1=Hn-p,where ζ(z) is Riemann zeta function [33, 34], ζ(z,a) is Hurvitz zeta function [35], and Hn(k) 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 Sp(N) as polynomial in N of degree p+1. In the case that the power p takes values in the set {0,1,…,4}, the corresponding polynomials are presented in Table 1.Table 1
Sum of powers of positive integers.
S 0 ( N ) N  S 1 ( N ) 1 2 ⋅ N 2 + N   S 2 ( N ) 1 6 ⋅ 2 N 3 + 3 N 2 + N   S 3 ( N ) 1 4 ⋅ N 4 + 2 N 3 + N 2   S 4 ( N ) 1 30 ⋅ 6 N 5 + 15 N 4 + 10 N 3 - NNow, let us consider the finite sums involving trigonometric functions(32)RpN,x=∑n=1Nnp·cos⁡2nx,IpN,x=∑n=1Nnp·sin⁡2nx.The following holds for their exponential form(33)EpN,x=∑n=1Nnp·ei2nx=12ip·∑n=1Ndpdxpei2nx=12ip·dpdxp∑n=1Nei2nx=12ip·dpdxpei2x·ei2Nx-1ei2x-1=12ip·dpfxdxp,where(34)fx=ei2x·ei2Nx-1ei2x-1.For arbitraryN∈N, let functions gp(z):C→C be defined in the following way:(35)g0z=z·zN-1z-1,gpz=2i·z·dgp-1zdz,p=1,2,….Clearly, forz=ei2x hold f(x)=g0(z) and (36)dpfxdxp=ddxdp-1fxdxp-1=ddzdp-1fxdxp-1·dzdx=ddzgp-1z·dzdx=gpz.Therefore, we have(37)EpN,x=12ip·gpei2x,and implicitly we obtain(38)RpN,x=Re12ip·gpei2x,IpN,x=Im12ip·gpei2x.The seriesAp(N,x), Bp(N,x), and Cp(N,x) can be represented, in the terms of Sp(N) and gp(x), in the following way:(39)ApN,x=12·SpN-12·Re12ip·gpei2x,(40)BpN,x=12·SpN+12·Re12ip·gpei2x,(41)CpN,x=12·Im12ip·gpei2x,Also, these series can be represented in the terms ofSp(N), Rp(N,x), and Ip(N,x) in the following simple notation:(42)ApN,x=12·SpN-12·RpN,x,(43)BpN,x=12·SpN+12·RpN,x,(44)CpN,x=12·IpN,x.In the case thatp takes values in the set {0,1,…,4}, the explicit forms of functions gp(z) obtained by p consecutive derivations according to formula (35) are presented in Table 2.Table 2
The functionsgp(z).
g 0 ( z ) z z - 1 ⋅ z N - 1   g 1 ( z ) 2 i z ( z - 1 ) 2 ⋅ N z N + 1 - ( N + 1 ) z N + 1   g 2 ( z ) ( 2 i ) 2 z ( z - 1 ) 3 ⋅ N 2 z N + 2 - ( 2 N 2 + 2 N - 1 ) z N + 1 + ( N 2 + 2 N + 1 ) z N - z - 1   g 3 ( z ) ( 2 i ) 3 z ( z - 1 ) 4 ⋅ N 3 z N + 3 - ( 3 N 3 + 3 N 2 - 3 N + 1 ) z N + 2 + ( 3 N 3 + 6 N 2 - 4 ) z N + 1 - ( N + 1 ) 3 z N + z 2 + 4 z + 1   g 4 ( z ) ( 2 i ) 4 z ( z - 1 ) 5 ⋅ N 4 z N + 4 - ( 4 N 4 + 4 N 3 - 6 N 2 + 4 N - 1 ) z N + 3 + ( 6 N 4 + 12 N 3 - 6 N 2 - 12 N + 11 ) z N + 2 - ( 4 N 4 + 12 N 3 + 6 N 2 - 12 N - 11 ) z N + 1 + ( N + 1 ) 4 z N - z 3 - 11 z 2 - 11 z - 1In the main section of this article, the synchronous model is observed, i.e., the numberN is chosen with respect to the relation Nx0=π. For p=0,1,…,4 and z0=ei2x0, this leads to the values Ep(N,x0) given in Table 3.Table 3
The values of functionsEp(N,x0).
E 0 ( N , x 0 ) 0  E 1 ( N , x 0 ) q 11 z 0 1 / 2 2 i ⋅ I m  z 0 1 / 2 q 11 = N  E 2 ( N , x 0 ) q 21 z 0 3 / 2 + q 22 z 0 1 / 2 + q 23 2 i ⋅ I m  z 0 1 / 2 3 q 21 = N 2 q 22 = - 2 N 2 - 2 N q 33 = N 2 + 2 N  E 3 ( N , x 0 ) q 31 z 0 1 / 2 + q 32 z 0 + q 33 + q 34 z 0 - 1 2 i ⋅ I m  z 0 1 / 2 4 q 31 = N 3 q 32 = - 3 N 3 - 3 N 2 + 3 N q 33 = 3 N 3 + 6 N 2 q 34 = - N 3 - 3 N 2 - 3 N  E 4 ( N , x 0 ) q 41 z 0 5 / 2 + q 42 z 0 3 / 2 + q 43 z 0 1 / 2 + q 44 z 0 - 1 / 2 + q 45 z 0 - 3 / 2 2 i ⋅ I m  z 0 1 / 2 5 q 41 = N 4 z N + 4 q 42 = - 4 N 4 - 4 N 3 + 6 N 2 - 4 N q 43 = 6 N 4 + 12 N 3 - 6 N 2 - 12 N q 44 = - 4 N 4 - 12 N 3 - 6 N 2 + 12 N q 45 = N 4 + 4 N 3 + 6 N 2 + 4 NIf we represent eachz0s in Ep(N,z0) in its expanded form as a polynomial in sin⁡x0 and cos⁡x0 of degree p (the degree of each term in this polynomial in two variables is the sum of the exponents in each term), then Rp(N,x0) can be transformed into the form given in Table 4.Table 4
The values of functionsRp(N,x0).
R 0 ( N , x 0 ) 0   R 1 ( N , x 0 ) N 2   R 2 ( N , x 0 ) 1 2 ⋅ N 2 + N + N ⋅ c t g 2 x 0   R 3 ( N , x 0 ) 1 4 ⋅ 2 N 3 + 3 N 2 + 3 N 2 ⋅ c t g 2 x 0   R 4 ( N , x 0 ) 1 4 ⋅ 2 N 4 + 4 N 3 - 2 N - 6 N ⋅ c t g 4 x 0 + ( 4 N 3 - 8 N ) ⋅ c t g 2 x 0In a similar way, by representing eachz0s in Ep(N,z0) in its expanded form as a polynomial in sin⁡x0 and cos⁡x0 of degree p, Ip(N,x0) can be transformed into the form given in Table 5.Table 5
The values of functionsIp(N,x0).
I 0 ( N , x 0 ) 0   I 1 ( N , x 0 ) - N 2 ⋅ c t g  x 0   I 2 ( N , x 0 ) - N 2 2 ⋅ c t g  x 0   I 3 ( N , x 0 ) 1 4 ⋅ 3 N ⋅ c t g 3 x 0 - ( 2 N 3 - 3 N ) ⋅ c t g  x 0   I 4 ( N , x 0 ) 1 2 ⋅ 3 N 2 ⋅ c t g 3  x 0 - ( N 4 - 3 N 2 ) ⋅ c t g  x 0Now, by formula (42), we obtain values Ap(N,x0), p=0,1,…,4, given in Table 6.Table 6
The values ofAp(N,x0).
A 0 ( N , x 0 ) N 2   A 1 ( N , x 0 ) N 2 4   A 2 ( N , x 0 ) 1 6 ( N 3 - N ) - 1 4 N ⋅ c t g 2 x 0   A 3 ( N , x 0 ) 1 8 ( N 4 - 2 N 2 ) - 3 8 N 2 ⋅ c t g 2 x 0   A 4 ( N , x 0 ) 1 30 ( 3 N 5 - 10 N 3 + 7 N ) + 3 4 N ⋅ c t g 4 x 0 - 1 2 ( N 3 - 2 N ) ⋅ c t g 2 x 0Using formula (43) we obtain values Bp(N,x0), for p=0,1,…,4, given in Table 7.Table 7
The valuesBp(N,x0).
B 0 ( N , x 0 ) N 2   B 1 ( N , x 0 ) 1 4 ( N 2 + 2 N )   B 2 ( N , x 0 ) 1 6 ( N 3 + 3 N 2 + 2 N ) + 1 4 N ⋅ c t g 2 x 0   B 3 ( N , x 0 ) 1 8 ( N 4 + 4 N 3 + 4 N 2 ) + 3 8 N 2 ⋅ c t g 2 x 0   B 4 ( N , x 0 ) 1 30 ( 3 N 5 + 15 N 4 + 20 N 2 - 8 N ) - 3 4 N ⋅ c t g 4 x 0 + 1 2 ( N 3 - 2 N ) ⋅ c t g 2 x 0Finally, by (44) we obtain values Cp(N,x0), for p=0,1,…,4, given in Table 8.Table 8
The valuesCp(N,x0).
C 0 ( N , x 0 ) 0   C 1 ( N , x 0 ) - N 4 ⋅ c t g  x 0   C 2 ( N , x 0 ) - N 2 4 ⋅ c t g  x 0   C 3 ( N , x 0 ) 1 8 ⋅ 3 N ⋅ c t g 3 x 0 - ( 2 N 3 - 3 N ) ⋅ c t g  x 0   C 4 ( N , x 0 ) 1 4 ⋅ 3 N 2  ⋅ c t g 3 x 0 - ( N 4 - 3 N 2 ) ⋅ c t g  x 0Now, by (27) we obtain the elements of Gram's matrix HT·H. Namely, for all i,j=1,2,3 and p=i+j-2, the elements aij are given by(45)a11=N2,(46)a12=a21=TN24,(47)a13=a22=a31=T2·16N3-N-14N·ctg2mωT,(48)a23=a32=T3·18N4-2N2-38N2·ctg2mωT,(49)a33=T4·1303N5-10N3+7N+34N·ctg4mωT-12N3-2N·ctg2mωT.Further, for alli,j=1,2,3 and p=i+j-2, the elements bij are given by(50)b11=N2,(51)b12=b21=T4·N2+2N,(52)b13=b22=b31=T2·16N3+3N2+2N+14N·ctg2mωT,(53)b23=b32=T3·18N4+4N3+4N2+38N2·ctg2mωT,(54)b33=T4·1303N5+15N4+20N2-8N-34N·ctg4mωT+12N3-2N·ctg2mωT.Finally, for alli,j=1,2,3 and p=i+j-2, the elements cij are given by(55)c11=0,(56)c12=c21=-TN4·ctgmωT,(57)c13=c22=c31=-T2N24·ctgmωT,(58)c23=c32=T38·3N·ctg3mωT-2N3-3N·ctgmωT,(59)c33=T44·3N2·ctg3mωT-N4-3N2·ctgmωT.The presented representation of Gram’s matrix shows that the computation of pseudoinverse(HT·H)-1·HT 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:(60)xt=atcos⁡2πf1t+ϕt,where a(t)=a0+a1sin⁡2πfat, ϕ(t)=ϕ0+ϕ1sin⁡2πfϕt, a0=ϕ0=1, a1=ϕ1=0.1, fa=fϕ=5 and N=24.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. Figures1(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.Figure 1
(a) Amplitude estimation, (b) phase estimation.
(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 as(61)TVE=Xr-XeXr,where Xr and Xe 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.Figure 2
Total Vector Error (TVE).


### 4.1. Amplitude Oscillation Case

The ability of the proposed approach to track an arbitrary oscillating amplitude is given in the following example where(62)xt=at∑h=0∞chcos⁡2πhf1tis the input signal to the estimator, with f1=50 Hz and nonzero Fourier coefficients (c1=1, c3=0.4, and c5=0.2), and where a(t) represents an amplitude oscillation, given by the following second-order polynomial(63)at=at2+bt+cwith coefficients a=-4, b=4, and c=0. The corresponding digital signal is sampled at 24 samples/cycle (f=1200Hz). 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.Figure 3
Measured and estimated waveforms in the case of a time-varying signal (62): (a) tracking the first phasors, (b) tracking the 3rd phasors, (c) tracking the 5th phasors, and (d) measured and estimated waveforms of main input signal.
(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)(64)xt=Xm1+kaut-tscos⁡ωt+kφut-tswhere u(t) is a unit step function, ka is the amplitude step size, kφ is the phase step size, and ts 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 Xm=2,ω=2π·50Hz, ka=+0.1p.u., kϕ=π/18, 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.Figure 4
Estimation of (a) amplitude, (b) phase, and (c) TVE in percentage.
(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 Figure4, 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 z 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 Table9 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.Table 9
Performance indices for phasor estimation under step changes in amplitude and phase using the proposed estimation algorithm.
Response time Delay time Overshoot Amplitude step 19.7ms 8.86ms 0.008 [p.u.] Phase step 9.67ms 0.02 [rad]


### 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]).(65)xn=cos⁡2πf1nΔt,0<n<10N1,xn=cos⁡2πf1+5nΔt,10N1≤n.The frequency was calculated as a ratio of roots of characteristic equation based on extracted unknown coefficientsAm0, Am1, Am2, Bm0, Bm1, Bm2 from (12) and 2π [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.Figure 5
Frequency estimation during frequency step test.


### 4.4. Harmonic Infiltration Test

Input signal defined with (66) was used to validate the proposed methods in harmonic condition:(66)xt=atcos⁡ω1t+ϕt+0.05cos⁡5ω1t+0.03cos⁡7ω1t,at=a0+a1sin⁡2πfat,ϕt=ϕ0+ϕ1sin⁡2πfϕt,a0=1,ϕ0=0.5,a1=0.1,ϕ1=0.05,fa=fϕ=5,N=24.Figure6 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 N=24 in this test, it is possible to estimate the dynamic phasor of the first N-1=23 harmonics. To increase the range of harmonics, the sampling number per cycle should be increased.Figure 6
Amplitude estimation in harmonic conditions.


### 4.5. Frequency Response

In order to evaluate the frequency characteristic of the proposed estimator, its transfer function is defined in the form(67)Hz=p^zXz,where X(z) and p^z are z 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].Figure 7
Frequency response: (a) magnitude response, (b) phase response.
(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:(68)at=1+0.5∑i=13et/τicos⁡2πf1t,where the time constants (τi) were generated by a uniform random process in the interval of [20,40] cycles. In a similar way, the three frequencies were randomly generated in the intervals of [1,3], [3,5], and [5,7]Hz. The error is calculated by(69)rms-error=∑an-a^n,where a(n) and a^(n) are real and estimated amplitudes, respectively. Figure 8 shows the histograms of the errors attained by proposed method.Figure 8
RMS error of amplitude estimation.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 by(70)tn=Sn-S^nTM=∑n=r-N1n=rtnwhere Sn is the real (measured) value of sample and S^n is the recomputed sample obtained by the estimated dynamic phasor. Estimation error tn, is the difference between these two quantities.TMis calculated by estimation error in every sample and can be used as a quality measure of phasor estimation. Figure 9 shows theTMas 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.Figure 9
Transient monitor versus noise variance.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 wasN=24.


### 4.8. Computational Complexity and Simulation Time

The inputs of the algorithm are matrixH of the type N×2(K+1) and matrix xmfil of the type N×1 given by (16) and (13), respectively. The output of the algorithm is column matrix y~m of the type 2(K+1)×1 obtained by (19). The proposed algorithm can be decomposed in the following steps:(A1)
Calculation of the elements of Gram’s matrixHT·H by (45)-(59).(A2)
Inversion(HT·H)-1 of Gram’s matrix HT·H.(A3)
Multiplication of(HT·H)-1 with HT.(A4)
Multiplication of(HT·H)-1·HT with xmfil.(A1) If the formulas (45)-(59) are used for calculation of the elements of Gram's matrix HT·H, then the task is reduced to 4(K+1)2 consecutive calculations of specific values based on polynomials. If c0 denotes the maximum of all computational times needed for executing (45)-(59), then the computational time for HT·H does not exceed 4(K+1)2c0.(A2) The solution to (19) includes the inversion of the Gram’s matrix. To find the inverse of 2(K+1)×2(K+1)-matrix HT·H one can use Lower-Upper (LU) decomposition method which takes 4/3·(2(K+1))3+3/2·(2(K+1))2-5/6·2(K+1) FLOPs.(A3) In this step, the matrix(HT·H)-1 of type 2(K+1)×2(K+1) is to be multiplied by the matrix HT of the type 2(K+1)×N which takes the 4(K+1)2N FLOPs.(A4) In this step, the matrix(HT·H)-1·HT of the type 2(K+1)×N is to be multiplied by the matrix xmfil of the type N×1 which takes the 2(K+1)N FLOPs.Based on the all above-mentioned procedures-steps, the proposed algorithm does not exceed(71)4K+12c0+43·2K+13+32·2K+12-56·2K+1+4K+12N+2K+1N=4K2+8K+4c0+323K3+38K2+1273K+15+4K2+10K+6Noperations-FLOPs, which practically defined its computational costs. If the step (A1) is replaced with(A’1) Multiplication ofHT with H,then the computational time of this step is increased with respect to the time needed for performing (A1). SinceH is N×2(K+1)-matrix, it follows that HT·H is 2(K+1)×2(K+1)-matrix and multiplication of matrices gives an algorithm that takes time 4(K+1)2N to multiply the matrix HT of type 2(K+1)×N with the matrix H of type N×2K. Clearly, the time needed for such an algorithm depends on the number of samples N. On the other hand, the time c0 in (A1) is independent of N, 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 c0 is the maximum of all computational times needed for executing (45)-(59), clearly, it is the time needed for realization of (54). If ctgx is regarded as ctgx≈1/x-x/3, then its computational time takes 3 FLOPs for basic arithmetic operations. Consequently, the computational time for (54) can be approximated with finite number c 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 Figure10. 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].Figure 10
Computation time versus sampling frequency.


## 4.1. Amplitude Oscillation Case

The ability of the proposed approach to track an arbitrary oscillating amplitude is given in the following example where(62)xt=at∑h=0∞chcos⁡2πhf1tis the input signal to the estimator, with f1=50 Hz and nonzero Fourier coefficients (c1=1, c3=0.4, and c5=0.2), and where a(t) represents an amplitude oscillation, given by the following second-order polynomial(63)at=at2+bt+cwith coefficients a=-4, b=4, and c=0. The corresponding digital signal is sampled at 24 samples/cycle (f=1200Hz). 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.Figure 3
Measured and estimated waveforms in the case of a time-varying signal (62): (a) tracking the first phasors, (b) tracking the 3rd phasors, (c) tracking the 5th phasors, and (d) measured and estimated waveforms of main input signal.
(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)(64)xt=Xm1+kaut-tscos⁡ωt+kφut-tswhere u(t) is a unit step function, ka is the amplitude step size, kφ is the phase step size, and ts 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 Xm=2,ω=2π·50Hz, ka=+0.1p.u., kϕ=π/18, 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.Figure 4
Estimation of (a) amplitude, (b) phase, and (c) TVE in percentage.
(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 Figure4, 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 z 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 Table9 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.Table 9
Performance indices for phasor estimation under step changes in amplitude and phase using the proposed estimation algorithm.
Response time Delay time Overshoot Amplitude step 19.7ms 8.86ms 0.008 [p.u.] Phase step 9.67ms 0.02 [rad]


## 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]).(65)xn=cos⁡2πf1nΔt,0<n<10N1,xn=cos⁡2πf1+5nΔt,10N1≤n.The frequency was calculated as a ratio of roots of characteristic equation based on extracted unknown coefficientsAm0, Am1, Am2, Bm0, Bm1, Bm2 from (12) and 2π [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.Figure 5
Frequency estimation during frequency step test.


## 4.4. Harmonic Infiltration Test

Input signal defined with (66) was used to validate the proposed methods in harmonic condition:(66)xt=atcos⁡ω1t+ϕt+0.05cos⁡5ω1t+0.03cos⁡7ω1t,at=a0+a1sin⁡2πfat,ϕt=ϕ0+ϕ1sin⁡2πfϕt,a0=1,ϕ0=0.5,a1=0.1,ϕ1=0.05,fa=fϕ=5,N=24.Figure6 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 N=24 in this test, it is possible to estimate the dynamic phasor of the first N-1=23 harmonics. To increase the range of harmonics, the sampling number per cycle should be increased.Figure 6
Amplitude estimation in harmonic conditions.


## 4.5. Frequency Response

In order to evaluate the frequency characteristic of the proposed estimator, its transfer function is defined in the form(67)Hz=p^zXz,where X(z) and p^z are z 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].Figure 7
Frequency response: (a) magnitude response, (b) phase response.
(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:(68)at=1+0.5∑i=13et/τicos⁡2πf1t,where the time constants (τi) were generated by a uniform random process in the interval of [20,40] cycles. In a similar way, the three frequencies were randomly generated in the intervals of [1,3], [3,5], and [5,7]Hz. The error is calculated by(69)rms-error=∑an-a^n,where a(n) and a^(n) are real and estimated amplitudes, respectively. Figure 8 shows the histograms of the errors attained by proposed method.Figure 8
RMS error of amplitude estimation.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 by(70)tn=Sn-S^nTM=∑n=r-N1n=rtnwhere Sn is the real (measured) value of sample and S^n is the recomputed sample obtained by the estimated dynamic phasor. Estimation error tn, is the difference between these two quantities.TMis calculated by estimation error in every sample and can be used as a quality measure of phasor estimation. Figure 9 shows theTMas 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.Figure 9
Transient monitor versus noise variance.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 wasN=24.


## 4.8. Computational Complexity and Simulation Time

The inputs of the algorithm are matrixH of the type N×2(K+1) and matrix xmfil of the type N×1 given by (16) and (13), respectively. The output of the algorithm is column matrix y~m of the type 2(K+1)×1 obtained by (19). The proposed algorithm can be decomposed in the following steps:(A1)
Calculation of the elements of Gram’s matrixHT·H by (45)-(59).(A2)
Inversion(HT·H)-1 of Gram’s matrix HT·H.(A3)
Multiplication of(HT·H)-1 with HT.(A4)
Multiplication of(HT·H)-1·HT with xmfil.(A1) If the formulas (45)-(59) are used for calculation of the elements of Gram's matrix HT·H, then the task is reduced to 4(K+1)2 consecutive calculations of specific values based on polynomials. If c0 denotes the maximum of all computational times needed for executing (45)-(59), then the computational time for HT·H does not exceed 4(K+1)2c0.(A2) The solution to (19) includes the inversion of the Gram’s matrix. To find the inverse of 2(K+1)×2(K+1)-matrix HT·H one can use Lower-Upper (LU) decomposition method which takes 4/3·(2(K+1))3+3/2·(2(K+1))2-5/6·2(K+1) FLOPs.(A3) In this step, the matrix(HT·H)-1 of type 2(K+1)×2(K+1) is to be multiplied by the matrix HT of the type 2(K+1)×N which takes the 4(K+1)2N FLOPs.(A4) In this step, the matrix(HT·H)-1·HT of the type 2(K+1)×N is to be multiplied by the matrix xmfil of the type N×1 which takes the 2(K+1)N FLOPs.Based on the all above-mentioned procedures-steps, the proposed algorithm does not exceed(71)4K+12c0+43·2K+13+32·2K+12-56·2K+1+4K+12N+2K+1N=4K2+8K+4c0+323K3+38K2+1273K+15+4K2+10K+6Noperations-FLOPs, which practically defined its computational costs. If the step (A1) is replaced with(A’1) Multiplication ofHT with H,then the computational time of this step is increased with respect to the time needed for performing (A1). SinceH is N×2(K+1)-matrix, it follows that HT·H is 2(K+1)×2(K+1)-matrix and multiplication of matrices gives an algorithm that takes time 4(K+1)2N to multiply the matrix HT of type 2(K+1)×N with the matrix H of type N×2K. Clearly, the time needed for such an algorithm depends on the number of samples N. On the other hand, the time c0 in (A1) is independent of N, 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 c0 is the maximum of all computational times needed for executing (45)-(59), clearly, it is the time needed for realization of (54). If ctgx is regarded as ctgx≈1/x-x/3, then its computational time takes 3 FLOPs for basic arithmetic operations. Consequently, the computational time for (54) can be approximated with finite number c 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 Figure10. 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].Figure 10
Computation time versus sampling frequency.


## 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.



---

*Source: 7613814-2018-11-26.xml*