Recursive SISO Impedance Modeling of Single-Phase Voltage Source Rectifiers

Due to the nonlinear time-periodic nature, precise impedance modeling of power-converter-based single-phase ac systems is complex. Small-signal modeling for the systems with a fixed operating point is easy, but the extension of the same procedure to the systems with time-periodic trajectory is still challenging. In this article, a recursive single-input single-output (SISO) impedance modeling framework for single-phase rectifiers is presented, which can greatly simplify both the modeling procedure and the resulting impedance model. A general linear-time periodic modeling method is proposed by extending the linear-time invariant modeling method, where frequency-coupling effects are modeled by the transfer function vector rather than the transfer function matrix, so the modeling complexity is reduced. Furthermore, based on the idea of mathematical induction, an analytical and accurate recursive SISO impedance model is derived. The resulting recursive SISO impedance model can characterize the frequency-coupling dynamics of arbitrary order with a low computation burden. Experiments validate the accuracy of the recursive SISO impedance modeling method.

Many research efforts have thus been made to deal with VSC-grid interactions. It has been reported that the frequency coupling effect (FCE) widely exists in the VSC systems, whose generation reasons include the phase-locked loop(PLL), asymmetric controllers, and so on [2], [4]. The FCE implies that there is coupling between different frequency components. Different system structures lead to different numbers of frequency coupling components. The frequency coupling components will form additional sideband loops, which poses serious challenges to accurate impedance modeling.
For three-phase balanced VSC systems, the number of frequency couplings is typically bounded by two [4]- [6]. In order to incorporate the coupling dynamics, the three-phase VSC is modeled as a 2 × 2 impedance matrix either in the dq synchronous rotating frame [7] or in the αβ stationary frame [8]. In dq frame, the VSC system is linearized at constant dc operating points, leading to a linear-time invariant (LTI) system. If performing linearization in the stationary frame, the resulting model is linear-time periodic (LTP). Linearizing in different reference frames yields different types of impedance (i.e., dq impedance, phasor impedance, and sequence impedance), and the correlation among different models are unveiled in [9] and [10]. On this basis, many detailed impedance models including different factors are formulated, e.g., constant power load [11] and PLL [12], [13]. Based on the complex vector notation [14] or algebraic calculation [15], the multiple-input multiple-output (MIMO) impedance matrix can be simplified into a single-input single-output (SISO) impedance model. The SISO impedance model not only facilitates the understanding of the consequence of impedance measurement under frequency couplings, but also provides an intuitive insight into the interactions among VSCs and power grids [16], [17].
However, in either three-phase unbalanced or single-phase systems, the time-periodic nature of the system introduces an infinite number of frequency couplings, which greatly complicates the modeling. Besides, the LTI modeling is no longer suitable because of the time-periodic operating trajectory. Based on the above considerations, the LTP framework-based modeling approaches are more feasible for the impedance modeling of such systems.
Generally speaking, there have been two general multifrequency modeling methods capable of dealing with time-periodic systems. The generalized averaging modeling (GAM) was earlier introduced to capture the cross-coupling dynamics in [18] by Sanders et al., which is also known as dynamic phasor modeling [19], [20]. The key of GAM is to decompose the actual time-periodic waveform into the phasors in different frequency domains based on the harmonic balance principle. On the other hand, the harmonic-state space (HSS) modeling method, which was first proposed by Wereley in [21], has been widely applied to address harmonic stability issues by Kwon et al. in [22] and [23] and by Salis et al. in [24] and [25]. Moreover, the counterpart of the HSS in the frequency domain is called as harmonic transfer function (HTF) approach [26], which directly establishes the harmonic impedance model in the frequency domain [2], [27], [28]. The core idea of HSS/HTF lies in the introduction of the exponential modulated periodic signal, which maps the LTP model as an infinite-dimensional LTI model. The models obtained from GAM or HSS/HTF are characterized by MIMO with infinite dimensions, bringing about significant complexity to stability analysis and physical measurement.
In order to simplify the MIMO model, many efforts have been devoted to SISO equivalent impedance modeling method. Harmonic linearization, proposed by Sun et al. in [29], gives the describing-function-based model. However, the neglect of cross-coupling dynamics limits its accuracy, especially in the low-frequency regions [30]. The harmonic signal-flow graph method is proposed to streamline and visualize the modeling process of harmonic linearization by Shah et al. in [31] and [32]. The additional loops closed through the grid impedance are not been considered. In [33], Qian et al. derive a singlefrequency impedance model, incorporating the effect of the other two frequency-coupling components. This incorporation does improve the accuracy of the model; it is still inadequate to predict the impact of certain factors, such as the modulation method. The latest work in [34], presented by Zhang et al., extends the model to include coupling components of arbitrary order by the HTF method. However, the developed SISO impedance model is represented by a series of large-dimensional matrices, which degrades model analyticity and increases computation time.
Besides the model simplification, another obstacle is the modeling complexity. In the describing function-based methods, the concept of perturbation injection is used for submodule modeling [29]- [32]. To solve the multiplication of injected perturbation and time-periodic gain, complex convolution is required. Moreover, the form of perturbations needs to be defined empirically. Considering more perturbations increase the modeling complexity while considering fewer perturbations degrade the model accuracy. As for the HTF modeling, each control loop of the VSC is modeled separately as a Toeplitz matrix, and then the models are combined based on their interconnections [26]- [28], [34]. In the modeling process, a series of matrix operations are required, especially the matrix inversion, which significantly increases the modeling complexity. Note that each row/column of the Toeplitz matrix can be obtained by frequency shifting another row/column. This observation implies that the information of the Toeplitz matrix is redundant in describing the FCEs.
In order to address the abovementioned challenges, a recursive SISO impedance modeling framework for single-phase voltage source rectifiers (VSRs) is proposed in this article. The main contributions of the article are summarized as follows.
1) A general and simple LTP modeling method is proposed by extending the LTI modeling. It does not require complex convolution and large-dimensional matrix operations for the submodule modeling. Therefore, the modeling complexity is reduced. 2) The recursive SISO impedance model is analytical and easy to include harmonics of arbitrary order based on the required model accuracy. Compared with the HTF-based SISO model [34], the recursive SISO impedance model avoids complicated matrix inversion. The rest of this article is organized as follows. In Section II, the basic procedures of the LTP modeling approach are introduced. In Section III, the single-phase VSR system configuration is described. Section IV presents the harmonic admittance modeling under no grid impedance. In Section V, the recursive SISO impedance model is derived and a detailed comparison with other existing models is given. In Section VI, the developed impedance model is validated by experiments in terms of impedance measurement and stability analysis. Finally, Section VII concludes this article.

