Abstract
In this paper, we explore the Riccatibased boundary feedback stabilization of the incompressible Navier–Stokes flow 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 optimal modelorder reduction scheme for reducedorder modeling, preserving the sparsity of the system. An extended form of the Krylov subspacebased twosided 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 reducedorder models with satisfying the Wilson conditions. Inverse projection approaches are applied to get the optimal feedback matrix from the reducedorder models. To validate the efficiency of the proposed techniques, transient behaviors of the target systems are observed incorporating the tabular and figurative comparisons with MATLAB simulations. Finally, to reveal the advancement of the proposed techniques, we compare our work with some existing works.
1. Introduction
In the present day, the incompressible Navier–Stokes flow is one of the most impressive contents of interest of the researchers. It has a prominent impact on fluid mechanics, naval engineering, oceanography, and water resource scientists. In the analysis of fluid properties and attributes of undersea atmospheres, the incompressible Navier–Stokes flow is one of the prominent issues for naval architects. It is an essential part of the investigation of searoute 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 finite element method without altering the time variable converts them into linear timeinvariant systems [3]. The incompressible Navier–Stokes flow can be written as the following differentialalgebraic equations:where is the nodal vector of the discretized velocity with , is the discretized pressure, and is the input with the output . The symmetric positive definite matrix represents the mass, whereas the matrices , , and are for system components including discretized gradient. and are the input matrices, whereas and are the output matrices. The inputoutput matrices quantify the velocity pattern at the corresponding nodes. is the direct transform matrix [4]. Considering is invertible, and have both full column ranks and is nonsingular. Hence, system (1) is of the index2 differentialalgebraic system having the dimension [5]. A system equivalent to system (1) can be written in the following form:where , , , , , and with appropriate dimensions. The matrix pencil corresponding to the system is defined as , where the matrix is singular and the matrix pencil has number of nonzero finite complex eigenvalues and infinite eigenvalues [6].
Reynolds number is one of the key features that instigate the system characteristics. For , the Navier–Stokes flow turns unstable and analysis of the flow attributes becomes troublesome [7]. Optimal feedback stabilization of the Navier–Stokes flow is indispensable and Riccatibased boundary feedback stabilization is the best apparatus in this situation [8]. The Riccati equation corresponding to the system (2) is of the form:
The solution of the Riccati equation (3) is symmetric positive definite and called stabilizing for a stable closedloop matrix . The optimal feedback matrix can be computed as and can be implemented for optimally stabilizing the target system by replacing by [9]. The stabilized system can be written as follows:
With the increasing number of flow components and finer meshes in the discretization process, the dimensions of the matrices in system (1) become very large. Due to the size of the system matrices, Riccatibased 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 modelorder reduction (MOR) technique needs to be applied [11]. A computationally convenient reducedorder model (ROM) of system (2) can be written as follows:where , , , , and . The projector and can be found in any compatible MOR technique [12]. Using the reducedorder matrices given in formula (5), the reducedorder 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). The reducedorder Riccati equation (6) is to be solved for a symmetric positivedefinite matrix by the MATLAB library command care. Then, the stabilizing feedback matrix corresponding to the ROM (5) can be estimated, and hence, the approximated optimal feedback matrix [13]. Finally, utilizing stabilization of system (2) can be carried out as system (4).
There are a number of MOR techniques that currently exist and are in use, namely, singular value decomposition (SVD)based balanced truncation (BT) and Krylov subspacebased iterative rational Krylov algorithm (IRKA). Those techniques are elaborately discussed in references [8, 14–21] and references therein. From the discussion of the BT method, some obstacles are identified, such as its demands to solve two largedimensional 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, lowrank alternative direction implicit (LRADI), integrated Newton–Kleinman (NK) method, and Krylov subspace associated rational Krylov subspace method (RKSM). Those methods are derived and analyzed in detail in references [22–25] and references therein. The Newton–Kleinman method is a very complex approach and at each Newton step, LRADI iterations need to be executed once, which is a very timelaborious task [26–29]. Also, the Newton–Kleinman method cannot ensure the definiteness 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 ineffective 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–33] and references therein. In this strategy, initially, we need to find a ROM of the target model implementing any conventional technique, then two generalized sparsedense Sylvester equations will be solved to find the system Gramian; hence, the required projection matrices constituted by their orthonormal columns constructed via QZdecomposition [34]. The ROM attained by the TSIA approach is stabilitypreserving and we will derive its sparsitypreserving form of it. Then, rest of the activities will follow the classical inverse projection scheme discussed previously. The validity of the proposed techniques will be investigated numerically through the transient behavior of the target models and accuracy will be justified by the norm optimality.
2. Background
Here, we discuss the basic approaches for the treatment of a generalized statespace system. Then, we briefly analyze the conversion of index2 descriptor systems into equivalent generalized systems with several cases. At the end of this section, we investigate sparsitypreserving reducedorder matrices and the corresponding reducedorder model. The generalized statespace system can be written as follows:where is nonsingular, and , , , and . The transfer function of system (7) is defined by , where . The CARE defined for system (7) is as follows:
2.1. Newton–Kleinman Method
The Newton–Kleinman method is used to convert the Riccati equation to an equivalent Lyapunov equation and to thereby solve this Lyapunov equation by a suitable strategy.
According to the Fre’chet derivative at and applying the Newton iteration with the condition , we obtain the following:
Now, we consider and , then formula (9) reduces to a generalized Lyapunov equation such as follows:
The generalized CALE (10) can be solved for by any conventional method, and the corresponding feedback matrix can be estimated [35]. The summary of the method is given in Algorithm 1.

