Sparsity-Preserving Two-Sided Iterative Algorithm for Riccati-Based Boundary Feedback Stabilization of the Incompressible Navier–Stokes Flow

. In this paper, we explore the Riccati-based boundary feedback stabilization of the incompressible Navier–Stokes fow via the Krylov subspace techniques. Since the volume of data derived from the original models is gigantic, the feedback stabilization process through the Riccati equation is always infeasible . We apply a H 2 optimal model-order reduction scheme for reduced-order modeling, preserving the sparsity of the system. An extended form of the Krylov subspace-based two-sided iterative algorithm (TSIA) is implemented, where the computation of an equivalent Sylvester equation is included for minimizing the computation time and enhancing the stability of the reduced-order models with satisfying the Wilson conditions. Inverse projection approaches are applied to get the optimal feedback matrix from the reduced-order models. To validate the efciency of the proposed techniques, transient behaviors of the target systems are observed incorporating the tabular and fgurative comparisons with MATLAB simulations. Finally, to reveal the advancement of the proposed techniques, we compare our work with some existing works.


Introduction
In the present day, the incompressible Navier-Stokes fow is one of the most impressive contents of interest of the researchers. It has a prominent impact on fuid mechanics, naval engineering, oceanography, and water resource scientists. In the analysis of fuid properties and attributes of undersea atmospheres, the incompressible Navier-Stokes fow is one of the prominent issues for naval architects. It is an essential part of the investigation of sea-route analysts.
Details of the Navier-Stokes equations with the discretization can be found in references [1,2]. Linearizing the Navier-Stokes equations in space and time variables by a mixed fnite element method without altering the time variable converts them into linear time-invariant systems [3]. Te incompressible Navier-Stokes fow can be written as the following diferential-algebraic equations: gradient. B 1 ∈ R n v ×n u and B 2 ∈ R n p ×n u are the input matrices, whereas C 1 ∈ R n y ×n v and C 2 ∈ R n y ×n p are the output matrices. Te input-output matrices quantify the velocity pattern at the corresponding nodes. D ∈ R n y ×n u is the direct transform matrix [4]. Considering M is invertible, A 2 and A 3 have both full column ranks and A 3 M − 1 A 2 is nonsingular. Hence, system (1) is of the index-2 diferential-algebraic system having the dimension n v + n p [5]. A system equivalent to system (1) can be written in the following form: where , and D � D with appropriate dimensions. Te matrix pencil corresponding to the system is defned as (A, E), where the matrix E is singular and the matrix pencil (A, E) has n v − n p number of nonzero fnite complex eigenvalues and 2n p infnite eigenvalues [6]. Reynolds number (Re) is one of the key features that instigate the system characteristics. For Re ≥ 300, the Navier-Stokes fow turns unstable and analysis of the fow attributes becomes troublesome [7]. Optimal feedback stabilization of the Navier-Stokes fow is indispensable and Riccati-based boundary feedback stabilization is the best apparatus in this situation [8]. Te Riccati equation corresponding to the system (2) is of the form: Te solution X of the Riccati equation (3) is symmetric positive defnite and called stabilizing for a stable closedloop matrix A − (BB T )XE. Te optimal feedback matrix K o can be computed as K o � B T XE and can be implemented for optimally stabilizing the target system by replacing A by A s � A − BK o [9]. Te stabilized system can be written as follows: With the increasing number of fow components and fner meshes in the discretization process, the dimensions of the matrices in system (1) become very large. Due to the size of the system matrices, Riccati-based boundary feedback stabilization of system (1) is infeasible as the associated Riccati equation cannot be solved by the usual matrix equation solvers [10]. To overcome this adversity, a suitable model-order reduction (MOR) technique needs to be applied [11]. A computationally convenient reduced-order model (ROM) of system (2) can be written as follows: where E � W T EV, A � W T AV, B � W T B, C � CV, and D � D. Te projector V and W can be found in any compatible MOR technique [12]. Using the reduced-order matrices given in formula (5), the reduced-order form of Riccati equation (3) can be obtained as follows: Now, the optimal feedback matrix for system (2) can be approximated by applying the ROM (5). Te reduced-order Riccati equation (6) is to be solved for a symmetric positivedefnite matrix X by the MATLAB library command care. Ten, the stabilizing feedback matrix K � B T XE corresponding to the ROM (5) can be estimated, and hence, the approximated optimal feedback matrix K o � KV T E [13]. Finally, utilizing K o stabilization of system (2) can be carried out as system (4).
Tere are a number of MOR techniques that currently exist and are in use, namely, singular value decomposition (SVD)-based balanced truncation (BT) and Krylov subspace-based iterative rational Krylov algorithm (IRKA). Tose techniques are elaborately discussed in references [8,[14][15][16][17][18][19][20][21] and references therein. From the discussion of the BT method, some obstacles are identifed, such as its demands to solve two large-dimensional Lyapunov equations for controllability and observability Gramian, which is very costly for computational time and memory requirements. On the other hand, IRKA is better in the sense of simulation time and needs less memory allocation but the stability of the ROM is not guaranteed. For solving the Riccati equation without explicit estimation of the ROM, some techniques are available in practice, for example, low-rank alternative direction implicit (LR-ADI), integrated Newton-Kleinman (NK) method, and Krylov subspace associated rational Krylov subspace method (RKSM). Tose methods are derived and analyzed in detail in references [22][23][24][25] and references therein. Te Newton-Kleinman method is a very complex approach and at each Newton step, LR-ADI iterations need to be executed once, which is a very time-laborious task [26][27][28][29]. Also, the Newton-Kleinman method cannot ensure the defniteness of the solution of the Riccati equation derived from an unstable system. On the contrary, RKSM is straightforward for simulation but lack of proper shift parameter selection sometimes makes this method inefective and adjustment of the stopping criteria may encounter the convergence of the method. An initial feedback matrix can be a remedy for unstable systems in the RKSM approach but this additional step makes the total process inconvenient in the sense of time management [30].
To overcome abovementioned troubles and complexities, we are proposing an extended form of Krylov subspacebased two sided iterative algorithm (TSIA) as discussed in references [31][32][33] and references therein. In this strategy, initially, we need to fnd a ROM of the target model implementing any conventional technique, then two generalized sparse-dense Sylvester equations will be solved to fnd the system Gramian; hence, the required projection matrices constituted by their orthonormal columns constructed via QZ-decomposition [34]. Te ROM attained by the TSIA approach is stability-preserving and we will derive its sparsity-preserving form of it. Ten, rest of the activities will follow the classical inverse projection scheme discussed previously. Te validity of the proposed techniques will be investigated numerically through the transient behavior of the target models and accuracy will be justifed by the H 2 norm optimality.

Background
Here, we discuss the basic approaches for the treatment of a generalized state-space system. Ten, we briefy analyze the conversion of index-2 descriptor systems into equivalent generalized systems with several cases. At the end of this section, we investigate sparsity-preserving reduced-order matrices and the corresponding reducedorder model. Te generalized state-space system can be written as follows: where E ∈ R n×n is nonsingular, and A ∈ R n×n , B ∈ R n×p , C ∈ R m×n , and D ∈ R m×p . Te transfer function of system (7) is defned by G Te CARE defned for system (7) is as follows: According to the Fre'chet derivative at X � X i and applying the Newton iteration with the condition ΔX i � X i+1 − X i , we obtain the following: then formula (9) reduces to a generalized Lyapunov equation such as follows: Te generalized CALE (10) can be solved for X i+1 by any conventional method, and the corresponding feedback matrix K i+1 � B T X i+1 E can be estimated [35]. Te summary of the method is given in Algorithm 1.
Te solution of the generalized CALE (10) can be found through the Cholesky factorization utilizing the low-rank Cholesky-factor alternative direction implicit (LRCF-ADI) techniques. Details of the LRCF-ADI techniques are provided in references [36][37][38] and references therein. In this approach, only the Cholesky factor Z of the solution X in needs to store. Te summary of LRCF-ADI techniques is given in Algorithm 2.

Rational Krylov Subspace Algorithm.
Te rational Krylov subspace algorithm is applicable to fnd the factored solution of CARE (8) without any conversion. Tis algorithm is defned for the LTI systems of symmetric cases. An orthogonal projector V ∈ R n×m spanned by the m-dimensional rational Krylov subspace needs to construct for fnding the an approximated low-rank CARE defned as follows: Te low-rank CARE (11) can be solved for X∈ R m×m by any conventional method or MATLAB care command, where X is taken as a low-rank approximation of X. Ten, the factored solution Z of the solution X is estimated to fnd the optimal feedback matrix K o . Te derivation and update of the RKSM techniques are given in references [39,40] and references therein. Te stepby-step RKSM approach is given in [41] and summarized in Algorithm 3.

Iterative Rational Krylov Algorithm. Iterative rational Krylov algorithm (IRKA) can be applied to construct an
such that its transfer function G(s) � C(sE − A) − 1 B + D interpolates the original one, G(s), at selected points in the complex plane along with selected directions. Te detailed procedure of IRKA is illustrated in reference [42]. Via the IRKA approach, two projection matrices V and W are to be constructed; hence, by the Petrov-Galerkin condition, we can construct the reduced-order matrices in formula (12) as follows: Using the reduced-order matrices defned in formula (13), the reduced-order CARE can be attained as formula (11) and can be solved for X as previously. Ten, the stabilizing feedback matrix K � B T XE corresponding to the ROM (12) can be estimated; hence, for stabilizing the full model (7), the approximated optimal feedback matrix K o can be retrieved by the scheme of reverse projection as follows: Te iterative rational Krylov algorithm (IRKA) introduced in reference [43] resolves the problem by iteratively correcting the interpolation points and the directions as summarized in Algorithm 4.

Conversion of the Index-2 Descriptor
System to Generalized System. Let us Assume M is a nonsymmetric positive defnite and A T 2 ≠ A 3 . Initially, we consider B 2 � 0. According to Section 3 of reference [44], the algebraic part of the system (1), p(t) can be expressed as v(t). Ten, by proper elimination and substitution, system (1) can be converted to an equivalent form: where the converted matrices are structured as follows: with Π T l v(0) � Π T l v 0 . Te projectors Π l and Π r are composed as follows: Input: E, A, B, C and X 0 (initial assumption). Output: approximate feedback matrix K ALGORITHM 2: LRCF-ADI for generalized systems Input: E, A, B, C, i max (number of iterations) and μ i (initial shifts).  Mathematical Problems in Engineering For the general case, we assume B 2 ≠ 0. According to Section of reference [45], we consider the velocity v(t) can be decomposed as follows: where Tus, the required converted system can be formed as follows: ) and the converted matrices are as follows: In addition, if both B 2 � 0 and C 2 � 0, the converted system will be the same as a system (15) but converted matrices will have the following form: According to the formulation derived in Section 4.1 of reference [46], sparsity-preserving reduced-order matrices for the index-2 descriptor systems can be written as follows: Here, the projection matrices V and W are composed of suitable sparsity-preserving matrix-vector formulation. Te ROM corresponding to the reduced-order matrices defned in (22) can be written as follows:

Two-Sided Iterative Algorithm for Index-2 Descriptor Systems
Tis section includes the main tasks of this work. We derive an improved version of the two-sided iterative algorithm (TSIA) for index-2 descriptor systems to stabilize the incompressible Navier-Stokes fow.

Formulation of the Generalized Sparse-Dense Sylvester Equation.
Initially, we assume a makeshift ROM with the desired dimension in the form (23) employing any classical iterative method accomplishing a few iterations. However, due to the use of a minimum number of iterations, the attained ROM cannot efciently approximate the original system. So, the quality of the ROM needs to be improved by further techniques such that it satisfes the Wilson conditions for H 2 optimality [42,47]. Now, we are to construct two generalized sparse-dense Sylvester equations. To do this, we need to compute the following matrices: (3) while (not converged) do (4) Compute the reduced-order matrices as (13).
Assuming P n � MB n B T and Q n � − MC n C, the required sparse-dense Sylvester equations can be formed as follows: AXE + EXA + P n � 0,

Solving the Generalized Sparse-Dense Sylvester Equation.
Sylvester equations formed in formula (25) can be solved by the techniques given in Section 3.2 of reference [48]. Since the simulations of the Sylvester equations in formula (25) are analogous, we will derive the method for computing X only.
To generate the Krylov subspace QZ-decomposition of E and A is to be computed, such that E � QSZ T and A � QTZ T , respectively. Here, S and T are the upper triangular matrices. As the techniques provided in reference [48] of the dense form, we have to rewrite the basis vector in the sparse form as follows: and the successive columns for the solution matrix X can be generated as follows: where P n (j) � P (j) Te updated sparse form of the desired techniques are provided in Algorithm 5.

