A Self-Adaptive Contractive Algorithm for Enhanced Dynamic Phasor Estimation

TL;DR

This paper introduces a self-adaptive contractive algorithm that improves dynamic phasor estimation in power systems by transforming signals, estimating parameters with filtering, and then recovering original signals with enhanced accuracy and harmonic rejection.

Contribution

The paper presents a novel self-adaptive contractive algorithm with convergence guarantees and superior harmonic rejection for dynamic phasor estimation in power systems.

Findings

Validated with IEEE standard tests under realistic noise conditions

Demonstrates convergence and robustness in diverse operating scenarios

Outperforms existing methods in harmonic rejection and accuracy

Abstract

In this paper, a self-adaptive contractive (SAC) algorithm is proposed for enhanced dynamic phasor estimation in the diverse operating conditions of modern power systems. At a high-level, the method is composed of three stages: parameter shifting, filtering and parameter unshifting. The goal of the first stage is to transform the input signal phasor so that it is approximately mapped to nominal conditions. The second stage provides estimates of the phasor, frequency, rate of change of frequency (ROCOF), damping and rate of change of damping (ROCOD) of the parameter shifted phasor by using a differentiator filter bank (DFB). The final stage recovers the original signal phasor parameters while rejecting misleading estimates. The most important features of the algorithm are that it offers convergence guarantees in a set of desired conditions, and also great harmonic rejection. Numerical examples, including the IEEE C37.118.1 standard tests with realistic noise levels, as well as fault conditions, validate the proposed algorithm.

Figures

Tables

Table 1. Table I: Worst-case error metrics for the different tests of Sections V-A , V-B and V-C

		AE (pu)	PE (rad)	TVE (%)	FE (Hz)	DE (1/s)	RFE (Hz/s)	RDE (1/s2)
A1	SAC algorithm	$8.9415\cdot {10}^{-5}$	$9.0944\cdot {10}^{-5}$	0.0096	0.0012	0.0078	0.1688	1.0641
	IEEE Std. algorithm	0.2193	$8.2336\cdot {10}^{-5}$	21.9288	0.0034	0.1423	14.7719	598.5252
A2	SAC algorithm	$9.7480\cdot {10}^{-5}$	$1.0019\cdot {10}^{-4}$	0.0102	0.0072	0.0078	0.1744	1.0670
	IEEE Std. algorithm	0.2231	0.0026	22.3070	0.1782	6.5213	554.8150	$2.5687\cdot {10}^4$
A3	SAC algorithm	$9.4300\cdot {10}^{-5}$	$1.0835\cdot {10}^{-4}$	0.0112	0.0084	0.0144	0.1555	1.0552
	IEEE Std. algorithm	0.3031	0.0487	30.3112	7.6028	321.4530	$1.7683\cdot {10}^4$	$1.1325\cdot {10}^6$
B1	SAC algorithm	$6.3672\cdot {10}^{-5}$	$6.5001\cdot {10}^{-5}$	0.0068	0.0010	0.0055	0.1004	0.4069
	IEEE Std. algorithm	0.2187	$3.5572\cdot {10}^{-4}$	21.8652	0.0036	0.1383	9.6899	450.7530
B2	SAC algorithm	$7.5743\cdot {10}^{-5}$	$7.4745\cdot {10}^{-5}$	0.0084	0.0016	0.0090	0.2663	1.4563
	IEEE Std. algorithm	$7.7686\cdot {10}^{-4}$	$7.8009\cdot {10}^{-4}$	0.1144	0.0032	0.0199	6.7630	45.1783
B3	SAC algorithm	$1.8363\cdot {10}^{-4}$	$1.7879\cdot {10}^{-4}$	0.0233	0.0121	0.0806	1.2611	8.6851
	IEEE Std. algorithm	0.0047	0.0047	0.6981	0.0261	0.1603	8.2548	55.2348
C1	SAC algorithm	$6.5510\cdot {10}^{-5}$	$6.4767\cdot {10}^{-5}$	0.0070	$8.8005\cdot {10}^{-4}$	0.0067	0.1372	1.0653
	IEEE Std. algorithm	0.0026	$4.7742\cdot {10}^{-4}$	0.2501	0.0023	0.0152	6.7818	48.3326
C2	SAC algorithm	$7.5912\cdot {10}^{-5}$	$6.6228\cdot {10}^{-5}$	0.0074	$8.5263\cdot {10}^{-4}$	0.0057	0.1290	0.7502
	IEEE Std. algorithm	$7.8080\cdot {10}^{-4}$	0.0010	0.1303	0.0023	0.0136	8.5739	50.0218
C3	SAC algorithm	$7.3535\cdot {10}^{-5}$	$6.6000\cdot {10}^{-5}$	0.0075	$8.3744\cdot {10}^{-4}$	0.0074	0.1367	1.2403
	IEEE Std. algorithm	0.0011	$7.0928\cdot {10}^{-4}$	0.1253	0.0025	0.0161	7.5543	52.7143

Table 2. Table II: Amplitude and phase step performance results

		Overshoot (%)	$t_{\text{res},\text{TVE}}$ (s)	$t_{\text{res},\text{FE}}$ (s)	$t_{\text{res},\text{DE}}$ (s)	$t_{\text{res},\text{RFE}}$ (s)	$t_{\text{res},\text{RDE}}$ (s)
Amplitude	SAC algorithm	1.5877	0.0141	0	0.0339	0	0.0385
	IEEE Std. algorithm	0.0282	0.0271	0	0.0349	$\infty$	$\infty$
Phase	SAC algorithm	0.7111	0.0161	0.0339	0.0339	0.0401	0.0380
	IEEE Std. algorithm	0.0191	0.0307	0.0349	0.0359	$\infty$	$\infty$

Equations

$$X(t_r+\tau )=X(t_r)e^{b(\tau )}{e}^{j\varphi (\tau )},\tau \in \mathcal{T}\triangleq [-\delta ,\delta ],$$

