Analysis of Grid-Connected Stability of VSG-Controlled PV Plant Integrated with Energy Storage System and Optimization of Control Parameters

: In the static stability analysis of the grid-connected photovoltaic (PV) generation and energy storage (ES) system, the grid-side is often simplified using an infinite busbar equivalent, which streamlines the analysis but neglects the dynamic characteristics of the grid, leading to certain inaccuracies in the results. Furthermore, the control parameter design does not consider the coupling relationships among parameters, resulting in arbitrary values and the inability to achieve overall optimality. To address these issues, this paper presents a comprehensive parameter optimization method for the oscillation characteristics of grid-connected PV generation and ES systems in various frequency ranges. Firstly, a detailed modeling of the grid-connected PV generation and ES system is conducted, resulting in the derivation of the system’s small-signal model. This study investigates the impact of parameters related to PV arrays, ES units, and the virtual synchronous generator (VSG) on the system’s characteristic roots using participation factors, sensitivity analysis, and eigenvalue root trajectories. Subsequently, based on the analysis of system root trajectories and the Particle Swarm Optimization (PSO) algorithm, a holistic parameter optimization design method for the system’s oscillation modes is proposed, and the parameter optimization results are obtained. Finally, the accuracy of the theoretical analysis is validated through perturbation testing using both Matlab/Simulink models and the small-signal model.


Introduction
Renewable energy and distributed generation have garnered widespread attention and rapid development as effective solutions to mitigate greenhouse gas emissions and address energy crises [1].To ensure the power quality and flexible control of distributed energy sources, power electronic converters play a crucial role as interfaces between distributed generators and the utility grid, and their optimized control is essential for stable system operation [2].Based on the constructiveness of grid frequency, power electronic converters can be categorized into Grid-Following (GFL) and Grid-Forming (GFM) synchronous control.Presently, most grid-connected distributed energy sources adopt GFL control based on phase-locked loops, which enable maximum power tracking [3] and exhibit high efficiency and stability in conventional strong power grids with high short circuit ratios.However, in weak power grids with a high penetration of renewable energy sources, GFL control is no longer appropriate [4,5].Virtual synchronous generator (VSG) control, as a representative GFM synchronous control method, emulates the rotor motion and excitation characteristics of traditional synchronous generators (SGs) by introducing inertia, damping, and droop control in the frequency and voltage amplitude domains of gridconnected devices.VSG control holds the potential to address the issues of voltage, inertia, and damping support in systems where renewable and distributed generation cannot provide these essential grid attributes.Due to its adaptability to weak power grids, VSG research has been widely applied in energy storage (ES) integration [6,7] and photovoltaic (PV) generation [8].However, VSG grid-connected systems encounter power stability issues similar to SGs [9] and inherent harmonic oscillations stemming from interactions in the Direct Current (DC)-side, voltage source converter, and grid-end system, leading to harmonic oscillations in the hundreds to thousands of hertz range.Therefore, effective modeling and analysis methods are needed to address the stability of VSG grid-connected systems with high renewable energy penetration.
To analyze the influence of key parameters on VSG stability, Yan et al. [10] conducted a small-signal modeling of a VSG-connected resistance-reactance load system under islanding conditions, investigating the effects of line and load parameter variations on system stability.Chen [11] established a small-signal model for a microgrid containing a VSG, analyzing the impact of control and circuit parameters on system stability.However, their models had high order, making further analysis and practical application challenging.Ref. [12] suggested that when a single VSG is connected to the grid, if the line resistance and reactance are of similar magnitudes, active and reactive power become coupled, leading to oscillations.Jiang et al. [13] used a state-space model to analyze the influence of inner loop proportional gain, integral coefficient, and virtual inductance parameters on harmonic oscillation characteristics.Xu et al. [14,15] improved dynamic characteristics by adding a differential loop in the VSG power loop but introduced high-frequency interference, affecting calculation accuracy.Liu et al. [16,17] established a state-space model for a multiinverter and interconnected system based on VSG algorithms, analyzing the contribution of various parameters to eigenvalue modes through sensitivity analysis.Zeng D et al. [18] analyzed the interaction between multiple VSGs and the grid, revealing that parallel oscillations originate from power oscillations but did not further analyze the power sources.The aforementioned references did not consider the dynamic characteristics of the gridend while conducting small-signal stability analysis, simplifying the analysis process but sacrificing accuracy.Additionally, they did not account for the coupling relationships among control parameters when analyzing them, relying solely on the impact of eigenvalue trajectories to design parameters, resulting in suboptimal parameter selection.
To address these issues, this paper proposes a global optimization method for control parameters related to system oscillation modes.First, based on the aforementioned literature, a detailed modeling of various components in a photovoltaic-storage power generation system controlled by a VSG is conducted, resulting in the derivation of the system's small-signal model.By combining system participation factors, sensitivity analysis, and eigenvalue root trajectory analysis, the properties of different oscillation modes of the system are determined, and the parameters to be optimized, along with their corresponding solution spaces, are identified.Subsequently, a global optimization method for controller parameters that enhance system stability is presented, and the accuracy of the theoretical analysis is validated through simulations.

