The topology of the GFM grid-connected system, including the virtual synchronous control and AC voltage outer loop-current inner loop control, as shown in Figure 1.

As shown in Figure 1, V _s is the terminal voltage vector of the converter, V and I are the voltage and current vectors of the point of common coupling (PCC), respectively; L _f and C _f are the filtering inductor and capacitor, respectively; L _g is the equivalent inductance of the power grid; I _g is the current vector of the grid, and E is voltage vector of the infinite bus. The above vectors are all established in the abc stationary coordinate frame.

The virtual synchronously controlled GFM converter is synchronized with the grid by simulating the second-order swing equation of the synchronous machine rotor, whose detailed expression is shown as follows: d θ d t = ω · ω b J d ω d t = P ref − P e − D ω (1) where θ is the output angle of virtual synchronous control, indicating the angle between the d-axis of the converter control coordinate system and the x-axis of the grid rotational coordinate system. ω _b is the fundamental frequency, ω is the output frequency of the GFM converter, J and D are the inertia and damping coefficients, respectively, P _ref is the reference value for active power, and P _e is the output active power of converter.

The converter with AC dual-loop control provides a reference signal to the current inner loop through the AC voltage outer loop, which in turn generates the pulse width modulation (PWM) control commands. The dynamic equations of the current inner loop and the AC voltage outer loop are shown as Equations 2, 3, respectively. V s d = V d + x 1 + k CCp I d ref − I d − ω 0 L f I q V s q = V q + x 2 + k CCp I q ref − I q + ω 0 L f I d (2) I d ref = x 3 + k VCp V d ref − V d + I g d − ω 0 C f V q I q ref = x 4 + k VCp V q ref − V q + I g q + ω 0 C f V d (3) where V _sd and V _sq are the d-axis and q-axis components of the converter terminal voltage, respectively; V _d and V _q are the d-axis and q-axis components of the voltage at the PCC point, respectively, with V _dref and V _qref as their references; I _dref and I _qref are the d-axis and q-axis reference currents, respectively; k _CCp and k _VCp are the proportional coefficient of the PI control in current inner loop and AC voltage loop, respectively; I _gd and I _gq are the d-axis and q-axis components of the grid current, respectively; ω ₀ is the nominal frequency, taken as 1. Furthermore, the integrator outputs of the two pairs of PI control for the current inner loop and voltage outer loop are selected as state variables, denoted as x _n (n = 1,2,3,4), whose expressions are shown as follows: d x 1 d t = k CCi I d ref − I d d x 2 d t = k CCi I q ref − I q (4) d x 3 d t = k VCi V d ref − V d d x 4 d t = k VCi V q ref − V q (5) where k _CCi and k _VCi are the integral coefficients of the PI control in the AC voltage outer loop and current inner loop.

In addition, V _dref in Equation 5 is controlled by reactive power, as shown in Equation 6: V d ref = V ref + K q Q ref − Q e (6) where K _q is the droop coefficient, Q _ref is the reactive power reference value, Q _e is the output reactive power of the converter.

According to Figure 1, the main circuit of the converter under abc coordinate frame can be obtained: L f d I d t = V s − V C f d V d t = I − I g (7)

Equation 7 can be expressed in the dq coordinate system as: L f ω b d I d d t = V s d − V d + ω L f I q L f ω b d I q d t = V s q − V q − ω L f I d C f ω b d V d d t = I d − I g d + ω C f V q C f ω b d V q d t = I q − I g q − ω C f V d (8)

Similarly, the dynamic function of the AC grid in Figure 1 can be described by the following equation: L g d I g d t = E − V (9)

Different from Equation 8, the dynamic equations of the AC grid should be discretized in the xy global coordinate frame, as shown in Equation 10: L g ω b d I g x d t = V x − E x + ω 0 L g I g y L g ω b d I g y d t = V y − E y − ω 0 L g I g x (10) where V _x and V _y are the x-axis and y-axis components of the voltage at the PCC, respectively; E _x and E _y are the x-axis and y-axis components of the grid voltage, respectively.

Combining Equations 1, 4, 5, 8, 10, the state space function for the converter grid-connected system can be obtained.

The GFM system model is a nonlinear system of state equations consisting of 12 and 5 input variables, as detailed in Supplementary Appendix A.

