Adaptive Polynomial Harmonic Distortion Compensation in Current and Voltage Transformers Through Iteratively Updated QR Factorization

Measuring current and voltage harmonics has paramount importance for improving the power quality of distribution grids. However, the achieved accuracy strongly depends on the adopted instrument transformer (IT). This article proposes an adaptive technique that enables an effective compensation of both the filtering behavior and the harmonic distortion (HD) introduced by current and voltage transformers (VTs), namely the strongest nonlinear effect at low-order harmonics. The approach is based on a flexible, linear in the parameters polynomial modeling of HD in the frequency domain. Model complexity can be different from one harmonic to the other, and it is selected through an automatic iterative process to suit the nonlinear behavior at each specific harmonic order, while avoiding overfitting. In particular, the number of parameters is increased by progressively updating the QR factorization of the regressor matrix trough Householder reflections until a convergence condition is reached. Experimental tests performed on an inductive VT and current transformer (CT) highlight the effectiveness of the approach.

Abstract-Measuring current and voltage harmonics has paramount importance for improving the power quality of distribution grids. However, the achieved accuracy strongly depends on the adopted instrument transformer (IT). This article proposes an adaptive technique that enables an effective compensation of both the filtering behavior and the harmonic distortion (HD) introduced by current and voltage transformers (VTs), namely the strongest nonlinear effect at low-order harmonics. The approach is based on a flexible, linear in the parameters polynomial modeling of HD in the frequency domain. Model complexity can be different from one harmonic to the other, and it is selected through an automatic iterative process to suit the nonlinear behavior at each specific harmonic order, while avoiding overfitting. In particular, the number of parameters is increased by progressively updating the QR factorization of the regressor matrix trough Householder reflections until a convergence condition is reached. Experimental tests performed on an inductive VT and current transformer (CT) highlight the effectiveness of the approach.

I. INTRODUCTION
M ODERN distribution grids (DGs) are characterized by a large penetration of nonlinear loads and generators, typically based on power electronics. This has noticeably raised the harmonic content of current and voltage waveforms [1], which increases heating, wear and stress of dielectric materials, it reduces overall efficiency, as well as it may trigger resonant modes.
Monitoring current and voltage harmonics [2] enables analyzing their propagation on the network, thus recognizing critical situations. In conjunction with proper signal processing techniques [3], [4], [5], [6], it would be possible to identify the users that are the most responsible for harmonic pollution. Moreover, in principle the large-scale availability of harmonic measurements would enable reconstructing the harmonic state of the grid [7], which helps planning actions aimed at improving power quality.
Obtained results are strongly related to the quality of the measurement data, expressed in terms of uncertainty. When estimating harmonics, two different uncertainty sources can be identified [8]. The first is the signal processing algorithm adopted to extract the harmonic components; the second is related to the hardware used to acquire those waveforms. In this respect, the most significant contribution comes from the instrument transformers (ITs) [9]. Among them, inductive voltage transformers (VTs) [10] and current transformers (CTs) [11] are still widely used, thanks to their favorable blend of ruggedness, metrological performance, and long-term stability.
Conventional inductive ITs are often employed for harmonic measurements [12], [13], [14], albeit their metrological performances are not guaranteed. They are typically regarded as affected by severe bandwidth limitations [15], [16], but the iron core makes that VTs [17], [18] and CTs also suffer from a nonlinear input-output relationship. In particular, the saturation of CTs produced by the dc transient component have been widely studied, and techniques aimed at mitigating its effect can be found in the literature [19]. However, the impact of core nonlinearity on harmonic measurements have been investigated only recently [20], [21], [22].
Several frequency-domain methods have been developed with the aim of reconstructing primary harmonics from the secondary side of an IT, while mitigating nonlinear effects [23], [24]. To reduce complexity, they exploit the typical spectral content of voltage and current waveforms, which consists of a strong fundamental and superimposed harmonics having significantly smaller amplitudes. It is worth noting that low-order harmonic measurements are the most heavily affected by nonlinearity. More specifically, in this case, the strongest nonlinear phenomenon is the harmonic distortion (HD) produced by the prevailing fundamental tone. For this reason, simplified methods that enable compensating only the HD (thus neglecting the intermodulation between different primary components) have been proposed [25], [26], [27].
In this respect, the authors of this article have developed a method for reconstructing primary harmonics through polynomial modeling and compensation of HD. It has been applied to This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ inductive VTs [25], [26] and CTs [28], but the method, based on behavioral modeling, is general and thus it abstracts from the physical operating principle. One of the main advantages is that identifying the parameters of the compensation formulas does not require a calibrator, capable of injecting a specific set of training waveforms with high accuracy. On the other hand, the underlying model of HD has a stiff structure, based on fixed-degree polynomials, the same for all the harmonics. The main contribution of [29] is overcoming this limitation: one may select, for each component, a different number of terms to compensate for HD. The improved method has been applied to a VT and CT through numerical simulations.
This article represents the technical extension of [29], thus the starting point is the same flexible structure for the HD compensation formulas. However, the number of terms at each harmonic order is not selected by the user, but it is the result of an adaptive identification procedure looking for the best tradeoff between accuracy of the reconstruction and complexity, while avoiding overfitting. This has key importance for a large-scale, automated implementation of HD compensation. Moreover, identification requires solving a least squares (LS) problem whose condition number becomes worse as complexity is increased; QR factorization is adopted to guarantee high numerical stability. In addition, the adaptive algorithm progressively updates the QR factorization through Householder reflections, thus improving computational efficiency since the identification problem is not solved from scratch at every iteration. Finally, the proposed technique has been experimentally applied to an inductive VT and CT. Obtained results show how the adaptive HD compensation approach provides an effective mitigation of nonlinearities by selecting the number of terms of the polynomial so that it reflects the impact of HD at each specific harmonic. In comparison, the former HD compensation method based on fixed-degree polynomials may result in lower accuracy, overparametrization, or numerical problems during identification.

