Derivative-Free Power Flow Solution for Bipolar DC Networks with Multiple Constant Power Terminals

This paper analyzes the power flow solution in bipolar direct current networks with radial structures considering multiple monopolar and bipolar constant power loads. The electrical configuration of the bipolar DC grid considers that the reference pole is non-grounded along the feeder, which produces important neutral currents and voltage imbalances along the DC grid. The power flow problem is formulated through the triangular-based representation of the grid topology, which generates a recursive formulation that allows determining the voltage values in the demand nodes through an iterative procedure. The linear convergence of the triangular-based power flow method is tested through multiple load variations with respect to the nominal grid operative condition. Numerical results in the 21- and the 85-bus grids reveal the relevant variations in the voltage profiles and total grid power losses when the neutral cable is solidly grounded or not.


Introduction
Recently, the direct current (DC) distribution network has gained attention in the research and industrial communities since they do not require generation synchronization, have better voltage profiles, present lower power and energy losses, and in these grids, reactive power and frequency concepts are nonexistent [1,2]. Additionally, the accelerated advances in power electronic converters make these grids attractive for integrating multiple distributed energy resources such as photovoltaic generation, battery energy storage systems, and fuel cells, since these operate directly with DC technology, reducing the number of converters required in comparison with classical alternating current (AC) networks [3][4][5].
Electrical distributions with DC technologies can be constructed with two possible configurations: monopolar and bipolar [6,7]. The former is the most typical structure, where a positive and a neutral pole are used to provide constant voltage to multiple loads connected in the feeder [8]. The latter, i.e., bipolar, is composed of three cables that work as positive, neutral, and negative poles, respectively [9,10]. The main characteristic of the bipolar connection is that it allows the transfer of double the power to the loads in comparison to the monopolar connection solely by including a new cable. Further, under a perfectly balanced load case, the neutral cable can be eliminated, which reduces about 33% of the total investment costs in conductors [11]. An additional advantage of the bipolar connection is that for certain special loads, it is possible to duplicate the voltage profile provided to the load by interconnection between the positive and negative poles [12].
To examine the DC distribution under steady-state conditions, even if this works with monopolar or bipolar configurations with multiple constant power loads, it is mandatory to implement the power flow methodologies owing to the nonlinearities introduced by the constant power terminals [13,14]. In the case of monopolar configurations in the current literature, multiple approaches based on derivative-free and derivative-based methods are found, such as backward/forward, successive approximations, triangularbased, Newton-Raphson, and Taylor-based power flow methods [15]. For bipolar DC grids, few works have recently been developed for power flow and optimal power flow analyses in these networks, some of which are discussed below. Authors of [9] proposed a mixed-integer linear multi-objective optimization model for phase-swapping in the bipolar DC grids; however, they do not consider the existence of constant power loads, which is an oversimplification of the power flow model. Authors of [16] presented an optimal power flow formulation for bipolar DC networks with multiple constant power loads, including dispersed generators; additionally, they incorporated the effect of the neutral pole and its resistance in their analysis. The objective of this work was to compute the locational marginal prices of the network; however, to obtain these, the authors have linearized the hyperbolic constraints regarding constant power loads, which implies that the power flow problem was again simplified, as in the previous work. Authors of [17] illustrated an admittance nodal formulation to solve the power flow problem in bipolar DC grids with constant power terminals. A numerical evaluation into a three-bus system is made using the PSCAD/EMTDC software; however, they do not present any analysis of the model complexity or convergence test under load variations, only altering the constant power terminals for controlled current sources and leaving the complications of the solution to a power system analyzer tool. In reference [18], the authors studied the optimal power flow problem in DC bipolar distribution networks considering high load unbalanced. The authors proposed a linearized optimal power flow model to determine the effect of the grid congestion on the local marginal prices per node. Numerical results in two example test feeders show the effectiveness of the proposed approach; however, the authors do not make any comparison with other alternative optimal power approaches. Authors of [10] presented an optimal power flow formulation with multiple load unbalances by applying the current injection method. Numerical results showed the effectiveness of the proposed formulation by using a sensitivity analysis based on the Jacobian matrix to determine the costs of the total grid unbalance. Their numerical achievements were verified with the help of the PSCAD/EMTDC software.
Owing to the existing gap in the current literature regarding numerical methodologies in solving the power flow problem in bipolar DC grids with constant power terminals, this research proposes the extension of the triangular-based power flow method expounded initially in [16] for bipolar configurations, considering that the neutral pole can or cannot be solidly grounded in all the nodes of the network. A recursive power flow formula is derived, and it can be easily applied for systems with multiple monopolar loads and pure bipolar loads even if high imbalances are present. The numerical results in the 21-and the 85-bus systems portray the negative effect of non-grounding the neutral cable, since important voltage values appear in the neutral pole added with important increments in the total power losses. Moreover, to show the linear convergence properties of the proposed power flow method, multiple load scenarios were considered.
This document is organized as follows: Section 2 presents the derivation of the recursive power flow formula based on the upper-triangular matrix representation; Section 3 unveils the iterative procedure to solve the power flow problem in bipolar DC grids; Section 4 reveals the main characteristics of the 21-bus and the 85-bus grid system. Section 5 details the numerical validations and their analyses and discussions. Finally, Section 6 describes the main concluding remarks and some possible future works derived from this research.

