Voltage Stability Analysis in Medium-Voltage Distribution Networks Using a Second-Order Cone Approximation

This paper addresses the voltage stability margin calculation in medium-voltage distribution networks in the context of exact mathematical modeling. This margin calculation is performed with a second-order cone (SOCP) reformulation of the classical nonlinear non-convex optimal power flow problems. The main idea around the SOCP approximation is to guarantee the global optimal solution via convex optimization, considering as the objective function the λ-coefficient associated with the maximum possible increment of the load consumption at all the nodes. Different simulation cases are considered in one test feeder, described as follows: (i) the distribution network without penetration of distributed generation; (ii) the distribution network with penetration of distributed generation; and (iii) the distribution grid with capacitive compensation. Numerical results in the test system demonstrated the effectiveness of the proposed SOCP approximation to determine the λ-coefficient. In addition, the proposed approximation is compared with nonlinear tools available in the literature. All the simulations are carried out in the MATLAB software with the CVX package and the Gurobi solver.


Introduction
Electrical distribution networks are a sub-component of the power system required for interconnecting end-users with the transmission and sub-transmission systems at the substation points [1]. They can cover thousands of kilometers in medium-voltage levels to provide electricity service to urban and rural areas [2,3]. Due to their extensions, they have higher power losses compared with the transmission system [4]. Furthermore, electrical distribution networks may also be more sensitive to voltage instabilities since they are typically operated in a radial structure [5,6]. In the scientific literature, the voltage stability problems in distribution networks are usually addressed via indices based on the voltage profile and load consumption. The authors in [7] studied the problem of voltage stability in AC distribution networks with large-scale photovoltaic generation developing a voltage stability index as a function of the line parameters, power demands, and voltages calculated with a classical power flow approach. They demonstrated via numerical simulations that depending on the renewable energy penetration on a particular node, the voltage stability index can be improved or worsened presenting guidelines for optimal location of these devices. In [8], a voltage stability index was proposed to identify the nodes close to the verge of voltage collapse. Its proposed voltage stability index value was calculated at each node of the distribution grid using a modified load flow method for voltage stability analysis. The modified load flow method incorporated composite load modeling and variations in load patterns at each node. Additionally, the results in [8] were found satisfactory compared with the other known methods of voltage sensitivity indices. The authors in [9] presented the classical Newton-Raphson method's extension based on the continuation point for voltage stability analysis in distribution networks. Their numerical results showed the efficiency in the 33-nodes test feeder with low-computational effort. Their results were also compared with specialized tools such as PSAT for MATLAB, reaching a minor error when the λ-coefficient is calculated. In [10], an uncertainty quantification approach for phasor measurement units (PMUs) in voltage stability evaluation was described. In [11], the calculation of a voltage stability margin in radial distribution system was presented considering a unique reactive loading index. In [12], a static voltage stability index in AC distribution systems was proposed using the network-load admittance ratio. This index showed a behavior highly linear with a load increase. In [13], a theoretical voltage stability analysis for AC distribution networks was presented. This analysis determined some essential practical points such as static power reserve until the voltage instability phenomenon, static voltage stability limit, and maximum active power. In [14], the impact on voltage stability margin in an islanded microgrid was studied considering photovoltaic generation. It is worth mentioning that in the the case of DC distribution networks, some approaches in relation with determining voltage stability margin have been reported in literature. In [15] the authors have proposed a semidefinite programming model calculation of the voltage stability margin in DC grids with constant power loads considering a monopolar grid connection. Numerical results demonstrated that the SDP approach reaches the global optimum comparing its results with multiple nonlinear solvers available in the general algebraic modeling system (GAMS). The authors of [16] have presented the problem of the λ-coefficient calculation in DC grids using interior point methods available in GAMS in a tutorial style. In reference [17] it has been proposed an approach to calculate the voltage stability index in DC networks by using the determinant of the Jacobian matrix as a function of load increments controlled by a linear function. Numerical results have demonstrated the possibility of reaching the global optimum when compared with nonlinear programming and semidefinite programming approaches. The main advantage lies mainly in its easy implementation at any programming language.
After the revision of the state-of-the-art, we identify a gap in the literature regarding the voltage stability analysis, which is related to the existence of convex mathematical formulation for voltage stability analysis in AC distribution networks in the optimization context. The convex analysis is a part of the exact mathematical optimization that exploits the structure of solution spaces and objective functions to solve complex problems that guarantee the global optimum [18]. For the voltage stability analysis, a second-order cone programming (SOCP) formulation of the nonlinear modified optimal power flow problem in the complex domain is proposed to compute the voltage stability margin (λ-coefficient) in distribution networks. Numerical results in the 33-nodes test feeder demonstrate that the proposed SOCP model reaches the global optimum of the problem compared with multiple nonlinear solvers available in the (GAMS) [16]. The proposed approach from convex optimization results comfortable to implement at any programming language and avoids the classical usage of the continuation method that requires an adequate parameterization to guarantee convergence.
This paper is organized as follows: Section 2 presents the exact nonlinear formulation of the voltage stability margin calculation using a modified version of the optimal power flow problem. Section 3 presents the SOCP formulation to calculate the λ-coefficient in distribution systems, highlighting the main advantage that consists of finding the optimum global. Section 4 provides the main characteristic of the 33-nodes test feeder and the simulation cases. Section 5 presents all the numerical simulations, analysis and discussion. Section 6 gives the main concluding remarks derived from this research and some guidelines for future avenues of research.

