Laplacian Matrix-Based Power Flow Formulation for LVDC Grids with Radial and Meshed Conﬁgurations

: There is an increasing interest in low voltage direct current (LVDC) distribution grids due to advancements in power electronics enabling efﬁcient and economical electrical networks in the DC paradigm. Power ﬂow equations in LVDC grids are non-linear and non-convex due to the presence of constant power nodes. Depending on the implementation, power ﬂow equations may lead to more than one solution and unrealistic solutions; therefore, the uniqueness of the solution should not be taken for granted. This paper proposes a new power ﬂow solver based on a graph theory for LVDC grids having radial or meshed conﬁgurations. The solver provides a unique solution. Two test feeders composed of 33 nodes and 69 nodes are considered to validate the effectiveness of the proposed method. The proposed method is compared with a ﬁxed-point methodology called direct load ﬂow (DLF) having a mathematical formulation equivalent to a backward forward sweep (BFS) class of solvers in the case of radial distribution networks but that can handle meshed networks more easily thanks to the use of connectivity matrices. In addition, the convergence and uniqueness of the solution is demonstrated using a Banach ﬁxed-point theorem. The performance of the proposed method is tested for different loading conditions. The results show that the proposed method is robust and has fast convergence characteristics even with high loading conditions. All simulations are carried out in MATLAB 2020b software.