A. Volterra Systems and ITs
ITs installed in ac power systems can be considered as (weakly) nonlinear time-invariant systems, where input and output are represented by the primary and secondary side quantities, x 1 (t) and x 2 (t). Furthermore, let us suppose that x 1 (t) is a periodic waveform at rated frequency f 0 , thus fully defined by the fundamental and the set of harmonics. Neglecting complex, chaotic nonlinearities, the corresponding steady-state response x 2 (t) is periodic with the same period. In this case, the IT behavior can be more effectively studied in the frequency domain, where it corresponds to a mapping rule between primary and secondary side spectral components. Under mild conditions, this mapping can be expressed using a frequency-domain Volterra [30] (or polynomial [31]) model. Adopting this approach, the mth order secondary-side harmonic (m ≥ 1) is where X 1 (n) is the nth-order component in the two-sided primary spectrum. Since x 1 (t) is real-valued, X 1 (n) is the complex conjugate of X 1 (−n). According to (1), X 2 (m) is the sum of infinite contributions X 2(d) (m), each produced by a dth degree homogeneous subsystem; in particular, that having d = 1 represents the underlying linear system. In turn, the output of each dth degree subsystem is a linear combination of all the products between sets of d input components (with repetition) whose sum of the harmonic indexes equals m. The input-output behavior is defined by the set of coefficients H (d) , which corresponds to the dth degree generalized frequency response function (GFRF) evaluated on a harmonic grid. The target is measuring primary harmonics, namely reconstructing their values from the secondary side of the IT: this requires somehow obtaining the post-inverse of (1). It can be proven that both the pre-and post-inverses of a Volterra system are still Volterra systems [32]. Therefore, separating the contribution produced by the underlying linear system where the coefficients K (d) are somehow related to H (d) in (1). Under the previous assumption, knowing the post-inverse would enable retrieving power systems harmonics from the IT secondary spectrum, thus compensating both its filtering and nonlinear behavior.

