Stability analysis of single-phase grid-feeding inverters with PLL using Harmonic Linearisation and Linear Time Periodic (LTP) theory

A general method for the stability analysis of power converters is presented in this paper. The method is based on Harmonic Linearisation and Linear Time Periodic (LTP) analysis techniques and a single-phase grid-feeding inverter with PLL is considered as case study. Although stability analysis has been developed using the average model of the converter, the obtained results can be extended to the switching model and it is possible to evaluate precisely the boundary between stability and instability.


I. INTRODUCTION
AC power systems based on power converters are complex and highly non-linear because both the converters and the interconnected loads often show a non-linear behaviour. A theory for the stability analysis of such systems is required and existing analysis techniques can be broadly categorized into frequency and time-domain methods. Frequency-domain methods are based on the impedance models of source and loads of the system and on the application of the Nyquist stability criterion [1]. For balanced and symmetric three-phase AC systems the analysis can be performed in the dq reference frame [2], but this is no longer the case for unbalanced or single-phase systems. In such cases the dynamic phasor approach allows one to determine the 2-dimensional source and load impedance, however a precise boundary between stable and unstable systems is not usually provided [3], [4]. Harmonic linearisation is a technique that develops a small-signal linear impedance model for a non-linear system along a periodically time-varying trajectory [5], [6]. These calculated small signal-impedances can be used along with the Nyquist stability criterion but potential limitations might affect this method, as stated in [1]. In this paper a time-domain method is developed and discussed. Starting from a non-linear state space model of the system, harmonic linearisation is applied in order to derive a linearised model, which is an LTP system. Subsequently LTP time-domain theory is applied in order to perform stability analysis [7], [8]; alternatively it is possible to obtain equivalent results in the frequency domain by defining a Harmonic State Space (HSS) small signal model [7], [9]. This method allows us to define accurately the stability boundaries of the system and no approximations are introduced apart from a truncation of the Fourier coefficients, which can be chosen according to the desired accuracy of the results. The single-phase inverter is supplied by an ideal DC source, V dc , and its output is connected through an L 2 C 1 L 1 low-pass filter to the grid, which is represented by an ideal sinusoidal voltage source, V g (t), in series with an R g L g impedance. Unity power factor operation is considered. A Phase Locked Loop (PLL) is used to measure the phase of the grid voltage and to generate the current reference signal for the inverter output, I ref (t). A PI control is then used to control the inverter. We will show that there is a threshold value for the parameter I ref , above which the system becomes unstable. That is, the instability is due to the fact that the PLL is no longer able to generate the correct phase reference for the controller. It is worth noticing that a damping resistor R c is used in the low pass filter in order to simplify the design of the current control and focus on the instability caused by the PLL. Using LTP techniques and theory we will also show that it is possible to evaluate precisely the instability threshold and perform a rigorous stability analysis. The analysis will be focused on the average model of the system. This model gives a good description of the switching system up to one third of the switching frequency. Since it is of interest to study the case where the instability is caused by the incorrect behaviour of the PLL, we know that such instability arises at frequencies far below one third of the switching one. So the switching model is replaced by the average one and the analysis is based on this approximation. Finally the results of the analysis and the time domain simulation of the switching model are compared and it is shown that both results agree.

MODEL
For single-phase applications the PLL needs to estimate the inquadrature component of the input voltage, V o (t), and among the several available topologies a linear filter that introduces a T g /4 delay at ω g is used: D(s) = ω 2 g /(s 2 + sω g + ω 2 g ). The whole system is equivalently represented by the state space model (1), which is an eighth order Non-Linear Time Periodic (NLTP) system, with all the state space variables T g -periodic and the non linearity due to the presence of the PLL. The only input signal is V g (t).