II. FROM LTI MODELING TO LTP MODELING
The information contained in the Toeplitz matrix is redundant, which increases the modeling complexity. In fact, one row/column of the Toeplitz matrix is enough to represent the complete system model. In the light of this observation, a simple LTP modeling is introduced as follows.

A. LTI Modeling Revisit
To better understand the LTP modeling, a brief revisit of LTI modeling is given first. Consider a one-dimensional nonlinear time-invariant system where xࢠR and uࢠR represent a state variable and an input variable, respectively, and f is a nonlinear function from R 2 to R. This system is autonomous, and its small-signal model can be derived as (2) where the superscripts "∼" and "-" represent the small-signal and steady-state quantities, respectively. Both α and β are constant values determined by the steady-state operating points.
Applying Laplace transform to (2), it gives that where s = jω and ω denotes the frequency variable. Then, the transfer function from input to state is obtained as By repeating the operations from (1) to (4), the small-signal model of a nonlinear time-invariant system in the frequency domain can be derived.

B. LTP Modeling Approach
The dynamics of a one-dimensional nonlinear time-periodic system is described as with its periodicity T 0 as The model of the single-phase grid-connected system is a typical example of this kind because of the steady-state timeperiodic operating trajectories.
The linearized model of such a nonautonomous system is where α(t) and β(t) are determined by the steady-state operating trajectoriesx(t) andū(t), satisfying As observed from (7), it is not able to derive the transfer function from the input to the state directly due to the timeperiodic coefficients.
According to the Fourier transform theory, periodic signals satisfying Dirichlet conditions could be represented by a complex exponential series of infinite dimensions [35]. Thus, (8) can be rewritten by the Fourier series as where ω 0 = 2π/T 0 . And a k and b k are complex Fourier coefficients, which are given by Substituting (9) and (10) into (7) yields Then, according to the displacement property of the Laplace transform, the frequency domain expression of (13) is It is found thatx at a given frequency ω is related tox andũ not only at frequency ω, but also multiple sideband frequencies ω ± kω 0 (kࢠZ). This proves the existence of the frequency-coupling effect in the time-periodic systems.
It is worth noting that (14) corresponds to one row of the Toeplitz matrix of the LTP system. Moreover, due to the arbitrariness of s, another row can be easily derived by frequency shifting as By applying the frequency-shifting operation on (14) repeatedly, the frequency response from input to any order of harmonic can be constructed. The total dynamic propagations of the LTP system in the frequency domain are shown in Fig. 1.
In this way, the use of matrix operation can be avoided and complex convolution is not required. Thus, the modeling complexity can be reduced. Also, the number of frequency-shifting operations required depends on the number of nonzero Fourier coefficients. In the practical analysis, truncation is usually required and the truncation order depends on the required accuracy [34]. By solving the obtained frequency-shifting equations, the SISO transfer function fromũ(s) tox(s) with consideration of frequency-coupling effects can be derived.

