Distributionally Robust Distributed Generation Hosting Capacity Assessment in Distribution Systems

Uncertainties associated with the loads and the output power of distributed generations create challenges in quantifying the integration limits of distributed generations in distribution networks, i.e., hosting capacity. To address this, we propose a distributionally robust optimizationbased method to determine the hosting capacity considering the voltage rise, thermal capacity of the feeders and short circuit level constraints. In the proposed method, the uncertain variables are modeled as stochastic variables following ambiguous distributions defined based on the historical data. The distributionally robust optimization model guarantees that the probability of the constraint violation does not exceed a given risk level, which can control robustness of the solution. To solve the distributionally robust optimization model of the hosting capacity, we reformulated it as a joint chance constrained problem, which is solved using the sample average approximation technique. To demonstrate the efficacy of the proposed method, a modified IEEE 33-bus distribution system is used as the test-bed. Simulation results demonstrate how the sample size of historical data affects the hosting capacity. Furthermore, using the proposed method, the impact of electric vehicles aggregated demand and charging stations are investigated on the hosting capacity of different distributed generation technologies.


Introduction
In recent years, integration of renewable-based distributed generations (DGs) has significantly increased in distribution networks.This is mainly due to governmental support, advancement in technology and economy of scale.Although financially DGs are one of the most viable options for end users, technically, high penetration of them may introduce many issues such as over-voltage, overloading of the feeders [1], degradation of power quality and even higher losses in some situations [2].These issues limit the integration of DGs in distribution systems.Thus, a hosting capacity (HC) assessment is required prior to the practical deployment of DGs.By doing so, distribution system operators (DSOs) can identify the maximum DG capacity that can be accommodated in their networks and make optimal decisions regarding the placement and sizing of DGs.Furthermore, DSOs can plan the future network augmentation and operation strategies to increase the HC of their system.
The HC problem has been studied considering different constraints [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].In [3,4], optimal power flow methods were proposed to solve the HC problem.Analytical [5,6], heuristic [7,8], and sequential power flow-based methods [9,10] were also developed to determine the HC.Nonetheless, the uncertainties related to loads, output power of DGs and the impact of electric vehicles (EVs) were not considered in those studies.The sensitivity of the HC to uncertainties is assessed in [16] using an optimization-based Monte Carlo framework.That study, however, used the time series simulations to address the uncertainties, which can be time consuming for analyzing large numbers of future grid scenarios.To address the uncertainties associated with DGs and loads, an optimization-based framework was presented in [11], which probabilistically assesses the HC.The study in [12] presented a framework to assess the effects of EVs and photovoltaic (PV) systems on the HC, separately.
Two dominant approaches that have been used to model the uncertainties in the HC assessment are (i) the robust optimization (RO) [2,13] and (ii) the stochastic optimization (SO) [14].In RO approach, the probability distribution functions (PDFs) of uncertainties are not required, and the only necessary data are the boundaries of the uncertain variables.Furthermore, the decisions are made based on the worst-case scenarios in the predefined uncertainty intervals, which result in a conservative solution.
In SO-based approach, the uncertainties are considered as random variables with given PDFs.Since obtaining the exact PDFs of uncertain variables is impossible, the robustness of the obtained solution cannot be guaranteed.
In practice, historical data can present more information than the boundaries of uncertain variables.However, historical data can be used to obtain an empirical PDF, not the exact PDF for the uncertain variable.To cope with such cases, distributionally robust optimization (DRO) has been recently developed.In DRO, it is supposed that the exact PDFs of uncertain variables are not available.However, it is supposed that the PDFs of uncertainties are in a confidence set.The confidence set can be constructed using the distributional information of uncertain variables, which is derived from the historical data.For instance, the confidence set can consist of all the PDFs with common mean and covariance matrix [15].In recent years, DRO have been applied to some power system optimization problems including unit commitment [17], and reserve scheduling [18].However, those studies often considered the first-and second-order moments of historical data to build the confidence set.Furthermore, all those studies neglected the empirical distribution of uncertainties.
In this paper, we propose a DRO-based method for evaluating the HC in distribution networks.The main contributions of this study are summarized as follows: 1. We propose a DRO-based method to evaluate the HC of distribution networks considering uncertainties associated with loads and DGs' output powers.As the statistic information (i.e., the first-and second-order moments) and the empirical distribution of uncertain variables are exploited in the assessment process, the results of the proposed method are practical and less conservative than those of the RO-based method.To the best of our knowledge, there is no such DRO application in power system that exploit both the statistic information and the empirical distribution simultaneously.2. In the proposed DRO-based method, a risk level is defined to control the robustness, enabling a trade-off between robustness and conservativeness.Unlike the robust optimization, tuning the risk level in the proposed method has a precise physical meaning, i.e., guaranteeing that the probability of operational constraints violation does not exceed a given risk threshold.3. The aggregated EVs and charging stations' loads are modeled in the HC problem and their effects on the HC have been assessed.To the best of our knowledge, there is no such study that assesses the effect of aggregated EVs and charging stations' loads on the HC of various types of DGs, i.e., PV, wind, and biomass.
The proposed method is examined on a modified IEEE 33-bus system and the sensitivity of the HC to the DG technology, the EVs aggregate load, the load of EVs charging station and historical data has been assessed.
The reminder of this paper is organized as follows: Section 2 formulates the problem.Section 3 presents the uncertainty modeling.In Section 4, the solution methodology is described.First, the DRO model is transformed into a joint chance constraint (JCC) problem.Then, a solution for the corresponding JCC is introduced.Section 5 presents the numerical results and discussion.Finally, the major conclusions are summarized in Section 6.

