A Multi-Domain Co-Simulation Method for Comprehensive Shifted-Frequency Phasor DC-Grid Models and EMT AC-Grid Models

To accurately capture the dynamics of a large-scale ac/dc system as a whole and the interactions between its individual components, a simulation method with high precision and efficiency is in great demand. For this purpose, we develop a multi-domain co-simulation method, in which the target system is partitioned into multiple dc and ac subsystems, represented by our proposed shifted-frequency phasor (SFP) models and the traditional EMT models, respectively. SFP models can be simulated with a much larger time step, leading to a significant improvement in simulation efficiency under a given requirement of precision. Further, a new interface model, namely, hybrid multi-domain transmission-line model (HMD-TLM), is developed to reflect the interactions between SFP models and EMT models, with the additional benefit of producing instantaneous and phasor waveforms simultaneously. Thus, the multi-domain co-simulation is implemented based on the efficient SFP models, the HMD-TLM, and a designed time sequences of simulation. The performance (efficiency and accuracy) of the proposed method has been validated on the well-known CIGRE ac/dc system as well as a practical system integrating large ac grids and modular-multilevel-converter-based multi-terminal dc grids.


I. INTRODUCTION
R ECENTLY, the modular multi-level converter (MMC)   has become a dominant topology for high-voltage direct current (HVdc) applications, such as the INELFE France-Spain HVdc project [1] and the Nanao MTDC Project in China [2].Compared to the conventional two or three-level voltagesource converter (VSC), an MMC allows higher power-handling capability with the additional benefits of reduced switching power losses and low harmonic distortion.As more power electronic devices, including MMC and the high-frequency-link (HFL)-based dc solid-state transformer (DCSST), are integrated into the power networks, system dynamics are re-shaped not only by the characteristics of the conventional ac components (e.g., synchronous generators, transformers, and transmission lines/cables) but also by the nonlinearities and fast-switching actions of power-electronic-based dc components.As a result, the frequency band of system dynamics becomes exceedingly higher and wider than before [3], [4].To precisely capture the dynamics of a large-scale ac/dc system as well as the detailed interactions between its individual components, an accurate and efficient electromagnetic transient (EMT) simulation method is highly necessary.
However, traditional EMT programs are usually inefficient to simulate a large ac/dc system containing MMC-based multiterminal dc (MTDC) grids because the computational burden is dramatically increased due to the following factors: 1) a significant increase in the number of nodes because of the expanding scale of the system; 2) a much smaller simulation time step (usually 20-50 μs or even less) to capture the switching actions of numerous semiconductors; 3) the time-consuming update of the system admittance matrix at each switching occasion; 4) much more additional calculations aroused by nonlinear components to solve the network equations.To improve the simulation efficiency from the aspect of MMC modeling, several accelerated models have been proposed in the previous work [8]- [11], for instance, the switch-function model [8], the Thevenin equivalent model [10], and the arm averaged model derived from state-space representations [11].However, the previous models have one or more of the following drawbacks.
Authorized licensed use limited to: UNIVERSITY OF ALBERTA.Downloaded on April 18,2022 at 18:54:22 UTC from IEEE Xplore.Restrictions apply.

