General Power Flow Calculation for Multi-Terminal HVDC System Based on Sensitivity Analysis and Extended AC Grid

The power flow calculation of the multiterminal high voltage direct current (MTDC) system is essential for planning sustainable energy sources and power flow analysis for the MTDC system. However, the traditional unified methods require a large system of non-linear equations leading to low calculational efficiency. Also, for a large DC grid, there is a concern about the convergence of sequential methods. This paper proposes a general power flow calculation method for voltage source converter (VSC) based MTDC systems. Based on an extended topology of an AC grid, a generalized calculation model of the MTDC power flow is proposed. Then, a novel sensitivity analysis-based power flow (SAPF) method is proposed, in which the state variables of the extended AC grid are calculated via sensitivity analysis. With a smaller system of equations, the proposed SAPF method has less computational burden than the traditional unified methods and causes no convergence problem compared to sequential methods. To further improve the calculational efficiency, the sensitivity-based variable updating is adopted to accelerate the iterative process. By comparing with the existing methods in calculating power flow for different MTDC systems, including large systems with multiple AC/DC grids, the effectiveness and scalability of the proposed methods are verified.


I. INTRODUCTION
W ith the advantages of large transmission capacity and low power loss [1], the high voltage direct current (HVDC) technology has been applied in offshore wind power systems [2] and power transmission between asynchronous grids [3]. However, the HVDC transmission is usually with the back-to-back and point-to-point configuration [4], which limits the application of the HVDC system. With the development of the direct current (DC) circuit breaker and voltage source converter (VSC), the interest of the multiterminal high voltage direct current (MTDC) system increases. As a promising power transmission network, the VSC based MTDC system can realize the connection of multiple asynchronous power grids, DC sources and loads, which shows flexible power control ability and wide applications [5].
The power flow calculation of the MTDC system is the foundation of the steady-state studies of the MTDC system, such as adaptive power control [6] and static security analysis [7]. To calculate the power flow of the MTDC system, the outputs of DC sources, the power control of multi-VSCs and the power loss of VSCs need to be considered. Generally, the existing methods for power flow calculation of the MTDC system can be mainly divided into two categories: the sequential methods [8]- [13] and the unified methods [14]- [19].
In sequential methods, the power flow for the alternating current (AC) and DC grids is calculated separately and repeatedly, where the state variables of each VSC are set as coupling variables for determining the consistency of power flow at the AC side and DC side of the MTDC system. The main advantage of sequential methods is that the existing AC power flow algorithms can be adopted without any modifications. As presented in [8], the traditional Newton-Raphson (NR) method is used to solve the AC power flow efficiently, and the DC state variables of the MTDC system are iterated until convergence. In [9] and [10], the sequential method is applied for solving the power flow of the MTDC system under droop control. A double-layer loop iterative method is proposed in [11] to improve the convergence of the sequential method where the approximation is introduced into the calculation. In [12], the MTDC system is divided into smaller sub-systems and analyzed separately to improve computational efficiency. Also, to further improve the computational speed, the iteration of the sequential method is simplified by approximating the VSC current in [13]. In unified methods, solving the power flow of the MTDC system is to solve a large system of nonlinear equations which can be built by combining all the power flow equations of AC grids and DC grids. The results of the unified methods are generally with high accuracy as no assumptions are included in the unified method. In [14], the system of equations of the power flow calculation for a common MTDC system is built. Also, considering the power loss of the VSC, a general unified method for MTDC is presented in [15] where the mismatch equations are built with a general model of the VSC station. The NR method is adopted in [16] and [17] to solve the large system of equations built in the unified method where the DC slack control and droop control are considered. To further accelerate the calculation, the smooth This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ approximation of the Fischer-Burmeister is used in [18], but the computational accuracy is reduced by the approximation. In [19], the performance of the typical unified method is tested based on simulations where the precision of the unified method is verified.
However, in certain situations, sequential methods may lead to the convergence problem [14]. For example, when the MTDC system contains multiple VSCs with the droop control or the system is with multiple AC and DC grids, multiple coupling variables need to be determined by the iterative process which does not converge until all coupling variables converge [11]. Therefore, with the increase of the number of VSCs, it would be a challenge to make a large number of coupling variables converge at the same time [12]. Although the approximation is adopted in [11] to improve the convergence of the sequential method, approximations will reduce the calculation accuracy. In unified methods, although there is no need to set coupling variables, it is computationally inefficient to solve a large system of equations, which limits the scalability of unified methods for optimal power flow calculation [20] and dynamic analysis [21] of a large system. Also, when the number of buses and VSCs increases, the computational burden of the unified method will increase nonlinearly and dramatically [10]. Besides, the modification of power flow equations in unified methods introduce more variables into the NR iteration [18], which could lead to the difficulty of finding suitable initial values for the convergence [22]. Therefore, to improve the previous methods, this paper proposes a novel power flow calculation method for the MTDC system. The main contributions of this paper are as the following, 1) Based on an extended topology of the AC grid, a general calculation model of the MTDC power flow is presented. One advantage of this topology is that the calculation equations of the AC grid do not need to be modified. Also, the power flow calculation model is built by combining the sub-model of each DC bus, which can be generalized to multiple large DC grids. 2) A novel sensitivity analysis-based power flow (SAPF) method is proposed to calculate the power flow for MTDC systems, in which the state variables of the extended AC grid are calculated via sensitivity analysis. Thus, the power flow calculation model can be represented by a small-scale system of equations, contributing to less computational burden and higher calculational efficiency compared to unified methods. In addition, the NR method is used in the SAPF avoiding the convergence problems occurred in the sequential method. 3) To further improve the calculation speed of the SAPF method, a fast SAPF (F-SAPF) method is proposed, in which the AC power flow is obtained by the sensitivity analysis instead of the traditional AC power flow calculation. The results obtained by the F-SAPF method are as accurate as that of the SAPF method because of the unchanged mismatch equations. 4) Different MTDC systems, including the scenarios of multiple DC grids and large AC networks, are studied to verify the effectiveness of the proposed method whose performance is compared with other existing methods.  Also, the power flow variations and contingencies in an offshore wind power system are investigated by the proposed method. The rest of the paper is organized as follows. Section II presents the extended AC grid and the power flow calculation model of the MTDC system. Then, the SAPF and F-SAPF methods are proposed in Section III. Case studies are implemented in Section IV, and Section V concludes the paper.