Exact Formulation
The problem of the voltage stability margin calculation in distribution networks can be represented with a nonlinear programming model, considering as the objective function the analysis of the λ coefficient [19]. This parameter allows to know the maximum simultaneous increment at all the loads before voltage-collapse [20]. The main complicated constraints of this model corresponds to the power balance equations, which are nonlinear due to the hyperbolic relation between voltages and powers at each node. The mathematical formulation of this problem is presented below.

Objective Function
The objective function of this problem corresponds to the loadability coefficient of the whole distribution grid, which is maximized and can be written mathematically, as follows where z is the objective function value, and λ represents the maximum chargeability coefficient.

Set of Constraints
The set of constraints associated with the problem of the voltage stability margin calculation are described below.
where V i ∈ C is the voltage value at node i; S cg i ∈ C is the apparent power generation at node i by the conventional generator, i.e., the slack node; S ds i ∈ C is the apparent power generation at node i by a distributed source, i.e., distributed generators or capacitor banks; S d i ∈ C is the apparent power consumption at node i; Y ij ∈ C corresponds to the component of the admittance matrix Y ∈ C n×n that relates nodes i and j. Note that N is the set that contains all the nodes which has a cardinality equal to n; and () is the complex conjugate operator.
where γ is a real positive constant that is related to the maximum deviation of the voltage profile at each node when the load increases constantly. For simulation purposes, it can be defined between 0.4 and 0.8 [9]. Due to the nature of the λ-coefficient the following constraint is added to the optimization model The interpretation of the NLP model (1)-(4) is presented below: Equation (1) corresponds to the objective function of the voltage stability problem in distribution networks which deals with computing the maximum loadability factor of the grid. This allows for understanding the margin of stable operation regarding load increments and possible expansion projects. The constraint (2) corresponds to the set of power balance equations regarding apparent power equilibrium based on the second Tellegen's theorem. Expression (3) corresponds to the voltage regulation constraint, which is relaxed to determine the maximum loadability factor of the grid. Finally, Expression (4) presents the positiveness definition of the λ-coefficient. Remark 1. The problem of voltage stability margin calculation in distribution networks is a nonlinear non-convex optimization problem due to Expression (2) that shows the hyperbolic relation between voltage and power at each node, which generates a set of non-affine quadratic constraints. To solve this problem, the power flow equations can be treated with a second-order cone equivalent model that transforms the exact NLP model defined from (1) to (4) into a SOCP equivalent convex reformulation.

