Dynamics Analysis Using Koopman Mode Decomposition of a Microgrid Including Virtual Synchronous Generator-Based Inverters

In the field of microgrids (MGs), steady-state power imbalances and frequency/voltage fluctuations in the transient state have been gaining prominence owing to the advancing distributed energy resources (DERs) connected to MGs via grid-connected inverters. Because a stable, safe power supply and demand must be maintained, accurate analyses of power system dynamics are crucial. However, the natural frequency components present in the dynamics make analyses complex. The nonlinearity and confidentiality of grid-connected inverters also hinder controllability. The MG considered in this study consisted of a synchronous generator (the main power source) and multiple grid-connected inverters with storage batteries and virtual synchronous generator (VSG) control. Although smart inverter controls such as VSG contribute to system stabilization, they induce system nonlinearity. Therefore, Koopman mode decomposition (KMD) was utilized in this study for consideration as a future method of data-driven analysis of the measured frequencies and voltages, and a frequency response analysis of the power system dynamics was performed. The Koopman operator is a linear operator on an infinite dimensional space, whereas the original dynamics is a nonlinear map on a finite state space. In other words, the proposed method can precisely analyze all the dynamics of the power system, which involve the complex nonlinearities caused by VSGs.


Introduction
To achieve environmental friendliness and sustainability, various types of power sources are freely connected to power systems. Moreover, distributed energy resources (DERs) such as batteries must be connected to use renewable energy, thereby avoiding fossil fuel depletion, and to shield the power system from accidental power outages caused by natural disasters such as hurricanes and earthquakes. It is a major challenge to maintain system stability in microgrids (MGs), which include small DERs with autonomous power, unlike traditional power systems, which are powered by large synchronous generators (SGs) with large inertial forces and short-circuit ratios (SCRs). Because the SCR is determined by the system voltage, impedance, and power supply capacity, the SCR of a small DER connected to a low-voltage distribution system cannot be increased. Furthermore, to protect the equipment from accidental overcurrent or overvoltage, a conventional grid-connected inverter, which cannot have inertial force, is controlled to respond immediately.
Numerous studies on new methods of grid-connected inverter control are underway to overcome the problem of inertial force. In particular, virtual synchronous generator (VSG) control of grid-connected inverters is attracting attention. VSG control can contribute to frequency stabilization by using battery energy as the inertial force of the SG [1][2][3][4]. VSG inverter controls, also called "smart inverters," were proposed at approximately the same time [1,2]. In the corresponding methods, the inertial control and voltage control of an SG supplement the conventional current control method installed in power conditioning systems. In contrast, Zhong et al. [3] based their research on the voltage control method and faithfully reproduced Park's generator theory with a controller. Meanwhile, Hirase et al. [4] analyzed the transient response of the VSG control and frequency response.
The grid-connected inverters that are used when DERs are connected to a grid are required to comply with the grid codes, although the methods of achieving compliance differ by manufacturer and are not disclosed. Because it is more complicated to analyze a system that includes a large number of black-boxed smart inverters, as described above, research is also underway to make the grid code applied to such inverters necessary and sufficient [5][6][7][8]. In [5,6], the correlations among frequency fluctuation, time, and control parameters were quantified by presenting analyses of the transient response of a VSG inverter based on the generator mechanism. Hashimoto et al. [7] presented a test method for obtaining comprehensive indices of the characteristics of various residential singlephase photovoltaic (PV) inverters. In addition, the use of power/control hardware-in-theloop (PHIL/CHIL) is being actively promoted for the testing of various power supplies and interoperable controllers [8]. Nevertheless, both types of analyses model and approximate the system dynamics within the controller and simulate its behavior.
In recent years, power semiconductor materials such as silicon carbide (SiC) and gallium nitride (GaN), and communication technologies such as the fifth-generation (5G) mobile communication system, have improved dramatically, and breakthroughs have also occurred in power system development. Previously, communication delays occurring in the supervisory control and data acquisition (SCADA) were fatal for controls requiring high speed, such as transient response controls. However, we can now analyze and control data measured in real time using the aforementioned innovative technologies. This method is called the "data-driven approach" [9][10][11][12][13][14]; it is characterized by calculation of the dynamics to be analyzed and controlled based on the acquired data. Park et al. [9] estimated the load consumption in a building by modeling based on an unsupervised learning method, as well as the demand response. Similarly, Gorjão et al. [10] calculated a stochastic model from the obtained data (market rules, consumer behavior, etc.) and performed a frequency analysis. Xinhui et al. [11] improved the calculation speed of the above methods by estimating the state from past data. Li et al. [12] estimated a dynamic model using synchronous phase measurements, whereas Almunif et al. [13] identified low-dimensional eigenvalues of power systems and enhanced noise immunity. Liu et al. [14] improved the robustness and calculation speed of the power system by evaluating the static voltage stabilization margin (VSM).
The key points of the preceding discussion can be summarized as follows. First, the analysis of a power system is complicated when numerous DERs are connected to it. Next, DERs are connected via grid-connected inverters, and conventional inverters lack stability owing to their responsiveness. Meanwhile, sophisticated autonomous controls for system stabilization, such as the VSG controls in smart inverters, are not publicly available. Nevertheless, thanks to 5G technology, the development of data-driven system analysis and controls that attempt to simultaneously perform system measurement and control is progressing. However, the dynamics of complex power systems are high-order, and the spectral characteristics of the dynamics can be compromised by reducing the order and approximating the dynamics in a data-driven approach.
To handle the high-order and nonlinear dynamics, a dynamic analysis method that uses the spectral characteristics of the linear Koopman operator (KO) has been proposed [15]. Although the introduction of the KO dates back to 1931 [16], the computational technology of the time hindered its application in industry. Dynamic mode decomposition (DMD) using the KO is called Koopman mode decomposition (KMD), and it calculates the eigenvalues (KEs) and modes (KMs) of KOs using empirical values obtained from nonlinear dynamics, even if the dynamic operator itself is not clear [17]. In recent years, increasingly many mathematical studies focusing on the spectral characteristics of the KO have appeared [17][18][19][20][21][22][23][24][25][26]. Cho and Yamazaki [18] classifies the properties of a composition operator, which is another name of KO. Williams et al. [19] presented a data-driven method for approximating the leading KEs and KMs. Surana [20] demonstrated forecasting and anomaly detection using nonlinear dynamic generative modeling of a time series based on the KO. Meanwhile, Lange et al. [21] described a spectral decomposition method that is applicable to nonlinear systems. The Fourier transform applied to a linear system was extended using KO, without assuming periodicity, thereby facilitating long-term predictions based on data measured at arbitrary sampling intervals. Arbabi and Mezic [22] demonstrated the spectral decomposition of KO in the DMD by assuming that the underlying dynamical system was ergodic. They also showed that the singular value decomposition, which is central to most DMD algorithms, converges to the appropriate orthogonal decomposition of the observable. Takeishi et al. [23] and Sinha et al. [24] described a KO approximation method used for observables with large measurement noise. Takeishi et al. [25] also introduced an application of KMD to a neural network and showed that KO has an invariant subspace even when an observable that spans it is not available. Further, Črnjarić-Žic et al. [26] discussed the characteristics of the KEs and KMs for the stochastic KO in various types of linear random dynamical systems (RDS).
Because of the high level of accuracy of its nonlinear dynamics estimations, the KO has been utilized for actual analyses. Among such analyses, those described in [27][28][29][30][31][32] are related to power systems. Susuki and Mezić [27] applied KMD to two major sets of actual accident data (the 2011 Arizona-Southern California power outage and the 2006 European interconnection grid system failure) to detect unstable power flow patterns. Susuki et al. [28] also described coherency identification regarding the swing of SGs that interfere with each other. Meanwhile, Jlassi et al. [29] determined the minimum number of calculations required for KMD convergence, and the results of the discrete Fourier transform (DFT) and KMD were compared to identify the nonlinear characteristics confirmed by KMD. Later, Korda et al. [30] introduced a method of combining KMD analysis results with model predictive control (MPC). Hernandez-Ortega and Messina [31] used KMD to improve the approximation accuracy when reducing the order of a nonlinear perturbation analysis model to an arbitrary order. Finally, Cassamo and van Wingerden [32] described a method of applying the KO to the aerodynamic analysis of a wind power generation system to reproduce a low-dimensional linear state-space model accurately, thereby reducing the design cost of the controller.
With the above research as a background, we focused on both the voltage stability and frequency stability of a power system in this study and attempted to extract the characteristics of their mutual interference using KMD. While prior literature [27,28] analyzes the steady-state power fluctuations of large-scale power systems in detail, this study focuses on the transient response in MGs including DERs and analyze the nonlinearity in transient response caused by their controls. The characteristic extraction results were visualized by widening the parameter range of the VSG control inverter. The test system was an Institution of Electrical and Electronic Engineering (IEEE) nine-bus model that included a main power supply and two distributed power supplies, which were inverters with storage batteries. The inverters were controlled by VSG control. The transient responses and DFTs of the system frequencies and terminal voltages were also calculated and compared with the results of KMD, thereby confirming the effect of cross-interference in the VSG control that appears only in KMD. The remainder of this paper is structured as follows. Section 2 reviews the KO and KMD. Section 3 explains the simulation test system used in this study. Section 4 presents analyses of the frequency response of the KMs when applying KMD. Finally, Section 5 summarizes the study.