A. The MTDC System and the VSC Station
The structure of the generalized MTDC system is presented in Fig. 1, where the DC grid includes DC transmission lines, VSCconnected DC buses (namely VSC buses) and pure DC buses (the DC buses without any VSCs connected). The DC sources denote DC power supplies such as photovoltaic (PV) power stations and battery storage systems. A complex MTDC system normally consists of multiple AC and DC grids, as shown in Fig. 1. Fig. 2 shows the steady state model of the VSC station connecting the k-th AC bus in the AC grid and i-th DC bus in the DC grid, where R T,i and X T,i are the equivalent resistance and reactance of the transformer, respectively. R c,i and X c,i are the equivalent resistance and reactance of the phase reactor respectively. Also, P k and Q k denote the active and reactive power injected to the gird at the k-th AC bus respectively. The magnitude and angle of the voltage at the k-th AC bus are V k and θ k respectively, and V f,i , θ f,i denote the magnitude and angle of the voltage at the node f i respectively.
According to [8], [23] and [24], the total loss of the VSC can be described as where a 0,i , a 1,i , a 2,i are three loss coefficients of the VSC connected to the i-th DC bus. I c,i is the magnitude of the converter current which can be calculated as where P f,i , Q f,i , S f,i are the active power, the reactive power and the magnitude of the apparent power, respectively, from the node f i to c i . And V f,i is the magnitude of the voltage at node f Fig. 2.

B. The Extended AC grid
The nodal balance for each AC bus in an N-bus AC grid can be described as where P z and Q z denote the active power injection and reactive power injection at the z-th bus (z = 1, 2, · · · , N) respectively. V z and θ z denote the magnitude and angle of the voltage at the z-th bus respectively. In (3) and (4), G zj and B zj are the real part and imaginary part of the admittance of the branch between the z-th and j-th AC bus respectively. So, by applying (3) and (4) to all buses in the N-bus AC grid, a system of 2N equations can be built to calculate the AC power flow. Then, we suppose this AC grid is connected to an MTDC system and one of the original AC buses, namely the k-th AC bus, is connected to the i-th DC bus via a VSC station, as shown in Fig. 2. In this paper, an extended AC grid is built by combining the transformer, filter, node f i and the N-bus AC grid. This extended model of the AC grid is different from the model presented in [11] and [19] in which both node f i and c i are considered as the additional AC buses. Also, in [11] and [19], Equations (1)-(4) are all needed to build the nodal balance model for node c i in each VSC station, which causes complexity to build the AC power flow model and reduces the computational efficiency of the AC power flow calculation. In the extended model proposed in this paper, only the node f i is regarded as an additional AC bus to the original AC grid, and the filter is regarded as an additional susceptance of the transformer branch between the node f i and the k-th AC bus. In this way, the node f i can be regarded as the (N+1)-th AC bus of this extended AC grid (which has N original AC buses). Without adding (1) and (2) to the AC power flow model, (3) and (4) are adequate for building the nodal balance model at the (N+1)-th AC bus. Note that, considering the direction of P f,i and Q f,i shown in Fig. 2, the active power injection at the (N+1)-th AC bus is P N+1 = −P f,i and the reactive power injection is Q N+1 = −Q f,i .
Normally, for each AC bus (including additional AC buses), the four variables P z , Q z , V z , θ z shown in (3) and (4) can be classified as input variables and state variables. Table I shows the input and state variables of the z-th AC bus according to its bus type, where two input variables are referred to as u 1,z , u 2,z and two state variables are referred to as y 1,z , y 2,z . Table I can be applied to each AC bus in the extended AC grid including additional AC buses. Note that the bus type of each additional AC bus depends on the reactive power control of the connected VSC, which will be presented in Part C of this Section.

