Modeling Techniques and Stability Analysis Tools for Grid-Connected Converters

Large-scale integration of renewable generation interfaced to the network through power electronic converters has led to drastic changes in power system dynamics. In islanded microgrids or weak grids, different control concepts for the synchronization of converters have been proposed to provide virtual inertia and improve their resilience against transient events, ensuring safe operation without heavy redundant design. The complexity of these power-related control algorithms and their interaction with the inner control loops causes problems in frequency components above the range of traditional studies which calls on modeling techniques with a wider bandwidth. This work aims to provide an outline of modeling methods for grid-connected converter dynamics from subsynchronous to switching sideband frequency range and relevant analyzing tools. The major contributions of this work are: 1. Theoretical foundations and the derivation processes are discussed for each of the modeling methods within the time domain, frequency domain and harmonic domain. 2. Similarities and differences between these methods are highlighted and recommendations are given regarding different grid situations. 3. A case study with an active front end converter is shown and the analysis results are validated by simulation and Hardware-in-the-Loop (HiL) test, illustrating the effectiveness of these methods.


I. INTRODUCTION
Conversion of electrical power from renewable power sources nowadays are accomplished mostly by power electronic systems [1].This trend introduced faster dynamics and a larger amount of harmonics.In particular, converters with high switching frequency and multiple nonlinear components inherently generate harmonics and could interact with each other, which are pretty common within electrified transportation systems.The case in Swiss locomotives shows how the coupling of harmonics within different frequency ranges could damage the system severely [2].As for the aircraft industry, the recent trend is to substitute hydraulic and pneumatic subsystems on conventional aircraft with their electric counterparts, namely More Electric Aircraft (MEA) [3].The varying frequency and multiple kinds of on-board electrical energy generation/utilization devices make it even more complex than locomotive systems.This has been motivating researchers for discovering an approach capable of characterizing the frequency coupling in both transient response and steady-state operation.Stability analyzing methods are also under great concern so as to couple with the ever-developing modeling methods.
Small-signal stability is one of the most important factors when designing the controller of a grid-connected converter.State-space-based time-domain (TD) modeling and eigenvalue analysis is a way to determine the oscillation mode and the parameters which could trigger resonances [4].However, building the state-space model requires knowledge of the structure and parameters of the system.Hence, when modeling a power electronics-based system, a change in the number of converters or the control loop may lead to a new state-space matrix.This problem is solved by the Component Connection Method (CCM) to some extent [5], yet the exact parameters of the whole system are still necessary.The focus of the frequency domain (FD) impedance modeling instead lies in the extrinsic properties of the converter output port and the power network input port [6].Such a model can either be built mathematically or measured in the real physical system [7]; hence grey-box systems can be analyzed accordingly.Further, according to [8], the linearized state-space transfer function is equivalent to the impedance model when the voltage state variable is considered as the input and the current state variable as the output.
However, in renewable power plants and microgrids, converters working in parallel could be using different switching frequencies, modulation techniques and topologies with bidirectional operation [2].Also, coupling between power-related controllers, inner voltage/current loops, and PLLs becomes much more severer.This tends to bring unexpected harmonics at frequencies ranging from hundreds of Hz to kHz due to the frequency coupling from the multiplication of several timevarying signals and bidirectional interaction [9].In this case, the nonlinear time-variant characteristic of converters are so obvious and the coupling between DC and AC components are so severe that they can no longer be approximated as a time-invariant system.Modeling approaches such as multifrequency averaging [10] and harmonic linearization [7] could solve this problem to some extent by introducing frequency shifting.However, the describing function and harmonic balance approaches enforce the single input to a single output map.Thus, the resultant model is actually part of the linear time-periodic system.To better analyze the effect of switching modulation and resultant frequency coupling, harmonicdomain (HD) models of converters need to be built [11].
The traditional phasor modeling approach was built on the assumption of "Sinusoidal-Steady-State" or "Quasi-Steady-State" whereby the dynamic process is relatively slow compared to the fundamental frequency.Therefore, it cannot predict precisely the fast dynamic response of switched power electronics-based systems.To address this issue, a method based on the generalized averaging method (GAM) called the dynamic phasor modeling (DPM) was proposed by S.R. Sanders and used to model a series resonant converter [12].Dynamic phasor (DP)-based methods are actually a timevariant Fourier decomposition with a sliding time window set according to the critical fundamental working period.Thus, it has a bandwidth close to twice the fundamental frequency and could predict the synchronous oscillation precisely with a significant reduction in calculation time compared to the electromagnetic transient (EMT) simulation, as the complex calculation of the switching process is avoided and the simulation step size can be increased [13].
Besides DPM, the Harmonic State-Space (HSS) modeling is also an effective method for modeling systems with multifrequency coupling.The HSS method was first developed for the analysis of helicopter blade dynamics [14], and was then utilized for motor drive systems in order to capture the frequency coupling of nonlinear and time-varying components [2].The DPM and HSS modeling methods are based on the same initial hypotheses and they both utilize the method of frequency segmentation to make the system time-invariant.However, HSS takes the form of the exponential Fourier series while the DPM appears to be trigonometric, and thus leading to different mathematical representations of electrical states.More importantly, for small-signal analysis, the HSS model was derived directly from the time series Linear Time-Periodic (LTP) state-space equation.In contrast, the DPs are derived by integrating the system variables over a moving time window [4].This leads to different properties and a different range of applications.
Stability analysis also develops with the ever-improving modeling methods.Basically, stability assessment can be categorized into signal frequency response-based [6] and energy-based [15].Energy-based methods have gained attention in renewable integration to the power grid as they give insight into the system from the perspective of power flow.The whole power system is always stable as long as each component stays in the passivity region [16].In contrast, a system oscillation mode is obtained from signal-based methods so that specific harmonics damping techniques are utilized to satisfy the grid code requirements.Moreover, with HD modeling methods, time-and frequency-domain analysis can be used in one unified model, and therefore more dedicated controllers can be designed accordingly [17], [18].Besides these numerical analysis methods, transient stability analysis with a large signal non-linear model in phase plane provides a visualized way to analyze the stability and performance under larger disturbance [19]- [21].
In this work, the derivation process and theoretical background of the above-mentioned modeling are discussed.Then, the relevant stability analysis methods are derived.Some representative work among recent research works are reviewed and a comparison is given on their properties and applicable scenarios.A case study on an active front end converter with these methods is presented in order to validate their accuracy and provide an example of the application scenario.Analysis results are validated by both PLECS simulation and HiL test.

