Initialisation of a hybrid AC/DC power system for harmonic stability analysis using a power flow formulation

: This study presents the initialisation of a hybrid alternating current/direct current (AC/DC) power system to determine the operation point around which the system is linearised as a starting point for an electromagnetic or harmonic stability analysis. To this end, each AC or DC power system component is modelled in the Fourier transform domain and further adjusted to a simplified representation compatible with a power flow formulation. The study also presents how voltage source converter (VSC) high-voltage DC-based systems can be analysed when the different VSC controls are applied. The result of the power flow solution is then used for the power converters’ initialisation.


Introduction
Our power system more and more takes the form of a hybrid alternating current/direct current (AC/DC) network. This is caused by increased penetration of new installations of high-voltage DC (HVDC), in addition to a rise in converter-interfaced renewable energy sources (e.g. large wind farms). These changes have introduced the need for more accurate modelling and simulation of the entire network. As the name suggests, the resulting hybrid network consists of both the AC and DC systems interfaced with power-electronic converters. These devices take the form of line commutated converters (LCCs) or the currently prevailing voltage source converters (VSCs). For HVDC applications, the voltage source technology can either be a simple two-level or three-level topology or the more recently developed modular multilevel converter (MMC) [1][2][3].
The increased deployment of power electronics in the transmission system has given rise to new types of stability problems. These are commonly referred to as 'electromagnetic stability' [4] or 'harmonic stability'. A detailed description of the frequency-dependency of traditional power system components along with a detailed model of the converter's circuit equations and control loop layout is needed to analyse the problem. At present, electromagnetic or harmonic stability analyses are typically carried out by considering either the AC or DC side [5] and by focusing on relatively small surrounding areas in the AC power system [6,7]. Different approaches can be used to perform such a stability assessment, but the impedance-based method is amongst the most commonly applied ones [6][7][8]. However, to enable an assessment of the impact of AC characteristics on the DC side stability and vice versa, it is necessary to consider the entire network in the analysis and not limit the analysis to one subsystem or a smaller part thereof, as it is commonly done. Since the actual stability analysis relies on Fourier transforms, a first step in the analysis is to determine the operating point of the entire hybrid AC/DC system to perform a linearisation of the non-linear system component models (e.g. the converters and their control loops).
The introduction of power electronics in the system complicates the computation of operating states of the power system in general. Over the past few decades, there has been a constant need for conceiving new power flow (PF) models and problem formulations that account for hybrid AC/DC systems; first for LCC HVDC [9] and later for VSC HVDC. As an example, the MATLAB-based tool Matpower [10] was improved by adding converters, their overall control capabilities, as well as a DC power system analysis in [11]. Similarly, optimal PF (OPF) formulations started to include hybrid AC/DC systems. For example, in [12], the hybrid AC/DC system description from [11] was incorporated into an OPF description for AC systems in the Julia programming language [13]. In addition, alternative formulations have been presented in the literature including shunt conductances within the HVDC grids [14].
The initialisation of standard AC systems for dynamic simulations using a PF has been the subject of study in the past [15,16] and is at present available in commercial software packages. On the contrary, the new PF and OPF algorithms for hybrid AC/DC systems have typically been targeted towards operational studies and contingency analyses. Consequently, their integration with other power system simulation tools has not received much attention before.
The presence of power converters as coupling elements between both AC and DC systems complicates the formulation of accurate initialisation procedures for a hybrid AC/DC network. Using a PF formulation for this purpose implies not only extending existing AC system PF formulations to take into account the peculiarities of DC systems. Indeed, an accurate initialisation requires that steady-state characteristics of the converter and its controls be accurately mapped to equivalent power-flow representations [17]. As a result, accurate initialisation routines are only starting to see the light of day as of recently. This is particularly true when the more complex MMC topology is concerned. For example, the design of an accurate initialisation procedure for an MMC was proposed in the work [18]. However, the presented algorithm focuses on determining the steady-state operating point of the MMC, while its interaction with the rest of the AC and DC power systems is left out of the study. An iterative AC/DC PF is used in [19] for estimating the MMC's operating point together with the system state to perform transient studies. Consequently, a steady-state phasor network model is used for the AC system [20]. The DC side network dynamics are considered, but only a simple equivalent π -section model composed of constant lumped R, L, and C parameters is used for the cable. Moreover, the approach does not consider meshed complex AC/DC networks.
This paper aims at filling this gap by introducing a detailed PFbased initialisation of complex hybrid AC/DC power systems to allow large-system electromagnetic or harmonic stability assessments using complex frequency-dependent component models. The power system contains both HVDC and AC subsystems, interconnected using state-of-the-art MMCs. The operating point analysis starts from defining equivalent PF models of frequency-dependent passive components (such as cables and overhead lines (OHLs)) together with a simplified equivalent steady-state model of the MMC incorporating the converter's outer controlling loops. After solving the hybrid AC/DC PF, the steadystate operating conditions of the network are used as a starting point for determining the operating conditions of each individual power electronic converter, by initialising inner and outer controlling loops. The developed procedure is described in the flow chart from Fig. 1. The practical implementation of this procedure ensures a short simulation time because (i) the system's equations are decoupled, (ii) complex component models are reduced to fit a PF formulation, and (iii) it is implemented in Julia, a C-related programming language. By making abstraction of inner converter control details, while solving the operating conditions of the network, the method also provides high flexibility in defining the converter's inner and outer controlling loops.
The methodology described in this paper is implemented in the Julia programming language, which offers a wide range of functionalities and fast computation. It includes the definition of each component's PF equivalent to fit the PF formulation from [13], using the package PowerModelsACDC.jl. The MMC is implemented based on the detailed model proposed in [21] since it is the most prevailing technology for VSC HVDC applications.
The paper is organised as follows. Section 2 reviews the basics for a multiport network implementation using ABCD and Y parameters and discusses the conversions between them. The section thereby explains the conversion between the internal modelling method used to construct equivalent grid impedances for the eventual stability assessment as in [22] and the format needed for the PF analysis. Section 3 describes the PF formulation with the corresponding models for the most commonly used power system components. The proper operation of the proposed PF initialisation algorithm is demonstrated in Section 4 on a detailed hybrid AC/DC power system consisting of a multi-terminal HVDC-based power system and small-scale adjacent AC networks. Section 5 presents the concluding remarks.

