Electromagnetic Oscillation Origin Location in Multiple-Inverter-Based Power Systems Using Components Impedance Frequency Responses

Existing impedance-based stability criterion (IBSC) for electromagnetic stability assessment of multiple-grid-connected-inverter (GCI)-based power systems suffers from several limitations. First, global stability feature is hard to be obtained if Nyquist-criterion-based IBSC is used. Second, heavy computational burdens caused by either right-half-plane (RHP) poles calculation of impedance ratios or nodal admittance matrix construction can be involved. Third, it's not easy to locate the oscillation origin, since the dynamics of individual components are missing in the aggregated load and source sub-modules. This article aims to overcome the aforementioned three limitations of the existing IBSC. First, frequency responses of the load impedance and source admittance defined at each node in a selected components aggregation path are obtained by aggregating individual components (e.g., GCIs and transmission lines), from which imaginary parts of RHP poles of these load impedances and source admittances are directly identified without knowing analytical expressions of these load impedances and source admittances. Then, based on the Nyquist plots of minor loop gains (defined as the ratios of the impedance frequency responses of these load and source sub-modules), stability features of these selected nodes are obtained. Finally, if some nodes are unstable, the oscillation origin is located based on numbers of the RHP poles of these load impedances and source admittances. Compared to the existing IBSC, the presented method can assess global stability and locate oscillation origin more efficiently. The local circulating current issue, as a main obstacle of the existing IBSC, can also be identified. Time-domain simulation results in Matlab/Simulink platform and real-time verification results in OPAL-RT platform of a four-GCI-based radial power plant validate the effectiveness of the presented electromagnetic oscillation origin location method.

The admittance and impedance frequency responses by seeing leftward and rightward at the right side of node N K . Y rl The admittance and impedance models by seeing leftward and rightward at the right side of node N K , and their frequency resposnes. Closed MLG between Z 2p−1 and GCI #p at node N p .