R E A D
O N L Y TABLE I CHARACTERISTICS OF DIFFERENT SIMULATION METHODS 1) They usually adopt a unanimous simulation time step, which, according to the numerical stability criterion [9], has to be limited to a very small value, for example, 10 μs or less, to simulate nonlinear components accurately.As a result, the advantages of the accelerated model of MMC cannot be fully utilized.
2) The simulation efficiency could be improved by these accelerated models; however, the accuracy is not always satisfactory, especially when the time step is larger than 20 μs [12].
3) The previous accelerated models mainly focus on the dynamics at the system level, failing to give the dynamics at the sub-module level simultaneously.
Besides, MMC-based MTDC grids generally contain emerging and sophisticated components, such as the HFL-based DC-SST.The accuracy and efficiency in modeling the individual component imposes a direct impact on the overall simulation accuracy and efficiency.Therefore, developing more efficient and accurate models for the components of an MTDC grid is still an open issue.Although the parallel processing technique based on transmission-line partitioning is applicable for EMT simulations such as PSCAD/EMTDC, the time-step generally has to be limited to 10-20 μs due to the existence of MMCs.Consequently, it is still inefficient to simulate very large-scale MMC-based ac/dc systems [18]- [20].Many efforts have also been made to counteract the skyrocketing scale of ac/dc grids, for instance, the work listed in Table I.In [13] and [14], a multiport Thevenin or Norton equivalent is used to represent the external ac system to improve the simulation efficiency.However, its drawback is obvious: it can hardly reflect the nonlinear and frequency-dependent EMTs as well as the interactions between MMCs and ac grids.The combination of transient stability (TS) and EMT simulations (so-called TS-EMT co-simulation) is another popular way for simulating large-scale ac/dc systems [15], [16].But in this method only low-frequency dynamics of ac grids can be captured and projected upon the dc subsystems.The interactions at high frequency cannot be accurately simulated.Multi-rate co-simulation was proposed in [9] and [17] to improve the accuracy of higher frequency responses.However, the numerical stability depends heavily on the strategy of network partitioning.Mostly the simulation of MTDC grids has to be carried out with a very small time step (e.g., 10-50 μs) to meet the accuracy requirement.As a result, the computation burden is so heavy that the simulation efficiency remains a big problem.
To address the aforementioned issues, a multi-domain cosimulation method is proposed in this paper.The contributions are threefold.
1) New shifted-frequency phasor (SFP) models based on nodal analysis and rotational matrix transformation are derived for typical components in dc grids, including MMC, HFL-based DCSST and long transmission lines.They use a much larger time step than the traditional EMT models, and thus, improve the simulation efficiency dramatically.
2) A new interface model between SFP and EMT subsystems, i.e., hybrid multi-domain transmission-line model (HMD-TLM), is proposed.It produces instantaneous and phasor waveforms simultaneously, and thus, either low-frequency or high-frequency dynamics or interactions can be captured and reflected accurately.
3) A new multi-domain co-simulation framework is developed based on SFP models, HMD-TLM, and a especially designed simulation time sequence.Case studies on a benchmark model and a practical system have fully demonstrated its capability of simulating very-large ac/dc grids with high efficiency.
The rest of this paper is organized as follows.Section II introduces the SFP modeling.Section III elaborates the HMD-TLM and multi-domain co-simulation.Section IV examines the performance of the proposed method on the modified CIGRE ac/dc system model and a practical system integrating ac and MTDC grids.Brief conclusions are finally drawn in Section V.

II. TRANSFORMATION-BASED SFP MODELING
Several modeling techniques were proposed to capture wide frequency-band interactions between different components in large-scale ac/dc grids.
1) State-space-based dynamic phasor (DP) model [21], which decomposes state variables into harmonics of different frequencies based on state-space equations.However, state-space-based modeling is usually more time consuming than the nodal analysis approach because the model order increases as each DP is included as a state to the model.Moreover, the system-level equations will be expanded by k times if the kth DPs are included for the MMC.2) Nodal-analysis-based DP model [22], [23], with which each variable is still represented by the corresponding DP.As a result, the interfaces between DP and EMT models should be carefully designed concerning the conversion between DP and EMT quantities.3) Frequency-adaptive model [24], [25], which, as a nodal analysis approach, is based on the Hilbert transformation.
For instance, the models in [24] and [25] focus on the frequency-adaptive model of traditional ac components, however, without giving frequency-adaptive models for dc components.In order to capture instantaneous values and wide frequency-band phasors in an accurate and efficient way, the nodal-analysis-based SFP model with a specific rotational matrix transformation is proposed and detailed hereafter.

A. Concept of Transformation-Based SFP Modeling
In a typical power system, the variables usually have bandpass characteristics with the frequency band centering around the fundamental frequency ω s , which can be represented by the SFP-derived phasor form [24], [25], [30] where S(t) is the complex envelope of the time-domain signal S(t) and reserves the low-frequency dynamics of S(t); and s I (t) and s Q (t) are the real and imaginary parts of S(t), respectively.Alternatively, the SFP-derived phasor can also be obtained by the Hilbert transformation, or with denoting the Hilbert transformation.

B. Derivation of SFP Modeling
Suppose the dynamic equation of each component is written as After some manipulations, ( 4) is converted into its SFP form, of which the differential equation is given by where û = [û x , ûy ] T , and ûx and ûy are the real and imaginary parts of û.
With the trapezoidal algorithm, (5) is discretized as where û(t) and F (t) denotes state variables and the differential term; and Δt is the time step.
Variables of the SFP form in (6) are then transformed back into their time-domain signals And ( 7) is rewritten in the real and imaginary parts separately as Theoretically, the nodal-analysis-based SFP model using the specific rotational matrix transformation has the following advantages.
1) They can use a much larger time step than the traditional EMT models so that simulation efficiency will be dramatically improved without the loss of accuracy.
2) The dynamics at both system level and sub-module level are preserved.
3) They produce instantaneous and wide frequency-band phasor waveforms simultaneously.
4) The implementation of SFP models is derived directly according to (9), in which the matrix transformation can be realized as an independent module.As a result, the coding of SFP models for practical utilization is very simple and efficient.
Thanks to the aforementioned distinct features, the SFP concept is adopted hereafter to develop the models of different components.