C. The DC Power Flow Calculation and Power controls
The active power P dc,i from the VSC to the i-th DC bus in an n-DC-bus MTDC system can be presented as where P G,i is the power injection from the DC sources to the i-th DC bus. Also, P L,i is the DC load and U dc,i is the DC voltage at the i-th DC bus. Y dc ij is the admittance between the i-th and j-th DC bus in this MTDC system.
The control of a VSC is with the corporation of an inner d-q current control and an outer power control, in which the active power and reactive power are controlled independently [25]. According to [3], [13], [25] and [26], the control objectives of the outer power control are the variables of the AC side, including P f,i , Q f,i , V f,i of the VSC station as shown in Fig. 2. From a steady state point-of-view, the outer power control (reactive and active power control) is the key to calculating the power flow of the MTDC system. 1) Reactive power control.
One of the following strategies is generally used to control the reactive power of a VSC: f,i and the node f i is regarded as an additional PQ bus of the extended AC grid); f,i and the node f i is regarded as an additional PV bus of the extended AC grid.
The superscript * denotes the reference value of control objectives. According to the model of the extended AC grid, the node f i is numbered as the (N+i)-th AC bus in the extended AC grid. As shown in Table I, with both the P-Q mode and P-V mode, the controlled variable (Q f,i or V f,i ) can be referred to as u 2,N+i . So, the reactive power control can be described as 1) Active power control (also known as DC voltage control). In this paper, the two most commonly used modes of DC voltage control are introduced: a) DC slack control: one of VSC buses is selected as the slack bus whose DC voltage is controlled to be constant, as shown in Fig. 3(a), and the other VSC buses are the non-slack DC buses whose P f,i is controlled to be constant, b) Droop control: for each VSC, either P f,i or I dc,i is controlled to be constant, or where K i is the droop coefficient of the VSC connected the i-th DC bus, and I dc,i is the DC current from the VSC to this DC bus shown in Fig. 2. The first type of droop control is called droop voltage-power (U-P) control, as shown in (9) and Fig. 3(b). The other one is called droop voltage-current (U-I) control, as shown in (10) and Fig. 3(c). According to (10), P dc,i at the i-th DC bus can be calculated by

D. Power Flow Calculation Model of the MTDC System
The power flow calculation model for an n-DC-bus MTDC system should be built on the basis of different control modes. Suppose the i-th DC bus is 1) a pure DC bus: P dc,i = 0 as the pure DC bus is not connected to the AC grid, so that an equation, namely f 1,i , for this pure DC bus can be obtained by substituting P dc,i = 0 into (5) as So, the following simultaneous equations F s,i can be built: where

2) a VSC bus with the DC slack control:
The active power balance at node f i can be described as By substituting (1) and (5) into (14), the f 1,i of the i-th DC bus can be built as Then, another equation, namely f 2,i , of the i-th DC bus can be obtained from (2): As the node f i is considered as the (N+i)-th bus (PQ or PV bus) of the extended AC grid, the set {V f,i , P f,i , Q f,i } can be referred to as the set {u 2,N+i , u 1,N+i , y 2,N+i } according to Table I. Also, according to (6), u 2,N+i is controlled to be the pre-set value. Therefore, the independent variables of (17) can be summarized as {I c,i , u 1,N+i , y 2,N+i }.
By combing the (15) and (17), the following simultaneous equations can be built for the i-th DC bus: 3) a VSC bus with the droop control: Two equations f 1,i and f 2,i for the i-th DC bus are as the same as (15) and (17) respectively. Then, another equation, namely f 3,i , is added. The f 3,i can be built according to the type of the droop control: if the connected VSC is under droop U-P control, by (9), or under droop U-I control, by substituting (11) into (5), Similarly, the following simultaneous equations can be built for the i-th DC bus with droop control: According to (13), (18) and (21) The sensitivity analysis is an effective tool to analyze the state of power systems, which has been used in the power flow coordination [27], line losses determination [28], convergency analysis [29], optimal power flow calculation [30], etc. In this Section, the sensitivity analysis is adopted to calculate the power flow for the MTDC system.

