Impedance-Based AC/DC Terminal Modeling and Analysis of the MMC-BTB System

Impedance-based small-signal stability analysis is widely applied in practical engineering with modular multilevel converters (MMCs). However, the deficiencies of existing impedance models (IMs) and the idealized extension for the single MMC influence the analyses in multiterminal systems. In this article, these gaps are filled by focusing on an MMC-based back-to-back system. To obtain the steady-state harmonics of the system, a numerical method is initially proposed based on Newton–Raphson iteration in the frequency domain. Then, by substituting the shared terminal dynamics with active or passive devices, theoretical AC/DC IMs based on typical control loops with pure time delays are directly established via multiharmonic linearization. Further aided by the derived IMs, two neglected aspects in the current literature are investigated, including the influence of power transformers on low-frequency impedance characteristics and the rationality of using simplified IMs for high-frequency resonance studies. The stability of interlinking systems needs to be comprehensively analyzed at both AC and DC terminals. The analyses help to position the instability source, obtain the stability margin, and guide the supplementary control strategy. All IMs and analyses are verified via frequency scans and simulations in PSCAD.


I. INTRODUCTION
MC-based multiterminal systems, featured by decoupled power control, good adaptability of various voltage levels, and high power quality [1], are expected to be widely applied in modern power systems.However, the recent real-world events, i.e., the 30 Hz subsynchronous oscillation [2] of the Nan-Ao wind power transmission project and the 1.2/0.7&1.8 Impedance-based AC/DC Terminal Modeling and Analysis of the MMC-BTB System Chongbin Zhao, Student Member, IEEE, and Qirong Jiang M kHz high-frequency oscillation [3], [4] of the Luxi/Yu-E grid interconnection project, attract broad concern for the system wideband small-signal stability, where the impedance-based method is welcomed for offering intuitive physical insights on the mechanism, protecting the privacy of vendors, and being able to adopt the field test.
Establishing the IM/AM that contains MMC is the basis of theoretical stability analysis, and extensive progress is reported in [2]- [20].IM describes the ratio of the small-signal response of voltage and current at each frequency of interest.Toward a typical linear time-periodic system, IMs can be indirectly and directly derived based on the 1) harmonic state-space in the time domain [2], [4]- [10], and 2) multiharmonic linearization in the frequency domain [11]- [16].The harmonic state-space-based method typically fails to approximate the pure time delay (cannot be omitted for MMC) above 100 μs using the rational polynomial at high frequency if starting from developing ordinary differential equations [5]- [7], [9], and even if great progress is made in [10] where the use of delay differential equations avoids the linearization of irrational transfer functions, the computation burden is high.Comparatively, the multiharmonic linearization sidesteps the flaws by eliminating intermediate variables and performing low-order matrix manipulation, and TIMs instead of AIMs that are common for TL-VSCs [17]- [19] are available (frequency responses are available for both TIM and AIM while explicit transfer functions are only available for AIM).Hence, multiharmonic linearization is more intuitive since IMs/AMs are frequency domain models and the theoretical basis of this paper.
In summary, efforts were mostly made on a single MMC in the earlier research [2]- [4], [6], [8]- [11], [16], [20], where ideal DC/AC voltage sources are deliberately assumed in general.The simplification diverges from the real world and leads to the low credibility of the obtained AC/DC IMs on stability identifying; thus urges a generalized IM framework that simultaneously considers both AC and DC terminal dynamics, which is not hitherto unified.It is also found that only the ideal transformer model is included in the current IMs and the mutual coupling is ignored, which is not considered for a long period.
In addition to the aforementioned imperfections of IMs for the single MMC, more issues are noteworthy when combining separate IMs for multiterminal stability analysis.Initially, the steady-state analytically computed in the time [21] or energy domain [22] with truncated harmonics considered, is doubtful with several simplifications, e.g., neglecting the active power loss or power flow constraints, while the measurement-based method [11] adopted in the frequency domain is invalid when systems are unstable, and the inaccuracy of steady-state impacts the IMs.Moreover, when the transformerless cases are focused in the distribution network for flexible interconnection [23], the ZS DM impedance of an MMC-BTB system is no longer infinite as same as the HVDC scenario, and the effect of virtual impedance on the stability induced by ZSCC is unclear.Finally, since the instability source can be any of the interconnections of AC and DC terminals [13], [18], [19], it is essential to analyze the stability with detailed IMs rather than simplified IMs (the complete control loops are modeled in the detailed IM while several control loops are omitted in the simplified IM) at each terminal to position the authentic unstable modes or verify some conclusions deduced from TL-VSCs, which can also help reveal the essence of some recent system-level stability analysis including MMCs [24], [25].In other words, executing only the DC stability analysis considering AC terminal dynamics [12] or vice versa [7] is not sufficient, After realizing the significance of facilitating the impedance-based stability analysis for complex networks with MMCs on the premise of satisfied system dynamic responses, basic works are completed by focusing on a three-phase symmetrical MMC-BTB system, which is a compromise between multiterminal and single MMC systems.This paper aims to contribute the following.
1) Since multiharmonic linearization is performed in the frequency domain, a practical frequency domain steady-state computing method is proposed based on Newton-Raphson iteration, which covers the cases of whether transformers are configured and is compared with the existing methods.
2) Based on the multiharmonic linearization, 1-d DC IMs are first obtained followed by the sequence domain 2-d AC IMs fully considering the frequency coupling, typical control loops, the pure time delay, and mutual coupling of transformers, whose accuracy is checked by frequency scans in PSCAD.
3) Aided by the derived accurate TIMs, sensitivity analyses are completed to explicate the influence on the stability of the nonideal transformer model and the virtual ZS DM impedance, and also to check the rationality of the deduced simplified IM with line-frequency NS current (suppressing) control.
4) Using the 1-d Nyquist criterion at the DC terminal and the 2-d Nyquist criterion at the AC terminal, the utility of IMs is verified by correctly positioning the instability source under multiple operating conditions.The adaptivity of a decoupled 1d Nyquist criterion [26] and the practical criterion proposed in [27] are also examined at the AC terminal.
The remainder of this paper is organized as follows.Section II presents the elementary knowledge of the investigated system.A novel steady-state computing method is introduced in Section III followed by the derivation and validation of TIMs in Section IV.Sensitivity and stability analyses are performed in Sections V and VI, respectively.Finally, Section VII concludes the work.

