Extended Nodal Admittance Matrix Based Stability Analysis of HVDC Connected AC Grids

Motivated by emerging oscillation issues of high-voltage direct-current (HVDC) transmission system resulted from the interaction with AC grids, impedance-based stability analysis becomes an attractive approach in the de-risking studies due to the capability of black-box modelling and clear physical meaning. However, impedance-based stability assessment of HVDC is typically conducted either from the AC side or from the DC side, without fully considering both AC and DC grid dynamics, thus impairs evaluation accuracy. In this paper, the nodal admittance matrix based resonance mode analysis method is extended to include both AC and DC dynamics through formulating a hybrid AC-DC admittance matrix. To achieve this goal, HVDC interfacing-converters are modelled as three-port admittance networks with one DC port and two AC ports (in dq-frame) integrating both the AC side frequency coupling and the AC-DC dynamic coupling. In contrast to the widely adopted Nyquist-based methods, the proposed method does not requires knowing the number of right half-plane poles of system minor-loop gain, thus simplifies stability analysis, and moreover, the participation factors of both AC and DC grid nodes in critical resonance modes can be identified through eigenvalue-decomposition analysis, which is in favor of the design of mitigation measures. The effectiveness of the proposed method is validated by EMT-simulations in MATLAB/Simulink.


I. INTRODUCTION
In order to achieve massive grid integration of renewable energy sources (RES) and long distance power transmission, multiple HVDC links have been constructed and integrated into German grids [1], [2]. Since both HVDC and RES are converter-interfaced components, the fast growth of such power electronics devices in modern power systems have resulted new dynamics and stability issues, and the most commonly encountered one is the wideband resonance stability [3]- [5]. Over the last decade, a series of oscillation incidents involving HVDC have been reported [6], [7]. The recorded oscillation frequencies range from several Hz to above 1 kHz. Typically they are triggered by the change of operating conditions in involved systems [3], [4], [8], such as the grid topology change induced by grid disturbances.
To reveal the mechanisms of such issues and prevent them happening again, e.g. through adopting proper measures in power system planning or operation stages, quite a few studies have been conducted, either through electromagnetic transient (EMT) analysis or analytical stability assessment [4], [7], [8]. Although EMT analysis can usually give plausible stability evaluation results through applying detailed nonlinear grid models, EMT simulation is in most cases time-consuming and incapable of capturing the underlying mechanism of oscillation phenomenon, thus will mainly be used for validation in the paper [9]. As for the analytical stability assessment, several widely used methods are: state-space based method, Nyquist based method and nodal-based Resonance Mode Analysis (RMA) method [5], [10]. For easy understanding, the nodalbased RMA method will be named as Ybus based method. Both Nyquist and Ybus based methods are impedance methods, which are more popular than state-space based method due to the strengths in natural association with physical circuits and the capability of black-box modelling. However, their application in the stability assessment of This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3177232, IEEE Access 2 VOLUME XX, 2017 HVDC and the resulted hybrid AC/DC grids is typically conducted either in AC side or DC side while partly neglecting or simplifying the dynamics of AC and DC side grids [11]- [14], which may fail to identify the potential arising stability problems when the AC-DC dynamic coupling is strong.
In order to get more complete model for hybrid AC/DC grids including detailed AC and DC grid dynamics that have been partly neglected or simplified in the past, multiport admittance model of HVDC converter that reserves the AC-DC dynamic coupling has been developed [15], [16] and furtherly applied in enhanced Nyquist based stability analysis [17], [18]. However, these methods fail to identify the critical dynamic modes and the participation factors (PFs) of grid components, which hinders the identification of risk sources. This issue was properly treated in the Ybus based method, which uses nodal admittance matrix for stability assessment [19]. But Ybus based method is typically applied to single sequence component system, which neglects the coupling between positive-and negativesequence impedances [20]. Moreover, the integration of AC-DC dynamic coupling into Ybus based method has not been found in literatures so far. To bridge the above gaps, this paper extends the traditional Ybus based method to integrate full dynamic couplings and analyzes various stability impacts in hybrid AC/DC grids, simultaneously contributing to accurate stability assessment and better support in mitigation design.
The rest of this paper is organized as follows: Section II compares existing stability analysis methods and introduces the proposed enhancement. Section III presents the test grid including the detailed nonlinear model and the small-signal modelling. Section IV shows the effectiveness of the proposed method through case studies. Section V concludes the work. Figure 1. (a) shows a simple hybrid AC-DC grid with embedded HVDC, synchronous machine (SM) and renewable energy source (RES). Two-level voltage-source converter (VSC) based HVDC is adopted in the study. The HVDC is parallel with a high-voltage AC (HVAC) branch. Figure 1. (b) shows the RES and the feeding route, e.g. for representing the application of wind park integration. An aggregated single-converter model was adopted. Multiple converter-interfaced components and controllers are involved in the system, with the potential for stability-critical interactions. For the stability check of such a test grid, the most widely used analytical methods including the state-space based method, Nyquist based method and Ybus based method are compared in (b) Aggregated RES and the feeding route  For the state-space based method, it requires detailed information on the circuit, control and parameters of a system to build system state-space model, which is naturally incapable of black-box modelling and lack of scalability [5], [21]. While for the Nyquist based methods, it is shown in [10] that bode plots can not predict system stability accurately under certain conditions and Nyquist diagrams are not straightforward in checking the right-half plane (RHP) poles and the encirclements of (-1, 0). Moreover, when checking the stability of strongly meshed AC/DC grids using Nyquist based methods, it is difficult to partition the system into source-load subsystems at any point and it may require the stability check at multiple partitioning points [14], which weakens the suitability of the method for complex systems. As for the Ybus based method, it has the best scalability but lacks the integration of asymmetrical control coupling and AC-DC dynamic coupling regarding the state-of-the-art. To simultaneously achieve accurate stability check and stability risk localization for a hybrid AC/DC grid, the Ybus based method will be enhanced to eliminate its limitations as listed in TABLE I.