B. Polynomial HD Compensation
Unfortunately, (2) cannot be implemented, being it defined by an infinite number of coefficients K (d) . Therefore, it is necessary to introduce simplifications that allow obtaining an expression having manageable complexity, but still providing an acceptable reconstruction of harmonics. The trivial idea is upper bounding to D the maximum degree of the nonlinear contributions, obtaining a truncated Volterra system. However, this is an ineffective approach: choosing small values of D undermines the capability to compensate for nonlinearities, but it still leads to a huge number of coefficients.
In addition, it is possible to exploit the typical spectral content of waveforms in power systems, which are quasi sinusoidal, namely made of a strong fundamental and significantly smaller harmonics, usually by at least an order of magnitude. Thanks to the weakly nonlinear behavior of ITs, this feature applies also to the secondary-side spectrum. Starting from this observation, one may neglect the contributions from products containing more than a harmonic component [23]. To further reduce complexity, it is possible to consider in (2) only the nonlinear contributions involving exclusively the fundamental and its image. In this case that is the input-output relationship of a subclass of pruned Volterra systems. As a result, (3) allows compensating the filtering behavior and the HD produced by the large fundamental tone that, under the aforementioned assumptions, represents the strongest nonlinear effect in ITs. On the contrary, (3) is not able to address other nonlinear phenomena, in particular the interaction (intermodulation) between different primary components. Equation (3) can be rewritten in a more compact form, which is Looking at (4) it is worth noting that considering the mth order harmonic, HD is compensated through a univariate polynomial function in the fundamental magnitude; the degree of each monomial is d, which is the degree of the corresponding nonlinear contribution. When considering an even order harmonic, the polynomial function just contains even powers greater or equal than m, up to D. Conversely, only odd powers between m and D appear in the polynomial functions that enable compensating HD at odd order harmonics. Equation (4) applies also for m = 1, thus mitigating the (typically very small) nonlinearity affecting the fundamental component measurement.

III. ADAPTIVE POLYNOMIAL HD COMPENSATION A. Improved Polynomial HD Compensation
According to Section II-B, (4) enables reconstructing the primary side harmonics in an IT. Both its filtering behavior and the strongest nonlinear effect are mitigated through univariate polynomial functions having degree D. The number of terms (powers of the fundamental magnitude) that are combined to compensate for HD affecting the mth order harmonic is Therefore, it reduces with the harmonic order, becoming zero if m > D. The first consequence is that mitigating nonlinearity at a given harmonic order demands a minimum value of D. From the opposite point of view, distortion at lowerorder harmonics is tackled with a progressively increasing number of terms, but some of them may be not significant. For example, this occurs when a predetermined accuracy target at that harmonic is reached also with a lower value of D. It may also happen that, for some harmonic orders, increasing D does not enable improving the quality of the reconstruction. The reason could be that HD has been reduced so that its impact has become negligible with respect to other uncertainty sources, such as noise or intermodulation between different spectral components in the primary waveform.
In order to increase the flexibility of HD compensation, it is possible to individually select L(m) [29], namely the number of monomials that allow modeling HD at each specific harmonic order. According to this approach, primary harmonics are obtained as . Before using (6) to retrieve a generic mth order primary harmonic, two passages are required. The first one, discussed in Section III-B, is presenting a numerically stable method for identifying the coefficients K m and B m for a predetermined number of terms L(m) (sometimes denoted as L for the sake of a lighter notation). Estimates are obtained by solving an LS problem through the QR factorization of the regressor matrix. The second one, covered by Section III-D, is defining an adaptive procedure to choose the value of L ensuring a proper compensation of HD. It is based on the computation of the QR factorization through Householder reflections, as explained in Section III-C.

B. Robust Identification
The compensation formula (6) is linear in the unknowns. Therefore, feeding the IT with a set of S quasi-sinusoidal periodic waveforms while measuring the corresponding steadystate response, permits writing a system of S linear equations in the form where ε (L) (m) is the column vector of the residuals, namely the deviations between the mth order harmonics in the applied signals (components of vector x 1 (m)) and their reconstructions from the secondary side using the parameters in the vector Finally, Y (L) m is an S × (L + 1) matrix that can be derived from the structure of (6) and the considered set of signals.
m has full column rank, an estimatê p (L) m of p (L) m is obtained in closed form by solving the ordinary LS problemp with || · || denoting the Euclidean norm. It is worth noting that the more the set of the identification signals forms a statistically significant sample of those found during regular operation, the more solving (9) corresponds to minimizing the second-order moment of the distance in the complex plane between the reconstructed and the actual primary side mth order harmonic phasor. The features of the LS problem depend on the regressor matrix, which can be partitioned as where x 2 (m) is the column vector of the secondary side mth order harmonics in the identification signals, F m is an S × S diagonal matrix whose sth entry is with X [s] 2 (m) being the mth order secondary component in the sth identification signal. Introducing • as the Hadamard (element-wise) power operator, U (L) m reads as that is the Hadamard square of a rectangular Vandermonde matrix. This is somehow expected, since estimating p (L) m involves polynomial fitting [33]. Unfortunately, because of the peculiar structure, the condition number rapidly increases with L. As a consequence, numerical instability may occur if (9) is solved through explicit inversion of the corresponding Gram matrix. A first expedient to mitigate this problem is adopting per unit values, so that the magnitude of the fundamental is not much larger than unity. Numerical robustness is highly enhanced if estimates are obtained via QR factorization, thus decomposing Y (L) m into the product between a unitary S × S matrix Q (L) m and an upper triangular S × (L + 1) matrix R (L) m . The reduced (or thin) QR factorization [34] is often adopted, thus leading to where m,1 consists of the first L + 1 rows of R (L) m . This enables writing andp (L) m is obtained through back-substitution (avoiding matrix inversion) thanks to the upper-triangular structure of R (L) m,1 .

