Solution Methods of Large Complex-Valued Nonlinear System of Equations

Nonlinear systems of equations in complex plane are frequently encountered in applied mathematics, e.g., power systems, signal processing, control theory, neural networks, and biomedicine, to name a few. The solution of these problems often requires a first- or second-order approximation of nonlinear functions to generate a new step or descent direction to meet the solution iteratively. However, such methods cannot be applied to functions of complex and complex conjugate variables because they are necessarily nonanalytic. To overcome this problem, the Wirtinger calculus allows an expansion of nonlinear functions in its original complex and complex conjugate variables once they are analytic in their argument as a whole. Thus, the goal is to apply this methodology for solving nonlinear systems of equations emerged from applications in the industry. For instances, the complex-valued Jacobian matrix emerged from the power flow analysis model which is solved by Newton-Raphson method can be exactly determined. Similarly, overdetermined Jacobian matrices can be dealt, e.g., through the Gauss-Newton method in complex plane aimed to solve power system state estimation problems. Finally, the factorization method of the aforementioned Jacobian matrices is addressed through the fast Givens transformation algorithm which means the square root-free Givens rotations method in complex plane.


Introduction
This work is a tribute to Steinmetz's contribution [1]. The reasons and motivations are stated throughout the whole document once the numerical solutions for solving power system applications are typically carried out in the real domain. For instance, the power flow analysis and power system state estimation are wellknown tools, among others. It turns out that these solutions are not well suited for modeling voltage and current phasor. To overcome this difficulty, the proposal described in this chapter aims to model the aforementioned applications in a unified system of coordinates, e.g., complex domain. Nonetheless, the solution methods of these problems often require a first-or second-order approximation of the set of power flow equations; such methods cannot be applied to nonlinear functions of complex variables because they are nonanalytic in their arguments. Consequently, for these functions Taylor series expansions do not exist. Hence, for many decades this problem has been solved redefining the nonlinear functions as separate functions of the real and imaginary parts of their complex arguments so that standard methods can be applied. Although not widely known, it is also possible to construct an extended nonlinear function that includes not only the original complex state variables but also their complex conjugates, and then the Wirtinger calculus can be applied [2]. This property lies on the fact that if a function is analytic in the space spanned by ℜ x f g and ℑ x f g in , it is also analytic in the space spanned by x and x* in . In complex analysis of one and several complex variables, Wirtinger operators are partial differential operators of the first order which behave in a very similar manner to the ordinary derivatives with respect to one real variable, when applied to holomorphic functions, non-holomorphic functions, or simply differentiable functions on complex domain. These operators allow the construction of a differential calculus for such functions that is entirely analogous to the ordinary differential calculus for functions of real variables [2,3]. Then, taken into account the Wirtinger calculus, this chapter shows how the Jacobian matrix patterns emerge in complex plane corresponding to the steady-state models of power flow analysis and power system state estimation, respectively.
In this chapter the classical Newton-Raphson and Gauss-Newton methods in complex plane aiming the numerical solution of the power flow analysis and power system state estimation are derived, respectively. Moreover, the factorization methods addressed to deal with the Jacobian matrices emerged from these approaches are included [4].
This chapter is organized as follows. The theoretical foundation which is based on Wirtinger calculus is summed up in Section 2. Section 3 describes two algorithms suggested to factorize Jacobian matrix in complex plane regardless if it is exactly determined or overdetermined. In Section 4, the complex-valued static model solution by using Newton-Raphson method is derived, whereas in Section 5, the Gauss-Newton method developed in complex plane is equally presented. Finally, in Section 6 some conclusions are gathered and stated the next issues to be investigated in the near future.

Complex differentiability
A complex function is defined as (1) are in general complex but may be real-valued in special cases, e.g., squared error cost function J e 2 À Á . The definition of complex differentiability requires that the derivatives defined as the limit be independent of the direction in which Δx approaches 0 in complex plane: This requires that the Cauchy-Riemann equations be satisfied, i.e., These conditions are necessary for f x ð Þ to be complex differentiable. If the partial derivatives of u a, b ð Þ and v a, b ð Þ are continuous on their entire domain, then they are sufficient as well. Therefore, the complex function f x ð Þ is called an analytic or holomorphic function [2]. As an example, let f which under differentiation rule leads to These results show that the Cauchy-Riemann equations hold, and hence f x ð Þ ¼ y ¼ x 2 is a holomorphic function.

