A Potential Issue of Using the MIMO Nyquist Criterion in Impedance-Based Stability Analysis

The impedance-based stability criteria are widely adopted for small-signal stability analysis of power systems. Among them, the determinant-based multi-input multi-output (MIMO) Nyquist criterion is particularly attractive since it simplifies the stability determination process without calculating the eigenvalues of MIMO systems. However, it is found in this letter that the determinant-based Nyquist criterion can potentially result in incorrect stability analysis results with pure inductance models in the MIMO models of the power network. This letter illustrates the issue and its root cause through mathematical derivation and numerical analysis. Then, the issue is demonstrated through an example system, and the corresponding solutions are provided as well.

With the increasing penetration of renewables in power systems, the small-signal stability issue across a wide frequency range becomes a concern [1], [2]. The impedance-based stability criterion offers an effective way to analyze the smallsignal stability of such systems since no detailed internal information about system components is required [3].
For three-phase ac systems, the impedance-based approach deals with multi-input multi-output (MIMO) systems, and hence the generalized Nyquist stability criterion (GNC) is used to determine the stability of the system [4]. Two main approaches, the eigenvalue-based GNC and the determinantbased GNC, have been used for stability analysis. Although the eigenvalue-based GNC could provide more insight like gain margin, phase margin, and the approximate resonance frequency, it could be tedious to examine each characteristic locus when using the eigenvalue-based GNC for large-scale systems. For the determinant based GNC, only one Nyquist plot needs to be checked to determine the stability of the whole system. Thus, the determinant-based GNC is often preferred in large-scale systems [5], [6], [7].
However, in this letter, it is found that when applying the determinant-based GNC for small-signal stability analysis of power systems, pure inductance models in MIMO system models may lead to wrong prediction results for the system's stability in some cases. The pure inductance models in system models can induce imaginary axis poles at 60 Hz and thus lead to incorrect results in both stable and unstable cases when using the determinant based GNC. This issue has not been reported in the literature so far, and therefore, this letter intends to: 1) identify the issue induced by pure inductance models when applying the determinant-based GNC in practical applications; 2) illustrate the issue based on the equivalence of eigenvalue-based and determinant-based GNC; and 3) provide solutions for this issue.
The rest of this letter is organized as follows. Section II illustrates the issue of the pure inductance models when applying the determinant based GNC, together with the derivation of the cause of this issue. A case study is provided in Section III to validate such analysis, together with the solutions for this issue. Section IV summarizes the conclusions.

A. THE EIGENVALUE-AND THE DETERMINANT-BASED GNC
The equivalent MIMO feedback system of a power system is shown in Fig. 1 [5], where G cl (s) is the closed-loop referenceto-output matrix, G cd (s) is the closed-loop disturbance-tooutput matrix, and G nw (s) is the MIMO impedance matrix model of the connection network.
The open-loop transfer function L(s) and the returndifference matrix F(s) can be expressed as: (1) G cl (s) should be stable since all the converters are designed to be stable individually with ideal external conditions. Thus, the stability can be analyzed by applying the eigenvalue-or the determinant-based GNC to L(s) or F(s), respectively.
For the eigenvalue based GNC, the stability of the system can be analyzed by examining the characteristic loci of L(s) according to (3): where Z () and P() represent the number of RHP zeros and RHP poles, respectively. Theoretically, L(s) should not have RHP poles. Thus, the system is stable if and only if the characteristic loci of L(s) do not encircle (−1, j0). Similarly, for the determinant-based GNC, the stability of the system can be analyzed by examining the Nyquist or Bode plot of the determinant of F(s) according to (4): The system is stable if and only if the Nyquist plot of det(F ) do not encircle (0, j0), or the phase angle change of the full frequency range is 0 in the Bode plot of det(F ).
The equivalence of these two GNCs is due to [8]: where λ n 's are the eigenvalues of L(s).

B. ISSUE WITH PURE INDUCTANCE MODELS
The determinant based GNC, however, can have issues when applied to systems with pure inductance models. Although there is no pure inductance in real systems, an inductor or a reactor is often simplified as a pure inductance when modeling the systems in many cases. The dq impedance model of a pure inductance L considering the coupling terms is: where ω 60 denotes angular frequency corresponding to the system frequency, i.e., 60 Hz. The dq admittance model of L can be derived as: As seen from (7), the dq admittance model of pure inductance L has imaginary axis poles at ±60 Hz. The MIMO network impedance matrix G nw (s), which is derived based on admittance models of passive components in the connection network and its network topology, can thus have such imaginary axis poles that are induced by these pure inductance components in corresponding matrix elements. Note that the Kron reduction and G nw (s) reformation when deriving G nw (s) [5] can still preserve such imaginary axis poles in many cases. Thus, when calculating the return-ratio matrix L(s) as the multiplication of diagonal matrix G cd (s) and network impedance matrix G nw (s) as in (1), the imaginary axis poles could be induced in L(s) because of these pure inductance components in the system, as will be shown later in an example.
In the GNC, the "Nyquist plot" means "the image as s goes clockwise around the Nyquist D-contour". When there are imaginary axis poles in the system, D-contour must avoid locations where L(s) has imaginary axis poles by making small indentations (semi-circles) around these points [9]. As pure inductances induce imaginary axis poles at ±60 Hz, the D-contour of the system should be like in Fig. 2. Finding of all eigenvalues of L(s) with frequency s varying along Nyquist D-contour in the complex plane becomes a major interest in coping with the stability problem for multi-variable systems.
To apply the MIMO Nyquist criterion, the eigenvalue loci should be calculated as a function of frequency. The eigenvalues can be computed by conventional techniques for a given frequency, and theoretically, one would have to repeat the entire computational procedure for all frequencies. However, it is impractical for such repetition in the actual stability analysis. For both eigenvalue-and determinant-based GNC, a more practical way of using software (e.g., MATLAB) is to calculate the eigenvalues at discrete frequencies in a wide frequency range, and then these points are connected in the software. An example case is used here to illustrate the impact of having imaginary axis poles when applying the MIMO Nyquist criterion, without loss of generality.

