Singular Perturbation and Small-signal Stability for Inverter Networks

TL;DR

This paper develops a singular perturbation-based framework to analyze small-signal stability in inverter-dominated electrical networks, providing efficient stability conditions that consider inverter dynamics and network topology.

Contribution

It introduces a systematic stability analysis method using singular perturbation theory, accounting for inverter dynamics and network structure, which improves upon full eigenvalue analysis.

Findings

Analytic sufficient stability condition reduces computational effort.

Network topology and inverter parameters significantly influence stability.

Numerical simulations validate the effectiveness of the proposed approach.

Abstract

This paper examines small-signal stability of electrical networks composed dominantly of three-phase grid-following inverters. We show that the mere existence of a high-voltage power flow solution does not necessarily imply small-signal stability; this motivates us to develop a framework for stability analysis that systematically acknowledges inverter dynamics. We identify a suitable time-scale decomposition for the inverter dynamics, and using singular perturbation theory, obtain an analytic sufficient condition to verify small-signal stability. Compared to the alternative of performing an eigenvalue analysis of the full-order network dynamics, our analytic sufficient condition reduces computational complexity and yields insights on the role of network topology and constitution as well as inverter-filter and control parameters on small-signal stability. Numerical simulations for a radial network validate the approach and illustrate the efficiency of our analytic conditions for designing and monitoring grid-tied inverter networks.

Figures

Tables

Table 1. TABLE I: Comparing the computation time and accuracy of conditions ( 18 ) and ( 15 ) for small-signal stability. The unit of the quantities T test subscript T test \mathrm{T}_{\textup{test}} and T lin subscript T lin \mathrm{T}_{\textup{lin}} are seconds, the quantity of p ^ ^ 𝑝 \widehat{p} is dimensionless and the quantities SPL SPL \mathrm{SPL} and SPL test subscript SPL test \mathrm{SPL}_{\mathrm{test}} and SPL static subscript SPL static \mathrm{SPL}_{\mathrm{static}} are integers.

$\hat{p}$	$T_{\text{lin}}\left(s\right)$	$T_{\text{test}}\left(s\right)$	$\mathrm{SPL}$	${\mathrm{SPL}}_{\mathrm{test}}$	${\mathrm{SPL}}_{\mathrm{static}}$
0.80	0.5853	0.0774	22	20	35
1.00	0.5636	0.0793	22	20	31
1.20	0.5546	0.0741	22	20	28
1.40	0.6525	0.0767	22	20	26
1.60	0.0767	0.5313	21	20	24
1.80	0.5063	0.0888	21	20	23
2.00	0.5241	0.0766	21	20	22

Table 2. TABLE II: Variables and their symbols for the inverter model.

Variable	Symbol	Variable	Symbol
PLL low-pass filter state	$v_{\text{PLL}}$	Current controller auxiliary state	${\gamma }_{\text{dq}}$
PLL PI controller state	${\varphi }_{\text{PLL}}$	Current controller output voltage	$v_{\text{idq}}$
PLL phase output	$\delta$	Output current	$i_{\text{odq}}$
Power controller low-pass filter state	$s_{\text{ave}}$	Power controller auxiliary state	${\varphi }_{\text{s}}$
Reference power injection	$s^{\mathrm{ref}}$	Current controller reference current	$i_{\text{ldq}}$
Power controller auxiliary state	${\varphi }_{\text{s}}$	PLL Frequency	${\omega }_{\text{PLL}}$
Current in the lines	${\xi }_{\text{DQ}}$	LC-filter voltage	$v_{\text{odq}}$
Grid voltage	$v_{\text{gDQ}}$	Grid frequency	${\omega }_{\text{nom}}$

Table 3. TABLE III: Dimensionless parameters of the inverter model.

Parameter	Symbol	Values from [24]	Values from  [17]
Grid frequency	${\omega }_{\text{nom}}$	$377$ rad/s	$377$ rad/s
Grid voltage amplitude	$V_{\text{g}}$	$169$ V	$169$ V
PLL time constant	${\tau }_{\text{PLL}}$	$1.27e-5$	$1.27e-5$
PLL time constant	${\tau }_{\text{PLL}}^{\prime }$	$2.36e-2$	$4.7e-3$
PLL time constant	$T_{\text{PLL}}$	$1.25e-1$	$1.25e-1$
Power controller time constant	${\tau }_{\text{s}}$	$1.99e-2$	$1.99e-2$
Power controller time constant	${\tau }_{\text{s}}^{\prime }$	$16.97$	$16.97$
Power controller time constant	$T_{\text{s}}$	$1.00e-1$	$1.00e-1$
Current controller time constant	${\tau }_{\text{c}}$	$1.70e-3$	$7.85e-4$
Current controller time constant	$T_{\text{c}}$	$1.00e-2$	$1.43e-3$
LC time constant	${\tau }_{\text{LC}}$	$2.37e-5$	$2.54e-5$
LC time constant	${\tau }_{\text{LC}}^{\prime }$	$2.37e-5$	$2.54e-5$
Line time constant	${\tau }_{\text{e}}$	$1.90e-3$	$2.70e-3$
Line time constant	${\tau }_{\text{e}}^{\prime }$	$7.40e-3$	$1.31e-0$

Equations

$$\mathcal{D}(\mathrm{u})=\left(\begin{array}{cc}{u}_1& u_2\\ u_2& -u_1\end{array}\right),{\mathcal{D}}^{\prime }(\mathrm{u})=\left(\begin{array}{cc}{u}_1& u_2\\ -u_2& u_1\end{array}\right).$$

$${\dot{v}}_{\mathrm{PLL}}$$

$${\dot{\varphi }}_{\mathrm{PLL}}$$

$$\dot{\delta }$$

$${\dot{s}}_{\mathrm{avg}}$$

$${\tilde{i}}_{l\mathrm{dq}}$$

$$s=\frac{3}{2}\left(\begin{array}{c}{v}_{\mathrm{od}}{i}_{\mathrm{od}}+v_{\mathrm{oq}}{i}_{\mathrm{oq}}\\ v_{\mathrm{oq}}{i}_{\mathrm{od}}-v_{\mathrm{od}}{i}_{\mathrm{oq}}\end{array}\right)=\frac{3}{2}\mathcal{D}(v_{\mathrm{odq}})i_{\mathrm{odq}}.$$

$$\begin{array}{rl}{\tilde{v}}_{\mathrm{idq}}& ={\mathrm{kp}}_{\text{c}}\left({\tilde{i}}_{\mathrm{ldq}}-i_{\mathrm{ldq}}\right)+{\mathrm{ki}}_{\text{c}}\int \left({\tilde{i}}_{\mathrm{ldq}}-i_{\mathrm{ldq}}\right)dt\\ & +{\omega }_{\mathrm{PLL}}{L}_{\text{f}}(\mathcal{J}{i}_{\text{ldq}}).\end{array}$$

$${\dot{i}}_{\mathrm{ldq}}$$

$${\dot{v}}_{\mathrm{odq}}$$

$${\gamma }_{\mathrm{DQ}}=\mathcal{R}(-\delta ){\gamma }_{\mathrm{dq}},i_{\mathrm{lDQ}}=\mathcal{R}(-\delta )i_{\mathrm{ldq}},v_{\mathrm{oDQ}}=\mathcal{R}(-\delta )v_{\mathrm{odq}}.$$

$${\dot{y}}_{\text{DQ}}-\mathcal{J}\mathcal{R}(-\delta )\dot{\delta }{y}_{\text{DQ}}=\mathcal{R}(-\delta ){\dot{y}}_{\text{dq}}.$$

$$\begin{array}{r}\dot{y}=f(y)+g(y)i_{\text{oDQ}}+C{s}^{\mathrm{ref}},\end{array}$$

$$\displaystyle Y=\scalebox{1.0}{$\begin{bmatrix}Y_{\textup{00}}&Y_{\textup{0L}}&Y_{\textup{0I}}\\ Y_{\textup{L0}}&Y_{\textup{LL}}&Y_{\textup{LI}}\\ Y_{\textup{I0}}&Y_{\textup{IL}}&Y_{\textup{II}}\end{bmatrix}$}.$$

$$\begin{array}{r}\dot{\mathbf{y}}=F(\mathbf{y})+G(\mathbf{y}){\mathrm{i}}_{\mathrm{oDQ}}+C{\mathbf{s}}^{\mathrm{ref}},\end{array}$$

$$\displaystyle\mathbf{v}_{\textup{L}}=\big{(}-[\mathbf{R}_{\textup{L}}]\otimes I_{2}\big{)}\mathbf{i}_{\textup{L}}.$$

$$([\mathbf{L}_{\mathcal{E}}]\otimes I_{2})\dot{\xi}_{\textup{DQ}}=\big{(}-[\mathbf{R}_{\mathcal{E}}]\otimes I_{2}-\omega_{\textup{nom}}[\mathbf{L}_{\mathcal{E}}]\otimes\mathcal{J}\big{)}\xi_{\textup{DQ}}\\ +(B^{\top}\otimes I_{2})\mathbf{v},$$

$$\dot{\mathbf{y}}$$

$$([{\mathbf{L}}_{\mathcal{E}}]\otimes I_2){\dot{\xi }}_{\text{DQ}}$$

$$+(B^{\top }\otimes I_2)\mathbf{v},$$

$${\hat{v}}_{\text{PLL}}$$

$${\hat{\varphi }}_{\text{s}}$$

$${\hat{i}}_{\text{lDQ}}$$

$${\hat{\xi }}_{\text{DQ}}:=V_{\text{g}}{s}_{\text{nom}}^{-1}{\xi }_{\text{DQ}},\hat{\mathbf{v}}:=V_{\text{g}}^{-1}\mathbf{v},\hat{\mathbf{i}}:=V_{\text{g}}{s}_{\text{nom}}^{-1}\mathbf{i}.$$

$${\dot{\hat{\mathbf{v}}}}_{\mathrm{PLL}}=[{\mathbf{\tau }}_{\text{PLL}}{]}^{-1}({\hat{\mathbf{v}}}_{\mathrm{od}}-{\hat{\mathbf{v}}}_{\mathrm{PLL}}),$$

$${\dot{\hat{\mathbf{\varphi }}}}_{\mathrm{PLL}}=-[{\mathbf{T}}_{\text{PLL}}{]}^{-1}{\hat{\mathbf{v}}}_{\mathrm{PLL}},$$