III. SYSTEM CONFIGURATION
The circuit configuration and control scheme of the studied single-phase VSR are depicted in Fig. 2. The main circuit is comprised of a power supply u g , a grid impedance L g /R g , an input inductance L f in series with a resistance R f , an H-bridge circuit, an output capacitance C dc , and a load resistance R dc .  In the control diagram, a second-order generalized integrator (SOGI)-PLL is used to obtain the phase θ of the point of common coupling (PCC) voltage u i . And the dc-link voltage u dc and input current i g are regulated by a voltage loop and a current loop, respectively.
The studied system is split into two subsystems at PCC: one is the source subsystem and the other is the load subsystem, as shown in Fig. 3. The analysis of source and load systems can separately bring to an impedance-based representation for both of them. Based on the input impedance Z(s) and the grid impedance Z g (s), the impedance-based stability criteria can be applied to determine the system stability.
In this system, the sampling delay is modeled as where ω i and ω v are the cutoff frequency of anti-aliasing filters for current and voltage measurement, respectively. Besides, an SOGI, a notch filter (NOTCH), two proportionalintegral (PI) regulators, and a proportional-resonant (PR) regulator are given as where ω 1 is the angular frequency of the grid, ξ is the damping ratio, σ is bandwidth coefficient, ω n is the center frequency, k pm and k im (m = 1, 2) are the PI controller parameters, and k pc and k rc are the PR controller parameters.
Zero-order hold effects and calculation delays are taken into consideration as where T s is the control period.

A. SOGI-PLL Modeling
The SOGI-PLL is used to track the phase θ of the PCC voltage u i , which is used for generating current reference. The SOGI generates the π/4-delay signal u iβ of u iα . Then, applying the Park transformation to u iα and u iβ , the q axis voltage u q is obtained Perturbing the variables in the time-domain, the linearized form of (23) is derived as where ω 1 t and V 1 are denoted as the phase and amplitude of the PCC voltage in steady-state, respectively. According to (13) and (14), (24) can be rewritten as Based on the SOGI-PLL diagram, the following relationships hold:ũ where H PLL (s) = PI 1 (s)/s. Substituting (26)- (28) into (25), it is deduced that where Applying the frequency shifting to (25)-(28) by ±kω 1 (kࢠZ), the harmonic signal-flow graph of the SOGI-PLL is drawn as Fig. 4, where s k is defined as s + jkω 1 (kࢠZ). It is observed that the perturbed phaseθ at a given frequency ω is related to the PCC voltageũ i at the two coupled frequencies ω ± ω 1 .