1) THE EIGENVALUE BASED GNC WITH IMAGINARY AXIS POLES
For the eigenvalue based GNC, assume there are two eigenvalues λ 1 (s) and λ 2 (s) for L(s): When L(s) has imaginary axis poles at 60 Hz due to pure inductance components, these eigenvalues become very large near 60 Hz in practical implementation when using software like MATLAB. Without loss of generality, assume the polarities of these eigenvalues are that the real part a 1 (s) and a 2 (s) of the eigenvalues are positive; the imaginary parts b 1 (s) and b 2 (s), quickly change from +∞ to −∞ at 60 − Hz and 60 + Hz. The characteristic loci of L(s) look like Fig. 3. Note that only the characteristic loci of the positive frequency range are shown in this figure. The characteristic loci of the negative frequency should be symmetrical to the positive frequency.

2) THE DETERMINANT BASED GNC WITH IMAGINARY AXIS POLES
According to (5), for the determinant based GNC, the determinant should be: Then, at 60 − Hz, the real part of the determinant is negative, and the imaginary part of the determinant is positive; at 60 + Hz, the real part is negative, and the imaginary part is also negative. When the points at 60 − Hz and 60 + Hz are connected in MATLAB, and the Nyquist plot of the determinant will look like Fig. 4. Similar to Fig. 3, only the positive frequency range of the Nyquist plot is shown here. Corresponding Bode plot of det(F ) will have a 360 • phase change in Bode plot since the Nyquist plot encircles the origin point, as shown in Fig. 4.  When applying the determinant-based GNC to the Nyquist or the Bode plots shown in Fig. 4, it shows that the system is unstable since the Nyquist plots encircle the origin point twice (considering the full frequency range from −∞ to +∞), and the phase angle change of the Bode plot in the full positive frequency range is 360 • . However, the eigenvalue based GNC shows that the system is stable since (−1, j0) is not encircled. Similarly, when in an unstable case, the analysis result predicted by the determinant-based GNC could be stable because of the wrong 360 • phase change, as shown in Case II of an example system in Section III. Thus, the stability prediction using the determinant-based GNC is in contradiction to the eigenvalue-based GNC in such cases which could provide wrong analytical results, and it is due to the imaginary axis poles induced by the pure inductances in the system.

III. CASE STUDY AND SOLUTIONS
To validate the above analysis, a simple case study is built in PSCAD simulation, and the analysis code using determinantand eigenvalue-based GNC is developed in MATLAB. Moreover, solutions for the issue when using the determinant based GNC are provided in this section.

A. SIMULATION INVESTIGATION
A 5-bus system shown in Fig. 5 is built in PSCAD as an example for the analysis of general multi-bus systems using the MIMO Nyquist criterion. For simplicity, bus 3 is connected with a grid-following (GFL) converter, while bus 1 and bus 2 are connected to ideal voltage sources. Parameters of the passive components are shown in Tables 1 and 2. The example system is stable in Case I, as can be seen from the steady-state voltage and current waveforms shown in Fig. 6. In Case II, the real power of the load at bus 3 is reduced to 200 MW, which leads to the unstable phenomenon in Case II, as the voltage and current waveforms at bus 3 shown in Fig. 7.
To apply the impedance-based stability analysis, the impedance or admittance models of all system components   should be derived first. The admittance models of the constant impedance loads, which are necessary to derive the connection network matrix G nw (s) in Fig. 1, are modeled as single S loads or RL loads in parallel connection depend on the parameters in TABLE 2. Note that when a load is modeled as RL in parallel connection, the admittance model of the load is calculated as the admittance of a pure resistance plus the admittance of a pure inductance. Therefore, the loads modeled as RL in parallel could also cause the imaginary-axis poles problem, and it can potentially lead to incorrect stability analysis results when applying the determinant-based GNC.
To analyze the small-signal stability of the 5-bus system, the closed-loop disturbance-to-output transfer function matrix and the impedance matrix model of the connection network are developed based on the system topology. The determinantbased GNC and the eigenvalue-based GNC are then applied to the return-difference matrix F(s) and the open-loop transfer function L(s), respectively. 1000 data points are calculated in the frequency range of 0.1 Hz to 10 5 Hz.    Fig. 8 shows the Nyquist plot of det(F ) in Case I, and Fig. 9 shows its Bode plot. Fig. 10 shows the characteristics loci of L(s). As shown in Fig. 9, when plotting the Nyquist plot, MATLAB will connect the two points at 59.8 Hz and 60.2 Hz since there are no other data points between the two frequencies. This leads to the encirclement of the origin point  (0, j0) in the Nyquist plot. Note that such encirclement is not due to the number of data points that were taken for this analysis. Even with more data points, due to the imaginary axis poles induced by pure inductance, the connection of the two frequency points that are closest to 60 Hz in the software will still encircle the origin point (0, j0) as analyzed above.