II. BASICS OF MMC-BTB SYSTEM
In this section, after introducing the configuration and control loops of the system in turn, the single-arm steady-state mathematical model is first derived in the time domain and subsequently transformed to the frequency domain, which lays the foundation of the subsequent computing method and IMs.Only the case of the ideal transformer model is focused on in this section for conciseness, and the cases of the nonideal transformer model or transformerless are comparatively discussed in Section III and Section V.

A. Configuration
Fig. 1 reveals the investigated MMC-BTB system based on the topology of real projects [3], [4].Current paths and definitions of variables are briefly marked in phase a of SE&RE.The adopted average model for each arm indicates that 1) the submodule capacitance-voltage that belongs to the single arm is balanced, and 2) the frequency range of concern is far below the equivalent switching frequency for stability analysis, and such assumptions are usually acceptable according to existing works.The lumped parameters are adopted for the grid side for simplicity since this paper focuses on the modeling of MMCs under various control modes.

B. Control Loops
Referring to the most common grid-following control, the loops are illustrated in Fig. 2 referring to [4], [15], where the SE regulates P and Q and the RE regulates vdc and vRMS for the outer loop control.NS current control and ZSCC are added in addition to the regular PS current control and CCSC for the inner loop control.Td represents the sum of control delay in the dq-frame while the sampling delay in the abc-frame is not considered, and the per-unit values are adopted for control.The key parameters are listed in Table I.