C. QR Factorization Through Householder Reflections
Several methods can be adopted to compute the QR factorization of matrix Y (L) m , but Householder reflections (or transformations) [35] are often employed thanks to the high numerical stability. A Householder transformation is defined by a column vector u (Householder vector), and it represents a (complex-valued) reflection in the hyperplane normal to this vector. It can be written as a premultiplication by the unitary matrix I−uu H (namely a rank-1 update of the identity I) that does not need to be explicitly computed or stored. Thanks to a proper choice of u, a Householder transformation applied to a vector v having the same size as u enables zeroing all the components but the first, which becomes equal to ν. Exploiting this property, R (L) m is obtained through L + 1 modifications of Y (L) m ; the kth consists in applying a Householder reflection to the last S − k + 1 rows and L − k + 2 columns, so that all the entries in the kth column below the main diagonal are zeroed. The ordered set of the associated Householder vectors {u k } k∈{1,...,L+1} is a factored representation of (Q (L) m ) H , which thus corresponds to a cascade of Householder reflections. Moreover, premultiplication by Q (L) m means applying the Householder reflections identified by the same vectors {u k } k∈{1,...,L+1} but in reversed order.

D. Adaptive Procedure
The next problem to be faced is selecting L(m) so that it enables reaching a satisfactory accuracy, without running into overfitting. In this respect, the overall accuracy in reconstructing the mth order harmonic can be quantified in terms of normalized root mean square error (NRMSE) is the vector of the reconstructed mth order primary harmonics of the S identification signals.
A straightforward idea to choose L(m) is progressively increasing its value starting from 0 (linear reconstruction), estimating the parameters of the corresponding compensation formula through (14) while computing the NRMSE from (15). The iterative process is stopped when one of the following conditions occurs.
1) NRMSE (L) (m) falls below a threshold NRMSE t (m). Its value should be set according to the optimal accuracy target, so that further improving the reconstruction of primary harmonics does not produce practical benefits.
2) The number of nonlinear terms reaches the maximum value L max (m), which bounds the complexity of the reconstruction formula.
3) The drop of NRMSE due to the added nonlinear term is below a predetermined threshold NRMSE t (m). Its value represents the minimum accuracy improvement that is considered as significant. The target of this stop condition is avoiding overfitting; when triggered, L(m) should be decreased by one. Considering the procedure summarized in Section III-C to perform the QR decomposition required by (14), it is possible to derive a computationally effective adaptive process for selecting the optimal value of L. Let us suppose thatL is the number of terms currently used to consider HD, thus corresponding to theL + 1th iteration step. From (10) and applying a further Householder reflection defined by vector uL +1 to zero its last S −L − 1 rows. It can be easily shown that (Q (L) m ) H is obtained by applying the Householder reflection defined by uL +1 to the last S−L rows of (Q (L−1) m ) H . Therefore, computing all the QR decompositions needed by theL + 1 iterative steps as progressive rank-1 updates basically requires performing the same number of Householder transformations that allows obtaining the QR decomposition of Y (L) m from scratch. Using (13) and (14) into (15), the NRMSE at theL + 1-th iteration step is where I S×S is the S × S identity matrix. Equation (16)  In conclusion, the adaptive procedure to find the optimal value of L(m) for each harmonic order is summarized by the pseudocode reported in Algorithm 1. Function build_y(x 1 , m, L) builds the last column of Y (L) m , house_vect(x) computes the Householder vector u that allows zeroing all the components in x but the first, which becomes equal to ν, also returned as output. house_refl(A, u) applies to A the Householder reflection defined by u, while compute_NRMSE(x 1 (m), Q (L) m,1 ) implements (16). Finally, function backward(A, b) solves a linear system of equations in the form Ax = b (with A upper triangular) via back substitution and 0 k×l denotes the k × l null matrix.

