A Specialized Long-Term Distribution System Expansion Planning Method With the Integration of Distributed Energy Resources

The electrical distribution system (EDS) has undergone major changes in the last decade due to the increasing integration of distributed generation (DG), particularly renewable energy DG. Since renewable energy resources have uncertain generation, energy storage systems (ESSs) in the EDS can reduce the impact of those uncertainties. Besides, electric vehicles (EVs) have been increasing in recent years leveraged by environmental concerns, bringing new challenges to the operation and planning of the EDS. In this context, new approaches for the distribution system expansion planning (DSEP) problem should consider the distributed energy resources (DG units, ESSs, and EVs) and address environmental impacts. This paper proposes a mixed-integer linear programming model for the DSEP problem considering DG units, ESSs, and EV charging stations, thus incorporating the environmental impact and uncertainties associated with demand (conventional and EVs) and renewable generation. In contrast to other approaches, the proposed model includes the simultaneous optimization of investments in substations, circuits, and distributed energy resources, including environmental aspects (CO2 emissions). The optimization method was developed in the modeling language AMPL and solved via CPLEX. Tests carried out with a 24-node system illustrate its effectiveness as a valuable tool that can assist EDS planners in the integration of distributed energy resources.


I. INTRODUCTION
The automotive industry has been going through a moment of transformations with the increase of electric vehicles (EVs). Worldwide projections indicate that by 2040 EVs should increase from 3 million to 66 million units circulating worldwide [1]. This increase in EVs is expected due to The associate editor coordinating the review of this manuscript and approving it for publication was Salvatore Favuzza .
incentives to promote their adoption combined with the growing cost-benefit ratio of this technology and its advantages related to environmental issues. Nevertheless, the main problems for the adoption of EVs are their high acquisition cost and the scarcity of infrastructure (i.e., the lack of EV charging stations (EVCSs)) prepared to meet their demand [2]. Since EVs charge on the electrical grid, this trend directly impacts the operation and planning of electrical distribution systems (EDSs). A large EV penetration can affect the EDS, causing VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ an increase in energy losses, degradation of voltage profiles, and overloads, thus compromising the quality of electricity supply [3], [4]. 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 EVCSs [5]. Thus, the distribution system expansion planning (DSEP) problem must be adapted to the demand requirements and the particularities of EVs. The DSEP problem determines the necessary expansion actions in the EDS to meet the growing demand, while maintaining a safe operation and guaranteeing the quality of the service for the users, all at the lowest possible cost [6]. The classic formulation for the DSEP problem includes investments in substations and circuits [7]. Subsequently, with the advancement of research and the growth of distributed generation (DG), there was a need to incorporate it into the DSEP proposals, adding DG units as problem decision variables [8].
The characteristics of the DSEP problem correspond to a mixed-integer non-linear programming (MINLP) model, which is highly complex due to its combinatorial nature and the discrete condition of the decision variables [9]. To solve it, several methods and models have been proposed in the specialized literature, whether using meta-heuristics such as Simulated annealing [10] and Ant colony algorithm [11], heuristic techniques such as branch exchange algorithm [12] or classic optimization methods such as linear programming [13], [14] and quadratically-constrained programming [15]. A review of the state of the art of this problem is presented in [16].
A few studies have addressed the DSEP problem including EVCSs. In some of them, the EV charging demand was included in the planning proposal, although without considering investments in EVCSs. On the other hand, most of the works address the allocation of EVCSs independently of the DSEP problem [17], [18]. A two-stage stochastic programming model to solve the allocation problem of EV parking lots has been proposed in [17], considering the uncertainties related to EV demand. The point estimation method has been used in [18] to consider the uncertainties associated with EV demand. Some approaches have included investment in EVCSs within the DSEP problem [19]- [23]. A mixed-integer linear programming (MILP) model was developed in [19] for the assessment of investments in substations, circuits, capacitor banks, non-renewable DG units, and EVCSs. Similarly, in [20], a MILP model identifies optimal investments in substations, circuits, battery banks, renewable DG units, energy storage systems (ESSs), and EVCSs. The stochastic behavior related to conventional demand was included in the models proposed in [19] and [20]; to deal with those uncertainties, robust optimization, and stochastic programming were adopted. Unlike [20], uncertainties related to EV demand were not included in the proposed model in [19]. A multiobjective DSEP model with high wind power penetration has been developed in [21], including the planning of EVCS. Optimal planning of EVCS and renewable DG was presented in [22] using k-means clustering technique to consider the uncertainties related to renewable power sources and EV demand. Finally, a genetic algorithm was implemented in [23] to solve the DSEP problem, which includes investments in EVCSs, ESS, and renewable GD units; in addition, the Monte Carlo method is used to model the uncertainties of renewable generation and EV demand.
In contrast to [17], [18], this paper proposes a MILP model for the joint planning problem of EVCS and EDS using an exact method guaranteeing that the optimal solution of the problem is found. Also, unlike [19], [20], the proposed model includes the EV demand uncertainty and the environmental issues associated with CO 2 emissions. Simultaneous planning of EVCS and EDS considering the presence of wind generation has been proposed in [21]. On the other hand, in contrast to [21], our proposal also includes investments in photovoltaic generation, non-renewable generation, and ESSs. Different from [22], this work takes into account investments in substations, circuits, non-renewable DG units, and ESSs. Furthermore, unlike [23], which focuses only on meeting EV demand, our paper also considers conventional demand in the planning problem.
The main contribution here is the simultaneous optimization of multi-period investments in substations, circuits, renewable and non-renewable DG units, EVCSs, and ESSs, in contrast to all the proposals mentioned previously. The stochastic behavior of conventional demand, EV demand, and renewable generation is modeled using two-stage Stochastic Programming. The following expansion alternatives are evaluated in the expansion plan: (1) construction of substations; (2) construction of circuits; (3) allocation of renewable and non-renewable DG; (4) allocation of EVCSs; and (5) allocation of ESSs. Furthermore, in this paper, the EDS is represented by an AC power flow model, different from most of the works found in the specialized literature, which uses the DC model [20]. DC approximation is adopted often as a compromise between solution accuracy and computational time, namely for transmission network systems [24]. Unfortunately, DC approximation assumptions do not hold well and become invalid for EDS [25]. Table 1 provides a comparison of previous related approaches with the proposed method, and the main contributions of this work are highlighted: 1) A mathematical model for the long-term stochastic planning model, including environmental aspects (CO 2 emissions). Such formulation considers the simultaneous optimization of multi-period investments in substations, circuits, renewable and non-renewable DG units, EVCSs, and ESSs. The joint optimization of these investments has not been addressed yet.
2) A MILP model for the multi-period DSEP problem to find the optimal solution of the problem. The operation of the EDS is represented by an AC power flow linear model.
3) A two-stage Stochastic Programming model for the multi-period DSEP that addresses uncertainties related to conventional/EV demand and renewable DG units (Wind turbine (WT) and photovoltaic (PV)) through a set of demand and generation scenarios.
The remaining part of this manuscript is organized as follows: section II presents the uncertainty modeling associated with conventional demand, EV demand, and renewable generation; section III shows the problem formulation; the application of the proposed method in a Case Study for a 24-node test system is presented in section IV and the conclusions are drawn in section V.