A. Sensitivity Analysis for the Extended AC Grid
Suppose an N-bus original AC grid is connected to an n-DC-bus MTDC system which contains m VSC buses. Thus, an extended AC grid containing N+m AC buses can be built by adding m additional AC buses (the node f i in each VSC station) to the original AC grid. By applying (3) and (4) to all buses in the extended AC grid, the following simultaneous equations can be built to solve the AC power flow: where u and y are the vectors consisted of input variables and state variables of all AC bus respectively. According to Table I, u = [u 1,1 , u 2,1 , u 1,2 , u 2,2 , · · · , u 1,N +m , u 2,N +m ] T and y = [y 1,1 , y 2,1 , y 1,2 , y 2,2 , · · · , y 1,N +m , y 2,N +m ] T . Normally, the AC power flow calculation is to solve y by inputting u into (22), which can be described as where G( · ) denotes the process of the AC power flow calculation. Suppose the extended AC grid is on a state with [u (k) , y (k) ], by the first-order Taylor polynomial, where Δu and Δy are the variation of input variables and corresponding state variables. By substituting (22) into (24), where S ac = Δy/Δu is called the sensitivity matrix whose entries describe the relationship between the variation of each state variable and input variable, which can be expressed as where B( · ) denotes the process to calculate S ac . Based on the data of the N-bus original AC grid, the input variables (u 1,1 , u 2,1 , u 1,2 , u 2,2 , · · · , u 1,N , u 2,N ) of N original AC buses are normally pre-known. Instead, for the m additional AC buses, the input variables (u 1,N +1 , u 2,N +1 , u 1,N +2 , u 2,N +2 , · · · , u 1,N +m , u 2,N +m ) are unknown. According to the reactive power control of each VSC as shown in (6), u 2,N +i = u * 2,N +i (i = 1, 2, 3, · · · , m) for these additional AC buses. Thus, the unknown variables in u can be summarized as u 1 = [u 1,N +1 , u 1,N +2 , · · · , u 1,N +m ] T , which means u can be determined once u 1 is given, namely u = g(u 1 ). By rewriting (23), the AC power flow calculation for the extended AC grid can be presented as Similarly, (27) can be rewritten as

