Stochastic Residential Harmonic Source Modeling for Grid Impact Studies

With the introduction of more non-linear loads, e.g., compact fluorescent lamps, electric vehicles, photovoltaics, etc., the need to determine the harmonic impact of the residential load is rising, illustrated by the many studies performed on their harmonic impact. Traditionally, these studies are performed for a single new device and single penetration level, neglecting the harmonic interaction between new types of devices, as well as giving little information at which moment in time possible problems may arise. A composite approach to access the impact of harmonic sources on the distribution network is therefore proposed. This method combines a bottom-up stochastic modeling of the residential load with harmonic measurement data and harmonic load-flows all based on a scenario analysis. The method is validated with measurement data and shows a good prediction of the current level of harmonics in a residential neighborhood for the current situation. To demonstrate the applicability of the proposed method, case studies are performed on the IEEE European Low Voltage Test Feeder. These case studies show a marked difference between applying individual device-based models and a composite modeling approach, demonstrating why the proposed approach is an adequate method for the determination of the impact of new devices on the harmonics.


Introduction
Household appliances based on power electronics (non-linear loads) are increasingly connected to the distribution network due to the rapid development of semiconductor technology. Although an individual device injects an insignificant amount of harmonic current in the network, the sum of a large number of such appliances could result in a large amount of harmonic distortion [1]. In the future, the penetration of power electronic-based loads is expected to increase due to the increase of new technologies and the changing of conventional appliances to power electronic-driven appliances, e.g., compact fluorescent lights (CFLs) and electric vehicle (EV) chargers. On the other hand, some devices also are producing less harmonics due to newer requirements of standards (e.g., PV) and reduced nominal powers. The injected harmonic current may cause serious problems, such as overheating and failure of equipment, errors in the operation of sensitive devices, false tripping of protection relays and interference with communication signals. To avoid the potential risks that affect the system operation and jeopardize the security of supply, harmonic impacts and propagation studies of large distribution network need to be performed for future scenarios. Distribution system operators (DSO) need to forecast how increased nonlinear loads will contribute to the harmonic distortion level and when the distortion level may exceed the regulation limit. DSO can make use of the result to take suitable measures to avoid deterioration and improve the power quality (PQ) level.
Conventionally harmonic impact studies have been performed with increased usage of a specific type of non-linear device. The harmonic impact due to CFLs was studied in [1][2][3]. In [4,5], the effects of residential photovoltaic installations were analyzed. With a high penetration of EV charges, the harmonic distortion is discussed as shown in [6,7]. Most of the research focuses on the harmonic impact by one single type of device. However, in real life, not only a single new harmonic source will be placed in the network, but many different types of harmonic sources will be added. The effect of these multiple sources is not simply the combined effect of the individual devices. A composite approach, which takes into account different new types of harmonic sources at the same time, is therefore required. By applying such an approach, the future harmonic distortion in the network can be better estimated.
A first attempt to evaluate the different types of new harmonic sources has been made [8]. Here, a detailed harmonic analysis has been conducted on actual circuits with mixed harmonic loads, which may be projected on the distribution network in 10-20 years. The base harmonic distortion is based on measurement data and remains constant throughout the analysis, while in reality, this would be subjected to changes, as well; moreover, the general framework of how these loads should be modeled and integrated with the conventional load is missing. Therefore, a composite framework for the harmonic load modeling needs to be defined. This framework should be able to include different types of harmonic sources whose penetration levels are determined by different scenarios; the aging, removing and replacement of the existed appliances (appliance development) are incorporated on the time series, and the efficiency improvement of the new connected appliances should be involved. In this paper, such a framework is proposed. The proposed framework builds on previous research on the modeling of the residential load [9,10] and the modeling of the harmonic sources [11]. With this framework, future harmonic impact studies can be performed; and not only the combined effect of many harmonic sources can be assessed, but also the time evolution of the level of harmonic distortion can be estimated.
This new framework for harmonic impact studies is presented in the following manner. In Section 2, the basics of the modeling of the harmonic load are discussed, starting with the bottom-up stochastic load modeling for both current, as well as future loads, followed by a discussion on the applied harmonic spectra and the harmonic aggregation. This framework is validated based on measurement results of a residential network in Section 3. Hereafter, the advantages of the proposed framework are shown in Section 4 based on case studies on an LV-feeder. The conclusions about the approach and applicability of the proposed framework are given in Section 5.