B. Voltage Loop Modeling
Assume that the power losses on the switches of the rectifier are neglectable, the power relationship between the ac side and the dc side in Fig. 2(a) can be described as Linearizing (33) gives (34) shown at the bottom of this page, and its corresponding frequency domain expression is obtained as (35) shown at the bottom of this page. I 1 and ϕ i1 correspond to the magnitude and initial phase of the fundamental current in steady-state, respectively. And G L (s), G x1 (s), and G x2 (s) are given by where Based on the above calculations, the harmonic signal-flow graph of the voltage loop can be drawn as Fig. 5. In the voltage loop,Ĩ * 1 at a given frequency ω is dependent onũ i andĩ g at two coupled frequencies ω ± ω 1 .

C. Current Loop Modeling
The grid current reference is The corresponding small-signal form of (41) is given as In the frequency domain, (42) is rewritten as Due to the current regulator, the perturbed converter voltage including the feedforward term is .
Similarly, cross-coupling dynamics exist in the current loop thatĩ * g at a given frequency ω is coupled withθ andĨ * 1 at two coupled frequencies ω±ω 1 . And the harmonic signal-flow graph of the current loop is depicted as Fig. 6.

D. Harmonic Admittance Modeling
The harmonic signal-flow graphs show that the single-phase VSR system is a MIMO system essentially, which is different from the usual time-invariant system. According to the above analysis, it is worth noting that different signal-flow paths are translated copies of each other with frequency shifting proportional to the fundamental frequency ω 1 . Therefore, the system model can be extended to include harmonics of arbitrary order by frequency-shifting any one signal-flow path without complicated calculations. The current dynamics in the average model is Substituting (29), (39), and (44) into (45), the dynamics between voltage and current at PCC is derived as (46) shown at the bottom of this page, where G i,-2 (s), G i,0 (s), G i,2 (s), G u,-2 (s), G u,0 (s), and G u,2 (s) are given in Appendix.
Neglecting the coupling terms in (44), the conventional input impedance is obtained as In order to include the interactions among different frequency domains, frequency-shifting (46) to s ± j2ω 1 gives where only the perturbed voltage termũ i (s) is retained due to no grid impedance and the harmonic currents with quadruple line frequency deviation are assumed to be zero. Substituting (48) and (49) into (46), the input impedance with capturing FCEs under no grid impedance is derived as (50) shown at the bottom of this page.
where Y op (s) = 1/Z op (s). The small-signal modeling illustrates that a perturbed PCC voltageũ i at a given frequency ω will generate three harmonic grid currentsĩ g at the frequencies ω and ω±ω 1 , as depicted in Fig. 7. The existence of the FCE described by (51) and (52) poses a challenge to the modeling of the VSC-grid interaction.

V. RECURSIVE SISO IMPEDANCE MODELING AND ANALYSIS
When there is a grid impedance on the grid side, the harmonic currents would generate the voltage at the same frequency and thus constitute additional feedback loops. The multifrequency diagram of a single-phase VSR with the additional closed loops is depicted in Fig. 8. The additional loops are the translated copies of the fundamental loop (denoted by the black line) with frequency shifting ±kω 1 . Moreover, the order of harmonic would extend to infinite at the presence of the grid impedance.