I. INTRODUCTION
Recently, renewable energies, such as wind power and solar power, have been increasingly penetrating into the existing utility grid [1]. Voltage source grid-connected inverters (GCIs), as efficient power electronic interfaces, are widely used to transmit the generated electricity into the utility grid [2]. However, impedance interactions between various control loops of the GCIs (e.g., outer power control loop, dc-link voltage control loop, inner current control loop, and phase-locked loop (PLL)) and passive components (e.g., transmission lines (TLs), underground cables, and transformers) may cause electromagnetic oscillation phenomena in various frequency ranges [3]- [6]. Though oscillation mechanisms of single-GCI-based power plants have been thoroughly explored in [4], [7]- [11], the oscillation mechanism and oscillation origin location of multiple-GCI-based power plants need to be further investigated.
Impedance-based stability criterion (IBSC) was originally proposed in [12] to assess the stability of DC power systems, and then modified in [13], [14] to cope with AC power systems. The whole power plant is partitioned into load and source sub-modules at an arbitrary node, based on which impedance ratio of the load and source sub-modules is calculated. Both encirclement number of Nyquist plot of the impedance ratio around (−1, j0) and number of right-halfplane (RHP) poles of the impedance ratio should be calculated [14], [15]. For single-GCI-based power plants, RHP poles calculation can be avoided, since both the GCI and the grid should be inherently stable [10], [14]. Furthermore, if the system is assessed as unstable, the oscillation origin can be located as the impedance interaction between the GCI and the grid. However, avoidance of RHP poles calculation and location of oscillation origin commonly cannot be achieved in multiple-GCI-based power plants by performing the Nyquist criterion (NC) one time at a selected node, since RHP poles may emerge during the components aggregation procedure, and the oscillation origin may be lost in the aggregated load and source sub-modules [16]- [18].
The RHP poles of the aggregated load and source submodules can be analytically derived based on system topology and impedance transfer functions of individual components [16], [17], [19]. However, this analytical approach sometimes cannot be implemented in practice, since internal structures and parameters of the GCIs can be confidential due to industry secrecy. In addition, the orders of the load and source sub-modules may be very high if a huge number of GCIs exist in the power plant, which will result in heavy computational burdens. The vector fitting (VF) and matrix fitting (MF) algorithms initially proposed in [20], [21] can fit a set of discrete frequency responses by continuous transfer function and transfer matrix, respectively, which have been employed in [18], [22]- [25] to facilitate the stability assessment of power electronic-dominated power systems. Specifically, the VF algorithm is used in [18], [23] to fit discrete load impedance and source admittance frequency responses as continuous transfer functions, so that extraction of the RHP poles from the analytically-derived transfer functions can be avoided. Discrete loop impedance models (LIMs) are established in [26]- [28] by further aggregating the impedance frequency responses of the load and source sub-modules, based on which the MF algorithm is used to generate continuous LIMs, and zeros of the determinant of the fitted LIMs are identified. Although the RHP poles identification can be more efficient or even avoided with the help of the VF and MF algorithms, oscillation origin still cannot be located.
According to [17], [29], nodal admittance model (NAM) and the aforementioned LIM are dual pairs. Therefore, the zeros of the determinants of system NAMs are calculated in [30]- [32] for stability analysis. Similar to the LIM method, the RHP poles calculation can also be avoided in the NAM method. In addition, a frequency-domain component connection method is presented in [16], [33]- [35], where the generalized NC is applied on the return ratio of the impedance matrices of the connection network and the composite model of all inverters. Since the two impedance matrices can commonly be guaranteed to not have RHP poles, the RHP poles calculation in the generalized NC can be avoided. According to the discussions in [16], [35], the derived NAMs in [30]- [32] are exactly the closed-loop transfer function matrices of multiple-input-multiple-output negative feedback systems with feed-forward gain being 1 and feed-back gain being the return-ratio matrices derived in [16], [33]- [35]. Since admittance information of each GCI is stored in the derived NAM and return-ratio matrix, global stability feature can be obtained, and oscillation origin can further be located based on participation factor analysis [35]. However, orders of both the NAM and the return-ratio matrix can be high, which may lead to heavy computational burdens. In addition, once GCIs are added/removed into/from the power plant or system topology changes, the high-order NAM and return-ratio matrix should be re-calculated.
The NC and generalized NC are sequentially implemented from the farthest node to point of common coupling (PCC) of radial power plants in [36], [37] and [19], [38], [39], respectively. No RHP poles calculation is needed, since both load and source sub-modules are inherently stable in each step. In addition, the node where the (generalized) NC is not satisfied is identified as the oscillation origin. Compared with the LIM, NAM, and component connection method, the multiple-step stability analysis method is more flexible, and can reduce the computational burdens. However, only stability of grid current instead of output currents of individual GCIs can be assessed, i.e., this analysis method cannot obtain global stability feature. Furthermore, the presented conditions in [19], [36]- [39] to guarantee the stability of grid current are sufficient yet not necessary.
To fill in the research gaps mentioned in the above literature review, a computationally efficient global stability analysis and oscillation origin location method based on only terminal impedance frequency responses of individual components is presented in this article. First, terminal impedance frequency responses of individual components are aggregated at both sides of each node in a selected components aggregation path. On its basis, RHP poles of these load and source sub-modules are directly identified from their impedance and admittance frequency responses, respectively. The stability of these selected nodes is then obtained with further help of the Nyquist plots of the impedance ratios. If a specific node is unstable, the oscillation origin is located based on numbers of RHP poles of these load and source sub-modules. Main contributions of this article can be highlighted as follows.
1) An NC-based sufficient and necessary condition for global stability is derived, where relation between global stability and local stability is discussed. 2) Based on the derived NC-based global stability condition, an efficient oscillation origin location method is presented. The problematic nodes where RHP poles occur during components aggregation can be identified. 3) RHP poles of load and source sub-modules for the NC are directly identified from load impedances and source admittances frequency responses, respectively, without knowing analytical expressions of load and source submodules, which can reduce the computational burdens and cope with black-box models. 4) Circulating current phenomenon is identified based on numbers of RHP poles of the load and source submodules along the selected components aggregation path. 5) A grid current stability enhancement method is initiated, i.e., the paralleled stable branch of an unstable branch is enforced to be the same as the unstable branch, which may be achieved by adjusting the transmission line impedance of the stable branch. The rest of this article is organized as follows. In Section II, an NC-based sufficient and necessary condition for global stability of a representative radial power plant is derived. On  its basis, the principle and implementation procedure of the presented global stability analysis and oscillation origin location method are explained in Section III. In Section IV, the presented method is implemented in a four-GCI-based radial power plant. The correctness of the theoretical analysis results in Section IV is validated in Section V. In Section VI, comparisons between the presented method and existing IBSCs are performed. Finally, Section VII draws the conclusion. Fig. 1 shows the single-line diagram of a radial gridconnected power plant, which consists of M GCIs. The detailed control structure of these GCIs are explained in Appendix A. If grid current I g and output currents of all GCIs (i.e., I 1 -I M ) are stable, the system is globally stable. The global stability can be assessed by further simplifying the equivalent circuit model in Fig. 1. The aggregation rules of two components/sub-modules during the circuit simplification procedure are stated here [18]. Each component/sub-module can be defined as either Y type or Z type which should be modeled as an Norton or a Thevenin equivalent circuit, respectively. In detail, individual GCIs and grid are defined as Y type and Z type, respectively. The transmission lines and underground cables are defined as Y type in parallel case and Z type in series case. Based on the definitions of individual components, the type of the aggregated sub-module of two components/sub-modules are further defined as follows. The aggregated sub-module is Z type, if two Z-type components/sub-modules are connected in parallel or in series, or one Z-type component/sub-module and one Y -type component/sub-module are connected in parallel. Furthermore, the aggregated sub-module is Y type, if two Y -type components/sub-modules are connected in parallel, or one Ztype component/sub-module and one Y -type component/submodule are connected in series.