Harmonic Injection Modeling
To model the harmonic current injected by a household, a closer look needs to be taken at the injections of the individual appliances. As the harmonic current injection of a household is dependent on the type of appliance that is activated, an accurate description of the harmonic current injected for the different appliances is a vital first step in the modeling of the household harmonic current injection. The main difficulty with defining these injected currents is the determination of when and for how long a certain appliance in the household is switched on. To tackle this problem, a bottom-up approach is adopted to model the households appliances [10]. Subsequently, measured harmonic spectra of household appliances are used to obtain the model of the current injection of a household.
As the behavior of the household members determines when an appliance will be switched on or off, the modeling of the current injection of a household starts with the household occupancy modeling. With a database of time use surveys (TUS) and a library of detailed appliance models, simulations are carried out to get the appliance usage patterns as a sequence of on/off states of each appliance. Harmonic spectra of various appliances are used in the harmonic load flow calculation to obtain the harmonic injection currents at each point of connection (PoC).
Since the occupancy plays a significant role in the electricity usage through the interactions between occupant and appliance, an improved version of techniques based on Markov chain Monte Carlo (MCMC) simulation is applied [10,12]. Allocating the appliances in a household can be completed with the database on the appliance ownership, the wealth level and the dwelling size [10,13]. The harmonic spectra of different household appliances are either measured in the laboratory or collected from the PANDA database [14], so the harmonic emission level of a specific appliance under laboratory conditions is obtained. To obtain the harmonic injection level of a specific household, the harmonic spectra are utilized as described in [15]. The scaling of magnitudes and shifting of phase angles have been done according to the recommended formulas.
With the knowledge of the occupancy and the set of allocated appliances, a more in-depth analysis of the appliance usage pattern can be obtained through simulation. Another Markov chain is applied to model the on/off status of each appliance. The weather variables and the time of the day are taken into consideration, as well.
In Figure 1, the flow chart of the methodology is illustrated. The detailed explanation of each of the steps is given in the following subsections.

Occupancy Modeling
The Dutch national statistics agency, Statistics Netherlands [16], provided a TUS of 2042 respondents' behavior with a 15-min resolution. Figure 2 shows the occupancy profile of three individuals of different households during two consecutive days. It can be seen that the behavior of an individual differs significantly from others. These differences in behavior will result in different usage of electrical appliances. Therefore, the model should incorporate these variations in the occupancy. A Markov chain Monte Carlo method is used for the occupancy modeling, which is explained in [17]. The Markov chain model refers to the fact that at each time instance depending on the current state and the associated probabilities, the model switches to another state or remains in the current state.
As shown in [10], the occupancy can be modeled for the purpose of determining the energy consumption within a household with just two states: active (and at home) or inactive (sleeping, away from home). By applying logistic regression to the time use survey data [12], a transition probability matrix is constructed for week and weekend days, the time of day, the number and the age of the household occupants. An individual Markov chain is assigned to each member of the household, and the probabilities within this Markov chain vary every 10 min within one week. The Markov chain modeling of the occupancy is defined as: where Oc(t) is the occupancy at time t, U (0, 1) is a random variable from a uniform distribution between zero and one and the probability P (Oc(t)|Oc(t − 1)) is the probability of switching the state of the occupant taken from the transition probability matrix, which only depends on the previous value. As described above, the probability of occupancy change in a specific group (particular age and household members) is estimated based on the data from TUS, which are used for the construction of the transition matrix. For weekdays, the occupancy of the same time interval of the day before is taken into account in the determination of the transition matrix, as it is found that the occupancy data between two weekdays that are consecutive are highly correlated.