Two-Sided Iterative Algorithm to Estimate the Optimal
Feedback Matrix. At frst, an iterative method with minimum iterations needs to be utilized to fnd an initial makeshift ROM with the desired dimension. With the required matrix-vector algebraic operations, the Sylvester equations defned in formula (25) need to be solved for X and Y, respectively. Ten, the approximated projector matrices V and W can be computed by the QR-decomposition of the matrices X and Y, respectively. Tose crude projector matrices need to be refned as V � V and W � W(V T M T W) − 1 , respectively. Implementing V and W the ROM (23) can be acquired with the reduced-order matrices defned in formula (22); hence, the reduced-order Riccati equation can be formed as follows: which can be solved for X by MATLAB care command, and then, the reduced-order feedback matrix can be estimated as K � B T XE. Finally, the optimal feedback matrix for system (1) can be approximated as follows: Te whole process is summed up in Algorithm 6.

Stabilization of the Index-2 Descriptor System.
Plugging the optimal feedback matrix K o , the unstable incompressible Navier-Stokes fow given by (1) can be optimally stabilized as follows:

H 2 Norm of the Error System. Let us consider G(s) and
G(s) are the transfer functions of the converted generalized system (15) and (23), respectively. Ten, the associated error system can be formulated as follows: where we have taken the reduced-order matrices of formula (22) into account and the system matrices are formed as follows: Authors in reference [49] derived the way of estimating H 2 norm of the error system (31) as follows: Here, H 2 norm of the full model is defned as ‖G(s)‖ H 2 and needs to be evaluated for once, which is inconvenient for the conventional simulation solvers. Assuming Z q is the low-rank Gramian factor of the Gramian Q that needs to be determined. In reference [50], a practically feasible technique is derived to overcome this situation. Ten, ‖G(s)‖ H 2 can be written as follows: Again, ‖G(s)‖ H 2 is the H 2 norm of the ROM that can be evaluated by the Gramian Q of the low-rank Lyapunov equation due to the involvement of the reduced-order matrices, the Lyapunov (35) is solvable by the MATLAB library command lyap.
Finally, trace(B T Q s B) can be computed by the Gramian Q s of the sparse-dense Sylvester equation that can be conveniently solved by reforming the techniques provided in Algorithm 5.