II. MODEL FOR THE UNCERTAINTIES RELATED TO DEMAND AND RENEWABLE GENERATION
A DSEP proposal that is more appropriate and committed to reality must address the uncertainties of the physical and operational parameters. In this work, the uncertainties related to consumer demand, EV demand, and renewable generation are taken into account. This section presents the proposed approach to model these uncertainties.

A. EV DEMAND UNCERTAINTY MODEL
Predicting the EV charging demand pattern is a complicated task due to the low number of real cases recorded, the uncertainties related to the driver's behavior, the level of EV penetration, the energy required by these vehicles, among others [26]. On the other hand, system operators are interested in forecasting the EV charging demand to assess the impacts and needs for updating the EDS infrastructure [27]. Therefore, to quantify the EV charging demand, a method to model its uncertainties is proposed here. The method presented here to estimate EV demand is adapted from [20]. In contrast to [20], which uses data from travel patterns of traditional vehicles [28], in this proposal, the electric vehicle infrastructure projection tool (EVI-Pro) Lite [29] was used to define the EV arrival times. Moreover, in [20], only one type of EV charger was considered (slow charger (3.3 kW)). On the other hand, in this proposal, the chargers used at the stations have powers of 7 kW (slow charger), 50 kW (quick charger), and 150 kW (super-fast charger). The proposed method comprises seven steps described as follows: 1) For each EV, randomly (uniform distribution) select an initial state of charge (SOC) value. 2) For an EV v and day d, randomly (uniform distribution) assign a total daily travel distance a (from 10 km to 200 km) [22] 3) If the SOC of the vehicle is enough to make the travel a, the SOC is updated, 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 electric vehicle infrastructure projection tool (EVI-Pro) Lite [29]. 5) If the last day of the year 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) Calculate the total demand (kW) for all EVs, for each hour and day (see equation (6)).
The simulation was performed considering 1,000 vehicles. 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 VOLUME 10, 2022 calculated for day 125; this day was chosen because it has the EV charging demand peak along the year: 2,430 kW at 19:00. Note that there is a variation in the demand profiles and the maximum demand registered on day 125 is about 33% higher than the average demand profile.
The SOC is separated into different categories (see Table 2) to facilitate the analysis. Figure 2(a) and 2(b) show the probability distribution function (PDF) of the initial and final SOC of the vehicle, respectively. Figure 2 (a) shows that the vehicle is more probable to be charged when the SOC is below 0.5 (Categories 1-4). On the other hand, after charging (see Figure 2 (b)), the SOC is more likely to be larger than 0.5 (Categories 5-9). The information provided by the calculated profiles will be used to represent the power related to EV charging in the DSEP model. Finally, the method proposed in this section is summarized in Figure 3.

