Methodology to Evaluate the Impact of Electric Vehicles on Electrical Networks Using Monte Carlo †

: In preparation for the electric mobility technological transition in Colombia, an impact assessment of the electric power system is required, considering the increasing loading that must be able to be managed in the future. In this paper, a plug-in electric vehicle (PEV) charging simulation methodology is developed in order to dimension the impact of this type of load on power grids. PEV electric properties, user charging behaviors, geographic location, trip distances, and other variables of interest are modeled from empirical or known probability distributions and later evaluated in different scenarios using Monte Carlo simulation and load ﬂow analysis. This methodology is later applied to the transmission network of Antioquia (a department of Colombia) resulting in load increases of up to 40% on transmission lines and 20% on transformers in a high PEV penetration scenario in 2030, increases that are well within the expected grid capacity for that year, avoiding the need for additional upgrades.


Introduction
In Colombia, the transportation sector is the largest energy consumer, consuming 39.8% of it and, at the same time, contributing 65% of the energy losses according to the 2015 Energy Balance. The consumed energy is taken mainly from diesel (37%) and gasoline (40%), with less than 1% of it being from electricity [1]. According to the International Energy Agency (IEA), sales of electric cars reached 2.1 million globally in 2019, accounting for 2.6% of global car sales that year and boosting the total plug-in electric vehicle (PEV) stock to 7.2 million. Colombia currently accounts for only 10,612 vehicles of that stock, according to the Colombian National Transportation Register (RUNT) [2], corresponding to a penetration level of 0.06%. In order to keep increasing those quantities, the EV30@30 campaign was launched at the 8th Clean Energy Ministerial meeting in 2017 with the endorsement of eleven countries and 29 companies, setting the goal of reaching a 30% market share for electric vehicles (EVs) of the total vehicle sales by 2030 [3].
The Colombian energy mix is known to be mainly based on hydropower [4], which being a renewable source, makes itself an attractive and clean power source for mobility. The future of urban transport systems will have an increasing portion of electric vehicles due to their sustainable and eco-friendly characteristics, especially when the energy used comes from renewable sources [5]. The EV penetration increase may have an impact on power systems' reliability, considering that the charging of a massive fleet of electric vehicles will increase the power load for the system itself and the local distribution network [6].
When looking at the state-of-the-art in power grids' asset response evaluation in the face of high PEV penetration, three main approaches can be found: static analysis, probabilistic analysis, and time series analysis.
Most static analyses consider that PEV charging processes start and stop always at the same hours of the day [7]. The studies carried out by Shafiee et al. [8] and by Pieltain Fernandez et al. [9] showed an increase in energy losses as one of the main concerns with widespread PEV usage. Other findings in these studies are related to a high peakto-average ratio (PAR) and high investment costs for distribution networks, both due to the coincidence between daily peak load and charging start times. The shortfall of the static approach is that only the worst-case scenario is considered; therefore, it tends to overestimate the impact of PEVs.
Improving on this approach, other studies build probabilistic models that may represent the charging load profile in a better way when compared to using deterministic charging patterns. Tehrani and Wang [10] proposed a probabilistic modeling approach, which makes use of statistical analysis, random simulation, and queuing theory in order to build the charging load distribution for PEVs. On the other hand, the simulation model proposed by Anand et al. [6] consists of two layers, where the first one represents the stochastic nature of the location and time of charging using traffic survey data as an input to a dynamic hidden Markov model (DHMM). The second layer adds the effects on the power distribution system, to later assess the impact of the PEV penetration on the distribution system's reliability performance using Monte Carlo simulation.
On the other hand, time series analyses make use of EV charging load profiles as inputs to a load flow analysis [7]. Some studies [11,12] make use of time series to measure the impact of PEVs under deterministic or stochastic configurations by simulating different PEV charging scenarios that take into account the charging start time, residual state-ofcharge (SOC), and final SOC, among other variables. In the study conducted in [11], one of the most interesting conclusions was that with a 20% PEV penetration, peak load would suffer a 35.8% increase in an uncontrolled scenario [7].
Other authors in Colombia have applied survey based methodologies and deterministic time series analysis to evaluate the impact of PEV massification. In some studies [13,14], the authors developed surveys to identify vehicle types and their characteristics, later estimating possible impacts on the Colombian electricity grid, while in others [15][16][17], the authors made use of a time series analysis with defined charging scenarios. However, a time series analysis with stochastic variables, using a method like Monte Carlo simulation, has yet to be implemented by any author in the local scope.
The contribution of this study is the application of Monte Carlo simulation as a statistical method to evaluate the impact of PEV penetration on the electricity grid of less developed countries, where currently there is not massive adoption of electric vehicles, and thus, real data from smart meters and EV chargers are not available. This approach uses origin-destination surveys, commonly available from governmental authorities, public EV sales data, and official PEV projection scenarios. Another novelty presented is the use of histogram data directly through quantile functions (inverse cumulative distribution method) and kernel density functions.
In this paper, a time series analysis with Monte Carlo simulation and charging scenarios is presented. This approach represents a novel study for Medellín city, and the presented methodology can be applied to other cities using location-specific origin-destination survey data and EV models. The content of this paper is as follows: In Section 2, a description of the methodology is shown. Section 2.2 describes the information required to apply the methodology. Section 3 shows the study case selected for this paper and its simulation results for the city of Medellín, and lastly, Section 5 resumes the conclusions reached and proposes identified related future work.

