Small-Signal Stability of Multi-Converter Infeed Power Grids with Symmetry

: Traditional power systems have been gradually shifting to power-electronic-based ones, with more power electronic devices (including converters) incorporated recently. Faced with much more complicated dynamics, it is a great challenge to uncover its physical mechanisms for system stability and/or instability (oscillation). In this paper, we ﬁrst establish a nonlinear model of a multi-converter power system within the DC-link voltage timescale, from the ﬁrst principle. Then, we obtain a linearized model with the associated characteristic matrix, whose eigenvalues determine the system stability, and ﬁnally get independent subsystems by using symmetry approximation conditions under the assumptions that all converters’ parameters and their susceptance to the inﬁnite bus ( B g ) are identical. Based on these mathematical analyses, we ﬁnd that the whole system can be decomposed into several equivalent single-converter systems and its small-signal stability is solely determined by a simple converter system connected to an inﬁnite bus under the same susceptance B g . These results of large-scale multi-converter analysis help to understand the power-electronic-based power system dynamics, such as renewable energy integration. As well, they are expected to stimulate broad interests among researchers in the ﬁelds of network dynamics theory and applications.


Introduction
The power system is generally believed to be one of the most complicated systems made by human beings, featured with some essential characteristics of complex systems, such as nonlinearity, large spatial scale, multiple timescale, interwinding between coupled elements, etc. [1][2][3]. In a traditional power system dominated by synchronous generators (SGs), the system of electromechanical dynamics is basically controlled by the rotor motion of SG under imbalance torque (or power) [4][5][6]. Its physical picture is clearly understood. Recently the traditional synchronous-generator-dominant power system has been gradually becoming a power-electronic-based one. In particular, with increasing penetration of renewable energy, a large number of devices, such as wind turbines and photovoltaic cells, have been connected to the power grid via power electronic converters, such as voltage source converter (VSC) as a major type of converter [7][8][9][10][11][12][13]. Different from the SG, the converter exhibits a greater complication, due to its diversified structures of the controller. Entangled interaction of controllers within and between converters also makes the system analysis more difficult [14]. Recently, a variety of wide-frequency-band oscillations of power-electronic-based power system have occurred frequently worldwide [2]. A relevant stability analysis has become a hot topic in the fields of both power electronics and power system [1][2][3][14][15][16][17]. Hence it is important and urgent to study the dynamics in a multiconverter system and uncover its stability mechanism.
Although the power-electronic-based power system is intrinsically nonlinear, for simplicity, the classical linear system theory has been extensively applied in the framework of small-signal stability under the condition that the operating point is subjected to a sufficiently small disturbance [4][5][6]. So far, three major methods have been developed, including the modal analysis method in the time domain, and complex torque and impedance methods in the frequency domain [18]. We will apply the modal analysis method as it is common and reliable. In addition, the multi-timescale decomposition technique has also been widely used in reduced-oder modeling and analysis of power-electronic-based power systems. There are several studies concentrating only on one single timescale. For instance, the slow-scale dynamics of VSC was studied within the DC-link voltage timescale by ignoring all fast-scale current controllers [19]. In contrast, it was reported that the fast oscillations can be put into the category of current-timescale by setting all slow-scale controller outputs as constants [20]. We will use the same strategy in our modeling.
On the other hand, the power system stability, as an important problem in realistic complex systems, has also attracted wide interests in physics in the past decade [21][22][23][24][25][26][27][28][29][30][31][32][33][34]. Most studies concentrated on cascade failure [26,27] and synchronization [28][29][30], based on the well-known second-order Kuromoto oscillators and mainly treating the problem how does the network structure shape power grid function. There are few papers on renewable integrations, except some on intermittent behaviors of renewable energies by studying stochastic differential equations [31]. Basically, their timescale is much slower than what we are interested in here. For recent progress on this direction, see, e.g., the focus issues [32][33][34]. Detailed studies on how renewable integration devices work still lack, to the best knowledge of the authors. We would like to stimulate the broad interests of researchers in physics.
Different from a single-converter system, the topologies of a multi-converter system can be diversified. However, due to the comparatively low energy density transferred by a converter, a symmetrical structure has been commonly observed, such as AC microgrid with distributed power sources [35], radial cluster topology of off-shore wind farms [36], grid-connected photovoltaic systems with a radial configuration [37], and multi-infeed VSC power systems [38], etc. The above systems have two common features. One is that the grid-connected devices have similar control structures and parameters, and the other is that the network topology has a certain symmetry. Taking the offshore wind farm as an example, the common topologies are as shown in Figure 1 [39]. Those features bring great convenience to the stability analysis of the multi-converter system. The main conclusions of our work will be based on such property of symmetry.
Based on the above literature review, we find that a generic model of a multi-converter infeed system hasn't been built yet and the stability mechanism research is still urgently needed. In this paper, relying on the above three major techniques: linearization, multitimescale decomposition and symmetry analysis, we will establish a nonlinear model of a multi-converter power system within the DC-link voltage timescale from the first principle. Then we will derive a linearized model and conduct linear stability analysis. Finally, by using the mentioned symmetry conditions, we will decompose the multi-converter system into several single-converter subsystems and identify the dominant one for the system's small-signal stability. We hope those step-by-step modeling and mathematical analysis help to uncover the internal mechanism of the stability of multi-converter infeed system and inspiring for future studies of complicated realistic systems. The offshore wind farm with wind turbines connecting each other is known as a collection system and then it is connected to the onshore grid usually by a high voltage transmission system. The collection system topologies can be radical connection (a), radical loop connection (b), and star connection (c).
The remainder of this paper will be arranged as follows. In Section 2, we establish a nonlinear model of multi-converter infeed power system within the DC-link voltage timescale. In Section 3, the nonlinear model is linearized around the system operation conditions. In Section 4, to obtain the eigenvalues of the characteristic matrix, we make some symmetry assumptions and decouple the whole system into several subsystems. Further analysis finds that only one single-VSC system is dominant in system stability. Section 5 is devoted to simulation results of root loci and time series. Finally, conclusions and discussions are given.In addition, more details about the converter parameters, approximate eigenvalues and the simulation data are represented in the Appendixes A-C. Figure 2a shows the representative and simplified topological structure for a multi-VSC infeed power system within the DC-link voltage control timescale, with each VSC connected to an infinitely strong bus U t,n+1 by an inductance. As only slow-scale dynamics are considered, all voltage and current can be expressed as phasors. The system has n + 1 buses, including n VSC buses and one infinite bus. The infinite bus is numbered as n + 1 and adopted as the reference bus. Both the voltage amplitude U t,n+1 (simplified as U g ) and the phase θ t,n+1 of the infinite bus are fixed (θ t,n+1 = 0).

