A Gray-Box Stability Analysis Method of Grid-Connected Inverter Considering Synchronization Dynamics

: The Grid-Connected Inverter (GCI) can be considered a gray box when circuit and controller parameters are missing due to intellectual property rights or parameter variations caused by aging, which poses an impediment to assessing the stability of the system. This paper presents a gray-box stability analysis method based on impedance identiﬁcation of GCI considering the synchronization dynamics. The impedance frequency responses of GCI are measured by the frequency scanning method on the dq-frame. Meanwhile, the inﬂuence of synchronization dynamics and background harmonics is theoretically investigated. A vector ﬁtting (VF) algorithm, co-designed with impedance identiﬁcation, is then applied to generate polynomial transfer functions. Based on the obtained transfer functions, the stability of the GCI can be judged by the distance relationship between the prohibited area boundary and the center of the gershgorin-circle through the distance formula. Finally, the experiments of both RT-LAB and experimental prototypes are conducted to verify the feasibility of the proposed method.


Introduction
Grid-connected inverters (GCIs) are commonly used with renewable energy resources, power transmission systems to transmit power into the utility grid [1,2]. The ever-increasing penetration of GCI is radically changing the dynamic operations of traditional power systems. And thus, instability phenomena triggered by the GCI under different grid conditions have been increasingly reported [3,4].
To tackle the stability problems introduced by the GCI, the impedance-based stability analysis method is applied in [5,6]. In [7], an impedance analysis based on the dq coordinate frame is presented. Further studies are reported in [8,9], which reveals that the interactions between inverter and grid are characterized by the generalized Nyquist stability criterion (GNSC), and the multiple-input and multiple-output (MIMO) impedance matrix based on the dq-frame can be decomposed into four single-input-single-output (SISO) transfer functions for the stability judgment. Yet, the detailed parameters of the system are needed when implementing the transfer function-based method, which is usually a challenging task in practice, caused by the confidentiality of control parameters [10], and the variation for hardware parameters due to temperature fluctuation, and aging [11,12]. Therefore, the impedance identification to obtain the system dynamic response for the stability analysis of the GCI is increasing demanded [13].
In [14,15], a vector fitting (VF) algorithm combining impedance identification based on a frequency scanning method is proposed, and the stability of the system can be measured by the system's eigenvalues. Furthermore, the instability sources in the system, component, and parameter levels of the multiple-inverter-fed power systems can be identified through VF and impedance frequency responses [16,17]. The VF algorithm was utilized to fit the impedance frequency response, from which the impedance transfer function was obtained. The control parameters are obtained by comparing with the actual control structure. Furthermore, the VF algorithm can be utilized for impedance remodeling in gray box systems [18] or reducing the model order of power cables [19]. In [18], a gray-box parameter identification method has been presented to identify the internal parameters of GCI. By adjusting the internal parameters of an inverter, the terminal output impedance can be reshaped to mitigate the oscillation issue that results from time-varying grid impedance. In [19], Prony analysis (PA) and VF are used to fit the frequency response and obtain the state space model of power cables, reducing the computational burden of power cable modeling and improving modeling and analysis efficiency. However, the identification of impedance by the aforementioned method does not consider the influence of perturbation, such as background harmonics, grid frequency deviation [20], and Phase Locked Loop (PLL) dynamics [21]. To improve the accuracy of the impedance measurement, more considerations, including grid harmonics, to improve the identification performances via twice-measurement can be found in [22], and the measurement method by the frequency selection principle has been investigated in [20,23]. During the impedance measurement, the voltage and current in the abc-frame need to be measured before being transformed into the dq-frame. The perturbation signals used for impedance measurement are usually designed in the dq-frame and then converted to the abc-frame. Therefore, PLL is typically used to track the system frequency online and estimate the synchronization phase angle for coordinate transformation. To reduce non-negligible phase errors caused by frequency deviation of the grid voltage, the PLL control bandwidth needs to be relatively high. However, controlling a PLL with an excessively high bandwidth may introduce undesired disturbances to impedance calculations [24,25]. And the errors of impedance identification increase with the bandwidth of PLL [21]. If a low-bandwidth PLL is employed, the tracking performance is poor under the circumstance of the time-varying frequency of grid voltage, which leads to non-negligible phase errors [26]. Hence, further research thinking about the impact of the synchronization dynamics on impedance identification should be carried out.
On the other side, the GNSC has been widely applied to the stability analysis of GCIs at the cost of complicated computation [27]. It is revealed that the GNSC is rigorous in mathematics, the eigenvalues of the impedance matrix of the MIMO systems need to be calculated at each frequency, which is a tedious process. Besides, the eigenvalue loci of MIMO systems exhibit implicit characteristics, which results in an unintuitive analysis process [28]. To facilitate the analysis, the non-diagonal elements are often neglected, which brings hidden problems to the stability analysis. In [27,29], there are two gershgorin circle bands around the diagonal elements loci when considering the non-diagonal elements, but it is still a cumbersome process on the premise of a reduced computation of the criterion. In [30,31], a G-Sum criterion has been presented, which is concisely compared to the GNSC, but it is highly conservative. In order to reduce the conservatism, the forbidden-region criterion has been proposed in [32], but the trouble still exists. In [33], the stability criterion expressions are derived by analyzing the distance between the center of gershgorin circle and the prohibited area in different cases, in which the computational procedure is relatively simple. However, four different cases need to be analyzed respectively, so that a cleaner approach is expected to evolve. For the moment, the balance between reduced computation and conservatism for the stability criterion is still an open issue.
To address the above challenges, this article presents an improved impedance calculation method, which makes it possible to accurately extract the phase difference in PLL during the frequency scanning, and the errors caused by synchronization dynamics are corrected. Besides, the frequencies of perturbations are selected to adapt background harmonics. Moreover, the distance relationship between the prohibited area and the gershgorin-circle in the complex field is divided into two different cases through the distance formula. Then, the stability of GCI can be judged. The proposed stability criterion has a much simpler computational procedure compared with the existing stability criterion.
Finally, based on the measured impedance and proposed stability criterion, the stability of the gray box GCI in different grid conditions can be analyzed.
The remainder of the paper proceeds as follows. Section 2 introduces the impedance measurement of the three-phase grid-connected inverter. The improved impedance calculation method and the proposed sufficient stability criterion and their corresponding simplified analysis are presented in Section 3. Section 4 investigates the proposed method by simulation and experiments. Finally, the conclusions are made in Section 5.