Notation
Uppercase bold identifiers are used for matrices, e.g. X, and lowercase bold identifiers are used for vectors, e.g. x. Matrix I n presents the n × n identity matrix. The following definition is used for a compact set of equations: For any matrix X, the notation X i, j presents an element of the matrix at the location i, j . The notation X i: j, k: l is used for a submatrix with j − i × l − k elements formed from a portion of matrix X. The imaginary unit is j = −1. The Fourier transform's variable s presents s = jω, for an angular frequency ω. ℜ x and ℑ x present real and imaginary part of the complex number x ∈ ℂ.
To transform three-phase voltages and currents from the stationary abc-frame to the rotating dqz-frame, the Park's transformation is used with ω 0 the angular frequency. The inverse Park's transformation is given as (2)

Power system components
For the system description at the basis of the electromagnetic stability assessment, each power system component is considered as a multiport network that can be described in the Fourier domain using ABCD, Y, Z, H, or S parameters. In this work, we have chosen to employ the definition using ABCD and Y parameters since they preserve the physical dimension of the component's voltages and currents. Assuming that power system components have the same number of input and output nodes-which is correct for all analysed three-phase AC and DC components as well as for both monopolar and bipolar DC power network configurations-each component can be represented using ABCD parameters as In (3), the current and voltage vectors at the component input (or primary nodes) have superscript 'p', while the output currents and voltages have superscript 's'. Matrices A, B, C, and D represent the interconnection between input and output variables (currents and voltages) specific for the power system component. If the number of input/output nodes of the component is n, then these matrices have sizes n × n. Similarly, the Y parameters for the same component can be defined as The Y parameters in the previous form can be obtained from ABCD parameters as follows [23]: Similarly, the ABCD parameters can be determined from Y parameters in a unified manner as