A. SFP Models of Single-Phase/Three-Phase Transformers
In developing the SFP model of a single-phase transformer, the couplings between primary and secondary sides of the transformer are represented as a controlled current source and voltage source [see Fig. 1(c)], respectively [20].The differential equations can be expressed in the following shifted phasor form: where ûxy m 1 and ûxy m 2 denote shifted phasors of primary and secondary voltages; and îxy m 1 and îxy m 2 denote shifted phasors of primary and secondary currents.According to the derivations of ( 4)-( 9), ( 10) is discretized by the Trapezoidal algorithm with where G T is the conductance matrix; J T (t − Δt) is the equivalent current; K u and K i are parameters for J T (t − Δt); and R(t) is the rotational matrix.
As shown in Fig. 2, the SFP model of a three-phase transformer can be easily developed by combining the Norton equivalents of all phase windings.Taking phase A as an example, the primary and secondary currents have the following relation: Since no external current is injected into node 5 or i xy 5 (t) = 0, there is ) Considering that i xy With all the equations of primary and secondary windings combined, the SFP model of the three-phase transformer is finally obtained as The saturation of the transformer is represented by additional nonlinear compensation current [26].Moreover, if there is external current, the equivalent current term can be added on the right side of (19).

B. SFP Model of Bergeron-Type Transmission Line
The Bergeron-type transmission line is adopted to reflect the effect of distributed parameters [27], [28].Here, the singlephase and three-phase transmission lines are adopted for dc and ac grids, respectively (see Fig. 3).Their SFP models will be derived as follows.

1) SFP Model of Single-Phase Transmission Line:
The equations of a single-phase transmission line, represented in the form of shifted phasors, are given by where û(x, t) and î(x, t) are shifted phasors of voltages and currents, respectively.Given the initial conditions of û(x, t), î(x, t), or û(x, 0) = 0, î(x, 0) = 0, the Laplace transformation of ( 20) is Applying the partial derivative for x on both sides of ( 21), there is Performing inverse Laplace transformation to (22) yields where the speed of wave is v = 1/ √ LC .By converting shifted phasors into their time-domain forms, the SFP model of a single-phase transmission line is obtained as with where τ = l/v, θ = ωl/v, and Z = L/C.The SFP model given in (24) can be further separated into their real and imaginary parts as follows: where "⊗" denotes the Kronecker product, or the element-byelement multiplication.Finally, the equivalent current vectors are calculated as 2) SFP Model of Three-Phase Transmission Line: For the SFP model of three-phase transmission line, (20)  where Û(x, t) and Î(x, t)are shifted phasors of three-phase voltages and currents; and L, C are three-phase inductance and capacitance matrix per length.Eigenvalue analysis is applied to produce diagonal matrices.Thereby, coupled equations in the phase domain can be transformed to decoupled equations in the modal domain.Each equation in the modal domain is solved independently by using modal travelling time τ and modal surge impedance Z.
The transformation matrices between phase and modal quantities are different for voltage and current, i.e., where T u and T i are the modal-phase transformation matrices for three-phase voltages and currents; and Ûm and Îm are shifted phasors of voltages and currents in the modal domain.Substituting ( 29) into (28) gives (30) As LC = CL, L, C ∈ M 3 (M 3 denotes the matrix of the dimension 3), L and C can be simultaneously diagonalized.In other words, there exists a matrix T satisfying T = T u = T i so that T −1 u LCT i and T −1 i CLT u are both diagonalizable.Consequently, the SFP model of the three-phase transmission line in the modal domain is formulated into where m 2 ; L 0 and L 1 are the three-phase inductance of zero and positive sequences, respectively; and C 0 and C 1 are the three-phase capacitance of zero and positive sequences, respectively.The SFP model of the frequency-dependent transmission line can be applied according to [28].