Allocation of Household Appliances
The allocation of household appliances is modeled with data on the penetration level of different types of appliances. It has been shown that the penetration levels of many appliances differ with respect to income and dwelling size [18]. Information on the penetration of different appliances in a household is estimated from Statistics Netherlands, which has for a number of appliances the penetration grade differentiated in regard to neighborhood wealth level and dwelling size. The logistic regression analysis is applied to model the dependency of the income and housing size on the appliances that are present in the household. The appliances for which the necessary data are not given are allocated according to the data on similar types of devices. The allocation of household appliances is performed by (2).
where A i is the i-th appliance, A H is the set of appliances present in the household, p A is the fraction of households with appliance A, r h is the household size adjustment factor, r w is the household wealth adjustment factor, and N A is the total number of appliances. The values of r h and r w can be adjusted between 0.5 and 2.

Appliances Modeling
To determine which appliances are switched on in the model, different modeling approaches are applied for different categories of appliances: base appliances, night appliances, temperature adjustment appliances, lighting appliances, behavioral appliances and other appliances. The list of appliances modeled in this paper can be found in Table A1 in Appendix A. The model should have a 10-min time step or shorter since it is the generically adopted time interval in PQ standards (e.g., EN50160); therefore, the modeling is done based on a 1-min time frame.
To cope with the development of the residential loads when the penetration level of technologies, like modern lighting (e.g., CFL), EV, PV, etc., increases, the modeling of these appliances is required. The usage pattern of EV is obtained with a mobility survey and the household occupancy. The generation output of PV is calculated with the weather information, the rooftop area database and the solar panel orientation. Figure 3 shows the approach to model the EV and PV. For the main other sources of harmonics, the load models from [10] can be used.

Base Appliances
Base appliances are those loads that are independent of occupancy. The presence of these appliances can be modeled with constant power.
where A pr (t) represents the presence of the appliance in the household at time t. t 0 is the initial time shift, obtained from U (−30, 30) randomly to consider the diversity of different appliances. p t,A represents the fraction of time appliance Ais on. The time unit is minutes, which are represented as min in the equation.

Night Appliances
Night appliances are assumed to be on at an exponentially-distributed time start after 21:00 when the off-peak tariff starts.
where X (λ) is the random variable from exponential distribution, t A,r is the average running time of appliance A and t A represents the time after which the appliance A became active. The 24 h equals 1440 min, and the 21 h is 1260 min.

Temperature Adjustment Appliances
Temperature adjustment appliances are switched on probabilistically with the information of the temperature and the wind speed, in which a heating degree approach is applied [19]. The heating systems are assumed to be switched on when either at least one occupant is active or the adjusted temperature of wind speed T w (t) falls below zero degrees and the previous operation of the heating installation was more than 2 h ago.
where T w is the wind speed adjusted temperature, T h represents the heating/cooling temperature, and t A,m is the minimum run time of appliance A.
The cooling system can be modeled with a similar approach.

Lighting Appliances
Lighting appliances can be modeled with the data of the solar irradiation and the household occupancy. The lighting system is assumed to be on when the irradiance is not more than 60 W/m 2 [20]. The usage scale of those appliances depends on the number of residents in the household at the moment [21].
where I(t) is the irradiance at time t.

Behavioral Appliances
Behavioral appliances are modeled based on statistical data. A Markov chain with a transition matrix that changes over the time is used to model the presence of these appliances.
where P (A pr (t) = 1|A pr (t − 1) = 0) is the probability that the appliance is switched on when it was off in the previous time instance.

Other Appliances
Other appliances are switched on when there is at least one active occupant. They are either used seldom or shorter than 10 min.