Theory of the Koopman Operator (KO) and Koopman Mode Decomposition (KMD)
We briefly review the concept of the KO and KMD calculation method, which yields approximate values of the KEs and KMs. Further details are presented in [19,27,28].

Koopman Operator (KO)
Let M be a state space and : a vector-valued function. Equation (1) presents the dynamics of the state variable vector t M  x , and Equation (2) presents the dynamics at time t (2)) of a continuous system into a discrete system. In a power system with various types of generators, F t is typically nonlinear.
To define the concept of the KO, let G be the space of the observable. The linear operator K is defined on G , which maps the observable : g M   to the new observable, as shown in Equation (4). As evident in the power system, any element, such as the frequency, phase, or voltage, at the terminal of each generator can be selected. : The observables are obtained in a discrete-time system, and the observable ( ) t g x with dynamics F can be rewritten as Equation (5) from Equations (3) and (4): k K is an operator that represents the time evolution of the observable g up to time t , and it is referred to as the KO. The KO K is a linear operator on an infinite dimensional space G even though the original dynamics F are a nonlinear map on a finite state space M . The expression of nonlinear dynamics using linear operators could be employed to analyze the factors contributing to system instability that are attributable nonlinear dynamics without linear approximation. In addition, it is useful to divert control library assets that are commonly applied to linear dynamics.

