Identiﬁcation of Grid Impedance by Broadband Signals in Power Systems with High Harmonics

: Grid impedance is an important parameter and is used to perform impedance-based stability analysis for the operation of grid-connected systems, such as power electronics-interfaced solar, wind and other distributed power generation systems. The identiﬁcation of grid impedance with the help of broadband signals is a popular method, but its robustness depends strongly on the harmonic disturbances caused by non-linear loads or power electronics. This paper provides an in-depth analysis of how harmonics affect the identiﬁcation of grid impedance while using broadband measurements. Furthermore, a compensation method is proposed to remove the disturbing inﬂuences of harmonics on broadband impedance identiﬁcation. This method is based on exploiting the properties of the used maximum-length binary sequence (MLBS). To explain the methodology of the proposed method, the design basis for the excitation signal is discussed in detail. The analysis from simulations and a real measurement in an industrial power grid shows the effectiveness of the proposed method in compensating the disturbing inﬂuences of harmonics on broadband impedance measurements.


Introduction
The increasing decentralization of the energy system through renewable energy sources and the electrification of all aspects of daily life are leading to more and more power electronics-based devices being connected to the power grid [1]. This increase in inverter-connected resources can lead to harmonic resonance, because of the interactions between the multiple-timescale of the converter control dynamics and grid impedance [2]. The harmonic resonance may first appear to be a power quality problem, but it is actually an indication of lack of system stability margin and may lead to instability and disruption of inverter operation if grid impedance or the share of inverter power further increases [3]. This is easily explained as small-signal stability depends on the source to load impedance ratio at any given interface, which must satisfy the Nyquist criterion, as has been demonstrated for both DC and AC systems [4,5]. Consequently, the identification of grid impedance is a critical point in the stability assessment and is often performed with broadband signals such as chirp signals [6], discrete-interval binary sequence [7] or maximum-length binary sequence [8]. The system response to the excitation signal can be used to calculate the non-parametric impedance, which represents the transfer relationship at discrete frequencies. Based on the non-parametric impedance the analytical parametric grid impedance of the system can be determined with the help of system identification algorithms [8]. However, these methods struggle in determining the impedance of power grids with a high share of harmonics due to non-linear loads as well as loads and sources connected to the grid via power electronic converters [9]. In [10] a method is proposed for compensating the disturbing influences of harmonics on the impedance measurement through a differential measurement by injecting harmonic currents with different phase angles. This approach, however, only works for mono-frequency excitation signals. Furthermore, no published work does an in-depth analysis of influence of harmonics on the identification of AC grid impedance using broadband signals and suggests how to solve these constraints.
This work attempts to fill this gap by proposing a method that compensates for the disturbing effects of harmonics on broadband measurements and enables system identification based on the non-parametric impedance. To understand why the identification of impedance is necessary, first the basics of impedance-based stability assessment is explained in Section 2 briefly. The standard procedure for an impedance measurement is then explained, starting with the design of the excitation signal. After the first step the excitation and measurement of the system, with subsequent calculation of the non-parametric impedance takes place. Finally, the system identification is performed. To compensate for the disturbing effects of harmonics, the properties of the excitation signal are exploited. Therefore, the selection, design and properties of the excitation signals are discussed in detail. Simulations are then used to show that the harmonic leakage effect is a major cause of the distorted impedance measurements. A method is proposed to compensate for these disturbing effects so that the system identification can be carried out based on the measurements. Based on simulations and a measurement in an industrial network with a high share of power electronics, it is shown in Section 3 that the proposed method can also be used in real scenarios to reduce the disturbing influence of harmonics on broadband impedance measurement, so that a reliable system identification is possible. The advantages, disadvantages and limitations of the proposed method are discussed in Section 4. In the end, the work is summarized in Section 5 and possible future fields of research are proposed.

Method for Compensation of Disturbing Influences on Impedance Identification
The impedance-based stability method was first used to study DC systems and was extended to 3-phase AC system [4]. The method is a type of small-signal stability analysis and can identify possible causes of instability problems. To apply impedance-based stability analysis, the system under investigation must be represented by a Thévenin and a Norton equivalent. In Figure 1, an inverter-grid system is analyzed, where the inverter is represented by a Norton equivalent and the power grid by a Thévenin equivalent. The current i g from the inverter into the grid can be calculated by (1), where Z g represents grid impedance and Z c the inverter impedance. v g is the grid voltage and G c is the transfer function from the current reference i * G c to the inverter output current i g , which must be stable when unloaded. The grid and inverter impedance forms the feedback path 1 + Z g Z c , where the ratio Z g Z c of grid impedance to the inverter output impedance must satisfy the Nyquist stability criterion in order for the inverter-grid system to be stable [4,11].