Power Flow Formulation
A bipolar DC distribution network with multiple constant power loads can be analyzed through a recursive formulation that allows the determining of the final voltage profile in the positive p, neutral o, and negative n poles from an initial voltage value provided in the first iteration [19]. To illustrate the formulation of the triangular-based bipolar power flow method, let us consider a small bipolar network composed of three nodes, depicted in Figure 1 [12].   Figure 1 shows that the nodes can be connected with multiple constant power terminals. Note that in an arbitrary node k, it is possible to have loads connected between the positive and neutral poles, i.e., p dk,po , loads between the negative and neutral poles, i.e., p dk,no , and also pure bipolar loads named p dk,pn . It is worth mentioning that the index d is used in referring to demand nodes. Now, if the first Kirchhoff's law is applied to each node, where the current through the line l is named j l,pon = [j l,p , j l,o , j l,n ] and the demanded current at node k is defined as i dk,pon = [i dk,p , i dk,o , i dk,n ] T , the following set of linear equations are formed: which can be easily compacted in a matricial structure with the form defined in Equation (2).
where J l,pon ∈ R 3l×1 is the vector that contains all the branch currents, i d,pon ∈ R 3(b−1)×1 is the vector that comprises all of the demanded currents (b being the total number of nodes/buses the network), and T ll,pon ∈ R 3l×3l is defined as the triangular matrix that relates branches with demand currents. Note that the triangular-based formulation defined in (2) is only applicable to DC bipolar networks with radial structures and a unique voltage controlled source (slack node). Now, if the second Kirchhoff's law is applied to each closed-loop trajectory from the slack node to each demand node, then the following set of linear equations is obtained: which can be compacted in its matricial form as follows: where V d,pon ∈ R 3(b−1)×1 is the vector that contains all the voltage profiles in the positive, neutral, and negatives poles for the demand nodes, V s,pon ∈ R 3s×1 is the vector with the voltage outputs in the slack source (known voltage values), E l,pon R 3l×1 is the vector that consists of all the voltage drops in the grid branches, and A ds ∈ R 3(b−1)×3 is a tridiagonal matrix composed of identity diagonal matrices that relate each demand node with the substation. To relate branch current and branch voltage drops, Ohm's law is applied to each branch l, which produces the following set of linear equations: which can be compacted as defined in Equation (6) E l,pon = R ll,pon J l,pon , where R ll,pon ∈ R 3l×3l is a diagonal matrix that contains all the resistive effects of the grid branches. Now, if Equation (2) is substituted in Equation (6) and its result is replaced in Equation (4), consequently, the triangular-based power flow formula for bipolar DC networks is obtained as defined in (7).
For simplicity in the formulation, the resistance-like nodal matrix, R dd,pon , as T ll,pon R ll,pon T ll,pon , is defined to convert Equation (7) into (8).
Note that Equation (8) is indeed a linear relation between demanded voltages and currents; however, due to the presence of multiple constant power terminals, this equation transformed into a general nonlinear non-convex equation. From Figure 1, it is possible to define the current demanded in an arbitrary node k for each pole as follows [17]: To arrive at a compact formula for the calculation of the demand current at each node, as defined in Equations (9)-(11), let us define the auxiliary variables that measure the voltage difference between poles as follows: v dk,po = v dk,p − v dk,o , v dk,pn = v dk,p − v dk,n , and v dk,no = v dk,n − v dk,o , which can be grouped into a vector as ∆v dk,pon = [v dk,po , v dk,pn , v dk,no ] . With these definitions, the demanded current in the k node can be calculated as presented below: where It is worth mentioning that in finding all the demanded voltages with Equation (8), it is necessary to use an iterative procedure since the demanded currents at each node k defined in (12) are a function of the demanded voltages, which implies that in order to know these voltages, it is necessary to have known them previously, which is only solvable through a recursive implementation, as will be presented in the next section [19].

