**5. Numerical Experiments**

In this section, several numerical experiments are carried out, and their results are analyzed. The studied Newton-SIP PF techniques are compared with NR in several test systems. All simulations have been done using Matpower v7.0 [38]. The studied systems have been taken from MATPOWER's database [39–41].

In all simulations, **g**(*x*)∞ ≤ 10−<sup>6</sup> has been taken as a convergence criterion, and a flat start has been considered for initializing the PF analysis. The reported execution times have been obtained under Windows 10 on a 3.4 GHz Intel Core i7-8750H CPU 2.2 GHz personal laptop (16.00 GB RAM) and calculated as the average value of 1000 simulations.

## *5.1. Well-Conditioned Cases*

Firstly, we have analyzed the performance of studied PF techniques in several wellconditioned systems, which range from 30 to 3120 buses. Figure 4 shows the obtained results of these systems. These results have been obtained for the maximum convergence rate of the studied Newton-SIP PF techniques, i.e., when the *s*-parameters are 1.

As expected, SIP1-J0 and SIP3-J0 are the fastest methods, which is strongly linked with the number of factorizations required (one should note that the LU decomposition is the heaviest part of any PF calculation [13]). Among all the studied Newton-SIP techniques, only SIP2-J and SIP4-J are occasionally slower than NR. These results may look not coherent with the analysis performed in Section 4; however, Figure 4c provides a clear explanation about this issue. In this figure, it can be appreciated that these two techniques frequently required more factorizations than NR, which is reflected in a higher computational burden and therefore less competitive execution times. Regarding the total iterations required to attain the solution, the results are expected since the highest convergence rate has the lowest

number of total iterations required for achieving the solution. Regarding those algorithms with linear convergence (Equations (4)–(7)), the following relations normally hold:

$$
\dot{\mathbf{t}}\_{\text{SIP4}-\text{J}\_0} \prec \dot{\mathbf{t}}\_{\text{SIP2}-\text{J}\_0} \prec \dot{\mathbf{t}}\_{\text{SIP3}-\text{J}\_0} \prec \dot{\mathbf{t}} \mathbf{t}\_{\text{SIP1}-\text{J}\_0} \tag{46}
$$

where *it*a indicates the total number of iterations of the method "a". One should note that relations (46) are coherent with the theoretical analysis performed in Section 3 (see Figure 2); nevertheless, there are exceptions such as case2746wop.

**Figure 4.** Comparison of the results obtained in well-conditioned systems. (**a**) Execution time (ms), (**b**) total iterations, and (**c**) total factorizations.

Now, let us consider the influence of the loading level. To do that, the active and reactive powers injected at PQ buses along the active powers injected at PV buses have been progressively increased in steps of *λ* = 0.0001 pu until all studied techniques have diverged. Figure 5 is analogous to Figure 4 for limit-loading cases. In this case, all techniques diverged for the same loading level; hence, all of them have been tested for the same conditions.

**Figure 5.** Comparison of the results obtained in well-conditioned systems in a limit load scenario. (**a**) Execution time (ms), (**b**) total iterations number and (**c**) total factorizations.

While similar conclusions can be extracted for SIP1-J, SIP2-J, SIP3-J, and SIP4-J, linear algorithms are not competitive in this scenario due to the huge amount of iterations required to reach the solution. Hence, although they frequently employed very few factorizations, their execution times are not competitive at all. In addition, relations (46) are not held in this situation for SIP2-J0 and SIP3-J0, since it can be observed in Figure 5 that clearly *it*SIP3−J0 < *it*SIP2−J0 .

In order to overcome the important drawbacks shown by linear techniques when the loading level is high, *s*-parameters can be occasionally set greater than one. In this case, these parameters have an accelerating effect on the convergence performance of the studied Newton-SIP techniques. Figure 6 shows the total number of iterations of different linear Newton-SIP techniques when their *s*-parameters are fixed greater than one and therefore achieve the fastest convergence. Under these settings, the results are compared with those results depicted in Figure 5. In this case, *s*-parameters have been obtained empirically. As can be seen, the total number of iterations can be drastically reduced; however, the criteria for determining the best values of s-parameters do not follow any specific pattern, as shown in Figure 6b. In light of the results obtained, this topic looks strongly case-dependent.