Modeling of PV Generation and ES System
The structure of the PV and ES system is depicted in Figure 1.The PV array and the battery are connected to the DC bus through a Boost converter and a bidirectional DC-DC converter, respectively.The inverter power source is connected to the Point of Common Coupling (PCC) with a pure power load and an SG.The PV array employs maximum power point tracking (MPPT) control to ensure it operates near its maximum power point (MPP), even in varying external conditions.The ES unit, through constant voltage bus control, serves to smooth out the energy fluctuations of the photovoltaic generation system, stabilize the DC bus voltage, and facilitate bidirectional energy flow.The inverter utilizes VSG control, endowing the grid-connected device with characteristics such as rotational inertia and damping components, significantly enhancing the power output and frequency disturbance resistance of PV and ES units.Compared to traditional photovoltaic and storage grid-connected systems, a photovoltaic-storage power generation system based on VSG control possesses rotational inertia and damping sharing features, greatly improving the output power and frequency disturbance resistance of PV and ES units.In Figure 1, when neglecting line losses, the power balance relationship in the system can be expressed as Equation (1): where ppv is the power output from the PV array, (positive); pbat is the power absorbed or released by the ES unit (positive when discharging, negative when charging); pinv is the power output from the inverter, which varies with the grid power, exporting or absorbing power; pg is the grid-connected power; Pload is the load and line loss power.All powers except for Pload are considered positive in the direction of output.
The PV and ES system based on VSG control, as shown in Figure 2, consists of the DC-side, including PV and ES converter hardware and control components, and the Alternating Current (AC)-side, which comprises VSG control and grid-equivalent components.The DC-side is modeled in the time domain, while the AC-side is modeled using dq coordinates.Finally, by coordinate transformation, a unified reference frame is achieved, allowing for the linearization of both DC and AC models and the construction of the small-signal model for the entire system.The following sections will analyze these components separately.In Figure 1, when neglecting line losses, the power balance relationship in the system can be expressed as Equation (1): where p pv is the power output from the PV array, (positive); p bat is the power absorbed or released by the ES unit (positive when discharging, negative when charging); p inv is the power output from the inverter, which varies with the grid power, exporting or absorbing power; p g is the grid-connected power; P load is the load and line loss power.All powers except for P load are considered positive in the direction of output.
The PV and ES system based on VSG control, as shown in Figure 2, consists of the DCside, including PV and ES converter hardware and control components, and the Alternating Current (AC)-side, which comprises VSG control and grid-equivalent components.The DC-side is modeled in the time domain, while the AC-side is modeled using dq coordinates.Finally, by coordinate transformation, a unified reference frame is achieved, allowing for the linearization of both DC and AC models and the construction of the small-signal model for the entire system.The following sections will analyze these components separately.
frequency disturbance resistance of PV and ES units.Compared to traditional photovoltaic and storage grid-connected systems, a photovoltaic-storage power generation system based on VSG control possesses rotational inertia and damping sharing features, greatly improving the output power and frequency disturbance resistance of PV and ES units.In Figure 1, when neglecting line losses, the power balance relationship in the system can be expressed as Equation (1): where ppv is the power output from the PV array, (positive); pbat is the power absorbed or released by the ES unit (positive when discharging, negative when charging); pinv is the power output from the inverter, which varies with the grid power, exporting or absorbing power; pg is the grid-connected power; Pload is the load and line loss power.All powers except for Pload are considered positive in the direction of output.
The PV and ES system based on VSG control, as shown in Figure 2, consists of the DC-side, including PV and ES converter hardware and control components, and the Alternating Current (AC)-side, which comprises VSG control and grid-equivalent components.The DC-side is modeled in the time domain, while the AC-side is modeled using dq coordinates.Finally, by coordinate transformation, a unified reference frame is achieved, allowing for the linearization of both DC and AC models and the construction of the small-signal model for the entire system.The following sections will analyze these components separately.

Modeling of DC-Side Power Generation System and Its Controller
The mathematical model for the main circuit of the PV generation system encompasses components such as the PV array output voltage, stabilizing capacitor, L-filter, and Boost step-up circuit.To describe the behavior of the PV side, this study employs an averagevalue model.This modeling approach yields the equation set presented in Equation (2), where d p is the duty cycle controlling the PV converter switch.
To improve the efficiency of solar energy utilization, the tracking of MPPT is achieved by controlling the on/off states of the Boost circuit's switching transistors.The control diagram is depicted in Figure 2, where the PV output voltage u pv is compared to the MPP voltage u * pv determined by the MPPT algorithm.Through the utilization of a Proportional-Integral (PI) controller, the duty cycle d p is calculated and subsequently compared to a carrier signal to generate Pulse Width Modulation (PWM) triggering pulses.As environmental factors such as sunlight intensity and temperature fluctuate, the error value also changes accordingly, resulting in adjustments to the duty cycle of the triggering pulses.This dynamic process enables the system to effectively track MPP, as demonstrated in Equation (3), where K p1 and K i1 represent the proportional and integral coefficients of MPPT control, and x 1 denotes the state variable within the PI controller.
ES and generation systems exhibit two topological configurations: In the Buck operating mode, energy flows from the DC bus to the ES, treating the ES as a load integrated into the system.Conversely, in the Boost operating mode, energy travels from the ES towards the DC bus, with the ES acting as a power source.To describe the behavior of the ES rectification side, this study utilizes an average-value model, resulting in the equation set outlined in Equation ( 4), where d b signifies the duty cycle governing the ES rectifier.
To enable the functionality of peak shaving and valley filling in the ES, a dual closedloop control strategy for voltage regulation is employed.The control diagram, as depicted in Figure 2, operates as follows: The deviation between the reference voltage u * dc and the actual voltage u dc on the DC bus is processed through a voltage loop PI controller, resulting in a reference current i * bat , which is subsequently compared to the ES output current.After passing through a current loop PI controller, it generates the desired duty cycle d b .Consequently, two sets of PWM pulse trains are generated by comparing with the carrier signal, controlling the switching transistors.This control process is represented by Equation (5), where K p2,3 and K i2,3 represent the proportional and integral coefficients governing constant voltage control on the bus, while x 2,3 denotes the state variable within the dual PI control loop.
Electronics 2024, 13, 1343 5 of 18 The DC port interfaces the PV generation and ES system with the converter's DCside, ensuring a constant voltage supply to achieve energy exchange.According to the topological structure depicted in Figure 2, the equation set can be outlined in Equation ( 6) as follows: where p dc is the instantaneous power output on the DC-side; p o is the instantaneous output power from the inverter; and u odq and i odq represent the voltage and current outputs from the inverter, respectively, in the dq coordinate system.