PF definition
To perform an electromagnetic or harmonic stability analysis, nonlinear descriptions of the components, such as power converters, must be linearised around an operating point calculated for the entire network. This operating point is determined in the initialisation phase by solving the PF equations representing the combined AC/DC system. To do so, the network is initialised using the OPF tool from [12] implemented as a Julia package, which can be found in the PowerModelsACDC.jl repository [13]. The package relies on the PF models developed for MatACDC [11], which extends the Matpower [10] AC power system models with DC system representations and with power converters. The mentioned PowerModelsACDC.jl library contains models for AC and DC power systems description provided in Table 1.
As a result, the constructed power system is divided into AC systems, DC systems, and converters. It contains AC and DC branches and buses, converters, generators, loads, and shunts. Storage elements can be represented as equivalent injections for the particular operation point. Components included for the electromagnetic stability assessment thus need to be represented using their equivalent models for the PF analysis to obtain the overall system description in terms of system voltages and currents. It should be noted that all models in this paper are developed under the assumption that the power system is balanced.
In this section, the fundamentals of the components used in the multiport power system are described. Only the most common components are considered, being impedance, transformer, transmission line, cable, AC and DC sources/grids, and given the focus on the electromagnetic instability assessment, the MMC converter. Components are presented using their physical equations, as well as by their model descriptions using ABCD parameters in the Fourier transform domain.

AC and DC branches
AC and DC branches typically represent three-phase AC and DC connections between AC and DC buses, respectively. Branches are grouped inside different AC or DC grids when several AC or DC networks are present. AC branches are typically defined in PF routines using their circuit parameters evaluated at 50 Hz, as described in [10], whereas DC branches can generally be reduced to their series resistance at DC [11].
The model of the AC branch as provided in [10] is depicted in Fig. 2. Beside the shunt susceptance j(b c /2), the Julia package PowerModelsACDC.jl [13] supports also a resistive shunt admittance which can also be unsymmetrical with respect to the sending and receiving side. In our formulation, we use a symmetrical resistive and capacitive shunts (g c /2) + j(b c /2). The full expression for the AC admittance parameters is given by the equation: It should be noted that the AC network is considered as balanced for the PF solution, and thus, all the components have diagonal matrix models, with equal elements on the matrix diagonal.
The DC branch is modelled with its equivalent series resistance [11], as depicted in Fig. 3.
For the initialisation using a PF calculation, a number of components need to be modelled either as AC or DC branches. Their models are described in detail in this subsection.
3.1.1 Impedance: Impedance can be formed between n input nodes and n output nodes. The impedance is defined with the n × n matrix Z, where (8) and each component Z i j represents the impedance between input node i ∈ 1, …, n and output node j ∈ 1, …, n .
The ABCD parameters of the impedance from matrix (8) are given by In the case of a DC resistance, all matrices are of size 1 × 1, while three-phase impedances are of size 3 × 3.
The corresponding DC branch model that is used for this component's representation is given as a Thevenin equivalent series impedance r = ℜ Z . The AC branch is modelled according to  It should be noted that this representation only refers to impedances interconnecting several nodes in the system (e.g. simplified grid representations by means of a Thevenin equivalent impedance). AC shunt impedances should be treated separately as shunt components. DC shunt resistors can only be added as DC loads.

Transformer:
The single-phase equivalent of a transformer is depicted in Fig. 4, which presents an accurate model for the transformer operating below the MHz range, as recommended in [24]. It relies on the model as detailed in [4,25] with the ABCD parametric description where The single-phase realisation can be extended into the YY and ΔY three-phase configuration [4].
Three-phase transformers are accounted for either in a YY or ΔY configuration, where each of the three single-phase transformers is represented by its equivalent from Fig. 4.
The YY configuration is derived from (10), such that The ΔY configuration is more complex and it is modelled using the following equations. The inner primary and secondary parameters of the transformer (i.e. all the components except the parasitic capacitances and the load impedance) are given by The ΔY configuration transforms voltages from the delta side v p Using the previous voltage/current relations and ABCD representation of the inner transfer function in (11), the inner impedance can be obtained as The transformer is now represented using ABCD parameters as Since the transformer model cannot easily be represented as the model in Fig. 2, Y parameters are extracted from the ABCD parameters as described in (5).
Under the assumption of a balanced system, the submatrices

