Statistical Modeling of Networked Solar Resources for Assessing and Mitigating Risk of Interdependent Inverter Tripping Events in Distribution Grids

SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 1 Statistical Modeling of Network ed Solar Resources for Assessing and Mitigating Risk of Interdependent In v erter T ripping Events in Distrib ution Grids Kav eh Dehghanpour , Member , IEEE, Y uxuan Y uan, Student Member , IEEE, Fankun Bu, Student Member , IEEE Zhaoyu W ang, Member , IEEE Abstract —It is speculated that higher penetration of in verter - based distributed photo-voltaic (PV) power generators can in- crease the risk of tripping events due to voltage fluctuations. T o quantify this risk utilities need to solve the interactive equations of tripping events f or networked PVs in r eal-time. However , these equations are non-differentiable, nonlinear , and exponen- tially complex, and thus, cannot be used as a tractable basis for solar curtailment prediction and mitigation. Furthermore, load/PV power values might not be a vailable in real-time due to limited grid observability , which further complicates tripping event pr ediction. T o address these challenges, we hav e employed Chebyshev’ s inequality to obtain an alternativ e probabilistic model for quantifying the risk of tripping for networked PVs. The proposed model enables operators to estimate the probability of interdependent inverter tripping ev ents using only PV/load statistics and in a scalable manner . Furthermore, by integrat- ing this probabilistic model into an optimization framework, countermeasures ar e designed to mitigate massive interdependent tripping events. Since the proposed model is parameterized using only the statistical characteristics of nodal active/r eactive powers, it is especially beneficial in practical systems, which hav e limited real-time observability . Numerical experiments have been performed employing real data and feeder models to verify the performance of the proposed technique. Index T erms —Probabilistic modeling; po wer statistics; risk assessment; tripping events; I . I N T RO D U C T I O N Increasing penetration of distributed energy resources (DERs), including in verter-based photo-voltaic (PV) power generators, in distribution grids represents opportunities for enhancing system resilience and customer self-sufficiency , as well as challenges in grid control and operation. One of these challenges is the potential increase in the risk of tripping of in verter -based resources due to undesirable fluctuations in the grid’ s voltage profile [1]. This can put a hard limit on the feasible capacity of operational PVs in distribution grids, re- duce the economic value of rene wable resources for customers, and cause loss of service in stand-alone systems [2], [3]. The possibility of DER po wer generation disruption due to v oltage- related vulnerabilities in unbalanced distrib ution grids has been discussed in the literature: in [4], [5], risk of interdependent tripping of PVs, with ON/OFF current interruption mechanism This work was supported by the Advanced Grid Modeling Program at the U.S. Department of Energy Office of Electricity under Grant de-oe0000875. ( Corr esponding author: Zhaoyu W ang ) K. Dehghanpour, Y . Y uan, F . Bu, and Z. W ang are with the Department of Electrical and Computer Engineering, Iowa State Univ ersity , Ames, IA 50011 USA (e-mail: kavehd@iastate.edu; wzy@iastate.edu). was demonstrated numerically in a distribution grid test case for the first time. It was shown that the unbalanced and re- sistiv e nature of networks can further exacerbate this problem by causing positi ve inter-phase voltage sensitivity terms that act as destabilizing positive feedback loops, leading to voltage deviations after tripping of an individual in verter . The impact of grid v oltage sensiti vity on DER curtailment was also studied and observed in [2]. Based on these insights, guidelines were provided in [6] to roughly estimate the impact of new DER capacity connections on the maximum v oltage deviations in the grid. It was sho wn in [7] that very large or small number of in verter-based resources in distribution systems can lead to interdependent failure ev ents that contribute to voltage collapse in transmission le vel. Detailed realistic numerical studies were performed on practical feeder models in [8]– [13] that corroborated the considerable impacts of extreme PV integration lev els, and in verter control modes on grid voltage fluctuations, which is the critical factor in causing massi ve solar curtailment scenarios. Most existing works relied on scenario-based simulations and numerical studies to capture the likelihood of inv erter tripping under high renewable penetration. While this has led to useful guidelines and inv aluable intuitions, it falls short of providing a generic theoretical foundation for predicting and containing tripping ev ents. Specifically , the dependencies between nodal solar power distributions, nodal voltage pro- files, and inv erter tripping e vents hav e not been explicitly analyzed in the literature thus far . These dependencies are influenced by in verter protection settings and gov erned by a set of networked power-flo w-based equations, which turn out to be non-differentiable and nonlinear . In this regard, sev eral fundamental challenges have not been addressed: (1) Lack of scalability: Solving the inv erter tripping equations directly in real-time requires a large-scale search process to explore almost all the joint combinations of “ON/OFF” configurations for the in verters. The computational complexity is due to the interactive and networked nature of tripping ev ents, meaning that the states of in verters influence each other and are not independent [14]. The source of interdependency in chances of in verter tripping is the dependencies in nodal voltages of power grid (i.e., disruption of power injection at one node impacts nodal voltages of other neighboring in verters, which in turn could influence their probability of tripping.) For example, tripping of an in verter (or a cluster of in verters) leads to a change in loading distributions, which SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 2 can most likely increase/decrease the chance of tripping for other inv erters during under/over -voltage scenarios, especially in weak grids. This interdependency pre vents the solver from decoupling the tripping equations into separate equations for individual in verters. Thus, the scale of search for finding the correct configuration increases exponentially ( 2 N ) with the number of in verters ( N ). Another factor that contributes to computational complexity is the volatility of PV power , which forces the solver to explore, not only various tripping configurations, but also numerous solar scenarios at granular time steps. (2) Limited tractability for mitigation: A direct solution strategy for tripping equations cannot be easily inte- grated into optimization-based decision models, since it has no predictiv e capability and cannot be use to answer what if queries, unless a thorough expensi ve search is performed over all possible future load/PV scenarios. Also, due to their non- differntiability , integrating the tripping equations into decision models complicates formulation by adding integer variables to the problem. (3) Limited access to online data: Practical distribution grids have lo w online observability , meaning that the values of real-time nodal power injections can be unknown in real-time for a large number of PVs/loads due to commu- nication time delays or limited number of sensors. Thus, we might not hav e access to sufficient online information to solve the tripping problem directly . T o tackle these challenges, we propose an alternative proba- bilistic modeling approach to quantify and mitigate the risk of voltage-dri ven tripping ev ents. Instead of complex scenario- based look-ahead search ov er numerous possible tripping configurations, our methodology is built upon probabilistic manipulation of power flow equations in radial networks to estimate the probability of inv erter tripping using only the av ailable statistical properties of loads/PVs. Interdependent Bernoulli random variables are used to model probabilities of in verter tripping and capture their mutua. These probabilities are voltage-dependent and serve as unknown micr o-states in the equations of tripping e vents. Then, Chebyshev’ s inequality [15] is applied to determine a stationary lo wer bound for the values that these micro-states can assume under any probable nodal power injection scenarios. This lower bound provides a conservati ve estimation of expected PV curtailment, and thus, represents a statistical risk metric for tripping e vents. Furthermore, due to its simple matrix-form and differentiable structure, the proposed probabilistic model can be con v eniently integrated into an optimization framework as a constraint, which enables mitigating unwanted solar curtailment events by designing optimal voltage regulation countermeasures. The proposed methodology is generic and can capture the behavior of arbitrary radial distribution feeders using only load/PV statistics and network topology/parameters. This implies that tripping ev ents can be conservati vely predicted using the proposed model and without the need for online access to granular PV/load data or expensi ve scenario-based search process, which makes our strategy specifically suitable for practical networks. Numerical experiments ha ve been performed using real ad- vanced metering infrastructure (AMI) data and feeder models from our utility partners to validate the developed probabilistic 𝑞 𝑖 𝑝 𝑗 𝑞 𝑗 𝑝 𝑘 𝑞 𝑘 𝑝 𝑖 𝑞 𝑖 𝑝 𝑗 𝑞 𝑗 𝑝 𝑘 𝑞 𝑘 𝑣 𝑗 𝑣 𝑘 𝑠 𝑖 ∈ { 1 , 0 } 𝑠 𝑗 ∈ { 1 , 0 } 𝑠 𝑘 = { 1 } 𝑝 𝑖 a b 𝒓 𝒊 𝒋 + 𝒋 𝒙 𝒊 𝒋 𝒓 𝒋 𝒌 + 𝒋 𝒙 𝒋 𝒌 𝑣 𝑖 c Fig. 1. Distribution feeder structure with PVs, loads, and voltage-sensiti ve current interruption mechanisms (i.e., switches). framew ork. The numerical validates the performance of the probabilistic model for both over - and under-v oltage scenarios, and show that ignoring the possibility of tripping in voltage regulation can exacerbate voltage deviations. I I . D E R I V I N G A C O N S E RV A T I V E P R O BA B I L I S T I C M O D E L O F P V T R I P P I N G E V E N T S In this section, we will dev elop and then parameterize a probabilistic model of networked in verter -based PVs to quantify the possibility of emergent tripping ev ents. T o do this, first, we begin with the original model of in verter tripping with ON/OFF voltage-dri ven current interruption mechanism, and then, we will show that by adopting a probabilistic approach tow ards the original model and using Chebyshev’ s inequality , tripping probabilities can be conserv ativ ely estimated using the statistical properties of nodal av ailable load/PV po wer . A. Original Interactive Switching Equations In this paper , it is assumed that PV resources are protected against voltage deviations using ON/OFF switching mecha- nisms. Note that here a “switch” can be a mechanical relay , as well as a non-physical in verter control function that stops current injection into the grid under abnormal voltage ev en if the inv erter is still physically connected to the grid [16]. The PV is tripped in case the nodal voltage deviates from a user- defined permissible range, [ V min , V max ] . In this paper, this range is adopted from the literature [4], as V min = 0 . 9 p.u. and V max = 1 . 1 p.u. . The switching mechanisms are simply modelled as binary micr o-state v ariables with the follo wing voltage-dependent function (see Fig. 1): s i ( t ) =      1 V min ≤ V i ( t ) ≤ V max 0 V i ( t ) < V min 0 V i ( t ) > V max (1) where, s i ( t ) is the micro-state assigned to the i ’th PV at time t as a function of the in verter node’ s v oltage magnitude V i . Here, s i ( t ) = 1 implies ON and s i ( t ) = 0 indicates OFF . The assumption in this switching model is that over long enough time intervals the impact of in verter dynamics, e.g., ride- through capabilities [17]–[20], can be conservati vely ignored. This assumption considerably enhances the tractability of the model at the expense of loss of accuracy . In this sense, the SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 3 switching model is a worst-case representation of inv erter tripping. Since the approximate power flow equations for distribution grids are linear with respect to the squared v alues of nodal voltage magnitudes [21], we re-write equation (1) using a v ariable transformation, v i = V 2 i , and employing unit step functions as follows: s i ( t ) = U ( v i ( t ) − v min ) − U ( v i ( t ) − v max ) (2) where, v min = V 2 min , v max = V 2 max , and the unit step function U ( · ) is defined as follows: U ( x ) = ( 1 x ≥ 0 0 x < 0 , (3) Note that in verters’ micro-states are influenced by nodal voltages and are thus highly interdependent on each other, as changes in the state of one switch will cause nodal power variations, which leads to a change of voltage at other nodes that can in turn influence probability of tripping e vents. T o obtain the overall governing equations of inv erter tripping, the mutual impacts of switch micro-states on each other are captured using an approximate unbalanced power flo w model for radial distribution grids [21], which determines voltage at node i as a function of activ e/reactive po wer injections of ev ery other node in a grid (with a total of N + 1 nodes): v i ( t ) = N X j =1 ˜ v ij + v 0 , ∀ i ∈ { 1 , ..., N } (4) where, v 0 = V 2 0 , with V 0 denoting the voltage magnitude at a grid reference bus, and the intermediary variable ˜ v ij represents the impact of activ e/reactive po wer injection at node j on v i , which is obtained as follo ws: ˜ v ij = R ij ˜ p j ( t ) + X ij ˜ q j ( t ) (5) where, R ij and X ij are the aggregated series resistance and reactance v alues corresponding to the intersecting branches in the paths connecting nodes i and j to the reference bus calculated as follo ws [21]: R ij = 2 X { n,m }∈ P a ( i,j ) r nm (6) X ij = 2 X { n,m }∈ P a ( i,j ) x nm (7) where, P a ( i, j ) represents the set of pairwise nodes consisting of the neighboring nodes that are on the intersection of the unique paths connecting nodes i and j to the reference b us; r nm and x nm denote the real series resistance and reactance of the branch connecting nodes n and m . Also, ˜ p j and ˜ q j denote the activ e and reactiv e po wer injections at bus j , which are in turn determined by the micro-state of the PV at node j (see Fig. 1): ˜ p j ( t ) = p j ( t ) s j ( t ) (8) ˜ q j ( t ) = q j ( t ) s j ( t ) (9) with p j and q j representing the av ailable load/PV power at node j , where p j > 0 implies generation. Equations (4)-(7) are obtained in vector form for all three phases of unbalanced distribution grids [21]. Equations (2)-(9) fully determine the states of networked PVs. The difficulty in solving these equations is due to three factors: (I) the size of solution space increases exponentially as the number of micro-states { s 1 , ..., s N } grows. Since these micro-states are not independent and influence each other in complex and non-trivial ways they cannot be obtained indi- vidually , and a thorough search process is needed to explore all possible switching configurations. This can be extremely expensi ve and impossible to scale to large systems with high population of inv erters. (II) Due to the discrete step functions in (2), tripping equations are nonlinear and non-differentiable. This contributes to problem dif ficulty since gradient-based methods cannot be applied. (III) p j and q j act as time-varying input parameters within the model. This implies that using the tripping equations for predicting probability of tripping ev ents requires extensiv e search process to cov er all probable PV/load time-series scenarios. This expensi ve search process hinders the tractability of optimization-based frameworks for designing tripping mitigation strategy . Not all the nodes in the tripping model are necessarily con- trolled by ON/OFF voltage-sensiti ve switching mechanisms. For examples, ordinary load nodes are generally not governed by equation (2). In this paper , for the sake of brevity , the switching equations are still written for all the nodes in the grid as presented, howe ver , we will simply assign constant values, s i ( t ) = 1 , ∀ t to the nodes without ON/OFF control and remov e their corresponding switching from the equations (see Fig. 1). B. Alternative Approximate Pr obabilistic Model W e adopt a probabilistic point of view towards tripping model. This allows us to obtain a stationary differentiable statistical model that has a simple matrix-form formulation. Accordingly , the ON/OFF current interruption mechanisms, s i ’ s, are modelled as random variables following Bernoulli probability distributions with parameters λ i , ∀ i ∈ { 1 , ..., N } : s i ∼ B ( λ i ) , where parameter λ i is defined as the probability of the i ’th in v erter switch being ON, λ i ( t ) = P r { s i ( t ) = 1 } . The goal is to transform micro-states from discontinuous binary variables ( s i ∈ { 0 , 1 } ) into continuous variables ( λ i ∈ [0 , 1] ). T o rewrite the equations in terms of new micro-states note that we have E { s i ( t ) } = λ i ( t ) for Bernoulli probability distributions, where E {·} represents the expectation operation. Thus, by performing an expectation operation ov er both sides of (2), probability of in verter tripping in terms of the new micro-states can be obtained as follows: λ i ( t ) = P r { v min ≤ v i ( t ) ≤ v max } (10) where, we hav e exploited E { U ( f ( x )) } = P r { f ( x ) ≥ 0 } . Note that the probability of tripping for an in verter is an implicit function of nodal voltage probability distribution, which in turn is influenced by the states of other inv erters. Due to the interconnected nature of the problem, no inde- pendency assumptions has been made on random variables λ i , ∀ i ∈ { 1 , ..., N } . Howe ver , the exact distributions of SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 4 nodal voltages are unknown and complex functions of nodal activ e/reactiv e injections, which implies that (10) cannot be determined analytically unless over -simplifying assumptions are made. Instead, we employ Chebyshe v’ s inequality [15] to provide a lower bound on micro-state as a function of nodal voltage statistics without making an y assumption on v oltage distributions, P r { v min ≤ v i ( t ) ≤ v max } ≥ 1 − σ 2 v i + ( µ v i − v max + v min 2 ) 2 ( v max − v min 2 ) 2 (11) where, σ 2 v i and µ v i are the variance and mean of v i , re- spectiv ely . Hence, the approximate probabilistic model can be formulated for each micro-state as follows: ˆ λ i ( t ) = 1 − σ 2 v i + ( µ v i − v max + v min 2 ) 2 ( v max − v min 2 ) 2 (12) This ne w tripping model has two features: (1) it is a conservati ve estimator of the original system since it ov er- estimates the probability of inv erter tripping, ˆ λ i ≤ λ i . (2) As will be shown in Section II-C, the approximate probabilistic model can be con veniently parameterized in terms of nodal av ailable activ e/reactiv e po wer statistics. Hence, as long as certain statistics are kno wn (or estimated), the model allows us to accurately track probability of inv erter tripping without running time-series simulations under numerous scenarios. C. Pr obabilistic Model P arameterization T o parameterize the alternati ve tripping model (12), nodal voltage statistics, σ 2 v i and µ v i , are obtained in terms of nodal av ailable activ e/reactive power statistics. T o do this, power flow/injection equations (4)-(9) are lev eraged. Stage 1: µ v i P arameterization - The expected value of voltage magnitude squared is determined using (4)-(5) as, µ v i = N X j =1 E { ˜ v ij } + v 0 = N X j =1 ( R ij E { ˜ p j } + X ij E { ˜ q j } ) + v 0 (13) T o calculate E { ˜ p j } and E { ˜ q j } , we will first obtain their cumulativ e distribution functions (CDFs) [15], F ˜ p j and F ˜ q j , respectiv ely . This process is shown for ˜ p j as follows ( F ˜ q j is obtained similarly): F ˜ p j ( P ) = P r { ˜ p j ( t ) ≤ P } = (1 − λ j ( t )) U ( P ) + λ j ( t ) F p j ( P ) (14) The rational behind (14) is that the distribution of po wer injection is determined by two functions: the distribution of PV switch (which is ON with probability λ j ( t ) ), and the CDF of av ailable PV power , F p j . Now , the probability density functions (PDF) of the realized active nodal power injection, f ˜ p j , can be calculated as a function of the av ailable active solar power , f p j (a similar operation is performed for reactive power): f ˜ p j ( P ) = d F ˜ p j ( P ) d P = (1 − λ j ( t )) δ ( P ) + λ j ( t ) f p j ( P ) (15) Then, using the activ e/reactiv e po wer injection PDFs, E { ˜ p j } and E { ˜ q j } , can be obtained through integration: E { ˜ p j } = Z + ∞ −∞ αf ˜ p j ( α )d α = λ j P j (16) E { ˜ q j } = Z + ∞ −∞ β f ˜ q j ( β )d β = λ j Q j (17) where, P j and Q j denote the mean values of the av ailable ac- tiv e and reactiv e powers at node j , respectiv ely ( P j = E { p j } and Q j = E { q j } ). Thus, the mean nodal voltage magnitude squared can be written in terms of in verter switch statistics and expected PV/load av ailable powers: µ v i = N X j =1 { R ij λ j ( t ) P j + X ij λ j ( t ) Q j } + v 0 (18) Stage 2: σ 2 v i P arameterization - Using (4), the v ariance of nodal voltage magnitude squared can be formulated as, σ 2 v i = N X j =1 σ 2 ˜ v ij + 2 X 1 ≤ k : ˆ λ ˆ λ ˆ λ ( t ) = a 0 a 0 a 0 + B ˆ λ ˆ λ ˆ λ ( t ) +    ˆ λ ˆ λ ˆ λ ( t ) > C 1 ˆ λ ˆ λ ˆ λ ( t ) . . . ˆ λ ˆ λ ˆ λ ( t ) > C N ˆ λ ˆ λ ˆ λ ( t )    (35) where, all the time-in variant parameters of the model are con- catenated into the vector a 0 a 0 a 0 and matrices B , and { C 1 , ..., C N } . The elements of these parameters are determined by or ganizing the previous deriv ations in Stages 1 and 2, as follows: a 0 a 0 a 0 ( i ) = 1 − ( 2 v 0 − v max − v min v max − v min ) 2 (36) B ( i, j ) = − 1 ( v max − v min 2 ) 2 Γ 1 ij − 2 v 0 − v max − v min ( v max − v min 2 ) 2 ( P j R ij + Q j X ij ) (37) C i ( j, k ) = ( − 1 ( v max − v min 2 ) 2 Γ 1 ij k j 6 = k 0 j = k , (38) where, a 0 a 0 a 0 ( i ) denotes the i ’th element of a 0 a 0 a 0 , and B ( i, j ) and C i ( j, k ) are the ( i, j ) ’th and ( j, k ) ’ th elements of B and C i , respectiv ely . The aggregate switching equation can be written as a function of approximate macr o-state , ˆ S = P N i =1 ˆ λ i , as follows: ˆ S ( t ) = [ N X i =1 a 0 a 0 a 0 ( i )] + [ N X i =1 B ( i, :)] · ˆ λ ˆ λ ˆ λ ( t ) + ˆ λ ˆ λ ˆ λ ( t ) > [ N X i =1 C i ] ˆ λ ˆ λ ˆ λ ( t ) (39) where, ˆ S is a conservati ve estimator of the real macro-state, S , which is the actual expected population of inv erter that are ON, i.e., ˆ S ( t ) ≤ S ( t ) . Also, B ( i, :) is the i ’th row of matrix B . T o summarize, the proposed approximate probabilistic model lev erages available load/PV power statistics shown in T able I. Previous works hav e used v arious data-dri ven and machine learning methods that can be applied for obtaining statistical properties of nodal load/PV powers in partially observable networks from limited av ailable data (for example see [22]– [24]). Also, although the micro-states in the probabilistic model are random variables, the model itself is governed by deterministic functions of load/PV statistics. A related problem in distribution grids that testifies to the dependent nature of tripping is known as sympathetic tripping of in verters in weak grids [14], [25]: ov erloading/faults on one feeder can trigger the voltage protection mechanism of inv ert- ers on a healthy neighboring feeder . The sympathetic tripping of inv erters is also caused by dependencies in nodal voltages within the distribution grid (i.e., excessiv e load/fault current on one feeder contributing to voltage drops on other nodes). While sympathetic tripping is not exactly what the proposed statistical model in this paper captures, it still provides further support that dependency in tripping is possible in practice. SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 6 T ABLE I N E ED E D S TA T I S TI C S F O R D E V E LO P I N G T H E P RO P O SE D M O DE L D. Discussion on Pr obabilistic T ripping Model Pr operties The probabilistic model (35) repres ents a set of self- consistent equations; in other words, any ˆ λ ˆ λ ˆ λ that satisfies these equations is a conservati ve estimator of probability of in verter tripping. Furthermore, this probabilistic model can be thought of as the asymptotic equilibrium of an abstract discr ete dynamic system : ˆ λ ˆ λ ˆ λ ( k + 1) = a 0 a 0 a 0 + B ˆ λ ˆ λ ˆ λ ( k ) +    ˆ λ ˆ λ ˆ λ ( k ) > C 1 ˆ λ ˆ λ ˆ λ ( k ) . . . ˆ λ ˆ λ ˆ λ ( k ) > C N ˆ λ ˆ λ ˆ λ ( k )    (40) where, the equilibrium is achieved at ˆ λ ˆ λ ˆ λ ( k + 1) = ˆ λ ˆ λ ˆ λ ( k ) and coincides with the solution of the proposed probabilistic model (35). This abstract dynamic system has an intuitiv e interpretation: matrix B represents the linear component of the dynamics, which as can be observed in (37) and (29), is de- termined only by each individual nodes’ activ e/reactive power statistics, including the expected values and self-correlation between active/reacti ve power at each node alone. Howe ver , matrices { C 1 , ..., C N } capture the nonlinear components of the dynamic system, where the element C i ( j, k ) determines the coef ficient assigned to the interacti ve nonlinear probability- product term ˆ λ j ( t ) · ˆ λ k ( t ) in dri ving ˆ λ i ( t + 1) . In other words, C i ( j, k ) quantifies the mutual impact of the j ’th and k ’ th PV micro-states on dynamics of the i ’th switch. Furthermore, as observed in (38) and (33) the elements of C i , unlike B , are determined by the mutual correlations in av ailable active/reacti ve powers of different PVs. The inherent nonlinearity of (40) hints at the possibility of stage transition and bifurcation at equilibrium of the abstract dynamic system as PV/load power statistics ev olve ov er time, which could potentially result into a cascading tripping ev ent, as pointed out in [4], [7], [18], [19]. A regime shift at the equilibrium of the abstract nonlinear dynamic system basically corresponds to qualitativ e changes in the solution of our probabilistic tripping model, potentially , leading to a sudden increase in the av erage chances of voltage-dri ven tripping events caused by the growing penetration of solar energy in the system. In this sense, the structure of the abstract dynamic model is similar to other complex interactive dynamic systems in the literature, including nonlinear combinatorial evolution models [26] and asymmetric Ising systems [27], which are also kno wn to demonstrate critical behavior and emergent non-trivial patterns at the macro-le vel under certain conditions. An important factor in tripping studies is the impact of setting of inv erter protection systems. This can be seen in equations (36), (37), and (38) that present the parameters of the proposed statistical model. Specifically , parameter C i ( j, k ) in (38), captures the joint impacts of inv erter j and in verter k on probability of tripping for in verter i . As can be seen, the absolute v alue of this parameter decreases with 1 ( v max − v min ) 2 . Thus, increasing the in verter upper protection threshold, v max , or decreasing the lower protection threshold, v min , (i.e., mak- ing the in verter less sensitive to voltage events) will result in a decline in mutual impacts of in verters on each other . In other words, relaxing the protection boundaries significantly weakens the dependencies among in verter tripping. If C i ( j, k ) is thought of as a measure of strength of interdependency among in verters, then our model suggests that loss of inter- dependency is approximately proportional to the inv erse of in verter protection dead-band width squared. E. Inte grating V oltage-Dependent Resour ces Into the Pr o- posed Pr obabilistic T ripping Model Note that so far we hav e assumed that the nodal activ e and reactiv e power injections, p j and q j , are external inputs to the model. Howe ver , active and reactiv e power injection of certain nodes can show high lev els of voltage-dependency and cannot be treated as external inputs. The voltage-dependency can be caused by reacti ve power support from the in verters or load power voltage-sensitivity . In this section, we will demonstrate that voltage-dependent resources can also be included in our probabilistic model. T o do this, the active/reacti ve power injections are linearized around the nominal squared voltage ( v n ): p j ( v j ) ≈ p j ( v n ) + d p j ( v j ) d v j     v j = v n × ( v j − v n ) (41) q j ( v j ) ≈ q j ( v n ) + d q j ( v j ) d v j     v j = v n × ( v j − v n ) (42) The activ e/reactiv e po wer injections in (41) and (42) consist of two terms: one is the voltage-independent term, and the second is caused by non-zero sensitivity to nodal voltage. Our model can conv eniently include the first term as outlined previously . The second term can also be integrated in the model if the operator has a rough estimation of acti ve/reactiv e power voltage-sensiti vity values. For example, this sensitivity can be obtained for ZIP loads [28] and inv erters that are capable of reacti ve power support [18], [19] as follows: d p j ( v j ) d v j     v j = v n = p j ( v n ) · ( B j + 2 C j 2 v n ) (43) d q j ( v j ) d v j     v j = v n = k j (44) where, B j and C j represent the ZIP coefficients corresponding to the fixed-current and fixed-impedance portions of ZIP load, respectiv ely , and k j < 0 is the local in verter droop coefficient. Giv en the voltage-sensiti vity values, the second terms in (41) and (42) simply serve as new additional nodal activ e/reactive power injections and can be treated in the model similar to other loads. For example, the surrogate nodal acti ve/reacti ve SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 7 injections for ZIP loads and inv erters with reactive support capability can be conservati vely estimated as follows: ∆ p j ≈ ( B j + 2 C j 2 v n )( ¯ v − v n ) p j ( v n ) (45) ∆ q j ≈ k j ( ¯ v − v n ) (46) where, ¯ v denotes a conservati ve user-defined value that can be used by the utilities to model worst-case tripping scenar- ios. Howe ver , note that (45) and (46) are still conservati ve estimations. Dev eloping more accurate models for integrating voltage-dependent power injection into tripping equations re- mains the subject of future research. I I I . S O L A R C U RTA I L M EN T Q U A N T I FI C A T I O N A N D M I T I G A T I O N Using (35) as a conserv ative probabilistic lower bound for the real system, an optimization problem is formulated to provide a realistic estimation of the actual values of the micro-states of the grid. This problem is solved at any given time-window at which av ailable nodal activ e/reactive po wer statistics are kno wn: min ˆ λ ˆ λ ˆ λ − ( P P P > · ˆ λ ˆ λ ˆ λ ) , s.t. ˆ λ ˆ λ ˆ λ = a 0 a 0 a 0 + B ˆ λ ˆ λ ˆ λ +    ˆ λ ˆ λ ˆ λ > C 1 ˆ λ ˆ λ ˆ λ . . . ˆ λ ˆ λ ˆ λ > C N ˆ λ ˆ λ ˆ λ    0 ≤ ˆ λ j ≤ 1 ∀ j ∈ { 1 , ..., N } (47) where, P P P = [ P 1 , ..., P N ] > . The objective of this optimization problem is to find the maximum achiev able expected solar power in the gird according to the conservati ve statistical model. While the solution to this problem is still a lower bound estimation of the real achiev able PV po wer, the esti- mation gap between ˆ λ ˆ λ ˆ λ and λ λ λ is minimized. In other words, the optimization searches for the most optimistic v alues for micro-states with respect to the conservati ve approximate probabilistic tripping model. The problem is constrained by the matrix equations that govern the probabilities of inv erter tripping. Furthermore, the physical characteristics of micro- states are constrained by v alid probability assignments within [0 , 1] interval. A similar problem can be formulated to provide counter- measures against massive tripping e vents at any giv en time window . In general, the proposed statistical tripping model can be integrated as a constraint into any volt-v ar optimiza- tion formulation [29]–[31] to represent the possibility of PV curtailment. For example, here we provide a formulation for minimizing solar curtailment by controlling the voltage magnitude at the system reference bus [29]: min ˆ λ ˆ λ ˆ λ,v 0 − ( P P P > · ˆ λ ˆ λ ˆ λ ) , s.t. ˆ λ ˆ λ ˆ λ = a 0 a 0 a 0 ( v 0 ) + B ( v 0 ) ˆ λ ˆ λ ˆ λ +    ˆ λ ˆ λ ˆ λ > C 1 ˆ λ ˆ λ ˆ λ . . . ˆ λ ˆ λ ˆ λ > C N ˆ λ ˆ λ ˆ λ    0 ≤ ˜ λ j ≤ 1 ∀ j ∈ { 1 , ..., N } v min ≤ v 0 ≤ v max v R min ≤ v 0 − v I 0 ≤ v R max v min ≤ µ v i ( ˆ λ ˆ λ ˆ λ, v 0 ) ≤ v max ∀ i ∈ { 1 , ..., N } (48) where, v 0 is integrated into the optimization problem as a decision v ariable. Constraints are added to ensure that the control action and the expected nodal voltage magnitudes remain within permissible boundaries [ v min , v max ] . Here, v I 0 represents the initial setpoint value for v 0 , and [ v R min , v R max ] is the permissible range of rate of change of voltage at the reference b us with respect to the initial voltage setpoint. T o integrate v 0 into the problem, the expected nodal voltage magnitude squared values are written as a function of network parameters, expected av ailable nodal activ e/reactiv e po wers, and the optimization decision v ariables using (18):    µ v 1 . . . µ v N    ≈    R 11 P 1 + X 11 Q 1 . . . R 1 N P N + X 1 N Q N . . . . . . . . . R N 1 P 1 + X N 1 Q 1 . . . R N N P N + X N N Q N    ˆ λ ˆ λ ˆ λ + v 0 v 0 v 0 (49) where, v 0 v 0 v 0 = [ v 0 , ..., v 0 ] > . Despite its conv enient differentiable matrix-form formula- tion, the probabilistic tripping model introduces quadratic non- con vex constraints into optimization problems. This challenge can be addressed using various relaxation techniques from the literature, such as semidefinite program (SDP) relaxation [32], second-order cone program (SOCP) relaxation [33], and parabolic relaxation [34]. T o handle the non-con vexity , these methods generally define an auxiliary matrix, Λ = ˆ λ ˆ λ ˆ λ ˆ λ ˆ λ ˆ λ > , which enables obtaining a con vex surrogate for the original problem. For example, by applying parabolic relaxation, the constraints defined by the model are replaced with the follo wing alterna- tiv e constraints: a 0 a 0 a 0 + ( B − I N ) ˆ λ ˆ λ ˆ λ +    C 1 • Λ . . . C N • Λ    −    + ≤ 0 0 0 (50) − a 0 a 0 a 0 − ( B − I N ) ˆ λ ˆ λ ˆ λ −    C 1 • Λ . . . C N • Λ    +    − ≤ 0 0 0 (51) ∀ i, j : ( Λ( i, i ) + Λ( j, j ) − 2Λ( i, j ) ≥ ( ˆ λ ( i ) − ˆ λ ( j )) 2 Λ( i, i ) + Λ( j, j ) + 2Λ( i, j ) ≥ ( ˆ λ ( i ) + ˆ λ ( j )) 2 (52) SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 8 UG 1 2 3 0 7 UG 1 1 9 9 7 UG 2 81 UG 2 85 OH 1 3 1 1 OH 1 1 3 6 OH 1 3 7 0 OH 5 1 6 1 OH 5 3 4 7 OH 5 2 6 5 OH 5 3 1 2 OH 5 1 7 0 OH 5 2 9 5 OH 1 3 2 2 OH 5 1 8 1 OH 5 2 8 1 OH 5 3 0 7 OH 5 3 0 3 OH 5 2 1 0 OH 5 2 5 9 OH 1 2 8 9 OH 5 1 6 4 OH 5 3 6 8 OH 5 2 4 1 OH 5 3 1 7 OH 5 1 8 1 OH 5 3 2 0 OH 5 3 3 1 OH 1 3 0 4 OH 1 1 0 3 OH 1 2 7 1 OH 1 1 3 6 OH 1 2 5 2 OH 1 2 3 5 OH 1 2 6 3 UG 2 2 8 0 UG 2 1 3 3 5 OH 1 2 4 3 OH 1 2 5 1 OH 1 3 8 8 OH 1 2 5 7 OH 1 3 3 2 OH 1 3 2 9 OH 5 1 4 4 OH 5 3 9 3 OH 1 2 1 9 OH 1 1 3 0 OH 1 2 2 6 OH 1 1 6 7 OH 1 3 6 0 OH 1 2 9 2 UG 2 1 2 8 UG 2 1 7 8 UG 2 1 7 8 UG 2 2 7 4 UG 2 77 UG 2 1 8 4 UG 2 4 4 3 UG 2 1 9 8 UG 2 1 5 7 UG 2 2 0 5 UG 2 1 2 3 UG 2 2 7 4 OH 5 1 2 1 OH 5 3 2 8 OH 5 4 1 8 OH 1 1 6 4 OH 1 3 2 9 UG 2 1 0 7 UG 2 1 4 1 UG 2 90 OH 1 20 OH 1 3 4 0 OH 1 8 7 7 OH 1 80 OH 1 1 9 1 OH 1 6 0 4 OH 1 1 9 5 OH 1 5 4 0 UG 2 1 2 5 UG 2 1 8 5 UG 1 1 4 6 1 UG 2 1 9 6 UG 2 5 1 5 UG 2 3 3 8 UG 2 2 3 9 UG 2 4 2 0 UG 2 1 1 2 UG 2 5 6 5 UG 2 2 2 2 UG 2 2 3 0 UG 2 1 2 5 UG 2 65 UG 2 2 9 6 UG 2 1 1 7 UG 2 2 8 4 UG 2 3 6 8 UG 2 1 8 1 UG 2 1 2 0 UG 2 2 2 3 UG 2 1 5 7 UG 2 2 2 5 UG 2 1 3 4 UG 2 4 6 6 UG 2 1 3 6 UG 2 1 5 2 UG 2 2 76 UG 2 3 2 8 UG 2 1 1 2 UG 2 4 5 6 UG 2 2 4 5 UG 2 1 4 9 UG 2 3 3 1 UG 2 2 5 5 UG 1 2 2 6 7 UG 2 4 6 4 UG 2 5 0 1 UG 2 1 7 7 UG 2 2 0 6 UG 2 3 5 0 UG 2 91 UG 2 67 UG 2 1 5 1 UG 2 3 4 3 UG 2 2 3 4 UG 2 2 3 8 UG 2 2 4 3 UG 2 2 3 0 UG 2 1 9 4 UG 2 1 4 5 UG 2 2 8 1 UG 2 2 4 0 UG 2 3 6 7 UG 2 2 5 2 UG 2 3 5 1 UG 2 2 8 9 UG 2 1 3 4 UG 2 1 9 3 UG 2 2 0 2 UG 2 1 4 7 UG 2 2 5 3 UG 2 1 3 9 UG 2 1 5 6 UG 2 2 6 0 UG 2 1 3 9 UG 2 1 1 9 UG 2 5 0 5 UG 2 3 1 9 UG 2 1 2 2 UG 2 1 7 8 UG 2 1 7 8 UG 2 1 7 8 UG 2 1 5 3 UG 2 2 9 5 UG 2 2 5 7 UG 2 2 4 3 UG 2 2 3 7 UG 2 1 7 6 1 3 0 0 2 3 0 0 4 3 0 0 6 3 0 0 5 3 0 0 9 30 10 30 11 30 12 30 14 30 13 30 16 30 17 30 18 30 19 30 20 30 21 30 23 30 24 30 25 30 26 30 27 30 28 30 29 30 31 30 32 30 33 30 34 30 35 30 36 30 37 30 38 30 41 30 42 30 43 30 44 30 45 30 47 30 48 30 49 30 50 30 5 1 30 5 2 30 5 3 30 5 4 30 5 6 3 0 5 7 3 0 5 8 3 0 5 9 3 0 6 0 3 0 6 1 3 0 6 2 3 0 6 3 3 0 6 4 3 0 6 5 3 0 6 6 3 0 6 7 3 0 7 0 3 0 7 1 3 0 7 2 3 0 7 3 3 0 7 4 3 0 7 7 3 0 7 8 3 0 7 9 3 0 8 1 3 0 8 3 3 0 8 4 3 0 8 5 3 0 8 6 3 0 8 7 3 0 8 8 3 0 8 9 3 0 9 0 3 0 9 1 UG 2 91 3 0 9 3 3 0 9 4 3 0 9 5 3 0 9 6 3 0 9 7 3 0 9 8 3 0 9 9 3 1 0 1 3 1 0 2 3 1 0 3 3 1 0 4 3 1 0 5 3 1 0 6 3 1 0 8 3 1 0 9 3 1 1 0 3 1 1 1 3 1 1 2 3 1 1 4 3 1 1 5 3 1 1 6 3 11 7 3 1 2 0 3 1 2 1 3 1 2 2 3 1 2 3 3 1 2 4 3 1 2 5 3 1 2 6 3 1 2 7 3 1 2 8 3 1 2 9 3 1 3 0 3 1 3 1 3 1 3 2 3 13 4 3 13 5 3 13 6 3 13 7 3 1 3 8 3 1 3 9 3 1 4 1 3 1 4 2 3 1 4 3 3 1 4 4 3 1 4 5 3 1 4 6 3 1 4 7 3 1 4 8 3 1 4 9 3 1 5 0 3 1 5 1 3 1 5 2 3 1 5 3 3 1 5 4 3 1 5 7 3 1 5 8 3 1 5 9 3 1 6 0 3 16 1 3 16 2 3 0 0 3 3 0 0 7 3 0 0 8 30 15 30 22 30 30 30 39 30 40 30 46 3 0 5 5 3 0 6 8 3 0 6 9 3 0 7 6 3 0 8 0 3 0 8 2 3 0 9 2 3 1 0 0 3 1 0 7 3 1 1 3 3 1 1 8 3 1 1 9 3 13 3 3 1 4 0 3 1 5 6 3 1 5 5 OH 1 1 8 8 B B B B B B B B B B B B C C C C C C C C C C C C C C C C C C C C C A A A A B B B B C C C C A A A A A A A B B B B B B B B B B B B C C C C C A A A A A C C C C C C A A A A A A A A A A A A A A B B B B A A A A A A A CB _ 302 R e g 1 0 0 2 1 0 0 3 UG 1 2 9 6 7 ft UG 2 3 7 2 ft 1 0 0 4 OH 1 6 3 8 ft 1 0 0 5 OH 1 3 9 4 ft OH 1 1 0 4 9 ft 1 0 0 6 1 0 0 7 OH 2 2 0 0 0 ft OH 1 4 5 4 ft 1 0 0 8 1 0 0 9 10 10 10 11 10 12 10 13 10 14 10 15 10 16 10 17 OH 1 1 0 8 2 ft OH 1 1 6 9 ft UG 2 6 0 ft UG 2 6 0 ft UG 2 3 0 6 ft UG 2 2 0 8 ft UG 2 1 2 6 ft UG 2 4 6 8 ft UG 2 7 1 ft AB B B C C Eq . S o ur c e UG 1 1 3 8 8 ft OH 1 1 2 8 1 ft OH 1 7 9 8 ft OH 1 3 1 7 ft OH 1 1 4 5 ft OH 5 2 4 5 ft OH 1 1 7 0 ft OH 1 2 7 5 ft OH 1 5 0 3 ft OH 1 2 6 8 ft UG 2 1 2 0 ft OH 1 2 5 7 ft UG 2 1 7 4 ft UG 2 4 4 7 ft OH 1 4 5 ft OH 1 4 3 ft UG 2 2 0 5 ft UG 2 4 0 8 ft UG 2 3 3 1 ft OH 1 3 6 9 ft UG 2 2 8 5 ft OH 1 1 4 9 ft OH 1 2 8 9 ft OH 1 1 6 0 ft OH 1 8 0 ft OH 1 1 9 2 ft UG 2 3 5 ft UG 2 1 5 8 ft UG 2 3 4 2 ft UG 2 3 2 3 ft UG 2 3 3 0 ft UG 2 1 8 9 ft UG 2 2 5 9 ft UG 2 1 3 0 ft UG 2 1 2 4 ft UG 2 1 3 9 ft UG 2 1 3 4 ft UG 2 4 6 4 ft UG 2 1 3 9 ft UG 2 1 5 6 ft UG 2 2 0 4 ft UG 2 4 0 7 ft UG 2 6 8 5 ft UG 2 3 0 4 ft 2 0 0 2 2 0 0 3 2 0 0 4 2 0 0 5 2 0 0 6 2 0 0 7 20 10 20 11 20 1 3 20 1 4 20 1 5 20 1 6 20 1 7 20 1 8 20 1 9 20 2 1 2 0 2 0 20 2 4 20 2 5 20 2 7 20 2 8 20 2 9 2 0 3 1 20 3 2 20 3 4 2 0 3 5 20 3 6 20 3 7 20 4 2 20 4 3 20 4 4 20 4 5 20 4 6 20 4 7 20 4 8 20 4 9 2 0 5 0 2 0 5 1 20 5 2 20 5 3 20 5 4 20 5 5 20 5 6 20 5 7 20 5 8 20 5 9 2 0 6 0 UG 2 7 6 1 ft OH 2 1 6 1 ft UG 2 1 2 5 ft UG 2 6 4 2 ft UG 2 1 5 0 ft UG 2 1 3 1 ft UG 2 1 6 5 ft UG 2 2 4 5 ft OH 5 1 8 5 ft OH 1 8 0 0 ft OH 1 3 0 9 ft OH 1 3 8 3 ft 2 0 0 8 20 2 2 20 2 3 2 0 3 0 20 3 3 20 3 8 20 3 9 2 0 4 0 2 0 4 1 2 0 0 9 OH 1 1 6 0 ft A A CB _ 202 CB _ 203 B B B B B B C C C C A A B B B B C C C A A A B B F e e de r A F e e de r B F e e de r C X f m r 20 12 20 2 6 3 0 7 5 C A P _ 201 C A P _ 301 A B B B CB _ 101 CB _ 201 CB _ 301 1 0 0 1 2 0 0 1 3 0 0 1 Fig. 2. Structure of the 240-node test system. Fig. 3. Nodal PV outputs in the test system. where, C i • Λ = P N n =1 P N m =1 { C i ( n, m )Λ( n, m ) } , I N is an N × N identity matrix, and    + /    + are positiv e/negati ve small-valued slack variables that are used for transforming equality constraints defined by the model into two equi valent inequality constraints. The obtained inequalities (50)-(52) are con vex constraints with respect to variables ˆ λ ˆ λ ˆ λ and Λ . I V . N U M E R I C A L E X P E R I M E N T S A N D V A L I D A T I O N Numerical experiments hav e been performed to validate the proposed probabilistic tripping model. In this, we ha ve used real feeder model of an Iowa distribution system from our utility partner as shown in Fig. 2. The network model in OpenDSS and detailed parameters are av ailable online [35]. T o perform simulations we have used real solar and load data with 1-second time resolution from [36]. Fig. 3 shows the PV outputs at different nodes in the system for one day . Fig. 4 demonstrates 15-minute average nodal demand. The load/PV data have been randomly distributed across the three phases of the unbalanced grid at each node. T o verify the performance of the proposed approximate statistical model, extensiv e time-series simulations were per - formed on the test system under various loading and solar generation scenarios o ver a course of day . Then, the real v alues of original micro-states, λ i , were determined empirically ov er time windows of length T = 60 minutes. Intuitiv ely , λ i serves as the gr ound truth and roughly represents the portion of time that s i is ON during each time window: λ i ( T ) ≈ P T t =1 s i ( t ) T (53) Fig. 4. A verage 15-minute nodal consumption in the test system. Thus, we have two distinct time windows throughout nu- merical studies: a 1-second time step is used to perform high-resolution simulations, and a 1-hour time window is employed to obtain tripping statistics and empirically verify the performance of the proposed probabilistic model. Fig. 5a demonstrates the empirical micro-states, λ i , at dif ferent time intervals, which are determined by applying (53) to simulation outcomes. Based on the values of these micro- states, the empirical macro-state v alue is calculated at all time intervals, which represents the expected percentage of PV switches in ON state, i.e., S p ( T ) ≈ P N i =1 λ i ( T ) N × 100 . Fig. 5b compares the empirical macro-state value and the lower bound value constructed using solutions of (47). As can be seen, the solution from the probabilistic model actually represents a lower bound to the empirical macro-state obtained from simulations at all time windows, which corroborates the per- formance of the method. This figure also shows another lower bound obtained by simply using maximum PV capacities and assuming zero nodal consumption. Howe ver , as can be seen, this lo wer bound gives fixed ov er-conserv ative outcomes that do not reflect the true conditions of the system and hav e no correlation with the time-series PV/load data. Fig. 5c depicts the aggregate maximum av ailable solar power (all switches ON at all time), empirical aggregate realized solar power from numerical simulations (53), and solar power corresponding to solution of (47). As observed, the lower bound solution still holds and provides a conservati ve yet close estimation for the empirical achie vable solar po wer outcome. Fig. 6 compares the empirical and model-based probabilities of in verter tripping in a heavy-loaded time interval. Unlike the previous case, these tripping probabilities are due to under-v oltages. As can be observed, the model still provides a conservati ve lower bound on the probability of tripping. Note that the reason for higher lev els of volatility in this figure is the shorter time window (15 minutes) used for assessing the empirical probability of tripping. The gap between the empirical macro-state obtained from numerical experiments and the proposed lo wer bound is an implicit function of PV penetration. Sensitivity analysis was performed to quantify the relationship between this gap and PV penetration percentage, as shown in Fig. 7. Here, PV penetration is defined as the mean value of peak nodal solar power o ver peak nodal demand. The maximum, minimum, and mean values of the gap between the provided lower bound and SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 9 (a) λ λ λ (b) S p vs ˆ S p (c) Aggregate expected PV power Fig. 5. Comparing the empirical and statistical lower bound solutions. Fig. 6. Model performance for a case of heavy-loaded system and 15-minute empirical tripping probability assessment time window . Fig. 7. Lower bound gap as a function of PV integration. Fig. 8. Solar curtailment sensitivity to inv erter control setpoints. the empirical macro-state is measured at various le vels of PV penetration. As is observed in the figure, the optimistic value of the gap drops and ev entually reaches 5% as PV penetration increases, which indicates that the lower bound approaches the true macro-state v alue in grids with higher PV penetration. On the other hand, the maximum value of the gap shows an increase after a certain PV penetration lev el which points out to higher v ariations in solutions obtained from the probabilistic model. Fig. 8 sho ws the overall daily solar curtailment levels, both empirical and the lower bound, as a function of changes in inv erter control parameter . The inv erters in the system are assumed to be controlled in constant power factor (PF) mode. As the reference PF setpoint increases and the system mov es towards unity PF the voltage fluctuations increase, which leads to higher solar curtailment. This confirms pre vious observations in the literature [9]. Furthermore, our proposed probabilistic lower bound always slightly over -estimates the curtailment lev el, as expected correctly from the conserv ative estimator . Further tests were performed to corroborate the performance of countermeasure design strategy introduced in (48). Fig. 9a shows the outcome of the optimization problem (48), compared to a base case without any voltage regulation. As observed, v 0 is optimally decreased during solar-rich intervals to compensate for the increased voltage fluctuation lev els. Fig. 9b compares the aggregate solar po wer injection values under the ne wly acquired v 0 values and the base case without voltage control. As can be seen, the obtained countermeasure has assisted significantly in mitigating the ov erall solar power curtailments during critical time interv als. SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 10 (a) Original and regulated v 0 (b) Macro-state under original and regulated v 0 Fig. 9. Solar curtailment countermeasure design verification W e ha ve performed another numerical experiment to analyse and verify the behavior of our tripping model during an under- voltage case study in a temporary heavy loading scenario in a weak grid under two strategies (see Fig. 10): (1) No voltage regulation is applied (baseline), and (2) V oltage regu- lation is applied with the objecti ve of minimizing the av erage squared voltage deviations across the whole system, subject to linearized power flow equations and the proposed statistical tripping model. As can be seen in Fig. 10a, under the baseline strategy (no voltage regulation) a portion of in verters (around 13%) have tripped due to under-v oltage protection during later hours of the day . This has resulted in a loss of rene wable po wer injection in the grid (Fig. 10b). Ho wever , by applying voltage regulation using the proposed tripping model we hav e been able to maintain the voltages much closer to their nominal values (see Fig. 10c) and prev ent tripping events and loss of solar generation resources altogether . Note that Fig. 10c shows the av erage value of nodal voltages across the whole system; thus, while most of the nodes maintain healthy voltage lev els (as they should), the excessi ve loading on weak system lines under the baseline has resulted to a temporary v oltage drop below in verters’ protection activ ation threshold, which has engaged their under-v oltage protection de vices. This issue was mitigated using the deployed voltage regulation strategy that leverages our proposed statistical tripping model. Fig. 11 demonstrates the average realized daily PV power ratio as a function of average PV penetration. As can be seen, the increasing penetration of solar has led to a regime shift after a certain threshold, from an initial state, in which the system shows almost no extensi ve tripping, to a new state, in which the av erage probability of solar curtailment steadily (a) Percentage of functional inv erters (b) Aggregate nodal solar power (c) A verage nodal voltage magnitude across the grid Fig. 10. An under-v oltage case study . increases and extended tripping ev ents can be expected. The existence of this threshold attests to a stage transition in the extent of switching ev ents, which has been observed in other nonlinear systems as well [26]. Above the PV integration threshold, which is around 30% for the test system, massive solar curtailment can be expected due to voltage fluctuations. It can be observ ed that the proposed statistical lower bound accurately tracks the behavior of the real system, and can be used to con vey information on the whereabouts of the transition. The exact value of the regime shift threshold depends on many factors, including network topology and spatial-temporal distribution of loads/generators. V . C O N C L U S I O N S In this paper , a probabilistic model of interdependent solar in verter tripping is presented to assess the risk of solar po wer curtailments due to v oltage fluctuations in distribution grids. This model is dev eloped using only the statistical properties of av ailable load/PV acti ve/reacti ve po wer . Numerical results on a SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 11 Fig. 11. Regime shift (stage transition) analysis. real distribution feeder using real data successfully validate the estimated conservati ve lo wer bounds on in verter micro-states. Furthermore, it is demonstrated that the proposed model can be used for identifying regime shifts in tripping ev ents and designing countermeasures to minimize risk of solar power curtailment. As a future research direction, we will explore integrating the more dynamic functions of in verter control and protection, including ride-through capabilities, [17]–[20] into the probabilistic tripping model. For example, the pro- posed statistical lower bound, which is based on Chebyshev’ s inequality , might become too conservati ve over short time windows if inv erters’ disturbance ride-through capabilities are activ ated. A less conserv ative lower bound that incorporates all aspects of inv erter behavior will enable operators to monitor the sequence and transitions of tripping ev ents, and mitigate potential cascading failure of resources. R E F E R E N C E S [1] R. T onkoski, D. T urcotte, and T . H. M. El-Fouly , “Impact of high PV penetration on voltage profiles in residential neighborhoods, ” IEEE T rans. Sustain. Energy , vol. 3, no. 3, pp. 518–527, Jul. 2012. [2] O. Gagrica, P . H. Nguyen, W . L. Kling, and T . Uhl, “Microinv erter curtailment strategy for increasing photovoltaic penetration in low- voltage networks, ” IEEE T rans. Sustain. Energy , vol. 6, no. 2, pp. 369– 379, Apr . 2015. [3] M. Hasheminamin, V . G. Agelidis, V . Salehi, R. T eodorescu, and B. Hredzak, “Index-based assessment of voltage rise and reverse power flow phenomena in a distribution feeder under high PV penetration, ” IEEE J. Photovolt. , vol. 5, no. 4, pp. 1158–1168, Jul. 2015. [4] V . Ferreira, P . M. S. Carvalho, L. A. F . M. Ferreira, and M. D. Ilic, “Distributed energy resources integration challenges in low-voltage networks: v oltage control limitations and risk of cascading, ” IEEE Tr ans. Sustain. Energy , vol. 4, no. 1, pp. 82–88, Jan. 2013. [5] P . M. S. Carvalho, L. A. Ferreira, J. C. Botas, M. D. Ilic, X. Miao, and K. D. Bachovchin, “Ultimate limits to the fully decentralized power in verter control in distribution grids, ” In P ower Systems Computation Confer ence (PSCC) , pp. 1–7, Jun. 2016. [6] P . M. S. Carvalho, L. A. Ferreira, and J. J. E. Santana, “Single-phase generation headroom in low-voltage distrib ution networks under reduced circuit characterization, ” IEEE Tr ans. P ower Syst. , vol. 30, no. 2, pp. 1006–1011, Mar . 2015. [7] C. Popiel and P . D. H. Hines, “Understanding factors that influence the risk of a cascade of outages due to inverter disconnection, ” In North American P ower Symposium (NAPS) , pp. 1–6, 2019. [8] A. Parchure, S. J. T yler , M. A. Peskin, K. Rahimi, R. P . Broadwater , and M. Dilek, “In vestigating PV generation induced voltage volatility for customers sharing a distribution service transformer , ” IEEE T rans. Ind. Appl. , vol. 53, no. 1, pp. 71–79, Jan. 2017. [9] D. Cheng, B. A. Mather , R. Seguin, J. Hambrick, and R. P . Broadwater, “Photov oltaic (PV) impact assessment for very high penetration levels, ” IEEE J. Photovolt. , vol. 6, no. 1, pp. 295–300, Jan. 2016. [10] P . Gupta, K. Rahimi, R. Broadwater, and M. Dilek, “Importance of detailed modeling of loads/PV systems connected to secondary of dis- tribution transformers, ” In North American P ower Symposium (NAPS) , pp. 1–6, Sep. 2017. [11] F . Ding and B. Mather, “On distributed PV hosting capacity estimation, sensitivity study , and impro vement, ” IEEE Tr ans. Sustain. Energy , v ol. 8, no. 3, pp. 1010–1020, Jul. 2017. [12] X. Xu, Z. Y an, M. Shahidehpour, H. W ang, and S. Chen, “Power system voltage stability e valuation considering rene wable energy with correlated variabilities, ” IEEE Tr ans. P ower Syst. , vol. 33, no. 3, pp. 3236–3245, May 2018. [13] H. Jiang, Y . Zhang, J. J. Zhang, D. W . Gao, and E. Muljadi, “Synchrophasor-based auxiliary controller to enhance the voltage sta- bility of a distribution system with high renewable energy penetration, ” IEEE T rans. Smart Grid , vol. 6, no. 4, pp. 2107–2115, Jul. 2015. [14] K. I. Jennett, C. D. Booth, F . Coffele, and A. J. Roscoe, “Investig ation of the sympathetic tripping problem in power systems with large penetrations of distributed generation, ” IET Gener . T ransm. Distrib. , vol. 9, no. 4, pp. 379–385, Dec. 2014. [15] S. M. Ross, Intr oduction to Probability Models . Academic Press, 2006. [16] “WECC PV power plant dynamic modeling guide, ” W estern Electricity Coordinating Council, T ech. Rep., Apr . 2014. [Online]. A vailable: https://www .wecc.org/Reliability/WECC [17] “NERC reliability guideline: Integrating in verter-based resources into weak power systems - Dec. 2017. ” [Online]. A v ailable: https://www .nerc.com/comm/PC Reliability Guidelines DL/Item 4a. Integrating%20 Inverter - Based Resources into Low Short Circuit Strength Systems - 2017- 11- 08- FIN AL.pdf [18] “Australian Energy Market Operator (AEMO) Information & Support Hub: Multiple voltage disturbance ride-through capability , justification of AEMO’ s proposal - Mar . 2018. ” [Online]. A v ailable: https://www . aemc.gov .au/sites/default/files/2018- 03/AEMO%20report%20updated% 20proposed%20multiple%20fault%20withstand%20obligation.pdf [19] “Australian Ener gy Market Operator (AEMO) Informa- tion & Support Hub: Electricity rule change proposal, generator technical requirement - Aug. 2017. ” [Online]. A v ailable: https://www .aemo.com.au/- /media/Files/Electricity/NEM/ Security and Reliability/Reports/2017/AEMO- GTR- RCP- 110817.pdf [20] M. G. Dozein, P . Mancarella, T . K. Saha, and R. Y an, “System strength and weak grids: Fundamentals, challenges, and mitigation strategies, ” In Austr alasian Universities P ower Engineering Conference (A UPEC) , pp. 1–7, 2018. [21] N. Li, G. Qu, and M. Dahleh, “Real-time decentralized voltage control in distribution networks, ” In 52nd Annual IEEE Conference on Commu- nication, Contr ol, and Computing (Allerton) , pp. 582–588, Oct. 2014. [22] R. Singh, B. C. Pal, and R. A. Jabr, “Statistical representation of distribution system loads using Gaussian mixture model, ” IEEE T rans. P ower Syst. , vol. 25, no. 1, pp. 29–37, Feb . 2010. [23] D. T . Nguyen, “Modeling load uncertainty in distribution network monitoring, ” IEEE T rans. P ower Syst. , vol. 30, no. 5, pp. 2321–2328, Sep. 2015. [24] L. Hatton, P . Charpentier, and E. Matzner-Lober , “Statistical estimation of the residential baseline, ” IEEE Tr ans. P ower Syst. , vol. 31, no. 3, pp. 1752–1759, May 2016. [25] H. L. van der W alt, R. C. Bansal, and R. Naidoo, “Pv based distributed generation power system protection: A review , ” Renewable Energy F ocus , vol. 24, pp. 33–40, Mar . 2018. [26] S. Thurner , R. Hanel, and P . Klimek, Intr oduction to the theory of complex systems . Oxford University Press, 2018. [27] M. Mezard and J. Sakellariou, “Exact mean-field inference in asymmet- ric kinetic Ising systems, ” Journal of Statistical Mechanics: Theory and Experiment , no. 07, p. L07001, Jul. 2011. [28] A. Arif, Z. W ang, J. W ang, B. Mather, H. Bashualdo, and D. Zhao, “Load modeling - a revie w , ” IEEE T rans. Smart Grid , vol. 9, no. 6, pp. 5986–5999, Nov . 2018. [29] J. O. Petinrin and M. Shaaban, “Impact of renewable generation on voltage control in distribution systems, ” Renewable and Sustainable Ener gy Revie ws , vol. 65, pp. 770–783, 2016. [30] D. Ranamuka, A. P . Agalgaonkar, and K. M. Muttaqi, “Online voltage control in distribution systems with multiple voltage regulating devices, ” IEEE T rans. Sustain. Energy , vol. 5, no. 2, pp. 617–628, Apr . 2014. [31] T . T . Ku, C. H. Lin, C. S. Chen, and C. T . Hsu, “Coordination of transformer on-load tap changer and PV smart in verters for voltage control of distribution feeders, ” IEEE T rans. Ind. Appl. , vol. 55, no. 1, pp. 256–264, Feb. 2019. SUBMITTED TO IEEE FOR POSSIBLE PUBLICA TION. COPYRIGHT MA Y BE TRANSFERRED WITHOUT NOTICE 12 [32] W . K. K. Ma, “Semidefinite relaxation of quadratic optimization prob- lems and applications, ” IEEE Signal Pr ocessing Magezine , vol. 1053, no. 5888/10, May 2010. [33] S. Kim and M. Kojima, “Second order cone programming relaxation of noncon vex quadratic optimization problems, ” Optimization methods and softwar e , vol. 15, no. 3–4, pp. 201–224, Jan. 2001. [34] M. Kheirandishfard, F . Zohrizadeh, and R. Madani, “Conve x relaxation of bilinear matrix inequalities part i: theoretical results, ” In IEEE Confer ence on Decision and Contr ol (CDC) , pp. 67–74, Dec. 2018. [35] F . Bu, Y . Y uan, Z. W ang, K. Dehghanpour, and A. Kimber . (2019) A time-series distribution test system based on real utility data. [Online]. A v ailable: https://arxiv .org/abs/1906.04078. [36] C. Holcomb, “Pecan Street Inc.: A test-bed for NILM, ” In International W orkshop on Non-Intrusive Load Monitoring , 2012. [Online]. A vailable: https://www .pecanstreet.org/