B. The SAPF Method
The basic idea of the SAPF method is to build a small-scale system of equations based on the equations of DC buses, and the mismatch equations of the NR method are modified by using the sensitivity matrix. Taking the n-DC-bus MTDC system mentioned in Section III Part A as an example, suppose the m VSC buses are with the DC slack control and the DC-1 (the 1st DC bus) is set as the DC slack bus. For each DC bus, simultaneous equations F s,i can be built based on its type, as shown in Section II Part D. Thus, by combining F s,i of each DC bus, the power flow calculation model of the MTDC system can be built, namely F dc , which is a system of equations shown as where I c = [I c,1 , I c,1 , · · · , I c,m ] T and y 2 = [y 2,N +1 , y 2,N +2 , · · · , y 2,N +m ] T . Then, a new matrix S MT DC =Δy 2 /Δu1 is set, whose entries can be extracted from S ac as shown in (31).  (30), the NR method is adopted, and the mismatch equations of (30) can be presented as Then, by substituting S MT DC =Δy 2 /Δu1 into (32), Fig. 4. The sensitivity analysis implemented in the SAPF method.
where X = [I T c , U T dc , u T 1 ] T , and J denotes the Jacobi matrix. According to (33), the sensitivity analysis is implemented between NR iterations for updating y 2 , as shown in Fig. 4. Note that, to reduce the computational burden, y obtained by each AC power flow calculation (represented by (28)) is adopted to be the initial value of the next AC power flow calculation in the sensitivity analysis, as the process represented by the green line in Fig. 4. The detailed steps of the SAPF method are presented in Fig. 5 where ε is the convergence tolerance.
Similarly, the calculation model presented in (30) can also be built for the MTDC system with the droop control, where the simultaneous equations F s,i of the VSC buses can be obtained by (21). Then, the SAPF method can also be implemented by following the steps shown in Fig. 5. For the MTDC system consisting of multiple DC grids, the calculation model F dc shown in (30) can be built by combining the simultaneous equations F s,i for all DC buses from different DC grids. To solve F dc , the SAPF method are also applicable by following the steps shown in Fig. 5. The multiple DC grids should be seen as a whole for building (5), where the disconnection of multiple DC grids can be represented by setting the admittance between disconnected DC buses to zero. For the MTDC system consisted by multiple AC grids (namely AC-1, 2, · · · , M), multiple S ac (S AC−1 ac , S AC−2 ac , · · · , S AC−M ac ) and corresponding S MTDC (S AC−1 MT DC , S AC−2 MT DC , · · · , S AC−M MT DC ) can be generated. Therefore, the SAPF method can be implemented as Δy 2 in mismatch (33) is calculated by are the corresponding u 1 in AC-1, 2, · · · , M. Thus, there are multiple sensitivity analysis processes for multiple AC grids in the SAPF method. Once the number of AC buses in one of the AC grids increases, it will only increase the computational burden in the sensitivity analysis for this AC grid, where the sensitivity analysis for other AC grids is not affected. While, in the unified method, the increased number of the AC buses will directly increase the scale of the large simultaneous equations for the whole system, which will significantly increase the computational burden of the unified method.
For the MTDC system with the DC slack control, the number of equations in (30) is m+n, and this number will be 2m+n for the MTDC system with the droop control. By introducing the control equations into the calculation model represented by (30), the power flow of the MTDC system with other control methods such as voltage margin control [31] and droop control with dead-band [32] can also be solved by the SAPF method. In addition, the model of the extended AC grid can also be built for the line-commuted converter (LCC) stations which is similar to the VSC stations except for the inner structure of the converter [33]. The sensitivity analysis shown in Fig. 4 can be used to determine the state variables of the extended AC models for both the LCC and VSC side. Therefore, the SAPF method can potentially solve the power flow problem for the hybrid LCC-VSC HVDC system, where the corresponding control strategies [33] of LCC stations will be applied to build the calculation model in (30).

C. The F-SAPF Method
In the iterative process of the SAPF method, the AC power flow calculation represented by (28) is implemented repeatedly to update y. To further improve the calculational efficiency, a fast SAPF (F-SAPF) method is proposed in this part. In the F-SAPF method, y is updated by sensitivity analysis instead of the AC power flow calculation, which can be described as y (k+1) = y (k) + Δy (k) = y (k) + S y · Δu (k) 1 (35) where S y = Δy/Δu1, whose entries can also be extracted from S ac : The detailed steps of the F-SAPF method are presented in Fig. 6. The AC power flow calculation, by (28), is implemented only twice in the F-SAPF method (once before the iteration and once after the iteration) and is not required in the iterative process, which greatly shortens the calculation time. Although y updated in the F-SAPF and SAPF methods will be a little different in each iteration, the accuracy of the F-SAPF method is high because the value of Δu1 in each iteration is small and the mismatch equations are the same as that of the SAPF method. Also, when the iteration converges, an AC power flow calculation will be implemented for correcting y as shown in Fig. 6. The number of iterations of the F-SAPF method might be different from the SAPF method, and this specific difference will be presented in case studies.