Modeling of AC-Side Virtual Synchronous Machine and Its Controller
The voltage-controlled ES VSG grid-connected system, as depicted in Figure 2, is subjected to small-signal modeling following procedures detailed in Yan et al. [10], Chen [11], and Xu et al. [14].For this modeling, the dq rotating coordinate system was selected as the reference coordinate system.This small-signal modeling encompasses five key components: power calculation, power control, virtual impedance control, dual-loop voltage-current control, and the grid interface circuit.The power calculation component employs a low-pass filter to compute the average output power of the inverter.In the power control module, electrical characteristics resembling those of an SG are emulated, enabling functions such as inertia support, primary frequency regulation, and reactive power-voltage regulation.Virtual impedance control introduces stator electrical characteristics of a synchronous machine into the control loop by simulating properties like synchronous reactance and armature resistance.The dual-loop voltage-current control module generates voltage reference signals and produces PWM drive signals to control various switching devices.Differential equations governing these five components are presented in Equations ( 7)- (11).
where p e and q e are the active and reactive power output of the VSG; ω c is the cutoff frequency of the low-pass filter; P ref is the active power command; J is the virtual rotational inertia; D p is the virtual damping factor of frequency; ω and ω 0 are the virtual angular frequency and rated angular frequency; K DP is the active frequency modulation coefficient; Q ref is the reactive power command value; U n is the rated voltage amplitude; K q is the reactive integral coefficient; D q is the coefficient governing reactive voltage regulation; e d and e q are the d-and q-axis components of the internal potential e, respectively; u od , u oq , i od , and i oq are the dq-axis components of the VSG output voltage and current, respectively; u * od and u * oq are the voltage reference values provided by the virtual impedance module for voltage and current closed-loop control; r v and L v are the virtual resistance and virtual inductance, respectively; i * ld and i * lq are the reference values from the voltage loop for the current loop; u * d and u * q are the modulation signals controlled by SPWM; K iv , K ic , K pv , and K pc are the integral and proportional coefficients employed in the dual closed-loop voltage and current control; Φ d and Φ q are the intermediate state variables in the PI control of the voltage loop; γ d and γ q are the intermediate state variables in the PI control of current loops; u gd and u gq are the components of the SG port voltage on the dq axis, respectively.

Network Terminal Equivalence
Both VSGs and SGs exhibit similar characteristics, resulting in stability challenges when small disturbances occur in grid-connected systems.Furthermore, as a consequence of the substantial integration of new energy sources into the grid, significant alterations are observed in the frequency characteristics at the grid interface.Unlike the scenario where the VSG is connected to an infinite system, as illustrated in Figure 2, the grid-connected system comprises a VSG, a constant power load, and an SG.The power load is integrated into the rotor motion equations, and the generator model's excitation system incorporates a first-order inertia excitation element.In the dq coordinate system, the expression is as depicted in Equation (12).
where ω sg and δ sg are the angular velocity of the SG and the angle between the q-axis of the SG's rotor and the x-axis of the xy coordinate system, respectively; u gd and u gq are the d-and q-axis voltage components of the SG, respectively; E ′ q and E f are the transient electromotive force and the excitation voltage output by the excitation system, respectively; X d and X ′ d are the steady-state and transient values of the direct-axis reactance; X q is the steady-state value of the quadrature-axis reactance; T sg J and D sg p are the generator's inertia time constant and deviation, respectively; T a and K a are the time constant and proportional coefficient of the excitation system; P msg and P esg are the mechanical and electromagnetic power of SG, respectively; U Nsg is the set values of the excitation controller output voltage of the SG.

Small-Signal Modeling of PV Generation and ES VSG Grid-Connected Systems
The analysis of interaction mechanisms in nonlinear systems is commonly approached through small-signal modeling methods [19].Among these methods, the state-space approach allows for a detailed examination of the changes in various state variables within the system.Through the incorporation of participation factors, sensitivity analysis, and root locus analysis, it becomes possible to determine the characteristics of different oscillatory modes within the system.Moreover, it facilitates the identification of parameters responsible for these modes, offering a holistic assessment of system stability.This study integrates the rectification systems of PV generation and ES, the VSG control system, and the grid interface equivalent portion to construct a full-order state-space model for GFM PV generation and ES systems.Equations ( 2)-( 11) collectively encompass 21 state variables and 21 differential-algebraic equations, effectively capturing the dynamic characteristics of the system.By linearizing this model around its steady-state operating point and organizing the equations, a comprehensive model is derived for transient stability analysis.
In Section 2, a model for the photovoltaic (PV) and energy storage grid-connected system was developed.This section will develop and validate the small-signal model of the PV generation and ES VSG grid-connected system.The analysis of interaction mechanisms in nonlinear systems is commonly approached through small-signal modeling methods [19].Among these methods, the state-space approach allows for a detailed examination of the changes in various state variables within the system.Through the incorporation of participation factors, sensitivity analysis, and root locus analysis, it becomes possible to determine the characteristics of different oscillatory modes within the system.Moreover, it facilitates the identification of parameters responsible for these modes, offering a holistic assessment of system stability.This study integrates the systems of PV generation and ES, the VSG control system, and the network terminal equivalence to construct a full-order state-space model for GFM PV generation and ES systems.Equations ( 2)-( 11) collectively encompass 21 state variables and 21 differential-algebraic equations, effectively capturing the dynamic characteristics of the system.By linearizing this model around its steadystate operating point and organizing the equations, a comprehensive model is derived for transient stability analysis: where ∆x = [∆x DC ∆x AC ] T is the system's state vector, as illustrated in Equation ( 14), which encompasses state vectors from both DC-and AC-sides; ∆u = [∆i pv ∆i bat ∆u gdq ] T is the system's input vector.The calculation of the eigenvalues of the state matrix A enables the analysis of the grid-connected system's static stability and frequency oscillation characteristics.

