Chance-Constrained Optimization Based PV Hosting Capacity Calculation Using General Polynomial Chaos

Increased penetration of renewable resources and new loads have increased the uncertainty levels in low voltage distribution systems (LVDS). This requires considering LVDS planning, such as computing photovoltaics (PV) hosting capacity (HC), as a stochastic problem. Traditionally, PV HC is computed using the iterative Monte Carlo method, which requires solving the power flow equations thousands of times. This paper proposes a chance-constrained optimization-based hosting capacity calculation technique, which eliminates the necessity of repetitive solving of power flow equations. The intrusive general polynomial chaos expansion is used to translate the input uncertainties defined by their probability density function to the hosting capacity of the network without the necessity of sampling, linearizing the power flow equations, or applying any relaxations. Chance constraints are applied for nodal voltages and thermal overload as per the norms, where the system can be congested for a certain time without affecting the power quality. Numerical illustrations show the computational time to compute the stochastic PV hosting capacity of a real large scale LVDS, and for the majority of the feeders, it falls within 100 sec. Furthermore, the estimated PV HC using this stochastic optimal power flow was, on average, 20% higher than its deterministic counterpart.

Abstract-Increased penetration of renewable resources and new loads have increased the uncertainty levels in low voltage distribution systems (LVDS).This requires considering LVDS planning, such as computing photovoltaics (PV) hosting capacity (HC), as a stochastic problem.Traditionally, PV HC is computed using the iterative Monte Carlo method, which requires solving the power flow equations thousands of times.This paper proposes a chanceconstrained optimization-based hosting capacity calculation technique, which eliminates the necessity of repetitive solving of power flow equations.The intrusive general polynomial chaos expansion is used to translate the input uncertainties defined by their probability density function to the hosting capacity of the network without the necessity of sampling, linearizing the power flow equations, or applying any relaxations.Chance constraints are applied for nodal voltages and thermal overload as per the norms, where the system can be congested for a certain time without affecting the power quality.Numerical illustrations show the computational time to compute the stochastic PV hosting capacity of a real large scale LVDS, and for the majority of the feeders, it falls within 100 sec.Furthermore, the estimated PV HC using this stochastic optimal power flow was, on average, 20% higher than its deterministic counterpart.
Index Terms-AC optimal power flow, hosting capacity, low voltage distribution system, general polynomial chaos, chanceconstrained.

E
Expectation operator.Set of the coefficient of orthogonal polynomials bases.lji ∈ T l← Tuple set of branch-l from bus-j to bus-i l ∈ L Set of branches in the feeder.lij ∈ T l→ Tuple set of branch-l from bus-i to bus-j p ∈ P Set of installed PV devices.pi ∈ T p Tuple set of PV installation-p connected to bus-i P rt  Set of available PV ratings.s ∈ S Set of source generator.si ∈ T s Tuple set of source-s and bus-i T l

M
Union Set branch tuples: T l← ∪ T l→ u ∈ U Set of units, load, PV or source.Series impedance.

I. INTRODUCTION
A. Background and Motivation D ETERMINING the hosting capacity (HC) of a Low Volt- age Distribution System (LVDS) is a planning problem where the additional generation and/or load is determined that can be connected to the system without violating its performance limits.Given the increased connection of weather-dependent photovoltaics (PV) generation and new types of load, e.g., electric vehicles and heat pumps, the level of uncertainty in LVDS has increased, and consequently determining its HC becomes a multidimensional stochastic problem.Furthermore, the performance limits of a LVDS are not absolute as they may be violated momentarily [1].Consequently, for LVDS, a risk-based stochastic HC model enables improved decision-making [2].
The uncertainties in LVDS are categorized into different types.In [3], three types of uncertainties are proposed: 1) planning, e.g., location, size and type of PV installations, 2) feeder, e.g., feeder length, consumer phase, and 3) operational, e.g., load and PV generation.
The former two types are generally described by discrete random variables, whereas the latter is described by continuous random variables.Despite the difference in these variables, these uncertainties are often sampled together in an iterative MCbased stochastic HC method for computational tractability.The analytical alternative, such as non-intrusive polynomial chaos expansion, can be used as an alternative to reduce the computation time of the Monte Carlo (MC) based HC computation if the uncertainties can be represented by continuous distributions [4].However, planning uncertainties are represented as unknowns rather than in terms of a Probability Density Function (PDF) [5].Therefore, the risks presented by the operational variables are different from the planning risk.A decoupled HC method can incorporate these two types of risks in the same model [3].However, such methods usually use a brute force approach where all planning scenarios are evaluated individually [3], [5].Consequently, there is a need for a novel method that allows the incorporation of continuous random variables, i.e., the operational uncertainties and discrete planning and feeder variables, in a single model, in a computationally tractable manner.