II. MODELING METHODS
Traditionally, there are two mainstream modeling routes, namely, time-domain modeling and frequency-domain modeling.Each methodology has its own strengths and is interchangeable due to the same physical theorem and circuit topology behind it.Recent studies have made the boundary of time and FD even vaguer with the combination of TD state equations and Fourier decomposition [17].In this section, three modeling methods are discussed regarding their derivation process and properties.A comparative review is given based on the application range.

A. TRANSFER FUNCTION BASED SMALL-SIGNAL MODEL
The fundamental concept of equivalent impedance modeling is to linearize the system around a stable operation point with transfer functions in the FD and represent each subsystem of a power system in the FD by an impedance connected to a voltage source (Thévenin equivalent).Based on this concept, any kind of power distribution network can be considered as a power source and load, as depicted in Fig. 1, where: The impedance of each part of the equivalent subsystem can be either calculated numerically by a transfer function, simulation, or measured in the real system [22].The advantages and basic properties of impedance-based modeling for AC and DC distribution networks have been reviewed in detail [23]- [25].Therefore, this paper mainly discusses recent progress and engineering practice, focusing on 3-phase AC distribution systems and other emerging application scenarios.Regarding different linearization coordinate systems, impedance models can be divided into three major types, namely: dq-domain impedance, sequence-domain impedance and phasor-domain impedance model.

1) dq-DOMAIN IMPEDANCE MODEL
The dq-domain impedance-based stability analysis applied with GNC was firstly carried out by Belkhayat in 1997 [26].Impedance models under dq-reference frame have the advantage of simple structure and effectiveness for modeling the electrical equipment with dq-reference controllers.Taking the grid-tied operation as an example, the major factor that could trigger synchronous instability lies in the negative-valued impedance introduced by the phase-locked loop (PLL) [27], the transfer function of which can be expressed as: with: Where V s d is the steady state operation point of d-axes voltage, k p,pll and k i,pll are the proportional and integeral parameter of PLL.Then the impedance model of a 3-phase VSC with current feedback controller can be derived: with: Where G i (s) and G ff (s) are the transfer function of the current controller and the feed forward compensator, L f is the inductance of the filter inductor on the inverter side and I d0 , I q0 denotes the d and q-axis current at equilibrium.It is also revealed that the impedance between I q and V q dominates the synchronization frequency due to the negative impedance within the PLL control bandwidth.The simulation results are compared with experimental results in [28], which offers a new approach to testify the accuracy of the AC impedance model regarding stability analysis.In [20], it is further testified that the synchronous reference frame (SRF)-PLL is threephase balanced under asymmetrical faults when a pre-filter is used.Thus dq-based model of the major components still holds in this case.The model has been further developed to be able to deal with a larger phase perturbation (up to 38 o ) by utilizing a second order linearized model of the PLL [29].Moreover, converters with Virtual Synchronized Generator (VSG)-based controller [30] and doubly-fed induction generator (DFIG)-based wind power generation systems [31] could also be modeled accordingly.

2) SEQUENCE-DOMAIN IMPEDANCE MODEL
Based on the harmonic linearization principle, the system impedance can be divided into positive and negative sequence components (Z p and Z n ) [7].TD perturbation can be written as: where V 1 is the fundamental voltage amplitude with frequency f 1 , V p , V n , φ vp and φ vn is the magnitude and phase of the sequence perturbation at frequency f p and f n , respectively.The sequence-domain impedance per phase can be expressed as: ϕ is the angle of the current vector I 1 , and f 1 denotes the fundamental frequency.When the PLL is considered, a negative-valued impedance is introduced, which will cause a low frequency oscillation when coupled with an inductive weak grid.More specifically, due to the asymmetric nature of the dq-axis control, a positive sequence disturbance at f p will trigger responses at both f p and (2 f 1 − f p ), which is defined as mirror frequency decoupled (MFD) [32].Such phenomena is also elaborated in [33] within the stationaryframe, which also reveals the influence of the steady-state voltage initial phase on the frequency coupling terms.Compared to the dq-domain model, the mathematical derivation of the sequence-domain impedance model is more complex.However, the impedance is not affected by reference frame or system asymmetries.By taking the entire power system as a closed-loop system, as depicted in Fig. 2, a single-input single-output (SISO) impedance model of each sequence is developed for 3-phase grid-tied VSCs in [34].Besides, each sequence is independent with real physical meaning and can be easily measured by commercial equipment.Based on the properties above, it is then used for modeling of MMC [35], DFIG power generation systems [36], and etc.In [37], the process of modeling a 3-phase VSC into self immittances and transfer immittances is discussed, which gives an insight of the coupling of different frequency between AC and DC ports.