II. LIMITATIONS IN EXISTING STABILITY ANALYSIS METHODS AND THE PROPOSED ENHANCEMENT
As the basis for enhancement, the nodal admittance matrix of the test system should be formulated, and before that, all components should be modelled as parallel or series admittance networks. To integrate the dynamic coupling induced by asymmetrical converter control in admittance modelling, the AC side of grid components should include positive-and negative-sequence ports when modelled in the natural reference frame, or d-and q-axis ports when modelled under dq-frame [7]. The sequence ports and dq ports can be equivalently transformed to each other through linear transformation and frequency shifting [7]. For the convenience of determining the PFs of d-and q-control loops in the critical dynamic modes, the dq-frame modelling is selected. With the objective of integrating AC-DC dynamic coupling in admittance modelling, each HVDC converter is modelled as a three port hybrid AC/DC admittance network, with two ports at AC side and one port at DC side [17]. Then the SM, RES, AC filters and external grid (as parallel AC grid components) are modelled as two port admittance network. The AC lines and transformers (as series AC grid components) are modelled as four port admittance network. The HVDC cable (as a series DC grid component) is modelled as a two port admittance network.
Since the admittance modelling of active components including SM, RES and HVDC converter is locally linearized, which depends on the angle of the local reference frame, the models developed under local reference frames will be rotated to a defined global reference frame [14]. Note that the modelling of DC components and DC ports of hybrid AC/DC components does not involve dq-frame rotations. In addition, the stationary frame in DC side can be mapped to the dq-frame of AC side without frequency shifting. On this basis, all components can be directly connected in circuits for system level analysis. Nodal admittance matrix of the studied system is formulated in (1) (1) In comparison to the commonly used nodal admittance matrix that is formulated using single phase or sequence component variables, the nodal current of each AC node in (1), e.g. AC1 i , is a column vector comprised of d-and q-axis currents under global reference frame, e.g. as denoted by d1 I and q1 I . Each AC admittance element in the nodal admittance matrix is a 2-order matrix instead of a single quantity. The integration of DC grid nodes in the node-voltage equation gives a hybrid AC-DC nodal admittance matrix. Having K L and K T as the left and right eigenvector matrices of KK Y obtained through matrix decomposition [22], the eigenvalue admittance matrix M Y of KK Y is also obtained, as illustrated in Figure 2. Through setting Performing frequency scanning to the frequency-dependent modal impedances m1 () Zs , i.e. the diagonal terms of M Z , the critical dynamic modes can be determined by the peaks of the impedance amplitude versus frequency curves, which is the basis concept of resonance mode analysis (RMA) [19]. The frequencies at curve peaks indicate the frequencies of dynamic modes. The stability of the k th dynamic mode can firstly be checked by the real part value of the corresponding complex-valued modal impedance at the oscillation  (4) where Qk denotes the quality factor of k th resonance circuit, as defined by the frequency-to-bandwidth ratio of the resonance circuit, and k f  denotes the resonance bandwidth at half maximum, which is the difference between the two neighbored frequencies (of k f ) at which the impedance amplitudes are equal to half of the peak amplitude value (at k f ) [23]. The other common nearly equivalent definition for Q is the ratio of the energy stored in the oscillating resonance circuit to the energy dissipated per cycle by damping processes [23]. After identifying critical dynamic modes, the system can be transformed from modal coordinate system back to nodal coordinate system, as illustrated in Figure 2. The eigenvector matrix K T or K L can be used to determine the PFs of both AC and DC grid nodes for the dynamic modes.
An alternative way of dynamic mode identification and stability assessment is to evaluate the poles of the closed loop system [24], [25]  In nodal coordinate system where k  is the attenuation factor of the mode and Note that the real part () k Rf of modal impedance in (3), in which the Laplace operator s is an imaginary value, and the attenuation factor k  of the closed system pole in (6), in which s is a complex value, are complete different quantities, which should not be mixed up in the formulation of damping ratio.
Moreover, finding the poles of closed loop system transfer function through solving the roots of is given in the form of continuous function of frequency [26], which is not the case as illustrated in Appendixes A-C, thus the point-by-point numerical manipulation over frequency and the approximation, e.g. by vector fitting [27], on each set of frequency response by a continuous function of frequency is needed, which is a major disadvantage of this approach. Furthermore, due to the better visualization and physical meaning of impedance versus frequency curves in dynamic mode identification than the poles of closed-loop transfer function, only the (modal impedance) frequency scanning method as illustrated in Figure 2 and Eq. (2)-(4) will be considered for further analysis.