C. Single-arm Time Domain Model
For the IM of a single MMC, by grasping the symmetric and periodic distribution of the internal multiharmonic variables, [11] introduces a focus on the single arm instead of the single phase [2], [6], [9].Thus, this idea is extended to the BTB system.The dynamics of the upper arm, i.e., phase a of RE, can be reflected in the time domain: where the subscript "uar" is neglected for simplification unless specified.v dz is the voltage between the grounding and neutral point and acts as the source of ZS DM circuit.An extra equation of the ideal transformer model is: (2) Substituting ( 2) to (1) eliminates vw but retains vv and i, which matches the demands of the voltage sampling in Fig. 2 and the arm current decomposition of deriving IMs.

D. Frequency Domain Interpretation
By replacing the time domain variable multiplication with the vector convolution and further matrix-vector multiplication both in the frequency domain [11], (1) and (2) are converted to: YLR and YC are diagonal transfer function matrices: Hence, vv, m, i, and v are the steady-state vectors to be solved and form the corresponding Toeplitz matrices in (3).v dz induces the coupling of both ends, which differs from the single MMC: Due to the transformer, h in (4) guarantees the admittance of ZS DM circuit to be zero, and v dz can be removed from (3), so ZS decouples from PS or NS in the DM circuit.Such coincidence should be assessed in the transformerless case.

A. Idea and Hypotheses
Solving the phase and amplitude of each order of harmonics for a time domain variable can be considered as solving ag⟨i⟩ and bg⟨i⟩ of frequency domain vectors with s=jω1 (2πf1), jω2, …, jωn.Several hypotheses are formulated as follows: 1) Supposing that the steady errors are neglected, the control objectives of both inner and outer loops are equal to their references, e.g., vdc=v * dc =vdcb, thus vdc becomes a known quantity.2) Due to the conjugate relation between the (n+1±i) th elements of each vector, v introduces (2n+1) unknowns; the use of CCSC and ZSCC/transformer lead i⟨2⟩ and i⟨3⟩ to be 0, hence i introduces (2n-3) unknowns; based on Fig. 2, m⟨0⟩=0.5 while m⟨1⟩ (led by PS/NS current control), m⟨2⟩ (led by CCSC), and m⟨3⟩ (led by ZSCC) can be nonzero; vv⟨1⟩ is related with the AC power flow and introduces two more unknowns.Hence, at most (4n+6) equations are required for each end.
3) The iteration between v and i in (3) requires the truncation of each vector.Considering m⟨3⟩ is specially focused, n is selected as 4 in this paper.

B. Process and Annotations
A flow chart is offered to explain the process by first computing coefficients of SE and subsequently RE, as shown in Fig. 3, and annotations are added as follows.1) With the given P * , Q * , Zg, and vg, 4 coefficients contributed by vv⟨1⟩ and i⟨1⟩ can be independently solved based on the AC power flow for SE.Suppose that the transformer is first configured and m⟨3⟩=0; then, the remaining 4n coefficients of SE contributed by m, v, and i can be independently solved with the same number of real equations featured by the balance of ag⟨i⟩ and bg⟨i⟩ extracted from the matrix form of (3).
2) The DC balance holds for SE and RE: Eq. ( 6) satisfies because of the negligible electrical distance between the DC terminals in the grid interconnection scenario and can be adjusted if the voltage drop of DC line impedance must be considered in other scenarios.It is also confirmed for RE that the remaining 4n+3 coefficients should be simultaneously solved using the 4n equations from (3), two equations to represent the voltage drop of Zgr, and an extra equation that represents the RMS voltage control.
3) The aforementioned nonlinear algebraic equations can be numerally solved by Newton-Raphson iteration, and part of the initial values are given as (the rest are default to 0): 4) When handling each coefficient for the calculation, the phase adjustment relative to vv⟨1⟩ is necessary: abs( ) a b , ang( )=arctan(b / a ).
5) If the transformerless case is focused on, these results can serve as the initials added with a 0 m⟨3⟩ =0, b 0 m⟨3⟩ =0.With identical ZS circuit sampled and ZSCC for both ends, the follow-up ( 9) should be combined with the aforementioned 8n+3 equations extracted from (3) for a simultaneous solution:

C. Validation
The steady-state of the transformerless case derived by three methods towards the operating condition in Table I is listed and compared in Table II and Fig. 4, respectively, as an example.The real and imaginary parts of each coefficient are listed in the upper and lower roles of each space in Table II.In [11], m is first obtained by simulation and then brought into (3) to calculate i and v for method a, but obvious errors are observed especially i⟨2⟩ (should be 0 theoretically due to CCSC) and i⟨0⟩ (as (6) indicates).Particularly, since a steady m cannot be measured when the system is unstable, neither method a nor method c is rigorous, while the proposed method b has high precision and low computational overhead (normally, 4-6 iterations are adequate with proper initial values as given ( 7)), and avoids the complicated but unnecessary transformation from the time or energy domain to the frequency domain [21], [22].The magnitude of m⟨3⟩ is nearly 70% of that of m⟨2⟩ due to the nonexcessive L and C in the transformerless case that is possible for distribution networks.Considering the critical impact of CCSC on the stability [11], [12], it can be referred ZSCC cannot be neglected in the IM if it works and will be elaborated in Section V.

IV. DERIVATION OF AC/DC IMS OF THE MMC-BTB SYSTEM
The core work of developing IMs is introduced in this section.The modeling idea is clarified using the power stage and efforts are added to the control stage, whose transfer function matrices are given in the Appendix for readability.DC and AC TIMs are then obtained and verified.

A. Basic Modeling Idea
There is a duality between AC and DC IMs but was not explicitly announced.Multiharmonic linearize on (3) gives: Δm is closely related to the sampling and control as illustrated in Fig. 2: (11) It is worth noting that Δv does not contribute to Δm as Fig. 2 reveals, but if energy-based control generates the reference of CCSC, an extra Bv should be developed in (11).Since IM describes the relation between selected elements of Δi with Δvv or Δvdc, Δm and Δv are the intermediate variables and can be eliminated by combining (10) and (11).In addition, an extra equation to eliminate Δvv or Δvdc for DC or AC IM is available since the dynamics can be described from the two sides of each partition point in Fig. 1 based on Kirchhoff's voltage law.Due to the simpler realization of DC IMs, they are first derived and subsequently used for deriving AC IMs.
Furthermore, two points are claimed: 1) AMs should be first developed followed by the inversion to IMs and 2) since s=jωp of the IM corresponds to s=jω1 of the steady-state, there is a position offset between two elements in (Δg)⟨i⟩ and (g)⟨i⟩ that hold the same CM-DM/sequence distribution law [11].

B. Control Stage
Eq. ( 11) is transformed to (12) for the modeling modularity: ) where J and K describe the influence of the inner loop control on Δm from Δi and Δvv with the fixed reference values; Δd and Δq cooperate with E and F to reflect the influence of the outer loop control on Δm without the inner loop control.
Since Δvv=0 ensures Δθ=0, a process of current sampling, (inverse) Park transformation, and linear control ensures Q a diagonal matrix.P is related to the feedforward and phaselocked loop (PLL), where the former is similar to Q and the latter introduces off-diagonal terms.The closed-loop PLL transfer function considering multiple delays reflected by ho (T1/4 phase-shift offset) and hd (aggregated control delay) is: Δd and Δq depend on the control modes, whose relations with the sampling are described by matrices U, R, S, T, W, and V.The cascade controllers are contained in E and F. Finally, Bi, Bv dc , and Bv v can be expressed with matrices in this subsection.

C. Derivation and Validation of DC IM
Since idc is 3 times i cz ua while iw is 2 times i dp ua , the extra equation to eliminate Δvv for the 1-d DC AM/IM Ydc/Zdc is: Bringing ( 11) and ( 14) into (10) gives: The accuracy of DC IMs of both ends is validated as shown in Fig. 5.The solid lines represent the theoretical IM by taking frequency of interest into (15) while the marks represent the frequency scans.It manifests the feature of introducing passivity only for RE at the high-frequency as the zoom-in shows due to the sampling delay of vdc.