3) PHASOR-DOMAIN IMPEDANCE MODEL
Phasors have been widely used in power system analysis so that magnitudes and angles of electrical states can be represented with a single variable.Considering the magnitude and angle of perturbations at the fundamental-frequency, phasors can be linearized as: where θ 0 and ϕ 0 denotes the angle of voltage U and current I at equilibrium, the subscripts g denotes the SRF of the power grid and P denotes the phasor domain transfer matrix of the output current.Then the terminal impedance model in the phasor domain can be described as [38]: where: It can be seen that the admittance matrix in this case takes the special form of: . Such matrices can be transferred to a symmetric similar matrix as: where: 1 j 1 -j .Then, according to the principles of the two-port network, the reciprocal two-port circuit of a symmetric admittance matrix can be built as shown in Fig. 3. Thus far, the multi-input multi-output (MIMO) problem introduced by the two-port equivalent circuit network can be transferred to a SISO problem.Different from two separate SISO paths in [34], all factors are compacted in one system so that the conventional FD analysis methods can be performed.Also, the output current limitation according to the gain margin could be calculated.
According to the discussion above, regardless of the reference coordinates being used, the major challenges of impedance modeling lie in the frequency coupling effect in AC/DC combined network, which still calls on further exploration when applied to large-scale, power electronic dominated network modeling.The interchangeability between phasor domain admittance and dq-domain, sequence-domain counterparts are discussed in [39].Also, a three-by-three transfer matrix is built so as to capture the coupling effect between AC and DC sides within the low-frequency range by separating the AC and DC component disturbance.Moreover, the impedance model of converters with multiple inner coupling effects is being extended to the switching frequency range with the HSS model [40], which will be further discussed in section II-C.The comparison of properties of small-signal impedance modeling methods is summarized as shown in Table 1.

B. DYNAMIC PHASOR MODELING
Different from the transfer function dominant methods mentioned above, the GAM to obtain the DPM is based on the property that any TD waveform x(t ) can be represented inside the interval τ ∈(τ − T, τ ] by the following Fourier series [41]: with ω s = 2π f s and X k (t ) being the complex Fourier coefficients, also named phasors.The dynamic k th -order phasor at time t, X k (t ), can be expressed as: Where x k (t ) denotes the average k th -order phasor over the period T .The major advantages of the DPM over traditional methods come from the three important properties.

1) CONJUGATIVE PROPERTY
Phasors can be expressed in the form of complex value, thus, except for zeroth-order DPs, a single DP is represented by real and imaginary variables.By separating the two parts of the dynamic equations, a nonlinear real-valued DP model has the property of: 2) DIFFERENTIAL PROPERTY The derivative of k th -order phasor can be expressed as: As the state-based dynamic models can generally be expressed as d dt x(t ) = f {x(t ), u(t )}, where u(t ) is the time-series input signal.The DP derivative can be expressed as: substituted with the electrical states of the basic circuit components, the state equation of which can be derived as shown in Table 2.

3) PRODUCT PROPERTY
For the two time series waveform x 1 (t ) and x 2 (t ), the product of their phasors is: That is, the time series waveform can be transformed to the convolution of their counterparts in phasor domain.According to (18), the voltage waveform at the switch terminals of a pulse-width modulated (PWM) system can be expressed as: where v s and q denote the voltage waveform of the power source and PWM switching function, respectively.Therefore, a link between the source and load component can be established accordingly.
As discussed above, the essence of DPM is to retain the dominant Fourier coefficients term to describe the major characteristics of the system.For instance, when analyzing high-frequency PWM switching circuits, the DC component should be addressed to capture the low-frequency behavior.Then the result would be the same as the state-space averaged model.For resonant converters with some predominantly fundamental frequency, first-order and the dominant harmonic DPs should be retained to capture the major sinusoidal response and the zeroth-order coefficients of those states that exhibit slowly varying behavior [12].As the DPM is a large signal model that varies relatively slowly over time, a larger simulation step size can be used compared to EMT models.In [42], the stability criterion for interfacing between DPM and EMT is given.However, as the amount of DPs being considered increases, the calculation complexity of the model raises significantly, which defeats the purpose of saving computing resources.The multi-frequency DPM is for small signal analysis when the system is linearized to form a Linear Time-Invariant (LTI) system.Also, it is worth mentioning that when the initial phase is extracted from the DP model and expressed in the form of a generalized dq (GDQ) transformation, a GDQ model is then derived, which puts an emphasis on the phase difference so that the power transmission is more convenient to calculate accordingly [43].Combined with the concept of sequence component, unbalance faults analysis can be performed [44].
It follows that the DPM could be used for both large-signal transient simulation [45]- [49] and steady-state stability assessment [50]- [52].In both cases, it can be coupled with other modeling methods by an appropriate interface.The tradeoff between model precision and system complexity can be achieved by proper selection of the phasor sequence.