The solution of the generalized CALE (10) can be found through the Cholesky factorization utilizing the lowrank Choleskyfactor alternative direction implicit (LRCFADI) techniques. Details of the LRCFADI techniques are provided in references [36–38] and references therein. In this approach, only the Cholesky factor of the solution in needs to store. The summary of LRCFADI techniques is given in Algorithm 2.

2.2. Rational Krylov Subspace Algorithm
The rational Krylov subspace algorithm is applicable to find the factored solution of CARE (8) without any conversion. This algorithm is defined for the LTI systems of symmetric cases. An orthogonal projector spanned by the dimensional rational Krylov subspace needs to construct for finding the an approximated lowrank CARE defined as follows:where and . The lowrank CARE (11) can be solved for by any conventional method or MATLAB care command, where is taken as a lowrank approximation of . Then, the factored solution of the solution is estimated to find the optimal feedback matrix .
The derivation and update of the RKSM techniques are given in references [39, 40] and references therein. The stepbystep RKSM approach is given in [41] and summarized in Algorithm 3.

2.3. Iterative Rational Krylov Algorithm
Iterative rational Krylov algorithm (IRKA) can be applied to construct an dimensional reducedorder modelsuch that its transfer function interpolates the original one, , at selected points in the complex plane along with selected directions. The detailed procedure of IRKA is illustrated in reference [42].
Via the IRKA approach, two projection matrices and are to be constructed; hence, by the Petrov–Galerkin condition, we can construct the reducedorder matrices in formula (12) as follows:
Using the reducedorder matrices defined in formula (13), the reducedorder CARE can be attained as formula (11) and can be solved for as previously. Then, the stabilizing feedback matrix corresponding to the ROM (12) can be estimated; hence, for stabilizing the full model (7), the approximated optimal feedback matrix can be retrieved by the scheme of reverse projection as follows:
The 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.

2.4. Conversion of the Index2 Descriptor System to Generalized System
Let us Assume is a nonsymmetric positive definite and . Initially, we consider . According to Section 3 of reference [44], the algebraic part of the system (1), can be expressed as . Then, by proper elimination and substitution, system (1) can be converted to an equivalent form:where the converted matrices are structured as follows:with . The projectors and are composed as follows:
For the general case, we assume . According to Section of reference [45], we consider the velocity can be decomposed as follows:where satisfies and Thus, the required converted system can be formed as follows:with and the converted matrices are as follows:
In addition, if both and , 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], sparsitypreserving reducedorder matrices for the index2 descriptor systems can be written as follows:
Here, the projection matrices and are composed of suitable sparsitypreserving matrixvector formulation. The ROM corresponding to the reducedorder matrices defined in (22) can be written as follows:
3. TwoSided Iterative Algorithm for Index2 Descriptor Systems
This section includes the main tasks of this work. We derive an improved version of the twosided iterative algorithm (TSIA) for index2 descriptor systems to stabilize the incompressible Navier–Stokes flow.
3.1. Formulation of the Generalized SparseDense 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 efficiently approximate the original system. So, the quality of the ROM needs to be improved by further techniques such that it satisfies the Wilson conditions for optimality [42, 47].
Now, we are to construct two generalized sparsedense Sylvester equations. To do this, we need to compute the following matrices:
Assuming and , the required sparsedense Sylvester equations can be formed as follows:
3.2. Solving the Generalized SparseDense 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 only. To generate the Krylov subspace QZdecomposition of and is to be computed, such that and , respectively. Here, and 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 can be generated as follows:where and . The updated sparse form of the desired techniques are provided in Algorithm 5.
3.3. TwoSided Iterative Algorithm to Estimate the Optimal Feedback Matrix
At first, an iterative method with minimum iterations needs to be utilized to find an initial makeshift ROM with the desired dimension. With the required matrixvector algebraic operations, the Sylvester equations defined in formula (25) need to be solved for and , respectively. Then, the approximated projector matrices and can be computed by the QRdecomposition of the matrices and , respectively. Those crude projector matrices need to be refined as and , respectively. Implementing and the ROM (23) can be acquired with the reducedorder matrices defined in formula (22); hence, the reducedorder Riccati equation can be formed as follows:
which can be solved for by MATLAB care command, and then, the reducedorder feedback matrix can be estimated as . Finally, the optimal feedback matrix for system (1) can be approximated as follows:
The whole process is summed up in Algorithm 6.

3.4. Stabilization of the Index2 Descriptor System
Plugging the optimal feedback matrix , the unstable incompressible Navier–Stokes flow given by (1) can be optimally stabilized as follows:
3.5. Norm of the Error System
Let us consider and are the transfer functions of the converted generalized system (15) and (23), respectively. Then, the associated error system can be formulated as follows:where we have taken the reducedorder matrices of formula (22) into account and the system matrices are formed as follows:
Authors in reference [49] derived the way of estimating norm of the error system (31) as follows:
Here, norm of the full model is defined as and needs to be evaluated for once, which is inconvenient for the conventional simulation solvers. Assuming is the lowrank Gramian factor of the Gramian that needs to be determined. In reference [50], a practically feasible technique is derived to overcome this situation. Then, can be written as follows:
Again, is the norm of the ROM that can be evaluated by the Gramian of the lowrank Lyapunov equation
due to the involvement of the reducedorder matrices, the Lyapunov (35) is solvable by the MATLAB library command lyap.
Finally, can be computed by the Gramian of the sparsedense Sylvester equationthat can be conveniently solved by reforming the techniques provided in Algorithm 5.
4. Numerical Results
In this section, the justification of the adaptability and efficiency 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. The 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 IntelXeon Silver 4114 CPU 2.20 GHz clock speed, 2 cores each, and 64 GB of total RAM.
4.1. Model Description
In this work, we have considered some unstable Navier–Stokes models with the Reynolds number, respectively. The 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 firstorder index2 descriptor system the submatrices and both are the sparse zero matrices. The full specifications of the target models are given in Table 1.
4.2. Approximation of the Full Models with the ReducedOrder Models
The 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 norm of the error systems for the objective models.
4.3. 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 3dimensional model with .
In Figure 1, the subfigure Figure 1(a) shows the comparison of the transfer function (sigma plot) of the full model and corresponding ROM, whereas subfigures Figures 1(b) and 1(c) evince the absolute and relative errors, respectively, of this reducedorder approximation. From the abovementioned figures, 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.
(a)
(b)
(c)
4.4. Norm of the Error System for the ROMs
Here, norm of the error system for the ROMs of the target models will be provided.
Table 2 conveys the desired error norm of the ROMs. From the measured error norm of the ROMs, it is evident that the level of the approximation process is suitably robust and confirms the compatible scale of the accuracy.
(a)
(b)
(c)
(a)
(b)
(c)
4.5. 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 3dimensional models.
4.6. Stabilization of the Eigenvalues
The subfigures of Figure 2 exhibit the eigenvalues of the original (unstable) 3dimensional systems with Reynolds number , respectively. In contrast, the subfigures of Figure 3 exhibit the stabilized eigenvalues of that of the systems. From those figures, it can be concluded that the eigenvalues of the desired models are properly stabilized.
4.7. Stabilization of the Step Responses
The subfigure of Figures 4 and 5 display the stepresponses of the 1stinput/1stoutput and 2ndinput/7thoutput relations, respectively, of the abovementioned original (unstable) systems. On the contrary, the subfigures of Figures 6 and 7 display the stabilized step responses of the 1stinput/1stoutput and 2ndinput/7thoutput relations, respectively, of that the systems. The scenario of the stepresponse stabilization duly confirms the efficiency of the proposed techniques.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
4.8. Comparison of the TSIA with IRKA
In this section, we are going to compare the proposed algorithm with the IRKA approach for reducedorder 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 stepresponses by the TSIA and IRKA via the reducedorder 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 error norm of the ROMs of selected 3dimensional models. Table 3 reveals the required information on computation time and error norm.
From Table 3, it is apparent that in the case of computation time, the TSIA algorithm requires about twothirds of that of IRKA. Also, a comparison of error norm of the ROMs manifests the betterment of TSIA over IRKA.
5. Conclusions
We have introduced a twosided projectionbased sparsitypreserving reducedorder modeling approach for the stabilization of nonsymmetric index2 descriptor systems explored from unstable Navier–Stokes models. The Wilson conditions satisfying Reducedorder models are derived by the twosided projection techniques implementing the sparsedense Sylvester equations. To solve the desired Sylvester equations, sparsitypreserving 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 proficiently applied for the stabilization of the target models through reducedorder modeling. We have efficiently approximated the full models by the corresponding ROMs with minimized error norm and stabilized the eigenvalues and stepresponses properly. A comparison of reducedorder 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 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. Thus, 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.
Data Availability
The computational data of the Navier–Stokes models are available at https://github.com/uddinmonir/systemDATA/find/main.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research work was funded by an NSUCTRG research grant under project no. : CTRG19/SEPS/05.