Quantitative Mapping of Modeling Methods for Stability Validation in Microgrids

Although power electronic converters have been the key enablers for the the integration of renewable generation, heterogeneous controllers and lower inertia increase the complexity of microgrids, more likely to cause instability. Hence, it is significant to study the stability of microgrids using prospective modeling tools such as, bode plots, Nyquist plots, eigenvalue loci, etc. The modeling of microgrids is an important step of stability analysis, where state-space-based and impedance-based modeling are the two commonly used stability evaluation approaches in microgrids. Miscellaneous modeling methods have respective pros and cons, which have been investigated in the existing literatures to some extent. However, it is still critical to quantify the modeling techniques so as to formalize the stability validation in an intuitive way, which is not addressed in existing literatures. Therefore, in this paper, modeling methods for stability validation are mapped based on an order-indicated complexity, and a quantitative framework for the mapping is provided. The proposed framework can be instructive when modeling microgrid systems with different sizes, topologies or control strategies, etc., and it turns out to be well applicable as demonstrated by exemplified simulations and experimental tests.


I. INTRODUCTION
Power electronics have been providing enormous possibilities for energy conversion and integration of new types of loads [1], [2]. In microgrids (MGs), a large number of renewable energy sources (RES) are integrated and utilized through power electronic converters, including photovoltaic arrays, batteries and fuel cells. The power electronic converters have fast dynamics due to the underlying controllers, and the power generation profiles from RES are also less steady than synchronous generators. These two aspects have led to lower inertia and less robustness in MGs [2], [3], [4]. On the other hand, the controllers will bring about interactions among converters leading to complicated harmonics and resonances [5]. Hence, a comprehensive framework to categorically assess the stability of MGs still remains a challenge in this field.
Many reported studies are focused on the modeling and evaluation of stability in MGs. In most cases, the controllers of converters are considered as the most critical part, which is not common in the modeling of conventional power systems [6]. The controllers can be modeled as, e.g., state spaces [7], [8], transfer functions [9] or equivalent impedances [10], [11], [12], [13], but due to approximations, different modeling methods exhibit distinct accuracy levels. Basically, there are two commonly used modeling mechanisms, namely statespace-based modeling [7], [8] and impedance-based modeling [9], [10]. MGs are modeled into state-space matrices or impedance networks, and the stability of the MG is then characterized through eigenvalue loci [7], [8], Nyquist plots [9], [11] or Bode diagrams [11], [12], [13], etc. In practice, state-space-based modeling is more suitable for mathematical programming with modular and reusable forms, while impedances have clearer physical interpretation of frequencydomain characteristics.
However, in [7], [8], [9], [10], [11], [12], [13], the modeled systems are similar in their scales, such that the applied modeling methods do not show much difference in simplicity. For example, in [7], a system consisting of three identical converters is selected as the study case. If complicated systems like the IEEE standard systems [14] (with larger scale or complicated architectures) or the CIGRE benchmark [15] (with RES or various controllers) are considered, the computational difference between the modeling methods will be of much importance during stability analysis, considering the following factors: heterogeneity of controllers and the size of system matrix which could increase in a non-linear manner with the increase in the number of nodes. Under such circumstances, the modeling methods will imply different computation burden. Consequently, there have been studies focusing on the modeling and performance validations of large-scale MG systems [16], requiring many approximations and much advanced solvers to reduce the computation time.
Under this scenario, the comparison of the modeling methods becomes a practical and general concern. There are many aspects under discussion, including the procedures of modeling and stability evaluation in [17], [18], and the availability of modeling methods in [19]. However, the literatures are mostly qualitative classifications, where the borderline between modeling methods are specified only in certain applications such as black-box systems, without intuitive mapping of system properties. The validations can be consequently ambiguous when the system grows complicated with large scale or multiple RES.
In light of the above, this paper aims at formalizing an explicit framework for stability validation in MGs, by mapping the modeling methods in a quantitative way. In our previous study [20], the complexity of modeling has been discussed, which is an important concern in terms of computation in practice. The mapping regarding different applications has also been demonstrated with study cases. In this paper, the mapping of modeling methods is further explored, and quantitative metrics based on order of elements in the system are proposed for the first time. The definitions of the metrics are introduced, which quantify the computational efforts for each modeling method for a given system in a handy way, formalizing stability validations and thereby serve as a clear baseline methodology for benchmarking the modeling methods. The efficacy is demonstrated by simulations and experimental results.