B. Literature Review
In [6], [7], existing PV HC methods for LVDS are reviewed, and both deemed considering uncertainties in calculating HC to be important.The available deterministic and stochastic PV HC calculation methods in the literature are evaluated and benchmarked in [2].Stochastic methods consider planning and operational uncertainties, while deterministic methods only take into account one possible operational scenario, e.g., peak generation PV with low load.The HC obtained from deterministic methods are based on a single snapshot of a number of possible combinations of stochastic variables and, therefore, provide a rough estimate of the HC.In contrast, the stochastic HC calculation provides a more comprehensive estimation of the HC.Another observation from [2] is that every feeder in a network had a different HC, and therefore there is no single number that defines the HC of all feeders.The underlying assumptions of every HC calculation method, e.g., size of PV, their allocation, and the increment rate of PV size and number, influence the HC of a feeder.Apart from the way of allocating PV, the assumption on the stochastic and operational constraints applied to the voltage and current directly influences the calculated HC.An important takeaway was the use of stochastic constraints, also known as chance-constraints (CC), in HC calculation which is deemed necessary to prevent the calculation of HC for the worst-case [2], [8].
The stochastic HC methods in literature usually rely on MC based methods [2], [9].Due to the high computation time associated with MC methods, these methods are usually demonstrated using representative feeders [10], reducing the input source of uncertainties [11] and checking only a few operational parameters [12].The calculated PV HC is a suboptimal solution, that is, the one with the highest total PV among the considered PV deployment scenarios.The Optimal Power Flow (OPF) and Stochastic Optimal Power Flow (SOPF) based methods to calculate PV HC in an optimization framework are scarce and in most approaches, many assumptions are made to linearize or relax the problem.In [13], [14], optimisation based methods were proposed for medium voltage (MV) feeders with a fixed location of PV.The authors proposed a worst-case analysis using robust optimization in [13] while linearizing the power flow (PF) equations and considering the load of consumers as the only source of uncertainty.The method improved the computational speed from 84 hrs in the traditional iterative method to 21 s.In [14], a two-stage optimization method was proposed that considers uncertainties in loads and distributed generation, where the operating point is approximated first before calculating the optimal PV HC.The method also used MC to formulate the scenarios of PV location and applied a worst-case analysis.In [15], an orthogonal-polynomial based system was proposed to calculate the PV HC when continuous probability density functions represent the uncertainties.This method resulted in the risk associated with a particular scenario of PV penetration without any optimization involved.In [16], an analytic-probabilistic approach to distributed generation HC evaluation for radial feeders using an analytical method called Herman Beta Extend transform was demonstrated.This is a moment-based transformation of the power flow equations, where the power flow equation are linearized using Taylor's expansion.The method works only when a Beta distribution represents the uncertainty of the inputs.The proposed method also assumes that the output voltage and current follow a Beta distribution, which is not necessarily true due to the nonlinear nature of the power flow equations.However, the method reduced the computational time compared to MC simulations, from days to a few minutes.A scenario-based CC-OPF was presented in [17], where the authors applied chance constraints in the curtailment risk, but an increasing number of scenarios had a huge toll on computation time.The authors of [8] used an actual LV feeder from the U.K. to demonstrate a stochastic OPF based HC of PV with storage, which calculated the minimum and maximum HC of each feeder.However, the method relies on linearizing the power flow equations, uses MC samples of load and PV generation, uses worst-case voltage and current limits instead of chance constraints, and assumes PV size as a constant rather than a variable.The only variable optimized was the number and the location of PV.In [18], a method for calculating stochastic HC for LV networks of the U.K. is proposed as an optimization problem.A bisection algorithm is used for obtaining optimal PV size and location where the maximum voltage in the feeder is equal to the voltage limit.This method uses linearized power flow and does not consider uncertainties in PV generation and consumer load.However, these methods clearly show the benefits of moving from the scenario-based method or iterative MC based methods to SOPF based calculations.The benefits are: 1) reduction in computational time from days to seconds, 2) similar or better accuracy, 3) easy scale-up for a large service area, and 4) minimum input and output data processing requirement.SOPF based methods also makes comparing different PV scenarios faster and easier compared to scenario based MC methods.So far, the methods proposed in literature do not include a full AC formulation but rely on relaxations of the power flow equations.The general Polynomial Chaos (gPC) expansion based CC-OPF introduced in [19] uses a full Rectangular Power-Voltage (ACR) formulation of the power flow.In [20], the authors use the Rectangular Current-Voltage (IVR) formulation of gPC-CC-OPF with auxiliary variables to make the problem scalable compared to the ACR formulation.It is to be noted that ACR is shown to have the best computational performance in solving OPF [21].However, the possibility of zero voltage in gPC-based stochastic reformulation might lead to a disconnected search path and convergence issues while using ACR formulation [22].Hence, a full IVR formulation of power flow is considered a suitable choice in this work.

