A novel iterative observer approach for real-time harmonic estimation in power distribution networks

In power distribution networks, the realisation of online harmonic state estimation puts forward a higher demand on the algorithm accuracy, stability as well as simplicity. To meet the requirements above, a novel iterative observer based method for harmonic state estimation is proposed. This method guarantees the stability of the harmonic estimation observer when output time delays are considered in the power distribution network. Meanwhile, signals with different frequencies can be estimated iteratively. System harmonics are regarded as the unknown exogenous disturbance, and each harmonic at multiple measurements is extracted and estimated iteratively. To guarantee the overall estimator stability, multiple transmission time delay is addressed in the system modelling for the design of observer gain. The proposed method achieves high estimation accuracy from multiple measurements while maintaining a low computation commitment. Comparing with traditional harmonic state estimation methods, the iterative observer is capable of estimating the unexpected harmonics. The IEEE 13-bus power distribution system built in the real-time digital simulator demonstrates the validity of the proposed method.


INTRODUCTION
Due to the fact that harmonics can cause communication errors, overheating, hardware damage etc., the harmonic problem has received a great deal of attention in power system to minimise the possible economic losses. Harmonic estimation is important to effectively identify and help to eliminate the harmonic sources in the system and thereby guarantees high levels of power quality. Many techniques have been proposed to estimate the harmonics generated by power electronic loads. These schemes usually use the phasor measurement and envelope tracking to extract harmonics. The most preferred method is the fast Fourier transform (FFT) [1][2][3][4] due to its simplicity. However, the FFT method has the limitation of aliasing, spectra leakage and picket fence effect caused by non-synchronous sampling [5][6][7][8]. To overcome these issues, parametric methods are proposed, such as singular value decomposition (SVD) [9], estimation of signal parameters via rotational invariance technique [10], multiple signal classification [11], Kalman filter (KF) [12], autoregressive model [13] and prony method  [14]. However, these methods are restricted the use in offline systems.
On the other hand, the development of real-time simulation platforms has prompted researchers to pay more attention on online harmonics estimation. In order to coordinate with the online implementation, algorithms for harmonic estimation should be more applicable while maintaining their high stability and high accuracy in the power distribution network. Some contributions have been made in the area of online harmonic research. In [15], phasor measurement units (PMUs) is utilised to obtain global positioning system (GPS) synchronised measurements in the harmonic monitoring system. But this approach is at the cost of high computational burden, which is not suitable for large-scale power distribution network. In [16], a lab-scale waveform generator with the realtime digital simulator (RTDS) and power amplifiers for testing the harmonic disturbances from physical equipment in a closed-loop environment is presented. However, this method cannot ensure estimation accuracy of the results due to the arbitrary selection of mother wavelet [17]. In [18], the proposed PMU-based HSE is able to determine harmonic states in non-fully observable power distribution networks with consideration of uncertainties over measurement, loads and harmonics. However, errors on large lateral feeders are increased, which indicates that optimisation problems with inequality constraints need to be considered. In [19], a method for assessment of harmonic levels along medium voltage distribution systems is presented. This method quickly assesses the harmonic performance of the network with low errors. However, it cannot be applied to assess harmonic voltages at the resonance frequencies. In [20], a time domain method for HSE is proposed in unbalanced three-phase power systems including non-linear loads. The optimal number of involved devices and the property of half-wave symmetry in voltage and current waveforms are considered in this method to reduce the computational burden.
As an alternative choice, observer-based HSE has proven its effectiveness and has been successfully applied in [21][22][23][24][25][26]. In [21], an approach based on synthesis of non-linear observer is proposed. However, bad input constitutes the singularity of the observer problem. In another method [22], SVD and pseudo-inverse are applied into HSE without the need of observability analysis while its computation progress is complex and has high computational burden. In [23], a composite observer based HSE is presented, which facilitates the application of dq transformation for a single-phase system. However, it has large fluctuation and cannot ensure the stability in the power distribution network. In [24], a method of output feedback control is presented for state estimation and disturbances estimation in electro-hydraulic servo actuators. This method reduces the online computational burden and guarantee the system stability in the unmeasured system. In [25], an adaptive output feedback robust control scheme is proposed for hydraulic servo systems, which handles parameter uncertainties, matched and mismatched disturbances simultaneously. Among these, the iterative observer avoids unnecessary computation on the non-existent harmonics by calculating the total harmonic distortion (THD) of residuals before each iteration, and thereby guarantees low computational burden [26].
In terms of long-distance transmission characteristics, how to deal with transmission time delay is significant for implementing real-time monitoring [27]. In [28], transmission time delays are discussed in power hardware in the loop simulations using discrete Fourier transform. In fact, time delay is one of the major sources causing inaccuracy in the dynamic model, which may intensively affect the system stability and degrade the system performance [29][30][31]. Hence, to test and validate the effects of transmission delays in real-time harmonic estimation, multiple time delays are non-negligible issues that have to be taken into account at the estimator design.
This paper proposes an iterative observer-based algorithm for HSE in the real-time environment, which is stable, accurate and does not need any parameter tuning to respond to harmonic changes. Compared with the traditional harmonic estimation methods, the proposed algorithm goes beyond the framework of fixed frequencies, which can estimate harmonics with iteratively frequencies. Additionally, the transmission time delays are taken into consideration in the observer design. Further, the distribution system is developed in RTDS to make it feasible towards online harmonic estimation in the power distribution network. The harmonics with different harmonic orders are extracted, and then their amplitudes and phases are obtained. The value of the THD is controlled below 5% according to IEEE Std-519. When THD of the output error is computed below 5%, the iteration will terminate, which indicates that only tolerable harmonics are left in the residual measurements.
The contributions of this paper are summarised as follows: (i) The iterative observer dynamically estimates multiple signals with various frequencies simultaneously, rather than estimating the fixed-frequency harmonic. (ii) Multiple transmission time delays are modelled for observer gain design. Furthermore, the dimension of the observer gain keeps low and constant irrespective of the increasing system dimension, which is well suited for the large-scale complex distribution system with multiple measurements. (iii) Evaluations on the accuracy, stability, and dynamic harmonic estimation of the proposed iterative observer are carried out online in the real-time environment using the RTDS.
In the remainder of the paper, the design of harmonic estimator is provided in Section 2, which include system and harmonic modelling, iterative observer derivation, and observer gain selection considering time delay. Section 3 is dedicated to describing the hardware and software configurations of realtime digital simulator, which implements real-time data acquisition and data transfer systems. In Section 4, a real-time case study of IEEE 13-node system is presented to demonstrate the flexibility and accuracy of the proposed approach, followed by conclusions in Section 5.

