Circular regression models of modern harmonic producing sources

: This study presents a novel measurement-based method for modelling harmonic current injections of low-voltage power electronics equipment. The aim of the developed methodology is to accurately represent harmonic current phase angle in response to the applied harmonic background voltage at the terminals of modelled device. In this method, the harmonic angles are analysed with the techniques of circular statistics. The performed study proves the adequacy of representing phase angles as sets of directional data rather than linear. After graphical and quantitative analysis, the regression models with underlying von Mises distribution are derived in the form of algebraic equation. This algebraic equation features sine and cosine terms and explains relationship between independent variable – voltage harmonic angle and dependent variable – current harmonic angle. The subsequent analysis of the residuals allows to conclude that models developed by applying this technique can be used in commercial simulation packages and allow to account on any value of background harmonic voltage angle.


Introduction
Voltage and current harmonic distortions are taking nowadays place among power quality issues with rising amount of complaints by customers. Increasing number of devices are being interconnected with the grid through power conditioning equipment. Recent detailed studies confirm that harmonics causing such well-known issues as overloading of power system components with tremendous thermal impact accompanied by poor power factor and additional power loss require multilateral research [1,2].
Accurate models of harmonic producing sources are essential for studying their impact on voltage distortion levels. The challenge in developing a model suitable for harmonic propagation studies in the distribution networks is to accurately represent the impact of background voltage distortion on harmonic current emission of certain device. This challenge is governed by a tradeoff between the complexity of a model and its level of details. Moreover, not only the modelling approach itself is important but also the convenience of applicability of developed model in the commercially available simulation packages.
The current modelling practice encompasses several approaches but among the most efficient is Norton-coupled frequency admittance matrices (FCM). This modelling method accounts on the sensitivity of the harmonic currents to the angle of background voltage distortion. Additionally, it allows to consider harmonic interaction between voltages and currents of different harmonic orders. A thorough attention to this method was given in [3] where coupled capacitor-smoothed bridge rectifiers were modelled as coupled admittances. The authors of this study proved that harmonic current is influenced significantly by applied harmonic voltage amplitudes and phase angles. Moreover, complex deviations of model output vector from the measurements were used to evaluate accuracy of the investigated models. The results showed 10-20% average deviations for 50% quantiles of data.
In [4], the authors emphasised the importance of metrological requirements to the measurement system used for obtaining initial data for Norton models. They demonstrated that derivation of FCM depends strongly on the accuracy of such a system.
The properties of admittance elements with respect to the topologies of power electronics equipment (PE) have been studied in [5,6]. The Fourier descriptors were proposed for estimating and analysing FCM. It was found out that depending on the topology of power electronic device, specifically on the type of power factor correction (PFC) some equipment types can indicate wide linear range in the coupled Norton approach -a brilliant observation leading to the prediction of harmonic behaviour of certain devices. Based on the measured values, the authors of these papers proposed the quantification of linearity -a non-linearity harmonic current response index.
To conclude, the Norton frequency-coupled matrices is a powerful modelling approach, which de facto requires significant amount of accurate measurement sets and ample mathematical derivations. Its implementation in commercial simulation packages can be somewhat limited since an iterative harmonic load flow would require a large number of admittance matrices comprised in advance in order to take into consideration all scenarios of interest.
As an alternative to Norton FCM, a component-based modelling of harmonic loads is capable of reproducing the instantaneous input current waveform of each individual power electronic device. Such models require the fundamental knowledge about topology of studied devices and numerical values of the components. In [7], a method for modelling residential aggregate loads for harmonic analysis was presented. Generic models based on the equivalent circuits of the power electronics-driven equipment have been developed. Based on the results, authors of this paper concluded that proposed aggregate models were suitable for preserving full information on the electrical characteristics of the modelled load mix. This approach was further developed in [8]. The main contribution of this work is a combination of detailed component-based modelling approach with the Markov-chain Monte Carlo routine allowing to account on user-load interaction. The case study involved modelling of groups of compact fluorescent lighting (CFLs) and important observation concerned diversity factors to quantify harmonic cancellation effects between different groups.
A rigorous component-based modelling of modern switchedmode power supply, photovoltaic (PV) inverter and electric vehicle charger has been shown in [9]. In this work, authors also introduced an influence of changing power operating mode on resulting harmonic emission levels as well as gave an important remark to test modern PE in presence of non-sinusoidal supply voltage waveform. This modelling approach in combination with probabilistic methods was also applied in [10] for studies of the commercial load sector.
As it can be deduced from these works, the component-based models introduce high level of details and associated complexity into the simulation scenario but they, nevertheless, require knowledge on generic underlying equipment topologies, type of the control system and values of the base components. Being able to accurately reproduce harmonic content of a certain device or a group of individual devices, they are bound to be used within specialised simulation software, for instance EMTP.
Finally, deterministic modelling methods are complemented by harmonic fingerprints. This modelling process was firstly described in [11] where harmonic current injections were comprised in lookup table according to the applied individual voltage background harmonics. The phase and magnitude of each individual harmonic voltage were changed stepwise providing at the end of procedure full harmonic current spectrum of equipment under study. The measurement approach is quite similar to the one used for FCM models, yet the post-processing part does not require complex mathematical equations. The idea was further developed in [12] for aggregate models of PV inverters.
Next, stochastic aggregate harmonic load models form a distinct family of methods focusing on application in stochastic harmonic propagation studies. Among research works involving development of such probabilistic models are [13,14]. These methods are best suited for simulating an impact of large bulk of aggregate lowvoltage load on medium-voltage networks with both composition and participation factor of loads usually modelled as random variables. As it was shown in these works, the accuracy of such models relies on the extent of measurement campaign typically lasting weeks at one particular location. On the other hand, the uncertainty of probabilistic models can be large if part of the measured data is substituted with historical or typical values.
As it can be seen, a large variety of modelling methods is currently available. Nevertheless, harmonic models described by certain deterministic function derived from the reasonable amount of field measurements are still widely sought. Evidently, development of such models suggest specific data-fitting procedures. A strong advantage of the algebraic equation representing behaviour of the harmonic source is its ease of application in simulation environment and reckoning on a large amount of simulation scenarios. Some work has been done in order to apply data-fitting routines to the measured sets of data. In [15], authors presented a method based on a linear least-square (LLS) estimator performed on the decoupled Norton admittances matrices of CFL. Due to the numerical instability of the algorithm it was not possible to obtain either accurate estimates of full Norton model or coupled admittance model. The results of this particular case study demonstrated reasonably good accuracy of the derived model in terms of current magnitudes but omitted angle derivation. The underlying reason in disregarding angle of harmonic is the nature of LLS, which requires modifications in order to resolve phase information.
In [16], an excellent method for modelling background harmonic voltage source has been presented. A non-linear leastsquare (NLSF) was applied to the sets of field data after statiscal processing which involved smoothing and filtering the data with Gaussian distribution. While the output of this algorithm in the form of normal distribution deemed to be appropriate for representing current magnitudes of combined operation of many harmonic sources, a different approach is required to model a response of the certain harmonic load.
To conclude, despite the research aiming to derive algebraic equations associated with harmonic sources, there is currently a lack of harmonic models of that type.
In this paper, we propose a novel method for modelling individual harmonic producing sources based on the foundations of circular statistics. A measurement procedure attributed to the harmonic fingerprinting [11] is used for obtaining the initial sets of data for three modern PE devices: PV inverter, battery charger and a group of CFLs. After that, coefficients of the models are estimated with the application of circular regression algorithm. The output of the modelling process is compared with the measured data to check for satisfactory performance.
The novelty of the proposed method lies in the procedure of deriving the relationships between applied harmonic background voltage angle and the phase angle response of the harmonic current emission of the device. The described modelling process brings an improvement over the fingerpint method [11]. Further, on contrary to another high-level details models, the proposed procedure is characterised by reduced level of complexity, yet providing stable output with regard to any input harmonic background voltage.
The expected accuracy of the models is in line with the most accurate state-of-the-art harmonic models available, e.g. as shown in [10,17,18]. The main advantage of the proposed approach is the possibility to reduce the model to an analytic dependency on the supplied voltage, which leads to the improvement of the iteration process in load-flow based solvers for harmonic analysis, as it avoids additional iterations in the simulations. This will allow to exclude possible convergence problems with the algorithm in case of large system studies with a significant number of sources present.
As it will be shown in this paper, an accurate modelling of harmonic angles cannot be performed on the basis of Gaussian distribution of linear data but rather needs to be tackled by approach of circular (directional) statistics. Naturally, circular statistics methods lay a foundation for explaining relationships between dependent and independent angular variables.
This paper is structured as follows. Section 2 gives fundamental information about directional statistics, Section 3 describes the analysis process of measured data and proves suitability of chosen methods, Section 4 provides knowledge about used circular regression model, Section 5 is focused on the obtained modelling results. Some additional implications about modelling process are given in Section 6 and summary of the research is presented in Section 7.