C. Contributions and Outline
This paper proposes a SOPF method to calculate the PV HC of an LVDS as an application based on gPC-CC-OPF, considering  operational uncertainties.As such, computing the PV HC of a LVDS is determined as a single-shot problem without the need for any sampling, relaxation, or linearization.The method is demonstrated using a single-phase equivalent of a radial network, where load and irradiance are the sources of uncertainty.
To the best of the authors' knowledge, this is the first work to propose a SOPF with chance-constraints to calculate the PV HC of a LVDS without sampling, relaxation or linearization.To this end, the main contributions of the paper are: 1) a gPC-CC-OPF-based method to calculate the PV HC of the single-phase equivalent of a LVDS, considering PV and load uncertainties, and 2) an empirical method to tune the moment-based reformulation of the probabilistic chance constraints for non-Gaussian uncertainty.It is to be noted that the HC calculated in this paper is the best planning scenario, which is always with three-phase PV with balanced generation.Hence, it can be solved using balanced power flow equations.
The remainder of the paper is organized as follows: Section II presents a deterministic OPF-based method to determine the PV HC.Subsequently, Section III presents a stochastic gPC-CC-OPF-based method to determine the PV HC.Section IV presents a numerical illustration where results from the stochastic HC method are compared with its deterministic counterpart in a realistic LVDS setting.In Section V, a case study regarding computational effort and accuracy on calculated HC in a real LVDS is presented.Lastly, Section VI concludes the paper.

II. DETERMINISTIC HC
In this section, the mathematical formulation of PV HC calculation in LVDS is framed using a simple feeder (Fig. 1) and the π-representation (Fig. 2) of a segment l of such feeder.The LV substation with a transformer is taken as the source (or slack) bus, with the active and reactive generation of the source Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
limited by the transformer capacity.Although each transformer supplies two to eight LV feeders, this assumption is justified by the ease of upgrading the transformers compared to underground cables.Consideration of transformer rating S rated t as the limit of the active power of the grid prevents unrealistic scenarios such as high rated PV at the beginning of the feeder.
LV distribution systems in Europe consist of a set of radial LV feeders.The LVDS HC is computed as the sum of the HC of the individual feeders.The PV HC of a feeder is the total PV that can be added to that feeder without violating its operational limits during any realistic operational scenario [2], [9].The usual operational limits are minimum voltage U min , maximum voltage U max , maximum voltage unbalance limit, cable current rating I rated or transformer rating S rated t .The operational scenarios are the result of the uncertainty in the irradiance and load of the consumer, while the planning scenarios are due to the uncertainty in the size of the PV, location of PV and type of PV [3].The feeder's HC is computed as the maximization of the sum of the additional PV capacity a subset of consumers on the feeder can accommodate, i.e., the planning scenario with the highest PV, without violating its operational limits: where, P is the set of new PV installations in the feeder, assuming each LV consumer d has the possibility to have a PV p of rating P kWp p .Fig. 2 shows the simple π-model of the single phase equivalent of a radial LV feeder segment to describe the OPF formulation.Buses i ∈ I are the vertices of the graph representing the feeder and branches l ∈ L are the edges of the graph.The bus voltage in rectangular coordinates is: In the deterministic formulation, the voltage magnitude The parameters y sh i , g sh i and b sh i denote the corresponding admittance, conductance and susceptance at a bus i.The current flowing through a branch l from bus i to j is: where, lij ∈ T l→ links a branch l to its from-i and to-bus j.
The union set T l = T l→ ∪ T l← has the reversed branch set T l← which links a branch l from its to-bus j and from-bus i.The current magnitude |I lij | is bounded by I rated l .z s l denotes the series impedance of a branch l, where r s l and x s l are the corresponding series resistance and reactance.y sh lij denotes the shunt admittance of a branch l at bus i, and g sh lij and b sh lij are the corresponding shunt conductance and susceptance.
Units u ∈ U generalize source, loads, and distributed generators in the feeder.A tuple ui ∈ T u links a unit u to a specific bus i.The current flowing from the bus i to the unit u and complex power are, of the order of the transformer rating.As previously mentioned, these limits are to avoid unrealistic PV scenarios and should not be confused with transformer loading (Fig. 3).To consider reverse power flow, the minimum real and imaginary powers of the source, P min s and Q min s , are considered to be negative of the maximum limits.In a deterministic HC problem, the assumption is made that each consumer can have up to P max p PV capacity.It is assumed that PV operates at a unity power factor such that Q min p = 0 and Q max p = 0. Peak irradiance is indicated by E in kilowatt/m 2 , so the actual PV generation is P p = E • P kWp , neglecting inverter efficiency, soiling loss, etc. PV systems smaller than 5 kWp are usually connected single-phase, while larger systems are connected three-phase.However, when solving as a single phase-equivalent, all PVs are considered balanced, and so is the load.The variable P kWp p is assumed to be a continuous variable and is a relaxed version of the discrete PV size.It is assumed that all consumers can have only one PV with a rating between P min p and P max p .Once the maximum PV size is computed, the individual P kWp p can be rounded to the lower value of PV size in PV ratings set (P rt ) if such a set is available.
Consequently, the feasible set of the HC formulation as IVR formulation of the standard deterministic OPF problem is presented below for more clarity: OPF HC Reference Bus Constraint: Source generator Constraint: Bus Constraints: Branch Constraints: Demand Constraints: PV Constraints: The objective of the optimization is to maximize the total PV capacity as given in (1).The uncertainties in both load and irradiance are neglected in this formulation.The final hosting capacity is the sum of the individual PV P kWp p .The output is the highest possible PV that the feeder can accommodate under network constraints for a particular snapshot of the operational scenarios.For PV HC with equal PV, the PV constraints in