SOCP Reformulation
The second-order cone programming approach is a branch of the convex optimization that allows to transform a type of nonlinear programming problems using conic constraints (i.e., representation using norms), that warranties the global optimum solution of the problem due to its convex structure [21,22]. In the case of the power flow problem, in the literature it has been demonstrated that SOCP approximations allow to have the zero-duality gap between the exact NLP model and its conic transformation [23]. Here we present the SOCP transformation of the voltage stability margin calculation model (1)-(4), which allows the global optimal solution, which is not possible with heuristic approaches. Note that the SOCP transformation of this model only requires to rewrite the product between voltage variables in (2) by using a new variable as presented below: with the new variable W ij , the power balance Equation (2) takes the following form: The power balance constraint defined in (6) is now at a set of an affine set of constraints in the complex domain, which implies that it is convex.
To guarantee that the new variable W ij allows recovering the original voltage variables V i , an additional set of constraints is incorporated as described below.
Let us pre-multiply both side of Expression (5) by W ij , which produces: Now, based on the definition (5) we can use an auxiliary variable as W ii = V i , which implies that Equation (7) can be redefined as Note that the right-hand-side of (8) can be represented using its hyperbolic representation as follows: which implies that if we substitute (10) in (8), then, we reach the following result Observe that Expression (10) can be shown using norms as presented below: (11) is still non-convex due to the equality structure; however, in power systems as demonstrated in [18], this can be relaxed by changing the equality symbol into the lower equal symbol as follows

Remark 3. Expression
where now it is a conic convex constraint.
To complete the SOCP representation of the voltage stability margin calculation, we rewrite (3) in its standard form as that can be rewritten by using W ii as follows Note that (14) is a box-type constraint which is convex. (1), (4), (6) and (12). The main advantage of the proposed SOCP model is the opportunity to guarantee the global optimum due to the convex structure of the solution space [18].

Test Systems and Simulation Scenarios
To validate the effectiveness and robustness of the proposed SOCP programming model to determine the voltage stability margin in AC medium voltage distribution test feeders with global optimum capabilities, we employ two classical AC medium-voltage networks composed of 33 and 69 nodes, respectively. The information of these test feeders and the simulation cases is presented below.

33-Nodes Test Feeder
The validation of the proposed SOCP model is made in a classical radial distribution test feeder composed of 33-nodes and 32 lines operated at a nominal voltage of 12.66 kV. The grid configuration of this test feeder is depicted in Figure 1, and its branches and load information are reported in Table 1.   In this feeder, the total active power demand is 3715 kW and total reactive power demand is 2300 kVAr. For simulation purposes, we consider 100 kVA and 12.66 kV as power and voltage base respectively.

69-Nodes Test Feeder
The 69-nodes test feeder is widely known in the literature as the Baran & Wu test feeder which is composed of 69 nodes and 68 branches, and operated at the voltage controlled node with 12.66 kV. The electrical configuration of this test feeder is reported in Figure 2 and all the parametric information regarding nodes and branches are presented in Table 2.   28 20 In this feeder the total active power is 3890.7 kW and the total reactive power is 2693.6 kVAr. For simulation purposes, we consider 100 kVA and 12.66 kV as power and voltage base respectively.

Simulation Scenarios
The validation of the proposed voltage stability margin calculation in AC distribution networks is performed in this research based on the following simulation cases.

Case 1 (C1):
The original configuration of the AC distribution network; i.e., without penetration of distributed generation or capacitor banks. Case 2 (C2): The operation of the distribution network considering the location of distributed generators reported in [24] which are operated with unity power factor. Case 3 (C3): The operation of the distribution network considering the location of fixed-step capacitor banks as recommended in [25].
For simulation cases 2 and 3, we consider the following distributed generators and fixed step capacitor banks: For the 33-nodes test feeder in the C2 the distributed generators are included at nodes 14, 24, and 30 with power injections of about 770.9 kW, 1096.9 kW and 1065.8 kW, respectively. On the other hand, for capacitor banks in the C3, we consider two banks of 450 kVAr located at nodes 13 and 24, and a bank of 900 kVAr positioned at node 30. For the 69-nodes test feeder in the C2, the distributed generators are included at nodes 12, 61, and 64 with power injections of about 813.1 kW, 1444.7 kW and 289.6 kW, respectively. For capacitor banks in the C3, we consider two banks of 300 kVAr located at nodes 11 and 18, and a bank of 1200 kVAr positioned at node 61.