II. NC-BASED SUFFICIENT AND NECESSARY CONDITION FOR GLOBAL STABILITY
Based on the components/sub-modules aggregation rules, Fig. 1 can be simplified as Fig. 2(a) by representing the left and right parts of PCC as an Norton and a Thevenin equivalent circuit, respectively. In addition, Fig. 1 can be simplified as Fig. 2(b) by representing the left and right parts of node N K (K ∈ [1, M]) as an Norton and a Thevenin equivalent circuit, respectively. Based on Fig. 2, I g and I K can be calculated as and respectively, where T clm1 Based on Appendix B, numbers of RHP poles of I g and I K can be calculated as and  [14]. However, Y le f t PCC and Z right N K may have RHP poles (This issue will be discussed in Section III-A.) [15].
A sufficient and necessary condition for the global stability of Fig. 1 can be summarized as follows.
Lemma 1: The radial grid-connected power plant in Fig. 1 is globally stable if and only if 1) I g is stable, i.e., P(T m1

III. PRESENTED GLOBAL STABILITY ANALYSIS AND OSCILLATION ORIGIN LOCATION METHOD
In this section, the NC-based sufficient and necessary condition for global stability derived in Section II is reformulated. On its basis, the principle of the presented global stability analysis and oscillation origin location method is explained.

A. REFORMULATION OF THE NC-BASED SUFFICIENT AND NECESSARY CONDITION FOR GLOBAL STABILITY
Detailed expressions of the aggregated Norton equivalent circuit in Fig. 2(a) and the aggregated Thevenin equivalent circuit in Fig. 2(b) should be derived to calculate the RHP poles of I L PCC , Y L PCC , Z R N1k , and V R N1k (k = 1, 2, . . ., M), which can be obtained by aggregating the components in Fig. 1(b) repeatedly.
Based on the components/sub-modules aggregation rules in Section II, the components aggregation procedure of Fig and respectively, where T clm Fig. 1 can then be derived based on Fig. 3(a) and (5)-(7) using recursive method, shown as Fig. 2(a) can then be calculated as I le f t N 0 and Y le f t N 0 , respectively. A sufficient yet not necessary condition for stability of I g can be derived based on item 1) of lemma 1 and (8), shown as follows.
and V lr respectively, where T clml Fig. 1 can then be derived based on Fig. 3(b), (5), and (9)-(10) using recursive method, shown as Based on the derived Norton and Thevenin equivalent circuits of the left and right parts of node N K in (8) and (12), respectively, Fig. 1 can be simplified as Fig. 4 at node N K , which can further be simplified as Fig. 2(b). The parameters of S right N K can thus be calculated as where T clmrl where T mrl A sufficient yet not necessary condition for stability of I K can be derived based on item 2) of lemma 1, (8), (12), and (14), shown as follows.
Lemma 3: To not lose the generality, Fig. 5(a) shows a complicated grid-connected power plant with multiple strings. Similar with the stability analysis of I g and I K (K ∈ [1, M]) in Fig. 1, which needs to calculate Y le f t PCC in Fig. 2(a) and Z right Fig. 2(b), respectively, the stability analysis of grid current and output currents of all GCIs in Fig. 5(a) also needs the corresponding aggregated admittance and impedance frequency responses, respectively. Take the stability analysis of I g and I 12 as an example. First, a main components aggregation path is selected, shown as the red line in Fig. 5(a), based on which the rest of the whole system can be divided into strings (i.e., strings #1, #3, #5, #6, and #7) and individual GCIs (i.e., GCIs #2, #4, and #8). Then, the equivalent circuit model of Fig. 5(a) can be established in Fig. 5(b) by applying the components can be calculated using Fig. 4. It can be seen that the presented components aggregation procedures in Fig. 3 are the basis for stability analysis of more complicated grid-connected power plants.

B. PRESENTED GLOBAL STABILITY ANALYSIS AND OSCILLATION ORIGIN LOCATION METHOD
The converse-negative propositions of lemmas 2 and 3 can be derived as lemmas 4 and 5, respectively. and Y le f t N K are needed for calculation of P(I K ), respectively. The derivation procedure of these impedance/admittance models in the conventional IBSCs can be explained as follows [14], [17], [19], [23]  According to the Kirchhoff's Law, I K in Fig. 4 can be expressed as which indicates that stability of I K can be revealed by stability of I (Z 2K−2 ) and I (Z 2K ), i.e., implementation point of the NC can be moved from node N K in the existing IBSCs to the right side of node N K . In this case, only Y rl Fig. 4 should be calculated without further calculating Z right N K , which can release computational burdens. On its basis, Fig. 6 shows the flowchart of the presented stability analysis and oscillation origin location method at node N K (K ∈ [1, M]).
In step 1, impedance frequency responses of individual components are obtained, which can be achieved using the following methods. First, the system planner can theoretically derive the impedance models of the GCIs if internal control structures and parameters are provided by the vendors [3], [8], [10], [16], [17], [19]. Second, when detailed information of the GCIs is confidential for the system planner due to the industry secrecy and intellectual property protection, the vendors tend to deliver lookup tables which provide discrete impedance/admittance frequency responses of the GCIs [18], [22], [23]. Third, when black-box simulation models are provided by the manufactures, the system planner can perform frequency scanning in commercial softwares, e.g., PSCAD/EMTDC and Matlab/Simulink [26], [27], [39], [41]. For example, the impedance frequency responses of the GCIs obtained by frequency scanning are plotted in Fig. 8(b). Finally, if both discrete impedance frequency responses and black-box simulation models are not provided by the vendors, the system planner can experimentally measure the impedance frequency responses of the GCIs by successively connecting these GCIs with a configurable strong grid [42]- [44]. It is not necessary to equip a lot of measurement devices in the network.
In step 2, frequency responses of Z rr , are calculated according to the components aggregation procedures in Figs. 3(a) and 3(b), respectively. On its basis, the stability of I (Z 2K−2 ) is assessed by , respectively. In detail, the magnitude-frequency curve's one peak whose phase angle is not within [−90 o , 90 o ] corresponds to a pair of RHP poles. If P(T clm N K ) is zero, I (Z 2K−2 ) is stable, and no oscillation origin should be located. Otherwise, step 3 should be performed to locate the oscillation origin. Similar with lemmas 4 and 5, the following lemma can be derived from (16).  1, j0) should be counted. In detail, if the encirclement number is zero, the oscillation phenomenon only results from the inherent instability of the two sub-modules. Otherwise, the oscillation phenomenon results from both unstable impedance interaction and inherent instability. The origin of the inherent instability can further be located based on Fig. 7.
Based on (5)-(7) and item 1) of lemma 6, Fig. 7(a) shows the flowchart of origin location of RHP poles of Y f rrl and Z 2K are aggregated.
Based on (5), (9), (10), and item 2) of lemma 6, P(Z f rrr N K ) can be calculated as According to (17), Fig. 7  using the recursive method. The RHP poles of the load and source sub-modules are commonly identified from the theoretically-derived impedance/admittance formulas, which may bring in heavy computational burdens if a large number of GCIs exist or even cannot be achieved if internal information of the GCIs is confidential [15], [16], [19], [27]. Although the RHP poles can be extracted from the impedance/admittance frequency responses using the VF and MF algorithms, selection of the optimal order is still challenging [18], [22], [23]. On the other hand, step 2 of the presented method directly identifies the RHP poles from the impedance/admittance frequency responses, which indicates that computational burden can be released, and the black-box issue can be solved. In addition, different from the conventional NC-based stability analysis methods which use the RHP poles of the load and source sub-modules only to assess the stability at the checking nodes, the RHP poles of the load and source sub-modules are used in step 3 of the presented method to also locate the oscillation origin.

