A Novel AC/DC Power Flow: HVDC-LCC/VSC Inclusion Into the PFPD Bus Admittance Matrix

In this paper, the matrix algorithm PFPD is generalized in order to compute the power flow solution of real and large AC/DC transmission networks. In particular, it is demonstrated that the HVDC-VSC/LCC links can be seen from the AC power systems as PV/PQ constraints, which englobe both the AC and DC characteristics of the HVDC links. The proposed analytical formulation to assess the PV/PQ constraints is valid for any other numerical methods (e.g., Newton-Raphson and derived, Gauss-Seidel, etc.). Furthermore, an iterative procedure for estimating the reactive power absorption of HVDC-LCC links from the power system is proposed. In order to validate the algorithm, solution comparisons with the commercial software DIgSILENT PowerFactoy are presented. This validation procedure shows that the algorithm can analyse large and real HVAC/HVDC networks (e.g., the Italian transmission one with its five HVDC links). Therefore, the conciseness, accuracy and performances of PFPD for studying real and large AC/DC power systems is confirmed.

In 2022, Roberto Benato published a novel AC matrix power flow algorithm named as PFPD [1]. This algorithm is based on an extensive use of the bus admittance matrix including generators, loads, and also the slack generator: the theoretical possibility of modelling the slack generator as a quasi-ideal current one gives a very easily-implementable, efficient and fast algorithm for AC networks. Notwithstanding, the paper [1] does not deal with the AC/DC power flow problem i.e., the compresence of both HVDC-VSC and HVDC-LCC point-to-point links inside an AC power system. This problem arises since the compresence of HVDC-LCC/VSC links in HV and EHV transmission networks is an increasingly widespread reality. In fact, in some contexts, several reasons make the HVDC technologies preferable to HVAC ones, such as: the greater interconnection capacity [2], [3]; the more economical, reliable, and sustainable operation [4]- [6]; and the possibility to integrate large scale of Renewable Energy Sources (RES) in power systems [7]. Therefore, it is firstly necessary to have powerful tools to compute the power flow solution of such AC/DC networks for steady-state operation evaluations, contingency analyses, and planning. Moreover, the power flow evaluation is the initial necessary step to make dynamic evaluations [8]. By considering all the implementation, computational, and research advantages of [1], a general formulation of PFPD can be a valid alternative to the Newton-Raphson formulations to compute AC/DC power flow.
This paper proposes a matrix power flow formulation for solving the AC/DC power flow problem efficiently by means of PFPD.

B. LITERATURE REVIEW
After the first HVDC industrial application in 1954, investigations on power flow formulations for networks containing HVDC links started in the 1960s [9]- [12]. All the methods in [9]- [12] use the Newton-Raphson approach to solve the non-linearity of the power flow problem. In such contributions, the DC and AC systems are treated separately and not solved simultaneously: this choice makes the approach inefficient. In [13], [14], the AC and DC equations are formulated simultaneously, therefore a unique Jacobian matrix is adopted.
Although the first power flow methods were matrix ones [15], [16], such approaches have been never adopted to solve the AC/DC power flow problem. The matrix approaches, in fact, are long abandoned. In [1], however, the strengths of matrix approaches for power flow purposes are highlighted.
Until the 1990s, the AC/DC power flow formulations in power systems with HVDC links consider only the HVDC-LCC technologies (which is the first-appeared HVDC technology). After the first experimental applications in the late 1990s of the HVDC-VSC model [17], [18] proposes a research modelling of such devices for power flow purposes [19]- [22]. In particular, the converters can control the active and reactive power independently. This greater controllability makes the power flow formulation more flexible.