Methodology
In order to present the methodology, the input data and steps required to implement it are described in detail in this section.

Description
The proposed methodology can be seen in detail in Figure 1 and can be summarized as follows. First task is to process available stochastic input data, then variables like charging start time, daily trip frequency, and daily trips are classified by density distribution estimation methods like kernel density estimation (KDE) and inverse cumulative distributions functions (ICDFs). These methods allow generating random values for specific variables. The second task requires creating copies of EV models initialized with specific characteristics according to the selected simulation scenario. The third task is sampling random values from estimated distributions, modeling the charging process during a week with hourly time steps. Each EV is connected to a random busbar of the electrical network and is considered to behave like a constant power load from the moment of plugging in and until fully charged, and thus, smart charging strategies are not considered. In the fourth task, a load flow analysis is performed, where each hour is a different study case and the base load is affected by charging processes occurring at different points of the system. Monte Carlo simulation is achieved repeating Tasks 2, 3, and 4 as needed and can be executed for as many scenarios as needed.

Model For EV Charging
To set a model to determine the behavior of EV charging and estimate the electric power demand, it is necessary to define and evaluate the availability of the following information.

Daily Mileage
To estimate the daily mileage traveled by a PEV per day (DM, daily mileage), it is necessary to determinate the number of trips in a single day (DT, daily trips) and obtain the mean distance traveled per trip in the city of study (DpT, distance per trip). Once this information is obtained, the daily mileage is calculated according to Equation (1).

Charging Behavior
In order to model the charging process and frequency of charging, it is necessary to know the preferences of users related to PEV charging. There are three parameters to take into account: • Start charging time: A parameter that indicates the time of day when the user begins to charge his/her PEV. • State-of-charge (SOC): This is the state of the battery at a given moment, assuming that the battery is fully charged after ending the charging processes; the SOC when the charging processes ends is 100%, and it decreases as the vehicle is used. Equation (2) is adapted from [18] and describes how to calculate SOC from daily mileage (DM), EV performance (P EV ), and battery capacity (C).
• Charge preferences (user SOC): This parameter indicates the user preference to charge his/her PEV, that is the battery charge level when the user chooses to connect the PEV to the charging station. From these parameters it is possible to determine the number of charging hours for each user and the starting time of charge. Charging time (H) depends on residual battery power (SOC), charger efficiency (η), battery capacity (C), and the rated power capacity of the charger (P charger ), as can be seen in Equation (3), which is adapted from [18].

EV's Electric Energy Consumption
To quantify the total power required from the distribution grid, the number of PEVs connected to a particular busbar in each zone needs to be determined. Total demand per busbar at a specific hour (P busbar ) is equal to the power demanded from all loads different from the PEVs at that hour (P base ) plus the sum of the power demanded at the busbar by the electric vehicles at the same time (P PEV ) [19,20], as can be seen in Equation (4).