Circular statistics for the analysis of harmonic phase angles
Owing its origins to medicine and astronomy and later finding applications in biology, geology and meteorology, circular or directional statistic methods have been mostly untouched by power system engineers and researchers. In the latter, the basic electrical parameters -currents and voltages are represented fundamentally by magnitudes and angles, making the linear approximations in analysis of such data possible, but not sufficient for routine data processing.
Approximate linearity may serve its purpose in large variety of power system analysis cases, but specific fields -for instance, analysis of harmonic distortions would require different approach. This is merely due to the fact that harmonic-producing sources exhibit non-linear behaviour and harmonic phase angles therefore can be analysed as two-dimensional orientations on a unit plane. The following subsections provide basic theoretical knowledge with respect to directional statistics and are mainly based on [19].

Fundamental parameters
The fundamental step in analysing directional sets of data is identifying whether or not measured set contains certain modal (directional) group or, perhaps, several modes (multimodal data). One of the most important quantities is a mean direction of circular data given as where IET Gener. Transm. Distrib., 2020, Vol. 14 Iss. 18, pp. 3826-3836 © The Institution of Engineering and Technology 2020 and θ i is individual angle measurement. Furthermore, the next descriptive quantity associated with the mean direction θ¯ is the mean resultant length R given by where n is a number of samples and R is a sum of variables of representative dataset defined on a cyclic interval [0, 2π]. Thus, the mean resultant length is the outcome of vectorial summation but not arithmetic. Its value lies in range (0,1) and R = 0 can potentially but not necessarily imply uniform distribution of samples around the circle. On the other hand, values close to R = 1 indicate rather large concentration of data points. In other words, R is an estimator measuring the concentration. These two quantities form angular and amplitude components of the first trigonometric moment ( 4) and trigonometric moment m 1 ′ is a basic population characteristic describing the underlying probability distribution function.

