Joint Optimal Allocation of Electric Vehicle Charging Stations and Renewable Energy Sources Including CO2 Emissions

In the coming years, several transformations in the transport sector are expected, associated with the increase in electric vehicles (EVs). These changes directly impact electrical distribution systems (EDSs), introducing new challenges in their planning and operation. One way to assist in the desired integration of this technology is to allocate EV charging stations (EVCSs). Efforts have been made towards the development of EVCSs, with the ability to recharge the vehicle at a similar time than conventional vehicle filling stations. Besides, EVs can bring environmental benefits by reducing greenhouse gas emissions. However, depending on the energy matrix of the country in which the EVs fleet circulates, there may be indirect emissions of polluting gases. Therefore, the development of this technology must be combined with the growth of renewable generation. Thus, this proposal aims to develop a mathematical model that includes EVs integration in the distribution system. To this end, a mixed-integer linear programming (MILP) model is proposed to solve the allocation problem of EVCSs including renewable energy sources. The model addresses the environmental impact and uncertainties associated with demand (conventional and EVs) and renewable generation. Moreover, an EV charging forecast method is proposed, subject to the uncertainties related to the driver's behavior, the energy required by these vehicles, and the state of charge of the EVs. The proposed model was implemented in the AMPL modelling language and solved via the commercial solver CPLEX. Tests with a 24-node system allow evaluating the proposed method application.


Introduction
Concerns about climate change and initiatives like the Paris agreement motivate many sectors of the economy to reduce CO2 emissions. CO2 emissions (the main greenhouse gas) in the transport sector are higher than those in other sectors. Still, these numbers are beginning to show a decline in regions with higher electric vehicle (EV) penetration rates (Canizes et al. 2019;Coignard and Macdougall 2019; Energy Information Administration (EIA) 2020). EV is the key to transport efficiency and sustainability. These technologies are being developed to reduce dependence on oil for transportation and limit CO 2 emissions. Similarly, renewable energy sources are being adopted to promote alternatives to fossil fuel-based electricity generation (Cai et al. 2017;Richardson 2013).
Mass adoption of electric vehicles will provide many benefits, such as the incentive to integrate renewable energy sources (RES) into the electrical distribution system (EDS). However, a large increase in EVs can harm the EDS, causing an increase in power losses, voltage degradation, and overloads, thus compromising the electric power quality (Etezadi-Amoli et al. 2010;Richardson 2013). Therefore, EDSs must be prepared for the demand requirements and the particularities of EVs.
The main obstacles for the integration of EVs are their high acquisition cost, limited cruising range and the lack of an infrastructure prepared to supply their demand. Thus, to contribute to the integration of this technology, it is necessary to invest in the system's infrastructure, which can be done in two ways: by expanding/reinforcing the grid and by installing EV charging stations (EVCSs) (Cai et al. 2017;Mauri and Valsecchi 2012).
Different methods have been used to solve the allocation problem of EVCSs such as classic optimization techniques (Neyestani et al. 2015) and metaheuristics (Amini and Islam 2014;Moradijoz et al. 2013). A mixed-integer linear programming (MILP) model to solve the allocation problem of EVCSs, aiming to minimize system costs, is proposed in (Neyestani et al. 2015). To this end, a two-stage model was developed. The first stage models the EVCS behavior considering market interactions (reserve and energy markets), aiming to maximize profit. In the second stage, the allocation problem of EVCSs is solved. Therefore, the method proposed in (Neyestani et al. 2015) is treated as an EVCS operation planning problem. Specialized methods using the genetic algorithm have been proposed in (Amini and Islam 2014;Moradijoz et al. 2013) to solve the allocation problem of EVCSs.
The uncertainties associated with EV demand must be addressed in the allocation problem of EVCSs, aiming to obtain a solution more committed to reality. A probabilistic approach based on the point estimate method has been proposed in (Mirzaei et al. 2016) to determine the optimal allocation of EVCSs, considering the driving patterns of EV owners. In that work, the allocation problem has been formulated as a two-stage stochastic programming model. Uncertainties based on the behavior of EV users are also considered in (Zeng et al. 2020). The allocation problem has also been formulated as a two-stage stochastic programming model. The first stage defines the planning actions, while the second stage evaluates the system operation. An important aspect that was not considered in any of the previously cited proposals is the allocation of EVCS integrated with RES investment.
The combination of EVs and RES offers a transformative impact on the world, making it possible to reduce the dependence on fossil fuels in the energy and transportation sectors and, consequently, reduces greenhouse gas emissions (Richardson 2013). Therefore, it is important to address EVCSs and RES allocation simultaneously (Bagherzadeh et al. 2019;Mozafar et al. 2017;Shaaban et al. 2019). In this context, a salp swarm algorithm has been proposed to solve joint planning of EVCSs and RES (Bagherzadeh et al. 2019). Multiobjective approaches have also been used to solve this optimization problem, such as (Mozafar et al. 2017;Shaaban et al. 2019). A non-dominated sorting Genetic algorithm (NSGA-II) was used in (Shaaban et al. 2019), while (Mozafar et al. 2017) proposed a genetic algorithm-particle swarm optimization hybrid to solve the allocation problem of EVCS and RES. Alternatively, some studies consider the presence of these technologies into the allocation problem of EVCS but without addressing the allocation of RES (Pashajavid and Golkar 2013).
Among the works cited, only (Shaaban et al. 2019) takes into account the environmental aspects associated with greenhouse gas emissions. Moreover, about the uncertainties addressed, only (Pashajavid and Golkar 2013) includes the uncertain behaviour of conventional demand. Thus, this work proposes a simultaneous approach for optimal allocation of EVCSs and RES, including environmental issues. A two-stage stochastic programming MILP model is used to solve the problem, which considers stochastic behavior of conventional demand, EV demand, and renewable generation. Table 1 provide a comparison of previous related approaches with the proposed method in which the main contributions are highlighted and described as follows: A two-stage stochastic programming model for the simultaneous allocation of EVCSs and RES that deals with uncertainties of renewable generation (wind turbine (WT) and photovoltaic (PV) generator) along with environmental issues related to CO 2 emissions. A mixed linear programming model (MILP) formulation for the addressed problem that can be optimally solved by efficient commercial solvers. A method for forecasting EV charging demand based on uncertain user behavior. An enhanced representation of active/reactive power limits of WT according to their capability curves and power factor limits.

