Stability Analysis of Droop-Controlled Inverter-Based Power Grids via Timescale Separation
Stefanos Baros Christoforos N. Hadjicostis Francis O'Sullivan
eess.SY cs.SY

	
	
	
	
	
	
	
	XBlock-EOS: Extracting and Exploring Blockchain Data From EOSIO
    Stefanos Baros, Christoforos N. Hadjicostis, and  Francis O’Sullivan
S. Baros is with the National Renewable Energy Laboratory, Golden, CO 80401, USA,
stefanos.baros@nrel.gov. This work was conducted while he was a postdoctoral researcher at MIT Energy Initiative, Cambridge, MA.
Christoforos Hadjicostis is with the ECE Department of University of Cyprus, Nicosia, Cyprus, chadjic@ucy.ac.cy.
Francis O' Sullivan was with the MIT Energy Initiative and is now with Orsted Onshore North America, frankie@mit.edu.
    Received: date / Accepted: date
==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

empty
	empty
	
	
	
	
		We consider the problem of stability analysis for distribution grids with droop-controlled inverters and dynamic distribution power lines. The inverters are modeled as voltage sources with controllable frequency and amplitude. This problem is very challenging for large networks as numerical simulations and detailed eigenvalue analysis are impactical. Motivated by the above limitations, we present in this paper a systematic and computationally efficient framework for stability analysis of inverter-based distribution grids. To design our framework, we use tools from singular perturbation and Lyapunov theories.  Interestingly, we show that stability of the fast dynamics of the power grid depends only on the voltage droop gains of the inverters while,  stability of the slow dynamics, depends on both voltage and frequency droop gains. Finally, by leveraging these timescale separation properties, we derive sufficient conditions on the frequency and voltage droop gains of the inverters that warrant stability of the full system. We illustrate our theoretical results through a numerical example on the IEEE 13-bus distribution grid.
		
		
		
	
	
	
	
	
	
	