C. HARMONIC STATE-SPACE
Similar to the DPM, the HSS is based on the principle that TD variables can be transferred to the FD by utilizing (13).Thus, when a small signal disturbance e st appears, the state matrix of the system becomes [53]: the differentiation of (20) could be expressed as: for an arbitrary LTP model: the HSS variable update equation can be derived as (23).
The coefficients of same frequency on each side of the equation should be equal, so the exponents can be eliminated, then the kth-order HSS equation is obtained: Consequently, the steady state equilibrium can be derived as: where: where [•] denotes the Toeplitz matrix or diagonal-constant matrix which is used for convolution calculation.According to (C), the HSS method provides a unified way to construct a multi-frequency model which is linearized around operating trajectories in the harmonic basis of the period T s .
In fact, harmonic linearization can be regarded as a combination of FD stabilization and perturbation linearization [25].Thus, power electronic converters of any topology, which can be expressed by a state-space model, can be easily transferred to an HSS model: single-phase converters [54], transformers [55], three-phase converters under both balanced and unbalanced situation [43], [56] and especially modular multilevel converters that have the nature of nonlinear inter-resonance between sub-modules [40], [57].Similar to the DPM, the selection of harmonics under consideration has been a critical concern for HSS regarding computational burden.The specific simulation time comparison of different orders of harmonics of a single-phase transformer has been made in [58], indicating that the computation duration of the model increases quadratically with the harmonic order.An effective solution to this problem is to select the critical harmonics according to the impedance/admittance map [59].By considering only the harmonic components with steadystate impedance larger than a certain threshold or significantly  larger than its adjacency, the order of the HSS model can be effectively reduced without causing an unacceptable error.
In [60], the HSS model of a single-phase Second-Order Generalized Integrator (SOGI)-PLL is built, improving the accuracy of the controller behavior compared to LTI models.Further, the nonlinear PLL dynamic is considered in statespace modeling in [61], which leads to a nonlinear modular state-space model.Time-domain simulations show comparable accuracy with the detailed switching model in both transient and steady-state.Though this method is still based on traditional state-space thus high-frequency ripples caused by switching dynamics are not present, the way of processing the low-frequency nonlinearity is also constructive to HSS modeling.
By reason of the foregoing, the derivation process of the two methods follows the same steps with different order and variable representation [43], thus leading to real states with the DPM while complex states with the HSS model [62].The flow charts of the derivation processes are shown in Fig. 4. The linearization process of DPM is encircled by a dashed line as it is not performed when modeling transient behavior.
Similarly, if the state-space model used for the HSS is not linearized for small-signal analysis, the consequent HD model is also accurate with the dynamic behavior, which is known as the dynamic harmonic (DHD) model [11], [63].Due to the difference in linearization principle, the HSS model is linearized around the time-periodic trajectories, whereas the DP model is based on linearization around the time-invariant points.Consequently, as the HSS model is the Fourier decomposition of an existing state-space model without changing the model topology, it is intrinsically more programmable and could be executed automatically by computer programming.In general, the DPM has the advantage that the TD simulation is either independently or interfacing with power system EMT models.However, the HSS is relatively less complex for impedance-based stability analysis and a larger amount of harmonics can be taken into account [64].

III. STABILITY ASSESSMENT METHODS
The above-discussed modeling methods serve as the foundation to evaluate the system dynamics and stability against different kinds of disturbance.The corresponding stability assessment methods will be discussed in this section.

A. EQUIVALENT IMPEDANCE-BASED ANALYSIS 1) IMPEDANCE RATIO CRITERIA
The equivalent impedance ratio-based analysis was introduced by Middlebrook to analyze DC-DC systems [6], which provides a sufficient criterion for the small-signal stability of SISO systems.Due to its over-conservative nature, modified criteria are proposed on its basis.Five of the most used methods ranked from high to low according to conservatism are: Middlebrook criterion, Gain Margin Phase Margin (GMPM) criterion [65], Opposing Argument criterion [66], Energy Source Analysis Consortium(ESAC) criterion [67] and perturbation based criterion [68].
However, the criterion used most for combined AC, DC networks is the Nyquist criterion in which the loop gains (Z S /Z L ) of the system described in (1) should satisfy the Nyquist stability while the converter system being studied is internally stable.Moreover, when analyzing a system with multiple converters, this stability assessment should be performed on each one with others considered as loads, so that all could satisfy the stability criterion.In [38], due to the special structured reciprocal SISO circuit model built in the polar-frame, MIMO inverter systems could also be analyzed with the Nyquist criterion.
In addition, the Generalized Nyquist Criterion (GNC) can be adopted in order to access the stability of the multi-port impedance system.Namely, the total numbers of clockwise and counterclockwise encirclements of the eigenvalues of the impedance ratio around the complex plane critical point (−1, j0) should be equal to the number of its open-loop righthalf-plane (RHP) poles.Consequently, for minimum phase systems, the stability is guaranteed as long as the eigenloci don't encircle (−1, j0).However, it is pointed out that for weak grid where non-minimum phase systems are becoming dominant, Bode plots of eigenvalues tend to fail for system stability analysis [69].This problem was addressed in [70] by identifying the open-loop RHP poles through magnitude slope and phase change rate variation.Moreover, a designoriented impedance specification procedure utilizing a Bode plot of individual impedances with a special focus on crossing direction is proposed to stabilize interconnect systems.Note that the determination of the crossing number at 0 Hz is very important in this method, which is discussed in [71].
Inspired by the participation factor (PF) analysis in the state-space model, the eigenvectors of the MIMO impedance ratio (or return ratio) matrix are used to identify the PFs of each converter in a multiconverter system when instability appears [72].The key converters that trigger the instability can be found accordingly.The installation locations of converters can also be chosen to ensure better stability margins.In [73], the residue of each node admittance is instead utilized to define the PF.Based on the chain rule, a three-layer gray-box method is developed to estimate the critical node and the exact passive component or control loop of it that affects the stability of the whole system in a perturbation and observation (P&O) manner.Further, the PF analysis is extended from the oscillation frequency points to critical frequency ranges where the eignloci are close to or tend to encircle (−1, j0) [74].Active damping techniques can be adapted accordingly to enlarge the PM and GM of the critical converters within critical frequency ranges.
The norm-based stability criterion has been proposed to simplify analyses and indicate the stability margin levels [75]- [77].Nevertheless, the resultant stability criteria are actually sufficient and unnecessary conditions because the system can be stable even if its stability requirements are not met.To further reduce the conservatism, the relevant stability margin was regulated according to Gerschgorin's Theorem when the infinity-norm is used [77].