C. SFP Model of HFL DCSST
To promote the efficient utilization of renewable energy, HFLbased DCSSTs are used to achieve voltage conversion and power transfer between HVdc and medium-voltage dc (MVdc) buses [29].As shown in Fig. 4, an HFL-based DCSST is composed of multiple identical dual-active-bridges (DABs) connected in parallel at the MVdc bus and in series at the HVdc bus.In a single DAB module, there are two full-bridge sub-modules on the high-frequency transformer.
To enable the flexible operations in the MVdc system, three control modes are used, namely, HVdc control mode, MVdc control mode, and power control mode [29].In the HVdc control mode, the voltage of the MVdc bus is fixed, and the voltage of the HVdc bus is controlled by the HFL-based DCSST.The amplitude and direction of the DCSST current is passively determined by the loads at the MVdc bus.In the MVdc control mode, the voltage of the HVdc bus is fixed, while the voltage of the MVdc bus is controlled by the DCSST.In the power The SFP model of an HFL-based DCSST is composed of SFP models of the DABs.Each DAB module contains a highfrequency transformer and two full-bridge sub-modules.The SFP model of the high-frequency transformer can be referred to as Part A. The SFP model of the full-bridge sub-module is modeled as a switch-dependent Thevenin equivalent circuit, comprising SFP models of capacitors and switches [30].For each full-bridge sub-module, it has two states in the normal condition, i.e., positive inserted and negative inserted [see Fig. 4(b)].Thus, the output voltage is the capacitor voltage or the opposite.
According to Eqs. ( 4)-( 9) [30], the SFP model of the capacitor, denoted as its Thevenin equivalents, is given by where R ceq is the Thevenin equivalent impedance, and v xy ceq (t − Δt) is the Thevenin equivalent voltage.
Parameters in (33) are calculated as If researchers would like to add frequency-dependent dynamics of high-frequency transformers, the SFP model of frequencydependent high-frequency transformers can be applied according to [24] and [28].

D. SFP Model of the MMC
The configuration of an MMC, which is composed of seriesconnected half-bridge sub-modules, is shown in Fig. 5.For each sub-module, the SFP model is the combination of the SFP model of its capacitor and of the switches.Represented in the Thevenin equivalent form, the SFP model of each sub-module is expressed as  The equivalent resistance and voltage of the sub-module in (35) are calculated by where R 1 and R 2 are the R on /R off resistance of the switches S 1 and S 2 ; v x i , v y i , i = 1, 2 are the corresponding ON/OFF-state threshold voltages in the xy coordinate; and R ceq and v xy ceq (t − Δt) are the equivalent resistance and voltage of the capacitor.
For each sub-module, the capacitor voltage is updated only when the sub-module is inserted.Specifically, the capacitor current for each sub-module is calculated by By connecting the Thevenin equivalents of upper and lower sub-modules in series, the overall switch-dependent Thevenin equivalent model of the arm, which are determined by the switch where "H" stands for the half-bridge topology.In Fig. 5, both the normal state [see Fig. 5(a)] and the blocking state [see Fig. 5(b)] are illustrated.Specifically, when the submodule is in its normal state, the capacitor voltage of each submodule can be charged or dis-charged according to the direction of the arm current.However, when the sub-module is blocked, the capacitor of each sub-module is only allowed to be charged due to the clamping diode.

IV. SYSTEM-LEVEL FORMULATION BASED ON SFP MODELS
With all the network components represented by SFP-derived Norton equivalents, the whole system is then formulated into nodal voltage equations, the solution to which prompts the efficient simulation of MTDC grids [27], [28].The system-level equation is written as ) where G is the complex-value-based conductance matrix; [v x (t), v y (t)] T is the vector of nodal voltages; [i x (t), i y (t)] T is the vector of external current sources; and [J x (t − Δt), J y (t − Δt)] T is the equivalent current vector.

V. MULTI-DOMAIN CO-SIMULATION OF AC/DC GRIDS
The SFP models derived in the previous section can adopt much larger time step than traditional EMT models without compromising the computation accuracy, and thus, improve the simulation efficiency dramatically.To take full advantages of SFP models in accelerating the efficiency of simulating large ac/dc systems, the concept of multi-domain co-simulation is proposed.The basic idea is that: the target system is partitioned into dc and ac subsystems, represented by our developed phasordomain SFP models and traditional time-domain EMT models, respectively.A novel interface model, namely HMD-TLM, is designed to exchange the boundary information between SFP and EMT subsystems, so as to fulfill the simulation of the whole system in a simultaneous and interactive way.

A. Framework of Multi-Domain Co-Simulation
As shown in Fig. 6, the target system is composed of the traditional generator-based ac grids, power electronic converterbased dc grids and transmission lines.The interactions among them are in different time scales.In the framework of multidomain co-simulation, the system is organized into SFP and EMT subsystems, with HMD-TLMs as the interface models among them.The MTDC grids are usually implemented into one or more SFP subsystems so as to accelerate the simulation.However, the ac grids have to been placed in the EMT subsystem, due to the difficulties in developing SFP models for the traditional generators and transmission lines that can represent their nonlinear and frequency-dependent characteristics.As the interface models, HMD-TLMs can reflect the couplings, especially those of high frequency, among EMT and SFP subsystems so that the interactions among different components, for instance, generators, converters, and networks [25]- [27], can be captured precisely.

