Quantifying Hosting Capacity for Rooftop PV System in LV Distribution Grids

1 Quantifying Hosting Capacity for Rooftop PV System in L V Distribution Grids Jingyi Y uan, Student Member , IEEE, Y ang W eng, Member , IEEE, and Chin-W oo T an, Member , IEEE Abstract —Power systems face increasing challenges on reliable operations due to the widespread distributed gener- ators (DGs), e.g., rooftop PV system in the low voltage (L V) distribution grids. Characterizing the hosting capacity (HC) is vital for assessing the total amount of distributed gener- ations that a grid can hold before upgrading. For analyzing HC, some methods conduct extensive simulations, lacking theoretical guarantees and can time-consuming. Therefore, there are also methods employing optimization over all necessary operation constraints. But, the complexity and in- herent non-convexity lead to non-optimal solutions. T o solve these problems, this paper provides a constructive model for HC determination. Based on geometrically obtained globally optimal HC, we construct HC solutions sequentially according to realistic constraints, so that we can obtain optimal solution even with non-convex model. For practical adaption, we also consider three-phase unbalance condition and parallel computation to speed-up. W e validated our method by extensive numerical results on IEEE standard systems, such as 8 -bus and 123 -bus radial distribution grids, where the global optimums are obtained and performance illustrations are demonstrated. I . I N T R O D U C T I O N I N recent years, an increasing number of renewable energy-based distributed generators (DGs) appears in power distribution grids. Such widespread use of DGs brings great benefits, e.g., voltage profile support, loss reduction, and lower capital cost. But, such change also raises operational challenges, like power quality problems and system control issues [1]. For example, high-level DG penetration can cause undesirable voltage flicker or equipment aging due to the excessive voltage regulating behaviors [2], [3]. T o avoid these problems, evaluating the hosting capacity of a system is vital to understand the impact of DG penetration level and conduct upgrades for reliable system operation. Hosting capacity (HC) is defined as the maximum amount of the power generation that a system can host without violating any operating standards [4]. W ith such a definition, Electric Power Research Institute (EPRI) proposes simulation-based methods, modeling each feeder and examine all the power quality and reliability issues with screening tools to determine HC at differ ent locations [5]–[7]. For example, we can keep on incr easing the PV (photovoltaics) penetration level for power flow analysis until constraint violation, where the last feasible PV generation is regarded as the HC. Such idea has already been deployed in commercial software. For example, our collaboration partner , CYME, focuses on the hosting capacity of one node per analysis [8]. Fig. 1 shows the three steps how CYME obtain the per-bus hosting capacity . However , no matter if one is focusing on one bus or for all buses, it is numerically impossible to exhaust all possible scenarios with simulations. Fig. 1: CYME for numerical hosting capacity analysis. Therefor e, there are also mathematical model-based methods. For example, [9] formulates HC problem as an AC optimal power flow (ACOPF) problem, extended by [10], [11] for a multi-period ACOPF . As solving ACOPF with many constraints are hard, [12], [13] propose to conduct sensitivity analysis as an approximation to char- acterize key factors, e.g., voltages, the network topology , load size, voltage limits, and generation locations [9], [14], [15]. There also recent study with a pr obabilistic model on key factors [16]. Based on these factors, additional control strategies can improve HC via 1 ) on-load tap-changers (OL TC) technology [17], 2 ) reactive power control [18]–[20], 3 ) active network management [10], 4 ) static and dynamic network re-configurations [11], and 5 ) adding soft open points from power-electr onic devices [21], and 6 ) incor- porating real-time information [9]. Besides, one can also model economic aspects to do a cost-benefit analysis for maximizing HC [22]. However , no matter how to approximate, e.g., select- ing key features, the aforementioned studies can not lo- cate the globally optimal solution. The objective in these methods is also inflexible to reduce optimizing variable numbers. Therefore, this paper focuses on obtaining exact solutions constructively with a flexible objective for unevenly deployed DGs. The first contribution is changing the objective in ACOPF to a flexible form and sequentially adjusting solution with constraints to preserve convexity as much as possible. For example, we transform into rectangular coordinates to convexify certain constraints. When con- 2 vexification is impossible, our second contribution lies in connecting geometric understanding for r eaching a global optimum even when the constraint is still non- convex. A detailed proof is provided. Our third con- tribution is for providing practical adjustments, where we extend our static model with sequence components for three-phase unbalance situation. W e also provide a method for paralleling computation for large systems. Extensive simulations verify the performance of the proposed method on IEEE distribution networks, e.g., 8 -bus and 123 -bus. The results show that our problem formulation and proposed theory can find exact HC in diversified scenarios. It also validates the integrated methods for multi-phase unbalanced situation and the distributed-optimization for computational reduction. The rest of the paper is organized as follows: Section II shows the problem formulation. Section III proposes a generalized theory of maximum HC solution based on geometric understanding and the proof. Section IV shows the algorithm adjustments for three-phase unbal- anced system and parallel computation. Section V tests the theory in various systems and Section VI concludes the paper . I I . M AT H E M AT I C A L R E M O D E L I N G The power flow equations are P inj i = n X k =1 | V i || V k | ( G ik cos θ ik + B ik sin θ ik ) , (1a) Q inj i = n X k =1 | V i || V k | ( G ik cos θ ik − B ik sin θ ik ) , (1b) where n is the bus number of the electrical system, P inj i is the active power injection at bus i , V i is the voltage magnitude at bus i , and G ik + j B ik is the element in the bus admittance matrix Y Bus ∈ C n × n . Specifically , G ik repr esents the line admittance between bus i and bus k . For flexible objective, we include a free parameter λ in front of each potential generator . H C = n X i =1 λ i · P inj i = n X i =1 λ i · ( P g en i − P load L ) , Here, λ can be a generating scaling factor , indicating how fast a linearly increas with “base” generation can be. For example, λ 1 = λ 2 = · · · = λ n are widely assumed in literature for simplification, wher e we ar e looking into an ideal and even growth for situational awareness. How- ever , practical setup won’t be homogeneous in general, therefor e, our model provide such a flexibility . • when λ i is a binary parameter , λ i = 1 indicates bus i has solor generation and λ i = 0 vise verse; • when λ i is a non-negative parameter , e.g., λ i ∈ R + , the weight parameter repr esents the importance of a bus, e.g., the relative price of installing solar panels. In this paper , we will first let λ i = 1 , ∀ i for simplification and focus on the variable P i for geometric insights, e.g., n X i =1 λ i · P inj i = n X i =1 P inj i . Then, we will remove such an assumption for generality . Operational constraints must be satisfied while calcu- lating hosting capacity . max n X i =1 λ i · P inj i (2a) s . t . Equ. (1) , (2b) V min ≤ | V i | ≤ V max , (2c) θ min ≤ θ ik ≤ θ max , (2d) I Lmin ≤ I ik ≤ I Lmax , (2e) pf i ≥ η , (2f) where (2b) is the power-flow constraint. (2c) and (2d) state the constraints over voltage magnitudes and the voltage angle differ ence. (2e) represents the thermal limit. (2f) is the power factor requirement of generator bus, where pf i = ( P inj i ) 2 ( P inj i ) 2 +( Q inj i ) 2 . I I I . F I N D I N G T H E G L O B A L O P T I M U M O F H O S T I N G C A PA C I TY B A S E D O N G E O M E T R I C A L U N D E R S TA N D I N G A. Geometrical Intuition for Obtaining Hosting Capacity It is challenging to obtain global optimum in 2, due to the nonlinear correlations of the optimization vari- ables in power flow equation, let alone the non-convex constraints. However , for a small system such as a 3 -bus system, hosting capacity can give geometrical results and intuition, shown in Fig. 2 . For example, we can let bus 0 be the refer ence bus with voltage 1 ∠ 0 ◦ . Let bus 1 , 2 be PV -type buses and let the branch resistance be R = 1 . For simplicity , let the voltage angle differences and line reactances be zeros. Based on the setup, we show below two geometric intuitions to find the global optimum. Fig. 2: A 3 D plot of P 1 + P 2 , V 1 , and V 2 that illustrates the P-V (total power-individual voltage) relevancy for the 3 - bus network. The value range of both V 1 and V 2 respect to the total power makes it a surface. The solid red point denotes the operation point with maximum generation where V 1 is at upper bound and V 2 is at lower bound. 3 For power system analysis, nodal P-V (power-voltage) curve is commonly used which is a 2-D plot. T o find the HC due to variable correlations, we extend to a 3- D P-V plot in Fig. 2a. In the 3 -bus visualization, the only voltage constraint ( 0 . 95 ≤ | V | ≤ 1 . 05 ) leads to the feasible solution region. If V 2 is fixed, we observe that the curve monotonically increases as V 1 increases, which leads to finding the possible solution at the top edge. Then we find the maximum generation capacity at one vertex where V 2 is at lower bound. As marked by a red point, the solution of HC is obtained when V 1 = 1 . 05 p.u. is the highest and V 2 = 0 . 95 p.u. is the lowest. Such a solution pattern indicates hosting capacity be in direct proportion to the differ ence of the voltage magnitudes over the neighbor buses. Fig. 3: The plot of P 1 v s. P 2 . The voltage constraints bound a feasible region under power flow equations. W e use a straight line: H C = max ( P 1 + P 2 ) to cut the maximum value. The same solution point is found considering pairwise power correlation. It denotes that bus 1 generates 0 . 1575 p.u. power while bus 2 consumes 0 . 095 p.u. power . In addition to voltage correlation, it is important to validate voltage solution more directly in the pairwise power domain to understand how voltage constraints impact the optimum solution. For this purpose, Fig. 2b is created by firstly plotting all possible ( P 1 , P 2 ) pairs as the gray dots according to power flow equations. Secondly , we add the voltage upper bounds and lower bounds of bus 2 and bus 3 to eliminate points outside the allowed region. Then, we push the line of P 1 + P 2 towards bigger values until we reach the global optimum. The red star in Fig. 2b confirms that the same optimal HC is reached when the neighbor buses’ voltages are with the pattern of one high and one low . The positive value of P 1 indicates generation while the negative value of P 2 means it absorbs power from bus 1 and hence the total is the maximum. Such hosting capacity take neighbour- bus correlation into consideration that allows all possible patterns for different buses. W e assign relative numbers to represent the bus order , e.g., bus 1 is odd and bus 2 is even. T o observe mathe- matically , the opposite but extreme values for neighbor buses attain individual maximization for the quadratic functional component in power flow equations. The con- clusion is feasible in both resistive and general inductive network. Based on the geometric summary , we can infer the general theory to the solution and the condition to achieve it next. B. Mathematical Theorem for Global Optimum Solution 1) HC due to V oltage Magnitude Constraints: Based on the geometric understanding on voltage magnitudes bounds, we can generalize the results from 3 -bus to n - bus radial network high-low-voltage pattern. Theorem 1. Let a radial distribution network be with resis- tive or inductive line impedance and let the voltage magnitude constraint be the only constraint on HC. The maximum power generation in ( 3 ), H C ∗ = H C max , is obtained when | V odd | = V max and | V ev en | = V min . For the proof, it is sufficient to show that the proposed hosting value H C ∗ is larger than arbitrary summation of power generations, e.g., H C ∗ = n X i =1 P i, solution ≥ n X i =1 P i , (3) Proof. See Appendix B. 2) HC due to V oltage Angle Constraints: W e add con- straints on the angle difference, leading to complex voltage phase. Theorem 2. Let θ min = 0 ≤ θ i ≤ θ max . The HC is given by θ ik = min[ π , θ max ] , | V | =      V max for all buses , if θ max > arccos( V max + V min 2 V max ); ( V max for bus 2 n − 1 , V min for bus 2 n, if 0 < θ max ≤ arccos( V max + V min 2 V max ) , where arccos( V max + V min 2 V max ) is the critical value of the voltage angle boundary that changes the solution pattern. Proof. When the angle constraint is added, we can still start with P N i =1 P i and obtain N X i =1 P i = N X i =1 N X k =1 ( − G ik ) · [( V i, re − V k, re ) 2 (4) + ( V i, im − V k, im ) 2 ] . The most valuable part is this quadratic form, although, sinusoidal components about angle is introduced. De- tailed proof of the derivation of (4) can be found in Appendix B. Subsequently , (4) can be written in the polar coordinates to quantify the impacts of voltage magnitude and phase angle separately . N X i =1 P i = N X i =1 N X k =1 ( − G ik ) · ( V i cos θ i − V k cos θ k ) 2 +( − G ik ) · ( V i sin θ i − V k sin θ k ) 2 . (5) As the multiplication is the key component in the quadratic optimization, we can simplify the objective 4 further by having | V i | = a, | V k | = b, cos θ i = a 1 , sin θ i = a 2 , cos θ k = b 1 , sin θ i = b 2 . Suppose bus i is the even bus and bus k is the odd bus. Such notation reshape the problem model into max n X i =1 ( − G ik ) · h ( aa 1 − bb 1 ) 2 + ( aa 2 − bb 2 ) 2 i (6) s . t . V max ≤ a ≤ V min , V min ≤ b ≤ V max , a 2 1 + a 2 2 = 1 , b 2 1 + b 2 2 = 1 . (7) Substitute the equality constraint (7) into the objective function (6) and we get renewed problem formulation, max n X i =1 ( − G ik ) · [ a 2 + b 2 − 2 ab · cos θ ik ] s . t . V max ≤ a ≤ V min , V min ≤ b ≤ V max . In the distribution grids, G ik is mostly negative, making ( − G ik ) positive. T o find the maximum objective value in a non-convex problem, we employ a piecewise analysis of 2 ab · cos θ ik according to θ . W e only show the reason- able upper and lower bounds of the angle differ ence in the following and the rest can be found in Appendix B. • θ ∈ [ − θ max , θ max ] , arccos( V max + V min 2 V max ) ≤ θ max ≤ π : W ith voltage angle constraint, cos θ max ≤ cos θ ik ≤ 1 . As a and b are positive, cos θ ik should be as small as possible to get a bigger value of the objective, which is cos θ max . T o maximize a 2 + b 2 − 2 ab · cos θ max , the voltage pattern should be at upper bounds ( a = b = V max ) when cos θ max is a negative or small positive number . Until θ max keeps decreasing to a specific number , the solution pattern changes to “high-low voltage” as in Theorem 1. For calculating this critical value, V 2 max + V 2 max − 2 · V max V max cos θ max ≤ V 2 max + V 2 min − 2 V max V min cos θ max , θ ≤ arccos( V max + V min 2 V max ) . If we plug in the boundary values V max = 1 . 05 and V min = 0 . 95 below , θ ≤ 0 . 3098( r ad ) . Therefor e, the solution for this angle range is a = b = V max , θ ik = θ max . • θ ∈ [ − θ max , θ max ] , θ max < arccos( V max + V min 2 V max ) : As we mention above, the solution pattern flips. a = V max , b = V min , θ ik = θ max , when voltage angle is limited below the critical value arccos( V max + V min 2 V max ) . When the voltage angle constraint narrows down, the variation of the results shows that the system acquires power generation prior to the difference of neighbor bus voltage angles. However , in the real systems, the angle differ ence between the two ends of a 100 km line is approximately 6 degrees. Therefor e, θ max can be set to be a small value below 0 . 1047( r ad ) in the distribution grids, under which the second piece is meaningful. 3) HC due to Other Constraints: In addition to the voltage constraints (for | V | and θ ), we obtain solution restricted to operational constraints. Solution Adjustment According to Thermal Limits The thermal limit plays a significant role in two-way power flow . Therefore, this constraint is represented by current flow limit that results from the heating effects of devices, I Lmin ≤ I ik ≤ I Lmax = − I Lmin = C . W e find the thermal limit boundary as inspired by the geometrical understanding and update the HC subject to both constraints. Theorem 3. Consider the following thermal limit in addition to the constraints on top of Theorem 3 : I Lmin ≤ I ik ≤ I Lmax , I Lmax = − I Lmin . Let I ∗ ik be the branch current flow of previous solution pattern without thermal limit boundary . • If I Lmin ≤ I ∗ ik ≤ I Lmax is satisfied, the HC is achieved by the same voltage solution in Theorem 3 ; • Else, the HC is achieved by the solution pattern calculated from a three variable cubic equation in Appendix B. W e will check if the solution point in Theorem 3 encounters the violation boundary or not. If not, the so- lution for HC stays the same. If yes, we will compute the new voltage solutions which is provided in Appendix B. Solution Adjustment According to Power Factor Limits PV (photovoltaics) generator produces DC power , meaning it injects power at unity power factor . Nonethe- less, the local load and power supplied by the grid include r eactive power , and PV inverter is also utilized to regulate the power factor by producing reactive power . Reactive power of all PV generators must lie within specific bounds, which in other words is to limit the power factor (Appendix B). Similar to the thermal limit correction, we substitute previous solution point to check if the inequality is true. If the reactive power is beyond the limit, we can switch bus type to fix reactive power at boundary value and compute the new voltage values and HC by Q-V sensitivity [23]. Such a constraint is non-convex, but we can still find the theoretically optimal solution. Other constraints can be analyzed following a similar way . I V . P R A C T I C A L S Y S T E M C O N D I T I O N S In real system operation, our proposed model needs to be accommodated to different scenarios. W e present two typical methodologies that integrate with other methods facing practical and operational difficulties. A. The Multi-phase system with Unbalanced Loads The single-phase power flow model has a strong assumption that ignores the unbalance due to some unbalanced loads and untransposed lines. T o provide a tool that is more analytical and practical, we transfer multi-phase model to sequence components and de- couple them into separate subproblems. After that, our proposed model can be applied and transfer the results back to the origin [24]. Let a , b , and c denote three phases, the scalar vari- ables of the single phase become vectors. | V abc i | = 5 [ | V a i | , | V b i | , | V c i | ] T and θ abc i = [ θ a i , θ b i , θ c i ] T denote voltage magnitude and angle at bus i . P abc i = [ P a i , P b i , P c i ] T denotes the active power and so is Q abc i = [ Q a i , Q b i , Q c i ] T . The self and mutual line admittance are considered, so a single element Y ik is replaced by a 3 × 3 matrix and then Y abc B us ∈ C 3 n × 3 n ( Y ik ∈ { Y s ik , Y z ik } if considering both series and shunt admittance). For single-phase or two-phase cases, simply set the quantities of the other phases to be zero and the model will degrade (e.g., | V abc i | = [ | V a i | , 0 , 0] for single phase.) W e simplify the analysis of multi-phase system by transferring the model into the the sequence form. In order to apply the unbalanced situation, we decouple the sequence model while using current injection from negative and zero sequences to repr esent the unbalance. According to [24], the positive-sequence voltage mag- nitude is much larger under unbalanced condition and the power flow equations are similar to single-phase formulation. The equations allow us to implement the proposed method and obtain voltage pattern solutions for positive sequence. Based on the decoupled models, we can compute the negative and zer o sequence voltages via nodal voltage equations. Combining with results of positive sequence power flow , it’s able to transfer back to the three-phase solution via (29) (we provide details in Appendix C). B. The Distributed-Optimization for Complexity Reduction For a large and complex system, the time on heuristi- cally solving the HC is much longer as too many vari- ables and constraints causing computational difficulties. In the following, we show that we can segment the tree structur e and compute distributedly . If a system is divided into m parts, the red seg- mentation point shown in Fig. 3 links 2 subsystems, it correlates to both parts geometrically . Mathematically in the problem model, the variables are divided into two sets. The first set includes local/private variables l t for subsystem t and the second set includes coupling variables c t which bridge two more subsystems. This turns the problem into the distributed formulation. max m X t =1 H C t ( l t , c t , c 0 t ) s . t . c 0 t = c t − 1 , V min ≤ | V | ≤ V max , θ min ≤ θ ≤ θ max , I Lmin ≤ I L ≤ I Lmax . (8) Fig. 4: The system is segmented into m parts for easier computation. s t is one segment bus linking 2 subsystems. In the model, H C t means the HC of subsystem t ; l t repr esents the local variables - bus voltages that only appear in one subsystem. c t and c 0 t repr esent the com- plicating variables - bus voltages that appear in more than one subsystems. The objective function looses the coupling as it changes from the summation of power injections for all buses to the summation of HC for all subsystems. The coupling shifts to the constraints where the segmentation is true only if c 0 t = c t − 1 . For each subsystem t , the objective is the same to analyze as (4) (the proof is shown in the Appendix D.) Thus, it is capable of segmenting the model without disrupting the solution pattern in theorems and parallel computing. V . N U M E R I C A L R E S U LT S For numerical verification, we use various distribu- tion systems, such as modified IEEE 8 -bus and 123 -bus systems. W e aim at 1 ) validating the hosting capacity solution theory , 2 ) evaluating the integrated method on multi-phase unbalanced system and 3 ) showing the effectiveness of parallel computation. W e use Matlab and its MatPower package for the demonstration [25], [26]. The optimization model is formed as non-linear programming, for which com- monly used solvers include Fmincon, Knitro, and Ipopt. By comparing the r esults fr om one of the classical solvers with our theoretical solution, we find that their perfor- mances are unstable. In other words, they are easily stocked at local optimal especially when the system size is large and the initialization is inappropriate. A. Hosting Capacity W e use the IEEE 8 -bus system for illustrating HC according to various constraint conditions. Specifically , each bus is modeled with DG as a PV -type bus and we apply the proposed method to calculate HC. (a) The comparison of HC subject to different constraints. (b) The voltage profiles with different λ i sets. Fig. 5: Numerical results for modified IEEE 8 -bus system. Fig. 5a presents hosting capacity subject to differ ent constraints and highlight the limiting bus node. The results of extensive simulation-based method are always smaller than ours. Considering differ ent constraints, the HC from our theory under voltage limits is 59 . 33 p.u. , which is relatively large. Resulted from the thermal limit 6 of the branch between node 3 and node 4, HC decreases to 22 . 63 p.u. Sometimes, only a few buses are interested or equipped with PV . For example, a binary λ i denotes if a bus is to be planned with PV or not. Here we set bus 2 , 5 , and 7 in the IEEE 8 -bus network, to have PV generators, e.g., λ i = 1 , i ∈ { 2 , 5 , 7 } and λ i = 0 , i ∈ { 1 , 3 , 4 , 6 , 8 } . The objective turns into P i =2 , 5 , 7 P i . As a result, Fig. 5b shows the corresponding voltage profile that compares to the original setup. As we can see that bus 2 , 5 , and 7 are automatically assigned high voltages to generate as much as possible due to their generator bus type. For comparison, we also plot the voltage values by another method in MatPower . It simply runs power flow heuristically and stops at a local optimum or flatten r egion of the objective. The calculated HC value is 15 . 38 p.u. , much less than 29 . 45 p.u. , show- ing the importance of the closed-form solution proposed in this paper . T o show the difference with respect to the first validation when all the λ 0 s are the same, we draw the voltage setup in gray in the figure as well. This shows how changing λ i to impact the way to achieve the hosting capacity . B. Multi-phase Unbalanced Condition Fig. 6: The Multi-phase setup in IEEE 8-bus network. Scenario Balance Condition Method All 1-phase W ell-balanced HC model All 2-phase W ell-balanced HC model Unbalanced load Integrate with at bus 4 sequence power-flow All 3-phase W ell-balanced HC model Untransposed lines Integrate with sequence line model Unbalanced load Integrate with at bus 4 sequence load current Multi-phase W ell-balanced HC model Untransposed lines Integrate with sequence line model Unbalanced load Integrate with at bus 4 , 5 , 8 sequence load current T ABLE I: The evaluation of different scenarios for multi- phase setups. T o evaluate the multi-phase modeling, we test differ- ent scenarios based on IEEE 8 -bus unbalanced radial network in Fig. 6. The summary of solving methods is shown in T able I. If the system is well-balanced, meaning transposed lines and balanced loads, we simply apply the original model to compute hosting capacity . When the lines are unbalanced, meaning the admittance matrix in phase variables is unsymmetrical, we use sequence decoupled line model before applying the proposed method to positive sequence. Besides, the loads can also be unbalanced, wher e the load for each phase is individually specified instead of averaging from total demand. After the HC calculation via positive sequence power flow , the specified current injections represent such unbalance and hence the unbalanced phase volt- ages can be computed by nodal voltage equations. C. Parallel Computation Fig. 7: Hosting capacity of IEEE 123 -bus system subject to differ ent constraints. Fig. 7 shows numerical results for IEEE 123 -bus sys- tem. Similar to the 8 -bus system, it shows HC subject to different constraints and compares the theoretical HC with the numerically found one by using NLP solver Fmincon for the same model, verifying our theorem. Fig. 8: The IEEE 123 -bus network structure shows the computational complexity . W e implement the distributed optimization model for speed up. The different segmen- tations are marked with different colors. The benchmark time on heuristically solving HC is much longer and the solver is unstable because of the much more complicated structure shown in Fig. 8. T oo many variables and constraints cause computational difficulties wher e the solver ’s results can be unstable and non-converge. Applying the distributed optimiza- tion model (8), the segmentations in Fig. 8 are marked with different colors. Key nodes bus- 16 and bus- 73 are chosen to partition the model. 7 Fig. 9: The comparison of solving time for differ ent systems with different segmentations. The results with parallel computation are also shown in Fig. 7 as the solid curve, which is almost the same as the bar plot, verifying the distributed algorithm. The dramatic time decrease (Fig. 9) indicates we are able to reduce computational complexity and impr ove ef ficiency as we calculate all the segments simultaneously . V I . C O N C L U S I O N In this paper , we propose a hosting capacity problem formulation for better solvability . A geometric intuition of solving HC motivates us to find the globally op- timal solution even with non-convex constraints, with which we reveal the conditions of voltages to achieve the hosting capacity under various operational limits. Practical conditions for HC calculation, such as multi- phase unbalance and solving-time complexity , are an- alyzed via integrated methodologies. A mathematical proof is provided along with numerical validations to support the theoretical insights. Future work includes introducing more operation limits in a similar fashion and considering more complex modeling. R E F E R E N C E S [1] E. K. Bawan, “Distributed generation impact on power system case study: Losses and voltage profile,” Australasian Universities Power Engineering Conference , pp. 1–6, Sep 2012. [2] P . P . Barker and R. W . D. Mello, “Determining the impact of dis- tributed generation on power systems: part 1 - radial distribution systems,” Power Engineering Society Summer Meeting , pp. 1645– 1656, Jul 2000. [3] G. K. Ari and Y . Baghzouz, “Impact of high PV penetration on voltage regulation in electrical distribution systems,” International Conference on Clean Electrical Power , pp. 744–748, Jun 2011. [4] J. Smith, “Stochastic analysis to determine feeder hosting capacity for distributed solar PV,” the Electric Power Research Institute (EPRI), T ech. Rep., Dec 2012. [5] M. R ylander , J. Smith, and W . Sunderman, “Streamlined method for determining distribution system hosting capacity ,” IEEE T rans- actions on Industry Applications , pp. 105–111, Jan 2016. [6] L. K. L. Estorque and M. A. A. Pedrasa, “Utility-scale DG planning using location-specific hosting capacity analysis,” IEEE Innovative Smart Grid T echnologies , pp. 984–989, Nov 2016. [7] N. Baldenko and S. Behzadirafi, “Determination of photovoltaic hosting capacity on radial electric distribution feeders,” IEEE International Conference on Power System T echnology , pp. 1–4, Sep 2016. [8] “T utorial integration capacity analysis CYME 8.0,” Eaton, CYME International , Mar 2018. [9] S. N. Liew and G. Strbac, “Maximising penetration of wind generation in existing distribution networks,” IEE Proceedings - Generation, T ransmission and Distribution , pp. 256–262, May 2002. [10] L. Ochoa, C. Dent, and G. Harrison, “Distribution network capacity assessment: V ariable DG and active networks,” IEEE T ransactions on Power Systems , pp. 87–95, Feb 2010. [11] F . Capitanescu, L. F . Ochoa, H. Margossian, and N. D. Hatziar- gyriou, “Assessing the potential of network reconfiguration to improve distributed generation hosting capacity in active distri- bution systems,” IEEE T ransactions on Power Systems , pp. 346–356, Jan 2015. [12] F . Ebe, B. Idlbi, J. Morris, G. Heilscher , and F . Meier , “Evaluation of PV hosting capacity of distribution grids considering a solar roof potential analysis — comparison of different algorithms,” IEEE PowerT ech , pp. 1–6, Jun 2017. [13] A. Dubey , S. Santoso, and A. Maitra, “Understanding photovoltaic hosting capacity of distribution cir cuits,” IEEE Power Energy Soci- ety General Meeting , pp. 1–5, Jul 2015. [14] K. Coogan, M. J. Reno, S. Grijalva, and R. J. Broderick, “Locational dependence of PV hosting capacity correlated with feeder load,” IEEE PES T & D Conference and Exposition , pp. 1–5, Apr 2014. [15] C. Schwaegerl, M. H. J. Bollen, K. Karoui, and A. Y agmur , “V oltage control in distribution systems as a limitation of the hosting capacity for distributed energy resources,” International Conference and Exhibition on Electricity Distribution , pp. 1–5, Jun 2005. [16] M. S. S. Abad, J. Ma, D. Zhang, A. S. Ahmadyar , and H. Mar- zooghi, “Probabilistic assessment of hosting capacity in radial distribution systems,” IEEE T ransactions on Sustainable Energy , pp. 1935–1947, Mar 2018. [17] D. A. Sarmiento, P . P . V ergara, L. C. P . da Silva, and M. C. de Almeida, “Increasing the PV hosting capacity with OL TC technology and PV V Ar absorption in a MV/L V rural brazilian distribution system,” International Conference on Harmonics and Quality of Power , pp. 395–399, Oct 2016. [18] N. Etherden and M. H. J. Bollen, “Increasing the hosting capacity of distribution networks by curtailment of r enewable energy resour ces,” IEEE PowerT ech , pp. 1–7, Jun 2011. [19] P . Hasanpor Divshali and L. Sder, “Improving hosting capacity of rooftop PVs by quadratic control of an L V-central BSS,” IEEE T ransactions on Smart Grid , pp. 919–927, Jan 2019. [20] J. Seuss, M. J. Reno, R. J. Broderick, and S. Grijalva, “Improving distribution network PV hosting capacity via smart inverter reac- tive power support,” IEEE Power Energy Society General Meeting , pp. 1–5, Jul 2015. [21] L. J. Thomas, A. Burchill, D. J. Rogers, M. Guest, and N. Jenkins, “Assessing distribution network hosting capacity with the addi- tion of soft open points,” IET International Conference on Renewable Power Generation , pp. 1–6, Sep 2016. [22] S. N. Salih, P . Chen, O. Carlson, and L. B. Tjernberg, “Optimizing wind power hosting capacity of distribution systems using cost benefit analysis,” IEEE T ransactions on Power Delivery , pp. 1436– 1445, Jun 2014. [23] J. Zhao, H.-D. Chiang, H. Li, and P . Ju, “On PV -PQ bus type switching logic in power flow computation,” Power Systems Com- putation Conference , Jan 2008. [24] M. Abdel-Akher , K. M. Nor , and A. A. Rashid, “Improved three- phase power-flow methods using sequence components,” IEEE T ransactions on power systems , pp. 1389–1397, Aug 2005. [25] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power systems resear ch and education,” IEEE T ransactions on Power Systems , pp. 12–19, Feb 2011. [26] R. Zimmerman, C. Murillo-Sanchez, and D. Gan, “Matpower-a matlab power system simulation package: User ’s manual,” Dec 1997. [27] Y . W eng, R. Rajagopal, and B. Zhang, “A geometric analysis of power system loadability regions,” IEEE T ransactions on Smart Grid , 2019. A P P E N D I X A I M P L E M E N T O N 3 - B U S T O Y E X A M P L E For the simple setup to obtain geometric understand- ing, voltage angles and line reactances are set to zeros. 8 In such a setup, the HC problem model simplifies to max n X i =1 P ing i s . t . P i = n X k =1 | V i || V k | G ik , V min ≤ | V i | ≤ V max . (9) A P P E N D I X B H O S T I N G C A PA C I T Y S U B J E C T T O C O N S T R A I N T S A. Pr oof of Theorem 1: Solution Adjustment According to V oltage Magnitude Constraints Proof. For the proof, it is sufficient to show that the proposed hosting value H C ∗ is larger than arbitrary summation of power generation. Specifically , we need to show that n X i =1 P i, solution ≥ n X i =1 P i , (10) where P n i =1 P i, solution repr esents the solution H C ∗ ac- cording to the theorem. P n i =1 P i repr esents the total power calculated by any setup of voltages within limits. Expanding the left-hand side of (10) with the voltages in the theorem n X i =1 P i, solution = n · | G ik | · (1 . 05 − 0 . 95) 2 (11a) ≥ n X i =1 | G ik | · ( V i − V k ) 2 (11b) = n X i =1 P i , (11c) where (11a) is the direct plug-in, (11b) is due to extreme numbers, and (11c) is due to power flow equations. B. Pr oof of Theorem 2: Solution Adjustment According to Angle Difference Constraint Proof. W ith the general complex impedance setup, there are two types of functionals, namely sinusoids poly- nomials in the polar coordinate. The multiplication of these two functionals in the power flow equation make the HC optimization difficult to analyze for the global optimum. Therefore, we transform the problem into the rectangular coordinates to reduce the function types into one, namely , the polynomial functions. max n X i =1 P i s . t . P i = − X k ∈ N i G ki · ( V i, real 2 + V i, imag 2 ) + X k ∈ N i ( G ki V k, real − B ki V k, imag ) · V i, real + X k ∈ N i ( V k, real B ki + V k, imag G ki ) · V i, imag , V min ≤ | V i | ≤ V max , θ min ≤ θ i ≤ θ max , (12) where bus k is the neighbor bus of bus i . Extend objective function for quadratic optimization [27], we obtain N X i =1 P i = N X i =1 N X k =1 ( − G ik ) · [( V i, real − V k, real ) 2 + ( V i, imag − V k, imag ) 2 ] = N X i =1 N X k =1 ( − G ik ) · ( V i cos θ i − V k cos θ k ) 2 + ( − G ik ) · ( V i sin θ i − V k sin θ k ) 2 . (13) As the multiplication is the key component in the quadratic optimization, we can simplify the objective further by having | V i | = a, | V k | = b, cos θ i = a 1 , sin θ i = a 2 , cos θ k = b 1 , sin θ i = b 2 . Suppose bus i is the even bus and bus k is the odd bus. Such notation reshape the objective into max n X i =1 ( − G ik ) · h ( aa 1 − bb 1 ) 2 + ( aa 2 − bb 2 ) 2 i (14) s . t . 0 . 95 ≤ a ≤ 1 . 05 , 0 . 95 ≤ b ≤ 1 . 05 , (15) a 2 1 + a 2 2 = 1 , b 2 1 + b 2 2 = 1 . (16) Substitute the equality constraint (16) into the objective function (14) and we get renewed problem formulation, max n X i =1 ( − G ik ) · [ a 2 + b 2 − 2 ab · cos θ ik ] s . t . 0 . 95 ≤ a ≤ 1 . 05 , 0 . 95 ≤ b ≤ 1 . 05 . (17) In the distribution grids, G ik is mostly negative, mak- ing − G ik positive. T o find the maximum HC, we show in the following that we can find the global optimum in a non-convex problem by piecewise analysis of 2 ab · cos θ ik according to θ . 1) When θ ∈ [0 , θ max ] , for a fixed θ max ∈ [ π , 2 π ] : θ ∈ [0 , θ max ] means that − 1 ≤ cos θ ik ≤ 1 . As a and b are positive voltage magnitudes, we obtain the maximal HC, when cos θ ik = − 1 . Therefore, one element of the objective function becomes a 2 + b 2 − 2 ab · cos θ ik = ( a + b ) 2 . The maximum is achieved when a = b = 1 . 05 , θ ik = π , meaning hosting capacity is under the condition that all the bus voltage magnitudes are the highest and the 9 differ ence between neighbor bus voltage angles is π . Therefore, we need to change the voltage angle of π within one branch interval, impossible in the real system. Meanwhile, the unconstrained voltage angle is unstable with observation from several attempts, so the next step is to narrow down the voltage angle. 2) When θ ∈ [0 , θ max ] , where arccos( V u + V l 2 V u ) ≤ θ max ≤ π : W ith voltage angle constraint, cos θ max ≤ cos θ ik ≤ 1 . As a and b are positive, cos θ ik should be as small as possible to get a bigger value of the objective, which is cos θ max . T o maximize a 2 + b 2 − 2 ab · cos θ max , the voltage pattern is the same as above ( a = b = 1 . 05 ), when cos θ max is a negative or small positive number . Until θ max keeps decr easing to a specific number , the solution pattern changes to “high-low voltage” as in Theor em 1. For calculating this critical value, V 2 u + V 2 u − 2 · V u V u cos θ max ≤ V 2 u + V 2 l − 2 V u V l cos θ max , θ ≤ arccos( V u + V l 2 V u ) , (18) where V u and V l are the upper and lower bounds. we plug in the boundary values of a and b below , θ ≤ 0 . 3098 . (19) Therefor e, the solution for this angle range is a = b = 1 . 05 , θ ik = θ max . 3) When θ ∈ [0 , θ max ] , where θ max < arccos( V u + V l 2 V u ) : As we mention above, the solution is a = 1 . 05 , b = 0 . 95 , θ ik = θ max , when voltage angle is limited below the critical value arccos( V u + V l 2 V u ) . The “high- low voltage” pattern appears again. C. Pr oof of Theorem 3: Solution Adjustment According to Thermal Limit In addition to the voltage constraints (for | V | and θ ), the thermal limit is rather important in two-way power flow . W e use current flow limit to represent this constraint that results from the heating effects of devices, I Lmin ≤ I ik ≤ I Lmax = − I Lmin = C . W e do a comparison of the previous solution and the thermal limit to estimate if it has the impact on the final solution. Proof. W e will plug in the voltage solutions from Theo- rem 3 to check if the new bounds are violated or not. If not, the solution for HC stays the same. If yes, we will find the new voltage solutions. Specifically , we plug in the voltage solutions into the current flow below I ik = | G ik + j B ik | · ( V i − V k ) . (20) Now , compar e the calculated current to the new boundary constraints. • If − C ≤ I ∗ ik ≤ C , meaning the thermal limit has no impact on the solution; • Otherwise, we plug in (20) into (2e) to find the solution. | G ik + j B ik | 2 · | ( V i, real − V k, real ) + j ( V i, imag − V k, imag ) | 2 ≤ C 2 , (21) equivalent to [( V i, real − V k, real ) 2 + j ( V i, imag − V k, imag ) 2 ] ≤ C 2 / ( G 2 ik + B 2 ik ) . (22) The quadratic functions on the left-hand side are exactly one component in (6) to represent the objec- tive P n i =1 P i . Therefor e, the power flow impacted by thermal limit is P ∗ ik = ( − G ik ) · [( V i, real − V k, real ) 2 + ( V i, imag − V k, imag ) 2 ] . (23) W e calculate the maximum value that this compo- nent can achieve via P ∗ ik max = ( − G ik ) · C 2 / ( G 2 ik + B 2 ik ) . (24) T o obtain the solution pattern (values of | V i | , | V k | and θ ik ), we solve this thr ee-variable-cubic equation. a 2 + b 2 − 2 ab · cos θ ik = C 2 / ( G 2 ik + B 2 ik ) , (25) where a and b represent | V i | and | V k | . Due to the coupling of variables in the power flow equations, the changing solution pattern of one branch can have the impact on the rest buses to achieve maximum generation. The voltage pattern might transpose from “high-low” to “low-high”. According to the previous theorems and proofs, the “high-low voltage” pattern in one branch is ob- tained because of the quadratic function type in (23). Therefor e, such transpose will not change the part of HC uncorrelated to this branch except the buses connected to it. In other words, only the quadratic functions containing the involved variables follow this change. D. Pr oof of Solution Adjustment According to Power Factor Constraint Limiting Reactive power of PV generator is to limit the power factor (2f). P | S | ≥ η , (26) Suppose we use η = 0 . 95 as the operation limit, P 2 P 2 + Q 2 ≥ η 2 , Q 2 ≤ 1 − η 2 η 2 P 2 , Q ≤ p 1 − η 2 p η 2 P f or Q > 0 , Q ≥ − p 1 − η 2 p η 2 P f or Q < 0 , (27) 10 A P P E N D I X C E X T E N D T O T H R E E - P H A S E U N B A L A N C E D N E T W O R K W e simplify the analysis of multi-phase system by transferring the model into sequence form via symmet- rical component transformation matrix. T =   1 1 1 1 α 2 α 1 α α 2   where α = 1 ∠ 120 ◦ , (28) a set of sequence variables is obtained, i.e. | V 012 i | = [ | V 0 i | , | V 1 i | , | V 2 i | ] T = T − 1 | V abc i | . (29) The sequence admittance matrix is Y 012 = T − 1 Y abc T . (30) In order to apply the unbalanced situation, we decouple the sequence model while using current injection from negative and zero sequence to represent the unbalance. I m =  S load V m  ∗ + X k ∈N ∆ I m k , m = 0 or 2 . (31) According to [24], the positive sequence voltage mag- nitude is much larger under unbalanced condition and the power flow is similar to single-phase formulation. Therefor e, the positive sequence power flow equations are given as P 1 i = n X k =1 | V 1 i || V 1 k | T ( G 1 ik cos θ 1 ik + B 1 ik sin θ 1 ik ) , (32a) Q 1 i = n X k =1 | V 1 i || V 1 k | T ( G 1 ik cos θ 1 ik − B 1 ik sin θ 1 ik ) . (32b) The equations allow us to implement the proposed method and obtain voltage pattern solutions for positive sequence. Based on the decoupled models, we can com- pute the negative and zero sequence voltages via nodal voltage equations. Y m V m = I m , m = 0 or 2 . (33) Combining with results of positive sequence power flow , it’s able to transfer back to three phase solution via (29). A P P E N D I X D E N A B L E P A R A L L E L C O M P U TAT I O N F O R L A R G E S C A L E Following the distributed optimization model in IV . B, we derive the objective HC t = s t X i = s t − 1 s t X k = s t − 1 ( − G ik ) · [( V i, real − V k, real ) 2 + ( V i, imag − V k, imag ) 2 ] ( s 0 = 1) = s j X i = s j − 1 s j X k = s j − 1 ( − G ik ) · ( V i cos θ i − V k cos θ k ) 2 + ( − G ik ) · ( V i sin θ i − V k sin θ k ) 2 , which shows the same pattern as (6). Thus, The capabil- ity of parallel computation with the proposed model is proved.