# Computationally Efficient Optimal Control for Unstable Power System Models

Authors: Mahtab Uddin; M. Monir Uddin; Md. Abdul Hakim Khan
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this article, the focus is mainly on gaining the optimal control for the unstable power system models and stabilizing them through the Riccati-based feedback stabilization process with sparsity-preserving techniques. We are to find the solution of the Continuous-time Algebraic Riccati Equations (CAREs) governed from the unstable power system models derived from the Brazilian Inter-Connected Power System (BIPS) models, which are large-scale sparse index-1 descriptor systems. We propose the projection-based Rational Krylov Subspace Method (RKSM) for the iterative computation of the solution of the CAREs. The novelties of RKSM are sparsity-preserving computations and the implementation of time-convenient adaptive shift parameters. We modify the Low-Rank Cholesky-Factor integrated Alternating Direction Implicit (LRCF-ADI) technique-based nested iterative Kleinman–Newton (KN) method to a sparse form and adjust this to solve the desired CAREs. We compare the results achieved by the Kleinman–Newton method with that of using the RKSM. The applicability and adaptability of the proposed techniques are justified numerically with MATLAB simulations. Transient behaviors of the target models are investigated for comparative analysis through the tabular and graphical approaches.


---

## Body

## 1. Introduction

For the practical purposes, optimal control is a vital part of the engineering interest, e.g., industrial control systems, system defense strategy, voltage stability of the power systems, and signal processing [1, 2]. Multitasking systems having various components arise in many fields of engineering applications, such as microelectronics, microelectromechanical systems, cybersecurity, computer control of industrial processes, communication systems, and the like. These systems are composed of branches of subsystems and are functioned by large mathematical models utilizing the interrelated inner mathematical system of higher dimensions.Power system models are one of the prime branches of the application of optimal controls. For the multiconnected power systems, avoiding cyber-attacks, ensuring the stability of the power connections, and maintaining the compatible frequency level are essential. For those inevitable issues, optimal controls have exigent roles [3].The dynamic of a large-scale power system model can be described by the Differential Algebraic Equations (DAEs) as(1)x˙t=fx1,x2,P,f:ℝn1+n2+n3⟶ℝn1,0=gx1,x2,P,g:ℝn1+n2+n3⟶ℝn2,where, x1∈X1⊂ℝn1 is the vector with differential variables, x2∈X2⊂ℝn2 is the vector with algebraic variables, and P⊂ℝn3 is the vector of parameters with n1+n2=n [4, 5]. In the state-space representation, the dynamic state variables x1 and instantaneous variables x2 are defined for the specific system, where the parameter P defines the configuration and the operation condition.The state variables inx1 are time-dependent generator voltages, and the parameter P is composed of the system parameters. The control devices together form fx1,x2,P, and the power flow balance forms gx1,x2,P. In case of voltage stability, some equations of the fx1,x2,P will not be considered. For a fixed parameter P0, linearizing system (1) around the equilibrium point will provide the following Linear Time-Invariant (LTI) continuous-time system with the input-output equations in the sparse form with the block matrices as(2)E1000︸Ex˙1tx˙2t︸x˙t=J1J2J3J4︸Ax1tx2t︸xt+B1B2︸But;xt0=x0,yt=C1C2︸Cx1tx2t︸xt+Dut,where E,A∈ℝn×n, B∈ℝn×p, C∈ℝm×n, and D∈ℝm×p with very large n and p,m≪n represent differential coefficient matrix, state matrix, control multiplier matrix, state multiplier matrix, and direct transmission map, respectively [6]. In system (2), xt∈ℝn is the state vector and ut∈ℝp is control (input), while yt∈ℝm is the output vector and considering xt0=x0 as the initial state. In most of the state-space representations, the direct transmission remains absent and because of that D=0. Because E is singular (i.e., detE=0), system (2) is called the descriptor system [7, 8].Here,x1∈ℝn1, x2∈ℝn2 with n1+n2=n are state vectors, and other submatrices are sparse in appropriate dimensions. If E1 and J1 have full rank, and J4 is nonsingular (i.e., detJ4≠0), the system is called the index-1 descriptor system [9]. In the current work, we will focus on the stabilization of index-1 descriptor system only. By proper substitution and elimination, descriptor system (2) can be converted to the generalized LTI continuous-time system(3)Ex˙t=Axt+ℬut,yt=Cxt+Dut,where we have considered the following relations:(4)x≔x1,E≔E1,A≔J1−J2J4−1J3,ℬ≔B1−J2J4−1B2,C≔C1−C2J4−1J3,D=D−C2J4−1B2.Lemma 1 (equivalence of transfer functions [10]).
Assume fors∈ℂ, the transfer functions Gs=CsE−A−1B+D and Gs=CsE−A−1ℬ+D are obtained from index-1 descriptor system (2) and converted generalized system (3), respectively. Then, the transfer functions Gs and Gs are identical, and hence, those systems are equivalent.Though systems (2) and (3) are equivalent, system (2) is sparse and system (3) is dense. This explicit conversion is contradictory with the aim of the work, and it will be bypassed by a more efficient way.LTI continuous-time systems are the pivot ingredient of the present control theory and many other areas of science and engineering [11]. Continuous-time Algebraic Riccati Equation (CARE) appears in many branches of engineering applications; especially in electrical fields [12, 13]. The CARE connected to system (3) is defined as(5)ATXE+ETXA−ETXℬℬTXE+CTC=0.If the Hamiltonian matrix corresponding to system (3) has no pure imaginary eigenvalues, then the solution X of the CARE (5) exists and unique [14]. The solution X of (5) is symmetric positive definite and called stabilizing for the stable closed-loop matrix A−ℬℬTXE. Riccati-based feedback matrix has a prime role in the stabilization approaches for unstable systems [15, 16]. To find an optimal feedback matrix Ko, the Linear Quadratic Regulator (LQR) problem technique can be applied, where the cost functional is defined as(6)Ju,x0=∫0∞xTtCTCxt+uTtutdt.Cost functional (6) can be optimized as Juo,x0=x0TXx0 by applying the optimal control uo=−Koxt generated by the optimal feedback matrix Ko=ℬTXE associated with the solution matrix X of the CARE (6). Using the optimal feedback matrix Ko, an unstable LTI continuous-time system can be optimally stabilized by replacing A with As=A−ℬKo. The stabilized system can be written as(7)Ex˙t=Asxt+ℬut,yt=Cxt+Dut.The eligibility of Rational Krylov Subspace Method (RKSM) for the large-scale LTI continuous-time systems has discussed by Simoncini [17], and the application of adaptive RKSM to solve large-scale CARE for finding optimal control of the LTI systems has narrated by Druskin and Simoncini [18]. Analysis of the basic properties of RKSM for solving large-scale CAREs subject to LTI systems was investigated by Simoncini [19], where the author briefed a new concept of shift parameters that is efficient for the perturbed systems. Very detailed discussion on the numerical solution of large-scale CAREs and LQR-based optimal control problems are given by Benner et al. [20], where the authors narrated the extensions Newton method by means of Alternating Direction Implicit (ADI) technique with convergence properties and the comparative analysis with the RKSM approach. For solving large-scale matrix equations, the Kleinman–Newton method based on Low-Rank Cholesky-Factor ADI (LRCF-ADI) technique is discussed by Kürschner [21].Because there is a scarcity of efficient computational solvers or feasible simulation tools for large-scale CAREs governed from the unstable power system models, the concentration of this work is to develop sparsity-preserving efficient techniques to find the solution of those CAREs. We are proposing a sparsity-preserving and rapid convergent form of the RKSM algorithm for finding solutions of CAREs associated with the unstable power system models and implementation of Riccati-based feedback stabilization. Also, a modified sparse form of LRCF-ADI-based Kleinman–Newton method for solving CAREs subject to unstable power system models and corresponding stabilization approach is proposed. The proposed techniques are applied for the stabilization of the transient behaviors of unstable Brazilian Inter-Connected Power System (BIPS) models [22]. Moreover, the comparison of the computational results will be provided in both tabular and graphical methods.


## 2. Preliminaries

In this section, we discuss the background of the proposed techniques, which are derived for generalized LTI continuous-time systems. It includes the derivation of the Rational Krylov Subspace Method (RKSM) technique for solving the Riccati equation, basic structure of the Kleinman–Newton method, and derivation of Low-Rank Cholesky-Factor Alternative Direction Implicit (LRCF-ADI) method for solving Lyapunov equations.


### 2.1. Generalized RKSM Technique for Solving Riccati Equations

In a study [19], Simoncini applied RKSM approach for solving the CARE in the form(8)ATXE+ETXA−ETXBBTXE+CTC=0,associated with the LTI continuous-time system(9)Ex˙t=Axt+But,yt=Cxt+Dut.If the eigenvalues of the matrix pairA,E satisfy λi+λ¯j≠0,∀i,j=1,2,…,m that ensures that the solution X of the CARE (8) exists and unique. The orthogonal projector V∈ℝn×m spanned by the m-dimensional rational Krylov subspace for a set of shift parameters μi∈ℂ+;i=1,2,…,m is defined as(10)Km=spanCT,AT−μ1ET−1CT,…,∏i=1mAT−μiET−1CT.Ifθj are the eigenvalues of VTAV,VTEV;Sm∈ℂ+ approximates the mirror eigenspace of A−BBTXE and δSm is its border, the shifts are computed from(11)μi+1=argmaxμ∈δSm∏j=1iμ−μj∏j=1iμ−θj.According to the Galerkin condition and after simplification by matrix algebra, a low-rank CARE can be obtained as(12)A^TX^E^+E^TX^A^−E^TX^B^B^TX^E^+C^TC^=0,where X^=VTXV, E^=VTEV, A^=VTAV, B^=VTB, and C^=CV. Equation (12) is an approximated low-rank CARE and can be solved by any conventional method or MATLAB care command. Here, X^∈ℝm×m is taken as low-rank approximation of X, corresponding to the low-rank CARE (12). Then, residual of the m-th iteration is(13)ℛF=SJSTF;J=0I0I0I0I0,where ⋅F denotes the Frobenius norm and S is a block upper triangular matrix in the QR factorization of the matrix U defined as(14)U=vm+1μm+1ETVmX^Hm−Temhm+1,mT−In−VmVmTATvm+1T,where Hm is a block upper Hessenberg matrix, and em is the matrix formed by the last p columns of the mp-order identity matrix. For CT=Q0R0 such that R0=β0, the relative residual can be estimated as(15)ℛFRelative=ℛFβ0Tβ0F.Through RKSM, the low-rank factorZ of the approximate solution X of the CARE (8) needs to be estimated, such that X≈ZZT. The low-rank solution X^ is symmetric positive definite, and the original solution X can be approximated as X=VX^VT. By the eigenvalue decomposition to the approximate solution X^ and truncating the negligible eigenvalues, the possible lowest order factor Z of X can be estimated as(16)X=VX^VT=VTΛTTVT=VT1T2Λ100Λ2T1TT2TVT=VT1Λ1T1TVT=VT1Λ11/2VT1Λ11/2T=ZZT.Here,Λ2 consists the negligible eigenvalues. Summary of the above process is given in Algorithm 1.Algorithm 1: Generalized RKSM.
Input:E,A,B,C, imax (number of iterations) and μi (initial shifts).Output: Low-rank factored solution Z such that X≈ZZT.(1)
ComputeCT=Q0R0 (QR factorization).(2)
ChooseV0=Q0.(3)
whilenot converged or i≤imaxdo(4)
Solvevm+1=AT−μm+1ET−1Vm.(5)
Compute shift for the next iteration.(6)
Using Arnoldi algorithm, orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vmv^m+1.(7)
AssumingE^=Vm+1TEVm+1, A^=Vm+1TAVm+1, B^=Vm+1TB, and C^=CVm+1, for X^, solve the reduced-order Riccati equation (12).(8)
ComputeℛFRelative for convergence.(9)
end while(10)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(11)
For negligible eigenvalues, truncateΛ2 and compute Z=Vm+1T1Λ11/2.


