Optimal Charging of Plug-In Electric Vehicle: Considering Travel Behavior Uncertainties and Battery Degradation

The negative environmental impacts of using fossil fuel-powered vehicles underlined the need for inventing an alternative eco-friendly transportation fleet. Plug-in electrical vehicles (PEVs) are introduced to cut the continuing increase in energy use and carbon emission of the urban mobility. However, the increased demand for mobility, and therefore energy, can create constraints on the power network which can reduce the benefits of electrification as a certain and reliable source. Thus, the rise in the use of electric vehicles needs electric grids to be able to feed the increased energy demand while the current infrastructure supports it. In this paper, we introduce a methodological framework for scheduling smart PEVs charging by considering the uncertainties and battery degradation. This framework includes an economic model for charging and discharging of PEVs which has been implemented in a 21-node sample distribution network with a wind turbine as a distributed generation (DG) unit. Our proposed approach indicates that the optimal charging of the PEVs has a high impact on the distribution network operation, particularly under the high market penetration of PEVs. Thus, the smart grid to vehicle (G2V) charging mode is a potential solution to maximize the PEV’s owner profit, while considering the battery degradation cost of the PEVs. The simulation result indicates that smart charging effectuation is economical.


Introduction
The rapid growth in urban mobility and the corresponding increase in energy use, greenhouse gas emissions, (GHG), air quality [1], and their externalities have driven more people to switch from their conventional cars to cleaner fleets, electric cars. The population of plug-in electric vehicles (PEVs) increased drastically over the past decade [2]. The PEVs are equipped with a large battery bank and the need for fast charging of these batteries and range anxiety issues impose an enormous load demand on the residential and distribution networks and this made many operation problems for the distribution system operator.
Two approaches might be considered when the charging infrastructure is needed. First is developing a charging infrastructure and increasing the capacity of the network which requires further investment in grid peak capacity. The second approach is managing the demand by considering the driving behaviors and charging schedules for cars and predicting the load pattern that the PEVs force on a daily basis. In this approach, if the pricing and rate structure are correctly planned, the demand for additional peak capacity investment will be declined and the cost of PEVs charging will be reduced, in one hand. It also improves the stability of the grid, on the other hand. Grid-to-vehicle (G2V), and vehicle-to-grid (V2G) are the two solutions suggested to avoid the loss of the power grid by supplying the network demand in a peak load period. G2V is the common type of PEV charging procedure, where power goes from grid to vehicle, but V2G, which has various applications in the power system auxiliary services, transfers power from vehicle to grid. In G2V and V2G, the aggregator can control and operate a group of PEVs and these vehicles act like storage systems.
There is a body of literature that proposed various managed PEV charging methods [3][4][5][6][7]. Prior studies in the area of PEVs and electric grid are based on the models, since the actual data of PEVs in a large scale are limited. However, the results of some of the studies are hampered because of much uncertainty surrounding the drivers' behaviors. For example, the current trend in the literature is to use a simplified PEV travel profile by mapping the travel behavior of conventional vehicles, which might introduce bias into the results. In addition, considering only a single PEV fleet in the simulation as a proof of concept or assuming that all vehicles have the same features are other shortcomings of prior studies [8][9][10][11][12]. The limitations of the distribution network are essential to be considered in the modeling framework, while they are neglected in some prior studies [13][14][15]. On-peak charging might cause serious problems such as voltage instability, reserve capacity, and voltage collapse, voltage insecurity, and power loss. In addition, while dynamic electricity price is critical in a smart charging system, the availability of accurate electricity price is an area of concern in a realistic optimization problem [16][17][18]. Dynamic pricing can encourage PEV owners to charge their vehicles at a lower rate.
In this study, we compare the uncoordinated, smart G2V versus V2G charging mode for a set of PEV fleets in a sample distribution network to understand which charging mode has supremacy over the others. We integrated wind energy into the power grid to consider using renewable energy resources in the grid. To fill the existing gap in the literature, the uncertainties of load, wind speed, electricity price, and travel behavior are studied in the PEV-grid model. A detailed driving pattern including the arrival and departure time and travel distance of PEVs drivers are taken into account. The Monte Carlo simulation with scenario reduction task is also conducted to model the stochastic nature of PEVs' travel behavior. The probabilistic economic dispatch is also applied to optimize the charging/discharging problem. This study provides insight into the effect of battery degradation in various charging modes (V2G and G2V) and presents the results in an economic and technical view.
The paper is structured as follows: Section 2 includes the problem statement and the method. In Section 3, the numerical results of this study are given in different approaches, which clarify the performance of the proposed method. Section 4 provides the numerical results of different charging procedures. Finally, Section 5 concludes the findings of this study.