IV. APPLICATION TO INDUCTIVE VTS
The proposed adaptive polynomial HD compensation algorithm has been first tested for the reconstruction of primary voltage harmonics from the secondary side of an inductive VT. The first step is defining the class E V of primary voltage multisine waveforms to be adopted, which should resemble those found in distribution grids, characterized by quasi-sinusoidal spectral content. The fundamental frequency is assumed to be equal to its rated value, while spectral components have been treated as random variables. Fundamental magnitude has been assumed to be uniformly distributed between 80% and 120% of the rated primary voltage, corresponding to the measuring range of VTs [10]. As far as relative harmonic magnitudes (orders up to 25), independent rectangular probability density functions (pdfs) between 0.2% and 2% have been considered. All the phases have been supposed to be independent and uniformly distributed between −π and π.

A. Experimental Setup
The developed method has been applied to a class 0.5 VT having ratio K n = 200 V/100 V, rated frequency f n = 50 Hz and 20 VA burden. This requires an experimental setup capable of supplying the VT under test with waveforms that belongs to the previously defined class, while measuring its secondary spectrum; the architecture reported in Fig. 1 (resembling that presented in [36]) has been adopted.
The secondary side of the VT under test has been connected to the rated burden, while its primary winding has been supplied by an AETechron 7548 power amplifier through a 100 V/400 V boost transformer to increase its output voltage capability. The input signal of the power amplifier has been provided by a National Instruments NI USB 6356 BNC board, featuring synchronous generation and acquisition. Primary and secondary voltages v 1 and v 2 have been acquired using calibrated resistive dividers connected to two input channels of the data acquisition board (sampling rate f s = 100 kHz) through Analog Devices AD215BY isolation amplifiers in voltage follower configuration; their frequency responses have been measured and calibrated. It is worth noting that v 1 is not just a scaled-up replica of the signal v g applied to the input of the amplifier, mostly because of the voltage boost transformer. For a better accuracy in voltage generation, the small-signal frequency response between v g and v 1 has been measured and used to pre-compensate the reference waveforms to be applied to the VT under test.