### 2.2. Modified LRCF-ADI Method for Solving Lyapunov Equations

The generalized Continuous-time Algebraic Lyapunov Equation (CALE) associated with LTI continuous-time system (9) is(17)ATXE+ETXA=−CTC.The shift parametersμi∈ℂ− are allowed, and the initial iteration is taken as X0=X0T∈ℝn×n. Assume Zi∈ℝn×ip as the low-rank Cholesky factor of Xi∈ℝn×n such that Xi=ZiZiT [23]. The ADI algorithm can be formed in terms of Cholesky factor Zi of Xi, and there will be no need to estimate or store Xi at each iteration as only Zi is required [24].Consideringγi=−2Reμi, the conventional low-rank Cholesky factor ADI iterations yield the form as(18)V1=AT+μ1ET−1CT,Z1=γ1V1=−2Reμ1AT+μ1ET−1CT,Vi=Vi−1−μi+μ¯i−1AT+μiET−1ETVi−1,Zi=Zi−1γiVi=Zi−1−2ReμiVi.Assume a set of adjustable shift parameters, for two subsequent block iteratesVi and Vi+1 of the ADI technique related to the pair of complex conjugated shifts μi,μi+1≔μ¯i, it holds(19)Vi+1=V¯i+2δiImVi,where ⋅¯ indicates the complex conjugate with δi=Reμi/Imμi and iterates associated to real shifts are always purely real [25].Then, the following matrix for the basis extension can be obtained:(20)ViVi+1=γiReVi+δiImViγiδi2+1ImVi.Then, for a pair of complex conjugate shifts at any iteration, the low-rank factorZi+1 can be computed as(21)Zi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.The modified techniques discussed above is summarized in Algorithm2.Algorithm 2: Generalized LRCF-ADI (modified).
Input:E,A,C, imax (number of iterations) and shift parameters μjj=1imax.Output: Low-rank Cholesky-factor Z such that X≈ZZT.(1)
ConsiderZ0=.(2)
fori←1 to imaxdo(3)
ifi=1then(4)
SolveV1=AT+μ1ET−1CT.(5)
else(6)
ComputeVi=Vi−1−μi+μ¯i−1AT+μiET−1ETVi−1.(7)
end if(8)
ifImμi=0then(9)
UpdateZi=Zi−1−2μiVi.(10)
else(11)
Assumeγi=−2Reμi and δi=Reμi/Imμi.(12)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.(13)
ComputeVi+1=V¯i+2δiImVi.(14)
i=i+1(15)
end if(16)
end for


### 2.3. Basic Kleinman–Newton Method

For system (9), the residual ℛX of the CARE (8) is(22)ℛX=ATXE+ETXA−ETXBBTXE+CTC.From the Fréchet derivative atX=Xi, we have(23)ℛ′Xi=A−BBTXETXiE+ETXiA−BBTXE.Consider the Newton iteration and apply (23), we have(24)ℛ′ΔXi+ℛXi=0.Now, putΔXi=Xi+1−Xi in (24). Then, after simplification, we get(25)A−BBTXiETXi+1E+ETXi+1A−BBTXiE=−CTC−ETXiBBTXiE.Then, assumeA˜i=A−BBTXiE and Wi=CTETXiB, then equation (25) reduces to a generalized Lyapunov equation such as(26)A˜iTXi+1E+ETXi+1A˜i=−WiWiT.Generalized CALE (26) can be solved for Xi+1 by any conventional method, such as Low-Rank Cholesky-Factor Alternative Direction Implicit (LRCF-ADI) method and the corresponding feedback matrix Ki+1=BTXi+1E can be estimated. The whole mechanism is called the Kleinman–Newton method for solving generalized CARE [26]. The summary of the method is given in Algorithm 3.Algorithm 3: Basic Kleinman–Newton.
Input:E,A,B,C, imax (number of iterations) and X0 (initial assumption).Output: Approximate solution X and feedback matrix K.(1)
whilei≤imaxdo(2)
ComputeA˜i=A−BBTXiE and Wi=CTETXiB.(3)
ForXi+1 solve A˜iTXi+1E+ETXi+1A˜i=−WiWiT.(4)
ComputeKi+1=BTXi+1E.(5)
end while


## 2.1. Generalized RKSM Technique for Solving Riccati Equations

In a study [19], Simoncini applied RKSM approach for solving the CARE in the form(8)ATXE+ETXA−ETXBBTXE+CTC=0,associated with the LTI continuous-time system(9)Ex˙t=Axt+But,yt=Cxt+Dut.If the eigenvalues of the matrix pairA,E satisfy λi+λ¯j≠0,∀i,j=1,2,…,m that ensures that the solution X of the CARE (8) exists and unique. The orthogonal projector V∈ℝn×m spanned by the m-dimensional rational Krylov subspace for a set of shift parameters μi∈ℂ+;i=1,2,…,m is defined as(10)Km=spanCT,AT−μ1ET−1CT,…,∏i=1mAT−μiET−1CT.Ifθj are the eigenvalues of VTAV,VTEV;Sm∈ℂ+ approximates the mirror eigenspace of A−BBTXE and δSm is its border, the shifts are computed from(11)μi+1=argmaxμ∈δSm∏j=1iμ−μj∏j=1iμ−θj.According to the Galerkin condition and after simplification by matrix algebra, a low-rank CARE can be obtained as(12)A^TX^E^+E^TX^A^−E^TX^B^B^TX^E^+C^TC^=0,where X^=VTXV, E^=VTEV, A^=VTAV, B^=VTB, and C^=CV. Equation (12) is an approximated low-rank CARE and can be solved by any conventional method or MATLAB care command. Here, X^∈ℝm×m is taken as low-rank approximation of X, corresponding to the low-rank CARE (12). Then, residual of the m-th iteration is(13)ℛF=SJSTF;J=0I0I0I0I0,where ⋅F denotes the Frobenius norm and S is a block upper triangular matrix in the QR factorization of the matrix U defined as(14)U=vm+1μm+1ETVmX^Hm−Temhm+1,mT−In−VmVmTATvm+1T,where Hm is a block upper Hessenberg matrix, and em is the matrix formed by the last p columns of the mp-order identity matrix. For CT=Q0R0 such that R0=β0, the relative residual can be estimated as(15)ℛFRelative=ℛFβ0Tβ0F.Through RKSM, the low-rank factorZ of the approximate solution X of the CARE (8) needs to be estimated, such that X≈ZZT. The low-rank solution X^ is symmetric positive definite, and the original solution X can be approximated as X=VX^VT. By the eigenvalue decomposition to the approximate solution X^ and truncating the negligible eigenvalues, the possible lowest order factor Z of X can be estimated as(16)X=VX^VT=VTΛTTVT=VT1T2Λ100Λ2T1TT2TVT=VT1Λ1T1TVT=VT1Λ11/2VT1Λ11/2T=ZZT.Here,Λ2 consists the negligible eigenvalues. Summary of the above process is given in Algorithm 1.Algorithm 1: Generalized RKSM.
Input:E,A,B,C, imax (number of iterations) and μi (initial shifts).Output: Low-rank factored solution Z such that X≈ZZT.(1)
ComputeCT=Q0R0 (QR factorization).(2)
ChooseV0=Q0.(3)
whilenot converged or i≤imaxdo(4)
Solvevm+1=AT−μm+1ET−1Vm.(5)
Compute shift for the next iteration.(6)
Using Arnoldi algorithm, orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vmv^m+1.(7)
AssumingE^=Vm+1TEVm+1, A^=Vm+1TAVm+1, B^=Vm+1TB, and C^=CVm+1, for X^, solve the reduced-order Riccati equation (12).(8)
ComputeℛFRelative for convergence.(9)
end while(10)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(11)
For negligible eigenvalues, truncateΛ2 and compute Z=Vm+1T1Λ11/2.


## 2.2. Modified LRCF-ADI Method for Solving Lyapunov Equations

