Small-signal analysis of multiple virtual synchronous machines to enhance frequency stability of grid-connected high renewables

The virtual synchronous machine technology is considered as an important technique to effectively control the shortcomings of renewable energy-based power electronics inter-faces, providing backup inertia and regulating grid stability. Conventionally, the virtual synchronous machine with a large capacity is responsible for controlling the entire grid stability against renewable penetration. It is usually operated as a centralised control system. But what if virtual synchronous machines with small capacities are independently operated by their additional droop control schemes, and will they present better performance than the single virtual synchronous machine? This study proposes the multiple virtual synchronous machine system with different active power-frequency ( P-f ) droop characteristics to improve inertia support regarding frequency stability improvement. The comprehensive small-signal modelling of the multiple virtual synchronous machine unit is designed to include the additional P-f droop characteristics. Then, the dynamic characteristics (steady-state and transient responses) and static stability of the multiple virtual synchronous machines are compared with the single virtual synchronous machine at the same rated capacity in both eigenvalue/sensitivity-domain and time-domain analysis. The obtained results reveal that the system with the presence of several virtual synchronous machines is more stable than the system with single virtual synchronous machine, maintaining stable and secure system operation during the contingency.

inertia, the rate of change of frequency (ROCOF) of the system could be high enough to trigger the load-shedding scheme, even at a small size of a disturbance. If the ROCOF increases beyond the permittable threshold, it can trigger the tripping of protection relays and generators. Severe frequency oscillations with large transients could happen during the disturbance, causing instability, cascaded failures and partial or complete blackouts [2]. As a result, a new challenge has appeared, that is, how to regulate the stable and secure operation of the system with RESs integration. With the growing share of RESs towards 100%, future power systems will be less stable and secure than ever.
To overcome this dilemma, the concept of virtual inertia control named virtual synchronous machine (VISMA) [3][4][5], virtual synchronous generator (VSG) [6][7][8] or synchronverter [9] is proposed to effectively emulate the inertial and primary responses through its droop control. Although the control techniques emulating virtual inertia in the mentioned literature are slightly different, the principle is similar in the aspect that all of them duplicate the inertial characteristic of the synchronous generator by imitating its basic swing equation. The concept of virtual inertia control can be formed by the short-term energy storage system (ESS), inverter and appropriate control method [10,11]. Thus, this concept can effectively emulate the additional inertia and damping supports, enhancing the system transient and steady-state characteristics. The further benefit of the VISMA over the synchronous generators is its capability to absorb power in cases of power excess. The processes of designing VISMA control parameters to meet the objective of maximising power system stability during RESs penetration have been discussed in [8,[12][13][14][15][16][17]. The authors in [13] discussed about the VISMA parameter tuning technique using trial and error. The authors in [8] proposed the VISMA parameter tuning method based on the transient energy analysis. The authors in [14] presented a method to select proper VISMA parameters regarding system frequency protection scheme. Then, the authors applied the robust control theory [17], model predictive control [12] and fuzzy logic [15,16] to determine optimal VISMA parameters regarding system stability enhancement. Consequently, it becomes a promising technique for integrating RESs into the existing power systems.
Focusing on the small-signal modelling and analysis of VISMA, frequency stability is considered as an important phenomenon that results in the system instability. By integrating high levels of RESs into the system, several research works have designed the small-signal models of VISMA to study the low inertia problem in both isolated and grid-connected systems. The small-signal models of VISMA were designed to control the stability of single area power systems with RESs [5,13,18]. Then, the small-signal VISMA models were utilised to improve the stability of interconnected HVDC systems [2,19,20] and interconnected HVAC system [21] with RESs. The small-signal models of the self-adaptive VISMA design were proposed in [22]. The dynamic impacts of the small-signal VISMA model were fully investigated in [13,23]. The effects of frequency measurement-based VISMA were analysed in [24]. However, the aforementioned references studied the small-signal VISMA models as a single (centralised) system with a large ESS capacity for responding to all disturbances and uncertainties.
Nowadays, RESs integration is more decentralised and distributed, creating more energy available with higher power quality for local communities [25]. Together with ESS development, this technology has been improved towards high energy density [26]. It is strongly believed that not only the amount of VISMA power but also the VISMA placement/location can significantly impact system efficiency [27,28]. The authors in [28] proposed a comprehensive analysis in an attempt to address the optimal inertia placement issues, whereas the robust H 2 optimisation technique is used, accounting for the worst-case disturbance location. The authors in [27] presented a techno-economic method to solve an optimal placement problem of VISMA regarding frequency stability point of view. The new framework is designed, which considers ESS features, lifetime, annual costs and state of charge into the optimum placement objective, improving frequency response with a minimum storage capacity. Subsequently, VISMA systems may own a smaller power rating compared with the entire system. Thus, it is expected that VISMA utilisation will be distributed around the system. The emulation of inertia power will require multiple-VISMA units to ensure the stable and secure performance of system operations. Compared with the single VISMA with a large capacity, the concept of multiple VISMAs with smaller sizes is that each VISMA system can independently operate based on its individual/additional control droops, helping each other in emulating inertia and damping supports in response to the whole disturbance of the system. In the near future, the multiple-VISMA control scheme will be implemented for a wide range of grid support. Until now, only a few studies have investigated this idea [29][30][31][32]. The authors in [30] proposed the fuzzy logic controller for regulating multiple VISMAs in an islanded system. The authors in [29] presented the game theory technique for controlling multiple VISMAs. In [31], an optimisation method for controlling multiple VISMAs has been designed. In [32], the authors proposed a stator reactance adjuster to improve the active power damping of parallel VISMAs. Although the mentioned references proposed the studies of multiple VISMAs, their multiple-VISMA systems have not analysed the dynamic effect of virtual control droop characteristics, which is an important factor for controlling the emulated active power set-point after the disturbance regarding frequency deviation. The dynamic behaviour greatly depends on such droop characteristic. Moreover, the mentioned references did not focus on the small-signal stability analysis of the multiple-VISMA system, which is essential for frequency study, control and design. The small-signal stability analysis can effectively identify the cause of power system deviations or oscillations, evaluating the solutions for system frequency stabilisation. To fill in this research gap, the comprehensive small-signal modelling of the multiple VISMAs should be designed and analysed. Further investigation on the dynamic stability of the multiple VISMAs among their individual control droops/parameters must be performed to effectively improve frequency stability and resiliency of the system. This work studies the control design of the multiple-VISMA system for improving inertia support regarding frequency stability. The small-signal model of the multiple-VISMA system is designed corresponding with additional active powerfrequency (P-f) droop characteristics for frequency stability analysis. To validate the efficacy of the multiple-VISMA control scheme, the transient and steady-state stability, including dynamic characteristics between the single VISMA and multiple VISMAs are compared under various system conditions. The main contribution of this research over the existing VISMA approaches are summarised as follows: (1) this is the first paper that designed the small-signal model of the multiple-VISMA unit with different control droop characteristics for frequency stability analysis; (2) transient and steady-state effects of the multiple-VISMA control scheme are validated by the timedomain and eigenvalues-domain analysis, determining a new trade-off between frequency control and virtual inertia emulation; (3) compared to the conventional methods (i.e. centralised VISMA) in [13,[18][19][20], the proposed multiple-VISMA technique provides superior frequency performance and stability to the system, which is expected to penetrate more RESs, avoiding instability and power blackouts.
The rest of this work is organised as follows. Section 2 describes a brief overview of frequency stability regarding fundamental inertia control based on the well-known swing equation. The small-signal model of the studied power system with RESs is also presented. Section 3 proposes the small-signal modelling of the multiple-VISMA system considering additional control droop characteristics. The practical overview of the multiple-VISMA system is also displayed. Section 4 reveals the simulation results of the multiple-VISMA approach under contrastive case studies. Lastly, it ends up with the conclusions in Section 5.

Inertia response for frequency stability control
Following a disturbance, the frequency response of the system is divided into three stages: (1) inertia (inertial) response, (2) primary response and (3) secondary response [33,34]. At the inertia response stage, the rotating kinetic energy accumulated in synchronous machines or generators is immediately released to regulate the mismatched power between generation and demand, effectively reducing the ROCOF of the system. If the ROCOF is higher than the permittable limit, it can activate the load shedding or trip the protection relays and generators. Thus, the inertia response is an indicator of the rotating kinetic energy that contributed to the system.
The inertia response regarding mechanical/electrical torque of a synchronous machine i can be expressed based on the swing equation relationship as where J s,i is the moment of inertia (kg/m 2 ) of the synchronous generator; ω i is the angular velocity of the synchronous rotor (rad/s); i = 2 f 0 , f 0 is the nominal frequency (Hz); T m,i, and T e,i are the mechanical and electrical torque of a synchronous generator; P m,i and P e,i are the mechanical and electrical power of a synchronous generator. The rotational kinetic energy is generally normalised to inertia constant (W•s) that represents the ratio of the accumulated kinetic energy at the rated rotating speed to the rated capacity of the synchronous machine. In real practices, the inertia constant indicates the time interval to provide energy for the demand, which equals to the rated capacity of the generator without any additional mechanical input. The dynamic equation or setting of the inertia constant (H) can be calculated as [1,33,34] where S i is the rated capacity of generator i, E i is the accumulated kinetic energy of generator i and H i is the inertia constant of generator i. For the multi-generator system, the inertia contributors from other synchronous generators can be considered. Thus, the equivalent inertia constant of the entire system is expressed as [33,34] where N is the number of the synchronous machine, H sys is the equivalent inertia constant of the system and S sys is the rated capacity of the system. The frequency stability problem is characterised by the power unbalance between generation and demand. The dynamic behaviours between active power and frequency of a synchronous machine following the disturbance can be effectively modelled by the swing equation as [33,34] where P m,i is the mechanical power of generator i, P e,i is the electrical power of generator i, which is characterised by the electrical load demand, f i is the rotor frequency of generator i and D i is the damping coefficient of generator i. It is noted that at the initial operating point, the mechanical input torque (or power) of a synchronous generator is assumed to be constant and equal to the electrical load power. As a result, the deviation of mechanical power input (ΔP m ) is zero. After the electrical power disturbance (ΔP e ) occurs, the P e will be adjusted by a turbine-governor system in order to minimise the frequency deviation (Δf). Following the disturbance, the ROCOF of the generator is arrested by the inertia constant and damping coefficient in the short time interval. Hence, the generator frequency cannot change immediately, enhancing frequency stability.
With an approximation, the dynamic behaviours of a control area can be defined as an equivalent synchronous generator, resulting in the aggregated swing equation model as [33,34] where P m,j is the total mechanical power of control area j, P e,j is the total electrical power of control area j, which is characterised by the total electrical load demand, f j is the aggregated frequency of control area j, H j is the equivalent inertia constant of control area j and D i is the damping coefficient of control area j.
During normal operation, all parameters in (5) can change around the steady-state operation point. Subsequently, the control law in (5) can be represented in terms of the small-signal model with the incremental formulation as By assuming Δ P m, j = 0, the Laplace transformation is then applied to both sides of (6). Consequently, the first-order transfer function characterised by H j and D j of the control area j can be derived as where ΔP e, j is the total electrical power deviation/change in control area j, and Δ f j is the frequency deviation of control area j.

Small-signal modelling of RESs-dominated power system for frequency analysis
The small-signal stability indicates the system potential to maintain synchronism when experienced with a small disturbance. When the disturbance is considered to be small, the dynamic equations that explain the system response can be linearised for the analysis objective. Due to the recent replacement of synchronous generators by RESs in the grid, the system frequency oscillations due to the absence of sufficient synchronising and damping torques can greatly occur, leading to instability and power failures. In real practices, the small-signal stability problem is related to the insufficient damping of the system oscillations [33].
In this study, the small-scale power system with RESs (i.e. microgrid) is constructed as depicted in Figure 1. It consists of the 12 MW of synchronous generator, 6 MW of photovoltaics, 7 MW of wind-turbine system, 10 MW of industrial load, 5 MW of residential load and VISMA. The study of VISMA systems is proposed by two control strategies: (1) single VISMA with the rated capacity of 6 MWh and (2) two-parallel VISMAs with each rated capacity of 3 MWh (total rated capacity is 6 MWh). The system base is 15 MW. The power plant-based synchronous generator offers inertia, primary and secondary frequency control services for the system [2,19]. To create more drastic operation, the photovoltaic and wind systems (RESs) have not controlled to penetrate their power into the system. The competence of such generations relies on solar irradiation and wind velocity. Subsequently, the RESs could not contribute to the frequency regulation of the system. In fact, the RES and load power generations are considered as the small-signal disturbances to the system. Due to RESs penetration, the system is lack of inertia and damping properties; thus, the VISMA system is expected to imitate the required inertia and damping, maintaining stable and secure system operation. To clearly reveal the effectiveness of the multiple VIS-MAs, the studied system has not connected to any large utility, performing the isolated operation.
To perform the small-signal stability analysis-based frequency control, the dynamic model of the studied system is presented in Figure 2. Technically, the turbine system receives the feedback mechanism from the governor unit and makes the generation to track the RES/load change (disturbance). The important physical limitation of generation rate constraint (GRC) is considered in this model, generating more actual dynamics. The GRC is represented by the limiter block, where V U and V L represent the upper and lower constant that limit the rate of turbine gate opening/closing speed. The small-signal (dynamic) equation of the turbine system can be derived as where ΔP m is the change in active power from the turbine unit, ΔP g is the change in speed-governor control action (signal) and T t is the time constant of the turbine unit. The speed-governor unit measures the variation in speed (frequency) through the primary and secondary control loops. The where ΔACE is the change in secondary control action (signal), Δ f is the change in system frequency, R is the droop constant of primary control and T g is the time constant of the speed governor. The secondary control unit (i.e. area control error: ACE) measures the feedback through frequency deviation and combines the signal to the primary control-based governor unit via the dynamic integral (K i ) controller. The effect of time delay (e −sT ) is added to perform the dynamic behaviour of filters and communication delays related to the secondary control process. In real practices, the dynamic controller is usually used as an integral controller or integral-proportional controller due to its simplicity. The small-signal (dynamic) equation of ACE for secondary control can be derived as where K i is the integral control constant, T is the time delay constant, and β is the bias factor. The wind-turbine system has not participated in frequency control. The ability in generating power depends on the wind speed for a specific time. The small-signal (dynamic) equation of the wind system can be derived as where ΔP W is the change in generated power from the windturbine, and ΔV wind is the initial wind power-based wind speed change.
Similarly, the photovoltaic (PV) system does not participate in frequency control. The ability in generating power depends on solar irradiation for a specific time. The small-signal (dynamic) equation of the photovoltaic system can be derived as where ΔP PV is the change in generated power from the PV system, and ΔG solar is the initial PV power-based solar irradiation change. According to Figure 2, the frequency experiences the dynamic change (Δ f ) after disturbances in the load power change (ΔP L ) and RES power change (ΔP W , ΔP PV ). Hence, the small-signal equation of frequency deviation that expresses the response of the studied system can be defined as where H is the system inertia constant, D is the damping coefficient of the system and ΔP VISMA is the change in net emulated (virtual) inertia power from VISMA system.
It is noted the small-signal (dynamic) model of the VISMA will be described in the following section. Based on [2, 17, 18-21, 24, 33, 35-38], these research works confirmed that the dynamic model Figure 2 is sufficient enough for frequency stability study.

Fundamental concept of VISMA
To resolve the low system inertia problem driven by the inverter-based RESs penetration, a virtual inertia control mechanism, such as VISMA, VSG or synchronverter, is proposed to emulate the inertia support based on the dynamic characteristics of synchronous machines, enhancing system stability and resiliency [2,4,17]. The virtual inertia system could be built by a short-term ESS, inverter and suitable inertia control mechanism. As described [4,11,13,39], the ESS unit is equivalent to a rotor, and the inverter unit is equivalent to a synchronous machine in imitating the required inertia power and damping property. Hence, the VISMA system can be operated as an actual synchronous machine in generating the additional inertia and damping for the short time intervals. Subsequently, this concept can offer a fundamental for controlling a high portion of RESs in the existing power system. Previous accomplishments with literature reviews on virtual inertia applications can be found in [1,10,11,40,41].

Small-signal modelling of VISMA for frequency analysis
To perform the small-signal stability analysis regarding frequency control, the dynamic model of the VISMA, which consisted of the virtual rotor control action and inverter-based ESS, is proposed in Figure 3. This dynamic model is constructed by considering the inertia emulation technology-based derivative technique [2,13,19,20] to determine the ROCOF of the system for modifying the additional active power with the required inertia and damping to the set point of the system during the disturbance. The blocks of virtual inertia (J vi ) and damping (D vi ) are used to compensate for the absence of system inertia and system damping driven by RESs integration. The low-pass filter block is applied to eradicate the noise problem and to perform the actual behaviour of fast response of the inverter-based ESS. The power limiter block is used to limit the ESS power output with respect to the capacity constraint. During normal operation, the emulated power variable (ΔP VIS ) in (14) can change around the steady-state operation point, tracking the system frequency deviation (Δ f ) caused by the disturbance. The incre-mental control equation of VISMA can be defined as [2] ΔP VISi (s) = where J vi is the virtual inertia constant of VISMA i, D vi is the damping coefficient of VISMA i, T ESSi is the time constant of the inverter-based ESS i and i is the number of VISMA in the control area. Nowadays, the ESS technology is developed towards high energy density. The VISMA unit can own a small power rating compared with the whole system due to economic constraints. Thus, the inertia power emulation will require multiple-VISMA units to ensure stable performance on system operation. Subsequently, in the future, multiple VISMAs with different control droops will be expected to widely integrate into the RESs-dominated power system. Accordingly, it is essential to further study the multiple VISMAs among their additional control droops/parameters, improving inertia support and damping supports regarding the frequency stability of the system.

P-f droop characteristics for multiple VISMAs
Previously, the RESs-dominated power system usually consists of a large VISMA unit with single P-f droop characteristics for emulating system inertia regulating frequency stability control. Considering the economic constraints and technology development from the ESS, this study proposes the multiple-VISMA concept to augment frequency stability of the system. The multiple-VISMA control technique provides a new structure of the system in which several inverter-based ESS units are independently controlled to help each other in providing inertia and damping supports. For understanding the VISMA droop characteristic concept, the dynamic equation of the VISMA control droop is proposed to establish a feedback loop responding to the control of emulated inertia power and frequency as where R vi is the virtual droop constant of VISMA i and P VISi_0 is the nominal value (reference) of emulated inertia power from VISMA i. As seen in Figure 4, two VISMA control schemes (i.e. single VISMA and multiple VISMAs) with different P-f droop characteristics are described. These VISMA units are operating at a unique nominal frequency with different emulated inertia powers. In this study, the parallel VISMAs with additional Pf droop characteristics are proposed instead of using a single VISMA with a single P-f droop at the same rated capacity. After the disturbance (e.g. RES/load power changes), each VISMA unit increases or decreases its emulated power until reaching a new common operating frequency. As shown in (15) and Figure 4, the amount of emulated power from each VISMA to compensate the disturbance depends on the VISMA's droop  characteristic. Then, the dynamic relationship in (14) can be modified to include the dynamic effect of the P-f droop control as Thus, the ratio of the frequency deviation to change in the emulated inertia power output of VISMA is known as the virtual droop and can be expressed as Lastly, the small-signal model of multiple VISMAs is proposed in Figure 5, which is modified to integrate the dynamic effect of different P-f droop control characteristics, providing additional inertia and damping supports to the system. The dynamic control equation of the multiple-VISMA system is defined as

SIMULATION RESULTS AND DISCUSSION
In this section, the small-signal stability assessment is focused on three essential control investigations; 1) the evaluation of VISMA control variables; 2) the effectiveness of multiple VIS-MAs under sudden RES and load changes; 3) the robustness of the multiple VISMAs under the critical lack of system inertia and damping triggered by high RES integration. The simulations are conducted using the MATLAB/Simulink® and examined in both eigenvalue/sensitivity domain and time domain. Two different control configurations of VISMA have been constructed; 1) single VISMA unit with a single P-f droop based on [13,[18][19][20]; and 2) multiple-VISMA unit with additional P-f droops. The single VISMA system has a rated capacity of 6 MWh. The proposed multiple-VISMA system consists of two small scales of VISMAs with the individual rated capacity of 3 MWh (i.e., the total rated capacity of 6 MWh), which are operated together as parallel P-f droops to uniquely modify the total emulated power set-point and emulate additional inertia/damping supports, resulting in better system response.