II. MAPPING OF STABILITY MODELING IN MICROGRIDS A. STATE-SPACE-BASED AND IMPEDANCE-BASED MODELING
State-space-based modeling is particularly suitable for multiple-in-multiple-out (MIMO) systems, as it is in the form of state matrix and state variable vectors [7]. Based on Lyapunov definitions, a MG is stable when all state variables describing the system converge towards or around a stable operation point with time, which co-aligns well with the small-signal stability of the MG. Accordingly, the stability of the MG is characterized by the eigenvalues of the state matrix or by the Lyapunov functions. The eigenvalue analysis is widely applicable for linear or linearized systems and easy to perform, but for systems that are not suitable for linearization, Lyapunov functions are preferred as a more general candidate approach.
In impedance-based modeling, the converters with filters are modeled into a Thevenin circuit or Norton circuit, where the control loops can also be modeled into impedances [11], [13]. Taking the Thevenin equivalence as an example, the MG is regarded as a network of voltage sources and impedances, and the stability of the MG is thereby translated by the stability of the impedance network. There are many approaches to conduct the impedance-based modeling, and two widely-used ones are: (1) developing the transfer functions by constructing the admittance matrix of the network [9], and (2) partitioning the MG into two parts and employing the Nyquist criteria [10]. Subsequently, Bode or Nyquist plots are often used for impedance-based stability analysis, which are quite straightforward to reflect the behavior of MGs at different frequencies, e.g., in response to inputs or disturbances.
A comparative evaluation of both the modeling methods considering different aspects has been summarized in Table 1 [19], [20]. Both the two modeling methods have respective pros and cons in particular applications, but their differences need to be comprehended in terms of complexity and applications when both can achieve high accuracy.

B. QUANTIFICATION OF COMPLEXITY OF MODELING METHODS
For stability characterization in MGs, the modeling methods can lead to differences in calculation complexity, fidelity and intuitiveness in terms of applications, etc. [19], [20]. The complexity is mostly related to the differentials or integrals inside, defining the order of MGs. Higher-order systems always have higher dimension in the state matrix or more complicated transfer functions derived by impedances, which are timeconsuming in computation.