§ INTRODUCTION

	Enviromental and economic reasons together with recent technological advances are the primary drivers for high penetration of renewable energy resources (RERs) in  power systems <cit.>, <cit.>, <cit.>, <cit.>. Today, renewable energy is mainly generated in distribution grids closer to the end-users by small distributed energy resources (DERs) that are interfaced with the main grid through alternating current (AC) inverters <cit.>.  When the penetration levels are relatively small, the power generation of RERs can be often treated as negative demand and the control objective for RERs is maximum power output generation <cit.>.  In such scenarios, the employable control techniques for AC inverters are limited to maximum power point tracking strategies <cit.>. However, as the penetration of RERs around the world increases, these control techniques are not offered anymore as viable options <cit.>, <cit.>.  It is imperative that RERs use control methods for their inverters that enable them to attain better regulation of their power outputs <cit.>, <cit.>, <cit.>. One such control method that emerged in the early `90s is droop control <cit.>. Droop control can allow inverter-interfaced RERs to satisfy a given load demand while achieving a certain power sharing distribution in steady-state.

	The stability analysis problem for droop-controlled inverter-based grids was first investigated in <cit.>, <cit.> and <cit.> using detailed small-signal analysis. This kind of analysis usually involves the repeated computation of the system's eigenvalues and extensive numerical simulations. In the same spirit,  stability of droop-controlled inverter-based power grids  was studied more recently in <cit.>, <cit.> and <cit.> via sensitivity and eigenvalue analysis. These approaches are easily implementable but they can be computationally very costly  as the network size increases. This is widely recognized by the scientific community and recently various researchers focused on examining stability of droop-controlled inverter-based grids via Lyapunov-like approaches. Representative publications along this line of work are <cit.>, <cit.> and <cit.>. In <cit.>, an approach for obtaining a reduced-order model for a droop-controlled inverter-based microgrid with electromagnetic network dynamics was first presented. A Lyapunov function for the reduced-order system was then constructed that led to decentralized sufficient stability conditions.   In <cit.>, sufficient conditions for stability of meshed microgrids with droop-controlled inverters were derived. As shown and explained in <cit.> however, it is absolutely necessary to consider the network dynamics in the stability analysis of droop-controlled inverter-based power grids as these may greatly impact the obtained stability regions. Finally, in <cit.>, stability of inverter-based power grids with grid-forming virtual oscillator control and line dynamics was studied, using a Lyapunov-like approach and singular perturbation theory. 
	Contributions. In this paper, we study the stability problem for distribution grids composed of constant impedance loads, droop-controlled inverters and dynamic distribution lines. We model inverters as voltage sources with controllable amplitude and frequency. Contrary to <cit.>, we do not require constant voltage amplitudes and incorporate line dynamics in our analysis.  Our main contribution is two-fold. We first present a systematic framework for conducting stability analysis of inverter-based power grids in a compositional and computationally efficient fashion. Subsequently, we apply our proposed framework and derive simple sufficient stability conditions.  Contrary to <cit.>, when our derived conditions hold they result in guaranteed stability of the original full system and not only of the low-order system approximation. 
	The remainder of the paper is structured as follows. In Section II, we review the dynamical model of a droop-controlled inverter-based power grid with line dynamics. In Section III, we formulate the stability problem. In Section IV we present our main results. In Section V, we verify our results through a numerical example. Finally in Section VI, we conclude this paper.
	

§ REVIEW OF INVERTER-BASED POWER GRID MODEL

	We consider a distribution grid comprised of N droop-controlled inverters, M transmission lines and L loads, described by the sets 𝒩, ℰ and ℒ, respectively.
	

 §.§ Droop-controlled Inverter-based Power Grid

	We model inverters as AC voltage sources with controllable amplitude and frequency, lines as “R-L"  circuits, and loads as constant impedances. Without loss of generality, we assume that loads exist only on inverter buses <cit.>, <cit.>.
	The state-variables associated with the inverters are the angles, frequencies and amplitudes of the output voltages denoted by, δ:=(δ_1,...,δ_N)^T∈ℝ^N, ω:=(ω_1,...,ω_N)^T∈ℝ^N and V:=(V_1,...,V_N)^T∈ℝ^N, respectively. The state-variables associated with the transmission lines are the currents I_D:=(I_D,1,...,I_D,M)^T∈ℝ^M and I_Q:=(I_Q,1,...,I_Q,M)^T∈ℝ^M, expressed in a dq reference frame. We note that, δ∈ℝ^N are measured in radians, ω∈ℝ^N in (r/s) while V∈ℝ^N, I_D∈ℝ^M, I_Q∈ℝ^M and P, Q∈ℝ^N are measured in per unit values. The matrices R=({R_i}_i=1^M), X=({X_i}_i=1^M) and L=({L_i}_i=1^M) are all diagonal matrices with the per unit resistances, reactances, and inductances on their diagonals, respectively.  The constant vectors ω_d:=(ω_d,1,...,ω_d,N)^T∈ℝ^N, V_d:=(V_d,1,...,V_d,N)^T∈ℝ^N, P_d:=(P_d,1,...,P_d,N)^T∈ℝ^N and Q_d:=(Q_d,1,...,Q_d,N)^T∈ℝ^N collect respectively, the desired frequencies, voltage amplitudes and the real and reactive power outputs of the inverters.  V_b and S_b  are the base voltage and base power, respectively. Finally, K_P:=K_P^puω_b=(N_PS_b), K_P:=({k_P,i}_i=1^N) are the scaled frequency droop gains that map per unit power output changes to frequency changes in (r/s) and K_Q:=K_Q^pu=K_Q(S_b/V_b),  K_Q:=({k_Q,i}_i=1^N) are the per unit voltage droop gains.   
	Given the above notation, the model of a distribution grid can be expressed compactly in state-space form <cit.>, <cit.>

    dδ/dt   =(ω-ω_b1_N)
    dω/dt   =T_P^-1[-ω+ω_d-K_P(P-P_d)]
    dV/dt   =T_Q^-1[-V+V_d-K_Q(Q-Q_d)]
    dI_D/dt   =ω_bL^-1[-RI_D+XI_Q+C^𝒯V_D] 
    dI_Q/dt   =ω_bL^-1[-RI_Q-XI_D+C^𝒯V_Q]

 The d-q components of the extended voltage vector V_D, V_Q are given by
	
    V_D=((({cos(δ_i)}_i=1^N)V)^T,     V_g,D)^T∈ℝ^N+1
    V_Q=((({sin(δ_i)}_i=1^N)V)^T,     V_g,Q)^T∈ℝ^N+1

	where V_g,D, V_g,Q∈ℝ are the d-q components of the feeder's voltage and ω_b is the grid's nominal frequency 2π60  (r/s). T_P=({T_P,i}_i=1^N) and T_Q=({T_Q,i}_i=1^N) are the time constants of the low-pass filters.   The real power of the inverters P∈ℝ^N can be expressed as
	
    P   =({cos(δ_i)}_i=1^N)({V_i}_i=1^N)[C^ℐI_D+C^ℒ𝒟V_DQ]
       +({sin(δ_i)}_i=1^N)({V_i}_i=1^N)[C^ℐI_Q+C^ℒ𝒬V_DQ]

	while the reactive power output Q∈ℝ^N as
	
    Q   =({sin(δ_i)}_i=1^N)({V_i}_i=1^N)[C^ℐI_D+C^ℒ𝒟V_DQ]
       -({cos(δ_i)}_i=1^N)({V_i}_i=1^N)[C^ℐI_Q+C^ℒ𝒬V_DQ]

 The vector V_DQ is defined as
	
    V_DQ=((({cos(δ_i)}_i=1^N)V)^T, (({sin(δ_i)}_i=1^N)V)^T )^T

	The load impedance matrices C^ℒ𝒟:={c^ℒ𝒟_ij}∈ℝ^N× 2N and C^ℒ𝒬:={c^ℒ𝒬_ij}∈ℝ^N× 2N are defined as
	
    c^ℒ𝒟_ij=R_L,k/R_L,k^2+X_L,k^2,    if load  k  lies at inv. bus  i  and  j=i 
    X_L,k/R_L,k^2+X_L,k^2,  if load  k  lies at inv. bus  i  and  j=i+N
    
    	0,  otherwise

	
    c^ℒ𝒬_ij=-X_L,k/R_L,k^2+X_L,k^2,  if load  k  lies at inv. bus  i  and  j=i  
    R_L,k/R_L,k^2+X_L,k^2,     if load  k  lies at inv. bus  i  and  j=i+N 
    
    	0,  otherwise

	Above, R_L,k and X_L,k represent the resistance and reactance respectively of load k, where k∈ℒ. Further,  C^ℐ:={c^ℐ_ij}∈ℝ^N× M is the incidence matrix, the entries of which are defined as
	
    c^ℐ_ij=-1,  if inverter  i  lies at the end of line  j
    
    	1,  if inverter  i  lies at the beginning of line  j
    
    	0,  otherwise

	In addition, C^𝒯:={c^𝒯_ij}∈ℝ^M× (N+1) is another incidence matrix, the entries of which can be defined as
	
    c^𝒯_ij=-1,  if inverter  j  lies at the end of line  i
    
    	1,    if inverter  j  lies at the beginning of line  i
    
    	0,             otherwise

	In defining C^ℐ, C^𝒯 we used the following convention. For a transmission line j:=(a,b)∈ℰ, we consider a to be its beginning and b its end when the current is flowing from node a to node b. 

	

§ PROBLEM FORMULATION

 Despite the wide use of droop-based control methods, today there is still a lack of systematic methodologies for stability analysis of power grids that accommodate numerous droop-controlled inverters. It is also not clear how we should best approach the stability analysis problem for such systems, in order to characterize the ranges for the droop control gains of the inverters that lead to guaranteed stability <cit.>, <cit.>. In light of this, the main research problem still remains open and can be formulated as follows.
	
		Consider the system (<ref>)–(<ref>) linearized around an equilibrium point w_0
		
    dw(t)/dt=Aw(t),  w(t)∈ℝ^3N+2M

		where w=[δ^⊤, ω^⊤, V^⊤,  I_D^⊤,  I_Q^⊤]^⊤ and f is the vector field of  (<ref>)–(<ref>) and A:=A(K_P,K_Q) is the Jacobian matrix of f. Develop a systematic and computationally efficient framework for tuning the frequency droop gains  K_P and voltage droop gains K_Q so that asymptotic stability of the equilibrium of (<ref>) is guaranteed.
		
	The rest of the paper is devoted to addressing Problem <ref>.
	

§ MAIN RESULTS

	In this section, we present our results on the stability problem for droop-controlled inverter-based power grids. Due to space limitation, we only provide (in the Appendix) the proofs of the main results in this section.
	

 §.§ Inverter-based Power Grid in Multi-parameter Singularly Perturbed Form

	We begin by making the reasonable assumption that the time-constants of the low-pass filters are the same.
	[Homogeneous filter time constants]
		
    T_P,1=...=T_P,N=T_Q,1=...=T_Q,N=ε_1∈ℝ_+

	


	Under Assumption <ref>, the linearized dynamics  of an inverter-based power grid described by (<ref>) can be written in the  multi-parameter singularly perturbed form <cit.>, <cit.>
	
    d x/dt   =A_xzz
    ε_1 dz/dt   = A_zzz+A_zxx+A_zyy
    
    	E_2 dy/dt   = A_yzz+A_yxx+A_yyy

	where the state-vector w can be broken down to:
	
    x   =[Δδ_1,...,Δδ_N]^T∈ℝ^N
    
    	z   =[Δω_1,...,Δω_N, Δ V_1,...,Δ V_N]^T∈ℝ^2N
    
    	y   =[Δ I_D,1,...,Δ I_D,M , Δ I_Q,1,...,Δ I_Q,M]^T∈ℝ^2M

We note that, all the matrices that appear in (<ref>)–(<ref>) can be obtained appropriately from the Jacobian matrix A.	The matrix E_2∈ℝ^2M× 2M can be compactly expressed as
	
    E_2:=[ ({ε_2,i}_i=1^M)          0_M× M;          0_M× M ({ε_2,i}_i=1^M) ]

	where ε_2,i:=(L_i/ω_b). 	It is evident that the dynamic behavior of (<ref>)-(<ref>) is dictated by (M+2) distinct time scales t, t_z:=(t/ε_1) and t_y,i:=(t/ε_2,i) where i∈ℰ.
	

 §.§ From Multi-parameter to Standard Singularly Perturbed Form

We assume that the values of the line inductances ε_2,i are all of the same order and this order is significantly different from the order of the time constant of the filters ε_1. This assumption allows us to specify ε_2:=(ε_2,1···ε_2,M)^1/M and bring the inverter-based power grid (<ref>)-(<ref>) to the more convenient standard singularly perturbed form <cit.>, <cit.>
	
    d x/dt   =A_xzz
    ε_1 dz/dt   = A_zzz+A_zxx+A_zyy
    ε_2 dy/dt   = D_2(A_yzz+A_yxx+A_yyy)

	where the new matrix D_2 is given by
	
    D_2:=[ ({ε_2/ε_2,i}_i=1^M)              0_M× M;              0_M× M ({ε_2/ε_2,i}_i=1^M) ]

	It should be now obvious that the dynamic behavior of the new system (<ref>)–(<ref>) is dictated by three timescales  t, t_z:=(t/ε_1) and t_y:=(t/ε_2) with the small parameters  ε_1 and ε_2 giving rise to these timescales.  We designate x,z and y to be the slow, fast and very-fast state-variables, respectively.
	

 §.§ Suppressing the Electromagnetic Network Dynamics

	Our goal here is to derive a reduced-order model for the inverter-based power grid (<ref>)–(<ref>) by suppressing the network dynamics. We start by expressing (<ref>)–(<ref>) with respect to t_z=(t/ε_1) in the form <cit.>
	
    d x/dt_z   =A_xzε_1 z
    dz/dt_z   = A_zzz+A_zxx+A_zyy
    ε_3 dy/dt_z   = D_2(A_yzz+A_yxx+A_yyy)

	where ε_3:=(ε_2/ε_1). By writing the system in this form, we uncover the timescale separation between the dynamics of the inverters' states (x,z) and the network's states y. 
	

 §.§ Very-fast Boundary-layer Subsystem

	Our next task is to derive the very-fast boundary-layer electromagnetic dynamics of the network and investigate their stability properties. To do that, we first compute the zeroth-order manifold of y,  y_0(x,z)=A_0zz+A_0xx	where A_0z:=(-A_yy^-1A_yz), A_0x:=(-A_yy^-1A_yx) and A_yy=[ -R  X; -X -R ]. It is easy to notice that A_yy is always invertible.  We perform change of variables ξ:=(y-y_0) to obtain from (<ref>)
	
    dξ/dτ=Ã_yyξ,ξ∈ℝ^2M

	where τ=(t_z/ε_3), Ã_yy:=(D_2A_yy). 	Practically, (<ref>) describes the trajectories of the network's currents when seen decoupled from the much slower dynamics of the output voltages of the inverters. 
	

 §.§ Stability of Very-fast Boundary-layer Subsystem

	At this point, we establish stability of the equilibrium of (<ref>) through the following lemma.
	
		The equilibrium ξ^*=0_2M of the very-fast boundary layer subsystem (<ref>) is asymptotically stable.
		
		
 Intuitively, Lemma <ref> says that the deviations of the network's currents ξ converge to  ξ^*=0_2M, as τ→∞, when the inverter's states  x,z are “frozen”.  The following corollary ensues from Lemma <ref>.
 
 	Let the matrix Ã_yy describing the decoupled transmission line dynamics be Hurwitz. Then, W:=ξ^TP_ξξ with P_ξ≻ 0  is a Lyapunov function for the very fast boundary-layer system (<ref>) satisfying
 	
    ∂ W/∂ξμ(x,η,ξ+y_0)   ≤ -α_4 ψ_4^2(ξ)
    
     	-∂ W/∂ξ[ ∂ y_0/∂ x  ∂ y_0/∂η ]h(x,η,ξ+y_0)   ≤
       β_4ψ_3(x,η)ψ_4(ξ)+γ_2ψ_4^2(ξ)

 	where  α_4:=λ_min(Q_ξ) with Q_ξ:=(Ã_yy^TP_ξ+P_ξÃ_yy). Further, ψ_4(ξ):=ξ_2 and ψ_3(x,η):=ψ_1(x), ψ_2(η)_2  with β_4:=max(σ_max(Θ_ξ x),σ_max(Θ_ξη))√(2) and γ_2:=-λ_min(Θ_ξξ), where
 	
    Θ_ξ x   :=	(P_ξ+P_ξ^T)[(A_0x+A_0zΓ_0)A_xzε_1Γ_0
       +A_0z(A_zzΓ_0+A_zx-Γ_0A_xzε_1Γ_0)
       +A_0zA_zy(A_0x+A_0zΓ_0))  ]
    Θ_ξη   :=(P_ξ+P_ξ^T)[(A_0x+A_0zΓ_0)A_xzε_1
       +A_0z(A_zz-Γ_0A_xzε_1+A_zyA_0z)]
    Θ_ξξ   :=(P_ξ+P_ξ^T)A_0zA_zy

 	
 	
In the sequel, we examine the stability properties of the reduced-order power grid, obtained by suppressing the electromagnetic network dynamics.
		

 §.§ Reduced-order Model of Inverter-based Power Grid

The reduced-order power grid model can be obtained from (<ref>)–(<ref>) upon substitution of the very-fast manifold y_0 <cit.>
	
    dx/dt   =A_xzz 
    ε_1dz/dt   =Ã_zzz+Ã_zxx
 
	where Ã_zz:=(A_zz+A_zyA_0z)∈ℝ^2N× 2N and Ã_zx:=(A_zx+A_zyA_0x)∈ℝ^2N× N. As the matrices Ã_zz and Ã_zx depend explicitly of the droop gains K_P and K_Q, stability of (<ref>)–(<ref>) will largely rely on their chosen values. Given that, we next focus on deriving conditions on the droop gains K_P and K_Q that assure stability of the equilibrium point (x^*, z^*)=(0_N,0_2N) of (<ref>)–(<ref>). 
	

 §.§ Fast Boundary-layer Subsystem

	It is easy to notice that the system (<ref>)–(<ref>) is in singularly perturbed form with two timescales t and t_z. The timescale separation here can be attributed to the time-constant of the inverters' low-pass filters. Once again, we employ singular perturbation <cit.> to analyze this system. We first compute the zeroth-order manifold z_0(x)=Γ_0x where Γ_0=(-Ã_zz^-1Ã_zx) and use a change of variables η:=(z-z_0) to recover the fast boundary-layer subsystem 
	
    dη/dt_z=Ã_zzη,η∈ℝ^2N

	where Ã_zz:=Ã_zz(K_P,K_Q)∈ℝ^2N× 2N.  This system characterizes the dynamics of the frequencies ω and amplitudes V of the inverters' output voltages when decoupled from the much slower dynamics of the voltage angles δ.
	

 §.§ Stability of Fast Boundary-layer Subsystem

	Through algebraic manipulations, we can obtain the matrix Ã_zz in the upper triangular form
	
    Ã_zz:=[   -I_N      ⋆; 0_N× N      E ]

	with the matrix E={e_ij}_i,j∈𝒩 being defined  as follows:
	
    e_ij=k_q,iν_i-1, i=j
    
    	k_q,iν_ij, i≠ jand  (j,i)∈ℰ
    
    	0, otherwise
 
	The terms ν_i and ν_ij can be expressed as:
	
    ν_i   = ∑_k∈ℰ_i[dQ_i/dI_D,k(R_k cos(δ_i)+X_ksin(δ_i)/R_k^2+X_k^2)
       -dQ_i/dI_Q,k(X_k cos(δ_i)-R_ksin(δ_i)/R_k^2+X_k^2)]θ_k-dQ_i/dV_i,  i∈𝒩

	
    ν_ij   =[-dQ_i/dI_D,k(R_k cos(δ_j)+X_ksin(δ_j)/R_k^2+X_k^2)
       +dQ_i/dI_Q,k(X_kcos(δ_j)-R_ksin(δ_j)/R_k^2+X_k^2)]θ_k,
        i,j∈𝒩,  k=(i,j)∈ℰ

	where
	
    θ_k=1,  if inverter  i  lies at the end of line  k
    
    	-1,  if inverter  i  lies in the beginning of line  k

	and ℰ_i denotes the set of distribution lines connected to inverter i.
	The next lemma affirms stability of  (<ref>).
	
		The equilibrium η^*=0_2N of the fast boundary-layer subsystem  (<ref>) is asymptotically stable when the voltage droop gains satisfy K_Q∈𝒜 where 𝒜:={K_Q∈ℝ^N  | E  is Hurwitz }.	
		
It is easy to see that the matrix E is Metzler. A strictly diagonally dominant Metzler matrix is stable; therefore, we can extract from Lemma <ref>, the following stability conditions. 
	
		Let ν_i<0, for every inverter i∈𝒩, hold. Then, the equilibrium η^*=0_2N of the subsystem (<ref>) is asymptotically stable when the voltage droop gains of the inverters satisfy:
		
    k_q,i   >0,        if  |ν_i|>=∑_j∈𝒩_i|ν_ij|
    
    		k_q,i   <1/∑_j∈𝒩_i|ν_ij|-|ν_i|,       if  |ν_i|<∑_j∈𝒩_i|ν_ij|

		where 𝒩_i is the set of inverters adjacent to inverter i.
		
		
	
	The main implication of the above lemma is that inverters can exploit decentralized stability criteria and tune their voltage droop gains properly using only information from the neighbors in order to assure stability of their fast dynamics. The following corollary is a byproduct of the above lemmas. Its proof can be found in the Appendix.
	
		Let f(x,z):=A_xzz,  	g(x,z):=Ã_zxz+Ã_zzz and the voltage droop gains K_Q∈ℝ^N be chosen so that E and thus Ã_zz are Hurwitz. Then, V_f:=η^TH η∈ℝ_+ with H≻ 0 is a Lyapunov function for the fast subsystem (<ref>) satisfying the following two inequalities
		
    ∂ V_f/∂ xg(x,η+z_0)   ≤ -α_2 ψ_2^2(η)
    
    		-∂ V_f/∂ηf(x,η+z_0)   ≤γ_1 ψ_2^2(η)+β_2ψ_1(x)ψ_2(η)

		with  α_2:=λ_min(Λ) where Λ:=-(Ã_zz^TH+HÃ_zz). Further, ψ_1(x):=x_2, ψ_2(η):=η_2, γ_1:=-λ_min(Z), Z:=(HΓ_0A_xz+H^TΓ_0A_xz), β_2:=σ_max(Θ), Θ:=ZΓ_0. 
		
	
	The key insight here is that stability of the fast dynamics of (<ref>) depends only on the voltage droop gains K_Q and is completely independent of the frequency gains K_P. 
	

 §.§ Slow  Subsystem

	The slow subsystem can be obtained as
	
    dx/dt=A^(s)x,  A^(s)∈ℝ^N× N

	where A^(s)=(-A_xzÃ_zz^-1Ã_zx). 
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	The following lemma puts forward a stability condition for A^(s).
	
		The equilibrium x^*=0_N of the slow dynamics is asymptotically stable when K_Q, K_P∈ℬ where ℬ:={K_Q, K_P∈ℝ^N  | A^(s) is Hurwitz }.	
		
	The following corollary naturally arises for system (<ref>). Its proof can be found in the Appendix.
	
		Let the voltage and frequency droop gains K_Q and K_P be chosen so that A^(s) is Hurwitz. Then, V_s:=x^TPx with P≻ 0 is a Lyapunov function for the slow  subsystem (<ref>) of the reduced inverter-based power grid  satisfying the following two inequalities
		
    ∂ V_s/∂ xf(x,z_0)   ≤ -α_1 ψ_1^2(x)
    ∂ V_s/∂ x[f(x,η+z_0)-f(x,z_0)]   ≤β_1 ψ_1(x)ψ_2(η)

		with α_1:=λ_min(Q^(s)), Q^(s):=-(A^(s)^TP+PA^(s)),  ψ_1(x):=x_2, ψ_2(η):=η_2, β_1:=σ_max((P+P^T)A_xz). 
		
	
	Next, we examine stability of the reduced power grid. 
	

 §.§ Stability of Reduced-order Inverter-based Power Grid

We will now explore what other additional condition is required, for stability of the original reduced-order system (<ref>)–(<ref>). We start by recasting system (<ref>)–(<ref>) using the manifold z_0 and the change of variables η:=(z-z_0)  as
	
    dx/dt   =A_xzΓ_0x+A_xzη
    ε_1dη/dt   =(-ε_1Γ_0A_xzΓ_0)x+(Ã_zz-ε_1Γ_0A_xz)η

 Next, we state Theorem<ref> whose proof can be found in the Appendix.
	
		Consider the reduced-order inverter-based power grid (<ref>)–(<ref>). Let Corollaries <ref> and <ref> hold with  V_f:=η^TH η and V_s:=x^TPx being Lyapunov functions for the fast and slow subsystems, (<ref>) and  (<ref>), respectively. Then, v(x,η):=(1-d_1)V_s(x)+d_1V_f(η) with d_1=d_1^*=β_1/(β_1+β_2) is a Lyapunov function for the system (<ref>)–(<ref>) and the equilibrium   (x^*,η^*)=(0_N,0_2N) of the system  (<ref>)–(<ref>) is asymptotically stable when 
		
    ε_1<ε_1^*:=α_1α_2/α_1γ_1+β_1β_2

		where α_1:=λ_min(Q^(s)) with Q^(s):=-(A^(s)^TP+PA^(s)) and α_2:=λ_min(Λ) with Λ:=-(Ã_zz^TH+HÃ_zz). Further, β_1:=σ_max((P+P^T)A_xz), γ_1:=-λ_min(Z), Z:=(HΓ_0A_xz+H^TΓ_0A_xz), β_2:=σ_max(Θ) and Θ:=ZΓ_0. 
		
	
So the reduced-order system (<ref>)–(<ref>) would be stable when ε_1 respects the upper bound in (<ref>) and the voltage droop gains K_Q and frequency droop gains K_P are chosen so that the fast and slow subsystems (<ref>) and (<ref>) are both stable. The following corollary (whose proof can be found in the Appendix) stems from Theorem <ref>.
	
		Let Theorem 1 hold with v(x,η)  being a Lyapunov function for the reduced inverter-based power grid (<ref>) and (<ref>). Then v(x,η) satisfies
		=0mu
		
    [ ∂ V/∂ x  ∂ V/∂η ]h(x,η,y_0)≤ -α_3 ψ_3^2(x,η)
       [ ∂ V/∂ x  ∂ V/∂η ][h(x,η,ξ+y_0)-h(x,η,ξ+y_0)]≤β_3 ψ_3(x,η)ψ_4(ξ)

		where ψ_3(x,η):=ψ_1(x), ψ_2(η)_2, ψ_4(ξ):=ξ_2, α_3:=ε_1λ_min(Q_v) and β_3:=σ_max((H+H^T)A_zy) with H:=d_1H. 
	
We are now ready to find conditions for stability of the overall system with the electromagnetic network dynamics.
	

 §.§ Stability of Full Inverter-based Power Grid

 We begin by specifying
 =0mu
	
    h(x,η,y):=
       [                                   A_xzε_1Γ_0x+A_xzε_1η;                                                        (A_xz+A_zzΓ_0-Γ_0A_xzε_1Γ_0)x+(A_zz-Γ_0A_xzε_1)η+A_zyy ]








	With the change of variables η=(z-z_0) and ξ:=(y-y_0), we can write the full system (<ref>)–(<ref>) with respect to time-scale t_z as
	
    dx/dt_z   =A_xzε_1Γ_0x+A_xzε_1η
    dη/dt_z   =(A_xz+A_zzΓ_0-Γ_0A_xzε_1Γ_0+A_zy(A_0zΓ_0
       +A_0x))x	+(A_zz-Γ_0A_xzε_1+A_zyA_0z)η+A_zyξ
    ε_3dξ/dt_z   =Ã_yyξ-ε_3A_0zdη/dt_z-ε_3(A_0x+A_0zΓ_0)dx/dt_z

	Clearly, our final goal is to find conditions for stability of (<ref>)–(<ref>). We designate
	
    μ(x,η,y):=D_2[(A_yx+A_yzΓ_0)x+A_yzη+A_yyy]∈ℝ^2M

	and compute μ(x,η,ξ+y_0):=Ã_yyξ. By combining the stability properties of the decoupled reduced-order power grid and the network dynamics, we arrive at our final result.
	
		Consider the full inverter-based power grid (<ref>)–(<ref>). Let Corollaries <ref> and <ref> hold with  W:=ξ^TP_ξξ being a Lyapunov function for the decoupled electromagnetic network dynamics (<ref>) and v(x,η)  a Lyapunov function for the reduced-order power grid (<ref>)–(<ref>). Then, 𝒱(x,η,ξ):=(1-d_2)v(x,η)+d_2W(ξ) with d_2:=β_3/(β_3+β_4) is a Lyapunov function for the system (<ref>)–(<ref>) and the equilibrium   
		(x^*,η^*,ξ^*)=(0_N,0_2N,0_2M) of the full inverter-based power grid (<ref>)–(<ref>) is asymptotically stable as long as
		
    ε_3:=(L_1/ω_b···L_M/ω_b)^1/M/T_P<	ε_3^*:=α_3α_4/α_3γ_2+β_3β_4

	 Further, α_3:=ε_1λ_min(Q_v), α_4:=λ_min(Q_ξ) where Q_ξ:=-(Ã_yy^TP_ξ+P_ξÃ_yy), β_3:=σ_max((H+H^T)A_zy), H:=d_1H, β_4:=max(σ_max(Θ_ξ x),σ_max(Θ_ξη)), γ_2:=-λ_min(Θ_ξξ) with Θ_ξξ given by (<ref>).
		
		
	Essentially, this theorem says that we have stability of the original inverter-based power grid (<ref>)–(<ref>) when: a) the reduced-order approximated system obtained by suppressing the network dynamics is stable (Theorem <ref>) and, b) the parameter ε_3 respects the upper bound ε_3^*. 
	Our analytical framework can be a useful tool for systematic and efficient stability analysis of  inverter-based power grids. The traditional approach to stability analysis usually involves adjusting the droop gains K_P and K_Q and repeatedly computing the eigenvalues of the (3N+2M)-dimensional full system matrix  until this becomes Hurwitz.   With our framework, one has to only assure that the N-dimensional matrices E and A^(s) are Hurwitz and the parameters ε_1 and ε_3 respect some well-defined bounds. 
	
	

§ NUMERICAL VALIDATION

	We corroborate our theoretical results numerically and illustrate how our framework can be practically implemented via an example on the IEEE 13-bus distribution grid.
	

 §.§ Set-up

	We use the standard IEEE 13-bus test feeder <cit.> model to validate our results which we modify by placing a single inverter at each bus. Further, we reduce four buses of the original system to two buses in the modified system. Each inverter i, where i∈𝒩, has a power rating S_n,i=10kVA and low-pass filter time-constant T_P,i=T_Q,i=0.0318s <cit.>.  Here, 𝒩 denotes the set of the 10 inverters.  In our example we choose V_b=381.58V as the base voltage, S_b=10kVA as the base power and compute the base impedance as Z_b=(V_b^2/S_b)=14.56Ω. 
	

 §.§ Implementation of the Proposed Framework

 By applying the conditions in Lemma <ref> we arrive at the following inequalities k_q,8<0.2,   k_q,i>0,∀ i∈𝒩∖{8}. Choosing k_q,i=0.05 yields the eigenvalues of E shown in Fig. <ref>. Although these inequalities are not very restrictive for the voltage droop gains K_Q, recall that, these droop gains have to still be chosen carefully, in conjuction with the frequency droop gains K_P, so that A^(s) is Hurwitz.








	

			
	



	

		
	



	Having chosen the voltage droop gains k_q,i, we then compute A^(s) symbolically in terms of the frequency droop gains K_P. We now have to choose frequency droop gains that yield a Hurwitz matrix A^(s). One suitable choice for these gains is k_p,i=0.6, ∀ i∈𝒩, as can be coroborrated by Fig. <ref>. The largest eigenvalue of A^(s) is -0.893. To proceed, we construct Lyapunov functions V_s=x^TPx and V_f=η^T Hη for the slow and fast subsystems by trivially choosing Q^(s)=I_N and Λ=I_2N where A^(s)^TP+PA^(s)=-Q^(s) and Ã_zz^TH+HÃ_zz=-Λ.  With these Lyapunov functions, we obtain through Theorem <ref> the following upper bound for the time-constant of the inverters' low-pass filters ε_1^*=0.0178· 10^-5. 	This bound reflects the maximum time constant of inverters' low-pass filters for which stability of E and A^(s) readily translates into stability of the reduced-order inverter-based power grid. For the particular Lyapunov functions that we chose here, this bound turns out to be quite convervative as the value of ε_1, which corresponds to the time-constants T_P,i, T_Q,i, is 0.0318. One has to experiment with different Lyapunov functions in order to arrive at a less conservative bound. This is beyond the scope of this work, but could be considered in future work.
	We now focus on the decoupled network dynamics which are always stable. We let Q_ξ=I_2M where	Ã_yy^T P_ξ+P_ξÃ_yy=-Q_ξ and construct a Lyapunov function W:=ξ^TP_ξξ for these dynamics. We let ε_1=ε_1^*/2 and use Theorem <ref> to obtain ε_3^*=7.84· 10^-10. As the actual value of ε_3, which can be computed using (<ref>), is 6.8705· 10^-4 we see that this bound is also quite conservative. 
		Overall, our numerical example shows that the bounds ε_1^* and ε_3^* can be quite conservative for certain  Lyapunov functions. One could try to come up with “better" Lyapunov functions that would lead to less conservative upper bounds.  As the conditions are only sufficient, the full system may still be stable even when they are not met. To illustrate this, we compute  the eigenvalues of the full system matrix A with the chosen droop gains, k_q,i=0.05 and k_p,i=0.6. As can be seen from Fig. <ref>,  all eigenvalues of A are negative despite the fact that the bounds ε_1^* and ε_3^* are not respected. 
	

		
	

	
	
	

§ CONCLUSIONS

	We studied the stability problem for distribution grids with droop-controlled inverters and electromagnetic network dynamics. We presented a systematic framework that builds on singular perturbation and Lyapunov theories for carrying out  stability analysis in a compositional and efficient manner. By deploying our framework, we derived sufficient stability conditions for the full system. Our theoretical results are numerically corroborated via an example on the IEEE 13-bus distribution grid.
	


	
	*
	unsrt
	
	Proof of Corollary <ref>
		
		When Ã_zz is Hurwitz, along the trajectories of the fast boundary-layer subsystem (<ref>) we have
		
    ∂ V_f/∂ηg(x,η+z_0)=η^T(Ã_zz^T H+HÃ_zz)η
       =-η^TΛη≤ -λ_min(Λ)η_2^2=-α_2ψ_2^2(η)

		where α_2:=λ_min(Λ), ψ_2(η):=η_2 and Λ≻ 0. With this, we conclude that V_f is a Lyapunov function for the fast system (<ref>). We are left to show that the perturbation term in (<ref>) is bounded appropriately. Expanding this term yields
		
    (   -∂ V_f/∂ y∂ z_0/∂ x)f(x,η+z_0)=-η^T Z η -η^TΘ x 
       ≤ -λ_min(Z)η_2^2+σ_max(Θ)η_2x_2
       =γ_1 ψ_2^2(η)+β_2ψ_1(x)ψ_2(η)

		where  ψ_1(x):=x_2, γ_1:=-λ_min(Z) and Z:=(HΓ_0A_xz+H^TΓ_0A_xz). Further, β_2:=σ_max(Θ) where Θ:=ZΓ_0 and with that we  complete the proof.
	
		Proof of Corollary <ref>
		
		With A^(s) being Hurwitz, along the trajectories of the slow subsystem (<ref>) we obtain
		
    ∂ V_s/∂ xf(x,z_0)   =x^T(A^(s)^TP+PA^(s))≤-α_1ψ_1^2(x)

		where α_1:=λ_min(Q^(s)), ψ_1(x):=x_2 and Q^(s)≻ 0. With this, we deduce that V_s is a Lyapunov function for the slow system (<ref>). Eventually, we can also bound the term:
		
    ∂ V_s/∂ x[f(x,η+z_0(x))-f(x,z_0(x))]=x^T(P+P^T)A_xzη
       ≤σ_max((P+P^T)A_xz)x_2η_2=β_1 ψ_1(x) ψ_2(η)

		where β_1:=σ_max((P+P^T)A_xz), ψ_2(η):=η_2 and conclude the proof.
		


		Proof of Theorem <ref>
	
		We construct a composite candidate Lyapunov function for the system (<ref>)-(<ref>) v(x,η):=(1-d_1)V_s(x)+d_1V_f(η).	Calculating the derivative of v(x,η) along the trajectories of the full system (<ref>)-(<ref>), we obtain
		
    dv/dt   =(1-d_1)∂ V_s/∂ xf(x,z_0)+(1-d_1)∂ V_s/∂ x[f(x,η+z_0)
       -f(x,z_0)]+d_1/ε_1∂ V_f/∂ηg(x, η+z_0)-d_1∂ V_f/∂η∂ z_0/∂ x f(x,η+z_0)

		Applying the inequalities  (<ref>), (<ref>), (<ref>), (<ref>) yields
		
    dv/dt   ≤ -(1-d_1)α_1 ψ_1^2(x)-d_1/ε_1α_2 ψ_2^2(η)
       +(1-d_1)β_1 ψ_1(x)ψ_2(η)+d_1γ_1ψ_2^2(η)+d_1β_2 ψ_1(x)ψ_2(η)

		which can be expressed in quadratic form as
		=0mu
		
    dv/dt≤ -[ ψ_1(x); ψ_2(η) ]^TQ_v[ ψ_1(x); ψ_2(η) ]

		where
		
    Q_v:=[               (1-d_1)α_1 -1/2(1-d_1)β_1-1/2d_1β_2; -1/2(1-d_1)β_1-1/2d_1β_2         d_1(α_2/ε_1-γ_1) ]

		Positive definiteness of Q_v is guaranteed when
		
    ε_1 ≤α_1α_2/α_1γ_1+1/4d_1(1-d_1)[(1-d_1)β_1+d_1β_2]^2_:=ε_1,d.

		Choosing d_1^*:=β_1/(β_1+β_2), yields the maximum value of ε_1^*=α_1α_2/α_1γ_1+β_1β_2. That, completes the proof.
	
	Proof of Corollary <ref>
		
		We specify ρ:=[ψ_1(x),  ψ_2(η)]^T, ψ_3(x,η):=ρ_2 and ψ_4(ξ):=ξ_2. Let W(ξ):=ξ^TP_ξξ be a candidate Lyapunov function for the very fast dynamics (<ref>). The derivative of W along the trajectories of (<ref>) is
		
    ∂ W/∂ξμ(x,η,ξ+y_0)   :=ξ^T(Ã_yy^TP_ξ+P_ξÃ_yy)ξ
       ≤ -λ_min(Q_ξ)ξ_2^2= -α_4 ψ_4^2(ξ)

		where, as Q_ξ≻ 0 (due to Ã_yy being Hurwitz) we conclude that W is a Lyapunov function for  (<ref>). Finally, we also have
		
    -∂ W/∂ξ[ ∂ y_0/∂ x  ∂ y_0/∂η ]h(x,η,ξ+y_0)
       =-ξ^T Θ_ξ xx-ξ^T Θ_ξηη-ξ^T Θ_ξξξ
       ≤max (σ_max(Θ_ξ x),σ_max(Θ_ξη))ψ_4(ξ) √(2)ψ_3(x,η)
        -λ_min(Θ_ξξ)ψ_4^2(ξ):= β_4ψ_3(x,η)ψ_4(ξ)+γ_2ψ_4^2(ξ)

		where β_4:=max (σ_max(Θ_ξ x),σ_max(Θ_ξη))√(2) and γ_2:=-λ_min(Θ_ξξ). We note that the matrices 
		Θ_ξ x, Θ_ξη, Θ_ξξ are as specified in (<ref>), (<ref>), (<ref>) and complete the proof.	
		
	Proof of Corollary <ref>
		
		From (<ref>), we readily have that
		
    [ ∂ v/∂ x ∂ v/∂ n ]h(x,η,y_0)≤-ε_1λ_min(Q_v)ψ_3^2(x,η)=-α_3ψ_3^2(x,η)

		where α_3:=ε_1λ_min(Q_v).		By letting P:=(1-d_1)P and H:=d_1H one can express v(x,η) as v(x,η):=x^TPx+η^T Hη and derive the following inequality
		
    [ ∂ V/∂ x  ∂ V/∂η ][h(x,η,ξ+y_0)-h(x,η,y_0)]
       =η^T (H+H^T)A_zyξ≤σ_max((H+H^T)A_zy)η_2ξ_2
       ≤β_3ψ_3(x,η)ψ_4(ξ)
	
		With that, we complete the proof.
		
		

	Proof of Theorem <ref>
		
	Similarly as before, we can employ a candidate Lyapunov  𝒱(x,η,ξ):=(1-d_2)v(x,η)+d_2W(ξ) and use the inequalities (<ref>), (<ref>), (<ref>), (<ref>) to finally obtain ε_3^*=α_3α_4/α_3γ_2+β_3β_4 where  d_2:=β_3/(β_3+β_4). That, completes the proof.
		
	
	
	