Impedance Measurement Three-Phase Grid-Connected Inverter
In order to simplify the analysis, an L-filtered GCI is used to discuss the proposed method. Figure 1 is the three-phase GCI and its control structure, which is based on PI controller in dq-frame. The U dc is DC bus voltage, L inv are filter inductors, U g is grid voltage, and Z g is equivalent grid impedance. I * dq is the reference value of the grid-connected current in the dq-frame, U abc , I abc is the grid-connected voltage and grid-connected current in the point of common coupling (PCC), respectively, I dq is the grid-connected current in the dq-frame. θ is the grid synchronization angle from PLL, PI is the proportional-integral current controller. The impedance analysis method shows that the GCI system is stable when the ratio of the grid impedance and the output impedance of the GCI satisfies the GNSC. In a weak grid, grid impedance exhibits resistance inductance, its expression in the dq-frame is shown in Equation (1). Here, ω represents the grid synchronous angular velocity.

PLL
The ratio between the output impedance of the GCI and the grid impedance in the dq-frame, which can be expressed in terms of the return rate matrix. The rate matrix is denoted as Equation (2).
In Equation (2), L dd (s), L dq (s), L qd (s), L qq (s) are the elements in the return rate matrix L(s). Y inv (s) is the output conduction of the GCI. This can be measured by injecting voltage or current disturbances at the PCC (Point of Common Coupling). The perturbation response obtained in the abc-frame is then transformed into the dq-frame. Subsequently, the amplitude and phase of the disturbance response can be calculated using Fast Fourier Transform (FFT). The mathematical model of the inverter voltage disturbance and current response can be expressed as In Equation (3), ∆ is the small disturbance of voltage and current at the equilibrium operating points. Z dd , Z dq , Z qd , Z qq are four unknown quantities to be measured, which represent the dq-axis impedance of the GCI and its cross-coupling impedance. In order to measure the four independent unknown quantities in Equation (3), two sets of linearly independent perturbations are injected into the system separately in the same frequency. The impedance model of the system in the dq-frame can be obtained as the following: In Equation (4), the subscript '1', '2' denotes it as two sets of linearly independent perturbations and their responses, U d , U q , I d , I q are the voltage and current in the dq-frame.