Monte Carlo Simulation
Monte Carlo simulation is a statistical method that uses random samples from probability distributions to model a phenomenon that contains stochastic elements. Due to the random nature of driving patterns, Monte Carlo simulation, as a stochastic modeling approach, is a common strategy used in electromobility research [21][22][23][24]. Monte Carlo simulation can be performed by following these steps [25]: • Model system variables as probability density functions (PDFs) or equivalent inverse cumulative distribution functions (ICDFs). • Repeatedly sample from the PDFs. • Compute the statistics and calculations of interest.

Study Case and Simulation Results
In order to test the methodology, a case study is carried out using an electrical network and mobility data corresponding to the city of Medellín, Antioquia, in Colombia for the year 2030.

Network Description
The power system considered in this study case has two voltage levels of 220 kV and 110 kV. The highest level network has a ring topology where generation units and network equivalents of larger regional systems are connected. For the current simulation objective and scope, external systems are supposed to deliver all power demanded from the base load and PEV charging process, only constrained by the specifications of 220 kV elements like lines and transformers. The 220 kV network is interconnected by different transformers to the inner 110 kV grid. This lower voltage grid is a mesh constituted by 20 substations. This electrical model is simulated using DIgSILENT Power Factory, representing an approximation of Medellín's distribution grid without including detailed distribution branches and loads. In Table 1, the main elements of the simulated electrical system are presented, and Figure 2 shows the single-line diagram of the network.
The input database is organized by study cases and operation scenarios. For each study case with a specific base load demand, there is an operation scenario that represents a single hour, starting with Monday for the period of 00:00 to 01:00 hours and ending Sunday for the period of 23:00 to 00:00. Base load demand in each study case represents projected consumption in the year 2030.

PEV Types and Characteristics
Defining the PEV model characteristics requires a series of steps prior to generating a simulation case with an important number of electric vehicles. The first step is to select which electric vehicles will be used to define the base models. The second step is obtaining the characteristics of those PEV models. The third and last step includes defining the base models, grouping the previously selected models by any criteria of interest, later averaging the specifications of each group, creating a unique set of specifications for each base model. Due to various power system structures and charging station characteristics, different charging levels are considered. In a real case, vehicles would have all kinds of different behaviors and power specifications, but still, by clustering vehicles into base groups, a time series simulation with stochastic distributions can be used to generate the scenarios of PEV fleet demand. For this study case, Colombian PEV sales data from 2016 to 2020 [26] were used to select the correct models. A total of 7 BEVs (BMW i3 60Ah, BMW i3 94Ah, Nissan Leaf, Mitsubishi i-MiEV, Renault Kangoo ZE, Renault Zoe, Kia Soul EV), 9 PHEVs (Mitsubishi Outlander PHEV, BMW X5 40e, Porsche Cayenne S E-Hybrid, BMW 330e, BMW 530e, BMW i8, Mini Cooper S E, Volvo XC60, BMW X5 45e), 2 taxis (BYD e5, BYD e6), and one bus (BYD K9GA) model were selected, and information about the battery capacity, performance, and charging power was collected or computed from each manufacturer's specifications; this information was later averaged for each PEV type. Table 2 shows the PEV types considered and their main characteristics.