Koopman Mode Decomposition
The KO K has eigenfunctions j G   that correspond to eigenvalues j   that satisfy Equation (6). j  and j  are called Koopman eigenvalues (KEs) and Koopman eigenfunctions (KEFs), respectively: Consider the m-dimensional vector  span the whole G , then any observable i g G  can be represented by the Fourier series as follows: Using Equations (3), (4) and (6), the time evolution of the i-th observable ( ) i k g x at the time step k derived from (7) can be expressed as Equation (8). Equation (9) demonstrates the decomposition of the observable vector g into an infinite number of oscillation modes, which is called KMD. Each oscillation mode is represented by the product of KEs, KEFs, and the complex vectors 1 2 , , , , which are called Koopman modes (KMs). Figure 1 illustrates the concepts of the KO and KMD.
where 0 x is the initial state variable vector. :  is defined as follows:

Arnoldi Algorithm
It is difficult to know the nonlinear dynamics F in advance, because DERs are repeatedly added and disconnected at arbitrary times. That is, it may not be possible to obtain the theoretical formula for KMD, because the KE and KEFs generally cannot be obtained. It is thus necessary to acquire the observables with minimal delay and to identify the system utilizing them. In this study, the Arnoldi-type algorithm was applied so that both the KEs and KMs could be calculated from time series observables without knowing the dynamics F or KO K . 1 N is a finite integer indicating the number of time samples, and we consider the The coefficients j c that satisfy Equation (10) can be calculated when the vector r is orthogonal to the space spanned by In particular, from the inner-product , ( ) 0 j c can be obtained by solving Equation (11). (11); that is, : Ac b holds, and by and can be considered equivalent to KEs ( j When the matrix  is obtained by using Equations (13) and (14) holds

the KEFs and
KMs) in Equation (9), where U represents the Vandermode matrix and j   are the empirical Ritz values: Because the differences between the empirical Ritz values and true values and between the finite and infinite series are small in the power system investigated in this study, we do not discuss how they differ here and rewrite j