C. CONTRIBUTIONS
In this paper, it is demonstrated that PFPD can be generalised for assessing the steady-state regime of modern AC/DC transmission networks. Therefore, PFPD is intended as the first iterative matrix method to solve the AC/DC power flow problem for transmission networks (different formulation from the well-known numerical methods e.g. Newton-Raphson, decoupled approaches, Gauss-Seidel, etc.). In fact, once the typical HVDC controlled variables are fixed, the classical PV/PQ constraint formulation is valid for modelling pointto-point HVDC links.
In particular, the HVDC-VSC links are modelled as PV constraints (they can also be modelled as PQ constraints, since the active and reactive power can be controlled independently), whereas HVDC-LCC links are modelled as PQ ones. Since, in LCC devices, the absorption of reactive power cannot be controlled independently from the active one, an iterative procedure (IRPE) estimating the reactive power absorption of the LCC converter for power flow purposes is presented. Therefore, the AC/DC power flow formulation can be reduced to an AC one, and the HVDC information is englobed in the ''all-inclusive'' admittance matrix of PFPD [1]. As a consequence, the DC variables quantities disappear from the iterative formulation, but their impact into the AC system is correctly assessed (e.g., global power loss computation). Notwithstanding, the actual DC values of voltages and currents can be computed normally. It is worth noting that such way of modelling an AC/DC network as an AC one is general, so it can be exploited for all the classical numerical methods (e.g. Newton-Raphson, decoupled approaches, Gauss-Seidel, etc.).
Eventually, some solution comparisons with the commercial software DGS are presented in order to validate the proposed method. In particular, the algorithm is tested by the Italian network, which has got five HVDC links (four with foreign nations), in order to ensure the applicability of PFPD to a real and large transmission network.

II. PFPD HVDC MODELLING BY MEANS OF THE DC ELIMINATION TECHNIQUE
In this study, a suitable modeling of both the HVDC-VSC (see Sect. III) and HVDC-LCC (see Sect. IV) technologies for PFPD is presented.
The HVDC links are modelled by means of PV or PQ constraints to be embedded inside the ''all-inclusive'' matrix Y SGL [1]. Obviously, such PV and PQ constraints do not represent the synchronous generators and the loads respectively, as in the classical AC power flow problem. However, they suitably assess the impact of the point-to-point HVDC technologies in the power system from a steady-state point of view. Moreover, the station converters and the DC link are not directly represented, but their presence is taken into account when the PV/PQ formulation (this fact is indicated as ''DC elimination technique'').
Although the converters introduce voltage/current harmonics in the power system, a ''fundamental frequency model'' approach is adopted, since the power flow problem is related to the power system fundamental frequency. In fact, the filters on both AC and DC sides cancel the harmonic power contribution.
However, the impact of the filters for harmonic compensations is taken into accont at the fundamental frequency only. In particular, the power losses due to current absorption in such devices is computable.
Eventually, HVDC converters are usually connected with OLTC transformers, which allows controlling the AC secondary voltage and so the DC link voltage (together with the angle α). Such elements are implemented and their modelling is described in [1].

III. HVDC-VSC LINK MODELLING A. THE TRANSMITTER STATION MODELLING (PV-CONSTRAINT MODELLING)
The VSC transmitter converter can control both the AC terminal voltage and the transmissible active power. Fig. 1 shows a monopolar HVDC-VSC link supplied by two-winding transformers.
The power flow direction conventionally goes from the transmitter converter (S) to the receiver one (R). The transmitter converter S is controlled to absorb, from the AC network, a scheduled active power to be transmitted i.e., P AC,S = P sched .
In fact, in VSC converters, the voltage phasor position can be controlled with respect to the network one. Moreover, the voltage value at its AC terminal can be scheduled: Thus, the AC terminal of the VSC transmitter converters can be set as PV constraint, where P = -P sched and V = V sched . The negative sign of P is due to the fact that the transmitter converter absorbs (and does not inject) active power from the network (active sign convention for PV constraints).