System Overview
For the i-th (i = 1, 2, · · · , n) VSC bus, the voltage amplitude and phase are denoted by U t,i and θ t,i , respectively. In our model, for the reactance, X = ω 0 L, where ω 0 is the fundamental frequency and ω 0 = 100π, as f 0 = 50 Hz. As shown in Figure 2a,b, X gi denotes the reactance between the i-th VSC and the infinite bus, termed as the convertergrid reactance, X gi = X ig . In contrast, X ij represents the reactance between the i-th and the j-th VSCs, referred to as inter-converter reactance, X ij = X ji . Figure 3 exhibits the detailed control diagram of a multi-VSC power system by utilizing a typical vector control scheme [40][41][42]. The control loops include the terminal voltage control loop (TVC) and DC voltage control loop (DVC) for the outer slow-scale control, and the AC current control loop (ACC) for the inner fast-scale control. In addition, the phase-locked control loop (PLL) for system synchronization belongs to the slow timescale. Different from the xy synchronous reference coordinate rotating at a fixed synchronous speed ω 0 , the cascaded control loops including the TVC, DVC and ACC should operate on the dq coordinate based on the information of the phase-locked angle θ i from the PLL. The coordinates and angle relations are illustrated in Figure 4a. The negative feedback controls of the TVC and DVC aim to maintain the terminal voltage and the DC voltage at normal values, respectively. The proportional-integral (PI) linear controllers are adopted in DVC, TVC, and PLL. After the cascade control, the output of the ACC should be transferred back from the dq coordinate to the xy one, serving as gating voltage for six insulated gate bipolar transistors (IGBTs). The waveform change is possible by using the well-known pulse-width modulation (PWM) technique; for more details, see, e.g., [40][41][42].

