Stability Assessment of Power-Converter-Based AC systems by LTP Theory: Eigenvalue Analysis and Harmonic Impedance Estimation

Stability analysis of power-converter-based ac systems poses serious challenges not only because of the nonlinear nature of power converters, but also because linearization is not generally applied around a steady-state operating point, as in the dc case, but around a time-periodic operating trajectory. Typical examples are single-phase and unbalanced three-phase systems. In this paper, two general methods to assess stability of the aforementioned systems are presented. Both are based on the linear time periodic (LTP) systems theory. The first is model-based and relies on the evaluation of the eigenvalues of the linearized model, assuming a complete knowledge of the parameters. By contrast, the second proposes a set of small-signal current injections to measure the harmonic impedances and applies the LTP Nyquist criterion, so that the stability of the system can be assessed with a black-box approach, without relying on the knowledge of the system parameters. The basic LTP systems theory is reviewed in order to provide a mathematical justification for the second method. As case study, a simple network, consisting of a source full-bridge converter in ac voltage-control mode and a load full-bridge converter in ac current-control mode including phase locked loop, is considered. Analytical results based on average modeling and simulations based on both average and switching models are presented, showing good accuracy in the identification of the stability thresholds for both the proposed methods.

Abstract-Stability analysis of power-converter-based AC systems poses serious challenges not only because of the non-linear nature of power converters, but also because linearisation is not generally applied around a steady-state operating point, as in the DC case, but around a time-periodic operating trajectory. Typical examples are single-phase and unbalanced three-phase systems. In this paper, two general methods to assess stability of the aforementioned systems are presented. Both are based on the Linear Time Periodic (LTP) systems theory. The first is model-based and relies on the evaluation of the eigenvalues of the linearised model, assuming a complete knowledge of the parameters. By contrast, the second proposes a set of smallsignal current injections to measure the Harmonic Impedances and applies the LTP Nyquist Criterion, so that stability of the system can be assessed with a black-box approach, without relying on knowledge of the system parameters. The basic LTP systems theory is reviewed in order to provide a mathematical justification for the second method. As case study, a simple network, consisting of a source full-bridge converter in AC voltage-control mode and a load full-bridge converter in AC current-control mode including PLL, is considered. Analytical results based on average modelling and simulations based on both average and switching models are presented, showing good accuracy in the identification of the stability thresholds for both the proposed methods. I N recent decades, improvements in Power Electronics have led to a significant spread of grid-connected converters, both for AC and DC applications. These systems, which are intrinsically non-linear due to the presence of switching elements, have been a challenging field of investigation for researchers. Power Electronics guarantees high conversion efficiency and high controllability of loads and sources. In fact, it is possible to achieve high dynamic performances using sophisticated controls and modulation schemes. However, these systems often show unstable behaviour, especially due to the harmonic couplings at the common connection points between sources and loads.
In general, there are two main possible approaches for stability analysis of non-linear systems: linearise the system in order to use well known linear analysis techniques, or use the Lyapunov stability criterion directly on the non-linear model of the system. The application of the latter in practical applications is usually complicated, so most stability analysis techniques rely on the linearised model of the system. One of the first investigations of these instabilities was presented by Middlebrook in [1], where the interaction between closed loop controlled DC-DC converters and their input filters was the cause of instability. This cascaded system was found to be stable only if the ratio of the source output impedance to the load input impedance satisfies the Nyquist stability criterion, with analysis based on the linearised model of the converter.
Other scenarios where similar instability issues arise are reported in [2] and [3] for different single phase railway distribution systems. In [4], the unstable behaviour of LCC HVDC converters is discussed. The analysis [5] shows that renewable energy systems involving grid-parallel inverters are affected by instability issues, and subsequently in [6] and [7] a report is provided of ongoing research, current harmonicrelated issues and future expected problems.
In case of AC systems, the Middlebrook approach is still applicable, but must be adapted slightly. In [8] an equivalent input-output impedance approach is proposed for the analysis of balanced three-phase AC power systems, in dualism with the DC case, based on the study of the synchronous dq frame return ratio matrix Z s,dq (s)Z −1 l,dq (s), exploiting multi-variable linear control theory. Later, Belkhayat [9] proposed a further improvement to this method by employing the Generalised Nyquist Criterion (GNC), which was originally introduced by MacFarlane and Postlethwaite, [10], and which extends the SISO frequency-response control methods to the MIMO case. Other examples of stability analysis of three-phase balanced AC systems are provided in [11], based on dq measurements of small-signal impedances. In fact, in the DC case the non-linear system is linearised about a fixed operating point, leading to an LTI linearised system whose stability is assessed with LTI techniques. For the AC balanced three-phase case, the system is first transformed into a MIMO system in the dq reference frame, then linearisation is performed as in the DC case and MIMO LTI techniques are exploited for the stability analysis. However, for single-phase and for unbalanced three-phase AC systems the methods described above are not suitable. The main reason is the linearisation process, that in this case must be performed around a steady-state time-periodic trajectory rather than around a fixed steady-state operating point.
There are several techniques available in literature to analyse single phase and generic unbalanced three phase systems. One of the most widely used is Harmonic Linearisation, where small-signal perturbations are used to determine the input-output small-signal impedances [12], [13], [14], and the Nyquist criterion is applied to their ratio. However, this method might be inaccurate in determining instabilities below the fundamental grid line frequency, as shown in [15]. An extension of this method is provided in [16], where a two-dimensional admittance for single-phase voltage source converters application is defined, able to capture the crosscoupling effects (up to twice the line frequency) and possibly overcome the aforementioned limitation.
Another method is the Dynamic Phasor approach, where the analysed system is described through a small-signal impedance-based model. In [17] an application example to a single-phase system is provided, with a Fourier series analysis truncated at the 1-phasor, which keeps the mathematical formulation closer to the dq transformation approach and therefore use the GNC to assess stability.
The method used in this paper is based on the Linear Time Periodic (LTP) systems theory, reported by Wereley and Hall in [18] and [19]. Single-phase and unbalanced three-phase systems are described by average models that are Non-Linear Time Periodic, where all the state-space variables are periodic with frequency equal to the fundamental line frequency or a multiple. Thus, linearising around the steady-state operating trajectories, the time-domain LTP model is derived, as well as the frequency-domain Harmonic State-Space (HSS) model. Stability analysis can then be performed on the latter, either evaluating the eigenvalue loci plot or exploiting the LTP Nyquist criterion. The Harmonic Transfer Function (HTF) matrix-operator describes the interaction of the harmonic coefficients of the input and output signals, giving a precise description of the cross-coupling effects. The formulation of the LTP analysis permits easy inclusion of an arbitrary number of harmonic coefficients in the stability assessment, and this is an advantage compared to Harmonic Linearisation and Dynamic Phasor approaches, where truncations are usually introduced to deal with the more complicated mathematical formulation. In [20] and [21] the HSS model and stability analysis are addressed for a single-phase and a three-phase grid connected converter, but without considering the dynamics due to the PLL and the impact of digital-computational delays. In [22] an impedance-based stability method is derived and a singlephase grid-connected converter with PLL is used as a test case; however, precise stability boundaries are not provided. In [23] the same circuit is analysed using an eigenvaluebased approach, with use of LTP theory. Precise stability boundaries can be evaluated, although the model does not include important elements such as digital computation delay, ZOH and PWM dynamics.
It is worth pointing out that the design of the fast control loops, i.e. inner current control, is normally carried out on the assumption that the system is linear, and ignoring potential cross-converter interactions caused by the non-linear loops such as PLLs, Power controllers, DC-link controllers. On the other hand, when modelling the impact of these non-linear loops on an interconnection of converters, the effect of the fast dynamics could often be neglected. To keep the analysis as general as possible, in this work stability assessment will be performed including the effect of all the dynamics, with the only approximation being that PWM and digital implementation delays will be modelled with their continuous-time estimate.
In this paper, LTP theory is exploited in an innovative way: based on the HSS model of the system, rewritten as an interconnection between a load and source system, the Harmonic Impedances are calculated for both subsystems using a small-perturbation current injection method. Thus a closed-loop transfer function is defined and the LTP Nyquist criterion is used to assess the stability of the overall system. The method is shown to be equivalent to the eigenvalue analysis, and good accuracy is obtained in the identification of the stability boundaries. The theoretical foundations of a method to measure these Harmonic Impedances is also provided, ensuring the practically relevant feature that this method is no longer model-based but rather can be applied to any AC network, provided that the harmonic impedances can be measured. Analytical results, as well as simulations based on both average and switching models of the two converter systems are given for both the eigenvalue analysis and the Harmonic Impedances method, showing good agreement and confirming the effectiveness of the proposed analysis. The paper is organized as follows: Section II provides a review of LTP theory; Section III describes the system used as case study and derives the non-linear average model. In Section IV, the steady-state time-periodic solution of the system is derived and in Section V the LTP model is calculated. In Section VI analytical and simulation results based on the average model are given and in Section VII simulation results based on a complete switching model are provided.