$$\dot{\mathbf{\delta }}=[{\mathbf{\tau }}_{\text{PLL}}^{\prime }{]}^{-1}\left({\hat{\mathbf{\varphi }}}_{\mathrm{PLL}}-{\hat{\mathbf{v}}}_{\mathrm{PLL}}\right),$$

$${\dot{\hat{\mathbf{s}}}}_{\mathrm{avg}}=[{\mathbf{\tau }}_{\text{s}}\otimes I_2{]}^{-1}(\hat{\mathbf{s}}-{\hat{\mathbf{s}}}_{\mathrm{avg}}),$$

$$\dot{{\hat{\mathbf{\varphi }}}_{\text{s}}}=[{\mathbf{\tau }}_{\text{s}}^{\prime }\otimes I_2{]}^{-1}({\hat{\mathbf{s}}}^{\mathrm{ref}}-{\hat{\mathbf{s}}}_{\mathrm{avg}}),$$

$${\dot{\hat{\mathbf{\gamma }}}}_{\text{DQ}}=[{\mathbf{\tau }}_{\text{c}}\otimes I_2{]}^{-1}\left({\tilde{\mathbf{i}}}_{\text{lDQ}}-{\hat{\mathbf{i}}}_{\text{lDQ}}\right)+\mathcal{J}[\dot{\mathbf{\delta }}]{\hat{\mathbf{\gamma }}}_{\text{DQ}},$$

$${\dot{\hat{\mathbf{i}}}}_{\mathrm{lDQ}}=([{\mathbf{\tau }}_{\text{LC}}{]}^{-1}[\mathcal{X}]\otimes I_2)\left({\hat{\mathbf{v}}}_{\mathrm{lDQ}}-{\hat{\mathbf{v}}}_{\mathrm{oDQ}}\right)+\mathcal{J}[\dot{\mathbf{\delta }}]{\hat{\mathbf{i}}}_{\mathrm{lDQ}},$$

$$\displaystyle\dot{\widehat{\mathbf{v}}}_{\mathrm{oDQ}}=([\boldsymbol{\tau}^{\prime}_{\textup{LC}}][\mathcal{X}]\otimes I_{2})^{-1}\big{(}\widehat{\mathbf{i}}_{\mathrm{lDQ}}-\widehat{\mathbf{i}}_{\mathrm{oDQ}}$$

$$\displaystyle\hskip 113.81102pt-[\boldsymbol{\tau}^{\prime\prime}_{\textup{LC}}\otimes I_{2}]\mathcal{J}\widehat{\mathbf{v}}_{\mathrm{oDQ}}\big{)},$$

$$\displaystyle\dot{\widehat{\xi}}_{\textup{DQ}}=\big{(}[\boldsymbol{\tau}_{\mathcal{E}}][\boldsymbol{\tau}^{\prime}_{\mathcal{E}}]\otimes I_{2}\big{)}^{-1}\big{(}(B^{\top}\otimes I_{2})\mathbf{\widehat{v}}-\mathcal{Z}\widehat{\xi}_{\textup{DQ}}\big{)}.$$

$${\tilde{\mathbf{i}}}_{\mathrm{lDQ}}$$

$${\hat{\mathbf{v}}}_{\mathrm{lDQ}}$$

$$\hat{\mathbf{i}}$$

$$\mathcal{Z}$$

$$\mathcal{X}$$

$$\hat{\mathbf{v}}$$

$$\begin{array}{rl}{Y}_{\text{red}}& :=Y_{\text{II}}-Y_{\text{IL}}{\left(Y_{\text{LL}}+[{\mathbf{R}}_{\text{L}}{]}^{-1}\otimes I_2\right)}^{-1}{Y}_{\text{LI}},\\ Y_{\text{C,red}}& :=Y_{\text{red}}+{\omega }_{\text{nom}}\left([{\mathbf{C}}_{\text{f}}]\otimes J\right),\\ Y_{\text{g}}& :=Y_{\text{I0}}-Y_{\text{IL}}{\left(Y_{\text{LL}}+[{\mathbf{R}}_{\text{L}}{]}^{-1}\otimes I_2\right)}^{-1}{Y}_{\text{L0}},\\ \mathbf{w}& :=-Y_{\text{red}}^{-1}{Y}_{\text{0g}}{v}_{\text{gDQ}},\end{array}$$

Full text

I Introduction

Problem description and motivation

The ongoing shift from fossil-fuel-driven synchronous generators to power-electronics-interfaced renewable energy is leading to changes in how power grids are modeled, analyzed, and controlled. While synchronous generators are generally rated at several hundreds of MVA and installed on the transmission backbone, power electronics inverters are distributed across both transmission and distribution subsystems and are generally much smaller in capacity. Furthermore, synchronous generators have large rotating masses that buffer supply-demand fluctuations and limit frequency excursions during transients, whereas inverters have very different dynamics (attributable dominantly to output filters and digital controllers [30]) and they possess no moving parts. Future grids will have large number of power-electronics-interfaced generations with highly distributed architecture as inverters assume a more prominent role, and this will necessitate the development of compatible models and computationally efficient analysis approaches to certify stability.

Most commercial inverters on the market for residential and utility-scale applications are grid-following. This means inverters inject currents while synchronized to the voltage at their terminals which is assumed to be set externally by the bulk grid. In effect, grid-following inverters act as voltage-following current sources. In recent years, there has been increased attention on grid-forming inverter technology, whereby—much like conventional synchronous generators—terminal voltage and frequency are modulated as a function of real- and reactive-power injections. Indeed, while grid-forming technology may very well be the solution for grids with penetration levels approaching $100\%$, in the foreseeable future, it is likely that the bulk grid will see conventional synchronous generators co-exist alongside a large number of grid-following inverters. For instance, in Oahu, Hawaii, at least $800,000$ micro-inverters interconnect photovoltaic panels to the grid, producing as much power as the state’s largest conventional power plant, Konkar [9]. This motivates the problem of stability analysis for large-scale networks composed dominantly of grid-following inverters that we examine in this work.

The inverter dynamical model that we investigate is composed of a current controller, a power controller, a PLL, and an LCL filter. This is prototypical and mirrors models published widely in the literature [16, 17, 18, 19, 20, 23, 24, 27]. We propose a framework that leverages singular perturbation methods to obtain an analytic, computationally light-weight condition for small-signal stability assessment of three-phase distribution networks with grid-following inverters. Our solution strategy yields an analytic condition for stability that is agnostic to system size and clearly highlights the role of the network topology and pertinent system parameters on system stability. Due to the complexity of the involved dynamics, most prior art on stability of inverter-based systems has typically focused on simplified models that neglect inner control loops which underpin fast dynamics [26] (see the discussions in [5]). Some exceptions are [1, 31], where stability of full-order inverter models are studied in grid-connected networks. However, these studies are restricted to parallel networks of inverters and are not applicable to networks with general topologies. In some cases, detailed inverter models have been considered in general networks, but system stability has only been studied for a single inverter or a small network of inverters using numerical eigenvalue analysis. For instance, small-signal stability is analyzed in the literature using eigenvalue analysis for a single grid-following inverter under unintentional islanding in [28], for the IEEE $37$-bus system with $7$ inverters in [23, 24], for a radial network consisting of $3$ inverters in [17], and for a single-machine-single-inverter network in [14]. Understandably, while numerical eigenvalue analysis of the linearized network is indeed a reasonable strategy, this approach comes with significant drawbacks. First, the large size of the network combined with the high dimensionality of the inverter model pose computational challenges to analysis. Furthermore, studying small-signal stability by computing eigenvalues of the linearized system does not reveal the role of critical network attributes on system stability, insights, which if formalized appropriately can facilitate analysis and design.

On a tangential note, it must be acknowledged that there is a growing body of work on stability assessment of synchronous generators and grid-forming inverter systems, and this includes approaches that have applied model-order reduction using singular perturbation analysis. For instance, in [6, 10], a model-reduction approach based on singular perturbation is proposed to study stability and control of grid-forming inverters. In [15], singular perturbation is applied to obtain a hierarchy of reduced-order models for inverters in the grid-forming mode. In [24], a suitable time-scale decomposition for a class of inverters is identified and an iterative scheme for model order reduction is proposed. We refer interested readers to [25] for a survey on singular perturbation methods and to [26] for a survey on application of singular perturbation in stability and control of inverters.

Contributions

We make several contributions to the study of small-signal stability of grid-following inverter networks. First, we show that adopting a static model for grid-following inverters (as fixed sources of active and reactive power) and neglecting fast dynamics induced by the inverters’ control loops may lead to erroneous conclusions regarding stability (see Example 3). This underscores the importance of acknowledging a full-order model for stability analysis. We start by introducing a model for grid-tied inverter networks that acknowledges line dynamics and where each inverter is modeled using a $13$th-order model. Next, we uncover a correspondence between the equilibrium points of the dynamics and the solutions of the algebraic power-flow equations. The main contribution of this paper, i.e., an analytic sufficient condition for small-signal stability, is derived in the context of a dimensionless transcription of the involved models which leads to the identification of a physically insightful parametrization of the inverters. We show that certain assumptions on the range of parameters result in a time-scale decomposition of the system. Using singular perturbation analysis, we propose an analytic sufficient condition which guarantees small-signal stability over a given parametric regime. As a unique contribution, we emphasize that the dimensionless form of the network equations as well as the regularity of the singular perturbation problem (i.e., existence of isolated quasi-steady state manifolds) are critical steps in a rigorous time-scale analysis. Over this specified parametric regime, our analytic sufficient condition can also be interpreted as a lower bound on the stability threshold of the network and allows us to check system stability with minimal computational complexity. Furthermore, in the special case of resistive networks, our sufficient condition reduces to checking Hurwitzness of a Metzler matrix, something that can be implemented efficiently via linear programming (see [22]).

In the literature, small-signal stability of systems is usually studied using eigenvalue analysis for the full-order models (e.g., see [17]) or for the reduced-order models (e.g., see the survey [26]). Compared to performing eigenvalue analysis for the full-order system, our proposed analytically driven sufficient condition reduces computational complexity by addressing the high dimensionality of the underlying dynamics (the inverter model we study has $13$ dynamical states) and it demarcates the role of the network (topology and constitution) and pertinent inverter dynamics (filter and controller parameters). Moreover, while many existing approaches in the literature only outline iterative schemes for model reduction (see [24]), our analysis provides an explicit reduced-order model. Finally, we provide several numerical case studies that validate the analysis.