B. Network Partitioning Based on HMD-TLM
There are two widely used network-partitioning methods, namely the transmission-line-based partitioning method [9] and the latency insertion-based partitioning method [34].The former is generally adopted for the partitioning of ac grids; while the latter is used in electronic circuits.However, both methods require that the subsystems to be partitioned should be modeled in the same way.Thus, neither of them is feasible for the multi-domain co-simulation, where SFP and EMT models are adopted.To address this problem, the HMD-TLM is proposed here to bridge the subsystems with different types of models.
As shown in Fig. 7, the HMD-TLM is a dual-port model, with each port represented by a Norton equivalent.The ports in EMT and SFP subsystem are of time-domain and phasor-domain representations, respectively.There is no direct connection between the terminals of the two ports (node k and node n).The influence from one subsystem on the other is reflected by a Norton equivalent that is composed of a paralleled impedance and a controlled current source with a time delay (traveling time) τ .Since the time delay is usually not an integer multiple of the time step, historical data of current sources on either side are extracted and interpolated to get the correct travelling time.Moreover, to utilize HMD-TLM for co-simulations, all its equivalent parameters and variables should be accurately and timely updated to reflect the interactions between EMT and SFP subsystems.

1) Update of HMD-TLM in the EMT Subsystem:
The HMD-TLM in the EMT subsystem is updated through the iterative calculation of the equivalent current source I k (t − τ ) using the following formula where u n (t − τ ), i n (t − τ ) are instantaneous interface voltages and currents of node n in EMT subsystems with the time delay τ , which are converted from the shifted phasors of interface voltages and currents in SFP subsystems.For interface voltages, the conversion is fulfilled by T are the interface voltages in the SFP subsystem, which are calculated using interpolations in phasor domain where the time delay τ is supposed to be within the interval For interface currents, the same conversion is used, or 2) Update of HMD-TLM in the SFP Subsystem: The update of HMD-TLM in the SFP subsystem is fulfilled in a similar way: first, the paralleled impedance matrix Z xy in the phasor domain is calculated by Z xy = Z ⊗ I 2×2 .Then, the equivalent current vector I xy n (t − τ ) in the SFP subsystem with the time delay of τ is calculated by where τ = l/v, θ = ωl/v, and Z = L/C.In ( 46), the phasor-domain interface voltage in the SFP subsystem, u xy k (t − τ ), is obtained by applying Hilbert transform on the time-domain interface voltages of the EMT subsystem, where R(t) is the rotational matrix; and H [.] is the Hilbert transformation.Specifically, when the given signal is an sinusoidal one, the Hilbert transformation can be simply obtained by adding a delay of π/2 to the signal [30].u k (t − τ ) is obtained by linear Authorized licensed use limited to: UNIVERSITY OF ALBERTA.Downloaded on April 18,2022 at 18:54:22 UTC from IEEE Xplore.Restrictions apply.

R E
Similarly, the interface currents are calculated as with Note that the HMD-TLM imposes no limits on the position of network partitioning and waveform distortion of interface variables.Moreover, its impedance parameters are calculated much more easily than the frequency-dependent network equivalent (FDNE) [35] does.Thus, HMD-TLM offers higher flexibility than the traditional methods and is suitable for different simulation scenarios.

C. Procedures of Multi-Domain Co-Simulation
Fig. 8 illustrates the overall procedure of multi-domain cosimulation.The first step is network partitioning.Then, EMT and SFP subsystems are initialized sequentially.Next, at each time step, they are simulated independently.Note that in our method, the same time step is adopted for both EMT and SFP subsystems and it can be extended to 100 μs, which, however, is infeasible for traditional methods.The interactions between different subsystems, especially those of high frequency, are reflected by the HMD-TLM.Its paralleled impedances and Norton equivalents are updated with the method presented in the previous subsection.The obtained interfacing variables are transmitted through the communication networks based on the shared memory technique [9] due to its capability of maximizing the simulation efficiency.Finally, as the total simulation time (T max ) is reached, the co-simulation comes to an end and all results are outputted to the users.