A. Case A: 3-Bus Small Test System
In this case, a small MTDC system with 3 DC buses (the total number of DC buses: n = 3) is used to illuminate the detailed calculation process of the SAPF and F-SAPF methods, in which the AC side is based on the CIGRÉ NORDIC 32-A grid [34]. This 32-bus AC grid (the total number of original AC buses: N = 32) can be considered as a scaled down version of the Swedish grid with three areas (North, Central, South), and the AC buses of this grid are renumbered as shown in Fig. 7. Three buses (No. 2,16,18) from different areas are connected to a 3-bus ±320kV MTDC system with DC slack control, and each DC bus is a VSC bus (the total number of VSC buses: m = 3). The control modes and corresponding references are presented in Table II. The resistance of DC cables is 0.012 Ω/km, and the parameters of VSC stations are shown in Table III [14] where B f,i denotes the susceptance of the filter.  Based on the extended AC model presented in Fig. 2, the node f 1 , f 2 , f 3 of VSC-1, 2, 3 are considered as the additional AC buses and numbered as the No. 33, 34, 35 bus of the extended AC grid respectively. Thus, the original 32-bus AC grid is extended to a 35-bus grid where the three additional AC buses are marked as green in Fig. 7. Accordingly, the sensitivity matrix S ac of the extended AC grid can be built by (26). The equations of each bus for building the calculation model represented by (30) are summarized in Table IV. Some independent variables in (30) are pre-set based on the power control of VSCs, which can be summarized as 1) For the DC-1 bus (i = 1), as the slack bus, U dc,1 = U * dc,1 by (7). Also, with the P-Q mode of the reactive power control, u 2,33 = u * 2,33 by (6), where u 2,33 denotes Q 33 (Q f,1 , as the AC bus 33 is renamed from the node f 1 of VSC-1) of the AC grid according to Table I. 2) For the DC-2 and DC-3 bus (i = 2 and 3), as the non-slack bus, P f,2 = P * f,2 and P f, 3 = P * f,3 by (8). Similarly, with the P-V mode, u 2,34 = u * 2,34 and u 2,35 = u * 2,35 by (6), where u 2,34 and u 2,35 denotes V 34 and V 35 (V f,2 and V f,3 ) respectively. So, considering the pre-set variables, (30) can be built and then simplified as where  Both S MT DC and S y can be extracted from S ac which can be calculated by inputting u and y into (29) in each iteration. The results are obtained after 5 NR iterations of both methods. Table V shows numerical results of the key independent variables in (33) with iterations. And the final result of y and y 2 of the F-SAPF method is obtained by inputting converged (28), as the steps shown in Fig. 6. From Table V, the results of the F-SAPF method are the same as that of the SAPF method except for the difference of 0.004% on P f,1 .

B. Case B: 6-Bus MTDC System
In Case B, the SAPF and F-SAPF methods are used to calculate the power flow of a 6-bus MTDC test system (n = 6). The results obtained by a typical unified method [19] are taken as the benchmark, and another sequential method proposed in [13] is adopted as the comparison. This MTDC system contains four VSC buses (m = 4), two pure DC buses (the DC buses without any VSCs connected), a PV power station, and a DC load of electric vehicle (EV) charging station. As shown in Fig. 8, the four VSC buses (DC-1, 2, 3, 4) are connected to the No. 8, 10, 30, 81 bus of an IEEE 118-bus system (N = 118) [35]. According to the extend AC model presented in Section II, the node f i in each VSC bus is considered as an additional AC bus, and these four additional AC buses are named as additional AC-119, 120, 121, 122 bus respectively, as shown in Fig. 8. Table VI shows the control references of VSCs. Also, two types VSCs, the half-bridge and full-bridge modular multilevel converter (MMC), are used in this case, of which the parameters and loss coefficients [13] are given in Table VII.
The calculation results of the four methods are presented in Table VIII. With DC slack control (DC-1 as the slack bus), the DC bus voltages obtained by the SAPF and F-SAPF methods are the same as the benchmark obtained by the method used in [19]. Instead, the sequential method in [13] has computational errors, which is because the converter current of each VSC is approximated. With droop control, the SAPF method and the method used in [19] also obtain the same results of DC bus voltages, while the results obtained by the F-SAPF method are very close to that of the SAPF method and the method used in [19] (only with an acceptable error of 0.001% on U dc,2 ). Meanwhile, the method used in [13] shows a relatively low accuracy. In Table VIII, the CPU time of the F-SAPF method is 63%, 70%, 72% less than that of the other three methods with the DC slack control, respectively, and 67%, 74%, 56% less with the droop control. This is because the calculation of the AC power flow in each iteration is replaced by the faster sensitivity analysis in the F-SAPF method. In addition, the numbers of iterations in the SAPF method, F-SAPF method and the method used in [19] are similar, which is because these three methods all iterate in a second-order speed based on the NR method with the same initial values. Although these three methods are with a similar number of iterations, the CPU time of each iteration of the method used in [19] is longer due to a larger system of nonlinear equations that needs to be solved [10]. Besides, the calculation speed of the method in [13] is faster than that of the SAPF method with droop control, but its calculation error is also larger, and more iterations are required. To summarize, in Case B, the SAPF and F-SAPF methods proposed in this paper are reliable in term of accuracy, and the F-SAPF method has a significantly higher computational efficiency than other methods.