The generalized Continuous-time Algebraic Lyapunov Equation (CALE) associated with LTI continuous-time system (9) is(17)ATXE+ETXA=−CTC.The shift parametersμi∈ℂ− are allowed, and the initial iteration is taken as X0=X0T∈ℝn×n. Assume Zi∈ℝn×ip as the low-rank Cholesky factor of Xi∈ℝn×n such that Xi=ZiZiT [23]. The ADI algorithm can be formed in terms of Cholesky factor Zi of Xi, and there will be no need to estimate or store Xi at each iteration as only Zi is required [24].Consideringγi=−2Reμi, the conventional low-rank Cholesky factor ADI iterations yield the form as(18)V1=AT+μ1ET−1CT,Z1=γ1V1=−2Reμ1AT+μ1ET−1CT,Vi=Vi−1−μi+μ¯i−1AT+μiET−1ETVi−1,Zi=Zi−1γiVi=Zi−1−2ReμiVi.Assume a set of adjustable shift parameters, for two subsequent block iteratesVi and Vi+1 of the ADI technique related to the pair of complex conjugated shifts μi,μi+1≔μ¯i, it holds(19)Vi+1=V¯i+2δiImVi,where ⋅¯ indicates the complex conjugate with δi=Reμi/Imμi and iterates associated to real shifts are always purely real [25].Then, the following matrix for the basis extension can be obtained:(20)ViVi+1=γiReVi+δiImViγiδi2+1ImVi.Then, for a pair of complex conjugate shifts at any iteration, the low-rank factorZi+1 can be computed as(21)Zi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.The modified techniques discussed above is summarized in Algorithm2.Algorithm 2: Generalized LRCF-ADI (modified).
Input:E,A,C, imax (number of iterations) and shift parameters μjj=1imax.Output: Low-rank Cholesky-factor Z such that X≈ZZT.(1)
ConsiderZ0=.(2)
fori←1 to imaxdo(3)
ifi=1then(4)
SolveV1=AT+μ1ET−1CT.(5)
else(6)
ComputeVi=Vi−1−μi+μ¯i−1AT+μiET−1ETVi−1.(7)
end if(8)
ifImμi=0then(9)
UpdateZi=Zi−1−2μiVi.(10)
else(11)
Assumeγi=−2Reμi and δi=Reμi/Imμi.(12)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.(13)
ComputeVi+1=V¯i+2δiImVi.(14)
i=i+1(15)
end if(16)
end for


## 2.3. Basic Kleinman–Newton Method

For system (9), the residual ℛX of the CARE (8) is(22)ℛX=ATXE+ETXA−ETXBBTXE+CTC.From the Fréchet derivative atX=Xi, we have(23)ℛ′Xi=A−BBTXETXiE+ETXiA−BBTXE.Consider the Newton iteration and apply (23), we have(24)ℛ′ΔXi+ℛXi=0.Now, putΔXi=Xi+1−Xi in (24). Then, after simplification, we get(25)A−BBTXiETXi+1E+ETXi+1A−BBTXiE=−CTC−ETXiBBTXiE.Then, assumeA˜i=A−BBTXiE and Wi=CTETXiB, then equation (25) reduces to a generalized Lyapunov equation such as(26)A˜iTXi+1E+ETXi+1A˜i=−WiWiT.Generalized CALE (26) can be solved for Xi+1 by any conventional method, such as Low-Rank Cholesky-Factor Alternative Direction Implicit (LRCF-ADI) method and the corresponding feedback matrix Ki+1=BTXi+1E can be estimated. The whole mechanism is called the Kleinman–Newton method for solving generalized CARE [26]. The summary of the method is given in Algorithm 3.Algorithm 3: Basic Kleinman–Newton.
Input:E,A,B,C, imax (number of iterations) and X0 (initial assumption).Output: Approximate solution X and feedback matrix K.(1)
whilei≤imaxdo(2)
ComputeA˜i=A−BBTXiE and Wi=CTETXiB.(3)
ForXi+1 solve A˜iTXi+1E+ETXi+1A˜i=−WiWiT.(4)
ComputeKi+1=BTXi+1E.(5)
end while


## 3. Solving Riccati Equations Arising from the Index-1 Descriptor Systems

In this section, we discuss the updated RKSM techniques for solving the Riccati equation derived from index-1 descriptor systems, it includes the stopping criteria, sparsity preservation, and estimation of the optimal feedback matrix. Then, the LRCF-ADI-based Kleinman–Newton method with the adjustment for index-1 descriptor systems and finding the optimal feedback matrix are discussed. Finally, the stabilized system is narrated accordingly.


### 3.1. Updated Rational Krylov Subspace Method

Let us consider LTI continuous-time system (3) and introduce an orthogonal projector V spanned by the m-dimensional rational Krylov subspace for a set of shift parameters μi∈ℂ+;i=1,2,…,m that is defined as(27)Km=spanCT,AT−μ1ET−1CT,…,∏i=1mAT−μiET−1CT.To find the shift parametersμi, the eigenvalues of the matrix pencil A,E are evaluated and sorted numerically. The required number of dominant eigenvalues are to store for the upcoming iterations. If some of the dominant eigenvalues have the negative real part, they are to shift to the positive complex plane by the mirror imaging process.Again, consider CARE (5) and apply the Galerkin condition on it. Then, after the simplification by matrix algebra, a low-rank CARE can be achieved as(28)A^TX^E^+E^TX^A^−E^TX^ℬ^ℬ^TX^E^+C^TC^=0,where X^=VTXV, E^=VTEV, A^=VTAV, ℬ^=VTℬ, and C^=CV. Equation (28) is a low-rank CARE and can be solved by MATLAB care command or any existing methods, such as Schur-decomposition method.For the quick and smooth convergence of the proposed algorithm, adjustable shift selection is crucial, and we are adopting the adaptive shift approach for index-1 descriptor systems [27]. This process required to be recursive, and in each step, the subspace to all the projector V generated with the current set of shifts will be extended.


#### 3.1.1. Stopping Criteria and Related Theorem

Arnoldi relation is a very essential tool for the computation of residual of the RKSM iterations. To avoid extra matrix-vector multiplies per iteration, the computation of projected matrixTm can be performed more efficiently than the explicit product Tm=VmTAVm. To find the stopping criteria, we have to consider the following lemmas.Lemma 2 (Arnoldi relation [28]).
LetKm=spanVm be the rational Krylov subspace with the shift parameters μi∈ℂ+;i=1,2,…,m. Then, Km satisfies the Arnoldi relation as follows:(29)ATVm=VmTm+v^m+1gmT;VmTVm=Im,where v^m+1β1=vm+1μm+1−In−VmVmTATvm+1 is the QR decomposition of the right-hand side matrix with gmT=β1hm+1,memTHm−1.Lemma 3 (building the projected matrix [29]).
Let the column vectors ofVm be an orthonormal basis of the rational Krylov subspace Km with the block diagonal matrix Dm=diagμ1ET,μ2ET,…,μmET, where μ1,μ2,…,μm is the set of shift parameters used in the algorithm. Then, for the projected matrix Tm=VmTAVm, the following relation holds:(30)Tm=Im+VmTETVmHmDm−VmTATvm+1hm+1,memTHm−1.Theorem 1 (residual of the RKSM iterations).
LetVm be the orthogonal projector spanned by the rational Krylov subspace Km and X=VmX^VmT is the solution of the CARE (5) using the low-rank solution X^. Then, the residual of m-th iteration can be computed as(31)ℛF=SJSTF;J=0I0I0I0I0,where ⋅F denotes the Frobenius norm, and S is a block upper triangular matrix in the QR factorization of the matrix U defined as(32)U=vm+1μm+1ETVmX^Hm−Temhm+1,mT−In−VmVmTATvm+1.Proof.
Assumef=In−VmVmTATvm+1 and consider the reduced QR factorization CT=V0β0. Then, by putting the relations in equation (29) of Lemma 2, the relation can be written as(33)ATVm=VmTm+vm+1μm+1−fhm+1,memTHm−1.
The residual of CARE (5) can be written as(34)ℛX=ATXE+ETXA−ETXℬℬTXE+CTC.
Consider the approximate solution using the low-rank solutionX^ as X=VmX^VmT and equation (30) in Lemma 3, then applying (33) in (34), we get(35)ℛXm=VmTmX^mVmTE+vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^TmTVmT+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT−ETVmX^VmTℬℬTVmX^VmTE+CTC=vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT+ℳTXE+ETXℳ−ETXℬℬTXE+CTC;ℳ=ATVmVmT=ETVmX^Hm−Temhm+1,mTμm+1Tvm+1T+vm+1μm+1−fhm+1,memTHm−1X^mVmTE−ETVmX^Hm−Temhm+1,mTfT+0=vm+1μm+1ETVmX^Hm−Temhm+1,mT−f0I0I0I0I0μm+1Tvm+1Thm+1,memTHm−1X^mVmTE−fT=SJST.
Thus, the proof follows from the Frobenius norm of equation (35).
Then, the relative residual with respect toβ0 can be estimated as follows:(36)ℛFRelative=ℛFβ0Tβ0F.


#### 3.1.2. Sparsity Preservation

The matrixA in (3) is in dense form, which is contradictory to the aim of the work and the rate of convergence of the converted system is slow enough. So, to bypass these drawbacks at each iteration, a shifted linear system needs to be solved for vi as(37)AT−μiETvi=Vi−1,orJ1T−μiE1TJ3TJ2TJ4TviΓ=Vi−10.Here,Γ=−J4−1J3vi is the truncated term. Linear system (37) is higher dimensional but sparse and can be solved by the conventional sparse-direct solvers very efficiently [30]. To improve the consistency of the RKSM approach, explicit form of the reduced-order matrices must not be used to construct reduced-order system. The sparsity-preserving reduced-order matrices can be attained by the following way:(38)E^=VTE1V,A^=VTJ1V−VTJ2J4−1J3V,ℬ^=VTB1−VTJ2J4−1B2,C^=C1V−C2J4−1J3V.


#### 3.1.3. Treatment for the Unstable Systems

If the system is unstable, a Bernoulli stabilization is required through an initial-feedback matrixK0 to estimate Af=A−ℬK0, and the matrix A needs to be replaced [31]. Then, system (3) and the CARE (5) need to be redefined as(39)Ex˙t=Afxt+ℬut,yt=Cxt+Dut,(40)AfTXE+ETXAf−ETXℬℬTXE+CTC=0.Then, for every iteration,K=ℬ^TX^VTE needs to be updated by the solution X^, which is the low-rank approximant of the solution X of (40). For this, the rational Krylov subspace for the projector X needs to be redefined as(41)Km=spanCT,AfT−μ1ET−1CT,…,∏i=1mAfT−μiET−1CT.For the stabilized system using the initial-feedback matrixK0, expressions (37) can be written as(42)AfT−μiETvi=Vi−1,orJ1−B1K0T−μiE1TJ3−B2K0TJ2TJ4TviΓ=Vi−10.To evaluate the shifted linear systems, explicit inversion ofA−ℬK should be avoided in practice, instead the Sherman–Morrison–Woodbury formula needs to be used as follows:(43)A−ℬK−1=A−1+A−1ℬI−KA−1ℬ−1KA−1.


#### 3.1.4. Estimation of the Optimal Feedback Matrix

