Measurement-Based Current-Harmonics Modeling of Aggregated Electric-Vehicle Loads Using Power-Exponential Functions

This paper presents an aggregate current-harmonic load model using power exponential functions and built from actual measurement data during the individual charging of four different fully electric vehicles. The model is based on individual emitted current harmonics as a function of state of charge (SOC), and was used to deterministically simulate the simultaneous charging of six vehicles fed from the same bus. The aggregation of current harmonics up to the 11th was simulated in order to find the circumstances when maximal current-harmonic magnitude occurs, and the phase-angle location. The number of possible identical vehicles was set to four, while battery SOC, the start of charging, and the kind of vehicle were randomized. The results are presented in tables, graphs, and polar plots. Even though simulations did not consider the surrounding harmonics, supply-voltage variation, or network impedance, this paper presents an innovative modeling approach that gives valuable information on the individual current-harmonic contribution of aggregated electric-vehicle loads. With the future implementation of vehicle-to-grid technology, this way of modeling presents new opportunities to predict the harmonic outcome of multiple electric vehicles charging.


Introduction
In recent years, energy policies and the desire to reduce air pollution have resulted in increased reception of electric vehicles (EVs). Nonetheless, the transport sector in 2016 still accounted for 24% of global emissions from fuel combustion, which is equivalent to more than 8 billion tons of emitted CO 2 [1]. In order to meet the goals of the Paris Declaration on Electromobility and Climate Change that set the requirements that at least 20% of all road transport vehicles globally need to be electrically driven by 2030, EVs play a major role in the future of sustainable transport systems, and the global proliferation of EVs is necessary [2].
In 2017, sales of new electric cars (battery EVs (BEVs), plug-in hybrid EVs (PHEVs), and fuel-cell EVs (FCEVs) exceeded 1 million vehicles, representing a growth in new electric car sales of 54% compared with 2016, and global stock surpassed 3 million vehicles [3].
The steady growth of EVs can, however, have adverse impact on distribution networks: increased peak demand and power losses, voltage violations, power-quality issues, and transformer heating and overloading. This impact was evaluated in multiple studies, for example, [4][5][6][7][8][9][10], and may vary on the basis of EV penetration level, battery characteristics, charging patterns, locations, modes, and times, battery state of charge (SOC) while charging, driving patterns, fleet charging profiles, driving distances, demand-response strategy, and tariffs [11].
The power-quality issue of harmonic distortion is mainly determined by the current-harmonic emissions of power electronic components and the network impedance of the grid. These current-harmonic emissions are normally quantified in terms of magnitude. However, phase angle is also of importance, as it allows for a better understanding of potential cancellation effects between devices and customers. Harmonic distortion is time-varying; current distortion varies with load levels, and harmonic voltages at the point of common coupling (PCC) aggravate the emission of harmonic currents. The authors in [12] stated that one can assume that harmonic currents are independent of bus voltages without much error. Other authors in [13,14] showed that, in some situations, when considering the sums of harmonics at different supply voltages, the variation can be significant. In [15], different voltage magnitudes and phase angles were applied while the harmonic currents of devices were measured. The results show that harmonic currents are mostly amplified by harmonic voltages.
In terms of EV charger models, there are two major approaches: circuit-based models in the time domain, and frequency-domain-based models. The first requires detailed a priori knowledge of the electrical equipment and its circuits, while the latter is based on excessive measurements that are matched, in an assumed or postulated model, with parameters fitted to the measured data. This work is based on the latter approach.
The first analysis of the effects of current harmonics generated by a cluster of EV battery chargers was theoretically done as early as the early 1980s. The authors concluded that consequences could be severe [4,16]. In those days, commercially available EVs were limited, and charging technology was in its early stages. In analysis made in 1996 [17], it was concluded that the improvement of EV battery-charging technology had affected the extent of emitted current harmonics, which were continuously decreasing. In [18], harmonic impact was investigated, and simulations showed that that the rise of total harmonic distortion (THD v ) would be very low. A method for predicting harmonic currents generated by a concentration of EVs was presented in [19]. The authors' analysis showed that the net produced harmonic current was not proportional to the number of EVs being charged, that is, the current-harmonics magnitude sum was far less than the total value of the current harmonics emitted by each charger. In [20], the integration of EVs in the Portuguese distribution grid was analyzed and showed that the simultaneous operation of multiple one-phase chargers in the same LV feeder could significantly raise current-harmonic levels. Simulations made in [21], where an EV charging station was modeled, showed that permissible current-harmonic limits were exceeded when six or more vehicles were simultaneously charging. In [22,23], it was shown that distortion limits would be exceeded when charging in relatively weak grids. The combination of four different EVs, and what their cancellation effect does to phase-angle diversity, was analyzed in [24]. Analysis was made during constant-charging mode, during which the phase angles were assumed to be constant. The results show significant harmonic cancellation. The same authors simulated a distribution network supplied by EV charging load, and the result indicate a minor increase in harmonic currents due to EV charging [13]. Kim et al. [25] showed that the low penetration and slow charging rate of EVs slightly impacted network harmonic distortion. However, the high penetration of EVs and fast charging rates may result in considerable voltage-and current-harmonic distortion [25,26]. Deilami et al. [27] established that randomly charging EVs could violate the standard level of voltage harmonics.
In some works, the interaction between EVs and other grid-connected equipment was examined. For example, in [28], the effect of combining EV charging and photovoltaics (PVs) was analyzed, and it was shown that the cancellation/amplification effect varied on the type of equipment and individual current harmonics. In [29], a harmonic decoupled power-flow model was developed to study the impact of multiple EVs with and without wind generators. The simulations demonstrated that low EV penetration had a significant impact on THD i (total harmonic distortion of current), but had low impact on THD v , while moderate-to-high EV penetration was shown to have an adverse impact on both THD i and THD v . With wind generation included in the simulations, THD could either be reduced or cause an increase; hence, the size and placement of wind generation is essential in order to manage THD levels. In [30], the interaction between EVs and other loads (foremost heat pumps) was analyzed. The study showed that the level of cancellation depends on the considered harmonic order, EV operating state, the phase angles of the background harmonic voltages, and the harmonic currents emitted by other loads.
In summation, the review of previous studies showed difficulties in predicting EV harmonic impact and the importance of further research within the field. The review also indicated that phase-angle fluctuation during charging is mostly neglected. In terms of hosting-capacity studies, there are no acceptable models for EV chargers [31].
The behavior of EV charging units is problematic for power systems, as each charger design differs from the next and usually contains some proprietary design. These design differences can result in different harmonic characteristics that can also fluctuate depending on the SOC. In other words, the magnitude and phase angle of each individual emitted current harmonic, at least in some of the measured vehicles and for some individual current harmonics, vary depending on vehicle and SOC. Since SOC is an important battery parameter that can be directly obtained and measured with high accuracy, its practical engineering value is great. It was shown that charging loads corresponding to different charging stages have significant difference; therefore, their impact on the power grid is also different [32]. In more recent works, Sheng et al. showed in their optimized charging model that the effect of a harmonic current can be reduced when EVs are charged at different initial SOC [33]. The authors in [34] analyzed the harmonic characteristics of multiple charging during asynchronous charging (charging of EVs with different initial SOC and starting at different points in time). They concluded that the harmonic current could be reduced by decentralizing EV charging start time and improving the charging state at the beginning of charging. When analyzing the effect on the distribution transformer, the authors in [35] concluded that the accuracy of harmonic assessment was improved if SOC was considered. In [36], the relationship between the amplitude and phase angle of each harmonic current and fundamental voltage, SOC, and the type of rectifier was nonlinearly mapped by the radial-basis-function network for an EV charging-station harmonic-modeling application. Similarly, the constructed EV charger models made in this work did not treat EVs as a constant power element, but also considered the SOC. Discrete SOC sets allow for aggregation analysis of multiple EVs at various states in the charging process. Even though harmonics can be reasonably statistic in the range of 20-80%, this is not always true, as presented in [37].
In this work, the magnitude and phase angles of current-harmonic emissions of four fully electric vehicles were measured several times during the mains part of the charging cycle. All calculations and analyses were based on nominal voltage magnitude (U n = 230 V); hence, on the assumption of nonvarying voltage supply. This was justified by the measured voltage levels and individual voltage harmonics. This is, of course, a simplification, since harmonics are affected by nearby loads and distortion in the supply voltage. It can, therefore, be argued that the modeling is only applicable in certain grid configurations. However, as an initial work that will be followed by a more exhaustive model (i.e., coupled Norton model) where these factors are considered, it is still justified because of the originality of the model construction.
Measurement results were checked for repeatability, i.e., maximal current-harmonic-magnitude deviation of 5%. On the basis of measurement data, power exponential functions were used to model individual current-harmonics behavior for each EV. In order to construct a continuous function with peaks and valleys where steepness and slopes were independent of each other, an analytically extended function was used that consisted of piecewise linear combinations of the power exponential function that was scaled and translated. Modeling the emitted harmonics as a function of SOC instead of using a look-up table presents new opportunities. One of these opportunities emerged as vehicle-to-grid (V2G) technology started to be implemented. This technology requires both new hardware (power electronics in the EV and the wallbox that enables bidirectional energy flow) and new software. The latter means that during the connection between the EV and complex wallbox, information, such as SOC information, is exchanged. When vehicles are connected, individual current-harmonic functions can be combined, and the aggregation can be made momentarily or be used to foresee aggregated effects to come. This method of harmonics modeling is a novelty in the research community. Technical standard ISO 15118, which handles the technical requirements on both hardware and software, is still under development.
Subsequently, the models were used to simulate the aggregation of multiple EV charging. This was followed by deterministic simulations to find situations (number of EV types, start time of charging, and level of SOC) when the maximal amplitude of individual current harmonics occurred.
The paper is organized as follows. Section 2 describes the theoretical background of harmonics and the aggregation of harmonic phasors. Section 3 describes the measurement setup. Section 4 describes the mathematical modeling, and Section 5 presents the results. In Section 6, results are discussed and conclusions are drawn.

Electric Vehicles and Harmonics
In a linear system, that is, a circuit with ideal resistors, capacitors, and/or inductors, the voltage and current waveforms are perfectly sinusoidal. In reality, such a system is highly unlikely, and there are always some nonlinearities. For example, certain types of loads (nonlinear loads) produce currents and voltages with frequencies that are whole multiples of the frequency at which the supply system was designed to operate (e.g., 50 or 60 Hz). These currents and voltages are called harmonics and range up to 2 kHz (40th harmonic).
A way to characterize the harmonic distortion of a signal is to use the so called THD that gives information on the harmonic content, and can be expressed as where h is the harmonic order, H is the highest order of the observed harmonic categorical, I h is the RMS value of the current h-th harmonic component, and I 1 is the current RMS value of the main frequency. THD v can be similarly calculated, but with the current values in Equation (1) replaced with respective voltage values. THD i does not always reflects the actual harmonic content. For example, if the fundamental current decreases while the harmonic content remains the same, the THD i value is higher, indicating higher emissions, but not reflecting the adverse impact on the grid. A more suitable way to express distortion in terms of current harmonics is the total harmonic current (THC): where I h is the RMS value of each harmonic order, and h is the harmonic order. In order to prevent harmonics from affecting the utility supply, standards such as IEC-61000-3-2 [38]/3-12 [39] and EN-50160 [40] were established with the goal of developing recommended practices and requirements for harmonic control in electrical power systems. IEC-61000-3-2 sets the emission limits for equipment connected to public low-voltage systems with input current ≤ 16 A per phase, and IEC-61000-3-12 for equipment with input current > 16 A and ≤ 75 A per phase. EN-50160 gives the main voltage characteristics at the customer's supply terminals in low-and medium-voltage networks; in terms of harmonics, it sets limits for harmonic voltage for orders up to the 24th. These standards concern both utilities and end users-the former because they are obliged to provide power of "good" quality, and the latter because they are responsible for not degrading voltage by drawing significant nonlinear or distorted currents. In this context, it becomes evident that power quality is of concern for both.
EV batteries require direct current (DC) for charging and power electronic converters for converting alternating current (AC) to DC, as well as charging control, presenting nonlinear load, and resulting in injected current harmonics. These harmonics are closely related to EV charger circuit topology that provides an interface with the AC network. While harmonic currents themselves do not directly affect customers, they can, as they are introduced into the distribution system, affect power losses in the entire low-voltage system, and cause nonlinear voltage drop that leads to distorted supply voltage, which can, in turn, affect electrical-distribution equipment (e.g., motors, transformers, lines, and capacitors) [41].
When a variety of vehicles with different SOCs and current variation during the charging cycle are connected to the electrical system, harmonic-frequency magnitude and individual harmonic phase-angle variation may be expected. The emissions of combined operations result in different cancellation and amplification effects, including mutual influence between voltage distortion and harmonic currents. As presented in the introduction, previous studies have not been univocal.

Harmonic Phasor Aggregation
Harmonic-aggregation studies are a complex topic, and research within the field is limited. For instance, up to now, little importance has been given to harmonic phase angles. The reason is the lack of standardization that has led to most instruments not providing that kind of information. Because of the difficulty in assessing phase angles from individual harmonics, most harmonic-aggregation studies used probabilistic approaches, for example, [9,12,[42][43][44][45][46].
An important aspect of realistic harmonic studies is an accurate representation of harmonic-current summations. When a number of loads are fed from the same bus, the total harmonic current injected at that bus is the sum of all individually generated harmonic currents. Since harmonic currents can be represented by vectors, the magnitudes and phase angles of all individual vectors are needed in order to make the vector summation. The diversification of devices with different circuit topologies can cause a diversity of current-harmonic phase angles, and may lead to a lower magnitude of vector sum than the arithmetical sum of harmonic currents [47].
The study of aggregated loads and the cancellation or amplification effects of harmonic currents needs the "absolute" harmonic phase angle; that is, the angle between current harmonics and fundamental voltage, as specified in IEC 61000-3-12 [39]. This differs from harmonic power-flow studies that require the "relative" harmonic phase angle (angle between harmonic voltage and harmonic current).

Measurement Setup
The power-quality (PQ) measurement equipment used in this work was a Yokogawa WT3004E precision power analyzer [48] compliant with IEC 61000-4-7 [49] measurement-instrumentation requirements that stores gapless 10/12-period values and each spectral component k at magnitude (rms) Y C,k and phase angle ϕ k . In a 50 Hz network, harmonic component Y H,h is identical to spectral component Y C,k , with k = h × 10, and is represented by harmonic phasor where Y is replaced by either I for current or U for voltage, and h indicates the harmonic order. The instrument has advanced computation installed that, in combination with the harmonic-/flicker-measurement software, allows us to perform harmonic measurements conforming to IEC 61000-3-2 and IEC 61000-3-12. This standard uses positive zero crossings as reference, thus calculating sine-based harmonic-phase angles (Figure 1). This differs from the Fourier definition, where cosine-based calculation is used.  Measurements were made on four fully electric vehicles in an environment within a local distribution grid; namely, the charging facilities of a distribution-system-operator, at 1 of the 8 outbound cables in an electrical enclosure; grid configuration can be seen in Figure 3. All vehicles were charged with one-phase charging, and, during measurements, only the charger included in the measurement was connected. Measurements were done several times and checked for repeatability. This feature was included in the software and compared, so that the measured magnitudes of the first 11 harmonics were within 5% deviation.
EVs indicate their SOC level in parts of full bar visualization on the instrument panel. The initial intention was to start charging and measure at SOC 25%. However, after measuring the amount of transferred energy from the grid to the vehicle, and comparing this to the battery capacity, it was concluded that the initial SOC varied from 20% to 36%.
Vehicles were chosen with respect to the most popular fully electric private vehicles in Sweden (Tesla models excluded). Technical and some measurement data of the EVs can be seen in Table 1. Due to sensitive issues, the names of these EVs are not disclosed.
The required measurement accuracy, referred to in the PQ instrument, is defined in IEC 61000-4-7 (exemplified in Table 2). The standard defines the required uncertainty for harmonic magnitudes, but excludes harmonic phase-angle uncertainty.
All vehicles were measured at least twice. As stated earlier, the measurements were checked for repeatability and examined so that the harmonic amplitude (up to the 11th harmonic) was lower than ±5%. Polar plots of the third current harmonic during the measurements of each vehicle can be seen in Figure 4. When measuring harmonics, the period of the fundamental signal (fundamental period) must be determined in order to analyze higher-order harmonics. The phase-locked loop (PLL) source is the signal that is used to determine the fundamental period. The data-update interval is automatically determined by the fundamental frequency of the PLL source and the number of periods of the PLL source used for analysis.
The measurement setup used in this work was a direct-input single-phase (two-wire) system. This setup can be seen in Figure 2.  Measurements were made on four fully electric vehicles in an environment within a local distribution grid; namely, the charging facilities of a distribution-system-operator, at 1 of the 8 outbound cables in an electrical enclosure; grid configuration can be seen in Figure 3. All vehicles were charged with one-phase charging, and, during measurements, only the charger included in the measurement was connected. Measurements were done several times and checked for repeatability. This feature was included in the software and compared, so that the measured magnitudes of the first 11 harmonics were within 5% deviation.
EVs indicate their SOC level in parts of full bar visualization on the instrument panel. The initial intention was to start charging and measure at SOC 25%. However, after measuring the amount of transferred energy from the grid to the vehicle, and comparing this to the battery capacity, it was concluded that the initial SOC varied from 20% to 36%.
Vehicles were chosen with respect to the most popular fully electric private vehicles in Sweden (Tesla models excluded). Technical and some measurement data of the EVs can be seen in Table 1. Due to sensitive issues, the names of these EVs are not disclosed.
The required measurement accuracy, referred to in the PQ instrument, is defined in IEC 61000-4-7 (exemplified in Table 2). The standard defines the required uncertainty for harmonic magnitudes, but excludes harmonic phase-angle uncertainty.
All vehicles were measured at least twice. As stated earlier, the measurements were checked for repeatability and examined so that the harmonic amplitude (up to the 11th harmonic) was lower than ±5%. Polar plots of the third current harmonic during the measurements of each vehicle can be seen Measurements were made on four fully electric vehicles in an environment within a local distribution grid; namely, the charging facilities of a distribution-system-operator, at 1 of the 8 outbound cables in an electrical enclosure; grid configuration can be seen in Figure 3. All vehicles were charged with one-phase charging, and, during measurements, only the charger included in the measurement was connected. Measurements were done several times and checked for repeatability. This feature was included in the software and compared, so that the measured magnitudes of the first 11 harmonics were within 5% deviation.
EVs indicate their SOC level in parts of full bar visualization on the instrument panel. The initial intention was to start charging and measure at SOC 25%. However, after measuring the amount of transferred energy from the grid to the vehicle, and comparing this to the battery capacity, it was concluded that the initial SOC varied from 20% to 36%.
Vehicles were chosen with respect to the most popular fully electric private vehicles in Sweden (Tesla models excluded). Technical and some measurement data of the EVs can be seen in Table 1. Due to sensitive issues, the names of these EVs are not disclosed. World Electric Vehicle Journal 2020, 11, x FOR PEER REVIEW 7 of 18  The required measurement accuracy, referred to in the PQ instrument, is defined in IEC 61000-4-7 (exemplified in Table 2). The standard defines the required uncertainty for harmonic magnitudes, but excludes harmonic phase-angle uncertainty. Table 2. Uncertainty requirements for current measurements (CLASS I instruments) [49]. I nom , nominal range of power-quality (PQ) instrument; I m , measured value.

Conditions Maximum Error (±ε)
Current I m ≥ 3% of I nom ±5% of I m I m < 3% of I nom ±0.15% of I nom All vehicles were measured at least twice. As stated earlier, the measurements were checked for repeatability and examined so that the harmonic amplitude (up to the 11th harmonic) was lower than ±5%. Polar plots of the third current harmonic during the measurements of each vehicle can be seen in Figure 4.    Table 2. Uncertainty requirements for current measurements (CLASS I instruments) [49]. Inom, nominal range of power-quality (PQ) instrument; Im, measured value.

Modeling Harmonic Currents Using Power Exponential Functions
In order to find a continuous function with multiple peaks, and where the steepness of the rise between peaks and the slope of the decaying parts are independent of each other, an analytically extended function was constructed. The function, with the same expression before and after maxima (or minima) time moments, but for different parameters, was used for approximating harmonic currents. This function consisted of a piecewise combination of the power exponential function that was scaled and translated, so that the resulting function was continuous.
First the power exponential function is, as described in [50], defined as

Modeling Harmonic Currents Using Power Exponential Functions
In order to find a continuous function with multiple peaks, and where the steepness of the rise between peaks and the slope of the decaying parts are independent of each other, an analytically extended function was constructed. The function, with the same expression before and after maxima (or minima) time moments, but for different parameters, was used for approximating harmonic currents. This function consisted of a piecewise combination of the power exponential function that was scaled and translated, so that the resulting function was continuous.
First the power exponential function is, as described in [50], defined as For non-negative values of t (time) and p (iterative obtained value for model fitting that could be negative), this function has a steeply rising part followed by a more slowly decaying part, exemplified in Figure 5. Then, a linear combination of a power exponential function, and a mirrored and translated power exponential function was constructed in order to obtain a function that best fit the measured values (see Appendix A for more details of the function): ( ; , , , , ) = ; Analysis showed that extending this function with more parameters pn marginally improves the model, but at the expense of computational effort.
where ( ) = | ( ) + ( )| = + , Since the harmonic curve can dramatically change shape at certain times due to an actual change to the charging circuit imposed by the control circuitry in the electric vehicle, a piecewise approximation of the harmonic currents using functions of the form g given in Equation (5) This approach is similar to the model used in [52,53] to describe electrostatic-discharge currents. The function's main advantages (as mentioned in [54]) are: simply adjustable derivative value, risetime value, time-to-peak value, and exact peak values chosen prior to adjusting other parameters. Changes in the current-harmonic curve could also reflect sudden changes on the supply side, e.g., Then, a linear combination of a power exponential function, and a mirrored and translated power exponential function was constructed in order to obtain a function that best fit the measured values (see Appendix A for more details of the function): Analysis showed that extending this function with more parameters p n marginally improves the model, but at the expense of computational effort.
In accordance with the IEC 61000-4-7 [49] standard, the following notation for frequency analysis of a current was used: where Since the harmonic curve can dramatically change shape at certain times due to an actual change to the charging circuit imposed by the control circuitry in the electric vehicle, a piecewise approximation of the harmonic currents using functions of the form g given in Equation (5) was used: 0, t < 0 g(t; t n , t n+1 , q m,n,1 , q m,n,2 , q m,n,3 ), t n ≤ t < t n+1 , 1 ≤ n ≤ N 0, t ≥ t N This approach is similar to the model used in [52,53] to describe electrostatic-discharge currents. The function's main advantages (as mentioned in [54]) are: simply adjustable derivative value, rise-time value, time-to-peak value, and exact peak values chosen prior to adjusting other parameters. Changes in the current-harmonic curve could also reflect sudden changes on the supply side, e.g., due to switching events. After examination of the voltage levels (and voltage harmonics), it was concluded that the curves in fact reflected the control circuit of the EVs reasonably well.
The first step in fitting the model is to choose intervals given by [t n , t n+1 ] in Expressions (7) and (8). This is done manually on the basis of the inspection of the measured data, that is, identifying peaks at certain times t. Expressions (7) and (8) show that values prior to t 1 are set to zero, that is, before the vehicle starts to charge. In many cases, choosing suitable intervals for one individual harmonic gives a set of suitable points for most or all other harmonics. The second step is then separately fitting the model in each interval. For this application, a least-squares fitting procedure was applied, and, since the model was nonlinear with respect to the parameters, a numerical procedure for finding the fitting was used. This was conducted with software [54] that uses an interior reflective algorithm based on [55,56].
The model construction can be summarized in the following steps: • inspection of measured value (imaginary and real part is done separately); • manually choosing moments in time when peaks occur; • choosing initial values for p 1 , p 2 , and p 3 (or q 1 , q 2 and q 3 ); • applying Function (7) or (8) between peaks; • least-squares fitting procedure (MATLAB) to find values of p 1 , p 2 and p 3 (or q 1 , q 2 and q 3 ) for a function that best fits measured values between peaks; • procedure is repeated until all peaks are handled. Functions between each peak are acquired, and a piecewise combination of the functions is scaled and translated, so that the resulting function is continuous.
Examples of functions over different time intervals are presented in Figure 6.
World Electric Vehicle Journal 2020, 11, x FOR PEER REVIEW 10 of 18 due to switching events. After examination of the voltage levels (and voltage harmonics), it was concluded that the curves in fact reflected the control circuit of the EVs reasonably well. The first step in fitting the model is to choose intervals given by , in Expressions (7) and (8). This is done manually on the basis of the inspection of the measured data, that is, identifying peaks at certain times t. Expressions (7) and (8) show that values prior to t1 are set to zero, that is, before the vehicle starts to charge. In many cases, choosing suitable intervals for one individual harmonic gives a set of suitable points for most or all other harmonics. The second step is then separately fitting the model in each interval. For this application, a least-squares fitting procedure was applied, and, since the model was nonlinear with respect to the parameters, a numerical procedure for finding the fitting was used. This was conducted with software [54] that uses an interior reflective algorithm based on [55,56].
The model construction can be summarized in the following steps: • inspection of measured value (imaginary and real part is done separately); • manually choosing moments in time when peaks occur; • choosing initial values for p1, p2, and p3 (or q1, q2 and q3); • applying Function (7) or (8) between peaks; • least-squares fitting procedure (MATLAB) to find values of p1, p2 and p3 (or q1, q2 and q3) for a function that best fits measured values between peaks; • procedure is repeated until all peaks are handled. Functions between each peak are acquired, and a piecewise combination of the functions is scaled and translated, so that the resulting function is continuous.
Examples of functions over different time intervals are presented in Figure 6. When performing the fitting, it is important to use initial parameters that at least qualitatively match the desired output; otherwise, the algorithm is likely to produce poor fitting. The result is graphically exemplified in Figure 7, where the real and imaginary part of the third current-harmonic during the charging of EV1 and the corresponding fitted function are presented in a complex plane. When performing the fitting, it is important to use initial parameters that at least qualitatively match the desired output; otherwise, the algorithm is likely to produce poor fitting. The result is graphically exemplified in Figure 7, where the real and imaginary part of the third current-harmonic during the charging of EV1 and the corresponding fitted function are presented in a complex plane.

Aggregation of Multiple Electric-Vehicle Loads
The simulation prerequisites are the following: we assumed that we had four vehicles of each model, and six charging points fed from the same bus. The numbers of charging points were based on an actual charging facility within a local distribution grid. The numbers of vehicles of each model were chosen so that at least two models would be simultaneously charged during the simulations. If the numbers of identical vehicles were equal to or more than the number of charging points, we likely had the largest current harmonics when identical vehicles were being charged, as the cancellation effect in this situation was limited. The measurement equipment stores gapless 10/12-period valuesamplitude and "absolute" phase angle for individual current harmonics. These values can be represented by real and imaginary values, as seen in Equation (6) The aggregated harmonic phasor can be expressed according to

Aggregation of Multiple Electric-Vehicle Loads
The simulation prerequisites are the following: we assumed that we had four vehicles of each model, and six charging points fed from the same bus. The numbers of charging points were based on an actual charging facility within a local distribution grid. The numbers of vehicles of each model were chosen so that at least two models would be simultaneously charged during the simulations. If the numbers of identical vehicles were equal to or more than the number of charging points, we likely had the largest current harmonics when identical vehicles were being charged, as the cancellation effect in this situation was limited. The measurement equipment stores gapless 10/12-period values-amplitude and "absolute" phase angle for individual current harmonics. These values can be represented by real and imaginary values, as seen in Equation (6), and the summation at each time instance can be made according to The aggregated harmonic phasor can be expressed according to The continuous SOC value of each model was calculated so that it corresponded to the time series of each model. For EV1, the time series started at SOC 20%; for EV2, the time series started at SOC 36%; for EV3, the time series started at SOC 32%; and for EV4, the time series started at SOC 25%.
The measurement series for all EVs and all individual current harmonics were randomly simulated (100,000 times) in terms of model, initial SOC, and charging start time (from time t = 1 min to t = 800 min). Lastly, a deterministic approach was applied to find circumstances when the maximal current-harmonic magnitude of individual harmonics occurred. The maximal current-harmonic magnitude of 6.89 A occurred when all four EV1 vehicles (and two EV2 vehicles) were being simultaneously charged. The polar plot visualizes the continuous function as charging started (from one end of the line at zero amps) until the charging stopped (at the other end of the line at zero amps). The phase angle during the entire time the vehicles were being charged was situated in the second and third quadrant. The measurement series for all EVs and all individual current harmonics were randomly simulated (100,000 times) in terms of model, initial SOC, and charging start time (from time t = 1 min to t = 800 min). Lastly, a deterministic approach was applied to find circumstances when the maximal current-harmonic magnitude of individual harmonics occurred. Figure 8 displays the result of the third current-harmonic simulation. The maximal currentharmonic magnitude of 6.89 A occurred when all four EV1 vehicles (and two EV2 vehicles) were being simultaneously charged. The polar plot visualizes the continuous function as charging started (from one end of the line at zero amps) until the charging stopped (at the other end of the line at zero amps). The phase angle during the entire time the vehicles were being charged was situated in the second and third quadrant. Table 3 shows the resulting situation for each individual current harmonic, the maximal magnitude and time when this occurred, the EV model being charged, the SOC level at the start, the time when the EVs started charging, and in which quadrant the phase angle was found.   Table 3 shows the resulting situation for each individual current harmonic, the maximal magnitude and time when this occurred, the EV model being charged, the SOC level at the start, the time when the EVs started charging, and in which quadrant the phase angle was found.

Discussion and Conclusions
This paper describes how power exponential functions can be used to model current harmonics emitted by the onboard charger, and thereafter simulate the aggregation of six fully electric vehicles.
On the basis of several measurements during the charging of each of the vehicles, fitted functions were created for all uneven individual current harmonics up to the 11th order. These fitted functions were then used to simulate aggregated EV loads fed from the same bus. The number of vehicles to be simultaneously charged was set to six, and the maximal quantity of each vehicle model was set to four. The SOC and start of charging were randomly selected, and each of the current harmonics were simulated 100,000 times. The simulation had a deterministic approach-to find the situation where the maximal current-harmonic magnitude could be found. Certainly, 100,000 simulations could not cover every single case, so the actual situation where the maximal magnitude occurred was not found. However, the results show an indication of the magnitude to be expected in which the quadrant the phase angle could be found, as well as the individual start and SOC of each of the vehicles. As expected, depending on the individual current-harmonic emissions, we mostly had all vehicles of the same model simultaneously charging because identical current-harmonic characteristics amplify the effects of aggregation and limit the cancellation effects. For all the harmonics except the ninth, a combination of the vehicles with the highest-rated charging power resulted in the highest harmonic magnitude.
The developed model is useful for determining the range of the current-distortion level of aggregated EV loads under different loading scenarios. Since the absolute phase angle varies during a charging cycle, cancellation and amplification vary too.
When technical standard ISO 15118 is implemented in both EVs and the wallbox, SOC information can be exchanged. If this information is collected, simulations can be continuously made, and the harmonic outcome can be predicted. This enables the more reliable operation and planning of the electrical grid.
Further improvements in the model can be achieved with load data from more EV models and measurements of a greater part of the battery's SOC. Most importantly, these simulations do not consider the surrounding environment, such as feeding voltage variations and other influencing nonlinear loads. Surrounding load characteristics, including prevailing phase angles and network impedance, affect the resulting influence of large-scale EV penetration. Modeling an entire network (e.g., using a harmonic coupled Norton-equivalent model) with included parameters largely improves the usability of the model. In future work, these parameters will be considered. Furthermore, the measurement equipment used in this work stored 10/12 period values, an unrealistic time resolution for permanent power-quality instruments located in the distribution grid. For this reason, future work will also be made on the impact of aggregation time intervals on harmonic phasors.