1) CASE I
Correspondingly, the Bode plot in Fig. 9 shows a 360 • increment in the positive frequency range. Both the Nyquist plot and Bode plot of det(F ) predict that the system should be unstable, which contradicts the stable simulation results shown in Fig. 6. As analyzed earlier, the mismatch between the analysis results and the simulation is due to the imaginary axis poles induced by the pure inductance models in the constant impedance loads. On the other hand, when judging the system's small-signal stability with eigenvalue-based GNC from Fig. 10, it appears that the system is stable in this case. The connection of two points at 59.8 Hz and 60.2 Hz do not encircle the critical point (−1, j0). Therefore, the issue of the imaginary axis poles does not affect the results using the eigenvalue based GNC in this case.  connects the two points at 59.8 Hz and 60.2 Hz, which results in no encirclement of the origin point (0, j0) in the Nyquist plot and zero phase change in the corresponding Bode plot. So, both the Nyquist and Bode plots using the determinantbased GNC predict that the system should be stable in this case, which contradicts the unstable simulation results shown in Fig. 7. In this case, the determinant-based GNC provides wrong analytical results when the system is unstable, and that is also due to the pure inductance models. When looking at the characteristic loci of L(s) using the eigenvalue-based GNC shown in Fig. 13, the system is predicted to be unstable, as the critical point (−1, j0) is encircled in this case. This is caused by the similar reason as in Case I, that the connection of the two points at 59.8 Hz and 60.2 Hz in the characteristic loci does not affect the encirclement of (−1, j0).

B. SOLUTIONS
In some system models, there are components (e.g., transmission lines, transformers, or loads) that have only L parameters in their models (as in Table 2). To solve the abovementioned issue, one practical way is to shift all imaginary axis poles into the LHP, for example, by replacing the 1/s with 1/(s+ ∈), where ∈ is a small positive number [9]. To correct the smallsignal stability analysis results of the example 5-bus system using the determinant-based GNC, a small ∈ is added for all pure inductance components. The updated stability analytical results of this example system with such modifications are shown in Figs. 14 and 15 for Case I and Case II, respectively. Analytical results are now the same as the simulations with this change. For Bode plots in both Figs. 14 and 15, the blue curves are the original results using the determinant-based GNC, while the orange curves are the results with such modifications using this solution. One can see from the Bode plots in both cases that the magnitude and phase plots show that such modifications almost don't change the characteristics of the system at other frequencies except around 60 Hz (where the pure inductance models cause the problem); moreover, such small modifications correct the error to provide the same analytical results as in simulations.  The second solution that do not change the original model is to modify the determinant based GNC by compensating the Nyquist and Bode plots of det(F ) with phase compensation. From the above analysis, it is concluded that the incorrect analytical results are due to the wrong phase change induced by the imaginary axis poles at 60 Hz. So, when such imaginary axis poles are found in systems, the Nyquist plot and the Bode plot should add a corresponding phase compensation by following the D-contour in Fig. 2 to get correct analytical results. In the example system, for Case I showed in Fig. 16, det(F ) should go to infinity from 59.8 Hz and then goes back to 60.2 Hz in the clockwise direction. With such compensation, in Case I, the Nyquist plot of det(F ) should not encircle (0, j0), which represents the system is stable in Case I. And correspondingly, the Bode plot of det(F ) should also have a 360 • phase compensation after 60 Hz, as shown in Fig. 16. The 0 • phase change in the Bode plot after the compensation is also showing that the system is stable in this Case. Similarly, in Case II, the det(F ) also goes to infinity from 59.8 Hz and then back to 60.2 Hz in the clockwise direction, as shown in the Nyquist plot in Fig. 17. One can see that the det(F ) encircles (0, j0) twice in the Nyquist plot, which means the system is unstable in Case II. Correspondingly, the phase change in the Bode plot in Fig. 17 also decreases 360 • . It also represents that the system is unstable in this case, making it the same as the simulation results.

IV. CONCLUSION
This letter identifies the issue caused by pure inductance models when applying the determinant-based GNC for smallsignal stability analysis. It is found that the pure inductance models in the system can induce the imaginary axis poles to the MIMO system model, and such an issue can potentially result in incorrect stability analysis results. A case study of an example system validates the analysis. Solutions to this issue include shifting the imaginary axis poles into the LHP or providing phase compensation for the determinant based GNC when the imaginary axis poles are involved due to pure inductance models.