**Figure 6.** (**a**) Comparison of the total number of iterations of different linear Newton-SIP solvers, taking the *s*-parameters equal to 1. (**b**) Considered *s*-parameters.

## *5.2. Ill-Conditioned Cases*

Now, we have considered case3012wp, case3375wp, and case13659pegase, which are available in MATPOWER's database. These cases correspond with snapshots of real cases, which demonstrates that ill-conditioned solvers may appear in real applications and they are even more frequent nowadays [42]. NR fails to solve these systems when a flat start is used, so that they can be categorized as ill-conditioned [13]. In previous simulations, we have taken the *s*-parameters to be equal to one, since the studied methods achieve their maximum convergence rate for these values. However, in ill-conditioned systems, this strategy may lead to divergence. Therefore, it is more suitable to study which values of the *s*-parameters in the considered Newton-SIP methods are reliable. Figures 7–9 show the areas of successful convergence for the studied ill-conditioned cases. Hence, it is considered to have failed since both divergence and convergence lead to inaccurate solutions.

Firstly, it can be easily appreciated that SIP1-J, SIP2-J, SIP3-J, and SIP4-J are typically less reliable than SIP1-J0, SIP2-J0, SIP3-J0, and SIP4-J0, since the latter normally showed wider convergence areas. There are some remarkable cases; for example, SIP1-J did not converge in the case3012wp and case3375wp; on the other hand, SIP1-J0 and SIP3-J0 frequently converged, regardless of the value of parameters. Finally, SIP4-J and SIP4-J0 look very sensitive to the values of *s*-parameters. Their convergence rate precisely explains the superior robustness features of linear methods. In [43], it is said that the methods with high convergence rates normally show narrow Regions of Attraction; in other words, the highest convergence rate has the most sensitivity with respect to the initial guess. This fact can also be appreciated for other PF techniques such as [10,14], which introduce a discrete step size to reduce the convergence rate and obtain robust techniques properly.

**Figure 7.** Convergence (green) and failure areas (red) in the s-parameter space for solving the case3012wp. Bold green areas indicate where the studied technique successfully converged, employing the least number of iterations.

**Figure 8.** Convergence (green) and failure areas (red) in the s-parameter space for solving the case3375wp. Bold green areas indicate where the studied technique successfully converged, employing the least number of iterations.

**Figure 9.** Convergence (green) and failure areas (red) in the s-parameter space for solving the case13659pegase. Bold green areas indicate where the studied technique successfully converged, employing the least number of iterations.

From Figures 7–9, it can also be seen that tuning the s-parameters is a strongly casedependent topic. For example, while the considered techniques performed very similar in case3012wp and case3375wp, the behavior of the techniques became erratic in the case13659pegase. In addition, the least number of iterations is not always achieved for the same s-parameters.

In order to properly compare the studied Newton-SIP techniques in these ill-conditioned cases, let us consider only those cases in which the considered techniques required the least number of iterations for successfully converging. Figure 10 shows the execution time along with the total number of iterations and factorizations employed by the Newton-SIP PF techniques in the studied ill-conditioned systems. As commented, NR failed in these cases.

As in well-conditioned systems, SIP1-J0, SIP2-J0, and SIP3-J0 are the fastest techniques, while SIP1-J, SIP2-J, SIP3-J, and SIP4-J typically reached the solution employing less iterations. It is worth mentioning that SIP3-J is occasionally faster than SIP4-J0 in the case13659pegase. This is because these two methods require the same number of factorizations in this system; however, SIP4-J0 computes more calculations per iteration. SIP1-J only converged in the case13659pegase.