$$b_{\text{pr}}(\tau )$$

$$Z(\tau )=X(t_r+\tau )e^{-b_{\text{pr}}(\tau )}{e}^{-j{\varphi }_{\text{pr}}(\tau )}=a_Z(\tau )e^{j{\varphi }_Z(\tau )},$$

$${\varphi }_Z(0)$$

$${\omega }_Z(0)$$

$${\alpha }_Z(0)$$

$$a_Z(0)$$

$${\sigma }_Z(0)$$

$${\gamma }_Z(0)$$

$${\varphi }_{\text{pr},k}[m]={\omega }_{\text{pr},k}Tm+\frac{1}{2}{\alpha }_{\text{pr},k}{T}^2{m}^2,m\in \mathcal{W},$$

$$b_{\text{pr},k}[m]=-{\sigma }_{\text{pr},k}Tm-\frac{1}{2}{\gamma }_{\text{pr},k}{T}^2{m}^2,m\in \mathcal{W},$$

$${\mathbf{Z}}_k[n_r]={\mathbf{D}}_{\text{pr},k}^{-1}\mathbf{Y}[n_r].$$

$${\hat{Z}}_k[n_r]$$

$${\hat{\omega }}_{Z_k}[n_r]$$

$${\hat{\sigma }}_{Z_k}[n_r]$$

$$\underset{{\mathbf{a}}_i}{\min }$$

$${\mathbf{a}}_i=\widetilde{{\mathbf{a}}_i}+{\mathbf{P}}_i^{-1}{\mathbf{G}}_i^T({\mathbf{G}}_i{\mathbf{P}}_i^{-1}{\mathbf{G}}_i^T{)}^{-1}({\mathbf{c}}_i-{\mathbf{G}}_i\widetilde{{\mathbf{a}}_i}).$$

$$f_i({\tilde{\mathbf{a}}}_i)=\frac{1}{2\pi }\Vert {\tilde{A}}_i{\Vert }^2={\tilde{\mathbf{a}}}_i^T{\mathbf{P}}_i{\tilde{\mathbf{a}}}_i,i=0,1,2,$$

$${\tilde{A}}_i(l{\omega }_{0}T)=0,i=0,1,2,l=1,\dots ,M.$$

$$\left|{\omega }_s-\frac{{\tilde{A}}_1({\omega }_{s}T)}{{\tilde{A}}_0({\omega }_{s}T)}\right|\le L_{\omega }|{\omega }_s|,\forall {\omega }_s\in {\Omega }_{\text{con}}.$$

$$X[n_r+m]$$

$$e_{\omega ,k+1}=e_{\omega ,k}-{\hat{\omega }}_{Z_{k+1}}[n_r]=e_{\omega ,k}-\frac{A_1(e_{\omega ,k}T)}{A_0(e_{\omega ,k}T)}\triangleq f_{\omega }(e_{\omega ,k}).$$

$$|f_{\omega }({\omega }_s)|=\left|{\omega }_s-\frac{A_1({\omega }_{s}T)}{A_0({\omega }_{s}T)}\right|\le L_{\omega }|{\omega }_s|,\forall {\omega }_s\in {\Omega }_{\text{con}}.$$

Full text

I Introduction

The proliferation of nonlinear loads and incorporation of distributed energy resources (DERs) are dramatically changing modern power systems [1]. In particular, their typically large inertia is steadily decreasing [2], which means that they are operating in more varied conditions, making them more prone to disturbances. Under this scenario, monitoring, protection and control of power systems are becoming much more challenging problems. However, the increasing deployment of phasor measurement units (PMUs) is allowing the development of wide-area measurement systems (WAMS), which in turn provides the opportunity to enhance the grid stability assessment, state estimation, fault detection, generator model validation, etc. [3].

Phasor estimation accuracy has an impact on all the aforementioned applications and is, therefore, a problem of fundamental importance in modern power systems. Although the IEEE Std. C37.118.1-2011 [4] on synchrophasors and its recent amendment C37.118.1a-2014 [5] (from now on simply denoted as the IEEE Std.) are firmly adopted for evaluation and certification of PMUs, it has received criticism, some of which are the following. First, given the importance of phase estimation accuracy in several applications, including active power flow computations and grid restoration [6], phase error (PE) may be a more useful metric than total vector error (TVE), a joint measure of the amplitude and phase errors [7]. On the other hand, amplitude error (AE) may be critical for other cases, such as high-impedance fault detection and voltage sag analysis [8]. Therefore, besides the TVE, it is important to characterize AE and PE separately. Another criticism arises from the modification and suspension of some of the RFE limits in the harmonic and interharmonic interference tests, since this could seriously limit the application of PMUs in island detection schemes and generation control for grid stabilization [9]. Similarly, since islanded grid operation may produce temporary frequency offsets even higher than those stipulated in the IEEE Std. [6], there is a current interest in wide frequency range algorithms [10, 11]. In addition, the IEEE Std. does not consider the impact of acquisition noise but this effect may be significant, particularly for frequency and rate of change of frequency (ROCOF) measurements [12]. On the other hand, in distribution systems the angular differences between nodes are quite small, so much more accuracy is typically required [8]. From the previous discussion, one could argue that the complete specification of PMUs performance requirements in different scenarios and applications is still an open problem. In any case, it is clear that next-generation PMUs should aim at providing enhanced performance.

In this context, a considerable amount of research was conducted in synchrophasor estimation algorithms [3], and interest in this area is still growing. These algorithms may be classified broadly in two main categories: non-adaptive (e.g., [13, 14]) and adaptive (e.g., [10, 15]). On the one hand, non-adaptive algorithms have fixed parameters so that they generally provide good performance only in a limited set of operating conditions. On the other hand, adaptive solutions are much more flexible and, in principle, allow for good performance on a larger set of conditions. For instance, the well-known Interpolated DFT (IpDFT) algorithms enhance the DFT performance by introducing an adaptive compensation of the phasor estimate based on an estimate of the frequency deviation from nominal [16]. Similarly, by including the harmonic components in a state-space model, a Kalman filter can provide not only better off-nominal behavior than a single DFT but also improve its off-nominal harmonic rejection [10]. Finally, the so-called polynomial phase-locked loop Taylor-Fourier filters use an adaptive phase polynomial to center the Taylor expansion around the actual operating condition [15], thus improving the approximation with respect to the plain Taylor-Fourier non-adaptive method [17]. However, adaptive systems are generally more computationally demanding. More seriously, convergence issues may appear in some operating conditions if this is neglected at the design stage.