IEEE Nine-Bus Model Test System
In this study, we considered the IEEE nine-bus model (three-machine system) ( Figure  2). One primary power supply was assumed to be a SG, and the other two generators were assumed to be distributed power sources of inverters with batteries which were equipped with the same VSG control. In the simulation tests, the SG model was that of the standard synchronous machine installed in the simulator. Figure 2 shows the rated capacities and voltages of the generators, bus voltages, line impedance, transformer capacities, and impedance. The total base load of 100 MVA was evenly distributed in the three locations of Buses 5, 6, and 8, and the simulation started with all the power being supplied from the main power supply (SG). A load of 100 MVA was added to Bus 8 in a stepwise manner, and the system fluctuation was analyzed immediately after adding the load. By adjusting the set point of the SG, the frequency and voltage just before the load was added were set to the rated values of 50 Hz and 13.8 kV. In the standard tests, the equipment constants and control parameters of all the generators were set to be nearly identical. In the other tests, the VSG control parameter was increased or decreased for comparative purposes. The main objective of this study was to investigate the instability caused by the interference between DERs (inverters) and the nonlinearity due to the interference of the control loop in each power supply.  Figure 3 shows a block diagram of VSG control. A Kawasaki VSG [2,4] was used in this study. The generator phase difference angle  [rad] was calculated from Equation (15), which corresponds to the block enclosed by the red rectangle in Figure 3 (17) is the impedance model, which represents the relationship among E f , V g , and the command value vector * I of the inverter current. The impedance r and x [pu] constitute a pseudo-impedance set in the VSG control software and 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 to the virtual resistance r was empirically set to 2  x r .The current command value was limited to avoid overcurrent, and was input into the minor control, which was composed of current control. The pulse width modulation (PWM) command value ( * abc V ) output from the minor control was input into the valve control. Because these controls are conventional, detailed explanations are not included here.

Simulation Tests
Numerical simulations were conducted using MATLAB/Simulink based on the connection diagram shown in Figure 2. The SG connected to Bus 1 and the VSGs connected to Buses 2 and 3 are referred to as SG1, VSG2, and VSG3, respectively. The parameters of each power supply ( PF T , etc.) are indicated by a subscript number _ n ( n = 1, 2, 3), and Table 1 lists the parameters for the 20 tests. In the standard tests (3,8,13,18), the VSG3 parameters are almost equivalent to those of the SG1 and VSG2. In all tests (1-20), p K and i K in the PI compensator of the AVR are 10 and 80, respectively. All analyses were conducted using unit capacitance with the rated capacitance of each power supply set to 100 MVA at 1 pu. In tests 1-10, an active load of 100 MW was added 15 s after the simulation onset. Similarly, in tests 11-20, a reactive load of 100 MVar was added.

