Constrained-Group-Delay-Optimized Multiple-Resonator-Based Harmonic Analysis

: Recently, a multiple-resonator structure was proposed as a robust and computationally efficient tool for harmonic analysis. Two trivial cases have been previously observed. The first case exhibits good out-of-band suppression and elimination of unwanted harmonics, but with a high latency. In the second case, a phase frequency response around the passband centre is zero-flat, that provides fast estimates. However, in this case, resonant peaks at the ends of the passband and high interharmonic gains cause large overshoots. In general, the basic algorithm performance requirements: selectivity and speed, are contradictory which makes it impossible to completely fulfil both of them. In this paper, the Constrained Linear Least-Squares (CLLS) optimization method is used to obtain a compromise solution. As a result, the resonant peaks in the passband are avoided and side lobes are mitigated, simultaneously minimizing the group delay in the middle of the passband. Performed simulations confirmed the effectiveness of the proposed algorithm. algorithm


INTRODUCTION
It is widely accepted that, in practice, the voltage and current harmonics are time variant due to the dynamic nature of linear and nonlinear loads and equipment in power systems and the continuous changes in system configurtion. The analysis of the time-varying characteristics of harmonic distortion is an important research topic for protection and control of power systems [1]. Based on this fact, there has been an increasing interest in signal processing techniques for detecting and estimating harmonic components of time parameters. The challenge is designing a harmonic estimator with high convergence ratio, high acuracy, low computational burden and immunity to the presence of interharmonics: conditions that are not easy to simultaneously deal with. In general, the design difficulty of these systems lies on the trade-off between a fast-dynamic response and filtering performance. In order to obtain a fast-dynamic response, a wide bandwidth is required. On the other hand, to filter interference signals and noise, bandwidth should be reduced as much as possible. This yields a tight constraint, which can be hardly fulfilled.
One of the new approaches to the estimation of harmonics of dynamic signals based on recursive filters realized using multiple resonators (MR) was proposed in [31]. This algorithm, in addition to the fundamental and harmonics phasors, can also estimate their first derivatives (the first, second and so on, depending on multiplicity of resonators). Two signals are present at the output of the parallel resonator structure for each harmonic. The first one is the output of the cascade, and it shows good attenuation of noise and harmonics, but it has a large latency in relation to the input signal. In the case of the second signal, which is a sum of outputs of all resonators belonging to the corresponding cascade, a resulted filter has a zero-flat phase response around the operation frequency and thus it provides a fast estimate, but with large overshoot caused by resonant peaks at the ends of the passband.
In [32], by expanding the estimation structure described in [31], a compromised solution is offered. This is made possible by using a linear combination of the outputs of the corresponding resonators in the cascade, instead of the sum. Both frequency responses, amplitude and phase, in the bandwidth range were taken into account. The group delay of the transfer function in the centre of the bandwidth was controlled simultaneously with a limitation of the peaks at the ends of the bandwidth. This is obtained through the shifting of the reference time point along the filters data window directly influencing the latency. The minimum value of the group delay for given conditions is determined heuristically.
A different design technique is proposed in this paper. The block diagram of the MR-based estimator is the same as in [32], but the proposed design method is different. Unlike a design method proposed in [32] which used the heuristic procedure, constrains on the amplitude frequency response and/or the group delay are prescribed in advance and an optimized frequency response is obtained corresponding to the chosen minimisation function. Design problems can be defined in different way for different purposes. Thus, on the basis of the basic MR structure, it is possible to obtain a whole bank of filters corresponding to a dozen of both stationary and dynamic conditions. Fig. 1 depicts a diagram of the triple-resonator-based estimation structure. A channel corresponding to the harmonic m consists of the cascade of triple resonators with pole z m = exp(jmω 1 ) at an angular frequency ω m = 2πmf 1 /f s (f 1 is the fundamental frequency, f s = 1/T is a sampling rate, T is a sampling period). Complex gains g m,k , m = −M, … , 0, …, M, k = 0, 1, 2, are assigned to the corresponding resonators in the cascade. All resonator cascades are embedded into the common feedback loop. Due to the infinite feedback loop gains at frequencies of the resonator poles, the transfer functions of each harmonic channel have unit values in own frequencies. This feature does not depend on other parameters of the system. Further, zero-flat gains around other harmonic frequencies are provided ensuring reinforcing of the required attenuations. The order of the overall closed system is 3(2M + 1). A maximum number of harmonics is M, that is equal to the number of resonator pairs, and depends on the fundamental component angular frequency the ω 1 = 2πf 1 /f s , because condition Mω 1 < π has to be satisfied. In [31], a closed form solution for the calculations of the resonator gains g m,1 , g m,2 and g m, 3 for a dead-beat estimator aregiven:

DESIGN METHOD
, and m, m, m, m, m, m, m m m Finally, it is:  It is supposed that a signal that will be analysed is defined by: The output of the k-th resonator (observed in the reverse order, starting from the end of the cascade and ending at its beginning) in the m-branch corresponding to the harmonic m is It consists of the complex envelope is also used.
The transfer functions corresponding to harmonic differentiator signals are as follows: Amplitude frequency responses of the transfer Eq. (6) to Eq. (8) for the fundamental component (m = 1) and the fifth harmonic (m = 5) are given in Fig. 2.
Two output signals belonging to the channel m are noleontable: is inherently formed and embedded in a common feedback signal. Corresponding transfer functions are , is a reason for huge gain peaks at the ends of the passband width. The first attempt to reduce overshoots is to weak this requirement for the flatness. It can be shown that for one of the particular cases, a signal  (for r m,1 = 1 and r m,2 = 1), the first derivative of the corresponding transfer function In this case, the narrower both amplitude and phase flatness around the estimating harmonic frequency are provided. A gain of the amplitude response at the ends of the passband is smaller than for it is still too high and large resonant frequency peaks are still presented, Fig. 3 (solid line).