Computational Validation
Implementing the SOCP approximation for voltage stability margin determination in distribution networks is carried out with the CVX tool using the MOSEK solver by using the MATLAB software version 2019b in on a desktop computer with an INTEL(R) Core(TM) i7 − 7700 2.8-GHz processor and 16.0 GB of RAM running on a 64-bit version of Microsoft Windows 10 Home.

33-Nodes Test Feeder
For this test feeder we consider two simulation conditions, the first one corresponds to the three cases simulation case, and the second one is associated to the effect that the renewables have in the loadability factor including daily load variations.

Evaluation of the Simulation Cases
In this simulation scenario it is considered the original formulation proposed in Equations (1), (4), (6) and (12) to determine the voltage stability margin for all the nodes considering all the simulation cases. Note that to verify that the solution provided by the SOCP model is optimal, we implement the exact nonlinear model (1)-(4) in the GAMS software with different NLP solvers. Table 3 presents the λ-coefficient for each simulation case. Results in Table 3 demonstrate that the proposed SOCP model finds the global optimum solution of the voltage stability margin determination in the case of the 33-nodes test feeder. Additionally, we can observe that: (i) the injection of active power by distributed generation (see C2) allows to increment the voltage stability margin about 27.92% respect to the base case (i.e., C1); and (ii) the injection of reactive power also allows to increase the voltage stability margin about 12.11%; however, its impact is lower when compared with the active power case. This behavior is explained since the total injection of active power is about 2933.6 kW, while the reactive power is 1800 kVAr. Table 3, neglecting their small differences, is the same reported by the continuation method reported in [9], which confirms the efficiency of our convex approximation in comparison with classical methods [7].  Note that throughout all the loadability cases, the minimum voltage at point of collapse is presented at node 18 with a magnitude of 0.3888 pu. However, we can observe that C2 allows lower voltage profiles in the case of nodes upper than 26. This behavior is expected since the load increments in this scenario are larger than C1 and C3, as presented in Table 3. Figure 4 presents the voltage behavior of the node 18 since, as depicted in Figure 3. This is the node that reaches the minimum voltage value at each simulation case.  Figure 4 confirms that the injection of the active and reactive power in some nodes of the distribution network increases the voltage stability margin in all the nodes of the network, as demonstrated in [26]. This implies that the problem of the optimal location of capacitor banks and/or distributed generators is a subject currently under study, since these devices have substantial impacts regarding voltage stability, power losses, and voltage regulation improvements [27]. Remark 6. The proposed SOCP approximation has allowed us to determine the maximum chargeability coefficient for the 33-nodes test feeder in different operation cases, guaranteeing its global optimal solution. In addition, numerical results demonstrate that the voltage collapse point is far from the nominal operation case for the 33-nodes test feeder. This means that to reach the instability point over all the loads, these must be increased at least 2.4 times, which in practice will not occur since the protection system will operate much before of this case. However, the proposed SOCP is general, useful to identify instability scenarios even if these are very closely to the nominal operation case.

Remark 5. The solution of the voltage stability margin in the 33-nodes test feeder via convex and nonlinear optimization methods reported in
Regarding processing times, we can mention that all the NLP solvers in GAMs and the proposed SOCP approach take less than 10 ms to find the solution of the voltage stability margin problem, which implies that from the computational effort perspective, we can confirm that all the methodologies are efficient; however, in terms of solution quality, the best methodology is the SOCP approach.