The main contribution of this work is to present the design of a novel adaptive dynamic phasor estimation algorithm. The proposed method offers a notable overall performance, a high harmonic rejection capability, and convergence guarantees in wide operating conditions that are specified at the design stage. It is based on a self-adaptive strategy working around a prototype differentiator filter bank (DFB) which allows to relax the design requirements of the filters and therefore leads to an enhanced performance. Convergence conditions are included in the design of the prototype system with the same convex semi-infinite programming (CSIP) framework as that used in [18]. Although the tool used for the design of the filters is the same as that in [18], it should be emphasized that the algorithms are quite different. In fact, note that the method in [18] is a fixed filtering algorithm (i.e., non-adaptive).

II Prior Model and Parameter Shifting

For convenience of presentation, the signal phasor or synchrophasor (see definitions in [4]) model is first presented in the continuous-time domain. Let $X(t)=a_{X}(t)e^{j\phi_{X}(t)}$, $t\in\mathbb{R}$, be the signal phasor, where $a_{X}(t)>0$ is its (instantaneous) amplitude and $\phi_{X}(t)$ its phase. Then, around a particular fixed time, say the $r$-th reporting time $t_{r}$, we can write

[TABLE]

where $b(\tau)=\ln a(\tau)$ is the log-amplitude modulation function, $a(\tau)=a_{X}(t_{r}+\tau)/a_{X}(t_{r})$ is the amplitude modulation function, $\phi(\tau)=\phi_{X}(t_{r}+\tau)-\phi_{X}(t_{r})$ is the phase modulation function and $\delta>0$ is a parameter that controls the width of the interval over which the signal phasor is to be approximated. In order to introduce the signal phasor prior model, the second-order log-amplitude and phase prior polynomials are defined as follows:

[TABLE]

where $\sigma_{\textrm{pr}},\gamma_{\textrm{pr}},\omega_{\textrm{pr}},\alpha_{\textrm{pr}}$ are, respectively, the prior damping, rate of change of damping (ROCOD), frequency and rate of change of frequency (ROCOF). For instance, during a fault condition, a damped sinusoidal model is appropriate for describing the synchrophasor [19] and damping gives important information about the instantaneous stability of the power system [20]. Moreover, second order information provided by the ROCOD may be valuable in order to predict critical grid operating conditions sooner. Finally, the signal phasor prior model may be written as $X_{\textrm{pr}}(\tau)=e^{b_{\textrm{pr}}(\tau)}e^{j\phi_{\textrm{pr}}(\tau)}$, $\tau\in\mathcal{T}$. Notice that this model is related to different proposals introduced in previous works. For instance, in [15] a Taylor polynomial is used to approximate the dynamic phasor with respect to an adaptive phase prior polynomial, but no adaptive prior model for the amplitude is proposed. On the other hand, Prony’s method in [20] corresponds to an adaptive signal phasor model $X(t_{r}+\tau)=X(t_{r})e^{-\sigma\tau}e^{j\omega\tau}$, $\tau\in\mathcal{T}$, where $\sigma$ and $\omega$ are the damping and frequency parameters, respectively. Therefore, $X_{\textrm{pr}}(\tau)$ gives a second-order generalization of this expression. Finally, in [17], a Taylor expansion is proposed for the dynamic signal phasor. Here, instead, the approximations are made separately for amplitude and phase.

In what follows, the parameter shifting process is defined and analyzed. Let

[TABLE]

be the parameter shifted signal phasor, where its amplitude and phase are, respectively, ${a_{Z}(\tau)=a_{X}(t_{r})e^{b(\tau)-b_{\textrm{pr}}(\tau)}}$ and ${\phi_{Z}(\tau)=\phi_{X}(t_{r})+\phi(\tau)-\phi_{\textrm{pr}}(\tau)}$. Consider now the effect of this transformation at the reporting time $t_{r}$ (i.e., at $\tau=0$). It is clear that for the phase, frequency (i.e., the first time derivative of phase) and ROCOF (i.e., the second time derivative of phase), the following relations hold:

[TABLE]

On the other hand, for amplitude, damping (i.e., minus the first time derivative of log-amplitude) and ROCOD (i.e., minus the second time derivative of log-amplitude) one has, respectively,

[TABLE]

In summary, the amplitude and phase of $Z(\tau)$ at $\tau=0$ are the same as those of $X(t_{r})$, while the frequency, damping, ROCOF, and ROCOD are shifted by the prior parameters. This means that if the prior parameters match the ones of $X(t_{r})$, the frequency, damping, ROCOF and ROCOD of $Z(\tau)$ will be zero at $\tau=0$. In general, if the match is not perfect but the prior parameters are not too different from the actual values, the magnitude of the parameters of $Z(\tau)$ at $\tau=0$ will be reduced with respect to those of $X(t_{r})$. Therefore, the signal $\{Z(\tau)\}_{\tau\in\mathcal{T}}$ is expected to be smoother than the signal ${\{X(t_{r}+\tau)\}_{\tau\in\mathcal{T}}}$ if $\delta$ is sufficiently small.

III Self-Adaptive Contractive Algorithm