Paper organization

In Section II, we present the dynamical model for a class of three-phase grid-following inverters. In Section III, we derive an equivalent dimensionless description for a grid-tied network of inverters and loads and we study the equilibrium points of the system. In Section IV, we provide a sufficient condition for existence of a locally exponentially stable equilibrium point. Finally, in Section V, we illustrate some applications of the theoretical results in network design and stability assessment.

Notation

Vectors and matrices. We denote the set of real numbers by $\mathbb{R}$, the set of complex numbers by $\mathbb{C}$, the set of complex numbers with negative real part by ${\mathbb{C}}_{-}$, the set of binary $n$-tuples by ${\mathbb{Z}}^{n}_{2}$, and the $n$-dimensional torus by $\mathbb{T}^{n}$. We define $\mathrm{i}=\sqrt{-1}$. We identify the complex plane $\mathbb{C}$ with the real plane $\mathbb{R}^{2}$. For a complex number $v=v_{1}+\mathrm{i}v_{2}\in\mathbb{C}$, the norm of $v$ is $|v|=\sqrt{v_{1}^{2}+v_{2}^{2}}$ and the argument of $v$, $\mathrm{arg}(v)$, is the angle between $v$ and the positive imaginary axis. We denote the identity matrix of dimension $n$ by $I_{n}$, the $n$-column vector of zeros with $\mathbbold{0}_{n}$, and the $n$-column vector of ones with $\mathbbold{1}_{n}$. For a matrix $A=\{a_{ij}\}\in\mathbb{C}^{n\times m}$, we denote the trace by $\mathrm{tr}(A)$, the determinant by $\det(A)$, and $\infty$-norm by $\|A\|_{\infty}=\max_{i}\sum_{j=1}^{n}|a_{ij}|$. For two real symmetric matrices $A,B\in\mathbb{R}^{n\times n}$ we write $A\succ B$ if $A-B$ is positive definite. A real square matrix $A\in\mathbb{R}^{n\times n}$ is Metzler if all its off-diagonal entries are non-negative. For two square matrices $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{m\times m}$, the tensor product is denoted by $A\otimes B$. For a vector $\mathrm{x}\in\mathbb{C}^{n}$, we denote $\mathrm{diag}(\mathrm{x})$ by $[\mathrm{x}]$.

From $n$-complex variables to $2n$-real variables. For every complex $Z=X+\mathrm{i}Y\in\mathbb{C}$, the associated real variable in the real plane is denoted by $z=(x,y)^{\top}\in\mathbb{R}^{2}$. We will frequently use matrix $\mathcal{J}=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}$, matrix $\mathcal{H}=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$, and the rotation matrix (parameterized by angle $\theta\in\mathbb{S}^{1}$) by $\mathcal{R}(\theta)=\begin{pmatrix}\cos(\theta)&\sin(\theta)\\ -\sin(\theta)&\cos(\theta)\end{pmatrix}$. Let $\mathrm{u}\in\mathbb{R}^{2}$. We define matrix-valued operators $\mathcal{D}:\mathbb{R}^{2}\to\mathbb{R}^{2\times 2}$ and $\mathcal{D}^{\prime}:\mathbb{R}^{2}\to\mathbb{R}^{2\times 2}$ by

[TABLE]

All the above matrices and operators can be extended to the $n$-dimensional complex and $2n$-dimensional real spaces using block diagonal structure. For the brevity of notation, we denote the extended $n$-dimensional complex ($2n$-dimensional real) matrix/operator with the same symbol as its $1$-dimensional complex ($2$-dimensional real) counterpart. Let $\mathbf{V}\in\mathbb{C}^{n}$ and $\mathbf{v}\in\mathbb{R}^{2n}$ be the associated real vector, then we define $\|\mathbf{v}\|_{\mathbb{C},\infty}=\|\mathbf{V}\|_{\infty}$.

Algebraic graph theory. We denote an undirected weighted graph by a triple $G=(\mathcal{N},\mathcal{E},A)$, where $\mathcal{N}=\{1,2,\ldots,n\}$ is the set of nodes and $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of edges. The matrix $A=\{a_{ij}\}\in\mathbb{R}^{n\times n}$ is the weighted adjacency matrix. For every node $i\in V$, the degree of the node $i$ is given by $d_{i}=\sum_{j=1}^{n}a_{ij}$. For a fixed orientation on $G$, the incidence matrix of the graph $G$ is denoted by $B\in\mathbb{R}^{n\times m}$. The Laplacian for the graph $G$ is defined by $L=D-A$, where $D=\mathrm{diag}(d_{1},d_{2},\ldots,d_{n})$.

Power systems. We consider two different reference frames. The first is the so-called global $\mathrm{DQ}$-frame, and it is a rotating reference frame tied to the nominal grid frequency $\omega_{\textup{nom}}$. The second frame, which is usually referred as the local $\mathrm{dq}$-frame, is with reference to each inverter’s terminal voltage vector. For a balanced three-phase signal $x:\mathbb{R}_{\geq 0}\to\mathbb{R}^{3}$, we denote the $\mathrm{dq}$-frame representation by $x_{\mathrm{dq}}=(x_{\mathrm{d}},x_{\mathrm{q}})^{\top}$ and the $\mathrm{DQ}$-frame representation by $x_{\mathrm{DQ}}=(x_{\mathrm{D}},x_{\mathrm{Q}})^{\top}$. These are related as follows: $x_{\mathrm{dq}}=\mathcal{R}(\delta)x_{\mathrm{DQ}}$, where $\delta$ is the angle between the $\mathrm{dq}$-frame and the $\mathrm{DQ}$-frame. We assume that all the electrical signals in the network are balanced [2, Chapter 2]. Therefore, the voltage of the grid is a three-phase AC signal given by the time-varying function $v_{\textup{g}}(t)=[v_{\textup{ga}}(t),v_{\textup{gb}}(t),v_{\textup{gc}}(t)]^{\top}$ and in the global $\mathrm{DQ}$-frame, it is represented as $v_{\textup{gDQ}}=[0,V_{\textup{g}}]^{\top}$, where $V_{\textup{g}}$ is the amplitude of the grid voltage.

II Model of Individual Inverter

We briefly overview the dynamics of the type of grid-following $3$-phase inverters examined in this work. For a more detailed description of the model, see [20, 23]. The model captures all relevant AC-side dynamics, and is composed of a: i) phase-locked loop (PLL), ii) power controller, iii) current controller, and iv) $LC$ output filter. An illustrative block diagram is given in Fig. 1.

The PLL consists of a low-pass filter with cut-off frequency $\omega_{c,\mathrm{PLL}}$ and a PI controller with gains $\mathrm{kp}_{\mathrm{PLL}}$ and $\mathrm{ki}_{\mathrm{PLL}}$. The PLL dynamics are

[TABLE]

where $v_{\mathrm{PLL}}$ and $\phi_{\mathrm{PLL}}$ denote states of the low-pass filter and PI controller, respectively, $\dot{\delta}$ is the output of the PI controller, and $v_{\textup{od}}$ is the $d$-component of the output voltage of the inverter. The frequency of the PLL loop is defined by $\omega_{\mathrm{PLL}}=\omega_{\mathrm{nom}}+\dot{\delta}$. The power controller consists of two low-pass filters with cut-off frequency $\omega_{\textup{s}}$ and two PI controllers with gains $\mathrm{kp}_{\textup{s}}$ and $\mathrm{ki}_{\textup{s}}$. The pertinent dynamics are given by:

[TABLE]

where $s_{\mathrm{avg}}=[p_{\mathrm{avg}},q_{\mathrm{avg}}]^{\top}$ collects the states of the low-pass filters, $\tilde{i}_{l\mathrm{dq}}$ capture the outputs of the PI controllers (these are the references for the current controller), $s^{\mathrm{ref}}=(p^{\mathrm{ref}},q^{\mathrm{ref}})^{\top}$ collects the active- and reactive-power references, and $s=(p,q)^{\top}$ collects instantaneous active- and reactive-power outputs (measured at the point of common coupling):

[TABLE]

The current controller consists of two PI controllers with gains $\mathrm{kp}_{\textup{c}}$ and $\mathrm{ki}_{\textup{c}}$, with outputs to be the references for the inverter voltage at the switching terminals $v_{\mathrm{idq}}$:

[TABLE]

Since the switching period is typically much shorter than the filter and controller time constants, we assume that $v_{\mathrm{idq}}=\tilde{v}_{\mathrm{idq}}$. The dynamics of the $LC$ filter are given by

[TABLE]

Finally, we introduce two new variables $\phi_{\textup{s}}:=\int(s^{\mathrm{ref}}-s_{\textup{avg}})dt$, and $\gamma_{\textup{dq}}:=\int(\tilde{i}_{\textup{ldq}}-i_{\textup{ldq}})dt$, that will aid in exposition. The inner- and outer-loop control architecture examined here is ubiquitous, see, e.g., [16, 17, 18, 19, 20, 23, 24, 27], where similar models are utilized.

Now, we write the dynamical system of the inverters in the global $\mathrm{DQ}$-frame. To this end, we introduce:

[TABLE]

For vector $y$ with time-varying entries, the time derivatives in the global $\mathrm{DQ}$-frame and local $\mathrm{dq}$-frame are related by

[TABLE]

Leveraging this identity, the dynamical model for the grid-following inverter in the $\mathrm{DQ}$-frame can be expressed as:

[TABLE]

where $y=(v_{\mathrm{PLL}},\delta,\phi_{\textup{s}},s_{\textup{avg}},\gamma_{\mathrm{DQ}},{i}_{\mathrm{lDQ}},{v}_{\mathrm{oDQ}})^{\top}\in\mathbb{R}^{13}$ captures states of the inverter in the global $DQ$ reference frame, $s^{\mathrm{ref}}=(p^{\mathrm{ref}},q^{\mathrm{ref}})^{\top}\in\mathbb{R}^{2}$ captures the references for active and reactive power, $f:\mathbb{R}^{13}\to\mathbb{R}^{13}$ is the drift vector field, and $g:\mathbb{R}^{13}\to\mathbb{R}^{13\times 2}$, and $C\in\mathbb{R}^{13\times 2}$ are control vector fields. The mappings $f,g$ and the matrix $C$ are obtained from the dynamics governing the current controller, the power controller, the PLL, and the $LC$ filter outlined previously.