Probability distribution on the circle
From Section 2.1 it follows that circular data cannot be described by normal Gaussian distribution since all basic statistical parameters are presented in vectorial rather than scalar form. The most often used model is von Mises symmetric unimodal distribution for which dispersion equals to a concentration parameter k on contrary to the normal distribution, where dispersion is represented by the variance σ 2 . It is worth noting that increasing k value indicates rising concentration around reference direction and typically if k ⩾ 2 it is reasonable to fit measured data to von Mises distribution and apply corresponding statistical methods.
The probability density function is given as where k is a concentration parameter and μ is a mean direction and is the Bessel function. Finally, the distribution function of von Mises distribution is represented by If k approaches 0, the distribution converges to the uniform distribution whilst if k = 1, the distribution tends to the point distribution concentrated in the direction μ.
The evaluation of these functions is complex and is best done by means of commercial software packages. In this paper, we use R software for analysing circular data.
It is worth of noting that in some cases the von Mises distribution can be closely approximated by wrapped normal distribution which can be obtained by wrapping normal distribution around the circle [19].

Statistical tests
In order to decide whether or not the particular set of angular data can be fit to unimodal distribution, a number of tests must be performed. Test of uniformness against a unimodal alternative or Rayleigh test is a suitable formal tool for that purpose. The Rayleigh test statistic is then given as where R is a mean resultant length, θ is a mean direction and μ 0 is the mean direction of the alternative unimodal model. It is typical to reject the null hypothesis of uniformity if R 0 is large and exceeds the significance probability of R 0 calculated as where Z 0 = (2n) 1 2 R 0 , n is a number of samples and p z = Φ(Z 0 ) is table values which define percentiles of the normal distribution.
Other tests exist to verify specific alternative assumptions, for instance Watson test aiming to check null hypothesis of randomness against von Mises distribution.