IV. IMPLEMENTATION OF THE PRESENTED OSCILLATION ORIGIN LOCATION METHOD
In this section, the presented global stability analysis and oscillation origin location method is implemented in a four-GCI-based radial power plant shown in Fig. 8(a), where global stability (case 1), local instability (cases 2-4), and global instability (cases 5-8) are investigated. Circuit and controller parameters of the four GCIs are listed in Table 3 in Appendix A. In addition, lengths of the transmission lines in the eight cases are listed in Table 1. The per-unit-length (p.u.l.) resistance and inductance of theses transmission lines are 10 μ /km and 10 μH/km, respectively.
According to step 1 of the proposed method in Fig. 6, the impedance frequency responses of individual GCIs in Fig. 8(a) should be obtained. 200 logarithmically-distributed frequency points between 400 Hz and 5 kHz are measured using the frequency scanning method (The detailed implementation procedure can be found in [39]). The Bode diagrams of theoretical impedance frequency responses derived by (21) and measured data are plotted as Z GCI theo and Z GCI mea , respectively, in Fig. 8(b), which shows that they highly agree with each other. Furthermore, phase angle of Z GCI theo is lower than −90 o in frequency range (1423 Hz, 1667 Hz). When length of the transmission line is between 7 km and 30 km, the magnitudes interaction point falls into the non-passive region, and corresponding phase angles difference is higher than 180 o , which indicates that the single GCI is unstable when length of the connected transmission line is between 7 km and 30 km.

A. CASE 1: GLOBAL STABILITY
In case 1, lengths of the transmission lines in Fig. 8(a) are l (Z 0 + Z g ) = 50 km, l (Z 1 ) = 2.5 km, l (Z 2 ) = 2 km, l (Z 3 ) = 2 km, l (Z 4 ) = 1 km, l (Z 5 ) = 1.5 km, and l (Z 6 ) = 1 km. Bode diagrams of the four load impedances and four source admittances, i.e., Z f rrr Fig. 8(a), are plotted in Fig. 9(a), which shows an unstable peak appears at 1633 Hz of the magnitude-frequency curve of Y f rrl has a pair of RHP poles whose imaginary parts are ±2π 1633 rad/s [18]. In addition, Nyquist diagrams of T m N K (K ∈ [1,4]) calculated by (16) at the four checking points, i.e., the right sides of N 4 -N 1 , are plotted in Fig. 9(a), which shows Nyquist plot of T m N 1 encircles (−1, j0) two times in anticlockwise direction. Information of the Bode and Nyquist diagrams is summarized in the table of  FIGURE 9. Implementation of the presented global stability analysis and oscillation origin location method in Fig. 8(a) for cases 1-4 Fig. 9(a), where Z and Y mean load and source sub-modules, respectively. According to (16), all the four checking points are stable, i.e., I (Z 6 ), I (Z 4 ), I (Z 2 ), and I g are stable. The system is globally stable, and no oscillation origin needs to be located.

. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
Besides the oscillation origin location method presented in this article, problematic components can also be identified using conventional eigenvalue-based method and participation factor analysis. However, weak nodes with low stability margins cannot be determined intuitively. On the other hand, the Nyquist plots drawn at different nodes in this article can easily tell the system planner the stability margins of these nodes, which can provide guidelines for further stability improvement. For example, the Nyquist plots in Fig. 9(a) show that nodes N 1 and N 2 are the weakest and strongest in nodes N 4 -N 1 , respectively, since the Nyquist plots of T m N 1 and T m N 2 are closest and farthest to (−1, j0).

