In this paper, we explore the Riccati-based 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 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 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 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 finite element method without altering the time variable converts them into linear time-invariant systems [3]. The incompressible Navier–Stokes flow can be written as the following differential-algebraic 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 input-output 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 index-2 differential-algebraic 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 Riccati-based 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 closed-loop 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, 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 , , , , and . The projector and 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). The reduced-order Riccati equation (6) is to be solved for a symmetric positive-definite 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 subspace-based iterative rational Krylov algorithm (IRKA). Those techniques are elaborately discussed in references [8, 1421] and references therein. From the discussion of the BT method, some obstacles are identified, 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). Those methods are derived and analyzed in detail in references [2225] and references therein. The 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 [2629]. 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 subspace-based two sided iterative algorithm (TSIA) as discussed in references [3133] and references therein. In this strategy, initially, we need to find a ROM of the target model implementing any conventional technique, then two generalized sparse-dense Sylvester equations will be solved to find the system Gramian; hence, the required projection matrices constituted by their orthonormal columns constructed via QZ-decomposition [34]. The ROM attained by the TSIA approach is stability-preserving and we will derive its sparsity-preserving 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 state-space system. Then, we briefly 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 reduced-order model. The generalized state-space 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.

Input: and (initial assumption).
Output: approximate feedback matrix
(1)While do
(2)Compute and .
(3)Solve the Lyapunov (10) for .
(4)Compute .
(5)end while.

The 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 [3638] and references therein. In this approach, only the Cholesky factor of the solution in needs to store. The summary of LRCF-ADI techniques is given in Algorithm 2.

Input: , (tolerance), (number of iterations) and shift parameters .
Output: low-rank Cholesky-factor such that .
(1)Consider , and .
(2)while or do
(3)Solve .
(4)if then
(5)Update ,
(6)Compute .
(8)Assume , ,
(9)Update ,
(10)Compute .
(12)end if
(14)end while.
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 low-rank CARE defined as follows:where and . The low-rank CARE (11) can be solved for by any conventional method or MATLAB care command, where is taken as a low-rank 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 step-by-step RKSM approach is given in [41] and summarized in Algorithm 3.

Input: (number of iterations) and (initial shifts).
Output: optimal feedback matrix
(1)Compute (QR factorization).
(2)Choose .
(3)while not converged or do
(4)Solve .
(5)Compute shift for the next iteration.
(6)Using the Arnoldi algorithm orthogonalize against to obtain , such that .
(7)Assuming , , and , for solve the Riccati (11).
(8)Compute for convergence.
(9)end while
(10)Compute eigenvalue decomposition .
(11)For negligible eigenvalues truncate and compute .
(12)Compute .
2.3. Iterative Rational Krylov Algorithm

Iterative rational Krylov algorithm (IRKA) can be applied to construct an -dimensional reduced-order 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 reduced-order matrices in formula (12) as follows:

Using the reduced-order matrices defined in formula (13), the reduced-order 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.

Output: optimal feedback matrix .
(1)Make the initial selection of the interpolation points and the tangential directions and .
(2)Construct the projection matrices , .
(3)while (not converged) do
(4)Compute the reduced-order matrices as (13).
(5)Compute and for , and , for .
(6)Repeat step 2.
(7)Repeat step 4.
(8)Solve the Riccati (11) for .
(9)Compute and hence .
(10)end while.
2.4. Conversion of the Index-2 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], sparsity-preserving reduced-order matrices for the index-2 descriptor systems can be written as follows:

Here, the projection matrices and are composed of suitable sparsity-preserving matrix-vector formulation. The ROM corresponding to the reduced-order matrices defined in (22) can be written as follows:

3. Two-Sided Iterative Algorithm for Index-2 Descriptor Systems

This 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 flow.

3.1. 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 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 sparse-dense Sylvester equations. To do this, we need to compute the following matrices:

Assuming and , the required sparse-dense Sylvester equations can be formed as follows:

3.2. 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 only. To generate the Krylov subspace QZ-decomposition 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.

Input: , (number of iterations).
Output: approximate solution of the Sylvester (25).
(1)Compute the QZ-decomposition and .
(2)Construct according to (24) and compute .
(3)Compute .
(4)while do
(5)Compute by formula (27).
(6)Solve the linear system (26) for .
(7)end while
(8)Compute .
3.3. Two-Sided 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 matrix-vector 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 QR-decomposition 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 reduced-order matrices defined in formula (22); hence, the reduced-order Riccati equation can be formed as follows:

which can be solved for by MATLAB care command, and then, the reduced-order 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.

Output: satisfying the Wilson conditions.
(1)Construct the Sylvester equations defined in (25) and solve for and , respectively.
(2)Compute the QR decomposition and .
(3)Compute and .
(4)Construct the Wilson conditions satisfying reduced-order matrices defined in (22).
(5)Form the reduced-order Riccati (28) and solve for .
(6)Approximate the optimal feedback matrix from (29).
3.4. Stabilization of the Index-2 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 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 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 low-rank 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 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, can be computed by the Gramian of the sparse-dense 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 real-world 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.

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 first-order index-2 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 Reduced-Order 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 3-dimensional 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 reduced-order 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.

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.

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 3-dimensional models.

4.6. Stabilization of the Eigenvalues

The subfigures of Figure 2 exhibit the eigenvalues of the original (unstable) 3-dimensional 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 step-responses of the 1st-input/1st-output and 2nd-input/7th-output 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 1st-input/1st-output and 2nd-input/7th-output relations, respectively, of that the systems. The scenario of the step-response stabilization duly confirms the efficiency of the proposed techniques.

4.8. 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 error norm of the ROMs of selected 3-dimensional 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 two-thirds 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 two-sided projection-based sparsity-preserving reduced-order modeling approach for the stabilization of nonsymmetric index-2 descriptor systems explored from unstable Navier–Stokes models. The 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 matrix-vector operations. From the numerical analyses, it is ascertained that the proposed techniques can be proficiently applied for the stabilization of the target models through reduced-order modeling. We have efficiently approximated the full models by the corresponding ROMs with minimized 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 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.


This research work was funded by an NSU-CTRG research grant under project no. : CTRG-19/SEPS/05.