Mathematical Modeling
The HC is defined as the maximum DG capacity that can be installed in a system without violating its technical constraints irrespective of DGs' locations.This implies that DGs' locations should not be defined as independent variables in the optimization model.This, however, raises a question on how to address the uncertainty associated with DGs' locations in the HC problem.This can be tackled using the HC calculation framework shown in Figure 1.As it can be seen, it is required to generate a large number of location combinations for DGs.Then, the maximum DG capacity for different location scenarios can be identified using the mathematical model.Finally, the minimum of the identified DG capacities for all location scenarios is defined as the actual HC.The focus of this paper is to develop a mathematical model to find the maximum DG capacity for any sets of DG locations.The objective function of this model is presented as follows: where, DG is the set of buses that have DGs, and Cap g j denotes the installed DG capacity at bus (j).The system operation constraints are presented in Sections 2.1 and 2.2.

Distribution System Model
Consider a radial network, where N = {0, . . ., N Bus } denotes the set of buses.Let B represents the set of all branches and (i, j) or i −→ j as a branch from bus (i) to bus (j) in the set B. Let Γ denotes the set of time periods.For every bus i ∈ N and for all t ∈ Γ, let V i,t denotes the complex voltage, where represents the load at bus (i) at time (t).For every line (i, j) ∈ B and for all t ∈ Γ, z ij = r ij + ix ij denotes the complex impedance, S ij,t = P ij,t + iQ ij,t denotes the sending-end complex power from bus (i) to bus (j), and I ij,t defines the sending-end complex current from bus (i) to bus (j).Let s g i,t = p g i,t + iq g i,t denotes the generation complex power at bus (i) at time (t), and s EV i,t = p EV i,t + iq EV i,t represents the EV aggregated complex demand at bus (i) at time (t).For every bus i ∈ DG, η g i,t represents the efficiency factor at time period (t).Assume T j as the subtree rooted at bus (j) (including (j)).We use k ∈ T j to refer to bus (k) in subtree T j and (k, l) ∈ T j to refer to line (k, l) in subtree T j .Let P j defines the set of lines on the path from substation to bus (j).The following equation can be written based on the Ohm's law: taking the magnitude squared from (2) (by taking the product of each side of (2) with its conjugate) results in: and according to the nodal power balance equations, we have: where, φ i,t is the power factor angle at bus (i) at time period (t).Considering the unique path from the substation to each bus, the voltage at bus (j) can be calculated as follows: The quadratic terms in ( 4), ( 5) and ( 7) are the sources of non-convexity in the optimization problem (1).These terms are small and ignoring them in the HC problem results in a sufficiently accurate linear model [15,19].Therefore, equations ( 4), ( 5) and ( 7) can be written as:

Technical Constraints
The following technical constraints are considered for determining the HC: where, υ i and υ i are the minimum and maximum voltage limits at bus (i), respectively.