VSC Modeling
U t,n+1 Figure 2. (a) Topological structure for a multi-voltage source converter (VSC) infeed system within the DC-link voltage timescale, with each VSC connected to an infinitely strong bus U t,n+1 . We use X ij and X ig to denote the inter-converter and converter-grid reactances, respectively. X ij = X ji and X ig = X gi . (b) Equivalent circuit topology of (a), where each VSC serves as a controlled current source and injects power to the network.  Figure 3. Control diagram of a multi-VSC power system with typical cascade vector controllers, including the terminal voltage control loop (TVC) and DC voltage control loop (DVC) for the outer slow-scale controls, and the AC current control loop (ACC) for the inner fast-scale control. The phase-locked control loop (PLL) is used for synchronization, by inputting the terminal voltage angle θ t,i and outputting the phase-locked angle θ i . The PLL also belongs to the DVC timescale and will be included in our model. All controllers including the TVC, DVC, and ACC are designed depending on the PLL reference frame. Figure 4. (a) Schematic show for the xy synchronous reference coordinate and the dq PLL coordinate used in the multi-converter system. The xy rotates at a fixed speed ω 0 . In contrast, the dq coordinate is based on the phase-locked angle θ i (from the terminal voltage phase angle θ t,i ) and implements the cascade vector control. Its speed ω pll,i is time-varying. (b) Similar to (a) but for a synchronous generator. Similarly, the xy frame rotates at a fixed speed ω 0 and the dq frame rotates at the rotor speed ω i , which is also time-varying. δ i denotes the power angle and θ t,i is the terminal voltage phase angle. In a traditional power system, the dynamics of δ i of the synchronous generator (SG) plays a key role in the system electromechanical behavior, similar to the phase-locked angle θ i in the PLL in (a).
Clearly the PLL, as a watchdog, plays an important role in the whole control. It provides the phase-locked angle θ i , keeps the converter synchronized with the external grid, and makes the cascade vector controls possible. Figure 4a schematically shows the dynamical relation in the xy synchronous reference coordinate and the dq coordinate, where the phase-locked angle of the PLL, θ i , determines its d-axis orient. In addition, θ t,i denotes the terminal voltage angle, and u td,i and u tq,i are the dand q-axis components of the terminal voltage, respectively. Finally, to make a comparison, we present some coordinates and angle relations in Figure 4b, but for the synchronous generator instead [4][5][6]. Similarly, the xy frame rotates at a fixed speed ω 0 and the dq frame for the rotor position rotates at the rotor speed ω i , which is also time-varying. Generally, δ i is termed as power angle and θ t,i is the terminal voltage angle. The dynamics of δ i in SG plays a key role in the system dynamical behavior and is similar to the phase-locked angle θ i in PLL.
According to the control design criterion, the dynamics of VSC can be roughly divided into the DC-link voltage and current control timescales [7]. For simplicity, in this paper we will only study the DC-link voltage timescale dynamics, by keeping all the outer loops (including the TVC, DVC, and PLL) and neglecting the inner ACC. Hence the current values, i d,i and i q,i , should simultaneously follow the current reference values, i dre f ,i and i qre f ,i . The simplified control diagram is shown in Figure 5. In addition, we can keep the active power input P in,i from the DC-link capacitor in Figure 3 as a constant. As a result, each VSC can be regarded as equivalent to a controlled current source [15], as schematically shown in the right of Figure 5. For the symbols in Figure 5, (k p1,i , k i1,i ), (k p2,i , k i2,i ), and (k p3,i , k i3,i ) represent the PI control parameters of the DVC, TVC, and PLL, respectively. C i stands for the DC capacitor. U t,i and U tre f ,i are the terminal voltage and the terminal reference voltage, respectively. P in,i denotes the active power input from the DC side and P i denotes the active power output via the VSC. θ i represents the PLL angle, while θ t,i is the terminal voltage angle. A fifth-order dynamic model of a single VSC has been developed by us recently [9], which will be developed for a detailed modeling of a multi-VSC system. Below we will start from each controller dynamics, network description and establish a whole system model gradually.
For the PI controller in the DVC, the output of the integrator of PI 1 is denoted as x 1,i . Let x 1,i be a state variable. Then the DVC can be represented by where the dot represents the time derivative, along with an algebraic equation, For the dynamics of the TVC, similarly let x 2,i be the output of the integrator of PI 2 , as a state variable. The TVC dynamics can be described by and where U t,i is the amplitude of the terminal voltage from the network, The PLL is composed of not only the PI 3 controller, but also another integrator 1/s, as shown in Figure 5. The output of PI 3 is the speed of the dq coordinate, ω pll,i , and it is integrated to obtain the phase of the PLL coordinate, θ i . Again we separate the PI 3 into one integrator part and the other proportional part, and denote its integrator output as x 3,i .
We choose x 3,i and θ i as state variables and yield the dynamical second-order differential equation of the PLL, From Figure 4a, we get the so-called phase-locked angle error: Note that in the steady-state, θ t,i = θ i and u tq,i = 0. Next for the DC capacitor, its dynamics is determined by the imbalance active power between the input and output of the VSC, i.e., where P in,i is assumed to be a constant in our model and P i denotes the output active power, determined by Finally, we can combine the dynamical Equations (1), (3), (6), and (9), and get the following fifth-order differential equations for each VSC: Clearly x 1,i , x 2,i , x 3,i , θ i , and U dc,i are five state variables, marked red in Figure 5, where x 1,i , x 2,i , and x 3,i are the outputs of the integrators of the corresponding PI controllers. The input variables from the network are u td,i and u tq,i . In contrast, as the VSC can be regarded as a controlled current, i d,i in (2) and i q,i in (4) are set as the VSC outputs. In addition, substituting (5), (7), and (10) into (3), (4), (6), and (8), U t,i , θ t,i , and P i as intermediate variables can be removed and further represented by u td and u tq .
In summary, the dynamic model for each VSC is composed of the state-space Equation (11), along with the input Equations (5), (7), and (10) and the output Equations (2) and (4), showing highly interwoven nonlinear relations between these variables.

