Participation Analysis in Impedance Models: The Grey-Box Approach for Power System Stability
CCopyright Statements
This work has been submitted to the IEEE for possible publication. Copyright may betransferred without notice, after which this version may no longer be accessible. a r X i v : . [ ee ss . S Y ] F e b his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 1 Participation Analysis in Impedance Models: TheGrey-Box Approach for Power System Stability
Yue Zhu,
Student Member, IEEE , Yunjie Gu,
Senior Member, IEEE , Yitong Li,
Member, IEEE ,Timothy C. Green,
Fellow, IEEE
Abstract —This paper develops a grey-box approach to stabilityanalysis of complex power systems that facilitates root-causetracing without requiring disclosure of the full details of theinternal control structure of apparatus connected to the system.The grey-box enhances the impedance transfer function approachpopular in power electronic systems and favoured by manufactur-ers of wind and solar systems for the limited disclosure required.Impedance participation factors are proposed and defined interms of the residues of the whole-system admittance matrix. Itis proved that, so defined, these impedance participation factorsequal the sensitivity of whole-system eigenvalues to apparatusimpedances. The classic state participation factors are related tothe impedance participation factors via a chain-rule. Based on thechain-rule, a three-layer grey-box approach, with three degreesof transparency, is proposed for root-cause tracing to differentdepths, i.e. apparatus, states, and parameters, according to theinformation available. The association of impedance participationfactors with eigenvalue sensitivity points to the re-tuning thatwould stabilize the system. The impedance participation factorscan be measured in the field or calculated from black-boximpedance spectra with little prior-knowledge required.
Index Terms —Impedance, Admittance, Participation Factor,Inverter-Based Resource, Eigenvalue Sensitivity N OMENCLATURE
K, k total nodes and node index, k ∈ { , , · · · , K } N total state number n, m state and mode index, n, m ∈ { , , · · · , N } ˆ Y whole-system admittance matrix Y net nodal admittance matrix for the network Z terminal impedance of grid-connected apparatus x, z state and mode vectors A, a mn state matrix and its elements Λ , λ eigen matrix and eigenvalue Ψ , ψ nm left-eigenvector matrix and its elements Φ , φ mn right-eigenvector matrix and its elements ρ system parameterRes λ G residue of G at λp mn state participation factor p λ,Z impedance participation factor p λ,ρ parameter participation factor (cid:104)· , ·(cid:105) Frobenius inner product (cid:107) · (cid:107)
Frobenius norm (cid:62) matrix transposecomplex conjugate ∗ conjugate transpose: ( · ) ∗ = ( · ) (cid:62) I. I
NTRODUCTION
Stability analysis for power systems must keep pace with thefast changing characteristics of the systems caused by emerg-ing inverter-based resources (IBRs) replacing synchronous generators (SGs) and becoming the dominant sources. Unsta-ble oscillations induced by IBRs are reported world-wide andthe characteristics of such oscillations are distinctly differentto the behaviour of a conventional SG-based grid [1]–[6]. Theunderstanding of the mechanisms of IBR induced instabilityis not yet comprehensive and nor are systematic solutionsavailable for ensuring system-wide stability.One of the impediments to stability analysis is the lackof standardised and precise analytical models of the dynamicbehaviour of IBRs. Unlike SGs, whose behaviour is largelydetermined by physics, IBR behaviour is shaped by its internalcontrol system which is extremely flexible and far from beingstandardized as yet [7], [8]. Manufacturers regard their controlalgorithms as critical proprietary technology and prefer to dis-close to system operators only black-box models, such as thefrequency-domain spectrum of the impedance (or equivalently,through duality, admittance) as seen from the terminal of theIBR [5]. For their part, system operators desire white-boxmodels for whole-system analysis. Such white-box models areusually state-space equations containing all of the physical andcontrol states of each component. An important advantage ofstate-space white-box models is the availability of participationfactor analysis that can identify the root cause of unstable orunder-damped oscillations throughout the grid [9]–[12].To bridge the gap between the black-box models divulgedby manufacturers and the white-box models desired by opera-tors, there have been efforts recently to incorporate impedance-based models into whole-system analysis. A whole-systemimpedance model is presented in [13], which formulates thedynamics of the grid as the whole-system admittance andimpedance matrix in the frequency domain. All of the elementsin the admittance and impedance matrices share the same poles(equivalent to the eigenvalues of the system) and thereforecontain the same information regarding system stability, butdifferent elements may have different resonant peaks at eachof the poles which reflect their differing participation in thecorresponding mode. Further to the whole-system model, [14]and [15] introduce nodal and branch participation factors basedon the diagonalization of the nodal admittance matrix andloop impedance matrix. These two approaches open up pathstowards participation analysis with impedance models, but therelationship between the participation factors in impedancemodels and those in state-space models is still unclear. As aresult, the participation factors in impedance models, definedin this way, cannot be used to look inside a black-box andso cannot determine exactly which state causes an instabilitynor indicate how internal parameters should be re-tuned tostabilize the system.1his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 2In this paper, we first prove that the sensitivities of state-space eigenvalues to apparatus impedances equal the residuesof the whole-system admittance seen at the correspondingnode. Here we use apparatus to refer to any device connectedin shunt at a node such as an SG, including a synchronouscondenser, an IBR or a FACTS device. Based on this finding,we define the residues as impedance participation factors ina way that is consistent with the classic state participationfactors. We further prove that the chain-rule for impedanceparticipation factors which links the state and parameter par-ticipation factors. This key step allows one to look insideblack-box models without disclosing the internal details. Tohighlight this feature, we name our method the grey-boxapproach . The proposed grey-box approach has not onlya rigorous mathematical basis but also good potential forpractical application. The residues can be measured in thefield or estimated from black-box impedance spectrum withvery little prior-knowledge and few presumptions.The paper is organized as follows. The whole-system ad-mittance model is briefly introduced in Section II. The theoryof impedance participation factor and the corresponding grey-box approach is presented in Section III. The major findingsare verified and applications demonstrated on an IEEE 14-bussystem in Section IV. The last section concludes the paper.II. W
HOLE -S YSTEM A DMITTANCE M ODEL
This section briefly describes the whole-system admittancemodel proposed in [13] based on which participation analysiswill be conducted in the next section. In a K -node meshednetwork, as illustrated in Fig. 1 (a), virtual voltage injec-tions ˆ v = [ˆ v , ˆ v , · · · , ˆ v K ] (cid:62) are introduced in series witheach item of apparatus in the system, such as a SG or anIBR, to create perturbations in the corresponding currents ˆ i = [ˆ i , ˆ i , · · · , ˆ i K ] (cid:62) . The transfer function matrix from ˆ v to ˆ i is called the whole-system admittance ˆ Y , that is, ˆ i ( s ) = ˆ Y ( s ) · ˆ v ( s ) .The whole-system admittance ˆ Y can be calculated from theequivalent feedback block diagram in Fig. 1 (b) ˆ Y = ( I + Y net Z ) − Y net (1)where Z = diag ( Z , Z , · · · , Z K ) (2)is a diagonal matrix containing the terminal impedance ofall apparatus, and Y net is the nodal admittance matrix of thenetwork. ˆ Y is a K × K transfer function matrix. All elements of ˆ Y share the same poles which is equivalent to the eigenvaluesof the system [13]. On the other hand, different elements mayhave different levels of excitation and thus different magnitudeof resonant peaks at each pole of the system which reflect thedifferent participation levels in each mode. This characteristicserves as the basis to extract precise participation factors from ˆ Y matrix which is the central theme of this paper.The diagonal elements of ˆ Y have a clear physical meaning.Taking the first diagonal element ˆ Y as an example, asillustrated in Fig. 1 (a), ˆ Y is essentially the admittance forthe loop containing Z and Z g , in which Z is the terminal (a)(b)Fig. 1. Illustration of the whole-system admittance model. (a) Virtual voltageinjection ˆ v to excite current perturbation ˆ i . (b) Equivalent feedback diagramfor ˆ v and ˆ i . impedance of the apparatus at the first node and Z g is theimpedance for the rest of the grid as seen by Z at that node.Therefore, we have ˆ Y = ( Z + Z g ) − and the same principleholds for any node k ˆ Y kk = ( Z k + Z g k ) − . (3)As we shall see in the next section, this property is very usefulfor impedance-based participation analysis.III. I MPEDANCE P ARTICIPATION F ACTOR AND THE G REY -B OX A PPROACH
Before proceeding to impedance-based participation analy-sis, we review the classic state-based participation analysis sothat the relationship between them is revealed to the reader.
A. State Participation Factor
The dynamics of a power system linearized around itsequilibrium point can be represented by a state equation ˙ x = Ax. (4)This state equation is usually very high-order but can be de-composed into a series of first-order equivalents via coordinatetransformations z = Ψ x and x = Φ z ( Φ = Ψ − ) such thatthe state matrix in the new coordinate z is diagonalized [9],that is, ˙ z = Λ z, Λ = Ψ A Φ = diag ( λ , λ , · · · , λ N ) (5)where λ n ( n = 1 , , · · · , N ) is the n -th eigenvalue of A ,and the rows and columns of Ψ and Φ correspond to theleft- and right-eigenvectors of A respectively. z is a vectorof modes of the system and determines stability according tothe corresponding eigenvalues. The transformation matrices Ψ and Φ describe the back-and-forth correlation between modes2his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 3 z and states x . Putting Ψ and Φ together, the participationfactor of the m -th state x m in the n -th mode z n is defined as p mn = ψ nm φ mn . (6)where ψ nm and φ mn are elements of Ψ and Φ respectively.In this paper, p mn is renamed state participation factor todistinguish it from the impedance participation factor to beintroduced in the next subsection.The state participation factor in (6) is defined heuristicallybut is endowed with a rigorous mathematical meaning due toits linkage to eigenvalue sensitivity [9], that is p mn = ∂λ n ∂a mm → ∆ λ n = p mn ∆ a mm . (7)In (7), a mm is the m -th diagonal element of the state matrix A and can be interpreted as the local dynamics for x m itself as ifdecoupled from the other states. In contrast, λ n describes theglobal dynamics for the whole system. Therefore, (7) impliesthat p mn determines how the local dynamics affect the globaldynamics. In the light of this interpretation, the state participa-tion factors can not only identify which local states participatein particular whole-system modes, but can also indicate howlocal parameters should be re-tuned to better damp modes.For this reason, state participation factors became an importanttool in stability analysis enabling tracing of root-causes, andtrouble-shooting of complex power systems [9]–[12].We now add inputs and outputs to the state equation (4)to see how the state participation factors are reflected in thetransfer functions of the system. For a system with input u and output y ˙ x = Ax + Buy = Cx (8)the transfer function from u ( s ) to y ( s ) is G ( s ) = C ( sI − A ) − B = C Φ( sI − Λ) − Ψ B. (9)For the sake of brevity, we take a case of a single-input andsingle-output (SISO) system, i.e., u ( s ) and y ( s ) are scalars.In such a case, G ( s ) can be expanded as G ( s ) = N (cid:88) n =1 Res λ n s − λ n (10)where Res λ n = N (cid:88) m =1 c m φ mn N (cid:88) m =1 ψ nm b m (11)is the residue of G ( s ) at λ n , and b m and c m are elements of B and C respectively. It is clear that the eigenvalues of thestate matrix appear as poles in the transfer function. If C and B are in such a form that their m -th elements equal 1 and allother elements equal 0, Res λ n can be simplified toRes λ n = ψ nm φ mn = p mn . (12)In this special case, the state participation factor is the same asthe residue of the transfer function, an observation which hintsat residues being useful for participation analysis. However,(12) is based on a strong assumption about B and C which may not hold for common cases. In the next subsection, weexplore the general residue-participation relationship whichyields the impedance participation factor. B. Impedance Participation Factor
In order to clarify the roles of residues in participationanalysis in a general form, we introduce Lemma 1.
Lemma 1.
For a square transfer function matrix G ρ depend-ing on parameters ρ , let H ρ be the inverse transfer functionof G ρ , i.e. H ρ = G − ρ , and λ be a non-repeated pole of G ρ .When the parameters ρ are perturbed infinitesimally by ∆ ρ , λ and H ρ are perturbed by ∆ λ and ∆ H ρ correspondingly andwe have the following relationship in between: ∆ λ = (cid:104)− Res ∗ λ G ρ , ∆ H ρ ( λ ) (cid:105) (13) in which (cid:104)· , ·(cid:105) is the Frobenius inner product of two matrices, Res λ G ρ is the residue matrix of G ρ at λ (the residue operateselement-wise on a matrix), and ∗ denotes the conjugatetranspose of the residue matrix.Proof. The proof of the Lemma is given in Appendix B, anda brief introduction to the mathematical preliminaries usedin this Lemma, including the definition of residue and theFrobenius inner product, are included in Appendix A.If we take G ρ to be the whole-system admittance seenat node k defined in Section II, that is, G ρ = ˆ Y kk , thecorresponding H ρ is H ρ = ˆ Y − kk = Z k + Z g k (14)according to (3). When Z k itself is subject to a perturbationwe can say that Z g k , representing the rest of the grid, isunchanged, that is, ∆ Z g k = 0 , so we have ∆ H ρ = ∆ Z k + ∆ Z g k = ∆ Z k (15)and hence ∆ λ = (cid:104)− Res ∗ λ ˆ Y kk , ∆ Z k ( λ ) (cid:105) . (16)As explained in the previous subsection, the poles of thetransfer function ˆ Y kk are exactly the eigenvalues of the system.Therefore, (16) implies that the sensitivity of an eigenvalue toan apparatus impedance is determined by the residue of thewhole-system admittance seen by the same apparatus. Sincesensitivity is equivalent to participation, we define the residueas the impedance participation factor p λ,Z k (cid:44) − Res ∗ λ ˆ Y kk (17)such that ∆ λ = (cid:104) p λ,Z k , ∆ Z k ( λ ) (cid:105) . (18)Due to impedance-admittance duality, we can similarlydefine the admittance participation factor as the residue ofthe whole-system impedance p λ,Y k (cid:44) − Res ∗ λ ˆ Z kk (19)such that ∆ λ = (cid:104) p λ,Y k , ∆ Y k ( λ ) (cid:105) (20)3his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 4 Fig. 2. The relationship between impedance participation factor and stateparticipation factor and the chain-rule for participation propagation. where Y k is the admittance of the k -th apparatus and ˆ Z kk is the whole-system impedance [13] at the k -th node. Theimpedance and admittance participation factors are equivalenttheoretically but each may better serve different applicationswhere impedance or admittance is more readily available. Wefocus on the impedance participation factor in this paper.If we know the sensitivity of the apparatus impedanceagainst its parameters ρ , that is, ∆ Z k ( λ ) = ∂Z k ( λ ) ∂ρ · ∆ ρ (21)we further define the parameter participation factor p λ,ρ = (cid:28) p λ,Z k , ∂Z k ( λ ) ∂ρ (cid:29) (22)such that ∆ λ = p λ,ρ ∆ ρ. (23)If ρ is selected as a mm in the state matrix A , the correspondingparameter participation factor is the state participation factor p mn = p λ n ,a mm = (cid:28) p λ n ,Z k , ∂Z k ( λ n ) ∂a mm (cid:29) . (24)Thus we establish the relationship between the three types ofparticipation factors and summarize this relationship in Fig. 2.Equations (22) and (24) are called the chain-rule of partici-pation factors which is of profound importance in participationanalysis. The impedance participation factor enables us toevaluate the participation of an apparatus in system oscillationsthrough using only black-box models. On top of this, thechain-rule yields parameter and state participation factorswhich enable us to look into the black-box and trace root-causes to detailed parameters and states without disclosing thestate equation. Based on this chain-rule of participation factors,we proposed the grey-box approach for power system stabilityanalysis, as will be described in the following subsection. Fig. 3. Illustration of the three-layer grey-box. In Layer-1, estimates of thepotential participants are created based on the upper bound of ∆ λ subjectto (cid:107) ∆ Z k (cid:107) = (cid:15) (cid:107) Z k (cid:107) . In Layer-2, the contribution of a participant to modedamping is estimated based on the real part of ∆ λ subject to ∆ Z k = (cid:15)Z k .In Layer-3, the root-cause of instability within the participating apparatusis identified, and parameter re-tuning facilitated, via the impedance-parametersensitivity. If the parameters are selected as values of the state-matrix, Layer-3has the same function as state-space participation analysis and is as transparentas a white-box. C. The Grey-Box Approach
The grey-box approach contains three layers with differenttransparencies according to the available prior-knowledge, asillustrated in Fig. 3. The higher the transparency, the moreprior-knowledge is needed but along with that comes moreuseful information for root-cause tracing and trouble-shootingin whole-system stability analysis. Now we describe in detaileach layer and the relationships between them.
1) Grey-Box Layer-1:
In this first layer, the only prior-knowledge available is the impedance Z k for all apparatusin the system, along with the nodal admittance matrix Y net determined by the topology and line impedances of the net-work. The impedance Z k can be provided either in the formof transfer functions Z k ( s ) or frequency spectra Z k ( jω ) . Fortransfer functions, the impedance participation factor p λ n ,Z k can be calculated directly, but for frequency spectra theimpedance participation factor needs to be estimated indirectly.Due to the three-wire, three-phase nature of power systems, Z k and p λ n ,Z k are × matrix blocks in the synchronous dq frame. We need a scalar index to represent the fourelements of p λ n ,Z k collectively so that different p λ n ,Z k canbe compared and the location of the dominant apparatus (fora given mode) determined. To this end, we assign a consistentmagnitude perturbation to each Z k and observe the effect onthe eigenvalue λ . The perturbation is normalized to (cid:107) Z k (cid:107) sothat it scales with the corresponding apparatus, that is, (cid:107) ∆ Z k (cid:107) = (cid:15) (cid:107) Z k (cid:107) (25)where (cid:15) is a small positive constant. According to the Cauchyinequality, we have | ∆ λ | = |(cid:104) p λ,Z k , ∆ Z k ( λ ) (cid:105)| ≤ (cid:107) p λ,Z k (cid:107) · (cid:107) ∆ Z k ( λ ) (cid:107) (26)4his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 5which yields | ∆ λ | max = (cid:107) p λ,Z k (cid:107) · (cid:107) ∆ Z k ( λ ) (cid:107) = (cid:15) (cid:107) p λ,Z k (cid:107) · (cid:107) Z k ( λ ) (cid:107) . (27)It is clear from (27) that (cid:107) p λ,Z k (cid:107) · (cid:107) Z k ( λ ) (cid:107) determines theupper bound of ∆ λ , | ∆ λ | max , which is the maximum possibleparticipation of the corresponding apparatus. Only apparatuswith relatively large | ∆ λ | max may possibly, but not necessarily,participate in the λ -mode. Thus, we use (cid:107) p λ,Z k (cid:107) · (cid:107) Z k ( λ ) (cid:107) asthe primary participation index in Layer-1 of the grey-boxapproach. Layer-1 roughly identifies potential participants ina mode and does so with very little prior-knowledge.
2) Grey-Box Layer-2:
Building on Layer-1, Layer-2 adds astipulation that the perturbation of apparatus impedance ∆ Z k is aligned to the original impedance Z k , that is, ∆ Z k = (cid:15) Z k (28)where (cid:15) is a very small positive real number. This is areasonable stipulation because it emulates the effect of discon-necting from the grid a portion of a power plant containingmultiple similar apparatus, e.g. disconnecting one turbine ina wind farm, and keeping the others in operation, in orderto investigate the impact of this plant on the system. This isin contrast to Layer-1 which emulates the effect of puttingpr removing an additional but not similar apparatus at thesame node thus inducing ∆ Z k with an arbitrary orientationnot aligned to Z k .Applying norm on (28) yields (cid:107) ∆ Z k (cid:107) = (cid:15) (cid:107) Z k (cid:107) which isconsistent with Layer-1, but the extra stipulation in Layer-2specifies the orientation of ∆ Z k and thus enables us to obtainthe detailed value of ∆ λ besides its upper bound (cid:107) ∆ λ (cid:107) max ∆ λ = (cid:104) p λ,Z k , ∆ Z k ( λ ) (cid:105) = (cid:15) (cid:104) p λ,Z k , Z k ( λ ) (cid:105) . (29)From (29) we see that ∆ λ is determined by (cid:104) p λ,Z k , Z k ( λ ) (cid:105) which we use as the new participation index for Layer-2.This index adds to Layer-1 extra information about eachparticipant’s contribution to the damping, positive or negative,and thus stability of the whole-system via the real part of thecorresponding ∆ λ .
3) Grey-Box Layer-3:
The final layer of the grey-boxapproach aims to look into the participating apparatus in orderto identify which physical component or control loop in theapparatus is the root-cause of instability, and thus provideinformation on which parameter should be re-tuned, and howto re-tune it, to stabilize the system.Layer-3 uses the sensitivity of an impedance to its inter-nal parameters, ∂Z k /∂ρ , so that a parameter perturbation ispropagated to an impedance perturbation and further into aneigenvalue perturbation via the chain-rule (see (21)-(23) inSection III-B). The chain-rule yields a parameter participationfactor that indicates which internal parameters can be re-tunedso that the eigenvalues move in the desired direction on thecomplex plane. If the internal parameters are selected as valuesin the state matrix, the corresponding parameter participationfactors are exactly the state participation factors. Importantly,the impedance-parameter sensitivity ∂Z k /∂ρ discloses littleinformation regarding the internal design or control of anapparatus and yet it enables root-cause tracing inside the Fig. 4. Modified IEEE 14-bus test system, with 3 extra IBRs connected atbus 12, 13 and 14. TABLE IP
ARAMETERS OF A PPARATUS (SG
S AND
IBR S ) Type Parameter Symbol ValueSynchronousGenerator Inertia J D X R f i
250 HzDC-link control bandwidth f dc
20 HzPLL bandwidth f PLL
20 HzFiltering Reactance X R C dc apparatus as effectively as a transparent white-box model. Thisis the great advantage of the Layer-3 grey-box approach.IV. C ASE S TUDY OF A C OMPOSITE P OWER S YSTEM
We now demonstrate the grey-box approach in a case studyof a composite power system with a high penetration of IBRs.The chosen system, Fig. 4, is based on the IEEE 14-bussystem [16] with three additional IBRs (Type-IV wind farms)connected to buses 11, 12, and 13. Parameters of the SGs andthe IBRs are listed in Table I. Each wind farm is aggregated asa single IBR with a scaled rating. To make the system prone tooscillation, the inertia and damping torque coefficient of M3, atbus-3, are deliberately de-tuned to 10% of normal values, andthe current control bandwidth of M7, at bus-12, is de-tuned tobe 140 Hz not 250 Hz. The data and codes used to generatethe simulation results can be found at: https://github.com/Future-Power-Networks/Power-System-Analysis-Toolbox.The whole-system admittance of the network, constructedfrom the impedances of all apparatus and admittances of allof the network lines, is displayed in the bode plot in Fig. 5.At each node, the whole-system admittance ˆ Y kk is a × matrix in the synchronous dq frame, and we only display oneof the four elements in the matrix since that is enough toillustrate the characteristics of the system. Only the nodes with5his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 6 Fig. 5. Bode diagram of whole system admittance ˆ Y kk at nodes with sources,presented in dd axis. sources (SGs or IBRs) connected are presented in the bode plotbecause other nodes are passive. A few resonant peaks appearon the bode plot, each representing an oscillation mode in thesystem. The mode at 60 Hz represents the flux dynamics ofwindings and lines which is well known [9]. We use the grey-box approach to trace the root-cause of the rest of the modesand find ways to damp the modes.Fig. 6 shows the results of applying grey-box Layer-1 (onthe left) and Layer-2 (on the right). For mode-1 (5.5 Hz), M3stands out in the Layer-1 pie chart. Further, the breakdowninto real and imaginary components for Layer-2 shows thatM3 affects both the damping and natural frequency of mode-1 as might be expected since M3 has low inertia and a lowdamping torque coefficient and is thus liable to create anunder-damped swing mode. The other modes in the dashedcircle at lower frequencies are similar to mode-1 but areparticipated by other SGs and therefore have lower magnitudeand frequency because of the high inertia. For mode-2, theLayer-1 pie chart indicates that M4 and M7 are the mainparticipants, and the Layer-2 decomposition shows that theSG (M3) and the IBR (M7) have different contributions tostability. The real part of ∆ λ corresponding to M7 is negative,and the imaginary part is positive. This implies that de-ratingof M7 (which is equivalent to increasing the impedance ofM7 by ∆ Z ) makes the eigenvalue move rightwards and thusbecome more stable. The role of M4 is opposite to that ofM7 since the components ∆ λ having opposite signs. We seethat M7 tends to destabilized the system whereas M4 helps tostabilize. This is an example of IBR-induced oscillation. Theexact cause of the destabilisation is not revealed until Layer-3of the grey-box which can point to particular components andcontrol loops. Nonetheless, Layer-1 and Layer-2 grey-boxesreveal rich information about the root-causes of modes 1 and2 with little prior-knowledge.After locating the participating apparatus and identifyingtheir role in system stability, the final step is to use the grey-box Layer-3 to re-tune the parameters in M3 and M7 to (a)(b)Fig. 6. Participation analysis of the two under-damped modes using Layer-1and Layer-2 of the grey-box. (a) Mode-1 (5.5 Hz), (b) Mode-2 (20.9 Hz).Fig. 7. Grey-box Layer-3: Parameter participation factors of M3 in mode-1and M7 in mode-2, showing the amplitude and direction of the eigenvaluemovement when subject to parameter variations in M3 and M7. stabilize mode-1 and mode-2 respectively. Based on the chain-rule, the parameter participation factors for mode-1 and mode-2 against the internal parameters of M3 and M7 are calculatedand illustrated in Fig. 7. It is clear that the damping ratio ofmode-1 is sensitive to the damping torque coefficient D , andthe natural frequency of mode-1 is sensitive to the reactance X and inertia J . This result agrees with the classic view of SGdynamics but is achieved here with the impedance-based grey-box model. For mode-2, it is clear that three control loops,i.e. the current control, dc-link voltage control, and phase-locked loop (PLL), have impact on this mode. Increasing thecurrent control bandwidth f i helps to stabilize the mode, butincreasing the dc-link control bandwidth f dc and the PLLbandwidth f PLL tend to destabilize the mode. Therefore, wecan remark that this mode is the result of coupling betweeninner-loops (current control) and outer-loops (dc-link controland PLL) in a relatively weak grid.6his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 7 (a)(b)(c)Fig. 8. Apparent power output of M3 and M7 during a 20% load stepchange at bus 4: (a) De-tuned system with obvious oscillations during transientprocess; (b) Increase f i _M7 by 30% and D M3 by 100%, giving significantimprovement in system stability; (c) Decrease f i _M7 by 30% leading to anunstable system. The major predictions of the three-layer grey-box analysisfor this case-study were verified via time-domain simulation.The apparent power output of M3 and M7 were measured. Thetransient response of the measured power subjected to a increase in load on bus-4 at t = 5 s is displayed in Fig. 8 (a).Two oscillation modes at 5.5 Hz (mode-1) and 20.9 Hz (mode-2) are visible in the power transient process of M3 and M7respectively, and these values agree with the resonant peaksin the admittance spectrum in Fig. 5. The parameters of M3and M7 were re-tuned to damp the modes according to thestabilization suggestions of the Layer-3 grey-box analysis inFig. 7, and the results are shown in Fig. 8 (b). In contrast,Fig. 8 (c) shows the results when the parameters are re-tunedin the opposite direction. It is clear that Fig. 8 (b) is properlydamped in both modes but Fig. 8 (c) becomes unstable.The grey-box-based participation analysis correctly locatedthe root-cause of oscillations and suggested an appropriateapproach for stabilization.V. C ONCLUSIONS
The grey-box approach establishes a systematic methodfor stability and participation analysis of complex powersystems with only impedance information. It has three layerswith different transparencies at each layer to facilitate root-cause tracing to different depths, i.e. apparatus, states, andparameters, according to the available knowledge. These grey-boxes provide a very useful tool to look inside a black-box(impedance) model to achieve almost the same transparencyas a white-box (state-space) model, but without the need formanufacturers of apparatus to disclose the internal details thatwould be required in the white-box approach. The proposed grey-box approach is based on rigorous mathematical analysiswith proof of the relationship between residues and impedanceparticipation factors, and elucidation of the chain-rule ofsensitivity propagation for internal states and parameters tobe carried forward to the impedance participation factor, thuspresenting a unified participation theory.A
PPENDIX AM ATHEMATICAL P RELIMINARIES
We summarise the mathematical preliminaries used in thispaper to assist the reader and to make the paper self-contained.
A. Residue
In complex analysis, the residue of a complex function G ( s ) is defined as the g − coefficient of the Laurent series of G ( s ) .This is, given the Laurent series of G ( s ) around λG ( s ) = ∞ (cid:88) h = −∞ g h · ( s − λ ) h (30)the residue of G at λ is defined asRes λ G = g − . (31)If λ is a non-repeated pole of G , the reside can be calculatedfrom Res λ G = lim s → λ ( s − λ ) G ( s ) . (32)This property is used in the proof of Lemma 1 in AppendixB. The residue can be applied element-wise on a matrix ofcomplex functions. B. Frobenius Inner Product
For two complex-valued matrices V and W with the samedimension, the Frobenius inner product of V and W is definedas (cid:104) V, W (cid:105) (cid:44) (cid:88) h,l V hl W hl (33)where h and l are the row and column indices of the matrices,and denotes complex conjugation. The complex conjugationin (33) ensures that the Frobenius inner product of a complexmatrix with itself is a non-negative real number, and thus isinduced the Frobenius norm (cid:107) · (cid:107)(cid:107) V (cid:107) (cid:44) (cid:112) (cid:104) V, V (cid:105) . (34)The Frobenius inner product and norm are derived from thecommon inner product in vector spaces, so the properties ofthe common inner product are naturally inherited. One of themost useful properties is the Cauchy inequality |(cid:104) V, W (cid:105)| ≤ (cid:107) V (cid:107) · (cid:107) W (cid:107) (35)where the equality holds if and only if V and W are alignedin orientation. This property is used in (26) in Section III-C.7his work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after whichthis version may no longer be accessible. 8A PPENDIX BP ROOF OF L EMMA G ρ is a scalar transferfunction and the pole λ is a zero of H ρ = G − ρ , that is, H ρ ( λ ) = 0 . (36)A perturbation on ρ induce a corresponding perturbation on λ , that is H ρ +∆ ρ ( λ + ∆ λ ) = 0 . (37)Since H ρ is analytic around its zero λ , we have the followingfirst-order Taylor expansion of (37) H ρ +∆ ρ ( λ ) + H (cid:48) ρ +∆ ρ ( λ )∆ λ = 0 (38)in which H (cid:48) represents the derivative of H . Combining (36)-(38) yields H ρ +∆ ρ ( λ ) − H ρ ( λ )+ H (cid:48) ρ ( λ )∆ λ + (cid:0) H (cid:48) ρ +∆ ρ ( λ ) − H (cid:48) ρ ( λ ) (cid:1) ∆ λ = 0 (39)and equivalently ∆ H ρ ( λ ) + H (cid:48) ρ ( λ )∆ λ + ∆ H (cid:48) ρ ( λ )∆ λ = 0 . (40)Suppressing the high-order infinitesimal ∆ H (cid:48) ρ ( λ )∆ λ in (40)yields ∆ λ = − H (cid:48) ρ ( λ ) − · ∆ H ρ ( λ ) . (41)As λ is a non-repeated pole for G ρ , the residue of G ρ at λ isgiven byRes λ G ρ = lim s → λ ( s − λ ) G ρ ( s ) = lim s → λ s − λH ρ ( s ) = 1 H (cid:48) ρ ( λ ) (42)in which the second equal sign results from L’Hôpital’s rule.Combining (41) and (42) yields ∆ λ = − Res λ G ρ · ∆ H ρ ( λ ) . (43)This is the reduced case of Lemma 1 with G ρ being a scalartransfer function.Now we prove the case where G ρ is a square matrix. Insuch a case, the pole λ is a zero for the determinant of H ρ ,that is det ( H ρ ( λ )) = 0 . (44)We take det ( H ρ ) (cid:44) H det as a scalar transfer function so asimilar result to (41) is obtained ∆ λ = − H (cid:48) det ( λ ) − ∆ H det ( λ ) . (45)Expanding H det along a column yields H det = (cid:88) h H ρhl F ρhl (46)in which F ρ is the cofactor matrix for H ρ and the subscript hl denotes the element in a matrix at the h -th row and l -thcolumn. It is clear to see from (46) that ∂H det ∂H ρhl = F ρhl (47) and hence ∆ H det = (cid:88) h,l ∂H det ∂H ρhl ∆ H ρhl = (cid:88) h,l F ρhl ∆ H ρhl = (cid:104) F ρ , ∆ H ρ (cid:105) (48)where the complex conjugate is associated to the Frobeniusinner product (cid:104)· , ·(cid:105) for complex matrices defined in (33).Since G ρ is now a matrix, its residue needs to be calculatedelement-wiseRes λ G ρ = lim s → λ ( s − λ ) G ρ ( s )= lim s → λ (cid:0) ( s − λ ) H ρ ( s ) − (cid:1) = lim s → λ (cid:18) s − λH det ( s ) F ρ ( s ) (cid:62) (cid:19) = F ρ ( λ ) (cid:62) H (cid:48) det ( λ ) (49)where we make use of the fact that H ρ ( s ) − = F ρ ( s ) (cid:62) /H det ( s ) . (50)Combining (45), (48) and (49) yields Lemma 1 ∆ λ = − H (cid:48) det ( λ ) − (cid:104) F ρ ( λ ) , ∆ H ρ ( λ ) (cid:105) = − (cid:42) F ρ ( λ ) H (cid:48) det ( λ ) , ∆ H ρ ( λ ) (cid:43) = (cid:104)− Res λ G ρ (cid:62) , ∆ H ρ ( λ ) (cid:105) = (cid:104)− Res ∗ λ G ρ , ∆ H ρ ( λ ) (cid:105) . (51)R EFERENCES[1] M. Amin and M. Molinas, “Small-signal stability assessment of powerelectronics based power systems: A discussion of impedance-andeigenvalue-based methods,”
IEEE Transactions on Industry Applica-tions , vol. 53, no. 5, pp. 5014–5030, 2017.[2] P. Kundur et al. , “Definition and classification of power system stabilityieee/cigre joint task force on stability terms and definitions,”
IEEEtransactions on Power Systems , vol. 19, no. 3, pp. 1387–1401, 2004.[3] F. Blaabjerg and K. Ma, “Wind energy systems,”
Proceedings of theIEEE , vol. 105, no. 11, pp. 2116–2131, 2017.[4] J. Bialek, “What does the GB power outage on 9 august 2019 tell usabout the current state of decarbonised power systems?”
Energy Policy ,vol. 146, p. 111821, 2020.[5] Y. Gu, J. Liu, T. C. Green, W. Li, and X. He, “Motion-inductioncompensation to mitigate sub-synchronous oscillation in wind farms,”
IEEE Transactions on Sustainable Energy , vol. 11, no. 3, pp. 1247–1256, 2019.[6] J. Shair, X. Xie, L. Wang, W. Liu, J. He, and H. Liu, “Overview ofemerging subsynchronous oscillations in practical wind power systems,”
Renewable and Sustainable Energy Reviews , vol. 99, pp. 159–168, 2019.[7] Y. Li, Y. Gu, Y. Zhu, A. Junyent-Ferré, X. Xiang, and T. C. Green,“Impedance circuit model of grid-forming inverter: Visualizing controlalgorithms as circuit elements,”
IEEE Transactions on Power Electron-ics , vol. 36, no. 3, pp. 3377–3395, 2021.[8] X. Wang, L. Harnefors, and F. Blaabjerg, “Unified impedance modelof grid-connected voltage-source converters,”
IEEE Trans. on PowerElectron. , vol. 33, no. 2, pp. 1775–1787, 2018.[9] P. Kundur,
Power system stability and control . McGraw-hill New York,1994, vol. 7.[10] J. Rommes and N. Martins, “Computing large-scale system eigenvaluesmost sensitive to parameter changes, with applications to power systemsmall-signal stability,”
IEEE Transactions on Power Systems , vol. 23,no. 2, pp. 434–442, May 2008.[11] K. W. Wang, C. Y. Chung, C. T. Tse, and K. M. Tsang, “Multimachineeigenvalue sensitivities of power system parameters,”
IEEE Transactionson Power Systems , vol. 15, no. 2, pp. 741–747, May 2000.[12] Y. Wang, X. Wang, F. Blaabjerg, and Z. Chen, “Harmonic instabilityassessment using state-space modeling and participation analysis ininverter-fed power systems,”
IEEE Transactions on Industrial Electron-ics , vol. 64, no. 1, pp. 806–816, Jan 2017. [13] Y. Gu, Y. Li, Y. Zhu, and T. Green, “Impedance-based whole-systemmodeling for a composite grid via embedding of frame dynamics,”
IEEETransactions on Power Systems , 2020.[14] E. Ebrahimzadeh, F. Blaabjerg, X. Wang, and C. L. Bak, “Bus par-ticipation factor analysis for harmonic instability in power electronicsbased power systems,”
IEEE Transactions on Power Electronics , vol. 33,no. 12, pp. 10 341–10 351, Dec 2018.[15] Y. Zhan, X. Xie, H. Liu, H. Liu, and Y. Li, “Frequency-domainmodal analysis of the oscillatory stability of power systems with high-penetration renewables,”
IEEE Transactions on Sustainable Energy ,vol. 10, no. 3, pp. 1534–1543, 2019.[16] R. C. Iraj Dabbagchi, “14 bus power flow test case,” Aug. 1993.[Online]. Available: https://labs.ece.uw.edu/pstca/pf14/pg_tca14bus.htm
Yue Zhu (S’21) received the B.Eng and M.Sc de-grees in electrical engineering, Zhejiang University,Hangzhou, China, in 2016 and 2019 respectively. Heis currently a Ph.D student at Department of Elec-trical and Electronic Engineering, Imperial CollegeLondon. His present research focuses on impedance-based stability analysis of power systems, and thenoise evaluation for impedance measurement.
Yunjie Gu (M’18-SM’20) received the B.Sc. andthe Ph.D. degree in Electrical Engineering fromZhejiang University, Hangzhou, China, in 2010 and2015 respectively. He was a Consulting Engineer atGeneral Electric Global Research Centre, Shanghai,from 2015 to 2016, and is now an EPSRC-funded In-novation Fellow at Imperial College London (awardEP/S000909/1). His research interests include powersystem control and stability, and the application ofpower electronics to power systems.
Yitong Li (S’17-M’21) received the B.Eng degreesin electrical engineering from Huazhong Universityof Science and Technology, China, and the Uni-versity of Birmingham, UK, in 2015. He receivedhis M.Sc degree in Future Power Networks fromImperial College London, UK, in 2016, where heis currently pursuing his Ph.D. degree. His currentresearch interests include control of power electronicconverters and analysis of power system dynamics.
Timothy C. Green (M’89-SM’02-F’19) receiveda B.Sc. (Eng) (first class honours) from ImperialCollege London, UK in 1986 and a Ph.D. fromHeriot-Watt University, Edinburgh, UK in 1990. Heis a Professor of Electrical Power Engineering atImperial College London, and Director of the EnergyFutures Lab with a role fostering interdisciplinaryenergy research. His research interest is in using theflexibility of power electronics to accommodate newgeneration patterns and new forms of load, such asEV charging, as part of the emerging smart grid.In HVDC he has contributed converter designs that reduce losses whilealso providing control functions assist AC system integration. In distributionsystems, he has pioneered the use of soft open points and the study of stabilityof grid connected inverters. Prof. Green is a Chartered Engineering the UKand a Fellow of the Royal Academy of Engineering.(M’89-SM’02-F’19) receiveda B.Sc. (Eng) (first class honours) from ImperialCollege London, UK in 1986 and a Ph.D. fromHeriot-Watt University, Edinburgh, UK in 1990. Heis a Professor of Electrical Power Engineering atImperial College London, and Director of the EnergyFutures Lab with a role fostering interdisciplinaryenergy research. His research interest is in using theflexibility of power electronics to accommodate newgeneration patterns and new forms of load, such asEV charging, as part of the emerging smart grid.In HVDC he has contributed converter designs that reduce losses whilealso providing control functions assist AC system integration. In distributionsystems, he has pioneered the use of soft open points and the study of stabilityof grid connected inverters. Prof. Green is a Chartered Engineering the UKand a Fellow of the Royal Academy of Engineering.