Impedance Measurement Error and Its Compensation Method
It is necessary to obtain the synchronous angle of the grid for disturbance injection and impedance calculation during impedance measurement. It is important to note that, the PLL for disturbance injection and impedance calculation is different from the PLL used in the control loop of the three-phase grid-connected inverter. Disturbance signals injected into the system and the system frequency variations affect the phase angle of the PLL, which in turn affects the calculation of impedance measurements. The GCI system exists in two different coordinate systems under the dq coordinate system in practice when the PLL dynamic is considered [30]. Considering the effect of the perturbation measurement signal injected into the GCI system on the phase angle of the PLL output during impedance measurement, the different dq-frame is represented as Figure 2. In Figure 2, the superscript 'p', 's' denotes the perturbation-voltage dq-frame and the PCC-voltage dq-frame, respectively. ∆θ 1 is the phase angle difference between the different dq-frames. The angular difference between different dq-frames has an impact on both the perturbation injection and impedance calculation processes, and the impact of PLL dynamics is smaller in the perturbation injection process; however, it has a greater impact on accuracy in the impedance calculation process [21]. The trigonometric function for conversion between different frames can be approximated as sin(∆θ 1 ) ≈ ∆θ 1 and cos(∆θ 1 ) ≈ 1. The Voltage components between different frames is expressed as In Equation (5), U s d0,q0 and U p d0,q0 are steady-state voltage components in each respective frame. U s d,q and U p d,q are minor disturbance voltages, θ 1 is phase angle in the perturbation-voltage dq-frame. Thus, the minor disturbance voltages between different frames is expressed as [30] The small signal transfer function of traditional PLL in the GCI can be represented as [26] ∆θ Therefore, the relationship of perturbation measurement signals in Equation (4) between different frames is Similarly, the relationship between the current in the PCC-voltage dq-frame and that in the perturbation-voltage dq-frame can be simplified as According to (8) and (9), the relationship between the measured impedance Z p dq and the actual impedance in PCC Z s dq is derived as The mismatch between the measured impedance and the actual impedance caused by the dynamic of the PLL can be mitigated when the term U s d0 G pll and I s d0 G pll are sufficiently small. From Equation (7), the terms U s d0 G pll and I s d0 G pll are equivalent to a low-pass filter. When the perturbation signal frequency is higher than the PLL bandwidth, it will be filtered, making it difficult to influence the measurement results. Yet, when the perturbation frequency is lower than the bandwidth of the phase-locked loop, the effect of the perturbation signal cannot be ignored. It is an effective way to attenuate the effect of disturbing signals on measurement results by reducing the PLL bandwidth when the frequency deviation in the grid has been ignored. If the bandwidth of the PLL is too low, there is a non-negligible angle difference ∆θ 1 when the system is subjected to harmonics and grid frequency deviation [26]. There is, thus, the impedance measurement error caused by synchronization dynamics, which cannot be completely ignored, even though the bandwidth of the PLL is designed suitably.
From the small signal model for PLL [31], the angle difference ∆θ 1 between the PLL output phase angle and the input voltage phase angle can be simplified as From (11), the angle difference ∆θ 1 can be extracted from PLL, and the voltage and current in the PCC-voltage dq-frame (for calculated actual impedance) can be acquired after the correction based on (8) and (9). The actual impedance of the GCI can be obtained consequently.
To calibrate the voltage and current in the perturbation-voltage dq-frame, it is crucial to accurately extract the phase differences ∆θ 1 when injecting the disturbed signal using the frequency scanning method. In this method, a single frequency excitation signal is injected at a time, and the resulting phase difference contains the corresponding frequency component. A phase extraction method is presented in Figure 3, taking into account both the filtering performance and the characteristics of the frequency scanning method. In this method, a second-order band-pass filter is utilized to extract the phase differences during impedance measurement. Its transfer function can be derived as follows: In Equation (12), ω n is the frequency of the disturbed signal in frequency scanning, and the ξ is chosen as 1.5.
In addition, the background harmonics of the inverter caused by the dead time and the background harmonics of the grid derived from nonlinear loads has the impact of impedance measurement [20,34]. Therefore, the frequencies of the injected perturbation signals must be carefully selected to minimize errors. Ideally, the frequencies of the perturbation signals should not overlap with the background harmonic frequencies (which are typically integral multiples of the fundamental frequency), and the signal should be sufficiently strong. By ensuring these conditions are met, the impedance can be accurately obtained. Figure 3. The improvement of PLL structure for impedance measurement.