B. CASES 2-4: CIRCULATING CURRENT-RELATED LOCAL INSTABILITY
In case 2, lengths of the transmission lines are l (Z 0 + Z g ) = 50 km, l (Z 1 ) = 2.5 km, l (Z 2 ) = 2 km, l (Z 3 ) = 2 km, l (Z 4 ) = 50 km, l (Z 5 ) = 20 km, and l (Z 6 ) = 20 km. Fig. 9(b) shows an unstable peak appears at 1498 Hz of the magnitudefrequency curve of Y f rrl N 3 , which indicates it has a pair of  Fig. 19(a) RHP poles whose imaginary parts are ±2π 1498 rad/s. In addition, Fig. 9(b) shows T m N 4 and T m N 3 encircle (−1, j0) two times in clockwise and anticlockwise directions, respectively. The analysis results summarized in the table of Fig. 9(b) show that the right sides of nodes N 3 -N 1 and node N 4 are stable and unstable, respectively. Therefore, I (Z 4 ), I (Z 2 ), and I g are stable, while I (Z 6 ) oscillates at one frequency, which means circulating current exists between GCIs #4 and #3. It can be derived from step 3 of Fig. 6 and Fig. 7 that aggregations between GCI #4 and Z 6 , and between GCI #3 and Z 5 bring in two RHP poles in Y le f t In case 3, lengths of the transmission lines are l (Z 0 + Z g ) = 100 km, l (Z 1 ) = 2.5 km, l (Z 2 ) = 50 km, l (Z 3 ) = 20 km, l (Z 4 ) = 10 km, l (Z 5 ) = 0 km, and l (Z 6 ) = 0 km, which indicates nodes N 3 , N 3 , N 4 , and N 4 are overlapped. Fig. 9(c) shows two unstable peaks appear at 1553 Hz and 1497 Hz of the magnitude-frequency curves of Z f rrr N 4 and Y f rrl N 2 , respectively, which indicates they have two pairs of RHP poles whose imaginary parts are ±2π 1553 rad/s and ±2π 1497 rad/s, respectively. In addition, Fig. 9(c) shows T m N 3 and T m N 2 encircle (−1, j0) two times in clockwise and anticlockwise directions, respectively. The analysis results summarized in the table of Fig. 9(c) show that the right sides of nodes N 2 -N 1 and nodes N 4 -N 3 are stable and unstable, respectively. Therefore, I (Z 2 ) and I g are stable, while I (Z 6 ) and I (Z 4 ) oscillate at one frequency, which means circulating current exits among GCIs #4-#2. It can be derived from step 3 of Fig. 6  In case 4, lengths of the transmission lines are l (Z 0 + Z g ) = 50 km, l (Z 1 ) = 21 km, l (Z 2 ) = 7 km, l (Z 3 ) = 0 km, l (Z 4 ) = 0 km, l (Z 5 ) = 0 km, and l (Z 6 ) = 0 km, which indicates nodes N 2 -N 4 and N 2 -N 4 are overlapped. Fig. 9(d) shows three unstable peaks appear at 1511 Hz, 1567 Hz, and 1488 Hz of the magnitude-frequency curves of Z f rrr , and Y f rrl N 1 , respectively, which indicates they have three pairs of RHP poles whose imaginary parts are ±2π 1511 rad/s, ±2π 1567 rad/s, and ±2π 1488 rad/s, respectively. In addition, Fig. 9(d) shows T m N 2 and T m N 1 encircle (−1, j0) two times in clockwise and anticlockwise directions, respectively. The analysis results summarized in the table of Fig. 9(d) show that the right sides of node N 1 and nodes N 4 -N 2 are stable and unstable, respectively. Therefore, I g is stable, while I (Z 6 ), I (Z 4 ), and I (Z 2 ) oscillate at one frequency, which means circulating current exists among GCIs #4-#1. It can be derived from step 3 of Fig. 6 and Fig. 7  are totally same. GCIs #4-#2 as a whole oscillate with GCI #1, i.e., sum of harmonic contents of I 4 -I 2 is equal to harmonic content of I 1 .

C. CASES 5-8: GLOBAL INSTABILITY
In case 5, l (Z 4 ) is reduced from 50 km in case 2 to 20 km. Compared with Fig. 9(b), Fig. 10(a) shows one additional unstable peak appears at 1455 Hz of the magnitude-frequency curve of Z f rrr . Therefore, I (Z 4 ), I (Z 2 ), and I g oscillate at one frequency, while I (Z 6 ) oscillates at two frequencies, which means a circulating current exists between GCIs #4 and #3, and a harmonic current flows through all nodes. Based on case 2, it can be concluded that the harmonic current flowing through all nodes originates from the unstable impedance interaction at the right side of node N 3 .
In case 6, l (Z 6 ) is increased from 20 km in case 5 to 25 km. Compared with Fig. 10(a), Fig. 10 . Therefore, I (Z 6 ), I (Z 4 ), I (Z 2 ), and I g oscillate at two frequencies. It can be concluded that, the discrepancy between l (Z 6 ) and l (Z 5 ) enforces the circulating current between GCIs #4 and #3 in case 5 to flow into the rest of the power plant.
In case 7, l (Z 6 ) is increased from 20 km in case 2 to 25 km. Compared with Fig. 9(b), Fig. 10 . Therefore, I (Z 6 ), I (Z 4 ), I (Z 2 ), and I g oscillate at one frequency. It can be concluded that, the discrepancy between l (Z 6 ) and l (Z 5 ) enforces the circulating current between GCIs #4 and #3 in case 2 to flow into the rest of the power plant.
In case 8, l (Z 0 + Z g ) is reduced from 100 km in case 3 to 50 km. Compared with Fig. 9(c), Fig. 10 . Therefore, I (Z 2 ) and I g oscillate at one frequency, while I (Z 6 ) and I (Z 4 ) oscillate at two frequencies, which means a circulating current exists between GCIs #4-#2, and a harmonic current flows through all nodes. Based on case 3, it can be concluded that the harmonic current flowing through all nodes originates from the unstable impedance interaction at the right side of node N 1 .
The Bode and Nyquist diagrams in Figs. 9 and 10 cannot identify oscillation frequencies of the unstable cases. On the other hand, Bode diagrams of the LIMs calculated at the four checking points for cases 1-8 are plotted in Fig. 11, based on which oscillation frequencies of I (Z 6 ), I (Z 4 ), I (Z 2 ), and I g are summarized in Table 1.