HARMONIC ESTIMATION USING AN ITERATIVE OBSERVER
In order to estimate and analyse the system harmonics in the real-time environment, a novel iterative observer based HSE is proposed. To enable the observer to estimate system states in unsynchronised measurements, the proposed approach is designed in time domain. Both continuous-time observer and discrete-time observer are considered.
As shown in Figure 1, the non-linear load, producing harmonics in the power distribution network, is designed as a linear load in parallel with several current sources at specified harmonic frequencies [32]. The direction of the current sources points to the voltage source of the system. Additionally, effects of renewable energy on harmonic disturbances in the power distribution network are taken into account, such as wind turbine.

Continuous-time observer
Consider a linear power system with harmonic disturbances consisting of multiple output delays. The system is given bẏ where u(t ) is the system input, x(t ) is the system state variable which consists of branch currents and node voltages in the power distribution system [33] and y(t ) is the measurement which is the combination of the branch current (node voltage) as well as the harmonics at the testing node. w(t ) represents the harmonic disturbance with the disturbance matrix E. E is unknown during the harmonic estimation. The matrices A, B, C j are constant in a given power distribution network, j ∈ [1, m] and j is treated as the output delay, m is the number of measurements. The delay profile j is defined as We have When time delay is ignored, w(t ) is a combination of harmonic content at each harmonic order, which can be written as where = 2 fh, f represents the fundamental frequency, v is the number of harmonics, h is the harmonic order and i is the phase angle at ith order. Since E cannot be directly measured, the harmonic disturbance can be regarded as an exogenous system, which is given byẇ where is the output of the exogenous system. Hence, with consideration of time delay, the harmonic disturbance can be expressed aṡw The distribution system dynamic model and exogenous harmonic disturbance dynamic model are augmented to yield: where Denoting and C ′ = [ C I 0 ], the Luenberger observer is constructed as .
The state variable can be written aŝx wherex is the actual value of system states and V iŵi is the exogenous states. V i is designed as a weighting matrix. LetS = I m+1 ⊗ S, g = I m+1 ⊗ g, A = I m+1 ⊗ A,B = I m+1 ⊗ B, and C = C I . Taking derivatives of Equation (9), we geṫ̂x Comparing (11) and (13), we have V i in Equation (15) can be seen as z in the Sylvester equation Fz + zG = H, with z calculated as Thus, the estimate at ith harmonic order iŝw Equation (17) can be set as the initial harmonic. In order to implement additional harmonic estimation, we define a set of integers P = {h i } for i = 1, … , n. Thus the estimate of additional harmonic order h i+1 can be written aṡ̂w Hence, the state estimate of the system can be expressed aŝx Taking derivatives with time, we havė̂x =Āx +Bū + (L Comparing Equations (11) and (21), the observer gain can be expressed as Consequently, when i = 1, and when i = 2, … , n, The expansion form is expressed as (27). The observer error dynamics are determined as Taking derivatives of Equation (27), we geṫ whereS h,i andḡ are, respectively, expressed as (30) and (31). Therefore, L h,i is designed such that (S h,i − L h,iḡ ) is Hurwitz, where the error approaches zero asymptotically and thus the stability of the observer can be guaranteed. In the proposed method, the only variable of the observer gain L h,i is the harmonic frequency and the observer gain remains 2 × 1 irrespective of the size of the power distribution network, thereby maintaining the accuracy.
For multiple measurements points, the state-space matrices remain unchanged, while the observer gain L h,i is designed to accommodate the number of measurements m and the number of harmonic orders i. Hence, L h,i ∈ R 2p×p , and p is the number of measurement points.

Discrete-time observer
The output signal in the exo-system is discretised to (t ) = (kT ) by the analog-to-digital converter for kT ≤ t < (k + 1)T, k = 0, 1, … , ∞ [34]. The exo-system can be modelled as where The sampling interval T is made implicit in the following progress. Consider a discrete-time plant model: The algorithm of the discrete-time iterative observer can be expressed as ALGORITHM 1 To determine TE based on VPSI and S SC Obtain node voltage or branch current measurements, y.
Design the discrete-time observer without consideration of harmonic disturbances: x(k + 1) =Āx(k) +Bū(k) +L(y(k) −Cx(k)) If THD of (y(k) −Cx(k)) > Threshold then Define a set of positive integers P = {h i } for i = 1, … , n Implementation of the multiple point approach is similar to the continuous-time case. Equations (30) and (31) can be directly applied to the discrete-time method.

Data acquisition and data transfer system specifications
Real-time implementation of developed harmonic estimation in a power distribution network is realised using the RTDS and the NI DAQ 6343 card. The simulation model data for real-time data acquisition is produced by the simulation of a power distribution network in RSCAD. RSCAD is the power system simulation software of the RTDS technology. The NI DAQ 6343 card is used for data acquisition in the working environment of Lab-VIEW. Figure 2 shows the progress of data acquisition between the RTDS and PC.

Configuration of data acquisition card in LabVIEW
To open the link between the RTDS and MATLAB, the key component NI DAQ 6353 card is used for processing the multiple output signals from the RTDS. The LabVIEW graphical program is shown in Figure 3. To acquire the real-time data from the RTDS, the physical channel of NI DAQ is configured based on the power distribution network in RSCAD, which means the voltage (current) ranges must cover the maximum and minimum value of the analogue data. The acquired data is retrieved in the while loop and saved as arrays.

LAB EXPERIMENT RESULTS
An unbalanced three-phase 13-node distribution network [38] is used as the test system. The one-line diagram of the IEEE 13-node system is shown in Figure 4. The end load at node  34 is replaced by a wind turbine. Minimum wind speed is limited to 5 m/s and the maximum is 25 m/s. The nominal speed of the wind turbine is temporarily set as 15 m/s. In practical, increased voltage of a wind turbine can cause burnout of components. Hence, the output voltage of the wind turbine is set as constant. The IEEE 13-node system is simulated in RSCAD, which is the power system simulation software by the RTDS  Technologies. Two balanced non-linear loads are added to the system at nodes 46 and 75. The spectra of harmonics are shown in Tables 1 and 2.
The matrices A, B, and C are determined by the algorithm presented in [35,36]. The optimum meter placement we apply to choose the measurement points is based on binary integer linear programming [37], which is able to ensure complete observability for power distribution networks.
For each phase, there are 12 measurements taken in the IEEE 13 node feeder, which include I 31−33 , I 33−34 , I 45−46 , I 32−71 , I 84−911 , I 92−75 , V 31 , V 33 , V 45 , V 92 , V 71 , and V 84 . I a−b represents the branch current measurement at branch a − b; V a represents the node voltage measurement at node a. All measurements are considered in the iterative observer. By calculation, the measurements, whose THD above the threshold, are I 45−46 , V 45 , I 92−75 , V 92 , and V 84 . The iterative observer is applied to these distorted measurements until the THD in iteration is lower than the threshold. Data acquisition card transfer the measurement data to MATLAB for harmonic states estimation. Figure 5 shows the measurement of branch current I 92−75 produced in RSCAD and collected in MATLAB. It can be  Table 3 shows the calculated THD of I 92−75 during harmonic estimation. At the beginning, all the value of THD are above the threshold, which is set as 0.5% according to Std-519. Apparently, the THD becomes smaller than 0.5% after estimating the seventh harmonic at each phase, thus the iteration goes to the termination. The iterative observer is carried out to all distorted measurements in the discrete-time form and Figure 6 presents the harmonic state estimation in branch current I 92−75 . This measurement point is picked due to its proximity to the harmonic source at node 75, which indicates that it has the most serious harmonic distortion.
To test the stability problem in transmission time delays, three different time delays are added to the three-phase output of the system dynamics, which are, respectively, 0.03, 0.04, and 0.01 s. Figure 7 shows the three-phase delayed measurements of branch current I 92−75 . The waveforms in Figure 8 shows the estimated harmonic contents in I 92−75 with the time delay  at each phase. Apparently, harmonics in the unbalanced threephase system are estimated with short period delays and quickly become stable.
In an actual power distribution system with wind turbines, the wind speed cannot keep constant due to natural laws. To test the effect of wind speed to harmonic estimation, the wind turbine is increased to 24 m/s in the power distribution network with same time delays on each phase. Figure 9 shows the measurements of I 92−75 with 24 m/s wind turbine and Figure 10 shows the corresponding harmonic estimation results. The amplitudes of the waveforms at each phase are increased when compared with the results of 15 m/s wind turbine. Hence, under a constant voltage output, the increased wind speed can raise the currents as well as the harmonic injections in non-linear loads in a power distribution network. Table 4 shows the comparison of actual and estimate values of all distorted measurement around node 75 in the system. The errors of all measurements are far below the threshold, which verifies the high accuracy of the iterative observer.
To further investigate the benefits of the iterative observer in the power distribution network, a comparison test between the proposed method and KF is established at node 46. A 0.01-s transmission time delay is considered when measuring the branch current. At 0.03 s, an extra fifth harmonic appears at node 46, acting as an unexpected harmonic, which is set as 14.8% of the fundamental and 0 phase angle. Figures 11 and 12 show a comparison between the iterative observer estimates and the KF estimates for the I 45−46 measurement. After a short transient, both KF and the iterative observer become stable within 0.02 s and identify with the actual value. Figure 13 shows the harmonic estimation results of the extra harmonic. The iterative observer estimates the harmonic after the time delay while KF fails. Due to the lack of setting in KF algorithm, the estimate value of KF approaches to zero approximately after a transient. Different with KF, the proposed method estimates harmonics with iteratively frequencies, which not only includes the fixed frequencies but also the unexpected ones.

CONCLUSION
This paper proposes a real-time iterative observer based harmonic estimation method. All the distorted measurements can be carried out simultaneously in the power distribution network. The presented algorithm is fast and flexible, which can measure both harmonic currents and harmonic voltages with the consideration of time delays. Case studies consider the effects of non-linear load and the unexpected  harmonics on the harmonic estimation. Simulation results show that the proposed approach is capable of estimating the harmonics in an unbalanced three-phase power distribution network with high accuracy. Besides, if the mutual influ-

FIGURE 11
Harmonic estimation for the fundamental frequency of I 45−46 measurement ence of harmonic sources among nodes is considered, the proposed method should be augmented with other harmonic source identification techniques for a better real-life application.

FIGURE 12
Harmonic estimation for the third harmonic of I 45−46 measurement