III. STOCHASTIC OPF HC
The formulation presented in Section II assumes deterministic load and irradiance set-points.However, in reality, these set points are random and time dependent.All variables in (10a) become stochastic variables for each time instance t under consideration.In [3], overvoltage was shown to be the most limiting factor for the calculation of PV HC.Furthermore, after a Probabilistic Power Flow (PPF) based evaluation, the maximum probability of overvoltage was deemed as a suitable index to calculate PV HC in an LVDS.Instead of running the OPF for HC calculation for all time period t ∈ T , the search space is reduced by only investigating the time-point τ where the probability that the maximum voltage in the feeder is higher than U max is highest, i.e., p ov t = P (max(U i,t ) > U max ) is the highest.This instance is found a priori by using PPF from [4] for all timestamps, assuming high PV penetration.The remainder of the CC-OPF formulation is described for single-period optimization for τ time where p ov τ = max(p ov t ).Stochastic variables are captured by a real-valued finitevariance stochastic germ ω = [ω 1 , . .., ω N ] T with N ∈ N, and corresponding set of possible realizations Ω ⊂ R N .The first (N-1) stochastic variables for the selected time τ are due to the uncertainties in load.The last variable ω N describes the solar irradiance E.
The gPC-CC-OPF used in this paper for stochastic HC is based on the Galerkin projection of the uncertain variables in the power flow equations.The projection is made in the form of the inner products of orthogonal polynomials for each uncertain variable considered.The gPC converts the infinite-dimensional space of the continuous uncertainty into a finite-dimensional vector K based on the orthogonal polynomials [23], [24].The value of dimension K depends on the number of uncertainty variables considered and polynomial degree considered, i.e., for second degree polynomials with N sources of uncertainty, the dimension of K is given by (2+N )! 2!N ! .The gPC based reformulation of the stochastic problem is mainly based on the properties of orthogonal polynomials summarized below, and interested readers are referred to [23], [24] for details: a) The continuous uncertain variable with known probability density function can be accurately represented by degree one orthogonal polynomials, i.e., with two coefficients.b) The tensor product of univariate orthogonal polynomials can be used to create multivariate orthogonal polynomials [Φ k ] k∈K that act as the base for all variables under consideration.c) The inner product of the orthogonal polynomial bases (also known as expectation operation) can be precomputed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and is given as: where, γ l is a positive scalar that can be precomputed and δ lk is the Kronecker delta.d) The first coefficient of the polynomial representation of the variable is the mean of that variable, and the square of the remaining coefficients is the variance of that variable.The approximate polynomial representation of a random variable x is: where, x k is the coefficient of Φ k , the orthogonal polynomial basis.The approximated mean and variance of x are as follows: The propagation of the uncertainty from the reference load and generation to the nodal voltages and current involves mainly two operations, summation and multiplication.For three approximated random variables x, ŷ, and ẑ which are defined on the same orthogonal basis, the Galerkin projection of the operations are as follows: where M denotes the multiplication tensor, and can be computed ahead of time.
In a stochastic context, enforcing deterministic bounds on continuous random variables is not possible.Alternatively, chance constraints ensure that the probability of a random variable x violating a deterministic bound x min or x max is below a certain stochastic bound ε: The choice of ε depends upon the nature of the solution required.
If ε = 0, the solution obtained is distributionally robust, which could lead to underestimating the actual HC of the network [2].This can be reformulated as a moment-based equation [19], [20], where λ(ε) > 0 is chosen based on knowledge of the random variable.In [19] and [20], λ(ε) was taken assuming the output as Gaussian for all nodes and lines, i.e., λ(ε) = 1.65 if ε = 0.05.In this paper, an empirical tuning method is used where instead of taking the same λ(ε) for all limits, each node is assigned a separate moment-based reformulation factor λ i (ε v ).Similarly, each line has a particular λ lij (ε i ).
The following auxiliary variables are used in this formulation.1) squared bus voltage magnitude: , 3) active source power: P s = U re i I re s + U im i I im s , and, 4) reactive source power: The reason was to eliminate the fourth-order tensor product in gPC-CC-OPF formulation on IVR formulation space.Interested readers are referred to [20] on how the fourth-order permutation of K presented in [25] is reduced to second-order permutations by the use of these auxiliary variables.
Finally, the feasible set of the IVR formulation of the gPC-CC-OPF for chance constrained PV HC is: gPC-CC-OPF HC Reference Bus Constraints -i ∈ I ref : Source generator Constraints -si ∈ T s Bus Constraints -∀i ∈ I : From Branch Constraints -∀lij ∈ T l← : Branch Constraints -∀lij ∈ T l : PV Constraints -∀pi ∈ T p : where M denotes the multiplication tensor for the Galerkin projection of uncertainties, and ] k∈K is employed for convenience to represent the vector of the coefficient variable of X (•) .Finally, an objective is introduced to maximize the total PV size : IV. NUMERICAL ILLUSTRATION In this section, the application of OPF and CC-OPF for HC calculation in a test feeder is shown.Then, the method is used to compute the PV HC of a set of actual LVDS feeders in [26] in Section V.The aim is to establish the method and discuss the limitations and applications of this tool.The code is available on github 1 as an independent project using StochasticPow-erModels.jl, along with the feederwise data of the network in [26].
All the simulations in this paper were done on a PC with a 2.11 GHz Intel i7-8650 U CPU processor with 16 GB memory, using JULIA 1.6.5.No parallelization of any sort was used in all of the calculations.The non-linear Ipopt solver v0.7.0 is used for the calculations.

