*2.1. Background*

The PF problem in polar coordinates can be established as a set of *n* nonlinear equations given by [34]:

$$\mathbf{g}(\mathbf{x}) = \begin{cases} \begin{array}{l} P\_i^{\text{sch}} - \sum |V\_i| \left| V\_j \right| \left| Y\_{ij} \right| \cos \left( \theta\_{ij} - \delta\_i + \delta\_j \right) = 0 \\\ Q\_i^{\text{sch}} - \sum |V\_i| \left| V\_j \right| \left| Y\_{ij} \right| \sin \left( \theta\_{ij} - \delta\_i + \delta\_j \right) = 0 \end{array} \tag{1}$$

where *Psc<sup>h</sup> i* and *Qsc<sup>h</sup> i* are the scheduled active and reactive power at *ith* bus, respectively; *Vi*∠*δi* is the complex voltage at the *ith* bus; *Yij*<sup>∠</sup>*θij* is the *ijth* element of the admittance matrix; and *x* ∈ R*n* is the PF vector of unknowns, which is defined in polar coordinates as follows:

$$\mathbf{x} = \begin{bmatrix} \delta \mathbf{p} \mathbf{v} \,|\, \delta \mathbf{p} \mathbf{Q} \,|\, \mathbf{V} \mathbf{p} \mathbf{Q} \end{bmatrix}^{T} \tag{2}$$

where *δ***PV** ∈ R*ng* is the vector of voltage angles at PV buses, *δ***PQ** ∈ R*nl* is the vector of voltage angles at PQ buses, *<sup>V</sup>***PQ** ∈ R*nl* is the vector of voltage magnitudes at PQ buses, and *nl* and *ng* are the total number of PQ and PV buses, respectively.

In the formulation above, only the well-known constant power loads have been considered, nevertheless, other type of consumers could be considered following the formulation described in the Appendix A. Explicit solutions of the system of nonlinear equations (1) cannot be directly obtained. In this regard, iterative methods are undoubtedly the most popular techniques for solving this kind of problem. Among them, the NR method has been the most widely used in PF analysis. The generic *kth* iteration of the NR for solving (1) is given by:

$$\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - \left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1} \mathbf{g}\left(\mathbf{x}^{(k)}\right) \tag{3}$$

where *J* ∈ R*n*×*n* is the Jacobian matrix, which is formed by the first partial derivatives of (1) with respect to (2). It is well known that the NR method has local quadratic convergence.

Alternatively, other robust and high-order Newton-like methods have been proposed to solve the PF and overcome the drawbacks posed by NR (See Section 1). This paper is focused on the family of Newton-SIP methods, which are described in the following subsections.

#### *2.2. Newton-SIP Methods (Type 1)*

In [32], two Newton-SIP methods have been developed. Firstly, we denote SIP1-J0 to the methodology whose generic *kth* iteration for solving the PF is carried out as follows:

$$\begin{cases} \mathbf{y}^{(k)} = (1 - a)\mathbf{x}^{(k)} + a \left( \mathbf{x}^{(k)} - \left[ I \left( \mathbf{x}^{(0)} \right) \right]^{-1} \mathbf{g} \left( \mathbf{x}^{(k)} \right) \right) \\ \qquad \mathbf{x}^{(k+1)} = \mathbf{y}^{(k)} - \left[ I \left( \mathbf{x}^{(0)} \right) \right]^{-1} \mathbf{g} \left( \mathbf{y}^{(k)} \right) \end{cases} \tag{4}$$

where *α*. On the other hand, another Newton-SIP method, namely SIP2-J0 has been proposed in [32]. In this case, the generic *kth* iteration of SIP2-J0 for solving the PF is given by: 

$$\begin{cases} \mathbf{y}^{(k)} = (1 - a)\mathbf{x}^{(k)} + a \left( \mathbf{x}^{(k)} - \left[ I \left( \mathbf{x}^{(0)} \right) \right]^{-1} \mathbf{g} \left( \mathbf{x}^{(k)} \right) \right) \\ \qquad \mathbf{x}^{(k+1)} = \mathbf{y}^{(k)} - \left[ I \left( \mathbf{y}^{(0)} \right) \right]^{-1} \mathbf{g} \left( \mathbf{y}^{(k)} \right) \end{cases} . \tag{5}$$

Both SIP1-J0 and SIP2-J0 are defined by only one *s*-parameter (namely *α*). The main difference between SIP1-J0 and SIP2-J0 lies in the latter requiring two Jacobian evaluations. It is noteworthy that the iterative algorithms defined by (4) and (5) only evaluate the Jacobian at *x*(0) and *y*(0); hence, they are a priori more efficient than NR.

#### *2.3. Newton-SIP Methods (Type 2)*

In [33], two other Newton-SIP algorithms are proposed. Firstly, let us denote SIP3-J0 to that method whose generic *kth* iteration for solving the PF is carried out as follows:

$$\begin{cases} \begin{aligned} \mathbf{z}^{(k)} &= (1-\theta)\mathbf{x}^{(k)} + \theta \left(\mathbf{x}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{x}^{(k)}\right)\right) \\\ \mathbf{y}^{(k)} &= (1-a)\left(\mathbf{x}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{x}^{(k)}\right)\right) + a\left(\mathbf{z}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{z}^{(k)}\right)\right) \\\ &\quad \mathbf{x}^{(k+1)} = \mathbf{y}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right) \end{aligned} \end{cases} \tag{6}$$

where *θ*. Finally, the methodology denoted SIP4-J0 proposed in [33] is carried out at its generic *kth* iteration for solving the PF problem as follows:

$$\begin{cases} \begin{aligned} \mathbf{z}^{(k)} &= (1-\theta)\mathbf{x}^{(k)} + \theta \left(\mathbf{x}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{x}^{(k)}\right)\right) \\\ \mathbf{y}^{(k)} &= (1-\mathbf{a})\left(\mathbf{x}^{(k)} - \left[I\left(\mathbf{x}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{x}^{(k)}\right)\right) + \mathbf{a}\left(\mathbf{z}^{(k)} - \left[I\left(\mathbf{z}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{z}^{(k)}\right)\right) \cdot \ . \end{aligned} \end{cases} \end{cases}$$
 
$$\begin{aligned} \mathbf{x}^{(k+1)} = \mathbf{y}^{(k)} - \left[I\left(\mathbf{y}^{(0)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right) \end{aligned}$$

As in the previous subsection, the main difference between SIP3-J0 and SIP4-J0 lies in the total number of Jacobian evaluations. While the latter requires three Jacobian evaluations, SIP3-J0 only requires one Jacobian evaluation. An important difference between the methodologies proposed in [32] and those developed in [33] is the number of *s*-parameters involved. While SIP1-J0 and SIP2-J0 are defined by only one *s*-parameter (*α*), SIP3-J0 and SIP4-J0 are characterized by a pair of *s*-parameters (*<sup>α</sup>*, *θ*). Finally, all studied Newton-SIP methods only evaluate the Jacobian matrix at the first iteration (in just one or various points); in this paper, we have considered alternative procedures in which the Jacobian matrices are updated each iteration (as in the standard NR). These alternative techniques have been called SIP1-J, SIP2-J, SIP3-J, and SIP4-J for the SIP1-J0, SIP2-J0, SIP3-J0, and SIP4-J0, respectively. With the aim to summarize, Table 1 collects the main characteristics of the studied PF solution techniques.


**Table 1.** Main features of the studied PF solution methods.

*K* → Total number of iterations.

It is also worth commenting that the *s*-parameters were only defined in the range [0,1] in [32,33]. However, we have not limited their analysis to this range. Therefore, we have avoided this definition in order to avoid misleading.

#### **3. Convergence Analysis of Studied Newton-SIP Methods**

In this section, the convergence rate of the studied Newton-SIP methods is derived. In this case, we consider that the Jacobian matrix is updated each iteration, since the convergence features for iterative procedures (4)–(7) can be derived from this analysis. For this study, the Taylor Expansion technique has been used (see [35] for details).

**Theorem 1.** Let **g** be sufficiently differentiable at each point of an open neighborhood D of *r* ∈ R*<sup>n</sup>*, this is a solution of the system **g**(*x*) = **0**. Let us suppose that **g**(*x*) is continuous and nonsingular in *x*. Then, the SIP1-J converges to *r* with the following error function.

$$\mathbf{e}^{(k+1)} = \left(1 - a^2\right) \mathbf{C} \mathbf{e}^{(k)^2} + \left(\left(a^3 - 3a^2 + 2\right) \mathbf{C}\_3 + \left(4a^2 - 2\right) \mathbf{C}\_2^2\right) \mathbf{e}^{(k)^3} + \mathbf{O}\left(\mathbf{e}^{(k)^4}\right) \tag{8}$$

$$\text{where } \mathbf{C}\_{j} = (1/j!) \left[ \mathbf{g}'(r) \right]^{-1} \mathbf{g}^{(j)}(r), j = 2, 3, \dots, \text{ and } \mathbf{e}^{(k)} = \mathbf{x}^{(k)} - r.$$

**Proof**. Taylor expansion of **g**(*x*) and *J*(*x*) about *r* yields:

$$\mathbf{g}\left(\mathbf{x}^{(k)}\right) = J(r)\left[\mathbf{e}^{(k)} + \mathbf{C}\_2 \mathbf{e}^{(k)^2} + \mathbf{C}\_3 \mathbf{e}^{(k)^3}\right] + O\left(\mathbf{e}^{(k)^4}\right) \tag{9}$$

$$J(\mathbf{x}^{(k)}) = J(r) \left[ I + 2\mathbf{C\_2}\mathbf{e}^{(k)} + 3\mathbf{C\_3}\mathbf{e}^{(k)^2} \right] + \mathbf{O}\left(\mathbf{e}^{(k)^3}\right) \tag{10}$$

where *I* ∈ R*n*×*n* is the identity matrix. Now, let us assume that:

$$\left[\boldsymbol{I}\left(\mathbf{x}^{(k)}\right)\right]^{-1} = \left[c\_1\boldsymbol{I} + c\_2\boldsymbol{\mathfrak{e}}^{(k)} + c\_3\boldsymbol{\mathfrak{e}}^{(k)^2}\right]\left[\boldsymbol{I}(\boldsymbol{r})\right]^{-1} + \boldsymbol{O}\left(\boldsymbol{\mathfrak{e}}^{(k)^3}\right) \tag{11}$$

where *c*s∈ R. Considering the following inverse definition:

$$\left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1}\mathbf{g}'\left(\mathbf{x}^{(k)}\right) = \mathbf{J}\left(\mathbf{x}^{(k)}\right)\left[\mathbf{g}'\left(\mathbf{x}^{(k)}\right)\right]^{-1} = \mathbf{I}.\tag{12}$$

By solving the resulting linear system, one can obtain:

$$\left[\boldsymbol{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1} = \left[\boldsymbol{I} - 2\mathbf{C\_2}\boldsymbol{\varepsilon}^{(k)} + \left(4\mathbf{C\_2}^2 - 3\mathbf{C\_3}\right)\boldsymbol{\varepsilon}^{(k)^2}\right] \left[\boldsymbol{J}(\boldsymbol{r})\right]^{-1} + \boldsymbol{O}\left(\boldsymbol{\varepsilon}^{(k)^4}\right). \tag{13}$$

Now, let us define "*e*(*k*) = *x*(*k*) − *<sup>J</sup>x*(*k*)−1**g***x*(*k*)− *r* and *e*(*k*) *y* = *y*(*k*) − *r*; thus,one can obtain:

$$\mathbf{e}\_y^{(k)} = (1 - a)\mathbf{e}^{(k)} + a\widetilde{\mathbf{e}}^{(k)} = (1 - a)\mathbf{e}^{(k)} + a\mathbf{C}\_2\mathbf{e}^{(k)^2} + \mathbf{O}\left(\mathbf{e}^{(k)^3}\right). \tag{14}$$

The Taylor expansion of **g**(*y*) about *r* yields:

$$\mathbf{g}\left(\mathbf{y}^{(k)}\right) = J(r)\left[\mathbf{e}\_{\mathcal{Y}}^{(k)} + \mathbf{C}\_2 \mathbf{e}\_{\mathcal{Y}}^{(k)^2} + \mathbf{C}\_3 \mathbf{e}\_{\mathcal{Y}}^{(k)^3}\right] + O\left(\mathbf{e}\_{\mathcal{Y}}^{(k)^4}\right). \tag{15}$$

Now, we can calculate the error vector at *k* + 1 as follows:

$$\mathbf{e}^{(k+1)} = \left(\mathbf{y}^{(k)} - \left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right)\right) - \mathbf{r}.\tag{16}$$

After some manipulations, one obtains:

$$\begin{aligned} \mathfrak{e}^{(k+1)} &= \left(1 - \mathfrak{a}^2\right) \mathfrak{C}\_2 \mathfrak{e}^{(k)^2} + \left(\left(\mathfrak{a}^3 - 3\mathfrak{a}^2 + 2\right) \mathfrak{C}\_3 + \left(4\mathfrak{a}^2 - 2\right) \mathfrak{C}\_2^2\right) \mathfrak{e}^{(k)^3} + \mathcal{O}\left(\mathfrak{e}^{(k)^4}\right) . \end{aligned} \tag{17}$$
  $\text{The proof is complete. } \Box$ 

**Theorem 2.** Let **g** be sufficiently differentiable at each point of an open neighborhood D of *r* ∈ R*<sup>n</sup>*; this is a solution of the system **g**(*x*) = **0**. Let us suppose that **g**(*x*) is continuous and nonsingular in *x*. Then, the SIP2-J converges to *r* with the following error function.

$$\begin{aligned} \mathbf{e}^{(k+1)} &= \begin{pmatrix} \mathbf{a}^2 - 2\mathbf{a} + 1 \end{pmatrix} \mathbf{C}\_2 \mathbf{e}^{(k)^2} + \left( \left( 6\mathbf{a}^2 - 2\mathbf{a}^3 - 6\mathbf{a} + 2 \right) \mathbf{C}\_3 + \left( 2\mathbf{a}^3 - 8\mathbf{a}^2 + 8\mathbf{a} - 2 \right) \mathbf{C}\_2^2 \right) \mathbf{e}^{(k)^3} \\ &+ \left( \left( 34\mathbf{a}^3 - 7\mathbf{a}^4 - 54\mathbf{a}^2 + 34\mathbf{a} - 7 \right) \mathbf{C}\_2 \mathbf{C}\_3 + \left( 4\mathbf{a}^4 - 22\mathbf{a}^3 + 37\mathbf{a}^2 - 22\mathbf{a} + 4 \right) \mathbf{C}\_2^3 \right) \mathbf{e}^{(k)^4} \\ &+ \mathbf{O} \left( \mathbf{e}^{(k)^5} \right) \end{aligned} \tag{18}$$

**Proof**. Taking *e*(*k*) *y* and the Taylor expansion of **g**(*y*) about *r* from (14) and (15), respectively, let us calculate the Taylor expansion of [*J*(*y*)]−<sup>1</sup> as (12). Thus, one can obtain:

$$\left[\boldsymbol{J}\left(\boldsymbol{y}^{(k)}\right)\right]^{-1} = \left[\boldsymbol{I} - 2\mathsf{C}\_{2}\mathsf{e}\_{\mathsf{Y}}^{(k)} + \left(4\mathsf{C}\_{2}^{2} - 3\mathsf{C}\_{3}\right)\mathsf{e}\_{\mathsf{Y}}^{(k)^{2}}\right] \left[\boldsymbol{J}(\boldsymbol{r})\right]^{-1} + \mathsf{O}\left(\mathsf{e}\_{\mathsf{Y}}^{(k)^{3}}\right). \tag{19}$$

Now, we can calculate the error vector at *k* + 1 as follows:

$$\mathbf{e}^{(k+1)} = \left(\mathbf{y}^{(k)} - \left[\mathbf{J}\left(\mathbf{y}^{(k)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right)\right) - \mathbf{r}.\tag{20}$$

After some manipulations, one obtains:

$$\begin{aligned} \mathbf{e}^{\left(k+1\right)} &= \left(\mathbf{a}^2 - 2\mathbf{a} + 1\right) \mathbf{C}\_2 \mathbf{e}^{\left(k\right)^2} + \left(\left(6\mathbf{a}^2 - 2\mathbf{a}^3 - 6\mathbf{a} + 2\right) \mathbf{C}\_3 + \left(2\mathbf{a}^3 - 8\mathbf{a}^2 + 8\mathbf{a} - 2\right) \mathbf{C}\_2^2\right) \mathbf{e}^{\left(k\right)^3} + \\ &\quad \left(\left(34\mathbf{a}^3 - 7\mathbf{a}^4 - 54\mathbf{a}^2 + 34\mathbf{a} - 7\right) \mathbf{C}\_2 \mathbf{C}\_3 + \left(4\mathbf{a}^4 - 22\mathbf{a}^3 + 37\mathbf{a}^2 - 22\mathbf{a} + 4\right) \mathbf{C}\_2^3\right) \mathbf{e}^{\left(k\right)^4} + \mathbf{O}\left(\mathbf{e}^{\left(k\right)^5}\right). \end{aligned} \tag{21}$$
  $\text{The proof is complete. } \square$ 

**Theorem 3.** Let **g** be sufficiently differentiable at each point of an open neighborhood D of *r*∈ R*<sup>n</sup>*; this is a solution of the system **g**(*x*) = **0**. Let us suppose that **g**(*x*) is continuous and nonsingular in *x*. Then, the SIP3-J converges to *r* with the following error function.

$$\mathbf{e}^{(k+1)} = \left(2 - 2a\theta^2\right)\mathbf{C}\_2^2 \mathbf{e}^{(k)^3} + \left(\left(3 + 2a\theta^3 - 9a\theta^2 + 4a\right)\mathbf{C}\_2\mathbf{C}\_3 + \left(-a^2\theta^4 + 14a\theta^2 - 4a - 5\right)\mathbf{C}\_2^3\right) \mathbf{e}^{(k)^4} + \mathbf{O}\left(\mathbf{e}^{(k)^5}\right) \tag{22}$$

**Proof**. Let us define "*e*(*k*) = *x*(*k*) − *<sup>J</sup>x*(*k*)−1**g***x*(*k*) − *r* and *e*(*k*) *z* = *z*(*k*) − *r*; then, one can obtain:

$$\mathfrak{e}\_{\mathbf{z}}^{(k)} = (1 - \theta)\mathfrak{e}^{(k)} + \theta\breve{\mathfrak{e}}^{(k)} = (1 - \theta)\mathfrak{e}^{(k)} + \theta\mathbf{C\_2}\mathfrak{e}^{(k)^2} + \mathcal{O}\left(\mathfrak{e}^{(k)^3}\right). \tag{23}$$

The Taylor expansion of **g**(*z*) about *r* yields:

$$\mathbf{g}\left(\mathbf{z}^{(k)}\right) = f(\mathbf{r}) \left[ \mathbf{e}\_z^{(k)} + \mathbf{C}\_2 \mathbf{e}\_z^{(k)^2} + \mathbf{C}\_3 \mathbf{e}\_z^{(k)^3} \right] + O\left(\mathbf{e}\_z^{(k)^4}\right). \tag{24}$$

Now, let us define *e*<sup>ˆ</sup>(*k*) = *z*(*k*) − *<sup>J</sup>x*(*k*)−1**g***z*(*k*)− *r* and *e*(*k*) *y* = *y*(*k*) − *r*; then, we can calculate:

$$\begin{split} \mathbf{e}\_{\mathbf{y}}^{(k)} &= (1 - a)\widetilde{\mathbf{e}}^{(k)} + a\mathbf{e}^{(k)} = \left(1 - a\theta^{2}\right)\mathbf{C}\_{2}\mathbf{e}^{(k)} + \left(\left(4a\theta^{2} - 2a\right)\mathbf{C}\_{2}^{2} + \left(a\theta^{3} - 3a\theta^{2} + 2a\right)\mathbf{C}\_{3}\right) \mathbf{e}^{(k)} + \\ & \quad \left(\left(4a - 9a\theta^{2}\right)\mathbf{C}\_{2}^{3} + \left(15a\theta^{2} - 5a\theta^{3} - 7a\right)\mathbf{C}\_{2}\mathbf{C}\_{3}\right) \mathbf{e}^{(k)} + \mathbf{O}\left(\mathbf{e}^{(k)}\right^{5}\right). \end{split} \tag{25}$$

The Taylor expansion of **g**(*y*) about *r* yields:

$$\mathbf{g}\left(\mathbf{y}^{(k)}\right) = J(r)\left[\mathbf{e}\_{\mathbf{y}}^{(k)} + \mathbf{C}\_{2}\mathbf{e}\_{\mathbf{y}}^{(k)^2} + \mathbf{C}\_{3}\mathbf{e}\_{\mathbf{y}}^{(k)^3}\right] + O\left(\mathbf{e}\_{\mathbf{y}}^{(k)^4}\right). \tag{26}$$

Now, we can calculate the error vector at *k* + 1 as follows:

$$\mathbf{e}^{(k+1)} = \left(\mathbf{y}^{(k)} - \left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right)\right) - \mathbf{r}.\tag{27}$$

of

After some manipulations, one can obtain:

$$\mathbf{e}^{(k+1)} = \left(2 - 2a\theta^2\right)\mathbf{C}\_2^2 \mathbf{e}^{(k)^3} + \left(\left(3 + 2a\theta^3 - 9a\theta^2 + 4a\right)\mathbf{C}\_2\mathbf{C}\_3 + \left(-a^2\theta^4 + 14a\theta^2 - 4a - 5\right)\mathbf{C}\_2^3\right) \mathbf{e}^{(k)^4} + \mathcal{O}\left(\mathbf{e}^{(k)^5}\right). \tag{28}$$
 
$$\text{The proof is complete. } \square$$

**Theorem 4.** Let **g** be sufficiently differentiable at each point of an open neighborhood D *r* ∈ R*<sup>n</sup>*; this is a solution of the system **g**(*x*) = 0. Let us suppose that **g**(*x*) is continuousandnonsingularin*x*.Then,theSIP4-Jconvergesto*r*withthefollowingerrorfunction.

$$\mathbf{e}^{(k+1)} = \mathbf{X\_4e^{(k)}} + \mathbf{X\_5e^{(k)}} + \mathbf{X\_6e^{(k)}} + \mathbf{X\_7e^{(k)}} + \mathbf{X\_8e^{(k)}} + \mathbf{O}\left(\mathbf{e^{(k)}}\right) \tag{29}$$

$$\mathbf{e}^{(k+1)} = \mathbf{X\_4e^{(k)}} + \mathbf{X\_5e^{(k)}} + \mathbf{X\_6e^{(k)}} + \mathbf{X\_7e^{(k)}} + \mathbf{X\_8e^{(k)}} + \mathbf{O}\left(\mathbf{e^{(k)}}\right) \tag{20}$$

where:

$$\mathbf{X\_4} = \left(a^2\theta^4 - 4a^2\theta^3 + 4a^2\theta^2 + 2a\theta^2 - 4a\theta\right)\mathbf{C\_2^3} \tag{30}$$

$$\begin{array}{ll} \mathbf{X\_{5}} = \left( 20a^{2}\theta^{4} \quad -4a^{2}\theta^{5} - 36a^{2}\theta^{3} + 28a^{2}\theta^{2} - 4a\theta^{3} - 8a^{2}\theta + 12a\theta^{2} - 12a\theta + 4a \right) \mathbf{C\_{2}^{2}} \mathbf{C\_{3}} \\ \quad + \left( 4a^{2}\theta^{5} - 24a^{2}\theta^{4} + 48a^{2}\theta^{3} - 36a^{2}\theta^{2} + 4a\theta^{3} + 8a^{2}\theta - 16a\theta^{2} + 16a\theta - 4a \right) \mathbf{C\_{2}^{4}} \end{array} \tag{31}$$

*X***6** = <sup>4</sup>*α*2*θ*<sup>6</sup> −24*α*2*θ*<sup>5</sup> + 60*α*2*θ*<sup>4</sup> − 80*α*2*θ*<sup>3</sup> + 60*α*2*θ*<sup>2</sup> − 24*α*2*θ* + <sup>4</sup>*α*<sup>2</sup>*<sup>C</sup>***2***C*2**3** <sup>+</sup>+2*α*3*θ*<sup>6</sup> − 12*α*3*θ*<sup>5</sup> − 8*α*2*θ*<sup>6</sup> + 24*α*3*θ*<sup>4</sup> + 56*α*2*θ*<sup>5</sup> − 16*α*3*θ*<sup>3</sup> − 146*α*2*θ*<sup>4</sup> + 184*α*2*θ*<sup>3</sup> −128*α*2*θ*<sup>2</sup> + 56*α*2*θ* + 6*αθ*<sup>2</sup> − 8*α*<sup>2</sup> − 12*αθ* + <sup>2</sup>*C*3**2***C***<sup>3</sup>** <sup>+</sup><sup>12</sup>*α*3*θ*<sup>5</sup> − 2*α*3*θ*<sup>6</sup> + 4*α*2*θ*<sup>6</sup> − 24*α*3*θ*<sup>4</sup> − 32*α*2*θ*<sup>5</sup> + 16*α*3*θ*<sup>3</sup> + 90*α*2*θ*<sup>4</sup> − 112*α*2*θ*<sup>3</sup> +72*α*2*θ*<sup>2</sup> − 32*α*2*θ* − 6*αθ*<sup>2</sup> + 4*α*<sup>2</sup> + 12*αθ* − <sup>2</sup>*C*5**2** (32) *X***7** = <sup>96</sup>*α*3*θ*<sup>6</sup> −12*α*3*θ*<sup>7</sup> − 288*α*3*θ*<sup>5</sup> + 396*α*3*θ*<sup>4</sup> − 24*α*2*θ*<sup>5</sup> − 240*α*3*θ*<sup>3</sup> + 144*α*2*θ*<sup>4</sup> + 48*α*3*θ*<sup>2</sup> − 288*α*2*θ*<sup>3</sup> +216*α*2*θ*<sup>2</sup> − 12*αθ*<sup>3</sup> − 48*α*2*θ* + 48*αθ*<sup>2</sup> − 48*αθ* + <sup>12</sup>*αC*6**2** <sup>+</sup><sup>24</sup>*α*3*θ*<sup>7</sup> − 180*α*3*θ*<sup>6</sup> + 516*α*3*θ*<sup>5</sup> − 696*α*3*θ*<sup>4</sup> + 48*α*2*θ*<sup>5</sup> + 432*α*3*θ*<sup>3</sup> − 264*α*2*θ*<sup>4</sup> − 96*α*3*θ*<sup>2</sup> +504*α*2*θ*<sup>3</sup> − 384*α*2*θ*<sup>2</sup> + 24*αθ*<sup>3</sup> + 96*α*2*θ* − 84*αθ*<sup>2</sup> + 84*αθ* − <sup>24</sup>*αC*4**2***C***<sup>3</sup>** <sup>+</sup>−12*α*3*θ*<sup>7</sup> + 84*α*3*θ*<sup>6</sup> − 228*α*3*θ*<sup>5</sup> + 300*α*3*θ*<sup>4</sup> − 24*α*2*θ*<sup>5</sup> − 192*α*3*θ*<sup>3</sup> + 120*α*2*θ*<sup>4</sup> + 48*α*3*θ*<sup>2</sup> −216*α*2*θ*<sup>3</sup> + 168*α*2*θ*<sup>2</sup> − 12*αθ*<sup>3</sup> − 48*α*2*θ* + 36*αθ*<sup>2</sup> − 36*αθ* + <sup>12</sup>*αC*2**2***C*2**3** (33) *X***8** = <sup>4</sup>*α*4*θ*<sup>8</sup> − 32*α*4*θ*<sup>7</sup> − 24*α*3*θ*<sup>8</sup> + 96*α*4*θ*<sup>6</sup> + 240*α*3*θ*<sup>7</sup> − 16*α*2*θ*<sup>8</sup> − 128*α*4*θ*<sup>5</sup> − 944*α*3*θ*<sup>6</sup> + 176*α*2*θ*7+ 64*α*4*θ*<sup>4</sup> + 1872*α*3*θ*<sup>5</sup> − 804*α*2*θ*<sup>6</sup> − 2016*α*3*θ*<sup>4</sup> + 1996*α*2*θ*<sup>5</sup> + 1216*α*3*θ*<sup>3</sup> − 2921*α*2*θ*<sup>4</sup> − 408*α*3*θ*2+ 2524*α*2*θ*<sup>3</sup> + 48*α*3*θ* − 1260*α*2*θ*<sup>2</sup> + 368*α*2*θ* + 16*αθ*<sup>2</sup> − 40*α*<sup>2</sup> − 32*θα* + <sup>4</sup>*C*7**2**<sup>+</sup>−7*α*4*θ*<sup>8</sup> + 56*α*4*θ*6+ 72*α*3*θ*<sup>8</sup> − 168*α*4*θ*<sup>6</sup> − 672*α*3*θ*<sup>7</sup> + 56*α*2*θ*<sup>8</sup> + 224*α*4*θ*<sup>5</sup> + 2516*α*3*θ*<sup>6</sup> − 580*α*2*θ*<sup>7</sup> − 112*α*4*θ*<sup>4</sup> − 4872*α*3*θ*5+ 2518*α*2*θ*<sup>6</sup> + 5280*α*3*θ*<sup>4</sup> − 6000*α*2*θ*<sup>5</sup> − 3280*α*3*θ*<sup>3</sup> + 8546*α*2*θ*<sup>4</sup> + 1128*α*3*θ*<sup>2</sup> − 7368*α*2*θ*<sup>3</sup> − 144*α*3*θ*+ 3766*α*2*θ*<sup>2</sup> − 1108*α*2*θ* − 28*αθ*<sup>2</sup> + 128*α*<sup>2</sup> + 56*αθ* − <sup>7</sup>*C*5**2***C***<sup>3</sup>** + <sup>624</sup>*α*3*θ*<sup>7</sup> − 49*α*2*θ*<sup>8</sup> − 2232*α*3*θ*<sup>6</sup> + 476*α*2*θ*7+ 4272*α*3*θ*<sup>5</sup> − 1984*α*2*θ*<sup>6</sup> − 4728*α*3*θ*<sup>4</sup> + 4628*α*2*θ*<sup>5</sup> − 72*α*3*θ*<sup>8</sup> + 3024*α*3*θ*<sup>3</sup> − 6598*α*2*θ*<sup>4</sup> − 1032*α*3*θ*2+ 5876*α*2*θ*<sup>3</sup> + 144*α*3*θ* − 3184*α*2*θ*<sup>2</sup> + 956*α*2*θ* − <sup>121</sup>*α*<sup>2</sup>*C*3**2***C*2**3** + <sup>24</sup>*α*3*θ*<sup>8</sup> − 192*α*3*θ*<sup>7</sup> + 648*α*3*θ*6− 1200*α*3*θ*<sup>5</sup> + 24*α*2*θ*<sup>6</sup> + 1320*α*3*θ*<sup>4</sup> − 144*α*2*θ*<sup>5</sup> − 864*α*3*θ*<sup>3</sup> + 360*α*2*θ*<sup>4</sup> + 312*α*3*θ* − 480*α*2*θ*<sup>3</sup> − 48*α*3*θ*+ 360*α*2*θ*<sup>2</sup> − 144*α*2*θ* + <sup>24</sup>*α*<sup>2</sup>*<sup>C</sup>***2***C*3**3** (34)

**Proof**. In this case, one can calculate the Taylor expansion of [*J*(*y*)]−<sup>1</sup> about *r* as (12). Thus, one can obtain:

$$\left[I\left(\mathbf{z}^{(k)}\right)\right]^{-1} = \left[I - 2\mathbf{C\_2}\mathbf{e\_z^{(k)}} + \left(4\mathbf{C\_2}^2 - 3\mathbf{C\_3}\right)\mathbf{e\_z^{(k)}}^2\right]\left[I(r)\right]^{-1} + O\left(\mathbf{e\_z^{(k)}}^{(k)^3}\right) \tag{35}$$

where *e*(*k*) *z* is already calculated in (23). Now, taking the Taylor expansion of **g**(*z*) about *r* from (24), let us define *e*<sup>ˆ</sup>(*k*) = *z*(*k*) − *<sup>J</sup>z*(*k*)−1**g***z*(*k*) − *r* and *e*(*k*) *y* = *y*(*k*) − *r*; then, we can calculate:

$$\mathfrak{e}\_y^{(k)} = (1 - \mathfrak{a})\widetilde{\mathfrak{e}}^{(k)} + \mathfrak{a}\mathfrak{e}^{(k)} \tag{36}$$

$$\text{where } \hat{\mathbf{e}}^{(k)} = \left(\mathbf{x}^{(k)} - \left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1} \mathbf{g}\left(\mathbf{x}^{(k)}\right)\right) - \text{r. Now, multiplying (31), one obtains: } \mathbf{x}^{(k)}$$

$$\mathbf{e}\_{\mathbf{y}}^{({\hat{k}})} = \left(a\theta^2 - 2a\theta + 1\right)\mathbf{C}\_2\mathbf{e}^{({\hat{k}})^2} + \left(\left(2a\theta^3 - 8a\theta^2 + 8a\theta - 2a\right)\mathbf{C}\_2^2 + \left(6a\theta^2 - 2a\theta^3 - 6a\theta + 2a\right)\mathbf{C}\_3\right)\mathbf{e}^{({\hat{k}})^3} + \mathbf{O}\left(\mathbf{e}^{({\hat{k}})^4}\right). \tag{37}$$

The Taylor expansion of **g**(*y*) and of [*J*(*y*)]−<sup>1</sup> about *r* yields:

$$\mathbf{g}\left(\boldsymbol{y}^{(k)}\right) = \boldsymbol{J}(\boldsymbol{r})\left[\boldsymbol{\varepsilon}\_{\mathcal{Y}}^{(k)} + \mathbf{C}\_{2}\boldsymbol{\varepsilon}\_{\mathcal{Y}}^{(k)^2} + \mathbf{C}\_{3}\boldsymbol{\varepsilon}\_{\mathcal{Y}}^{(k)^3}\right] + O\left(\boldsymbol{\varepsilon}\_{\mathcal{Y}}^{(k)^4}\right) \tag{38}$$

$$\left[I\left(\mathbf{y}^{(k)}\right)\right]^{-1} = \left[I - 2\mathbf{C\_2}\mathbf{e}\_y^{(k)} + \left(4\mathbf{C\_2}^2 - 3\mathbf{C\_3}\right)\mathbf{e}\_y^{(k)^2}\right] \left[I(r)\right]^{-1} + O\left(\mathbf{e}\_y^{(k)^3}\right). \tag{39}$$

At this point, we can calculate the error vector at *k* + 1 as follows:

$$\mathbf{e}^{(k+1)} = \left(\mathbf{y}^{(k)} - \left[\mathbf{J}\left(\mathbf{x}^{(k)}\right)\right]^{-1}\mathbf{g}\left(\mathbf{y}^{(k)}\right)\right) - \mathbf{r}.\tag{40}$$

After some manipulations, one can obtain:

$$\mathbf{e}^{(k+1)} = \mathbf{X\_{4}e^{(k)}}^{4} + \mathbf{X\_{5}e^{(k)}}^{5} + \mathbf{X\_{6}e^{(k)}}^{5} + \mathbf{X\_{7}e^{(k)}}^{5} + \mathbf{X\_{8}e^{(k)}}^{5} + \mathbf{O}\left(\mathbf{e^{(k)}}^{\circ \text{s}}\right) \tag{41}$$

where the *X*s are defined in (30)–(34); hence, the proof is complete. 

At this point, one can easily check that the studied techniques achieve their highest convergence rates when the *s*-parameters are equal to 1. Thus, SIP1-J, SIP2-J, SIP3-J, and SIP4-J can show third, fourth, fourth, and eighth convergence rates, respectively. In this regard, it is important to note that the mentioned techniques present a higher convergence rate compared with NR, which presents quadratic convergence. On the other hand, it is worth mentioning that the studied Newton-SIP methods have the characteristic lowest convergence rate even when *α*, *θ* -= 1. By observing the error functions of these techniques, one can deduce that the order of convergence of SIP1-J and SIP2-J is at least two, while it is three in the case of SIP3-J and four for SIP4-J.

In the case of the standard forms of studied Newton-SIP methods (Equations (4)–(7)), since their convergence rates do not remain constant during the iterative procedure, it is more suitable to study them by successive substitution and Taylor expansion in their respective algorithms (details of this technique can be found in [36]). Figure 1 shows the convergence rates of the algorithms SIP1-J0, SIP2-J0, SIP3-J0, and SIP4-J0 as a function of the iteration number. In this figure, the convergence rate of the conventional NR has been also included for comparison. From this figure, it can be deduced that SIP4-J0 and SIP1-J0 show the highest and the lowest convergence rate, respectively, while both SIP2-J0 and SIP3-J0 have the same convergence order. Anyway, for these techniques, the convergence rate is linear after the first iteration. Therefore, their convergence orders are always less than two, being so overcome by NR.

**Figure 1.** Convergence rates of the algorithms (4)–(7).

To compare the algorithms (4)–(7) and their respective counterparts (in which the Jacobian is updated each iteration), let us refer to the convergence rate of the initial error vector. It means, let us suppose that the error vector evolves as *e*(0) *J*. Thus, as *j* grows, the solution is assumed to be more closely approached. Therefore, it is assumed that an algorithm will converge faster as *j* grows rapidly. Figure 2 plots the value of *j* of the studied Newton-SIP methods at different iterations. From this figure, it can be seen that *j* exponentially grows in the case of SIP1-J, SIP2-J, SIP3-J, and SIP4-J, while it grows linearly with SIP1-J0 SIP2-J0, SIP3-J0, and SIP4-J0. To complete the section, Table 2 summarizes the convergence analysis of the considered techniques and the NR. As commented, the studied techniques achieve their maximum convergence rate when the *s*-parameters are equal to 1.

**Figure 2.** Convergence degree of the error vector *e*(0) at different iteration counters.

**Table 2.** Comparison of the convergence rates of studied PF solution techniques.


#### **4. Comparison of the Efficiency of Different Iterative Algorithms**

In this section, we compare the efficiency of the studied PF solution methods. To do that, let us consider the following efficiency index [37].

$$FEI = p^{\frac{1}{\Box \Box}} \tag{42}$$

where *p* ∈ R<sup>+</sup> is the order of convergence, and *CO* stands for the total computational cost of an iteration. In this sense, the following theorem is generally used to estimate the cost of a LU decomposition [35].

**Theorem 5.** The number of products and quotients required for solving *q* linear systems of equations with the same matrix of coefficients, using LU factorization, is:

$$
\rho(n,q) = \frac{1}{3}n^3 + qn^2 - \frac{1}{3}n.\tag{43}
$$

It is also suitable to consider the computational cost of each function evaluation *o*(**g**) and Jacobian evaluation *o***g**; hence, the total computational cost of an iteration can be estimated as follows: 

$$\mathbf{CO} = o(\mathbf{g}) + o\left(\mathbf{g'}\right) + o(n, q). \tag{44}$$

The value of the index (42) for studied Newton-SIP methods and NR is depicted in Figure 3 for different sizes of the PF state vector (*n*). From this figure, it can be appreciated that Newton-SIP techniques are more efficient than NR. In addition, the following relations hold:

$$FEI\_{\rm NR} < FEI\_{\rm SIP2-J} = FEI\_{\rm SIP4-J} < FEI\_{\rm SIP1-J} < FEI\_{\rm SIP3-J}.\tag{45}$$

**Figure 3.** Efficiency index (42) for different PF solvers.