To conclude this analysis, the influence of the initial guess *x*0 on the convergence features of each studied technique has been analyzed. To this end, two cases have been compared. On the one hand, we assume a flat start as in previous simulations, taking the best value of the *s*-parameters for each case. On the other hand, we took the default starter provided in Matpower, which is normally closer to the solution than the flat initialization. Figure 11 shows the number of iterations for these two cases in the case3012wp. As observed, the quality of the initialization directly affects the convergence of the studied methods, normally employing more iterations when a flat start is used. These results are coherent, since the closer the solution is to the start point, the faster the convergence. This same conclusion was attained in other recent papers such as [44].

**Figure 10.** Comparison of the results obtained in ill-conditioned systems. (**a**) Execution time (ms), (**b**) total iterations number, and (**c**) total factorizations. *S*-parameters have been tuned as all studied techniques employed the least number of iterations.

**Figure 11.** Total iterations employed in case3012wp with different solvers considering a flat start and the default starter provided in Matpower.

It is also worth noting how the starter point affects the robustness of the solvers. Particularly, in the case of SIP1-J, the solution is successfully achieves starting from the default point, while this solver failed from the flat starter. This is coherent with the definition of ill-conditioned systems provided in [13], which is strongly related with the quality of the starting point.

#### *5.3. Influence of the R/X Ratio*

It is well-known that the R/X ratio may negatively impact on the convergence features of PF solvers [21]. In this section, we analyze how the studied techniques are affected by this parameter. To this end, we have considered the small-scale case\_ieee30, since high R/X ratios are more frequently in this kind of network. To properly analyze this aspect, we have considered that the branch resistances are multiplied by a real factor *μ*, and the number of iterations for various values of *μ* is compared in Figure 12. In this case, the *s*-parameters were fixed equal to 1 in order to attain the highest convergence rate. As expected, the higher the R/X ratio, the higher the number of iterations employed to converge. This fact is especially remarkable in those solvers with fixed Jacobian, provoking divergence in some cases.

**Figure 12.** Total iterations employed in case\_ieee30 with different solvers and various R/X ratios.

Lastly, we analyze the performance of the developed solvers on a real radial distribution system. To this end, we have considered case18 (18-bus radial distribution system from Grady, Samotyj, and Noyola) and case141 (141-bus radial distribution system from Khodr, Olsina, De Jesus, and Yusta) from Matpower's database. In this case, the results obtained by the studied solvers have been compared with the Forward–Backward sweep algorithm (FBS) [45], which is frequently considered the most conventional solver for radial distribution systems. Figure 13 shows the total number of iterations in the studied radial systems. As observed, the Newton-SIP methods normally outperformed FBS, thus proving their efficacy to handle a wide variety of networks with different features and topologies.

**Figure 13.** Total iterations employed in the various radial distributions test networks.

#### **6. Conclusions and Future Works**

In this paper, the applicability of different Newton-SIP methods for solving the PF problems has been comprehensively studied. Four techniques [32,33] have been considered, and two different iterative schemes have been analyzed. Firstly, we have considered the standard form of Newton-SIP techniques, i.e., those algorithms in which the Jacobian matrix is only updated at the first iteration. Secondly, a fully updated iterative scheme has been considered, in which the Jacobian matrix is updated each iteration as NR.

The convergence characteristics of the considered PF techniques have been studied. In case a fully updated iterative scheme is considered, the techniques can achieve up to the eighth order of convergence. It has been also demonstrated that the highest convergence rate is achieved when the *s*-parameters are set equal to one. On the other hand, the convergence rate of the standard form of Newton-SIP methods is always linear.

The efficiency of the studied techniques has been compared using a well-known efficiency index. The results indicate that the SIP3-J is the most efficient Newton-SIP technique, since it is able to achieve the 4th order of convergence by only factorizing one Jacobian matrix for each iteration. Nevertheless, all the studied techniques showed higher efficiency indices than NR.

Various numerical experiments have been carried out for several well and ill-conditioned systems with different sizes and topologies. The most remarkable conclusion is the Newton-SIP techniques' ability to manage both well and ill-conditioned systems efficiently. They typically outperformed NR in well-conditioned systems. In addition, they are more robust than NR in ill-conditioned cases. These features make Newton-SIP techniques very suitable for widespread industrial applications, as it was pointed out in [31]. Drawbacks showed by the linear Newton-SIP techniques in heavy loading cases can be overcome using their fully updated schemes. Comparing the studied Newton-SIP techniques, the best trade-off between robustness and efficiency is normally obtained with SIP3-J.