Electric Vehicle
For the EV, the most important modeling from this perspective is the state of charge and whether the EV is at home. If the household has an EV based on the analysis of the penetration level, the EV can either be at home or away from home. When the EV is at home, the charging process may start if the battery is not full. As shown in Figure 3, a mobility survey is used to derive these variables, including home arriving time, departure time and daily driving distance [22,23]. These inputs, as well as household occupancy are used to determine when the EV is present at home. When the EV has arrived home, the charging strategy needs to be modeled. This charging is based on a plug-in and charge strategy as described in the second part of this section. Presence of EV at home: assuming the probability that which occupant drives the EV is equal, the probability that the EV is at home is defined as: where p EV,in (t) is the probability that the EV is at home at t, N po (t) represents the number of occupants in the household at time t, N p is the total number of occupants in the household and P av (t) is the probability of EV arrival at time t. Charging mode of EV: the battery state of charge (SOC) at the time of arrival can be determined with the driven distance, as shown in (10).
where L t is the distance driven with a full EV battery and L d represents the driven distance. The required charging duration (t ev,r ) can be calculated as: where C EV is the EV battery capacity, DoD represents the EV battery depth of discharge, η c is the EV battery charger efficiency, and P R,EV is the charging power. Fully discharging of the battery may cause serious and irreparable damage to the battery. Therefore, the depth of discharge, which refers to the fraction of power that can be withdrawn from the battery, needs to be specified. With the information of t ev,r , charging strategies can be implemented [23,24]. In this study, it is assumed that all of the EVs just start charging within an hour after the EV returns to the household, as shown in (12).
where t 0 ∼ U (0, 1) is implemented to model a random start in the next one hour after the driver with the EV has arrived home.

Photovoltaic System
The PV system of a household is modeled based on the roof area of a household from which the amount of PV panels can be placed is known. Further, the output of the PV system is obtained with the data of solar irradiance (weather conditions) [25]. The greatest potential of PV installation can be assessed with this approach.
The level of solar irradiation is the most significant environmental condition for the energy production of the PV system. Besides, it is also affected by the tilt and azimuth (orientation) angle [26,27] of the PV panel. The radiation on a tilted surface (I T ) can be obtained by: where I b is the beam irradiance on a horizontal surface, R b is the ratio of beam radiation on the tilted surface to that of the horizontal surface, I d the diffuse horizontal irradiance (DHI), β is the tilt angle, I g is the global horizontal irradiance (GHI), and ρ is the ground reflectance. R d is the ratio of diffuse radiation on the tilted surface to that of the horizontal surface. Multiple models are available to get R d and obtain I T [27]. In this paper, the irradiance model called the Liu and Jordan model, which assumes an isotropic diffuse sky, is applied. The detailed calculation process is shown in [28], and the radiation level can be defined as a function of two angles (β and the orientation angle γ).
The optimal angle should be generated for the location of the PV system, and a band around these angles can be employed to generate the expected variability.
With the measured GHI and DHI, the tilted irradiance is obtained with (14). The amount of irradiance in a day is mainly dependent on climate conditions. A sampling approach is employed to determine the irradiance in the future, as the climate conditions do not change much over time.
To generate an adequate database to sample from, solar irradiance data of GHI and DHI of each day, as well as the five days before and five days after a day in the past 40 years are regarded as samples of that day in the year. Therefore, a database of 440 samples of time series solar irradiance is available for each day in the future (i.e., the solar irradiance of 1 September 2016 is randomly selected from that of 27 August-6 September from 1976-2016).

Harmonic Current Spectra of Appliances
The harmonic spectra of over 20 typical household appliances have been obtained. For the appliances that have high levels of the current harmonic distortion, e.g., CFL, TV, appliances of different brands are measured so that more than one harmonic spectrum is available. For the appliances whose operating conditions can be considerably different (e.g., vacuum cleaner), the harmonic spectra with both high and low power output are obtained, as well as the running and the standby states. With the solution of the fundamental frequency power flow, I 1 ∠θ 1 , which is the fundamental current injection, are obtained. The harmonic magnitudes and phase angles (I h and Θ h ) of h order are calculated according to [15].
where I h−spec is the magnitude of h order from the spectrum and Θ h−spec is the phase angle of h order from the spectrum. Formula (15) scales the measured harmonic current spectrum of a typical appliance with the fundamental current (I 1 ) of the used appliance. Formula (16) adjusts the spectrum according to Θ 1 . The h-times shifting of h-th order harmonics to the fundamental harmonic is also taken into account.

