Branch and bus data.
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.
- large nonlinear system of equation solution in complex plane
- complex-valued Newton-Raphson and gauss-Newton iterative algorithms
- Cartesian coordinates
This work is a tribute to Steinmetz’s contribution . 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 well-known 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 . This property lies on the fact that if a function is analytic in the space spanned by and in , it is also analytic in the space spanned by and * 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, 4]. 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 .
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.
2. Theoretical foundation
2.1 Complex differentiability
A complex function is defined as
where and , are real functions, , : . Functions like Eq. (1) are in general complex but may be real-valued in special cases, e.g., squared error cost function . The definition of complex differentiability requires that the derivatives defined as the limit be independent of the direction in which approaches in complex plane:
This requires that the Cauchy-Riemann equations be satisfied, i.e.,
These conditions are necessary for to be complex differentiable. If the partial derivatives of and are continuous on their entire domain, then they are sufficient as well. Therefore, the complex function is called an analytic or holomorphic function . As an example, let be a complex function with . Then,
which under differentiation rule leads to
These results show that the Cauchy-Riemann equations hold, and hence is a holomorphic function.
2.2 CR-Calculus or Wirtinger calculus
Introduced by Wilhelm Wirtinger in 1927 , 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 given by Eq. (1) if and have continuous partial derivatives with respect to and , yielding
Since we have
and by setting to zero, it follows that
Note that the Cauchy-Riemann conditions for to be analytic in can be expressed compactly using the gradient as , i.e., is a function of only .
Similarly, if we take the derivative of with respect to , that is,
By setting to zero, we get
Again, the Cauchy-Riemann conditions for to be analytic in can be expressed compactly using the gradient as , i.e., is a function only of .
In other words, the gradient (respectively conjugate gradient) operator acts as a partial derivative with respect to (respectively to ), treating (respectively ) as a constant. Formally, we have
As an example, let be a real function of complex variable which is the squared Euclidean distance to the origin, with . Then,
as ; clearly the Cauchy-Riemann equations do not hold, and hence 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., . 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.
3. 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 .
3.1 Three-angle complex rotation algorithm
The first studied algorithm is the three-angle complex rotations (TACR), which is derived in polar coordinates . Nonetheless, the key idea behind 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.
3.2 Complex-valued fast givens rotations
On the other hand, the fast plane rotation  that is derived in complex plane and rectangular coordinates is a square root- and division-free Givens rotations . In this sense, the fast plane rotation which is also referred as the complex-valued fast Givens rotations () is a very efficient algorithm, aiming a QR-decomposition 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 performance will be compared to the well-known approaches which were successfully applied to the power system state estimation [12, 13, 14], 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.,  and , to cite a few.
The complex fast Givens transformation is computed using Algorithm 1 such that the second component of is zero and is a diagonal matrix, as shown below:
Notice the superscript denotes the Hermitian operation, i.e., complex conjugate transpose.
Algorithm 1. Complex fast Givens transform.
[,, r, type,,] = fast.givens (f, g,,);
type; ; ;
type; ; ;
type; ; ;
type; ; ;
The matrices , , and 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.
[,, type,]=fast.givens_QRD (A);
Step 1: Getandusing Algorithm 1:
[,, R, type,,]=fast.givens (R, new,,);
Step 2: Update elements based ontype
ifand type=1 then
else ifand type=2 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.
4. 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 that satisfies
where in Eq. (21), is a vector of specified quantities, i.e., constant term; is a vector of calculated quantities at each iteration. Consequently, the linearization of Eq. (21) from one step to the sequel leads to
where is the complex-valued Jacobian matrix which the dimension is . 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, 4], 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 . Thereby, the Jacobian matrix in expanded form may be represented through four partitions matrix, yielding
4.1 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 , resulting
where is an upper triangular matrix of dimension () and comprises the corresponding rows in the updated rhs vector of dimension (). 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 , 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 , 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.
4.2 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 , i.e., , is
where 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), is the general off-nominal tap transformer model which is composed by an ideal transformer with complex turns ratio in series with its admittance or impedance. Thus, if the corresponding branch is referred to.
1. Off-nominal tap transformer: and .
2. Pure-shifter: and .
3. Phase-shifter: .
4. transmission line: and .
4.3 Small example
The power flow model in complex plane as detailed in  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 (), are presented in Table 1.
|Bus||Specified quantities in pu|
The nodal admittance matrix, , 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.
4.4 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 is assumed.
|Algorithms||Number of iterations||Time/iteration (s)||Total time (s)|
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.
5. Complex-valued weighted-least-squares method (CV-WLS)
As shown in , the complex-valued WLS state estimator minimizes an objective function defined as
where = is a vector of specified complex values of dimension (), = is a vector of complex-valued state variables of dimension (), () is a vector of nonlinear functions of dimension () that maps to , is a vector of random errors in complex plane which dimension is (), and is a Hermitian positive-definite covariance matrix of which dimension is (). The superscript stands for Hermitian operator, i.e., the transpose complex conjugate operation. Thus, the necessary condition of optimality is given by
By applying a first-order Taylor series expansion of about , we get
yielding the updated estimated state vector expressed as
Notice and . Thus, the iterations are stopped when
where is the infinity norm and is the iteration counter.
Note that in Eq. (40), is the Jacobian matrix of dimension () defined in the complex domain, yielding
Let and be Jacobian submatrices of dimension (). 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 is a real-valued function of complex variables, the following property holds
where is a swap operator that permutes blocks of rows or blocks of columns depending upon whether pre-multiples or post-multiples a matrix, respectively. Moreover, this operator is an isomorphism from to the dual space , which obeys the properties . It shows that is symmetric and is equal to its own inverse, that is, . For instance, as shown in , this matrix is defined as
where is the ()-identity matrix.
Now, the complex-valued gain matrix in expanded form can be expressed as
where and , all of dimension (). Then, it follows from Eq. (47) that
Similarly, we get
The investigation of the sparsity structure of the Jacobian matrix in complex plane given by Eq. (49), e.g., Figure 4, reveals that the Jacobian matrices in the -domain are sparser than its counterpart in the -domain . Figure 4 is referred to the Brazilian equivalent 730-bus system which the Jacobian matrix in the -domain has 10,396 nonzero elements, while in the -domain, it has 18,403 nonzero elements; it is about 45% sparser.
Nonetheless, the recommended numerical procedure aiming to solve the complex-valued power system state estimation problem is addressed by solving the weighted form of the right-hand-side (rhs) of Eq. (40) instead of Eq. (44), yielding
where is of dimension () and is of dimension (). Thus, the incremental changes in the state vector are calculated via
where the operator is defined as the Moore-Penrose pseudoinverse .
Aiming to avoid explicitly store, the -matrix, we apply the -transformation to the augmented matrix, , given by
By storing the rotations in compact form, the complex-valued Jacobian matrix can be kept constant and only the right-hand-side vector is updated throughout the final iterations. The solution of the state vector increment given by Eq. (55) is found by executing a simple back-substitution of Eq. (56) after performing unitary transformations to the latter matrix, resulting in
Here, is an upper triangular matrix of dimension (), and comprises the corresponding rows in the updated rhs vector of dimension (). Finally, Eq. (55) is solved by performing a back-substitution via
5.1 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., and , is equal to the corresponding power flow measurement, i.e., and , respectively.
As observed in the above expressions of , the sub-matrix is more sparse than the sub-matrix , 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.
|0.0000 − j 0.3131||+1.4747 + j 0.1745||+1.8873 + j 0.4815||+1.8827 + j 0.4248|
|0.0000 − j 0.3131||−1.4014 − j 0.0707||−1.8036 − j 0.5094||−1.7997 − j 0.4500|
|0.0000 − j 0.3131||+1.4747 + j 0.1745||+1.8873 + j 0.4815||+1.8827 + j 0.4248|
|0.0000 − j 0.3131||−1.4014 − j 0.0707||−1.8036 − j 0.5094||−1.7997 − j 0.4500|
5.2 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 is assumed.
|CV-Jacobian matrix dimension||Algorithms||Number of iterations||Time/iteration (s)||Total time (s)|
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.
6. Conclusions and future developments
The chapter’s prime goal is to present the advances aiming to solve nonlinear system of equations in complex plane, regardless if it is exactly or overdetermined. Firstly, it develops the former Newton-Raphson method addressed to solve exactly determined problem. Thereafter, the weighted-least-squares (WLS) is developed through the Gauss-Newton method. In both cases the applications are carried out on small examples and larger systems. All results demonstrate the advantages of the algorithms developed in complex plane over the former procedure, although the computational burden is still a bottleneck for larger systems. Highlight that there are many enhancements that can be addressed to mitigate this problem. To cite a few, they are shown as follows:
Vector of the state variables in the conjugate coordinate system
Complex and complex conjugate state variables
Complex and complex conjugate tap position
Real and imaginary part of a complex variable
Complex-valued Jacobian matrix for the exactly determined system of equations
Complex-valued Jacobian matrix for the overdetermined system of equations
Complex-valued mismatch vector
Quantity in the conjugate coordinate system
Squared Euclidean norm