VI. CASE STUDIES
In this section, the proposed method as well as the traditional method is applied to two test systems, i.e., the modified CIGRE ac/dc system and a practical ac/dc system, to evaluate their performances in terms of efficiency and accuracy.Case studies on the first system focus on dynamic interactions caused by transmission lines, and thus, synchronous generators are simplified into ideal voltage sources.But in the second system, each generator is represented by its detailed EMT model containing governor and exciter so as to investigate dynamic interactions among generators.Different time steps (from 20 to 100 μs) are adopted to compare the proposed method with the traditional one.Meanwhile, the simulation results obtained with a unanimous time step of 20 μs are used as the high-fidelity reference values.

A. Modified CIGR AC/DC System
The modified CIGR ac/dc system is shown in Fig. 9. Its ac grid is the IEEE 39-node system and the dc grid is adapted from CIGR's benchmark model.In the ac grid, the generator is modeled as an ideal voltage source and the Bergeron model is used for each transmission line.In doing so, the electro-mechanical and electro-magnetic dynamics of generators are neglected.But the EMTs of the network, mainly determined by transmission lines, are reserved.The dc grid is divided into two MTDC subsystems with one in mesh and the other in radial connections, respectively.In the MTDC grid #1, converter no. 1 is responsible for maintaining the constant dc voltage; and converters no.2-4 control the power flow.Similarly, in the MTDC grid #2, converter no. 5 adopts constant dc voltage control; and converters no.6-8 adopts constant power control.The ac grid and the MTDC grids belong to EMT and SFP subsystems, respectively.
One of the challenges confronted by the co-simulation is to meet the accuracy requirement even when the MTDC grids contain various power-electronics-based dc converters.Different scenarios are simulated later to show how this challenge is overcome by our proposed method.In all scenarios, the time steps of SFP and EMT subsystems are set as 20 μs.

1) Scenario 1-1. Step Change in the d-axis Current Reference of a Converter:
In this scenario, the d-axis current reference (I dref ) of converter no. 4 is changed step wise, which stimulates dynamics in the dc voltages of converters as well as sub-module capacitors.Figs. 10 and 11 show the dynamic curves obtained with our proposed method side by side with the reference curves.The differences are also displayed.The results match each other very well, even during the fast EMTs.The maximum relative errors of dc voltages, capacitor voltages, and consumed reactive power are less than 0.003%, 0.005%, and 0.001%, respectively.Thus, the very high accuracy of the proposed method is verified in this scenario.
2) Scenario 1-2.A Symmetric Resistive Fault in the AC System: A symmetric resistive fault is triggered at Bus 29 from t = 3.0 s and is cleared 0.05 s later.DC quantities, capacitor voltages, as well as the simulation errors are shown in Fig. 12.All of plots have demonstrated that the results obtained with the proposed method matches the reference values very well.The maximum errors of dc voltage, dc current, and capacitor voltage are no more than 0.04 kV, 1 A, and 0.01 kV, about 0.02%, 0.1%, and 0.005% of the rated values, respectively.Fig. 13 shows the ac voltages of Bus 29, including instantaneous values in the xy coordinate and their corresponding SFP curves.Clearly, the latter matches the envelopes of the former exactly.Of particular note is that the SFP can perfectly reflect interactions of both high-and low-frequency EMTs, even around the occasions of the occurrence and clearance of the fault.Since the generators in the ac system are simplified into ideal voltage sources, the said high-frequency EMTs are mainly aroused by the transmission lines.