Introduction
The load flow analysis in AC grids has been one of the most studied topics since the 1960s. Load flow analysis in DC grids has received attention rather recently due to the revitalization of low voltage direct current (LVDC) grids [1][2][3][4][5][6][7][8][9][10][11][12]. The DC load flow in conventional power grids should not be confused with the DC load flow in LVDC grids [1]. Load flow in an LVDC grid is that in which there is actual DC present and there are constant power terminals while the term DC load flow is a linearization of the load flow equations in AC grids, which is named this way because of the analogy. Moreover, the load flow in LVDC distribution systems (DS) is also different from its AC counterpart. In LVDC grids, the load flow equations remain in the real domain. The Y bus matrix in LVDC grids is always diagonally dominant and the Z bus is monotone. The unknown variables in LVDC grids are only real power and voltage magnitudes and, finally, load flow methods that are valid for DS with a low Resistance/Reactance (R/X) ratio and weakly meshed structures [13][14][15][16] are not necessarily valid for an LVDC load flow analysis. Thus, the load flow problem in LVDC grids cannot be interpreted with the simplified conventional methods used for its AC counterpart [16]. In an LVDC grid, a power electronic converter can be operated in different modes such as constant current and constant power. In constant power mode, the system of network equations becomes non-linear and requires an iterative power flow for a steady state analysis [17]. Due to the non-linear and non-convex characteristics of 2 of 21 load flow equations, it is often difficult to determine whether a solution for these nonlinear equations is unique and it is not possible to solve these equations analytically with classical circuit methods [5]. The convergence is not always guaranteed with load flow algorithms when the load flow equations are non-linear [15]. An algorithm could even diverge or converge to a non-realistic solution or a solution may not even exist so the uniqueness of the solution cannot be taken for granted. In the literature, different load flow methods have been proposed for LVDC grids based on well-known and studied methodologies for AC grids such as Gauss-Seidel (GS), Newton-Raphson (NR) [17], successive approximation [5,6], loop analysis-based methods [18,19], the backward/forward sweep (BFS) class of methods [20], iterative and linear methods based on the series of Taylor and Laurent [2,5,8] and the triangular matrix-based method [12]. Reference [21] implemented the convex relaxation and linear approximation formulation for optimal power flow in an AC/DC transmission system. A confined exact power flow formulation was proposed in [22] for a DC distribution system based on a model predictive control methodology incorporating the dynamics of energy storage systems and power output from renewables. For a gentle introduction to the relaxation and approximation of the power flow equations readers are referred to the survey presented in [23]. These methods can solve load flow problem via interior point methods but they generate a quadratic increment of variables to be analyzed that increase the processing time. The linear methods proposed in the literature sacrifice precision over speed [2,5,8]. In [10], a review of classical and emerging methods for a load flow solution in DC grids is presented. Reference [3] reports a comparative study on the different methods of power flow solutions in LVDC grids considering processing time and numerical convergence. Reference [4] proposes an approximate power flow solution method for LVDC grids based on a logarithmic transform of voltage magnitudes. In general, the studies presented in the literature are only for radial DS and the convergence is taken for granted. The convergence of the GS method and the successive approximation [6,11] were proven with the Banach space theorem while the convergence of the NR method [7,17] was proven with the Kantorovich theorem. The method proposed in [6] requires more iterations as opposed to NR [5]. The conditions presented in [17] are not only sufficient but also necessary so the presented results are conservative. If any of the conditions are not met, the solution may diverge. The convergence in heavily loaded conditions is not guaranteed. The BFS class of solvers for power flow solutions are widely used for radial and weakly meshed configurations [14,24]. Although BFS-based methods have been favored for power flow solutions in DS due to their simplicity, these methods present challenges in handling highly meshed distribution networks and convergence characteristics [25].
A generic methodology based on the application of the Newton method to networks formulated with a modified nodal analysis that is expanded to accommodate various component models is presented in [26]. This methodology can handle arbitrary DS topologies and devices with arbitrary constraints and variables such as load tap changers, regulators [27] and induction machines [28]. It converges much faster than the BFS class of solvers and does not have their topology limitations. On the other hand, it was only demonstrated for AC DS. A load flow solution using graph theory-based techniques is also possible [25,29,30]. The graph-based method considered in these references requires the recalculation of the matrices proposed for the load flow solution during the iterative process, which increases the processing time. References [10,11] demonstrated that improved graph-based methods for radial networks had the best processing time compared with classical methods especially when the number of nodes were increased in the grid. Reference [19] proposed a method based on a loop analysis and its linear formulation for radial and meshed LVDC grids hosting different load types where a formulation based on the injected power instead of injected currents was derived. However, the approximations made to speed up the convergence of the proposed method could lead to errors in the solution with heavy load conditions. The proof of convergence was not provided for the proposed method. The convergence proof of a graph theory-based method [12] is Energies 2021, 14, 1866 3 of 21 addressed in [9,11]; however, this method is only valid for radial configurations. Meshed networks have many advantages over radial networks in terms of voltage profile, reliability and losses [31]. Most of the studies in the literature for DC grids focus on radial networks and the convergence is taken for granted. Moreover, AC networks are often transformed to their DC counterparts by simply considering the resistance of lines without taking into account the skin effect [2,5,8,12,32]. Note that DC resistance is less than AC resistance for a given cross-sectional area of a conductor. This paper proposes a load flow solution for LVDC grids based on a graph theory for both radial and meshed configurations. Although graph-based methods have been previously reported in the literature for AC networks [24,[33][34][35][36][37], their application to LVDC grids considering a meshed configuration has not been studied. In the proposed method, first a formulation is developed for radial networks then it is extended to meshed networks. The main contributions of this paper are as follows:

1.
A novel graph-based method for a load flow solution in LVDC grids hosting different load types including a constant power load (CPL) for both radial and meshed configurations is proposed. Any network configuration can be modelled with the proposed method with only two incidence matrices by taking advantage of the DS structure. The only input required by the proposed algorithm is the conventional node to branch oriented data. The incidence matrices are constant matrices and totally depend on the network configuration, i.e., there is no need to recalculate these matrices in the iterative process. This enhances the performance of the method in terms of the processing time.

2.
As line parameters and load characteristics are modelled in separate matrices, modifications in the network can be performed easily without the need to recalculate the bus matrix. In addition, Distributed Generation (DG) integration is also accounted for in the formulation. The equivalent DC resistance of lines is used to achieve better accuracy levels.

3.
The uniqueness of the solution of the proposed method is proven by contraction mapping using the Banach fixed-point theorem.

4.
Different loading conditions and load types are considered to validate the robustness of the proposed method. In addition, the results are compared with the direct load flow (DLF) method proposed in [24] as a benchmark. Other than benchmark results readers are referred to references [25,33,38], which provide information about the convergence of different power flow algorithms and their processing time.
The rest of the paper is organized as follows. In Section 2 the formulation of the proposed method is derived and in Section 3 the proof for the uniqueness of the solution of the proposed method is provided. Simulations and results are presented in Section 4 while the convergence analysis of the proposed method and the potential applications of the LVDC grids are discussed in Section 5. Finally, conclusions are drawn in Section 6. Note all of the simulations were carried out in MATLAB 2020b using a PC with the following specifications: CPU: Intel core i7 @ 3.21 and 3.19 GHz, 16 GB RAM, 64-bit operating system with Windows 10.

Formulation of the Proposed Method
The formulation of the load flow analysis for the LVDC grid presented in this paper is based on a graph theory with following conditions. C-1: The graph is connected and there is no islanding in the feeder, which guarantees that the bus matrix is a non-singular matrix.
C-2: There is at least one CPL and there is one constant voltage node present in the system. Moreover, the constant voltage node has the capability to provide the combined power demand of the loads and network losses. C-3: In the steady state operation of the LVDC grid, voltages remain within the boundary of (0 < v min ≤ v ≤ v max ). This condition is required for voltage regulation and stability.
Energies 2021, 14, 1866 4 of 21 C-4: Short circuit currents > normal operating currents (this condition is a useful observation for any electrical power system and helpful in the proof of convergence).

Graph Theory
Let us consider an LVDC grid with n number of nodes, b number of branches and number of loops; for this grid, let us also define the primitive resistance matrix as follows: where R ij is the resistance of the branch connected between node i and node j and the p ∈ R bxb space. The DC resistance can be calculated from the AC resistance as follows [39]: Here R dc is DC resistance, R ac is AC resistance, δ is the skin depth of the conductor and r is the conductor radius. The value of the skin depth is frequency dependent. The studied systems in this work are 50 Hz systems. Figure 1 shows the graph for an arbitrary simple LVDC grid with six nodes, five branches and two loops. By taking advantage of a graph theory for a DC resistive network, a relationship can be developed between branch currents B ∈ R bx1 and nodal injected currents I ∈ R (n−1)x1 with an incidence matrix H ∈ R bx(n−1) . The algorithm to construct the incidence matrix is as follows.
Step 1. Construct the matrix H ∈ R bxn and fill it with zeros.
Step 2. Bool (H n (i,j)) = 1 if there is a path of length 1 ≤ n from node i to node j. Bool (H n (i,j)) = 0 otherwise where H n is the number of walks of length n from node i to node j.
Step 3. Remove the first column from matrix H (i.e., the slack node), which constitutes H ∈ R bx(n−1) . C-3: In the steady state operation of the LVDC grid, voltages remain within the boundary of (0 < vmin ≤ v ≤ vmax). This condition is required for voltage regulation and stability.
C-4: Short circuit currents > normal operating currents (this condition is a useful observation for any electrical power system and helpful in the proof of convergence).

Graph Theory
Let us consider an LVDC grid with n number of nodes, b number of branches and l number of loops; for this grid, let us also define the primitive resistance matrix as follows: where Rij is the resistance of the branch connected between node i and node j and the space. The DC resistance can be calculated from the AC resistance as follows [39]: Here Rdc is DC resistance, Rac is AC resistance, δ is the skin depth of the conductor and r is the conductor radius. The value of the skin depth is frequency dependent. The studied systems in this work are 50 Hz systems. Figure 1 shows the graph for an arbitrary simple LVDC grid with six nodes, five branches and two loops. By taking advantage of a graph theory for a DC resistive network, a relationship can be developed between branch currents x1 B b   and nodal injected currents . The algorithm to construct the incidence matrix is as follows.
Step 1. Construct the matrix and fill it with zeros.
Step 2. Bool (H n (i,j)) = 1 if there is a path of length 1 ≤ n from node i to node j. Bool (H n (i,j)) = 0 otherwise where H n is the number of walks of length n from node i to node j.
Step 3. Remove the first column from matrix H (i.e., the slack node), which constitutes

Formulation for Radial Networks
Removing the dashed lines in Figure 1 will leave the system with a radial configuration. The construction of an incidence matrix by applying Kirchhoffʹs Current Law (KCL) can be formed as: = + + + + B I I I I I .
(3) Figure 1. Graph of a simple low voltage direct current (LVDC) distribution system.

Formulation for Radial Networks
Removing the dashed lines in Figure 1 will leave the system with a radial configuration. The construction of an incidence matrix by applying Kirchhoff's Current Law (KCL) can be formed as: B 1 = I 2 + I 3 + I 4 + I 5 + I 6 . (3) Energies 2021, 14, 1866

Load Model
For node k its demand current can be written as follows: where r, i and p stand for constant resistance, constant current and constant power load, respectively, and α, β and γ are the coefficients of the constant resistance, constant current and constant power load current percentage of consumption at node k, respectively. It is important to note that the sum of α + β + γ = 1; e.g., for a constant power load α = β = 0 and γ = 1.
Equation (11) in matrix form: We considered DGs with active and reactive power (PQ) generation constraints in our study.
The net current for node k for a Resistance (R), Current (I), Power (P) (RIP) load model can be calculated as follows: where P k is the total injected power at node k. Including the generated power in the formulation turns the network from passive to active and it is possible to include DGs in the system.

Power Flow Formulation
The voltage drops on each branch are calculated as follows: where ∆V ∈ R bx1 is the vector that contains the voltage drop of all branches. The slack node voltage (V s ) is known (by default it equals to 1 per unit (p.u)). For the remaining nodes, the voltage drop from a given node to the slack node can be calculated by applying the Kirchhoff's Voltage Law KVL on the closed paths in between the node and the slack node. [ Note that matrix H is an upper triangular matrix so its transpose will be a lower triangular matrix. A combination of (9), (12), (15) and (16) gives the following equation. Let us define: Here, Φ is a Laplacian matrix that is weighted by the branch resistances of the network. For a connected network, Φ is always non-singular whether the incidence matrix H is a square matrix or not. [ Equation (19) is a non-linear expression for the power flow solution of a radial LVDC grid due to having a CPL in the system according to C-2, which means it can only be solved by an iterative numerical method. The proposed iterative solution updates the V vector as follows: where t is the iteration count.

Formulation for Meshed Networks
Now let us consider Figure 1 with the dashed lines (tie-lines) connected, which transforms the network architecture from radial to meshed. A new incidence matrix is required to model the loops created by the tie-lines. Loops do not affect the node current injection but new branch currents need to be added into the system. Taking the new branches into account, the current injections of nodes 3, 5 and 6 can be written as: By incorporating the loops in the formulation, the following expression is achieved: where L ∈ R bx is the second incidence matrix that contains the information of loops, B ∈ R x1 is a vector of the tie-line currents and is the number of tie-lines. Matrix L can be constructed with the following algorithm: Step 1. Construct a matrix of dimension (L ∈ R bx ) and fill it with zeros.
Step 2. For a tie-line between node i and j, subtract column i in H from column j and place the result in column of L.
Step 3. Repeat step 2 for all tie-lines.
For a system shown in Figure 1 the matrix L is: The voltage drop in the tie-lines is calculated by applying KCL in each fundamental loop as follows: By using the transpose of the second incidence matrix, a relationship is built for tie-line voltage drops as follows:  (29) where R is the diagonal matrix of the tie-line resistances. It is possible to determine tie-line currents based on branch currents by rearranging (29) as follows: . (30) [S] [H] where U b is the identity matrix of dimension b. C-1 guarantees that the primitive resistive matrix is a non-singular matrix so its inverse exists. The relationship vector B between branch currents to the nodal injection will take the following form: Finally, the load flow equations for a meshed configuration can be written as follows: The load flow expression in (33) is valid for both radial and meshed configurations, i.e., for a radial network, matrix S will be equal to matrix H. Hence, the proposed formulation can handle a variety of distribution network configurations. If we observe (33), matrix H and matrix S are constant matrices that will remain constant during the iterative process ultimately reducing the processing time. Both matrices depend on the network configuration and only need to be built once.
The total power loss of the network can be calculated as follows: Table 1 shows the pseudo-code of the proposed power flow algorithm. Algorithm 1 starts with a flat start while Algorithm 2 starts with a random initial guess. As one of the reasons for the divergence of the algorithms is that when they are provided with a bad initial guess, to check the robustness of the proposed algorithm we tested the algorithm with a random initial guess. In both algorithms l is the number of loops, i.e., for a radial network the value of l will be zero and for a meshed network its value will be equal to the total number of loops.  Table 1 shows the pseudo-code of the proposed power flow algorithm. Algorithm 1 starts with a flat start while Algorithm 2 starts with a random initial guess. As one of the reasons for the divergence of the algorithms is that when they are provided with a bad initial guess, to check the robustness of the proposed algorithm we tested the algorithm with a random initial guess. In both algorithms l is the number of loops, i.e., for a radial network the value of l will be zero and for a meshed network its value will be equal to the total number of loops. Table 1. Proposed load flow algorithms (with a flat start and with a random initial guess). Calculate t I from (34);

Proposed Power Flow Algorithm
starts with a flat start while Algorithm 2 starts with a random initial guess. As one of the reasons for the divergence of the algorithms is that when they are provided with a bad initial guess, to check the robustness of the proposed algorithm we tested the algorithm with a random initial guess. In both algorithms l is the number of loops, i.e., for a radial network the value of l will be zero and for a meshed network its value will be equal to the total number of loops. Table 1. Proposed load flow algorithms (with a flat start and with a random initial guess). Calculate t I from (34);

Convergence Proof
A classical reference for the convergence that includes theorems such as contraction mapping is [40], which was used in our analysis. We used the Banach fixed-point theorem to prove the existence and uniqueness of the solution. Although the Banach fixed-point theorem is a well-explored theorem in functional analysis, it has not been fully explored in a power system analysis yet. It is noteworthy that contraction mapping has only a fixed-point and according to the Banach fixed-point theorem any contraction mapping on a non-empty metric space has a unique fixed-point. Because of that, a non-linear function converges to that unique point with an iterative process.

Remark 1.
The matrix Φ is a Laplacian matrix that is weighted by network resistances and has the characteristics of a diagonal dominant and a positive definite so we are sure that |Φ ii | ≥ Φ ij , ∀i = j.
Equation (36) in terms of node voltages can be written as: Equation (38) is a recursive formulation of the load flow problem; we took advantage of the Banach fixed-point theorem representation [41][42][43] to prove the existence and uniqueness of the solution. Theorem 1. The recursive load flow formulation in (38) is stable and contraction mapping can be formed.
For any V that is contained in closed box (0 < v min ≤ v ≤ v max ) according to C-3, regardless of the initial condition V 0 , let X ∈ R (n−1)x1 be the solution of the load flow such that: Here ψ is a contraction constant also known as the Lipschitz constant whose value is (0 < ψ < 1) and ≤ is called the Lipschitz inequality. (40) is a recursive function like the Banach fixed-point theorem, which is a mapping of V to itself. According to the Banach fixed-point theorem, the solution X of the load flow problem satisfying X = f (X) exists and is unique if and only if f (V) is contraction mapping on V. The fixed-point theorem is valid for any Banach algebra so from expressions (39) and (40) the following result can be deduced:

Proof. The expression in
where: As mentioned before, Φ is a diagonal dominant matrix so we can guarantee that Φ ≤ max{φ ii }, which allows us to rewrite (42) as follows: where φ ii is the Thevenin equivalent resistance at each node except the slack node. The Lipschitz constant can then be expressed in terms of the maximum allowable load and short circuit current, which has a physical sense for power systems.
Remark 2. The proof is completed here. From expression (44) we can conclude that the value of the contraction constant will be less than one as from the earlier defined condition C-4 (short circuit current > normal operation current). (41) and (43) one can conclude that the system will collapse near ψ = 1 . In addition, a smaller (Φ * P) means a smaller ψ, which implies a stronger network that allows a larger permissible region for the proposed method with a fast convergence. Therefore, the Lipschitz constant has a physical meaning that is directly related to the network configuration, loading conditions and voltage limitations defined for a given network.

Method Validation
There is no (IEEE) benchmark feeder available for LVDC grids. Therefore, we modified two AC IEEE test feeders to validate the performance of the proposed method, i.e., 33-and 69-node test feeders. These test feeders have been used in the literature for LVDC grids in different aspects. For instance, in [3,4] these feeders were modified to test a load flow solution in LVDC grids. These test feeders have been transformed in the literature to study solvers proposed for LVDC grids. In the literature, these test feeders have been considered as hosting only CPLs but in our study we considered different load types, i.e., Case 1: all nodes having a CPL; Case 2: a RIP load with different percentages of α, β and γ with at least one node having a CPL. We also tested the robustness and processing performance of the proposed method with heavy loading conditions. The proposed method was compared with the DLF approach [24] applied to DC radial and meshed networks. The DLF approach is a fixed-point approach like the BFS class of solvers but uses connectivity matrices that enables the handling of meshed networks more easily. The main difference between the DLF algorithm and the proposed algorithm is the way (15) and (33) were obtained. In the DLF method, the power flow solution is achieved in two sweeps, i.e., in forward sweep the relationship between the bus injection to the branch current is achieved while in the backward sweep voltage drops are computed using the branch current to branch voltage relationship. In the case of the proposed algorithm this was done only in one sweep with the help of the Laplacian matrix.

Modified 33-Node Test Feeder
This test feeder was basically an adoption of an AC test feeder used for power loss reduction with an optimal location of DGs in [44]. We transformed this feeder from AC to DC by neglecting the reactive power and reactance of all branches. We calculated the resistive elements of the branches with (2). The base values for this test feeder were as follows: V base = 12.66 kV, P base 10 MW, R base = 16.02756 Ω. We considered three DG units at node 12, 15 and 31; the optimal location and values of the DGs was considered according to the study conducted in [44]. The test feeder is shown in Figure 2 with a radial and meshed configuration and its parameters are given in the Appendix A (Table A1). Note for Case 1 the values of α = β = 0 and γ = 1 for all demand nodes while the values of α, β and γ for Case 2 are given in Appendix A (Table A1). to the study conducted in [44]. The test feeder is shown in Figure 2 with a radial and meshed configuration and its parameters are given in the Appendix A (Table A1). Note for Case 1 the values of α = β = 0 and γ = 1 for all demand nodes while the values of α, β and γ for Case 2 are given in Appendix A (Table A1).

Modified 69-Node Test Feeder
This test feeder was also an adoption of the AC test feeder used for power loss reduction with an optimal location of DGs in [44]. We transformed this test feeder from AC to DC in a similar way as we transformed the 33-node test feeder with the same base values. The test feeder is shown in Figure 3 with a radial and meshed configuration. We created four loops for the 69-node test feeder. The relevant data are given in the Appendix A (Table A2). The tie-line data of both test feeders are given in the Appendix A (Table A3).

Modified 69-Node Test Feeder
This test feeder was also an adoption of the AC test feeder used for power loss reduction with an optimal location of DGs in [44]. We transformed this test feeder from AC to DC in a similar way as we transformed the 33-node test feeder with the same base values. The test feeder is shown in Figure 3 with a radial and meshed configuration. We created four loops for the 69-node test feeder. The relevant data are given in the Appendix A (Table A2). The tie-line data of both test feeders are given in the Appendix A (Table A3).

Modified 69-Node Test Feeder
This test feeder was also an adoption of the AC test feeder used for power loss reduction with an optimal location of DGs in [44]. We transformed this test feeder from AC to DC in a similar way as we transformed the 33-node test feeder with the same base values. The test feeder is shown in Figure 3 with a radial and meshed configuration. We created four loops for the 69-node test feeder. The relevant data are given in the Appendix A (Table A2). The tie-line data of both test feeders are given in the Appendix A (Table A3).  Figure 4 shows the voltage profile of the modified 33-node test feeder for both the radial and meshed configurations. As expected by adding loops in the network, the voltage profile was improved compared with the radial configuration. We also plotted the results by using the AC resistance of the network as is considered in the literature and with actual equivalent DC resistances of the branches. As mentioned earlier, classical methods such as NR, fast decoupled, GS and BFS class solvers for load flow solutions used for transmission networks have convergence issues when used for DS. DS with special characteristics, i.e., a shorter length of lines, a high R/X ratio and tight voltage restrictions require more accurate study for a load flow solution. By considering the actual resistance of the lines for the DC current makes this study more accurate.  Figure 4 shows the voltage profile of the modified 33-node test feeder for both the radial and meshed configurations. As expected by adding loops in the network, the voltage profile was improved compared with the radial configuration. We also plotted the results by using the AC resistance of the network as is considered in the literature and with actual equivalent DC resistances of the branches. As mentioned earlier, classical methods such as NR, fast decoupled, GS and BFS class solvers for load flow solutions used for transmission networks have convergence issues when used for DS. DS with special characteristics, i.e., a shorter length of lines, a high R/X ratio and tight voltage restrictions require more accurate study for a load flow solution. By considering the actual resistance of the lines for the DC current makes this study more accurate.   Figure 5 shows the voltage profile of the 33-node test feeder with different load types both for radial and meshed configurations. The results showed that the proposed method was as accurate as DLF but with a faster convergence. The convergence performance is discussed in the next section. Note that the DLF method considered here for comparison has been proven to be the fastest in the literature among other classical methods for a power flow solution in DC grids [3,10] so the comparison with classical methods was out of scope in this study. Figure 6 shows the voltage profile of the 33-node test feeder with different loading conditions, i.e., 300% increase in the nominal load with 100% increase at each step.  Figure 5 shows the voltage profile of the 33-node test feeder with different load types both for radial and meshed configurations. The results showed that the proposed method was as accurate as DLF but with a faster convergence. The convergence performance is discussed in the next section. Note that the DLF method considered here for comparison has been proven to be the fastest in the literature among other classical methods for a power flow solution in DC grids [3,10] so the comparison with classical methods was out of scope in this study. Figure 6 shows the voltage profile of the 33-node test feeder with different loading conditions, i.e., 300% increase in the nominal load with 100% increase at each step. was as accurate as DLF but with a faster convergence. The convergence performance is discussed in the next section. Note that the DLF method considered here for comparison has been proven to be the fastest in the literature among other classical methods for a power flow solution in DC grids [3,10] so the comparison with classical methods was out of scope in this study. Figure 6 shows the voltage profile of the 33-node test feeder with different loading conditions, i.e., 300% increase in the nominal load with 100% increase at each step.  With a 300% increase in load, the minimum voltage for the weakest node (node 18) in the network was 0.8440 p.u. for a radial system while for the meshed configuration the value was 0.9200 p.u., which was still within the limit. Figure 7 shows the voltage profile comparison of the 33-node test feeder with and without DG penetration for both radial and meshed configurations. Here, one thing worth mentioning is that in the case of the radial configuration, the voltage at node 15 was greater than 1 p.u. while in the case of the meshed configuration it was less than 1 p.u. The reason is that in the meshed configuration power was more equally distributed in the system as it had multiple paths of current flow. With a 300% increase in load, the minimum voltage for the weakest node (node 18) in the network was 0.8440 p.u. for a radial system while for the meshed configuration the value was 0.9200 p.u., which was still within the limit. Figure 7 shows the voltage profile comparison of the 33-node test feeder with and without DG penetration for both radial and meshed configurations. Here, one thing worth mentioning is that in the case of the radial configuration, the voltage at node 15 was greater than 1 p.u. while in the case of the meshed configuration it was less than 1 p.u. The reason is that in the meshed configuration power was more equally distributed in the system as it had multiple paths of current flow.

Simulation Results
We observed similar results by increasing the load and DG penetration from a nominal up to 300% with a 100% increase at each step, as seen in in Figure 8. Figure 9 shows the voltage profile of the 69-node test feeder at a nominal load both for radial and meshed configurations considering R ac and R dc .
With a 300% increase in load, the minimum voltage for the weakest node (node 18) in the network was 0.8440 p.u. for a radial system while for the meshed configuration the value was 0.9200 p.u., which was still within the limit. Figure 7 shows the voltage profile comparison of the 33-node test feeder with and without DG penetration for both radial and meshed configurations. Here, one thing worth mentioning is that in the case of the radial configuration, the voltage at node 15 was greater than 1 p.u. while in the case of the meshed configuration it was less than 1 p.u. The reason is that in the meshed configuration power was more equally distributed in the system as it had multiple paths of current flow. We observed similar results by increasing the load and DG penetration from a nominal up to 300% with a 100% increase at each step, as seen in in Figure 8.

Performance and Convergence
To test the performance and convergence of the proposed method, simulations at

Performance and Convergence
To test the performance and convergence of the proposed method, simulations at different loading conditions were performed. We did a comparative analysis of the pro-

Performance and Convergence
To test the performance and convergence of the proposed method, simulations at different loading conditions were performed. We did a comparative analysis of the proposed method with DLF in terms of processing time and the number of iterations required to find the final solution. For the convergence test, we considered different random initial points (i.e., between 0.5 and 1.5 p.u.) and observed the behavior of the proposed method in terms of convergence. Figure 10 shows the comparison of the solution time of the proposed method relative to the DLF methodology for different loading conditions. The number of iterations required to find the final solution is shown in Figure 11. The results suggested that the proposed method was faster than the DLF approach for the same level of accuracy. The proposed method did not sacrifice accuracy over speed even in high loading conditions. We performed 1000 simulations for each case and calculated the average value to be more precise in the evaluation of the solution time. quired to find the final solution is shown in Figure 11. The results suggested that the proposed method was faster than the DLF approach for the same level of accuracy. The proposed method did not sacrifice accuracy over speed even in high loading conditions. We performed 1000 simulations for each case and calculated the average value to be more precise in the evaluation of the solution time.

Convergence
To test the convergence performance of the proposed methods we considered 20 random initial points (between 0.5 to 1.5 p.u.) and observed the maximum error in each case at a nominal load as well as at 300% loading, as shown in Figure 12. The results showed that the solution converged (voltage mismatch between two consecutive iterations up to four decimal places, i.e., 10 −4 ) in the three iterations for all random initial values even un- quired to find the final solution is shown in Figure 11. The results suggested that the proposed method was faster than the DLF approach for the same level of accuracy. The proposed method did not sacrifice accuracy over speed even in high loading conditions. We performed 1000 simulations for each case and calculated the average value to be more precise in the evaluation of the solution time.

Convergence
To test the convergence performance of the proposed methods we considered 20 random initial points (between 0.5 to 1.5 p.u.) and observed the maximum error in each case at a nominal load as well as at 300% loading, as shown in Figure 12. The results showed that the solution converged (voltage mismatch between two consecutive iterations up to four decimal places, i.e., 10 −4 ) in the three iterations for all random initial values even un-

Convergence
To test the convergence performance of the proposed methods we considered 20 random initial points (between 0.5 to 1.5 p.u.) and observed the maximum error in each case at a nominal load as well as at 300% loading, as shown in Figure 12. The results showed that the solution converged (voltage mismatch between two consecutive iterations up to four decimal places, i.e., 10 −4 ) in the three iterations for all random initial values even under high loading conditions, which meant that the proposed method was very robust, i.e., it guaranteed a fast convergence even in poor conditions. In addition, as aforementioned in Section 3, we provided the convergence proof and uniqueness of the solution using contraction mapping, which has a real physical meaning in an electrical system. From the value of the contraction constant, we analyzed the behavior of the system. in Section 3, we provided the convergence proof and uniqueness of the solution using contraction mapping, which has a real physical meaning in an electrical system. From the value of the contraction constant, we analyzed the behavior of the system.  Figure 13 shows the allowable loading range for 33 and 69 test feeders at vmin of 0.85 p.u. The results showed that for the 33-node test feeder, the maximum allowable loading was 300% of the nominal value while in the case of the 69-node test feeder the value was 180% of the nominal value. It is noteworthy that the above values were for systems with AC branch resistances (Rac). By considering Rdc we achieved more accurate results, as shown in Figure 14. The value of the contraction constant for Rdc was less than from the corresponding value for Rac, which meant a stronger network that could ultimately bear a greater load without violating the voltage restrictions. Figure 14 shows the value of the contraction constant for each node for the 33-node test feeder at different loading conditions. Nodes 24 and 25 had a maximum value at 600%  Figure 13 shows the allowable loading range for 33 and 69 test feeders at v min of 0.85 p.u. The results showed that for the 33-node test feeder, the maximum allowable loading was 300% of the nominal value while in the case of the 69-node test feeder the value was 180% of the nominal value. It is noteworthy that the above values were for systems with AC branch resistances (R ac ).
Energies 2021, 14, x FOR PEER REVIEW 16 of 24 in Section 3, we provided the convergence proof and uniqueness of the solution using contraction mapping, which has a real physical meaning in an electrical system. From the value of the contraction constant, we analyzed the behavior of the system.  Figure 13 shows the allowable loading range for 33 and 69 test feeders at vmin of 0.85 p.u. The results showed that for the 33-node test feeder, the maximum allowable loading was 300% of the nominal value while in the case of the 69-node test feeder the value was 180% of the nominal value. It is noteworthy that the above values were for systems with AC branch resistances (Rac). By considering Rdc we achieved more accurate results, as shown in Figure 14. The value of the contraction constant for Rdc was less than from the corresponding value for Rac, which meant a stronger network that could ultimately bear a greater load without violating the voltage restrictions. By considering R dc we achieved more accurate results, as shown in Figure 14. The value of the contraction constant for R dc was less than from the corresponding value for R ac , which meant a stronger network that could ultimately bear a greater load without violating the voltage restrictions. loading, i.e., above this loading condition there would be a voltage violation at these two nodes. Similarly, for the 69-node feeder, node 61 had a maximum value at 300% loading, as shown in Figure 15. There are different approaches to finding the maximum loadability of the network; for instance, the authors in [45] proposed a sequential convex optimization method for the maximum loadability problem for meshed networks based on a semi-definite relaxation approach but the proposed approach is not valid for large scale systems and also the performance was sensitive to the selection of the penalty parameter. To stress the convergence of the proposed method, we increased and decreased the branch length five times, respectively. The results in Figure 16 show that the proposed method converged for both cases in less than three iterations.  Figure 14 shows the value of the contraction constant for each node for the 33-node test feeder at different loading conditions. Nodes 24 and 25 had a maximum value at 600% loading, i.e., above this loading condition there would be a voltage violation at these two nodes.
Similarly, for the 69-node feeder, node 61 had a maximum value at 300% loading, as shown in Figure 15. There are different approaches to finding the maximum loadability of the network; for instance, the authors in [45] proposed a sequential convex optimization method for the maximum loadability problem for meshed networks based on a semidefinite relaxation approach but the proposed approach is not valid for large scale systems and also the performance was sensitive to the selection of the penalty parameter.
To stress the convergence of the proposed method, we increased and decreased the branch length five times, respectively. The results in Figure 16 show that the proposed method converged for both cases in less than three iterations.  Similarly, for the 69-node feeder, node 61 had a maximum value at 300% loading, as shown in Figure 15. There are different approaches to finding the maximum loadability of the network; for instance, the authors in [45] proposed a sequential convex optimization method for the maximum loadability problem for meshed networks based on a semi-definite relaxation approach but the proposed approach is not valid for large scale systems and also the performance was sensitive to the selection of the penalty parameter. To stress the convergence of the proposed method, we increased and decreased the branch length five times, respectively. The results in Figure 16 show that the proposed method converged for both cases in less than three iterations.

Potential Applications of LVDC Grids and Challenges
The International Electrotechnical Commission (IEC) conducted a survey in 2016 on the potential application of such grids considering a market assessment from different organizations and industries worldwide [46]. Considering the various aspects, LVDC grids are likely to have a profound impact not only on developing economies by providing the framework for enabling electricity access even in the remotest locations but also in developed economies where LVDC is already seen as a solution toward greener and more sustainable energy. In developed countries, the main driver for moving towards LVDC grids is energy efficiency. Although the concept of High Voltage Direct Current (HVDC) for transmission systems is matured enough, LVDC for DS is a new concept and still there are many challenges such as a lack of standardization, safety issues, a fault analysis, grid structures and reactive power support in case the system has AC loads that require reactive power support. These challenges are holding back the implementation of LVDC grids; further research is needed in this area before the implementation of such grids.

Conclusions
This paper proposed a power flow solution based on the use of a Laplacian matrix for LVDC grids both for radial and meshed configurations. The proposed method was faster than the DLF methodology for the same level of tolerance. The robustness and performance of the proposed method was tested for different types and loading conditions. The proposed method presented fast convergence characteristics without sacrificing the accuracy even at heavy loading conditions. In addition, the uniqueness of the solution was proven with a Banach fixed-point theorem. The non-linear mapping (contraction constant) used in the proof has a physical meaning in electrical systems, i.e., it gives information about the system state at which the proposed method guarantees that the solution will converge and it will be unique. The proposed method is useful for LVDC grids and for DC microgrids as they have DGs that require real-time monitoring and a faster power flow algorithm. The formulation is suitable for an extension to study LVDC networks with multiple slack nodes and with only one adjacency matrix for both radial and meshed networks. Moreover, it can be also integrated into existing solvers to develop a unified power flow method for a hybrid AC/DC system. Moreover, it will be possible to expand this formulation to accommodate various DS component models in the power flow solution such as a solid-state transformer in the case of hybrid AC/DC systems.

Potential Applications of LVDC Grids and Challenges
The International Electrotechnical Commission (IEC) conducted a survey in 2016 on the potential application of such grids considering a market assessment from different organizations and industries worldwide [46]. Considering the various aspects, LVDC grids are likely to have a profound impact not only on developing economies by providing the framework for enabling electricity access even in the remotest locations but also in developed economies where LVDC is already seen as a solution toward greener and more sustainable energy. In developed countries, the main driver for moving towards LVDC grids is energy efficiency. Although the concept of High Voltage Direct Current (HVDC) for transmission systems is matured enough, LVDC for DS is a new concept and still there are many challenges such as a lack of standardization, safety issues, a fault analysis, grid structures and reactive power support in case the system has AC loads that require reactive power support. These challenges are holding back the implementation of LVDC grids; further research is needed in this area before the implementation of such grids.

Conclusions
This paper proposed a power flow solution based on the use of a Laplacian matrix for LVDC grids both for radial and meshed configurations. The proposed method was faster than the DLF methodology for the same level of tolerance. The robustness and performance of the proposed method was tested for different types and loading conditions. The proposed method presented fast convergence characteristics without sacrificing the accuracy even at heavy loading conditions. In addition, the uniqueness of the solution was proven with a Banach fixed-point theorem. The non-linear mapping (contraction constant) used in the proof has a physical meaning in electrical systems, i.e., it gives information about the system state at which the proposed method guarantees that the solution will converge and it will be unique. The proposed method is useful for LVDC grids and for DC microgrids as they have DGs that require real-time monitoring and a faster power flow algorithm. The formulation is suitable for an extension to study LVDC networks with multiple slack nodes and with only one adjacency matrix for both radial and meshed networks. Moreover, it can be also integrated into existing solvers to develop a unified power flow method for a hybrid AC/DC system. Moreover, it will be possible to expand this formulation to accommodate various DS component models in the power flow solution such as a solid-state transformer in the case of hybrid AC/DC systems.