Network Modeling
For simplicity, we only consider inductance and neglect resistance and conductance. Based on the network topology in Figure 2 and applying the electric circuit theory [43], the mutual-susceptance between the device buses i and j, B ij (or called inter-converter susceptance), and the self-susceptance of the i-th device bus, B ii , can be expressed as, respectively, and Similarly, the mutual-susceptance between the i-th VSC bus and the infinite bus, B gi (or called converter-grid susceptance), and the self-susceptance of the infinite bus B gg can be expressed as, respectively, and With the voltage and current phasors, the bus admittance equation of the network [43] where We can further derive the network admittance equation in the xy synchronous coordinate, as shown in (19).
As mentioned before, the i-th VSC accepts voltage as inputs (u td,i and u tq,i ) from the network, and meanwhile provides current as outputs (i d,i and i q,i ) to the network. The coordinate transformations of these voltage and current between the dq coordinate and the xy coordinate are needed; based on Figure 4, we get From (16), we obtain Combining the above (17), (18) and (19) yields where (21) denotes the phase-locked angle mismatch between i and j. Note that we will always use bold to represent a matrix in the paper.

Whole System Modeling
Clearly Equations (11) and (20), along with the input Equations (5), (7), and (10) and the output Equations (2) and (4) as the interface between VSCs and network, have constituted a nonlinear model for the whole multi-VSC infeed system within the DVC timescale. In this picture, each node serves as a controlled current source, whereas the network works as a reservoir for the interaction. In addition, we like to emphasize that differing from the usual treatment in power system analysis in textbooks [4], where the whole system is analyzed within the common xy synchronous coordinate, here oppositely we keep each VSC within its own PLL coordinate and deal with the interaction on the network within the dq coordinate of each VSC. Actually, it keeps the converter inner structure and makes system modeling and analysis more convenient.

Linear Modeling for Multi-Converter Infeed System
Now we like to step further for a small perturbation analysis, in the fashion of power electrical engineering [4]. Therefore, the above nonlinear equations will be linearized around the operating point (or called fixed point in the nonlinear dynamics field) with the linear system theory employed.

Operating Point
By setting the right side of (11) equal zero, the operating point can be solved: U t0,i = U tre f ,i , θ i0 = θ t0,i , and thus u d0,i = U t0,i and u q0,i = 0. U dc0,i = U dcre f ,i , P i0 = P in,i . We will always use the subscript 0 to represent the operating point under the steady-state. Let U t0,i = U dc0,i = 1 in the per unit system, to simplify the calculation. In addition, U g = 1 for the infinite bus. In this way, the output power P i0 and the terminal voltage amplitude U t0,i of each VSC can be obtained. By treating these VSC buses as PV buses, the system steady-state can be analyzed by solving the so-called power flow [4][5][6]. As a result, we obtain the current at the operating point i d0,i and i q0,i , and the steady-state phase-lock angle θ i0 . All parameters are listed in Appendix A.

VSC Linearization
The state variable of VSC is denoted by The state-space equations for VSCs can be inferred by linearizing (11). Here we use the prefix ∆ to represent a small-signal variable. The output of VSC (also the network input) is denoted by For convenience of matrix calculation later, here we also add ∆θ i in the input variables of VSC, As a result, by linearizing and combining (5), (7), (10), and (11), we get the small-signal state-space equations of VSC in the matrix form where Note that the last column of H i is all zero, indicating that the ∆θ i in ∆U i has no any effect according to (25).
Correspondingly, the VSC output matrix can be obtained by linearizing and combining (2), (4), and (5), where Note also that among the four matrices F i , H i , J i , and K i , only H i depends on the steady-state values of i d0,i and i q0,i , whereas all other matrices depend on fixed system parameters only.