II. REVIEW OF LINEAR TIME PERIODIC SYSTEMS THEORY
Given a general Non-Linear Time Periodic (NLTP) system, with T -periodic state-space variables, T = 2π/ω T = 1/f T : (1) and given a steady-state T -periodic input,ū(t), the system equations can be solved either applying Harmonic Balance, [18], or numerically (in Matlab for example) and the periodic steady-state solutions,x(t),ȳ(t), can be obtained. To proceed with stability analysis, the LTP model is derived by linearising the NLTP system around the calculated steady-state solutions. In contrast to DC systems, where the steady-state operating point is constant, AC systems have a time-periodic steady-state operating trajectory. Thus, linearisation is applied by adding small-signal perturbation terms around the periodic steadystate operating trajectory: with |ũ(t)| |ū(t)|, |ỹ(t)| |ȳ(t)| and |x(t)| |x(t)|. Then (2) is substituted into (1): and the LTP system is then obtained by taking into account only first-order terms, ignoring steady-state and second-order terms: with the matrices A(t), B(t), C(t) and D(t) being T -periodic. The test input signal for the frequency analysis of LTP systems is the Exponentially Modulated Periodic (EMP) signal [18], [19], which is defined as: When such a signal is given as input to the LTP system, all the state-space variables and the output are EMP signals. Thus input and output spaces are consistent and it is possible to define a transfer function operator. Substituting (5) into (4) and applying the Fourier matrix expansion to the matrices A(t), B(t), C(t) and D(t) allows the system to be rearranged as: This set of equations is not convenient for analytical manipulation because of the convolution products. A simpler form can be derived by applying the Toeplitz transform to the periodic matrices, defined as follows for the generic matrix P (t): which is a doubly infinite block Toeplitz matrix and the blocks P i are the Fourier matrix coefficients of the T -periodic matrix P (t). This definition applies to the matrices The Toeplitz form of (6) provides a more compact description of the system, enabling the definition of the Harmonic State-Space (HSS) model associated with the LTP system (4): . . , N n , . . . ) and N n being a square diagonal matrix of the same dimensions as A n with diagonal coefficients equal to jnω T . Rearranging (8), the relationship between the harmonic coefficients of the input and output signals is represented by the Harmonic Transfer Function (HTF) operator of the system, defined as: where I is the Toeplitz form of the identity matrix. Stability analysis can now be performed in two equivalent ways.