C. Case C: MTDC System With 2 AC Grids
In Case C, a New England 39-bus AC grid [36] is connected to the IEEE 118-bus system by extending the existed 6-bus MTDC system in Case B, where three scenarios with different topologies of the MTDC grid are presented in Fig. 9. In scenario (a), a VSC bus (DC-7) is added to build a 7-DC-bus system (n = 7, m = 5) with 2 AC grids (N = 118+39). In scenario (b), two VSC buses (DC-8, 9) and a pure DC bus (DC-10) with another PV power station are added to build another MTDC gird which is connected to the New England 39-bus AC grid. Thus, this power network contains two AC grids (N = 118+39) and two separated DC grids (n = 7+3, m = 5+2). In scenario (c), another two lines (yellow lines shown in Fig. 9) are added to connect the two separated DC grids into one DC grid (N = 118+39, n = 10, m = 7). In Case C, the newly added VSC-5, 6 and 7 are connected to bus No. 39, 14, 23 of the New England 39-bus AC grid respectively, and all VSCs are with the droop control. Detailed control modes and reference values of VSC-5, 6, 7 are  Table IX, and the parameters of them are consistent with those of the half-bridge MMC shown in Table VII. Table X shows the power flow calculation results of the test system in these three scenarios, where U dc , V f and P loss dc denote the average voltage of DC buses, average voltage of additional AC buses, and power loss of the MTDC grid (including the power loss of VSCs) respectively. Compared to scenario (a), the CPU time of the SAPF, F-SAPF and the method used in [19] in scenario (b)/(c) increases by 9%/12%, 8%/10% and 16%/17% respectively. In terms of calculation accuracy, for all three scenarios, the calculation results of the SAPF method and the method used in [19] are still consistent except for the 0.1% differences in P loss dc . Also, the difference between the results obtained by the F-SAPF method and the method used in [19] is less than 0.3%.
In order to study the influence of the initial values X (0) of the power flow calculation on the CPU time and convergence of given methods, the initial values of the power flow calculation are selected in two ways and tested based on scenario (b). The first way to select initial values is by where X * is the accurate result of the power flow obtained by the previous simulation, and c is the proportion coefficient. Fig. 10(a) shows the CPU time cost by the SAPF, F-SAPF and the method used in [19] with c = 1, 2, 3, · · · , 20, where the calculation of the method used in [19] does not converge with c = 1 and c = 2. Besides, the CPU time of the three methods is less when the initial values are closer to the accurate results X * . Also, the CPU time of the F-SAPF method is the least, as shown in Fig. 10(a). The initial values selected by another way are the random number in the interval of (0.5, 1.5). Fig. 10(b)  shows the CPU time with the random initial values. Similarly, the F-SAPF method costs the least time to calculate the power flow in the 10-bus DC system compared to other two methods. In Fig. 10, it can be concluded that the small initial values will increase the calculation time and may cause the calculation of the NR method to be non-convergent.

D. Case D: Offshore Wind Power System
In this case, a fraction of a developing offshore wind power grid called the North Seas Countries Offshore Grid [37] is introduced, in which offshore wind farms of three countries (the United Kingdom: U.K., The Netherlands: NL and Germany: DE) are connected to onshore grids by DC submarine power cables (with resistance 0.0195 Ω/km) and VSC stations. Fig. 11 shows the topology of this system (n = 9, m = 6). The proposed methods are implemented to study the impact of the varying wind farm outputs and presumptive contingencies on the power flow of the offshore wind power system. Three independent IEEE 300-bus systems [35] are used to represent the U.K., NL  Table XI. Fig. 12 shows the variations of power flow indicators with the changing output of the U.K. offshore wind farm, while the other two wind farms are with the rated power output (1000MW). In Fig. 12(a), the average voltage of DC buses drops when the output of the U.K. wind farm decreases. Compared to the droop control, the DC slack control maintains the average DC voltage at a higher level. Also, a higher average voltage of DC buses can be observed under droop control with a higher droop coefficient. The opposite trend can be observed in Fig. 12(d) that the maximum current of DC branches rises as the decreasing output of the U.K. wind farm. Besides, the DC slack control leads to more power losses in the offshore wind power system than the droop control as shown in Fig. 12(b). In Fig. 12(c), the inverter power of DE VSC with DC slack control shows a more is the maximum current among the DC cables. As shown in Table XII, the CPU time of the F-SAPF method is over 70% less than that of the other three methods, showing the high efficiency of the F-SAPF method. The power flow results obtained by the SAPF, F-SAPF and the method used in [19] are almost the same while a relatively large difference can be seen between the results obtained by the method in [13] with other methods, which is similar to the observations in Case B and C. In addition, the CPU time of all methods for the calculation with droop control is larger than that of the DC slack control. This is because the number of equations for the droop control is more than that of the DC slack control.
According to the power flow results shown in Table XII, when the U.K. wind farm is forced to an outage, the total power loss of the DC grid drops by 21% with droop control but increases by 77% with DC slack control. Also, with DC slack control, the inverter mode of the DE VSC needs to be adjusted to the rectifier mode after the outage of the U.K. wind farm, which is difficult to implement in a short time. On the contrary, with the droop control, three onshore VSCs in this MTDC system are maintained in stable operation status. When the NL or DE onshore VSC is offline, the offshore wind power system can keep stability only by the droop control.