Network Linearization
The network Equation (20) should also be linearized. As we know, there are n + 1 buses with three input or output variables for each bus, as shown in (23) and (24). The dimension of the linear nodal admittance matrix should be 3(n + 1) × 3(n + 1). However, as the voltage amplitude and phase on the infinite (n + 1)-th bus are constant (U g = 1 and θ t,n+1 = 0), their perturbed variables can be eliminated completely. We obtain the final linear nodal admittance matrix Y c with the dimension of 3n × 3n: where Y c can be divided into the n × n partitioned matrix with the diagonal block of Y c expressed as and the off-diagonal blocks of Y c (i = j) expressed as Here θ ij0 = θ i0 − θ j0 denotes the phase-locked angle difference between i and j under the steady-state.

Characteristic Equation for Linear System Modeling
Finally, by combining (25), (27), and (29), and eliminating ∆U i and ∆I i , we get the linearized model, where the small-signal stability should be determined by A; A is the characteristic matrix with the dimension of 5n × 5n, with F, H, J, and K being all n-block diagonal partitioned matrices, Now we have established the complete state-space model of the n-VSC system after knowing all system parameters and the steady-state operating points through power flow calculation. We get the unified mathematical form of the characteristic matrix which can be directly extended into the system of any number of VSCs. If all eigenvalues of A are located in the left-half-plane, the system would be locally stable. In this way, the eigenvalue calculation of a huge matrix A should be conducted. Below, however, we like to step ahead for a possible order reduction and search for a physical result.

Small-Signal Stability Analysis for a Single VSC System
First, we start from the simplest case of a single VSC infeed system in Figure 6, where the VSC-grid reactance X g represents the reactance between the VSC and the infinite grid. Under this situation, the characteristic matrix can be easily obtained from (34), where F s (F s = F 1 ), H s (H s = H 1 ), J s (J s = J 1 ), and K s (K s = K 1 ) correspond to the only one block matrix determined by the single VSC in (26) and (28), with the corresponding network matrix derived from Y c,ii in (31), Here, the VSC-grid susceptance B g = 1 X g and θ 0 denotes the steady-state value of the phase-locked angle of the single VSC. VSC X g Figure 6. Schematic show for a single VSC infeed system.
For such a small system, the five eigenvalues of A single in (36) can be easily calculated. We will see that this simple-converter system plays a crucial role in our multi-converter system stability analysis.

Decoupling of Multi-VSC Infeed System with Symmetry
Now to analyze the multi-VSC infeed system, we consider some symmetry assumptions. As the first assumption, all VSCs should have the same control loops, control parameters, and DC capacity. Hence F i , H i , J i , and K i are all identical and can be simplified as F s , H s , J s , and K s , respectively. Therefore, we get F = diag(F s , . . . , F s ), H = diag(H s , . . . , H s ), J = diag(J s , . . . , J s ), K = diag(K s , . . . , K s ), (38) corresponding to the original (35).
In addition, Y c in (34) should also be decoupled and then the characteristic matrix A is solvable. Observing that both Y c,ii in (31) and Y c,ij (32) highly depend on the steady-state phase-locked angle difference θ ij0 (θ ij0 = θ i0 − θ j0 ), we can make the second assumption that the VSC-grid susceptance B ig should be identical. It is also notable that the inter-VSC susceptance B ij can be different. Under this assumption, we define for all i's (i = 1, ..., n). Hence we have θ i0 = θ 0 for all i's (i = 1, . . . , n). The purpose of the second assumption is to let all VSCs have the same operating points under the steady-state, including their state variables and the phase-locked angle. The current through their inter-VSC conductance is zero, just the same as the corresponding single-VSC infeed system.
Meanwhile, based on or from (12) and (13), and θ ij0 = 0, the expression of Y c can be simplified, with its selfadmittance and mutual-admittance matrices expressed as and from (31) and (32), respectively. According to (30), we have (44) under the above assumptions. Based on Y c,ij = Y c,ji , Y c is a partitioned symmetric matrix. Since symmetric matrices can be diagonalized in matrix theory [44], we perform the block-diagonalization of the block symmetric matrix as follows.
Observing that Y c,ji (i, j = 1, . . . , n) has similar structure with each other, we simply exchange the rows and columns of Y c according to the following rules: move the (3k − 2)-th row to the k-th one, move the (3k − 1)-th row to (k + n)-th one, move the 3k-th row to (k + 2n)-th one for all k's (k = 1, . . . , n), and applying the same rules to all columns, we obtain (45). Clearly, the above switching manipulation on the rows and columns of Y c is a similarity transformation, which does not change the eigenvalues of Y c .
where the symbol ∼ is used to denote the similarity diagonal transformation, I n is the identity matrix with the dimension of n × n, 0 n is the zero matrix with the dimension of n × n, and the matrix Y e with the dimension of n × n can be expressed as Then supposing B ei 's (i = 1, · · · , n) are the eigenvalues of Y e , we have Y e is diagonalized first and then let the rows and columns of Y c transformed in the reverse order of the previous sequence. Thus we realize the block diagonalization of Y c .
Substituting (48) into (34), we get the decoupled form of the characteristic matrix A and Till now, we achieve the decoupling of the whole n-VSC system. The eigenvalues of the sub-matrix A i (i = 1, . . . , n) together constitute the eigenvalues of A matrix. Clearly, A i has the exact same form as that for the single VSC system in (36), expect their specific matrices Y d,i in (49) for the decoupled subsystems and Y c,single in (37) for the single-converter system are slightly different. Below we would further analyze their relation between Y d,i and Y c,single .