A. Eigenvalue Loci
This method requires a full model of the system, i.e. knowledge of all the component values and controller parameters. Stability analysis is performed by evaluating the eigenvalues of the matrix (A − N ). If all the eigenvalues have Re[λ i ] ≤ 0, where those with Re[λ i ] = 0 have algebraic multiplicity equal to 1, then the system is stable, otherwise the system is unstable.
In real applications, only a finite number of harmonic coefficients can be considered, for computational reasons. A truncation of the Fourier series in (6) is applied, defining the truncation order M , which refers to the number of harmonics taken into account. If M = 2, for example, the DC-component and the first and second harmonics are considered. The corresponding truncated Fourier expansion involves the Fourier coefficients n = −2, −1, 0, 1, 2 and yields the approximation and the following associated truncated Toeplitz form: with Z A being a zero matrix of the same dimension as the A n . For a general LTP system of order k and truncation order M , where k refers to the total number of state-space variables, A − N is a square (2M + 1)k × (2M + 1)k matrix, with the total number of associated eigenvalues equal to (2M + 1)k. Among these, the number of significant eigenvalues equals the order k, while the remaining 2M k are translated copies of the significant ones, with translation equal to jnω T , n = ±1, . . . , ±M . Fig.1 reports an example of LTP eigenvalue loci with k = 4 and M = 2. There are four vertical lines of eigenvalues and a total number of eigenvalues equal to (2 × 2+1)×4 = 20. The significant ones are depicted in red, while their translated copies are in black. Increasing the truncation order M results in longer vertical lines of eigenvalues. For example, incrementing the truncation order to M = 3 will result in eight more eigenvalues, since now there is a total number of (2 × 3 + 1) × 4 = 28 eigenvalues. Thus each line will be incremented with two more eigenvalues, one located on the top and the other on the bottom of the eigenvalue line, giving rise to two additional translated copies of the significant eigenvalue. The system associated with the eigenvalue loci plot in Fig.1 is clearly unstable, since there is one significant eigenvalue in the right-half plane. Remark: the locations of the eigenvalues depend on the chosen truncation order. HSS analysis is defined on the assumption that an infinite truncation order is used, so that all the significant eigenvalues (plus their translated copies) are correctly located in the complex plane. When a finite truncation order is used, the eigenvalues are shifted, and to ensure the robustness of the numerical results a truncation order M = M * is determined, such that for M > M * , the eigenvalue loci plot exhibits no noticeable shift compared to ones evaluated for M M * , which provide the correct eigenvalue locations; for M < M * there is a noticeable shift of the eigenvalue locations in the complex plane. The evaluation of M * has been addressed in the literature [24], but an explicit formula has not yet been provided. Thus, M * is usually evaluated iteratively. In the work presented below, M * = 23 has been found. Hence, to guarantee an accurate analysis, in section VI and VII an order M = 40 is used.