V. SIMULATION AND REAL-TIME VERIFICATION
In this section, the correctness of the global stability analysis and oscillation origin location results in Section IV based on the presented method is validated by simulation results in Matlab/Simulink platform and real-time verification results in OPAL-RT platform. Fig. 12 shows time-domain simulation results of output currents of the four GCIs and grid current in Fig. 8(a) for cases 1-5. Fig. 12(a) shows simulation results when case 1 changes to case 2 at 3 s, which indicates the system is globally stable before 3 s and locally unstable after 3 s. In detail, the 1500-Hz frequency components of both I 4 and I 3 in case 2 are 244 A, which indicates this harmonic current circulates between GCIs #4 and #3. The analysis results in Figs. 9(a), 9(b), 11(a), 11(b), and Table 2 are thus validated. Fig. 12(b) shows simulation results of case 3, which indicates the system is locally unstable. In detail, the 1500-Hz frequency components of both I 4 and I 3 are 682 A and half of that of I 2 , which indicates this harmonic current circulates between GCI #2 and the sub-module composed of GCIs #4    Table 1 are thus validated. Fig. 12(c) shows simulation results of case 4, which indicates the system is locally unstable. In detail, the 1491-Hz frequency components of I 4 -I 2 are 572 A and one third of that of I 1 , which indicates this harmonic current circulates between GCI #1 and the sub-module which consists of GCIs #4-#2. The analysis results in Figs. 9(d), 11(d), and Table 1 are thus validated. Figs. 12(a)-(c) show grid current is stable even output currents of some GCIs are unstable, which indicates stability of grid current can be enforced by installing a parallel branch which is totally the same as the unstable branch. Fig. 12(d) shows simulation results of case 5, which indicates the system is globally unstable. In detail, the 1500-Hz frequency components of both I 4 and I 3 are 1065 A, which indicates this harmonic current circulates between GCIs #4 and #3. Furthermore, the 1439-Hz frequency components satisfy that F 1439 (I 4 + I 3 + I g ) = F 1439 (I 2 + I 1 ), which indicates this harmonic current circulates between the sub-module composed of GCIs #4, #3 and grid and the sub-module composed of GCIs #2 and #1 (F indicates the FFT operator.). The analysis results in Figs. 10(a), 11(e), and Table 1 are thus validated. Fig. 13 shows picture of the real-time OPAL-RT digital simulator platform in laboratory, where OP5600 hardware platform and OPAL-RT's RT-LAB software platform are combined to validate the model developed in Matlab/Simulink platform in a real-time manner. The real-time verification results are then post-processed in Matlab software. Fig. 14(a) shows real-time verification results of case 6, which indicates the system is globally unstable. In detail, the 1433-Hz frequency components of I 4 -I 1 and I g satisfy that F 1433 (I 4 + I 3 + I g ) = F 1433 (I 2 + I 1 ), which indicates this harmonic current circulates between the sub-module composed of GCIs #4, #3, and grid and the sub-module composed of GCIs #2 and #1. In addition, the 1484-Hz frequency components of I 4 -I 1 and I g satisfy that F 1484 (I 4 + I 2 + I 1 ) = F 1484 (I 3 + I g ), which indicates this harmonic current circulates between the sub-module composed of GCIs #4, #2, and #1 and the sub-module composed GCIs #3 and grid. The analysis results in Figs. 10(b), 11(f), and Table 1 are thus validated. Since the 1484-Hz harmonic current oscillates much faster than the 1433-Hz harmonic current, only the frequency contents at initial oscillation stage from 28.5 s to 28.7 s are extracted. Fig. 14(b) shows real-time verification results of case 7, which indicates the system is globally unstable. In detail, the 1481-Hz frequency components of I 4 -I 1 and I g satisfy that F 1481 (I 4 + I 2 + I 1 ) = F 1481 (I 3 + I g ), which indicates this harmonic current circulates between the sub-module composed of GCIs #4, #2, and #1 and the sub-module composed of GCI #3 and grid. The analysis results in Figs. 10(c), 11(g), and Table 1 are thus validated. Fig. 14(c) shows real-time verification results of case 8, which indicates the system is globally unstable. In detail, the 1441-Hz harmonic currents satisfy that F 1441 (I 4 + I 3 + I 2 + I g ) = F 1441 (I 1 ), which indicates this harmonic current circulates between the sub-module composed of GCIs #4, #3, #2, and grid and the sub-module composed of GCI #1. Furthermore, the 1504-Hz frequency components of both I 4 and I 3 are 4.679 A and half of that of I 2 , which indicates this harmonic current circulates between GCI #2 and the submodule composed of GCIs #4 and #3. The analysis results in Figs. 10(d), 11(h), and Table 1 are thus validated.

VI. DISCUSSIONS
In this section, the oscillation origin location method presented in this article is compared with the existing NC-based IBSCs. In addition, the computational burdens of the presented method when network complexity and GCIs number increase are investigated. Some other issues related with the presented method are also discussed.