Thermal Capacity Constraints
The apparent power flow of lines and substation transformer is limited by a higher bound as follows: where, |S ij | is the maximum apparent power of line (i, j).Equation ( 12) is a quadratic nonlinear constraint, which needs to be linearized to simplify the final problem formulation.For doing so, we rotate the tangent (13) around the original circular constraint (12) using the counter clockwise rotation matrix (14).
where, θ is the rotation angle.Applying ( 14) to ( 13) results in: One of the important characteristics of distribution systems is SCL, which is defined as the maximum acceptable fault level.A basic requirement for connecting DGs to the networks is that the SCL in presence of DGs should remain below the designed SCL [20].It is shown in [21] that the SCL can be modeled as follows: where, a j is the dependency of SCL to the capacity of the DG at bus (j), SCL sub is the SCL at substation without DGs, and SCL Rated is the rated SCL at substation.It is worth mentioning that a j depends on the structure of the system and the DG type; e.g., a j is very small for PVs [22].

Deterministic Problem Formulation Summary
The deterministic HC problem is formulated as follows: In practice, the output power of DGs and the loads are uncertain variables and modeling them as deterministic variables is not realistic.Therefore, the HC optimization model should include such uncertainties, as explained in the next subsection.

Distributionally Robust Optimization (DRO) HC Model
The uncertainties of the DG output (p g i,t ), load (p L i,t ), and EVs aggregated demand (p EV i,t ) at bus (i) are modeled as follows: where,  18) in ( 6), ( 8)- (10), and denoting ξ = {ξ g i,t , ξ L i,t , ξ EV i,t } to the vector of uncertain variables, the power flow state variables (P ij,t , Q ij,t , υ j,t ) can be written as functions of Cap g and ξ.Thus, all constraints can be expressed as follows: where, C eq and C are the sets of equality and inequality constraints for the HC problem.Hence, the uncertain HC model is presented as follows: To immunize the HC model ( 21) against the uncertainty vector, ξ, we adapt a distributionally robust approach.In this approach, a risk level is defined to adjust the conservatism of the solution.From a modeling prospective, the distributionally robust HC problem is defined as follows: where, τ is the risk level, D φ is the confidence set, and T (Cap The next step is modeling the uncertainty vector, ξ, and the confidence set, D φ , which is detailed in the next section.

Uncertainty Modeling
In this section, the uncertainties associated with the outputs of DGs, loads and aggregated demands of EVs as well as the confidence set are modeled based on the historical data.Herein, the generation technologies are PV, wind, and biomass.

Uncertainty Modeling of PV Generation
The stochastic variations of PVs from their predicted output values follow a Beta distribution [23].This distribution, defined by two shape parameters, α and β, enables us to represent the prediction error of a predicted output power, pg i,t , with the normalized predicted output power and a standard deviation, which varies with ηg i,t .The Beta function for modeling the occurrence of an output power, x, if a prediction value, pg i,t , has been forecast, is as follows: The shape parameters α and β are related to the normalized predicted output power and variance as follows:

Uncertainty Modeling of Wind Generation
A Beta function is justified to model the prediction error of the wind power [24].The relation between the parameters of Beta distribution (i.e., α and β) and the normalized predicted output power and variance are presented in (24) and (25).

Uncertainty Modeling of Biomass Generation
The biomass DGs are considered as firm generation.Thus, the output powers of such DGs are considered constant at their rated capacities [25].

Uncertainty Modeling of Load
Load uncertainty is usually modeled by a normal distribution, in which, the mean value is the forecast load and the standard deviation is set to be 2% of the mean value [26].

Uncertainty Modeling of EV Demand
The demand of EVs depends on the number of EVs, their state of charge (SOC), charging start time and its duration.These variables are uncertain, so the overall charging demand of EVs is uncertain as well.Considering the charging behavior of EVs, the overall demand can be categorized as follows:

Overall Charging Demand of EVs in a Local Residential Community
It was shown in [27,28] that the aggregated demand of EVs in a residential community follows a normal distribution at each hour, and the average and standard deviation of the distribution depend on the number of EVs and the transportation data of the area.

Overall Charging Demand of an EV Charging Station
The demand of a charging station depends on the number of EVs that arrive in different time intervals, the duration of a charging process and the charging power profile [27,29].It was shown in [27] that the demand uncertainty of a charging station can be modeled by using a Weibull distribution.

Modeling the Confidence Set
A challenging difficulty of handling the uncertainties in HC optimization is the accessibility to the exact PDFs.Since the historical data is limited or there might be not much trust in it, assuming the perfect knowledge about the PDFs of uncertainties is unrealistic.To address this, an ambiguity set of distributions can be used.A common method to create an ambiguity set exploiting the empirical distributions is φ−divergence, which is defined as follows: where, f and f 0 are the actual and estimated density functions, respectively; ξ ∈ R K represents a K-dimensional random vector defined on a probability space, Ω, and φ : R −→ R is a convex function on R + .We refer the readers to [30] for more details on characteristics of φ−divergences.One of the most common used members of the φ−divergence families is Variation Distance, which is defined as: Based on the φ−divergences, a confidence set can be built as follows: where, M + represents the set of all cumulative density functions (CDFs), and ψ denotes the risk-aversion level.The higher the risk-aversion level (ψ) is, the bigger the ambiguity set and the more conservative the result of the optimization model would be.However, as compared to uncertainty sets in the RO or the moment-based ambiguity sets, the confidence set, D φ , can more accurately depict the profile of PDFs, and so provides a less conservative result.

Solution Methodology
To solve the HC problem ( 22), we convert it to an equivalent JCC optimization by employing a theory presented in [30].Then, the sample average approximation is used to solve the equivalent JCC [31,32].Before presenting the equivalent JCC problem, we need to review the definition of conjugate duality.For a given function g : R −→ R, the conjugate g * : R −→ R ∪ {+∞} is defined as follows: g
Therefore, if we have a function τ(τ + , φ, ψ) that maps the τ + to the original risk level, τ, for the given tolerance, ψ, and φ-divergence, it can be stated that if (30) is satisfied by some Cap g , then inf [T (Cap g , ξ)] is equal to the optimal value of the following optimization problem: minimize r,s,tt>0 λr + P 0 s s.t.
µr + P 0 s + (1 − P 0 )tt = 1 where, we let P 0 = P 0 (T (Cap g , ξ)) for notation brevity, and this follows that the mapping function τ(τ + , φ, ψ) can be defined as the optimal objective value of ( 31) with setting P 0 to be 1 − τ + .Thus, the new risk level, τ + , can be evaluated using a bisection line search algorithm, as given in Algorithm 1.

Sample Average Approximation (SAA)
Due to the intractability of JCC programs, we use the sample approximation in which the distribution of the random ξ vector is replaced with the empirical distribution obtained from the historical data.In the SAA method, an approximation problem, which is based on an independent Monte Carlo samples of the uncertainty vector, is solved.To simplify the presentation, we replace the constraints set T (Cap g , ξ) with a real-valued function G as follows: Thus, the JCC HC problem can be represented as: where, X where, X ⊂ R n represents a deterministic feasible region.Let ξ 1 , ξ 2 , • • • , ξ N be N independent Monte Carlo samples of the uncertainty vector ξ.Then, for γ ∈ [1, 0), the sample approximation problem is defined as follows: where, where, I(•) is the indicator function, which takes value one when the argument is true and zero otherwise; P τ + , and P N γ are referred as the true and SAA problems, respectively.It can be proved has an optimal solution, then: where, If γ > τ + , then: Using M set of N independent samples of ξ, the lower bound that is valid with confidence 1 − δ can be obtained using Theory 1.
Theory 1.If δ ∈ (0, 1), γ ∈ [0, 1), N, H, and M be positive integer variables such that H ≤ M and, where, ẑN γ, [H] is the optimal value of the set H such that ẑN γ, [1] [31].Furthermore, it is proven that if γ < τ + , the feasible region of the P N γ will be a subset of the feasible region of the P ξ) is a Caratheodory function, and ( 5) there is an optimal solution Cap g of the true problem such that for all ε > 0 there is The sample approximation problem (35) is reformulated as a mix-integer problem as follows: where, Z j is a binary variable and M j is a large positive number.