E. Case E: the CIGRÉ B4 DC Grid With Large AC Systems
To evaluate the performance of the proposed methods for large networks, the CIGRÉ B4 DC grid [38] is adopted in this case, which contains 15 DC buses (n = 15), including 11 VSC buses (m = 11) and 4 pure DC buses, as shown in Fig. 13. This DC grid shows the ability of the MTDC system to integrate renewable generation from offshore grids into different synchronous networks. In this case, the DC-DC converters in the MTDC system are considered as DC sources with constant power flow to the connected DC buses. Also, the five AC offshore VSCs (Cm-C1, Cb-C2, Cb-D1, Cm-E1, Cm-F1) are set to be with DC slack control, and corresponding DC buses are non-slack DC buses. The other VSCs are with the droop U-P control where K i = 40 and U * dc,i = 1 p.u. The VSC Cm-A1 and Cm-B2 are set to be working in PQ mode with Q * dc,i = 0, and the other VSCs are working in PV mode where V f,i = 1 p.u. The setpoint of P * dc,i and parameters of the transmission line can be found in [38]. The VSCs in this case are all the half-bridge MMC. Note that, there is a situation where two VSCs are connected to the same AC bus, where two additional AC buses will be added and connected to this AC bus. In this case, the CIGRÉ NORDIC 32-A system [34], a 3012-bus AC test grid [35], and a scaled down version of the European system model [39] containing 2869 buses are used for representing the real AC grids. Different AC grids are connected to the bus Ba-A0 and Ba-B0 shown in Fig. 13, in which the connection buses of these large AC grids are shown in Table XIII. The power injections from AC bus Ba-A1, Ba-B1, Ba-B2, Ba-B3 to the DC grid are denoted as P AC Ba−A1 , P AC Ba−B1 , P AC Ba−B2 , P AC Ba−B3 respectively, which are obtained by three methods as shown in Table XIII. Based on the results shown in Table XIII, it can be observed that the impact of the scale of the AC grid on the accuracy of the proposed SAPF and F-SAPF methods is limited, as the difference between the results by the SAPF method and the method used in [19] is less than 0.15%. Besides, with the significantly larger AC systems, three methods cost a longer CPU time to calculate the power flow than previous cases. In the calculation with N = (32+2869), the CPU time of the F-SAPF method is less by 69% and 74% than that of the SAPF method and method used in [19], respectively, which shows the high efficiency of the F-SAPF to calculate the power flow for the MTDC system with large AC grids. In the calculation with N = (3012+2869), compared with the results with N = (32+2869), the CPU time of the three methods is increased by 91%, 78%, 208% respectively, which shows the calculational burden caused by introducing the very large AC system into the calculation. Because the proposed SAPF and F-SAPF methods can separately generate the sensitivity matrix for each AC grid instead of combing the two AC grids into a large system as the method used in [19], the increase in the CPU time of these two methods is significantly less than that of the method used in [19].

V. CONCLUSION
This paper proposes a novel power flow calculation method for the MTDC system based on the sensitivity analysis and the extended AC grid. In the proposed power flow calculation model, the mismatch equations of the DC side have been built according to the control mode of the VSCs, while the AC grid has been extended without any change to the AC power flow calculation. The sensitivity analysis has been adopted in the proposed SAPF method which has higher calculation efficiency than the unified method and causes no convergence problem compared to the sequential method. To further improve the calculation efficiency, the F-SAPF method has been proposed, where the time-consuming AC power flow calculation in each iteration of the SAPF method has been replaced by the faster sensitivity-based calculation. A small test system has shown in Case A for presented the detailed calculation process with numerical results. In Case B, the accuracy and higher computational speed of the proposed SAPF and F-SAPF methods have been shown by comparing with the other two existing methods. Case C has presented the performance of the proposed methods with different topologies of the MTDC system. With high convergence, the wide applicability and high robustness of the proposed method have also been illuminated in Case D, where the proposed methods have been applied to determine the power flow variations of an offshore wind power system with the changing output of offshore wind farms. It also shows the advantage of droop control over DC slack control in maintaining power flow stability in given contingencies. The accuracy of the proposed methods to calculate power flow for large power networks has been presented in Case E, where the comparison of the results obtained by different methods has also shown the advantage of the F-SAPF method in computational efficiency with the scenarios of large power networks.