A. COMPARATIVE ANALYSIS OF SEVERAL BREEDS OF NC-BASED IBSCS
Based on the components aggregation procedures in Fig. 3, several impedance aggregation methods have been reported in the existing literature, which can mainly be grouped into the following five breeds [13], [16]- [19], [23], [36]- [39], [45]- [49]. In IBSC #1, the NC is performed one time at a certain node, e.g., PCC, with calculating RHP poles of the minor loop gain (MLG) and encirclement number of its Nyquist plot around (−1, j0), where the RHP poles of the MLG are identified from the theoretically-derived formulas of the load impedance and source admittance models [13], [16], [19], [45]- [47]. Different from IBSC #1, to better cope with the black-box models, the RHP poles of the MLG are identified in IBSC #2 from impedance frequency responses (IFRs) of load and source sub-modules based on the VF algorithm [18], [23]. Furthermore, to avoid calculating the RHP poles of the MLG, stability of grid current is assessed in IBSC #3 by performing the NC sequentially from farthest node to PCC [19], [36]- [39]. The grid current is assessed as stable if the Nyquist plots at all the checking nodes do not encircle (−1, j0). However, only local stability, i.e., stability of grid current, can be obtained in IBSCs #1, #2, and #3. In IBSC #4, the NC is performed at all nodes to assess stability of grid current and output currents of all GCIs with calculating RHP poles of the MLGs and encirclement numbers of their Nyquist plots around (−1, j0), where the RHP poles of the MLGs are identified from the theoretically-derived formulas of the load impedance and source admittance models [17], [48], [49]. In IBSC #5, the NC is performed sequentially at all nodes to assess stability of grid current and output currents of  Fig. 8(a).
all GCIs without calculating RHP poles of the MLGs. Similar with IBSC #3, the system is assessed as globally stable if the Nyquist plots at all the checking nodes do not encircle (−1, j0).
Performances of the aforementioned five breeds of IBSCs for Fig. 1 are summarized in Table 2. It can be seen that only local stability feature can be obtained by IBSCs #1, #2, and #3. Global stability feature can be obtained by IBSCs #4 and #5. However, white-box models are required in IBSC #4, and oscillation origin cannot be located. In addition, the stability analysis result obtained by IBSC #5 is conservative, and the implementation procedure is relatively complicated. Take case 2 of Fig. 8(a) as an example. When IBSCs #1-#3 are used at PCC, although stability of I g can be successfully assessed as stable, the circulating current-related local instability of I 4 and I 3 is ignored. Furthermore, Fig. 15 shows the analysis results based on IBSC #5. The stability of I g can be assessed based on T m N 4 -T m N 1 and T mr Although T m N 4 , T m N 3 , and T mr N 3 encircle (−1, j0) two times in clockwise, clockwise, and anticlockwise directions, respectively, I g in Fig. 12(a) is stable. Similar analysis can be performed for I 2 and I 1 . Therefore, IBSC #5 is a sufficient yet not necessary condition for global stability.
In the presented IBSC, the NC is performed sequentially at all nodes of the checking path to assess stability of grid current and output currents of all GCIs with calculating RHP poles of the MLGs and encirclement numbers of their Nyquist plots around (−1, j0), where RHP poles of the MLGs are directly identified from the IFRs. The theoretical analysis in Section IV and simulation verification in Section V show that, compared with the aforementioned five breeds of IBSCs, the presented IBSC provides a sufficient and necessary condition for global stability and can locate the oscillation origin only based on IFRs of the black-box models.