2) PASSIVITY-BASED CRITERIA
Passivity is mathematically defined when a one port network system only absorbs energy from the outside physical world, namely: The stability of a passive system is inherently guaranteed because the energy dissipation process always drives the system back to its equilibrium point [78].Consequently, the passivity criterion of the one port system can be defined as [79]: r Z (s) does not have right half plane (RHP) poles and the entire Nyquist curve is within the RHP r Z (s) has a non-negative real part for every ω Compared with impedance-ratio-based methods, passivitybased methods are more suitable for multi-converter systems.When the impedance of the converters within the system all satisfy the passivity margin criterion, their summation will also be in between the passivity region with an even larger passivity margin.For example, as depicted in Fig. 5, a passivity margin angle of α is required for the system to make sure possible resonance is well damped.Z −1 inv1 and Z −1 inv2 are the output admittance of two converters in parallel.It can be seen that the parallel admittance Z −1 inv = Z −1 inv1 + Z −1 inv2 still lies within the expected passivity region.This provides a decentralized way of identifying the overall system stability [16].Moreover, the source and load concept is no longer necessary to be defined.Consequently, with the purpose of increasing the overall system passivity [80], the control structure could be simplified while maintaining a high control bandwidth.However, as proven in [16], the passivity margin criterion is a sufficient and unnecessary condition.Thus, a trade-off has to be made between the stability objective and the conservation of the controller design.

B. STATE-SPACE EIGENVALUE ANALYSIS
The eigenvalue analysis has been a widely used method for accessing the stability and oscillation mode of the traditional power systems.The state matrix A of a small-signal linearized LTI model can be extracted to form the characteristic equation: the resultant eigenvalues could be an indication of the dynamic modes of the power system.And the participation rate calculated by the corresponding left and right eigenvectors can also reveal the impact of the grid-connected converters on the instability [4].However, the PFs of the state-space model are all scalars instead of transfer functions with frequencydomain impedance ratio matrices.Thus, they can only predict the PFs of each node under oscillation modes instead of port dynamic over a frequency range [74].That is, the PF analysis with state-space models offers limited guidance in controller design to optimize the oscillation damping.This problem is solved by nonlinear modular state-space modeling of power-electronics systems [61], in which the VSCs models are based on reduced-order impedance model and nonlinear PLL dynamic is considered.By tuning different parameters, the nonlinear model can be linearized at different equilibrium points with the resultant Jacobian matrix.The eigenvalues distribution trend gives information on the system damping ability of different oscillation modes in response to the parameter change.
With the increase of control complexity and power electronic converter penetration rate, the order of the system becomes larger.Thus, the difficulty of building and deriving the system model increases significantly.This problem was solved by the CCM, which defines each critical part as a subsystem and connects them according to their linear algebraic relation.Its modular structure also improves the extensibility of the model [5].

C. FLOQUET-THEORY-BASED METHOD
The Floquet theory illustrates the relationship between LTP and LTI systems [81].By introducing the time-periodic matrix P(t, t 0 ), any LTP system can be transformed into an LTI model.According to [82], the solution of an LTP system (22) can be expressed as: where the matrix (t, t 0 ) is a state transition matrix and is a periodic time-varying matrix with period T s .If P(t, t 0 ) can be numerical obtained, the original time-varying state transition matrix (STM) A(t ) can be transferred into a time-invariant matrix Q by: The eigenvalues of Q are called Floquet exponents (FEs) and can be used to illustrate the system dynamics [83].If P(t, t 0 ) can not be numerical obtained (which is true in most cases), the stability of an LTP system is instead analyzed mainly by (T A , 0), where T A is the minimum period [84].
The eigenvalues of (T A , 0) are called Floquet multipliers (FMs).In [85], a complete model is derived for an asymmetric single-phase cascaded H-bridge inverter.Different steadystate operating trajectories are considered, loci and moduli of resultant FMs are used to evaluate the system stability.Compared to FEs, FMs can only indicate the system stability instead of characterizing the oscillation frequency in the fundamental strip [14].The relationships between different analysis methods are summarized in [86], where an LTP-MMC model is built and a numerical solution of the constant matrix Q in the TD is used to calculate the FEs and analyze the interaction between different controllers under unbalanced grid conditions.Further, it is proved that FEs provide a physically sound way for truncation harmonic number selection of HSS models for eigenvalue analysis [87].