General Procedure for Impedance Measurement
To perform impedance-based stability analysis, grid impedance must be identified, which can be accomplished in different ways. Figure 2 tries to illustrate the following general procedure: Design of Excitation Signal: The excitation signal must be selected and designed very carefully, as it has a very large influence on the accuracy of the impedance measurement. With respect to the amplitude of the excitation signal, on one hand it must be small to ensure the system stays in close vicinity of the equilibrium point. On the other hand, the amplitude must be sufficiently large to reject noise disturbances [12]. As a general rule, the amplitude of the excitation signal is chosen between 5% and 10% of steady-state values [8].
Excitation and Measurement: To determine grid impedance, the system must be excited with 2 linearly independent signals, which can be realized in the phase sequence or in the DQ domain [13,14]. The excitation signal is then converted back into the ABC-natural reference frame where it is injected into the system by a series voltage excitation [15] or shunt current excitation [16,17]. Meanwhile the voltage response of the system during the shunt current excitation or the current response during the series voltage excitation are measured and recorded in the ABC-natural reference frame.
Calculation of Non-Parametric Impedance: The measurement of current and voltage is transferred to the desired domain where the non-parametric grid impedance is calculated using a Fast Fourier Transform [17].
Identification of Parametric Impedance: The parametric grid impedance is identified based on the non-parametric grid impedance. There are commercial tools available to achieve this in an easy way, such as the System Identification Toolbox of MATLAB [18], which is used in this work.

Design of Excitation Signal
Excitation and Measurement

Identification of parametric impedance
Calculation of non-parametric impedance MLBS Figure 2. General procedure for impedance measurement.

Design of an Excitation Sequence Based on the MLBS
There are different types of excitation signals to determine the response of a system. The simplest is a single sine wave at a frequency of interest. The single sin excitation provides the highest possible SNR and hence offers the highest accuracy in determining grid impedance [17]. Since each frequency of interest must be measured separately and an accurate measurement is only possible after the transient has decayed after each frequency injection step, this approach is very time consuming [19]. This drawback of long measurement times can be overcome using broadband signals that excite several frequencies at once.
Excitation signals can be divided into general-purpose signals, optimized test signals and advanced test signals. General-purpose signals do not require optimization and should be able to excite a system with a flat spectrum in a user-specified frequency band [20]. The periodic chirp signal [6] or the impulse [21] are examples for general-purpose excitation signals used for grid impedance measurement. Next an example for an optimized test signals is the optimized multisines excitation signal, which is a sum of user-specified sinusoids with varying amplitudes and phases [20]. By optimizing the phases of the frequency components, the time-domain amplitude can be reduced and the amplitude in the frequency domain maximized [22]. This property can be described by the crest factor (CF) which is defined as where T is the period, x peak is the peak value and x rms is the root mean square (rms) value of the signal [23]. Advanced test signals, such as the multilevel excitation signals, should be used only in critical conditions, where the special shape of the excitation gives a significant advantage [20]. The Maximum-Length Binary Sequence (MLBS) used in this work is a general-purpose signal from the Pseudo Random Binary Sequence group. It has the smallest possible CF making it suitable for sensitive systems [7,24]. It can be generated very easily by feedback shift registers as shown in Figure 3 and offers several degrees of freedom such as selecting the number of bits n of the shift register. MLBS is a popular choice because it can be easily designed using the parameters illustrated in Figure 4. The length N of the MLBS exist for N = 2 n − 1 and every bit of the MLBS is generated with a frequency of f clock = 1 T clock , which leads to a MLBS period duration of T MLBS = NT clock . The MLBS is repeated k j times, resulting in a period length for the excitation sequence of T seq = k j T MLBS . In practice, the values 0 and a are mapped to −a and +a to produce a symmetrical MLBS resulting in an average close to zero [20]. The generated excitation sequence, as shown in Figure 5a has the following properties [9]:   With MLBS different frequencies can be excited depending on the design of the MLBS. The frequency bins q n excited by the MLBS are illustrated in Figure 5b as blue dots and the characteristic frequency resolution ∆f exct between these bins can be calculated by From Equation (4) it can be derived that the resolution ∆f exct of the excited frequency bins depends only on the number of bits of the shift register n and the maximum frequency f EFB of the EFB, under the assumption of optimal choice of the clock frequency. By repeatedly exciting a system with the same MLBS, the frequency resolution ∆f seq can be achieved, which can be calculated by However, the repetition of the MLBS does not lead to further increase in excited frequency bins q n . It results in not excited frequency bins added by the repetition of the MLBS are shown in Figure 5b as blue circles p. The MLBS is repeated k j = 4 times to create the corresponding frequency spectrum. Each additional repetition of the MLBS results in another nonexcited frequency bin. According to (5), f EFB can be considered another degree of freedom in the design of ∆f seq . Unlike n, there is no logarithmic relationship between ∆f seq and the MLBS repetitions k j . This degree of freedom can be used to increase the frequency resolution by decreasing ∆f seq without losing amplitude of the effective frequency band Φ EFB . The power spectrum of a time series describes the distribution of power into frequency components composing that signal [25]. The power spectrum Φ MLBS of a MLBS can be described as (6).
with sinc(x) = sin(x)/x and q is the frequency bin under consideration. The MLBS is mapped from a to −a resulting in a peak-to peak amplitude of 2a [20]. Since the sinc function for q = 1 approaches 1, Equation (6) can be simplified to (7). As the slope of the effective frequency band can be considered flat, since it only drops by −3 dB up to f EFB , as illustrated in Figure 5a. Since the power spectrum of the EFB is flat up to f EFB , its amplitude Φ EFB can be calculated by (7), which equals the power of the first harmonic |Φ MLBS (q 1 )| of the MLBS.
From (7) it can be seen that the effective power spectrum Φ EFB of an MLBS excitation signal is only based on the amplitude a and the number of bits of the shift register n. It can be deduced that Φ EFB decreases by the ratio of (8) when increasing the size of the shift registers by one.
The only way to compensate for this decrease in the magnitude Φ EFB of the EFB is to increase the amplitude a of the MLBS.