B. THE RECEIVER STATION MODELLING (PV-CONSTRAINT MODELLING)
The receiver converter typically controls both the AC and the DC terminal voltages: The voltage V DC,R , in fact, must be controlled and kept constant for two reasons: ensuring both the proper operation (V DC smooth voltage) of the VSC converters and the active power balancing from the transmitter to the receiver converter.
The receiver converter R converts the DC-bus power P DC,R into the AC active power P AC,R , which can be computed by knowing the active power losses P loss occurring in the VSC converters and in the DC line: The expression of P loss is the following: P loss = P loss_conv,S + P loss_conv,R + P loss_line (2) where P loss_conv,S and P loss_conv,R represent the losses due to converter S and R respectively, whereas P loss_line is due to the Joule losses in the DC line. The value of the Joule losses in the DC line can be computed by considering the current I DC : (3) VOLUME 10, 2022  By considering the ''fundamental frequency model'' approach, the losses of each converter are given by the sum of three terms [23]: P loss_conv = P loss_noload + P loss_switch + P loss_cond (4) where P loss_noload represents the no load loss contribution, P loss_switch represents the switching loss contribution and P loss_cond represents the resistive loss contribution. As shown in Table 1, for each term of the right-hand side of (4), an input loss factor is defined: P no,load for the no load losses, k for the switching losses and k for the resistive losses. Thus the addends of (4) can be computed by knowing such loss factors characteristic of the converter [24].
For the no-load losses P loss_noload , the conductance coefficient G no_load must be computed by knowing the its input loss factor P no,load [MW] and the DC rated voltage [kV] as follows [24]: Therefore, the no-load losses, P loss_noload , can be computed: The switching losses, P loss_switch , depend on the current I DC circulating in the DC line as follows: where the coefficient V drop is: and k is the loss factor characterising the switching losses per unit current. The resistive losses P loss_cond can be computed according to the following relation: where k is the resistive loss input factor. In order to compute I DC , the following system is considered: which leads to: Therefore, once the controlled variables are fixed, I DC is a term computable a priori (i.e., before the power flow computation) by means of (11). By considering the set of equations (2)÷(9), the value I DC allows determining P loss . Therefore the PV active power constraint of the receiver converter AC terminal can be computed by means of (1).
It is worth underlying that modelling the transmitter and receiver converters by means of conventional PV constraints allows finding the AC-side converter phase angles as the solution of the power flow problem.
By extending the above-mentioned procedure, the following HVDC-VSC typical configurations [2] can be immediately modelled: • Type 1: Monopolar HVDC-VSC supplied by two winding transformers (see Fig. 1), • Type 2: Bipolar HVDC-VSC supplied by two winding transformers (see Fig. 2a)), • Type 3: Monopolar HVDC-VSC supplied by three winding transformers (see Fig. 2b)), • Type 4: Bipolar HVDC-VSC supplied by three winding transformers (see Fig. 2c)). As it can be seen in Fig. 2, every AC converter terminal is considered as a PV constraint. Therefore, multibridge converter configurations can be studied, by using one PV constraint for each converter. However, in industrial reality, for high power transmission the double-bridge solution is the most used in order not to complicate the DC control [25].
The schedule active power of the sending converter station S can be equally considered shared among the different n S PV constraints. Therefore, the power of the sending i-th PV constraint which models the i-th converter is: each of which is characterised by the corresponding scheduled AC voltage: With regard to the constraints of the receiving converter, (2) can be generalised in the following: (12) Therefore, by subdividing the global losses of (12) for all the receiving PV constraints, the active power of the i-th PV receiving constraint is: each of which is characterised by the AC voltage controlled by the corresponding converter: C. HVDC-VSC LINKS (PQ CONSTRAINTS) Since active and reactive power can be independently controlled, the VSC converters can be also modelled as PQ nodes, where P and Q are the scheduled power values. In this way of modelling, the AC-bus converter voltage magnitude cannot be controlled and depends on the power flow solution. Therefore, this control mode can be used only if it is assured that the active and reactive power assume values not causing high voltage drops on the VSC converter AC bus node (strong AC network as seen from the AC-bus converter). Hence, this control mode has got effects which are different from the VSC converter modelling meant as PV constraints, where the voltage magnitude is set (but not the reactive power).
The active power loss computation is the same of the one described in sub-section B. The transmitter converter controls the transmittable DC power. In fact, under the hypothesis of constant current I DC , the value of the DC voltage can be controlled through the delay firing angle α (0 • < α < 90 • ) as it follows [23]: (15) where the coefficients k , V drop assume the same meaning already explained for VSC converters [23]. Eq. (15) is valid for both inverter and rectifier operation modes and take into account the loss contribution. Therefore, it is possible to transmitt the scheduled DC power from the trasmitter convert station, i.e., In order to compute I DC , the following system is considered: which leads to: Thus, it is possible to find I DC .
Similarly to the VSC case, once the controlled variables are fixed, I DC is a term computable a priori (i.e., before the power flow computation).
The transmitter converter absorbs both active and reactive power from the network, so it is possible to define a PQ constraint at its AC terminal. The AC active power absorption of the transmitter converter depends on the active power losses in the converter S: According to the loss computing approach of [23], for the thyristor converters, the relations (5)÷(9) are still valid. Therefore, the active power absorbed of (19) as seen from the AC LCC rectifier terminal can be computed.
For the reactive power computation, see subsection C) of the present section.