Small-Signal Analysis of VISMA Control Variables
To examine the small-signal stability of the system equipped with two VISMAs, trajectories of dominant eigenvalues are derived by changing the virtual droop (R vi ), virtual inertia (J vi ) and virtual damping (D vi ). The arrow specifies the increasing/decreasing of corresponding control parameters. Then, the result of small-signal stability evaluation of the multiple VIS-MAs regarding the eigenvalue analysis is validated in the presence of a step load/RES change of ±0.08 p.u. as a disturbance under the normal condition of high system inertia and damping (100%) [34]. The control parameters for this scenario are presented in Table A.1. It is noted that the ranges of R vi , J vi and D vi in this study are collected from [13,[17][18][19][20]. Figure 6 shows the influence of the decreasing virtual droops (R vi ) of parallel VISMAs at the same value against the system frequency response. To clearly observe the actual behaviour of R vi , other control VISMA variables are set as small values (J vi = 0.6, D vi = 0.3). Considering the eigenvalue analysis with the decreasing of R vi , the dominant imaginary (conjugate) eigenvalues of the system move towards the left side of the s-plane, indicating the enhancement of small-signal stability (see Figure 6(a) and (b)). Obviously, the dominant conjugate eigenvalues are distinctly further from the imaginary axis, indicating better damping ratios in the dominant oscillation modes of the system. This phenomenon is validated by applying the step disturbance of 0.08 p.u. against the reduction of R vi (see Figure 6(c) and (d)). As R vi decreases, the frequency response of the system reaches more stable performance with a gradual decrease in the frequency drop/overshoot. However, if the R vi is too small, the system will be very sensitive to the disturbance. The VISMA will greatly respond to the disturbance even with a small deviation, whereas the VISMA power output could extensively change, vulnerating the stability margin of the system. Figure 7 displays the influence of the increasing virtual inertia constants (J vi ) of parallel VISMAs at the same value against system frequency response. To clearly observe the actual behaviour of J vi , other control VISMA variables are fixed at low values (R vi = 0.8, D vi = 0.3). Considering the eigenvalue analysis with the increasing of J vi , the dominant conjugate eigenvalues travel to the left side of the s-plane, indicating superior damping ratios (see Figure 7(a) and (b)). However, increasing J vi too high will lead the real eigenvalues to move towards the right side of the s-plane, resulting in the reduction of the damping. Hence, while a high value of J vi is beneficial in terms of frequency stability, it should not be set too high since it would then reduce the damping of the system and might lead to system instability. A higher J vi value will lead to a smaller frequency deviation, but it will also make the system requires a longer time to settle subject to the disturbances. Then, the result of the eigenvalue analysis is validated by applying the step disturbance against the reduction of J vi (see Figure 7(c) and (d)). With the increasing J vi , the system performance becomes small-signal stable, and the frequency nadir/overshoot of the system is greatly arrested, resulting in a vast reduction in frequency drop/dip. Apparently, by applying the larger J vi (obviously over 2.4), the system response   Figure 8(a) and (b)). This phenomenon is validated by applying the step disturbance of 0.08 p.u. against the increase of D vi (see Figure 8(c) and (d)). It is obvious that increasing D vi could enhance the stabilising performance (lower settling time) with a slight reduction in frequency nadir/overshoot (see Figure 8(c) and (d)). Clearly, the stabilising effect/time is alleviated by increasing the suitable D vi . However, if a larger value of D vi (obviously over 2.4) is implemented, the dominant imaginary (conjugate) eigenvalues will move far away from the real axis towards the right side (i.e. higher oscillation), indicating lower damping ratios. This side effect can trigger the large second overshoot to the system as validated in Figure 8(c) and (d). Thus, when the D vi is applied in VISMA, it is better to keep a small value of D vi (less than 2) to suitably maintain the stabilising effect with less oscillation and give the main respon-sibility of arresting frequency nadir/overshoot to the R vi and J vi units.
Finally, following the above analysis, the R vi , J vi and D vi of multiple VISMAs are designed taking initial loadings of the RES/load and maximising the system stability into account. To mitigate the overloading of all RESs/loads, the initial loadings of RESs and loads have been set at +0.08 and −0.08 p.u. (respectively) of the rated power. Under the initial loadings of RES/load, the analysis results of frequency deviations in Section 4.1 by varying the R vi , J vi and D vi are carried out. Obviously, the smaller of R vi , the lower the frequency nadir/overshoot can be achieved. On the contrary, the higher of J vi and D vi , the lower the frequency nadir/overshoot can be achieved, while the longer stabilising time of the system can be significantly observed. As a result, the R vi , J vi and D vi of multiple VISMAs in this study are appropriately and equally set at 2.70, 1.60 and 1.30, respectively.

