Synthesis of fast and superfast solvers of large systems of linear algebraic equations in Power Flow Calculation

Algorithms for fast and superfast solvers of large systems of linear algebraic equations, used in solving mathematical problems of electric power industry, are proposed. These algorithms are constructed based on a novel method of multistep decomposition of a multidimensional linear dynamic system. Examples of the analytical synthesis of iterative solvers for matrices of the general form and for large numerical systems of linear algebraic equations are presented. 1 Problem description One of the central objects in linear algebra is a large system of linear algebraic equations (SLAEs) Ax = b, (1) where ∈ × – Jacobi matrix of size n×n, ∈ – previously given vector of dimension n. Here and below, F denotes either the field of real numbers R or/and the field of complex numbers C. There is a large number of iterative methods and software packages of iterative solvers designed for system of linear equations Eq (1). However, these methods are effective under fairly strong conditions on the properties of the matrix А, among which there are cyclicity, symmetry, and positive definiteness [1–3]. Furthermore, no iterative solvers that could solve the equations with the general matrix А. are known. For example, the methods based on Krylov’s subspaces, such as the conjugate gradient method or the generalized minimal residual method of order r GMRES(r), demonstrate good global convergence for positive definite and normal matrices [1]. In the absence of the normality property of the convergence condition become fairly complex, maybe even chaotic, and they are poorly understood. For example, in [4] one can find examples of simple equations (1), for which the iterations of GMRES(1) converge rapidly reaching an accuracy 10–15, while the error of the process GMRES(2) demonstrates chaotic behavior within the range 0,1...1. The aim of this work is to develop fast and superfast iterative solvers based on a novel method of multistep decomposition and synthesis of control of a multidimensional dynamic system [5] for SLAEs. These solvers must ensure rapidly converging iterative processes in numerical computations. In other words, we want these iterative solvers to be stable and converging within several first iterations independently of the properties of the matrices in the SLAEs. The minimization of algorithmic complexity is not considered here. 2 Analysis of iterative solvers based on control theory Parameters of iterative solvers, such as shift, relaxation, etc., can be considered as controls, and the solvers themselves can be considered as controlled dynamic systems. One of the first approaches in this direction can be found in the early book by Bellman [6]: «...We have already referred to the fact that computing can be considered to be a control process in which we want to blend complexity, time, and accuracy in some appropriate way. It should, in addition, be treated in many cases as an adaptive control process in which the results of the previous calculation are used to guide the subsequent calculations, producing not only a choice of parameters but a choice of algorithms. » In recent times a number of numerical algorithms have been studied from the viewpoint of stabilization of dynamic systems using feedback control synthesis. Within this direction of research, [1] cites the papers [7, 8], where the choice of the step size in ODE-solvers was improved by finding an appropriate control. In [9–11] , the concept of controllability was applied to finding eigenvalues matrix pencils and solvers of linear matrix equations (1). In [12,13] methods of the robust stability theory were used to optimize an SLAE solver. Let us consider some popular iterative solvers for the linear matrix equation (1), based on the simple iteration , 0 Web of Conferences https://doi.org/10.1051/e3sconf/20191390 E3S 139 (2019) 1


Problem description
One of the central objects in linear algebra is a large system of linear algebraic equations (SLAEs) where ‫ܣ‬ ∈ ‫ܨ‬ × -Jacobi matrix of size n×n, ܾ ∈ ‫ܨ‬previously given vector of dimension n. Here and below, F denotes either the field of real numbers R or/and the field of complex numbers C. There is a large number of iterative methods and software packages of iterative solvers designed for system of linear equations Eq (1). However, these methods are effective under fairly strong conditions on the properties of the matrix А, among which there are cyclicity, symmetry, and positive definiteness [1][2][3]. Furthermore, no iterative solvers that could solve the equations with the general matrix А. are known. For example, the methods based on Krylov's subspaces, such as the conjugate gradient method or the generalized minimal residual method of order r GMRES(r), demonstrate good global convergence for positive definite and normal matrices [1]. In the absence of the normality property of the convergence condition become fairly complex, maybe even chaotic, and they are poorly understood. For example, in [4] one can find examples of simple equations (1), for which the iterations of GMRES(1) converge rapidly reaching an accuracy 10-15, while the error of the process GMRES(2) demonstrates chaotic behavior within the range 0,1…1.
The aim of this work is to develop fast and superfast iterative solvers based on a novel method of multistep decomposition and synthesis of control of a multidimensional dynamic system [5] for SLAEs. These solvers must ensure rapidly converging iterative processes in numerical computations. In other words, we want these iterative solvers to be stable and converging within several first iterations independently of the properties of the matrices in the SLAEs. The minimization of algorithmic complexity is not considered here.