B. Experimental Results
The first step to be performed is collecting the data required for identifying the parameters of the compensation formulas. For this purpose, a set of S = 100 identification signals belonging to class E V has been obtained by sampling the corresponding pdfs. They have been applied to the VT under test using the previously described setup; for each identification signal, once steady state operation has been reached, P = 100 periods of the primary and secondary voltages have been acquired. Thanks to synchronized sampling and generation, since f s / f n is an integer and high enough to avoid aliasing, it is possible to extract the harmonics in each period through discrete Fourier transform. Their values have been averaged over the P periods to mitigate the effect of noise and disturbances.
The measured primary and secondary voltage harmonics have been used, together with (13) and (14), to estimate the parameters of the HD compensation formulas, with L ranging from 0 (corresponding to the best linear approximation, BLA) to 5 and for harmonic orders up to the 11th. The proposed adaptive technique introduced in Section III-D has been implemented in MATLAB and executed with L max = 5, NRMSE t = 2·10 −4 , NRMSE t = 3·10 −5 , identical for all the harmonic orders. As a rough quantification of the computational burden, the adaptive process runs in about 40 ms on a core of a highend PC.
For evaluating performance, an independent set of 500 validation signals belonging to class E V have been synthesized and applied to the VT under test. Considering the sth signal, the resulting steady-state primary and secondary side harmonics V [s] 1 (m) and V [s] 2 (m) have been measured with the same method previously adopted for the identification signals. Accuracy in harmonic measurements is often quantified in terms of ratio and phase errors. However, their distributions over the validation signals are expected to be extremely similar, as a result of the combination of the nonlinear behavior of the IT and the independent rectangular pdfs (between −π and π ) of the harmonic phases, as discussed in [13] and highlighted by the results reported in [25] and [28]. Therefore, the total vector error (TVE) can be adopted as a concise indicator, capable of taking into account both ratio and phase errors simultaneously. Considering the sth validation waveform and the mth order harmonic it is whereV [s] 1 (m) is the reconstructed primary harmonic obtained with (6) and the previously estimated parameters for different values of L, or using the nominal ratio K n .
Obtained results are summarized in Fig. 2. For each harmonic order and reconstruction method, it reports both the rms (TVE rms , solid bar) and the 95th percentile value (TVE 95 , error bar) of the TVE computed over the set validation signals; arrows highlight the number L of terms selected by the adaptive procedure.
When considering the fundamental, using the nominal ratio produces a TVE 95 value of about 0.5%, with the corresponding rms value being just slightly smaller. As aforementioned, polynomial HD compensation with L = 0 corresponds to the BLA, namely to the complex-valued ratio ensuring the best reconstruction of a spectral component (in the LS sense) for a given class of signals. In this case, the 95th percentile value of the TVE drops below 0.03%. Adding nonlinear terms, thus L ≥ 1, produces a further accuracy improvement (TVE 95 becomes smaller than 5·10 −3 %), albeit not significant in most of the practical applications. As from Fig. 2, the adaptive identification procedure selects L = 0: the iterative process stops because the NRMSE (strongly related to the accuracy) is already below the threshold value.
As far as the second-order harmonic, using the BLA produces a minor accuracy improvement with respect to the nominal ratio. This means that nonlinearity is the main performance bottleneck in this case. Including up to two nonlinear terms leads to a progressive decrease of TVE 95 , from 1.5% to about 0.2%. Further increasing the complexity of the reconstruction formula does not improve performance. On the contrary, a small increase of TVE is noticeable. This is clearly due to overfitting: the compensation formulas include non-significant terms that lead to unnecessary complexity. Even worse, they may also reduce performance when moving outside the set of identification signals. During adaptive identification, the iterative procedure stops after trying L = 3: the reduction of the NRMSE is below the threshold, therefore L = 2 is selected, which represents the optimal choice.
It is well-known that an inductive VT mainly introduces odd nonlinearity: considering its Volterra representation (1), the most significant output contributions come from odddegree homogeneous subsystems. The consequence is that the impact of HD is substantially weaker at even-order harmonics. This is clearly noticeable from the results: at the fourth-order harmonic, HD compensation allows a small error reduction with respect to the BLA, while being virtually negligible for higher order, odd harmonics. The adaptive method selects L = 1 for the fourth and sixth-order harmonics, L = 0 for the eighth and the tenth.
The third-order harmonic is typically the most heavily affected by nonlinearity. In this case, adopting K n leads to TVE 95 above 13%, but almost the same value is obtained by using the BLA. Representing HD with just one term reduces TVE 95 by an order of magnitude. With L = 4, the 95th percentile TVE drops to 0.25%, while L = 5 produces basically the same result: HD has been already made negligible with respect to other nonlinear phenomena, such as intermodulation. Also in this case, the adaptive method selects the best value L = 4, the iterative identification procedure is terminated because the decrease of NRMSE is below the corresponding threshold.
Considering the fifth, seventh, and ninth-order harmonics, using nonlinear reconstruction formulas produces significant error reductions; however, accuracy improvement decreases with the harmonic order, since also the impact of HD becomes progressively weaker. In all the cases, TVE stops decreasing before L reaches 5. Nevertheless, TVE 95 is reduced by more than a factor 10, 3 and 2. The adaptive procedure selects L = 3, 2 and 3, respectively, always exiting the iterative process for the small decrease of NRMSE.
At the 11th-order harmonic, polynomial HD compensation still enables a better performance: TVE 95 passes from 0.85% (K n ) to 0.35% (L = 5), while 0.5% is obtained with the BLA. The adaptive method selects L = 3 (TVE 95 = 0.37%) in this case; including this contribution with the conventional HD compensation approach would have resulted in a 15th degree model, thus leading to many unnecessary nonlinear terms in the compensation formulas, especially at even-order harmonics. This highlights the benefits of the new flexible structure of the compensation formula (6), in particular when combined with the adaptive identification: it enables automatically finding the best tradeoff between accuracy and complexity, without stepping into overparametrization. Finally, it is worth highlighting that choosing the thresholds NRMSE t and NRMSE t does not represent a critical task. Having selected reasonable values, small variations may just slightly change the number of terms selected by the adaptive procedure.