Influence of Harmonics on Broadband AC Grid Impedance Identification
The increase in non-linear loads, as well as loads and sources connected to the grid via power electronic converters cause harmonics. In [9] it was shown that the harmonic currents emitted by these devices have a highly disturbing effect on the measurement of the response of the excitation signal and strongly corrupts the calculated impedance. One reason is the spectral leakage, which occurs when the periodic extension of a signal not commensurate with its natural period exhibits discontinuities at the boundaries of the observation [26]. This can be a problem for broadband excitation with an MLBS, since the harmonics in the power grid are not necessarily multiples of the frequency resolution ∆f seq . This fact can be better illustrated by the measurement window x wdw shown in Figure 6, which duration can be calculated by (9).
Since impedance is calculated based on the system response to the spectrum of the excitation signal, only the duration T wdw of the MLBS is measured using the measurement window x wdw . For this example, a shift register with n= 9 bits and a repetition of the excitation signal of k j = 1 is used. The EFB is intended to go up to f EFB = 9 kHz, which leads to a measurement duration of T wdw = 22.7 ms. According to (5) the frequency resolution of the MLBS is ∆f seq = 44.12 Hz, which is illustrated in Figure 6. In this figure, the harmonics are represented by the blue sine x 50Hz with a frequency of 50 Hz as an example.  Figure 6 shows that exactly one period of the excitation signal x exct is captured with the rectangular measuring window x wdw . The resulting spectrumx exct is shown in orange in the figure below. The signal x exct is completely mapped by the frequency bin at 44.12 Hz and there is no harmonic leakage to the frequency bins at 0 Hz and 88.24 Hz. In contrast, spectral leakage occurs in the power spectrumx harm , since no integer multiple of the harmonic signals x 50Hz is sampled through the window x wdw . This spectral leakage is reflected in the frequency components of the power spectrumx harm at 0 Hz and 88.24 Hz, which are not present in the time-domain signal x 50Hz .

Harmonic Leakage from Harmonics Due to Improper Window Size
The harmonic leakage effect due to different window sizes T wdw in the measurement of the grid voltage and the harmonics is investigated using the base scenario shown in Figure 7. The power grid is represented by a Thévenin equivalent with a voltage v g and an impedance Z g . The analytic expression for impedance in the dq-frame is a multi-input and multi-output system [27], which is derived in (10), with the inductance L g and the resistance R g of impedance as well as angular frequency ω g of the dq transformation.
The harmonic current caused by power electronic components and non-linear loads, is represented using a current source, which injects the current i harm into the grid. The harmonic current i harm is generated by a three-phase diode rectifier shown in Figure 8a and has the corresponding characteristic curve illustrated in Figure 8b. The injected harmonic current i harm has the typical spectrum ( Figure 9a) of a six-pulse rectifier, which was sampled with a proper measurement window of 20 ms, resulting in a frequency resolution of the spectrum of 50 Hz. Since the grid frequency f abc is 50 Hz, the fundamental and all higher order harmonics can be mapped by the FFT into the frequency domain without spectral leakage. According to [28] the corresponding power spectrum has the following characteristics.