D. Derivation and Validation of AC IM
To derive the sequence domain 2-d AC AM/IM Yv/Zv, the following dual expression of ( 14) to eliminate Δvdc is: (16) where "^" indicates that only the subscript s or r of Zdc is inconsistent with that of the other matrices in (17) With the conjugate symmetry of 2-d Zac [5], only the two elements that belong to the same row are plotted for each end in Fig. 6, and the modeling accuracy is confirmed similar to Fig. 5. Since matrices A, B, Bi, Bv dc , and Bv v are shared in ( 15) and ( 17), there must be a generalized modeling framework for AC and DC IMs.

A. Issues of Nonideal Transformer Model
The classical modeling approach of the nonideal transformer is based on the following frequency domain formula in PSCAD: where Lvv, Lww, and Lvw are the self-inductances of winding v/w and the mutual inductance.Thus, (2) of the ideal transformer model turns to (19) of the nonideal transformer model: By calculating with ka, Ima, Xv, and Sb in Table I, kc is 0.9995 while both (Lvv/Lvw) and (Lvw/Lww) are close to ka.Generally, the coupling between iv and vv in the second equation of (19) induces an extra negative impedance Zvv shunt connected with the deduced Zac in (17), which mainly affects the low-frequency impedance characteristics of the AC IM and doubts the simple handling in [12].However, DC IM is almost unaffected since Zvv is far greater than Zg.Such inferences are supported as shown in Fig. 7, where the steady-state is recalculated considering (19), and the appropriate adjustments are given on control stages, e.g., replacing ka to (Lvv/Lvw) and updating Bv v for both AC and DC IMs while further revising (17) for DC IMs.

B. Issues of Transformerless Case
Here the doubt of the transformerless case on decoupling ZS with PS/NS DM stability analysis is removed, where ZSCC is activated and an extra submatrix is introduced for J: With the computation of m⟨3⟩ in ( 9), the principle of i dz ua =0 and v dz elimination is clarified in Fig. 8 (a).The sufficient virtual impedance balance Mv⟨3⟩ with its small-signal (MΔv+VΔm)⟨3⟩ of two ends over the wideband by hZSCC(s) and the idea of h in ( 4) can be moved from YLR to its product terms in A and B. Hence, two more effects on IMs induced by ZSCC are: 1) the change of M as Table II shows and 2) the release of YC(x,x), x=n+1±3 in Av, Bv, Adc, and Bdc.
A 1-d positive sequence Z p v is transformed from 2-d Zv and Zg for a more intuitive analysis [26]: By comparing the case transformer configured (ka=1, Xv=0) for isolation and the transformerless case in Fig. 8 (b), a larger kpZSCC damps the resonant peaks over 3f1 to 5f1 and has a cumulative effect at low-frequency according to (3).In addition, the influence on the stability of 3 rd harmonic injection to suppress the fluctuation of v [28] remains an open issue based on the conclusion of this subsection and is worth investigating.C. Issues of High-frequency Simplified IMs PS AIMs are welcomed to design notch-filter-based or selfadaptive high-frequency oscillation suppressing control [4], [16], [20].The proposed simplified AIM is given in (22) at the bottom of this page with the following annotations: 1) With the limited bandwidth of outer-loop control, only the inner-loop current control and voltage feedforward are included, which is a cognition for high-frequency simplified IM.The coefficients k should be elaborately managed based on the control loops in Fig. 2.
2) Distinguishing from the PS control that an fp/fn (fn=fp-2f1) perturbation in the sequence domain results in an fp-f1/fn+f1 perturbation in the dq domain through park transformation and then yielding an fp/fn response through inverse park transformation, the effect of NS control is that an fp/fn perturbation in the sequence domain results in an fp+f1/fn-f1 perturbation in the dq domain and finally an fp/fn response in the sequence domain, so the NS control should be fully modeled in (22) but is not considered for real project analysis in [4] or the related works in [15]- [20].
3) Another critical observation is that v⟨0⟩ replaces vdc⟨0⟩ since MMC uses v rather than vdc for modulation, and it proves the necessity of steady-state calculation again in Section III.
Rationality check of the simplified AC IMs over the highfrequency is reflected in Fig. 9. Whether NS current control is modeled or not in (22), nearly 20-Hz passivity region deviation or unacceptable errors due to the T1/4 phase offset are observed, as the blue or red solid lines show.It is also confirmed that the 4th-order Padé approximation on a 300 μs pure time delay at multiple sampling paths is most effective below 200 Hz, which casts doubts on the rational approximation when establishing TIMs considering time delay.Since the frequency responses are of practical requirements for stability analysis, keeping the preliminary irrational from e -sT is the key point of correctly treating the time delay, as shown in this paper and [10].