Numerical Results
In this section, simulations are carried out to assess the performance of the proposed methodology.Firstly, the modified IEEE 33-bus network is presented in Section 5.1.Then, the efficacy of the proposed method is evaluated in Section 5.2.

Test System
The performance of the proposed method is evaluated on the modified IEEE 33-bus system, shown in Figure 2. The bus 1 is connected to the 33 kV grid using a Dy transformer with 8.5% reactance.The SCL and the rated SCL at the 33kV side are 200 MVA and 250 MVA, respectively.The system nominal voltage is 12.66 kV.There are four candidate sites for DGs, i.e., DG = {15, 22, 25, 32}.There are two charging stations at buses 26 and 30.Additionally, the charging demand of EVs in the residential area is aggregated at buses 6 and 13.The detail description of the test system is presented in [33].To have a real situation, three load models (i.e., industrial, commercial, and residential) are considered, as shown in Figure 2. The PV and wind generation profiles, the aggregated EVs charging demand in residential area, the EV charging station demand curve and the load profiles are provided in Figure 3.The PV output profile and the residential, commercial, and industrial load profiles are general profiles that are derived from [13].The aggregated EVs demand profile is derived from [34,35], which are based on the National Household Travel Survey (NHTS).As for the demand of charging stations, the obtained profile using the NHTS data provided in [29] is used.Since the NHTS data presents the driving behavior of American public, the wind generation profile is derived from the data provided for the continental United States in [36].The presented curves in Figure 3 are normalized multipliers.Please note that basic values of DGs and loads are the capacity of DGs and the nominal loads, respectively.The basic value for charging stations' demand is 224 kW and for aggregated demand of residential EV is 285 kW.