System under Study Measurement
• The absence of third harmonics, which corresponds to zeros sequence components. • The presence of harmonics of orders 6k ± 1 for integer values of k.

•
The harmonics of orders 6k + 1 are of positive sequence and the harmonics of orders 6k − 1 are of negative sequence. Since there are no zero-sequence components, the three-phase quantities are modeled as equivalent two-phase quantities in the sequence domain by the symmetric component transform where I pos is the positive sequence current phasor and I neg the negative sequence current phasor at an arbitrary frequency. a = e j2π/3 corresponds to an 120°phase shift applied to the current phasors in the abc-domain. The sequence domain relates to dq domain according to Equation (12) [14].
Since the positive sequence components rotate in the same direction as the synchronous reference frame, the resulting frequency ω dq after the dq transformation is ω g less than the frequency ω pos of the positive sequence components (pos), with ω g as the angular frequency of the grid. The negative sequence components (neg) appear faster relative to the dq axes, because they rotate in the opposite direction as the synchronous reference frame. The resulting frequency ω dq is ω g higher than the frequency ω neg of the negative sequence components. As a result, the positive and negative sequence components in the spectrum of the harmonic current i harm in the abc-domain (Figure 9a) are superimposed on the spectrum in the dq-domain shown in Figure 9b. It becomes clear that the harmonic spectrum of a three-phase diode rectifier in the dq domain are limited to a multiple of 300 Hz. The amplitude of the harmonics decreases with increasing order. Figure 10a shows the effect of the mismatched size of measurement window caused by the duration T wdw of the excitation signal on the power spectrum of the harmonic current i harm in the frequency range from 100 Hz to 10 kHz. It is calculated according to (9) with the parameters Table 1. As defined in Figure 5b, the unfilled squares and diamonds indicate that the frequency bins are not excited by the MLBS. A measurement window of T wdw = 327.9 ms corresponding to a frequency resolution ∆f seq = 3.05 Hz used in this figure is the result of the MLBS designed with the parameters in Table 1. Table 2 shows the nearest multiples of the frequency resolution ∆f seq to the harmonics with their order of the multiple. As can be seen, the harmonics from Figure 9b are not integer multiples of the frequency resolution ∆f seq which leads to a leakage of the power spectrum of the harmonics into the neighboring frequencies. Neglecting the voltage v g , which is the fundamental component, the measured power spectrumv meas represents the voltage harmonics caused by the voltage drop due to the current harmonicsî harm at the grid impedance Z g . The corresponding relationship can be expressed byv meas = Z gîharm (13)

MLBS Impedance Measurement in Power Grid without System Harmonics
Based on Figure 11, the scenario of an ideal impedance measurement is investigated, in which the power grid is free of harmonics. The power grid is represented by a Thévenin equivalent with a voltage v g and an impedance Z g . The excitation current i exct is designed using the parameters in Table 1 and is equal to the measured current i meas . Since impedance Z g in the dq domain is a 2 × 2 matrix, two linear independent measurements are needed to determine the components [13]. Based on the two linear independent measurements M1 and M2 grid impedance can be calculated with the following equation: In the first measurement M1 the system is excited with the excitation current in the d-axis, whose power spectrum is shown in Figure 12a in the frequency range from 100 Hz to 10 kHz. As defined in Figure 5b, the unfilled circle and triangle indicate not excited frequency bins. Filled circle and triangle represents frequency bins which are excited by the MLBS. The measurement window of T wdw = 245.6 ms corresponding to a frequency resolution ∆f seq = 4.07 Hz used in this figure is the result of the MLBS designed with the parameters in Table 1. The power spectrum of the excited d-axisî d M1 has the characteristic flat curve of the MLBS. Since there is no excitation in the q-axis, the frequency component ofî q M1 is negligible. The voltage response of the system to the excitation in the d-axis can be seen in Figure 12b. The spectrum of the d componentv d M1 has a typical ohmic inductive curve as described by Z dd g in (10). The spectrum of thev q M1 describes the coupling between the d and q-axis indicated by Z dq g . Figure 13a shows the excitation with an MLBS excitation current in the q-axis and Figure 13b illustrates the corresponding voltage response of the system. Since the same MLBS is used, the same basic relationships apply as in Figure 12.
System under Study Excitation Measurement Figure 11. Ideal impedance measurement scenario in which the power grid is free of harmonics.  The relationship between the spectrum of the measured excitation currentî meas and the spectrum of the voltage responsev exct is given by (15), which shows that the voltage drop of the harmonics current at grid impedance is measured.