Analysis for Subsystems
Now the focus is on B ei as it is the main variable of the matrix Y d,i in (49). Based on B ij = B ji in (12), clearly Y e is a real symmetric matrix, whose eigenvalues are real according to the matrix theory [44]. In addition, inputting (41) into (46) we get where Again on the basis of the matrix theory, we know that Y e is a Laplacian matrix, whose minimal eigenvalue is zero [44]. Therefore, a direct product is that the minimal eigenvalue of Y e is B g , namely, if we arrange B ei (i = 1, . . . , n) in an increasing order, we have and for all i ≥ 2. Therefore, inputting (54) into Y d,i in (49), we immediately find that This indicates that the first subsystem B e1 is just a single-VSC system under B g . Next let us study the cases for i ≥ 2. By a careful comparing of Y d,i in (49) and Y c,single in (37), we can see that these two matrices have a similar structure. If we define an equivalent susceptance B g,e f f for the subsystem i, from (37), we have its corresponding equivalent matrix, Note that here θ 0 corresponding to B g,e f f can be different with θ 0 under B g . The remaining point is to uncover the relation between Y d,i for the subsystem i in (49) and Y c,single for an equivalent single-VSC system in (58). After some cumbersome calculations, we can prove that and under certain conditions which are usually fulfilled. The details will be presented and the applicable conditions for the approximate equivalence will be discussed in Appendix B. Therefore, we find that each subsystem can be regarded as a single-VSC system under an equivalent grid susceptance B g,e f f , which is identical to the corresponding eigenvalue B ei of the matrix Y e . Meanwhile, based on the mature experience in power systems, a smaller line susceptance (or a larger line inductance) could ruin its stability in the single-VSC system. Hence, we conclude that the first subsystem under B e1 (B e1 = B g ) should determine the stability of the whole system, based on the fact that B e1 is minimal among all eigenvalues of B ei . In this respect, we obtain that the converter-grid susceptance B g (or reactance X g ) is the most important parameter, whereas all inter-converter susceptance B ij 's (or reactance X ij 's) could only influence the stability of subsystem i (i ≥ 2), but not that of the whole system.
The whole derivation process is clearly illustrated in the flow chart in Figure 7.
Nonlinear model of n-VSC system Linearized model of n-VSC system Unified characteristic matrix A of n-VSC system Characteristic matrix A single of single-VSC infeed system, whose network parameter is B g .
Characteristic matrix A of n-VSC indeed system. Under symmetry conditions n characteristic matrics A i (i =1,…,n) Decouple network matrix first and then decouple characteristic matrix .

One characteristic matrix
A 1 , whose network parameter is B e1 =B g .
n-1 characteristic matrics A i (i=2,…,n),whose network parameter is B ei , respectively. n-1 characteristic matrics A i (i=2,…,n) of single-VSC infeed systems, whose network parameter is B ei , respectively. approximately equal equal Conclusion 1: A n-VSC system under symmetry conditions can be decoupled into n single-VSC systems B g is always minimal among B ei (i=1,…,n) Conclusion 2: The stability of a n-VSC system is determined by the single-VSC system corresponding to A 1 .