A. Recursive SISO Impedance Modeling
According to Fig. 8, the relationship betweenĩ g andũ i in different frequency domains can be consistently established as  where KࢠZ. In addition, due to the existence of grid impedance, it is established that where KࢠZ\0.
In order to obtain a model that is suitable for analysis, truncation is required. Assume that the multiloop system is truncated to P (PࢠN) positive loops and N (NࢠN) negative loops as depicted in Fig. 9, the truncation conditions are expressed as where XࢠN and YࢠN. The positive loop truncation is taken as an example to illustrate how to obtain a truncated model. Based on the idea of mathematical induction, assume first that where 1 ≤ X ≤ P + 1, F P (P − (X − 1)) is a mapping function fromũ i (s 2(X-1) ) toũ i (s 2X ) that needs to be determined, and F P (0) = 0.
In the same way, a recursive relationship for the negative loop truncation can be derived as is given as (61), shown at the bottom of this page. Finally, substituting X = 1 into (56) and substituting Y = 1 into (60) yield Substituting K = 0 into (53) gives Combining (62)-(64), a recursive SISO impedance model can be derived as .
(65) The recursive SISO impedance model has an analytical expression. Moreover, the recursive SISO impedance is easy to incorporate harmonics of arbitrary order by only four transfer functions (Y op (s), Y p (s), Y n (s), and Y g (s)) with their frequencyshifting copies. In addition, the matrix inversion is avoided compared with [34]. When no grid impedance exists (Z g (s) = 0), it yields F N (•) = 0 and F P (•) = 0 such that the recursive impedance Z(s) degrades to Z op (s).
Overall, the recursive SISO impedance modeling method is summarized as Fig. 10, which includes two main steps. The first step is to establish the harmonic admittance model under no grid impedance condition, which can be deemed as the fundamental frequency coupling. The next step is to construct the complete system dynamics by frequency-shifting the fundamental frequency coupling, and then the recursive SISO impedance model can be derived based on the idea of mathematical induction.
Remark: Considering a multiphase system under linear grid impedance, the relationship between the voltage vector and the grid current vector in an M-phase system (MࢠN + ) can be expressed as N 1 ࢠN, N 2 ࢠN, and KࢠZ. And Y pK , Y 0 , and Y nK are all definite M-by-M matrices. The additional loop resulted from the grid admittance can be expressed as The truncation conditions are expressed as (68) Fig. 11. Comparative analysis of SISO impedance with different truncation orders under L g = 5.5 mH and R g = 1.5 Ω.

TABLE I MAIN CIRCUIT AND CONTROL PARAMETERS
According to the idea of mathematical induction, two recursive relationships can be derived as where F pm and F nl are two recursive functions. Combining (66), (69), and (70) and performing the frequencyshifting, the relationship between i(s 0 ) and u(s 0 ) can be solved. Therefore, the developed recursive SISO modeling method has a great potential of being promoted to other time-periodic systems.

B. Frequency Response Analysis of Recursive Impedance
The model accuracy depends on the chosen truncation order. When finite-order truncation is performed, frequency response deviation may occur, as shown in Fig. 11 (main parameters are listed in Table I). The selection of truncation order mainly influences the frequency response of Z(s) around the fundamental frequency. And the frequency response deviation decreases  as the truncation order increases. The evaluation of truncation order has also been discussed in the literature [34], but an explicit formula has not yet been provided. Thus, the truncation parameters P and N are usually evaluated iteratively. In the work presented in the following, N = 3 and P = 3 have been chosen.
In order to better visualize the influence of grid impedance, the impedance model Z(s) under different grid impedances are plotted in Fig. 12. The symbol Z X (s) (X = A, B, C, D) are used to represent Z(s) under different grid impedances (Z A (s): L g = 0 mH, R g = 0 Ω; Z B (s): L g = 1.5 mH, R g = 1 Ω; Z C (s): L g = 3.5 mH, R g = 1 Ω; Z D (s): L g = 5.5 mH, R g = 1 Ω). From Fig. 12, it is observed that evident discrepancies occur around the fundamental frequency as the grid impedance changes, which shows that the rectifier impedance depends on the grid impedance. There exist two valley points in the low-frequency regions, whose frequencies are symmetric about the fundamental frequency. And, with the increase of the grid impedance, the magnitude gain of the rectifier impedance decreases, implying that the ability to suppress voltage disturbances near the valley frequencies degrades.
Based on the MATLAB platform, the computation time required to calculate the impedance response of the HTF matrixbased model [34] and the proposed recursive SISO model at one thousand frequency points is shown in Fig. 13. As seen, the computation time required by the recursive SISO model is always less than that of the HTF matrix-based model, verifying the superiority of the recursive SISO model in the computational burden. The reason for the superiority lies in the elimination of large-dimensional matrix operations. The recursive SISO model can save approximately 16.7% time when the truncation order is 1. Moreover, the polynomial fitting curve shows that the computation complexity of the HTF matrix-based model will increases exponentially as the truncation order increases. In contrast, the computation complexity of the recursive SISO model increases almost linearly as the truncation order increases.