The self-adaptive contractive (SAC) algorithm proposed here is based on three stages: parameter shifting, filtering performed by an DFB, and parameter unshifting. In order to perform the system adaptation, the frequency, damping, ROCOF and ROCOD estimates of the signal phasor are fed back into the algorithm. This information is used to perform a parameter shifting in order to attract the input signal phasor to nominal conditions and, if necessary, to perform an adaptation of the filters coefficients in order to match their zeros with the off-nominal harmonic frequencies. After filtering the parameter shifted signal, a parameter unshifting of the estimates is performed to obtain the ones corresponding to the original signal phasor. The algorithm is iterative so it will produce estimates depending on an iteration index $k$, being the reported estimates the ones corresponding to the final iteration. The proposed method is shown schematically in Fig. 1 and summarized as Algorithm 1 at the end of this section for convenience and reference. In the following, a detailed description of each stage of the algorithm is provided.

III-A Input Signal Vector

Consider the discrete-time signal phasor $X[n]$, $n\in\mathbb{Z}$. From (1), within the measurement window at the $r$-th reporting time, $X[n]$ takes the general form $X[n_{r}+m]=X[n_{r}]e^{b[m]}e^{j\phi[m]}$, $m\in\mathcal{W}$, where $n_{r}=t_{r}/T$ is assumed to be an integer, $b[m]$ and $\phi[m]$ are, respectively, the discrete-time log-amplitude and phase modulation functions, and $\mathcal{W}=\{-N/2,\ldots,N/2\}$ is the window index set, being $N+1$ the window length, which is assumed to be odd. The input signal to the SAC algorithm is a complex baseband signal obtained from a real three-phase signal sampled with a sampling period $T$ and by performing an abc-dq transformation. As in [18], the input signal is denoted as $Y[n]=X[n]+I[n]$, where $I[n]$ is a general disturbance signal which includes the harmonics, the negative-sequence component, noise, etc. The samples around each reporting time index $n_{r}$ are collected in the input signal vector ${\boldsymbol{Y}}[n_{r}]=\left[Y[n_{r}-N/2],\ldots,Y[n_{r}],\ldots,Y[n_{r}+N/2]\right]^{T}$.

III-B Parameter Shifting

From (2), the discrete-time phase prior polynomial at iteration $k$ is defined as follows

[TABLE]

where $\omega_{\textrm{pr},k}$ and $\alpha_{\textrm{pr},k}$ are, respectively, the prior frequency and ROCOF parameters at iteration $k$. These parameters are initialized using the previously reported estimates as ${\omega_{\textrm{pr},0}=\widehat{\omega}_{X}[n_{r-1}]+\Delta t\,\widehat{\alpha}_{X}[n_{r-1}]}$ and $\alpha_{\textrm{pr},0}=\widehat{\alpha}_{X}[n_{r-1}]$, where $\Delta t=1/F_{s}$, being $F_{s}$ the reporting rate of the PMU. For the first window of samples, we use the nominal values of the parameters for initialization. Starting from these values, $\omega_{\textrm{pr},k}$ and $\alpha_{\textrm{pr},k}$ are adapted to match the signal phase dynamics at the current reporting time index $n_{r}$. Similarly, the log-amplitude polynomial is defined as follows:

[TABLE]

where $\sigma_{\textrm{pr},k}$ and $\gamma_{\textrm{pr},k}$ are the prior damping and ROCOD parameters at iteration $k$. The initialization of these parameters is completely analogous to that of the parameters $\omega_{\textrm{pr},k}$ and $\alpha_{\textrm{pr},k}$. The adaptation scheme for the parameters $\omega_{\textrm{pr},k}$, $\alpha_{\textrm{pr},k}$, $\sigma_{\textrm{pr},k}$ and $\gamma_{\textrm{pr},k}$ is presented in Section III-E.

Now, the phasor prior matrix $\boldsymbol{D}_{\textrm{pr},k}\in\mathbb{C}^{(N+1)\times(N+1)}$ is defined as a diagonal matrix with elements $(\boldsymbol{D}_{\textrm{pr},k})_{m,m}=e^{b_{\textrm{pr},k}[m]}\;e^{j\phi_{\textrm{pr},k}[m]}$, where $m\in\mathcal{W}$. The parameter shifted input signal vector is then obtained as follows:

[TABLE]

Notice that (12) is the discrete-time equivalent of (3).

III-C Filtering

The estimates (respectively, the parameter shifted signal phasor estimate and its derivatives, frequency, damping, ROCOF and ROCOD) $\widehat{Z}_{k}[n_{r}]$, $\widehat{Z}^{\prime}_{k}[n_{r}]$, $\widehat{Z}^{\prime\prime}_{k}[n_{r}]$, $\widehat{\omega}_{{Z}_{k}}[n_{r}]$, $\widehat{\sigma}_{{Z}_{k}}[n_{r}]$, $\widehat{\alpha}_{{Z}_{k}}[n_{r}]$ and $\widehat{\gamma}_{Z_{k}}[n_{r}]$ are obtained by filtering ${\boldsymbol{Z}}_{k}[n_{r}]$ with the current DFB, which is represented by the linear-phase coefficient vectors $\mathbf{a}_{0}\in\mathbb{R}^{N/2+1}$, $\mathbf{a}_{1}\in\mathbb{R}^{N/2}$ and $\mathbf{a}_{2}\in\mathbb{R}^{N/2+1}$ [18]. The expressions are shown in (III-C) at the bottom of the next page (the rationale behind these computations can be found in [17]). For simplicity of notation, three constant matrices are defined, $\boldsymbol{Q}_{0}\in\mathbb{R}^{(N/2+1)\times(N+1)}$, $\boldsymbol{Q}_{1}\in\mathbb{R}^{N/2\times(N+1)}$ and $\boldsymbol{Q}_{2}\in\mathbb{R}^{(N/2+1)\times(N+1)}$, to express the impulse response of the linear filters in terms of $\mathbf{a}_{0}$, $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$, respectively. The expressions of these matrices are omitted for space reasons.

III-D Parameter Unshifting