Analysis of the measured data
In this paper, three distinct types of low-voltage PE are present, namely 1 kW single-phase PV inverter, Li-ion battery charger with rated input current of 16 A and a group of CFL lamps comprised of 36 pieces of different brands.
The measurements were performed according to [11] for every equipment type. Firstly, a harmonic source is connected to the undistorted voltage supply and the emission level at this condition is measured. Further, individual harmonics with magnitudes between 0.5 and 5% depending on the harmonic order are added one by one. Next, within each level of harmonic voltage magnitude, a phase shift is being varied between 0 and 360° with standard step of 30° and the harmonic current response of equipment under test is recorded. Finally, all the obtained data are organised into convenient look-up tables.
In order to investigate the shape of distributions of measured values, the harmonic currents were split into the magnitude and phase angle components. After that, harmonic phase angles were organised on a unit circle pane. Fig. 1 presents measured harmonic current angles of PV inverter in response to voltage background angles varied as 0, 30°, 60°, 90°, 120°, 150°, 180°, 210°, 240°, 270°, 300° and 330°. The composition of these plots is typical for representing angular data by means of circular statistics. On the shown unit panes every black dot corresponds to the one specific value of an angle of emitted harmonic current. Whenever the phase of harmonic current falls within the same sector multiple times, the dots are plotted on a 'stacked' manner. At the later stages this technique allows to assess the modes of distribution. The magnitude of voltages was kept at the same level. Furthermore, the recorded data were analysed with circular statistics methods as described in Sections 3.1 and 3.2.

Exploratory analysis
At this stage, the purpose of the initial graphical analysis is to merely assess the fundamental harmonic behaviour of equipment under test and to conclude if there are certain patterns in angle distribution. The black filled circles in Fig. 1 designate directions of harmonic phase angles. The prevailing direction is clearly visible for fifth harmonic, where points falling under the same group are drawn in a stacked manner. Furthermore, the spread of angles is increased for 7th and 9th harmonic and becomes nearly evenly distributed for 13th harmonic component.
For visual evaluation of the presence of certain directional modes, another non-parametric statistical instrument is applied. Fig. 2 shows kernel density estimates of the measured data. In general, the influence of each data point is stretched over reasonably small arc containing this data point. Next, the final density estimate for prevailing direction is then represented as the sum of all the contributions of smoothed points. The unimodal groups falling between 0 and 75° are clearly observable in Figs. 2a-c. Some irregularities in the form of additional bumps can be visible in Fig. 2d corresponding to 13th harmonic order.
Based on this initial analysis, it can be concluded there is an evidence of prevailing directions of harmonic phase angles as consequence of changing background voltage distortion angle.
Furthermore, once the presence of unimodal groups in measured data is established, a graphical assessment of goodnessof-fit for the von Mises model is performed. One way to check whether samples of data belong to a specified distribution is quantile-quantile (Q-Q) plots. Along the y-axis the input data is plotted versus the x-axis which is a collection of theoretical values of expected distribution. In case when the resulting plot shows signs of linearity it can be concluded that the data sample comes from the sought distribution. Fig. 3 demonstrates that studied circular distribution model is suitable for the representing measured set of data. There is little evidence of the departure from the model for 5th (Fig. 3a) and 13th (Fig. 3d) harmonic orders, however, Q-Q plots can be fluctuating and therefore formal statistical tests are required.
The exploratory phase of analysis was performed in a similar fashion for the other two equipment types and results are shown in Figs. 4-7. It is interesting to note the spread of harmonic current angles of battery charger ( Fig. 4). From these directional plots it is observable that for seventh and ninth harmonic the angles are distributed visually evenly, and this is confirmed by Figs. 5b and d. The latter kernel density estimate graphs demonstrate the presence of several modes of distribution. This in theory can point to the lesser goodness-of-fit, however, at this stage the data samples cannot be discarded based only on the visual analysis. It will be shown in Section 3.2 that formal mathematical tests identify the presence of one large mode. Moreover, the effect of even spread of the majority of current harmonic angles can possibly be referred to the notion of linearity studied in [5,20] and can be attributed to the  On the other hand, a group of CFLs demonstrates highly nonlinear behaviour with angles being concentrated strictly in certain quadrants (Fig. 6). As it was mentioned before, a current angle response belonging to the same sector is being stacked on the response of previous test scenarios of changing voltage background angle with a step of 30°. The fact is confirmed by kernel density plots (Fig. 7) showing quite narrow modal groups for selective harmonics.