C. Comparison Between the Existing Impedance Models and the Proposed Recursive Impedance Model
Many modeling methods can lead to accurate models, but their modeling procedures and the resulting models are different. Table II summarizes the comparison results of the existing impedance modeling methods and the proposed modeling method for the single-phase converters. As shown in Table II, the modeling methods can be mainly categorized into two types: describing-function-based methods and linearizationbased methods.
The harmonic linearization modeling method and the harmonic signal-flow graph modeling method can be directly applied to the nonlinear time-periodic systems without linearization. However, the additional closed-loops caused by the grid impedance are neglected and complex convolutions are needed for the submodule modeling. The multifrequency modeling method, the HTF modeling method, and the proposed modeling method all directly linearize the nonlinear time-periodic system around ac steady-state trajectories. The multifrequency modeling method considers a limited number of frequency coupling loops, which may be a potential limitation for model accuracy. The use of large-dimensional matrices degrades model analyticity and increases the computation effort of the HTF model. The proposed recursive SISO impedance modeling method gives an analytical model, which can include the harmonics of arbitrary order with low computation burdens. Thus, it is a promising choice for single-phase converter modeling.

VI. EXPERIMENTAL VALIDATION
In this section, a prototype of the single-phase VSR is built to verify the correctness of the proposed models, as shown in Fig. 14. The system specifications are the same as Table I. The control platform is based on a floating-point digital signal processor (DSP) TMS320F28335. The grid voltage and the injected small perturbation from 10 Hz to 1 kHz are generated by a regenerative grid simulator (Chroma 61830). And the measured signals are sent to the frequency response analyzer Bode 100.

A. Frequency Scan Validation
The implementation of the impedance measurement scheme is depicted in Fig. 15. A small-signal perturbation voltage u p is injected, which excites the perturbations of corresponding frequency on the PCC voltage and input current. The PCC voltage signal u i sends to the CH1 port, and the grid current signal i g transmits to the CH2 port. The measured impedance is obtained by plotting the bode diagram of CH1/CH2.
The frequency response characteristics of the built impedance models and the corresponding measured results under different  grid impedances are plotted in Fig. 16. As shown, the conventional impedance Z c (s), which only considers the fundamental perturbed components, has a large deviation in the regions around the fundamental frequency. The Z op (s) achieves higher accuracy compared with Z c (s), but the accuracy degrades as the grid impedance increases. In contrast, the measured results are in good accordance with Z(s), which validates the accuracy of the  proposed impedance modeling approach. Since some unmodeled perturbations (e.g., the nonlinearity of the passive components, precision of measurement device, parasitic parameters) exist in the actual experimental setup, there still exists a small deviation between the measured results and calculated results of Z(s), which is a common problem in impedance measurement [30], [33].
The experimental measurement shows that the couplings between the grid impedance and the rectifier play an important role in the model accuracy. Thus, it is suggested that the interactions resulted from the grid impedance should be considered in the input impedance modeling of the single-phase VSR system.