Effect of Renewables in the Stability Margin
To present the effect of the renewable energy penetration in the voltage stability index for distribution networks, consider the 33-nodes test feeder as recommended in [28,29] with the penetration of photovoltaic and wind turbine sources. The connection of the generators for each test system is presented as follows: at node 13, it is connected a photovoltaic generator PV 1 and a wind turbine WT 1 with nominal rates of 450 kW and 825 kW, respectively. At node 25, it is connected a second photovoltaic source PV 2 with a nominal power rate of 1500 kW while at node 30, it is connected the second wind power source named WT 2 with a rate capability of 1200 kW. The hourly behavior of these generators is reported in Table 4. Table 4. Behavior of the photovoltaic and wind power sources during a sunny day [28]. In the case of the demand power's daily behavior, the load curve reported in [30] is considered as depicted in Figure 5.  Figure 6 presents the behavior of the voltage stability margin in the 33-nodes test feeder considering different levels of power generation penetration, which is defined from 0 to 100% of the power generation reported in Table 6 in steps of 25%. From Figure 6, we can observe that: In the time periods when the demand is low (see the period between 0 h and 8 h in Figure 5), the λ-coefficient is higher taking values upper than 6 with a maximum between 18 and 22 in the period between 3 and 5 h. This behavior is explained by the fact that the chargeability coefficient is a factor that multiplies the load at each operating condition. This implies that for low demands, this will increase until the point of the voltage-collapse. For the time periods with higher demand (upper than 10 h in Figure 5) it is possible to see that the λ-coefficient oscillates between 2 and 4, and we also observe that depending on the level of distributed generation penetration, this increases respect to the base case (0% of renewable generation availability), i.e., the renewable generation has positive effects regarding voltage stability margin since for all the penetration cases the loadability coefficient is enlarged. Note that the minimum λ-coefficient is presented at 19.5 h. At this point when the renewable energy penetration is 0%, the λ-coefficient is 2.407 (see Table 3 for the C1), and when the renewable energy penetration is 100%, this factor is 2.883. Note that this modest increment is because the renewable generation based on photovoltaic is zero at this time and the wind power is also in low values; nevertheless, this increment contributes to the voltage stability margin improvement.

69-Nodes Test Feeder
In this section it is computed the voltage stability margin for the 69-nodes test feeder considering all simulation cases 1 to 3. Table 5 presents the numerical results reached by the proposed SOCP and the NLP solvers available in the GAMS optimization package. Numerical results in Table 5 confirms that the proposed SOCP reformulation can deal with the global optimum of the studied optimization problem which matches with the solutions reported by the GAMS solvers. Additionally, the numerical performance of the cases C2 and C3 show that the injection of the active and reactive power allows enlarging the voltage stability margin in 32.84% and 12.03% respectively. These are expected results, since these local injections reduce the total consumption observed in the distribution network, which is directly connected to the rate of power increments allowed on them and measured by the λ-coefficient.
Regarding processing times as reported for the 33-nodes test feeder, the proposed SOCP reformulation and the GAMS solvers take less than 10 ms to find the optimal solution of the problem. This low-computational effort shows the efficiency of all the studied methods to solve the voltage stability margin in radial distribution networks with considerable number of nodes, i.e., 33 and 69, respectively.
When the voltage profiles are analyzed in the 69-nodes test feeder, in the case C1 the minimum voltage is experienced at node 65 with 0.4700 pu; in addition, for the cases C2 and C3 the minimum values result at the same node (i.e., node 65) with magnitudes of 0.4740 pu and 0.4660 pu respectively. The interpretation of these values is as follows: (i) these correspond to the maximum voltage drop on the distribution network under extreme load conditions; however, these values will never occur under real operation, since protection schemes regarding low-voltages will disconnect the distribution network to protect all the components; and (ii) the λ-coefficients reported in Table 5 can be understood as grid efficiency indicators, since the farther from zero they are, the more robust the distribution network is unser load increments, due to that the λ-coefficient measures the extreme load condition where the distribution grid will present a voltage collapse.

