Computationally efficient optimal control for unstable power system models

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.


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]. Multi-tasking systems having various components arise in many fields of engineering applications, such as microelectronics, micro-electro-mechanical systems, cybersecurity, computer control of industrial processes, communication systems, etc. These systems are composed of branches of sub-systems and are functioned by very 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 multi-connected 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) aṡ where x 1 ∈ X 1 ⊂ R n1 is the vector with differential variables, x 2 ∈ X 2 ⊂ R n2 is the vector with algebraic variables and P ∈ R n3 is the vector of parameters with n 1 + n 2 = n [4,5]. In the state-space representation, the dynamic state variables x 1 and instantaneous variables x 2 are defined for the specific system, where the parameter P defines the configuration and the operation condition. The state variables in x 1 are time-dependent generator voltages and the parameter P is composed of the system parameters. The control devices together form f (x 1 , x 2 , P ) and the power flow balance form g(x 1 , x 2 , P ). In case of voltage stability, some equations of the f (x 1 , x 2 , P ) will not be considered. For a fixed parameter P 0 , linearizing the 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 where E, A ∈ R n×n , B ∈ R n×p , C ∈ R m×n and D ∈ R 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 the system (2), x(t) ∈ R n is the state vector and u(t) ∈ R p is control (input), while y(t) ∈ R m is the output vector and considering x(t 0 ) = x 0 as the initial state. In most of the state-space representations the direct transmission remains absent and because of that D = 0. Since E is singular (i.e. det(E) = 0), the system (2) is called the descriptor system [7,8].
Here x 1 ∈ R n1 , x 2 ∈ R n2 with n 1 + n 2 = n are state vectors and other submatrices are sparse in appropriate dimensions. If E 1 and J 1 have full rank, and J 4 is non-singular (i.e. det(A) = 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 the descriptor system (2) can be converted to the generalized LTI continuous-time system where, we have considered the following relations Lemma 1.1 (Equivalence of transfer functions [10]). Assume the transfer functions G(s) = C(sE −A) −1 B +D and G(s) = C(sE −A) −1 B +D are obtained from the index-1 descriptor system (2) and the converted generalized system (3), respectively. Then, the transfer functions G(s) and G(s) are identical and hence those systems are equivalent.
Though the systems (2) and (3) are equivalent, the system (2) is sparse and the 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 the system (3) is defined as If the Hamiltonian matrix corresponding to the 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 − (BB T )XE. Riccati-based feedback matrix has a prime role in the stabilization approaches for unstable systems [15,16]. To find an optimal feedback matrix K o , the Linear Quadratic Regulator (LQR) problem technique can be applied, where the cost functional is defined as The cost functional (6) can be optimized as J(u o , x 0 ) = x 0 T Xx 0 by applying the optimal control u o = −K o x(t) generated by the optimal feedback matrix K o = B T XE associated with the solution matrix X of the CARE (5). Using the optimal feedback matrix K o , an unstable LTI continuous-time system can be optimally stabilized by replacing A by A s = A −BK o . The stabilized system can be written as The eligibility of Rational Krylov Subspace Method (RKSM) for the largescale LTI continuous-time systems has discussed by Simoncini in [17] and the application of adaptive RKSM to solve large-scale CARE for finding optimal control of the LTI systems has narrated by Druskin et al. in [18]. Analysis of the basic properties of RKSM for solving large-scale CAREs subject to LTI systems investigated by Simoncini in [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. in [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 Kuerschner in [21].
Since 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.

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.

Generalized RKSM technique for solving Riccati equations
In [19] Simoncini applied RKSM approach for solving the CARE in the form associated with the LTI continuous-time system If the eigenvalues of the matrix pair (A, E) satisfy λ i +λ j = 0, ∀i, j = 1, 2, · · · , m, that ensures the solution X of the CARE (8) exists and unique.
The orthogonal projector V ∈ R n×m spanned by the m-dimensional rational Krylov subspace for a set of shift parameters µ i ∈ C + ; i = 1, 2, · · · , m is defined as If θj are the eigenvalues of (V T AV, V T EV ); Sm ∈ C + approximates the mirror eigenspace of A − B(B T XE) and δSm is its border, the shifts are computed from .
According to the Galerkin condition and after simplification by matrix algebra, a low-rank CARE can be obtained aŝ (10) is an approximated low-rank CARE and can be solved by any conventional method or MATLAB care command. HereX ∈ R m×m is taken as low-rank approximation of X, corresponding to the low-rank CARE (10). Then residual of the m-th iteration is where . F denotes the Frobenius norm and S is a block upper triangular matrix in the QR factorization of the matrix U derived as 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 C T = Q0R0 such that R0 = β0, the relative-residual can be estimated as Through RKSM the low-rank factor Z of the approximate solution X of the CARE (8) needs to be estimated, such that X ≈ ZZ T . The low-rank solutionX is symmetric positive definite and the original solution X can be approximated as X = VXV T . By the eigenvalue decomposition to the approximate solutionX and truncating the negligible eigenvalues, the possible lowest order factor Z of X can be estimated as Here, Λ2 consists the negligible eigenvalues. Summary of the above process is given in the Algorithm-1.
Compute shift for the next iteration.

Modified LRCF-ADI method for solving Lyapunov equations
The generalized Continuous-time Algebraic Lyapunov Equation (CALE) associated with the LTI continuous-time system (9) is The shift parameters µi ∈ C − are allowed and the initial iteration is taken as X0 = X0 T ∈ R n×n . Assume Zi ∈ R n×ip as the low-rank Cholesky factor of Xi ∈ R n×n such that Xi = ZiZ T i [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 = −2 Re(µi), the conventional low-rank Cholesky factor ADI iterations yield the form as Assume a set of adjustable shift parameters, for two subsequent block iterates Vi, Vi+1 of the ADI technique related to the pair of complex conjugated shifts {µi, µi+1 := µi} it holds 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 Algorithm 2: Modified LRCF-ADI.
Input : E, A, C, i max (number of iterations) and shift parameters Then, for a pair of complex conjugate shifts at any iteration, the low-rank factor Zi can be computed as The modified techniques discussed above is summarized in Algorithm-(2).

Basic Kleinman-Newton method
For the system (9), the residual R(X) of the CARE (8) is Form the Fre'chet derivative at X = Xi, we have Consider the Newton iteration and apply (17), we have Now, put ∆Xi = Xi+1 − Xi in (18). Then, after simplification we get Then, assumeÃi = A − BB T XiE and Wi = C T E T XiB , then equation (19) reduces to a generalized Lyapunov equation such as Algorithm 3: Basic Kleinman-Newton.
Input : E, A, B, C and X 0 (initial assumption).
Output: Approximate solution X and feedback matrix K.
The generalized CALE (20) 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 = B T Xi+1E can be estimated. The whole mechanism is called the Kleinman-Newton method [26] for solving generalized CARE. The summary of the method is given in Algorithm-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, 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.

Updated Rational Krylov subspace method
Let us consider the LTI continuous-time system (3) and introduce an orthogonal projector V spanned by the m dimensional rational Krylov subspace for a set of given shift parameters µi ∈ C + ; i = 1, 2, ..., m is defined as Again, consider the CARE (5) and apply the Galerkin condition on it. Then, after the simplification by matrix algebra, a low-rank CARE can be achieved aŝ where, (21) 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.

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 matrix Tm can be performed more efficiently than the explicit product Tm = V T m AVm. To find the stopping criteria, we have to consider following lemmas. Lemma 3.1 (Arnoldi relation [28]). Let Km = span(Vm) be the rational Krylov subspace with the shift parameters µi ∈ C + ; i = 1, 2, · · · , m. Then Km satisfies the Arnoldi relation as follows is the QR decomposition of the right hand side matrix with g T m = β1hm+1,me T m H −1 m . Lemma 3.2 (Building the projected matrix [29]). Let the column vectors of Vm be an orthonormal basis of the rational Krylov subspace Km with the block diagonal matrix Dm = diag(µ1E T , µ2E T , · · · , µmE T ), where {µ1, µ2, · · · , µm} is the set of shift parameters used in the algorithm. Then for the projected matrix Tm = V T m AVm the following relation holds Theorem 3.1 (Residual of the RKSM iterations). Let Vm be the orthogonal projector spanned by the rational Krylov subspace Km and X ≈ VXV T is the solution of the CARE (5) using the low-rank solutionX. Then, the residual of m-th iteration can be computed as where . F denotes the Frobenius norm and S is a block upper triangular matrix in the QR factorization of the matrix U is defined as Proof. Assume f = (In − VmV T m )A T vm+1 and consider the reduced QR factorization C T = V0β0. Then, by putting the relations in equation (22) of Lemma-(3.1), the relation can be written as The residual of the CARE (5) can be written as Consider the approximate solution using the low-rank solutionX as X = VmXV T m and equation- (23) in Lemma-(3.2), then applying (26) in (27), we get Thus, the proof follows from the Frobenius norm of the equation (28).
Then the relative-residual with respect to β0 can be estimated as follows

Sparsity preservation
The matrix A 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 Here Γ = −J −1 4 J3vi is the truncated term. The linear system (30) 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 following waŷ

Treatment for the unstable systems
If the system is unstable, a Bernoulli stabilization is required through an initialfeedback matrix K0 to estimate A f = A − BK0 and the matrix A needs to be replaced [31]. Then, the system (3) and CARE (5) need to be re-defined as Then, for every iterations K =B TX V T E needs to be updated by the solutionX of (10) and the rational Krylov subspace for the projector V needs to be redefined as For the stabilized system using the initial-feedback matrix K0, the expressions (30) can be written as To evaluate the shifted linear systems, explicit inversion of A − BK should be avoided in practice, instead the Sherman-Morrison-Woodbury formula needs to be used as follows
Solve the linear system (34) for v m+1 . 6 Compute adaptive shifts for the next iterations (if store is empty). 7 Using Arnoldi algorithm orthogonalize v m+1 against V m to obtain v m+1 , such that V m+1 = V m ,v m+1 .   Compute R m (relative) for convergence. 11 end while 12 Compute eigenvalue decomposition 13 For negligible eigenvalues truncate Λ 2 and construct Z = V m+1 T 1 Λ 1 2 1 . 14 Compute the optimal feedback matrix

Estimation of the optimal feedback matrix
The low-rank solutionX is symmetric and positive definite and can be factorized aŝ X = Y Y T . The original solution can be approximated as Finally, the desired low-rank factored solution Z = V Y of the CARE (5) will be stored and the optimal feedback matrix K o = B T XE = B T (ZZ T )E can be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in the Algorithm-4.

Updated LRCF-ADI Based modified Kleinman-Newton method
In each iteration of Algorithm-3, 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 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.

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 3.3 (Convergence of ADI iterates [32]). Let Xi be the solution of the CALE (11) and consider X (k) i be an iterate of the ADI method. Then for all µi ∈ C − the relation holds where Lemma 3.4 (Residual of the ADI iterates [32]). Let X (k) i be an iterate of the CALE (11) by the ADI method. Then considering A k, Theorem 3.2 (Residual factor of the ADI iterations). The residual for CALE (11) at i-th iteration of ADI method is of rank at most m and is given by Proof. Consider X is the solution of the CALE (11) 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 For all µi ∈ C − , the Stein's equation is equivalent to the CALE (11) and it can be written as Using the Stein's equation (39) in Lemma-3.3, we can find some initial iterates X
Thus, for WiW T i ≤ τ 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 Then, using (41) the residual factor Wi can be derived in a recursive form as In case of real setting, µi+1 :=μi needs to be considered to find the following form The summary of above techniques is given in the Algorithm-(5).

Adjustment for the unstable systems
For the unstable index-1 descriptor system, initial feedback matrix K0 needs to be introduced and instead of (Ã (i) , E), corresponding shift parameters are needed to be computed from eigen-pair (Ã (i) − BK0, E). The sparse form of the eigen-pair can be structured as To find V (i) j , in each ADI (inner) iteration a shifted linear system needs to be solved as Thus, V (i) j can be obtained from the sparse form of the shifted linear system structured as

Estimation of the optimal feedback matrix
The feedback matrix T )E needs to be computed in each ADI (inner) iteration and the optimal feedback matrix K o = K (imax) needs to be stored after the final Newton (outer) iteration. The summary of the modified method is given in the Algorithm-(6).

Computation of the optimally stabilized system and the optimal control
The optimal feedback matrix K o = B T XE can be achieved by the feasible solution X of the Riccati equation (5). Then, applying As = A − BK o , optimally stabilized LTI continuous-time system can be written as (7). To preserve the structure of the system, it needs to back to the oginial structure (2), and for this the submatrices J1, J3 are replaced by J1 − B1K o , J3 − B2K o , respectively. Finally, the desired optimal control u o = −K o x(t) of the targeted power system models can be computed.

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.

Brazilian inter-connected 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 Inter-connected Power Systems (BIPS) is one of the most convenient examples of the power system models with various test systems [34]. The following Table 1 provides the details about the models. The detailed structure of those will be found at https://sites.google.com/site/rommes/software, 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 − 1963 Algorithm 6: Modified KN-LRCF-ADI (sparsity-preserving).
Compute adaptive shifts µ j . 22 end for have stable eigenvalues [35]. Here the names of the models are considered according to their number of states.   Table 2 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 both of the tables, Table 2 and Table 3 we have analyzed the same features of the stabilized BIPS models, we can easily compare the efficiency and robustness of the proposed methods.
From the above tables 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.   From the sub-figures in the above mentioned figures, it has been observed that the eigenvalues of the models mod − 606 and mod − 2476 are stabilized very well but its marginal for the model mod − 1998. Thus, it can be concluded that both the RKSM and LRCF-ADI based Kleinamn-Newton techniques have adequate efficiency to stabilize the unstable descriptor systems by closed-loop structure via Riccati based feedback stabilization. Fig.4 displays the applicability of the RKSM for the semi-stable descriptor system, whereas LRCF-ADI based methods are ineffective in this case. Since the eigenvalues of the semi-stable 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.

Stabilization of step-responses
The investigation of the figures from Fig.5-Fig.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. Since the power system models are of 4 × 4 input-output relations, there are 16 step-responses for each models.
For the effective comparison we have investigated only some graphically significant step-responses.  Fig.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 indicates by RKSM is suitably robust. On the other hand, though sometimes the Kleinman-Newton approach provides very good accuracy but it has some scattered behaviors.
Moreover, Fig.11 and Fig.12 shows the applicability of the RKSM technique for the Riccati based feedback stabilization for the semi-stable index-1 descriptor systems.

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 semi-stable 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, while 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 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 and apply the "Restarted Hessenberg Method" for solving shifted sparse linear systems. Since the LRCF-ADI techniques are inefficient for the semi-stable systems, we will work for it as well.