The low-rank solutionX^ is symmetric and positive definite and can be factorized as X^=YYT. The original solution can be approximated as X=VX^VT=VYVYT. Finally, the desired low-rank factored solution Z=VY of the CARE (5) will be stored, and the optimal feedback matrix Ko=ℬTXE=ℬTZZTE can be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in Algorithm 4.Algorithm 4: Updated RKSM (sparsity-preserving).
Input:E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), imax (number of iterations), and μi (initial shifts).Output: Low-rank factored solution Z such that X≈ZZT optimal feedback matrix Ko.(1)
ComputeC1V−C2J4−1J3VT=Q0R0 (QR factorization).(2)
ChooseV0=Q0.(3)
ChooseK=K0.(4)
whilenot converged or i≤imaxdo(5)
Solve the linear system (42) for vm+1.(6)
Compute adaptive shift parameters for the next iterations (if store is empty).(7)
Using Arnoldi algorithm, orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vmv^m+1.(8)
AssumingE^, A^, ℬ^, and C^, for X^ solve the reduced-order Riccati equation (28).(9)
UpdateK=ℬ^TX^Vm+1E1.(10)
ComputeℛFRelative for convergence.(11)
end while(12)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(13)
For negligible eigenvalues, truncateΛ2 and compute Z=Vm+1T1Λ11/2.(14)
Compute the optimal feedback matrixKo=B1−J2J4−1B2TZZTE1.


### 3.2. Updated LRCF-ADI-Based Modified Kleinman–Newton Method

In each iteration of Algorithm3, the generalized CALE needs to be solved for once, and there are several techniques available to do it. In his Ph.D. thesis, Kuerschner discussed low-rank ADI approaches (Algorithm 3.2 in Chapter 3 and Algorithm 6.2 in Chapter 6) for solving generalized CALE derived from generalized CARE in the iterative loops of Kleinman–Newton algorithm [21]. Now, we need to implement above mechanisms for the index-1 descriptor system in the sparse form. For the adjustment, some modifications are required as given below.Here, the initial stage of the computation of the shift parametersμi is the same as narrated in previous section. But due to the structure of the subspace for the ADI technique, the dominant eigenvalues with the negative real part are to store, and those having the positive real part are to shift to the negative complex plane by the mirror imaging process.


#### 3.2.1. Convergence Criteria and Recurrence Relations

The computation of residuals of the ADI iterations can be achieved using the simplified technique implementing the adjustable shift parameters. This approach will be efficient for time dealing and memory allocation. Following lemmas represent some important properties of the ADI iterates.Lemma 4 (convergence of ADI iterates [32]).
LetXi be the solution at the i-th iteration of the CALE as follows:(44)ATXE+ETXA=−CTC.ConsiderXik and Xik+1 are the successive iterates of the ADI method. Then, for all μi∈ℂ−, the relation holds(45)Xik+1−Xi=∏j=1kAk,μjXi0−Xi∏j=1kAk,μjT,where Ak,μj=AT−μ¯jETAT+μjET−1.Lemma 5 (residual of the ADI iterates [32]).
LetXik an iterate of the CALE (44) by the ADI method. Then, considering Ak,μj=AT−μ¯jETAT+μjET−1 for all μi∈ℂ−, the residuals at Xik have the form(46)Rik+1=ATXikE+ETXikA+CTC=∏j=1kAk,μjATXi0E+ETXi0A+CTC∏j=1kAk,μjT.Theorem 2 (residual factor of the ADI iterations).
The residual for CALE (44) at i-th iteration of ADI method is of rank at most m and is given by(47)ℛXi=ATXiE+ETXiA+CTC=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT, for all μi∈ℂ−.Proof.
ConsiderX is the solution of CALE (44) and Xi is its i-th iterate by ADI method. Then, the residual of the i-th iteration in terms of Xi can be written as(48)ℛXi=ATXiE+ETXiA+CTC.
For allμi∈ℂ−, the Stein’s equation is equivalent to CALE (44), and it can be written as(49)X=AT−μ¯iETAT+μiET−1XiAT−μ¯iETTAT+μiET−T−2μiAT+μiET−TCTCAT+μiET−1.
Using the Stein’s equation (49) in Lemma 4, we can find some initial iterates Xi0=0 for which the residuals ℛik+1≥0. Then, according to Lemma 5, residual (48) can be defined as(50)ℛXi=∏i=1mAT−μ¯iETAT+μiET−1CTC∏i=1mAT−μ¯iETAT+μiET−1T=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT.
Thus, forWiWiT≤τ, the ADI iterations in the LRCF-ADI algorithm needs to be stopped, where τ is a given margin of tolerance. Then, with the residual factor relation for Vi can be derived as(51)Vi=AT+μiET−1Wi−1.
Then, using (51) the residual factor Wi can be derived in a recursive form as(52)Wi=Wi−1−2ReμiETVi.
In case of real setting,μi+1=μ¯i needs to be considered to find the following form:(53)Wi+1=Wi−1+2γi2ETReVi+δiImVi.
The summary of above techniques is given in Algorithm5.Algorithm 5: Updated LRCF-ADI (real version).
Input:E,A,C, τ (tolerance), imax (number of iterations), and shift parameters μjj=1imax.Output: Low-rank Cholesky-factor Z such that X≈ZZT.(1)
ConsiderW0=CT, Z0= and i=1.(2)
whileWiWiT≥τ or i≤imaxdo(3)
SolveVi=AT+μiET−1Wi−1.(4)
ifImμi=0then(5)
UpdateZi=Zi−1−2μiVi.(6)
ComputeWi=Wi−1−2μiETVi(7)
else(8)
Assumeγi=−2Reμi and δi=Reμi/Imμi.(9)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.(10)
ComputeWi+1=Wi−1+2γi2ETReVi+δiImVi.(11)
i=i+1(12)
end if(13)
i=i+1(14)
end for


#### 3.2.2. Adjustment for the Unstable Sparse Systems

For the unstable index-1 descriptor system, initial feedback matrixK0 needs to be introduced, and instead of A˜i,E, corresponding shift parameters are needed to be computed from eigenpair A˜i−ℬK0,E. The sparse form of the eigenpair can be structured as(54)A˜i−ℬK0,E=A−ℬK0−ℬℬTZiZiTE,E=J1−B1K0−B1B1TZiZiTE1J2J3−B2K0−B2B1TZiZiTE1J4,E1000.To findVji, in each ADI (inner) iteration, a shifted linear system needs to be solved as(55)A˜i−ℬK0+μjiETVji=Wj−1i,orA−ℬK0−ℬℬTZiZiTE+μjiETVji=Wj−1i.Thus,Vji can be obtained from the sparse form of the shifted linear system structured as(56)J1−B1K0−B1B1TZiZiTE1T+μjiE1TJ3−B1K0−B2B1TZiZiTE1TJ2TJ4TVjiΓ=C1TE1TZiZiTB1C2T0.


#### 3.2.3. Estimation of the Optimal Feedback Matrix

The feedback matrixKi=ℬTXiE=B1−J2J4−1B2TZiZiTE1 needs to be computed in each ADI (inner) iteration, and the optimal feedback matrix Ko=Kimax needs to be stored after the final Newton (outer) iteration. The summary of the modified method is given in Algorithm 6.Algorithm 6: Modified KN-LRCF-ADI (sparsity-preserving).
Input: E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), τ (tolerance), and imax (number of iterations).Output: Low-rank Cholesky-factor Z such that X≈ZZT and optimal feedback matrix Ko.(1)
fori←1 to imaxdo(2)
ChooseZ0i=, K0i=K0 and j=0.(3)
AssumeW0i=C1TKi−1TC2T0.(4)
Compute adaptive shift parametersμ1i,…,μji from the eigenpair defined in (54).(5)
whileWji2>τW0i2do(6)
j=j+1(7)
Solve the linear system (56) for Vji.(8)
ifImμi=0then(9)
UpdateZji=Zj−1i−2μjiVji.(10)
ComputeWji=Wj−1i−2μjiE1TVji.(11)
ComputeKji=Kj−1i−2μjiB1−J2J4−1B2TVjiVjiTE1.(12)
else(13)
Assumeγji=−2Reμji, δji=Reμji/Imμji, and βji=ReVji+δjiImVji.(14)
ComputeZd=γjiβjiγjiδji2+1ImVji.(15)
UpdateZj+1i=Zj−1iZd.(16)
ComputeWj+1i=Wj−1i+2γji2E1Tβji.(17)
ComputeKj+1i=Kj−1i+B1−J2J4−1B2TZdZdTE1.(18)
j=j+1(19)
end if(20)
end while(21)
UpdateZi=Zji and Ki=Kji(22)
end for


### 3.3. Computation of the Optimally Stabilized System and the Optimal Control

The optimal feedback matrixKo=ℬTXE can be achieved by the feasible solution X of Riccati equation (5). Then, applying As=A−ℬKo, optimally stabilized LTI continuous-time system can be written as (7). To preserve the structure of the system, it needs to back to the original structure (2), and for this, the submatrices J1, J3 are replaced by J1−B1Ko, J3−B2Ko, respectively.Finally, the desired optimal controluo=−Koxt of the targeted power system models can be computed.


## 3.1. Updated Rational Krylov Subspace Method

Let us consider LTI continuous-time system (3) and introduce an orthogonal projector V spanned by the m-dimensional rational Krylov subspace for a set of shift parameters μi∈ℂ+;i=1,2,…,m that is defined as(27)Km=spanCT,AT−μ1ET−1CT,…,∏i=1mAT−μiET−1CT.To find the shift parametersμi, the eigenvalues of the matrix pencil A,E are evaluated and sorted numerically. The required number of dominant eigenvalues are to store for the upcoming iterations. If some of the dominant eigenvalues have the negative real part, they are to shift to the positive complex plane by the mirror imaging process.Again, consider CARE (5) and apply the Galerkin condition on it. Then, after the simplification by matrix algebra, a low-rank CARE can be achieved as(28)A^TX^E^+E^TX^A^−E^TX^ℬ^ℬ^TX^E^+C^TC^=0,where X^=VTXV, E^=VTEV, A^=VTAV, ℬ^=VTℬ, and C^=CV. Equation (28) is a low-rank CARE and can be solved by MATLAB care command or any existing methods, such as Schur-decomposition method.For the quick and smooth convergence of the proposed algorithm, adjustable shift selection is crucial, and we are adopting the adaptive shift approach for index-1 descriptor systems [27]. This process required to be recursive, and in each step, the subspace to all the projector V generated with the current set of shifts will be extended.


### 3.1.1. Stopping Criteria and Related Theorem