V. APPLICATION TO INDUCTIVE CTS
The effectiveness of the presented adaptive polynomial HD compensation technique has been assessed also when employed for reconstructing the primary current harmonics from the secondary of a CT. The steps resemble those followed during the application to the VT, but there are some important differences. A new class E I of periodic multisine primary current waveforms has been introduced: it should reflect both the wider measurement range of CTs with respect to VTs, as well as the stronger harmonic content of currents. Therefore, a uniform distribution between 20% and 120% of the nominal value has been assumed for the fundamental magnitude. Relative harmonic amplitudes (orders up to 25 also in this case) have independent and identical rectangular pdfs ranging from 0.5% to 5%. It is worth noting that according to the relevant standards, some current components absorbed by loads may exceed 5%; on the other hand, measuring stronger harmonics is notably easier. Finally, phases are characterized by independent and uniform distributions in the range [−π, π ).

A. Experimental Setup
The developed technique has been applied to a class 0.5 CT with 10 VA burden, 50 Hz rated frequency, and nominal ratio K n = 100 A/5 A. An experimental setup conceptually similar to that described in Section IV-A has been adopted to apply the periodic multisine current waveforms to the CT under test [37]; its diagram is shown in Fig. 3.
In this case, the AETechron 7548 power amplifier feeds the primary of the CT under test (loaded with rated burden) through a 5 A/200 A current boost transformer. Primary and secondary currents (i 1 and i 2 ) have been measured with calibrated class 0.2 coaxial shunts having 20 kHz bandwidth and rated currents of 100 and 10 A, respectively. Their outputs have been acquired using calibrated Analog Devices AD215BY isolation amplifiers in noninverting configuration (100 nominal gain) connected to two input channels of the NI USB 6356 BNC board, operating with f s = 100 kHz sampling rate. The small-signal frequency response between v g and i 1 has been measured and employed to pre-compensate the reference waveforms to be applied to the CT under test.