III Model for Grid-tied Network of Inverters

In this section, we derive the dynamical system model governing the grid-tied network of inverters and loads and study the equilibrium points of the system. We model the network using an undirected, connected, complex-weighted graph $G$ with node set (buses) $\mathcal{N}$, edge set (branches) $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$, and the symmetric matrix-valued edge weights (admittances) $a_{kj}=a_{jk}=(R_{kj}I_{2}+\omega_{\textup{nom}}L_{kj}\mathcal{J})^{-1}$, where $R_{kj}$ is the resistance and $L_{kj}$ is the inductance of the line $(k,j)$, for every $(k,j)\in\mathcal{E}$. Suppose $B$ is the incidence matrix of $G$. There are three types of nodes in the network: we have one grid bus with voltage $v_{\textup{gDQ}}=[0,V_{\textup{g}}]^{\top}$ denoted by [math]. We have $n\geq 1$ inverter buses collected in the set $\mathcal{N}_{\textup{I}}$, and $\ell$ load buses collected in the set $\mathcal{N}_{\textup{L}}$. Without loss of generality, we assume that $\mathcal{N}_{\textup{L}}=\{1,\ldots,\ell\}$ and $\mathcal{N}_{\textup{I}}=\{\ell+1,\ldots,n+\ell\}$ such that $\mathcal{N}=\{0\}\cup\mathcal{N}_{\textup{L}}\cup\mathcal{N}_{\textup{I}}$. Therefore $|\mathcal{N}|=1+\ell+n$ and we assume that $|\mathcal{E}|=m$. Associated to the matrix-weighted graph, $G$, we define the nodal admittance matrix by $Y=(B\otimes I_{2})\mathcal{A}(B\otimes I_{2})^{\top}\in\mathbb{R}^{2(1+\ell+n)\times 2(1+\ell+n)}$, where $\mathcal{A}\in\mathbb{R}^{2m\times 2m}$ is given by $\mathcal{A}=\mathrm{blkd}(a_{jk})$. The partition $\mathcal{N}=\{0\}\cup\mathcal{N}_{\textup{L}}\cup\mathcal{N}_{\textup{I}}$ induces the following decomposition of incidence matrix $B$: $B^{\top}=\begin{pmatrix}B^{\top}_{\textup{0}}&B^{\top}_{\textup{L}}&B^{\top}_{\textup{I}}\end{pmatrix}$, where $B_{\textup{0}}\in\mathbb{R}^{1\times m}$, $B_{\textup{L}}\in\mathbb{R}^{\ell\times m}$, and $B_{\textup{I}}\in\mathbb{R}^{n\times m}$, and the following partition for the admittance matrix $Y$:

[TABLE]

We also establish the following convention: for a given variable (parameter) $y$ corresponding to the inverter, we define vector $\mathbf{y}=(y_{1}^{\top},\ldots,y_{n}^{\top})^{\top}$, where $y_{k}$ is the associated variable (parameter) for the $k$th inverter.

Inverter model

Using (5), the governing dynamics for all inverters in the network can be expressed as:

[TABLE]

where, in $\mathbf{y}=(y_{1}^{\top},\ldots,y_{n}^{\top})^{\top}\in\mathbb{R}^{13n}$, $y_{k}$ captures all the dynamic states for the $k$th inverter, $F(\mathbf{y})=(f^{\top}_{1}(y_{1}),\ldots,f^{\top}_{n}(y_{n}))^{\top}$, $G(\mathbf{y})=\mathrm{diag}\left(g_{1}(y_{1}),\ldots,g_{n}(y_{n})\right)$, and $C=\mathrm{diag}\left(C_{1},\ldots,C_{n}\right)$, where $f_{k}$ is the drift vector field and $g_{k}$ and $C_{k}$ are control vector fields of inverter $k$.

Load and line models

Let $v_{k}\in\mathbb{R}^{2}$ be the $k$th load voltage (in the $\mathrm{DQ}$-frame) and $i_{k}\in\mathbb{R}^{2}$ be the current demand (in the $\mathrm{DQ}$-frame) of the $k$th load. We collect the nodal voltages and current demands for the loads in $\mathbf{v}_{\textup{L}}=(v_{1}^{\top},\ldots,v_{\ell}^{\top})^{\top}\in\mathbb{R}^{2\ell}$ and $\mathbf{i}_{\textup{L}}=(i_{1}^{\top},\ldots,i_{\ell}^{\top})^{\top}\in\mathbb{R}^{2\ell}$, respectively. We assume loads are purely resistive; suppose $\mathbf{R}_{\textup{L}}\in\mathbb{R}^{\ell}$ is the vector of load resistances, then the loads can be described by

[TABLE]

Suppose that the vector of line resistances and line inductances are denoted by $\mathbf{R}_{\mathcal{E}}\in\mathbb{R}^{m}$ and $\mathbf{L}_{\mathcal{E}}\in\mathbb{R}^{m}$, respectively; the nodal current injections by $\mathbf{i}=(\mathrm{i}_{\textup{g}},-\mathbf{i}_{\textup{L}},\mathbf{i}_{\textup{oDQ}})^{\top}$ and nodal voltages by $\mathbf{v}=(v_{\textup{gDQ}},\mathbf{\mathbf{v}}_{\textup{L}},\mathrm{\mathbf{v}}_{\textup{oDQ}})^{\top}$. The governing dynamics for the transmission lines are [26, Equation 4.10]:

[TABLE]

where $\xi_{\textup{DQ}}\in\mathbb{R}^{2m}$ is the vector of current flows in the lines. Thus we have $\mathbf{i}=(B\otimes I_{2})\xi_{\textup{DQ}}$.

III-A Network Model and Dimensionless Transcription

From (6)–(8), the grid-tied inverter-network dynamics are:

[TABLE]

where $\mathbf{v}=(v_{\textup{gDQ}},\mathbf{\mathbf{v}}_{\textup{L}},\mathrm{\mathbf{v}}_{\textup{oDQ}})^{\top}$ and \mathbf{v}_{\textup{L}}=\big{(}-B_{\textup{L}}[\mathbf{R}_{\textup{L}}]\otimes I_{2}\big{)}\xi_{\textup{DQ}}. We transcribe the differential equations (III-A) in a dimensionless format. We assume that $s_{\textup{nom}}$ is the nominal power generation/consumption in the network. For each inverter, we introduce the following dimensionless variables:

[TABLE]

For the network, we introduce the dimensionless parameters

[TABLE]

We also isolate time-constants of different sub-systems:

• •

for the PLL low-pass filter: $\tau_{\textup{PLL}}=\omega^{-1}_{\textup{c,PLL}}$ and $\tau^{\prime}_{\textup{PLL}}=(V_{\textup{g}}\mathrm{kp}_{\textup{PLL}})^{-1}$;

• •

for the PLL PI controller: $\mathrm{T}_{\textup{PLL}}=\frac{\mathrm{kp}_{\textup{PLL}}}{\mathrm{ki}_{\textup{PLL}}}$;

• •

for the low-pass filter of the power controller: $\tau_{\textup{s}}=\omega^{-1}_{\textup{s}}$;

• •

for tracking in the power controller: $\tau^{\prime}_{\textup{s}}=(V_{\textup{g}}\mathrm{ki}_{\textup{s}})^{-1}$;

• •

for the PI controller in the power controller: $\mathrm{T}_{\textup{s}}=\frac{\mathrm{kp}_{\textup{s}}}{\mathrm{ki}_{\textup{s}}}$;

• •

for the current controller: $\tau_{\textup{c}}=V_{\textup{g}}^{2}(\mathrm{ki}_{\textup{c}}s_{\textup{nom}})^{-1}$;

• •

for the PI controller in the current controller: $\mathrm{T}_{\textup{c}}=\frac{\mathrm{kp}_{\textup{c}}}{\mathrm{ki}_{\textup{c}}}$;

• •

for the LC filter: $\tau_{\textup{LC}}=\frac{\mathrm{L}_{\textup{f}}}{\sqrt{\omega^{-2}_{\textup{nom}}\mathrm{C}^{-2}_{\mathrm{f}}+\omega^{2}_{\textup{nom}}\mathrm{L}^{2}_{\textup{f}}}}$ and $\tau^{\prime}_{\textup{LC}}=\frac{\mathrm{C}_{\textup{f}}}{\sqrt{\omega^{2}_{\textup{nom}}\mathrm{C}^{2}_{\mathrm{f}}+\omega^{-2}_{\textup{nom}}\mathrm{L}^{-2}_{\textup{f}}}}$ and $\tau^{\prime\prime}_{\textup{LC}}=\mathrm{C}_{\textup{f}}\omega_{\textup{nom}}s^{-1}_{\textup{nom}}V_{\textup{g}}^{2}$;

• •

for line $e\in\mathcal{E}$ in the network: $\tau_{\textup{e}}=\frac{L_{e}}{\sqrt{R^{2}_{e}+\omega^{2}_{\textup{nom}}L_{e}^{2}}}$ and $\tau^{\prime}_{\textup{e}}=s_{\textup{nom}}V_{\textup{g}}^{-2}\sqrt{R^{2}_{e}+\omega^{2}_{\textup{nom}}L_{e}^{2}}$.

With these preliminaries in place, the dimensionless grid-tied inverter-network dynamics can be expressed as:

[TABLE]

Above,

[TABLE]

The dimensionless grid-tied inverter network dynamics (10a)–(10i) is $(13n+2m)$-dimensional. Our first goal is to find the equilibrium points of the dynamical systems (10a)–(10i).

III-B *Equilibrium points of the Dimensionless Grid-tied

Inverter-network Dynamics*

We start by introducing some notation. Let $Y\in\mathbb{R}^{(2n+2)}$ be the admittance matrix of the network. Then

[TABLE]

and the dimensionless parameters:

[TABLE]

We show that the equilibrium points of the dimensionless grid-tied inverter network dynamics (10) are in correspondence with the solutions of the following power-flow equations:

[TABLE]

Lemma 1**.**