A. Test Feeder
A simplified test feeder, as shown in Fig. 3, is considered for a detailed study.The feeder has seven consumers, with each consumer having the possibility to have a PV installation between 0-15 kWp.The single time-period optimization is done on the timestamp with the maximum probability of overvoltage on a typical high irradiance spring day.The time-stamp chosen was 15:00-15:15 hrs as it was the time period with the highest p ov t in a priori PPF with all consumers having 10 kWp.The PV generation was considered perfectly correlated as the span of the feeder is 300 meters only.The consumers were classified into four categories based on their annual consumption and connected kVA.The consumers in the same category were considered to be linearly correlated.All the load and PV distributions obtained from the historical profiles were fitted as a beta distribution [3], [27].In Table I, the details of the uncertainty terms considered for the test case are given.In this work, all consumers in the same group are assumed to be perfectly correlated.However, if available, the correlation factor between the consumers can be readily endorsed in the model.The deterministic calculations were done assuming the maximum irradiance and minimal load of 0.1 kW for all consumers.

B. Deterministic OPF HC
This section compares the result obtained from OPF HC with the deterministic scenario based HC method.In the scenario based method, several planning scenarios were taken, and the best planning scenario within the operational limits was considered to be the PV HC of that feeder.The main drawback of this method is the number of scenarios to be considered.For this test case, for seven consumers with the possibility of having PV from 0 to 15 kWp in steps of 1 kWp the number of scenarios is 268 435 456 (16 7 ).These scenarios are only due to PV size, as other uncertainties, e.g., PV type and connected phase, are not even considered.A brute force deterministic evaluation of all these scenarios might take a few months in computation time, if not years.If the scenario based method is stochastic as proposed in [3], the time for computation will increase by multiple orders.For this practical reason, the OPF HC and deterministic scenario based HC methods are compared, assuming that a similar relationship will hold between stochastic OPF HC and stochastic scenario based HC methods.
A histogram is plotted for feasible and infeasible PV scenarios using only two million randomly selected possible scenarios (only 3% of the possible scenarios) using deterministic limits along with the OPF HC (Fig. 4).The time taken by the OPF HC was 0.241 sec, while it was ≈ 2 days for deterministic scenario based HC (Table II).The OPF HC calculated was 50.4 kWp, while the scenario-based HC was 50 kWp, the closest possible integer value below OPF-HC.There are 41 621 scenarios where the sum of PV installed is equal to 50 kWp, of which 3252 were feasible combinations (≈8% of total scenarios with 50 kWp).In Fig. 4, all the feasible solutions are less than (on the left-hand side) the OPF HC.While the unfeasible solutions are on both Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.sides of the OPF HC.The OPF HC is hence the theoretically possible upper limit as there are unfeasible scenarios where the total PV is less than this limit when the location of PV installations changes.An iterative algorithm was applied for scenario reduction where the scenarios with a sum of PV below a feasible scenario were considered to be feasible automatically.For e.g., if the scenario with 30 kWp PV is feasible, all the PV scenarios below that value were not evaluated.Still 3 000 scenarios were required for the test case to reach the accuracy of the OPF HC method.The computation time of the method using scenario reduction was still 20 mins for deterministic scenario based HC calculation, compared to 0.241 sec for OPF HC calculation.Another possible scenario reduction method is to assume all consumers to have equal PV as proposed in [18].In such case, the maximum number of PV deployment scenarios required is 25, if the maximum PV is considered to be 15 kWp in a step of 0.1 kWp.Using this method, the deterministic HC calculated is 10% less than obtained from the OPF HC method, while the computation time is ten times the time needed for the OPF HC method.
The scenario reduction schemes for operational random variables (load and irradiance samples) was dealt with in detail in [2], where 1 000 samples from Sobol or Halton sequences were found to be enough for the probability of congestion calculation.Hence, the computation time for MC based stochastic HC was calculated considering those quasi Monte Carlo schemes and hence no other scenario reductions were considered.