MLBS Impedance Measurement in Power Grid with Harmonics
The influence of harmonics in the power grid on the impedance measurement with MLBS is investigated for the realistic scenario as shown in Figure 14. The power grid is represented by a Thévenin equivalent with a voltage v g and an impedance Z g . The harmonic current i harm originates from a three-phase diode rectifier, which is represented as current source and has the characteristic curve shown in Figure 8b. The excitation current i exct is designed using the parameters in Table 1 and is equal to the measured current i meas . The measured excitation current i meas leads to the voltage response v exct of the system. As already derived in Equation (14), two linear independent measurements M1 and M2 are necessary to determine impedance Z g . In the first measurement M1 the system is excited with the excitation current in the d-axis. Since the power spectrum in Figure 15a is the same as in Figure 12a, the same observations apply. The spectrum of the voltage response in Figure 15b is the superposition of the harmonics, distorted by the leakage effect, shown in Figure 10b and the voltage response from the ideal scenario shown in Figure 12b. The voltage response can be expressed using the following equation:

System under Study Excitation Measurement
v meas =v harm +v exct = Z g (î harm +î meas ) The harmonics strongly dominate the measurement around their respective frequency, which are given in Table 2. One effect of this is that the frequency binsv d M1,Exct excited in the d-axis by i exct are difficult to separate from the nonexcited frequency binsv d M1,NotExct in the neighborhood of the harmonics. In contrast, the system provides a good voltage response in the excited d-axis at the frequency between the harmonics. The obvious cross coupling in Figure 15b between the excitation current in the d-axis and the voltage response in the q-axis is difficult to identify. It is superimposed by the harmonics, which is apparent in that the excited frequenciesv  A similar condition can be observed in the second measurement M2 when the q-axis is excited as shown in Figure 16a. Since the power spectrum is the same as in Figure 12a, the same observations apply. The spectrum of the voltage response in Figure 16b is the superposition of the harmonics, distorted by the leakage effect, in Figure 10b and the voltage response from the ideal scenario in Figure 13b. The harmonics do not distort the measurement M2 as much as the measurement M1, because the q component of the harmonic currentî harm is smaller than the d component. This results in a clearer voltage response of the system, which is represented by the strong deviation of the excited frequency binsv q M2,Exct compared to the nonexcited frequenciesv q M2,NotExct . The obvious cross coupling in Figure 13b between the excitation current in the q-axis and the voltage response in the d-axis is difficult to identify. It is overlaid by the harmonics, which can be seen by the fact that the excited frequenciesv d M2,Exct do not lift off from the neighboring nonexcited frequenciesv d M2,NotExct .

Distortion of Estimated Impedance
If Equation (14) is applied to the measurements M1 and M2, the dotted red nonparametric impedance Z g,NP in Figure 17 is obtained. Here the blue line is the reference bode plot of the analytical transfer function Z g,A according to (10). The green dots represent the calculated grid non-parametric impedance Z g,NP based on scenario in Figures 15 and 16 with harmonic distortions. The harmonic distortions of the three-phase diode rectifier at the multiples of 300 Hz can be seen in the magnitude and phase of the non-parametric impedance Z g,NP . Since the system has given a good voltage responsev d M1 to the excitation currentî d M1 in the d-axis and a good responsev q M2 to the excitation currentî q M2 in the q-axis, the disturbances in the diagonal elements Z dd g,NP and Z qq g,NP of impedance are not so strongly present and there is a good correspondence with the analytical functions Z dd g,A and Z qq g,A . In contrast, the harmonic distortions dominate the non-parametric coupling impedances Z dq g,NP and Z qd g,NP , which becomes apparent in the comparison with the analytical functions Z dq g,A and Z qd g,A . This is not surprising since the system gives a low voltage responsev in the q-axis.
Using the MATLAB System Identification Toolbox, the effects of harmonic disturbances during impedance measurements on a system identification algorithm is investigated. A dotted green parametric impedance Z g,P with two poles and one zero is fitted to every non-parametric impedance Z g,NP . To demonstrate the worst-case identification scenario the algorithm is used with default settings and without the use of a compensating filter. The identification algorithm can estimate impedances Z dd g and Z qq g well, as can be seen from the superposition of the analytical functions Z dd g,A and Z qq g,A with the estimated impedances Z dd g,P and Z qq g,P . In contrast, the identification algorithm is not able to estimate the coupling impedances Z dq g,P and Z qd g,P correctly. For impedance Z dq g,P the algorithm fits the impedance to the harmonic distortion peak at 300 Hz. The impedance Z qd g,P was fitted to the harmonic distortion peak at 2.7 kHz. This shows that the identification algorithm is not able to estimate the correct parametric impedance Z g,P based on the non-parametric impedance Z g,NP if the disturbance by harmonics is too strong.