Single-VSC Infeed System
As a comparison, we first study the stability of the single-VSC system in Figure 6 with the parameters summarized in Appendix A. The root locus results for different network reactance X g 's are presented in Figure 8, with an increase of X g from 0.1 (circle) to 1.1 (star) representing a weaker system. These calculations show a clear critical reactance X g,c : X g,c ≈ 0.98, for a pair of λ 1 and λ 2 coming across the imaginary axis. For X g > X g,c , the system would become unstable.
Furthermore, the participating factor analysis is employed to identify the relations between eigenvalues and control loop parameters [4][5][6]. We find that λ 1 and λ 2 as a pair are mainly related to the DVC, λ 3 and λ 4 as a pair are connected with the PLL, and λ 5 is associated with the TVC. With X g gradually increasing for a weaker grid, λ 1 and λ 2 become the dominant modes and make the system unstable, but meanwhile λ 5 oppositely moves from right to left, seemingly enhancing the system stability. Imaginary Part Figure 8. Root loci of the single-VSC infeed system with the converter-grid reactance X g varying from 0.1 (circle) to 1.1 (star). In particular, a critical reactance X g,c exists; X g,c ≈ 0.98. For X g > X g,c , the system becomes unstable.

Multi-VSC Infeed System
Without losing generality, a three-VSC infeed system with symmetry is studied under the condition that all VSC-grid reactances are identical X gi = X ig = X g , for i = 1, . . . , n. The VSC parameters remain unchanged in Appendix A. Extensive numerical simulations have been conducted and found to support the above theoretical analysis.
In the first test, the VSC-grid reactance is fixed but the inter-VSC reactance changes. We fix X g1 = X g2 = X g3 = X g = 0.5, X 13 = 0.2, and X 23 = 0.8, but change X 12 from 0.1 to 1.1. Again the root loci are calculated and illustrated in Figure 9a. In addition, we calculate the three eigenvalues of Y e : B e1 , B e2 , and B e3 , and treat them as the effective network conductances in three independent single-VSC systems. For example, under X 12 = 0.1 (circle), B e1 = 2 (or X g1,e f f = 1/B e1 = 0.5), B e2 ≈ 10.7 (or X g2,e f f = 1/B e2 ≈ 0.09), and B e3 ≈ 25.9 (or X g3,e f f = 1/B e3 ≈ 0.04). Additively under X 12 = 1.1 (circle), B e1 = 2 (or X g1,e f f = 1/B e1 = 0.5), B e2 ≈ 5.2 (or X g2,e f f = 1/B e2 ≈ 0.2), and B e3 ≈ 13.1 (or X g3,e f f = 1/B e3 ≈ 0.08). Then, the root loci of all these three decoupled single-VSC systems are plotted in Figure 9b-d, respectively Note that Figure 9b is actually the decoupled system corresponding to X g . Comparing these three panels (Figure 9b-d) with Figure 9a, one can clearly see that the whole three-VSC system indeed has been divided into three single-VSC systems, with nearly unchanged eigenvalues. This perfectly demonstrates the validity of our theory. Finally, recall that though λ 5 's in Figure 9c,d are near to the imaginary axis compared to Figure 9b, it does not mean that these two subsystems in Figure 9c,d are more likely to cause instability than the dominant subsystem in Figure 9b. Because λ 5 always moves from right to left for a larger X g and won't cause instability, as we can see in Figure 8.
For the second test, we fix the inter-converter reactance, X 12 = X 23 = X 13 = 0.5, but vary the converter-grid reactance. The root loci for an increasing X g from 0.1 (circle) to 1.1 (star) are illustrated in Figure 10a. Similarly, Figure 10b-d show the corresponding results for the three decoupled single-VSC systems. By comparing these sub-figures, again the whole system has been well divided into three single-VSC systems and clearly the dominant subsystem is in Figure 10b. Some small visible deviations in Figure 10 might come from the approximation induced by Equation (60).   Figure 9, but for another parameter test, instead. Now we fix X 12 = X 23 = X 13 = 0.5 and change X g from 0.1 (circle) to 1.1 (star).
To be clearer, we present the corresponding eigenvalues data of Figures 8-10 in Appendix C.

