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

TL;DR

This paper introduces a probabilistic modeling approach using Chebyshev's inequality to assess and mitigate the risk of interdependent inverter tripping events in distribution grids with high PV penetration, addressing real-time observability limitations.

Contribution

It develops a scalable probabilistic model based on PV/load statistics for predicting inverter tripping risks and integrates it into an optimization framework for mitigation strategies.

Findings

Model accurately estimates tripping probabilities using only statistical data.

The approach is scalable and suitable for systems with limited real-time observability.

Numerical experiments validate the effectiveness of the proposed method.

Abstract

It is speculated that higher penetration of inverter-based distributed photo-voltaic (PV) power generators can increase the risk of tripping events due to voltage fluctuations. To quantify this risk utilities need to solve the interactive equations of tripping events for networked PVs in real-time. However, these equations are non-differentiable, nonlinear, and exponentially complex, and thus, cannot be used as a tractable basis for solar curtailment prediction and mitigation. Furthermore, load/PV power values might not be available in real-time due to limited grid observability, which further complicates tripping event prediction. To address these challenges, we have employed Chebyshev's inequality to obtain an alternative probabilistic model for quantifying the risk of tripping for networked PVs. The proposed model enables operators to estimate the probability of interdependent inverter tripping events using only PV/load statistics and in a scalable manner. Furthermore, by integrating this probabilistic model into an optimization framework, countermeasures are designed to mitigate massive interdependent tripping events. Since the proposed model is parameterized using only the statistical characteristics of nodal active/reactive powers, it is especially beneficial in practical systems, which have limited real-time observability. Numerical experiments have been performed employing real data and feeder models to verify the performance of the proposed technique.

Equations

$$s_i(t)=\{\begin{array}{l}1V_{min}\le V_i(t)\le V_{max}\\ 0V_i(t)<V_{min}\\ 0V_i(t)>V_{max}\end{array}$$

$$s_i(t)=U(v_i(t)-v_{min})-U(v_i(t)-v_{max})$$

$$U(x)=\{\begin{array}{l}1x\ge 0\\ 0x<0,\end{array}$$

$$v_i(t)=\sum _{j=1}^N{\tilde{v}}_{ij}+v_0,\forall i\in \{1,...,N\}$$

$${\tilde{v}}_{ij}=R_{ij}{\tilde{p}}_j(t)+X_{ij}{\tilde{q}}_j(t)$$

$$R_{ij}=2\sum _{\{n,m\}\in Pa(i,j)}{r}_{nm}$$

$$X_{ij}=2\sum _{\{n,m\}\in Pa(i,j)}{x}_{nm}$$

$${\tilde{p}}_j(t)=p_j(t)s_j(t)$$

$${\tilde{q}}_j(t)=q_j(t)s_j(t)$$

$${\lambda }_i(t)=Pr\{v_{min}\le v_i(t)\le v_{max}\}$$

$$Pr\{v_{min}\le v_i(t)\le v_{max}\}\ge 1-\frac{{\sigma }_{v_i}^2+({\mu }_{v_i}-\frac{v_{max}+v_{min}}{2}{)}^2}{(\frac{v_{max}-v_{min}}{2}{)}^2}$$

$${\hat{\lambda }}_i(t)=1-\frac{{\sigma }_{v_i}^2+({\mu }_{v_i}-\frac{v_{max}+v_{min}}{2}{)}^2}{(\frac{v_{max}-v_{min}}{2}{)}^2}$$

$${\mu }_{v_i}=$$

$$=$$

$$F_{{\tilde{p}}_j}(P)=Pr\{{\tilde{p}}_j(t)\le P\}=(1-{\lambda }_j(t))U(P)+{\lambda }_j(t)F_{p_j}(P)$$

$$f_{{\tilde{p}}_j}(P)=\text{derivative}{F}_{{\tilde{p}}_j}(P)P=(1-{\lambda }_j(t))\delta (P)+{\lambda }_j(t)f_{p_j}(P)$$

$$E\{{\tilde{p}}_j\}={\int }_{-\infty }^{+\infty }\alpha f_{{\tilde{p}}_j}(\alpha )\mathrm{d}\alpha ={\lambda }_j{P}_j$$

$$E\{{\tilde{q}}_j\}={\int }_{-\infty }^{+\infty }\beta f_{{\tilde{q}}_j}(\beta )\mathrm{d}\beta ={\lambda }_j{Q}_j$$

$${\mu }_{v_i}=\sum _{j=1}^N\{R_{ij}{\lambda }_j(t)P_j+X_{ij}{\lambda }_j(t)Q_j\}+v_0$$

$${\sigma }_{v_i}^2=\sum _{j=1}^N{\sigma }_{{\tilde{v}}_{ij}}^2+2\sum _{1\le k<j\le N}\Omega \{{\tilde{v}}_{ij},{\tilde{v}}_{ik}\}$$

$${\sigma }_{{\tilde{v}}_{ij}}^2=R_{ij}^2{\sigma }_{{\tilde{p}}_j}^2+X_{ij}^2{\sigma }_{{\tilde{q}}_j}^2+2R_{ij}{X}_{ij}\Omega \{{\tilde{p}}_j,{\tilde{q}}_j\}$$

$${\sigma }_{{\tilde{p}}_j}^2=E\{s_j^2{p}_j^2\}-E\{{\tilde{p}}_j{\}}^2$$

$$f_{s_j^2{p}_j^2}(\zeta )=(1-{\lambda }_j(t))\delta (\zeta )+\frac{{\lambda }_j(t)}{2\sqrt{\zeta }}(f_{p_j}(\sqrt{\zeta })+f_{p_j}(-\sqrt{\zeta }))$$

$${\sigma }_{{\tilde{p}}_j}^2={\lambda }_j(P_j^{+}+P_j^{-})-{\lambda }_j^2{P}_j^2$$

$${\sigma }_{{\tilde{q}}_j}^2={\lambda }_j(Q_j^{+}+Q_j^{-})-{\lambda }_j^2{Q}_j^2$$

$$\Omega \{{\tilde{p}}_j,{\tilde{q}}_j\}=E\{{\tilde{p}}_j{\tilde{q}}_j\}-E\{{\tilde{p}}_j\}E\{{\tilde{q}}_j\}$$

$$\Omega \{{\tilde{p}}_j,{\tilde{q}}_j\}={\lambda }_j{P}_j{Q}_j-{\lambda }_j^2{P}_j{Q}_j+{\lambda }_j\Omega \{p_j,q_j\}$$

$${\sigma }_{{\tilde{v}}_{ij}}^2={\lambda }_j{\Gamma }_{ij}^1-{\lambda }_j^2{\Gamma }_{ij}^2$$

$${\Gamma }_{ij}^1=$$

$$+2R_{ij}{X}_{ij}(P_j{Q}_j+\Omega \{P_j,Q_j\})$$

$${\Gamma }_{ij}^2=2R_{ij}{X}_{ij}{P}_j{Q}_j+P_j^2{R}_{ij}^2+Q_j^2{X}_{ij}^2$$

Full text

I Introduction

Increasing penetration of distributed energy resources (DERs), including inverter-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 inverter-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, reduce the economic value of renewable resources for customers, and cause loss of service in stand-alone systems [2, 3]. The possibility of DER power generation disruption due to voltage-related vulnerabilities in unbalanced distribution grids has been discussed in the literature: in [4, 5], risk of interdependent tripping of PVs, with ON/OFF current interruption mechanism was demonstrated numerically in a distribution grid test case for the first time. It was shown that the unbalanced and resistive nature of networks can further exacerbate this problem by causing positive inter-phase voltage sensitivity terms that act as destabilizing positive feedback loops, leading to voltage deviations after tripping of an individual inverter. The impact of grid voltage sensitivity 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 voltage deviations in the grid. It was shown in [7] that very large or small number of inverter-based resources in distribution systems can lead to interdependent failure events that contribute to voltage collapse in transmission level. Detailed realistic numerical studies were performed on practical feeder models in [8, 9, 10, 11, 12, 13] that corroborated the considerable impacts of extreme PV integration levels, and inverter control modes on grid voltage fluctuations, which is the critical factor in causing massive solar curtailment scenarios.

Most existing works relied on scenario-based simulations and numerical studies to capture the likelihood of inverter tripping under high renewable penetration. While this has led to useful guidelines and invaluable intuitions, it falls short of providing a generic theoretical foundation for predicting and containing tripping events. Specifically, the dependencies between nodal solar power distributions, nodal voltage profiles, and inverter tripping events have not been explicitly analyzed in the literature thus far. These dependencies are influenced by inverter protection settings and governed by a set of networked power-flow-based equations, which turn out to be non-differentiable and nonlinear. In this regard, several fundamental challenges have not been addressed: (1) Lack of scalability: Solving the inverter 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 inverters. The computational complexity is due to the interactive and networked nature of tripping events, meaning that the states of inverters influence each other and are not independent [14]. The source of interdependency in chances of inverter 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 inverters, which in turn could influence their probability of tripping.) For example, tripping of an inverter (or a cluster of inverters) leads to a change in loading distributions, which can most likely increase/decrease the chance of tripping for other inverters during under/over-voltage scenarios, especially in weak grids. This interdependency prevents the solver from decoupling the tripping equations into separate equations for individual inverters. Thus, the scale of search for finding the correct configuration increases exponentially ($2^{N}$) with the number of inverters ($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 integrated into optimization-based decision models, since it has no predictive capability and cannot be use to answer what if queries, unless a thorough expensive 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 low 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 communication time delays or limited number of sensors. Thus, we might not have access to sufficient online information to solve the tripping problem directly.

To tackle these challenges, we propose an alternative probabilistic modeling approach to quantify and mitigate the risk of voltage-driven tripping events. Instead of complex scenario-based look-ahead search over numerous possible tripping configurations, our methodology is built upon probabilistic manipulation of power flow equations in radial networks to estimate the probability of inverter tripping using only the available statistical properties of loads/PVs. Interdependent Bernoulli random variables are used to model probabilities of inverter tripping and capture their mutua. These probabilities are voltage-dependent and serve as unknown micro-states in the equations of tripping events. Then, Chebyshev’s inequality [15] is applied to determine a stationary lower bound for the values that these micro-states can assume under any probable nodal power injection scenarios. This lower bound provides a conservative estimation of expected PV curtailment, and thus, represents a statistical risk metric for tripping events. Furthermore, due to its simple matrix-form and differentiable structure, the proposed probabilistic model can be conveniently 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 events can be conservatively predicted using the proposed model and without the need for online access to granular PV/load data or expensive scenario-based search process, which makes our strategy specifically suitable for practical networks.

Numerical experiments have been performed using real advanced metering infrastructure (AMI) data and feeder models from our utility partners to validate the developed probabilistic framework. The numerical validates the performance of the probabilistic model for both over- and under-voltage scenarios, and show that ignoring the possibility of tripping in voltage regulation can exacerbate voltage deviations.

II Deriving A Conservative Probabilistic Model of PV Tripping Events

In this section, we will develop and then parameterize a probabilistic model of networked inverter-based PVs to quantify the possibility of emergent tripping events. To do this, first, we begin with the original model of inverter tripping with ON/OFF voltage-driven current interruption mechanism, and then, we will show that by adopting a probabilistic approach towards the original model and using Chebyshev’s inequality, tripping probabilities can be conservatively estimated using the statistical properties of nodal available load/PV power.

II-A Original Interactive Switching Equations

In this paper, it is assumed that PV resources are protected against voltage deviations using ON/OFF switching mechanisms. Note that here a “switch” can be a mechanical relay, as well as a non-physical inverter control function that stops current injection into the grid under abnormal voltage even if the inverter 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 micro-state variables with the following voltage-dependent function (see Fig. 1):

[TABLE]

where, $s_{i}(t)$ is the micro-state assigned to the $i$’th PV at time $t$ as a function of the inverter node’s voltage 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 inverter dynamics, e.g., ride-through capabilities [17, 18, 19, 20], can be conservatively ignored. This assumption considerably enhances the tractability of the model at the expense of loss of accuracy. In this sense, the switching model is a worst-case representation of inverter tripping. Since the approximate power flow equations for distribution grids are linear with respect to the squared values of nodal voltage magnitudes [21], we re-write equation (1) using a variable transformation, $v_{i}=V_{i}^{2}$, and employing unit step functions as follows:

[TABLE]

where, $v_{min}=V_{min}^{2}$, $v_{max}=V_{max}^{2}$, and the unit step function $U(\cdot)$ is defined as follows:

[TABLE]

Note that inverters’ 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 events. To obtain the overall governing equations of inverter tripping, the mutual impacts of switch micro-states on each other are captured using an approximate unbalanced power flow model for radial distribution grids [21], which determines voltage at node $i$ as a function of active/reactive power injections of every other node in a grid (with a total of $N+1$ nodes):

[TABLE]

where, $v_{0}=V_{0}^{2}$, with $V_{0}$ denoting the voltage magnitude at a grid reference bus, and the intermediary variable $\tilde{v}_{ij}$ represents the impact of active/reactive power injection at node $j$ on $v_{i}$, which is obtained as follows:

[TABLE]

where, $R_{ij}$ and $X_{ij}$ are the aggregated series resistance and reactance values corresponding to the intersecting branches in the paths connecting nodes $i$ and $j$ to the reference bus calculated as follows [21]:

[TABLE]

[TABLE]

where, $Pa(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 bus; $r_{nm}$ and $x_{nm}$ denote the real series resistance and reactance of the branch connecting nodes $n$ and $m$. Also, $\tilde{p}_{j}$ and $\tilde{q}_{j}$ denote the active and reactive power injections at bus $j$, which are in turn determined by the micro-state of the PV at node $j$ (see Fig. 1):

[TABLE]

[TABLE]

with $p_{j}$ and $q_{j}$ representing the available 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 individually, and a thorough search process is needed to explore all possible switching configurations. This can be extremely expensive and impossible to scale to large systems with high population of inverters. (II) Due to the discrete step functions in (2), tripping equations are nonlinear and non-differentiable. This contributes to problem difficulty 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 events requires extensive search process to cover all probable PV/load time-series scenarios. This expensive search process hinders the tractability of optimization-based frameworks for designing tripping mitigation strategy.

Not all the nodes in the tripping model are necessarily controlled by ON/OFF voltage-sensitive 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, however, we will simply assign constant values, $s_{i}(t)=1,\ \forall t$ to the nodes without ON/OFF control and remove their corresponding switching from the equations (see Fig. 1).

II-B Alternative Approximate Probabilistic Model

We 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 $\lambda_{i},\forall i\in\{1,...,N\}$: $s_{i}\sim\mathpzc{B}(\lambda_{i})$, where parameter $\lambda_{i}$ is defined as the probability of the $i$’th inverter switch being ON, $\lambda_{i}(t)=Pr\{s_{i}(t)=1\}$. The goal is to transform micro-states from discontinuous binary variables ($s_{i}\in\{0,1\}$) into continuous variables ($\lambda_{i}\in[0,1]$). To rewrite the equations in terms of new micro-states note that we have $E\{s_{i}(t)\}=\lambda_{i}(t)$ for Bernoulli probability distributions, where $E\{\cdot\}$ represents the expectation operation. Thus, by performing an expectation operation over both sides of (2), probability of inverter tripping in terms of the new micro-states can be obtained as follows:

[TABLE]

where, we have exploited $E\{U(f(x))\}=Pr\{f(x)\geq 0\}$. Note that the probability of tripping for an inverter is an implicit function of nodal voltage probability distribution, which in turn is influenced by the states of other inverters. Due to the interconnected nature of the problem, no independency assumptions has been made on random variables $\lambda_{i},\forall i\in\{1,...,N\}$. However, the exact distributions of nodal voltages are unknown and complex functions of nodal active/reactive injections, which implies that (10) cannot be determined analytically unless over-simplifying assumptions are made. Instead, we employ Chebyshev’s inequality [15] to provide a lower bound on micro-state as a function of nodal voltage statistics without making any assumption on voltage distributions,

[TABLE]

where, $\sigma_{v_{i}}^{2}$ and $\mu_{v_{i}}$ are the variance and mean of $v_{i}$, respectively. Hence, the approximate probabilistic model can be formulated for each micro-state as follows:

[TABLE]

This new tripping model has two features: (1) it is a conservative estimator of the original system since it over-estimates the probability of inverter tripping, $\hat{\lambda}_{i}\leq\lambda_{i}$. (2) As will be shown in Section II-C, the approximate probabilistic model can be conveniently parameterized in terms of nodal available active/reactive power statistics. Hence, as long as certain statistics are known (or estimated), the model allows us to accurately track probability of inverter tripping without running time-series simulations under numerous scenarios.

II-C Probabilistic Model Parameterization

To parameterize the alternative tripping model (12), nodal voltage statistics, $\sigma_{v_{i}}^{2}$ and $\mu_{v_{i}}$, are obtained in terms of nodal available active/reactive power statistics. To do this, power flow/injection equations (4)-(9) are leveraged.

Stage 1: $\mu_{v_{i}}$ Parameterization - The expected value of voltage magnitude squared is determined using (4)-(5) as,

[TABLE]

To calculate $E\{\tilde{p}_{j}\}$ and $E\{\tilde{q}_{j}\}$, we will first obtain their cumulative distribution functions (CDFs) [15], $F_{\tilde{p}_{j}}$ and $F_{\tilde{q}_{j}}$, respectively. This process is shown for $\tilde{p}_{j}$ as follows ($F_{\tilde{q}_{j}}$ is obtained similarly):

[TABLE]

The rational behind (14) is that the distribution of power injection is determined by two functions: the distribution of PV switch (which is ON with probability $\lambda_{j}(t)$), and the CDF of available PV power, $F_{p_{j}}$. Now, the probability density functions (PDF) of the realized active nodal power injection, $f_{\tilde{p}_{j}}$, can be calculated as a function of the available active solar power, $f_{p_{j}}$ (a similar operation is performed for reactive power):

[TABLE]

Then, using the active/reactive power injection PDFs, $E\{\tilde{p}_{j}\}$ and $E\{\tilde{q}_{j}\}$, can be obtained through integration:

[TABLE]

[TABLE]

where, $P_{j}$ and $Q_{j}$ denote the mean values of the available active and reactive powers at node $j$, respectively ($P_{j}=E\{p_{j}\}$ and $Q_{j}=E\{q_{j}\}$). Thus, the mean nodal voltage magnitude squared can be written in terms of inverter switch statistics and expected PV/load available powers:

[TABLE]

Stage 2: $\sigma_{v_{i}}^{2}$ Parameterization - Using (4), the variance of nodal voltage magnitude squared can be formulated as,

[TABLE]

where, $\sigma_{\tilde{v}_{ij}}^{2}$ is the variance of $\tilde{v}_{ij}$, and the operator $\Omega\{x_{1},x_{2}\}$ denotes the covariance of the two random variables $x_{1}$ and $x_{2}$, which itself can be written in terms of their correlation, $\rho_{x_{1},x_{2}}$, and standard deviations, $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$, as $\Omega\{x_{1},x_{2}\}=\rho_{x_{1},x_{2}}\sigma_{x_{1}}\sigma_{x_{2}}$. To fully parameterize $\sigma_{v_{i}}^{2}$ using available load/PV power statistics, $\sigma_{\tilde{v}_{ij}}^{2}$ and $\Omega\{\tilde{v}_{ij},\tilde{v}_{ik}\}$ have to be determined separately.

Stage 2-1: $\sigma_{\tilde{v}_{ij}}^{2}$ Parameterization - Using (5), $\sigma_{\tilde{v}_{ij}}^{2}$ is formulated as a function of $\tilde{p}_{j}$ and $\tilde{q}_{j}$ statistics:

[TABLE]

where, $\sigma_{\tilde{p}_{j}}^{2}$ and $\sigma_{\tilde{q}_{j}}^{2}$ are the active/reactive power injection variances, which can in turn be determined as follows:

[TABLE]

where, $E\{s_{j}^{2}p_{j}^{2}\}$ is calculated through a similar process involved in (14)-(17) (i.e., obtain the CDF, determine the PDF, and integrate). Noting that in our case $s_{j}^{2}=s_{j}$, the PDF of $s_{j}^{2}p_{j}^{2}$ is derived as follows (similar derivation applies to $s_{j}^{2}q_{j}^{2}$):

[TABLE]

By integrating (22) and using (16)-(17) to substitute for $E\{\tilde{p}_{j}\}$ and $E\{\tilde{q}_{j}\}$, the following results are obtained to parameterize the variances of nodal active/reactive power injections:

[TABLE]

[TABLE]

where, $P_{j}^{+}=E\{p_{j}^{2}|p_{j}\geq 0\}$ and $P_{j}^{-}=E\{p_{j}^{2}|p_{j}<0\}$; similar definitions apply to $Q_{j}^{+}$ and $Q_{j}^{-}$. Note that given that $p_{j}\geq 0$ for PVs, $P_{j}^{+}=\sigma_{p_{j}}^{2}+P_{j}^{2}$ and $P_{j}^{-}=0$. Employing an analogous logic, $P_{j}^{+}=0$ and $P_{j}^{-}=\sigma_{p_{j}}^{2}+P_{j}^{2}$ for loads.

To obtain $\Omega\{\tilde{p}_{j},\tilde{q}_{j}\}$ in (20), we leverage the fact that $\Omega\{x_{1},x_{2}\}=E\{x_{1}x_{2}\}-E\{x_{1}\}E\{x_{2}\}$ as follows:

[TABLE]

where, the term $E\{\tilde{p}_{j}\tilde{q}_{j}\}$ is calculated similar to previous derivations (i.e., CDF$\rightarrow$PDF$\rightarrow$integration), which combined with (16) and (17) yields the following result:

[TABLE]

where, $\Omega\{p_{j},q_{j}\}$ can be determined in terms of available active/reactive power statistics, including correlation and standard deviations as $\Omega\{p_{j},q_{j}\}=\rho_{p_{j},q_{j}}\sigma_{p_{j}}\sigma_{q_{j}}$.

Thus, using (23), (24), and (26), $\sigma_{\tilde{v}_{ij}}^{2}$ can be parameterized in terms of the available active/reactive power statistics, and with respect to micro-states:

[TABLE]

where, the time-invariant parameters $\Gamma_{ij}^{1}$ and $\Gamma_{ij}^{2}$ are given below:

[TABLE]

[TABLE]

Stage 2-2: $\Omega\{\tilde{v}_{ij},\tilde{v}_{ik}\}$ Parameterization - Similar to (26), $\Omega\{\tilde{v}_{ij},\tilde{v}_{ik}\}$, is broken down to its components:

[TABLE]

By adopting a CDF$\rightarrow$PDF$\rightarrow$integration strategy, $E\{\tilde{v}_{ij}\tilde{v}_{ik}\}$ is determined in terms of active/reactive power injection statistics as follows:

[TABLE]

where, using previous derivations and through algebraic manipulations, the following parameterization is obtained in terms of available active/reactive power statistics for $\Omega\{\tilde{v}_{ij},\tilde{v}_{ik}\}$:

[TABLE]

where, the parameters $\Gamma_{ijk}^{1}$ and $\Gamma_{ijk}^{2}$ are determined as:

[TABLE]

[TABLE]

By substituting (32) and (27) into (19), $\sigma_{v_{i}}^{2}$ is now fully determined as a function of micro-states and in terms of available nodal active/reactive power statistics.

Stage 3. Probabilistic Inverter Tripping Model Representation: Finally, using the parameterized $\sigma_{v_{i}}^{2}$ and $\mu_{v_{i}}$, the probabilistic model (12) yields a the following bilinear matrix-form representation for the approximate micro-state vector $\boldsymbol{\hat{\lambda}}=[\hat{\lambda}_{1},...,\hat{\lambda}_{N}]^{\top}$:

[TABLE]

where, all the time-invariant parameters of the model are concatenated into the vector $\boldsymbol{a_{0}}$ and matrices $B$, and $\{C_{1},...,C_{N}\}$. The elements of these parameters are determined by organizing the previous derivations in Stages 1 and 2, as follows:

[TABLE]

[TABLE]

[TABLE]

where, $\boldsymbol{a_{0}}(i)$ denotes the $i$’th element of $\boldsymbol{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}$, respectively. The aggregate switching equation can be written as a function of approximate macro-state, $\hat{S}=\sum_{i=1}^{N}\hat{\lambda}_{i}$, as follows:

[TABLE]

where, $\hat{S}$ is a conservative estimator of the real macro-state, $S$, which is the actual expected population of inverter that are ON, i.e., $\hat{S}(t)\leq S(t)$. Also, $B(i,:)$ is the $i$’th row of matrix $B$. To summarize, the proposed approximate probabilistic model leverages available load/PV power statistics shown in Table I. Previous works have used various data-driven and machine learning methods that can be applied for obtaining statistical properties of nodal load/PV powers in partially observable networks from limited available data (for example see [22, 23, 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 inverters in weak grids [14, 25]: overloading/faults on one feeder can trigger the voltage protection mechanism of inverters on a healthy neighboring feeder. The sympathetic tripping of inverters is also caused by dependencies in nodal voltages within the distribution grid (i.e., excessive 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.

II-D Discussion on Probabilistic Tripping Model Properties

The probabilistic model (35) represents a set of self-consistent equations; in other words, any $\boldsymbol{\hat{\lambda}}$ that satisfies these equations is a conservative estimator of probability of inverter tripping. Furthermore, this probabilistic model can be thought of as the asymptotic equilibrium of an abstract discrete dynamic system:

[TABLE]

where, the equilibrium is achieved at $\boldsymbol{\hat{\lambda}}(k+1)=\boldsymbol{\hat{\lambda}}(k)$ and coincides with the solution of the proposed probabilistic model (35). This abstract dynamic system has an intuitive interpretation: matrix $B$ represents the linear component of the dynamics, which as can be observed in (II-C) and (29), is determined only by each individual nodes’ active/reactive power statistics, including the expected values and self-correlation between active/reactive power at each node alone. However, matrices $\{C_{1},...,C_{N}\}$ capture the nonlinear components of the dynamic system, where the element $C_{i}(j,k)$ determines the coefficient assigned to the interactive nonlinear probability-product term $\hat{\lambda_{j}}(t)\cdot\hat{\lambda}_{k}(t)$ in driving $\hat{\lambda_{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 (II-C) the elements of $C_{i}$, unlike $B$, are determined by the mutual correlations in available active/reactive 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 evolve over time, which could potentially result into a cascading tripping event, as pointed out in [4, 7, 18, 19]. A regime shift at the equilibrium of the abstract nonlinear dynamic system basically corresponds to qualitative changes in the solution of our probabilistic tripping model, potentially, leading to a sudden increase in the average chances of voltage-driven 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 known to demonstrate critical behavior and emergent non-trivial patterns at the macro-level under certain conditions.

An important factor in tripping studies is the impact of setting of inverter protection systems. This can be seen in equations (36), (II-C), and (38) that present the parameters of the proposed statistical model. Specifically, parameter $C_{i}(j,k)$ in (38), captures the joint impacts of inverter $j$ and inverter $k$ on probability of tripping for inverter $i$. As can be seen, the absolute value of this parameter decreases with $\frac{1}{(v_{max}-v_{min})^{2}}$. Thus, increasing the inverter upper protection threshold, $v_{max}$, or decreasing the lower protection threshold, $v_{min}$, (i.e., making the inverter less sensitive to voltage events) will result in a decline in mutual impacts of inverters on each other. In other words, relaxing the protection boundaries significantly weakens the dependencies among inverter tripping. If $C_{i}(j,k)$ is thought of as a measure of strength of interdependency among inverters, then our model suggests that loss of interdependency is approximately proportional to the inverse of inverter protection dead-band width squared.

II-E Integrating Voltage-Dependent Resources Into the Proposed Probabilistic Tripping Model

Note that so far we have assumed that the nodal active and reactive power injections, $p_{j}$ and $q_{j}$, are external inputs to the model. However, active and reactive power injection of certain nodes can show high levels of voltage-dependency and cannot be treated as external inputs. The voltage-dependency can be caused by reactive power support from the inverters or load power voltage-sensitivity. In this section, we will demonstrate that voltage-dependent resources can also be included in our probabilistic model. To do this, the active/reactive power injections are linearized around the nominal squared voltage ($v_{n}$):

[TABLE]

[TABLE]

The active/reactive power 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 conveniently 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 active/reactive power voltage-sensitivity values. For example, this sensitivity can be obtained for ZIP loads [28] and inverters that are capable of reactive power support [18, 19] as follows:

[TABLE]

[TABLE]

where, $B_{j}$ and $C_{j}$ represent the ZIP coefficients corresponding to the fixed-current and fixed-impedance portions of ZIP load, respectively, and $k_{j}<0$ is the local inverter droop coefficient. Given the voltage-sensitivity values, the second terms in (41) and (42) simply serve as new additional nodal active/reactive power injections and can be treated in the model similar to other loads. For example, the surrogate nodal active/reactive injections for ZIP loads and inverters with reactive support capability can be conservatively estimated as follows:

[TABLE]

[TABLE]

where, $\bar{v}$ denotes a conservative user-defined value that can be used by the utilities to model worst-case tripping scenarios. However, note that (45) and (46) are still conservative estimations. Developing more accurate models for integrating voltage-dependent power injection into tripping equations remains the subject of future research.

III Solar Curtailment Quantification and Mitigation

Using (35) as a conservative 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 available nodal active/reactive power statistics are known:

[TABLE]

where, $\boldsymbol{P}=[P_{1},...,P_{N}]^{\top}$. The objective of this optimization problem is to find the maximum achievable expected solar power in the gird according to the conservative statistical model. While the solution to this problem is still a lower bound estimation of the real achievable PV power, the estimation gap between $\boldsymbol{\hat{\lambda}}$ and $\boldsymbol{\lambda}$ is minimized. In other words, the optimization searches for the most optimistic values for micro-states with respect to the conservative approximate probabilistic tripping model. The problem is constrained by the matrix equations that govern the probabilities of inverter tripping. Furthermore, the physical characteristics of micro-states are constrained by valid probability assignments within $[0,1]$ interval.

A similar problem can be formulated to provide counter-measures against massive tripping events at any given time window. In general, the proposed statistical tripping model can be integrated as a constraint into any volt-var optimization formulation [29, 30, 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]:

[TABLE]

where, $v_{0}$ is integrated into the optimization problem as a decision variable. 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_{0}^{I}$ 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 bus with respect to the initial voltage setpoint. To integrate $v_{0}$ into the problem, the expected nodal voltage magnitude squared values are written as a function of network parameters, expected available nodal active/reactive powers, and the optimization decision variables using (18):

[TABLE]

where, $\boldsymbol{v_{0}}=[v_{0},...,v_{0}]^{\top}$.

Despite its convenient differentiable matrix-form formulation, the probabilistic tripping model introduces quadratic non-convex 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]. To handle the non-convexity, these methods generally define an auxiliary matrix, $\Lambda=\boldsymbol{\hat{\lambda}}\boldsymbol{\hat{\lambda}}^{\top}$, which enables obtaining a convex surrogate for the original problem. For example, by applying parabolic relaxation, the constraints defined by the model are replaced with the following alternative constraints:

[TABLE]

[TABLE]

[TABLE]

where, $C_{i}\bullet\Lambda=\sum_{n=1}^{N}\sum_{m=1}^{N}\{C_{i}(n,m)\Lambda(n,m)\}$, $I_{N}$ is an $N\times N$ identity matrix, and $\boldsymbol{\epsilon}^{+}$/$\boldsymbol{\epsilon}^{+}$ are positive/negative small-valued slack variables that are used for transforming equality constraints defined by the model into two equivalent inequality constraints. The obtained inequalities (57)-(59) are convex constraints with respect to variables $\boldsymbol{\hat{\lambda}}$ and $\Lambda$.

IV Numerical Experiments and Validation

Numerical experiments have been performed to validate the proposed probabilistic tripping model. In this, we have 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 available online [35]. To 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.

To verify the performance of the proposed approximate statistical model, extensive time-series simulations were performed on the test system under various loading and solar generation scenarios over a course of day. Then, the real values of original micro-states, $\lambda_{i}$, were determined empirically over time windows of length $T=60$ minutes. Intuitively, $\lambda_{i}$ serves as the ground truth and roughly represents the portion of time that $s_{i}$ is ON during each time window:

[TABLE]

Thus, we have two distinct time windows throughout numerical 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, $\lambda_{i}$, at different time intervals, which are determined by applying (60) to simulation outcomes. Based on the values of these micro-states, the empirical macro-state value is calculated at all time intervals, which represents the expected percentage of PV switches in ON state, i.e., $S_{p}(T)\approx\frac{\sum_{i=1}^{N}\lambda_{i}(T)}{N}\times 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 performance of the method. This figure also shows another lower bound obtained by simply using maximum PV capacities and assuming zero nodal consumption. However, as can be seen, this lower bound gives fixed over-conservative outcomes that do not reflect the true conditions of the system and have no correlation with the time-series PV/load data. Fig. 5c depicts the aggregate maximum available solar power (all switches ON at all time), empirical aggregate realized solar power from numerical simulations (60), and solar power corresponding to solution of (47). As observed, the lower bound solution still holds and provides a conservative yet close estimation for the empirical achievable solar power outcome. Fig. 6 compares the empirical and model-based probabilities of inverter tripping in a heavy-loaded time interval. Unlike the previous case, these tripping probabilities are due to under-voltages. As can be observed, the model still provides a conservative lower bound on the probability of tripping. Note that the reason for higher levels 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 lower 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 over peak nodal demand. The maximum, minimum, and mean values of the gap between the provided lower bound and the empirical macro-state is measured at various levels of PV penetration. As is observed in the figure, the optimistic value of the gap drops and eventually reaches $5\%$ as PV penetration increases, which indicates that the lower bound approaches the true macro-state value in grids with higher PV penetration. On the other hand, the maximum value of the gap shows an increase after a certain PV penetration level which points out to higher variations in solutions obtained from the probabilistic model.

Fig. 8 shows the overall daily solar curtailment levels, both empirical and the lower bound, as a function of changes in inverter control parameter. The inverters in the system are assumed to be controlled in constant power factor (PF) mode. As the reference PF setpoint increases and the system moves towards unity PF the voltage fluctuations increase, which leads to higher solar curtailment. This confirms previous observations in the literature [9]. Furthermore, our proposed probabilistic lower bound always slightly over-estimates the curtailment level, as expected correctly from the conservative 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 levels. Fig. 9b compares the aggregate solar power injection values under the newly acquired $v_{0}$ values and the base case without voltage control. As can be seen, the obtained countermeasure has assisted significantly in mitigating the overall solar power curtailments during critical time intervals.

We have 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) Voltage regulation is applied with the objective of minimizing the average 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 inverters (around 13%) have tripped due to under-voltage protection during later hours of the day. This has resulted in a loss of renewable power injection in the grid (Fig. 10b). However, by applying voltage regulation using the proposed tripping model we have been able to maintain the voltages much closer to their nominal values (see Fig. 10c) and prevent tripping events and loss of solar generation resources altogether. Note that Fig. 10c shows the average value of nodal voltages across the whole system; thus, while most of the nodes maintain healthy voltage levels (as they should), the excessive loading on weak system lines under the baseline has resulted to a temporary voltage drop below inverters’ protection activation threshold, which has engaged their under-voltage protection devices. 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 extensive tripping, to a new state, in which the average probability of solar curtailment steadily increases and extended tripping events can be expected. The existence of this threshold attests to a stage transition in the extent of switching events, 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 observed that the proposed statistical lower bound accurately tracks the behavior of the real system, and can be used to convey 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 Conclusions

In this paper, a probabilistic model of interdependent solar inverter tripping is presented to assess the risk of solar power curtailments due to voltage fluctuations in distribution grids. This model is developed using only the statistical properties of available load/PV active/reactive power. Numerical results on a real distribution feeder using real data successfully validate the estimated conservative lower bounds on inverter micro-states. Furthermore, it is demonstrated that the proposed model can be used for identifying regime shifts in tripping events and designing countermeasures to minimize risk of solar power curtailment. As a future research direction, we will explore integrating the more dynamic functions of inverter control and protection, including ride-through capabilities, [17, 18, 19, 20] into the probabilistic tripping model. For example, the proposed statistical lower bound, which is based on Chebyshev’s inequality, might become too conservative over short time windows if inverters’ disturbance ride-through capabilities are activated. A less conservative lower bound that incorporates all aspects of inverter behavior will enable operators to monitor the sequence and transitions of tripping events, and mitigate potential cascading failure of resources.