The parameter unshifting stage is a transformation required to counteract the effects of the parameter shifting stage to obtain the desired parameter estimates. Note that, as in Section II, the parameter unshifted estimates are denoted analogously to the parameter shifted ones by replacing $Z$ with $X$. From (4) and (7), the transformation required for the phasor is trivial: $\widehat{X}_{k}[n_{r}]=\widehat{Z}_{k}[n_{r}]$. For the other parameters, the strategy is more complex for the sake of enhanced performance and robustness. Note that in an actual power system the range of values of the parameters is limited. However, some transient scenarios may produce very large misleading parameter estimates within a short time window. This observation leads to the following strategy to address the issue, which will be exemplified with the frequency estimate, but it applies in an exactly analogous manner to the other parameters. If $|\widehat{\omega}_{Z_{k}}[n_{r}]+\omega_{\textrm{pr},k}|\leq\Delta\omega_{\max}$, where $\Delta\omega_{\max}$ is an upper frequency threshold specified by the algorithm designer, a normal regime operation is considered. In such case, using the relation (5), the transformation is simply $\widehat{\omega}_{X_{k}}[n_{r}]=\widehat{\omega}_{Z_{k}}[n_{r}]+\omega_{\textrm{pr},k}$. Otherwise, the estimate $\widehat{\omega}_{Z_{k}}[n_{r}]$ is disregarded and the transformation reduces to using the prior value: $\widehat{\omega}_{X_{k}}[n_{r}]=\omega_{\textrm{pr},k}$. Note that to determine the parameter unshifting process completely we need the additional upper thresholds of damping, ROCOF and ROCOD (respectively, $\Delta\sigma_{\max},\Delta\alpha_{\max},\Delta\gamma_{\max}$).

III-E Prior Parameters Adaptation

Note that at a particular reporting time $n_{r}$, there may be some parameters that require adaptation and others that do not. For instance, during the onset of a frequency ramp, only the ROCOF prior should be modified. In addition, a large and fast frequency change may indicate a phase step, in which case the algorithm should not modify wildly its prior parameter and filter coefficients but rather perform as prescribed by the prototype filters in that transient condition. Concretely, the algorithm should avoid very large changes in the prior parameters for the sake of reliability as well as very small changes for the sake of computational efficiency. With this purpose in mind, we introduce the algorithm lower thresholds: $\Delta\omega_{\min}$, $\Delta\alpha_{\min}$, $\Delta\sigma_{\min}$, and $\Delta\gamma_{\min}$. Now, the parameter adaptation set (i.e., the prior parameters for which the adaptation will be made at iteration $k$) $P_{k}\subseteq\{\omega_{\textrm{pr},k},\alpha_{\textrm{pr},k},\sigma_{\textrm{pr},k},\gamma_{\textrm{pr},k}\}$ is defined according to whether the quantities $|\widehat{\omega}_{Z_{k}}[n_{r}]|,|\widehat{\alpha}_{Z_{k}}[n_{r}]|,|\widehat{\sigma}_{Z_{k}}[n_{r}]|$, and $|\widehat{\gamma}_{Z_{k}}[n_{r}]|$ belong to its corresponding ranges defined by the lower and upper algorithm thresholds. For example, $\omega_{\textrm{pr},k}\in P_{k}$ if and only if ${\Delta\omega_{\min}\leq|\widehat{\omega}_{Z_{k}}[n_{r}]|\leq\Delta\omega_{\max}}$. These criteria are illustrated in Fig. 2 for the frequency. For each parameter in $P_{k}$, its value is updated according to the following simple rules: $\omega_{\textrm{pr},k+1}=\widehat{\omega}_{X_{k}}[n_{r}]$, $\alpha_{\textrm{pr},k+1}=\widehat{\alpha}_{X_{k}}[n_{r}]$, $\sigma_{\textrm{pr},k+1}=\widehat{\sigma}_{X_{k}}[n_{r}]$, and $\gamma_{\textrm{pr},k+1}=\widehat{\gamma}_{X_{k}}[n_{r}]$. The algorithm will then iterate until either $P_{k}=\emptyset$ or $k$ becomes greater than some predefined limit $k_{\max}$.

III-F Filters Coefficients Adaptation

In order to guarantee that the DFB behaves properly with large off-nominal and high total harmonic distortion levels, it is required to adapt the filters to have zeros as close as possible to the off-nominal harmonic frequencies. In addition, it is desirable that the filters have ideal gains at the fundamental frequency. It is possible to express these conditions as linear constraints in the filter coefficients $\mathbf{a}_{i}$ and therefore they can be written in matrix form as $\boldsymbol{G}_{i}\mathbf{a}_{i}=\mathbf{c}_{i}$. The matrices $\boldsymbol{G}_{i}$ depend on $\omega_{h,l}=l(\omega_{0}+\widehat{\omega}_{X_{k-1}}[n_{r}])$, $l=1,\ldots,M$, the off-nominal harmonic frequency estimates available at iteration $k$ prior to filtering, where $\omega_{0}$ is the angular nominal power system frequency, and $M$ is the number of harmonics considered. The expressions of $\boldsymbol{G}_{i}$ are omitted for space reasons, while $\mathbf{c}_{0}=[1,0,\ldots,0]^{T}\in\mathbb{R}^{N/2+1},\mathbf{c}_{1}=\mathbf{0}\in\mathbb{R}^{N/2},\mathbf{c}_{2}=\mathbf{0}\in\mathbb{R}^{N/2+1}$.

On the other hand, the goal of the filters adaptation is to modify the filters coefficients as little as possible in order to preserve the prototype filter features, but also enforcing the aforementioned constraints. The $L^{2}[-\pi,\pi]$ squared distance between the adapted filters and the prototype filters are, by Parseval’s theorem, given by ${\|\boldsymbol{P}_{i}^{1/2}(\mathbf{a}_{i}-\widetilde{\mathbf{a}}_{i})\|^{2}}$, where $\boldsymbol{P}_{0}=\textrm{diag}(1,\frac{1}{2},\ldots,\frac{1}{2})=\boldsymbol{P}_{2}$, and $\boldsymbol{P}_{1}=\textrm{diag}(\frac{1}{2},\ldots,\frac{1}{2})$, are diagonal matrices [21]. Therefore, the adapted filters coefficients are obtained by solving the following linearly constrained least-squares optimization problem:

[TABLE]

The solution to this problem can be readily found by using the method of Lagrange multipliers, yielding

[TABLE]

Note that the adaptation of the filters coefficients will occur if and only if $\omega_{\textrm{pr},k}\in P_{k}$.

III-G Computational Cost

The online computational cost of each stage of the SAC algorithm for one iteration is as follows. The computation of the prior polynomials, the parameter shifting operation, and the filtering stage have a cost of $O(N)$ operations. The parameter unshifting and prior parameters adaptation take only $O(1)$ operations. Finally, the filters coefficients adaptation has a cost of $O(M^{2}N)$ operations, where it is assumed that $N\gg M$. Therefore, the adaptation of the filters coefficients is clearly the dominant cost component of the algorithm. Thus, a proper selection of the frequency thresholds $\Delta\omega_{\min}$ and $\Delta\omega_{\max}$ should be made to minimize it as much as possible for a given desired performance in a practical setting. Note also that the overall cost of the algorithm is similar or lower than that of comparable self-adaptive solutions [22, 15].

IV Prototype Differentiator Filter Bank Design

The DFB consists of a parallel architecture of three linear-phase filters to estimate the signal phasor and its first two derivatives [18]. The order of the filters are all assumed to be equal to $N$, and the expressions for their amplitude responses are given by $\widetilde{A}_{i}(\omega T)=\mathbf{g}_{i}^{T}(\omega T)\;\widetilde{\mathbf{a}}_{i},\ i=0,1,2,$ where $\mathbf{g}_{0}(\omega T)=[1,\cos(\omega T),\ldots,\cos(\omega TN/2)]^{T}$, $\mathbf{g}_{1}(\omega T)=[\sin(\omega T),\ldots,\sin(\omega TN/2)]^{T}$ and $\mathbf{g}_{2}(\omega T)=\mathbf{g}_{0}(\omega T)$. The analysis and design of an DFB structure using the CSIP framework is thoroughly presented in [18], so the reader is referred to that paper for the details of this approach. Here, we only discuss the differences that arise in the design of the DFB when it is to be used in the SAC algorithm. As will be shown, the self-adaptation scheme relaxes most of the performance constraints of the problem. Moreover, in many cases, the constraints imply that, under such cases, the different errors, which are the total vector error (TVE), frequency error (FE), rate of change of frequency error (RFE), damping error (DE) and rate of change of damping error (RDE), converge to zero (see Appendix), which is a much stronger result than a performance bound on them as in the formulation of [18].

IV-A Cost Functions

The cost functions are chosen in order to obtain minimum norm filters, yielding

[TABLE]

It can be proved that this choice minimizes the error contribution of a small model mismatch between the prior phasor and the signal phasor and, in fact, it is interesting to note that the DFT and Taylor-Fourier filters can also be interpreted as minimum norm filters. However, note that the design constraints and the structure of these algorithms are very different.

IV-B Interference Rejection Constraints

As discussed in Section III-F, in off-nominal frequency conditions, the filters adaptation strategy forces zeros at the corresponding frequencies. Therefore, only linear constraints on $\widetilde{\mathbf{a}}_{i}$ to ensure zeros at nominal harmonics are required:

[TABLE]

This relaxation of the original harmonic rejection constraints is of critical importance. In fact, if one includes the joint effect of off-nominal frequency conditions and multiple high-level harmonics as constraints in the original CSIP filter design problem, the resulting optimization problems may result unfeasible [18].

IV-C Contractive Constraints

In steady-state off-nominal frequency conditions, the signal phasor around a fixed reporting time $n_{r}$ can be modeled as $X[n_{r}+m]=X[n_{r}]e^{j\omega_{s}Tm}$, where $\omega_{s}\in\Omega_{\textrm{con}}$. A sufficient condition for the convergence of $\widehat{\omega}_{X_{k}}[n_{r}]$ to $\omega_{s}$ as $k\to\infty$ is given by the following convex semi-infinite (CSI) contractive constraint on $\widetilde{\mathbf{a}}_{1}$ (see Appendix for details)

[TABLE]

The set $\Omega_{\textrm{con}}$ is called the contractive frequency set and is defined as $\Omega_{\textrm{con}}=[-\omega_{\textrm{con}}\,T,\omega_{\textrm{con}}\,T]$, where $\omega_{\textrm{con}}$ is the frequency contractive range (in rad/s). Also, $L_{\omega}\in(0,1)$ is the frequency contraction parameter, which controls the rate of convergence of the frequency estimate. This constraint, imposed on the design of the filter $\widetilde{A}_{1}$ given $\widetilde{A}_{0}$, is convex. It is interesting to compare this constraint with the one obtained in [18]. Note that the bound $L_{\omega}|\omega_{s}|$ is greater than $2\pi\textrm{FE}_{\textrm{STA}}$ for all frequencies $f_{s}$ such that $|f_{s}|>\textrm{FE}_{\textrm{STA}}/L_{\omega}$, which in general is in the order of a few mHz. This shows that a relaxation of the original constraint is achieved for significant off-nominal frequencies.

Similar contractive constraints are posed for damping, ROCOF, and ROCOD, by considering, respectively, the following signal phasor models:

[TABLE]

Note that, once again, the constraints are completely specified by their contractive ranges ($\sigma_{\textrm{con}},\alpha_{\textrm{con}},\gamma_{\textrm{con}}$) and contraction parameters ($L_{\sigma},L_{\alpha},L_{\gamma}$). The constraint expressions are omitted for the sake of brevity.

IV-D Other Constraints