D. PHASE PORTRAIT
Unlike the above-mentioned methods, which access the smallsignal performance with linearized models, phase portrait methods focus on the large-signal stability (transient stability).Namely, the ability of grid-connected converters to maintain synchronism when severe grid disturbances occur or when the fundamental frequency varies [88].Simultaneously, a large signal model needs to be built in order to capture the dynamics under serious phase jumps.Taking the model of the SRF-PLL as an example, different from the small-signal model as shown in (2), which is based on the assumption that the phase angle difference is sufficiently small, the largesignal model is expressed as: Where θ and ω denotes the difference between the real value of the phase angle and the frequency and their estimation, V is the amplitude of the grid voltage.Based on (33), the global stability of the SRF-PLL was studied in [19].
The resultant phase portrait is shown as Fig. 6, which separates the operation range into multiple small zones around  stable equilibrium points.Disturbances within each small zone will converge to the unique stable point.By that, the phase portrait offers a visualized way to access the dynamic performance with both frequency and phase jump, while oscillation can also be confirmed when under specific disturbance.In [21], based on the DPM, the time-domain simulation and phase plane trajectory of a grid-forming converter is given for different control strategies.Combined with the impedance model and eigenvalue analysis, the coupling effect between the power-related control and inner closed-loop control is further revealed.

IV. DISCUSSION
Thus far, the above discussed modeling and stability analysis methods can be categorized based on the modeling domain as shown in Fig. 7, and the corresponding application range [89] is depicted in Fig. 8.A comparison of some main features of the modeling approaches is given in table.3.For stiff grids with large inertia provided by synchronous generators, PWM modulation of converters usually won't affect the whole TABLE 3. Comparison of the Modeling Approaches system stability.Also, power-related controller or synchronization unit and inner controller can be considered separately because the coupling is neglectable.So the common approach is to analyze the large-signal synchronous stability (frequency dynamics) with TD EMT equations [20], [21] and simplified symmetrical small-signal models for aiding the inner controller design [90].
However, in weak grids like renewable power plants and microgrids, the large impedance and frequency multiplication of different converters enhance the coupling between different control loops and synchronization units.In this case, different control loops should be considered comprehensively when building the model.In particular, the transfer function models can be used to identify the frequency response of the terminal currents and voltages usually caused by the PLL, the inner voltage and current control loop and passive filter components (voltage dynamics).Thus, the bode plots of different ports can be used as guidance for controller design.Also, both impedance ratio and passivity-based methods allow black-box operation, so a detailed system configuration is not necessary.The major difference between these two methods is that passivity-based methods are decentralized, which means the system is stable as long as all units are passive.However, the passivity-based criteria are over-conservative and thus may lead to poor dynamic performance.Impedance-based methods instead could identify the critical node that may cause instability and give constrictive information on parameter tuning [73], [74].
Time-domain state-space model instead can also provide insights into the coupling path based on the participation factors [21].Note that the TD state-space model is actually a series of differential equations.Thus, it can be both linear [48] and nonlinear [61] based on different modeling assumptions.The nonlinear modular state-space models can also be linearized around different operating points, the elements of the resultant Jacobian matrix can also be transformed to transfer function and FD methods can be used.Also, the transient stability of the synchronization units or power-related controllers can be analyzed with the help of the nonlinear model in the phase plane.However, a major drawback of TD modeling is that detailed data is needed which is not feasible in a largescale interconnected network.
The harmonic modeling is a combination of TD and FD.Thus, it can be analyzed with some of the existing methods in both tracks.More specifically, it can be transferred to TD with

TABLE 4. Model Parameters of the Grid-Connected Converter
Fourier inversion and FD with Harmonic Transfer Function (HTF).Besides, Floquet theory for analyzing LTP systems can also be used for stability analysis and oscillation mode identification [85], [87].Compared to classic TD and FD methods, it gives an insight into the switching behavior and harmonic coupling.The frequency of interest is then expended to switching dynamics which is extremely important in power electronics dominate microgrids with different and a synchronized switching frequencies, making it a solid supplement to traditional modeling.However, the above-mentioned methods for stability assessment still need transform or indirect execution.In [17], based on Harmonic Riccati Equation, globally and asymptotically stabilizing periodic feedback control laws are proposed which directly come from HSS models without invoking the Floquet theory.

V. CASE STUDY
In this section, the above-mentioned modeling and stability assessment methods will be utilized to analyze a two level 3phase grid-connected converter with proportional-integral (PI) controller for the DC-link voltage control and the dq-domain current controller as shown in Fig. 9. Theoretical analysis will be compared against the simulation results.The parameter is shown in Table 4 .The design reference of the controller can be found in [91].