Equations 1–10 illustrate that the complete state-space expression contains 12 variables, whose exact form is given in Supplementary Appendix A. However, the grid-side dynamic equations are modelled under the xy coordinate frame, while the rest are modelled under the dq coordinate frame. To facilitate the analysis, the grid-side dynamic equations need to be transformed under the dq-axis, and such a process is realized by the transformation matrix T : T = cos ⁡ θ sin ⁡ θ − sin ⁡ θ cos ⁡ θ (11)

The relationship between the voltage vector at the PCC point V , the grid current vector I _g, and the infinite bus voltage vector E in the xy coordinate frame and dq coordinate system can be described as follows: V x y = T − 1 V d q E x y = T − 1 E d q I g x y = T − 1 I g d q T − 1 = cos ⁡ θ − sin ⁡ θ sin ⁡ θ cos ⁡ θ (12) where subscripts xy and dq represent the forms of the corresponding variables in the xy coordinate frame and dq coordinate frame, and T ^-1 represents the inverse of T . Therefore, linearizing Equations 1, 4, 5, 8, 10 with Equation 11, the small-signal linearization equation of the AC power grid in the dq coordinate system can be obtained as shown in Supplementary Appendix A.

By combining (A2) and (A6) in Supplementary Appendix A, the system state space equation can be obtained as follows: E d d t Δ x t = A Δ x t (13) where E is a 12th-order nonsingular matrix, A is the constant matrix of the closed-loop system at the equilibrium point, and its block matrix expression is shown in Supplementary Appendix A, (A7)-(A8). The small-signal state space model of the GFM converter grid-connected system can be described by Equation 12.

To clarify the impact of grid strength on the stability of GFM converter grid-connected systems, the SCR is commonly used as a metric. SCR at a certain point of the grid is usually defined as the ratio of the short-circuit capacity to the rated capacity of the system, as shown as follows: SCR = S ac S B = 1 S B Z (14) where S _ac represents the short-circuit capacity of the PCC point in the AC system, S _B is the rated capacity, Z = ω ₀ L _g is the equivalent impedance of the AC power grid.

The SCR at the critical stability condition is defined as CSCR, which can be employed to characterize the stable boundary of the system (Badran et al., 2025). By calculating the difference between CSCR and SCR, the small-signal stability margin of the GFM converter grid-connected system can be determined.

The impedance method is one of the commonly employed methods in the stability analysis of systems with small disturbances. Figure 2 illustrates a multivariate negative feedback system in the xy coordinate system, with the device-side conductance matrix Z _GFM-xy (s) as the forward transfer function and the inverse matrix of the grid-side conductance matrix Z _sys-xy (s) as the feedback transfer function.

The grid-side conductance matrix in the xy coordinate system is given directly by Wang et al. (2023): Z sys − x y s = Z β s − α s α s β s (15) where α s = ω 0 2 / s 2 + ω 0 2 , β s = s ω 0 / s 2 + ω 0 2

According to Section 1, the state space equations containing the device-side parameters are built in dq coordinate frame. Therefore, to obtain the device-side impedance matrix in the xy coordinate frame, it is necessary to derive the impedance matrix from the state space equations, and then a coordinate transformation will be employed to gain the impedance matrix under the xy coordinate system.

A subsystem state space containing only the converter state variables can be obtained by performing the Laplace transform on the equations and variables of the converter control and primary circuit.

The state space equations containing only the dynamics of the converter can be obtained by transforming Equation 13, as shown in Equation 16: s E 10 × 10 Δ x = A 10 × 10 Δ x + B 10 × 2 u y = C 2 × 10 Δ x (16) where s is complex variable, A _10×10 matrix composed of the remaining elements after removing A (3:4,1:12) and A (1:12,3:4) from A , y = [ Δ V_d, Δ Vq] T , u = [ Δ Igd, Δ I _gq ]^T, E _10×10 is the matrix composed of the remaining elements after removing E (3:4,1:12) and E (1:12,3:4) from E , B 3 : 4 , 1 : 2 = − 1 C f 1 0 0 1 , C 1 : 2 , 3 : 4 = 1 0 0 1 . The other elements of B and C are 0. Equation 14 can derive the device-side impedance matrix in the dq-axis as shown in Equation 17: Δ V d Δ V q = Z GFM − d q s Δ I g d Δ I g q (17) where Z _GFM-dq (s) = C _2×10(s E _10×10- A _10×10)⁻¹ B _10×2 represents the impedance matrix of the converter.