A contractive constraint in frequency, damping, ROCOF and ROCOD is sufficient to guarantee good performance in amplitude and phase modulation conditions, so no additional constraints are introduced for this scenario. However, it should be noted that in these cases there is a small model mismatch between the prior phasor and the true signal phasor, even with perfect prior parameters, but its order is $O(T^{3})$. Moreover, as discussed in Section IV-A, the impact of this mismatch is minimized by the choice of the cost functions.

On the other hand, as explained in Section III-E, during a step transient, the SAC algorithm should behave like a fixed DFB if the thresholds are correctly chosen. Therefore, the required constraints are the same as those presented in [18].

V Numerical Results

For the sake of a better assessment of the results, the algorithm was implemented in a sliding-window fashion. The following parameter values are used in all simulations. The nominal frequency is set to $f_{0}=60$ Hz and the sampling period is $T=1/(32\cdot f_{0})\approx 520.83$ $\mu$s. The reporting rate is set to $F_{s}=60$ frames-per-second (fps). After a parameter tuning process, the thresholds are defined to be as follows: $\Delta\omega_{\min}=2\pi\cdot 10^{-3}$ rad/s, $\Delta\omega_{\max}=2\pi\cdot 15$ rad/s, $\Delta\sigma_{\min}=6\cdot 10^{-3}$ 1/s, $\Delta\sigma_{\max}=4$ 1/s, $\Delta\alpha_{\min}=2\pi\cdot 0.1$ rad/s2, $\Delta\alpha_{\max}=2\pi\cdot 16$ rad/s2, $\Delta\gamma_{\min}=0.6$ 1/s2, and $\Delta\gamma_{\max}=110$ 1/s2. The number of iterations is limited to $k_{\max}=5$. Due to aliasing effects, only $M=11$ harmonics are considered instead of 50 as in the IEEE Std. In any case, one should consider that in practice the anti-aliasing filter limits the harmonic content of the signal fed to the algorithm. To properly account for noise, note that the signal-to-noise ratio (SNR) was measured at the distribution level to be around 60 dB [23], while SNR at the transmission level is higher. However, it should be noted that this value was obtained with a sampling rate of 1024 samples per cycle. Assuming that noise is white its power may be reduced approximately by the factor $1024/32=32$ by using a proper decimator filter. Thus, we use $\textrm{SNR}=75$ dB, which is the expected SNR at the sampling rate considered.

The prototype filters length is set to two nominal cycles, which is particularly well suited for protection and control applications (i.e., for use in a P class PMU). Parameter values for their design are $\omega_{\textrm{con}}=2\pi\cdot 15$ rad/s, $\alpha_{\textrm{con}}=2\pi\cdot 16$ rad/s2, $\sigma_{\textrm{con}}=4$ 1/s, $\gamma_{\textrm{con}}=110$ 1/s2, $L_{\omega}=L_{\sigma}=0.3$ and $L_{\alpha}=L_{\gamma}=0.9$. The contractive ranges were chosen to ensure convergence over the most stringent range of parameters of the benchmark signals considered in the IEEE Std. However, $\omega_{\textrm{con}}$ is chosen so that the algorithm has a wide frequency range basin of attraction, motivated by other works [10, 11]. Note that the contractive constraints on the second-order variables are weaker, which gives looser constraints for the prototype CSIP filters design problem and therefore allows to obtain a smaller value for their norms. This translates to better performance in noisy conditions, which is a critical issue for ROCOF and ROCOD estimates. For the step constraints, the uniform scaling parameters are set to $\rho_{i}=0.5$ (see [18]).

The worst-case errors of the SAC algorithm for the different tests performed are summarized in Table I. In addition, as a benchmark, the same values are reported for the method proposed in the IEEE Std. for a P class PMU (from now on simply the IEEE Std. algorithm) [4, 5]. The step test results are presented in Table II as they involve different metrics of performance. For the sake of completeness, damping and ROCOD estimates were incorporated to the IEEE Std. algorithm using central finite differences analogous to the ones used for frequency and ROCOF estimation but applied to the estimate of the log-amplitude signal.

V-A Stationary Tests

V-A1 Off-nominal frequency

The system is tested with a sinusoidal input signal whose frequency is swept in the off-nominal range defined by the contractive set $\Omega_{\textrm{con}}$ with steps of 0.1 Hz (referred to as test A1). The results show a good performance of the SAC algorithm over the entire contractive frequency range. Moreover, it is interesting to observe that the values of $2\pi\cdot\textrm{FE}$ and DE are comparable, as well as those of $2\pi\cdot\textrm{RFE}$ and RDE. This could be useful to set tentative limits for these quantities.

V-A2 Harmonic distortion

The harmonic distortion rejection test is considered under two different scenarios, jointly with an off-nominal frequency condition in the entire frequency contractive range. In the first scenario, each harmonic is introduced separately with an amplitude equal to $1\%$ of the fundamental (test A2). In the second scenario, all harmonic amplitudes are simultaneously set to $10\%$ of the fundamental (test A3). Although this is an unrealistic condition, it is nevertheless useful to verify the SAC algorithm robustness to high harmonic distortion levels. In addition, it is interesting to note that in all of the above tests, the number of iterations $k$ was only nonzero at the initialization. Therefore, in steady-state conditions, the cost of the algorithm is virtually the same as that of a non-adaptive DFB.

V-B Dynamic Tests

V-B1 Frequency Ramp

A frequency ramp test is performed starting at $t=1$ s and ranging from 45 to 75 Hz with a ramp rate of 1 Hz/s (test B1). Maximum errors on this time interval (discarding exclusion intervals [5]) are presented in Table I. Again, due to the wide frequency range covered in the test, the performance difference of the methods is quite large.

V-B2 Amplitude and Phase Modulations

A joint amplitude and phase modulation test is performed using a modulation frequency of 2 Hz and modulation factors of 0.1 (test B2). The results of the SAC algorithm shown in Table I indicate that the model mismatch impact is not severe for such low-frequency modulations. In fact, once again, no adaptations occur after the initialization. On the other hand, the same test but for a 5 Hz modulation frequency is also considered (test B3). In this case, the model mismatch impact on the SAC algorithm is much more noticeable. In fact, adaptations occur on more than $67\%$ of the times, with 0.94 average iterations. Nevertheless, results are well within the IEEE Std. established limits [5].