Small-Signal Model Validation
In the state-space model and the Simulink time domain dynamic model, step disturbances were applied to the grid frequency, grid voltage, and power commands within the 2-4 s timeframe.∆ω g = 0.6 π•rad/s, ∆u g = 0.3 V, and ∆P ref = 500 W illustrate the comparison of the active power P o and reactive power Q o output curves of the VSG between both models.In Figure 3, the waveforms of the small-signal model closely match those of the practical simulation model, thereby validating the accuracy of the small-signal model developed in this research.
ances were applied to the grid frequency, grid voltage, and power commands within the 2-4 s timeframe.Δωg = 0.6 π•rad/s, Δug = 0.3 V, and ΔPref = 500 W illustrate the comparison of the active power Po and reactive power Qo output curves of the VSG between both models.In Figure 3, the waveforms of the small-signal model closely match those of the practical simulation model, thereby validating the accuracy of the small-signal model developed in this research.

Small-Signal Stability Analysis
The small-signal model established in Section 3 has been validated.This section, as a case study, performs a static stability analysis of the PV generation and ES system, as illustrated in Figure 2. The parameters for the numerical model are detailed in Table 1.Subsequently, these parameters and initial values are input into the system's state matrix A for the computation of all eigenvalues, as displayed in Table 2. Eigenvalues λ1-4 and λ9-14 are complex, corresponding to 10 oscillatory modes.The real components reflect the system's damping concerning oscillations, while the imaginary components indicate the oscillation frequencies.All the system's eigenvalues are in the complex plane's left half, signifying that the system exhibits small disturbance stability.

State Matrix Eigenvalue Calculation and Stability Analysis 4.1. Small-Signal Stability Analysis
The small-signal model established in Section 3 has been validated.This section, as a case study, performs a static stability analysis of the PV generation and ES system, as illustrated in Figure 2. The parameters for the numerical model are detailed in Table 1.Subsequently, these parameters and initial values are input into the system's state matrix A for the computation of all eigenvalues, as displayed in Table 2. Eigenvalues λ 1-4 and λ 9-14 are complex, corresponding to 10 oscillatory modes.The real components reflect the system's damping concerning oscillations, while the imaginary components indicate the oscillation frequencies.All the system's eigenvalues are in the complex plane's left half, signifying that the system exhibits small disturbance stability.

Parameter
Value Parameter Value  To explore the relationship between system state variables and the associated oscillatory modes of eigenvalues, participation factors for each mode were calculated, and the results are depicted in Figure 4. To explore the relationship between system state variables and the associated oscillatory modes of eigenvalues, participation factors for each mode were calculated, and the results are depicted in Figure 4.

Eigenvalue Sensitivity Analysis
Several factors can introduce errors in computed eigenvalues, including practical data inaccuracies and limitations in calculation precision.Additionally, due to the interdependencies among various parameters within the system, it becomes challenging to isolate the influence of individual parameter variations on oscillatory modes.Therefore, conducting eigenvalue sensitivity analysis becomes highly valuable.On the DC-side, six parameters were selected, denoted as K p1-3 and K i1-3 .Meanwhile, on the AC-side, eight parameters were chosen, including ω c , J, D p , K pv , K iv , K pc , K ic , and L v , to investigate the sensitivity of eigenvalues concerning system parameters.The outcomes of this sensitivity analysis are illustrated in Figure 5.
The analysis of the 21 characteristic eigenvalues, as depicted in Figures 4 and 5, offers significant insights into the stability and sensitivity of the PV generation and ES system.Specifically, the eigenvalues are divided into two categories, with λ 1-8 representing characteristics on the DC-side and λ 9-21 pertaining to the AC-side.
On the DC-side, λ 1,2 are closely associated with i Lpv and u pv , with substantial impacts from parameters K p1 and K i1 .λ 3,4 predominantly relate to u dc and i Lbat , primarily influenced by parameters K p2 and K i2 .Furthermore, λ 5 demonstrates a strong correlation with ∆x 2 , indicating its sensitivity to the outer loop PI control parameters governing the ES system.Lastly, λ 7,8 are associated with state variables ∆x 3 and ∆x 1 , emphasizing the significance of designing three appropriate sets of PI parameters for influencing the oscillatory behavior on the DC-side.
late the influence of individual parameter variations on oscillatory modes.Therefore, conducting eigenvalue sensitivity analysis becomes highly valuable.On the DC-side, six parameters were selected, denoted as Kp1-3 and Ki1-3.Meanwhile, on the AC-side, eight parameters were chosen, including ωc, J, Dp, Kpv, Kiv, Kpc, Kic, and Lv, to investigate the sensitivity of eigenvalues concerning system parameters.The outcomes of this sensitivity analysis are illustrated in Figure 5.The analysis of the 21 characteristic eigenvalues, as depicted in Figures 4 and 5, offers significant insights into the stability and sensitivity of the PV generation and ES system.Specifically, the eigenvalues are divided into two categories, with λ1-8 representing characteristics on the DC-side and λ9-21 pertaining to the AC-side.
On the DC-side, λ1,2 are closely associated with iLpv and upv, with substantial impacts from parameters Kp1 and Ki1.λ3,4 predominantly relate to udc and iLbat, primarily influenced by parameters Kp2 and Ki2.Furthermore, λ5 demonstrates a strong correlation with Δx2, indicating its sensitivity to the outer loop PI control parameters governing the ES system.Lastly, λ7,8 are associated with state variables Δx3 and Δx1, emphasizing the significance of designing three appropriate sets of PI parameters for influencing the oscillatory behavior on the DC-side.
Conversely, on the AC-side, λ9-12 are primarily linked to the inductance current dqaxis components Δildq and capacitance voltage dq-axis components Δuodq, revealing their dependence on the selection of LC filter parameters.λ13,14 correspond to grid-connected current dq-axis components Δiodq and are predominantly influenced by LC filter parameters, virtual impedance, and line impedance.Additionally, λ15,16 are related to active power ΔP and reactive power ΔQ, exhibiting sensitivity to the ωc parameter in the low-pass filter.Notably, λ17 is mainly associated with Δω, ΔP, and ΔQ, with its value influenced by rotor inertia and damping coefficients (Kf/Jωn+Dp/J).Furthermore, λ18,19 are linked to the state variables of PI regulators' Δγdq in the current control loop, with values approximating Kiv/Kpv.Finally, λ20,21 are primarily associated with the state variables of PI regulators' Δxdq in the voltage control loop, with values approximating Kic/Kpc.