C. gPC-CC-OPF PV HC Assuming Gaussian Output
In Fig. 5, the OPF HC method is compared with the gPC-CC-OPF HC method.The gPC-CC-OPF considers that the voltage or current violation probability is not strict, but can be risk based, i.e., can be violated for some scenarios, which is denoted by the stochastic limit .The gPC-CC-OPF HC of the test network is calculated to be 52.4 kWp, 4% more than the result obtained with the deterministic limit.The computation time for gPC-CC-OPF HC method was 4.04 sec, almost 16 times more than the OPF HC method.In Table II, the computation time for all methods is given.
The PV size of each consumer is a variable in the above calculation.The advantage is that each PV installed can have different constraints for PV size; for example, all consumers can have equal PV.This means that the proposed CC-OPF HC calculations can be used when there is limited knowledge of the scenario or to evaluate different fairness schemes in the feeder.The individual limits of the consumer PV size is an input for the method.In Fig. 6, four different scenarios of PV size limits are depicted, and the respective gPC-CC-OPF HC are as follows:

D. Tuning of Moment-Based Chance-Constraints
One of the limitations of gPC-CC-OPF is that the momentbased reformulation in (20), wherein we assumed the λ i (ε v ) to be based on the Gaussian distribution where its value is 1.65 for v = 0.05 [19], [20].But the input we assume is not Gaussian, and the non-Gaussian uncertainties propagated through nonlinear power flow equations resulting in an output that is certainly not Gaussian.In Fig. 7, two distributions are taken, i.e., Normal and Beta, with the same expectation and variance.The figure shows that the Beta distribution has a more skewed tail compared to the Normal distribution.The cyan region corresponds to the area for ε v = 0.05 for the Normal PDF, which is exact.For the Beta distribution, the area covered by the PDF after the point given by Gaussian reformulation is higher than 0.05 (the space between the two PDFs).For the moment-based reformulation of a such probability distribution, the value of λ i (ε v ) needs to be increased to get ε v = 0.05.An iterative method is proposed to correct the value of λ i (ε v ).In the first run of gPC-CC-OPF, the dual variable is recorded for the W i constraint in (25f).The λ i (ε v ) is tuned only for the bus with a binding constraint, i.e., the constraints with a non-zero dual variable.For the test case, this is bus 5.The voltage is sampled at that bus and based on the 95% quantile of U i , a new λ i (ε v ) is decided.The second run of gPC-CC-OPF is done after the tuning of moment-based formulations.The empirical tuning method is proposed, where a new λ i (ε v ) is assigned to the node with binding chance constraints to reduce the discrepancy created due to assuming the output as Gaussian.
a) If 95% quantile of U i > 1.001U max : After tuning, the probability P (U i < U max ) ≥ 0.95 is satisfied.The new computation time is 9.26 sec, which is 2.3 times the computation time required for the original gPC-CC-OPF problem.The gPC-CC-OPF is solved twice, and the remaining time is due to the sampling of the voltage, finding quantiles and tuning λ.The initial and after-tuning probability of overvoltage   for each node is shown in Fig. 8.The new HC calculated is 52.7 kWp which is 0.38% higher than before tuning.As the name suggests, the tuning is done to improve the accuracy of the gPC-CC-OPF HC method and is an additional feature to the original solution.