PEV Penetration Scenarios
In order to conduct the assessment of PEV impact on the Colombian power system, three scenarios are considered regarding the PEV penetration level. These scenarios were constructed by Colombia's National Planning Department (DNP) aiming to fulfill the greenhouse emissions reduction, and they are based on government announcements regarding the penetration of 400,000 electric vehicles, equivalent to 5% of the vehicle fleet, in 2030 [1]. Considering the penetration goal of the countries belonging to the Electric Vehicles Initiative (EVI) is 30%, the government scenario is a low penetration scenario with a limited entrance of PEVs in the market. The medium scenario considers a faster increase of electric vehicles, so the penetration of this scenario is 10% of the vehicle fleet in 2030. The last scenario considers a high rate replacement by PEVs toward 2030, achieving a penetration of 15% of the vehicle fleet in 2030. It is estimated that the sum of private vehicles, taxis, and buses in 2030 will be 7,682,755 [1]. This study was simulated on a regional power distribution system, making it necessary to scale the number of vehicles in 2030 to a region or a city. Currently, the Antioquia department is the region with the largest amount of PEVs in Colombia, reaching 24.5% of Colombia's fleet [27]. The number of PEVs of each type was calculated using information from the DNP PEV penetration scenarios, while the distribution of private cars between BEVs and PHEVs was estimated making use of Colombian PEV sales data from 2016 to 2020 [26]. Table 3 presents the estimated number of PEVs in 2030 for the Antioquia department.

Daily Mileage
After processing the data from the metropolitan area of Medellín city [28], the dataset of daily trips and the histogram of the distance per trip were obtained. The probability density function of daily trips and distance per trip was generated from this source with the purpose of getting random values for daily mileage using Equation (1). For taxis, an average value of 148.9 km traveled per day was considered and 300 km traveled per day for buses [29].

Daily Trips
Kernel density estimation (KDE) was used to represent the probability density function of the number of daily trips using the dataset obtained from [28]. Figure 3 shows the probability density curve obtained. This approach considers every observation in the dataset represented as a normal distribution, and all of these normal distributions were added to calculate the final probability density function, as explained in [30]. The KDE algorithm from the Scikit Learn library [31] for Python was used with the Gaussian kernel and a bandwidth of 0.1. This value acts as a smoothing parameter, controlling the tradeoff between bias and variance in the estimation [32]. The KDE method is explained by general Equation (5).
where kernel K(x; h) is a positive function controlled by bandwidth h. K(x; h) ∝ e −x 2 2h 2 when the Gaussian kernel is selected.

Distance Per Trip
Data for the distance per trip are available through histograms from [28]. One method to use these resources for Monte Carlo simulation is using quantile functions (ICDF) [25]. For the ease of sampling purposes, ICDFs were used. In order to calculate the ICDFs, an equivalent histogram with bins of equal width is necessary. The ICDF can be defined as the inverse function of a CDF. Equation (6) shows the CDF for a discrete random variable, and the definition of the ICDF can be seen in Equation (7).
The ICDF allows knowing the probability that a number sampled randomly from the PDF will be less than or equal to x. Figure 4 shows the histogram from which the ICDF of distance per trip was calculated. This histogram is a modification of the original presented by the source, as it had unequally sized bins. The modification consisted of establishing a standard bin size of 1 km, so that narrower bins (short distances) were joined into a single one, while wider bins (long distances) were split into multiple bins of the same width and frequency.

Charging Behavior
In order to model the charging process and frequency of charging, it is necessary to know the preferences of the users for the start charging time and the preference of the users for choosing when to charge the vehicle's battery. Once these values are obtained, it is possible to determine the charging start time and the hours of charge needed to reach 100%.

Start Charging Time
The histograms shown in Figure 5a,b were used to obtain user start charging time preferences, which were used as inputs for the obtaining quantile functions (ICDF). Using the same method described before resulted in probability densities equivalent to the histograms. Buses in Medellín city follow a constant schedule for the charging process. The maximum end time for bus service is at 22 h, and a range between 21 and 23 h was assumed for start charging time [33].

Charge Preferences (user SOC)
There are two patterns to charge PEVs in accordance with the Electric Power Research Institute (EPRI)study in [34]. In Figure 6a, the histogram of the BEV users' charging preference related to the SOC level of the battery can be appreciated, and Figure 6b shows the histogram of the PHEV users' charging preference related to the SOC level of the battery. For public service vehicles (E-buses and E-taxis), it was assumed that they have a daily charging process, with the purpose of guaranteeing the availability of these to provide service and mobilize passengers.