Attenuation and Diversity Effect
In [11], it is concluded that the attenuation effect becomes significant when the total harmonic distortion of voltage (THDV) as calculated in (17) is very high (exceeding 10%), which is unlikely to be observed in the LV network [29]. The diversity effect is due to the differences of harmonic current phase angles. Applying (16), the variation of phase angles due to the variation of the supply voltage phase angles and the phase angle differences because of the different positions on the feeder are incorporated. In the future scenario, the voltage may be too distorted to ignore the attenuation effect. As shown in [30], the attenuation and diversity effects are modeled as: where V h is the h-order harmonic voltage and V 1 is the fundamental voltage, α h (THDV) gives the magnitude of the hth harmonic current as a function of the THDV, and β h (THDV) gives the phase angles. The functions are determined with curve fitting technique based on measurement data. Thus, the injected harmonic current by the harmonic source at bus k is defined as: An iterative frequency domain harmonic analysis method needs to be applied to obtain the solution of (20) and (21). The current spectra of the harmonic sources are revised according to (18) and (19) with an iterative process, which is illustrated step by step in [30].

Appliances Power and Harmonic Injection
Once the usage patterns of the appliances are derived as described above, the power and harmonic injection need to be assigned when the appliance is on.
In (22), when the appliance is on (A pr (t) = 1), P R,A , the running power of the appliance, and H R,A , the harmonic spectrum of switched-on appliance, are assigned. When the appliance is off (A pr (t) = 0), P S,A , the standby power of appliance, and H S,A , the harmonic spectrum of the switched-off appliance, are assigned. For the appliances without standby state, P S,A and H S,A are considered to be zero. Figure 4 shows the equivalent circuit for a household with the calculated impedance Z house and the harmonic injection current I sh . Z house is calculated with the fundamental current of load, and I sh is determined with the procedure described in this paper.  With the measured harmonic spectrum as shown in Figure 5, a TV is taken as an example of the harmonic current modeling in a household with three persons. Figure 6a shows the occupancy profile of the household obtained from the occupancy modeling process. The switching on/off probability of a TV is generated from TUS data and is shown in Figure 6b. P (TV pr (t) = 1|TV pr (t − 1) = 0) is the probability that TV is switched on when it was off in the previous time instance, and vice versa, P (TV pr (t) = 0|TV pr (t − 1) = 1) is the probability that TV is switched off when it was previously on. When no active occupant is present in the household, the probability is shown in a transparent manner, because there is no one who could switch on the TV. A random number is generated from U (0, 1) and compared to the probability value to determine the on/off state, which is shown in Figure 6c. The power of a TV (P R,A ) is obtained from a uniform distribution to consider the different TV types, and then, the fundamental, 3rd, 5th and 7th harmonic currents of TV in this household are shown in Figure 6d. For other appliances, a similar procedure is applied. Although the harmonic currents are only presented up to seventh order in Figure 5, the modeling has been done up to 40th to be in accordance with the standard [31].