Iterative Solution
To solve the power flow problem in bipolar DC grids with multiple constant power loads, the iterative procedure implemented is presented in Algorithm 1.
In Algorithm 1, the parameter ζ is the convergence error between two consecutive voltages, which is defined as recommended in [19] as 1 × 10 −10 . Note that once the power flow problem is solved by implementing the Algorithm 1, the total grid power losses in all the branches of the network can be determined as follows: Algorithm 1: Power flow solution using the triangular-based formulation for multiple constant power loads. Data: Select the bipolar DC network under analysis. 1 Transform the data in their per-unique equivalent; 2 Calculate the diagonal resistance matrix R ll,pon ; 3 Calculate the upper-triangular matrix T ll,pon ; 4 Calculate the resistance-like nodal matrix R dd,pon ; 5 Define the voltage input matrix A ds,pon ; 6 Select the maximum iterations' number t max ; 7 Chose the convergence error ζ; 8 Define the slack voltages: V s,pon = [1, 0, −1] ; 9 Make t = 0; 10 Define the initial voltage as V t d,pon = A ds V s,pon ; 11 for t ≤ t max do 12 for k = 2 : n do 13 Calculated the demanded current I t dk,pon using Equation (12); 14 Actualize the demanded voltages V t+1 d,pon using Equation (7); < ζ then 16 for k = 2 : n do 17 Calculate the demanded current I t dk,pon using Equation (12); 18 Calculate and report the branch currents using Equation (2); 19 Report the nodal voltages as V pon = V s,pon ; V t+1 d,pon ;

Test Feeders
To evaluate the efficiency of the proposed triangular-based power flow approach for bipolar DC grids with multiple constant power loads, two test feeders are employed. The characteristic of each one of these test feeders is presented below.

21-Bus System
The 21-bus system was adapted from [19] to include loads connected to the positive and negative poles with respect to the neutral as well as pure bipolar loads. The electrical configuration of this grid is depicted in Figure 2.
For this test feeder, the slack bus works in the positive and negative poles with ±1 kV. All the electrical parameters regarding branches and loads are reported in Table 1. Notably, in this research, we assume that all the poles, including the neutral, were constructed with the same caliber. 22 submitted to Sensors 6 of 13 flow problem is solved by implementing the Algorithm 1, the total grid power losses in all the branches of the network can be determined as follows:

Test feeders
To evaluate the efficiency of the proposed triangular-based power flow approach for 120 bipolar DC grids with multiple constant power loads, two test feeders are employed. The characteristic of each one of these test feeders is presented below.

21-bus system
The 21-bus system was adapted from [19] to include loads connected to the positive and negative poles with respect to the neutral as well as pure bipolar loads. The electrical 125 configuration of this grid is depicted in Figure 2 For this test feeder, the slack bus works in the positive and negative poles with ±1 kV. All the electrical parameters regarding branches and loads are reported in Table 1. Notably, in this research, we assume that all the poles including the neutral were constructed with the same caliber.