CR-Calculus or Wirtinger calculus
Introduced by Wilhelm Wirtinger in 1927 [2], the CR-Calculus, also known as the Wirtinger calculus, provides a way to differentiate nonanalytic functions of complex variables. Specifically, this calculus is applicable to a function f x ð Þ given by ð Þ and v a, b ð Þ have continuous partial derivatives with respect to a and b, yielding Since we have and by setting ∂x * ∂x to zero, it follows that Note that the Cauchy-Riemann conditions for f Á ð Þ to be analytic in x can be expressed compactly using the gradient as ∂f ∂x * ¼ 0, i.e., f Á ð Þ is a function of only x. Similarly, if we take the derivative of f Á ð Þ with respect to x * , that is, By setting ∂x ∂x * to zero, we get Again, the Cauchy-Riemann conditions for f Á ð Þ to be analytic in x * can be expressed compactly using the gradient as ∂f ∂x ¼ 0, i.e., f Á ð Þ is a function only of x * . In other words, the gradient (respectively conjugate gradient) operator acts as a partial derivative with respect to x (respectively to x * ), treating x * (respectively x) as a constant. Formally, we have As be a real function of complex variable which is the squared Euclidean distance to the origin, with as v ¼ 0; clearly the Cauchy-Riemann equations do not hold, and hence Þ¼x * x is not analytic and thus is non-holomorphic function. To overcome this apparent difficult, by applying the CR-Calculus leads to which suggests the geometric interpretation shown in Figure 1. Its analysis allows us to infer that the direction of maximum rate of change of the objective function is given by the conjugate gradient defined in Eq. (13). Observe that its positive direction is referred to a maximization problem (dot arrow), whereas the opposite direction concerns to the cost function minimization. Hereafter, a real-valued or complex-valued function and its argument are provided with a subscript c if it is a function in the complex conjugate coordinates, i.e., x, x * ð Þ. Moreover, when the CR-Calculus is extended to the vector case, it is denoted that the multivariate CR-Calculus and the basic rules for the scalar case remain unchanged.

Solution of the problem:
As well in the real domain, there are two classes of methods for the numerical solution of large linear system of equations in complex plane: • Direct methods: Which produce the exact solution assuming the absence of truncation and round-off errors, by performing a finite number of flops in a finite known number of steps. These methods are usually recommended when most of the entries in the coefficient matrix are nonzero and the dimension of the system is not too large, for instance, the Gaussian elimination, the LU decomposition and QR factorization, to cite a few.
• Iterative methods: This class of methods is beyond of this work. Notice the solution provided by these methods is approximated and the accuracy is imposed by the user. The number of ν accomplished iterations depends on the given precision or convergence criterion. This class of methods is preferred when the majority of the coefficients are equal to zero and the number of unknowns is very large. The methods of Jacobi and Gauss-Seidel are good examples, besides the classical Conjugate gradient method [5].

Three-angle complex rotation algorithm
The first studied algorithm is the three-angle complex rotations (TACR), which is derived in polar coordinates [6]. Nonetheless, the key idea behind QR decomposition is to eliminate the square roots needed for the computation of the cosine and sine which represent a bottleneck in real-time applications. Consequently, this algorithm is devoid of interest for the solution of large linear systems of equations.

Complex-valued fast givens rotations
On the other hand, the fast plane rotation [7] that is derived in complex plane and rectangular coordinates is a square root-and division-free Givens rotations [8]. In this sense, the fast plane rotation which is also referred as the complex-valued fast Givens rotations (CV À FGR) is a very efficient algorithm, aiming a QRdecomposition of matrices once the computations are performed incrementally, i.e., as the data arrives sequentially in time. Thus, it allows us to reduce the overall latency and hardware resources drastically. In the forthcoming contribution, the CV À FGR performance will be compared to the well-known approaches which were successfully applied to the power system state estimation [9][10][11], once they are accordingly converted from real to complex domain. Further proposals in the updated state of the art will be equally considered, e.g., [12] and [13], to cite a few.
The complex fast Givens transformation M is computed using Algorithm 1 such that the second component of M H x is zero and M H DM is a diagonal matrix, as shown below: Notice the superscript H denotes the Hermitian operation, i.e., complex conjugate transpose.
The matrices Q , M, and D are connected to the following equation: In the sequence the QR-decomposition using complex fast Givens transformations is presented as Algorithm 2. Algorithm 2. Sequential fast Givens QRD decomposition.
Step 1: Get alpha and beta using Algorithm 1: Step 2: Update elements based on type if j < n and type=1 then Note that the complex fast Givens QRD does not require any square root operation, and during each incremental QRD-update step, the incoming input data row-vector is stored, i.e., vector new. In the sequence, the input data row-vector elements are zero-out (inner for loop) in order to update upper triangular matrix R.
The new vector is overwritten each time till the QRD-algorithm has exhausted all the input data, i.e., the upper triangular matrix R is entirely updated.