B. SCENARIO GENERATION METHOD
The method presented previously in section II-A aims to estimate the total charging demand for all EVs. Thus, the EV demand profiles obtained are correlated with the other uncertain parameters (conventional demand, solar irradiation, and wind speed). To represent the uncertainties, a set of annual scenarios is obtained by using historical demand data (Brazil) [31], solar irradiation data [32], wind speed data [32], and the EV charging demand provided by the method presented in section II-A. Since a high number of scenarios can be obtained, the k-means scenario reduction method is applied to achieve computational tractability. For medium and long-term planning, the k-means clustering technique is widely accepted to model the uncertainties [33]. Such clustering technique maintains correlation among uncertain data [34]. Similar to [8], 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; that 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 (8 clusters × 2 seasons × 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. Finally, the scenario generation process is illustrated in Figure 4. Thus, the representative scenarios obtained here are applied to the long-term stochastic planning model (section III). Renewable generation profiles are obtained following the method in [35]. The power profiles of the WT units are obtained, for each scenario, through linearization of the VOLUME 10, 2022 power curve, which is a function of the wind speed (7). Moreover, power profiles of the PV units are obtained in (8), depending on the cell junction temperature (9).

III. MATHEMATICAL MODEL FOR THE DSEP PROBLEM WITH INTEGRATION OF DISTRIBUTED ENERGY RESOURCES
This section describes the proposed mathematical formulation that aims the optimal expansion plan for the multi-period DSEP problem. The problem is formulated first as an MINLP model and then, to reduce the complexity of the problem and ensure that the optimal solution is found, the original model is approximated by a MILP model. The model proposed here is based on the formulations presented in [19]. Similar to [19], herein, a MILP model for the integrated planning of EDS and EVCSs is proposed. However, different from [19], this proposal considers investments in DG units (renewable WT and PV and non-renewable gas turbine (GT)), uncertainties related to generation and EV charging demand, and investment in ESS. The model proposed here for the long-term stochastic DSEP model assumes that: 1) the operation of the EDS is represented by an AC power flow linear model; 2) The expansion decisions are defined through a multi-period model; 3) the EDS is operated using a radial topology; 4) CO 2 emissions are penalized with a given emission cost; 5) the planning of EVCSs and EDS is centralized. These assumptions are usually adopted for the long-term planning problem [19], [20].
Uncertainties related to conventional demand, EV charging demand, and renewable generation are represented in the model through two-stage Stochastic Programming. In the first stage, decisions about investments in substations, allocation of DG units, EVCSs, and ESS are executed before the uncertainties are known (here and now decisions). On the other hand, the second stage evaluates the expected operation cost of the EDS after to the uncertain parameters assume a specific value (wait-and-see decisions).

A. OBJECTIVE FUNCTION
The objective function of the problem is presented in (10) and minimizes the present value of the expected total cost, being composed by the following costs: 1) Investments (I p ) (11) associated with substation, circuits, DG units, EVCSs, and ESSs; 2) Operational costs (O p ) (12) related to the energy supplied by substations, jointly with distributed energy resources maintenance costs, and CO 2 emissions costs. These costs are calculated for each period p. The function f (τ, λ) = 1 − (1 + τ ) −λ /τ allows the calculation of the present value of the annualized cost.