Improved Stability Criteria Based on Gershgorin Circle Theorem
This paper considers the stability of GCI connected to the grid in weak grid conditions, but does not consider the situation when the grid is unstable (Z g (s) exists the unstable pole) or when the new energy-generating units are unstable (Y inv (s) exists the unstable pole). The GCI system is stable when the eigenvalue locus of the return rate matrix L(s) is unencircled (−1, j0). In GNSC, calculating system eigenvalues is a tedious process. The circle theorem uses the diagonal elements of the matrix to replace the eigenvalues of the system, ignoring the errors caused by the non-diagonal elements of the matrix, and estimates the range of the eigenvalues of the system by the diagonal and non-diagonal elements of the matrix. For the matrix of system A = a ij , the magnitude of the offdiagonal elements in row i is r i (A) = ∑ i =j a ij . All eigenvalues distribute in a gershgorin circle, with centers at the diagonal elements a ii and radii r i (A). Thus, relationship between the eigenvalues of the return rate matrix L(s) and its four elements can be represented as follows: The forbidden region of the forbidden region based criterion (FRBC) is on display in Figure 4.
Eigenvalue trajectories of L(s) satisfy the GNSC when any frequency point of two gershgorin circles lies to the right of the parallel line Re = −1 (all gershgorin circles are not intersecting with the grey areas in Figure 4). Expressions for the stability of FRBC are as follows: From the foregoing, it is clear that the discussion on the diagonal elements L dd and L qq is basically the same. Thus, the proposed method is discussed using the element L dd as an example.
In order to establish the exclusion zone, the stability margin parameters h (0 < h < 1) and p (0 < p < π/2) are introduced. The forbidden region of the proposed method is displayed in Figure 5.  In Figure 5, the complex field is divided into two cases (case1 and case2) by two straight lines (l 1 : Im = tan(π/2 − p)(Re + h) and l 2 : Im = − tan(π/2 − p)(Re + h)), and the grey area is the forbidden region in the proposed method. When L dd is located in case1, it is represented as follows: The distance between L dd and the forbidden region is indicated in Figure 5 by a green line, and its equation is as follows: When L dd is located in case2, that is The distance between L dd and the forbidden region is indicated in Figure 5 by a red line, and the equation can be simplified as follows: Then, the stability judgment can be obtained by comparing D(s) and L dq .
When the stability criterion condition W(s) is less than 0, the system is unstable. Apparently, the conservative condition is related to the stability margin parameters (h and p).