III. GRID MODEL
The 8-bus grid in Figure 1 is adopted for analysis. The component models and their main parameters are given in TABLE II. For the SM, the detailed round rotor synchronous machine model from MATLAB / Simscape library is adopted [28]. For the RES, the virtual synchronous machine (VSM) control, as one of the grid-forming (GFM) control methods, is adopted, as shown in Figure 3. The readers are referred to [29] for the details on the design of the VSM. For the HVDC, the generally used grid-following (GFL) control in real applications is adopted. Detailed circuit and control diagrams of the HVDC converters are shown in Figure 4.
GFL converters use phase-locked loop (PLL) to synchronize with grid, while the GFM converter automatically synchronize with grid through the swing blocks in VSM core, as shown in Figure 3 and Figure 4 individually. For the transformers, the typically used T model is considered, the series branch is described by the short-circuit voltage expressed in percentage (ur as resistive part and ux as inductive part), and the shunt magnetizing branch is described by paralleled resistance and inductance, both of which are set as 500 pu [7].
With the above configuration, the interaction of GFL converter, GFM converter and SM is involved in the test system. To prepare for the stability analysis, the dq-frame admittance models of SM and RES are derived referring to [30], as presented in Appendix A and Appendix B. Similarly, the three-port admittance modelling of HVDC converters through combining the control and power stage equations (including the AC-DC power balance) is conducted and validated through EMT simulations, as presented in Appendix C. Note that two DC voltage compensation modes for the modulation of HVDC converters are commonly seen in literatures [18], [31], as distinguished by constant DC voltage Udc0 compensation and instantaneous DC voltage udc compensation, as shown in in Figure 4. Both modes will be considered to show their impact on stability assessment.

A. STABILITY IMPACT OF WEAK GRID CONNECTION
Firstly, the simplified grid as shown in Figure 5 is adopted for analysis. The AC Grid-A and AC Grid-B are both modelled as ideal voltage sources with series RL impedances, and they are characterized by the short circuit ratios of SCR-A and SCR-B. The X/R ratio of their RL impedances is fixed to 50. The left side HVDC converter is operated in the active power control mode, and its setpoint is set as P=0.

B. STABILITY IMPACT OF DC VOLTAGE COMPENSATION MODE
DC voltage compensation using either the constant value Udc0 or the instantaneous value udc is widely used in converter modulation [18], [31]. For the PQ-controlled HVDC converter (VSC-A), replacing instantaneous udc compensation with constant Udc0 compensation mainly influences the 7 Hz resonance mode, as indicated by the modal impedances shown in Figure 6. For the Udc0 compensation mode, the damping of the 7 Hz resonance (red dashed line) is positive and the damping ratio is around 1 Hz / (2 7 Hz) 0.07 k     according to (4), which is slightly larger than the damping ratio of the udc compensation mode around 0.05, as calculated in Section IV.C. For simplicity, this statement can be qualitatively obtained by looking at the sharpness of the amplitude curve peaks. A less sharped amplitude peak indicates larger damping ratio, which can be validated by comparing the EMT simulation results in Figure 7