B. THE RECEIVER STATION MODELLING
The receiver converter R converts the DC bus power P DC,R into the AC active power P AC,R : the converter operates in the inverter mode (90 • < α < 180 • and α < 180 • -γ to avoid commmutation failure [26]). In the inverter operation, the active power transmitted in the AC network depends on the DC line and on the converter R: where the formulations (3)÷(9) are still valid. It is worth noting the negative sign of P AC,R , since the positive power is injected into the network. As it is explained in sub-section C), the inverter reactive power computation is the same of the rectifier one. By extending the above-mentioned procedure, four different HVDC-LCC configurations are modelled (see Fig. 4): • Type 1: HVDC-LCC monopolar, supplied by two winding transformers (see Fig. 3) • Type 2: HVDC-LCC bipolar, supplied by two winding transformers (see Fig. 4

b)
• Type 4: HVDC-LCC bipolar, supplied by three winding transformers (see Fig. 4 c) As it can be seen in Fig. 4, every AC converter terminal is considered as a PQ constraint. As for VSC, multibridge converter configurations can be studied, by using one PQ constraint for each converter. The active power of each sending converter station S can be computed by equally dividing P DC,sched into n S converters: P DC,S(i) = P DC,sched /n S , and then by applying (19) for each i-th converter. In order to model the active power injected by the receiving converters, the sum of the losses in all the lines and receiving converters must be taken into account: Therefore, by subdividing the global losses of (21) for all the receiving PQ constraints, the active power of the i-th PV receiving constraint is: The negative sign of (22) is due to the fact that the PQ constraint is injecting active power. The reactive power computation for both the rectifier/inverter converters is estimated by means of the IRPE method, described in the following sub-section.

C. THE IRPE METHOD: THE LCC REACTIVE POWER ABSORPTION COMPUTATION
In this subsection, a method to estimate the LCC reactive power absorptions for power flow purposes is presented.
Firstly, it is assumed that the commutation reactance and the converter losses are negligible. The reactive power absorbed by an LCC converter depends on its control angle ϑ. This angle corresponds to the delay firing angle α for a rectifier and to the leading firing angle β for an inverter (hence, 0 • < ϑ < 90 • ).
By considering an LCC converter, the following relation is valid: where Q is the reactive power set absorbed by the converter and P is alternativelly the active power absorbed from the AC system (for rectifiers) or injected into the AC system (for inverters).
Differently from the VSC case, the LCC converter AC side voltage cannot be controlled, since it depends on the power flow solution, and the value of ϑ is set to control the active power. Therefore, the value of the reactive power Q computable by means of (20), cannot be controlled a priori, but it is a function of the power flow problem.
In order to compute the reactive power Q, it is possible to proceed iteratively.
By considering an LCC converter (a rectifier or an inverter), its reactive power absorption is and by substituting the relation between the DC current and the AC one [26] I AC = √ 6 π I DC (26) in (25), the expression (27) for the reactive power computation can be derived: which basically links the reactive power Q conv with the unknown AC voltage module V AC . The values of I DC and P sched are constant (I DC is computable a priori by solving (18) and P sched is the controlled quantity). It is worth noting that Q conv is always positive for both rectifier and inverter.  Formulation (27) is suitable to be computed iterativelly as a cycle external to the PFPD.
The initial guess for the AC voltage module of the converter is set to 1 p.u. This value is near to the power flow solution since, for transmission purposes, LCC converters ought to be connected to strong AC nodes.
By substituting this initial guess in (27), an estimation of the reactive power Q conv is obtained. Therefore, the LCC converter can be modelled as PQ nodes, where Q = Q conv . Eventually the PFPD computation update the voltage module of the converter. Therefore (27) can be applied again, and the procedure is iterativelly repeated (see Fig. 5) until the mismatch between the voltage module of two consecutive iterations is less than a predifinied tolerance tol: In this study, a typical value for this voltage tolerance tol is equal to 100 V. The present procedure can be summarized and visualized by considering the flow chart of Fig. 5. Such iterative procedure can be written in matrix formulation (vectorization), thus it is possible to consider all the HVDC link reactive power updates simultenously.
This reactive power estimation is named as IRPE (Iterative Reactive Power Estimation). The method presented is valid for LCC converters operating both in inverter and rectifier mode (both the operation modes need to absorb positive reactive power from the AC system) [27].