1) STEADY-STATE OPERATION
Expressions (13)-(16) set the steady-state operation of the EDS and are based on [15]. The active and reactive power balance is expressed by (13) and (14). The voltage drop is determined by (15), while (16) 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. Auxiliary equations (17)

2) OPERATIONAL LIMITS
Expressions (20)-(29) define the operating limits for the different assets that may be installed in the EDS. Thus, voltage limits are guaranteed by (20), while (21)-(23) limit the current, active, and reactive power flows through circuit ij according to the conductor type. Besides that, constraints (24)-(26) establish the limits for the current, active and reactive power flows, respectively, in terms of the circuit operating conditions. The square of the apparent power supplied in each substation is determined by (27) and limited by (28). Finally, constraint (29) limits variable b ij,ω,p , which is used in (15) to maintain the feasibility of the problem when the circuit ij is not installed; if the circuit is operated then b ij,ω,p = 0, otherwise, b ij,ω,p is free to take any value limited by b.

3) CONSTRAINTS OF INVESTMENTS AND OPERATION
Constraints (30)- (34) guarantee that only one investment decision is carried out for substations, circuits, EVCSs, ESSs, and WT/GT units, in each node. On the other hand, (35)- (37) limit the number of PV units, EV chargers, and battery banks that can be installed at each node during the planning horizon. The operating conditions of circuits is addressed by (38). In this regard, y L+ ij,t + y L− ij,t is equal to one if the circuit ij is operated; otherwise, it is zero. Expression (39) relates the operational variable z ij,a,p with its corresponding investment variable x L ij,a,t . Finally, expressions (40) and (41) refer to the binary characteristic of the investment and operation variables of substations, circuits, EVCSs, and WT/GT units and integer nature of the components that can be added to the EDS (PV units, EV chargers, and battery banks).
z ij,a,p = P t=1 x L ij,a,t ; ∀ij, a, p are defined to represent the feasible region for power injection. For the modeling of WT units, the capability curve of a double-fed induction generator was considered. 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. Capability curves are linearized using constraints presented in (42)-(51). The details of this formulation can be found in [36].

4) DG UNIT MODEL
The operational limits of active and reactive power by DG units are shown in (52) and (53) for GT units, (54) and (55) for WT units, and (56) and (57) for PV units.

5) EVCS MODEL
It is assumed a centralized approach for the planning of EVCSs and EDS, which is usually adopted for the long-term planning problem [19], [20], aiming computational tractability since the DSEP problem is highly complex due to its combinatorial nature, containing many binary and integer variables. Constraint (58) ensures that the charging demand in EVCSs does not surpass the capacity of the stations, while (59) determines that the charging demand in EVCS corresponds to the EV demand obtained with the method presented in section II-A.

6) ESS MODEL
To mitigate the impacts of the uncertain behavior of renewable DG units, investment in ESSs is considered an expansion alternative. It is assumed that this technology will store power during the day and will inject power into the system at night when there is no photovoltaic generation. The charge and discharge processes of the ESS will have a pre-established daily operational time (d ESC and d ESD ). Expressions (60)-(65) represent the ESS model based on [37] with adaptations to include the variables associated with the energy stored (charged) and injected (discharged) by the ESS. Constraints (60) and (61) limit the ESS charge and discharge power, while expressions (62) and (63) determine the energy stored and injected by the ESS. Moreover, the energy stored in the ESS is limited by (64).
There is an additional difficulty in medium and long-term EDS planning since k-means is used to reduce the operational scenarios. When the data is organized into scenarios, the chronological information is lost. Therefore, it is not possible to maintain the exact chronological sequencing of ESS charge and discharge. Thus, similar to [20], the ESS operation does not consider the chronological sequence of charge and discharge in this model. Besides, expression (65) determines that the net energy at an ESS is zero along a block of time. This constraint has been formulated for each scenario within the block of time bl. This simplification makes sense for a long-term planning problem. A detailed inclusion of the sequencing of ESS charge and discharge can be computationally infeasible and is unnecessary from the point of view of long-term planning since it is not possible to have detailed and accurate information for all the parameters involved in a long-term planning horizon [37].