Small-signal analysis under sudden RES and load changes
This section investigates the effectiveness of the multiple-VISMA system against various step RES and load variations, as Renewable energy source (RES) and load power disturbance depicted in Figure 9. With the introduction of small RES penetration, the system inertia (H) and damping (D) are slightly disturbed by reducing 30% from its nominal values. The control parameters are given in Table A1. The VISMA control parameters are selected based on the small-signal analysis in Section 4.1. Figure 10 shows the system response-based time-domain analysis under sudden RES/load changes. Following the 30% reductions of system inertia and damping, the system is more sensitive to the disturbances. When the RES power increment of 0.2 p.u. is immediately integrated at 20 s, it greatly increases the power mismatch of the system and causes the sudden frequency overshoot. Consequently, the system without VISMA becomes unstable and could not maintain frequency stability, leading to the higher oscillation, system instability and collapse. As seen in Figure 10(a), the systems with single VISMA and multiple VISMAs could properly maintain frequency stability after the applied disturbance, preventing instability and system collapse. By applying the multiple VISMAs, the transient stability is significantly improved by the deployment of parallel P-f droop characteristics, which successfully provides the additional damping torque regarding inertia and damping properties. The frequency nadir and overshoot are effectively arrested, resulting in the lowest ROCOF (see Figure 10(b)). This advantage leads to the lower frequency drop and dip during the contingency. Thus, after using the multiple-VISMA control scheme, the significant reduction in frequency drop (nadir), overshoot and ROCOF can be achieved, leading to more stable and secure stability and resiliency of the system. Figure 10(c) displays the power responses of synchronous generator and VISMA units, respectively. It is obvious that the multiple VISMAs with additional droops could extract more emulated inertia power than the single VISMA with one droop. Subsequently, the synchronous generator is less stressed and demands lower power generation, improving the stability margin of the system. Specifically, with the additional derivative terms to imitate inertial response, the power response of the multiple VISMAs will be faster than the single VISMA, which  contributes to the smoother change of system frequency. This also confirms that the multiple-VISMA system has a better dynamic response than the single-VISMA system under the same (total) rated capacity of inverter-based ESS.
The results for small-signal stability are presented in Figure 11 and Table 1. Compared with the case of no VISMA, after applying the single VISMA, the eigenvalue trajectories travel to the left side of the s-plane, and the damping ratio increases from 0.2891 to 0.7071, indicating small-signal stable. By applying the multiple VISMAs, the eigenvalues further move towards to the real axis, indicating the improvement of small-signal stability. The damping ratio finally reaches 0.9624, which is more than 25% that achieved in the case of single VISMA. This is confirmed that the multiple-VISMA system has a wider stable region than the single-VISMA system against sudden RES/load changes, including system inertia and damping reductions.