Complex-valued Newton-Raphson method
Aiming the solution of any set of exactly determined equations in complex plane, the vector of unknowns is regularly taken into account in the iterative algorithm as follows: and the residual vector hereafter termed as "mismatches" vector leads to Nonetheless, here the goal is to calculate x c that satisfies where in Eq. (21), Y s is a vector of specified quantities, i.e., constant term; Y e x c ð Þ ¼ J c Δx c is a vector of calculated quantities at each iteration. Consequently, the linearization of Eq. (21) from one step to the sequel leads to where J is the complex-valued Jacobian matrix which the dimension is 2 N À 1 ð ÞÂ2 N À 1 ð Þ. It means that at least one complex-valued state variable have to be specified, i.e., is known.
As a further advantage provided by the Wirtinger calculus [2,3], the Jacobian matrix which emerged in Cartesian coordinates needs lesser algebra task as well as minor implementation effort (encoding) than the former procedure in real domain [14]. Thereby, the Jacobian matrix in expanded form may be represented through four partitions matrix, yielding For instance, Figure 2 displays the pattern for the IEEE-57 bus system in domain and complex plane. This latter is given by Eq. (24).

Jacobian matrix factorization
In order to factorize the Jacobian matrix required in Eq. (23), the recommended procedure is to operate the factorization on the augmented Jacobian matrix, yielding which the dimension is 2n Â 2n þ 1 ð Þ, resulting where T c is an upper triangular matrix of dimension (2n Â 2n) andM x νÀ1 ð Þ c comprises the corresponding rows in the updated rhs vector of dimension (2n Â 1).

Finally, Eq. (23) is solved by performing a back-substitution via
On the other hand, it is recommendable to perform the convergence checking over the infinity norm of two vectors. Firstly, as the former, it occurs over the corrections to be applied to the state variables and simultaneously over the mismatches vector. This latter is included, aiming to be aware against ill-conditioned systems [14], yielding If Eq. (28) is satisfied, stop and print out the results. Otherwise, the state vector is updated as shown below: and the iteration counter is increased followed by the updating of the mismatch vector and the Jacobian matrix factorization. This latter task can be mandatory or not once the Jacobian matrix may be kept constant throughout the iterative process (approximate, instead of full gain) which is a decision very often adopted after the second iteration aiming to lighten the computational burden. Further details can be found in [14], but in the sequence, a small example is presented forwarded of simulations carried out on large systems. However, as any other application in the industry, the power flow model lies in the solution of a system of linear algebraic equations which is summarized thereafter.

Nodal equation
This approach requires the nodal admittance matrix building, e.g., thus the complex nodal power can be expressed as Then, the nodal complex power at bus À k, i.e., S k , is where N is the number of network nodes. Therefore, the unknowns to be determined are the voltages at each node or bus into the system and the general power flow equations that model any type of branch in an electrical network, i.e., transmission lines and phase-and phase-shifting-transformers can be written, yielding and their complex conjugate counterpart are In Eqs. (34)-(37), t km ¼ a km e Àjφ km is the general off-nominal tap transformer model which is composed by an ideal transformer with complex turns ratio t km : 1 in series with its admittance or impedance. Thus, if the corresponding branch is referred to.