Arnoldi relation is a very essential tool for the computation of residual of the RKSM iterations. To avoid extra matrix-vector multiplies per iteration, the computation of projected matrixTm can be performed more efficiently than the explicit product Tm=VmTAVm. To find the stopping criteria, we have to consider the following lemmas.Lemma 2 (Arnoldi relation [28]).
LetKm=spanVm be the rational Krylov subspace with the shift parameters μi∈ℂ+;i=1,2,…,m. Then, Km satisfies the Arnoldi relation as follows:(29)ATVm=VmTm+v^m+1gmT;VmTVm=Im,where v^m+1β1=vm+1μm+1−In−VmVmTATvm+1 is the QR decomposition of the right-hand side matrix with gmT=β1hm+1,memTHm−1.Lemma 3 (building the projected matrix [29]).
Let the column vectors ofVm be an orthonormal basis of the rational Krylov subspace Km with the block diagonal matrix Dm=diagμ1ET,μ2ET,…,μmET, where μ1,μ2,…,μm is the set of shift parameters used in the algorithm. Then, for the projected matrix Tm=VmTAVm, the following relation holds:(30)Tm=Im+VmTETVmHmDm−VmTATvm+1hm+1,memTHm−1.Theorem 1 (residual of the RKSM iterations).
LetVm be the orthogonal projector spanned by the rational Krylov subspace Km and X=VmX^VmT is the solution of the CARE (5) using the low-rank solution X^. Then, the residual of m-th iteration can be computed as(31)ℛF=SJSTF;J=0I0I0I0I0,where ⋅F denotes the Frobenius norm, and S is a block upper triangular matrix in the QR factorization of the matrix U defined as(32)U=vm+1μm+1ETVmX^Hm−Temhm+1,mT−In−VmVmTATvm+1.Proof.
Assumef=In−VmVmTATvm+1 and consider the reduced QR factorization CT=V0β0. Then, by putting the relations in equation (29) of Lemma 2, the relation can be written as(33)ATVm=VmTm+vm+1μm+1−fhm+1,memTHm−1.
The residual of CARE (5) can be written as(34)ℛX=ATXE+ETXA−ETXℬℬTXE+CTC.
Consider the approximate solution using the low-rank solutionX^ as X=VmX^VmT and equation (30) in Lemma 3, then applying (33) in (34), we get(35)ℛXm=VmTmX^mVmTE+vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^TmTVmT+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT−ETVmX^VmTℬℬTVmX^VmTE+CTC=vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT+ℳTXE+ETXℳ−ETXℬℬTXE+CTC;ℳ=ATVmVmT=ETVmX^Hm−Temhm+1,mTμm+1Tvm+1T+vm+1μm+1−fhm+1,memTHm−1X^mVmTE−ETVmX^Hm−Temhm+1,mTfT+0=vm+1μm+1ETVmX^Hm−Temhm+1,mT−f0I0I0I0I0μm+1Tvm+1Thm+1,memTHm−1X^mVmTE−fT=SJST.
Thus, the proof follows from the Frobenius norm of equation (35).
Then, the relative residual with respect toβ0 can be estimated as follows:(36)ℛFRelative=ℛFβ0Tβ0F.


### 3.1.2. Sparsity Preservation

The matrixA in (3) is in dense form, which is contradictory to the aim of the work and the rate of convergence of the converted system is slow enough. So, to bypass these drawbacks at each iteration, a shifted linear system needs to be solved for vi as(37)AT−μiETvi=Vi−1,orJ1T−μiE1TJ3TJ2TJ4TviΓ=Vi−10.Here,Γ=−J4−1J3vi is the truncated term. Linear system (37) is higher dimensional but sparse and can be solved by the conventional sparse-direct solvers very efficiently [30]. To improve the consistency of the RKSM approach, explicit form of the reduced-order matrices must not be used to construct reduced-order system. The sparsity-preserving reduced-order matrices can be attained by the following way:(38)E^=VTE1V,A^=VTJ1V−VTJ2J4−1J3V,ℬ^=VTB1−VTJ2J4−1B2,C^=C1V−C2J4−1J3V.


### 3.1.3. Treatment for the Unstable Systems

If the system is unstable, a Bernoulli stabilization is required through an initial-feedback matrixK0 to estimate Af=A−ℬK0, and the matrix A needs to be replaced [31]. Then, system (3) and the CARE (5) need to be redefined as(39)Ex˙t=Afxt+ℬut,yt=Cxt+Dut,(40)AfTXE+ETXAf−ETXℬℬTXE+CTC=0.Then, for every iteration,K=ℬ^TX^VTE needs to be updated by the solution X^, which is the low-rank approximant of the solution X of (40). For this, the rational Krylov subspace for the projector X needs to be redefined as(41)Km=spanCT,AfT−μ1ET−1CT,…,∏i=1mAfT−μiET−1CT.For the stabilized system using the initial-feedback matrixK0, expressions (37) can be written as(42)AfT−μiETvi=Vi−1,orJ1−B1K0T−μiE1TJ3−B2K0TJ2TJ4TviΓ=Vi−10.To evaluate the shifted linear systems, explicit inversion ofA−ℬK should be avoided in practice, instead the Sherman–Morrison–Woodbury formula needs to be used as follows:(43)A−ℬK−1=A−1+A−1ℬI−KA−1ℬ−1KA−1.


### 3.1.4. Estimation of the Optimal Feedback Matrix

The low-rank solutionX^ is symmetric and positive definite and can be factorized as X^=YYT. The original solution can be approximated as X=VX^VT=VYVYT. Finally, the desired low-rank factored solution Z=VY of the CARE (5) will be stored, and the optimal feedback matrix Ko=ℬTXE=ℬTZZTE can be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in Algorithm 4.Algorithm 4: Updated RKSM (sparsity-preserving).
Input:E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), imax (number of iterations), and μi (initial shifts).Output: Low-rank factored solution Z such that X≈ZZT optimal feedback matrix Ko.(1)
ComputeC1V−C2J4−1J3VT=Q0R0 (QR factorization).(2)
ChooseV0=Q0.(3)
ChooseK=K0.(4)
whilenot converged or i≤imaxdo(5)
Solve the linear system (42) for vm+1.(6)
Compute adaptive shift parameters for the next iterations (if store is empty).(7)
Using Arnoldi algorithm, orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vmv^m+1.(8)
AssumingE^, A^, ℬ^, and C^, for X^ solve the reduced-order Riccati equation (28).(9)
UpdateK=ℬ^TX^Vm+1E1.(10)
ComputeℛFRelative for convergence.(11)
end while(12)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(13)
For negligible eigenvalues, truncateΛ2 and compute Z=Vm+1T1Λ11/2.(14)
Compute the optimal feedback matrixKo=B1−J2J4−1B2TZZTE1.


## 3.1.1. Stopping Criteria and Related Theorem

Arnoldi relation is a very essential tool for the computation of residual of the RKSM iterations. To avoid extra matrix-vector multiplies per iteration, the computation of projected matrixTm can be performed more efficiently than the explicit product Tm=VmTAVm. To find the stopping criteria, we have to consider the following lemmas.Lemma 2 (Arnoldi relation [28]).
LetKm=spanVm be the rational Krylov subspace with the shift parameters μi∈ℂ+;i=1,2,…,m. Then, Km satisfies the Arnoldi relation as follows:(29)ATVm=VmTm+v^m+1gmT;VmTVm=Im,where v^m+1β1=vm+1μm+1−In−VmVmTATvm+1 is the QR decomposition of the right-hand side matrix with gmT=β1hm+1,memTHm−1.Lemma 3 (building the projected matrix [29]).
Let the column vectors ofVm be an orthonormal basis of the rational Krylov subspace Km with the block diagonal matrix Dm=diagμ1ET,μ2ET,…,μmET, where μ1,μ2,…,μm is the set of shift parameters used in the algorithm. Then, for the projected matrix Tm=VmTAVm, the following relation holds:(30)Tm=Im+VmTETVmHmDm−VmTATvm+1hm+1,memTHm−1.Theorem 1 (residual of the RKSM iterations).
LetVm be the orthogonal projector spanned by the rational Krylov subspace Km and X=VmX^VmT is the solution of the CARE (5) using the low-rank solution X^. Then, the residual of m-th iteration can be computed as(31)ℛF=SJSTF;J=0I0I0I0I0,where ⋅F denotes the Frobenius norm, and S is a block upper triangular matrix in the QR factorization of the matrix U defined as(32)U=vm+1μm+1ETVmX^Hm−Temhm+1,mT−In−VmVmTATvm+1.Proof.
Assumef=In−VmVmTATvm+1 and consider the reduced QR factorization CT=V0β0. Then, by putting the relations in equation (29) of Lemma 2, the relation can be written as(33)ATVm=VmTm+vm+1μm+1−fhm+1,memTHm−1.
The residual of CARE (5) can be written as(34)ℛX=ATXE+ETXA−ETXℬℬTXE+CTC.
Consider the approximate solution using the low-rank solutionX^ as X=VmX^VmT and equation (30) in Lemma 3, then applying (33) in (34), we get(35)ℛXm=VmTmX^mVmTE+vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^TmTVmT+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT−ETVmX^VmTℬℬTVmX^VmTE+CTC=vm+1μm+1−fhm+1,memTHm−1X^mVmTE+ETVmX^Hm−Temhm+1,mTvm+1μm+1−fT+ℳTXE+ETXℳ−ETXℬℬTXE+CTC;ℳ=ATVmVmT=ETVmX^Hm−Temhm+1,mTμm+1Tvm+1T+vm+1μm+1−fhm+1,memTHm−1X^mVmTE−ETVmX^Hm−Temhm+1,mTfT+0=vm+1μm+1ETVmX^Hm−Temhm+1,mT−f0I0I0I0I0μm+1Tvm+1Thm+1,memTHm−1X^mVmTE−fT=SJST.
Thus, the proof follows from the Frobenius norm of equation (35).
Then, the relative residual with respect toβ0 can be estimated as follows:(36)ℛFRelative=ℛFβ0Tβ0F.


## 3.1.2. Sparsity Preservation

The matrixA in (3) is in dense form, which is contradictory to the aim of the work and the rate of convergence of the converted system is slow enough. So, to bypass these drawbacks at each iteration, a shifted linear system needs to be solved for vi as(37)AT−μiETvi=Vi−1,orJ1T−μiE1TJ3TJ2TJ4TviΓ=Vi−10.Here,Γ=−J4−1J3vi is the truncated term. Linear system (37) is higher dimensional but sparse and can be solved by the conventional sparse-direct solvers very efficiently [30]. To improve the consistency of the RKSM approach, explicit form of the reduced-order matrices must not be used to construct reduced-order system. The sparsity-preserving reduced-order matrices can be attained by the following way:(38)E^=VTE1V,A^=VTJ1V−VTJ2J4−1J3V,ℬ^=VTB1−VTJ2J4−1B2,C^=C1V−C2J4−1J3V.