For a given reference-power injection $\widehat{s}^{\mathrm{ref}}$ to the inverters, the following statements are equivalent:

1. (i)

$(\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}},\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}})^{\top}\in\mathbb{R}^{4n}$* is a solution for the power-flow equations (12) and  (13);* 2. (ii)

for every $\alpha=(\alpha_{1},\ldots,\alpha_{n})^{\top}\in\mathbb{Z}^{n}_{2}$, $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}\in\mathbb{R}^{13n+2m}$ is an equilibrium point of the dimensionless grid-tied inverter-network dynamics (10) given by

[TABLE]

with

[TABLE]

Proof.

Regarding $(ii)\Longrightarrow(i)$, from the power-controllers’ dynamics in (10), we can conclude that if $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$ is an equilibrium point, then we have $\widehat{\mathbf{s}}=\widehat{\mathbf{s}}^{\mathrm{ref}}$. This implies $(\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}},\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}})^{\top}$ satisfies (12). Moreover, at the equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$, the line dynamics (10i) will simplify to $\widehat{\mathbf{i}}^{\mathrm{ref}}=\widehat{Y}\widehat{\mathbf{v}}^{\mathrm{ref}}$. Using Kron reduction [8], this implies that $(\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}},\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}})^{\top}$ satisfies (13).

Regarding $(i)\Longrightarrow(ii)$, suppose $(\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}},\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}})^{\top}$ is a solution to the power flow equations (12) and (13). Then from the PLL dynamics in (10), we have $\widehat{\mathbf{v}}_{\mathrm{PLL}}=\mathbbold{0}_{n}$ and $\widehat{\mathbf{v}}_{\mathrm{od}}=\mathbbold{0}_{n}$. Note that, $\widehat{\mathbf{v}}_{\mathrm{od}}=\mathbbold{0}_{n}$ can be written in the trigonometric form

[TABLE]

This implies that, for every $k\in\{1,\ldots,n\}$, there exists $\alpha_{k}\in\mathbb{Z}_{2}$ such that $\delta^{k,{\mathrm{ref}}}=-\mathrm{arg}(v^{k,{\mathrm{ref}}}_{\mathrm{oDQ}})+\alpha_{k}\pi$. From the power-controller dynamics in (10), we have $\widehat{\mathbf{s}}_{\textup{avg}}=\widehat{\mathbf{s}}^{\mathrm{ref}}$. Finally, one can find the value of $\widehat{\boldsymbol{\phi}}^{\mathrm{ref}}_{\textup{s}}$, $\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{lDQ}}$, $\widehat{\boldsymbol{\gamma}}^{\mathrm{ref}}_{\textup{DQ}}$, and $\widehat{\xi}^{\mathrm{ref}}_{\textup{DQ}}$ by solving the remaining equations in (10). ∎

Remark 2**.**

The power-flow equations (12) and (13) have been studied extensively in the literature and many sufficient conditions for existence and uniqueness of solutions have been developed; see, e.g., [3, 7, 29]. Lemma 15 in Appendix C restates a result from [29] pertaining to uniqueness that is leveraged in subsequent results.

IV Stability Analysis of the Dimensionless Grid-tied Inverter-network Dynamics

In bulk power-systems dynamics literature, it is commonplace to assume that the dynamics of the grid-following inverters are much faster than the dynamics of grid-forming inverters and synchronous machines [26]. This assumption justifies the use of a static model for grid-following inverters such that the inverter nodes are considered to be sources of constant (active and reactive) power. Subsequently, the network operation is described by the following power-flow equations:

[TABLE]

Quite obviously, the static representation (14) does not capture stability. The following example shows that the internal dynamics of the inverters can induce instabilities, even if the power-flow equations in (14) admit a high-voltage solution.

Example 3**.**

(Instabilities Induced by Inverter Dynamics)* Consider the radial grid-connected network consisting of $25$ identical inverters (a sketch is provided in Fig. 2). Suppose the inverters have uniform reference power injections $\mathbf{p}^{\mathrm{ref}}=p\mathbbold{1}_{25}$, each line has resistance $R=10^{-2}\ \Omega$ and inductance $L=10^{-5}\ \mathrm{H}$, and the grid voltage is $v_{\textup{gDQ}}=[0,120\sqrt{2}]^{\top}\ \mathrm{V(peak)}$ with constant frequency $\omega_{\mathrm{nom}}=120\pi\ \mathrm{rad/s}$.*

1. (i)

Static model:* If inverter power injections satisfy*

[TABLE]

then, there exists a unique high-voltage, low-current solution for the power-flow equations (14) (see Lemma 15). 2. (ii)

Dynamic model:* We use the dynamic model (5) for the inverters (parameters are given in the fourth column of Table III). The governing equations for the network are in (10), and condition (15) via Lemma 1 guarantees the existence of a family of equilibrium points $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$, $\alpha\in\mathbb{Z}^{n}_{2}$ for the system (10). Linearizing the system (10), we study local stability of the equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\mathrm{0}}$.*

Figure 3 plots the maximum real part of the eigenvalues of the linearized system (10) around the equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\mathrm{0}}$ as a function of active-power injection. The red vertical line is the threshold of the power injection for which the power-flow equations admit a unique solution (obtained from (15)). Notice that there are power injections for which a high-voltage solution of the power-flow equations exists, however, the corresponding equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\mathrm{0}}$ is not stable.

IV-A Small-signal Stability via Time-scale Separation

We now focus on the small-signal stability of the dimensionless grid-tied inverter-network dynamics (10). Due to high dimensionality and nonlinearity of the dynamic model, studying small-signal stability is not analytically tractable. Therefore, it is a reasonable goal to reduce the model order. To that end, we identify a physically meaningful parametrization of the inverters, and we show that under suitable assumptions, this parametrization leads to a time-scale decomposition of the system (10) which simplifies analysis. We begin by uncovering time constants of different sub-systems.

Definition 4** (Singular perturbation parameter).**

For dynamics (10), we define:

[TABLE]

Using these, we define the singular perturbation parameter:

[TABLE]

Remark 5**.**

1. (i)

Definition 4 establishes a physically meaningful time-scale separation of the components of the inverter network when $\epsilon\ll 1$. In this case, $\|\boldsymbol{\tau}^{\prime}_{\mathrm{LC}}\|_{\infty},\|\boldsymbol{\tau}_{\mathrm{LC}}\|_{\infty},\|\boldsymbol{\tau}_{\mathcal{E}}\|_{\infty},\|\boldsymbol{\tau}_{\textup{c}}\|_{\infty},\|\mathbf{T}_{\textup{c}}\|_{\infty}\leq\epsilon^{2}$ which implies that the line dynamics, the $LC$ filter, and the current controller of the inverter are the fastest components of the network. Moreover, $\|\boldsymbol{\tau}_{\textup{PLL}}\|_{\infty},\|\boldsymbol{\tau}^{\prime}_{\textup{PLL}}\|_{\infty},\|\boldsymbol{\tau}_{\textup{s}}\|_{\infty}\leq\epsilon$ which implies that the PLL, and averaging-part of the power controller (i.e., $\widehat{\mathbf{s}}_{\mathrm{avg}}$) are slower than the current controller, the line dynamics, and the $LC$ filter but they are faster than the steady-state power-tracking controller (i.e., $\widehat{\boldsymbol{\phi}}_{\mathrm{s}}$) dynamics. 2. (ii)

The assumption $\epsilon\ll 1$ is equivalent to $\epsilon_{\textup{I}},\epsilon_{\textup{L}}\ll 1$, which is realistic in practice. For instance, in the inverter models with the parameters in Table III and with the parameters used in **[17, 23, 24]**, it holds that $\epsilon_{\textup{I}},\epsilon_{\textup{L}}\leq 0.1$. 3. (iii)

$\tau_{e}=\frac{L_{\textup{e}}}{\sqrt{R^{2}_{\textup{e}}+\omega^{2}_{\textup{nom}}L^{2}_{\textup{e}}}}\leq\omega^{-1}_{\textup{nom}}$, for every $e\in\mathcal{E}$. Similarly, one can show that $\|\boldsymbol{\tau}^{\prime}_{\textup{LC}}\|_{\infty},\|\boldsymbol{\tau}_{\textup{LC}}\|_{\infty}\leq\omega^{-\frac{1}{2}}_{\textup{nom}}$ and this gives an upper bound for the LC and line parameter $\epsilon_{\mathcal{E}}\leq\omega^{-\frac{1}{2}}_{\textup{nom}}$. 4. (iv)

For a different set of parameters and control architectures, one can conceivably identify a different singular perturbation parameter and time-scale decomposition. However, we expect the general nature of the stability result that follows to be similar.

Theorem 6** (Small-signal Stability).**

Consider the dimensionless grid-tied inverter-network dynamics (10) with states in $\mathbb{R}^{13n+2m}$ and references $\widehat{\mathbf{s}}^{\mathrm{ref}}\in\mathbb{R}^{2n}$. The following hold:

1. (i)

if,

[TABLE]

then, there is a unique solution $(\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}},\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}})^{\top}$ for the power-flow equations (12) and (13) and a family of equilibrium points $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$, $\alpha\in\mathbb{Z}_{2}^{n}$, for the grid-tied inverter network dynamics (10) satisfying:

[TABLE] 2. (ii)

if additionally, $[\mathbf{T}_{\textup{PLL}}]\succ[\boldsymbol{\tau}_{\textup{PLL}}]$ and the $2n\times 2n$ matrix

[TABLE]

is Hurwitz, then, there exists an $\epsilon^{*}>0$ such that, for every $\epsilon\leq\epsilon^{*}$, the equilibrium point $\widehat{\mathbf{x}}_{\mathrm{0}}^{\mathrm{ref}}$ is locally exponentially stable.

Proof.

Regarding part (i), the proof follows from combining Lemma 15 and Lemma 1. Regarding part (ii), consider the dimensionless grid-tied inverter network dynamics (10). Using (16) and defining the variable $\Delta\mathbf{x}\in\mathbb{R}^{13n+2m}$ by $\Delta\mathbf{x}=\mathbf{x}-\mathbf{x}^{\mathrm{ref}}_{0}$, we get a three-time-scale decomposition of the system. We define the states $\mathbf{z}_{1}\in\mathbb{R}^{6n+2m}$, $\mathbf{z}_{2}\in\mathbb{R}^{5n}$, and $\mathbf{z}_{3}\in\mathbb{R}^{2n}$ as follows:

[TABLE]

where $\mathbf{z}_{1}$ is faster than $\mathbf{z}_{2}$ and $\mathbf{z}_{2}$ are faster than $\mathbf{z}_{3}$. The corresponding time-scales are given by $\tau=t/\epsilon^{2}$ and $\tau^{\prime}=t/\epsilon$. Using these time scales, the dimensionless grid-tied inverter network dynamics (10) can be written as follows:

[TABLE]

where $g_{1},g_{2},g_{3}$ are suitably defined functions. The quasi-steady-state manifold for the time-scale $\tau$ is the manifold $\mathbf{z}_{1}=h_{1}(\mathbf{z}_{2},\mathbf{z}_{3})$ obtained by solving the algebraic equations $g_{3}(\mathbf{z}_{1},\mathbf{z}_{2},\mathbf{z}_{3},0)=\mathbbold{0}_{(6n+2m)}$ [12, §11.2]. After some algebraic computations, the quasi-steady-state manifold for the time-scale $\tau$ is time-invariant and is given by:

[TABLE]

Since the quasi-steady state manifold for the time-scale $\tau$ is an isolated manifold, the singular perturbation problem is in the standard form. The boundary-layer dynamics are:

[TABLE]

We first show that for the boundary-layer dynamics (20), the origin is the exponentially stable equilibrium point. It is easy to see that since the graph is connected, we have $\mathrm{Ker}(B_{\textup{I}}^{\top})=\{\mathbbold{0}_{n}\}$, the matrix $([\boldsymbol{\tau}^{\prime\prime}_{\textup{LC}}\otimes I_{2}])\mathcal{J}$ is skew-symmetric and the matrix $\mathcal{Z}_{\mathrm{L}}+\mathcal{Z}^{\top}_{\mathrm{L}}$ is negative definite. Hence, using Lemma 14, the origin is the locally exponentially stable point of the boundary-layer dynamics (20). Similarly, the quasi-steady-state manifold for time-scale $\tau^{\prime}$ is $\mathbf{z}_{2}=h_{2}(\mathbf{z}_{3})$ obtained by solving the algebraic equations $g_{2}(h_{1}(\mathbf{z}_{2},\mathbf{z}_{3}),\mathbf{z}_{2},\mathbf{z}_{3},0)=\mathbbold{0}_{5n}$. Thus, the quasi-steady-state manifold for time-scale $\tau^{\prime}$ is time-invariant and is given by:

[TABLE]

where $\boldsymbol{\delta}^{BL}=\boldsymbol{\delta}^{BL}(\Delta\widehat{\boldsymbol{\phi}}_{\textup{s}})$ is the solution to:

[TABLE]

Since the quasi-steady state manifold for the time-scale $\tau^{\prime}$ is an isolated manifold, the singular perturbation problem is in the standard form with boundary-layer dynamics:

[TABLE]

We show that the boundary-layer dynamics (21) are stable around the origin. Note that the linearized boundary-layer equation dynamics around the origin have the form $\Delta{\mathbf{\dot{z}}_{2}}=S\Delta\mathbf{z}_{2}$, where $S\in\mathbb{R}^{5n\times 5n}$ has the upper block triangular form $S=\begin{pmatrix}S_{11}&S_{12}\\ \mathbbold{0}_{2n\times 3n}&S_{22}\end{pmatrix}$ with:

[TABLE]

Since $[\mathbf{T}_{\textup{PLL}}]\succ[\boldsymbol{\tau}_{\textup{PLL}}]$, by Lemma 14, matrix $S$ is Hurwitz. Therefore, the origin is the locally exponentially stable point of the boundary-layer equations (21). Now, we consider the reduced-order dynamics. The reduced-order model for the multi-time-scale analysis is given by:

[TABLE]

In order to study the stability of the reduced-order model (22), we introduce the new variable $\boldsymbol{\eta}=\mathcal{R}(-\boldsymbol{\delta}^{BL})\mathcal{H}\Delta\widehat{\boldsymbol{\phi}}_{\textup{s}}$. First note that, for $\epsilon=0$, we have $\Delta\widehat{\boldsymbol{\phi}}_{\mathrm{s}}=\mathbbold{0}_{2n}$ and $\boldsymbol{\delta}^{BL}=\boldsymbol{\delta}^{\mathrm{ref}}$. This implies that $\dot{\boldsymbol{\eta}}=\mathcal{R}(-\boldsymbol{\delta}^{\mathrm{ref}})\mathcal{H}\Delta\dot{\widehat{\boldsymbol{\phi}}}_{\textup{s}}$ and the reduced order model (22) is given by:

[TABLE]

Linearizing the above equation, we get $\Delta\dot{\boldsymbol{\eta}}=M^{\prime}\Delta\boldsymbol{\eta}$, where

[TABLE]

Note that, for any scalar $c\in\mathbb{R}_{>0}$ and any matrix $A\in\mathbb{R}^{2n\times 2n}$, the matrix $A$ is Hurwitz if and only if the matrix $cA$ is Hurwitz. Moreover, by [11, Theorem 1.3.22], Hurwitzness of the matrix $M^{\prime}$ is equivalent to the Hurwitzness of matrix $M$ defined in (18). This means that $M^{\prime}$ is Hurwitz and the origin is a locally exponentially stable equilibrium point of the reduced-order model (22). Stability of the equilibrium point of the dimensionless grid-tied inverter network dynamics follows from condition (ii) [12, Theorem 11.4].∎

We provide a few contextualizing and clarifying remarks.

Remark 7**.**

1. (1)

It is well-known in singular perturbation analysis that the largest value of $\epsilon^{*}$ for which Theorem 6(ii)* holds is hard to compute. A standard lower bound on $\epsilon^{*}$ can be obtained using techniques in [13, Lemma 2.2].* 2. (2)

Small-signal stability analysis of the grid-tied inverter network dynamics (10) is computationally complicated for large networks and, in general, it requires linearizing the system (10) and checking stability of a $(13n+2m)$-dimensional matrix. Theorem 6(ii)* uses a time-scale analysis to eliminate the line dynamics, the dynamics of the current controller, PLL, and $LC$ filter from small-signal stability analysis. The sufficient condition in Theorem 6(ii)** significantly reduces computational complexity by reducing the problem to checking that a $2n$-dimensional matrix is Hurwitz;* 3. (3)

By part (i), there exists a family of equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$ for the grid-tied inverter network dynamics (10). However, in part (ii), we focus on the equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\mathrm{0}}$. The reason is that this trajectory is locally exponentially stable when there is no power injection for the inverters; 4. (4)

Theorem 6(ii)* brings out the role of network topology, the power injections/demands, and inverter parameters in small-signal stability of inverter networks. For the matrix $M$, the term $\widehat{Y}_{\textup{red}}$ reveals the role of network topology. The terms $\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{oDQ}}$, $\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oDQ}}$, and $\delta^{\mathrm{ref}}$ (obtained by solving the power-flow equations (12)–(13)) reveal the role of network topology and power injections/demands. Finally, $\boldsymbol{\tau}^{\prime}_{\textup{s}}$ reveals the role of inverters’ internal dynamics.*

IV-B Corollaries

We now present two corollaries that may be applicable in different settings and shed more light on to the main result.

Corollary 8** (Single inverter connected to the grid).**

Consider a single inverter with the reference power injection $s^{\mathrm{ref}}=(p^{\mathrm{ref}},q^{\mathrm{ref}})^{\top}$ connected to the grid with voltage $v_{\textup{gDQ}}=\left(\begin{smallmatrix}0\\ V_{\textup{g}}\end{smallmatrix}\right)$ through a line with resistance $R$ and inductance $L$. If

[TABLE]

then the following statements hold:

1. (i)

there exist two equilibrium points $\widehat{\mathrm{x}}^{\mathrm{ref}}_{0}$ and $\widehat{\mathrm{x}}^{\mathrm{ref}}_{1}$ satisfying

[TABLE] 2. (ii)

if we have $\mathrm{T}_{\textup{PLL}}>\mathrm{\tau}_{\textup{PLL}}$, then there exists a $\epsilon^{*}>0$ such that, for every $\epsilon\leq\epsilon^{*}$, $\widehat{\mathrm{x}}^{\mathrm{ref}}_{0}$ is locally exponentially stable.

Proof.

Regarding part (i), we define the dimensionless resistance $\widehat{R}=V_{\textup{g}}^{-2}s_{\textup{nom}}R$ and dimensionless inductance $\widehat{L}=V_{\textup{g}}^{-2}s_{\textup{nom}}L$. Therefore, we get $\widehat{Y}^{-1}_{\textup{red}}=\begin{pmatrix}\widehat{R}&-\omega_{\mathrm{nom}}\widehat{L}\\ \omega_{\mathrm{nom}}\widehat{L}&\widehat{R}\end{pmatrix}$. Then, (23) is equivalent to (17) and the result follows from Theorem 6(i). Regarding part (ii), note that we have $\widehat{\mathrm{v}}^{\mathrm{ref}}_{\textup{oDQ}}-(\widehat{R}I_{2}+\omega_{\textup{nom}}\widehat{L}\mathcal{J})\widehat{\mathrm{i}}^{\mathrm{ref}}_{\textup{oDQ}}=\left(\begin{smallmatrix}0\\ 1\end{smallmatrix}\right)$. Then, matrix $M$ in (18) becomes:

[TABLE]

Using simple algebraic manipulations, we get $\det(M)=(2\widehat{v}^{\mathrm{ref}}_{\textup{oQ}}-1)(\tau^{\prime}_{\textup{s}})^{-1}$ and $\mathrm{tr}(M)=-2\widehat{v}^{\mathrm{ref}}_{\textup{oq}}(\tau^{\prime}_{\textup{s}})^{-1}$. Since $\|\widehat{\mathrm{v}}^{\mathrm{ref}}_{\textup{oDQ}}-\left(\begin{smallmatrix}0\\ 1\end{smallmatrix}\right)\|_{2}\leq\tfrac{1}{2}$, we have $2\widehat{v}^{\mathrm{ref}}_{\textup{oQ}}\geq 1$. This implies that $\det(M)\geq 0$. Moreover, for the equilibrium points $\widehat{\mathrm{x}}^{\mathrm{ref}}_{0}$, we have $\widehat{v}^{\mathrm{ref}}_{\textup{oq}}>0$. This implies that $\mathrm{tr}(M)<0$ and $M$ is Hurwitz. Thus, by Theorem 6(ii), $\widehat{\mathrm{x}}^{\mathrm{ref}}_{0}$ is locally exponentially stable. ∎

In the next corollary, we study the special case of resistive networks and provide a computationally efficient numerical method for checking sufficient condition (18).

Corollary 9** (Resistive network of inverters).**

Consider the dynamics (10). Suppose the lines are purely resistive and the reference power injections and demands are purely active. The following statements hold:

1. (i)

if $\left\|[\widehat{\mathbf{u}}]\widehat{L}^{-1}_{\textup{red}}[\widehat{\mathbf{u}}]^{-1}[\widehat{\mathbf{p}}^{\mathrm{ref}}]\right\|_{\mathbb{C},\infty}\leq\tfrac{3}{8}$, where $\widehat{\mathbf{u}}=\widehat{L}^{-1}_{\textup{red}}\widehat{L}_{\textup{0g}}v_{\textup{gDQ}}$, then there exists a family of equilibrium points $\widehat{\mathbf{x}}^{\mathrm{ref}}_{\alpha}$ for (10) with the property that

[TABLE] 2. (ii)

if additionally $[\mathbf{T}_{\textup{PLL}}]\succ[\boldsymbol{\tau}_{\textup{PLL}}]$ and the Metzler matrix

[TABLE]

is Hurwitz, then there exists $\epsilon^{*}>0$ such that for every $\epsilon\leq\epsilon^{*}$, equilibrium $\widehat{\mathbf{x}}_{\mathrm{0}}^{\mathrm{ref}}$ is locally exponentially stable.

Proof.

Regarding part (i), since the network is resistive and power injections/demands are purely active, we know that

[TABLE]

Therefore: $\left\|\mathcal{D}^{\prime}(\widehat{\mathbf{w}})\widehat{Y}^{-1}_{\textup{red}}(\mathcal{D}^{\prime}(\widehat{\mathbf{w}}))^{-1}\mathrm{D^{\prime}}({\mathbf{s}}^{\mathrm{ref}})\right\|_{\mathbb{C},\infty}=\left\|[\widehat{\mathbf{u}}]\widehat{L}^{-1}_{\textup{red}}[\widehat{\mathbf{u}}]^{-1}[\widehat{\mathbf{p}}^{*}]\right\|_{\mathbb{C},\infty}$. The result then follows from Theorem 6(i). Regarding part (ii), since there are no reactive-power injections from the inverters, we have $\widehat{\mathbf{i}}^{\mathrm{ref}}_{\textup{od}}=\mathbbold{0}_{n}$. Since the input voltage for the PLL and power controller is the output voltage of the $LC$ filter (i.e., $\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{odq}}$), and the network is purely resistive, we have $\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oD}}=\mathbbold{0}_{n}$ and as a result $\boldsymbol{\delta}^{\mathrm{ref}}=\mathbbold{0}_{n}$. Alternatively, this observation can be proved rigorously as follows. Since the power injections are purely active, the power injection vector $\mathbf{s}^{\mathrm{ref}}$ has the form $\mathbf{s}^{\mathrm{ref}}=\mathbf{p}^{\mathrm{ref}}\otimes\left(\begin{smallmatrix}1\\ 0\end{smallmatrix}\right)$. Therefore, starting from the initial condition $\mathbf{v}^{(0)}=\widehat{\mathbf{w}}=\widehat{\mathbf{u}}\otimes\left(\begin{smallmatrix}0\\ 1\end{smallmatrix}\right)$, the $k$th iteration in Lemma 15(ii) has the form $\mathbf{v}^{(k)}=\mathbf{u}^{(k)}\otimes\left(\begin{smallmatrix}0\\ 1\end{smallmatrix}\right)$, for every integer $k\in\mathbb{Z}_{\geq 0}$. This implies that, in the limit, we have $\widehat{\mathbf{v}}^{\mathrm{ref}}_{\textup{oD}}=\mathbbold{0}_{n}$ and as a result $\boldsymbol{\delta}^{\mathrm{ref}}=\mathbbold{0}_{n}$. Therefore, the matrix $M$ in condition (18) simplifies as shown below:

[TABLE]

Using Lemma 13, matrix $M$ is Hurwitz if and only if matrices

[TABLE]

are Hurwitz. Note that the active-power injections from the inverters to the grid are non-negative in steady state. This implies that $\widehat{\mathbf{i}}^{\mathrm{ref}}_{\mathrm{oq}}\geq\mathbbold{0}_{n}$. Thus, the matrices (25) are similar to the following matrices:

[TABLE]

Since the matrix $[\widehat{\mathbf{v}}^{\mathrm{ref}}_{\mathrm{oq}}]$ is positive definite and the matrix $[\widehat{\mathbf{i}}^{\mathrm{ref}}_{\mathrm{oq}}]^{\frac{1}{2}}\widehat{L}^{-1}_{\textup{red}}[\widehat{\mathbf{i}}^{\mathrm{ref}}_{\mathrm{oq}}]^{\frac{1}{2}}$ is positive semidefinite, the matrix

[TABLE]

is Hurwitz. Moreover, note that the matrix $L_{\textup{red}}$ is a grounded Laplacian matrix and by [4, E 9.10] its inverse $L_{\textup{red}}^{-1}$ is non-negative. Also, the matrices $[\boldsymbol{\tau}^{\prime}_{\textup{s}}]$, $[\widehat{\mathbf{i}}^{\mathrm{ref}}_{\mathrm{oq}}]$, and $[\widehat{\mathbf{v}}^{\mathrm{ref}}_{\mathrm{oq}}]$ are all diagonal with non-negative diagonal elements. This implies that the matrix $-[\widehat{\mathbf{v}}^{\mathrm{ref}}_{\mathrm{oq}}][\boldsymbol{\tau}^{\prime}_{\textup{s}}]^{-1}+[\widehat{\mathbf{i}}^{\mathrm{ref}}_{\mathrm{oq}}]\widehat{L}^{-1}_{\textup{red}}[\boldsymbol{\tau}^{\prime}_{\textup{s}}]^{-1}$ has non-negative off-diagonal elements and therefore, it is Metzler. The proof of Corollary (9) then follows from Theorem 6. ∎

Remark 10** (Computational complexity).**

There are computationally efficient methods for checking Hurwitzness of Metzler matrices. In particular, one can reformulate the Hurwitzness of (24) as the following feasibility problem:

[TABLE]

The feasibility problem (26) is a linear program and can be checked using distributed methods whose computational time scales linearly with the number of non-zero elements in $N$ [22].

V Numerical Simulations

In this section, we present numerical simulation results for radial networks (see Fig. 2) with identical inverters and reference-power setpoints with the goal of answering the following questions:111All the numerical simulations are performed in MATLAB R2016a on a computer with Intel Core i5 processor $@\ 1.6\ \mathrm{GHZ}$ CPU and $4\ \mathrm{GB}$ RAM.

• •

For which values of the inverter parameter $\epsilon_{\textup{I}}$ does Theorem 6(ii) guarantee small-signal stability?

• •

How efficient is the condition (18) in Theorem 6(ii) to analyze small-signal stability?

The first question can be interpreted as an inverter design problem, while the second question can be interpreted as a network monitoring problem. Solutions to these problems are provided in Sections V-A and V-B, respectively. Our test case is a family of radial networks $\{G(n)\}_{n\in\mathbb{N}}$, where $G(n)$ is the weighted undirected connected graph with the node set (buses) $\mathcal{N}_{n}=\{0,\ldots,n\}$ and the edge set (branches) $\mathcal{E}_{n}=\{(i,i+1)\mid i=0,\ldots,n-1\}$ (see Fig. 2). For every radial network $G(n)$, the node [math] is the slack bus connected to the grid with voltage $v_{\textup{g}}=\mathrm{i}(120\sqrt{2})\ \mathrm{V(peak)}$ and frequency $\omega_{\textup{nom}}=120\pi\ \mathrm{rad/s}$. Nodes $\{1,\ldots,n\}$ are the inverters with $s_{\textup{nom}}=1000\ \mathrm{VA}$, and

[TABLE]

For every $(j,k)\in\mathcal{E}_{n}$, the admittance of line $(i,j)$ is given by $a_{jk}=a_{kj}=(RI_{2}+\omega_{\mathrm{nom}}L\mathcal{J})^{-1}$, where the line resistance is $R=0.02\ \Omega$ and line inductance is $L=2\times 10^{-5}\ \mathrm{H}$. We first define the notion of a safe penetration level.

Definition 11** (Safe Penetration Level).**

Given a family of networks $\{G(n)\}_{n\in\mathbb{N}}$ and uniform reference power injection $s^{\mathrm{ref}}\in\mathbb{R}^{2}$, the Safe Penetration Level (SPL) for $\{G(n)\}_{n\in\mathbb{N}}$ is the maximum $n\in\mathbb{N}$ such that the dimensionless grid-tied inverter network dynamics (10) with underlying graph $G(n)$ and reference powers $\mathbf{s}^{\mathrm{ref}}=s^{\mathrm{ref}}\otimes\mathbbold{1}_{n}$ has a locally stable equilibrium point.

Remark 12**.**

1. (1)

For our test case network, SPL of the network depends on the parameter $\epsilon_{\textup{I}}$; 2. (2)

One can use the matrix $M$ defined in (18) to estimate the SPL of the network; from Theorem 6(ii)**, there exists $\epsilon^{*}>0$ such that, for every $\epsilon\leq\epsilon^{*}$, SPL of the network is larger than this estimate.

V-A Designing Grid-tied Inverter Networks

In this part, we examine the efficiency of our analytic results in Theorem 6 to design grid-following inverter networks focusing on the small-signal stability of the grid. In particular, we focus on estimating the SPL for the radial network shown in Fig 2 and numerically computing the largest range of parameter $\epsilon_{\textup{I}}$ for which Theorem 6(ii) holds. In order to carry out this task, we study the effect of parameter $\epsilon_{\textup{I}}$ for different active power injections on the SPL of the system. The result is shown in Fig. 4.

Note that the overlapping dashed lines in Fig. 4 are the estimates of SPL based on the sufficient condition in Theorem 6(ii). Therefore, according to the data in Fig. 4, the largest domain of the inverter parameter $\epsilon_{\textup{I}}$ for which Theorem 6(ii) holds is $(0,0.0025]$, for the active power injections $p^{\mathrm{ref}}=1000\ W$, $p^{\mathrm{ref}}=1500\ W$, and $p^{\mathrm{ref}}=2000\ W$, respectively. It is interesting that, as $\epsilon_{\textup{I}}$ becomes smaller, the dashed line and the solid lines get closer to each other and Theorem 6(ii) can be used to find the exact SPL of the network.

V-B Monitoring Grid-tied Inverter Networks

In this section, we study the accuracy and computational efficiency of condition (18) in Theorem 6(ii) for monitoring the small-signal stability of the network. For our test case, we pick the parameter $\epsilon_{\textup{I}}=0.001$. Therefore based on the discussion in Section V-A, we are in the range of applicability of Theorem 6(ii). Recall that SPL is the largest number of inverters in the network for which the equilibrium point $\widehat{\mathbf{x}}^{\mathrm{ref}}_{0}$ of full-order system (10) is locally asymptotically stable. We denote the computational time of finding SPL, using the eigenvalue analysis for the linearized system, by $\mathrm{T}_{\textup{lin}}$. We denote the largest number of inverters in the grid for which the matrix $M$ in (18) is Hurwitz by $\mathrm{SPL}_{\mathrm{test}}$. Similarly, we denote the computational time for checking the Hurwitzness of matrix $M$ in (18) by $\mathrm{T}_{\textup{test}}$. Finally, we denote the largest number of inverters in the grid for which condition (15) holds, by $\mathrm{SPL}_{\mathrm{static}}$. We start with different uniform reference active-power injections $\widehat{\mathrm{s}}^{\mathrm{ref}}=(\widehat{p},0)^{\top}\otimes\mathbbold{1}_{n}$ and compute the thresholds $\mathrm{SPL}$, $\mathrm{SPL}_{\mathrm{test}}$, and $\mathrm{SPL}_{\mathrm{static}}$ together with the computational times for $\mathrm{SPL}$ and $\mathrm{SPL}_{\mathrm{test}}$, for each active-power injections. The results are shown in Table I. From Table I, one can see that the static condition (15) overestimates the value of $\mathrm{SPL}$. Moreover, it is clear that condition (18) gives an accurate lower bound for safe penetration level and its corresponding computation time (i.e., $\mathrm{T}_{\textup{test}}$) is almost one order of magnitude less that the computation time to perform eigenvalue analysis for the full-order system (i.e., $\mathrm{T}_{\textup{lin}}$).

VI Conclusion

We studied small-signal stability of grid-tied networks of grid-following inverters and loads. Using a time-scale analysis and a suitable choice of a family of parameters for the inverters, we presented an analytic sufficient condition for local exponential stability. We showed that, compared to the direct eigenvalue analysis of the full-order system, this sufficient condition has the advantages of reducing the computational complexity of checking small-signal stability as well as providing insights about the role of the network topology and inverter parameters on stability.

Appendix A Table of variables and parameters

Table (II) collects the variables and their symbols for grid-following inverter model and Table (III) collects the parameter values for this class of inverters.

Appendix B Pertinent Results from Linear Algebra

Lemma 13**.**

Let $\eta_{1},\ldots,\eta_{m}$ be the eigenvalues of a matrix $C\in\mathbb{C}^{m\times m}$. For $A,B\in\mathbb{C}^{n\times n}$, $\lambda$ is an eigenvalue of $A\otimes I_{m}+B\otimes C$ if and only if it is an eigenvalue of $A+\eta_{k}B$, for some $k\in\{1,\ldots,m\}$.

We omit the proof of this elementary result.

Lemma 14**.**

Let $\Gamma,\Pi,\Xi,\Upsilon,\Sigma,\in\mathbb{R}^{n\times n}$ and $\Theta\in\mathbb{R}^{m\times m}$ be diagonal matrices with positive diagonal entries such that $\Gamma\succ\Sigma$. Suppose that $K\in\mathbb{R}^{m\times n}$ is an arbitrary matrix with $\mathrm{Ker}(K)=\{\mathbbold{0}_{n}\}$, $P\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix, and $Z\in\mathbb{R}^{m\times m}$ is such that $Z+Z^{\top}$ is negative definite. Then the following matrices are Hurwitz:

[TABLE]

Proof.

The characteristic polynomial of matrix $A$ is: $\lambda^{3}I_{n}+\lambda^{2}\Gamma+\lambda\Upsilon\Xi+\Upsilon\Xi\Sigma=\mathbbold{0}_{n}$. Since all matrices $\Gamma$, $\Sigma$, $\Xi$, and $\Upsilon$ are diagonal, $\lambda_{i}$ is an eigenvalue of $A$ if and only if $\lambda_{i}^{3}I_{n}+\lambda_{i}^{2}(\Gamma)_{i}+\lambda_{i}(\Upsilon)_{i}(\Xi)_{i}+(\Upsilon)_{i}(\Xi)_{i}(\Sigma)_{i}=0$. The diagonal elements of the matrices $\Gamma$, $\Sigma$, $\Xi$, and $\Upsilon$ are positive. Therefore, using the Routh–Hurwitz criteria, $\lambda_{i}\in\mathbb{C}_{-}$ if and only if $(\Gamma)_{i}(\Upsilon)_{i}(\Xi)_{i}>(\Sigma)_{i}(\Upsilon)_{i}(\Xi)_{i}$. Thus, $A$ is Hurwitz if and only if $\Gamma\succ\Sigma$. To show that the matrix $B$ is Hurwitz, we use LaSalle’s invariance principle [12, Theorem 4.4]. Consider the dynamical system $\dot{\mathrm{x}}=B\mathrm{x}$, where $\mathrm{x}=(\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3},\mathrm{x}_{4})^{\top}\in\mathbb{R}^{(6n+m)}$. We define the Lyapunov function $V:\mathbb{R}^{(6n+m)}\to\mathbb{R}$ by: V(\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3},\mathrm{x}_{4})=\tfrac{1}{2}\big{(}\mathrm{x}_{1}^{\top}\Gamma^{-1}\mathrm{x}_{1}+\mathrm{x}_{2}^{\top}\Xi^{-1}\mathrm{x}_{2}+\mathrm{x}_{3}^{\top}\Pi^{-1}\mathrm{x}_{3}+\mathrm{x}_{4}^{\top}\Theta^{-1}\mathrm{x}_{4}\big{)} Then, it is easy to check that $\dot{V}(\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3})=-\mathrm{x}_{2}^{\top}\Upsilon\Xi^{-1}\mathrm{x}_{2}+\frac{1}{2}\mathrm{x}_{4}^{\top}(Z+Z^{\top})\mathrm{x}_{4}$. Therefore, by LaSalle’s invariance principle, the trajectories of the system $\dot{\mathrm{x}}=B\mathrm{x}$ converges to the largest invariant set inside $S=\{\mathrm{x}\in\mathbb{R}^{(6n+m)}\;|\;\dot{V}(\mathrm{x})=0\}$. It is easy to see that $S=\{\mathrm{x}\in\mathbb{R}^{6n}\;|\;\mathrm{x}_{4}=\mathrm{x}_{2}=\mathbbold{0}_{n}\}$. Let us denote the largest invariant set inside $S$ by $L$. Our goal is to show that $L=\{\mathbbold{0}_{(6n+m)}\}$. Suppose that $\gamma:t\mapsto(\gamma_{1}(t),\gamma_{2}(t),\gamma_{3}(t),\gamma_{4}(t))$ is a trajectory which belongs identically to $S$. Then we have $\gamma_{4}(t)=\gamma_{2}(t)=\mathbbold{0}_{n}$. First note that $\dot{\gamma}_{4}(t)=0$ implies that $\Theta K\gamma_{3}(t)+\Theta Z\gamma_{4}(t)=\Theta K\gamma_{3}(t)=\mathbbold{0}_{m}$. Since $\mathrm{Ker}(K)=\{\mathbbold{0}_{m}\}$, we deduce that $\gamma_{3}(t)=0$. Moreover, we see that $\dot{\gamma}_{2}(t)=\mathbbold{0}_{n}\ \Longrightarrow\ \Xi\gamma_{1}(t)=\mathbbold{0}_{n}$. This implies that $\gamma_{1}(t)=\gamma_{2}(t)=\gamma_{3}(t)=\gamma_{4}(t)=0$. Thus, the only invariant set inside $S$ is $\{\mathbbold{0}_{(6n+m)}\}$ and thus $B$ is Hurwitz. ∎

Appendix C Solutions to power flow equation

Lemma 15**.**

Consider (12) and (13) with $\widehat{\mathbf{w}}=-\widehat{Y}^{-1}_{\textup{red}}\widehat{Y}_{\textup{0g}}$. Suppose

[TABLE]

Then the following statements hold:

1. (i)

the power flow equations (12) and (13) has a unique solution $(\widehat{\mathbf{v}}^{*}_{\mathrm{oDQ}},\widehat{\mathbf{i}}^{*}_{\mathrm{oDQ}})$ with $\widehat{\mathbf{v}}^{*}_{\mathrm{oDQ}}\in\Omega$, where

[TABLE] 2. (ii)

for every $\mathbf{v}^{0}\in\Omega$, the iteration procedure

[TABLE]

converges to $\widehat{\mathbf{v}}^{*}_{\mathrm{oDQ}}$, where $(\widehat{\mathbf{v}}^{*}_{\mathrm{oDQ}},\widehat{Y}_{\textup{red}}(\widehat{\mathbf{v}}^{*}_{\mathrm{oDQ}}-\widehat{\mathbf{w}}))$ is the unique solution to (12) and (13).

Proof.

By considering $\mathbb{R}^{2n}\simeq\mathbb{C}^{n}$, part (i) and (ii) are straightforward generalizations of [29, Theorem 1]. One should note the fact that the nodal variables in [29, Theorem 1] are average power injections/demands and therefore the power flow equations have the form $S=V\overline{I}$. However, in this paper, the nodal variables are instantaneous power injections/demands and the power flow equations read $s=\frac{3}{2}\mathcal{D}(v_{\mathrm{odq}})i_{\mathrm{odq}}$. ∎