Remark I:
It is true that operations like matrix inversions and eigenvalue calculations can be time-consuming in large systems or some advanced modeling methods, but the differentials are generally playing the major role in complexity considering its large number of appearances in the modeling. But if there are repetitive or nested advanced operations, the complexity should be the total effect of all operations instead.
A MG normally consists of multiple renewable generations and various loads, which are mostly interfaced by converters and passive filters. When modeled within control bandwidth (several kHz), the behaviors of converters can be averaged by switching period and then linearized. The passive filters and controllers are then described into transfer functions or differential equations with respective orders.
Under this circumstance, the contributions of units in a MG (controllers or passive components, etc.) are listed in Table 2. In this paper, two-level DC-AC converters are considered, while multi-level converters could also contribute much to the total order and could be studied in details in the future. In Table 2, when a unit is synthesized as a transfer function, the order is the maximum of the numerator and denominator, or the maximum order of inputs and outputs in the corresponding differential equation.
To introduce the metrics for quantifying the modeling methods, it is assumed that there is a total of n nodes in the MG system, and there is one converter and up to one passive load at each node.
where, x MG is the vector of observed state variables, u MG is the input vector consisting of the reference of controllers, initial states or small-signal disturbances, y MG is the output vector indicating the behavior or performance of the MG, and A MG , B MG , C MG , D MG are the state matrices. If the total number of state variables is denoted as λ, then the vector of state variables x MG is in the dimension of λ×1.
where, V conv and I conv are respectively the voltage and current of converters at the point of grid integration. The converters are divided into two types dependent on the control references: Type I with voltage reference (e.g., grid-forming converters) and Type II with current reference (e.g., gridfollowing converters). G 1 and G 2 are the respective transfer functions regarding the references or disturbances, which are diagonal for local controllers without considering crossconverter coupling [9]. Y MG is the admittance matrix considering both lines and loads. μ is the dimension of the admittance matrix, which can be equal to m (e.g., in single-phase or symmetric single-in-single-out (SISO) cases), 2m (e.g., in the dq or αβ frame) or 3m (e.g., in the dq0 or abc frame).
In both state-space-based and impedance-based modeling, there are variables under observation and matrices indicating the relationship between the variables. The vector containing the variables under observation is the state vector, i.e., x MG in (1) or [V conv,I , I conv,II ] T in (3). One of the matrices can be regarded as the representative matrix, which describes the overall architecture and is normally the largest one, for example, A MG in (1) or Y MG in (2). The size of the representative matrix determines the space complexity of the modeling method in computation, as listed in [20].
Based on the basics and the assumptions above, two metrics are thereby defined to evaluate the modeling complexity, maximum-order complexity (MOC) and apparent-order complexity (AOC).
MOC is aimed at estimating a lower bound of the complexity of modeling methods when all couplings are neglected: Definition I: The MOC is defined as the sum of: a) Vector Part (VP): the total order of calculation directly applied to the state vector; b) Matrix Part (MP): the sum of maximum orders of each row of the representative matrix. Furthermore, AOC can be used to include the non-diagonal elements in state matrices or admittance matrices as well and to obtain the overall computational burden: Definition II: The AOC is defined as the sum of: in the representative matrix. The calculations of AOC and MOC for state-space-based (denoted as the subscript SS) and impedance-based (denoted as the subscript Imp) modeling are listed in Table 3. The MOC can also be regarded as the size of systems measured by orders instead of the number of converters, where the complexities of controllers are also included. When either advanced controllers are employed or the system consists of more nodes, the MOC will increase in a linear way according to the total order. On the other hand, the AOC for state-space-based modeling is actually equal to the MOC, since the state matrix only consists of scalars. In contrast, for impedance-based modeling, the order of impedances (both lines and loads) across every two converters are additionally included, which means AOC ≥ MOC. In this case, the VP only describes the complexity of the state vector, so it is not necessarily influenced by the coupling relationship between states. The MP, however, is used to look into the representative matrix, where the non-diagonal elements are the coupling. For example, in state-space-based modeling, the MP is always 0, indicating that the coupling among states will not likely influence the modeling complexity. Therefore, the MOC can be regarded as the case where the representative matrix is the sparsest, or MOC can normally be the lower bound of AOC with the same modeling method.
An example is given based on a multi-converter MG to illustrate the metrics. In Fig. 1, all m converters in parallel are controlled by droop controllers and double-loop voltage regulation. It is assumed that the control framework is uniformly deployed with proportional-integral (PI) controllers in the dq frame. Transmission lines between the converters can all be regarded as resistor-inductor (RL) in series, and the loads are RL loads.
In each converter, the total order involves: r The first-order power filters for droop controllers (including both active and reactive power), r The integration of ω in P-f droop, r The double-loop PI controllers, and r The LC filters.
Therefore, the total order of a single converter is obtained as: LPF and ∫ωdt in droop controllers Voltage and current PI controllers (dq) while, the order is obtained as (5) instead if only one of the dq components is considered (or neglecting the dq coupling, which is denoted as "single axis" in this paper). In this case, the admittance matrix is m×m instead of 2m×2m.   Accordingly, the MOC and AOC of the two modeling methods are calculated and listed in Table 4. Apart from the given modeling methods, there are advanced modeling methods such as harmonic-domain Toeplitz matrix in [24], wherein the complexity can also be evaluated similarly. It can be concluded that, the MOC is in the increasing rate of O(m) (i.e., being bounded by m multiplied by finite constants), or linearly related to the size of system. The AOC, which is related to the total number of differentials in the mathematical model, increases faster in impedance-based modeling and shows the major differences of the modeling methods. Technically, the AOC could be used for evaluating the relative computational burden of modeling methods, while the MOC could be used as a benchmark for e.g., improving the modeling for sparse matrices.