Simulation Results on the Modified IEEE 33-Bus System
The proposed DRO-based method is examined on the IEEE 33-bus system for the four cases as follows: Figure 4 shows the estimated HC of different generation technologies for different risk levels.Observe that the HC grows by increasing the risk level, which indicates that looser security requirement leads to a higher HC.Furthermore, as can be seen in Figure 4, the HC curve slope depends on the DG technology.For instance, accepting 20% security risk, increases the HC by 2.48%, 20.3% and 1.17% for PV, wind, and biomass technologies, respectively.This implies that accepting a higher risk to increase the HC of the system is not a good idea if the DGs technologies are PV and biomass.An important factor that affects the estimated HC of the system is the existence of the historical data.Intuitively, as the sample size of historical data increases, the empirical distribution becomes closer to the actual distribution.Thus, the conservativeness of the HC estimation decreases.In DRO-based method, ψ is the parameter to present the value of data.The higher the number of samples, the more accurate the empirical distribution, the smaller the risk-aversion level, ψ, and the smaller the confidence set, D φ .Figure 5 demonstrates the effect of risk-aversion level (hence historical data) on the HC for different technologies.As it can be seen in Figure 5, the shortage of historical data increases the conservativeness of the HC estimation.For instance, the shortage of data decreases the estimated HC for PV, wind, and biomass technologies by 3.09%, 14.52% and 2.06%, respectively.Observe that the value of data in the test system for wind technology is much higher than that of PV or biomass technologies.Furthermore, as the slope of HC curves for small values of risk-aversion level (i.e., ψ ∈ [0%, 5%]) is small, increasing the size of historical data to get a more accurate empirical distribution does not have a tangible effect on the estimated HC.As shown in Figure 5, 5% error in empirical distributions only decreases the HC by 0.21%, 1.41%, and 0.11% for PV, wind, and biomass technologies, respectively.A prospective technology that may affect the HC of a system is EVs. Figure 6 indicates the effect of EVs' aggregated demand on the HC.Observe that increasing the peak value of EVs' aggregated demand to 640 kW can only increase the HC of the system by 0.29%, 0.25%, and 0.13% for PV, wind and biomass technologies, respectively.Thus, the aggregated impact of EV loads on the HC of the system is not tangible.This is because the HC reaches its maximum value during time period t = 11 for PV, t = 5 for wind, t = 5 for biomass and t = 9 for combined generation.We will call these periods as the critical time periods.During critical time periods, the normalized aggregated demand of residential EVs is below 10%, which means that increasing the EVs cannot effectively increase the load.To demonstrate that the concluded result does not depend on the aggregation points of EVs, the sensitivity of the HC to the location of aggregated EVs demand is also presented using error bars in Figure 6.Observe that changing the location of aggregated demand of residential EVs deviates the obtained HC by 0.068%, 0.072%, and 0.032% for PV, wind, and biomass technologies, respectively.In other words, irrespective of locations that demands of residential EVs are aggregated, increasing the peak value of aggregated EVs demand does not increase the HC, effectively.
Finally, the effect of charging stations' demands on the HC is assessed.To do so, the peak value of charging stations' demands is varied from 0 kW to 480 kW, as shown in Figure 7. Observe that increasing the demand of charging stations can linearly increase the HC of the system.However, the rate of increase in the HC depends on the generation technology.Increasing the peak value of charging stations' demands from 0 to 480 kW grows the HC by 1.81%, 0.15%, and 0.1% for PV, wind, and biomass technologies, respectively.This is because the critical time period depends on the DGs technologies.The charging station demand profile at critical time period for PV, wind, biomass, and combined generations are 58.26%,1.83%, 1.83%, and 27.06%, respectively.Since the demand profile of charging stations during the critical time period of PV generation is higher than that of other technologies, the charging stations' demands have a higher effect on the HC of the system for PV generation.