V-B3 Amplitude and Phase Steps

Firstly, an amplitude step test with a step size of $0.1$ pu is performed. It should be noted that the DE and RDE response times have been defined by using the rough empirical correspondences $\textrm{DE}\sim 2\pi\cdot\textrm{FE}$ and $\textrm{RDE}\sim 2\pi\cdot\textrm{RFE}$ in order to assess their speed of response. Response times are denoted as $t_{\textrm{res},\textrm{TVE}},t_{\textrm{res},\textrm{FE}}$, etc. Secondly, a phase step test with a step size of $\pi/18$ rad is considered. The results in Table II show that the SAC algorithm adaptation scheme works successfully under fast changes. In fact, the performance is compliant with the prototype specifications. Note that the IEEE Std. algorithm RFE and RDE response times are reported as infinity since these error metrics are always greater than the limits established in the IEEE Std.

V-C Faults

In all fault conditions, the performance is reported in Table I after the response time allowed by the IEEE Std. (i.e., $2/f_{0}\approx 33.3$ ms for AE, PE and TVE, $4.5/f_{0}=75$ ms for DE and FE, and $6/f_{0}=100$ ms for RFE and RDE) in a step test.

V-C1 Small Power System

The two-area system of Fig. 3 is considered [24], where it is assumed that two PMUs are installed at buses 1 and 13. The simulation is performed with the Power System Toolbox (PST) [25] (using the d2asbeg.m file with default parameters) and numerical results are reported for the PMU at bus 1. Two situations are considered. First, a multiple loss of lines is applied at $t=0.1$ s between the buses 3 and 101 (test C1). Results in Table I show that in such severe case the algorithm will allow for a smooth transition between the two grid topologies. Second, a three-phase fault is applied at the same time and between the same buses, and it is cleared after 5 cycles (test C2). The evolution of the estimates of all the parameters is shown in Fig. 4, showing a close agreement between the true values and the estimates after a fast transient.

V-C2 Large Power System

A three-phase fault between buses 6 and 7 is simulated for the well-known IEEE 50-machine 149-bus system, and measurements at bus 6 are assessed (test C3). In this case, we observed that the dynamics are significantly more pronounced but correctly tracked by the SAC algorithm.

VI Conclusion

A novel algorithm for enhanced dynamic phasor estimation with a self-adaptive capability has been presented to cope with the stringent performance requirements found in different applications of PMUs measurements. There are several interesting features that make this algorithm a powerful and flexible tool. Firstly, it has a guaranteed convergence in a prescribed set of conditions that can be defined at the design stage. Secondly, it offers great harmonic rejection even with large frequency deviations. Thirdly, it offers a simple control on the trade-off between accuracy, speed and computational cost. In addition, it is proposed to estimate (and report) damping and ROCOD, which are more physically meaningful than plain amplitude derivatives. Thus, the knowledge of these quantities could lead to the design of novel amplitude control schemes. Simulations of the IEEE Std. tests under realistic noisy conditions as well as different faults have shown that remarkably accurate and fast estimates can be obtained with this approach. These results show that algorithm is promising for use in both transmission and distribution networks.

[Convergence Analysis]

For brevity, only the evolution of frequency is considered. However, the other quantities are treated in an exactly analogous manner. Assuming for simplicity that the adaptation is performed always, the prior frequency at iteration $k+1$ is simply $\omega_{\textrm{pr},k+1}=\omega_{\textrm{pr},k}+\widehat{\omega}_{Z_{k}}[n_{r}]$. To derive the convergence conditions, consider the signal phasor model ${X[n_{r}+m]}=X[n_{r}]e^{j\omega_{s}Tm}$, where $\omega_{s}\in\Omega_{\textrm{con}}$. Neglecting the influence on the other parameters for simplicity (i.e., considering a one-dimensional error flow), it is found that $e_{\omega,k}=\omega_{s}-\omega_{\textrm{pr},k}$ satisfies the recursive relation

[TABLE]

It is clear that $f_{\omega}$ has a fixed point at the origin (i.e., $f_{\omega}(0)=0$) since $A_{1}(0)=0$ for any $\mathbf{a}_{1}\in\mathbb{R}^{N/2}$. Therefore, convergence of the prior frequency $\omega_{\textrm{pr},k}$ to $\omega_{s}$, for any $\omega_{s}\in\Omega_{\textrm{con}}$, can be guaranteed if the following contractive constraint holds:

[TABLE]

This is actually an application of the celebrated contraction mapping theorem, which guarantees both uniqueness of the fixed point zero and convergence of the fixed-point iteration $e_{\omega,k+1}=f_{\omega}(e_{\omega,k})$ for all the desired conditions [26]. By design, from (18), condition (23) holds exactly for the prototype filters, that is, for the function $\widetilde{f}_{\omega}(\omega_{s})=\omega_{s}-\widetilde{A}_{1}(\omega_{s}T)/\widetilde{A}_{0}(\omega_{s}T)$. In addition, the filters adaptation strategy (15) changes the filter amplitude responses mainly around harmonic frequencies. Thus, the approximation $f_{\omega}(\omega_{s})\approx\widetilde{f}_{\omega}(\omega_{s})$ for all $\omega_{s}\in\Omega_{\textrm{con}}$ is reasonable for our analysis. This implies that $Z_{k}[n_{r}+m]\to X[n_{r}]$ for all $m\in\mathcal{W}$ as $k\to\infty$, that is, ${\boldsymbol{Z}}_{k}[n_{r}]$ converges to a constant signal. Therefore, the DFB ideal gain nominal frequency conditions $A_{0}(0)=1$ and $A_{1}(0)=0$ guarantee that $\widehat{\omega}_{X_{k}}[n_{r}]\to\omega_{s}$ as $k\to\infty$.