Eigenvalue Trajectory Analysis
In the context of GFM PV generation and ES systems, it is essential to analyze and optimize control parameters since circuit parameters are typically pre-designed based on practical requirements and are not frequently modified due to component constraints.Conversely, on the AC-side, λ 9-12 are primarily linked to the inductance current dqaxis components ∆i ldq and capacitance voltage dq-axis components ∆u odq , revealing their dependence on the selection of LC filter parameters.λ 13,14 correspond to grid-connected current dq-axis components ∆i odq and are predominantly influenced by LC filter ters, virtual impedance, and line impedance.Additionally, λ 15,16 are related to active power ∆P and reactive power ∆Q, exhibiting sensitivity to the ω c parameter in the low-pass filter.Notably, λ 17 is mainly associated with ∆ω, ∆P, and ∆Q, with its value influenced by rotor inertia and damping coefficients (K f /Jω n +D p /J).Furthermore, λ 18,19 are linked to the state variables of PI regulators' ∆γ dq in the current control loop, with values approximating K iv /K pv .Finally, λ 20,21 are primarily associated with the state variables of PI regulators' ∆x dq in the voltage control loop, with values approximating K ic /K pc .

Eigenvalue Trajectory Analysis
In the context of GFM PV generation and ES systems, it is essential to analyze and optimize control parameters since circuit parameters are typically pre-designed based on practical requirements and are not frequently modified due to component constraints.
In this regard, the key control parameters with a significant influence on system stability are categorized into two main aspects: the DC-side, encompassing control coefficients K p1 and K i1 for the PV array and the dual-loop control coefficients K p2,3 and K i2,3 for the ES constant voltage bus; and the AC-side, including the cutoff angular frequency ω c , rotor inertia J, active power frequency control coefficient K Dp , damping coefficient D p , reactive power voltage control coefficient D q , virtual impedance r v , virtual inductance L v , and dual-loop PI control loop coefficients K pv , K iv , K pc , and K ic .It is noteworthy that the values of K Dp and D q should satisfy the requirements of virtual synchronous machine frequency regulation and voltage droop characteristics, typically requiring no modifications [20].This study primarily focuses on the analysis and optimization of the following parameters: To investigate the range of parameter values that ensure system stability and to understand the patterns of eigenvalue variations as parameters change, trajectory plots of relevant eigenvalues for both DC and AC systems are generated, as shown in Figures 6 and 7.It is important to note that due to the large span of root trajectories, eigenvalues with small variations are not displayed in the figures for clarity.
Figure 6a,b depict the eigenvalue trajectories with varying parameters K p1 and K i1 .As K p1 ranges from 0 to 100, eigenvalues λ 1,2 and λ 6 move towards the unstable region.When K i1 varies from 0 to 0.03, eigenvalues λ 1,2 shift to the right, entering the positive half-plane, indicating increased system instability.This suggests that setting PI control coefficients for the PV part too high is unfavorable for system stability.
In Figure 6c,d, the eigenvalue trajectories are shown for parameters K p2 and K i2 .As K p2 changes from 0.05 to 20, eigenvalues λ 3,4 transition from the positive half-plane to the negative half-plane, and λ 5 moves to the right.When K i2 varies from 0 to 100, eigenvalues λ 3,4 shift to the right, entering the positive half-plane, indicating reduced system stability.
Figure 6d,e illustrate the eigenvalue trajectories with changes in parameters K p3 and K i3 .When K p3 varies from 0 to 11, eigenvalues λ 3,4 move leftward after entering the positive half-plane, and λ 7 shifts towards the unstable region.Changes in K i3 from 0 to 80 do not significantly affect the system.
Figure 7a reveals that when the parameter ω c varies from 1 to 22 π•rad/s, it results in simultaneous changes in four eigenvalues.Specifically, eigenvalues λ 13,14 exhibit an increasing real part and decreasing damping ratio, which leads them into the unstable region.Meanwhile, λ 16,17 display a decreasing real part as they gradually move towards the stable region.In this regard, the key control parameters with a significant influence on system stability are categorized into two main aspects: the DC-side, encompassing control coefficients Kp1 and Ki1 for the PV array and the dual-loop control coefficients Kp2,3 and Ki2,3 for the ES constant voltage bus; and the AC-side, including the cutoff angular frequency ωc, rotor inertia J, active power frequency control coefficient KDp, damping coefficient Dp, reactive power voltage control coefficient Dq, virtual impedance rv, virtual inductance Lv, and dual-loop PI control loop coefficients Kpv, Kiv, Kpc, and Kic.It is noteworthy that the values of KDp and Dq should satisfy the requirements of virtual synchronous machine frequency regulation and voltage droop characteristics, typically requiring no modifications [20].This study primarily focuses on the analysis and optimization of the following parameters: Kp1-3, Ki1-3, ωc, J, Dp, rv, Lv, Kpv, Kiv, Kpc, and Kic.
To investigate the range of parameter values that ensure system stability and to understand the patterns of eigenvalue variations as parameters change, trajectory plots of relevant eigenvalues for both DC and AC systems are generated, as shown in Figures 6  and 7.It is important to note that due to the large span of root trajectories, eigenvalues with small variations are not displayed in the figures for clarity.In Figure 7b, parameter J undergoes variations from 0.01 kg•m 2 to 1.5 kg•m 2 , causing eigenvalue λ 17 to increase in the real part and gradually move towards the unstable region.While increasing virtual inertia can enhance the system's frequency support capability in an inertial environment [21], excessively high virtual inertia can hinder system stability [22].
Figure 7c illustrates the impact of parameter D p changing from 0.001 N•m•s/rad to 1 N•m•s/rad.It results in simultaneous changes in eigenvalues λ 15,16 .Both eigenvalues move towards the stable region, with an increasing damping effect.The oscillation frequency decreases as damping increases, contributing to a more stable system.Typically, increasing inertia J strengthens the system's frequency support ability but leads to larger overshoot and a longer settling time.Increasing damping reduces overshoot but may increase system error.Therefore, selecting appropriate values for virtual inertia and damping should consider practical factors [23][24][25].
Figure 7d shows that parameter L v varying from 0.2 mH to 20 mH leads to simultaneous changes in two eigenvalues.Both eigenvalue λ 13,14 exhibit a decreasing real part and an increasing damping ratio, moving towards the negative half-plane, which enhances system stability.Figure 6a,b depict the eigenvalue trajectories with varying parameters Kp1 and Ki1.As Kp1 ranges from 0 to 100, eigenvalues λ1,2 and λ6 move towards the unstable region.When Ki1 varies from 0 to 0.03, eigenvalues λ1,2 shift to the right, entering the positive half-plane, indicating increased system instability.This suggests that setting PI control coefficients for the PV part too high is unfavorable for system stability.
In Figure 6c,d, the eigenvalue trajectories are shown for parameters Kp2 and Ki2.As Kp2 changes from 0.05 to 20, eigenvalues λ3,4 transition from the positive half-plane to the negative half-plane, and λ5 moves to the right.When Ki2 varies from 0 to 100, eigenvalues λ3,4 shift to the right, entering the positive half-plane, indicating reduced system stability.
Figure 6d,e illustrate the eigenvalue trajectories with changes in parameters Kp3 and Ki3.When Kp3 varies from 0 to 11, eigenvalues λ3,4 move leftward after entering the positive half-plane, and λ7 shifts towards the unstable region.Changes in Ki3 from 0 to 80 do not significantly affect the system.In Figure 7e, when parameter K pv varies from 1 to 100, multiple eigenvalues undergo simultaneous changes.Eigenvalues λ 9-12 move towards the negative half-plane, indicating improved system stability.Eigenvalues λ 20,21 exhibit an increasing real part, moving towards the positive half-plane, reducing damping ratios, corresponding to increased oscillation frequencies and decreased system stability.
Figure 7f shows the effect of parameter K iv varying from 20 to 2500, resulting in simultaneous changes in multiple eigenvalues.Real parts of eigenvalues λ 9-12 decrease, gradually moving towards the unstable region, while eigenvalues λ 20,21 increase in the real part, moving towards the stable region.
Figure 7g illustrates that parameter K pc ranging from 14 to 100 causes simultaneous changes in four eigenvalues.Real parts of eigenvalues λ 9-12 decrease, gradually moving towards the negative half-plane.Specifically, λ 8,9 exhibit an increase in the absolute value of their imaginary parts, leading to increased oscillation frequencies, while λ 10,11 show less pronounced changes.When initial parameter values are too small, eigenvalues λ 8,9 fall into the positive half-plane, resulting in system instability.
In Figure 7h, parameter K ic varies from 20 to 2500, causing simultaneous changes in multiple eigenvalues.Eigenvalues λ 15,16 both move towards the stable region.However, real parts of eigenvalues λ 9-12 and λ 18,19 increase, all moving towards the unstable region, accompanied by reduced damping ratios, which are unfavorable for improving system stability.