7) RADIALITY CONSTRAINTS
Expressions (66)-(68), jointly with (16) and (17), guarantee the radial operation of the EDS. For this purpose, (66) ensures that circuits connected to substation nodes are always operating in the forward direction. In addition, (67) guarantees that each load node is connected to only one circuit operating in the forward direction and (68) allows the use of transfer nodes (i.e., nodes without demand that may connect two loads nodes) [38].

C. LINEARIZATION OF EXPRESSIONS (16) AND (27)
The MINLP model for the DSEP problem defined by the set of equations (10)-(65) contains nonlinear constraints in (16) and (27). Thus, to reduce the complexity of the problem and ensure that the optimal solution is found, the original model is approximated by a MILP model, as proposed in [9], [39]. Linearizations and approximations are described as follows:

1) LINEARIZATION OF (16)
The linearization of the product V sqr i,ω,p ·Î sqr ij,ω,p is performed by using a constant value for the voltage (V i,ω,p ). The value of this constant is defined using the maximum and minimum voltage limits [9], as shown in (69). Finally, the sum of the square of the power is linearized using the piecewise linear approximation that is given by (70)-(79). The variableŝ P 2 ijω,p andQ 2 ij,ω,p are approximated by the product of slope of the γ th of the piecewise discretization(m G y ) and the discretization variable ( P ij,y,ω,p ) for active power or ( Q ij,y,ω,p ) for reactive power as shown in (70). Equations (71) and (72) representP ijω,p andQ ij,ω,p by non-negative auxiliary variables (P + ij,ω,p , P − ij,ω,p , Q + ij,ω,p , Q − ij,ω,p ). Additionally, (73) and (74) determine that the termsP ijω,p andQ ij,ω,p are equal to the sum of all values in each discrete block. Expressions (75) and (76) limit the values of the discretization blocks. Finally, expressions (77) and (78) calculate the parameter values used in the discretization.
∀i, ω, p (69) (79) 2) LINEARIZATION OF (27) The same technique presented above is used to linearize (27). Variable Sg sqr i,ω,p (square of the apparent power supplied by substation) can be approximated from the linear expression (80). Equations (81) and (82) determine that the terms P s i,ω,p and Q s i,ω,p are equal to the sum of the values of each discretization block. Besides, equations (83) and (84) limit the values of the discretization blocks. Finally, expressions (85) and (86) define the parameter values used in the discretization.

D. MILP MODEL FOR THE DSEP PROBLEM INTEGRATED WITH THE EXPANSION OF EVCSs
The MINLP model is transformed into an equivalent linear model using the piecewise similar to [9], [39]. The proposed linear model is represented by the following expressions: min CT: Equation (10) subject to :