Formal tests
Having concluded graphical analysis of measured angular data, it is reasonable to calculate all fundamental statistic values and perform tests of null hypothesis. In this section, the statistic quantities are calculated as it is explained in Section 2.1 and Rayleigh tests are performed as per Section 2.3. The results for every studied device are summarised in Table 1.
The calculated statistical quantities are in a good agreement with theoretical information of Section 2.1 stating that mean resultant vector length R close to 0 does not necessarily imply uniform distribution of phase angles. Looking at the rows related to battery charger in Table 1, it is visible that mean resultant lengths equal to 0.11 and 0.13 for 7th and 13th harmonic orders, respectively, correspond to large values of Rayleigh statistics. Moreover, the outcome of the Rayleigh test shows that significant probability values P 0 are exceeded by R 0 by several orders of magnitudes for these harmonics. It can be concluded therefore that unimodal directional groups are indeed present in the measured set of data, however, initial graphical analysis can be insufficient for drawing this conclusion.
On the other hand, R = 0.11 of 13th harmonic PV current angles corresponds to R 0 = 0.18 and is less than calculated significance probability P 0 = 0.29. This confirms the visual analysis presented in Fig. 2d and can indicate ambiguity in determining prevailing phase angle directions.
Finally, the results presented for a group of CFL in Table 1 implicate high level of data concentration in relation to prevailing directions for every harmonic order.

Circular regression model
The analysis of measured data which was carried out in Section 3 proves that current harmonic angles measured with fingerprint approach can be used as the input for circular modelling process. In fact, one of the advantages of treating data as a sample of von Mises distribution is the possibility to model a circular response variable with respect to the explanatory variable. In the case of harmonic modelling it means that within fixed background voltage magnitude level there will be a distinct response of equipment current angle to the change of applied voltage angle. In other words, a circular regression model can be derived which explains the systematic variability of harmonic current angles. This means that the response of estimated model for any input value of background voltage angle can be predicted.
A fundamental assumption in the modelling process is that all data points have been taken from von Mises distribution with specified prevailing mean direction. The general purpose of the model therefore is to estimate the mean direction with relation to where X i is the explanatory variable, β = (β 1 , …β k ) is the vector of coefficients to be estimated and g is the link function, subject to selection in advance. This model explains the behaviour of response variable when explanatory variables are in the form of linear covariates. Given the case when both response and independent variables, u and ν, are circular, the general regression model was modified by Jammalamadaka and Sarma [21]. The following model was proposed: E(e iν u) = ρ(u)e iμ(u) = g 1 (u) + ig 2 (u) (11) with e iν = cosν + isinν being the conditional expectation of estimated mean direction, μ(u) is the conditional mean direction of ν with respect to u and ρ(u) is the conditional concentration parameter for periodic functions g 1 (u) and g 2 (u). It is interesting to note that model proposed in [21] estimates not only mean direction but also the concentration parameter. The predicted value of mean direction is derived as The trigonometric polynominals are then utilised as the approximations of g 1 (u) and g 2 (u) where k is the trigonometric order and A k , B k , C k , D k are the parameters to be estimated.
Finally, the generalised least squares estimation algorithm is used for estimating parameters of (13). If ε = (ε 1 , ε 2 ) is the vector of random errors, then the following is in order: As it can be deduced from (14), the final circular regression model is represented by both cosine and sine terms. According to (12), the estimated mean direction is given as This model is conveniently integrated in the 'Circular' package of R software and is used in this paper. The proposed algorithm is summarised in the flow chart of Fig. 8.