Small-signal analysis under low system inertia and damping triggered by high RESs penetration
In this section, the robustness of the multiple-VISMA system is examined under the critical operating condition of various changes in primary control, secondary control, system inertia and damping. From the frequency stability point of view, the critical operating condition is when the system operates under high RESs penetration (see Figure 12(a)) with low load  System response under low system inertia and damping triggered by high renewable energy sources (RESs) integration consumption (see Figure 12(b)), resulting in a huge mismatch of power unbalance. In this condition, the system inertia and damping are greatly decreased to 20% of its nominal values, performing actual behaviour of the RES/inverter-dominated system. In real practices, the operating condition of the system also changes with time, which may cause the uncertainty, degrading system performance and stability. Thus, the parameter variations of primary control (T g ≈ +55% and T t ≈ +25%) and sec-ondary control (K i ≈ −50%), including communication delay (T ≈ +120%) are applied as shown in Table A.1 to perform more drastic operating condition. The VISMA control parameters are selected based on the small-signal analysis in Section 4.1. Figure 13 exhibits the frequency response of the system under the severe operating condition. In this condition, the frequency response widely oscillates with larger transients, indicating the degradation of system stability and performance. In

FIGURE 14
Virtual synchronous machine (VISMA) power deviation under low system inertia and damping triggered by high renewable energy sources (RESs) integration the case of no VISMA, the system is totally unstable, leading to instability and power blackout. Obviously, during the wind integration (at 400 s), the system with the single VISMA yields the largest frequency overshoot/transient of 0.38 Hz, which may trigger the protection schemes or generation tripping. By applying the multiple-VISMA system, the frequency overshoot reduces to 0.17 Hz, which is more than 50% reduction in frequency excursion. Cleary, the frequency oscillation is smoother with lower frequency nadir and overshoot. It is obvious that the frequency drop and dip are reduced by a half when compared with the single VISMA at the same (total) rated capacity of inverter-based ESS. This is because of the contribution from the multiple P-f droops, which could virtually create additional inertia and damping components, improving small-signal stability. Therefore, the higher frequency drop/overshoot and ROCOF can be alleviated by the additional P-f droops in the multiple-VISMA system. The power response of each VISMA unit under this condition is also shown in Figure 14. Following the severe operating condition, it is obvious that the multiple-VISMA system still maintains a better dynamic performance in emulating inertia property than the single-VISMA system under the same rated capacity of inverter-based ESS.
The small-signal stability of this scenario is illustrated in Figure 15. Following the severe operating condition, the system FIGURE 15 Eigenvalue trajectory under critical system inertia and damping reduction triggered by high renewable energy sources (RESs) integration without VISMA becomes small-signal unstable as there exists a pair of conjugate eigenvalues (0.0580 ± j5.0252) on the right side of the s-plane, whereas the systems with the single VISMA and multiple VISMAs are small-signal stable against a wide range of RES/load changes and parameter variations. Clearly, the dominant conjugate eigenvalues in the multiple VISMAs have lower natural oscillation frequency and higher damping ratio compared with the single VISMA, which indicates the multiple VISMAs have better dynamic stability and performance than the single VISMA under the same (total) rated capacity of inverter-based ESS.
In summary, the small-signal stability of the system is examined against the reduction of system inertia and damping driven by RESs integration as shown in Table 2. In the case of no VISMA, by reducing the inertia and damping, the damping ratio decreases from 0.36888 to −0.01154, leading to smallsignal unstable. On the contrary, by applying the single VISMA and multiple VISMAs, the damping ratios increase against the reduction of system inertia and damping, indicating the stable and robust system operation. Compared with the single VISMA, the multiple VISMAs could improve the damping ratios in each case around 20%, which makes the system become less sensitive to the lack of inertia and damping triggered by RESs integration. Thus, it is confirmed that the multiple-VISMA system is more robust against parameter changes and RES/load uncertainties than the single-VISMA system.