However, the performance of Newton-SIP techniques is notably influenced by the values of the involved *s*-parameters. For example, in heavy loading systems, it has been shown that the convergence characteristics can be notably improved by taking advantage of the accelerating effect of *s*-parameters. On the other hand, the robustness properties of the Newton-SIP techniques are strongly affected by the parameters involved. However, the analysis carried out in this work shows that tuning the *s*-parameters is a strong casedependent topic. It also looks very difficult to be tackled, since any common pattern has been observed.

In radial distribution systems, the studied techniques outperformed FBS, thus demonstrating their ability to handle a wide variety of networks. These results along those obtained in large-scale well and ill-conditioned cases manifest the suitability of the Newton-SIP methods to be applied in real industry tools and even for voltage stability analysis and optimization problems. Further results will be obtained in future works to confirm that point.

Consequently, future works should be focused on further tackling alternative schemes for optimally tuning the *s*-parameters in order to ge<sup>t</sup> a good trade-off between efficiency and robustness and avoid the necessity to be initially set by the user. In this sense, optimal conditions should be derived, thus allowing implementing auxiliary routines to properly set those parameters.

**Author Contributions:** Conceptualization, M.T.-V.; Data curation, S.K. and I.B.M.T.; Formal analysis, M.T.-V., S.K., I.B.M.T. and F.J.; Funding acquisition, S.K. and I.B.M.T.; Investigation, M.T.-V., S.K. and F.J.; Methodology, M.T.-V., S.K. and F.J.; Project administration, S.K., I.B.M.T. and F.J.; Resources, F.J.; Software, F.J.; Supervision, S.K., I.B.M.T. and F.J.; Validation, S.K.; Visualization, S.K., I.B.M.T. and F.J.; Writing—original draft, M.T.-V. and S.K.; Writing—review and editing, I.B.M.T. and F.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Taif University Researchers Supporting Project number (TURSP-2020/61), Taif University, Taif, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors would like to acknowledge the financial support received from Taif University Researchers Supporting Project Number (TURSP-2020/61), Taif University, Taif, Saudi Arabia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Consideration of Composite Loads**

In this paper, the loads have been considered as constant power models, which is the approach frequently taken in similar papers. However, constant current, constant impedance, or a polynomial combination of the three models can be easily considered as follows:

$$P\_i^{sch}(V\_i) = P\_{G\_i} - P\_{D\_i}(V\_i) = P\_{G\_i} - P\_{0\_i} \left[ a\_{0\_i} \left(\frac{V\_i}{V\_{0\_i}}\right)^2 + a\_{1\_i} \left(\frac{V\_i}{V\_{0\_i}}\right) + a\_{2\_i} \right] \tag{A1}$$

$$Q\_i^{\text{sch}}(V\_i) = Q\_{\text{G}\_i} - Q\_{D\_i}(V\_i) = Q\_{\text{G}\_i} - Q\_{0\_i} \left[ b\_{0\_i} \left( \frac{V\_i}{V\_{0\_i}} \right)^2 + b\_{1\_i} \left( \frac{V\_i}{V\_{0\_i}} \right) + b\_{2\_i} \right] \tag{A2}$$

where *PGi* and *QGi* are the active and reactive power generation at bus *i*, respectively. *PDi* and *QDi* are the active and reactive power demand at bus *i*, respectively, *<sup>P</sup>*0*i* and *Q*<sup>0</sup>*i* are the active and reactive power consumption at rated voltage *<sup>V</sup>*0*i* at bus *i*, respectively. Parameters *<sup>a</sup>*0*i* , *b*0*i* , *<sup>a</sup>*1*i* , *b*1*i* , and *<sup>a</sup>*2*i* , *b*2 are used to model the constant impedance, constant current, and constant power loads, respectively.

Models (A1) and (A2) can be easily incorporated in the proposed methods by substituting these expressions in the system of Equation (1).