Modelling results
This section presents the outcome of the modelling process described in details in Section 4. Only selected harmonics are discussed in the subsequent subsections. Refer Appendix for the complete tables of model coefficients. Fig. 9 shows the results of deriving circular regression model from the measured sets of data for selective harmonics produced by PV inverter. The presented scatter plots provide convenient way of visualising the dependence of one variable on another, i.e. the response of harmonic current angle to the changes of background voltage phase angle. In Fig. 9, the measured values (designated with coloured bullet •) are plotted together with fitted data (designated with coloured plus +).

PV inverter model
The regression models are estimated separately for different levels of magnitudes of voltage background distortion. It can be seen that harmonic injection patterns are significantly different for every studied frequency. Additionally, at fifth harmonic, the magnitude of the applied background voltage impacts to the great extent the pattern of harmonic current angle evolution.
On the other hand, nearly linear trend is observable for 13th harmonic, with all points being quite close to each other regardless of the magnitude of voltage.
All in all, the fitted model demonstrates good agreement with measured data. Moreover, in certain cases of manifested linearity (13th harmonic) the estimated harmonic angles overlap the measured ones.   Fig. 10 shows the residual analysis corresponding to the derived regression models, where fitted data are plotted against estimated residuals. Fifth harmonic plot (Fig. 10a) suggests there are some outlying points corresponding to −40 and 40°, but the majority of residuals is enclosed between −10 and 10°. A similar behaviour is observed for seventh ( Fig. 10b) and ninth (Fig. 10c) harmonics. Fig. 10d, on contrary, is characterised by the residuals mostly concentrated between −3 and 3° indicating good performance of the regression model.

Battery charger model
Figs. 11 and 12 demonstrate the fitted data together with measured values and model residuals for the battery charger. It is observable that for every plotted harmonic order the linear trends are explicit, with the exception of ninth harmonic. Furthermore some outlier points are clearly present in the 13th harmonic. Once again, the nearly linear behaviour of harmonic current angle response is in agreement with the analysis presented in Section 3.1.
The spread of residuals is mostly enclosed between −10 and 10° with some outlier points visible in Fig. 10d (13th harmonic).

CFL model
The studied group of energy-efficient lights is characterised by strong non-linear and asymmetric harmonic patterns. Fig. 13 allows to conclude that for 5th and 7th harmonic order the current angle response alternates drastically under the influence of different harmonic voltage magnitudes. On the other hand, this is not the case for 9th and 13th harmonic components. For the former it can be observed that all data points are roughly concentrated within one prevailing direction regardless of the magnitude of applied voltage. Furthermore, in agreement with Fig. 7d the dispersion of the current angles for 13th harmonic order is larger than for other frequencies, nevertheless, different voltage magnitudes levels do not influence significantly the response of harmonic current phase angles.
An analysis of the residuals presented in Fig. 14 shows that most of the points lie between −2.5 and 2.5° suggesting that measured set of data is represented well by derived circular regression model. An exception to this is however, 13th harmonic with its residuals being enclosed between −10 and 10°.