Problem formulation
The joint allocation problem of EVCSs and RES is formulated in this section as a MILP model in which the uncertainties are handled through a two-stage stochastic programming model. In the first stage, the investment decisions related to EVCSs and RES are carried out before the uncertainties are realized (here and now decisions). In the second  (Tabares et al. 2016) and assumes that: the operation of the EDS is represented by a balanced AC linearized power flow model and the loads are modeled as constant powers; CO 2 emissions are penalized in the objective function with an emission cost; and the RES that will be allocated belong to the distribution system operator. Similar to (Tabares et al. 2016), the original non-linear model is transformed into an equivalent linear model using the piecewise g-function for the square calculation of a given value (Wolsey 1998).

Objective function
The objective function aims the minimization of the present value of investment and operational costs. Thus, the total cost (1) considers the investment costs in EVCSs (2) and RES (3), jointly with the operational costs related to the energy imported by the substation (4), RES maintenance costs (5), and CO 2 emissions costs (6). The function f(τ, λ) = 1 − (1 + τ) −λ /τ allows the calculation of present value of the annualized cost. The description of sets, indices, parameters, and variables used in the model is available in the Nomenclature at the end of this paper.
min TC : where:

Constraints
The proposed model has the following types of constraints: Steady state operation, operational limits, investment limits, RES model, and EVCS model.

