Time-Domain Aggregation of Interharmonics from Parallel Operation of Multiple Sustainable Sources and Electric Vehicles

Abstract

This paper examines the random nature of interharmonics generated by power converters connected to sustainable energy sources and loads, such as wind turbines, photovoltaic (PV) panels, and electric vehicles (EVs). Current research often overlooks the stochastic behavior of interharmonics and their impact on power system reliability and resilience, leading to gaps in effective modeling and mitigation strategies. Thus, this study examines a low-voltage installation with a PV panel, an EV and a microwave operating simultaneously, providing practical insights into real-world scenarios of interharmonic related disruptions and solutions for enhancing the reliability and resilience of sustainable energy grids. By leveraging real-time measurements of interharmonics, suitable probability distribution functions (PDFs) are initialized to develop a probabilistic model using Monte Carlo simulation. This enables the derivation of a time-domain aggregation model of interharmonics from multiple sources operating together at the point of common coupling (PCC). The findings reveal that the peak values of voltage or current fluctuations at the PCC are influenced by the randomness in the number of devices connected and the frequency components originating from different sources. Through multiple case studies, the dependency of these fluctuations on stochastic parameters is systematically established. Empirical relationships are formulated to predict aggregated interharmonic values under varying scenarios, enhancing the accuracy and applicability of the model. The results demonstrate that higher interharmonic frequencies and fewer randomly connected devices significantly increase the probability of elevated aggregated peak values. These insights can serve as benchmarks for grid operators and policymakers in mitigating interharmonic related issues in modern power systems.

1. Introduction

Interharmonics are frequencies that are non-integer multiples of a power system frequency and can exist as discrete frequency components or as wideband spectrums [1]. Estimation of interharmonics is important and has gained significant attention recently, due to their impacts on power systems, like power transformer saturation, interference with data acquisition, unwanted tripping of relay protection circuits, instabilities in PLL-based converter control systems even with very low amplitudes, light flickering and monitor image fluctuations, aging and thermal effects, sub-synchronous oscillations of systems, ability to excite dominant resonances, interferences with harmonic compensation and harmonic filtering mechanisms, erroneous firing of thyristors [2,3,4,5,6], etc. The occurrence of extremely low-frequency sub-synchronous interharmonic voltages, amounting to less than 1% of rated voltage in the supply, has led to an increase in the peak value of the no-load current of a single phase transformer to about nine times the normal no-load value, causing hysteresis loop asymmetry and saturation [6].

Besides the classical interharmonic sources like cycloconverters, static frequency converters, AC/DC electrical drives, arc furnaces, induction furnaces and pulsating loads asynchronous with fundamental frequency, etc., inverter-based sources like PV (photovoltaic) and wind farms, as well as HVDC links, are also sources of interharmonics and their penetration in power systems is rapidly increasing [2]. It is shown in [7] that even for interharmonic values of only 5% amplitude, the modulation amplitude reaches as high as 10% (up to the seventh harmonic), thus significantly impacting the peak amplitude under very low interharmonic levels. The interharmonics’ aggregation in grid current leads to cycle-by-cycle variations in peak values, resulting in uneven charging and discharging of DC capacitors in rectifier loads. These fluctuations, if severe, can cause equipment malfunctions and shorten the lifespan of end-user devices [7]. It is therefore important to monitor the instantaneous peak value variations in the time-domain, with increased interharmonic emissions in the grid due to the increasing proliferation of power electronic devices. Differences between harmonics and interharmonics are in the aggregation of the emission from different sources. The emission of interharmonics varies strongly between devices and is without a common reference (e.g., for harmonics, it is normally the zero crossing of fundamental voltage) so that it will not be likely that many devices emit asynchronised interharmonics at about the same frequency. The interharmonic distortion subgroup, established within the framework of the IEEE PES Transmission & Distribution Committee’s Harmonics Working Group (519) [2], concluded that while modern devices currently in use may not face significant issues related to interharmonic interference, this does not eliminate the need to anticipate potential interharmonic problems in the future. It emphasizes the importance of verifying equipment compatibility with interharmonics and implementing stricter limits to address these challenges proactively.