IV. TESTS AND RESULTS
The proposed MILP model was implemented in the modeling language AMPL [40] and solved using the commercial solver CPLEX [41] in a computer with an Intel Xeon E5-2650 processor and 64GB of RAM. The proposed model was validated using the 24-node EDS adapted from [9], which has a nominal voltage of 20 kV. The 10-year planning horizon is divided in two periods (5 years each). The upper and lower voltage limits are defined as 1.10 and 0.90 p.u., respectively. Furthermore, it is assumed that the EV demand grows 11% in the second period of the planning [19]. The energy supplied cost in the substation is 0.1 USD/kWh and the interest rate is 10% [9]. The chargers that can be installed at the EVCS have powers of 50 kW (quick charger) and 150 kW (super-fast charger). The initial system topology, investment cost data (substations, circuits, DG units, ESS, and EVCS), emission rate data (DG units, grid power), and other case study data are available in [42].
The proposed model was analyzed according to the following case studies: I) The EV demand and conventional demand feed only from the grid. DSEP with investments only in circuits, substations, and EVCSs; II) DSEP including investments in non-renewable DG (GT units); III) DSEP including investments in renewable DG (PV/WT units). In this case, the EV demand and conventional demand feed from a mix of renewable/non-renewable energy and the grid; IV) In this case, all investment alternatives mentioned in this work are considered (substations, circuits, renewable/non-renewable DG, EVCSs, and ESS); V) A sensitivity analysis with different ESS investment costs; VI) A sensitivity analysis was performed with different CO 2 emission rates; VII) In this case study, different candidate nodes for allocation of EVCSs are analyzed.
Case I: Investments in substations, circuits and EVCSs are considered. The solution found has a cost of 141.70 millions of USD with a total of 609.61 kTons of CO 2 emissions. The expansion plan includes the following actions: (1) First period: construction of substations located at nodes 23 and 24, allocation of EVCSs with the installation of 2 quick EV chargers and 2 super-fast EV chargers for the following nodes 4, 10, 13, and 14, installation of 2 quick EV chargers and 1 super-fast EV charger at node 15, and construction of 15 new circuits. (2) Second period: installation of 1 superfast EV charger at node 15.
Case II: In addition to investments in substations, circuits, and ECVSs, DSEP includes allocation of GT units. The total cost of the expansion plan is 124.28 millions of USD with a total CO 2 emission of 609.39 kTons. The investment plan provides the following actions: (1) First period: construction of substations located at nodes 23 and 24, allocation of EVCSs with the installation of 2 quick EV chargers and 2 super-fast EV chargers at nodes 4, 10, 13, and 14, installation of 2 quick EV chargers and 1 super-fast EV charger at node 15, allocation of two GT unit at nodes 3 and 11, and construction of 15 new circuits. (2) Second period: allocation of 1 superfast EV charger at node 15.
Case III: Investments in substations, circuits, allocation of renewable/non-renewable DG units, and allocation of EVCSs. The total cost of the expansion plan is 88.10 with a total CO 2 emission of 337.02 kTons. The expansion plan includes the following actions: (1) First period: construction of substations located at nodes 23 and 24, allocation of EVCSs with the installation of 2 quick EV chargers and 2 super-fast EV chargers at nodes 4, 10, and 13, installation of 2 quick EV chargers and 1 super-fast EV charger at nodes 14 and 15, installation of two GT units at nodes 3 and 5, six WT units at nodes 3,5,9,11,16, and 19, one hundred thirty-two PV units (thirteen at node 3, five at node 4, nine at node 6, nineteen at node 8, forty three at node 12, sixteen at node 13, twenty-five at node 14, and two at node 15), and construction of 16 new circuits. (2) Second period: allocation of 1 super-fast EV charger at nodes 14 and 15, and installation of twenty-five PV units (seventeen at node 8, eight at node 12) Case IV: Similar to Case III, but including the possibility to invest in ESS. Li-Ion batteries have been considered in the ESS with an energy to power E/P ratio of 4 hours, capacity of 250kW/1000kWh, and an associated cost equal to USD 343,000 ($ 271/kWh and $ 288/kW) [43]. Due to the high cost of investment in ESS, this case study chose not to install any battery bank. Thus, the results for this case are the same as for Case III.
The computational times to solve Cases I-IV are 108.50 s, 1050.03 s, 1347.36 s, and 3838.95 s, respectively. A summary of the main DSEP results for the first four cases is shown in Table 3. Note that in the first case (without investments in DG units) there was a higher total cost, a difference of approximately 12.29% and 37.83% related to Cases II and III/IV. Investment in renewable DG units to meet the EV demand and conventional demand benefits the system and contributes to total cost reduction (Case III/IV). Moreover, Case III, where there is a mix of investments in substation, circuits, EVCSs, and renewable/non-renewable DG units, obtained the more economical expansion plan.
To demonstrate the quality of the solutions found using the proposed model, the objective functions of Cases I-IV are compared to the original model (nonlinear formulation). It is important to highlight that both models obtained the same investment plans. Moreover, Table 4 shows that the linearization errors are negligible. Therefore, the proposed model provides optimal solutions within a good accuracy. Figure 5 shows the active power profile injected by the substation in all the periods and scenarios for cases I, II, and III (see Figure 5 (a), (b) and (c)). It can be highlighted that in Case III (Figure 5 (a)), with the support of the renewable DG units, there was a reduction in the energy supplied by the grid, thus lowering the total cost  of the expansion plan (see Table 3). Case III also obtained the best results in environmental aspects, with a significant reduction in CO 2 emissions of approximately 44.71% when compared to Case I, and 44.69% compared to Case II. The best expansion plan (Case III) is illustrated in Figure 6.
Case IV, unlike Case III, could carry out investments in ESSs. However, the expansion plan for this case chose not to install any ESS due to the high cost of the battery bank. Therefore, cases III and IV obtained the same results (see Table 3). Because of this, a sensitivity analysis is performed with different ESS prices (Case V).
Case V: Price of ESSs is expected to be annually reduced by 8% [44]. From that perspective, a sensitivity analysis with different ESS prices is performed to evaluate the corresponding changes of the expansion plan ( Table 5). The available battery banks have a capacity of 250 kW/1000 kWh. Note that when the investment price is lower ($96/kWh and $107.2/kW), the expansion plan chooses to install more battery banks and PV units. On the other hand, when the investment cost is higher ($271/kWh and $288/kW), the model does not install any ESS. Therefore, in the coming years with the reduction in the price of batteries, it will become more feasible to invest in that technology, contributing to deferring reinforcements in the EDS and integrating distributed energy resources.
Case VI: A sensitivity analysis was performed with different CO 2 emission rates. In this case, one more constraint that limits CO 2 emission rates is included in the model. The results of this analysis are summarized in Tqble 6. The VOLUME 10, 2022  first and the last solutions (1 and 6) represent the extreme solutions, with solution 1 presenting the lowest cost and maximum CO 2 emission value and solution 6 the opposite. The first solution represents Cases III and IV, 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 6 that, to reduce CO 2 emissions, investment in renewable DG units increases. Figure 7 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, due to the increase in investment in the renewable DG units, which contributes to the reduction of CO 2 emissions.
Case VII: Similar to [19], this work proposes a centralized approach for planning EVCSs. Thus, simulations were carried out varying the installation locations of EVCSs, according to the interests of different investors. Table 7 presents the expansion plan results for each location of EVCSs: solution 1 (4,10,13,14,15), solution 2 (3,6,8,12,19), solution 3 (5,9,11,16,18). The set of nodes for allocation of EVCSs used in solution 1 was the same used for the previous case studies. The results of solution 1 are very close to solutions 2 and 3, with a difference of 0.00135% and 0.0027%, respectively. The expansion plan of solutions 1 and 2 are identical. The only difference between solution 3 and the other solutions is the installation of a PV unit.
The results of the cases highlight the advantage of an integrated planning including different expansion alternatives. Moreover, the centralized approach for planning EVCSs can  serve as an auxiliary tool so that different agents, investors, network owners can define, according to their needs, the locations where the EVCSs will be installed, to maximize the social benefit.