Problem Statement
In this research, the objective function in the optimization of PEVs smart charging and discharging is to minimize the cost function subject to the network constraint: where Ra s (t) is the active and Rr s (t) is the reactive part of the power supplied in the scenario s at time t. Indices s and t are presenting scenario number and time respectively, where s ranges from 1 to 10 to consider all 10 scenarios in this study, and t ranges from 1 to 24, to include all hours in a single day. OCa s (t) and OCr s (t) are the price for active and reactive power, and DEC s (t) is PEVs fleet total battery degradation cost in the scenario s, respectively. The probability of each scenario is presented by Pr(s). The load and generation equation can be written as follows [16]: where Wa s (t) is the active power that the wind turbines produce, La s (t) is the real part of load demand, and PEVa s (t) is the absorbed or injected active power by PEVs in time t and scenario s. Wr s (t) is the reactive power produced by the wind turbines, Lr s (t) is the reactive part of the demand and PEVr s (t) is the absorbed reactive power that is consumed by the charging equipment. Plossa s (t) and Plossr s (t) are power loss in a transmission line that can be calculated using Equations (4) and (5) [19]: where i and j are indices of buses and n is the number of buses. G and B are the real and imaginary part of the admittance matrix, and θ is the voltage angle in this equation. V s (i, j, t) is the voltage difference between buses i and j at time t in scenario s.
In this framework, power and constraint of all nodes such as a voltage of nodes and frequency of system should not surpass the boundaries limitation. The nodes' constraints are introduced in Equations (6) and (7) [20].
. . , n eq , n eq = total number o f buses . . , n LD t = 1, 2, . . . , 24 n LD = total number o f load nodes P u (t) and Q u (t) represent active and reactive consumed power by the equipment u, which must not pass the equipment's maximum active and reactive power that are presented by P u max and Q u max . Voltage magnitude of all nodes V l (t) must stay between the minimum and maximum value of them, which are presented by V min , V max . The most significant restraint in PEVs charging and discharging problem is the state of charge (SOC) of the batteries which should fulfill the minimum boundaries at the moment of departure. This constraint is presented in Equation (8) [19]: In this study, for customers convenience, it is assumed that the SOC dep is at least 80 percent of the total capacity [19]. Here, SOC p,s (t dep ) shows SOC of p-th PEV battery at the departure time (t dep ) in scenario s. The undermentioned constraint must be satisfied for all PEVs' battery.
where, P chr p (t), P dichr p (t) are charging and discharging rate of each PEVs in time interval t, and η chr , η dichr are the total efficiency of battery charging/discharging and converter, respectively [16].
In Equation (10), battery degradation cost has been calculated based on the PEVs total battery degradation at each time for each scenario [21].
where C Batt p is the p-th PEV battery cost and D TOT p,s represents summation of calendar (D CAL p,s ) and cyclic (D CYC p,s ) degradation of PEVs' battery in the cycle of T cyc c,s in scenario s and cycle c. Also, Q 0 p , Q Use p stand for the nominal and useful capacity of p-th PEV, respectively. Q Use p is extracted based on Equation (11).
Stem from the United States Advanced Battery Consortium (USABC) PEVs' battery should be changed if the battery capacity is reduced to 80% of its total capacity [21], thus ρ = 0.8.
Equations (13) and (14) illustrate the mechanism of the calculating calendar (D CAL p,s ) and cyclic (D CYC p,s ) degradation [21].
To model the degradation of the PEVs' batteries, two main criteria including calendar degradation (D CAL p,s ) and cycling degradation are considered. Employing calendar degradation (D CAL p,s ), which is a function of the SOC of p-th PEV over period T (SOC p,s (T cyc c,s )), nominal SOC of p-th PEV (SOC), ambient temperature (ψ B ), nominal temperature (ψ 0 ), and charging time (t s ), is a vital part in battery degradation task. Also, a and b are calendar degradation's fitting parameter. The cycling degradation, which is a function of DOD over period T (DOD p,s (T cyc c,s )), has a significant effect on the total battery life. a and b in Equation (13) are fitting parameters for calendar degradation. α 1 , α 2 , and α 3 are fitting parameters for cycling degradation related to DOD. In Equation (15), the benefits of PEV aggregator has been introduced.
where the welfare is calculated as the difference between power consumed by charging and discharging of PEVs' batteries in two charging modes: uncoordinated and smart charging [21].
The power loss in the transmission line in the smart charging mode is different from the uncoordinated mode, and the difference between this two is considered as the network benefit.
The total benefit of using smart charging and discharging of PEV batteries in smart G2V and V2G mode is calculated using Equation (17). (17) where the C inf is installing smart charging and equipment cost that is extracted from [16]. Figure 1 shows the methodological framework of the study. In general, the driving pattern of PEVs drivers, the effects of wind energy, electricity price, and load demand are considered in the simulation. Then the probabilistic economic dispatch is applied to optimize the charging/discharging problem.

Modeling Load Demand, Price, and Wind Speed Uncertainties
The uncertainties in the load demand, electricity price, and wind speed are simulated by using a Monte-Carlo simulation which produces 1000 scenarios of possible probabilities [24]. The normal and Weibull distribution functions are suitable solutions for modeling system objective uncertainties [25].
The input data for the Monte-Carlo simulation is extracted from the PDF created by normal distribution function for the load demand and electricity price. The PDF created by the Weibull distribution function for the wind speed is also extracted for the Monte-Carlo simulation [24]. Normal and Weibull distribution functions are presented in Equations (19) and (20).
The generated scenarios for load demand, price, and wind are shown in Figures 5-7, respectively. The relationship between wind speed and wind power is derived from [26] which is shown in Equation (21). This function transforms the wind speed to electricity power for a wind turbine. It is important to specify the initial amount of PEVs' SOCs. Monte-Carlo approach is conducted to model the PEVs' travel behaviors under uncertainties [22]. Monte-Carlo is applied to produce different scenarios based on the input variables such as arrival and departure time, initial SOCs derived from travel distance, and the type of vehicle in order to get battery characteristics. The initial SOC of each PEVs is derived from Equation (18).
where, d n p is the traveled distance of p-th PEV in q-th iteration by Monte-Carlo simulation, C e f f (km/kWh) presents the efficiency coefficient of the battery during the driving mode, and Cap bat,p (kWh) is the capacity of p-th PEV's battery. The parameters of this equation are selected according to the PEVs model. C e f f and Cap bat,p are obtained from [16,19,23]. The input data are extracted from the probability distribution function (PDF) of travel data to produce the necessary information for PEVs fleet. The distribution function is used to create the PDF of arrival time and travel distance is the generalized extreme value (GEV) distribution which are shown in Figures 2-4. Also, for departure time Weibull distribution function is taken into account as shown in Figure 3.

Scenario Reduction
To handle the uncertainty of various stochastic parameters such as electricity price, wind speed, and load demand, the Monte-Carlo simulation task is employed in this study. Monte-Carlo simulation generates various scenarios based on the probability distribution functions. In this study, to increase the accuracy of the forecasting results, 1000 scenarios are generated in the first step. This large volume of data makes high computation cost because calculating economic dispatch for all 1000 scenarios is time-consuming. To handle this problem a scenario reduction task is employed based on the forwarding scenario selection and the backward scenario reduction approaches that are elaborated in [20]. This method eliminates the similar scenarios to 10 final scenarios which have the Density Figure 4. PDF of travel distance.

Modeling Load Demand, Price, and Wind Speed Uncertainties
The uncertainties in the load demand, electricity price, and wind speed are simulated by using a Monte-Carlo simulation which produces 1000 scenarios of possible probabilities [24]. The normal and Weibull distribution functions are suitable solutions for modeling system objective uncertainties [25].
The input data for the Monte-Carlo simulation is extracted from the PDF created by normal distribution function for the load demand and electricity price. The PDF created by the Weibull distribution function for the wind speed is also extracted for the Monte-Carlo simulation [24]. Normal and Weibull distribution functions are presented in Equations (19) and (20).
The generated scenarios for load demand, price, and wind are shown in Figures 5-7, respectively. The relationship between wind speed and wind power is derived from [26] which is shown in Equation (21). This function transforms the wind speed to electricity power for a wind turbine.
where P wd n w ,t is the power generated by the wind turbine. In general, generated power equals to zero when the wind speed is lower than the cut-in speed (v ci ) or is greater than the cut-out speed (v co ). When the wind speed is between rated and cut-out speed, the turbine produces its nominal power; also if it is between a cut-in and rated speed (v r ), the generated power is calculated by related function that is shown in Equation (21). Table 1 shows the parameters of the distribution functions for the arrival and departure time and travel distance.

Scenario Reduction
To handle the uncertainty of various stochastic parameters such as electricity price, wind speed, and load demand, the Monte-Carlo simulation task is employed in this study. Monte-Carlo simulation generates various scenarios based on the probability distribution functions. In this study, to increase the accuracy of the forecasting results, 1000 scenarios are generated in the first step. This large volume of data makes high computation cost because calculating economic dispatch for all 1000 scenarios is time-consuming. To handle this problem a scenario reduction task is employed based on the forwarding scenario selection and the backward scenario reduction approaches that are elaborated in [20]. This method eliminates the similar scenarios to 10 final scenarios which have the most probability. This procedure is applied for handling the system uncertainties. For PEV's modeling, 1000 scenarios which are generated by the Monte-Carlo simulation are reduced to 420 scenarios where 30 PEVs are considered in each bus in distribution network. Then, the economic dispatch is calculated for these scenarios and the vehicle owners' savings are derived by multiplying the optimized cost function by their relative probability.
The boxplot in Figure  5 shows the load profile of the modeling scenarios after the scenario reduction. Each whisker shows the variability of the load profile of scenarios for each hour of the day. On average, the load profile is higher between 12 pm to 5 pm. Electricity price is at the highest point at 5 pm, when the load profile is also the highest, on average, as shown in Figure 6. Figure 7 shows the scenarios for six wind turbines which generate active and reactive power in the distribution network.

Optimization Procedure
A mixed-integer non-linear problem (MINLP) is presented to optimize the cost function by considering the PEVs batteries' SOC at the departure time as the optimization constraints. MINLP is considered because of the complexity of the economic dispatch, power loss considerations and network constraints compounded with the charging and discharging of PEVs' batteries. The main goal is to minimize the Equation (1) by considering the other mentioned constraints such as SOC values.
In the optimization procedure, the maximum value of the PEV's SOC is not considered as a constraint for the problem formulation. The minimum value of batteries' SOC is employed as a constraint, so the operator can charge the batteries as much as possible if this procedure is beneficial (i.e., depends on the electricity price and power system load). For solving the MINLP problem, General Algebraic Modelling System (GAMS) software has been used [27]. The "Conopt-4" solver is employed which has the capability to solve the convex nonlinear problem.

Results
For the numerical study, a 21-node sample distribution network has been employed which is 13.8 Kv and has 20 feeder section which is shown in Figure 8. In this sample, six wind turbines which are able to produce 500 kVA are set up in buses 14, 15, and 18 to 21. The network characteristics and constraints such as voltage deviation and HV/MV substation capacity are derived from [16]. The load demand, electricity price, and wind speed data are simulated by the Monte-Carlo simulation. The normal distribution function for the load demand and electricity price and the Weibull distribution

Scenario Reduction
To handle the uncertainty of various stochastic parameters such as electricity price, wind speed, and load demand, the Monte-Carlo simulation task is employed in this study. Monte-Carlo simulation generates various scenarios based on the probability distribution functions. In this study, to increase the accuracy of the forecasting results, 1000 scenarios are generated in the first step. This large volume of data makes high computation cost because calculating economic dispatch for all 1000 scenarios is time-consuming. To handle this problem a scenario reduction task is employed based on the forwarding scenario selection and the backward scenario reduction approaches that are elaborated in [20]. This method eliminates the similar scenarios to 10 final scenarios which have the most probability. This procedure is applied for handling the system uncertainties. For PEV's modeling, 1000 scenarios which are generated by the Monte-Carlo simulation are reduced to 420 scenarios where 30 PEVs are considered in each bus in distribution network. Then, the economic dispatch is calculated for these scenarios and the vehicle owners' savings are derived by multiplying the optimized cost function by their relative probability.
The boxplot in Figure 5 shows the load profile of the modeling scenarios after the scenario reduction. Each whisker shows the variability of the load profile of scenarios for each hour of the day. On average, the load profile is higher between 12 pm to 5 pm. Electricity price is at the highest point at 5 pm, when the load profile is also the highest, on average, as shown in Figure 6. Figure 7 shows the scenarios for six wind turbines which generate active and reactive power in the distribution network.

Optimization Procedure
A mixed-integer non-linear problem (MINLP) is presented to optimize the cost function by considering the PEVs batteries' SOC at the departure time as the optimization constraints. MINLP is considered because of the complexity of the economic dispatch, power loss considerations and network constraints compounded with the charging and discharging of PEVs' batteries. The main goal is to minimize the Equation (1) by considering the other mentioned constraints such as SOC values.
In the optimization procedure, the maximum value of the PEV's SOC is not considered as a constraint for the problem formulation. The minimum value of batteries' SOC is employed as a constraint, so the operator can charge the batteries as much as possible if this procedure is beneficial (i.e., depends on the electricity price and power system load). For solving the MINLP problem, General Algebraic Modelling System (GAMS) software has been used [27]. The "Conopt-4" solver is employed which has the capability to solve the convex nonlinear problem.

Results
For the numerical study, a 21-node sample distribution network has been employed which is 13.8 Kv and has 20 feeder section which is shown in Figure 8. In this sample, six wind turbines which are able to produce 500 kVA are set up in buses 14, 15, and 18 to 21. The network characteristics and constraints such as voltage deviation and HV/MV substation capacity are derived from [16]. The load demand, electricity price, and wind speed data are simulated by the Monte-Carlo simulation. The normal distribution function for the load demand and electricity price and the Weibull distribution function for the wind speed data are taken into account. The load, electricity price, and wind speed data are for Ontario, Canada [28]. The parameters of the sample network are shown in Table 2. It is assumed that 30 PEVs are connected to each demand node. The information on arrival and departure time, travel distances, and type of vehicles are based on the NHTS database [29]. To show the performance of the proposed charging method, three different approaches including uncoordinated, smart G2V, and smart V2G are considered.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 10 of 17 function for the wind speed data are taken into account. The load, electricity price, and wind speed data are for Ontario, Canada [28]. The parameters of the sample network are shown in Table 2. It is assumed that 30 PEVs are connected to each demand node. The information on arrival and departure time, travel distances, and type of vehicles are based on the NHTS database [29]. To show the performance of the proposed charging method, three different approaches including uncoordinated, smart G2V, and smart V2G are considered.

Uncoordinated Charging
In the uncoordinated charging approach, PEVs drivers plug in their vehicle as soon as they arrive home. The cost of charging is equal to the summation of the consumed power by the PEVs battery to charge the vehicle. There is a considerable peak in the total load demand of charging PEVs because the residual load demand and PEVs load demand occur semi simultaneously. This peak intensifies the importance of using smart charging during a day. In uncoordinated charging during the night, PEVs charging rate is close to zero. The costs and power loss of the uncoordinated charging approach are shown in Table 3. Table 2. Parameters of the Sample Network.

Smart G2V Charging
In the G2V approach of charging PEVs, the aggregator can control the rate of power and time of charging. In this way, the amount of money that PEV owners should pay is much lower than the uncoordinated charging modes since the cars can be charged at the off-peak time. Electricity price is negative in several hours a day in Ontario, so it is beneficial for PEVs owners and operators to charge the PEVs at this time. PEV owners will be paid by the network operator due to generation constraints. As Table 4 shows, purchased energy cost in the smart G2V mode is lower than the uncoordinated charging mode, likewise the degradation cost. The reduction in degradation cost is caused by spreading the charging time interval.

Smart V2G Charging
In smart V2G charging, the PEVs can operate in both vehicles to the grid (V2G) and grid to vehicle (G2V) modes. The charging rate in every hour is varied to minimize the operation cost and maximize the SOCs of the PEVs. The results are illustrated in Table 5. At this point, the same as G2V charging mode, the purchased energy cost is reduced in comparison with the uncoordinated and smart G2V charging mode. Degradation cost in V2G mode is higher than the G2V charging approach.

Discussion
The benefits of various smart charging mode in comparison with uncoordinated mode are shown in Table 6. The PEVs battery degradation cost in smart G2V charging mode is higher than the one in V2G charging mode. This difference exists because in V2G charging mode, PEVs are used as the electricity storage. Using V2G charging mode reduces the purchased energy cost in comparison with smart G2V charging mode, while the degradation cost is higher than G2V. It is because batteries are always in interaction with the network by trading their power, in V2G mode. The difference between the degradation cost in three charging approaches is shown in Figure 9. As illustrated by the numerical results, the proposed method has a significant effect on the total charging cost of the PEVs' owners. Thus, implementing smart charging is recommended with high penetration of PEVs.

Conclusions
In this study, we compared three ways of PEVs charging including uncoordinated charging, smart charging, and V2G charging in a sample distribution network with/without considering the wind generation. The findings suggest that considering degradation cost of PEV's batteries is a vital factor which should be calculated for different charging approaches. While it is almost ignored in common modeling frameworks, modeling the battery degradation of PEV's and its non-linear behavior needs to be further studied. In this study, we showed how considering the behaviors of the battery can change the results. An algorithm based on Conopt-4 solver in GAMS was used to minimize the operation and battery degradation cost. Based on the optimization results, a G2V method is much beneficial than two other methods in terms of owner saving. However, the smart The simulation task is done by a PC with an Intel Core i7, 3.4 GHz CPU, and 32 GB DDR4 RAM, and computational cost for different approaches is presented in Table 7.

Conclusions
In this study, we compared three ways of PEVs charging including uncoordinated charging, smart charging, and V2G charging in a sample distribution network with/without considering the wind generation. The findings suggest that considering degradation cost of PEV's batteries is a vital factor which should be calculated for different charging approaches. While it is almost ignored in common modeling frameworks, modeling the battery degradation of PEV's and its non-linear behavior needs to be further studied. In this study, we showed how considering the behaviors of the battery can change the results. An algorithm based on Conopt-4 solver in GAMS was used to minimize the operation and battery degradation cost. Based on the optimization results, a G2V method is much beneficial than two other methods in terms of owner saving. However, the smart V2G method is more efficient than the uncoordinated method. The results show that in smart G2V charging mode, the net profit is equal to 287.85 ($/day) in the presence of wind power while in V2G charging mode, the net profit is much less than smart G2V Charging. The huge difference between the charging cost in the uncoordinated method and the two other methods is because of the time of charging. Charging in the off-peak time is much cheaper than the peak time. In addition, the PEV owner can save more if the vehicle is charged during the night when the energy price is negative.   Normal distribution function f (x|λ, κ)

Conflicts of
Weibull distribution function La s (t) Active load at time t in scenario s Lr s (t) Reactive load at time t in scenario s Pr(s) Probability of the scenario s P u (t) Active power of equipment u at time t P chr p (t) Charging rate of p-th PEV ad time t P dichr p (t) Discharging rate of p-th PEV at time t PEVa s (t) PEV's active power at time t in scenario s PEVr s (t) PEV's reactive power at time t in scenario s Plossa s (t) Active power loss at time t in scenario s Plossr s (t) Reactive power loss at time t in scenario s P wd n w ,t Power generated by the wind turbine P Wmax Nominal generated power by the wind turbine Q u (t) Rective power of equipment u at time t in scenario s Q USE p Useful capacity of p-th PEV Q 0 p Nominal capacity of p-th PEV Ra s (t) Purchesed active power at time t in scenario s Rr s (t) Purchesed reactive power at time t in scenario s SOC n init,p SOC at initial time t for p-th PEV SOC p,s (t dep ) SOC at departure time t for p-th PEV Wind active power at time t in scenario s Wr s (t) Wind reactive power at time t in scenario s V s (i, j, t) Line voltage between bus i and j at time t in scenario s V l (t) Voltage value for bus p at time t v t Wind speed v ci Wind turbine cut-in speed v co Wind turbine cut-out speed v r Wind turbine rated speed θ s,i,t Voltage's angel of node p at time t θ s, j,t Voltage's angel of node q at time t D.