One of the research gaps identified with respect to interharmonics for establishing limits [2] is the lack of understanding of how interharmonics aggregate under parallel operation of multiple interharmonic sources, to which this work contributes. This finds application and relevance in current power systems where parallel operations of multiple inverters in a PV farm, as well as in a wind farm, and multiple EVs’ charging (which are sources of interharmonics) at the point of common connection exist. The aggregation of interharmonics from multiple sources has not been extensively analyzed. Such analyses could be of interest to standardization committees in establishing standards. Hence, the data obtained from multiple experimentations and real-time measurements, along with the Stochastic Model, mainly contribute to determining the probability of obtaining peak values of interharmonic aggregated current from multiple parallel sources. Thus, in the future, long-term measurements will not be required for interharmonic aggregation analysis. By analyzing short-term measurements and estimating the dominant interharmonic frequencies, the aggregated peak value can be predicted using the developed formula or deterministic model and compared against the established standards.

In this article, the division of sections is structured to present a clear and logical flow of the research. Section 2 outlines the Proposed Methodology, detailing the approach for investigating the time-domain aggregation of interharmonics. Section 3 delves into the Characteristics of Interharmonics, focusing on observations from field measurements in photovoltaic (PV) systems, wind parks, and electric vehicle (EV) charging. Section 4 presents the Stochastic Model Results, exploring the likelihood of interharmonic aggregation. This is followed by Section 5, which provides Verification with Real Data, validating the findings against empirical evidence. Finally, Section 6 concludes this paper, summarizing key insights and implications.

2. Proposed Methodology

The proposed methodology shown in Figure 1 for investigating the time-domain aggregation of interharmonics will address the following aspects:

For identification of interharmonic sources and impacts in real-world scenarios, a review of the existing literature and field data to map the origins and impacts of interharmonics in systems like PV systems, wind turbines, LED lamps, and electric vehicles (EVs) was performed. Further, interharmonics observed in both controlled laboratory settings and field measurements were analyzed to determine interharmonic characteristics such as peak magnitude value of aggregated interharmonics, phase angle, etc. Various signal-processing techniques such as Discrete Fourier Transform (DFT), Short-Time Fourier Transform (STFT), Modified Sliding Window ESPRIT (SWESPRIT), Modified Desynchronized Processing (DP) technique, High-Accuracy Frequency Interpolation Technique were used in MATLAB 2020 to assess the realistic values of interharmonics in measurements of both steady and time-varying waveform distortions in a real scenario. The graphical tool spectrogram was used for visualization.

Deterministic and Stochastic Models were used to simulate interharmonics under various scenarios. Further, a Monte Carlo simulation was used to examine the aggregation of interharmonics from multiple sources under various scenarios at the point of common coupling (PCC). Mathematical expressions for two cases of aggregated interharmonic current peak, (a) randomness in the number of interharmonic frequencies aggregated at PCC and (b) constant number of interharmonic frequencies at PCC, were determined and further verified using real-time data.

This methodology will provide a systematic framework for understanding, estimating, and setting standards for interharmonics in modern power systems. The prediction model developed will also help to understand the worst-case scenario due to the aggregation of interharmonics and to check whether it will exceed the limits in standards, thereby addressing critical gaps in current knowledge and practices.

3. Characteristics of Interharmonics Observed in Field Measurements from PV, Wind Park, and EV Charging

A simple technique implemented widely in PV inverters for active islanding detection is grid impedance estimation through the injection of interharmonic frequencies into the grid. Under unfavorable circumstances of asynchronous operation of multiple inverters with such control in a PV farm, cases have been reported of phase angle cancellation and dilution of the interharmonic peak value, leading to a higher grid impedance estimation and false tripping of protection circuits [8,9,10]. On the contrary, it can also lead to phase angle aggregation of the interharmonics and an increased peak value can lead to a lower grid impedance estimation, causing a failure to trip, in addition to the upstream grid being disconnected and thus forming islands, which is a serious issue. In the worst case, the aggregated interharmonic peak value can be the number of PV inverters times the interharmonics from a single inverter, as shown in [8,9,10]. As reported in [11,12], cascaded H-bridge (CHB) converters, commonly used in PV applications due to their modularity and high efficiency, can generate interharmonics that are more pronounced than those produced by individual PV inverters. This occurs because of the potential superposition of MPPT perturbations on the DC voltages of individual CHB cells. In the worst-case scenario, the interharmonics can be up to n times larger than those from a single PV inverter (n representing the number of cascaded cells) if the DC-side oscillations of the CHB cells are in phase. Consequently, interharmonic issues in CHB PV inverters are significantly more severe and complex compared to two-level PV inverters.

Figure 2a represents the interharmonic current (after filtering out fundamental and harmonics using a desynchronized processing technique [13]) from a single PV inverter and emulated interharmonic current from 10 such PV inverters, all with the same active islanding technique of grid impedance estimation using periodical interharmonic injections of 75 Hz into the grid when the interharmonic currents are in phase. It is visible that the peak value of the aggregated interharmonics has increased 10 times and has a symmetrical variation.

Figure 2b presents the emulated aggregated interharmonic current at the output of 11 PV inverters connected in parallel, where 10 inverters were emulated as before. Measurements from another inverter with MPPT control that introduced low-frequency interharmonics of less than 100 Hz due to MPPT perturbations [14] were considered as the 11th inverter. It is visible from Figure 2 that the peak value of the aggregated interharmonics has a wide range of asymmetrical variation depending on phase angle cancellation or aggregation.

Figure 3 presents a spectrogram of grid current (fundamental removed) for a day at the output of a 20 kW three-phase rooftop PV installation where two 10 kW PV inverters operate in parallel with fixed axis tracking. Figure 4 and Figure 5 represent the same measurements zoomed in between 5:30 a.m. and 08:30 a.m. and 18.42 p.m. and 22:00 p.m., respectively for easy visibility.

The observed interharmonics are discrete components and appear early in the morning during sunrise and during the evening hours of sunset when the power production is smaller than 30%, as shown in Figure 3. The reasons explained in [14,15] are the global horizontal irradiance dominated by the direct horizontal irradiance at those instances and the reactive power control activated by the converters operating in parallel. The interharmonic frequencies are subharmonics (10 Hz, 16.8 Hz, 25 Hz) and their multiples, which are time-varying. These interharmonics co-exist at the same time at many instants; for example, the marked points A, B and C in Figure 4 and Figure 5 lead to aggregation and peak value variations in the time-domain waveform. From an extensive analysis of measurements, it was identified that three interharmonic frequencies with the highest amplitude (can be any multiple of these subharmonics) exist pre-dominantly and can be contributed by each of the two paralleling operating inverters [15]. It was verified that the random behavior of these interharmonic frequencies contributed to the peak value variations in the time-domain aggregated interharmonic current.

Figure 6a represents an example of interharmonic current waveform of the same installation where multiple frequencies aggregate and the peak value of interharmonics varies with respect to time. In Figure 6b, the 3 days’ data of a photovoltaic system with its instantaneous individual cycle peak values of interharmonics are shown. The increase in the peak value of interharmonic current is visible during the morning hours of the day and the evening hours of the day for all three days.

Another example of issues caused by aggregation and amplification of interharmonics in wind farms and HVDC transmission systems is reported in [16,17,18]. The sub-synchronous interactions (SSIs) leading to instability issues are caused by interharmonics generated by the operation of multiple doubly fed induction generators, as discussed in [15]. If the inter-harmonics in the grid side contain a frequency component that coincides with the system’s sub-synchronous oscillation frequency, resonance may occur and, in turn, provoke a more severe SSI. If the amplitude of the aggregated inter-harmonics frequency components exceeds 1%, the rotor shaft torsional vibration is excited, which would cause great damage to the rotor shaft system. The other is the frequency of inter-harmonics component close to the frequency of generator rotor shaft, which would not excite serious torsional vibration that could cause damage to the rotor shaft in a short time.

An example of measurements performed in a 14-turbine wind park with each turbine rated at 2500 kW is illustrated in Figure 7. It represents a 2D spectrum of 4-day measurements plotted one over the other at the turbine terminal for easy identification of the time variations in the amplitude and frequency range of interharmonics. The fundamental component and harmonics are removed and residues of filtered away harmonics are visible. From the spectral analysis of measurements corresponding to 4 days, it was identified that interharmonic frequency appears as a time-varying narrow band spectrum close to the 13th harmonic of which 635 Hz was dominant, i.e., with the highest amplitude. The same frequency is observed at each of the turbine terminals and becomes aggregated at the collection grid.

To understand the phase angle characteristics of the interharmonics for model initializations, the complex plot of the observed interharmonics was plotted and an example is shown in Figure 7. From the complex plot, the random phase angles of interharmonics compared to harmonics were inferred, as also stated in [19].

Figure 8 presents another example of a 2D spectrum of 60 measurements at the output terminal of a single-phase 6 kW AC EV charger connected in a network, as shown in Figure 1 in [20], plotted one over the other. The observed interharmonics (residues of harmonics are also visible) were time-varying (both amplitude and frequencies) in the narrow band range from 1000 to 1400 Hz.

The effect of inrush transient interharmonics caused by the operation of a vacuum cleaner together with a PV inverter, causing false tripping of the PV inverter (with implemented active islanding injecting interharmonics into the grid), under different scenarios is experimented with and explored in [21]. The transient interharmonics lasting for less than a second can momentarily increase or decrease the peak value of the aggregated interharmonic current. Figure 9 represents an example of a measurement carried out at the PCC of a low-voltage installation with a series of events, as PV is in continuous operation, the microwave starts operating at 4.5 s and the EV starts to charge at 6.5 s. The impact of the PV operation is distinctly seen in the spectrogram every 1 s. The increase in a broadband of frequencies (transient interharmonics) due to the start of the microwave’s operation and EV charging is also visible.

From the experience with interharmonics, it can be concluded that interharmonics have a random behavior and are time-varying, appearing either as discrete components or as a narrow band or as broadband/variable frequencies [22] (transient nature) with random phase angle characteristics [23]. The amplitude of an interharmonic in most cases varies randomly between 0 and 1 A. With this background, a Stochastic Model was developed to understand the probability distribution of the aggregated interharmonics and bring forth an approximate mathematical expression for predicting the aggregated peak value of interharmonics in Section 4.

In addition to very few studies on harmonic aggregation [24,25] and supraharmonic aggregation in a black box model [26] being reported in the literature, no studies on interharmonic aggregation are available. It was mentioned in [2] that most of the current inferences about interharmonics are based on limited data and experiences and little consideration has been given to important factors like the severity of effects and likelihood of occurrence. To address these issues, probabilistic studies offer more insight (a range of possibilities) to evaluate the extremities due to interharmonics and their aggregation effects. This paper presents a generic Stochastic Model of the interharmonic aggregation from multiple sources. The statistics of interharmonic aggregation are derived using Monte Carlo simulations and empirical relations are formulated.

4. Stochastic Model Results for Likelihood of Interharmonic Aggregation

The summation of the individual interharmonic current waveform at the point of common coupling (PCC) from multiple devices can be represented as follows (1) and (2):

$$I_{ih}={\overline{I}}_{1ih}+{\overline{I}}_{2ih}+{\overline{I}}_{3ih}+\dots +{\overline{I}}_{nih}$$

(1)

$$I_{ih}={\sum }_{k=0}^n{\sum }_{l=0}^{nf}{A}_l\cos (2\pi f_{l}t+{\varnothing }_l)$$

(2)

$A_l$ is the amplitude, ${\varnothing }_l$ is the phase angle corresponding to the interharmonic frequency $f_l$, $nf$ is the number of interharmonic frequencies from one device and $n$ is the total number of devices.

The field measurements indicate that interharmonic current amplitude and phase angles have a random nature and are difficult to correlate with any system parameters. Hence, random variables are used to represent the different parameters, such as amplitude, phase angle and the uncertainty of the number of devices with interharmonics.

$$P\left[I_{ih}\right]=max\{{\sum }_{k=0}^n{{\sum }_{l=0}^{nf}\left(R\right(A}_{k1})\cos (2\pi f_{l1}t+{R(\varnothing }_{k1})))+{R(A}_{k2})(\cos (2\pi f_{l2}t+{R(\varnothing }_{k2})))+\cdots \}$$

(3)

$P\left[I_{ih}\right]$ is the probability of the time-domain aggregated interharmonic peak value, ${R(A}_k$) is the random amplitude chosen from a normal distribution within the interval [(0, 1)], and ${R(\varnothing }_k$) is the random phase angle chosen from a uniform distribution within the interval [(−π, π)]. n corresponds to the random number of sources of interharmonics connected to PCC at the same instant, and $f_l$ is a random interharmonic frequency from a uniform distribution within the interval [(1–2000 Hz)].

A Monte Carlo simulation is utilized to analyze the probability distribution of the time-domain aggregated interharmonic peak value over 10 s windows, based on the parameters specified in model (3): n = a random number between 1 and 100 (for each iteration), $nf$ = 4, simulation cases = 1000, 10,000, 20,000. CP95, which is the 95th percentile value [27], is considered for evaluating the aggregated peak in each iteration of a 10 sec window. It is important to highlight that the total interharmonic frequencies aggregating at the PCC are critical in the time-domain. These frequencies exhibit randomness, either due to the unpredictable connection or disconnection of multiple interharmonic sources (random n) or the variability in multiple frequencies originating from a single source (random $nf$).

Multiple simulations were carried out in order to evaluate the evolution of the aggregated interharmonic peak value. A left skewed Type I extreme value distribution (Gumbel distribution) [27] with the maximum likelihood estimates (MLEs) of [25.5346 6.0861] and confidence intervals of [25.5346 6.0861; 25.9323 6.3949] was obtained, as can be observed in Figure 10. The skewness is defined by the limits of perfect phase angle cancellation (zero peak) and perfect phase angle aggregation (max. peak) of the aggregating interharmonic frequencies at the PCC.

A Type I extreme value distribution (Gumbel distribution) is defined as follows:

$$f\left(x\right)={\displaystyle \frac{1}{\sigma }}{e}^{-{\displaystyle \frac{x-\mu }{\sigma }}}{e}^{{-e}^{-{\displaystyle \frac{x-\mu }{\sigma }}}}$$

(4)

The cumulative distribution function, CDF is defined as follows:

$$F\left(x\right)=\exp \left[{-e}^{-{\displaystyle \frac{x-\mu }{\sigma }}}\right]$$

(5)

μ is the location parameter (corresponds to the 95th percentile peak value with the highest probability, Iagg_peak95), and σ the scale parameter (corresponds to the population density of the samples, tabulated using interquartile percentile, iqr).

It was confirmed that the left-skewed distribution results from the variability in the number of interharmonic frequencies aggregating at the PCC, whereas if it is not random, the left-skewed extreme value distribution in all the cases converged to a normal distribution (see model case study 2 in the forthcoming section).

From Figure 10, it is clear that the highest probability of occurrence of 95th percentile aggregated peak value lies between 20 and 25 A.

Some of the significant observations (refer to

Table 1. Comparison of parameters.

Cases	μ (A)	iqr (%)	$\mathit{P}[{\mathit{I}}_{\mathit{i}\mathit{h}}$] (%)
1000	19.99	3.45	7.26
10,000	21.02	4.07	7.37
20,000	20.93	4.14	7.73

) are:

(1)

μ is independent and remains almost constant for all the simulation cases (1000, 10,000, 20,000).

(2)

The interquartile (iqr) range, the difference between the 25th and 75th percentiles in the data—that is, the width of the region that contains the middle 50% of the data—increased slightly with the simulation cases.

(3)

The probability of the occurrence of the 95th percentile peak value of time-domain aggregated interharmonic remained almost constant with the simulation cases.

Two simulation model case studies with different outcomes are detailed in Section 4.1 and Section 4.2. Due to the above inferences from 1 to 3, the simulation cases in the model case studies discussed further are limited to 1000.

4.1. Model Case Study 1—Randomness in the Number of Interharmonic Frequencies Aggregated at PCC

The model parameters considered are as follows: $n$ = (1–10), (1–20), (1–30), …, (1–100), which is selected randomly within the defined limits, $nf$ = 2, 4, 8, and Simulation case = 1000.

Some of the significant observations from Figure 11a,b and

Table 2. Comparison of parameters, case study 1.

nf = 2
n	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
1–10	4.12	17.32	25.84
1–50	9.63	8.48	14.57
1–100	13.31	5.17	10.21
nf = 4
n	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
1–10	6.02	13	21.1
1–50	13.91	5.91	10.07
1–100	20.8	3.81	7.33
nf = 8
n	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
1–10	8.55	9.29	16.1
1–50	19.41	3.81	7.51
1–100	27.9	2.5	4.96

are:

(1)

As the randomness in the number of interharmonic frequencies aggregated at PCC increases, the skewness toward the left increases (wider range for zero peak to max. peak) but the likelihood of achieving an aggregated 95th percentile peak value decreases exponentially.

(2)

μ increases with the increase in the number of interharmonic frequencies.

(3)

The interquartile range, the distance between the 25th and 75th percentile in the data—that is, the width of the region that contains the middle 50% of the data values—decreases with the number of interharmonic frequencies.

(4)

The probability of the occurrence of the 95th percentile peak value of time-domain aggregated interharmonic decreases with the increase in the number of interharmonic frequencies.

The location parameter μ represents the time-domain aggregated 95th percentile peak value of interharmonic current with the highest probability of occurrence (Iagg_peak95). An empirical relationship between μ and the total number of interharmonic frequencies randomly aggregating at the PCC is approximated in Equation (6). This relationship is established through extensive Monte Carlo simulations and curve fitting using MATLAB tools:

$$Iag{g}_{peak95}~=\sqrt{N}$$

(6)

Here, n represents the total number of randomly aggregated interharmonic frequencies at the PCC. This is also inferred from the same values tabulated under different conditions in rows 5 and 9 in Table 2. A goodness-of-fit test of (6) is performed and tabulated in

Table 3. Comparison of goodness of fit of (6) for different simulated cases, case study 1.

nf	R2	Linear Model Poly1: y = p1 × x + p2 Coefficients (with 95% Confidence Bounds):
2	0.9924	p1 = 1.07 (0.9931, 1.146) p2 = −1.042 (−1.844, −0.2412)
4	0.9960	p1 = 1.041 (0.9878, 1.095) p2 = −0.715 (−1.508, 0.07856)
8	0.9918	p1 = 0.9876 (0.9143, 1.061) p2 = 0.01888 (−1.519, 1.557)

and presented in Figure 12. The analysis, combining Monte Carlo simulations and goodness-of-fit tests, suggests that the peak value of time-domain aggregated interharmonic current with the highest likelihood of occurrence can be effectively modeled by the value (6).

4.2. Model Case Study 2—Constant Number of Interharmonic Frequencies at PCC

The model parameters considered are as follows: n = 10, 20, 30, …, 100 (constant), nf = 2, 4, 8, Simulation case = 1000.

Some of the significant observations from Figure 13a,b and

Table 4. Comparison of parameters, case study 2.

nf = 2
nc	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
10	5.127	31.11	49.67
50	11.87	32.2	45.87
100	16.37	31.14	45.98
nf = 4
nc	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
10	7.2	29.69	48.46
50	16.45	33.71	48
100	23.13	32.51	45.98
nf = 8
nc	$Iag{g}_{peak95}$ (A)	iqr (%)	$P[I_{ih}$] (%)
10	10.12	34.38	44.56
50	23.45	31.09	48.52
100	32.89	30.31	46.47

are:

(1)

As the constant number of interharmonic frequencies becomes aggregated at the PCC, the skewness toward the left disappears and the likelihood of obtaining an aggregated 95th percentile peak value remains almost constant with the increase in the number of interharmonic frequencies.

(2)

μ increases with the increase in the number of interharmonic frequencies.

(3)

The interquartile range tends to remain almost a constant independently of the number of interharmonic frequencies.

(4)

The probability of the occurrence of the 95th percentile peak value of time-domain aggregated interharmonic tends to remain almost constant independently of the number of interharmonic frequencies.

An empirical relationship between μ and the constant number of aggregating interharmonic frequencies can be approximated as in (7).

$$Iag{g}_{peak95}~=1.2\times \sqrt{N}$$

(7)

n denotes the total constant number of interharmonic frequencies observed at the PCC.

A goodness-of-fit test of (7) is performed, tabulated in

Table 5. Comparison of goodness of fit of (7) for different simulated cases, case study 2.

nf	R2	Linear Model Poly1: y = p1 × x + p2 Coefficients (with 95% Confidence Bounds):
2	0.9995	p1 = 1.007 (0.9881, 1.026) p2 = −0.1965 (−0.4312, 0.03814)
4	0.9995	p1 = 0.9875 (0.9698, 1.005) p2 = 0.03835 (−0.2692, 0.3459)
8	0.9999	p1 = 0.9996 (0.9916, 1.008) p2 = −0.1526 (−0.3499, 0.04469)

and represented in Figure 14. Monte Carlo simulations and goodness-of-fit analysis indicate that the peak value of time-domain aggregated interharmonic current with the highest probability of occurrence can be approximated as (7). This value differs from (6) by 20% more, i.e., a factor of 1.2. This 1.2 factor in (7) is attributed to the constant number of interharmonic frequencies that are aggregated at the point of common coupling (PCC).

The difference with randomness (model case study 1) and without randomness (model case study 2) can be observed with the distribution of 95th percentile peak value plotted w.r.t simulation cases, as shown in Figure 15.

From Figure 15, it is observed that the aggregated peak value varies between 0 and 25 A in the first case, and it varies between 20 and 25 A in the second case.

A simple logic developed with an illustrative example to identify the frequency of occurrence of peak value in case of random interharmonic aggregation:

Using the lowest interharmonic frequency as the reference, 100 samples were collected per half-cycle of its positive peak. With four interharmonic frequencies aggregating randomly, the probability of observing a maximum value within these 100 samples is estimated to be 1/100⁴. Consequently, the estimated time (T) for 100% probability of encountering the maximum aggregated peak can be determined as per Equation (8). For the lowest frequency of 1 Hz, the expected occurrence of the maximum aggregated peak with 100% probability is once every 1.58 years.

$$T={\displaystyle \frac{\left(\mathrm{T}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{y}\right)}{\begin{array}{c}(\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\mathrm{o}\mathrm{f}\mathrm{g}\mathrm{e}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{f}\\ \mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{y})\end{array}}}$$

(8)

The preceding relationship suggests that the maximum peak occurs more frequently when a greater number of lower-frequency interharmonics are aggregated. Conversely, the likelihood of observing the maximum peak value diminishes with fewer aggregating interharmonics, particularly those residing within higher frequency ranges. This is also in line with Figure 11, where, because the number of randomly aggregating interharmonic frequencies is low, the probability of occurrence of a 95th percentile aggregated peak value is high.

5. Verification with Real Data

To verify the mathematical expressions derived in (6) and (7), the PDFs are plotted for the interharmonic measurements in PV installation shown in Figure 4, Figure 5 and Figure 6 and wind park shown in Figure 7.

5.1. Data Set 1

A left-skewed probability distribution indicates a min aggregated peak and a max. aggregated peak are obtained for the aggregated interharmonic of the PV installation with two PV inverters operating in parallel and injecting three dominant interharmonic frequencies in the grid current in a random way. The data are collected for a defined time every day after sunrise for a long-term period of 6 months. According to (6), the following expression is obtained:

$$Iag{g}_{peak95}~=\sqrt{2\times 3}=2.44\mathrm{A}.$$

(9)

From Figure 16, it is observed that the 95th percentile aggregated interharmonic current value with the highest probability of occurrence is 2.33 A, which is close to the estimated value using (6) with a deviation of 0.1 A.

Figure 17 represents the cycle-by-cycle peak value variations in interharmonic current in the time-domain of a phase of the PV inverter for a period of 4 months extracted directly from the installed power-quality meter (Figure 6b extended). It can be inferred that for 95% of the observed period, the aggregated interharmonic peak value is less than 2.4 A, which agrees with the mathematical expression in (6).

It was verified that even though subharmonics and multiples of subharmonics were observed in the spectrogram of current measurements, it is the three dominant interharmonic frequencies (with the highest magnitude) that contributed to the peak value.

5.2. Data Set 2

A normal distribution is obtained for the data measured at the collection grid of the wind park with 14 turbines, where each turbine injects a single interharmonic frequency of 635 Hz into the grid, as shown in Figure 7. The collected data are of 10 min resolution for 4 days. According to (7), if a constant number of interharmonic frequencies aggregate at the PCC, the following expression can be obtained:

$$Iag{g}_{peak95}~=1.2\times \sqrt{14\times 1}=4.5\mathrm{A}$$

From Figure 18a, it is observed that 95% of the time the percentile aggregated interharmonic current value with the highest probability of occurrence is 4.48 A, which closely agrees with the estimated value using (7) with a deviation of 0.02 A. From Figure 19, it can be observed that for 95% of the observed cases, the peak value of 635 Hz is less than 4.48 A. It was verified that, for some measurements, the value of 635 Hz is occasionally amplified and approaches values greater than 10 A due to the closeness to the parallel resonance frequency of the network estimated as the 13th harmonic.

Figure 18b represents the peak value of the aggregated interharmonic current for 4 days, where it is seen that, in 95% of the cases, the peak value is less than the estimated value of 4.5 A.

From the two analyzed data sets, it could be inferred that the long-term interharmonic data from PV installation follow model case study 1 and the short-term data from wind park follow case study 2. Thus, it can be concluded that in the analyzed data, there exists randomness in the number of interharmonic frequencies aggregating for the PV installation at the PCC, whereas for the wind park, it is in the constant number of interharmonic frequencies aggregating at the collection grid.

6. Conclusions

This paper presents approximate mathematical expressions to estimate the aggregated 95th percentile peak value of time-domain interharmonic current with the highest probability of occurrence. It was systematically established that the randomness and the constancy in the number of dominant frequency components that aggregate at the PCC is what matters for the aggregated interharmonic peak value in time-domain. This is besides the general random characteristics of interharmonics in terms of phase angles, amplitudes and frequency ranges, and is the main difference between the two model case studies. This is evident in the skewness observed in model case study 1, which arises from the inherent randomness, and the absence of skewness in model case study 2. Furthermore, it was logically determined that the probability of observing a linear summation of peak values increases within shorter time periods when aggregating interharmonics with higher frequencies and fewer components.

The findings confirm that the possibility of encountering a worst-case scenario of high aggregated interharmonic peak values at the PCC cannot be discounted, especially in systems where multiple sources, such as PV farms, wind farms, EV charging stations, and low-voltage installations, operate concurrently. From the two case studies, it is evident that the number of interharmonic frequencies, whether random or constant, significantly impacts the aggregated peak values at the PCC. A key benchmark quantitative result is that under constant conditions, the aggregated peak value is 20% higher than under random conditions, as shown by the factor of 1.2 between the models. This highlights the importance of considering the variability in interharmonic frequency aggregation when designing and operating sustainable energy systems.

The developed time-domain aggregation model provides a scientific foundation for developing interharmonic-specific standards. Current standards like IEEE 519 [28] or IEC 61000 [29] do not comprehensively address interharmonic aggregation limits, and this model can help committees establish new guidelines tailored to the growing penetration of sustainable energy and EVs.

While harmonic limits are well established in standards like IEEE 519, interharmonic limits remain underexplored. The model allows standardization committees to:

• Understand the interharmonic aggregation patterns at PCC due to multiple sources and defining acceptable aggregation limits for interharmonic levels, thus recommending dynamic or time-domain thresholds instead of fixed steady-state limits.
• Develop test methodologies for interharmonic assessment in labs.
• Establish compliance criteria for equipment manufacturers to ensure that devices like EV chargers and inverters operate within permissible interharmonic levels.

This ensures that interharmonics from multiple sources do not exceed grid thresholds. Hence, this model and the methodology can support a standardization committee to set new aggregation limits and revisions to standards to accommodate modern grid challenges.