Transforming the Equation 15 to the xy coordinate frame, as detailed in Supplementary Appendix B. Finally, it can be obtained by Equation 18: Δ V x Δ V y = Z GFM − x y s Δ I g x Δ I g y (18) where Z _GFM-xy is the impedance matrix in the xy coordinate system. Further partitioning the matrix yields can be obtained by Equation 19: Z GFM − x y s = − 1 S B Z x x s Z x y s Z y x s Z y y s (19) where Z _xx (s), Z _xy (s), Z _yx (s) and Z _yy (s) are block matrices for Z _GFM-xy (s).

Thus, the characteristic equation of the closed-loop system shown in Figure 2 can be described as follows: det I + Z GFM − x y Z sys − x y − 1 = 0 (20) where I is a 2 × 2 identity matrix, det (∙) represents finding the determinant of a matrix.

Substituting Equation 14 into Equation 20 yields: c s = det Z x x s Z x y s Z y x s Z y y s + 1 SCR β s − α s α s β s (21)

Equation 21 can be further represented with the explicit function of SCR: SCR 2 + a s SCR + b s = 0 (22) where, a s = β s Z x x s + Z y y s + α s Z x y s − Z y x s Z x x s Z y y s − Z x y s Z y x s b s = α 2 s + β 2 s Z x x s Z y y s − Z x y s Z y x s (23)

The solution of Equation 22 is the eigenvalue of the system. The small-signal stability of the system can be judged by whether the eigenvalues of the system exist in the right half-plane. When a right half-plane eigenvalue exists, the system is small destabilized. For the GFM device, increasing SCR leads to a gradual move of the dominant eigenvalue from the left half-plane to the right half-plane, and the SCR when the dominant eigenvalue is near the imaginary axis is defined as CSCR.

In addition, the coefficients a(s) and b(s) in Equation 22 contain the elements of the device-side impedance matrix with the converter control dynamics, as shown in Equation 23. This indicates that the converter control dynamics affect the solution of Equation 22 and thus the CSCR of the system.

Figure 3a shows the variation of the dominant eigenvalue with increasing SCR for a specific case. The dominant eigenvalues gradually move to the right half plane as the SCR increases, eventually reaching a critical steady state at SCR = 21.5. Besides, the relationship between SCR and damping ratio is shown in Figure 3b. As SCR increases from 10 to 50, the damp ratio gradually decreases from 0.16 to a negative value, indicating a continuous decrease in the system stability margin.

Based on the theory described in Section 2, the data set is constructed with the damping coefficient, inertia coefficient of GFM, and PI control parameters of inner and outer loops as the input, and the system CSCR as the output. The specific method is as follows:

Eigenvalue extraction: Selecting data that enables the system to operate in steady conditions to construct an input dataset. Selecting one of the input data sets and computing the eigenvalue s _m according to Equation 22.

After completing the establishment of the dataset, it is also necessary to select appropriate models and algorithms to process the data. The principles of the support vector regression model and the particle swarm optimization algorithm will be introduced next.

The SVR model is an extended form of the SVM in regression problems, having a good effect in dealing with small sample, nonlinear, and high-dimensional problems (Abo-Khalil and Lee, 2008).

Let the training set be x k , y k k = 1 n , where x k ∈ R z represents the input vector of the kth training sample and y k ∈ R represents the corresponding output vector. Transforming the input vector to a high-dimensional feature space by nonlinear mapping, the following nonlinear regression function can be established as follows: f x = μ T φ x + b (24) where φ x is the nonlinear mapping function, μ is the weight vector, and b is the bias term. By introducing relaxation variables ξ k and ξ k * , the original optimization problem can be expressed as: min 1 2 μ 2 + C ∑ k = 1 n ξ k + ξ k * s . t . y k − μ T φ x k − b ≤ ε + ξ k μ T φ x k + b − y k ≤ ε + ξ k * ξ k , ξ k * ≥ 0 (25) where C is the penalty factor, and ε represents the error of the insensitive loss function.

Kernel functions can directly compute the result of the inner product of two vectors in a high-dimensional space, without performing the mapping process explicitly. The radial basis kernel function is commonly used due to its wide applicability, whose specific expression is shown in Equation 26: K x k , x j = exp − γ x k − x j 2 (26) where γ represents the kernel function parameter.