Method for Compensating the Disturbing Effects of Harmonics on Broadband Measurement
The compensation method proposed in this paper is based on the evaluation of the system response and extends the measurement procedure as shown in Figure 18. Therefore, the system response is first evaluated in relation to the harmonics distorted by the spectral leakage effect. If the system response does not meet the requirements in certain frequency ranges, these are ignored. Afterwards the frequency components in the neighborhood of integers of the grid frequency are removed, because these can corrupt the results [29]. The last step in the proposed compensation method is the thinning of frequency response data. The thinning enforces equal weighting over the full frequency window of the nonparametric data [8]. Figure 18. Extended procedure for impedance measurement by proposed compensation method.

Response Quality Evaluation
The response evaluation is the first step in the preparation of the measurement for the calculation of the non-parametric impedance. It consists of several steps which are shown in Figure 19. In the first step, the system response in the frequency domain is calculated based on the measurements, through which we obtain the frequency spectrum q n with the blue frequency bins excited by the MLBS. Based on the repetitions k j of the MLBS, the frequency spectrum also consists of additional nonexcited frequency bins marked in red. For the proposed method, an MLBS repetition of k j ≥ 2 is necessary, since in step two the quality of the system response Ψ n is evaluated based on the nonexcited frequency bins. The Quality assessment involves determining the strength of the system response relative to an estimated baseline measurement Θ n according to following equation.
The estimated baseline measurement Θ n for an excited frequency bin q n is calculated based on the κ neighboring nonexcited frequency bins p n . Since no excited frequency bins q n are to be included in the estimation, κ ≤ kj − 1 must be satisfied. The estimated quality of the system response Ψ n is calculated according to equation.
For this purpose, the ratio of the magnitude of the frequency bin q n and the estimated base measurement Θ n is calculated. If the quality of the system response Ψ n does not meet the requirements Ψ tresh , these frequency bins Γ n are set to zero and thus not taken into account in the calculation of the non-parametric impedance. The threshold Ψ tresh must be set individually for each system based on the intensity of the harmonics in the power system. The threshold must be set individually for each system based on the strength of the harmonics in the power grid. In general, for a case with a high percentage of harmonics, it is recommended to go up from a small value, as the Ψ tresh can be very sensitive when removing frequencies.

Harmonic Neighborhood Cancelation
The frequency components in the neighborhood of integers multiples of the grid frequency can corrupt the non-parametric impedance [29]. These distortions are not reliably removed by the compensation step, which can be seen at frequency bin q 8 in Figure 19 step 3). The gray curve represents the trajectory of harmonics that have been distorted by the spectral leakage effect. The system response q 8 at this frequency is directly on the gray trajectory, which indicates that it is dominated by harmonics. Since the neighboring nonexcited frequency bins p 7 and p 9 are far below q 8 , the estimated frequency bin Θ 8 also turns out to be low. This results in a good quality of the system response Ψ 8 , which means that this frequency bin q 8 is not set to zero, even though the system response is very poor in reality. To avoid this problem, the frequency components in the neighborhood of integers of the grid frequency in the range of f HNC are completely removed.

Data Thinning
A thinning imposes uniform weighting over the entire frequency range of the nonparametric data. The data thinning procedure is used to obtain a logarithmic spaced subset of the data points. For this purpose, on the one hand, the lower and upper frequency limits must be specified and, on the other hand, the number σ of data points for the data thinning routine must be defined. All data points located outside this frequency range will not be considered for further processing [8]. Figure 19. Procedure for determining the quality of the system response.

Simulative Validation of Compensation Method
Applying the proposed compensation method to the scenario Figure 14, we obtain the compensated measurement M1 shown in Figure 20 and the compensated measurement M2 shown in Figure 21 can be obtained. These were determined under the conditions in Table 3. Comparing the uncompensated case in Figure 15b, where the d-axis gets excited, with the compensated case Figure 20b, it is noticeable that a clear curve of the voltage responsev d M1,Exct in the d-axis can be seen. It follows the progression described by Z dd g in (10). The outliers at the harmonics of order 6k with k as integer value were reliably removed by the compensation method. This also applies to the coupling response in the q-axis which is described by Z dq g . Comparing the uncompensated response of M2 in Figure 16b with the compensated response Figure 21b, the dominant d component of the harmonic distortion current i harm becomes apparent. The system response in the q-axis is clear, as can be seen from the well-defined voltage responsev q M2,Exct , which describes the progression of Z qq g in (10). In contrast the coupling described by Z qd g in (10) is too low, which leads to the fact that only a few frequency binsv d M2,Exct fulfill the requirements Ψ tresh . The total measured system responsev d M2,Exct is dominated by the d component of the current i harm , which is apparent in the uncompensated case shown in Figure 16b. To weight the individual frequency components equally during system identification, 30 frequency bins were selected, which are logarithmically distributed across the frequency range.   Based on these measurements, the compensated non-parametric impedance Z g,CNP shown in Figure 22 can be calculated according to (14). In comparison with the uncompensated non-parametric impedance Z g,NP shown in Figure 17, the much lower proportion of outliers in the coupling terms is noticeable.

Identification of Industrial Power Grid Impedance
To demonstrate the effectiveness of the proposed method, a system identification is performed in an industrial grid with lots of power electronics. The measurement setup is shown in Figure 23. The injection of the excitation current into the power grid is realized by an Spitzenberger & Spies APS15000 linear amplifier. This linear amplifier operates in current-controlled mode and feeds the reference excitation current i exct into the grid. The amplitude of the reference excitation current is set to i exct = 10A. The measurement of the injected current i meas and the voltage response v meas of the system as well as the generation of the reference excitation current i exct is performed by an Opal OP5600 real time simulator. The linear amplifier allows an almost ideal injection of the excitation current and the Opal OP5600 ensures synchronized measurement of i exct and v exct . To align the excitation signal to the d-axis and the q-axis to obtain the two required linearly independent measurements, the phase angle of the power grid is needed. This phase angle is provided by a DSOGI-FLL implemented on the Opal OP5600, which provides a fast dynamic response with low overshooting [30].  The MLBS used to excite the industrial power grid was designed according to the parameters in Table 1. A shift register with a length of n = 11 bits is used. The effective frequency band of the MLBS is f EFB = 10 kHz, which leads to an excited frequency resolution ∆f exct = 12.21 Hz. Since the MLBS is repeated k j = 4 times, the resulting measurement window is T wdw = 327.9 ms corresponding to a frequency resolution of ∆f seq = 3.05 Hz.
The frequency response M1 of the industrial grid after excitation of the d-axis is shown in Figure 24. The frequency spectrum of the current in the d-axisî d M1 depicts the characteristic curve of the MLBS shown in Figure 24a. In contrast to the simulation in Figure 15, there is a small coupling betweenî d M1,Exct andî q M1,Exct in this case. The voltage responsev d M1,Exct of the system illustrated in Figure 24b is good for the excited d-axis. The coupling response in the q-axis, on the other hand, is dominated by harmonic disturbances. The excited frequency binsv q M1,Exct do not stand out from the nonexcited frequency binŝ v q M1,NotExct , which makes compensation necessary that is shown in Figure 25.  The compensation of the measurement was performed with the parameters in Table 4. The compensated current frequency spectrum Figure 25a shows the logarithmic thinning on 150 bins and the cancelation in the neighborhood of multiples of the grid frequency. Since the current response is good and meets the requirements for the quality of the system response Ψ tresh , no frequency bins are compensated. In contrast, a large part of the frequency bins of the voltage response Figure 25b was compensated by the compensation method. The high proportion of harmonic disturbances at 300 Hz leads to the complete compensation of the system responsev d M1,Exct in the d-axis up to 400 Hz. From then on, the requirements for Ψ tresh are met forv d M1,Exct . The coupling responsev q M1,Exct is compensated to a large extent due to the high proportion of harmonics.v q M1,Exct meets the requirements for the quality of the system response Ψ tresh only in the range from 500 Hz to 2 kHz. The same applies to measurement M2 with the excitation of the q-axis shown in Figure 26. The frequency spectrum of the current in the q-axisî q M2 depicts the characteristic curve of the MLBS shown in Figure 26a. In contrast to the simulation in Figure 15, there is a small coupling betweenî q M2,Exct andî d M2,Exct in this case. The voltage responsev q M2,Exct of the system illustrated in Figure 16b is good for the excited q-axis. The coupling response in the d-axis, on the other hand, is dominated by harmonic disturbances. The excited frequency binsv d M2,Exct do not stand out from the nonexcited frequency binsv d M2,NotExct , which makes compensation necessary.
The compensated current frequency spectrum Figure 27a shows the logarithmic thinning on 150 bins and the cancelation in the neighborhood of multiples of the grid frequency. Since the current response is good and meets the requirements for the quality of the system response Ψ tresh , no frequency bins are compensated. In contrast, a large part of the frequency bins of the voltage response Figure 27b was compensated by the compensation method. The high proportion of harmonic disturbances at 300 Hz leads to the complete compensation of the system responsev q M2,Exct in the q-axis up to 400 Hz. From then on, the requirements for Ψ tresh are met forv d M1,Exct . The coupling responsev d