Numerical Results
In this section, the justifcation of the adaptability and efciency of the proposed techniques will be discussed. Computations are performed through MATLAB simulations to observe the betterment of the proposed techniques in comparison to the present methods. Te investigation involves both the graphical and tabular approaches. Computation time and accuracy of the approximation, and robustness of the stabilization of the transient behaviors will be the prime concern of this discussion. We have implemented the proposed techniques to the realworld Navier-Stokes models of unstable type. All the results have been achieved using 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.

Model Description.
In this work, we have considered some unstable Navier-Stokes models with the Reynolds numberRe � 300, 400, 500, respectively. Te dimension of the models varies from 1 to 5. Each of the models is nonsymmetric having 2 inputs with 7 outputs. For the target models, in the conventional frst-order index-2 descriptor system the submatrices B 2 and C 2 both are the sparse zero matrices. Te full specifcations of the target models are given in Table 1.

Approximation of the Full Models with the Reduced-Order
Models. Te accuracy of the approximation of the full models with the ROMs will be validated here. We will verify the level of approximation graphically and then evaluate the H 2 -norm of the error systems for the objective models.

Comparison of the Transfer Functions.
A graphical comparison between the full models and ROMs will be depicted. Since the Navier-Stokes models of various dimensions have the same system structure and transitional properties in the graphical analysis, we will only show the properties of the 3-dimensional model with Re � 500.
In Figure 1, the subfgure Figure 1(a) shows the comparison of the transfer function (sigma plot) of the full model and corresponding ROM, whereas subfgures Figures 1(b) and 1(c) evince the absolute and relative errors, respectively, of this reduced-order approximation. From the abovementioned fgures, it can evident that by the proposed techniques the objective full models can be properly approximated by the ROMs with a reasonable level of accuracy.

H 2 Norm of the Error System for the ROMs.
Here, H 2 norm of the error system for the ROMs of the target models will be provided.
(2) Construct B n according to (24) and compute P n � MB n B T .
ALGORITHM 6: TSIA for the optimal feedback matrix of index-2 descriptor systems.
Mathematical Problems in Engineering Table 2 conveys the desired H 2 error norm of the ROMs. From the measured H 2 error norm of the ROMs, it is evident that the level of the approximation process is suitably robust and confrms the compatible scale of the accuracy.

Graphical Validation of Stabilization of the Unstable
Systems. In this section, the stabilization of the unstable Navier-Stokes models will be graphically illustrated. For the compactness of this article, we will only demonstrate the stabilization of the transient behaviors of the 3-dimensional models. Figure 2 exhibit the eigenvalues of the original (unstable) 3-dimensional systems with Reynolds number Re � 300, 400, 500, respectively. In contrast, the subfgures of

Comparison of the TSIA with IRKA.
In this section, we are going to compare the proposed algorithm with the IRKA approach for reduced-order modelling. Initially, we aimed to compare with the RKSM approach as well. However, due to the nonsymmetric structure of the Navier-Stokes models, the RKSM approach is applicable for here. Since the stabilization of the eigenvalues and step-responses by the TSIA and IRKA via the reduced-order modeling is very identical, we have ignored their graphical comparisons and chosen the tabular comparisons. In this comparative analysis, we will investigate both of the computation time and the H 2 error norm of the ROMs of selected 3-dimensional models. Table 3 reveals the required information on computation time and H 2 error norm.   From Table 3, it is apparent that in the case of computation time, the TSIA algorithm requires about two-thirds of that of IRKA. Also, a comparison of H 2 error norm of the ROMs manifests the betterment of TSIA over IRKA.

Conclusions
We have introduced a two-sided projection-based sparsitypreserving reduced-order modeling approach for the stabilization of nonsymmetric index-2 descriptor systems explored from unstable Navier-Stokes models. Te Wilson conditions satisfying Reduced-order models are derived by the two-sided projection techniques implementing the sparse-dense Sylvester equations. To solve the desired Sylvester equations, sparsity-preserving Krylov subspaces are structured via the system of linear equations with a compact form of matrixvector operations. From the numerical analyses, it is ascertained that the proposed techniques can be profciently applied for the stabilization of the target models through reducedorder modeling. We have efciently approximated the full models by the corresponding ROMs with minimized H 2 error norm and stabilized the eigenvalues and step-responses properly. A comparison of reduced-order modelling by TSIA and IRKA is performed. From the comparative discussion, it is manifested that the techniques involved in TSIA outplayed that of the IRKA in both the saving computation time and H 2 error norm. Also, we intended to include RKSM in the comparative analysis but found that RKSM is not applicable to the target models due to the nonsymmetric structure. Tus, the bottom line of the article is that the proposed techniques can be utilized to stabilize the unstable Navier-Stokes models with better accuracy and less computing time.