Analysis of iterative solvers based on control theory
Parameters of iterative solvers, such as shift, relaxation, etc., can be considered as controls, and the solvers themselves can be considered as controlled dynamic systems. One of the first approaches in this direction can be found in the early book by Bellman [6]: «…We have already referred to the fact that computing can be considered to be a control process in which we want to blend complexity, time, and accuracy in some appropriate way. It should, in addition, be treated in many cases as an adaptive control process in which the results of the previous calculation are used to guide the subsequent calculations, producing not only a choice of parameters but a choice of algorithms. » In recent times a number of numerical algorithms have been studied from the viewpoint of stabilization of dynamic systems using feedback control synthesis. Within this direction of research, [1] cites the papers [7,8], where the choice of the step size in ODE-solvers was improved by finding an appropriate control. In [9][10][11] , the concept of controllability was applied to finding eigenvalues matrix pencils and solvers of linear matrix equations (1). In [12,13] methods of the robust stability theory were used to optimize an SLAE solver. Let us consider some popular iterative solvers for the linear matrix equation (1), based on the simple iteration (Richardson's) method and the methods of polynomial and linear iteration [1]. For invertible matrix A Eq.(1) has the well-known solution The iterative solver based on Richardson's method can be classified as a discrete bilinear control system with a single input and multiple outputs (SIMO), defined over Fn: Here ‫ݑ‬ -is the control parameter (scalar control) and x -is the initial condition. The fixed (steady) point in n-dimensional space of this iterative process is solution (2) to the linear equation (1). Various strategies for ‫ݑ‬ and certain families of matrices were proposed [14,15]. In particular, for ‫ݑ‬ = ‫ݑ‬ = const the discrete system (3) has a steadystate value (is asymptotically stable), if the set of eigenvalues (denoted by) of the matrix E − ‫ݑ‬A, where E -is the identity matrix of order n, is within the unit circle on the complex plane, i.e.∀ߣ ∈ eig(E − ‫ݑ‬A), |ߣ| < 1.
If ‫ݑ‬ is a feedback control, then ‫ݑ‬ = ୣ ೖ * ୣ ೖ ‖ୣ ೖ ‖ మ , e = b − Ax , where e -is the current residual, then This approach is used in the method GMRES(1). It is known that the dynamic system(3) with this control -is asymptotically stable and the iterative process converges if the matrixA + A * ≻ 0 is positive definite. Clearly, this condition is not satisfied for arbitrary matrices. In [13], an alternative to the discrete SIMO-system (3) was proposed in the form of a multi input multi output (MIMO) linear system: where is a given matrix that generally is not directly related to The asymptotic stability and, therefore, the convergence of the iterative process (4) can be ensured using the feedback control that can be designed by linear-quadratic optimization. This yields the known LQRES [1] solver.

Controllability of iterative solvers
It is known that a necessary and sufficient condition for the efficiency of feedback control (5) is the complete controllability condition of the dynamic system [6]. In [1], the controllability of the iterative solver based on Richardson's method (3) was analyzed. It was shown that, for every input sequence ‫ݑ(‬ ), the trajectory of the state (x ) converges to the solution of Eq. (2), if and only if the sequence of residuals converges to zero. Thus, the dynamics of the discrete system (3) is equivalent to the dynamics of (6). The iterative solver (3) based on Richardson's method is completely controllable, when A has n different real eigenvalues [1]. Otherwise, controllability is not guaranteed.
The controllability of the iterative solver based on the discrete MIMO-system (4) can be determined using the well-known tests [5]. It is clear that (4) is completely controllable if and only if the Kalman rank condition is satisfied 4. Linear control of iterative solvers based on the multistep decomposition method Consider a method of linear control of linear multidimensional dynamic systems that can be used for designing fast and superfast iterative solvers of SLAEs (1) [17]. Under the condition that the matrix A, is invertible, solution (2) can be written in the form of the power series [3] x Thus, the inverse matrix A ିଵ can be determined using the MIMO system (4). The behavior of residual e is described by the linear system [1] An analysis of (8) proves the following propositions. 4) The discrete MIMO system (8) is completely controllable if and only if the pair (E − A, −AG) is completely controllable in the sense of criterion (7). 5) If the MIMO system (8) is completely controllable, then there exists a feedback control (5) (which is not unique in the case m>1), such that the asymptotic stability condition holds: Let, as before, A ෩ = E − A, G ෩ = −AG. Denote by G ෩ ୄ the matrix that is a left divisor of the matrix G and, therefore, satisfies the conditions [5]: where G ෩ ୄ ା is a pseudoinverse matrix for G ෩ ୄ . Determine the controller matrix K in the formula for the feedback control (5) by , Here, / is a matrix with the prescribed (desired) eigenvalues. It can be shown that A drawback of this feedback control (10) is that the stability of the matrix G ෩ ୄ A ෩ G ෩ ୄ ା cannot be guaranteed in the general case. For this reason, we consider the multistep decomposition of system determined as follows: initial state the first step the kth (intermediate) step the Sth (final) step (14) Here ܵ=floor ቀ ቁ − 1; floor -the operation of rounding a number ݊/݉ to the nearest integer in the smaller side, e.g., floor(0,4) = 0, floor(2,9) = 2, etc.
If we also determine in the step-by-step manner the matrices Then, after closing the system by the feedback control (5), with the matrix K given by we obtain the equality of the sets of eigenvalues: Here Λ ିଵ (݅ − 1 = 0,1, . . . , ܵ) -are the matrices with the prescribed eigenvalues at each level of decomposition and ⋃ eig(Λ ିଵ ) = eig(Λ) ௌାଵ ୀଵ . Thus, the multistep decomposition of system (11) -(14) and the procedure of determining the matrix K (19) based on (15) -(18) ensure the prescribed location of the eigenvalues ߉ = ൛ߣ̑ଵ, ߣ̑ଶ, . . . ߣ̑ൟ, as specified in (20). Consider the problem of designing control (5), that guarantees that the transient process of the discrete MIMO system takes a finite time. In this case, the matrix A ෩ − G ෩ K has only zeros as its eigenvalues [18], i.e., eig൫A ෩ − G ෩ K൯ = {0,0, . . . ,0} (21) This requirement implies that we can take any nilpotent matrices with the index of nilpotency not exceeding ‫ݎ‬ as the matrices Λ [2,5]. In this case, the convergence of the designed solvers considered as discrete systems will take place not more than n steps. We call such solvers fast. If we use the zero matrices as Λ , then this case can be assigned the discrete system x ାଵ = b, in which the convergence is obviously already at the second step, i.e., x = 0, x ଵ = b. Such solvers can be called superfast.
For the zero matrices Λ of the corresponding size, the final stage of algorithm (15) - (18) is It follows from (22) -(23), that, if the number of decomposition steps is ܵ = 1, then the formula for the matrix K in control (5), that guarantees not only the finite duration of the transient process but also the superfast convergence (during the first two iterations independently of the size of the discrete MIMO system) takes the form: For ܵ = 2 и ܵ = 3, the formulas for K are somewhat more complicated: The analysis of the final formulas (24) -(26), shows that the matrices K obtained using the multistep decomposition contain enclosures whose form is determined both by the matrix A of system (1), and by an arbitrarily chosen matrix G. and by an arbitrarily chosen matrix G , we have a wide choice of the matrices K. For example, if we require in advance that G has not less than n/2, columns, then control with the matrix K (24) guarantees that the iterative process has fast convergence rate:  where r, w, O are given numbers. The nonciclicity of matrix (28) is directly indicated by the diagonal 2×2 blocks the first of which is a Jordan block and the second one is an upper triangular matrix. The fact that this matrix is not normal is verified straightforwardly:‫ܣ‬ ் ‫ܣ‬ ≠ ‫ܣܣ‬ ் . Let the vector b be also given in the general form