V. CASE STUDY: EUROPEAN TEST NETWORK
In this section, the OPF HC method is compared with gPC-CC-OPF HC calculation with tuning for the large-scale network presented in [26] (Fig. 9).
The OPF HC is calculated assuming maximum irradiance in the day and a load of 0.1 kW for all consumers.The analysis is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.done feederwise instead of the whole network at once.For gPC-CC-OPF HC computation, a clustering was performed to identify similar loads according to the annual energy consumption, and the loads in the same clusters are assumed to be perfectly correlated.The load of each type of consumer is fitted to a Beta distribution per 15-minute time period for a typical spring day.Similarly, the irradiance of the typical spring day is fitted to a Beta distribution using historical measurements in EnergyVille premises of Genk.It should be noted that a Beta distribution is taken in this study as it is mostly used to represent both irradiance and consumer load at household level [27].However, the method is valid if the input uncertainties are represented by any other standard distribution, such as Normal and Uniform distributions, which are usually referred to under Askey schemes [23], [28].For details on supported distributions, readers are pointed to the book Hypergeometric Orthogonal Polynomials and their q-Analagous by Koekoek, Lesky & Swarttouw [29], which has the current extension of the Askey scheme.The reactive power is considered to be 5% of the active power for all loads.
Out of 147 feeders with consumers in the network, gPC-CC-OPF does not converge to an optimum value for 8 feeders within 500 secs or 3000 iterations.The results shown are based on 139 feeders where the OPF HC and gPC-CC-OPF HC calculation both converged to an optimum value.For gPC-CC-OPF HC, the time-stamp of maximum congestion probability is calculated beforehand using probabilistic power flow when all consumers have 10 kWp PV installed.Single time period gPC-CC-OPF solutions with tuning are obtained for each feeder with uncertainties in PV generation and consumer load.

A. OPF HC vs gPC-CC-OPF HC in the Actual Network
In Fig. 10 The histogram shows that the deterministic OPF HC is below the HC of CC-OPF HC, concretely on 126 feeders.The lower value from OPF HC is attributed to the assumption of a constant load of all consumers and peak irradiance, which is only one snapshot of all possible operational scenarios due to load and irradiance.It shows that the assumption of low load and high irradiance gives us the solution in a more distributionally robust space (worst-case).
The computation time for gPC-CC-OPF HC was compared as a function of the number of nodes and uncertainty sources.The computation time of gPC-CC-OPF HC for 139 LVDS feeders with respect to the number of sources of uncertainty is shown in Fig. 11.The sources of uncertainty in the LV feeders range from two to eleven.The computation time increases with the number of uncertainties.Moreover, the number of feeders where a solution is not obtained within 500 sec and 3000 iterations is high (three out of five feeders) for the case with ten uncertainties.In Fig. 12, it is seen that the number of nodes in the LVDS feeder also influences the computation time.However, the feeders with fewer nodes can also have higher computation time if the number of sources of uncertainty is higher.In this study, most feeders had less than ten uncertainty sources.Furthermore, the stochastic HC of 90% of the LV feeders could be computed within 100 secs.While gPC-CC-OPF method has some room for algorithmic improvement, in particular for systems with a higher number of uncertainties, the results show that the method has a strong potential.