Small example
The power flow model in complex plane as detailed in [14] is applied to a small example system in which the diagram is shown in Figure 3, while the corresponding branch parameters and bus data, both in pu (V base ¼ 230 kV; S base ¼ 100 MVA), are presented in Table 1.
The nodal admittance matrix, Y bus , leads to The whole set of intermediary results throughout the iterative process is presented in the sequence.
As in the real domain, the elements of the complex-valued Jacobian matrix remain practically unchanged after the second iteration, which suggest that we may keep them constant thereafter. Moreover, the computation of some entries can be avoided because they are complex conjugates of other entries; it turns out that these elements are PV-nodes.  Interestingly, the numerical values of the state variables corrections, state variables, and mismatches vectors calculated in the complex plane are displayed in the Tables 2-4, respectively.

Performance in larger systems
The algorithm described earlier was encoded in MATLAB by using sparsity technique and column approximate minimum degree (colamd) ordering scheme. The numerical tests were executed by using an Intel® Core™ i5-4200 CPU at 1.60 Hz 2.30 GHz, 6GB of RAM, and 64-bit operating system. A flat start condition is assigned to the state variables in all simulations.
Thereby, Table 5 presents the performance on larger systems of the Newton-Raphson method in complex plane which is highlighted in bold. The corresponding performance is compared with those of the former Newton-Raphson method developed in polar and rectangular coordinates, both in real domain. In all simulations the convergence criterion of 1 Â 10 À12 is assumed.     Table 3.

Unknown vector.
The results presented in the aforementioned table allows us to infer that the Newton-Raphson method in complex plane has very good performance. Except for the SIN-1916 bus system, the time consuming required to achieve the solution is lower than the remainder approaches.
In the next section, the Gauss-Newton method is presented. The goal is to solve overdetermined systems of equations, i.e., the former nonlinear least-squares method in complex plane.

Complex-valued weighted-least-squares method (CV-WLS)
As shown in [3], the complex-valued WLS state estimator minimizes an objective function defined as where z c = z, z * ð Þis a vector of specified complex values of dimension (2m Â 1), Þis a vector of complex-valued state variables of dimension (2n Â 1), h c (x c ) is a vector of nonlinear functions of dimension (2m Â 1) that maps z c to x c , ω c is a vector of random errors in complex plane which dimension is (2m Â 1), and ð Þ ∂x * be Jacobian submatrices of dimension (m Â n). They are obtained through the Wirtinger partial derivatives with respect to the complex and the complex conjugate state variables using the chain differentiation rule. Now, let us define the Jacobian matrix as In the important special case given by Eq. (39) where J x c ð Þ is a real-valued function of complex variables, the following property holds Here, T c is an upper triangular matrix of dimension (2n Â 2n), and Δz c comprises the corresponding rows in the updated rhs vector of dimension (2n Â 1). Finally, Eq. (55) is solved by performing a back-substitution via

Small example
The one-line diagram of a 2-bus system is depicted in Figure 5. The system is provided with two PMU measurements that meter the nodal voltage magnitudes and phase angles and two real and reactive power flow measurements, which are identified by means of black bullets and red triangles, respectively. In Table 6, the transmission line parameters are given in pu.
From Figure 5, the complex power injections at sending and receiving end, respectively, are written yielding Applying the Wirtinger calculus to Eqs. (59) and (60) leads to the Jacobian matrix given by Eq. (47) at each iteration as shown in the sequence. Here, the complex power injection measurement at each end, i.e., S 1 and S 2 , is equal to the corresponding power flow measurement, i.e., S 12 and S 21 , respectively.  As observed in the above expressions of J c , the sub-matrix J h is more sparse than the sub-matrix J d h , which affects the sparsity structure of the complex-valued Jacobian matrix as shown earlier in Figure 4. The state variables at each iteration are provided in Table 7, while the estimated measured quantities are given in Table 8.  The residual vector and the chi-squared index throughout the iterations are presented in Table 9. Notice that the estimated values obtained in the -domain are exactly equal to those extracted from the power flow report in the -domain. Furthermore, the Gauss-Newton iterative algorithm converges in three iterations in both vector spaces, i.e., real and complex domain.

Performance in larger systems
In the sequel, Table 10 presents the performance on larger systems of the Gauss-Newton method in complex plane which is highlighted in bold. In all simulations the convergence criterion of 1 Â 10 À3 is assumed.
The comparative analysis of the results presented in the aforementioned table allows us to infer that the Newton-Raphson method in complex plane has very good performance. Except for the SIN-1916 bus system, the time consuming required to achieve the solution is lower than the remainder approaches.