## 3.1.3. Treatment for the Unstable Systems

If the system is unstable, a Bernoulli stabilization is required through an initial-feedback matrixK0 to estimate Af=A−ℬK0, and the matrix A needs to be replaced [31]. Then, system (3) and the CARE (5) need to be redefined as(39)Ex˙t=Afxt+ℬut,yt=Cxt+Dut,(40)AfTXE+ETXAf−ETXℬℬTXE+CTC=0.Then, for every iteration,K=ℬ^TX^VTE needs to be updated by the solution X^, which is the low-rank approximant of the solution X of (40). For this, the rational Krylov subspace for the projector X needs to be redefined as(41)Km=spanCT,AfT−μ1ET−1CT,…,∏i=1mAfT−μiET−1CT.For the stabilized system using the initial-feedback matrixK0, expressions (37) can be written as(42)AfT−μiETvi=Vi−1,orJ1−B1K0T−μiE1TJ3−B2K0TJ2TJ4TviΓ=Vi−10.To evaluate the shifted linear systems, explicit inversion ofA−ℬK should be avoided in practice, instead the Sherman–Morrison–Woodbury formula needs to be used as follows:(43)A−ℬK−1=A−1+A−1ℬI−KA−1ℬ−1KA−1.


## 3.1.4. Estimation of the Optimal Feedback Matrix

The low-rank solutionX^ is symmetric and positive definite and can be factorized as X^=YYT. The original solution can be approximated as X=VX^VT=VYVYT. Finally, the desired low-rank factored solution Z=VY of the CARE (5) will be stored, and the optimal feedback matrix Ko=ℬTXE=ℬTZZTE can be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in Algorithm 4.Algorithm 4: Updated RKSM (sparsity-preserving).
Input:E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), imax (number of iterations), and μi (initial shifts).Output: Low-rank factored solution Z such that X≈ZZT optimal feedback matrix Ko.(1)
ComputeC1V−C2J4−1J3VT=Q0R0 (QR factorization).(2)
ChooseV0=Q0.(3)
ChooseK=K0.(4)
whilenot converged or i≤imaxdo(5)
Solve the linear system (42) for vm+1.(6)
Compute adaptive shift parameters for the next iterations (if store is empty).(7)
Using Arnoldi algorithm, orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vmv^m+1.(8)
AssumingE^, A^, ℬ^, and C^, for X^ solve the reduced-order Riccati equation (28).(9)
UpdateK=ℬ^TX^Vm+1E1.(10)
ComputeℛFRelative for convergence.(11)
end while(12)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(13)
For negligible eigenvalues, truncateΛ2 and compute Z=Vm+1T1Λ11/2.(14)
Compute the optimal feedback matrixKo=B1−J2J4−1B2TZZTE1.


## 3.2. Updated LRCF-ADI-Based Modified Kleinman–Newton Method

In each iteration of Algorithm3, the generalized CALE needs to be solved for once, and there are several techniques available to do it. In his Ph.D. thesis, Kuerschner discussed low-rank ADI approaches (Algorithm 3.2 in Chapter 3 and Algorithm 6.2 in Chapter 6) for solving generalized CALE derived from generalized CARE in the iterative loops of Kleinman–Newton algorithm [21]. Now, we need to implement above mechanisms for the index-1 descriptor system in the sparse form. For the adjustment, some modifications are required as given below.Here, the initial stage of the computation of the shift parametersμi is the same as narrated in previous section. But due to the structure of the subspace for the ADI technique, the dominant eigenvalues with the negative real part are to store, and those having the positive real part are to shift to the negative complex plane by the mirror imaging process.


### 3.2.1. Convergence Criteria and Recurrence Relations

The computation of residuals of the ADI iterations can be achieved using the simplified technique implementing the adjustable shift parameters. This approach will be efficient for time dealing and memory allocation. Following lemmas represent some important properties of the ADI iterates.Lemma 4 (convergence of ADI iterates [32]).
LetXi be the solution at the i-th iteration of the CALE as follows:(44)ATXE+ETXA=−CTC.ConsiderXik and Xik+1 are the successive iterates of the ADI method. Then, for all μi∈ℂ−, the relation holds(45)Xik+1−Xi=∏j=1kAk,μjXi0−Xi∏j=1kAk,μjT,where Ak,μj=AT−μ¯jETAT+μjET−1.Lemma 5 (residual of the ADI iterates [32]).
LetXik an iterate of the CALE (44) by the ADI method. Then, considering Ak,μj=AT−μ¯jETAT+μjET−1 for all μi∈ℂ−, the residuals at Xik have the form(46)Rik+1=ATXikE+ETXikA+CTC=∏j=1kAk,μjATXi0E+ETXi0A+CTC∏j=1kAk,μjT.Theorem 2 (residual factor of the ADI iterations).
The residual for CALE (44) at i-th iteration of ADI method is of rank at most m and is given by(47)ℛXi=ATXiE+ETXiA+CTC=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT, for all μi∈ℂ−.Proof.
ConsiderX is the solution of CALE (44) and Xi is its i-th iterate by ADI method. Then, the residual of the i-th iteration in terms of Xi can be written as(48)ℛXi=ATXiE+ETXiA+CTC.
For allμi∈ℂ−, the Stein’s equation is equivalent to CALE (44), and it can be written as(49)X=AT−μ¯iETAT+μiET−1XiAT−μ¯iETTAT+μiET−T−2μiAT+μiET−TCTCAT+μiET−1.
Using the Stein’s equation (49) in Lemma 4, we can find some initial iterates Xi0=0 for which the residuals ℛik+1≥0. Then, according to Lemma 5, residual (48) can be defined as(50)ℛXi=∏i=1mAT−μ¯iETAT+μiET−1CTC∏i=1mAT−μ¯iETAT+μiET−1T=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT.
Thus, forWiWiT≤τ, the ADI iterations in the LRCF-ADI algorithm needs to be stopped, where τ is a given margin of tolerance. Then, with the residual factor relation for Vi can be derived as(51)Vi=AT+μiET−1Wi−1.
Then, using (51) the residual factor Wi can be derived in a recursive form as(52)Wi=Wi−1−2ReμiETVi.
In case of real setting,μi+1=μ¯i needs to be considered to find the following form:(53)Wi+1=Wi−1+2γi2ETReVi+δiImVi.
The summary of above techniques is given in Algorithm5.Algorithm 5: Updated LRCF-ADI (real version).
Input:E,A,C, τ (tolerance), imax (number of iterations), and shift parameters μjj=1imax.Output: Low-rank Cholesky-factor Z such that X≈ZZT.(1)
ConsiderW0=CT, Z0= and i=1.(2)
whileWiWiT≥τ or i≤imaxdo(3)
SolveVi=AT+μiET−1Wi−1.(4)
ifImμi=0then(5)
UpdateZi=Zi−1−2μiVi.(6)
ComputeWi=Wi−1−2μiETVi(7)
else(8)
Assumeγi=−2Reμi and δi=Reμi/Imμi.(9)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.(10)
ComputeWi+1=Wi−1+2γi2ETReVi+δiImVi.(11)
i=i+1(12)
end if(13)
i=i+1(14)
end for


### 3.2.2. Adjustment for the Unstable Sparse Systems

For the unstable index-1 descriptor system, initial feedback matrixK0 needs to be introduced, and instead of A˜i,E, corresponding shift parameters are needed to be computed from eigenpair A˜i−ℬK0,E. The sparse form of the eigenpair can be structured as(54)A˜i−ℬK0,E=A−ℬK0−ℬℬTZiZiTE,E=J1−B1K0−B1B1TZiZiTE1J2J3−B2K0−B2B1TZiZiTE1J4,E1000.To findVji, in each ADI (inner) iteration, a shifted linear system needs to be solved as(55)A˜i−ℬK0+μjiETVji=Wj−1i,orA−ℬK0−ℬℬTZiZiTE+μjiETVji=Wj−1i.Thus,Vji can be obtained from the sparse form of the shifted linear system structured as(56)J1−B1K0−B1B1TZiZiTE1T+μjiE1TJ3−B1K0−B2B1TZiZiTE1TJ2TJ4TVjiΓ=C1TE1TZiZiTB1C2T0.


### 3.2.3. Estimation of the Optimal Feedback Matrix

The feedback matrixKi=ℬTXiE=B1−J2J4−1B2TZiZiTE1 needs to be computed in each ADI (inner) iteration, and the optimal feedback matrix Ko=Kimax needs to be stored after the final Newton (outer) iteration. The summary of the modified method is given in Algorithm 6.Algorithm 6: Modified KN-LRCF-ADI (sparsity-preserving).
Input: E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), τ (tolerance), and imax (number of iterations).Output: Low-rank Cholesky-factor Z such that X≈ZZT and optimal feedback matrix Ko.(1)
fori←1 to imaxdo(2)
ChooseZ0i=, K0i=K0 and j=0.(3)
AssumeW0i=C1TKi−1TC2T0.(4)
Compute adaptive shift parametersμ1i,…,μji from the eigenpair defined in (54).(5)
whileWji2>τW0i2do(6)
j=j+1(7)
Solve the linear system (56) for Vji.(8)
ifImμi=0then(9)
UpdateZji=Zj−1i−2μjiVji.(10)
ComputeWji=Wj−1i−2μjiE1TVji.(11)
ComputeKji=Kj−1i−2μjiB1−J2J4−1B2TVjiVjiTE1.(12)
else(13)
Assumeγji=−2Reμji, δji=Reμji/Imμji, and βji=ReVji+δjiImVji.(14)
ComputeZd=γjiβjiγjiδji2+1ImVji.(15)
UpdateZj+1i=Zj−1iZd.(16)
ComputeWj+1i=Wj−1i+2γji2E1Tβji.(17)
ComputeKj+1i=Kj−1i+B1−J2J4−1B2TZdZdTE1.(18)
j=j+1(19)
end if(20)
end while(21)
UpdateZi=Zji and Ki=Kji(22)
end for


## 3.2.1. Convergence Criteria and Recurrence Relations