B. Practical AC/DC System
The target system is adapted from a practical ac/dc system that has a four-terminal MMC-based ±200-kV MTDC grid.As shown in Fig. 14, the MTDC grid is connected to the neighboring ac grid at four nodes, namely, SJ, SD, GN, and MM.To operate the MTDC grid, converter no. 1 maintains the dc voltage and the other converters regulate power flows.An HFL-based DCSST is connected to the dc bus at converter no. 1 to transform the high dc voltage down to 15-40 kV so as to supply an MVdc system.As compared with the first test system, this practical system is much more complicated: the number of generators and transmissions are dramatically increased and each generator is represented by the detailed EMT model with governor and exciter models; the dc grid becomes much more sophisticated due to the inclusion of multiple MMCs and the HFL-based DC-SST.The DCSST example is modified according to a practical project in Shenzhen, China.The project is now in the design stage.The parameters are given in Appendix.In China, the DC-SSTs will be used in distribution network.So, this simulation case will be helpful for the design and implementation of these projects.
Four simulation scenarios, i.e., asymmetric/symmetric ac system faults, dc pole-to-pole fault, and step change of the DCSST voltage, are investigated to check the performance of our proposed method in terms of efficiency and accuracy.For the scenarios of the ac system fault and step change of the DCSST voltage, the time steps of SFP and EMT subsystems are set as 20 μs; while for the scenario of dc pole-to-pole fault, the time steps are extended to 100 μs.
1) Scenario 2-1.A Three-Phase Resistive Fault: The fault is triggered at Bus SJ at t = 3.0 s and is cleared 0.05 s later.Bus SJ is the interface bus between the dc grid and the adjacent ac grid.DC quantities as well as capacitor voltages of converter no. 1 are shown in Fig. 15, which indicates that the results obtained with the proposed method match the reference values very well.The maximum errors of the dc voltage, dc current and capacitor voltage are no more than 0.0023, 0.04, and 0.001 p.u., respectively.Comparing the results in Fig. 15   in Fig. 12 has also revealed that the recovery of the dc current from the fault takes more time (about 0.8 s) when the generator is represented by its nonlinear model instead of the ideal voltage source, due to the response of the generator's governor control.
Figs. 16 and 17 display instantaneous and SFP curves of the ac voltages and currents.Clearly, SFP curves match the envelopes of instantaneous values exactly.The results also indicate that the system has two types of oscillation dynamics.One is the highfrequency electro-magnetic oscillation, lasting for one or two cycles following the disturbances and sometimes accompanied by significant current or voltage overshoots.This mainly originates from the frequency-dependent dynamics of the networks,    which is dominated by transmission lines.The other is the sustained low-frequency oscillations of small amplitude after the fault is cleared, as shown in Fig. 18.A deep investigation reveals that such low-frequency oscillation is caused by the dynamics of generators' and converters' control in response to power swings.Both types of dynamics can be accurately captured by SFP models.It should be noted that if researchers would like to study the frequency-dependent dynamics of transmission lines, their SFP models can be applied according to [28], which is regarded as a future research work.
2) Scenario 2-2.An Asymmetric Resistive Fault: An asymmetric double line (phase A and phase B of Bus SJ) to ground fault is trigged at t = 3.0 s and the fault lasts for 0.05 s.DC currents and capacitor voltages of the proposed model and PSCAD results are compared in Fig. 19.The maximum errors of dc cur-rents and capacitor voltages are less than 0.0052 and 0.003 p.u., respectively.
Figs. 20 and 21 demonstrate SFP-derived instantaneous and phasor curves of the ac voltages and currents.The following observations can be made: 1) SFP phasor curves match the envelopes of instantaneous values very closely under unbalanced fault conditions; and ii) the comparisons between Figs. 16 and 17,and Figs. 20 and 21 indicate that the system experiences sustained low-frequency oscillations of small magnitude after the clearance of the three-phase fault.However, the low-frequency oscillations are damped out after the double line to ground fault is cleared.Consequently, the proposed SFP model has satisfied the accuracy expectations under both balanced and unbalanced fault conditions, with additional benefit of producing wide frequency-band phasors.

3) Scenario 2-3. A DC Pole-to-Pole Fault:
A dc pole-topole fault is triggered at the dc bus between converter T3 and T4 from t = 2.0 s.The dc voltage drops immediately.The fault current grows rapidly.When it reaches the pre-set threshold, the protection actions to block all MMCs.As a result, all MMCs operate in the rectifier mode and the capacitor voltage of each sub-module remains unchanged [see Fig. 23(c)].At t = 2.1 s, the ac-side circuit breaker (CB) opens to isolate MMCs.After another 0.1 s, the dc fault is cleared.Next, the CBs are reclosed at t = 2.3 s and the ac voltages is restored.All MMCs are restarted sequentially from t = 2.5 s.
Fig. 22 displays the dc quantities during the dc fault.Other quantities obtained with different time steps are depicted in Fig. 23.For the proposed method, the maximum and mean errors under different time steps are listed in Table II, where the reference values are obtained by PSCAD/EMTDC with the same time step.As can be seen, the errors are reduced to less than 1e-3 p.u.Thus, it is well demonstrated that the proposed method can achieve satisfied accuracy even when the time step is extended to 100 μs.(15-40 kV).In this simulation scenario, the voltage of the MVdc bus is increased from 15 to 30 kV.Fig. 24 displays the primary and secondary voltages of DABs as well as their inductor currents.Fig. 25 shows the dc voltage of the MVdc bus.Again, the waveforms from the proposed method match very well with the reference curves.The maximum error of the MVdc voltage is below 0.03 p.u., indicating the high preciseness of the proposed method.