B. LTP Nyquist Criterion
LTP Nyquist Criterion does not necessarily require full knowledge of the model; this can be extracted through smallsignal current injections, thus making the method suitable for a generic "black-box" system, provided that the AC terminals are accessible. Fig.2 illustrates the case of a grid-connected power converter, where a small-current injection is performed at the interface between the AC grid and the power converter. Such an injection is small enough to not affect the steady-state operation of the system. The system in Fig.2

is linearised
and the LTP model is derived as in (4), and for convenience is reproduced here: where the outputs areỹ(t) = [ỹ 1 (t);ỹ 2 (t);ỹ 3 (t)] = [ĩ g (t);ĩ c (t);ṽ o (t)] and the only input isũ(t) =ĩ x (t). The associated HTF satisfies: where the Toeplitz transform has been applied also to the input and output signals: Rearranging these equations gives us the Harmonic Impedances: The equivalent block diagram is shown in Fig.3, referring to which we have V o = Z g (jω)I g = Z c (jω)I c and I x = I g + I c , thus: The equations in (15) are analogous to the input-output closed-loop transfer functions of a MIMO system; in this case the inputs are the harmonics associated with the injected perturbation current and the outputs are the harmonics associated with the grid and converter perturbation currents, respectively. These equations can be represented by the simple block diagram shown in Fig.4. The LTP Nyquist Criterion, described by Hall and Wereley, is here exploited to assess the stability of the system. For a detailed derivation of the LTP Nyquist Criterion, the reader is encouraged to refer to the original work, [19]. The main result useful for our purposes can be stated as follows. Theorem: consider the open-loop harmonic transfer function Z −1 g (s)Z c (s) (or equivalently Z −1 c (s)Z g (s)) and assume that it has no unobservable or uncontrollable right-half plane poles. Then, the LTP system (4) , included by the LTP Nyquist contour plot, as shown in Fig.6 [19] .