Introducing the Lagrange multiplier method as well as the Karush Kuhn Tucker (KKT) condition and combining it with the kernel function shown in Equations 23, 24 can be transformed to construct the following dyadic problem, which is shown as Equation 27: max − 1 2 ∑ k = 1 n ∑ j = 1 n α k − α k * α j − α j * K x k , x j + ∑ k = 1 n α k * − α k y k − ε ∑ k = 1 n α k * + α k s . t . ∑ k = 1 n α k − α k * = 0 0 ≤ α k , α k * ≤ C (27) where α k , α k * , α j , a n d α j * are the Lagrange coefficients.

After solving the above optimization problem, the final regression function can be obtained as Equation 28: f x = ∑ k = 1 n α k − α k * K x k , x + b (28)

The performance of the SVR model is affected by the penalty factor C, allowable error ε, and radial basis kernel function parameters γ. To improve the prediction accuracy of the model, the Particle Swarm Optimization algorithm is employed to optimize the combination of each parameter (Xu et al., 2025).

PSO was first proposed by Kennedy and Eberhart in 1995 (Kennedy and Eberhart, 1995). Each potential solution is represented as a ‘particle’ in the search space by simulating the behavior of bird flocks searching for food together, which dynamically updates its position and speed according to its own historical optimal experience and the flock’s optimal experience. The fitness function serves as the evaluation index to finally converge to the global optimal solution.

During the iteration process, the velocity and position of each particle will be updated and adjusted based on its historical optimal position and the global optimal position of the entire population. The updated rules for particles are shown in Equation 29: v m o h + 1 = β v m o h + c 1 r 1 p m o − l m o h + c 2 r 2 g o − l m o h l m o h + 1 = l m o h + v m o h + 1 (29) where m represents the mth particle, o is the dimension of the particle, h is the number of iterations, β is the inertia weight, v m o h and l m o h are the velocity and position of the mth particle, respectively, c _1, and c ₂ are the learning factors, r ₁ and r ₂ are the random numbers within the range of [0,1], p _mo and g _o are the optimal position of the mth particle and the population, respectively. Additionally, to avoid overfitting and improve the model’s generalization ability when optimizing parameter selection, this paper adopts a 5-fold cross-validation method.

In summary, the overall flowchart of the PSO-SVR prediction model is as Figure 4.

To verify the effectiveness of establishing a GFM converter model based on the state space method in this paper, a time-domain simulation model was built on the MATLAB/Simulink platform. The model parameters are shown in Supplementary Table S1.

The active power output of the converter obtained from the Simulink time-domain model is shown in Figure 6, and the dominant oscillation frequency obtained from this waveform is 3.731 Hz. Based on the established state space model, the system characteristic value is obtained, and the dominant characteristic value is 0.581 + j23.9, corresponding to an oscillation frequency of 23.9/(2π) = 3.803 Hz. The two oscillation frequencies are the same, which demonstrates the effectiveness of the state space model established in this paper.

The accuracy of the model in this paper can be verified by comparing the time-domain simulation results under the same parameters with the CSCR results calculated by Equation 19. Adjusting the parameters of the time-domain simulation model to the critical steady state, the oscillation frequency of the output active power P is 3.14 Hz.

The dominant eigenvalue is −0.328 + j19.81 by solving (Equation 20), corresponding to an oscillation frequency of 3.15 Hz, which is basically the same as that of the time-domain simulation. This proves the validity of the GFM state-space model described in this paper.

To verify the feasibility of the evaluation method described in this paper, two cases are given: case 1: SCR = 5; case 2: SCR = 10.

The GFM control parameters are the same as in 5.1. Applying a small disturbance of 0.1 p.u. to the grid voltage E at 0.5 s, and the active power P under the two cases is obtained as shown in Figure 10. The system under case 1 is small-signal stabilized, while the system under case 2 is small-signal destabilized. According to the small-signal stability assessment process in Section 4, combined with the PSO-SVR prediction model proposed in this paper, the corresponding CSCR of the control parameters is 6.26. Thus, σ = 1.26 > 0 in case 1 and σ = −3.74 < 0 in case 2, indicating that case 1 is small-signal stabilized while case 2 is small-signal destabilized, which is in accordance with the phenomenon shown in Figure 10.