CONCLUSION
This work presents the small-signal modelling of the multiple-VISMA system with additional P-f control droops to effectively provide inertia and damping supports regarding frequency stability enhancement. 1. Under the same control variables of multiple VISMAs, by decreasing the R vi , the frequency response of the system reaches more stable performance with a gradual decrease in the frequency drop/overshoot. However, if the R vi is too small, the system will be very sensitive to the disturbance, vulnerating stability margin of the system. When the J vi is increased, the system performance becomes smallsignal stable, and the frequency nadir of the system is greatly arrested, leading to a huge reduction in frequency drop/overshoot. However, the increasing J vi introduces the longer stabilising time due to excessive inertia property, which makes the system slower in responding the disturbance. This effect is prone to cause instability, deteriorating the stabilising/settling performance of the system. By increasing D vi , it could enhance the stabilising performance (lower settling time) with a slight reduction in frequency drop/overshoot. However, if a large value of D vi is used, it can trigger the large second overshoot, causing higher oscillation and degrading system stability. Thus, it is recommended to use a small value of D vi to suitably maintain the stabilising effect with less oscillation and give the main responsibility of arresting frequency drop/overshoot to R vi and J vi units. 2. The results of eigenvalue/sensitivity analysis reveal that with the deployment of the multiple-VISMA system, the dominant conjugate eigenvalues move to the left side of the s-plant towards the real axis, indicating the higher damping ratio and small-signal stable. As a result, the multiple-VISMA system has a wider stable region than the single-VISMA system against various RES/load changes, system inertia and damping reductions, including parameter uncertainties. 3. The results of time-domain analysis show that with the deployment of additional P-f droops and suitable inertia and damping control variables, the multiple-VISMA system provides better stability and performance than the single-VISMA system. The significant enhancement in frequency nadir, overshoot and ROCOF reductions can be achieved, resulting in more robust system operation under a wide range of disturbances and uncertainties.