V. INSERTION OF THE HVDC-VSC/LCC POWER FLOW CONSTRAINTS INSIDE Y SGL
In PFPD the PV constraints modelling the HVDC-VSC links are modelled as passive admittances [1] as in the following: where i ∈ a ÷ g. The admittance y i is positioned in the PV submatrix, P i and V i are fixed in accordance with the procedure of Sect. III A and III B, and Q i is the unknown reactive power absorbed by the HVDC-VSC link.
Once the PV constraint admittances are computed, they can be embedded in the i-th diagonal position of the square Y SGL admittance matrix (see Fig. 6a). It is worth noting that in such admittances the reactive power estimation is the same as that for estimating reactive power of the conventional AC generators in PFPD [1].
Similarly, the PQ constraints modeling the HVDC-LCC links are modelled as passive admittances by means of (30): where i ∈ h ÷ m. y i is positioned in the PQ submatrix P i and Q i are fixed and computed according to the procedure of Sect. IV A/IV B and Sect. IV C. It is worth noting the opposite signs of the expressions (29), (30). Q i is the unknown reactive power absorbed by the HVDC link.  Once the PQ constraint admittances are computed, they are embedded in the i-th diagonal position of the square Y SGL admittance matrix (see Fig. 6b).
Therefore, by considering a network with n buses, the dimension of the AC/DC problem is always described by an (n × n) Y SGL matrix, which contains the HVDC information also by means of PV/PQ constraints.
This fact is valid independently from the number of HVDC links. Therefore, the insertion of further HVDC in an n-bus network does not increase the dimension of the problem.
It is of note that the elimination of the DC part can create different separated subsystems if the two connected systems by the HVDC link have no other connections between them [25]. However, this fact does not introduce any computational critical issue.
Hence, once the matrix Y SGL containing HVDC information is built, the iterative procedure is the same of the AC situation.
Moreover, the modelling of the HVDC links by means of shunt admittances (and not branch admittances) makes the matrix further well-conditioned [1], [28], since such admittances are in the matrix diagonal position.

VI. THE DC QUANTITIES
By considering the typical converter control schemes adopted in the HVDC links and described in this paper, the DC quantities can be computed directly, without the need of using any iterative procedure.
For both the HVDC-LCC and VSC configurations, in fact, the transmitted active power and the DC voltage at the sending side are constant. As it is explained in Sect. II, the DC voltage of the receiver side converter is constant, whereas the current I DC can be computed a priori by means of (11) and (18).
Therefore, the value of the sending DC voltage can be computed by considering the voltage drop due to the DC current circulation. In general, by considering the series of the converters at both side: where r i and i are the resistances and the lengths of the DC connection sections.
Eventually, the value of α for rectifiers and β for inverters can be computed by the knowledge of the AC active/reactive power absorbed by the converters after the power flow calculation.

VII. AC/DC POWER FLOW SIMULATIONS BY MEANS OF PFPD
PFPD is implemented in Matlab environment. In order to test its effectiveness in including HVDC links, solution comparisons with the commercial software DGS are shown.
All the HVDC connections of this paper (see Fig. 1, 2, 3, and 4) are considered one at a time.
After the computations, the power flow solutions in PFPD and DGS are compared in terms of voltage magnitude and angle. For the voltage magnitude, the following vector of the relative mismatch between PFPD and DGS solutions is computed: where V PFPD and V DGS are the voltage magnitude vectors, and V is the voltage magnitude mismatch vector. For the voltage angles the following mismatch vector between the PFPD and DGS solutions is computed: where δ PFPD and δ DGS are the angle magnitude vectors, and δ is the angle magnitude mismatch vector.  Table 2 reports the order of magnitude of the maximum solution mismatches between PFPD and DGS computed by means of (32) and (33). Such values confirm the very good agreement between the two computational methods.

B. VALIDATION OF PFPD: HVDC-LCC LINKS
The four HVDC-LCC configurations are inserted between buses 15 and 16 of the AC 18-bus test network (see App. I). Therefore, the following four test cases are analysed: 1. Case E: Network containing monopolar HVDC-LCC link supplied by two-winding transformers.   The tolerance of PFPD is set to 10 −8 p.u. Since in such test networks, there are HVDC-LCC links, a tolerance for the IRPE method equal to 100 V is fixed. Table 3 reports the order of magnitude of the maximum solution mismatches between PFPD and DGS computed by means of (32) and (33). Also in this case, the values confirm the very good agreement between the two computational methods.