C. MAPPING OF THE MODELING METHODS BY COMPLEXITY
After the calculation of respective MOC and AOC, the modeling methods can be subsequently mapped to MGs with different sizes. Taking the case in Fig. 1 as an example, the increase of complexity can be plotted in Fig. 2. Another modeling method introduced in [10] is also compared together, namely the impedance-based modeling by partitioning a MG into two parts (the grid part, a voltage source with serial impedance Z g , and the load part, a current source with parallel impedance Z 0 ). For this method, it should be noted that [Z g , Z 0 ] T is regarded as the representative matrix, and the modeling is always under the single-axis condition.  The modeling methods show significant differences in large-scale systems. For example, if both the d and q components are modeled when there are 10 converters, the AOC of impedance-based modeling (510) is almost three times to that of the state-space-based modeling (168), which indicates that the state-space-based modeling is more appropriate. However, if the d and q components are symmetric, impedance-based modeling is simpler in return. More generally, the impedancebased method is more suitable for smaller systems, while the state-space-based method for larger ones, and the boundary can be defined by the proposed metrics.
Remark II: a) In this case, all converters are assumed to be with the same PI control scheme and filters, but this framework is   also applicable for heterogeneous controllers. If e.g., advanced controllers are employed, the pattern of the curve may vary. Based on Table 3, the number of converters and the number of the states will both contribute to the total computational burden as revealed by the proposed metrics.  b) The metrics are an estimation of the modeling complexity, while the actual complexity can also be decreased by e.g., the cancellation of zeros and poles in impedancebased modeling, which requires deeper inspection with specific parameters. Methods like model-order reduction by clustering of nodes [25], [26] can also reduce the complexity. c) If there are nested advanced operations, e.g., taking the inversion of impedance matrix to get the admittance matrix, the proposed metrics are still taking the major part of modeling complexity, but may not be linearly related to the computation time in practice. To reflect the differences of the modeling methods more clearly, the number of independent state variables under observation can be chosen as a base value of the complexity, denoted as the freedom or rank of the state variables. The AOC is thereby normalized into relative AOC. For example, in Fig. 1, if the output voltage and current of the m converters are the focus, the base value will be 2m. But when the load current is included, the base value will be 3m-1, where the Kirchhoff's equations will reduce the freedom of currents by 1. In principle, this base value reflects the relative scale of a MG system similar to the MOC, but it is only determined by the number of nodes (system architecture) and the observed variables, without focusing on the internal complexity (controllers, filters, etc.) of each node.
In the case where the base value is 2m, the relative AOC is plotted in Fig. 3.
In this way, the intersection points can be illustrated more clearly by the base value, and the increasing rate of relative AOC can show the scalability of modeling methods: The lower the increasing rate is, the higher this method allows adaptation for larger-scale systems.

D. MAPPING OF THE MODELING METHODS BY APPLICATIONS
Stability modeling methods may also vary by applications. The tools employed in the modeling procedure can be mapped by the respective accuracy in certain applications. Based on Table 1 and our previous study [20], a comprehensive mapping can be concluded in Fig. 4. The frequency under study and the linearization in the modeling are illustrated.
The modeling methods can thereby be properly linked with the scale of the MG system and the applications under study. The accuracy of a modeling method is the rationality of  possible modeling results, while the proposed metrics basically sort the eligible methods by complexity, measuring the amenity of the modeling methods. For example, the resonance in single-or double-converter systems can easily be revealed by impedance-based modeling with Bode plots. But eigenvalue loci in state-space-based modeling is more suitable for multi-converter system or cases regarding dq coupling where the number of states is double of that in symmetric cases.
It has to be noted that, the mapping is not unique in most cases, but it could serve as a clear and practical guideline for optimizing the modeling methods. There also exist cases, where certain modeling methods can best fulfill the requirement, e.g., phasor portraits for transient stability [27].