Appliance Development Modeling
The development in appliances needs to be implemented, which deals with aging, removing, replacing and gaining new household appliances. The overview of the modeling approach is shown in Figure 7. The changes in the composition of the household needs to be incorporated at first. Since the probability matrices for the occupancy are categorized into different groups according to occupants' age and household size, implementing future scenarios can be achieved by changing the used probability matrices to the matrices of another group. The changing of the household composition is reflected by the change of occupants' age distribution versus years, and the probability metrics are updated each half year based on the new age groups in the model.
The lifetime of each appliance (l t,A ) is assumed to be given by N (L t,A , L t,A /3), where L t,A is its average reported lifetime. At a predetermined number of intervals, the appliances' remaining lifetime is calculated. The list of available appliances will be updated since appliance ratings will change (e.g., CFLs replace the existing incandescent lamps, and the size of televisions is increasing). If an appliance remaining lifetime is zero, the appliance should be removed from the household, and there is a chance that a new appliance of the same type or with more advanced technology is allocated to the household with the consideration of the penetration level.
Simultaneously, an additional appliance is possible to be obtained by a household based on the household characteristics and penetration level. If the penetration level of an appliance changes, an additional chance of either adding or removing the appliance from the household is performed. Figure 7 shows an example that a new TV will be assigned in the house since the old one is broken (l t,TV N (L t,TV , L t,TV /3)). At the same time, an EV as a new appliance is gained by the owner since U (0, 1) ≤ p EV × r h × r w . However, the development of power supplies of devices cannot be known for sure.

Validation of the Harmonic Injection Model
To assess the performance of the modeling approach, field measurements at the secondary side of the MV/LV transformer from a residential network in the Netherlands have been recorded, and the simulation results are compared to these measurement data. The household size adjustment factor r h and the household wealth adjustment factor r w are set to be the average value 1.2.

Residential Network Modeling and Simulation Approach
A residential LV network in the Netherlands is selected for the validation of the harmonic modeling approach. The network provides the electricity to 383 residential customers through a 630-kVA (10/0.4 kV) transformer via seven outgoing feeders of underground cables. The network topology is shown in Figure 8, and the network parameters can be found in [32]. As shown in the figure, the measurement campaign has been implemented at the secondary side of the transformer. Each household is represented as the current injection model since the bottom-up procedure has given the harmonic injection value at each order. The background harmonics are represented by a voltage harmonic source in Figure 4. The parameters of the Thevenin equivalent can be estimated with two operating conditions according to [33]. Assuming no change in operating conditions in the grid-side LV network at the same time of different days, the parameters for the Thevenin equivalent circuit can be obtained with the procedure shown in [33]. The comparison between the measurement and simulation results has been done as presented in Figure 9.

Measurement and Simulation Results
A whole week of measurement data has been compared with the mean values of the results from the Monte Carlo simulations (N MC = 1000). Figure 10 shows the magnitudes of the measured (| I m h |), the mean and the 5-95 percentile range of the simulation results (| I s h |) for fifth and ninth harmonic currents and the current calculated up to 40th order during two days (Sunday and Monday). The differences between | I m h | and the mean of | I s h | are generally lower during the night than during the day since more appliances are used during the day. Harmonic currents of higher orders are found to be mostly very low (less than 1 A). From Figure 10c, it can be seen that the model gives results of the total harmonic distortion of current (THDI) close to the measured values. Next to the magnitudes, the performance of the model with respect to the harmonic phase angle is also investigated. Figure 11 shows the phase angles of fifth and ninth harmonic currents relative to the phase angle of the fundamental current. The mean of the simulated phase angle fits well with the measured one. However, a lower variation is shown for the simulated curves than the measured one because of the aggregation effect of the 1000 Monte Carlo simulations. The THDI is calculated using (17) when V h is replaced by I h , which is the harmonic current of order h, and V 1 is replaced by I 1 , which is the fundamental current. The RMSE is calculated between | I m h | and the mean of | I s h | on time series of the week for both magnitude and phase angle from fundamental to 40th harmonic. Due to the decrease of the injection level with the increase of harmonic order, the RMSE values of magnitudes are normalized by applying (23) in order to investigate the performance difference among different harmonic orders. The results in Table 1 indicate that the error becomes more significant when the harmonic order is getting higher for both the magnitude and phase angle. The value of NRMSE for the THDI is about 22.1%, and the value of RMSE for THDI is about 5.7%.