Transmission line:
Transmission lines can take the form of cables, OHLs and mixed OHL-cable systems. Since cross-bonded cables [25] and transposed OHLs do shift the resonance spectrum and hence the electromagnetic stability assessment, they are considered here as well. All transmission lines have ABCD model parameters, which can be estimated as in [26,27] using known values for the series impedance Z′ and the shunt admittance Y′ per unit length as where Γ = Z′Y′, Y c = Z′ −1 Γ, and l is the line length. The series impedance Z and shunt admittance Y are derived differently depending on the actual OHL or a cable realisation. For the electromagnetic stability assessment, the models are derived as frequency-dependent distributed parameter models.
A transmission line (OHL, cable, cross-bonded cable, or mixed OHL-cable) is represented using its nominal π -model depicted in Fig. 5, where For the DC lines and cables, the shunt admittance is not considered, while the branch resistance is equal to r = ℜ Z 0 . In the symmetrical monopolar configuration, the transmission line (usually underground cable) is represented with two resistances r. This model can easily be extended to employ shunt conductance as in [14]. For the balanced AC transmission line, the impedance and admittance matrices are diagonal such that Z jω = Z jω 1, 1 and Y jω = Y jω 1, 1 can be chosen. Then, the AC branch PF model is given by

Shunt components
Shunt reactors and capacitors are defined with their admittance value as y = g s + jb s [10]. Shunt elements are modelled according to [25].

Generators
Models of the generators are represented with the desired values for the voltage magnitude V m and phase θ, and with an estimation of the active power P and reactive power Q injections. In this work, they are used for modelling both grid equivalents and AC power sources (hence a simplified representation of synchronous machines), implying a high degree of idealisation. Generators are thus defined as reference buses. Grid equivalents are thus represented in the PF formulation as generators behind complex impedances, using the AC branch model introduced previously.

Power electronic converters
A power electronic converter is modelled in accordance with [11,28], including its phase reactor, potentially a filter, and a transformer. To match it with the MMC model, only the reactor is considered as part of the converter model. Internal losses of the converter for the PF formulation are calculated in the form of P loss = a + bI c + cI c 2 . The phase reactor, however, contains a series resistance and reactance as z pr = R r + jX r , where X r = ωL r .
Depending on the actual realisation of the converter's controls, the parameters of the converter can control DC voltage, the converter active power or use a DC side voltage-power droop. The model of the power converter with the corresponding default PF directions is depicted in Fig. 6, where the AC node and parameters are presented in black and DC node is presented in blue.

MMC power converter
As mentioned before, the focus of this work is the integration of the MMCs inside a hybrid AC/DC power system. The model incorporated relies on the work presented in [21,29] and is depicted in Fig. 7, where each submodule is represented as a halfbridge module. The switches inside submodules are considered as ideal in this model and thus, losses of the converter are given only with the constant c = R arm /2 .
Submodules are represented with their averaged equivalent, and thus, the following equations for voltage and current can be written for the upper and lower arms: where m j U, L are the corresponding upper and lower arm insertion indices.
The converter model is developed by following the methodology reported in [21,29,30]. Using the Σ − Δ nomenclature, the variables in the upper and lower converter arms can be determined as The MMC is described using differential equations where N being the number of submodules, C the capacitance of the halfbridge submodule. L arm and R arm are the equivalent inductance and resistance of the arm, respectively. The expressions for the equivalent inductance and resistance are L eq ac = L r + (L arm /2) and R eq ac = R r + (R arm /2). The model incorporates dqz-frames at the angular frequency −2ω components for the Σ variables and components at the frequencies ω and 3ω for Δ state variables, in accordance with [21]. The insertion indices m dqz Σ T are obtained by using the following algebraic relation: The voltage references in (20) are obtained as an output of the controllers, which accompanies the set of 13 differential (18)  , which is set to (v dc /2) by default, the other voltage references are set to zero unless they are specified by a certain controller.
This MMC model can easily be incorporated with different controlling loops. Supported controlling loops can be grouped into inner and outer controlling loops, see Fig. 8. Each control loop is implemented as proportional-integral (PI) controller. Controllers for the output current, circulating, and zero-sequence current are tuned using a pole-zero placement method [31].