VI. CONCLUSION
This paper introduces a SOPF method for stochastic HC calculation using gPC-CC-OPF.This method uses stochastic load and irradiance input for a peak summer day to get the best PV placement scenario for low voltage distribution feeders.Stochastic chance constraints in the operational variables, i.e., nodal voltage and branch current, are considered.The use of gPC-CC-OPF enables us to solve the problem as a single-shot problem without any linearization, sampling, or relaxation.This method gives the stochastic HC in a few seconds while solving a similar problem using MC-based methods would take a few Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
days.The comparison between HC obtained from deterministic OPF and stochastic gPC-CC-OPF in actual LVDS feeders shows that the OPF HC calculated are an underestimate of the CC-OPF HC.Finally, an iterative empirical tuning approach for the moment-based reformulation of the chance constraint is proposed to obtain a correct chance constraint limit when the inputs are not Gaussian distributed and the distribution of nodal voltage and current are unknown.The method is made openly available as a branch of STOCHASTICPOWERMODELS.JL along with the network data and the uncertainties definitions.
For a planning application, such as the calculation of the PV HC, a single-time period optimization for a balanced system is sufficient.However, for computing the PV-battery systems and electric vehicles, a multi-time period SOPF is recommended.Similarly, if the interest is in obtaining the HC of single-phase PV systems, an unbalanced SOPF is recommended.Hence, an unbalanced, multi-period SOPF is identified as future work.Furthermore, improvement of the computation time of the SOPF using gPC by leveraging the sparsity of the coefficients will be the future work.
Arpan Koirala (Graduate Student Member, IEEE) received the master's degree under the Erasmus Mundus Master Course in sustainable transportation and electrical power systems from the consortium of the University of Oviedo, Gijón, Spain, La Sapienza, Rome, Italy, University of Nottingham, Nottingham, U.K., and Polytechnic Institute of Coimbra, Coimbra, Portugal, in 2018, and the Ph.D. degree from the KU Leuven, Leuven, Belgium, in 2023.His research interests include stochastic programming, optimization, distribution system modelling, energy market, and high-scale integration of renewable resources.
Tom Van Acker born on 1990 in Roeselare, Belgium.He received the M.Eng., M.Sc.and Ph.D. degrees in electrical engineering from KU Leuven, Leuven, Belgium, in 2012, 2014, and 2020, respectively.In 2020, he was a Postdoctoral Researcher with KU Leuven.Since 2021, he has been a Power System Expert with BASF Antwerp.His main research interest include stochastic processes, optimization and uncertainty quantification applied in a power system context, specifically in the topics: availability, harmonics, and state estimation.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 2 .
Fig.2.π-model of a branch lij, with series impedance z s l and from-and to-side shunt admittances y sh lij and y sh lji .

Fig. 3 .
Fig. 3. Single phase representation of an LV feeder with all consumers with PV for OPF HC calculation.

Fig. 5 .
Fig.5.HC scenarios obtained from OPF and Scenario based method (PV HC for each scenario is given in the legend).

Fig. 6 .
Fig.6.PV HC scenarios obtained from gPC-CC-OPF with different PV size limits with HC in the legend.

r
S1:The PV installations of consumers are in between 0 to 15 kWp.The corresponding CC-OPF HC is 52.5 kWp.

r
S2: The PV installations of consumers are equal, and between 0 to 15 kWp.The corresponding CC-OPF HC of the test feeder is 44.8 kWp.

r
S3: Consumer 1 and 2 in the test feeder can have PV installations between 0 and 10 kWp, consumer 3 has a fixed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 7 .
Fig. 7. Moment based reformulation of stochastic limits.The cyan region is the area where ε v = 0.05 for Normal PDF.

r:
S4All consumers have PV installations between 5 and 10 kWp.The corresponding CC-OPF HC is 47.4 kWp.

Fig. 10 .
Fig. 10.Histogram of mismatch between OPF HC and gPC-CC-OPF HC normalized by gPC-CC-OPF HC in percentage.The number of feeders with no mismatch or positive mismatch are plotted as dot.

Fig. 11 .
Fig.11.Boxplot of log 10 of the computation time for gPC-CC-OPF HC in real LV feeders with respect to the number of uncertainties considered.The numbers in the red brackets represent the number of LV feeders whose HC solution did not converge in 500 seconds or 3000 iterations.

Fig. 12 .
Fig. 12. Log 10 of computational time for gPC-CC-OPF HC in real LV feeders with respect to the number of nodes in the LV.The red dots represent the HC solutions not converged in 500 seconds or 3000 iterations.

TABLE II COMPUTATION
TIME OF SCENARIO BASED MC, OPF HC, AND GPC-CC-OPF HC FOR TEST FEEDER (7 CONSUMERS, EACH CONSUMER CAN GET 0-15 KWP IN THE STEP OF 1 KWH)