C. Stability Impact of the parallel AC line for HVDC
As already shown in above subsections, the simplified test grid in Figure 5 is stable at SCR-A=4 and SCR-B=4, and becomes unstable when SCR-B is reduced to 2. If increasing SCR-A to 6 and keeping SCR-B as 2, the overall SCR of the HVDC connections is still 8. Under this condition, it will be checked if the 18 Hz resonance mode can be stabilized by connecting the parallel AC line for HVDC. Then the new test grid is illustrated by Figure 8.  Figure 9 shows the modal impedances considering the following Scenarios: (a) no AC line connection; (b) single AC line connection; (c) double AC line connection. All the scenarios have SCR-A=6, SCR-B=2 and Udc0 compensation mode for VSC-A. Note that all the SCR values in the paper refer to the HVDC capacity. In Scenario (a), the 18 Hz resonance mode is still unstable due to the negative damping as indicated in the impedance real part plot, while in Scenarios (b) and (c), the system is stabilized by connecting the parallel AC line(s). In comparison to Scenario (c), the resonance peak in Scenario (b) is much sharper, which indicates a smaller damping ratio according to (4). The 7 Hz resonance mode observed in Section IV.A and IV.B disappears because the AC Grid-A is strengthened. Above analytical results are then validated by the EMT simulation results as shown in Figure 10  coupling in HVDC is weakened by the parallel AC line in this case.

D. STABILITY IMPACT OF SM AND RES
In above analysis, the left and right side external grids are both modelled as Thevenin circuits. To create a more practical impedance profile, the AC Grid-A is replaced by SM and RES (see Figure 1). And the influence of replacing the Thevenin circuit (with SCR-A=2.8) with SM and RES (in total of 2.8 GVA capacity) is analyzed for the scenarios of varying SCR-B from 4 to 2. Note that the parallel double AC lines of the HVDC link are kept connected in the system. Figure 11 shows the analytical results. When AC Grid-A is modelled as Thevenin circuit, varying SCR-B from 4 to 2 distabilize the system at the 18 Hz dynamic mode. In comparison to Section IV.C, the strength of AC Grid-A is set as SCR-A=2.8 (unstable) instead of SCR-A=6 (stable). After replacing Thevenin circuit with SM / RES, the variation of SCR-B from 4 to 2 does not destabilize the system. And this stability assessment result is furtherly validated by the EMT simulations as shown in Figure 12.
Note that the 7 Hz resonance mode observed in Section IV.A and IV.B has not been observed again as the strength of AC Grid-A has not been decreased down to SCR-A=2 again and the parallel AC lines of HVDC are connected in the system. For the scenarios with SCR-B=2, the PFs of grid nodes are presented in TABLE VI. After replacing the AC Grid-A (Thevenin circuit) with SM and RES, the left side AC grid has more participation in the 18 Hz dynamic mode as the summed PF of AC1-AC4 becomes larger, while the AC-DC dynamic coupling becomes weaker as indicated by the smaller PFs of the DC grid nodes. Another finding from the PF results is, the SM branch participates several times more than the RES branch in the 18 Hz dynamic mode. If furtherly increasing the RES penetration, the stabilizing effect from SM and RES in this case may become weakened. Further analysis on this effect will be presented in a separate paper.

E. DISCUSSION OF OSCILLATION MITIGATION DESIGN
In above subsections, participation factor analysis for the dynamic modes of a simple AC/DC grid was conducted considering different stability impacts. According to the results in TABLE IV -TABLE VI, the grid nodes AC3 and AC4 have the dominating PFs for the around 7 Hz dynamic mode and the grid nodes AC7 and AC8 have the dominating PFs for the around 18 Hz dynamic mode in most grid conditions as shown in Figure 1, Figure 5 and Figure 8, thus an oscillation mitigation scheme might be preferably configured at these grid nodes [9]. A very effective way is to add active damping control in the HVDC converters at AC3 and AC7, so that no additional passive or active damping devices should be installed there [7], [20]. Moreover, proper damping design in HVDC converters can de-risk the potential oscillation stability for a wide frequency range until the limits of their PWM modulation capability are reached [8], [32], [33]. Apart from the mitigation design in HVDC converters, the damping capability of the RES at AC3, i.e. the active damping in the VSM control scheme here, could also be fully deployed to furtherly increase the overall damping level of the system [29]. However, careful design should be conducted to prevent potential control interactions.