85-bus system
The 85-bus system is a DC bipolar adaptation of the IEEE 85-bus system presented in [20] to locate and size fixed-step capacitor banks. This system is operated with 11 kV per pole, and it has radial configuration, i.e., 84 distribution lines. The electrical configuration of this test feeder is depicted in Figure 3.

85-Bus System
The 85-bus system is a DC bipolar adaptation of the IEEE 85-bus system presented in [20] to locate and size fixed-step capacitor banks. This system is operated with 11 kV per pole, and it has radial configuration, i.e., 84 distribution lines. The electrical configuration of this test feeder is depicted in Figure 3.

IEEE-85 bus system
Single line diagram of the system is shown in Fig. 9 and the system data are listed in Table A2, in appendix. The base val-ues are considered as 100 MVA and 11 kV. The total real and reactive power loss for the base case is computed using MATLAB and losses are found to be 316.12 kW and 198.6 kVAr respectively, Operating cost is 1,93,845 $/kWh  The complete parametric information for the 85-bus grid with bipolar structure is presented in Table 2.

Numerical Validations
The numerical implementation of the proposed triangular-based bipolar power flow method was executed in the MATLAB programming environment using the researcher's own scripts. For this implementation, 2021b was used on a PC with an AMD Ryzen 7 3700 2.3 GHz processor and 16.0 GB RAM, running on a 64-bit version of Microsoft Windows 10 Single language.

Results in the 21-Bus System
To demonstrate the effectiveness of the proposed method on addressing the power flow solution in bipolar DC grids with multiple constant power loads, we consider two simulation cases regarding the neutral pole: (i) this pole is solidly grounded at each one of the nodes of the network; and (ii) this pole is only solidly grounded at the substation bus. Table 3 exhibits the voltage in the positive and negative poles when both cases are simulated.  The numerical data provided in Table 3 permit noting the following: (i) in the first simulation case when the neutral is solidly grounded, the positive pole presents some nodes with voltage regulations more than 10%, presenting a minimum voltage at node 17 with a magnitude of 890.1027 V, which implies that the grid regulation is about 10.99%. Conversely, for the negative pole, the worst regulation voltage appears at node 18 with a value of 9.14%. The voltage unbalance between the positive and negative poles is caused by the load disequilibrium, since the positive pole has about 554 kW of constant power consumption while the negative pole has only 445 kW of the same. (ii) In the second case, when the neutral pole is not connected to the electrical ground at each node, it is possible to observe that nodes 15 to 18 have voltages higher than 18 V, which can cause the misoperation of some devices in the adjacent areas, especially when sensitive to voltage variations on the reference point. Additionally, the non-grounded neutral pole can worsen the voltage profile in some nodes of the network, as seen for the positive voltage profile at node 17, where the voltage has decreased 1.8433 V with respect to the first simulation case.
The main issue when the neutral pole is not solidly grounded in all the nodes corresponds to the total grid power losses since these are increased with respect to the solidly grounded case. For example, in the 21-bus system for the first simulation case, the total power losses are 91.2701 kW, while for the second, this number increased to 95.4237 kW. These values indicate that there is about 4.1536 kW that is transformed into heat owing to the currents flowing through the neutral cable.
As evidence of the convergence properties of the studied triangular-based bipolar power flow method, we plot the logarithm value of the convergence error, i.e., log( ) ), for the worst operative case (i.e., non-grounded simulation scenario) considering that the load varies from 60% to 180% of their nominal values listed in Table 1. The convergence rate is depicted in Figure 4. Note that the 100% corresponds to the benchmark case of the network where the load takes its nominal value.  The numerical performance of the triangular-based method for bipolar DC grids in the 21-bus system depicted in Figure 4 detail the following: (i) the load variations directly impact the total number of iterations to reach the desired convergence. For example, when the load is 60% of the nominal value, 10 iterations are required, whereas in the benchmark case, 13 iterations are used, and when the load is 1.8 times the nominal value, 26 iterations are employed to solve the power flow problem. (ii) For all the load increments, it is observed that the convergence rate of the triangular-based bipolar power flow is linear, implying that under normal operating conditions, this method can ensure convergence to the power flow solution. However, this is only possible if the system is operated far from the voltage collapse point.