B. Experimental Results
In order to identify the reconstruction formulas, S = 100 random current waveforms belonging to class E I have been applied to the CT under test. For each signal, P = 100 periods of the primary and secondary currents (i 1 and i 2 ) have been acquired under steady-state conditions. Primary and secondary harmonics have been computed using the procedure described in Section IV-B. Starting from these data, the parameters of the compensation formulas (L from 0 to 15, harmonics up to the 11th order) have been estimated. It is worth highlighting that such large values of L could be selected thanks to the robust procedure based on QR decomposition. The proposed adaptive identification method has been executed with L max = 15, NRMSE t = 1·10 −4 , NRMSE t = 8·10 −6 for all the considered harmonics.
Validation data have been gathered by applying another set of 500 current waveforms to the CT under test (belonging to class E I ), while measuring the corresponding primary and secondary harmonics, also in this case extracted by averaging over the P acquired periods. Accuracy in reconstructing primary current components has been evaluated in terms of TVE, defined by (17) with the straightforward substitution of voltage phasors with current phasors. Its rms and 95th percentile values have been computed over the set of validation signals; results are summarized in Fig. 4.
Starting from the fundamental, using the nominal ratio TVE 95 exceeds 0.5%, but it drops below 0.3% thanks to the BLA. This value is substantially larger if compared with that obtained with the VT: the reason is related to the wider variability of the current amplitude, which magnifies the impact of nonlinearity at the fundamental. When adding nonlinear terms to the reconstruction formula, the TVE progressively decreases, until L = 8 is reached. This behavior, caused by the large measuring range of CTs, is different with respect to that observed with the VT, where passing to L = 1 results in a huge decrease of the TVE, but adding more than two nonlinear terms does not improve accuracy. As far as the CT, the adaptive algorithm selects L = 3 (TVE 95 = 0.05%), exiting the iterative process because the NRMSE is below the threshold value.
When reconstructing the second-order harmonic using K n , TVE 95 is around 1%, and it does not reduce noticeably when adopting the BLA, thus confirming that nonlinearity is the major uncertainty source in this case. Accuracy in terms of TVE rms improves noticeably up to L = 5 (the corresponding TVE 95 is about 0.4%), while further increasing complexity produces negligible benefits; moreover, L ≥ 11 results in slightly higher errors, which is an indication that overfitting has occurred. The adaptive identification method correctly selects L = 5, the optimal choice in this case. The iterative process ends because adding another nonlinear term produces a drop of NRMSE below the threshold, thus detecting possible overparametrization.
As typical for inductive ITs, the strongest nonlinearity occurs at the third-order harmonic. Reconstruction with the nominal ratio produces TVE 95 above 3%, not mitigated by the BLA. On the contrary, adding more and more nonlinear terms leads to a slow decrease of the error, with L = 12 corresponding to TVE 95 below 0.4%; adding further terms, does not improve the estimates. The behavior differs from that observed with the VT, where L = 1 produces an abrupt decrease of the TVE; the reason is once again related to the wider variability of the fundamental current with respect to the voltage. The adaptive identification algorithm selects L = 12. It is important to highlight that including all these terms at the third-order harmonic requires at least a 25th degree compensation model (L = 15 would have required a 31st degree model) with the old truncation approach. Its inherently stiff structure would have resulted in many useless nonlinear terms at the other harmonics, thus producing overfitting.
To better understand how L affects performance, the range of variation of the fundamental magnitude has been divided into ten uniform intervals. The validation signals have been categorized accordingly, the rms and 95th percentile values of the TVE have been computed for each class. Fig. 5 reports the obtained values for the third-order harmonic. It is evident that introducing up to five terms to model HD enables a consistent accuracy improvement when the fundamental magnitude is above 40% of its rated value; errors are significantly larger at smaller fundamental amplitudes. A further increase of L enables improving performance in this range, while accuracy at larger fundamental magnitudes is barely affected. However, going beyond L = 12 increases errors in the largest magnitude class.
When considering higher-order, even harmonics, as expected HD is rather small and decreases with the harmonic order. For example, at the eighth-and tenth-order components, error is mitigated by passing from K n to the BLA, but adding nonlinear terms does not significantly refine the estimates. It is worth noting that for high values of L, TVE rms may increase noticeably: it is the clear indication that the compensation formulas suffer from overparametrization. The adaptive method selects a progressively smaller number of terms to model and compensate HD: L = 3 at the fourth-order harmonic, 2 at the sixth, 1 at the eighth, and 0 at the tenth, thus corresponding to a linear reconstruction. In all cases, the iterative procedure stops because the decrease of the NRMSE is smaller than NRMSE t .
As expected, HD is much more significant at odd-order components, and it reduces as the harmonic order increases. Estimates of the seventh, ninth, and 11th order harmonics may be jeopardized by overfitting if too large values of L are selected. In agreement with the previous considerations, the number of nonlinear terms in the compensation formulas chosen by adaptive identification becomes progressively smaller when increasing the harmonic order. In particular, it selects 11 terms at the fifth-order harmonic (not far from L max ), 3 at the seventh, 1 at the ninth, and 0 at the 11th, thus meaning that HD has become negligible.

VI. CONCLUSION
Nonlinearity, and in particular, HD introduced by ITs has a detrimental effect on the accuracy of low-order harmonic measurements. However, its impact could be very different according to the harmonic order. This article proposes a flexible, adaptive method aimed at compensating HD using polynomials in the frequency domain. The number of terms appearing in the reconstruction formula for a specific harmonic is tailored by an iterative, LS-based identification procedure, so that it reflects the nonlinear behavior of the IT at that harmonic. Computational efficiency and robustness are achieved by progressively updating the QR decomposition of the regressor matrix. The method has been experimentally validated on both a VT and CT. Results show that the number of terms required for a proper compensation of HD at different components may be very different. In this respect, the new technique can effectively cope with the wide measurement range of CTs and it allows reaching unprecedented performance in terms of tradeoff between accuracy and complexity, while avoiding overfitting. In fact, including non-significant terms in the compensation formulas may also degrade performance.
Future developments may be addressed by investigating (and possibly improving) accuracy under more complex operating conditions, for example, considering a significant offnominal frequency. Moreover, it would be interesting to scout the possibility of adopting the same compensation coefficients to different, nominally identical ITs.