Problem Statement
In addition to inherent signals formed as a linear combination of the zero, first and second derivative estimates with coefficients 1, r m,1 and r m,2 , respectively. The transfer function corresponding to the signal where: Eq. (10) has the following matrix form: where:

Linearization of Constraints
Eq. (12) are nonlinear. These constraints can be linearized using an appropriate approximation [33], allowing to solve the problem through a constrained linear least-squares (CLLS) optimization method [34]. In [33], an octagon has been used to approximate a circle. In this pa-per, an approximation of the circle is generalized to a polygon of any order providing an enough high approximation accuracy. It is easily proved to be valid: Figure 4 Approximation of a cycle with a square and an octagon Since i a is unknown, the following system of the linear inequalities can approximate nonlinear one (12): Using matrix notation, Eq. (15) becomes the following linear form: where matrix m A and vector m b are given by: At the end, the design problem can be cast to the typical optimization one: It should be mentioned that the linear programming (LP) technique can also be used.

DESIGN EXAMPLE
This section shows two design cases. The first case considers only limitation of the gain in the passband, in the absence of constrains in the stopbands. The value of the maximum gain in the passband is settled to 1.    T z filter yields an excessive bandwidth at the expense of a weaker stopband attenuation, which is detrimental to the estimator interference rejection capability. In the next section, simulations with two performance tests will be performed to analyse the influence of the chosen constrains, and consequently designed parameters, on tracking performance. Due to space savings, only frequency responses for the basic component (m = 1) are depicted. The frequency responses for other harmonics have the same shapes. It is important to note that, for asynchronous sampling, the parameters designed for one harmonic also apply to the others.

SIMULATION RESULTS
The aim of the following examples is only to demonstrate the ability of the algorithm to track the fundamental and harmonic components, but not to demonstrate a compliance with any standard. Tests have been performed by computer simulations. A sampling frequency of f s = 800 Hz was adopted. A multiple harmonic content is injected to the signal. Odd harmonics (up to seventh) were added simultaneously to the signal waveform, with the individual harmonic levels 30%, 10% and 10%, respectively. Even harmonics are set to zero. Parameters (coefficients r m,1 and r m,2 corresponding group delay) of all considered cases are summarized in Tab. 1.
The estimates of the amplitude, phase, and TVE in per- with f 1 = 50 Hz and f s = 800 Hz, for the phase angle modulated and step changes, are shown in Fig. 6 and Fig. 7. Similar performances were obtained for changes of the phasor amplitude too. Given the brevity, estimates for the basic component are shown only. The waveforms of estimates for harmonic phasors have similar shapes.
A total vector error (TVE) is defined by Eq. (18), and it is used as an index of cumulative information, which aggregates three possible sources of error: amplitude, phase angle, and timing [2]. . The TVE is thus a convenient mean to express the accuracy of measurement. This representation offers the advantage of being extremely compact, but it is often more useful to refer back to individual errors components, to understand better the impact of measurement accuracy on the specific application [2]. The first test is referred to the phase modulation of the input signal, reflecting phase angle variations that can occur in the voltage waveforms during a swing of the stable power. The standard [35] proposes the following test signal: where X m is an amplitude of the fundamental component,  T z need to be improved in the step test, but they are better in terms of reducing the group delay and improving tracking of the modulated signals.
It is obvious that a wider passband provides a lower latency, but resulted overshoots for inputs with step changes are much larger. And vice versa, a reduction in bandwidth provides lower overshoots for step-changed input responses and longer response times and group delays. The bandwidth modification allows a trade-off between these opposite requirements.

CONCLUSION
This paper studies the spectral estimation based on the MR filters presented in [31]. In this paper, the CLLS optimization technique is applied to find a compromised solution between the conflicting requirements arising from the need for instantaneous estimation, which is achieved by sufficient both amplitude and phase frequency responses flatness in the bandwidth centre, and admissible levels of peaks at the bandwidth ends. The results of the conducted computer simulations showed that the obtained performances of this estimation algorithm are in the middle compared to the best and the worst case obtained for two extreme inherent cases. Thus, the obtained estimators can be at the same time sufficiently accurate and resistant to the harmonics and/or interharmonics and the presence of noise, and yet fast enough to enable monitoring of the content of nonstationary signals. It is worth noting the possibility of obtaining a set of filters based on the core MR-structure, to simultaneously estimate different harmonic phasor with different estimation performances. For example, the fundamental component and dominant low order harmonics can be estimated with a smaller group delay and faster response, since latency may be more important than accuracy for feedback control in many applications such as wide-area control and protection, while low level harmonic phasors can be estimated with filters with a narrow bandwidth and a slow response. Even more, the same harmonic phasor can be estimated with different dynamics and so obtained estimates further be used for different purposes.
Further works will be focused to the compliance of this estimation technique with the relevant standards.