Conclusions
In this paper, a DRO-based method is proposed to assess the HC of distribution networks.In the proposed method, the empirical distributions of uncertain variables are used to define a confidence set.The more accurate the empirical distribution is, the smaller the confidence set and the less conservative the obtained result would be.The effectiveness of the proposed method was examined on the modified IEEE 33-bus system and the effects of risk level and size of historical data on the HC of the system were assessed.It was demonstrated that the proposed method yields in a less conservative solution by increasing the risk level.However, the increase in the HC by accepting a higher risk level depends on the DG technology.For instance, accepting 20% risk level can increase the HC of the system by 2.48%, 20.3% and 1.17% for PV, wind and biomass technologies, respectively.Thus, accepting a higher risk level to increase the HC of the test system may not be a good idea if the DGs technologies are PV and biomass.As for the historical data, it was observed that shortage of data can exponentially increase the conservativeness of the estimated HC.For instance, the shortage of historical data decreases the estimated HC for wind technology up to 14.52%.
The proposed method is also used to assess the effects of the aggregated residential EVs and charging stations' demand on the HC for different DG technologies.Following are the conclusions derived from the simulation results: • Although aggregated demand of residential EVs increases the peak load of the distribution system, it does not affect the HC significantly.This is because the HC of the system reaches its maximum value during the time periods that the aggregated demand of EVs is below 10% of its peak value.
• The effect of charging stations' demand on the HC depends on the DG technology.If DGs are PV units, increasing the charging stations' demand to its maximum increases the HC up to 1.81%.Nevertheless, it cannot effectively increase the HC of the system if DGs are wind and biomass generators.
A future extension would be considering the uncertainty associated with the location of DGs.To do so, the developed mathematical model should be integrated in a Monte Carlo-based framework.
is held for φ and ψ.However, we still need to know how to obtain the τ + and the mapping function.It can be shown that the worst-case probability bound of inf f ∈D φ Pr ξ∼ f

Figure 3 .
Figure 3.The daily PV, wind, aggregated EVs, charging stations and load normalized profiles.

•
PV-HC: All DGs are PV systems.• Wind-HC: All DGs are wind generators.• Biomass-HC: All DGs are biomass generators.• Combined-HC: The DGs at buses 15 and 32 are PV units, and the DGs at buses 22 and 25 are wind and biomass generators, respectively.

Figure 4 .
Figure 4.The HC obtained via the proposed DRO-HC method with different risk levels for different technologies including: (a) PV; (b) wind; (c) biomass; and (d) combination of PV, wind and biomass.

Figure 5 .
Figure 5.The value of data for the HC obtained via the proposed DRO-HC method with τ = 10% for different technologies including: (a) PV; (b) wind; (c) biomass; and (d) combination of PV, wind and biomass.

Figure 6 .
Figure 6.The effect of aggregate demand of residential EVs on the HC obtained via the proposed DRO-HC method with τ = 1% for different technologies including: (a) PV; (b) wind; (c) biomass; and (d) combination of PV, wind and biomass.

Figure 7 .
Figure 7.The effect of charging stations' demands on the HC obtained via the proposed DRO-HC method with τ = 1% for different technologies including: (a) PV; (b) wind; (c) biomass; and (d) combination of PV, wind and biomass.
and ξ EV i,t model the prediction error of the DG output, load and aggregated demand of j , substituting (