V. CONCLUSION
A mixed-integer linear programming (MILP) model for the distribution system expansion planning (DSEP) problem has been proposed. Uncertainties related to conventional demand, EV demand, and renewable generation along with investments in distributed generation (DG) units and energy storage systems (ESSs) were considered.
In contrast to previous works, the proposed model includes the simultaneous optimization of investments in substations, circuits, and distributed energy resources, including environmental aspects (CO 2 emissions). Moreover, this paper presented an algorithm to estimate EV charging demand. The EV demand uncertainty model showed that more than 90% of the EVs were charged when they had an initial state of charge of up to 50%.
Tests carried out in a 24-node system suggest that considering a simultaneous investment plan in substation, circuits, EV charging stations (EVCSs), DG units, and ESSs within the DSEP results in expansion plans at lower costs since the operational costs of the system is reduced. The coordinated investment plan proved to be more promising, contributing to making the system less dependent on the support of substations. Besides, this expansion plan brings environmental benefits contributing to the reduction of CO 2 emissions. Another important aspect is that some expansion plan do not propose installation of battery banks due to the high cost of this technology; however, in the coming years with the reduction in the price of batteries, it will become more feasible to invest in that technology, contributing to deferring reinforcements in the EDS and integrating distributed energy resources.
Future work would address the robust optimization that has been little explored in the specialized literature compared to stochastic programming. Investment variable for a PV unit at node u and period p. P ij,a,ω,p Active power flow through circuit ij for conductor a in scenario ω and period p. P ij,ω,p Active power flow through circuit ij in scenario ω and period p. P ESC b,ω,p Active power stored of an ESS at node b, scenario ω, and period p. P ESD b,ω,p Active power supplied by the ESS at node b, scenario ω, and period p. P gt f ,ω,p Active power injected by the GT units at node f , 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,a,ω,p Reactive power flow through circuit ij for conductor a in scenario ω and period p. Q ij,ω,p Reactive power flow through circuit ij in scenario ω and period p. Q Operational variable related to the backward/forward direction of circuit ij and period p. z ij,a,p Operational variable related to circuit ij using conductor type a, at period p.