C. VALIDATION OF PFPD: COOPERATIOON OF HVDC-VSC AND HVDC-LCC LINKS
In this section, the possibility of computing the power flow solution by considering more HVDC links in the same network is assessed. A 60-bus fictitious network is considered and the three following cases are studied: 1. Case I: Network containing a monopolar HVDC-VSC link supplied by two-winding transformers and a bipolar HVDC-VSC link supplied by three-winding transformers. 2. Case L: Network containing a monopolar HVDC-LCC link supplied by two-winding transformers and a bipolar HVDC-LCC link supplied by three-winding transformers. 3. Case M: Network containing a bipolar HVDC-VSC link supplied by three-winding transformers and a bipolar HVDC-LCC link supplied by three-winding transformers.   The tolerance of PFPD is set to 10 −8 p.u.; a tolerance of the IRPE method equal to 100 V is chosen. Table 4 reports the maximum solution mismatches between PFPD and DGS computed by means of (32) and (33).
Such values confirm the very good agreement between the two computational methods and the capability of PFPD to treat different HVDC links inside the same network.

VIII. A REAL-WORLD APPLICATION OF PFPD: THE ITALIAN TRANSMISSION NETWORK
Due to its geographical collocation inside the synchronous ENTSO-E transmission system, the Italian network is connected with six nations by means of twenty-five interties [29]. VOLUME 10, 2022     Four of these links are HVDC ones, whose characteristics are summarized in App. II. Another HVDC link connects the peninsular part of Italy with the Sardinia region (see in App. II the SA.PE.I link).
In order to test the industrial applicability of PFPD, a final validation is shown by considering the real-world Italian transmission network data.
Two different configurations of the Italian network are considered: a high load scenario and a low load scenario one. Their data are implemented both in Matlab environment and in DGS. The tolerance of the algorithm is set to 10 −8 p.u. in both environments. In PFPD a tolerance for the IRPE method equal to 100 V is fixed. In Table 5 solution displacements in terms of voltage angles and magnitudes are presented. It is of note the good agreement between the two solutions: by considering both the scenarios, the maximum angle displacement computed for the voltage angle is 0.53 • (found in the low load scenario) and 0.57 • (found in high load scenario), whereas the order of magnitude of the maximum voltage displacement computed in both the scenarios is about 0.1%.

IX. HVDC-LCC RECTIVE POWER ABSORPTION: THE MON.ITA LINK CASE STUDY
In this section, an example of the applicability of PFPD to assess the impact of HVDC-LCC links in real power systems is shown. An experimental validation of PFPD by considering real steady-state measurements is assessed. It is known that, in HVDC-LCC links, the higher the transmitted active power, the higher the converter reactive power absorption [30] as inferable from (34):    In order to maintain the reactive power absorption from the AC network inside a specified range, RPC logic is used [31], [32]. RPC logic basically insert/disconnect compensation devices in the PCC node to limit the reactive power absorption from the AC network. Fig. 7 shows the experimental measurements of an RPC test carried out during the commissioning of the MON.ITA HVDC link in 2019 [32]. These measurements are related to the rectifier side. To perform this test, the transmittable power is changed by step of 50 MW, and the RPC compensators are inserted/disinserted according to the active power transmitted. There are four reactive power compensators [32], whose combinations define different compensation stages depending on the transmitted active power.
It can be observed that the reactive power injected by the AC network (red line) is inside a ±50 Mvar range, so complying with the RPC specifications. The blue line represents the VOLUME 10, 2022 reactive power absorbed by the rectifier side: its requirement is guaranteed by both the AC networks and the AC compensation devices. The rated power of the adopted compensation stages is indicated by the black line.
The test configuration during the RPC test is modelled in PFPD, by considering the real Italian network data. The active power transmittable by the MON.ITA link is changed by considering 50 MW steps and by inserting the corresponding compensation stages. Fig. 8 illustrates the correlation between active and reactive power absorbed by the rectifier side of the MON.ITA link. A comparison between the experimental active/reactive measurements [32] and PFPD is shown. This case study confirms that PFPD results are consistent with the experimental measurements.