VI. CASE STUDY
Three cases are studied in this section: the change in short-circuit ratio (a), the parameter tuning of hv dc (s) (b), and the step of P (c), are discussed.Since some cases may not be close to those of the real world and the overhead lines are replaced by a lumped Zg in particular, the studies mostly verify the necessity and wideband accuracy of the TIMs.The simulation parameters are based on Table I without specification.
An alternative method to determine whether there is a righthalf plane pole in the transfer function matrix LR is to examine whether the subsystems in Fig. 10 can stably operate through simulation.The LR of each terminal in Fig. 1 is defined as: This method is acceptable in engineering since TIMs cannot directly offer system zeros and poles of MMC, which is emphasized to correctly apply Nyquist criteria with frequency responses.Considering the conjugate symmetry of eigp and eign for AC subsystems, and the self-symmetry of eigdc for DC subsystem, the symmetries are both about f1 for the sequence domain model, -2900 Hz to 3000 Hz of eigp and 50 Hz to 3000 Hz of eigdc with a step of 0.1 Hz are selected for the eigenvalue calculations.

A. Case A: Change of Short-Circuit Ratio
If the nonideal instead of the ideal transformer model is adopted and Lgr is increased from 22 mH to 25 mH (short-circuit ratio from 1.45 to 1.27), a diverged oscillation is first observed at t=4 s and subsequently becomes a sustained oscillation for all of ivra, ivsa (the harmonic amplitude is very small), and idc in Fig. 11.The FFT shows that |fp-f1|=46.5 Hz.Only eigpr of the nonideal transformer model judges an instability mode as shown in Fig. 12, whose model frequency (the initial oscillating frequency, a deviation with the sustained oscillation frequency is reasonable) and PM are approximately 98.2 Hz (|fp-f1|=48.2Hz) and -0.3°, respectively, aided by the equivalent 1-d PS Bode plot, as Fig. 13 shows.
The above theoretical analyses are consistent with the simulations with two additional remarks on case a: 1) using the ideal transformer model to deduce IM may omit the sub/super synchronous instability mode; 2) since the instability mode is only correctly identified at RE AC instead of SE AC or the DC terminal, the supplementary control should be added in the control loops of RE related to its pure AC dynamics for priority such as PLL or inner-loop current control instead of the DC voltage control, even if DC voltage control also influence the AC sub-super synchronous dynamic responses [7].