Case Study
Case studies for the future are performed on a benchmark of the LV network as depicted in Figure 12 to show the possible application of the approach and investigate the impact of the development of residential loads. The benchmark network is known as the IEEE European LV Test Feeder, which is operated in a radial structure [34]. There are 55 household loads in total on the feeder. For the EV, a mobility survey carried out in the Netherlands [35] is used to generate the synthetic driver profiles as described in [23]. Figure 13 shows the probability of the arrival times of drivers in one day. For the PV generation, since the roof area of the houses in the European LV Test Feeder is unknown, an estimate is generated using the statistics from a Dutch distribution system operator (DSO). The dataset includes the rated power values of PV installations of many households. With the information of solar irradiance and the installed rated power of PV, the PV generation curve of each household is obtained. In the Netherlands, the PV module should face south to get higher solar irradiance. Therefore, the orientation angle (γ) is assumed to be between southwest (SW) and southeast and β is between 30 • and 60 • [27]. The solar irradiance database of last 40 years is obtained from the KNMI [36].

Scenarios
Future scenarios need to be implemented. The scenario of appliance efficiency (η A ) can be obtained from [37], which studied the energy savings potentials in European Union countries. The expectation of Dutch economic growth (Π E ) is in [38], which shows the scenario up to the year 2040. In [39], the scenarios on the penetration level of EV and PV (ρ EV and ρ PV ) are discussed, where three drivers have been concerned: a more nationally-or internationally-focused policy, a more economicallyor environmentally-oriented society and the economic growth. According to those, three typical penetration level scenarios are defined as the low, medium and high penetration level scenarios. For the low scenario, it is the output that the number of households does not increase and the electricity demand grows slowly, the investments in new technologies and the renewable energy production growth are low. For the high scenario, the number of households grows with 1% per year, the investment in sustainable energy technologies increases much and the share of distributed electricity generation grows due to the high economic growth. The medium scenario is in the middle of the high and low scenarios. Therefore, the penetration levels of PV and EV of the three scenarios up to 2050 are shown in Figure 14.
The parameters of (2) are updated for the future analysis, the factor f A changes based on the scenarios on the penetration of technologies; r w is assigned based on the income level of the neighborhood and subsequently based on the economic growth, as well; the factor increases the chance of a household having non-essential appliances; and the values P R,A and P S,A change with the efficiency scenarios. In the case study, the methodology is applied to a scenario with a medium level of η A , a high level of Π E and high values of ρ EV and ρ PV .

Harmonic Analysis for Compliance Assessment
With the defined scenario, the simulation of one week has been done at the benchmark network in both winter and summer. Figure 15 shows the 95 percentile value of voltage harmonic distortion at the secondary side of the MV/LV transformer in one week up to 40 years. It can be seen that the voltage is more distorted in winter than summer due to the seasonal differences in appliances usage; lighting is used more often, and variable speed drives are used more in pumped central heating. With the defined scenario, the values of individual harmonic distortion (third and fifth) and THDV all grow with time. It also can be predicted that the distortion level could exceed the regulated limits at a certain point in time. In future studies, different scenarios can be implemented from which the harmonic impact can be analyzed with regard to the aspect of regulatory compliance.
Year EN IEC IEEE Limit

Harmonic Distortion throughout the Network
Since the harmonic distortion levels are distinct from each other at different nodes on the grid, investigating the harmonic distortion value at one specific location may not be sufficient. The approach stimulates the harmonic current injection of each load, which implies that the harmonic sources over the entire network are known; the harmonic distortion levels over the entire circuit can be calculated rather than only the secondary side of the MV/LV transformer. Figure 16 shows the 95 percentile value of THDV at each node of the network and the root mean square (RMS) value of total harmonic current as defined in (24). The results are obtained from the simulation on the year 2016. The increment of the harmonic distortion along the feeder is better observed through the color plot of the simulated harmonics. It can be seen from the example in Figure 16 that THDV at the end of the feeder is 12.5% more than that of the feeder head. Therefore, it may lead to an underestimation if only looking at the secondary side of the transformer for regulatory compliance. The visualization tool is also applicable to identify the significant harmonic source and the location exceeding the higher distortion level. This is important when it comes to the design of mitigation measures, such as the filter locating and the isolation of sensitive loads [40,41].