Discussion
One of the fundamental assumptions allowing to make progress in modelling method described in this paper is that all measured values fit to the underlying von Mises distribution for circular data. While this distribution is a natural choice for many theoretical problems, it can be difficult to use in practice, in particular when concentration parameter k associated with measured data set is less or approaching the value of 2 (see Section 2.2 for the details about von Mises distribution). However, since circular statistics includes vast varieties of methods and analysis capabilities, it can be reasonable to estimate confidence intervals associated with estimated prevailing direction of current harmonic angles and uncertainty of concentration parameter. For the specific case studied in this paper, i.e. when the amount of samples in every dataset is small, the so-called parametric bootstrap method is in order.
The aforementioned is closely related to the explicit symmetrical patterns visible on some of the residual plots, for instance, in Figs. 10b-d and 14d. This indicates that underlying von Mises distribution and consequently the chosen regression model does not fully describe the behaviour of measured set of data. Besides improving the model, one of the ways to increase accuracy of the modelling method is an application of outlier filtering technique. Calculating confidence intervals, improving the model and analysing the outlier problem is left for a future work.
Furthermore, since proposed modelling method concerns derivation of algebraic equations with respect to the harmonic phase angles, a distinct procedure for handling harmonic current magnitudes must be applied. An influence of circular variable (voltage phase angle) on a linear variable (current magnitude) This allows to approximate current magnitude as a function of background voltage by simply applying interpolation or appropriate curve-fitting procedure. As an example, Fig. 15 shows an evolution of fifth harmonic current magnitude with respect to the voltage angle as measured at the terminals of battery charger. Moreover, system studies which analyse large networks often require the use of a relatively large number of harmonic source models. Power-flow based calculation engines use iterations to reach a solution, which can be further complicated, and have their convergence jeopardised, in case of additional iterations of each source model due to e.g. frequency coupling [4]. For this reason, an analytic relation between the supplied voltage and the injected harmonic current is a very useful simplification to use in studies on large networks. This paper proposes one solution to this problem, with an approach to develop a source model for any type of load which meets specific preconditions described in this paper. Another example conveying the idea of modelling simplification while retaining high accuracy of non-linear model was recently reported by Laurano et al. [22].

Conclusion
In this paper, a new method for modelling harmonic producing sources is proposed. Three distinct power electronics devices have been modelled: single-phase PV inverter, battery charger and CFL. For this method harmonic fingerprint measurements are used for characterising current emission of modelled equipment as a function of applied harmonic background voltage.
The measured sets of data are separated into magnitude and phase angle components. The current phase components are then analysed with the methods of circular statistics, proving that they can be described by prevailing angle directions when the full range of background voltage angles has been applied at the terminals of equipment.
Theoretically, angular datasets characterised by prevailing angle directions cannot be analysed using linear distribution functions, for instance Gaussian distribution. Therefore, a unique von Mises probability distribution for circular data is proposed as underlying model. The regression technique based on von Mises distribution provided algebraic relationship between explanatory variablevoltage phase angle and response variable -harmonic current angle. This relationship is valid within certain level of harmonic voltage magnitude to which regression technique is applied.
The performance of the circular regression models was checked against the measured values in Figs. 9, 11 and 13. The graphs show satisfactory matching with the experimental data, despite the low number of parameters required by the model. Based on the available scientific works, e.g. models shown in [3][4][5][6][7][8][9][10], it can be concluded that the reduction of model complexity does not lead to poor performance in comparison with the state-of-the-art models used in the literature, as long as the devices in question satisfy the given assumptions of angle distributions.
Generally, the accuracy of these models depends on how well theoretical von Mises distribution fits to the measured data. In most of the cases this is attributed to the underlying equipment topology which in turn affects harmonic behaviour. The results of this work demonstrated that derived models are characterised by fair performance with residuals falling between −2.5 and 2.5° in best case and between −10 and 10° for the worst situation. The implications about improving the quality of the models were briefly discussed together with a method for handling harmonic current magnitude during modelling process.

Acknowledgments
This work has received funding from the European Union's Horizon 2020 research and innovation programme MEAN4SG under the Marie Sklodowska-Curie grant agreement 676042. The authors thank Prof. Jan Desmet and the research staff of LEMCKO laboratory, Ghent University, for providing access to the laboratory facility and supporting experiments.