Parameter Optimization Design for GFM PV ES System Parameters
To improve the stability and stability margin of the PV generation and ES system, Milano et al. [9] and Jiang et al. [13] employed the root locus method to analyze the impact of control parameter variations on system stability.However, root locus-based parameter optimization methods often exhibit subjectivity and lack well-defined, scientifically quantifiable design criteria.They struggled to address the intricate interactions among parameters and their combined effects on multiple eigenvalue characteristics within complex systems.Therefore, in the process of parameter tuning and selection, achieving an optimal position for system eigenvalues and identifying better parameter combinations require the application of alternative optimization methods.
The Particle Swarm Optimization (PSO) algorithm has several distinct advantages over other optimization methods such as genetic algorithms and neural networks.Firstly, PSO does not require complex operations like crossover and mutation, making it easier to code and implement.Secondly, it requires fewer computational steps, and the computational cost of updating particle positions is comparatively low, which improves the efficiency of PSO in solving certain optimization problems.Furthermore, compared to genetic algorithms, PSO requires fewer adjustable parameters, simplifying the application of the algorithm and the process of parameter tuning.Based on these advantages, building upon the general principles of parameter design, this paper utilizes a parameter optimization design approach that combines the root locus method with the PSO algorithm.Taking the root locus analysis in Section 4 as the value range for the parameters to be optimized, through iterations of PSO, the overall characteristic values move towards the negative half-axis and the damping ratio of each oscillation mode approaches the optimal damping ratio of 0.707, in order to increase the damping of system modes and enhance the transient stability of the system.The detailed steps of this approach are illustrated in Figure 8.The Particle Swarm Optimization (PSO) algorithm has several distinct advantages over other optimization methods such as genetic algorithms and neural networks.Firstly, PSO does not require complex operations like crossover and mutation, making it easier to code and implement.Secondly, it requires fewer computational steps, and the computational cost of updating particle positions is comparatively low, which improves the efficiency of PSO in solving certain optimization problems.Furthermore, compared to genetic algorithms, PSO requires fewer adjustable parameters, simplifying the application of the algorithm and the process of parameter tuning.Based on these advantages, building upon the general principles of parameter design, this paper utilizes a parameter optimization design approach that combines the root locus method with the PSO algorithm.Taking the root locus analysis in Section 4 as the value range for the parameters to be optimized, through iterations of PSO, the overall characteristic values move towards the negative half-axis and the damping ratio of each oscillation mode approaches the optimal damping ratio of 0.707, in order to increase the damping of system modes and enhance the transient stability of the system.The detailed steps of this approach are illustrated in Figure 8.