Charging Schedule
To obtain the charging schedule of a PEV, the number of charging hours is calculated from Equation (3), considering charge preferences as the SOC of the equation. After charging time is calculated, the days of the week when the vehicle is charging must be determined. The daily SOC of the vehicle is calculated using Equation (2) and compared with the user SOC. If the SOC of the day is less than or equal to the user SOC, this day is marked as the charging day; otherwise, the remaining SOC of the day is saved for the next day. For days marked as charging days, the vehicle is charged for the number of hours calculated from the start charging time obtained previously until the charging process is complete.

Monte Carlo Simulation
In order to carry out a statistical approach to the behavior of PEVs in the city of Medellín for 2030, a Monte Carlo simulation was created under three penetration scenarios and using the deterministic parameters for the base vehicles described in Table 2 A total of 500 iterations were used in the simulation. In each iteration, the PEV quantities specified in the Table 3 were created according to the penetration scenarios mentioned before. Every PEV was connected to a system busbar with the corresponding charging schedule. The number of simulations can be determined according to the desired precision of the results.

Results' Analysis
In order to observe the voltage profile in the busbars and the load of the equipment due to the penetration of electric vehicles in the city of Medellín for the year 2030, load flows were simulated using Python scripts integrated with DIgSILENT Power Factory for all the iterations generated by the Monte Carlo simulation, first without the EV connected to the network (base case) and subsequently under the three penetration scenarios (see Table 3).
To analyze the results obtained, the average of the variables of interest resulting from every load flow simulated for each of the iterations of the Monte Carlo simulation was calculated.

Voltage Results
The voltage results show that the impact of PEVs on network voltages is negligible; the voltages in the base case are nearly equal to the voltages for the low, medium, and high penetration scenarios; in other words, voltages are invariant across all simulations. Figure 7a shows the five minimum voltages present in the busbars for all simulation cases, allowing one to appreciate that none of the busbar voltages drop under the lower voltage limit in Colombia (0.9 p.u.). On the other hand, Figure 7b shows the five maximum voltages present in the busbars for all simulation cases, where none of the busbar voltages exceeds the upper voltage limit in Colombia (1.1 p.u.).

Lines' Load Results
The most significant results of the lines' load are presented in Figure 8, and from this figure, a significant increase in load due to PEV penetration can be seen between 18 and 24 h. Figure 8a shows a line corresponding to a residential area where the impact of PEVs is very low. Figure 8b shows the load of a line corresponding to a mixed area (residential-commercial) where a large impact of PEVs can be seen, which implies a load increase of 10%, 30%, and 40% for the three simulated scenarios. Figure 8c shows the load of a line in an industrial zone where the load profile of the zone is flattened thanks to the penetration of PEV. Finally, Figure 8d shows the load of a line in a residential area where the increase in load in the last hours is very high, which could lead to load problems at some point.

Transformers' Load Results
The most significant results of the transformers' load are presented in Figure 9, and from this figure, a significant increase in load due to PEV penetration can be seen between 18 and 24 h. Figure 9a shows a transformer that mainly supplies industrial users where the impact of PEVs increases the load of the transformer between 18 and 24 h. Figure 9b shows a transformer that mainly supplies residential users where a large impact of PEVs can be seen, which implies a load increase of 4.8%, 12.9%, and 21.8% for the three simulated scenarios. Figure 9c shows the load of a transformer that supplies mixed users (residentialcommercial) where the load profile again presents an increase between 18 and 24 h. Finally, Figure 9d shows the load of a transformer that supplies commercial users where the load profile is flattened thanks to the penetration of PEVs.