Transient and Frequency Response Analyses of Simulation Results
Figures 4-7 show the simulated transient responses of the three power supplies obtained from tests 1-5, 6-10, 11-15, and 16-20, respectively. In each figure, (a-f) depict the results for SG1, VSG2, and VSG3, respectively. The left column (a,c,e) and the right column (b,d,f) present the results of frequencies and voltages measured at the terminal of each generator. Each observable (system frequency or terminal voltage of N = 3000 points) was measured for 10 s from the onset time (15 s). A sample rate capable of measuring a rated system frequency of 50 Hz was required, and in this study, the sample rate was set to 250 Hz (sampling every 4 ms).
In tests 1-10, the system frequency changed according to the drooping characteristics owing to the addition of the active power load, and transient fluctuations due to the swing characteristics of the generators were confirmed. The swing characteristics of the generators reflect the relationship between the rotor angular velocities of the generators and the input/output torques. Thus, the addition of active power significantly affects the system frequency (a,c,e). However, the terminal voltage (b,d,f) is also affected. The high-frequency vibration of the terminal voltages and sudden fluctuations at 15 s are due to the load switch, the voltage controls, and the phase-locked loops (PLLs) in the VSG controls that were used to convert the control signals to a dq reference frame.
Conversely, the results of tests 11-20 reveal the voltage drooping characteristic resulting from the addition of reactive power. Unlike the system frequency, the terminal voltage differs for each measurement point. Therefore, the parameter differences among tests 11-20 considerably affected the terminal voltage of VSG3. Simultaneously, transient fluctuations and accelerations in the system frequencies due to the addition of reactive power were also confirmed. In tests 16-20, Table 1.  Table 1.  Table 1.  Table 1.
As evident in the transient responses in Figures 4-7, natural vibrations occurred in the low-frequency bands (0.1-10 Hz) and high-frequency bands near the system fre-

Koopman Mode Decompositions (KMDs) of Power System Observables
The system frequencies and terminal voltages (Figures 4-7) measured at the output terminals of each power supply shown in Figure 2 were used as the observables , to which the Arnoldi algorithm described in Section 2.3 was applied to perform KMD. Subsequent analyses were performed assuming that the terminal voltages and system frequencies, which are the observables of the target power system, can be represented as shown in Equation (7). There were N = 3000 measurement points (sampled every 4 ms for 12 s), and at each measurement point, the observable vector ( ) j g x ( 0 j N   ) consisted of the components 1 ( ) j g x , 2 ( ) j g x , and 3 ( ) j g x , which were the measurement values of SG1, VSG2, and VSG3, respectively. Using Equations (13) and (14), we calculated KE j  , KM j v , and the vibration frequency corresponding to each KE j  at each measurement point ( 0 j N   ). KEs j  were complex scalars from which the arguments and vibration frequencies could be obtained, whereas KMs j v were complex vectors, where the KM components 1j v , 2 j v , and 3 j v corresponded to 1 ( ) j g x , 2 ( ) j g x , and 3 ( ) j g x , respectively. Hereafter, by specifying "SG1", "VSG2", and "VSG3", the subscripts 1, 2, and 3 will be omitted.  Figure 12A(a,c,e) do not appear in Figure 8a,c,e. It was also confirmed in other tests that the power spectrum belonging to Group 2 was caused by the discontinuous changes of the observables due to the disturbance (load power addition) that occurred 15 s after the start of the simulation. As described in the previous section, the difference can be attributed to the fact that Figure  8 represents a system that is linearly approximated around the equilibrium point using the DFT, whereas in using KMD, Figures 12A and B represent the exactly equal (not approximate) system containing both linear and nonlinear elements.
In the controls of the SG and VSGs, the active power and system frequency exhibit a control loop with the same mechanism as in Equation (15), and the control loops between the reactive power and terminal voltage can be expressed by Equation (16). Furthermore, these control loops cross-interfere, as shown in Equation (17). In tests 1-20, the effects of the active power change on the terminal voltage and of the reactive power change on the system frequency, which correspond to the cross-interference, were also confirmed in Group 4. The cross-interference in Group 4 does not appear when using the DFT, but appears when using the KMD; this can be reasonably confirmed in the right columns (b,d,f) for tests 1-10 ( Figures 12A and 13A) and in the left columns (a,c,e) for tests 11-20 ( Figures  14A and 15A). These tendencies indicate that KMD can express the interference between cross-loops relatively clearly.  When the values of the main control parameters change, the amounts of vibration of the observables change slightly. In particular, for tests [16][17][18][19][20], no significant differences in the characteristics were observed within the range of parameters listed in Table 1. Furthermore, the system operated stably in all tests and did not collapse. In future research, it will be important to increase the parameter ranges and to perform a sensitivity analysis.