Results in the 85-Bus Grid
After the implementation of the proposed triangular-based method for bipolar DC grids in the 85-bus system under the nominal operating conditions, we found that (i) when the neutral wire is solidly grounded in all the nodes of the grid, then the total power losses is 452.2981 kW, i.e., 6.76% of the total power consumption, while when the neutral wire is non-grounded, the total power increases to 489.5759 kW, i.e., 7.32% of the net power consumption; (ii) the total number of iterations was about 10 for the solidly-grounded case and 13 for the non-grounded case.
On the other hand, in Figure 5 are presented the voltage profiles in the 85-bus system for the positive and negative poles considering solidly-grounded and non-grounded operation cases.
Voltage profiles in Figure 5 show that (i) for the solidly-grounded operative costs, the minimum voltage in the positive pole is 0.9189 pu at node 54, and for the negative pole it is 0.8950 pu at the same node; and (ii) in the non-grounded scenario, these values are 0.9204 pu for the positive pole, and 0.8925 pu for the negative pole.
Note that when positive and negative voltage profiles are compared (see Figure 5), it is observed that the positive pole has similar behavior to the negative pole, which is very similar to the reflection of the voltage profile with respect to the axis of nodes; even if the system is perfectly balanced, v p must be equal to the absolute value of v n .

Complementary Analysis
The following results were also obtained for both test feeders: In the solidly grounded scenario after 1000 consecutive evaluations of the power flow methodology, the average processing times for the 21-and 85-bus grids were 0.6498 ms and 6.4858 ms, respectively. In the non-grounded case, these times were 0.7318 and 6.7226 ms. Note that the increments in these times are expected since the number of iterations increases in the non-grounded neutral wire simulation scenario. Comparative simulations with the classical backward/forward power flow method adapted for bipolar DC grids demonstrated that both methodologies reach the same voltage profiles and power losses calculations. However, the main advantage of the triangular-based power flow method is the processing time required, since for both tests feeders were 20 % faster than the backward/forward power flow method, which confirmed the results presented by authors of [15] in the case of monopolar DC grid equivalents.

Conclusions and Future Works
The power flow problem in bipolar DC grids with multiple constant power loads was addressed in this research through the application of the triangular-based power flow formulation. The proposed formulation is applicable in radial configurations, considering the possibility that the neutral cable is or is not grounded in all the nodes of the system. The numerical results in the 21-bus system demonstrate that when the neutral pole is solidly grounded, the maximum voltage regulation bound under normal operative conditions was about 10.99% and total grid power losses about 91.2701 kW; conversely, for the nongrounded simulation case, the voltage regulation of the network was 11.17% and total grid losses increased to 95.4237 kW. Numerical evaluations where the total load was incremented from 60% to 180% of the nominal value illustrated that the convergence of the triangular-based power flow method is linear, and the effect of the load increments is transferred to the number of iterations.
In the case of the 85-bus grid, the minimum voltage profile occurs in the negative pole for the solidly-grounded and non-grounded operative scenarios with magnitudes of 0.8950 pu and 0.8925 pu, respectively; while in the case of the power losses, the solidlygrounded case has a total power loss of 452.2981 kW, which increased to 489.5759 kW in the non-grounded scenario.
For future research, it will be possible to execute the following works: (i) propose a recursive power solver for bipolar DC grids by linearizing the demanded currents around the initial voltage values and then updating this point until the desired convergence error is reached; (ii) determine via metaheuristic or convex optimization the optimal subset of nodes where the neutral cable must be solidly grounded to reduce the total grid power losses; and (iii) design an experimental scenario where the power flow methodologies for DC bipolar networks can be validated. Funding: This work was supported in part by the Centro de Investigación y Desarrollo Científico de la Universidad Distrital Francisco José de Caldas under grant 1643-12-2020 associated with the project: "Desarrollo de una metodología de optimización para la gestión óptima de recursos energéticos distribuidos en redes de distribución de energía eléctrica".