The computation of residuals of the ADI iterations can be achieved using the simplified technique implementing the adjustable shift parameters. This approach will be efficient for time dealing and memory allocation. Following lemmas represent some important properties of the ADI iterates.Lemma 4 (convergence of ADI iterates [32]).
LetXi be the solution at the i-th iteration of the CALE as follows:(44)ATXE+ETXA=−CTC.ConsiderXik and Xik+1 are the successive iterates of the ADI method. Then, for all μi∈ℂ−, the relation holds(45)Xik+1−Xi=∏j=1kAk,μjXi0−Xi∏j=1kAk,μjT,where Ak,μj=AT−μ¯jETAT+μjET−1.Lemma 5 (residual of the ADI iterates [32]).
LetXik an iterate of the CALE (44) by the ADI method. Then, considering Ak,μj=AT−μ¯jETAT+μjET−1 for all μi∈ℂ−, the residuals at Xik have the form(46)Rik+1=ATXikE+ETXikA+CTC=∏j=1kAk,μjATXi0E+ETXi0A+CTC∏j=1kAk,μjT.Theorem 2 (residual factor of the ADI iterations).
The residual for CALE (44) at i-th iteration of ADI method is of rank at most m and is given by(47)ℛXi=ATXiE+ETXiA+CTC=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT, for all μi∈ℂ−.Proof.
ConsiderX is the solution of CALE (44) and Xi is its i-th iterate by ADI method. Then, the residual of the i-th iteration in terms of Xi can be written as(48)ℛXi=ATXiE+ETXiA+CTC.
For allμi∈ℂ−, the Stein’s equation is equivalent to CALE (44), and it can be written as(49)X=AT−μ¯iETAT+μiET−1XiAT−μ¯iETTAT+μiET−T−2μiAT+μiET−TCTCAT+μiET−1.
Using the Stein’s equation (49) in Lemma 4, we can find some initial iterates Xi0=0 for which the residuals ℛik+1≥0. Then, according to Lemma 5, residual (48) can be defined as(50)ℛXi=∏i=1mAT−μ¯iETAT+μiET−1CTC∏i=1mAT−μ¯iETAT+μiET−1T=WiWiT,with Wi=∏i=1mAT−μ¯iETAT+μiET−1CT.
Thus, forWiWiT≤τ, the ADI iterations in the LRCF-ADI algorithm needs to be stopped, where τ is a given margin of tolerance. Then, with the residual factor relation for Vi can be derived as(51)Vi=AT+μiET−1Wi−1.
Then, using (51) the residual factor Wi can be derived in a recursive form as(52)Wi=Wi−1−2ReμiETVi.
In case of real setting,μi+1=μ¯i needs to be considered to find the following form:(53)Wi+1=Wi−1+2γi2ETReVi+δiImVi.
The summary of above techniques is given in Algorithm5.Algorithm 5: Updated LRCF-ADI (real version).
Input:E,A,C, τ (tolerance), imax (number of iterations), and shift parameters μjj=1imax.Output: Low-rank Cholesky-factor Z such that X≈ZZT.(1)
ConsiderW0=CT, Z0= and i=1.(2)
whileWiWiT≥τ or i≤imaxdo(3)
SolveVi=AT+μiET−1Wi−1.(4)
ifImμi=0then(5)
UpdateZi=Zi−1−2μiVi.(6)
ComputeWi=Wi−1−2μiETVi(7)
else(8)
Assumeγi=−2Reμi and δi=Reμi/Imμi.(9)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi.(10)
ComputeWi+1=Wi−1+2γi2ETReVi+δiImVi.(11)
i=i+1(12)
end if(13)
i=i+1(14)
end for


## 3.2.2. Adjustment for the Unstable Sparse Systems

For the unstable index-1 descriptor system, initial feedback matrixK0 needs to be introduced, and instead of A˜i,E, corresponding shift parameters are needed to be computed from eigenpair A˜i−ℬK0,E. The sparse form of the eigenpair can be structured as(54)A˜i−ℬK0,E=A−ℬK0−ℬℬTZiZiTE,E=J1−B1K0−B1B1TZiZiTE1J2J3−B2K0−B2B1TZiZiTE1J4,E1000.To findVji, in each ADI (inner) iteration, a shifted linear system needs to be solved as(55)A˜i−ℬK0+μjiETVji=Wj−1i,orA−ℬK0−ℬℬTZiZiTE+μjiETVji=Wj−1i.Thus,Vji can be obtained from the sparse form of the shifted linear system structured as(56)J1−B1K0−B1B1TZiZiTE1T+μjiE1TJ3−B1K0−B2B1TZiZiTE1TJ2TJ4TVjiΓ=C1TE1TZiZiTB1C2T0.


## 3.2.3. Estimation of the Optimal Feedback Matrix

The feedback matrixKi=ℬTXiE=B1−J2J4−1B2TZiZiTE1 needs to be computed in each ADI (inner) iteration, and the optimal feedback matrix Ko=Kimax needs to be stored after the final Newton (outer) iteration. The summary of the modified method is given in Algorithm 6.Algorithm 6: Modified KN-LRCF-ADI (sparsity-preserving).
Input: E1,J1,J2,J3,J4,B1,B2,C1,C2,K0 (initial feedback matrix), τ (tolerance), and imax (number of iterations).Output: Low-rank Cholesky-factor Z such that X≈ZZT and optimal feedback matrix Ko.(1)
fori←1 to imaxdo(2)
ChooseZ0i=, K0i=K0 and j=0.(3)
AssumeW0i=C1TKi−1TC2T0.(4)
Compute adaptive shift parametersμ1i,…,μji from the eigenpair defined in (54).(5)
whileWji2>τW0i2do(6)
j=j+1(7)
Solve the linear system (56) for Vji.(8)
ifImμi=0then(9)
UpdateZji=Zj−1i−2μjiVji.(10)
ComputeWji=Wj−1i−2μjiE1TVji.(11)
ComputeKji=Kj−1i−2μjiB1−J2J4−1B2TVjiVjiTE1.(12)
else(13)
Assumeγji=−2Reμji, δji=Reμji/Imμji, and βji=ReVji+δjiImVji.(14)
ComputeZd=γjiβjiγjiδji2+1ImVji.(15)
UpdateZj+1i=Zj−1iZd.(16)
ComputeWj+1i=Wj−1i+2γji2E1Tβji.(17)
ComputeKj+1i=Kj−1i+B1−J2J4−1B2TZdZdTE1.(18)
j=j+1(19)
end if(20)
end while(21)
UpdateZi=Zji and Ki=Kji(22)
end for


## 3.3. Computation of the Optimally Stabilized System and the Optimal Control

The optimal feedback matrixKo=ℬTXE can be achieved by the feasible solution X of Riccati equation (5). Then, applying As=A−ℬKo, optimally stabilized LTI continuous-time system can be written as (7). To preserve the structure of the system, it needs to back to the original structure (2), and for this, the submatrices J1, J3 are replaced by J1−B1Ko, J3−B2Ko, respectively.Finally, the desired optimal controluo=−Koxt of the targeted power system models can be computed.


## 4. Numerical Results

The stability of the target models is investigated, and the unstable models are stabilized through the Riccati-based feedback stabilization process. The proposed methods are employed to find the solution of the Riccati equation arising from the BIPS models, and corresponding feedback matrices are generated for system stabilization. Also, initial Bernoulli feedback stabilization is implemented for the convenient rate of convergence. All the results have been achieved using the MATLAB 8.5.0 (R2015a) on a Windows machine having Intel-Xeon Silver 4114 CPU 2.20 GHz clock speed, 2 cores each, and 64 GB of total RAM.


### 4.1. Brazilian Interconnected Power System Models

Power system models are an essential part of engineering fields that consists of simulations based on power generations and grid networks. The computation required to analyze electrical power systems employing mathematical models utilizing real-time data. There are several applications of the power system model, i.e., electric power generation, utility transmission and distribution, railway power systems, and industrial power generation [33].The Brazilian Interconnected Power Systems (BIPS) is one of the most convenient examples of the power system models with various test systems [34]. Table 1 provides the details about the models. The detailed structure of those will be found at the webpage of Joost Rommes, where all of them are index-1 descriptor system. The models mod-606, mod-1998, mod-2476, and mod-3078 have the unstable eigenvalues, whereas the models mod-1142, mod-1450, and mod-1693 have stable eigenvalues [35]. Here, the names of the models are considered according to their number of states.Table 1
Structure of the models derived from BIPS test systems.
Test systemsDimensionsStatesAlgebraic variablesInputsOutputsBIPS987135606652944973511428593441126514509815441327516931158244BIPS07150661998130684416861247614385442112830781805044


### 4.2. Comparison of the Results Found by RKSM and KN-LRCF-ADI

The CAREs arising from the modelsmod-606, mod-1998, and mod-2476 are efficiently solved and stabilized the corresponding models by both RKSM and KN-LRCF-ADI techniques. As model mod-3078 is semistable, the computation of CARE derived from this model is not possible by LRCF-ADI techniques, but by the RKSM approach, model mod-3078 successfully stabilized, and the numerical result for model mod-3078 is investigated for RKSM only.Table2 depicts the numerical results of the stabilization process via RKSM for the unstable BIPS models and various properties of the stabilized systems are illustrated, whereas Table 3 displays the several modes of ADI techniques in KN-LRCF-ADI method for stabilizing the unstable BIPS models, including characteristics of the stabilized models. In Tables 2 and 3, we have analyzed the same features of the stabilized BIPS models, we can easily compare the efficiency and robustness of the proposed methods.Table 2
Results of RKSM-applied BIPS models.
ModelNo. of iterationsToleranceSpace dimensionNumerical rankCPU time (second)mod-60610010−104001951.81×102mod-199820010−108002651.41×103mod-247624810−109922673.06×103mod-307825710−510282953.01×103Table 3
Results of KN-LRCF-ADI-applied BIPS models.
ModelNo. of newton iterationsTotal iterationsToleranceNumerical rankCPU time (second)mod-606531110−54812.74×10254410−109537.21×1021050810−54735.12×10285310−109697.92×102mod-1998527710−56632.30×10348510−1012015.83×1031051410−54973.27×103100310−1014171.17×104mod-2476525410−54733.13×10336310−109375.75×1031036610−54413.02×10369810−108499.50×103From Tables2 and 3, it can be said that the proposed RKSM approach has quick convergence ability and occupies very small solution space to provide the efficient solution of the CAREs. In contrast, LRCF-ADI-based Kleinman–Newton has several approaches for finding the solution of CAREs, whereas most of the approaches required higher computation time. Also, there are deviations in the numerical ranks of the factored solution of CAREs in the Kleinman–Newton approaches, and in all of the cases, RKSM provides significantly better result.


### 4.3. Stabilization of Eigenvalues