B. Case B: Tuning of Hv dc (s)
To ensure the universality of conclusions in subsection A, case b is studied, where the leakage inductance of transformer is deliberately omitted to examine the necessity of Section V, subsection B. A group of heuristic simulations is given for each subsystem, which has the time logic as follows: 1) At t=3 s, Lgr is set to 0 and kpv dc is set to 30; 2) At t=4 s, Lgr becomes its rated value L * gr (10 mH); 3) At t=6 s, kpv dc becomes 4; 4) At t=7 s, kpv dc becomes 0.01, Tiv dc becomes 0.005, and Lgr becomes 0 again; 5) At t=8 s, Lgr becomes L * g .AC/DC terminal and arm capacitor voltages as well as the phase/circulating currents are shown in Fig. 14, and further simulations prove that subsystem ② can stably operate whether the transformers are configured when hv dc (s)=30+1/0.05sinstead of hv dc (s)=0.01+1/0.005s,while subsystem ① cannot stably operate with any controllers and subsystem ③ can stably operate with both controllers (the simulation is neglected).The remarks on case b are: 1) the instabilities induced by the too large or too small kpv dc regard to different system modes, and detailed IMs should be used to correctly identify them; 2) the discussion in Section V, subsection B is necessary, and the effect of ZSCC on the stability over 3f1 to 5f1 should be faithfully modeled; 3) in subsection A, the instability mode is only identified at the RE AC terminal, whereas part of cases of kpv dc tuning can be analyzed at both AC and DC terminals of the same end (RE) but is no longer feasible for the AC SE terminal; thus, the independent black-box analyses are invalid for the (inter)harmonics of ivsa, which should be explained as a phenomenon of harmonic balance, thus urging an "impedance interaction" between multiterminal for stability modeling and instability positioning in future privacy-protected analyses.Certainly, the instability source of the studied hv dc (s) tuning cases is positioned at RE with the DC control loop to be improved in priority to avoid the unstable mode.The power control can also be selected to tune for hv dc (s)=30+1/0.05swhen transformers are configured since it acts on the coupling between AC and DC terminals.

C. Case C: Step of P
Finally, as a practical concern, the stability of P steps from 0.5 pu to 1 pu, is focused on, where L * gs is changed to 10 mH to emulate a weak grid.A group of heuristic simulations is given for each subsystem, which has the time logic as follows: 1) At t=3 s, Lg remains at L * g , kpP is set to 2, and P * is set to 0.5 pu; 2) At t=4 s, kpP becomes 2.5; 3) At t=5 s, P * steps to 1 pu; 4) At t=7 s, Lg becomes 0; 5) At t=8 s, kpP becomes 2; 6) At t=9 s, Lg turns to L * g .Fig. 17 claims that kpP=2 or 2.5 will make the system unstable when P * =1 pu on the premise that only subsystem ② cannot stably operate, where the frequency of sustained oscillation increases with a larger kpP (|fp-f1|=230 Hz and 236 Hz when kpP=2 and 2.5).Nevertheless, the frequency domain analyses using Fig. 18 and 19 illustrate that the same unstable mode can be identified at both DC and AC SE terminals when kpP=2 but can only be identified at the AC SE terminal when kpP=2.5, while the instabilities cannot be determined at the AC RE terminal (the simulation is neglected).The model frequency increases from 296 Hz to 299 Hz while PM decreases from -4.9° to -21.7°, which satisfies Fig.    Case c verifies the previous finding that supposing a specified terminal as the instability source for analysis may leave out the authentic unstable mode.Regarding the adaptivity of the simplified AC stability criteria for MMC from TL-VSC, more interesting findings are noted.First, the extra unstable modes appear in the 1-d PS loop criterion [26] as the red solid lines in the zoom-in of Fig. 18

VII. CONCLUSION
By focusing on an MMC-BTB system, this paper discusses several issues in extending the impedance-based analysis from a single MMC to multiterminal systems.With theoretical analyses and validation, the following conclusions are drawn: 1) A numerical method is proposed to rapidly and accurately obtain steady-state of the system regardless of whether transformers are configured, which helps unify the complete small-signal stability analysis in the frequency domain.
2) By clarifying the modeling idea, a generalized framework is offered to uniformly and succinctly develop AC/DC TIMs, which considers the typical control loops and the pure time delay in particular, and is verified by frequency scans.
3) Sensitivity analyses are adopted to emphasize the essential effects of the non-ideal transformer and ZSCC on the impedance characteristics, while the rationality of a novel simplified AIM considering the NS current control is comparatively studied.4) Multiple case studies collectively confirm the necessity of executing AC/DC stability analyses with detailed IMs.The unstable mode can be identified at the AC or DC terminal of one MMC and the harmonics spread to the entire system, and it may not be suitable to replace Nyquist criteria with the simplified criteria in wideband stability analyses from TL-VSC to MMC.
Deeply understanding the modeling idea of this paper opens access to more complicated device-level and system-level impedance-based modeling and analysis.Future studies will focus on the couplings between multiple physical paths and the stability identification with multiple unstable subsystems in the multiterminal system.