B. Stability Analysis
This section will further analyze the validity of the proposed SISO models in terms of small-signal stability. The system stability is assessed by analyzing the zero effects of 1 + Z Y (s)/Z g (s) [28], [36], where Z Y = Z c , Z op , and Z. First, with a fixed voltage loop bandwidth f b = 21 Hz, different SISO impedance models under Z g = 3 mH/1 Ω (i.e., Case 1) and Z g = 4.5 mH/1 Ω (i.e., Case 2) are compared in Fig. 17. The frequency response curves concerning Z c and Z op remain almost the same under different Z g , and the small deviation is caused by the increased steady-state voltage drop across the grid impedance.
By analyzing Fig. 17(a), no right-half plane (RHP) zero effect is observed, indicating a stable system in both cases. On the contrary, both cases are predicted to be unstable by Z op as shown in Fig. 17(b), where the complex RHP zero effect is found at 66 Hz. From 17(c), only complex left-half plane (LHP) zero effects are found in Case 1, indicating a stable system. While for Case 2, the complex RHP zero effects are found at 34.6 and 65.4 Hz, meaning that the system is unstable. The harmonic and stability identified by different models are summarized in Table III.
Experimental waveforms under Cases 1 and 2 are shown in Figs. 18 and 19. When Z g = 3 mH/1 Ω, no unstable oscillation exists in both the grid current i g and dc-link voltage u dc . However, when the grid impedance Z g increases to 4.5 mH/1 Ω, the grid current i g and the dc-link voltage u dc exhibit low-frequency oscillation, indicating an unstable system. By analyzing the spectrum of i g , the oscillation frequencies of 34 and 66 Hz are all identified as shown in Fig. 19(b), which closely agrees with the analyzed results in Fig. 17(c). Based on the above discussion, the neglect of frequency couplings will lead to inaccurate stability results.
Further, the influence of voltage loop bandwidth f b on the system stability is analyzed. Fig. 20 shows the frequency response curves of 1 + Z(s)/Z g (s) under a fixed grid impedance  (L g = 5.5 mH, R g = 1 Ω) when the voltage loop bandwidth f b is configured as 19 and 21 Hz, respectively. According to Fig. 20(a), only complex LHP zero effects are found when f b = 19 Hz, indicating a stable system. While for f b = 21 Hz as shown in Fig. 20(b), the complex RHP zero effects are found at 35 and 65 Hz, implying that the system is unstable. The   Table IV.
Experiments are carried out with f b = 19 Hz and f b = 21 Hz in Figs. 21 and 22, respectively. It is seen that for f b = 19 Hz, the grid current i g remains sinusoidal and the average value of dc-link voltage u dc tracks its reference, meaning that this system is stable. Whereas for f b = 21 Hz, instability occurs that dc-link voltage u dc and the grid current i g oscillate with low frequency. The FFT analysis of i g is shown in Fig. 22(b), where the oscillation frequencies of 34.5 and 65.5 Hz can be identified. Similarly, the observed oscillation frequencies closely agree with the analyzed results in Fig. 20.
The stability analysis demonstrates that the proposed impedance model is accurate for stability prediction. The impedance models Z c and Z op have limitations in stability analysis due to the neglect of FCEs. Additionally, the small-signal stability analysis shows that smaller grid impedance Z g and lower voltage loop bandwidth f b can achieve better stability in terms of the VSC-grid interaction. Increasing grid impedance Z g or voltage loop bandwidth f b will move two complex LHP zeros to the RHP, and thus the low-frequency oscillation occurs.

VII. CONCLUSION
In this article, a recursive SISO impedance modeling framework for the single-phase converters is proposed. Through the displacement property of the Laplace transform and the frequency-shifting operation, the traditional LTI modeling is extended to LTP modeling. The LTP modeling method transferred the time-varying model into a time-invariant transfer function vector model. The complex convolution and the Toeplitz matrix are eliminated in the submodule modeling, so the modeling complexity is reduced. Furthermore, based on the idea of mathematical induction, a recursive SISO impedance model can be obtained by solving the additional loops closed by the grid impedance. The recursive SISO impedance model is analytical and accurate. Moreover, compared with the HTF-matrix-based SISO impedance model, the recursive SISO impedance model has lower computation burdens.
On the other hand, it is found that both the increase of voltage loop bandwidth and grid impedance lead to low-frequency oscillations. Also, the typical impedance without considering any FCE is over-optimistic on the stability assessment, whereas the impedance neglecting the grid impedance effect is overconservative in terms of stability. An accurate stability analysis must take all the frequency coupling components into account.
The proposed modeling method has a great potential of being promoted to other time-periodic systems. The detailed modeling and analysis of the multiphase system under complex grid conditions will be the topic of future research.