A. TRANSFER FUNCTION-BASED SMALL-SIGNAL ANALYSIS
The αβ-domain model is widely used in the modeling of VSCs [33] and has been proved to be interchangeable with dq-domain model [92].In this paper, dq-domain model is carried out as the controller is designed within SRF.The dqdomain small-signal model of the converter system is shown in Fig. 10.G dc and G i are the PI controller of the DC-link voltage and dq-domain current, with G dei denotes the decoupling term.G del and G PW M are the first-order delay and PWM drive coefficient, respectively, the subscript of electrical state s denotes the real signal with the physical system and c denotes the signal detected by the controller.G pv is the linearized transfer function between the active power disturbance and the DC-link voltage response [92], which can be expressed as: The dashed line parts illustrate the small-signal disturbance introduced by the PLL, where the small-signal transfer function is shown in (2).The d denotes the small-signal disturbance, and subscript 0 denotes the steady-state value.As the system is a SISO system but has two state variables in d and q-channel respectively, expanding matrices are used as depicted in the figure in order to generate the small-signal variable on q-channel.The derivation process of each block and the impedance measurement algorithm has expatiated in [28], [92].
The bode plot of the derived admittance by this model is shown in Fig. 11 and compared with the measurement results obtained by a simulation conducted with the software Plexim PLECS.It can be observed that the d-d channel admittance Y dd of the converter system has a positive increment in the low-frequency range, which is known as the constant power load (CPL) behavior.Such phenomena are introduced by the DC-link voltage control loop and could cause lowfrequency oscillation and voltage overshoot.In contrast, the low-frequency phase jump in q-axis is mainly caused by the PLL as the input current reference from G v is always 0. Note that the disturbance signal within the bandwidth of the PLL may trigger a nonlinear response in the error signal output.Fig. 12 illustrates the variation in the PLL frequency error signals, with the disturbance of the same amplitude and waveform, when the bandwidth is 3 Hz and 10 Hz.Thus, the measured impedance in the coupling admittance may have a relatively larger bias.
Combined with the LC-type impedance of the transmission line, GNC-based analysis can be performed as shown in Fig. 13, where the two eigenvalues of the system are depicted in both Nyquist and Bode plots.The zoom-in plot in Fig. 13(a) shows that the characteristic root locus plots of the system do not encircle (−1, 0) point, the system is thus stable.However, combined with Fig. 13(b) the root locus of λ2 intersects the unit circle at approximately 715 Hz, which could trigger convergent oscillation around this frequency under transient event.

B. TRANSIENT STABILITY ANALYSIS WITH PHASE PORTRAIT
As the bandwidth of the inner current loop is much larger than the bandwidth of the PLL, the dynamic of which is always neglected when studying the effect of the PLL.Considering the grid-PLL interaction in the above-mentioned system and assuming that the impedance is static, a second-order nonlinear model of the dynamic around the fundamental frequency can be derived [20]: Where V amp is the amplitude of the grid voltage, R line and X line are the line resistance and reactance.According to the global stability zone discussed in [19], the phase domain stability zone of the converter system is shown in Fig. 6.An 8 Hz drop in PLL reference frequency is added to the system so as to get the trajectory in the phase plane.Compared with the PLECS simulation shown in Fig. 14, the large-signal model provides an envelope curve of the PLL transient response with 10.5% larger predicted amplitude within the second oscillation.Due to the absence of the inner current loop dynamic, the model fails to predict the high-frequency oscillation which can be seen in the characteristic loci shown in Fig. 13 and the switching component that will be discussed in V-C.In contrast, a more dedicated model considering the inner loop with dynamic phasor is derived in [21].However, the model is so complicated that it cannot be described by a first-or second-order non-linear differential equation.Therefore, the stability analysis is based on a reduced-order time-domain simulation framework rather than a design-oriented graphical analysis with the trajectory in the phase plane.

C. HARMONIC DOMAIN MODELING AND STABILITY ANALYSIS
In order to access the HD dynamics of the converter system, the HSS model, having flexibility in harmonics order, is used to develop the system model in this paper.The open-loop HSS model is first built for the converter, which is also referred to as the topology model [93].Time-domain discrete Forward-Euler iteration of the model is performed so as to illustrate the accuracy in time-domain simulation.The simulation result and comparison with PLECS are shown in Fig. 15.It can be seen that the major harmonic components in both models agree with each other during the steady-state of the open-loop operation.
To include the controller behavior, the model needs to be linearized around the time-periodic trajectories calculated by the open-loop topology model.Assuming the switching instant of the controller output signal is very small, the relation between AC and DC electrical states can be expressed as: Where d denotes small variation of signals.This principle also applies for HSS models.and outputs Y c = [dS a dS b dS c ] T , where d denotes the small-signal disturbance and the physical meaning of the variables are shown in Fig. 9. Assuming that the PLL can track the grid frequency perfectly, ideal Park and Inverse Park transformation are used in the controller model to link the static and synchronized reference frames.Consequently, the HSS model of the converter topology and the controller can be derived as: (37) shown at the bottom of this page. Where: ) .I m and Z m denote the m = (2n + 1)-order identity matrix and zero matrix as the Toeplitz matrix of these two kinds are themselves.n denotes the number of harmonics being considered in the model.
For symmetrical Sinusoidal PWM modulation (SPWM), the Fourier decomposition of the switching signal in upper arm of phase A can be expressed as [94]: (43) shown at the bottom of the next page.Where M r is the modulation rate, ω sw and ω 0 are the switching and fundamental angular frequency, m and n denote the order of the carrier wave and the adjacent harmonic wave, respectively.The expression besselj(•) denotes the Bessel function of the first kind, which is: besselj Other signals can be derived by changing the initial phases of the mathematical expression.The resultant time-series waveform is shown in Fig. 16, the modulation rate is 15 and the harmonics calculated is 255 th .The switching function provides a link between DC and AC components and reproduces the harmonics introduces by the switching behavior.Thus, the more harmonics being calculated, the more it appears like a square wave.
Fig. 17 shows the schematic diagram of the combined topology and controller model.As illustrated by the dashed blue arrow, the state variable of the topology model is the input of the control model, and the outputs of the control model act as inputs of the topology model.This coupling relationship between the two systems requires further regulation of the state matrix to form the compacted linearized system, but the deviation process is not shown in [93].One way to combine the two parts and eliminate the duplicate states is shown below.
In order to derive the compacted model, first consider the variation of the DC-link voltage reference signal as controller input set U c1 , the DC-link voltage and converter output current disturbance as controller input set U c2 , the disturbance in the PCC voltage as topology input set U p1 and the switching signal as topology input set U p2 and grid side current disturbance.Separate the input matrix of the two systems accordingly and merge the state matrices of the two systems based on the relation in Fig. 17, the close-loop HSS state iteration can be derived as: Further, to obtain the I/O relation between different electrical states, set the combined states as dX a , voltage disturbance at the point of common connection (PCC) and DC-link voltage reference as input dU a , the state variables of the topology as the outputs dY a , and the resultant state matrices as A a , B a , C a and D a .The HTF can be established as [95]: By setting s as zero and considering up to the 35th order harmonic components, a steady-state linear mapping between the Fourier coefficients of the converter electrical states and voltage disturbance can be built as shown in Fig. C. The figure can be interpreted as HD coupling admittance since the converter input current is seen as the output of the HTF.It can be seen that voltage disturbance at a fundamental frequency can trigger a current disturbance at both fundamental frequency and around the switching frequency and vice versa.
Besides, the resonant zone can also be predicted as encircled in red, which is consistent with the major components of the THD shown in Fig. 14(b).