Results and Discussion
Based on the topology and the control block diagram shown in Figure 1. The impedance identification simulation model of the grid-connected inverter system is built, and the specific parameters of the model are shown in Table 1. Due to the lack of an adequate hardware platform, the proposed impedance measurement method is verified by RT-LAB. Using the RT-LAB real-time simulation platform from OPAL-RT Technologies, the GCI simulation model is built according to the system block diagram in Figure 1 and the parameters in Table 1, and the real-time simulation experimental setup is shown in Figure 6. The grid-connected inverter control and disturbance injection are implemented in the OP5700 real-time digital simulator, and the data obtained in the OP5700 are swept for impedance calculation in the Matlab software of the host computer. The discrete impedance data obtained are fitted with VF to obtain the analytical formula of the inverter output impedance for the stability analysis of the inverter in the subsequent section. The VF algorithm is a rational function approximation algorithm that can fit the discrete data obtained from impedance measurements to obtain the impedance transfer function of the GCI [13].
Based on the GCI model be built in the OP5700 real-time digital simulator, the voltage disturbance sources with amplitudes ranging from 5∼10% of PCC voltage are sequentially connected in series at the PCC. In addition, the bandwidth of the PLL in the impedance measurement needs to be high to track the frequency deviation, which is k ppll = 0.5 and k ipll = 50.
In Figure 7, it shows the impedance measurement results based on different measurement methods considering both the influence of the varying system frequencies, background harmonics and injected perturbations. The blue line and green dots indicate that the high bandwidth of the traditional method significantly influences the results of Z dq , Z qd and Z qq . Nevertheless, the measured results match the analytical model to a large extent when the method proposed in this article is adopted, as shown by the red points. However, the dq coupling of the model utilized in this paper is relatively small, particularly within the low-frequency range. As a result, the perturbation voltage used for measuring Z dq and Z qd is unable to elicit a sufficient response across the cross axis, leading to significant discrepancies in Z dq and Z qd at low frequencies. Furthermore, the Z dq and Z qd of the model utilized in this paper is relatively small. Thus, it does not affect the stability judgment through the method proposed in this paper. The VF algorithm is able to fit the frequency response of the discrete impedance series in the s-domain and derive a polynomial transfer function [15,16], and the polynomial transfer function can be shown as follows: In Equation (20), m represents the order of the polynomial transfer function, E is nonzero when fitting the impedance model of the L-filtered GCI. During the VF fitting process, the impedance matrix based on the dq-frame, which is multiple-input and multipleoutput (MIMO), can be broken down into four single-input and single-output (SISO) transfer functions. Before applying the VF algorithm to fit the discrete frequency response data, it is necessary to determine the order and initial poles of the function. Taking the fitting process of Z dd as an example, the VF procedure can be explained in detail.
According to [35], the order of the model is an important parameter for VF, which can be chosen based on some information extracted from the system frequency response samples. In short, the number of initial poles in the VF algorithm can be quickly determined based on the resonance peaks found in the sample data, and one resonance peak corresponds to a conjugate pole pair. In addition, the frequency response of real poles is smooth. The major resonance peaks were detected using the Equation (21) when the frequency response samples Y(si) were obtained. y 11 s p−1 < y 11 s p > y 11 s p+1 (21) According to [35] and Equation (21), the order of the transfer function about Z dd can be selected as 3. However, increasing the order can reduce the fitting error, and in practice, selecting an order of 5 for the GCI results in a small enough fitting error [17].
After determining the initial poles, the frequency response data in Figure 7 were fitted to obtain a polynomial transfer function, as shown in Equation (20), and the coefficients of the fitted 5-order polynomial transfer functions about Z dd are shown in Table 2. Similarly, we can use the same method to fit the data of Z dq , Z qd and Z qq , as shown in Figure 8. Then, the return rate matrix can be obtained based on the fitting results and grid impedance. The curves of h1 and h2 (h1, h2 is the eigenvalue curves based on the stability criterion proposed in this article) under different SCRs are obtained in Figures 9 and 10.    In Figures 9 and 10, the measurement error in Z dq and Z qd will not affect the stability judgment. According to the stability criterion proposed in this paper, when h1 or h2 is less than zero at any frequency, the system can be considered unstable. As shown, as the SCR drops to 3.5, the curve of h2 has two intersection points with the zero plane. Consequently, the GCI system is unstable. According to Figure 10, the frequency of the intersection of λ 2 with the unit circle less than 215 Hz.
To enhance the depiction of the instability phenomenon witnessed during the experiment, the Nyquist curves of the eigenvalue of the return rate matrix L(s) is introduced, as shown in Figure 11.
According to Figure 11, the characteristic value λ 2 of L(s) is surround the point (−1, j0), the GCI system is unstable. And the frequency of the intersection of λ 2 with the unit circle is 173 s.
The experimental results of the GCI connected into a weak grid with different SCRs are based on the experimental platform in Figure 12.
Before the stability experiment, the dynamic experiment in Figure 13 was used to demonstrate the dynamic performance of the inverter system.
The experimental results under different SCRs are given in Figures 14 and 15. As shown, when the SCR decreases from 5 (shown in Figure 14) to 3.8 (shown in Figure 15), the system changes from a stable state to an unstable one. According to Figure 14a, the output current of the GCI is stable when SCR = 5. Yet, in Figure 15a, the output current of the GCI suffers from harmonic instability [2]. By importing data from an oscilloscope into matlab, THD analysis results of the gridconnected current are obtained, as shown in Figures 14b and 15b. According to Figure 15b, the frequency of harmonic instability is 195 Hz, and the coupling frequency is 95 Hz. In [2], the frequency of harmonic instability is almost consistent with the result of adding or subtracting the fundamental frequency at the frequency of the intersection of λ 2 with the unit circle. Therefore, the frequency of harmonic instability is 173 ± 50 Hz in the theoretical analysis. Clearly, discrepancies can be observed between the theoretical analysis results and the experimental findings. However, based on the aforementioned analysis, the proposed method remains feasible even when considering the errors between the actual circuit and the simulation model.

Conclusions
This paper has analyzed the impact of synchronization dynamics and background harmonics on impedance measurement. An improved impedance calculation method is proposed, incorporating a frequency selection principle for injected voltages during frequency scanning. This enables accurate identification of perturbations injected into the measurement system. The proposed method ensures a sufficient bandwidth for the PLL to track grid frequency deviations and mitigates the influence of synchronization dynamics caused by perturbations. Consequently, the impedance can be accurately determined using this method.
To assess stability based on the measured impedance, a stability criterion based on the Gershgorin circle theorem is introduced. This criterion reduces the computational burden compared to existing methods. Simulation and experimental results of the GCI system validate the effectiveness of the proposed method.

Conflicts of Interest:
The authors declare no conflict of interest.