Examples of the synthesis of superfast iterative solvers
Specify the matrix G by: We construct the iterative solver using algorithm (11) -(18), (22) -(23). It is seen from the analysis of matrices (28) и (30), that in this case ܵ=floor(݊/݉) − 1 = floor(4/2) − 1 = 2 − 1 = 1 therefore, after performing the two-step decomposition consisting of the zero and the first step, we may use formula for the matrix K (24). Thus, according to (24), we obtain In this case eig൫A ෩ − G ෩ K൯ = {0,0,0,0}, which was to be proved. By running the iterative solver (27) with the parameters for the zero initial conditions indicated above, we obtain the following sequence of residuals > @ We see that the iterative process converged to the exact solution already at the second iteration step. An analysis of the proposed superfast solvers of linear equations with the matrices and vectors the elements of which are distributed normally and the eigenvalues of the matrices A re distributed by Girko's circular law [19]: A = randn(1000 × 1000), G = randn(1000 × 500), b = randn(1000 × 1), A = randn(5000 × 5000), G = randn(5000 × 2500), b = randn(5000 × 1), showed that in all the cases the iterative process converged at the third or fourth iteration step independently of the size of the problem. For each value of the problem size, at least n+1 examples were, n -dimension of the vector in the right part (1). Superfast iterative solver for power flow calculation of large electric power system Let's consider the example the calculation of superfast iterative solver for a model of a real power system [20].
It's scheme includes several concentrated subsystems connected by relatively weak connections. Electric power system considers 286 buses, 531 branches, 129 synchronous generators. There are several variants of the Jacobi matrix A for this system, including size 258×258. «Portrait» of the specified matrix for one of the regimes of electric power system is shown in Fig.  1.
Computing the determinant of the Jacobi matrix yields a number that close to infinite ‫ݐ݁݀‬ ‫ܣ‬ = ܽ = ߣ ଵ ߣ ଶ ⋯ ߣ ଶହ଼ = 6,310810258, i.e. it's very difficult to solving SLA (1) by the standard method by calculating the attached matrix. In addition, the Jacobi matrix demonstrates weak scalability, which means large differences in the values of elements (in this case items vary from -500 to 1500).  It is necessary to adjust the matrix A in such a way that fast convergence of iterative processes in multivariate calculations of electric power systems regimes is carried out. The «Portrait» of the corrected Jacobian matrix of electric power system and the distribution of its eigenvalues are shown in Fig. 1 and 2, respectively. As we see at the Fig. 1 the scaling of the matrix is quite satisfactory (the values of the elements vary from -3 to 2). From Fig. 2 it follows that all eigenvalues of the corrected matrix with high accuracy (10-7…10-6) are located at the beginning of the coordinates of the com-Plex plane. This provides the highest rate of convergence of the iterative process with the error rate of 10-13... 10-12. In this case, the convergence of the iterative process in the analyzed case occurs on the second iteration independently of the right part of the SLA (1).

Conclusion
Algorithms for fast and superfast solvers of large SLAEs based on a novel method of the multistep decomposition of a linear multidimensional dynamic system are proposed. Examples of the analytical design of iterative solvers for matrices of the general form demonstrated that the iterative processes converge already at the second iteration step. The statistical analysis of the constructed superfast solvers of linear equations with matrices and vectors consisting of normally distributed elements showed that in all the cases the iterative process converged at the third or fourth iteration independently of the size of the problem. The difficulty of the implementation of the proposed iterative algorithm could be caused by the absence of efficient procedures for computing pseudoinverse matrices of rank deficient matrices. However, this difficulty can be overcome as described in [21]. Application of the fast solver in the practical problem of calculation of the real EPS mode it is shown that the convergence of the iterative process with the residual norm 10-13…10-12 starts on the second iteration independently of the right side of the SLA (1).