IV. LINEAR TIME PERIODIC APPROACH
The stability analysis that we are going to describe involves an LTP system. However, most real AC system average models are NLTP, as is the one we want to analyse (1). Then the first step to perform is to linearise the system. In DC systems, the linearisation is applied around a constant steady state operating point, while in AC systems the steady state operating point is time-periodic rather than constant. So in this case traditional linearisation techniques cannot be applied, but Harmonic Linearisation must ne applied , which means a linearisation around the steady state operating trajectory.
Given a general NLTP system, T -periodic: and given a steady state inputū(t) we solve numerically (in Matlab for example) this NLTP system and we obtain the steady state solutionx(t). Now Harmonic Linearisation is applied, so a small signal perturbation is added to the steady state input, output and solution: and making this substitution in the NLTP system: Taking into account only the first order terms, those proportional tox(t),ỹ(t) orũ(t), and ignoring all the others gives us the linearised model, which results in an LTP system: with the matrices A(t), B(t), C(t) and D(t) being T -periodic. In Fig.2 the steps to obtain the LTP model are summarised.

Fig. 2: From NLTP to LTP
Before proceeding with the stability analysis we want first to review some properties of Linear Time Invariant (LTI) systems which will be useful in the analysis of LTP systems. Consider the LTI system:ẋ The test signal of fundamental interest in the study of LTI systems is the complex exponential u(t) = u 0 e jΩt = u 0 e st , s ∈ C, u 0 ∈ C m . Substituting this signal in the LTI state space model and evaluating the steady state output response gives us: So the transfer function operator G(s) is defined and this is possible only because both input and output spaces are the same; the input signal has one harmonic component and the output signal has the same harmonic component, with possibly different amplitude and phase. Now the stability analysis can be performed by evaluating the eigenvalues of the A matrix.
We are now interested in obtaining the transfer function operator for the LTP system. Consider first the case where the same signal used for LTI systems is now used for LTP systems. The important result that comes out is that the output space is different from the input space, as shown in Fig.3. In this case the output signal is composed of the same input harmonic Ω plus infinite harmonics that are Ω + qω T , with q ∈ Z and ω T = 2π/T = 1/f T . So the main result is that with a complex exponential as test signal to an LTP system we are unable to define the transfer function operator and proceed with stability analysis. A different test signal is now defined, which will be referred to as Exponentially Modulated Periodic (EMP) signal [7], [8], When an EMP signal is injected as test signal to an LTP system, we find the significant result that both input and output spaces are the same and they possess the same set of harmonic components, as shown in Fig.4. So now it is possible to define the transfer function operator and apply stability analysis.
Consider the LTP model (5). The dynamic matrix is expanded in complex Fourier series as: and similarly for B(t), C(t) and D(t). An EMP test signal implies that the steady state and output responses are also EMP signals: We now make these substitutions in (5), which gives us: This system of equations has to hold for ∀n ∈ Z, which means: This is a concise representation of the input-output relationship between the Fourier coefficients of the input and output signals. However the manipulation of Fourier series usually leads to complicated calculation and for this reason the Toeplitz transform is introduced in order to make the analysis simpler.
A Toeplitz transformation is defined as follows:

is a doubly infinite block Toeplitz matrix and the matrices A i are the Fourier matrix coefficients of the Tperiodic matrix A(t). A similar definition applies to the matrices T [B(t)] = B, T [C(t)] = C, T [D(t)] = D and to the vectors T [x(t)] = X , T [u(t)] = U, T [y(t)] = Y.
If we examine the system of equations (15) for ∀n, we notice that the Toeplitz transforms can be used in order to obtain a clearer and more compact notation. So we are able to define the harmonic state space model (HSSM) of the LTP system: with N = blkdiag{jnω T I}. With simple steps the Harmonic Transfer Function (HTF) of the system is defined: Stability analysis is now performed 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. So far an infinite number of coefficients has been considered in the Fourier series expansion. To implement this analysis on a computer we introduce a truncation order N which refers to the number of harmonics taken into account. For N = 2, for example, this means that we are taking into account the DC-component and the first and second harmonic.
The corresponding Fourier expansion involves the Fourier coefficients for n = −2, −1, 0, 1, 2: and we consider the following associated truncated Toeplitz form: with Z a zero matrix of proper dimension. So when the truncation order is increased a larger number of Fourier coefficients of the Fourier expansion is taken into account and at the same time the dimension of the associated Toeplitz form increases . Given an LTP system of order p, which means there are p state space variables, and a truncation order N , the number of eigenvalues associated with the matrix A − N will be (2N + 1) × p. However, only p of these eigenvalues are the important ones for the analysis; all the others are translated copies of the original ones, with translation equal to jnω T , n = ±1, . . . , ±N . Fig.5 shows an example of LTP pole loci with q = 6 and N = 1. In red are depicted the important poles and in black their translated copies. It can be noticed that for a large truncation order the eigenvalue loci result in long vertical lines of eigenvalues. In order to apply Harmonic Linearisation the steady state solutions of the system (1) must be evaluated first. We will discuss some mathematical considerations and then solve numerically the system equations. The set of steady state solutions will be of the form: x i (t) = |x i | cos(ω g t + arg(x i )) =x i e jωgt + c.c. 2 for i = 1, 2, 5, 6, 7, 8 ;x 3 (t) = ω g t +x 03 (21) Substituting these expressions in the non-linear state space model gives us: Ignoring the complex conjugate parts and making some simplifications gives: (24) A similar approach is applied to the other equations: These 6 equations are now solved numerically in Matlab and the solutionsx 5 ,x 6 ,x 7 ,x 8 ,V o ,V conv are obtained. We find: so the solutionsx 1 ,x 2 are given by: The last two quantities to be defined arex 3 andx 4 : Applying trigonometric simplification gives us: So it follows thatẋ 4 (t) = 0. But from the state space model, using again trigonometric simplifications, we have: and soẋ which impliesV o −x 03 = 0 or ±π. In our case,V o =x 03 , which gives the last two solutions: VI. LINEARISED MODEL Following the steps illustrated in the previous paragraph we now derive the LTP Model of the NLTP system (1), which is of the formẋ(t) = A(t)x(t), with A(t) being a T g -periodic matrix (42). It is now possible to derive the Toeplitz form A and evaluate the eigenvalues of A − N to determine whether the system is stable or not.

VII. SIMULATION RESULTS
In this section simulation results are presented to validate the proposed stability analysis method. From the study of the eigenvalue loci it can be seen that the maximum value for the reference input I ref for which the system is stable is around 22.5 A. Two simulations have then been performed: the first with I ref = 22 A, which gives a stable system and the second with I ref = 23 A, for the unstable case. In Fig.6 are shown the eigenvalue loci for the two cases: all the eigenvalues are in the left half-plane for the stable case, while for the unstable case some of them are in the right-half plane. Some spurious eigenvalues arise due to the fact that a truncation is introduced in evaluating the Toeplitz form (in this case a truncation order N = 40 is considered) but they have no physical meaning. In Fig.7 we show the time evolution of the relevant currents, voltages and phases of the switching model in the two cases. It is worth noticing that in Fig.7(d) the inductor current I l (t) correctly tracks the reference current I ref (t), which confirms that the instability is not due to the current control but to the PLL. A good accuracy is provided in finding the stability boundaries of the system.
8 (t) − k p2 cos(x 3 (t))V o (t)x 3 (t) + k p2 cos(x 3 (t))x 1 (t) − k p2 sin(x 3 (t))x 1 (t)x 3 (t) x 4 (t) = −k i2 sin(x 3 (t)) L 1 R g − L g R c L g + L 1 x 6 (t) − k i2 sin(x 3 (t)) L g R c L g + L 1x 7 (t) − k i2 sin(x 3 (t)) L g L g + L 1x 8 (t) − k i2 cos(x 3 (t))V o (t)x 3 (t) + k i2 cos(x 3 (t))x 1 (t) − k i2 sin(x 3 (t))x 1 (t)x 3 (t) x 5 (t) = −I ref sin(x 3 (t))x 3 (t) −x 7 (t) x 6 (t) = − R c + R g L g + L 1x 6 (t) + VIII. CONCLUSION In this paper a general method is presented, based on harmonic linearisation and LTP theory, to perform stability analysis of complex non-linear time periodic power systems. A case study of a single-phase grid-feeding inverter with PLL has been used to show the practical application of the method. Numerical simulations have been performed to validate the proposed technique, showing good accuracy in the identification of the instability boundary.