Additional Results
The proposed SOCP model corresponds to an adequate optimization methodology to deal with finding the global optimum, ensuring convergence with low-computational effort. Even if the GAMS package can deal with the global optimum there are two main problems when using it for large-scale nonlinear optimization as follows: (i) it is not possible to ensure mathematically the global optimum due to the non-convexities of the solution space given by the power balance equations, and (ii) it requires a commercial license to deal problems with thousands of variables, while the proposed SOCP model is suitable for being implemented in free software such as Phyton with available tools for solving convex optimization problems [31].
In addition, here, we present a possible application where the usage of the SOCP allows to reach the global optimum of mixed-integer programming problems which is not possible with mixed-integer NLP models solvable in the GAMS software. For doing so, consider that in the 33-nodes test feeder it is possible to locate three distributed generators with capabilities of generation from 0 kW to 1200 kW considering a penetration equivalent to 60% of the total demand under normal operative conditions. Table 6 presents the solution of the optimal sizing and location of distributed generators in AC distribution networks, to improve the voltage stability margin with different GAMS solvers as well as using a branch MI-SOCP solver via branch & bound methods [32]. Once the numerical results reported in Table 6 are obtained, we can mention that: All the MINLP solvers available in GAMS, different optimal solutions are reached, which oscillate between 3.3010 and 3.3184. In addition, these solvers identify different nodes and sizes of the distributed generators being node 32 the most recurrent location in all the solutions. The MI-SOCP strategy allows finding the global optimal solution of this problem with a λ-coefficient of 3.3187. We can ensure that this is indeed the global optimum since an exhaustive search has been implemented to evaluate all the possible combinations. After 2 h of simulations, the combination of nodes 15, 18, and 32 is the best possible scheme to enlarge the voltage stability margin in the 33-nodes test feeder. Regarding processing times it is worth mentioning that the GAMS solvers and the proposed MI-SOCP have faster performance to obtain optimal solutions since running time of all the simulation cases was lower than 10 s. This is considered in literature a negligible time in planning purposes, since physical installations of these distributed generators can take several weeks or months.
The most important fact from the results reported in Table 6 is that the optimization approach that uses the SOCP formulation to solve the continuous optimization part of the optimal location of distributed generators in distribution networks, allow to find the global optimum, if the branch & bound method is combined with it as demonstrated in [33]. This result confirms the importance of using convex optimization tools when possible in distribution network analysis.

Conclusions and Future Works
The problem of determining the voltage stability margin in AC distribution networks has been addressed in this paper under the context of the exact mathematical optimization, proposing a second-order cone programming reformulation for the exact nonlinear model. Numerical validations in the 33-and 69-nodes test feeders have demonstrated that the proposed SOCP approach converges to the global optimum compared with different NLP solvers available in the GAMS optimization package. The injection of active and reactive power in some nodes allowed to demonstrate that the λ-coefficient increases respect to the base case, which implies that these power injections had positive effects on the network's voltage stability margin. The 24-h operative scenario confirmed these results with renewable power source and load variation. Depending on the percentage of renewable, the loadability coefficient always increase respect to the base case.
The hybridization of the proposed focus of SOCP with a branch & bound approach demonstrates the efficiency of using convex optimization in complex optimization problems that involves binary variables, since it (i.e., MI-SOCP) has demonstrated global optimization capabilities, that are not possible with MINLP solvers available in GAMS. These latter have stuck in local optimal solutions due to the non-convexities of the continuous optimization part of the problem. This result confirms the effectiveness and robustness of using convex optimization in AC distribution analysis when possible.
For future work, it will be possible to perform: (i) extend the proposed SOCP programming formulation to an MI-SOCP (mixed-integer-SOCP) to optimal locating and sizing distributed generators and capacitor banks for voltage stability improvement in AC distribution networks; (ii) apply the proposed analysis to the case of the operation of distribution networks with battery energy storage systems; and (iii) make a comparison between voltage stability indexes reported in the literature and the proposed SOCP model to identify if this voltage stability analysis is compatible.