Remark:
1) It is worth noting that the LTP Nyquist contour plot comprises the four segments N A , N B , N C , N D , as shown in Fig.6, and it is used to determine the number of righthalf plane poles of the open-loop HTF. Considering the example in Fig.6, there are three right-half plane poles, thus the LTP Nyquist plot should have three counterclockwise encirclements around the (−1, 0) point, for a stable system. In practice, however, the Nyquist plot is evaluated only for s ∈ N A , since the other segments of the LTP Nyquist contour do not contribute to generate encirclements around the critical point [3]. 2) It follows by calculation that: III. CASE STUDY AND AVERAGE MODEL The system considered as case study is a single-phase network made up of a voltage-controlled full-bridge inverter (source converter) and a current-controlled full-bridge inverter (load converter), connected through an L 2 C 1 L 1 low-pass filter. The voltage control of the source converter ensures that the voltage v o (t) follows the reference v g (t) = V g sin(ω g t), while the load converter transfers power from the DC source V dc2 towards the source converter. To obtain unity-power-factor operation mode, the current-controlled inverter is synchronized  (19) with the source voltage-converter through a Phase Locked Loop (PLL), where the in-quadrature component of the input voltage v o (t) is estimated using a linear filter that introduces a T g /4 delay at ω g : D(s) = ω 2 g /(s 2 + sω g + ω 2 g ). The amount of power transferred by the load converter is controlled through the parameter I ref , which sets the reference amplitude for the current controller. It will be shown that there is a threshold value for this parameter, I * ref , below which the system operates stably, and above which the system shows unstable operation. Such instability is caused by the fact that the PLL is no longer able to generate the correct phase reference. It is worth noting that, in the unstable operation mode, if the current controller is driven with the correct phase rather than the one provided by the PLL, the system does not show unstable behaviour, confirming that the instability is in fact caused by the PLL and not by the other controllers or by limitations in the power transfer capability [22].
The analysis of the network is based on the average model of the system, considering that PLL instability is expected at relatively low frequencies, where the switching behaviour does not have significant impact on the stability boundaries. The average model is given in (20), which is a fourteenth-order Non-Linear Time Periodic system, as in (1), with all the statespace variables being T g -periodic. To maximise the accuracy of the model while retaining low complexity and to generalise the analysis, the elements modelling the digital control, such as the computational delay (T s = 1/f s ), along with the zeroorder hold (ZOH) delay of the pulse-width modulator (PWM), 0.5T s , are included in the analysis through their continuoustime equivalents [25], [26]. In fact, these digital blocks are represented in the continuous-time domain by the following transfer function: The complex exponential is replaced with a first-order Padé approximation of the form: 6 where the coefficients n 1 , n 0 , d 1 and d 0 have been calculated with the Matlab command pade, with specified time delay equal to T s and order equal to one. Substituting (18) in (17) gives an approximated transfer function of the form: The state-space variables of the non-linear model of the system in Fig.5 have the following physical meaning: x 1 -x 4 describe the internal dynamics of the PLL; x 5 is associated with the PI 1 current controller; x 6 and x 7 are related to the voltage controller PI 3,4 ; x 8 , x 9 and x 10 , x 11 represent the internal dynamics of the computational delay, ZOH and PWM for the current and voltage inverters, respectively; x 12 represents the inductor current, i L1 ; x 13 the inductor current, i L2 ; x 14 the voltage across the capacitor. The outputs are: y 1 the current towards the current-source inverter; y 2 the current towards the voltage-source inverter and y 3 the voltage at the common connection point, v o . Please refer to Fig.5 to identify the location of the different state variables. Explicit indication of the time dependency of the variables , v conv2 andĩ x has been suppressed in the following, in order to reduce the quantity of mathematical notation. v All the theoretical analysis is based on the parameters summarised in Table I. Based on (20), in the following sections the steady-state solution is evaluated, linearisation is applied and stability analysis is performed on the linearised model.

IV. STEADY-STATE SOLUTION
In this section the steady-state solutions of the NLTP system (20) are mathematically evaluated by applying Harmonic Balance. The results are then used to perform linearisation and stability analysis. Since the converter is designed for AC systems, the solutions take the form: x i = |x i | cos(ω g t + arg(x i )) = x i e jωgt + c.c. /2 for i = 1, 2, 5, ..., 14 ;x 3 = ω g t +x 03 ;x 4 =x 04 (21) v q = |v q | cos(ω g t + arg(v q )) = v q e jωgt + c.c. /2 where c.c. stands for complex conjugate of the term preceding it within the square brackets. Substituting these expressions in the eighth equation of model (20) gives: A similar approach is applied to the other equations: These equations, (23)-(35), can now be solved numerically in Matlab and the steady-state solutionsx 5 , ... ,x 14 ,v o ,v conv1 , v conv2 obtained. Next,x 3 is first written in the formx 3 = ω g t +x 03 , theṅ Applying trigonometric simplification gives: so it follows thatẋ 4 (t) = 0. But from the state-space model, again using trigonometric simplifications, it follows that: and soẋ which impliesv o −x 03 = 0 or ±π. In our case,v o =x 03 , which gives the solutions: Finally, the last two equations arex 2 = jω gx1 andx 1 = −jv o , for which the solutionsx 1 ,x 2 are given by: The steady-state solution has been calculated by applying Harmonic Balance. Please note that in [27], a similar method is developed for steady-state analysis of a Modular Multilevel Converter. Regardless the specific application, Harmonic Balance is a convenient technique that can be exploited when a steady-state solution of a system of non-linear differential equations is needed.