5) Comparisons of Simulation Efficiency:
In this paper, the simulation efficiency is measured in terms of execution time and   speedup.The speedup is defined as the ratio of the execution time consumed by PSCAD (base scenario) and the proposed method.The time step of PSCAD models is restricted to 20μs to guarantee the numerical stability of the simulation [9].To evaluate the simulation efficiency under different system scale, the number of MMCs is gradually increased from 1 to 8. All simulations are carried out on a computer with 2.67-GHz Intel i7 CPU, 8-GB RAM, and 64-b Windows 7 operating system.To make it fair to compare different simulation methods, the parallel processing technique is not adopted by any simulation method.The results are shown in Fig. 26.As can be seen, the execution time of the proposed method is much lower than PSCAD.The speedup is evident, especially when the time step of the proposed method is extended to 100 μs.As the scale of system or number of MMCs increases, the speedup grows as a whole.Note that when the number of MMCs increases to 8, the proposed method with the time step of 100 μs achieves a speedup of 90 times.
The reason why our proposed method can achieve such a high efficiency is given as follows.The calculation of the MMC by our method is actually a two-stage model [30].The first stage is to solve the global node voltage equation, where each arm of the MMC is represented by one Thevenin equivalent circuit according to (39).The number of nodes will not expand as the number of sub-modules increases, and only one node will be added in the node voltage equations.This node is the connecting node between upper and lower arms.In the second stage, the capacitor voltages are refreshed in a parallel and independent way [see (33)], which, taking the real and imaginary parts into account, is only a two-dimensional equation, and thus, can be very efficiently worked out.Similarly, the simulations of the DCSST can be efficiently calculated.Therefore, the proposed method has improved the efficiency significantly in simulating large-scale ac/dc grids.

VII. CONCLUSION
To accurately simulate large-scale ac/dc systems, especially to capture the interactions between the various components, this paper proposes a multi-domain co-simulation method.It can improve simulation efficiency considerably without sacrificing the accuracy.The co-simulation method starts by partitioning the target system into dc and ac subsystems, which are represented by our developed phasor-domain SFP models and traditional time-domain EMT models, respectively.SFP models accelerate the speed of simulating dc grids by using a larger time step.To reflect the wide-band dynamics between SFP and EMT models, a new interface model, namely HMD-TLM, is developed.It has an additional benefit of producing instantaneous and phasor waveforms simultaneously.Thus, the multi-domain co-simulation is implemented in a coordinated and decoupled The performance (efficiency and accuracy) of the proposed method has been validated on the modified CIGR ac/dc system and a practical ac/dc system.The results have demonstrated the following conclusions.
1) The curves obtained with SFP models match the envelopes of instantaneous values exactly.In other words, SFP models reflect both low-and high-frequency interactions perfectly.
2) With the proposed method, simulation errors are reduced to less than 1e−3 p.u. under different time-steps.
3) As compared with the traditional method, the proposed method has achieved a speedup of 90 times in simulating a practical system with a time-step of 100 μs.
Therefore, the proposed method offers an efficient and accurate option for analyzing very-large-scale ac/dc grids.
Authorized licensed use limited to: UNIVERSITY OF ALBERTA.Downloaded on April 18,2022 at 18:54:22 UTC from IEEE Xplore.Restrictions apply.mode, voltages of both the MVdc bus and the HVdc bus are fixed, where the amplitude and direction of the whole power of the DCSST are passively determined by the load at the MVdc bus.

Fig. 16 .
Fig. 16.Instantaneous and SFP curves of the ac voltages at Bus SJ.

Fig. 17 .
Fig. 17.Instantaneous and SFP curves of the ac currents at Bus SJ.

Fig. 19 .
Fig. 19.Simulation results/errors of dc currents and capacitor voltages (phase A) of Station 1.

Fig. 20 .
Fig. 20.Instantaneous and SFP curves of the ac voltages at Bus SJ.

Fig. 21 .
Fig. 21.Instantaneous and SFP curves of the ac currents at Bus SJ.

Fig. 22 .
Fig. 22. Simulation results of dc voltages and dc currents of Station 1.

Fig. 23 .
Fig. 23.Simulation results of arm currents (phase A), ac currents (phase A), and capacitor voltages (phase A) of Station 1 under different time steps.

4 )
Step Change in the DC Voltage of the HFL-Based DCSST: The HFL-based DCSST, composed of three independent DAB modules, transforms high ac voltage down to

Fig. 26 .
Fig. 26.Comparisons of execution time and speedup for different number of MMCs.
Authorized licensed use limited to: UNIVERSITY OF ALBERTA.Downloaded on April 18,2022 at 18:54:22 UTC from IEEE Xplore.Restrictions apply.on SFP models, HMD-TLM, and a designed time sequences of simulation.