APPENDIX DETAILS OF CONTROL STAGE
The basic J can be decomposed according to Fig. 3, The elements not mentioned default 0 and the same below:

Fig. 3 .
Fig. 3. Flow chart of the frequency domain steady-state computing method.

Fig. 4 .
Fig. 4. Amplitude error of methods a and b compared with method c.

Fig. 6 .
Fig. 6.Validation of the AC IMs.(a) RE.(b) SE.V. SENSITIVITY ANALYSIS Since the derived IMs are based on (3) of the ideal transformer model in Section IV, the influence on IM of the nonideal transformer model and the case without transformers remain unclear.Detailed IMs can also serve as standards to check the errors of simplified IMs at high frequencies.Due to the pure time delay and the inverse of Av or Adc in AMs, AIMs are impractical, and the sensitivity analyses are completed by adjusting key parameters and replotting TIMs point-by-point.

Fig. 8 .
Fig. 8. Influence on AC/DC IM of whether any transformer configured.(a) Diagrams of the principle of i dz ua =0 for various cases.(b) 1-d PS Bode plots.

Fig. 10 .
Fig. 10.Definition of subsystems for stability analyses at various terminals.
Furthermore, two instability cases are identified by eigdc when L * gr is accessed to RE with transformers configured in Fig 15 and 16, whose corresponding Bode plots match the FFTs in Fig. 14 (|fp-f1|=206 Hz & 12 Hz).eigdc also confirms the stable operation of the transformerless case with a PM of 1.15° when hv dc (s)=30+1/0.05s.However, AC IMs cannot determine the instabilities that are induced by the improper hv dc except eigpr in

Fig. 15 (
Fig. 15 (two eigps are omitted and they cannot determine the instabilities).
and 19.Since |fp-f1|=243.6Hz or 238.2 Hz deviates from the real 246 Hz or 249 Hz in two cases, it may affect the estimation of sustained oscillation frequency when kpP increases.In fact, 1-d and 2-d Nyquist criteria are only equivalent in marginal stability judging[26].Furthermore, a practical determinant(Det(s))-based criterion[27] is examined since the zeros and poles are difficult to obtain for IMs, as previously mentioned.By extracting a pair of conjugate right-half plane poles from Det(s), a stable mode satisfies: ), Im(Det(s)) and Re(Det(s)) separate the imaginary and real parts of Det(s), and szc is the zero-crossing angular frequency of Im(Det(s)).In Fig.20when kpP=2.5, the 1540.6Hz mode is misjudged to be instability except for the expected 299 Hz, which is induced by incorrectly handling the time delay, as shown by the purple triangles shown in Fig.12 (e).Overall, the 2-d Nyquist criterion is essential AC stability analyses of MMC-BTB systems.

Fig. 20 .
Fig. 20.Validation of the practical criterion [27] in the MMC-BTB systems considering the nonignorable time delay.

E
Δm contributed by PLL are divided into: x n 2, y n 1, w 0.5 (x n) sign[(x n) (y n)], Δd and Δq are different for each end as (x=n, y=n±1): and F are given for each end as (x=n): dc , and Bv v are expressed for each end as: (single) apostrophe means to integrate the matrices with hd(s) ([ho .hd](s)) of Δd or Δq into E or F.

TABLE II COMPARISON
BETWEEN CALCULATING METHODS OF THE SYSTEM STEADY-STATE (UPPER ARM, PHASE A, TRANSFORMERLESS CASE)Variable RE SE method a method b method c method a method b method c