Steady state operation
Expressions (7)-(10) set the steady-state operation of the EDS and are based on (Franco et al. 2014). The active and reactive power balances are expressed by (7) and (8). The voltage drop is determined by (9), while (10) represents the relationship between the voltage at node i, the current and the active/reactive power flow through circuit ij in scenario ω and period p. Function g(P ij, ω, p , Q ij, ω, p , Γ) in (10) approximates the sum of the squares of variables P ij, ω, p and Q ij, ω, p , using Γ blocks in a piecewise representation (Wolsey 1998).
Operational limits Expressions (11)- (17) define the operating limits in the EDS. Thus, voltage limits are guaranteed by (11), while (12)-(14) limit the current, active and reactive power flows through circuit ij. The square of the apparent power supplied in each substation is determined by (15) and limited by (16). Finally, (17) ensures that there is only a single direction of flow through circuit ij.

RES model
The operation of the WT units is represented by limiting the injected active/reactive power according to the capability curves and power factor restrictions (Rueda-medina et al. 2013). Thus, WT units are modeled using the capacity curve of a double-fed induction generator. To establish the generation limits, the points ðP wt k;1 ; Q wt k;1 Þ; ðP wt k;2 ; Q wt k;2 Þ; ðP wt k;3 ; Q wt k;3 Þ, and ðP wt k;4 ; Q wt k;4 Þ were defined. Capacity curve is linearized using constraints presented in (24)-(28). The details of this formulation can be found in (Rueda-medina et al. 2013). The operational limits of active and reactive power by RES are shown in (29) and (30) for WT units, and (31) and (32) for PV units.
jQ wt k;ω;p j ≤ P wt k;ω;p tanðcos −1 ðφ wt ÞÞ ∀k; ω; p ð30Þ  (34) determines that the charging demand in EVCS corresponds to the EV demand obtained with the method presented in the following section "Uncertainty modelling".
Uncertainty modelling The uncertainties brought by intermittent renewable generation, future demand growth, and EV charging demand represent challenges to the reliable operation and planning of power systems (Zhou et al. 2017). Thus, a proposal for the allocation problem of EVCSs and RES that is more appropriate and committed to reality must include these uncertainties. In this work, the uncertainties related to consumer demand, EV demand, and renewable generation are considered.

Method for forecasting EV charging demand
Electrical system operators are interested in forecasting the EV charging demand to assess the impacts and needs for updating the EDS infrastructure (Knezović et al. 2017).
To quantify the EV charging demand, a method to model its uncertainties is proposed here, based on (De Quevedo et al. 2019), that comprises seven steps described as follows: 1. For each EV, randomly select an initial state of charge (SOC) value.
2. For each EV v and day d, randomly (uniform distribution) assign a total daily travel a (from 10 km to 200 km). 3. If the SOC of the EV is enough to make the travel a, the SOC is updated and we go to step 5. Otherwise, go to step 4. 4. The EV should be charged with at least the SOC necessary to carry out the travel.
After charging the vehicle, the SOC and the charging demand (kW) are updated. EV arrival times are selected based on the user behavior of the fuel vehicles presented in (Federal Highway Administration 2009). 5. If the last day was evaluated, go to step 6. Otherwise, evaluate the next day and go to step 2. 6. If the last EV was evaluated, go to step 7. Otherwise, evaluate the next EV and go to step 1. 7. Repeat steps 1 to 6 until each EV completes all daily trips. Finally, calculate the total demand (kW) for all EVs, for each hour and day.
The algorithm above was implemented in MATLAB (Mathworks 2017) considering two types of EVs with 40 kWh (Nissan leaf) and 60 kWh (Chevrolet Bolt) batteries. An analysis to obtain the charging profiles is carried out for the planning period of 1 year. As mentioned earlier, for each vehicle, a random SOC within the interval [0.2, 1.0] is initially selected; note that this happens only on the first trip of each vehicle. Eq. (35) determines the SOC of the vehicle v on day d, while (36) calculates the remaining battery capacity. Finally, (37) is used to determine the energy used by vehicle v on day d.
The simulation was performed considering 1000 vehicles (half with 40 kWh and half with 60 kWh). After its execution, the EV charging demand for each hour of the year was obtained. Figure 1 illustrates the obtained EV average daily demand profile as well as the demand profile calculated for day 235; this day was chosen because it has the EV charging peak demand along the year: 4650 kW at 19:00. Note that there is a variation in the demand profiles, and the maximum demand registered on day 235 is about 20% higher than the average demand profile.
The SOC is separated into ten different categories (see Table 2) to facilitate the analysis. Figure 2(a) and (b) show the probability distribution function (PDF) of the initial and final SOC of the vehicle, respectively. Figure 4(a) shows that the vehicle is more probable to be charged when the SOC is below 0.5 (Categories 1-5). On the other hand, after charging (see Fig. 2(b)), the SOC is more likely to be larger than 0.5 (Categories 6-10). Figure 2(c) and (d) show the cumulative distribution function (CDF) of the initial and final SOC, respectively. The probability of the initial SOC to be between 0 and 0.5 (Categories 1-5) is 85.6% (Fig. 2(c)), while the probability of the final SOC being above 0.5 (Categories 6-10) is 91.4% (Fig. 2(d)). The information provided by the calculated profiles will be used to represent the power related to EV charging in the proposed model. Finally, the proposed method to forecast EVs demand is summarized in Fig. 3.

