Characteristic Analysis and Indexing of Multimachine Transient Stabilization Using Virtual Synchronous Generator Control

: As distributed power sources via grid-connected inverters equipped with functions to support system stabilization are being rapidly introduced, individual systems are becoming more complex, making the quantiﬁcation and evaluation of the stabilizing functions difﬁcult. Therefore, to introduce distributed power sources and achieve stable system operation, a system should be reduced to a necessary but sufﬁcient size in order to enable the quantiﬁcation of its behavior supported by transient theory. In this study, a system in which multiple distributed power supplies equipped with virtual synchronous generator control are connected is contracted to a two-machine system: a main power supply and all other power supplies. The mechanical torque of each power supply is mathematically decomposed into inertia, damping, synchronization torques, and the governor effect. The system frequency deviations determined by these elements are quantitatively indexed using MATLAB/Simulink. The quantiﬁcation index displayed in three-dimensioned graphs illustrates the relationships between the various equipment constants of the main power supply, the control variables of the grid-connected inverter control, and the transient time series. Moreover, a stability analysis is performed in both the time and frequency domains.


Introduction
Conventional power systems use synchronous generators as their main power source, the swing characteristics of which suppress the power system swing owing to the load power and generated power fluctuations [1]. However, due to recent advancements in distributed power sources, such as renewable energy, a large number of grid-connected inverters are being connected to power systems. The inertial power of the entire system, which generally suppresses the fluctuations of the power system, is insufficient when the capacity of the grid-connected inverters exceeds the capacity of the main power supply, making it challenging to maintain stability [2,3]. There is thus a need for the grid-connected inverters of distributed power sources to have power generation functionality and to suppress the fluctuations of the power system [4]. Obviously, it is necessary to address other existing technical issues such as the intermittentness of solar energy [5], limitations on the data capacity of smart grids [6], the influence of power electronics equipment on power quality [7], and the effectiveness of energy trading; various studies have investigated these issues [8]. Although the problem of the capacity of grid-connected inverters exceeding that of the main power supply was previously regarded as being characteristic of microgrids (MGs), a similar problem has recently been identified in large-scale smart grids. Field research using actual data is therefore being actively conducted (in Korea [9], Canada [10], Oahu, HI [11][12][13], and the USA in general [14]). To leverage open-source system validation platforms (SVPs) for interoperable test procedures, research using power/controller hardware-in-the-loop (PHIL/CHIL) is also being actively promoted [15][16][17][18]; Ref. [19] summarizes these trends.
So-called "smart inverters" for photovoltaic power generation in detached houses have a functionality that can contribute to maintaining the balance between electricity supply and demand at steady state [20,21]. Some smart inverters realize grid-support functions (GSFs) based on power control according to the communication of commands at a higher level of integrated control, in which case it is relatively easy to estimate the effects and formulate the required specifications. In contrast, some inverters autonomously calculate the power imbalance of the system without higher-level control commands, in a manner similar to that of synchronous generators, immediately after fluctuations in the load/generated power. In this way, high-performance grid-connected inverters can control their own inputs and outputs using high responsiveness, which has been attracting increased attention. The autonomous functions are generally installed in the primary control unit of the inverter, and the details are rarely disclosed by manufacturers. It is therefore difficult to estimate their effects and formulate the required specifications [22]. A protocol study was conducted to comprehensively test the capability of smart inverters to stabilize MGs [23,24]. Reference [24] examines various testing methods using three commercially available, single-phase, residential-scale photovoltaic (PV) inverters from three different manufacturers. In [23], a test protocol for smart inverters that can utilize battery energy storage systems is studied, and a summary is provided of issues related to four interoperability function tests defined in IEC TR 61850-90-7: the specified active power test (INV4), the var-priority volt-var test (VV), the specified power factor test (INV3), and frequency-watt control (FW) from storage.
Research on optimizing parameter settings to improve the performance of distribution systems by making smart inverters compatible with volt-var and volt-watt functions is in progress [9,14,[25][26][27][28][29]. It has been proposed that if the volt-watt function is prioritized over the volt-var function, the system voltage will vibrate continuously [27]. Conversely, when suppressing the volt-watt function, fairness is required to reduce the amount of power generated by distributed power sources (PVs, etc.) [29]. To establish a method for handling these problems, the simulation analysis of an actual field is required: an analysis using actual data in an actual feeder [9,14], quasi-static time series simulation using an open distribution system simulator (OpenDSS) [9], and HIL simulation [29] have been reported thus far. Moreover, in [26], the effects of improved volt-var control using the particle swarm optimization (PSO) algorithm on line loss and operation cost reduction are demonstrated, and in [25,28], the optimum utilization of PV inverter capacity is described.
Some smart inverters obtain the same characteristics by simultaneously mimicking the swing and excitation characteristics of a synchronous generator [30][31][32][33], while others obtain similar characteristics without focusing on the excitation characteristics [34,35]. Unfortunately, there is little research on volt-var, volt-watt, and frequency-watt control that focuses on transient responses related to the swing characteristics of synchronous generators. Although the GSFs of frequency-watt controls between grid-forming and gridfollowing PV inverters are discussed in [12], the transient response, which corresponds to the inertial response, has not been discussed. However, research on improving transient responses with inverters is ongoing [35][36][37][38][39][40]. References [38][39][40] propose an improved GSF by using the direct current (DC) energy buffer of a grid-connected inverter. In [40], a new distributed consensus-based control that reduces the circulating power between parallel running inverters is demonstrated. In [39], experimental results that achieve an inertial response by connecting only a small capacitor to the DC link part are presented. Reference [37] presents the experimental results for a decentralized control that controls multiple PV inverters as master/slave without using a storage battery. In [38], a control that achieves both a droop response and an inertial response using a hybrid energy storage system (ESS) and that reduces the deterioration of the battery storage by using a super capacitor is proposed. A method that strictly simulates the swing characteristics of a synchronous generator using battery power is called virtual synchronous generator (VSG) control, and many studies have been conducted to improve this method [30,31]. For example, reference [35] presents the prediction model of a VSG to enhance suppression performance of the voltage and frequency variations. In [36], a method for adding a phase adjustment function to multi-loop VSG control was introduced, and it was reported that renewable energy-based inverters were used in a plug-and-play manner to reduce the harmonics and oscillation and to improve the dynamic response and power quality.
The swing characteristics of the synchronous generator to be mimicked in this study are determined to be the characteristics of the rotor, expressed by both the inertial and damping responses. The response attributed to the synchronization force is added when multiple synchronous machines are interconnected [41]. After some time, the control of the governor becomes effective, and a governor-free response is added to the inertial and damping responses [42]. The deviation of the system frequency represents the difference between the electrical angular velocity and rated angular velocity of the rotor reference, which is equal to the deviation of the internal phase difference angle (angular velocity).
That is, to accurately analyze the transient fluctuation of the system frequency, it is essential to evaluate the inertial effect, damping effect, synchronization effect, governor effect, etc. separately from the swing equation of the generator rotor angle (mechanical angle). Nevertheless, few studies have analyzed the system frequency by separating and analyzing the components of the swing equation of the generator rotor in detail. Therefore, in this study, the swing characteristics of the generator mimicked by the grid-connected inverter are accurately separated by effect, and the correlation between each effect and the values of the control parameters is clarified. The mechanical torque of each power supply, which is calculated mathematically, is decomposed into inertia, damping, synchronization torque, and governor effect, and the decomposition components were analyzed and indexed from both time and frequency domains. This provides an indicator of the parameter settings required for grid-connected inverters not only in MGs but also large smart grids. To acquire the time-series data that is used to decompose the above-mentioned mechanical torque, the Institution of Electrical and Electronic Engineers (IEEE) 9 bus model is adopted, and two distributed power sources are connected to the existing power system through a grid-connected inverter. Although simulation results are to be verified separately, it should be noted that it is not theoretically necessary to be limited to IEEE 9 bus from the viewpoint of reducing multiple machines to a two-machine model. The transient characteristics are analyzed using numerical simulation results from MATLAB/Simulink according to a mathematical model, and the correlation between each control parameter and each response is clarified based on the concept of the center of inertial frequency. In this system, a synchronous generator is the main power source. The two inverters are reduced to one, and each is equipped with the same VSG control. The total capacities of the grid-connected inverters are assumed to be almost the same as that of the main power supply, thereby making it difficult to maintain system stability. That is, the multimachine large-scale system was reduced to a two-machine model. This not only reduces the mathematical analysis associated with a multi-machine model but also allows analysis that can be focused on various parameter ratios between the conventional main power supply and other distributed power sources.
The structure of this study is as follows: Section 2 reviews the interconnection principle of multiple synchronous machines, and Section 3 analyzes the transient response characteristics from simulation test results in response to changes in each parameter of the VSG control inverter. The frequency of the center of inertia is also calculated in Section 3, while Section 4 discusses indicators for the VSG inverter based on the outcomes of the torque response and the center of inertia in both the time and frequency domains. Section 5 summarizes this study. Figure 1 is based on the swing equation for a synchronous generator (SG) (Equation (1)) [41]:

Synchronous Generator Model
where s is the Laplacian differential operator, ∆P  In the case of a conventional SG, the governor time constant T [s] is several hundreds of milliseconds, and its effect appears after the delay from the onset of the fluctuation in the input/output power. The synchronization torque term Kδ is relevant when the SG is connected to another synchronous machine, and the synchronization coefficient K [pu/rad] is determined by the interconnection style (bus model). Therefore, assuming that the governor response is not on time, the transfer function G sg (s) in a stand-alone SG immediately after the onset of input/output power fluctuations is only determined by the inertial torque term (M/ω 0 )s 2 δ and the dumping torque term (D/ω 0 )sδ.
The synchronous machine block, a device in the standard library of Simscape Electrical in the MATLAB/Simulink environment, is used for the governor model H sg (s), for which Equation (1) holds. A model of a general excitation system is also included, but details are omitted here because the swing equation is the focus of this study.

Virtual Synchronous Generator Control Model
. T Q [s] is the filter time constant, and the conventional proportional-integral (PI) compensator is used as the AVR.
The vector diagram enclosed in the red frame in Figure 2 is the impedance model, which illustrates the relationship between the generator voltage vector E f , terminal voltage vector V g , command value vector I * of the inverter current, and impedance r, x [pu]. The impedance r, x [pu] is a pseudo-impedance set in the VSG control software, and it becomes dominant over the alternating current (AC) output filter impedance in low-frequency regions, such as commercial frequencies. The ratio of the virtual inductance x and virtual resistance r is empirically set to x = 2r. The current command vector value I * is expressed by Equation (2) [30,31], as a vector of the dq coordinates:

System Model
In this study, N synchronous machines are connected in a system: one primary power supply is assumed to be an SG, and (N = 1) generators are assumed to be distributed power sources. The IEEE 9 bus model (N = 3 machine systems) is adopted, and two inverters equipped with the same VSG control are connected. The two VSG inverters are contracted, and a test in which the total capacity, inertia, reactance, and other variables are converted based on the combination of the two inverters being approximately equal to the SG is used as the standard test. In other tests, the VSG control parameter was increased or decreased for comparative purposes. In this study, we do not focus on the DC side of the inverter: we assume a constant DC power supply, and the capacity of DC side is not limited. Figure 3 shows a connection diagram for the system. The base loads are evenly connected to the three buses A, B, and C between pairs of power supplies. For the purpose of running the simulation smoothly, the total base-load capacity is supplied from the SG, and an additional load is stepped onto bus B between the two VSG inverters after the power supply for the base loads has settled. The frequency deviation for the base load is compensated by the SG torque command value, and the phase difference angle after settling is subtracted as an initial offset value. The voltage deviation is considered to be small, and the terminal voltage of each power supply where the additional load is connected is assumed to be equal to the standard voltage at V 0 ∠0.  Figure 4 presents a power-frequency control block diagram of two synchronous machines [31]. G sg (s) and G vsg (s) indicate the same responses as those in Figures 1 and 2, and the governor models H sg (s) and H vsg (s) are considered along with the inertial and dumping characteristics. Herein, the subscripts "sg" and "vsg" represent the SG and the contracted devices of VSG2 and VSG3 in Figure 3, respectively, and subscript "0" represents the initial value. The block diagram in Figure 4 shows the equation of motion (the swing equation) of the difference in rotor phase angles (δ sg and δ vsg ) in the two-machine system. ∆P sg and ∆P vsg denote the difference when the load is added between the mechanical input and the electrical output of the generator for the SG and VSG, respectively. The mechanical input from the prime mover is received by the rotor of each generator. Equation (1) and each order term of the phase difference angles δ sg and δ vsg in Figure 4 represent the torque types, which can be decomposed as shown in Equation (3). Hereafter, the subscript "n" refers to either "sg" or "vsg": where K = K sg K vsg /(K sg + K vsg ) is the synchronization force generated when the two synchronous machines are operated in parallel, which can be represented by the feedback of an opposite force determined by multiplying the difference in the phase angle between the two machines δ sg − δ vsg by K. K sg and K vsg can be expressed by Equation (4): (4) These expressions can be simplified to K sg = 1/x sg and K vsg = 1/x vsg by assuming that the voltage deviations are small and that they offset the initial value of the phase difference angle. In that case, x sg and x vsg correspond to the transient reactance of the SG (x d ) and the virtual inductance of the VSG control (x), respectively.
∆P sg and ∆P vsg can be expressed as the sum of each type of torque, as shown in Equation (5); they are input by proportionally distributing the variable load power ∆P LD that corresponds to K sg and K vsg , as shown in Equation (6): The first-order differential of the phase difference angle of the rotor is the angular velocity of the system (at the reference voltage), and the second-order differential is the inertial torque, which has a relatively large value. Therefore, in the following sections, the mechanism responsible for the frequency deviations of the system during load fluctuation is analyzed using the differential equations of the phase difference angles (Equations (3)-(6)).

Analysis of Mechanical Torque
Numerical simulations were conducted with MATLAB/Simulink based on the connection diagram in Figure 3, and the resulting phase difference angles δ sg and δ vsg were used to calculate Equations (3)-(6), described in Section 2.3. However, because differentiation of the simulation log data would greatly amplify the errors in the simulation calculation, the data was decimated for approximately 30 ms. Even though there were few samples due to the limitations of the simulation algorithms and the lack of computer processing power, the trend of data deviations was not suppressed. Table 1 lists the main machine constants and control variables of the SG and VSG for the 13 tests that were conducted. In the standard test (#1), the total capacity, inertia, reactance, and droop coefficient of the VSG are approximately the same as for the SG. All analyses were conducted using unit capacitance. The rated capacitance of the SG and VSG on the simulation circuit in standard test #1 was set to 200 MVA as 1 pu, and a line impedance of approximately 300 km was added to buses A, B, and C. The governor time constant of the VSG (T vsg ) was fixed at 0 s according to the Kawasaki topology, and the droop coefficient of the SG (R sg ) was fixed at 2.5%. The relationship between the dumping coefficient D n , unit inertia constant M n , and damping time constant T dn was established with the equation T dn = 2M n /D n . Thus, the phase difference angle δ n (or sδ n ) estimated in the simulations was used to read the damping time constant T dn and to calculate the dumping coefficient D n . We also conducted a test in which the unit inertia constants of VSG2 and VSG3 were set to different values: 2.50 s and 5.00 s, respectively (results not included in table). Compared with test #1, in which the inertia constants of the two VSG machines were equal, only a slight difference was observed. Therefore, all subsequent analyses were performed using the results from the test in which the inertia constants of the two VSG machines were equal. Figure 5 shows the transient response of the internal phase difference angle between the SG and two VSGs in test #1. A load of 60 MW was added 45 s after the onset of the simulations. Figure 5a, in which an additional load was connected to bus B between the two VSGs, can be compared with Figure 5b, in which the loads connected to bus C between the SG and one VSG within a 0.5 s interval. The fluctuation in the onset time slightly before 45 s can be attributed to the necessary interpolations after the decimation of the log data. As shown in Figure 5, the phase difference angle of the SG (δ sg ) increases relatively uniformly at any load connection position (Figure 5a,b), and positive and negative vibrations of sδ n and s 2 δ n occur for approximately 200 ms after the load is added.  Figure 5. Transient responses of phase differences: δ sg and δ vsg are the phase differences of the SG and VSG, respectively; sδ sg and sδ vsg are respectively the first derivatives of δ sg and δ vsg ; and s 2 δ sg and s 2 δ vsg are respectively the second derivatives of δ sg and δ vsg . Figure 6 shows the classification of the mechanical torque calculated using the results shown in Figure 5 and Equations (3)-(6): the synchronization torque ∆P n K (red); the sum of the dumping and synchronization torques (∆P n D and ∆P n K ) (blue); the sum of the governor effect (∆P n R ) and the dumping and synchronization torques (∆P n D and ∆P n K ) (orange); and the sum of the inertial torque (∆P n M ), governor effect (∆P n R ), and dumping and synchronization torques (∆P n D and ∆P n K ) (green). The lighter lines are for the SG, the darker lines are for the VSG, and the thick lines indicate the total value of the SG and VSG.
The area between the red and blue lines is the dumping characteristic: because the lines are nearly overlapping, the dumping characteristic can be considered to be minimal. The area between the blue and yellow lines represents the governor effect, with governor time constants for the SG and VSG of 0.25 s and 0 s, respectively. Only the VSG governor effect ∆P vsg R is expressed immediately after the load is added. The dumping and governor effects are minimal because sδ n is minimal relative to s 2 δ n , and the inertial torque acting on s 2 δ n becomes predominant in the time domain.
The inertial torque ∆P n M is represented by the area between the yellow and green lines. Both inertial torques ∆P sg M and ∆P vsg M oscillate between positive and negative for approximately 200 ms immediately after the load is added, and SG generates a larger inertial torque than VSG (the lighter green lines are larger than the darker green lines). The effect of the connection position of the additional load is reflected in the inertial torque, which is larger when the synchronous machine is closer to the additional load and smaller when it is farther away. This may be due to the increase and decrease in the phase difference angles δ sg and δ vsg resulting from the line impedance. Figure 7 provides a comparison of the total mechanical torque ∆P LD in tests #1 through #9. The parameters for tests #10-#13 (the droop coefficient R n and governor time constant T n ) are omitted because they cause minor changes in the time range for which the transient torque response is considered. Test #1 is used as a reference (purple). The results of changing the VSG unit inertia constant M vsg are shown in Figure 7a (red), and the results of changing the virtual impedance X vsg are shown in Figure 7b (blue). The dumping coefficient D vsg depends on the unit inertia constant M vsg ; this is not discussed here because the ratio of dumping effects is small, as mentioned above. The virtual impedance X n , which is composed of the reactance and half the value of its resistance, was also changed. In Table 1, X n refers to the transient reactance of the SG x d or virtual inductance x of the VSG, and the resistance is not listed. Figure 7. Comparison of the total mechanical torque: M vsg is the VSG unit inertia constant, and X vsg is the virtual impedance of the VSG. Tests #1 to #9 are the test cases listed in Table 1. Figures 5 and 6, the inertial torque becomes predominant at approximately 200 ms, after load fluctuations. Thus, Figure 7a,b indicate that changes in the VSG unit inertia constant M vsg (Figure 7a) have a greater impact on the mechanical torque than do changes in the virtual impedance X vsg (Figure 7b). An increase in the VSG unit inertia constant M vsg increases the total deviations of the mechanical torque ∆P LD (Figure 7a). Increasing the VSG virtual impedance value has a similar effect (Figure 7b). According to Equation (3), X vsg affects the coefficient that determines the synchronization force ∆P K , but the contribution of ∆P K to the total torque is minimal, as indicated by the red line in Figure  6. That is, in the change in X vsg , it is presumed that the change of the output distribution in Equation (6) has a dominant effect on the mechanical torque ∆P LD rather than the change of the synchronization force ∆P K .

As shown in
The results of test #6 in Figure 7b also indicate that the mechanical torque begins to oscillate. In test #6, the VSG virtual impedance X vsg is set to its smallest value. Whether the oscillation can be attributed to the control stability limit of the Kawasaki VSG or to the resonance of impedance between the SG and VSG is debatable; this can constitute a topic for a future study. Figure 8 shows a comparison of the frequency deviations at the terminal voltages of each power supply. Although the droop and governor effects appear at the steady-state frequency, transitions from the transient period cannot be ignored. Hence, subsequent frequency evaluations were conducted to compare all the tests (#1 to #13). The evaluation period was set to 1.0 s after the load-power fluctuation; this exceeds the torque evaluation period of 0.5 s. Figure 8a shows the results for the cases in which the VSG unit inertia constant M vsg was changed (red); Figure 8b shows the results when the virtual impedance X vsg was changed (blue); Figure 8c shows the results when the droop coefficient R vsg was changed (green); and Figure 8d shows the results when the SG governor time constant T sg was changed (orange). The upper and lower graphs in Figure 8a-d show the frequency deviations ∆F sg and ∆F vsg [Hz] calculated from the terminal voltages of the SG and VSG, respectively. The positive and negative signs are inverted in Figure 8 to match the plots of mechanical torque in Figures 5-7 in the previous section.  Figure 8a shows that for approximately 50 ms after the load is added, there is no difference in test conditions #1-#5 for the frequencies of both the SG and VSG. During the subsequent 100 ms, the frequency deviations of the SG and VSG can be recovered by setting a large unit inertia constant M vsg . This finding is consistent with Figure 7a. Meanwhile, as shown in Figure 8b, there is a minor difference in the frequencies of the test cases immediately and approximately 150 ms after the load is added. Although this result does not appear in the order of magnitude X vsg , the frequency deviations of the SG and VSG can be recovered by setting a larger VSG virtual impedance X vsg , which is also consistent with the results in Figure 7b. These results show that the impedance has a profound influence on the frequency deviation at the output terminals of the SG and VSG immediately after the load fluctuation and at intervals of 50 ms in subsequent periods. It may be that difference in the frequencies of the test cases are not observed 50 ms after the load is added (particularly for Figure 8a), due to the influence of the actual impedance (which affected all test cases equally) rather than that of the inertial torque shown in Figure 7a.

Analysis of System Frequency
In test case #10 in Figure 8c, the VSG supplies most of the load power in the steady state, and the frequency deviations during this period are smallest for both the SG and VSG. Meanwhile, the frequency deviations are large for both the SG and VSG when the proportion of the load power in the steady state is large in the SG, as in test case #11. The differences, in the effects of changing the VSG droop coefficient (Figure 8c) and SG governor time constant (Figure 8d), are not confirmed immediately and 100 ms after the load is added and gradually increase during the transition period. Finally, they can be clearly confirmed in their steady states. These tendencies are similar to Figure 8a.
The frequency deviations in both generators diminish as the SG load sharing decreases (Figure 8c) and the SG governor delay time decreases (Figure 8d). This is because the swing (slack) bus is set such that the imbalance between the generated power and the load power of the entire system is compensated by the SG. It can be said that the mechanical torque lost from the swing node in the transient period after the load-power fluctuation determines the frequency deviation in the entire system. Thus, compensating for the lost mechanical torque at an early stage leads to an improvement in the frequency deviation.
We therefore introduce the concept of the frequency deviation of the center of inertia ∆F c [Hz] [42], which can be considered to be the frequency deviation that exists as a common component of the entire power system. The frequency deviation of the center of the inertia ∆F c [Hz] can be obtained from Equation (7) [42]; the results are shown in Figure 9. Herein, ∆F sg and ∆F vsg [Hz], W sg and W vsg [VA], and M sg and M vsg [s] are the frequency deviations by the terminal voltages, the rated capacities, and the unit inertia constants of SG and VSG, respectively, F 0 [Hz] is the nominal frequency: As in Figure 8, the results associated with changes in the VSG unit inertia constant M vsg are shown in Figure 9a (red); those associated with changes in the virtual impedance X vsg are shown in Figure 9b (blue); those associated with changes in the droop coefficient R vsg are shown in Figure 9c (green); and those associated with changes in the SG governor time constant T sg are shown in Figure 9d (orange). The positive and negative signs are inverted, as in Figures 7 and 8. Figure 9 can be regarded as depicting the center of gravity of the frequency deviations of both the SG and SVG in Figure 8. This can be understood based on the fact that ∆F c [Hz] is calculated from the rotor kinetic energy by considering the size of each power source. As seen in Figure 8, because the frequency measured by the terminal voltage varies with respect to the capacity and parameters of each power supply, it is difficult to obtain a unique index of the frequency of the entire system. With regard to the aforementioned factors, we propose utilization of the concept of the center of inertia frequency and use it as a representative of the frequency of the entire system. In the following section, the indices of the deviation and the rate of change of the center of inertia frequency are measured for the purpose of stabilizing the frequency of the entire system. The correlations between the indices and the generator parameters are also provided.

Indices of the System Frequency Deviation Based on the Parameters and Time Range
We calculated the correlation between the generator control variables and equipment constants and the frequency deviation and rate of change of frequency. Figures 10 and 11 show the frequency deviation and rate of change of frequency (RoCoF), respectively.  Figure 10 shows the maximum frequency deviation in each interval calculated from Figure 9. Similarly, in Figure 11, the RoCoF in each interval is calculated from Figure 9. The calculated maximum frequency deviation and RoCoF are approximated as those in the middle of the intervals. Tables A1-A4 in Appendix A list the maximum frequency deviations of the intervals for all parameters used in the calculations for Figures 10 and 11, and Tables A5-A8 in Appendix A list the RoCoFs of the intervals. In Figures 10 and 11, (a) shows the result of changing the VSG unit inertia constants M vsg , (b) shows the result of changing the virtual impedance X vsg , (c) shows the result of changing the droop coefficients R vsg , and (d) shows the result of changing the SG governor time constants T sg . Figure 9 shows the results of changing the parameters centered on the standard test (#1, in which the two machines have almost the same capability). Therefore, in Figures 10 and 11, the VSG to SG ratio of each parameter is displayed on the Y-axis. However, because the VSG governor time constant of the Kawasaki topology is zero, the SG governor time constant T sg is displayed on the Y-axis in Figures 10d and 11d.
Both the time response of the frequency deviations and their response to changes in the main parameters can be understood simultaneously in Figures 10 and 11. It can be seen from Figures 10a and 11a that the inertia constant ratio and the reactance ratio strongly affect both the frequency deviation and RoCoF in this time domain. For a small inertial constant ratio, a small frequency recovery is observed immediately after the load power fluctuates, and the absolute value of RoCoF is large. However, after the frequency nadir, the recovery effect is large and the absolute value of RoCoF is small compared to the case where the inertial constant ratio is large. The opposite tendency of the reactance ratio can be observed in Figures 10b and 11b. After the load power fluctuation, the smaller the reactance ratio, the smaller the frequency fluctuation and the larger the absolute value of RoCoF. However, this difference disappears over time. It is worth noting that both the inertia constant and the impedance ratios show tendencies of mountains and valleys when the ratio value is between 1 and 1.5, especially in the period immediately after the load power fluctuation. Alternatively, Figures 10c and 11c show that when the droop ratio is small, the frequency deviation reduces, and the RoCoF fluctuation increases. From  Figures 10d and 11d, it can be inferred that the smaller the governor time constant of the SG, the smaller the frequency deviation. However, there is no significant difference in the RoCoF.
Because both the frequency deviation and the RoCoF are approximated by the values at the midpoint time of each interval, the results do not exactly match the actual measured values. The accuracy might be improved by lengthening the time period and/or by increasing the number of intervals. We plan to re-test with an increased accuracy in the future.

Indices of the System Frequency Deviation Based on the Frequency Characteristics
The frequency deviation index of each parameter in the time domain, calculated in the previous section, is examined from another angle in this section. The frequency characteristic from the load power fluctuation to the frequency fluctuation measured by the terminal voltage of each generator (G vsg (s) and G sg (s)) can be obtained from Equation (8) [30,31]. Here, as in Figure 1, G vsg (s) = −1/(sM vsg + 1/R vsg ) and G sg (s) = −(sT sg + 1)/(s 2 M sg T sg + sM sg + 1/R sg ) are stand-alone responses of the VSG and SG, respectively, and as explained above, Equation (4) can be simplified to K vsg = 1/X vsg and K sg = 1/X sg . ω 0 is the rated angular frequency: To achieve an index for controlling the entire system, we introduce the center of inertial frequency in place of the frequencies by terminal voltages. The frequency characteristic of the center of inertial frequency can be derived from Equations (7) and (8); it is expressed in Equation (9) characteristic of the center of inertial frequency can be derived from Equations (7) and (8); it is expressed in Equation (9): Using Equations (10)- (12), G vsg (s) and G sg (s) in Equation (8) can be rewritten as Equation (13). Here,Ĝ vsg (s) andĜ sg (s) can be considered as the transfer functions of the frequency deviations from the load power fluctuation when the VSG and SG are connecting to an infinite bus. It is clear thatĜ vsg (s) andĜ sg (s) are the key factors that determine the characteristics of G vsg (s) and G sg (s): The Bode diagram in Equation (9) is shown in Figure 12. Similar to Figures 8-11, Figure 12a shows the results of changing the VSG unit inertia constant M vsg ; Figure 12b shows the results of changing the virtual impedance X vsg ; Figure 12c shows the results of changing the droop coefficient R vsg ; and Figure 12d shows the results of changing the SG governor time constant T sg .  Figure 12a shows that the frequency deviation can be suppressed by setting a large unit inertia constant for the VSG, as in Figure 10a. At the same time, it is evident that its phase is delayed in the band of approximately 1.0 to 10 rad/s (0.16 to 1.6 Hz). Furthermore, although interference between the two machines is seen in the band of approximately 10 rad/s (1.6 Hz), the gain in that band is sufficiently attenuated, so this does not pose a major problem. Conversely, in test cases #2 and #4, where (M vsg /M sg ) < 1, the parallel operation of two machines in the frequency band where the gain is amplified must be undertaken carefully. As evident in Figure 12b, the effect of the virtual impedance of the VSG is small. It is well known that frequency deviation can be suppressed in the low frequency band (steady state) when the droop coefficient of the VSG is extremely small and the load sharing ratio is large. However, in test case #10 in Figure 12c, the gain is greatly amplified in the frequency band of approximately 10 rad/s (1.6 Hz) where the interference of the two machines occurs. Figure 12d shows that in the low frequency band (approximately 1.0 rad/s (0.16 Hz)), setting the governor time constant of the main power supply causes the gain to be amplified and the phase to lead. Figure 13 shows the poles and zeros of G c (s), G vsg (s), and G sg (s). As in previous figures, Figure 13a-d respectively show the results for the cases in which the VSG unit inertia constant M vsg , the virtual impedance X vsg , the droop coefficient R vsg , and the SG governor time constant T sg are changed. The red, blue, and green symbols denote the poles and zeros of G c (s), G vsg (s), and G sg (s), respectively. Figure 13. Poles and zeros of the transfer functions: the light blue, yellow and pink arrows/circles indicate that the parameters are increasing; red represents the deviation of the center of inertial frequency (G c (s)); blue denotes the frequency deviation ratio of the SG to the VSG ( G vsg (s)); and green denotes the frequency deviation ratio of the VSG to the SG ( G sg (s)).
The three significant types of oscillating poles of G c (s) can be confirmed, and they are highlighted by light blue, yellow, and pink arrows/circles. The arrows indicate the direction in which the parameters are increasing. In the tests in which M vsg , X vsg , and R vsg change, G c (s) has a natural frequency band (the poles in Figure 13a-c marked by light blue arrows) that is strongly influenced by G vsg (s), and conversely, in the test in which T sg changes, it has a natural frequency band (the poles in Figure 13d marked by pink arrows) that is affected by G sg (s). The poles in Figure 13a-d marked by yellow arrows/circles are related to both the VSG and SG. As shown by the arrows in this figure, when the parameters are increased, the poles approach the origin. This means that even though the oscillation frequency reduces, the oscillation tends to continue. The tendency of oscillation due to the specified parameters is difficult to judge from the transient response in Figure 9, and it is also difficult to judge the time required for damping from the Bode diagram in Figure 12. Based on the aforementioned factors, it can be inferred that the pole zero map of Figure 13 effectively shows the indexed graphs of Figures 10 and 11 from another angle.

Conclusions
In this study, a synchronous generator was used as the main power source, and multiple grid-connected inverters equipped with VSG controls were connected to the IEEE 9 bus system. The mechanical torque characteristics of each generator were decomposed and analyzed to clarify the mechanism through which each type of torque causes system frequency deviations. The correlation between the time, each parameter (control variables and equipment constants) of the main generator and grid-connected inverters, and the frequency responses was also analyzed. The findings of this study are summarized as follows: - In systems that include smart inverters (e.g., VSG inverters) in addition to the conventional synchronous generator, various mechanical torques in each power source cause fluctuations in the system frequency. -Because the frequency at the terminal voltage of each power supply of the multimachine system is strongly influenced by the equipment constants and control variables of the power supply, we introduced the concept of "center of inertial frequency." -Both the time and frequency domains of the system frequency, which are affected, differ with respect to various mechanical torque types (damping, inertia, synchronization torques, and governor effect). -Frequency deviations that strongly depend on both device constants/control variables and time were indexed using 3D graphs. Furthermore, to evaluate the oscillation intensity and the damping time, the pole arrangement of the system was studied. -A multi-machine system can be reduced to a system with only two machines consisting of a main power supply and other power supplies. The frequency fluctuation system is determined through the ratios of the equipment constants/control variables of those two machines.
It has been previously confirmed empirically that a grid-connected inverter with the characteristics of a synchronous machine can stabilize a system as the main power source in a small-scale system or an isolated MG. However, a similar problem has recently been identified in large-scale smart grids and the magnitude of faults and their contributions in the grids were not quantified. Thus, that uncertainty has impeded the introduction of distributed power sources. In this study, therefore, a quantification index that applies the two-machine contraction model is utilized to suitably set the parameters of multiple grid-connected inverters, making it possible to predict and control their effects on large smart grids. It is expected that this knowledge will contribute to accelerating the introduction of distributed power sources originating from renewable energies, thereby reducing carbon emissions.