Time-Domain Simulation Results
To verify the above small-signal stability analysis results, sufficient time-domain simulations on the original nonlinear system were conducted in the MATLAB/Simulink. The controller parameters used are listed in the Appendix A. First, the single-VSC system in Figure 6 is studied with the identical device parameters. Three tests for different VSC-grid reactances, including X g = 0.9, X g = 0.98, and X g = 1.0 were conducted. A small step disturbance of amplitude, 0.01 p.u., was added on i d0 at the beginning, which indicates that there is a small node injection current disturbance. The corresponding three terminal voltages of VSC are shown in Figure 11a, where the three curves are for stable, critical stable, and unstable behaviors, respectively. The occurrence of a critical VSC-grid reactance at X g,c ≈ 0.98 in the plots is in accordance with the eigenvalue analysis in Figure 8 before.  Figure 11. Time-domain simulation results for the single-VSC infeed system (a) and the three-VSC infeed system (b). Both cases show that the critical reactance X g is 0.98. In (b), for other parameters, X 12 = 0.5, X 13 = 0.2, and X 23 = 0.8 are arbitrarily chosen.
For comparison, a three-VSC system has been studied. As the condition for identical VSC-grid reactances is necessary in our theory, we consider X 1g = X 2g = X 3g = 0.9, X 1g = X 2g = X 3g = 0.98, and X 1g = X 2g = X 3g = 1.0, corresponding to the above three cases of a single-VSC system. The other parameters for the inter-VSC reactances X 12 = 0.5, X 13 = 0.2, and X 23 = 0.8 has been arbitrarily chosen. The same disturbance is given to one VSC, with its corresponding terminal voltage waveform shown in Figure 11b. Comparing these two panels in Figure 11, one can clearly see that the critical VSC-grid reactance in the multi-converter system is nearly unchanged; X g,c ≈ 0.98. This again demonstrates that the whole large-scale system dynamics is indeed solely described by only one dominant subsystem with the key parameter X g .

Conclusions
In conclusion, we have developed an appropriate theory for the small-signal stability in multi-converter systems, from nonlinear modeling to linear modeling and its approximation, based on three important techniques: multi-timescale decomposition, linearization, and symmetry decoupling. We have found that, under the approximation assumptions, the whole huge system can be well-decomposed into several subsystems and each subsystem can be treated as a single-converter system with an effective network susceptance. The stability of the large-scale multi-converter system is determined by only one dominant subsystem, whose network parameter is equal to the original converter-grid susceptance.
Essentially, the key to decouple the multi-machine power system into several singlemachine systems is to decouple the power network, which is often represented by the admittance matrix in the mathematical model. The matrix itself has a natural symmetry. Only when the matrix of each device connected in the network does not destroy this symmetry, the characteristic matrix of the whole system keeps symmetric and then we can decouple the power system. Therefore, we put forward strict conditions for the device and network in this paper to meet the requirement of precise decoupling. However, when the parameters of the device or network appear little deviation or the power flow in the network is not obviously asymmetric, the decoupling condition may be relaxed, so that it can be better applied and popularized in real scenes. We believe that this theory provides a theoretical tool for analyzing integrated renewable power systems and helps to find the key transmission lines.
On the other hand, as is known, the multi-converter system is also a typical networked dynamical system consisting of dynamical nodes and a stationary network. Possible methods and techniques from relevant disciplines, such as complex systems, nonlinear dynamics, and statistical mechanics as well, are highly appreciated. Different from the well studied synchronization process in traditional power grids, the physical mechanism of the power-electronic-based power grid remains obscure. We hope that our paper could fill the gap between power electrical engineering and contemporary physics by a clear description of modeling details of power electronic devices. It could arouse extensive interests among physicists with network dynamics background for further investigations. g,e f f − P 2 in , (A8) whose squared deviation is expressed as We find that if P in = 0 or B g,e f f = B g , it will vanish. From this, one understands the following two conditions: (1) when P in,i is relatively small, i.e., the power delivered by transmission line is far from the steady-state-stability-limited power, or (2) B g,e f f should not be far away B g , or B ij should be relatively small compared to B g , which means that the inter-converter connection should be weak compared to the converter-grid connection. Under these two conditions, we can have the approximative equivalence of Y d,i (1,3) and Y c,single (1,3) in (60). Clearly the above conditions also fit with our intuition.
Finally, there is one slight point that remains to be addressed. As we have found that among the four matrices F i , H i , J i , and K i in (26) and (28), only H i depends on the steadystate values of i d0,i and i q0,i , whereas all other matrices depend on system parameters only, we should also prove that H i (or i d0,i and i q0,i ) is roughly unchanged under the condition of B g,e f f . By using the phasor relation in Figure A1, we have found that it is right under the same two conditions above. The detailed derivation is not shown here.  1   Table A3. The detailed eigenvalues data of Figure 10.