Scenario generation method
The set of scenarios is created from annual historical data of stochastic parameters using historical demand data (ONS 2019), solar irradiation data (Renewables 2021), wind speed data (Renewables 2021), and the EV charging demand previously provided. Since a high number of scenarios can be obtained, the k-means scenario reduction method is applied to achieve computational tractability. The hourly data are divided  into 2 seasons; each season has 2 sub-blocks (day and night). Moreover, each sub-block is reduced to 10 clusters using k-means (the number of clusters was arbitrarily defined). The adopted process can be described in the following steps: 1 Historical data on demand, solar irradiance, and wind speed are normalized by their corresponding maximum values. 2 The annual duration curve for each stochastic parameter is classified into two seasons (blocks of time): winter and summer. The data contained in the time blocks are classified into two sub-blocks (night and day). 3 The number of required clusters is defined. The k-means method is applied to each sub-block. Thus, the data contained in each sub-block is categorized according to the number of clusters previously established. At the end of the process, centroids of each cluster are determined. 4 The set of scenarios is stored in a matrix that represents 40 operating conditions with four columns for the uncertain parameters (10 clusters x 2 seasons x 2 sub-blocks). 5 The probabilities for each scenario are calculated by dividing the number of hours in the respective scenario by the sum of hours within a block of time.
Renewable generation profiles are obtained following the method in (Montoya-Bueno et al. 2015). The power profiles of the WT units are obtained, for each scenario, through linearization of the power curve, which is a function of the wind speed (38). Moreover, power profiles of the PV units are obtained in (39), depending on the cell junction temperature (40).

