This article generalizes a recently introduced procedure to solve nonlinear systems of equations, radically departing from the conventional Newton–Raphson scheme. The original nonlinear system is first unfolded into three simpler components: (i) an underdetermined linear system, (ii) a one-to-one nonlinear mapping with explicit inverse, and (iii) an overdetermined linear system. Then, instead of solving such an augmented system at once, a two-step procedure is proposed in which two equation systems are solved at each iteration, one of them involving the same symmetric matrix throughout the solution process. The resulting factored algorithm converges faster than Newton's method and, if carefully implemented, can be computationally competitive for large-scale systems. It can also be advantageous for dealing with infeasible cases.
Solving nonlinear equations efficiently and reliably is of paramount importance in physics, engineering, operational research and many other disciplines in which the need arises to build detailed mathematical models of real-world systems, all of them being nonlinear in nature to a certain extent [1–4]. Moreover, large-scale systems are always sparse, which means that the total number of additive terms in a coupled set of n nonlinear functions in n variables is roughly O(n).
According to John Rice, who coined the term mathematical software in 1969, ‘solving systems of nonlinear equations is perhaps the most difficult problem in all of numerical computations’ [5, p 355]. This surely explains why so many people have devoted so much effort on this problem for so long.
Since the seventeenth century, the reference method for the solution of nonlinear systems is Newton–Raphson's (NR) iterative procedure, which locally approximates the nonlinear functions by their first-order Taylor-series expansions. Its terrific success stems from its proven quadratic convergence (when a sufficiently good initial guess is available), moderate computational expense, provided the Jacobian sparsity is fully exploited when dealing with large systems, and broad applicability to most cases in which the nonlinear functions can be analytically expressed .
For large-scale systems, by far the most time-consuming step lies in the computation of the Jacobian and its LU factors . This has motivated the development of quasi-Newton methods, which make use of approximate Jacobians usually at the expense of more iterations . Well-known examples in this category are the chord method, where the Jacobian is computed only once, and the secant method , which approximates the Jacobian through finite differences (no explicit derivatives are required). Broyden's method is a generalization of the secant method which carries out rank-one updates to the initial Jacobian .
Another category of related methods, denoted as inexact Newton methods, performs approximate computations of Newton steps, frequently in combination with pre-conditioners .
More sophisticated higher order iterative methods also exist, achieving super-quadratic convergence near the solution at the expense of more computational effort. These include Halley's cubic method.
When the initial guess is not close enough to the solution or the functions to be solved exhibit acute non-convexities in the region of interest, the NR method works slowly or may diverge quickly. This has motivated the development of the so-called globally convergent algorithms which, for any initial guess, either converge to a root or fail to do so in a small number of ways [6,11]. Examples of globally convergent algorithms are line search methods, continuation/homotopy methods [12,13], such as Levenberg–Marquardt methods, which circumvent Newton's method failure caused by a singular or near singular Jacobian matrix, and trust region methods.
Other miscellaneous methods, more or less related to Newton's method, are based on polynomial approximations (e.g. Chebyshev's method), solution of ordinary differential equations (e.g. Davidenko's method, a differential form of Newton's method) or decomposition schemes (e.g. Adomian's method ).
While each individual in such a plethora of algorithms may be most appropriate for a particular niche application, the standard NR method still remains the best general-purpose candidate trading off simplicity and reliability, particularly when reasonable initial guesses can be made [15,16].
A feature shared by most existing methods for solving nonlinear equations is that the structure of the original system is kept intact throughout the solution process. In other words, the algorithms are applied without any kind of preliminary transformation, even though it has long been known that by previously rearranging nonlinear systems better convergence can be achieved. According to [3, pp. 174–176], ‘the general idea is that a global nonlinear transformation may create an algebraically equivalent system on which Newton's method does better because the new system is more linear. Unfortunately, there is no general way to apply this idea; its application will be problem-specific’.
This paper explores such an idea through a new promising perspective, based on unfolding the nonlinear system to be solved by identifying distinct nonlinear terms, each deemed a new variable. This leads to an augmented system composed of two sets of linear equations, which are coupled through a one-to-one nonlinear mapping with diagonal Jacobian. The resulting procedure involves two steps at each iteration: (i) solution of a linear system with symmetric coefficient matrix and (ii) computation of a Newton-like step.
The basic idea of factoring the solution of complex nonlinear equations into two simpler problems, linear or nonlinear, was originally introduced in the context of hierarchical state estimation , where the main goal was to geographically decompose large-scale least-squares problems. Later on, it has been successfully applied to nonlinear equations with a very particular network structure, such as those arising in power systems analysis and operation [18,19].
In this paper, the factored solution method, initially derived for overdetermined equation systems, is conceptually re-elaborated from scratch, and generalized, so that it can be applied to efficiently solving broader classes of nonlinear systems of equations.
This paper is structured as follows: in the next section the proposed two-stage solution approach is presented. Then, §§3 and 4 introduce canonical and extended forms of nonlinear functions to which the proposed method can be applied. Next, §5 discusses how, with minor modifications, the proposed method can reach different solution points from the same initial guess. Section 6 is devoted to the analysis of cases which are infeasible in the real domain, where convergence to complex values takes place. Section 7 briefly considers the possibility of safely reaching singular points, whereas §8 closes the paper by showing the behaviour of the proposed method when applied to large-scale test cases.
2. Factored solution of nonlinear systems
Consider a general nonlinear system, written in a compact form as follows: 2.1 where is a specified vector and is the unknown vector.1
By applying the NR iterative scheme, a solution can be obtained from an initial guess, x0, by successively solving 2.2 where subindex k denotes the iteration counter, Δxk=xk+1−xk, Hk is the Jacobian of h(⋅), computed at the current point xk, and Δpk=p−h(xk) is the mismatch or residual vector.
In essence, the new method assumes that suitable auxiliary vectors, y and , with m>n, can be introduced so that the original system (2.1) can be unfolded into the following simpler problems: 2.3 2.4 2.5 where E and C are rectangular matrices of sizes n×m and m×n, respectively, and vector f(⋅) comprises a set of one-to-one nonlinear functions, also known as a diagonal nonlinear mapping  2.6 each with closed-form inverse 2.7
By eliminating vector u the above augmented system can also be written in more compact form: 2.8 and 2.9 Note that (2.8) is a linear underdetermined system, whereas (2.9) is an overdetermined one. Among the infinite solutions to (2.8) only those exactly satisfying (2.9) constitute solutions to the original nonlinear system (2.1). As explained in §3, many systems can be found in practice where such a factorization is possible, the aim being to reduce m as much as possible. In the limit, if a vector y of size m=n can be found, then solutions to the original nonlinear system will be obtained directly (i.e. without iterating) by sequentially solving the unfolded system of equations (2.3)–(2.5). But this case will arise only when the set of n equations comprises just n nonlinear distinct terms. Therefore, the need to iterate in the factored method stems from an excess of nonlinear terms (m>n) rather than from the nonlinear nature of the original problem.
Obviously, when the remaining auxiliary vector y is eliminated the original ‘folded’ system is obtained in factored form: 2.10 This leads to an equivalent expression for the Newton step (2.2) 2.11 where F is the trivial Jacobian of f(⋅). Whether (2.11) offers any computational advantage over (2.2) will mainly depend on the complexity involved in the computation of h(⋅) and its Jacobian H, compared to their much simpler counterparts f(⋅) and F (plus the required matrix products), but the convergence pattern will be exactly the same in both cases.
Yet the augmented system (2.8) and (2.9), arising when the factored form is considered, opens the door to alternative solution schemes which may outperform the ‘vanilla’ NR approach. The scheme developed in the sequel begins by considering the solution of the factored model as a linearly constrained least-squares (LS) problem. Then, an auxiliary least-distance problem is formulated aimed at providing improved linearization points at each iteration.
(a) Solution of the factored model as an optimization problem
Finding a solution to the unfolded problem (2.8) and (2.9) can be shown to be equivalent to solving the following equality-constrained LS problem: 2.12 which reduces to minimizing the associated Lagrangian function 2.13
The first-order optimality conditions (FOOCs) give rise to the following system: 2.14
Given an estimate of the solution point, xk, we can choose yk in such a way that (2.9) is satisfied, i.e. 2.15 Then, linearizing f−1(⋅) around xk allows (2.14) to be written in an incremental form 2.16 where is a diagonal matrix with positive elements.
From the above symmetric system, Δyk can be easily eliminated, yielding: 2.17 or, in more compact form 2.18 If the Jacobian Hk remains non-singular at the solution point, then μ=0 and the above system reduces to (2.11), as expected. Therefore, the Lagrangian-based augmented formulation has the same convergence pattern and converges to the same point as the conventional NR approach.
(b) Auxiliary least-distance problem
The driving idea behind the proposed procedure is to linearize (2.14) around a point which, being closer to the solution than yk, can be obtained in a straightforward manner. The resulting scheme will be advantageous over the NR approach if the extra cost of computing the improved linearization point is offset by the convergence enhancement obtained, if any.
For this purpose, given the latest solution estimate, xk, we consider a simpler auxiliary optimization problem, as follows: 2.19 with yk=f−1(Cxk). The associated Lagrangian function is 2.20 which leads to the following linear FOOCs: 2.21 where . Besides being linear, the unknown x is missing in this simpler problem, which can be solved in two phases. First, λ is computed by solving: 2.22 Then, is simply obtained from 2.23 By definition, is as close as possible to yk while satisfying (2.8).
Next, the FOOCs of the original problem (2.14) are linearized around , hopefully closer to the solution than yk, leading to 2.24
Once again, eliminating yields 2.25 or, in a more compact form 2.26 The above linear system, to be compared with (2.18), provides the next value xk+1, which replaces xk in the next iteration. Moreover, as happened with the solution of (2.18), so long as the Jacobian remains non-singular, μ=0 and xk+1 can be obtained with less computational effort from 2.27
(c) Two-step solution procedure
The two-step algorithm, arising from the proposed factored representation of nonlinear systems, can then be formally stated as follows:
Step 0: Initialization. Initialize the iteration counter (k=0). Let xk be an initial guess and yk=f−1(Cxk).
Step 1: Auxiliary least-distance problem. First, obtain vector λ by solving the system 2.28 and then compute from, 2.29
Step 2: Non-incremental computation of x. Solve for xk+1 the system, 2.30 where is the factored Jacobian computed at . Then update yk+1=f−1(Cxk+1). If ∥xk+1−xk∥ (or, alternatively ∥p−Eyk+1∥) is small enough, then stop. Otherwise set k=k+1 and go back to step 1.
As explained above, step 1 constitutes a linearly constrained least-distance problem, yielding a vector which, being as close as possible to yk, satisfies (2.8). As a feasible solution is approached, and . Note that the sparse symmetric matrix EET needs to be factorized only once, by means of the efficient Cholesky algorithm (in fact, only its upper triangle is needed). Therefore, the computational effort involved in step 1 is very low if Cholesky triangular factors are saved during the first iteration. Moreover, the vector p−Eyk is available from the previous iteration if it is computed to check for convergence.
It is worth noting that step 2 is expressed in a non-incremental form. It directly provides improved values for x with less computational effort than that of a conventional Newton step (2.11), which may well offset the moderate extra burden of step 1.
In [18,19], the two-step procedure based on the factored model was originally developed from a different perspective, related to the solution of overdetermined systems in state estimation (i.e. filtering out noise from a redundant set of measurements).
(d) Factored versus conventional Newton's method
The factored scheme can be more easily compared with NR method if step 2 is written in incremental form. For this purpose, the term is subtracted from both terms of (2.30). Considering that f(yk)=Cxk, this leads to 2.31 Note that the above expression, being equivalent to (2.30), is less convenient from the computational point of view owing to the extra operations involved.
Next, the Taylor's expansion around of the set of nonlinear functions f(y) is rearranged as follows, for y=yk: 2.32 where the remainder comprises the second- and higher order terms of the polynomial expansion 2.33 where is the diagonal matrix containing the jth derivatives of f(y), computed at , and e is a vector of 1's.
(i) If the current point, xk, is so close to the solution that , then will be negligible compared to Δpk. In these cases, the only difference of the factored method with respect to NR lies in the use of the updated Jacobian, , rather than Hk, and the local convergence rate of the factored algorithm will be of the same order as that of Newton's method.
(ii) If the current point, xk, is still so far away from the solution that , then will dominate the right-hand side in (2.34). In these cases, the behaviour of the factored scheme can be totally different from that of Newton's method, which fully ignores when computing xk+1.
(iii) As expected, if step 1 is omitted (i.e. ), the factored scheme reduces to the standard Newton's method.
Let and denote the solution point. Further insights into remark (i), regarding local convergence rate, can be gained by subtracting from both sides of (2.30) and rearranging, which leads for the factored method to 2.35 Manipulating (2.11) in a similar manner yields, for the NR method 2.36 The last two expressions confirm that both the NR and the factored schemes converge quadratically near the solution (third- and higher order terms are negligible in ). Moreover, if is closer than yk to the solution , then the factored method converges faster than Newton's.
In summary, from the above discussion it can be concluded that the local behaviour of the proposed factored algorithm is comparable to Newton's method (quadratic convergence rate), while the global behaviour (shape of basins of attractions) will be in general quite different. The fact that step 1 tries to minimize the distance between and yk surely explains the much wider basins of attractions found in practice for the factored method, but there is no way to predict beforehand the global behaviour of any iterative algorithm in the general case (the reader is referred to §3, ‘Aspects of Convergence Analysis’, in , for a succinct but clarifying discussion on this issue).
3. Application to canonical forms
The factored model and, consequently, the two-stage procedure presented above can be applied in a straightforward manner to a wide range of nonlinear equation systems, which are or can be directly derived from the following canonical forms. Small tutorial examples will be used to illustrate the ideas. None of those examples are intended to assess the potential computational benefits of the factored solution process, but rather to show the improved convergence pattern over the NR method (the reader is referred to §8, where very large equation systems arising in the steady-state solution of electrical power systems are solved by both methods and solution times are compared).
(a) Sums of single-variable nonlinear elementary functions
The simplest canonical form of nonlinear systems arises when they are composed of sums of nonlinear terms, each being an elementary invertible function of just a single variable. Mathematically, each component of h(x) should have the form: where cij is any real number and can be analytically obtained. In this case, for each different nonlinear term, hij(xk), a new variable is added to the auxiliary vector y, keeping in mind that repeated definitions should be avoided. This leads to the desired underdetermined linear system to be solved at the first step: while matrix C of the overdetermined linear system u=Cx is trivially obtained from the definition of u:
Consider the simple nonlinear function represented in figure 1. Introducing two auxiliary variables leads to
Therefore, with the above notation, the involved matrices are
For p=1, there are two solutions (points A and B in figure 1). Table 1 shows, for different starting points, the number of iterations required to converge (convergence tolerance in all examples: ∥Δx∥1<0.00001) by both the NR and the proposed factored procedure. Note that the factored procedure converges much faster than the NR method. The farther x0 from the solution point, the better the behaviour of the proposed procedure compared to Newton's method.
For values of x0<0.75 (the minimum of the function, where the slope changes its sign), the NR method converges to point A. On the other hand, the factored scheme always converges to point B no matter which initial guess is chosen. This issue is discussed in more detail in §5, where it is also explained how to reach other solution points. As expected, the NR method fails for x0=0 (null initial slope), whereas the factored solution does not suffer from this limitation (the updated Jacobian corresponding to the first value is not singular).
Finally, the components of vector λ (and obviously μ) are always null at the solution point, indicating that a feasible solution has been found in all cases (see §6).
In this example, the following periodic function will be considered: for which two auxiliary variables are needed with u1=u2=x. The relevant matrices arising in this case are:
For p=1.4, the two solutions closest to the origin are x=0.6435 and x=0.9273 (points A and B in figure 2).
Table 2 shows, for different starting points, the number of iterations required to converge by both methods and the final solution reached. In this particular example (p=1.4), the NR method sometimes outperforms the factored scheme in terms of number of iterations. However, note that, for starting points far away from the origin, the NR method may converge to ‘remote’ periodic solutions (shown in bold), whereas the factored scheme always converges to point A or B (this issue if further discussed in §5). Other advantages of the factored method will become apparent for infeasible values of p (see the same example, with p=1.5, in §6).
In all cases, the final solution provided by the factored scheme is real (or, more precisely, its imaginary component is well below the convergence tolerance). However, complex intermediate solutions are obtained in some cases, owing to the fact that and return complex values for |y|>1 (this is discussed in §6).
(b) Sums of products of single-variable power functions
Another canonical form to which the factored method can be easily applied arises when the nonlinear equations to be solved are the sum of products of nonlinear terms, each being an arbitrary power of a single variable. Mathematically, each component of h(x) has the form: where cij and qk are arbitrary real numbers.
This case can be trivially handled if the original set of variables, x, is replaced by its log counterpart
Then the auxiliary vector y is defined as in the previous case, avoiding duplicated entries which leads to the desired underdetermined linear system
The second key point is to embed the function in the nonlinear relationship u=f(y), as follows: which leads to the overdetermined linear system to be solved at the second step:
Once vector α is obtained by the factored procedure, the original unknown vector x can be recovered from
Consider the following nonlinear system in two unknowns: which, upon introduction of four y variables is converted into an underdetermined linear system
The nonlinear relationships u=f(y) are which lead to the following Jacobian: and the final overdetermined system in the log variables
Once the problem is solved in the log variables, the original variables are recovered from
For p1=24 and p2=20, there is a known solution at x1=2 and x2=3. Table 3 provides the number of iterations for different starting points. Apart from taking a larger number of iterations, the NR method sometimes converges to the alternative point x1=31.1392 and x2=0.5103 (marked with an ‘A’ in the table), unlike the factored scheme, which seems to converge invariably to the same point. For the farthest starting point tested (last row), the NR method does not converge in 50 iterations. On the other hand, the factored scheme is much less sensitive to the initial guess, since it converges systematically to the same point.
4. Transformation to canonical forms
There are many systems of equations which can be easily transformed into one of the canonical forms discussed above, by strategically introducing extra variables aimed at simplifying the complicating factors. Then, the nonlinear expressions defining the artificial variables are added to the original equation system, leading to an augmented one.
Consider the following nonlinear system: and
By defining a new variable it can be rewritten in augmented canonical form, as follows:
Note that in this case only the original unknowns, x1 and x2, need to be initialized, since x3 can be set according to its definition.
5. Extending the range of reachable solutions
When there are multiple solutions, the NR method converges to a point determined by the basin of attraction in which the initial guess is located. Some heuristics have been developed enabling the NR method to reach alternative solutions in a more or less systematic manner, such as the differential equation approach proposed in , which modifies the equation structure according to the sign of the Jacobian determinant.
In this regard, the behaviour of the proposed factored algorithm is mainly affected by the computable range of f(y) components. When this range does not span the whole real axis, the two-stage procedure may not be able to reach some solution points, irrespective of the starting point chosen. For instance, if y=xq, with q even, then during the solution process, the expression will always return a positive u, provided y is positive (or a complex with positive real component if y is negative). This will be determinant of the computed x values. In those cases, an easy way of allowing the factored procedure to reach alternative solution points is by considering negative ranges for f(y), i.e.
In example 3.1, it was noted that the factored procedure always converges to point B (the positive root in figure 1). However, if the original expression for u1 is replaced by2 then the factored procedure always converges to point A (negative root) in fewer iterations than in example 3.1.
A similar situation arises when the nonlinear system to be solved contains periodic functions, since their inverses will have a limited range related to the period of the original function. For instance, if , then will naturally lie in the range −π/2<u<π/2. If we want to obtain u values in the extended range for any integer q, then we should replace by
In a similar manner, , which naturally lies in the range 0<u<π, should be replaced by to obtain u values in the same extended range.
If, in example 3.2, the expression is replaced by then, the factored method converges, in a similar number of iterations, to the points 6.9267 and 7.2105, at a distance 2π from points A and B in figure 2 (in this simple example, involving a purely periodic function, this is totally an expected result, but the reader is referred to the not so trivial example 5.3).
The fact that the range of f(y) adopted during the iterative process determines the solution point which is finally reached, irrespective of the starting point, is a nice feature worth investigating in the future, since by fully exploring all alternative values u=f(y) one might be able to systematically reach solution points in given regions of interest, without having to test a huge number of starting points. The following examples illustrate this idea.
Consider the extremely nonlinear non-periodic function which is graphically represented in figure 3.
As explained in the previous examples, the following augmented system is to be factored: and
The relevant components and relationships of the factored solution approach are where variables have been introduced
The inverse of the Jacobian is in this case
For p=5, no real solution exists when x<3π/2. In order to obtain real solutions for x>3π/2 (points A, B, C and D in figure 3) as explained before, we have to replace by the more general expression for q>1.
Table 4 shows, for increasing values of q, the points to which the factored solution converges, as well as the number of iterations (starting with x0=qπ). The complex solutions provided for q=0 and 1 (whose real component is very close to the local maximum), indicate that no real solution exists for x<3π/2.
One more example will be worked out to illustrate the capability of the factored procedure to reach different solutions, just by extending the computed range of certain f(y) components.
Consider the 2×2 system proposed by Boggs  and which is known to have only three solutions: (0,1), and (−1,2).
In this case, the following relationships are involved in the factored method:
When u1 and u4 are defined as above, the factored method always converges to (0,1), irrespective of the starting point x0. What is as important, the number of iterations is only slightly affected by x0, no matter how arbitrary it is.
If we extend the computed range of u1 to negative values, by replacing the original definition with then the factored method always converges to , irrespective of the starting point. Finally, when the ranges of both u1 and u4 are extended, as follows: and the third solution point, (−1,2), is reached for arbitrary values of x0.
This kind of controlled convergence cannot be easily implemented in the NR method, characterized by rather complex attraction basins in this particular example. In fact, Newton's method may behave irregularly even when started near the solutions .
Interestingly enough, if the fourth combination of ranges is selected and then the factored procedure converges to the complex point (1.7174+0.2131i,3.9041+0.7320i) showing that no real solution exists in that region. This issue is discussed in §6.
6. Infeasibility and complex solutions
A nonlinear system h(x)=p with real coefficients will have in general an undetermined number of real and complex solutions. It can be easily shown that complex solutions constitute conjugate pairs if and only if f(y*)=[ f(y)]*, which happens for many functions of practical interest. When no real solutions exist we say that vector p is infeasible.
The basin of attraction associated with each solution (real or complex) is characteristic of the iterative method adopted. As discussed in §2d and illustrated by the previous examples, the basins of attraction associated with the factored method are quite different from those of Newton's algorithm. Depending on which basin of attraction the initial guess x0 lies in, the factored method will converge to a real or complex solution.
A summary of the possible outcomes of the factored algorithm is as follows:
— Convergence to a real solution: for feasible cases, starting from a real x0 close enough to a solution, the factored scheme eventually converges to the nearby real solution. This is favoured by the introduction of the first step, aimed at minimizing the distance between successive intermediate points. Note however that, even if both x0 and the solution are real, depending on the domain of existence of the nonlinear elementary functions f(y), the two-stage procedure may give rise to complex intermediate values xk during the iterative process. A real solution (or, more precisely, a solution with negligible imaginary component) can also be reached starting from a complex x0.
— Convergence to a complex solution: for infeasible cases the factored method can only converge to a complex solution, provided the right x0 is chosen. Note that complex solutions cannot be reached if x0 is real and the domains of f(y) and F−1 span the entire real axis. In those cases, selecting a complex initial guess is usually helpful to allow the algorithm to reach a complex solution. For feasible cases, if x0 is real but far from a real solution, the factored iterative procedure might converge to a complex point.
— Inability to converge: as any iterative algorithm, the factored scheme may numerically break down owing to unbounded values provided by the elemental functions f(⋅) or their inverse f−1(⋅), which is the case, for instance, when an asymptote is approached. It can also remain oscillating, either periodically or chaotically, which typically happens when x0 is real, no real solution exists and complex intermediate values cannot pop-up.
The following examples illustrate the behaviour of the factored scheme as well as the nature of the complex solutions obtained when solving infeasible cases.
Let us reconsider example 3.2 with p=1.5, which constitutes an infeasible case for which the NR fails to converge to a meaningful point. The two-stage procedure always converges in a reduced number of iterations to the same complex value (or its conjugate), as shown in table 5 for the same real starting points as in example 3.2 (the same happens with other real starting points).
However, if we prevent f(y) from taking complex values by restricting y to its real domain, then the two-stage procedure remains oscillating around real values, much like the NR method.
We can use this example to see what happens if p is further increased. Table 6 presents the number of iterations and the solution points for increasing values of p (starting from x0=0). The maximum value for which there is a feasible (real) solution is p=1.4142 (critical point). For larger p values, the factored approach converges to complex values with the same real component and increasing absolute values. Eventually, for p=4.204 the two-step procedure breaks down, which means that x0=0 is not in the basin of attraction for this infeasible value of p.
Experimental results show that the complex values to which the factored method converges in infeasible cases are not arbitrary. In this example, evaluating the nonlinear function at the real component of the solution point (0.7854) yields pr=1.4142 (point C in figure 2). This lies almost exactly on the maximum of the function .
In example 3.1, if we set p=−0.2, which is infeasible in the real domain, the factored method converges from nearly arbitrary real starting points to the complex value 0.8090+0.2629i. Taking the real part of the solution yields pr=−0.1011 (point C in figure 1), which is indeed very close to the minimum of the function x4−x3 (pm=−0.1055 at x=0.75).
Let us consider the period function which is infeasible for −2<p<2 and has an infinite number of solutions otherwise (figure 4).
This example is very similar to example 3.2, except for the way the auxiliary variables are defined
In this case, both the elementary nonlinear functions f(y) and the Jacobian functions are well defined all over the real axis (except for the asymptotes of the functions).
Therefore, starting from a real value x0, the factored procedure (and also the NR scheme) will always remain in the real domain, unlike in example 3.2 containing the functions and .
With p=3, starting from x0=1, the two-stage algorithm converges in five iterations to 1.2059, whereas the solution point 0.3649 is reached, also in five iterations, when x0=−1.
For the infeasible value p=1.9, the algorithm remains oscillating for any real starting point. However, with the complex starting value x0=1+i, it converges to the complex solutions shown in table 7 for different infeasible values of p.
Note that the real component of the solution is always the same (0.7854). This value is exactly π/4, a local minimum of the nonlinear function.
The above examples suggest that, when the two-stage method fails to find a real solution, it tries to reach a complex point whose real component is in the neighbourhood of an extreme point of the function, which is indeed a nice feature. This is surely a consequence of the strategy adopted in step 1, which attempts to minimize the distance from new solution points to the previous ones, but this conjecture deserves a closer examination which is out of the scope of this preliminary work.
In any case, it is advisable to adopt complex starting points, particularly when infeasibility is suspected, in order to facilitate complex solutions being smoothly reached.
7. Critical points
Feasible (real) and infeasible (complex) solution points can be obtained as long as the Jacobian matrix () remains non-singular throughout the iterative process. This applies to both the NR and the factored methods. Even though critical points are theoretically associated with singular Jacobians, in practice they can also be safely approached so long as the condition number of the Jacobian is acceptable for the available computing precision and convergence threshold adopted. However, the convergence rate tends to significantly slow down near critical points.
It is worth noting that, in the surroundings of a critical point, it may be preferable to deal with the augmented equation (2.26) rather than solving the more compact one (2.27), which assumes that μ=0.
Let us consider again the periodic nonlinear function of example 6.2: which has an infinite number of critical points for p=2. Table 8 provides the number of iterations and the solution points reached by both the NR and factored methods for different starting points (as in previous examples the convergence threshold is ∥Δx∥1<0.00001). As can be seen, the factored method is much more robust against the starting point choice, and converges always to the same critical point (π/4).
In the neighbourhood of the critical point, the number of iterations is significantly reduced, but the factored procedure still performs much better than the NR method. For the sake of comparison, Table 8 also shows the values corresponding to p=2.1.
In addition, numerical problems may arise during the factored solution procedure if, by a numerical coincidence, any element of the diagonal matrix F−1 becomes null or undefined. This risk can be easily circumvented by preventing diagonal elements of F−1 from being smaller than a certain threshold or abnormally large.
8. Large-scale test cases
This section provides test results corresponding to the nonlinear systems arising in the steady-state analysis of electrical power systems, known in the specialized literature as the ‘power flow’ problem . The results presented herein are representative of many other engineering problems characterized by very sparse equations with network or graph-based structure, in which state variables associated with nodes are linked directly only to their neighbours in the graph .
In the power flow problem, the active and reactive power injections specified at each bus i, and , respectively, are linearly related with unknown powers flowing through branches linking i with its neighbours j∈i, as follows: 8.1 In turn, branch power flows are nonlinearly related to nodal state variables (voltage magnitudes, V i, and phase angles θi): 8.2 where gij, bij and bsh represent branch parameters . For a network with N buses the phase angle of a given node is arbitrarily taken as a reference (θs=0), while the active power injection at this same node (termed slack bus) is not specified to account for unknown network losses. Replacing the pair (8.2) into (8.1) leads therefore to a set of n=2N−1 nonlinear equations in the same number of unknowns, which are customarily solved by means of NR method. In practice, the voltage magnitude rather than the reactive power is specified at generator buses, accordingly reducing the size of the equation system. Using the so-called flat-start profile (V i=1, θi=0), the nonlinear power flow equations are solved typically in three to five iterations by Newton's method, depending on the network loading level and ill-conditioning.
As explained in , the natural choice for the intermediate vector y in this particular application is 8.3 where for every node i and, for each branch connecting nodes i and j, the following pair of variables is adopted: 8.4 and 8.5 Vector y then comprises m=2b+N variables (b being the number of branches). It can be easily checked that, when expressed in terms of y, the power flow equations become linear (the reader is referred to  for the details).
Table 9 compares the number of iterations and solution times (milliseconds) for both the proposed factored procedure and Newton's method, when applied to several benchmark systems ranging from 30 to over 3000 nodes (i.e. some 6000 nonlinear equations) available in (http://www.pserc.cornell.edu/matpower/). The convergence threshold adopted is .
Solution times have been obtained on an Intel Core i7 processor (2.8 GHz, 8 GB of RAM) by running a prototype Matlab implementation. Every effort has been made to optimize the code in both cases, by exploiting available sparsity features and avoiding duplicated operations (e.g. the Jacobian Hk in Newton's method is obtained almost free in this particular application, as a by-product of the computations involved in the mismatch vector Δpk).
As can be seen, the factored method always converges and takes one iteration less than NR in most cases. In scenarios representing highly loaded conditions (i.e. more nonlinearity), where the customary flat start is not so close to the solution point, the factored algorithm tends to save two or more iterations with respect to the NR method .
Further insights into the local convergence rate of the factored method can be gained from figure 5, where the evolution of when solving the 300-bus system from flat start is represented, both for Newton's and factored methods.
In this paper, a new iterative procedure to reliably solve nonlinear systems of equations is presented and illustrated with several examples. The basic idea is to unfold the original system with the help of auxiliary variables so that elementary nonlinearities with explicit inverse can be individually handled. At each iteration, the algorithm first solves a trivial least-distance problem; then a Newton-like step is computed. For realistic problems, the overall computational effort required by the two steps should be of the same order or even lower than that of the standard NR scheme, provided suitable sparsity techniques are adopted.
Practical experience so far, including large sets of nonlinear equations arising in power systems problems, shows that the proposed method clearly outperforms the NR method near the solution point, while its behaviour, far away from NR typical basins of attraction, is less erratic and more controllable. Future efforts should be directed to more thoroughly investigating its behaviour in terms of global convergence and capability to deal with infeasible cases.
Applying the factored scheme to the development of globally convergent and/or higher order iterative methods, as well as other related problems such as nonlinear programming and the solution of nonlinear ordinary differential equations, also seems to be a promising venue of research.
↵1 Mathematicians usually write nonlinear equations as h(x)=0, the problem being that of finding the roots of the nonlinear function. In engineering and other disciplines dealing with physical systems, however, it is usually preferable to explicitly use a vector of specified parameters, p, since the evolution of x (output) as a function of p (input) is of most interest. When no real solutions exist we say that p is infeasible.
↵2 When using Matlab, the function nthroot should be used to obtain real qth roots of negative numbers, with q odd. Alternatively, the use of the equivalent fractional power will provide complex roots, if needed.
- Received March 22, 2014.
- Accepted June 3, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.