Numerical simulations
This section describes a numerical example for the case of a hybrid AC/DC power system including an MMC-based HVDC grid, shown in Fig. 9. The figure shows the hybrid AC/DC system, with the DC system indicated in blue, while the AC is depicted in black. The example is considered in the publication [32] to analyse the design of a passivity-based control. The analysed example contains eight AC nodes and four DC nodes. The converters' parameters are provided in Table 2. Power electronic converters have control loops for circulating current, zero-sequence current and output current control, together with the zero-energy balancing control. MMC 1 controls the DC voltage, while MMCs 2-4 have control loops for active and reactive power with the PI coefficients as given in Table 2.
AC branches, denoted as tl i for i ∈ 1, 4 , are all implemented as OHLs with the parameters provided in Table 3, while the DC branches, denoted as c i for i ∈ 1, 5 , are a pair of two underground cables ( Table 3). The OHLs are considered as single-pole equivalents using (16), while the underground cables are in a symmetrical monopolar configuration.
It is desired that the MMC 1 receives an active power of P 1 = 120 MW and reactive power Q 1 = 0, respectively. The power controlling converters are set to control the active power injected/ absorbed to/from the AC grid: P 2 = 60 MW, P 3 = 20 MW, and P 4 = 40 MW, and to have no reactive power injection, Q 1 − 4 = 0. Power converter MMC 1 ensures that the HVDC network's voltage at DC node 1 is V dc = 640 kV, while it is desired that the voltage of AC buses 1, 4, 5 and 8 have a magnitude of V m = 380 kV and a phase angle of zero, thus the wind farm and the AC grid connections are represented as reference nodes of the four asynchronous areas.
The determined operating point for the converters is given in Table 4. These values are further used as a reference for solving the equations for the converter's equilibrium.
The following equations are used to initialise non-linear equations with the result of the initialisation as input The operating points (equilibrium) of the converters are presented in Table 5 for all 13 state variables. By comparing the values from Tables 4 and 5 we can see the great correspondence between them, which demonstrates that the operating point of an MMC can be estimated in two steps: first by solving the PF problem with the desired goals and then by solving the set of converter's non-linear equations for determining the equilibrium for all state variables. The solution time for this example is 29.141 s, including the formulation of the PF problem, solving it, and determining the operating point of each converter on an Intel(R) Core(TM) i7-7700HQ on 2.80 GHz CPU and with 16 GB RAM. The solution of the PF problem itself takes 54 ms, while the remaining time goes to solving the non-linear equations for the four converters' steadystate estimation and the complete system initialisation. The PF is solved as an AC/DC PF problem using PowerModelsACDC.jl [13] with Ipopt as a solver. Solving the non-linear equations of the MMC is done using the NLsolve package written in the Julia programming language [33].

Conclusion
This paper presented the use of a PF formulation to determine the operating point of the hybrid AC/DC power system, as a starting point for an electromagnetic or harmonic stability assessment. The      simulation is performed in two steps since it couples procedures for the AC/DC power system's PF and the solution of non-linear differential equations for the operating point of an MMC. In the first step, the PF problem of the combined power system is solved, in which power converters are modelled only with their losses and controlling loops for DC voltage and active and reactive power. The PF solution is used in the second step to set the initial values of the power converters, and then the complete solution for the converter's equilibrium is found with respect to the initial values. This procedure shortens the time for solving the PF problem, as it avoids including the equations belonging to the converter's internal dynamics and controls at the beginning. At the same time, it provides sufficient accuracy, which is demonstrated by the fact that the full order converter's equilibrium can easily be determined. The proposed procedure ensures a high degree of flexibility in the definition of the power system as well as for the inner and outer controlling loops of the MMCs.

Acknowledgment
This work is part of the Neptune project, supported by the Energy Transition Fund, FOD Economy, Belgium. This work was supported by the Research Foundation Flanders (FWO) under grant no. G0D2319N. The research of J. Beerten was funded by the FWO under grant no. 12D1117N.