Figures1–3 illustrate the eigenvalue stabilization of the models mod-606, mod-1998, and mod-2476 by both the RKSM and LRCF-ADI-based Kleinman–Newton techniques.Figure 1
Comparisons of the eigenvalues for the model mod-606. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)Figure 2
Comparisons of the eigenvalues for the model mod-1998. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)Figure 3
Comparisons of the eigenvalues for the model mod-2476. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)From Figures1–3, it has been observed that the eigenvalues of the models mod-606 and mod-2476 are stabilized very well, but it is marginal for the model mod-1998. Thus, it can be concluded that both the RKSM and LRCF-ADI-based Kleinman–Newton techniques have adequate efficiency to stabilize the unstable descriptor systems by closed-loop structure via Riccati-based feedback stabilization.Figure4 displays the applicability of the RKSM for the semistable descriptor system, whereas LRCF-ADI-based methods are ineffective in this case. Because the eigenvalues of the semistable model mod-3078 are very close to the imaginary axis, it is not possible to realize in the normal view. So, for the simulation tool capacity and visual convenience, we have considered the a very magnified view of the eigenspaces.Figure 4
Comparisons of the eigenvalues (stabilized by RKSM) for the model mod-3078.


### 4.4. Stabilization of Step Responses

The investigation of Figures5 to 12 consist the step responses for some dominant input/output relations to compare the RKSM and the LRCF-ADI-based Kleinman–Newton approaches via the system stabilization. Because the power system models are of 4×4 input-output relations, there are 16 step responses for each model. For the effective comparison, we have investigated only some graphically significant step responses.Figure 5
Comparisons of step responses for the modelmod-606 for the first-input/third-output.Figure 6
Comparisons of step responses for the modelmod-606 for the second-input/first-output.Figure 7
Comparisons of step responses for the modelmod-1998 for the third-input/second-output.Figure 8
Comparisons of step responses for the modelmod-1998 for the fourth-input/third-output.Figure 9
Comparisons of step responses for the modelmod-2476 for the first-input/third-output.Figure 10
Comparisons of step responses for the modelmod-2476 for the third-input/second-output.Figure 11
Comparisons of step responses for the modelmod-3078 for the second-input/first-output.Figure 12
Comparisons of step responses for the modelmod-3078 for the fourth-input/third-output.From Figures5 to 10, it is evident that both RKSM and LRCF-ADI-integrated Kleinman–Newton techniques are applicable for the Riccati-based feedback stabilization of unstable power system models. The graphical comparisons indicate by RKSM is suitably robust. On the other hand, though sometimes the Kleinman–Newton approach provides very good accuracy, it has some scattered behaviors.Moreover, Figures11 and 12 signify the applicability of the RKSM technique for the Riccati-based feedback stabilization for the semistable index-1 descriptor systems.


## 4.1. Brazilian Interconnected Power System Models

Power system models are an essential part of engineering fields that consists of simulations based on power generations and grid networks. The computation required to analyze electrical power systems employing mathematical models utilizing real-time data. There are several applications of the power system model, i.e., electric power generation, utility transmission and distribution, railway power systems, and industrial power generation [33].The Brazilian Interconnected Power Systems (BIPS) is one of the most convenient examples of the power system models with various test systems [34]. Table 1 provides the details about the models. The detailed structure of those will be found at the webpage of Joost Rommes, where all of them are index-1 descriptor system. The models mod-606, mod-1998, mod-2476, and mod-3078 have the unstable eigenvalues, whereas the models mod-1142, mod-1450, and mod-1693 have stable eigenvalues [35]. Here, the names of the models are considered according to their number of states.Table 1
Structure of the models derived from BIPS test systems.
Test systemsDimensionsStatesAlgebraic variablesInputsOutputsBIPS987135606652944973511428593441126514509815441327516931158244BIPS07150661998130684416861247614385442112830781805044


## 4.2. Comparison of the Results Found by RKSM and KN-LRCF-ADI

The CAREs arising from the modelsmod-606, mod-1998, and mod-2476 are efficiently solved and stabilized the corresponding models by both RKSM and KN-LRCF-ADI techniques. As model mod-3078 is semistable, the computation of CARE derived from this model is not possible by LRCF-ADI techniques, but by the RKSM approach, model mod-3078 successfully stabilized, and the numerical result for model mod-3078 is investigated for RKSM only.Table2 depicts the numerical results of the stabilization process via RKSM for the unstable BIPS models and various properties of the stabilized systems are illustrated, whereas Table 3 displays the several modes of ADI techniques in KN-LRCF-ADI method for stabilizing the unstable BIPS models, including characteristics of the stabilized models. In Tables 2 and 3, we have analyzed the same features of the stabilized BIPS models, we can easily compare the efficiency and robustness of the proposed methods.Table 2
Results of RKSM-applied BIPS models.
ModelNo. of iterationsToleranceSpace dimensionNumerical rankCPU time (second)mod-60610010−104001951.81×102mod-199820010−108002651.41×103mod-247624810−109922673.06×103mod-307825710−510282953.01×103Table 3
Results of KN-LRCF-ADI-applied BIPS models.
ModelNo. of newton iterationsTotal iterationsToleranceNumerical rankCPU time (second)mod-606531110−54812.74×10254410−109537.21×1021050810−54735.12×10285310−109697.92×102mod-1998527710−56632.30×10348510−1012015.83×1031051410−54973.27×103100310−1014171.17×104mod-2476525410−54733.13×10336310−109375.75×1031036610−54413.02×10369810−108499.50×103From Tables2 and 3, it can be said that the proposed RKSM approach has quick convergence ability and occupies very small solution space to provide the efficient solution of the CAREs. In contrast, LRCF-ADI-based Kleinman–Newton has several approaches for finding the solution of CAREs, whereas most of the approaches required higher computation time. Also, there are deviations in the numerical ranks of the factored solution of CAREs in the Kleinman–Newton approaches, and in all of the cases, RKSM provides significantly better result.


## 4.3. Stabilization of Eigenvalues

Figures1–3 illustrate the eigenvalue stabilization of the models mod-606, mod-1998, and mod-2476 by both the RKSM and LRCF-ADI-based Kleinman–Newton techniques.Figure 1
Comparisons of the eigenvalues for the model mod-606. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)Figure 2
Comparisons of the eigenvalues for the model mod-1998. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)Figure 3
Comparisons of the eigenvalues for the model mod-2476. (a) Stabilized by RKSM. (b) Stabilized by KN-LRCF-ADI.
(a)(b)From Figures1–3, it has been observed that the eigenvalues of the models mod-606 and mod-2476 are stabilized very well, but it is marginal for the model mod-1998. Thus, it can be concluded that both the RKSM and LRCF-ADI-based Kleinman–Newton techniques have adequate efficiency to stabilize the unstable descriptor systems by closed-loop structure via Riccati-based feedback stabilization.Figure4 displays the applicability of the RKSM for the semistable descriptor system, whereas LRCF-ADI-based methods are ineffective in this case. Because the eigenvalues of the semistable model mod-3078 are very close to the imaginary axis, it is not possible to realize in the normal view. So, for the simulation tool capacity and visual convenience, we have considered the a very magnified view of the eigenspaces.Figure 4
Comparisons of the eigenvalues (stabilized by RKSM) for the model mod-3078.


## 4.4. Stabilization of Step Responses

The investigation of Figures5 to 12 consist the step responses for some dominant input/output relations to compare the RKSM and the LRCF-ADI-based Kleinman–Newton approaches via the system stabilization. Because the power system models are of 4×4 input-output relations, there are 16 step responses for each model. For the effective comparison, we have investigated only some graphically significant step responses.Figure 5
Comparisons of step responses for the modelmod-606 for the first-input/third-output.Figure 6
Comparisons of step responses for the modelmod-606 for the second-input/first-output.Figure 7
Comparisons of step responses for the modelmod-1998 for the third-input/second-output.Figure 8
Comparisons of step responses for the modelmod-1998 for the fourth-input/third-output.Figure 9
Comparisons of step responses for the modelmod-2476 for the first-input/third-output.Figure 10
Comparisons of step responses for the modelmod-2476 for the third-input/second-output.Figure 11
Comparisons of step responses for the modelmod-3078 for the second-input/first-output.Figure 12
Comparisons of step responses for the modelmod-3078 for the fourth-input/third-output.From Figures5 to 10, it is evident that both RKSM and LRCF-ADI-integrated Kleinman–Newton techniques are applicable for the Riccati-based feedback stabilization of unstable power system models. The graphical comparisons indicate by RKSM is suitably robust. On the other hand, though sometimes the Kleinman–Newton approach provides very good accuracy, it has some scattered behaviors.Moreover, Figures11 and 12 signify the applicability of the RKSM technique for the Riccati-based feedback stabilization for the semistable index-1 descriptor systems.


## 5. Conclusion

From the tabular and graphical comparisons of the results of numerical computations, we have observed that by both RKSM and KN-LRCF-ADI techniques, the CAREs arising from the unstable index-1 descriptor systems are efficiently solved, and the corresponding models are properly stabilized. The semistable index-1 descriptor system is successfully stabilized through Riccati-based feedback stabilization by RKSM, whereas KN-LRCF-ADI is still not suitable for it. There are deviations of the numerical ranks of the factored solutions of CAREs in the Kleinman–Newton approaches, whereas RKSM provides significantly better results for all the cases. RKSM approach has quick convergence ability and occupies very small solution spaces to provide efficient solutions to the CAREs. In contrast, LRCF-ADI-based Kleinman–Newton method has several approaches for finding the solutions of CAREs, where almost all of the approaches required higher computation time. Riccati-based feedback stabilization for the index-1 descriptor systems by the RKSM approach is very effective and robust. Contrariwise, LRCF-ADI-based Kleinman–Newton method is slightly scattered in case of the stabilization of step responses. Thus, it can be concluded that the RKSM is suitably applicable to the unstable index-1 descriptor systems for Riccati-based feedback stabilization, and this method is more preferable to the Kleinman–Newton method in the sense of computation time and memory allocation.In this work, the MATLAB library command care is used in RKSM to find the solution of the CAREs governed from the reduced-order models and used the direct backward inversion technique for solving shifted sparse linear systems. In the future, we will try to find the self-sufficient RKSM algorithm for solving CAREs by applying the “Restarted Hessenberg Method” with a quick rate of convergence and implement the preconditioning techniques for solving shifted sparse linear systems (see [36–38] and the references therein). Because the LRCF-ADI techniques are inefficient for the semistable systems, we will work for it as well.



---

*Source: 8071869-2021-10-19.xml*