M2,Exct
is compensated to a large extent due to the high proportion of harmonics.v d M2,Exct meets the requirements for the quality of the system response Ψ tresh only in a narrow frequency range from 1 kHz to 2 kHz.  Based on Figures 24 and 26 measurements, the blue non-parametric impedance Z g,NP shown in Figure 28 can be calculated. As in the analytical analysis, harmonic disturbances are very dominant at the frequencies 300 Hz and 600 Hz, which leads to the characteristic peaks in the non-parametric impedance Z dd g,NP due to the spectral leakage effect. The coupling between the d and the q-axis is also very weak, which leads to the strong scattering of the non-parametric impedance in the coupling terms Z dq g,NP and Z qd g,NP . After applying the proposed compensation method, the compensated red non-parametric impedance Z g,CNP can be calculated from the compensated measurements Figures 25 and 27. The non-parametric impedances Z dd g,CNP and Z qq g,CNP have a clear curve, which allows a very good fit of the green parametric impedances Z dd g,CP and Z qq g,CP . Despite the weak coupling and harmonic disturbances, a good fit of the parametric impedance Z qd g,CP was possible. However, this does not apply to the coupling impedance Z dq g,CP where only a limited frequency bins meet the system response requirements Ψ tresh .

Discussion
The proposed compensation method is reliably able to remove frequency components that do not meet the system response requirements. In Section 3 it was possible to compensate the harmonic disturbances on the broadband measurement in a simulative scenario and a real measurement in an industrial power grid. Based on these measurements, it was possible to determine the transfer function of the system under investigation. The determination of the coupling impedances Z dq g and Z qd g proved to be difficult. Due to the weak coupling in the systems studied, the quality of the system response in the q-axis during an excitation of the d-axis and the response in the d-axis during an excitation of the q-axis was not sufficient. This leads to a cancelation of the large frequency components in the non-parametric coupling impedances, which makes an accurate system identification difficult or even impossible. If this is the case, it is reasonable to neglect impedances Z dq g and Z qd g due to the low coupling. The system under consideration thus becomes two decoupled SISO systems, which simplifies subsequent stability considerations.

Conclusions
In this paper, a method is presented to compensate harmonic disturbances during the identification of grid impedances using broadband excitation signals. The design of MLBS broadband excitation signal has been addressed in depth, as its inherent characteristics are the basis for the developed compensation method. Based on simulations, it is shown that the leakage effect due to the harmonics caused by power electronic components and an inappropriate window size in the FFT transformation are a major cause of the disturbances. The characteristics of the MLBS excitation signal to excite only certain frequency components are used to estimate the quality of the system response. The characteristic of the MLBS to excite only certain frequency components is used to estimate the quality of the system response, allowing frequency components that do not meet the requirements to be compensated. Based on simulations and a real measurement in an industrial power grid, it was shown that harmonic disturbances in broadband impedance measurements can be compensated with the proposed method. Using the compensated system response, the diagonal elements of the power grid impedances could be reliably determined with the help of a system identification algorithm. With the compensation method, the determination of the coupling impedances has proven to be challenging. Due to the weak coupling, the system response is dominated by the harmonics in the corresponding axes. If the proposed method is used, large frequency components are canceled, which makes system identification problematic.
Future work will investigate the dynamic adjustment of the frequency resolution of the MLBS. This will be achieved by measuring the grid frequency with a PLL and adjusting the generation of the MLBS accordingly. This should reduce the leakage effect due to an inappropriate window size of the FFT and the harmonics caused by power electronics.