Conclusions
As the introduction of DERs using REs is accelerating, the analysis of power system dynamics is becoming increasingly complicated due to the following factors:


The inverters that connect DERs to power grids comprise control loops for active/reactive power and frequency/voltage, and thus interfere with each other and exhibit nonlinear characteristics.


The system stability contribution function for the synchronization power of a smart inverter is an interference factor between power supplies, and it is difficult to perform mathematical analyses for complicated connection configurations.  Functions such as the fault ride-through (FRT) in the grid code require fast responsiveness from the grid-connected inverters, which depends on the detailed data acquisition and accurate numerical analysis of these operations.
Expectations have been increasing for data-driven analyses and control using data measured in real time. Therefore, to develop future data-driven control, we developed a new analysis method for the nonlinear dynamics of a power system using offline data. Frequencies and voltages, which are the observables of the main control loops, were obtained from simulation results obtained using MATLAB/Simulink and the IEEE nine-bus system. By applying linear KOs to these observables, we could analyze the nonlinear dynamics of the power system. The simulation system included two VSG inverters, which are classified as smart inverters, and one conventional synchronous generator as the main power source. The vibration frequencies were calculated from the KEs, and the frequency responses were confirmed through comparison with the DFT, which is a conventional analysis method that linearly approximates a system. The following conclusions are drawn based on this comparison:


The KMs calculated by using the Arnoldi algorithm in this study are approximations of the eigenfunctions of the KO. That is, when a system is series-expanded using the KEs and KEFs of the KO as the base, there are series terms represented by specific KEs and KEFs.  Even if the frequency bands of the KEs are almost identical, the absolute values of the KM components tend to be divided into large and small groups. The absolute values of both differ by several orders of magnitude.  The results obtained from the KMD method have KMs with large absolute values that cannot be confirmed by the DFT. Such power spectrum that does not appear in the DFT but appears only in the KMD can be considered to represent the nonlinearity in the dynamics. The simulation results indicate the power spectrum contained in the observables (frequency and voltage) that fluctuate discontinuously when the load power added.  The AC power system is generally controlled in the dq reference frame. Therefore, the cross-couple effect due to the asymmetric term, which causes interference between the d-axis and q-axis control, must be determined. Although it is generally difficult to formulate the nonlinearity in the dynamics and analyze the control crossloops, KMD can clearly reveal these effects in the frequency bands at particular KEs.
In the simulation tests of this study, it was confirmed that changes in active power affect not only system frequency but also terminal voltage by changing the control parameters of the inverter. Similarly, changes of the reactive power effect the system frequency as well as the terminal voltage.
When KMD is utilized in a power system, it can be applied for fault prediction and control stability by focusing on new phenomena that are not addressed by conventional approximation methods. However, the following tasks remain:


Confirmation of how KMs with large absolute values behave in actual systems;  Indexing of the magnitudes of the absolute values of the KM components;  Quantification of the errors caused by observable measurement delays and sampling cycles;  Sensitivity analysis with KMs comprising a wider range of control parameters, compared to this study; and  Comparison of the computer load and calculation time with those of conventional methods such as DFT.
We will address these challenges in our future work.