V. LINEARISED MODEL
Linearisation around the steady-state solution calculated in the previous section is applied to the model (20), according to (3). Each variable in (20) is replaced by its steady-state solution plus a small-signal perturbation. After substitution, the product of perturbations are neglected, as well as the purely steady-state terms (which balance with one another). Doing so, the LTP small-signal model (43) is derived, in the forṁ with A(t), B(t), C(t) and D(t) being T g -periodic matrices, u(t) =ĩ x (t) being the injected current perturbation andỹ(t) the output vector. Applying Fourier series expansion to these matrices gives the result that only A i , B i for i = −1, 0, 1, and C j , D j for j = 0 are different from zero, with A −1 = A * 1 and B −1 = B * 1 (complex conjugate). All the other coefficients of the Fourier series are zero. The Toeplitz transform is applied and the HSS model is calculated. Below, in section VI, the stability of the system is assessed based on the average model (20), comparing time-domain simulations of the average model with the analytical results from the two LTP analysis methods, i.e. Eigenvalue and LTP Nyquist based on the return-ratio of the Harmonic Impedances, showing a very good accuracy in the identification of the stability boundaries of the system. In section VII, time-domain simulations are performed based on the switching model, showing that the stability boundaries are still well predicted by the LTP analysis. Subsequently, a possible method to measure the Harmonic Impedances from the switching model is presented, allowing one to treat the circuit as a "black-box", providing a practical solution to assess stability using the LTP Nyquist Criterion. It is worth pointing out that in this case the full analysis discussed above is not required, as stability assessment would be purely based on measurements.

VI. AVERAGE MODEL RESULTS
To limit the complexity of the results, and to focus on the effect of a single parameter, the current reference for the load converter, I ref , has been selected as the only variable quantity, while components and controllers are assumed to be constant. In the system under study, the stability boundary, I * ref , lies between 11.3 and 11.4A. For I ref < I * ref , the system operates stably, while for I ref > I * ref , the PLL is no longer able to generate the correct phase reference and the system exhibits unstable operation. These results are supported by the following analysis.

A. Time-domain simulation
Time-domain simulations have been performed in the Matlab-Simulink environment, based on the NLTP system (20). In Fig.7 the time-domain evolution of currents, i L1 (t), i L2 (t), voltage, v o (t), and phase, θ(t) are shown, with the reference current being equal to I ref = 11.3A for t < 0.2s and to I ref = 11.4A for t ≥ 0.2s. It can be observed that the system starts to strongly oscillate for t > 0.5s, which reflects the fact that the two unstable significant eigenvalues, as reported in the next Section, are very close to the imaginary axis, providing a long transient. In Fig.7 (e) − (l) a zoom of both stable and unstable operation regions is reported. For I ref = 11.3A the system is stable, Fig.7 (e) − (h), whereas for I ref = 11.4A the system starts to oscillate until currents and voltages exceed the rated values and the control switches off, Fig.7 (i) − (l).

B. Eigenvalue analysis
Based on the HSS model, obtained from the LTP system (43), the eigenvalue loci plot of the matrix A − N is evaluated with I ref = 11.3, 11.4A. A truncation order M = 40 has been chosen, in order to guarantee a good accuracy of the results and at the same time to keep low the computational complexity. It can be seen from Fig.8 that for I ref = 11.3A all the significant eigenvalues lie in the left half-plane, as required by a stable system, whereas for I ref = 11.4A some significant eigenvalues lie in the right half-plane, confirming that in this case the system is unstable. The significant eigenvalues that cause instability are located at s unst = 1.175 ± j5238. Thus, the natural frequency associated with these modes is ω n = Re[s unst ] 2 + Im[s unst ] 2 = 5240 rad/s, which corresponds to an oscillating frequency f n = 830 Hz, which is consistent with the time-domain simulations, where the frequency of the oscillations is measured to be f n = 820 Hz (see Fig.7(g)). Some spurious eigenvalues arise because of the applied truncation, but they do not have a physical meaning and are therefore neglected.