D. HARDWARE-IN-THE-LOOP VALIDATION
In order to validate previous analysis, HiL tests with the aid of PLECS RT-Box have been done.The setup of HiL test is shown in Fig. 19.The system parameters are the same as shown in Table 4. Similar to the simulation, an 8 Hz frequency drop is added to the PLL reference frequency is triggered at 0.45 s.The resultant converter output voltage and Fast Fourier transform results are shown in Fig. 20 and the PLL error signal is shown in Fig. 21.The medium frequency at around 700 Hz can be observed in Fig. 20(b), which agrees with the result of the GNC analysis.The PLL frequency error signal conforms to the phase plane trajectory, the dynamic is better damped than simulation which is caused by sampling delay in real-time system.Also, the magnitudes harmonic components are consistent with the major components shown in the HSS admittance map.

VI. CONCLUSION
In this paper, three mainstream modeling methods and relevant analysis tools are discussed for the study of gridconnected converters in the following aspects: derivation process, key feature, and application scenarios.The small-signal impedance model deals with low-frequency and synchronize stability issues.The TD nonlinear large-signal model can also be used to aid the controller design with the visualized phase plane portrait.Besides, HD modeling methods could extend the range of interest to the switching frequency and the wide-band controller dynamic.More specifically, the DPM could be used for both large-signal TD simulation and stability analysis with the system order properly selected, while HSS is more suitable for steady-state harmonic stability assessment due to its LTP nature and complex-valued state matrices.A case study with a three-phase VSC further validates these methods and gives a paradigm on how they are utilized with PWM converters.This critical evaluation of the modeling and analysis tools gives the engineers theoretical guidance for the design and evaluation of power electronics-dominated power systems.

FIGURE 1 .
FIGURE 1. Thevenein equivalent representation of a subsystem with connected load.

FIGURE 2 .
FIGURE 2. Closed-loop representation of a grid-tied VSC based power system.

FIGURE 3 .
FIGURE 3. Reciprocal SISO circuit in the phasor domain.

FIGURE 4 .
FIGURE 4. Flow chart of the DPM and HSS derivation process.

FIGURE 5 .
FIGURE 5. Schematic diagram of the passivity region for two converters in parallel.

FIGURE 6 .FIGURE 7 .
FIGURE 6. Phase plane global stable region and disturbance trajectory of the SRF-PLL.

FIGURE 8 .
FIGURE 8. Application range of different modeling methods.

FIGURE 9 .
FIGURE 9. Structure diagram of the grid-connected converter and its control.

FIGURE 10 .
FIGURE 10.Small-signal model of the grid-connected converter system.

FIGURE 13 .FIGURE 14 .
FIGURE 13.Eigenvalue of the simulated impedance from the small-signal model.

FIGURE 15 .
FIGURE 15.Simulated A capacitor voltage of the HSS model.
A p dX p + B p1 dU p1 + B p2 dY c = A p dX p + B p1 dU p1 + B p2 C c dX c + D c1 dU c1 + D c2 dX p1 = A p dX p + B p2 D c2 dX p1 + B p2 C c dX c+ B p1 dU p1 + B p2 D c1 dU c1 , (45)d Ẋ c = A c dX c + B c1 dU c1 + B c2 dU c2 = A c dX c + B c1 dU c1 + B c2 dX p1 .(46)Subsequently, fill the empty elements with zero matrices, the close-loop state matrix can be derived as:[B p2 D c2 |Z ] B p2 C c [B c2 |Z ]

FIGURE 20 .
FIGURE 20.Converter output voltage under phase jump.

FIGURE 21 .
FIGURE 21.PLL frequency error signal under frequency drop.

TABLE 1 . Comparison of Impedance Modeling in Different Domains
• stands for the advantages and stands for the drawbacks.

ga dI gb dI gc dV Ca dV Cb dV Cc dI ca dI cb dI cc dV dc
Thus, the linearized close-loop HSS model can be built by considering the small-signal instant of the topology electrical states and the controller.Set the topology state variables as X p = dI T, inputsU p = dV

loop input admittance magnitude map of VSC HSS model. FIGURE 19. HiL test setup.
FIGURE 18. Close-