Discussion
In the course of the development of the presented methodology, some limitations emerged, either from the design decision or external factors. The first limitation is related to the survey data, which are considered an approximation of real behavior and may be altered by poorly proposed questions and bad sample selection, among other reasons. The second limitation is the lack of information about the real location of charging stations, which prevents the methodology from estimating precisely the effect on the grid equipment.
Colombia has no infrastructure to develop smart charging, as there are not enough smart meters nor an hourly tariff scheme for households. Smart charging strategies require infrastructure investments and legislation changes in order to be viable and attractive to consumers. For these reasons, smart charging was not considered in the presented methodology; however, it could be later explored as a method to ease the adverse effects of unmanaged charging. Applying EV demand side management strategies [36][37][38][39] could enable users of shared charging stations to initiate the charging process taking in mind parameters like the battery SOC or charging power, even allowing charge authorization auctions in some stations. These strategies could avoid increasing the load in specific electrical network areas.
The proposed methodology was applied to the Antioquia department, while previous studies conducted in Colombia only considered the impacts in Bogota city. This involves different charging scenarios with random start charging times and connection busbars, while other studies in Colombia consider deterministic one-condition scenarios like office charging hours or household charging hours.
Other studies in Colombia have reached more alarming results for the grid impact of PEV charging, like the one carried out by Tellez Gutierrez et al. [16], in which one of the findings was that the distribution network would not be able to supply the power demand required by semi-fast charging systems if the EV introduction rates were higher than 8% of the total projected car fleet. This result differs from the one of the presented study case due to various reasons: the physical size of the studied electrical network is significantly lower; charging schedules and locations are static (causing concurrent charging); and a single charging station's power is considered. This serves also as an example of the difference in the results between static and stochastic/time series analysis.
Another research topic that can be explored is analyzing the effect in electrical networks of scenarios with bus rapid transit (BRT) systems [40] using in-motion charging (IMC) [41]. This could allow a reduction of the battery capacity required for buses, therefore causing the charging load to be distributed throughout the day at different points of the network.

Conclusions
This paper applies a Monte Carlo based methodology to assess the EV impact in a reallife-inspired electrical network, considering official Colombian government information such as penetration scenarios, reports of vehicles sold, user preference,s and the expected growth of the public and private vehicle fleets. The development of the methodology and the case study carried out allow reaching the following conclusions: • The quantile function (ICDF) is a useful tools to produce non-uniform random values in Monte Carlo simulations. These can be used as an alternative when using histograms (raw data are scarce) and regional probability density functions to model stochastic variables required in EV simulations, like mileage driven and start charging time, are not known. • In a high PEV penetration scenario, a load increase of up to 40% in lines between Hours 18 and 24 is expected. Lower penetration scenarios may carry load increases between 10% and 30% in lines. It is necessary to promote the creation of public access charging stations in commercial areas, universities, and companies, in order to encourage different schedules of vehicle charging and to try to balance loadability and avoid bottlenecks in networks. • Due to the stochastic nature of the location of residential EV charging stations, it is necessary to encourage the charging of EVs in different areas of Colombian cities and at different times of the day, in order to prevent distribution system equipment overloading in certain areas of a city at load peak hours. In a high PEV penetration scenario, the maximum load increase of the transformers is close to 22%; however, this equipment does not reach an overload state. • It is necessary to classify the total load into residential, industrial, and commercial categories in order to correctly evaluate the impact of EV charging on different distribution systems. It is possible that if the correct category is not assigned to an area, the results would not show problems that may arise due to the nature of the load. • The proposed methodology can be applied to transmission, as well as distribution networks. A model that represents the PEV charging station distribution across the network in a deterministic way could allow measuring more precisely the impact on the electrical networks. The application of this kind of charging model on distribution networks, where the impact of PEVs is expected to be higher, could allow performing a more successful system planning analysis. • The presented methodology was applied to the Antioquia transmission network. According to the results in the high PEV penetration scenario, the capacity of the grid in 2030 can withstand the effect of EV charging and does not need an upgrade. • Electrical systems with electrical urban rail transit systems, like the study case, can make use of the existing traction substations to feed future charging points for buses and taxis at valley hours, allowing the use of the existing robust infrastructure without increasing the load on other distribution networks.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the data use agreement.

Acknowledgments:
The authors would like to thank the support of the project: "Estrategia de transformación del sector energético colombiano en el horizonte de 2030" funded by the research call 778 of MinCiencias Ecosistema Cientifico, Contract FP44842-210-2018, and the research group of Transmission and Distribution of Electric Power (T&D) at Pontificia Bolivariana University. Thanks are also given to professor Ernesto Pérez from National University and Solenium for providing key information to build the study case.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: