The first-order index-1 descriptor system consisting of sparse sub-matrices is formed in a system of matrix equation having input-output relations

E_1 x ̇_1 (t)=A_1 x_1+A_2 x_2+B_1 u(t),

0=A_3 x_1+A_4 x_2+B_2 u(t)                                   (1)

y(t)=C_1 x_1+C_2 x_2+Du(t)

x(t_0 )=x_0,   t≥t_0.

In the system (1), E_1∈R^(n_1×n_1 ) is the differential coefficient matrix and A_1∈R^(n_1×n_1 ), A_2∈R^(n_1×n_2 ), A_3∈R^(n_2×n_1 ), A_4∈R^(n_2×n_2 )are the state sub- matrices, respectively, whereas B_1∈R^(n_1×p),B_2∈R^(n_2×p)are the control multiplier sub-matrices, C_1∈R^(m×n_1 ), C_2∈R^(m×n_2 )are the state multiplier sub-matrices, and D∈R^(m×p)is the direct gain for n_1+n_2=n with p,m≪n.  In the power system models, the matrix D is zero or absent.  Also, x_1 (t)∈R^(n_1 ),x_2 (t)∈R^(n_2 )are the state vectors, u(t)∈R^p is the control (input), and y(t)∈R^m is the output. The mentioned sub-matricesE_1 and A_1 are invertible (Hossain and Uddin, 2019).

Over the years, computational techniques are modified for the feasibility of massive systems. But compu-tational intricacy and memory extravagant simulation tolls keep the approaches impractical. Thus, large dimensional empirical systems are needed to be reduced to smaller dimensional systems. The strategy of finding a reduced-order model (ROM) of a given system is called model order reduction (MOR), currently earning vital attention to the researchers. ROM-based techniques are extensively adopted in the technological fields of the modern scientific world (Flagg et al., 2010). An adjustable ROM size and keeping the structural design invariant of the test system is the earnest sought (Rahman et al., 2020). The mechanism of the MOR techniques needs to essentially robust maintaining minimized norm of the error. 

To duly synchronize the iterative steps of the manipulation, the index-1 descriptor system (1) is to be restructured by implementing the well-known Schur-complements as

x∶=x_1, 

E≔E_1,

A≔A_1-A_2 A_4^(-1) A_3,

B≔B_1-A_2 A_4^(-1) B_2,                                           (2)

C≔C_1-C_2 A_4^(-1) A_3,

D≔D-C_2 A_4^(-1) B_2.

Then, the index-1 descriptor system (1) can be converted to an equivalent generalized LTI con-tinuous-time system as

Ex ̇(t)=Ax(t)+Bu(t),

y(t)=Cx(t)+Du(t).                                   (3)

The essence of LTI continuous-time systems is undeniable in the phenomena of the real-world engin-eering models, such as mechatronics, aeronautics, and system and control theory (Rahman et al., 2021). In the analysis of system stability and its appliances, the Continuous-time Algebraic Riccati Equation (CARE) is the pivot factor (Bänsch et al., 2015). The CARE explored from (3) can be written as

A^T XE+E^T XA-E^T X〖BB〗^T XE+C^T C=0.     (4)

A finite solution Xof the CARE (4) is obtainable if no eigenvalue of the corresponding Hamiltonian matrix is purely imaginary (Abou-Kandil et al., 2012). The symmetric positive-definite solution X is called stabilizing for the stable closed-loop matrix 〖A-BB〗^T XE.

Riccati-based feedback stabilization is obligate for an unstable type of system (3), and the opti-mal feedback matrixK^o=B^T XEis the premiering radiant in this process (Chen and Qui, 2015). To stabilize the target system optimally, the matrix A is to be replaced byA_s=A-BK^o. Then, the stabilized system can be formed as

Ex ̇(t)=A_s x(t)+Bu(t),

y(t)=Cx(t)+Du(t).                                      (5)

Computationally proficient solvers or analytical tools for large-scale CARE arising from the index-1 descriptors systems are still not technically cheap. Some Kleinman-Newton methods exist, which are very complicated, time laborious, and predefined structures are required (Mena and Saak, 2008). Alternating direction implicit (ADI) based approaches is computationally extravagant for the prerequisite large solution space, preconditioned initial parameters, and time-costly matrix factorization (Hasan and Uddin, 2018). 

The linearization ability and enforcement of initial priories boost up the convergence rate of the simulations by the projection-based approach RKSM, which is efficiently applicable in systems with perturbations (Simoncini et al., 2013). We are extending the techniques of the RKSM approach for the standard systems that is discussed in (Simoncini, 2016). The extension is including the modification of the existing work in such a way that it can apply for index-1 descriptor systems. We will investigate the effect of the shift parameters in the RKSM algorithm and compare this for the adaptive ADI shifts and heuristic shifts, respectively.

Also, to optimally stabilize the unstable index-1 descriptor system, we will apply the Riccati-based feedback matrix and analyze the stabilized system through eigen space and transient behaviors.

2. Rational Krylov subspace method for solving index-1 descriptor systems

If the eigenvalues of the matrix pair(A,E) sati-sfyλ_i+λ ̅_j≠0; for the all i,j=1,2,…,m that ensures the solutionXof the CARE (4) exists and unique. For a given set of competent shift para- meters〖 μ〗_i∈C^+; i=1,2,•••,m, the requisite m- dimensional rational Krylov subspace can be generated to pursue the orthogonal projectorV∈R^(m×m). The desired rational Krylov subspace can be constructed as

K_m=span(█(C^T,(A^T-μ_1 E^T )^(-1) C^T,@……,∏_(i=1)^m▒〖(A^T-μ_i E^T )^(-1) C^T 〗)).

Again, consider the CARE (4) and apply the Galerkin condition to it. After the simplification by matrix algebra, a low-rank CARE can be achieved as

A ̂^T X ̂E ̂+E ̂^TX ̂A ̂-E ̂^T X ̂〖B ̂B ̂〗^T X ̂E ̂+C ̂^T C ̂=0,       (6)

Here, X ̂=V^T XV, E ̂=V^T EV, A ̂=V^T AV, B ̂=V^T B, and C ̂=CV.

The equation (6) is a low-rank CARE and can be solved by MATLAB care command or any existing methods, such as the Schur-decomposition method. 

For the quick and smooth convergence of the proposed approach, adjustable shift selection is crucial (Druskin et al., 2010). There are many schemes to find suitable shift parameters, several approaches are given in (Moret and Popolizio, 2014; Benner et al., 2010) and references there in. In the present work, we are adopting the adaptive ADI shift and the heuristic shift parameters approach for index-1 descriptor systems.

In the iterative simulation process, to generate the feasible orthogonal projector V the existing set of shift parameters must be extended within the solution space. By the solutionX ̂ of the low-rank CARE (6), the solutionXof the CARE (4) will be approximated as X≈VX ̂V^T.

To have a stopping condition, a suitable approximation of the residual is required. The residual of m-th iteration can be executed in terms of the Frobenius norm ‖.‖_(F ) as

‖R‖_F=‖SJS^T ‖_F  ;   J=[■(0&I&0@I&0&I@0&I&0)].                 (7)    

Here the block upper triangular matrixSis found from the QR-factorization of the matrix U is defined as

U=[v_(m+1) μ_(m+1) E^T V_m X ̂H_m^(-T) e_m h_(m+1,m)^T-(I_m-V_m V_m^T)A^T v_(m+1)].                                                 (8)

The factor β_0=R_0 needs to be computed from the QR-factorizationC^T=Q_0 R_0. Then, the relative residual can be estimated as

‖R‖_F^((Relative))=‖R‖_F/‖〖β_0^T β〗_0 ‖_F .

The discussion about the terms and notations of the equation (8) are given in (Uddin, 2020) and the references there in. 

3. Sparsity preservation techniques

In system (3), the matrix A  is dense and for this reason, the rate of convergence of their structured system will be disturbed and accuracy will be affected. To avoid those inconveniences, instead of the explicit conversion of A, a more convenient approach will be done. To findv_i, at each of the iterations, a simple linear system will be solved as

〖(A〗^T-μ_i E^T)v_i=V_(i-1),

[■(A_1^T-μ_i E_1^T&A_3^T@A_2^T&A_4^T )][■(v_i@Γ)]=[■(V_(i-1)@0)].                               (9)

HereΓis the truncated term. Due to sparse system structure, the higher dimensional linear system (9) is conveniently solvable by the conventional direct solvers (Uddin, 2020). For the improvement of the consistency of the RKSM approach, explicit forms of the reduced-order matrices will not be used to construct the reduced-order system. The reduced-order matrices in sparse form can be attained in the following way

E ̂≔V^T E_1 V,

A ̂≔〖V^T A〗_1 V-(〖V^T A〗_2)A_4^(-1) 〖(A〗_3 V),

B ̂≔〖V^T B〗_1-(〖V^T A〗_2)A_4^(-1) B_2,                            (10)

C ̂≔C_1 V-C_2 A_4^(-1) 〖(A〗_3 V).

The low-rank solution X ̂ is symmetric and posi-tive-definite and can be factorized asX ̂=YY^T. Consider the desired low-rank factored solution of the CARE (4) asZ=VY. Thus, the solutionXof the (4) can be attained as X≈ZZ^T such that X≈VX ̂V^T=(VY) (VY)^T=ZZ^T.

For the future steps, the factored solutionZ=VYwill be stored and hence the optimal feedback matrixK^o=B^T XE=B^T (ZZ^T)Ecan be estimated. This process is iterative and will continue until the desired convergence is achieved. The whole process is summarized in Algorithm-1.

Finally, applying A_s=A-BK^o, the optimally stabilized LTI continuous-time system can be written as (5). To preserve the structure of the system, it needs to back to the original structure (1), and for this, the target system needs to be re-written as

E_1 x ̇_1 (t)=〖(A〗_1-B_1 〖K^o)x〗_1+A_2 x_2+B_1 u(t),   

0=〖(A〗_3-B_2 K^o)x_1+A_4 x_2+B_2 u(t),             (11)

y(t)=C_1 x_1+C_2 x_2+Du(t). 

4. Numerical results

The eligibility and efficiency of the proposed tech-niques for desired shift parameters are investigated by subjecting a Brazilian Interconnected Power System (BIPS) model, namely BIPS-606 with 606 differential and 6529 algebraic variables, which is of the type unstable index-1 descriptor system (Rommes et al., 2019; Freitas et al., 2008). The proposed techni-ques are applied with a tolerance level 〖10〗^(-10)to find low-rank approximate solutions for both adaptive ADI and heuristic shifts.

To stabilize the target models the optimal feedback matrix K^o isused. The comparisons of several compu-tational aspects are given in Table 1.

Table 1: Comparative analysis

From the above table, it is observed that by using adaptive ADI shifts a smaller solution space can be achieved, which is feasible for the memory allocation of the computational tool. In contrast, heuristic shifts are better for minimizing computational time. 

The comparison of the relative residuals and eigen-value stabilization for the RKSM approach exerting the adaptive ADI and heuristic shifts are shown in Fig 1, and Fig 2.

From the above figures, it is evident that that adaptive ADI shifts provide better relative residuals than heuristic shifts. But the stabilized eigenvalues for both of the shifts are very identical. 

The sub-figures of Fig 3 depict the stabilization of the step-responses of dominant input-output relations of the target model.

From the above sub-figures, it can be said that by utilizing both of the adaptive and heuristic shifts the unstable index-1 descriptor systems can be efficiently stabilized.

In this work, the techniques to optimally stabilize the unstable index-1 descriptor system by the optimal feedback matrix attained by the solution of the Riccati equation, are introduced and embedded. To serve the purpose, the CARE corresponding to a model derived from the Brazilian Interconnected Power System has been efficiently solved by the iterative RKSM app-roach utilizing adaptive ADI shift and heuristic shift parameters. The robustness and time-dealing efficiency of the computation is justified by applying the techniques to the target model. From the tabular and graphical comparisons, it can be concluded that the RKSM approach utilized both the adaptive ADI and heuristic shifts, stabilized the target model with the desired efficiency. But adaptive ADI shifts are feasible for the number of iterations and space dimensions, whereas heuristic shifts are suitable for the numerical rank and computation time. 

Thus, adaptive ADI shifts are applicable for the feasibility of memory allocation and heuristic shifts are comfortable for quick convergence.

We would like to express our profound gratitude to Md. Motlubar Rahman, Associate Professor in Mathematics, Directorate of Secondary and Higher Edu-cation, Ministry of Education, Bangladesh for his in-valuable suggestions and support in writing this article successfully. 

The authors declared that there is no conflict of interest in this article.