Losses Due to Harmonic Distortion
Applying the methodology, the network losses due to harmonic distortion can be assessed. In the Netherlands, the three-phase four-wire system is used in the LV network. The losses on the neutral conduct should be calculated since it is mainly due to the PQ phenomenon (unbalance and harmonic). Figure 17 shows the 95 percentile fundamental and harmonic losses both in the phase conductors and neutral conductor where the values are expressed as the percentages of consumption power.

Impact of Specific Appliances
The approach can be used to gain more insight into the harmonic impact of specific appliances. For example, the increased penetration of EV for the sake of less environment pollution and the replacement of incandescent lamps with energy saving lighting system can be assessed with the proposed methodology. Figure 18 shows the impact of the EV, CFL, EV + CFL to THDV when the penetration level moves from 0%-100%. The penetration level of CFL means that a certain percentage of incandescent lamps is replaced by CFL lamps. The superposition of EV and CFL is also given in the figure, which has higher values than EV + CFL. This is due to the cancellation effect between EV and CFL, which becomes more and more significant with the increase of penetration levels. More exactly, the value of THDV when the penetration level of CFL is 100% is about 4.35%, and the value when the penetration level of EV is 100% is about 4.9%. However, the CFL and EV can be both connected to the network. If only the individual analysis result is obtained, the combined effect of EV and CFL has to be obtained by the superposition of the two values, which is about 5.3%. Applying the composite approach, the EV and CFL are both simulated, so the harmonic interaction between different appliances are incorporated. The distortion line in the figure shows that the THDV level when the penetration levels of two appliances are 100% is about 5.1%. The integration of two devices has made a difference of about 0.2%. If more appliances are simulated, more significant differences are expected. 20 40 60 80 100 The penetration rate of appliances (%) 3.5

Conclusions
The proposed approach can be applied to model the harmonic pollution from all kinds of residential loads through a bottom-up procedure, e.g., CFL, EV and PV. By applying the method, the combined effect of different harmonic sources in LV networks can be assessed for future load conditions. In the modeling of the future household harmonic sources, not only new devices are added, but also the evolution of the current harmonic sources within a household are modeled. The replacement of old devices with more efficient power electronics-based devices can in this way be taken into account for the harmonic injection modeling. The future penetration level of a particular appliance is determined with other relevant statistics instead of an assumption, which leads to more realistic results, as well as it gives the ability to estimate when problems with harmonics can be expected. The three most important characteristics of the framework for the future harmonic load modeling are: multiple harmonic sources, scenario implementation and the appliance development. The advantages of the proposed methodology are illustrated through a case study with certain scenarios. It shows that the harmonic distortion level grows by the year. Different applications, e.g., regulatory compliance, the harmonic distortion through the whole network, the harmonic losses and the influence of a specific appliance, are demonstrated. The circumstances that lead to the exceeding of the harmonic distortion limits can be predicted. The differences in the harmonic distortion within the network can be calculated, which simplifies the installation of mitigation equipment. The losses due to the deterioration of PQ are assessed, together with the regulation compliance analysis; the necessities to include the PQ issue in the future network planning procedure can be shown. The comparison of simply adding the effects of multiple sources, compared to an integrated analysis of the development of the household harmonic injection, shows clear differences in the resulting harmonic distortion, due to the harmonic cancellation. The harmonic distortion level increases with the growth of the penetration level of CFL and EV; indicating that the application of a composite approach as presented in the paper would give a better assessment of the expected levels of harmonics in the network.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. The list of appliances.

Category Appliances
Base Appliances Charger, cordless phone, freezer, modem, refrigerator, water bed.
Night Appliances Dishwasher, electric Boiler, electric blanket.
Temperature Adjustment Appliance air conditioner, electric heater, electric radiator, fan, mechanical ventilation, pumped central heating system, underfloor heating.
Other Appliances Coffee machine, drill, electric kettle, hairdryer.