Tests and results
The proposed approach was implemented in the mathematical language AMPL (Fourer et al. 2003) and solved using the commercial solver CPLEX (IBM 2019), in a computer with an Intel Xeon E5-2650 processor and 64GB of RAM. The proposed model is validated using the 24-node EDS adapted from (Tabares et al. 2016), which has 20 load nodes, four substations, a nominal voltage of 20 kV, and the horizon of 10 years is divided in two periods (5 years for each one). The upper and lower voltage magnitude limits are 1.05 and 0.95 pu. The EVCS operational cost is considered to be 10% of the installation cost and is included in the investment cost (Banol Arias et al. 2017) Furthermore, it is assumed that the EVs demand grows 10% in the second period of the planning. The energy supplied cost in the substation is 0.1 USD/kWh and the interest rate is 10%. Data associated with demand, substation capacity, scenarios used, EVCSs, and RES are available in (De Lima et al. 2021). The results for the planning were analyzed according to the following case studies: Case I) The allocation problem of RES is optimized without investments in EVCSs; Case II) Joint optimal allocation of EVCSs and RES; Case III) Sensitivity analysis with different CO 2 emission rates; and Case IV) Sensitivity analysis with EV penetration levels.
The investment plan for Case I determined the following actions: (1) First period: allocation of eight WT units at nodes 3, 4, 5, 9, 11, 14, 16 and 19; and installation of one hundred one PV units. (2) Second period: Allocation of six PV units at node 10. On the other hand, the investment plan for Case II defined the following proposals for each period: (1) First period: allocation of eight WT units at nodes 3, 4, 5, 9, 11, 14, 16 and 19; installation of eighty nine PV units; and allocation of five EVCSs 12 at nodes 3, 6, 8, 14, and 15. (2) Second period: Allocation of sixteen PV units (three at node 6, eight at node 10, and five at node 13). A summary of the main results for Cases I and II is shown in Table 3. It is possible to note, that the total cost for Case I (without investments in EVCSs) is higher than the cost related to Case II, with a difference of approximately 3.62 %. Moreover, Case II, which simultaneously optimizes the allocation of EVCSs and RES, resulted in lower emissions than Case I (a difference of 5.53%). The best expansion plan (Case II) is illustrated in Fig. 4.
A sensitivity analysis was performed with different CO 2 emission rates (Case III). The results of this analysis are summarized in Table 4. The first and the last solutions (1 and 9) represent the extreme solutions, with solution 1 presenting the lowest cost and maximum CO 2 emission value and solution 9 the opposite. The first solution represents Case II, in which the problem is optimized without directly restricting CO 2 emissions, aiming to minimize investment and operational costs (including emissions costs). Moreover, it is noted in Table 4 that to reduce CO 2 emissions, investment in RES increases. Figure 5 shows the conflict between minimizing costs and reducing emissions. It is noted that the decrease in emission values leads to an increase in the total planning cost.
A sensitivity analysis was performed varying the EV penetration level. Figure 6 illustrates how the total cost increases as the EVs demand increases, with vehicle penetration ranging from 60% to 140%, emphasizing the importance of adequately forecasting the EVs charging demand.

Conclusions
A two-stage stochastic programming model for the simultaneous allocation problem of EV charging stations (EVCSs) and renewable energy sources (RES) in electrical distribution system (EDS) has been proposed in this paper. The proposed mixed-integer linear programming model aims to minimize investment and operational costs (including CO 2 emission costs). Additionally, uncertainties related to conventional demand, EV demand, and renewable generation were considered.
Results showed that the installation of EVCSs has a positive impact on the EDS, reducing the expansion plan costs, despite the investment costs in this technology. Moreover, the integration of EVCS with RES contributes to the reduction of CO 2 emissions, showing that the simultaneous allocation of these technologies can bring both economic and environmental benefits. A sensitivity analysis has been performed with different CO 2 emission rates. Results demonstrate that an increase in RES investments leads to CO 2 emission reduction. However, the reduction of CO 2 emissions results in increasing costs, evidencing the conflict between cost and emission reduction. I sqr ij;ω;p Square of the current through circuit ij in scenario ω and period p. N EV r;c;p Number of EV chargers at node r, type c and period p. ORES p RES maintenance costs at period p. P ij, ω, p Active power flow through circuit ij in scenario ω and period p. P pv u;ω;p Active power injected by the PV units at node u, scenario ω, and period p. P s s;ω;p Active power supplied by the substation at node s, scenario ω, and period p. P wt k;ω;p Active power injected by the WT units at node k, scenario ω, and period p. P R Rated electrical power. Q ij, ω, p Reactive power flow through circuit ij in scenario ω and period p. Q pv u;ω;p Reactive power injected by the PV units at node u, scenario ω, and period p. Q s s;ω;p Reactive power supplied by the substation i at node s, scenario ω, and period p Q wt k;ω;p Reactive power injected by the WT units at node k, scenario ω, and period p. Sg sqr i;ω;p Square of the apparent power supplied by substation at node i, scenario ω, and period p. V sqr i;ω;p Square of the voltage at node i, scenario ω, and period p. x EV r;p Investment variable for installing an EVCS at node r and period p. x wt k;p Investment variable for installing a WT unit at node k and period p. y − ij;p Operational variable related to the forward direction of circuit ij and period p y þ ij;p Operational variable related to the backward direction of circuit ij and period p.