X. OPEN QUESTIONS
The present power flow procedure considers the two-terminal HVDC modelling by means of PV/PQ constraints. No HVDC multi-terminal configuration is considered.
Notwithstanding, some existing HVDC links could be connected to other ones: in the future, there could be more and more DC portions of radial/meshed transmission networks [2]. These DC network portions must be located inside the AC/DC boundaries which connect them to the AC network, also having integrated DC sources and DC loads. The authors hope that, in such future networks, new DC power flow constraints due to the converters can be fixed to set the power flow equations solving the problem. By starting from the DC network solutions, AC/DC boundaries could be set as PV/PQ constraints in order to study the AC network, similarly to the approach explained in this paper. Therefore, authors think that the presented power flow procedure can be further generalised to model HVDC multi-terminal links by means of PV/PQ constraints.

XI. CONCLUSION
The paper presents how PFPD can be suitably generalized to include HVDC-LCC/VSC links in the power flow of real and large AC/DC power systems. This matrix algorithm keeps all its performance in terms of CPU-time and convergence efficiency. HVDC-LCC/VSC are simultaneously considered by means of PQ and PV constraints: these ones follow the typical controls with the setpoints used in the different HVDC typologies. For HVDC-LCC, the reactive power to be included into PQ constraint are pre-estimated by an iterative procedure named as IRPE (which could be implemented also in N-R methods). Even if the DC networks are eliminated, after the convergence, all the DC quantities can be computed: transmission line and converter power losses, LCC converter firing angles etc. Filter behaviours are also considered at power frequency with their power losses.
Extensive comparisons with DIgSILENT PowerFactory also demonstrate the high level of accuracy of PFPD power flow solutions. Eventually, PFPD is efficiently applied to the real Italian grid with four HVDC-LCC and one HVDC-VSC links: some comparisons with measurements carried out during the real-time operation of Italy-Montenegro HVDC-LCC shows the power and accuracy of PFPD. Such open-access algorithm also plays and will play a key role for the planning and operation of the future Italian HVDC systems i.e., the Adriatic and Tyrrhenian links.

APPENDIX I: THE 18-BUS TEST NETWORK
The AC 18-bus test network used in this study is represented in Fig. 9. The different HVDC link typologies are inserted between nodes 15 and 16.

APPENDIX II: THE ITALIAN HVDC LINKS
In the Italian transmission network there are five HVDC links. Table 6 reports the main characteristics of these connections. For the Italian transmission network, two different scenarios are studied: a low load scenario and a high load scenario. In PFPD, the AC side terminals of the HVDC converters are modelled by considering the foreign networks as PV nodes, i.e., nodes in which the active power and voltage is scheduled. In DGS, instead, such nodes are modelled as external grids, which are specific blocks modelling entire networks as seen from one terminal.
There are several ways to set these external grids, however, in order to perform comparisons at the same conditions, they are modelled as PV nodes.

APPENDIX III: NUMERICAL RESULTS
In this appendix, some numerical results of the simulations described in Sect. VII and VIII are reported. The quantities represented in all the tables are rounded to five decimal places.
As an example, in Table 7 the comparisons between PFPD and DGS for the Case E (see Sect. VII) in terms of the AC nodal voltages are reported. Table 8, 9, and 10 report the computed quantities associated with the HVDC link between bus 15 and 16 (see Fig. 9). These tables show the significant quantities related to the HVDC link, and are subdivided into the rectifier, inverter, and DC line parts. It is worth noting from Table 7 to Table 10 the very good agreement between PFPD and DGS of the AC/DC results. In particular, the DC quantities assume the same value in both the environments: as it is said in Sect. VI, their computation is analytical. In particular, by means of (11) and (18), the value of the DC currents I DC can be computed analytically, and so the computation of the DC voltages, and power losses. Therefore, the evaluation of the power losses introduced by the HVDC links can be estimated accurately by means of the approach described by (2)÷(9).
The above-mentioned consideration can be extended to the simulations of the Italian transmission network. By considering, as an example the low load scenario (see Sect. VIII for its description), Table 11, 12, and 13 report the computed quantities associated with the five HVDC links existing in the Italian network. Fig. 11 and 12 synthetize the results of Table 7, and Tables 11/12 respectively. Fig. 11 represents the AC voltage magnitudes and angles of all the AC nodes of the Case E network (see Fig. 9). Fig. 12 shows the AC voltage magnitudes and angles, and the DC voltages of the five Italian HVDC links.
The p.u. DC voltages are computed by considering the DC base voltages of each HVDC link (which correspond to the inverter-side controlled voltages).