III. CASE STUDY FOR STABILITY VALIDATIONS
To illustrate the stability validation based on the proposed quantitative mapping, a simple case is selected for simulations and experimental tests. As shown in Fig. 5, the system consists of two parallel converters with droop controllers and LC filters. A resistive load is connected to the point of common coupling (PCC). The DC link is created by rectifying from a 230 V from the AC end. The objective is to verify the role of dq decoupling in the current loop (the feedforward terms related to ω 0 L f ), which is consequently enabled or disabled for comparing the differences in stability performances Parameters of the validations are listed in Table 5. The parameters of simulations and experiments are different due to available hardware configurations. In experiments, Case I is aimed to verify the steady-state performance, and Case II is focused on dynamics.
Based on the proposed complexity metrics, the AOC of state-space-based modeling is: The AOC of impedance-based modeling is: Therefore, state-space-based modeling is more appropriate in this case with less complexity, and eigenvalue loci is capable of evaluating the stability of the system. Similar to [7], the MG can be modeled into the following small-signal state space: where, x conv consists of the state variables of converters: x conv, k,dq = where, for the k-th converter, i s is the current of filter inductors, v o and i o are, respectively, the voltage and current at the output of LC filters, and E v and E i are the integral error of v o and i s , respectively. The subscript dq indicates that the MG is modeled in a dq frame synchronized at the PCC.

A. SIMULATION RESULTS
Based on the state space and the given parameters, the eigenvalue loci of the system can be plotted in Fig. 6. According to the eigenvalue loci, there are right-half-plane (RHP) poles when the dq decoupling is disabled, which indicates the potential instability of the system. Simulations based on PLECS is then conducted, as shown in Fig. 7.
In the simulations, when the dq decoupling is disabled, the waveforms are distorted with harmonics with the frequency around 1.1 kHz, which accords with the eigenvalue loci in Fig. 6, as there are RHP poles. It preliminarily illustrates the application of the proposed metrics for complexity evaluation.

B. EXPERIMENTAL RESULTS
The configuration of the experimental setup is shown in Fig. 8. Two converter racks are used, each consisting of two threephase 7.5-kW DC-AC converters. Different from the simulations, only one inductor is employed as the transmission line L line2 in Fig. 6 due to available hardware configurations.
Similar to the simulations, the eigenvalue loci of the experimental parameters are plotted in Fig. 9. The control parameters are modified to ensure the performance in the hardware setup. As evident from Fig. 9, there are RHP poles in both cases, but when the dq decoupling is enabled, the poles are closed to the imaginary axis, so the system can be critically stable. In contrast, the system is much more unstable when the dq decoupling is disabled, or more likely to diverge under large-signal disturbances.

1) CASE I: STEADY-STATE PERFORMANCE
Experimental results in steady-state case are shown in Fig. 10, where the load is 57.5 . The voltage and current waveforms are sinusoidal in normal operation, but when the dq decoupling is disabled, harmonics appear with the frequency of around 550 Hz, indicating the instability of the system.

2) CASE II: DYNAMIC PERFORMANCE
Further results are presented in Fig. 11. In this case, the load is changed from 0.6 kW (57.5 ) to 0.9 kW (another 115 in parallel). The system can be critically stable when the operation point is carefully selected. However, it will lose its stability in response to large-signal disturbances, when the coupling between dq components will boost the disturbances.
Therefore, it can be concluded that the selected modeling methods mapped by the quantitative metrics is eligible in this case for stability validation. Generally, the mapping of modeling methods could be well employed in stability validations especially for MGs with a large number of converters or high penetration of RES, where the difference between methods will be even more significant.

IV. CONCLUSION
In this paper, the modeling methods for stability validation of microgrids are mapped in terms of the properties of microgrids and the focus of study in practice. The modeling methods may differ in complexity and accuracy regarding different applications. Quantitative metrics for evaluating the complexity of modeling methods are proposed, which is obtained based on the order of system, and the metrics have been proved to be feasible and able to clearly reflect the complexity difference between the methods. A study case with simulation and experimental results is presented, demonstrating the validation of stability based on the proposed quantitative mapping policies.
The mapping of modeling methods can serve as a practical step in the stability validation of microgrid systems, and can effectively formalize the validation of stability in microgrids. However, the mapping in this paper is aimed at optimize the validations with the particular objective of complexity reduction, which could vary and accordingly determine the definition of the metrics. Besides, the topologies of converters in microgrids like multi-level converters remain unconsidered in this paper, which could be a future extension of this topic.