Parameter Optimization
The parameter ranges determined through root locus analysis in Section 4.3 serve as the solution space for the optimization algorithm.Initial velocities and positions for each particle are randomly generated, as shown in Equation ( 15

Parameter Optimization
The parameter ranges determined through root locus analysis in Section 4.3 serve as the solution space for the optimization algorithm.Initial velocities and positions for each particle are randomly generated, as shown in Equation (15), where x k n and v k n represent the positions and velocities of particles after the k-th iteration, and n represents the number of parameters to be optimized.


The fitness function p i for the i-th particle is defined as the basis for determining the target direction of iterations, evaluating individual performance, and assessing the global optimum, as shown in Equation (16).Understanding the impact of eigenvalues on the transient characteristics of the system, the overall objective of parameter optimization is to move the system's eigenvalues towards the negative half-plane to increase damping and shorten transient responses.Simultaneously, the aim is to achieve optimal damping ratios for all modes.With this goal in mind, two fitness functions p i1 and p i2 are designed separately.
where α and β are the weights assigned to the two functions, with values of 0.3 and 0.7, respectively.The expressions for p i1 and p i2 are detailed in Equation ( 17), where N denotes the number of eigenvalues to be optimized, and the real parts are arranged in descending order as δ 1 ~δN .The coefficient (N + 1 − j) represents the weight of the j-th real part.M stands for the number of modes, and ξ k represents the damping ratio of the k-th mode.
The objective of the optimization algorithm is to maximize the fitness function p i .The update formulas for velocity and position are shown in Equation (18), where ω represents the inertia weight, c 1 = c 2 = 2 are the cognitive factors, and r is a random number ranging between 0 and 1.During each iteration, each particle updates its position and velocity based on its individual optimal P b and the optimal population G b .
The parameters of the case study model are listed in Table 1.The overall optimization was performed on PI control parameters K p1-3 and K i1-3 , cutoff frequency ω c , virtual inertia J, damping coefficient D p , virtual inductance L v , and double-loop PI control parameters K pv , K iv , K pc , and K ic .The initial values and ranges for the parameters to be optimized are provided in Appendix A. A total of 500 iterations were set for the optimization process, and the iteration progress is illustrated in Figure 9.
After 500 iterations, the optimal population G b reached its maximum value, and the parameters to be optimized converged to their best positions.The results of the parameter optimization are presented in Table 3.
The data from Table 3 were used to calculate the optimized eigenvalues λ * i , which were then compared to the original eigenvalues λ i from Table 2, as shown in Figure 10.The eigenvalues before optimization and after optimization are represented by solid orange and black circles, respectively.In Figure 10, after optimization, the system eigenvalues have moved to varying degrees toward the negative half-plane.The damping ratios of the system modes have all improved.This enhancement in damping ratios signifies an improvement in the transient stability of the system.It also demonstrates the feasibility and effectiveness of the comprehensive parameter optimization approach.

Simulation Validation
To validate the effectiveness of the optimization algorithm, this study employed Matlab/Simulink to establish the simulation model depicted in Figure 1.Three operational scenarios were selected to analyze the output power of the PV array, ES, and VSG before and after parameter optimization.The power waveforms after optimization for these three scenarios are compared uniformly, as shown in Figure 11  In Figure 10, after optimization, the system eigenvalues have moved to varying degrees toward the negative half-plane.The damping ratios of the system modes have all improved.This enhancement in damping ratios signifies an improvement in the transient stability of the system.It also demonstrates the feasibility and effectiveness of the comprehensive parameter optimization approach.

Simulation Validation
To validate the effectiveness of the optimization algorithm, this study employed Matlab/Simulink to establish the simulation model depicted in Figure 1.Three operational scenarios were selected to analyze the output power of the PV array, ES, and VSG before and after parameter optimization.The power waveforms after optimization for these three scenarios are compared uniformly, as shown in Figure 11.Specifically, P pv , P * pv ; P bat , P * bat ; P vsg , and P * vsg represent the pre-and post-optimization output powers of the PV array, ES, and VSG.In Figure 11a, a disturbance is introduced by raising the active power command from 3 kW to 4 kW for 2-4 s.This increase in the command requires additional output power from the virtual machine and ES.After 4 s, the command value is restored, resulting in a reduced ES output, returning to the initial steady state.
In Figure 11b, the disturbance involves an increase in solar irradiance from 200 W/m 2 to 300 W/m 2 for 2-4 s.This leads to an increase in the PV array output power from 1000 W to 1500 W. To maintain power balance, the ES reduces its output from 2000 W to 1500 W. After 4 s, solar irradiance returns to its original level, causing the ES output to increase back to the initial steady state.
In Figure 11c, a minor disturbance is applied by decreasing the grid frequency from 50 Hz to 49.7 Hz for 2-4 s.This reduction in grid frequency requires an increase in the virtual machine output and additional power from the ES.After 4 s, the grid-side frequency is restored, leading to a decrease in the ES output back to the initial steady state.
The simulation results demonstrate that the controller parameters obtained through the optimization method result in smaller overshoot, a shorter settling time, and smaller steady-state error in the output powers of the PV array, ES, and VSG.This indicates that the GMF PV generation and ES system exhibits improved transient performance, an enhanced stability margin, and remains stable under minor disturbances, consistent with the dynamic analysis conducted earlier in this study.

Conclusions
This paper has developed a comprehensive small-signal model for the transient sta- In Figure 11a, a disturbance is introduced by raising the active power command from 3 kW to 4 kW for 2-4 s.This increase in the command requires additional output power from the virtual machine and ES.After 4 s, the command value is restored, resulting in a reduced ES output, returning to the initial steady state.
In Figure 11b, the disturbance involves an increase in solar irradiance from 200 W/m 2 to 300 W/m 2 for 2-4 s.This leads to an increase in the PV array output power from 1000 W to 1500 W. To maintain power balance, the ES reduces its output from 2000 W to 1500 W. After 4 s, solar irradiance returns to its original level, causing the ES output to increase back to the initial steady state.
In Figure 11c, a minor disturbance is applied by decreasing the grid frequency from 50 Hz to 49.7 Hz for 2-4 s.This reduction in grid frequency requires an increase in the virtual machine output and additional power from the ES.After 4 s, the grid-side frequency is restored, leading to a decrease in the ES output back to the initial steady state.
The simulation results demonstrate that the controller parameters obtained through the optimization method result in smaller overshoot, a shorter settling time, and smaller steady-state error in the output powers of the PV array, ES, and VSG.This indicates that the GMF PV generation and ES system exhibits improved transient performance, an enhanced stability margin, and remains stable under minor disturbances, consistent with the dynamic analysis conducted earlier in this study.

Conclusions
This paper has developed a comprehensive small-signal model for the transient stability analysis of the PV generation and ES system for VSG control.Combining participation factors, system sensitivity, and root locus analysis, a global optimization method for controller parameters has been proposed.The following conclusions have been drawn: (1) In the DC generation system, the PI control of the PV converter and the dual-loop control of the ES converter have a significant and comprehensive impact on the system's eigenvalues.Specifically, K i1 , K i2 , and K i3 can affect the system's damping if their values are too large; in the AC inverter system, when D p and L v remain constant or decrease, an increase in J will lead to a decrease in system stability.The selection of inner loop control parameters has a more complex impact on system eigenvalues.
(2) This study revealed the impact of various parameters on the system's stability from a small-signal stability perspective.The overall optimization method effectively overcomes the coupling characteristics between control parameters and oscillatory modes.It provides parameter optimization results and serves as a reference for parameter tuning in the integrated application of PV generation and ES with VSG control.These findings hold practical value in real-world applications.

Figure 2 .Figure 1 .
Figure 2. Block diagram of PV generation and ES system based on VSG control.

Figure 1 .
Figure 1.Topology of PV generation and ES system.

Figure 2 .Figure 2 .
Figure 2. Block diagram of PV generation and ES system based on VSG control.

Figure 3 .
Figure 3.Comparison between small-signal model and Simulink simulation model.

Figure 3 .
Figure 3.Comparison between small-signal model and Simulink simulation model.
positions and velocities of particles after the k-th iteration, and n represents the num-

Figure 10 .
Figure 10.Eigenvalue analysis of the system under global optimal controller parameters.
. Specifically, Ppv, P * pv ; Pbat, P * bat ; Pvsg, and P * vsg represent the pre-and post-optimization output powers of the PV array, ES, and VSG.

Figure 10 .
Figure 10.Eigenvalue analysis of the system under global optimal controller parameters.

Figure 11 .
Figure 11.Small disturbance stability simulation results: (a) power setpoint increased by 1 kW from 2 to 4 s; (b) irradiance increased by 100 W/m 2 from 2 to 4 s; and (c) grid frequency dropped by 0.3 Hz from 2 to 4 s.

Figure 11 .
Figure 11.Small disturbance stability simulation results: (a) power setpoint increased by 1 kW from 2 to 4 s; (b) irradiance increased by 100 W/m 2 from 2 to 4 s; and (c) grid frequency dropped by 0.3 Hz from 2 to 4 s.
Electronics 2024, 13, x FOR PEER REVIEW 3 of 19 frequency disturbance resistance of PV and ES units.Compared to traditional photovoltaic and storage grid-connected systems, a photovoltaic-storage power generation system based on VSG control possesses rotational inertia and damping sharing features, greatly improving the output power and frequency disturbance resistance of PV and ES units.
load Figure 1.Topology of PV generation and ES system.

Table 1 .
Key parameters of PV generation and ES system.

Table 2 .
Eigenvalues of PV generation and ES system.