V. CONCLUSION
A straightforward nodal admittance matrix based method, that integrates the detailed dynamics of AC and DC grids as well as AC-DC dynamic coupling, is proposed to assess the stability of converter-interfaced hybrid AC/DC grids. The proposed method not only prevents knowing the number of RHP poles of system minor loop gain as required by Nyquist based method, but also helps identify dynamic modes and PFs of grid components. It is especially advantageous compared to state-space based method, because it does not strictly rely on white-box models, the partial admittance models can also be measured or obtained from product vendors.
The effectiveness of the proposed method is illustrated by case studies with a simple hybrid AC/DC grid consisting embedded HVDC, SM and RES. The stability impacts of varying grid strengths, converter DC voltage compensation mode and switching states of parallel AC line as well as the model variation of AC Grid-A are analyzed. Influences on the AC-DC dynamic coupling are analyzed through grid node PF analysis. It is shown, the GFL controlled HVDC faces potential resonance risk from weak grid connection, and the observed AC-DC dynamic coupling from HVDC converters can introduce the interaction of the connected AC grids. Besides, modelling the connected AC grid of an HVDC converter as simplified Thevenin circuit greatly changes the stability assessment result compared to detailed AC grid model. Moreover, the GFL and GFM converters and SM can operate stably in the considered conditions, and the SM branch participates more than the VSM-controlled RES branch in stabilizing the system in the case study.
In the future, highly meshed AC/DC grids are expected, and the AC-DC dynamic coupling effect is expected to get stronger when HVDC is configured with GFM control (due to the conflict between AC side power availability and DC voltage control) or in a multi-terminal HVDC system with drooped DC voltage control. Moreover, the ever-growing complexity of DC grid may introduce new dynamic modes and/or potentially have larger impact on the dynamic modes which are presently dominated by AC grid, i.e. some dynamic mode could be simultaneouly dominated by both AC and DC grid nodes. The proposed method could be used to analyze the stability risks from such configurations. Furthermore, as the study case is a small system, whether there is any limitation of the proposed method from system size or not would be addressed in future works.

A. SM MODEL
A round rotor SM is considered. The machine has one damping winding in the direct-axis (subscript D) and two damping windings in the quadrature-axis (subscript Q1 and Q2). The field winding is oriented in the direct-axis. When parameterized with fundamental parameters and written in per-unit form, the flux equations of the SM are formulated in (8) L  denoting stator leakage inductance and fd L denoting rotor field circuit inductance.
Swing and electrical torque equations in per-unit form are: where e m and T T are the electrical and mechanical torques,  is the rotor angle, and 0  is the per-unit steady-state frequency. On the basis that the prime-mover's speed governor is slow compared to the fast transients under consideration, m T is assumed to be constant. Moreover, the exciter and power system stabilizer (PSS) dynamics are disregarded, resulting a constant field voltage f u , if desired, these dynamics can be included in the model [34]. Then performing small-signal analysis to (8) and (9) (12) where Zdq is the SM output impedance under system dqframe. The corresponding admittance can be simply obtained by matrix inversion. The linearized impedance model is then validated using detailed time-domain model from MATLAB / Simscape, the validation results are omitted for saving sapce. The adopted standard parameters, as given in TABLE VII, can be transformed to the fundamental parameters as used in the derived model according to [35].

B. RES MODEL
The VSM model in [29] is applied for the RES generation, as shown in Figure 3. The VSM is composed of a VSM core, a virtual circuit and an active damping. The VSM core is used to provide a reference for the inner source voltage. The virtual circuit introduces a virtual impedance in series to the RES output filter, which increases the damping of the dynamic modes around nominal frequency. The active damping employs a lead-lag filter to counteract the highfrequency LC oscillation from the AC filter.
The  (17) where the superscript c indicates the variables under control dq-frame. The transformation between control dq-frame and system dq-frame can be illustrated by

C. HVDC Converter Model
Circuit and control diagrams of HVDC converters are shown in Figure 4.