B. VARIATION OF COMPUTATIONAL BURDENS OF THE PRESENTED METHOD WITH NETWORK COMPLEXITY
To validate effectiveness of the presented oscillation origin location method for more complicated network structure, the radial power plant in Fig. 1 with eight GCIs (i.e.,  Fig. 16 show that the right sides of nodes N 7 -N 1 and node N 8 are stable and unstable, respectively. Therefore, I (Z 12 ), I (Z 10 ), I (Z 8 ), I (Z 6 ), I (Z 4 ), I (Z 2 ), and I g are stable, while I (Z 14 ) oscillates at one frequency, which means circulating current exists between GCIs #8 and #7. It can be derived from step 3 of Fig. 6 and Fig. 7  are totally same. Fig. 17 shows the Matlab/Simulink-based simulation verification results of Fig. 16, which indicates I 8 and I 7 are unstable, while I (Z 12 ) and I g are stable (I 6 -I 1 are also stable, and their simulation results are omitted for simplicity.). The simulation results indicate that one harmonic current circulates between GCIs #8 and #7. The correctness of the analysis results in Fig. 16 is thus validated.
Furthermore, Fig. 18 shows the variations of the running times of the presented IBSC and IBSC #5 when number of the GCIs in Fig. 1 increases from 1 to 1000 with step size 10 (Intel  i7-7600 U processor and 8 GB of RAM). It can be seen that the running time of the presented IBSC increases linearly as number of checking nodes increases, and the average execution time for each checking node is 23.5 s/1000 = 23.5 ms. On the other hand, the average execution time for each checking node using IBSC #5 is 105.3 s/1000 = 105.3 ms, which is about four times slower than the presented IBSC.

C. OTHER ISSUES RELATED WITH THE PRESENTED METHOD
It is possible to further identify the problematic controller or even controller parameters based on the identified problematic components. The discrete impedance frequency responses of the GCIs are fitted as continuous transfer functions in [18], [41], and [43] using vector fitting algorithm, system identification technique, and least-squares fitting routine, respectively, based on which the circuit and controller parameters of the GCIs are identified by equalizing the fitted transfer functions and theoretical impedance models. Note that internal structure should be known for further circuit and controller parameters identification.
This paper aims to present a global stability analysis method to identify the oscillation origin in the multiple-GCIbased power systems, where the effects of various control loops on system stability are not focused. Therefore, the digital time delay-related high-frequency stability issue, as a representative, is used to verify the validity of the proposed method in this paper. The DC-link voltage dynamics are, thus, neglected by assuming a large DC-link capacitor. In fact, similar simplification has been widely adopted in existing works. For example, the dynamics of DC-link voltage control and PLL are neglected in [31], [40], where high-frequency harmonic instability is focused. In addition, the DC-link voltage is assumed as constant in [45] to focus on the PLL-related low-frequency instability issue. The detailed impedance modeling procedure of PMSG considering DC-link voltage dynamics can be found in [50].
When outer controller dynamics cannot be ignored, the presented three-step-based oscillation origin location method in Fig. 6 can be modified for further low-frequency stability analysis. In step 1, the MIMO instead of SISO impedance/admittance frequency responses of GCIs are measured, as shown in the authors' recent work [39]. On its basis, the MIMO impedance/admittance frequency responses of the aggregated load and source sub-modules are derived still using Fig. 3, except that the SISO impedance/admittance models in (5)-(7) and (9)-(10) are updated as the corresponding MIMO impedance/admittance models. In step 2, similar with the identification method in this article, the RHP poles of the load and source sub-modules are directly identified from their MIMO impedance and admittance frequency responses, respectively. In addition, the generalized NC is applied on the ratio of the MIMO impedance frequency responses of the load and source sub-modules. In step 3, the oscillation origin location based on the MIMO impedance/admittance frequency responses can be achieved in a similar way as Fig. 7, since the RHP poles generation mechanism during components aggregation procedure using the SISO and MIMO impedance models are the same.

VII. CONCLUSION
This article presents a global stability analysis and oscillation origin location method of multiple-GCI-based power systems. Based on terminal impedance frequency responses of individual components, frequency responses of equivalent load and source sub-modules are calculated at all nodes along the selected components aggregation path, based on which RHP poles are directly identified from Bode diagrams of the load and source frequency responses without knowing analytical expressions. Stability features of these nodes and oscillation origin can further be obtained based on the improved IBSC. Effectiveness of the presented oscillation origin location method is validated by theoretical analysis, simulation results in Matlab/Simulink platform, and real-time verification in OPAL-RT platform. The presented method can assess global stability and locate oscillation origin location more efficiently. Furthermore, local instability issues induced by impedance/admittance order reduction phenomena when paralleled branches are totally the same and unstable can also be identified. Inspired by this, one future work can be to improve stability of grid current by enforcing the paralleled branch to be the same as the existing unstable branch. Fig. 19(a) shows the control strategy of an LCL-filtered GCI under grid current control mode with PLL synchronization. The circuit and controller parameters of the GCI are listed in Table 3. The impedance/admittance model of the GCI is an MIMO matrix when considering PLL dynamics [3], [10], [16]. However, PLL only shapes low-frequency impedance/admittance characteristics of the GCI [3], [10]. In addition, the influence of PLL can be ignored if its bandwidth is narrow enough. In this case, d-axis and q-axis impedance components are decoupled, and d-axis impedance component suffices for stability analysis. In fact, similar simplification has been widely adopted in existing work [4], [14], [18], [22], [23], [30], [31], [37], [40]. The d-axis control block diagram can be derived as Fig. 19(b), where G i = K p + K i s is the current controller, and G del = e −1.5T s s represents the digital time delay. In addition, G x1 and G x2 are expressed as G x1 = 1 L f 1 C f s 2 + C f K cp G del s + 1 (18) and G x2 = L f 1 C f s 2 + C f K cp G del s + 1 L f 1 L f 2 C f s 3 + L f 2 C f K cp G del s 2 + (L f 1 + L f 2 )s ,

APPENDIX A RELATION BETWEEN RHP POLES OF CURRENT SOURCE AND ADMITTANCE OF NORTON EQUIVALENT CIRCUIT OF GCIS
respectively. On its basis, the closed-loop transfer function G cl and output admittance Y o can be calculated as (20) and respectively. Based on (21), I g,d can be derived as Therefore, the Norton equivalent circuit of the GCI can be established as Fig. 19(c). It can be seen from (20) and (21) that G cl = G i G del G x1 Y o . Since G i , G del , and G x1 are designed to not have RHP poles, the current source G cl I re f g and admittance Y o in Fig. 19(c) share the same stability characteristics. There- Since D 2 (I le f t PCC ) does not have RHP zeros, P(I g ) can be calculated as Since D 2 (I le f t N K ) and D 2 (V right N K ) do not have RHP zeros, P(I K ) can be calculated as Therefore, (4) is proved. Assume that Z right and Y bot ) . According to (10), Z lr can be calculated as Z lr Therefore, P(Z lr N K−1 ) can be calculated as P(Z lr Therefore, (14) is proved. It can be concluded from (24), (26), and (28) that I * g , I * K , and Z right N K−1 need not be considered for the calculation of P(I g ), P(I K ), and P(Z lr N K−1 ), respectively, due to the poles cancellation phenomena, as shown in (23), (25), and (27), respectively.