C. LTP Nyquist Criterion
Using the theory developed in Section II, Harmonic Impedances Z g and Z c are evaluated using the LTP model (43). The estimation of the LTP Nyquist Criterion based on small-signal perturbations will be discussed in Section VII. Considering the truncation order M , these matrices will have the form (44). It can be noted that, for s ∈ (−j∞, +j∞), the elements on the same diagonal of the matrix impedance are shifted copies of each other, with shifting proportional to the fundamental period of the grid T g = 2π/ω g . In fact, Z 0 g (s) and Z 0 g (s + jω g ), evaluated for s ∈ (−j∞, +j∞), give the same amplitude and phase Bode plot, except for a frequency shift of ω g between the two plots. Thus, the Harmonic Impedance associated with the 0-component (which is related to the ω g -component) can be obtained from any Z 0 g (s + jkω g ), k ∈ [−M/2, .., M/2], and the same applies to the Harmonic Impedances of other orders. However, in the application of the LTP Nyquist Criterion, the LTP Nyquist contour is defined by s ∈ (−jω g /2, +jω g /2). Hence, the same range is used to derive the plots of the impedances. The Bode plot of the Harmonic Impedance Z g is obtained as follows: gives the Bode plot for ω ∈ (−ω g /2, +ω g /2), Z 0 g (s+jω g ) gives the Bode plot for ω ∈ (ω g /2, 3ω g /2), and so on; 3) these individual Bode plots are merged together, giving the overall 0-Harmonic Impedance Bode plot.
The same procedure applies to Z c and to the Harmonic Impedances of other orders. It is worth noting that the Harmonic Impedances are different from zero only for the even components, i.e. Z ±q g,c (s) = 0 for q = 1, 3, 5.. . Fig.9 shows these impedances for the system under study: for the grid impedance the only relevant component is the 0 one, since the amplitude of all the other components is significantly lower, while for the converter impedance components are relevant up to the third harmonic. The evaluation of the impedances provides a check on the quality of the selected truncation order by observing the relative amplitudes of the different impedances. It should also be noted that Z q g,c (s) = [Z −q g,c (−s)] * . LTP Nyquist Criterion is now applied, with the return-ratio matrix Z −1 g Z c considered in the following analysis (note that the LTP Nyquist Criterion can also be applied by considering . Two cases are considered: a robustly stable system, with I ref = 10A, and a strongly unstable system, with I ref = 13A. This choice has been made in order to obtain clearer results, which are more easily read from the Nyquist plots. In any event, the same stability boundaries discussed before can be calculated with this approach. The Nyquist plot is obtained by evaluating the eigenvalues of this return-ratio matrix for s ∈ (−jω g /2, +jω g /2). Then, the system is stable if the number of counter-clockwise encirclements of the point (−1, 0) is equal to the number of righthand plane poles of Z −1 g Z c included within the LTP Nyquist  contour. Using (16), the poles of H 1 H −1 2 are evaluated using a contour plot. This is because an explicit expression for these matrices in terms of s is not easily obtainable analytically, but evaluation for single values of s is more feasible. Note that the contour plot has been used here only for the sake of theoretical validation of the method, and would not be feasible in practice, as discussed below. In fact, the contour plot allows one to chose a grid of values for s, such that s = p + jq, p, q ∈ R in the domain of interest and for each of these values H 1 H −1 2 is evaluated (recall equations (13), (14), with jω substituted by p + jq), the modulus of each element of the matrix is calculated and the largest of these is stored and used in the contour plot. The idea is that as s approaches the location of a pole, the modulus of one of the elements of the matrix H 1 H −1 2 will diverge and it will be shown in red-yellow colour in the contour plot, whereas if s is away from a pole, the greatest modulus will be small. Thus, with an appropriate set of values of s, the contour plot allows a good estimation of the pole locations, as shown in Fig.10. Since the return-ratio matrix Z −1 g Z c depends on the steady-state solution, which is a function of I ref , the poles in both stable and unstable cases are reported. Fig.11 shows a zoom of the region of interest in the Nyquist LTP contour. In the stable case, two poles are included inside this contour, at s = 0 and s = 903.3, while in the unstable case again two poles are included inside the contour, now at s = 0 and s = 1139.2, confirming the dependence upon I ref . Thus, in both cases, two poles are present inside the Nyquist LTP contour. This is a particularly interesting feature of the system under analysis. In fact, most of the impedance-based stability approaches take for granted that for a stable system, the open-loop transfer function does not have any poles with Re[s] ≥ 0. Hence the system is assumed to be stable if the Nyquist plot does not encircle the (−1, 0) point, and it is unstable when there are encirclements of this point. However, as demonstrated in this work, this assumption is not generally true, and in the system under analysis it is actually incorrect. Therefore, the identification of the number of poles with Re[s] ≥ 0 becomes a crucial task to perform in the stability analysis of a system using the impedance method. From a practical perspective, when the analytical model of the system is known one possibility is to apply the procedure based on the contour plot, as above discussed in detail, whereas when the analytical model is not known, it is possible to fit the measured Harmonic Impedances Z g and Z c using one of the several methods for curve fitting from Bode plots, in order to obtain an estimation of the transfer functions as function of s.
The Nyquist plots are reported in Fig.12, with the stable case shown in Fig.12(b), where it can be seen that there are two counter-clockwise encirclements of the point (−1, 0), ensuring that the system is stable; by contrast, the unstable case is depicted in Fig.12(c), where the external counter-clockwise encirclement is compensated by the internal clockwise encirclement, thus no net encirclements of the point (−1, 0) are present, leading to an unstable system. It is worth noting that the external encirclement is obtained with an infiniteclosure, which provides a counter-clockwise encirclement, in accordance with the theory for infinite-closures, [19], which states that when two segments of the Nyquist diagram are diverging towards plus and minus infinity, respectively, then these two segments must be joined by n semicircles (with infinite radius) in a counter-clockwise fashion, where n is the number of eigenvalues of the open-loop HTF in the origin of the complex plane. In Fig.10 it can be seen the presence of one pole in the origin, thus the Nyquist plot has an enclosure of one counter-clockwise semicircle.

VII. SWITCHING MODEL RESULTS
In this Section simulations are reported on the full switching model of the system. The control has been discretised and implemented in C-language and the circuit has been simulated in the Plecs-Matlab environment, including digital computation delays as is usual in real systems. The stability boundary is slightly different to the one evaluated based on the average model and I * ref lies in 11.4-11.5A. Such a difference, quantified to be less than 2%, is entirely expected, because of the linear modelling of the digital delays. The predicted stability boundaries have been validated through a time-domain simulation and implementation of an impedance measurement technique based on current injections to practically evaluate the LTP Nyquist Criterion with a "black-box" approach, more oriented to practical scenarios.

A. Time-Domain Simulation
Stable and unstable operation modes are reported in Fig.13, showing relevant currents and voltages. A good agreement with the time-domain simulations and with the analytical results based on the average model is obtained.

B. Current-injection measurement
In the switching model, the Harmonic Impedances can be obtained by injecting a small-perturbation current at the common-connection point and measuring the relative perturbation currents and voltages, as shown in Fig.2. From Section II, it is known that the LTP Nyquist contour for which it is relevant to calculate the LTP Nyquist plot is s ∈ [−jω g /2, +jω g /2], thus for Ω ∈ [−ω g /2, +ω g /2]. A succession of values for Ω is chosen and the LTP Nyquist plot is evaluated separately for each value. A total number of 2M + 1 separate current injections is required, which are of the form (M is the truncation order of the system): i 2 (t) =Ĩ 2 cos((Ω + (−M + 1)ω g )t) . . .
These grid small-signal currents are measured and an FFT is applied to extract the harmonic components, which are then stored in the following matrix: from which the Harmonic Impedances are obtained by inverting the respective current matrices in (49). This procedure must be repeated for each value of Ω ∈ [−ω g /2, +ω g /2], remembering that the LTP Nyquist plot for Ω ∈ [−ω g /2, 0] is the complex conjugate of the one for Ω ∈ [0, +ω g /2], so only the evaluation of the latter is required. Fig.14 shows the Nyquist plot obtained from the measured harmonic impedances compared with the one from analytical calculations (section VI-c), for the stable operation mode at I ref = 10A. The measured LTP Nyquist plot is obtained injecting 2 sets of perturbation currents: the first with f Ω = 1Hz, f Ω = Ω/(2π) (blue dots in Fig.14), and the second one with f Ω = 17Hz (red dots in Fig.14). The calculated LTP Nyquist plot is obtained by plotting the eigenvalue loci of the open-loop HTF Z g (jΩ) −1 Z c (jΩ), for Ω ∈ [1, 3, .., 49] ( (15) and Theorem of Section II.B), and is indicated by the black asterisk in Fig.14. The number of the set of injections can be decided according to the required resolution. In this case, two sets are enough to describe the LTP Nyquist plot around the critical point (-1,0). Good agreement is achieved, showing the potential of the current-injection method as a means of estimating harmonic impedances in practical applications.

VIII. CONCLUSION
In this paper two general methods are presented, based on LTP theory, to perform stability analysis of complex nonlinear time-periodic power systems. A case study of a singlephase network, comprising a voltage-source and a currentsource converter is considered for the development of the stability analysis. One method is based on the eigenvalueloci of the linearised system and requires a full knowledge of the model. The other, based on the LTP Nyquist Criterion, has the advantage that it allows the stability analysis to be performed using small-signal current injections: by measuring the relevant perturbed quantities, Harmonic Impedances can be calculated and stability assessed, considering the system as a black-box. Simulation results are provided based on both average and switching models of the system, showing a good accuracy in the identification of the stability boundaries.