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

Authors: Md. Toriqul Islam; Mahtab Uddin; M. Monir Uddin; Md. Abdul Hakim Khan; Md. Tanzim Hossain
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

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 ℋ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 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.


---

## Body

## 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:(1)Mv˙t=A1vt+A2pt+B1ut,0=A3vt+B2ut,yt=C1vt+C2pt+Dut,where vt∈ℝnv is the nodal vector of the discretized velocity with v0=v0, pt∈ℝnp is the discretized pressure, and ut∈ℝnu is the input with the output yt∈ℝny. The symmetric positive definite matrix M∈ℝnv×nv represents the mass, whereas the matrices A1∈ℝnv×nv, A2∈ℝnv×np, and A3∈ℝnp×nv are for system components including discretized gradient. B1∈ℝnv×nu and B2∈ℝnp×nu are the input matrices, whereas C1∈ℝny×nv and C2∈ℝny×np are the output matrices. The input-output matrices quantify the velocity pattern at the corresponding nodes. D∈ℝny×nu is the direct transform matrix [4]. Considering M is invertible, A2 and A3 have both full column ranks and A3M−1A2 is nonsingular. Hence, system (1) is of the index-2 differential-algebraic system having the dimension nv+np [5]. A system equivalent to system (1) can be written in the following form:(2)E¯x˙t=A¯xt+B¯ut,yt=C¯xt+D¯ut,where xt=vtpt, E¯=M000, A¯=A1A2A30, B¯=B1B2, C¯=C1C2, and D¯=D with appropriate dimensions. The matrix pencil corresponding to the system is defined as A¯,E¯, where the matrix E¯ is singular and the matrix pencil A¯,E¯ has nv−np number of nonzero finite complex eigenvalues and 2np infinite eigenvalues [6].Reynolds numberRe is one of the key features that instigate the system characteristics. For Re≥300, 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:(3)A¯TXE¯+E¯TXA¯−E¯TXB¯B¯TXE¯+C¯TC¯=0.The solutionX of the Riccati equation (3) is symmetric positive definite and called stabilizing for a stable closed-loop matrix A¯−B¯B¯TXE¯. The optimal feedback matrix Ko can be computed as Ko=B¯TXE¯ and can be implemented for optimally stabilizing the target system by replacing A¯ by A¯s=A¯−B¯Ko [9]. The stabilized system can be written as follows:(4)E¯x˙t=A¯sxt+B¯ut,yt=C¯xt+D¯ut.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:(5)E¯^x˙^t=A¯^x^t+B¯^u^t,y^t=C¯^x^t+D¯^u^t.where E¯^=WTE¯V, A¯^=WTA¯V, B¯^=WTB¯, C¯^=C¯V, and D¯^=D¯. The 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:(6)A¯^TX^E¯^+E¯^TX^A¯^−E¯^TX^B¯^B¯^TX^E¯^+C¯^TC¯^=0.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 X^ by the MATLAB library command care. Then, the stabilizing feedback matrix K^=B¯^TX^E¯^ corresponding to the ROM (5) can be estimated, and hence, the approximated optimal feedback matrix Ko=K^VTE¯ [13]. Finally, utilizing Ko 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, 14–21] 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 [22–25] 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 [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 subspace-based 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 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 ℋ2 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:(7)Ex˙t=Axt+But,yt=Cxt+Dut,where E∈ℝn×n is nonsingular, and A∈ℝn×n, B∈ℝn×p, C∈ℝm×n, and D∈ℝm×p. The transfer function of system (7) is defined by Gs=CsE−A−1B+D, where s∈ℂ. The CARE defined for system (7) is as follows:(8)ATXE+ETXA−ETXBBTXE+CTC=0.


### 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 atX=Xi and applying the Newton iteration with the condition ΔXi=Xi+1−Xi, we obtain the following:(9)A−BBTXiETXi+1E+ETXi+1A−BBTXiE=−CTC−ETXiBBTXiE.Now, we considerA˜i=A−BBTXiE and Wi=CTETXiB, then formula (9) reduces to a generalized Lyapunov equation such as follows:(10)A˜iTXi+1E+ETXi+1A˜i=−WiWiT.The generalized CALE (10) can be solved forXi+1 by any conventional method, and the corresponding feedback matrix Ki+1=BTXi+1E can be estimated [35]. The summary of the method is given in Algorithm 1.Algorithm 1:Basic Kleinman–Newton method
Input:E,A,B,C and X0 (initial assumption).Output: approximate feedback matrixK(1)
Whilei≤imax do(2)
ComputeA˜i=A−BBTXiE and Wi=CTETXiB.(3)
Solve the Lyapunov (10) for Xi+1.(4)
ComputeKi+1=BTXi+1E.(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 [36–38] and references therein. In this approach, only the Cholesky factor Z of the solution X in needs to store. The summary of LRCF-ADI techniques is given in Algorithm 2.Algorithm 2: LRCF-ADI for generalized systems
Input:E,A,C, τ (tolerance), imax (number of iterations) and shift parameters μjj=1imax.Output: low-rank Cholesky-factorZ such that X≈ZZT.(1)
ConsiderW0=CT,Z0=, and i=1.(2)
whileWi−1Wi−1T≥τ or i≤imax do(3)
SolveVi=A˜T+μiET−1Wi−1.(4)
ifImμi=0 then(5)
UpdateZi=Zi−1−2μiVi,(6)
ComputeWi=Wi−1−2μiETVi.(7)
else(8)
Assumeγi=−2Reμi, δi=Reμi/Imμi,(9)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi,(10)
ComputeWi+1=Wi−1+2γi2ETReVi+δiImVi.(11)
i=i+1(12)
end if(13)
i=i+1(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 projectorV∈ℝn×m spanned by the m-dimensional rational Krylov subspace needs to construct for finding the an approximated low-rank CARE defined as follows:(11)A^TX^E^+E^TX^A^−E^TX^B^B^TX^E^+C^TC^=0,where X^=VTXV,E^=VTEV,A^=VTAV,B^=VTB and C^=CV. The low-rank CARE (11) can be solved for X^∈ℝm×m by any conventional method or MATLAB care command, where X^ is taken as a low-rank approximation of X. Then, the factored solution Z of the solution X is estimated to find the optimal feedback matrix Ko.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.Algorithm 3: RKSM for generalized systems.
Input:E,A,B,C,imax (number of iterations) and μi (initial shifts).Output: optimal feedback matrixKo(1)
ComputeQ0R0=CT (QR factorization).(2)
ChooseV0=Q0.(3)
while not converged orm≤imax do(4)
Solvevm+1=AT−μm+1ET−1Vm.(5)
Compute shift for the next iteration.(6)
Using the Arnoldi algorithm orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vm,v^m+1.(7)
AssumingE^=Vm+1TEVm+1, A^=Vm+1TAVm+1, B^=Vm+1TB and C^=CVm+1, for X^ solve the Riccati (11).(8)
Computeℛmrelative for convergence.(9)
end while(10)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(11)
For negligible eigenvalues truncateΛ2 and compute Z=Vm+1T1Λ11/2.(12)
ComputeKo=BTZZTE.


### 2.3. Iterative Rational Krylov Algorithm

Iterative rational Krylov algorithm (IRKA) can be applied to construct anr-dimensional r≪n reduced-order model(12)E^x˙^t=A^x^t+B^u^t,y^t=C^x^t+D^u^t,such that its transfer function G^s=C^sE^−A^−1B^+D^ interpolates the original one, Gs, 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 matricesV and W are to be constructed; hence, by the Petrov–Galerkin condition, we can construct the reduced-order matrices in formula (12) as follows:(13)E^=WTEV,A^=WTAV,B^=WTB,C^=CV,D^=D.Using the reduced-order matrices defined in formula (13), the reduced-order CARE can be attained as formula (11) and can be solved for X^ as previously. Then, the stabilizing feedback matrix K^=B^TX^E^ corresponding to the ROM (12) can be estimated; hence, for stabilizing the full model (7), the approximated optimal feedback matrix Ko can be retrieved by the scheme of reverse projection as follows:(14)Ko=B^TX^E^VTE=K^VTE.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.Algorithm 4: IRKA for generalized systems.
Input:E,A,B,C,DOutput: optimal feedback matrixKo.(1)
Make the initial selection of the interpolation pointsαii=1r and the tangential directions bii=1r and cii=1r.(2)
Construct the projection matricesV=α1E−A−1Bb1,…,αrE−A−1Bbr, W=α1E−A−TCTc1,…,αrE−A−TCTcr.(3)
while (not converged) do(4)
Compute the reduced-order matrices as (13).(5)
ComputeA^zi=λiE^zi and yi∗A^=λiyi∗E^ for αi←−λi, bi∗←−yi∗B^ and ci∗←C^zi∗, for i=1,…,r.(6)
Repeat step 2.(7)
Repeat step 4.(8)
Solve the Riccati (11) for X^.(9)
ComputeK^=B^TX^E^ and hence Ko=K^VTE.(10)
end while.


### 2.4. Conversion of the Index-2 Descriptor System to Generalized System

Let us AssumeM is a nonsymmetric positive definite and A2T≠A3. Initially, we consider B2=0. According to Section 3 of reference [44], the algebraic part of the system (1), pt can be expressed as vt. Then, by proper elimination and substitution, system (1) can be converted to an equivalent form:(15)Ex˙t=Axt+ℬut,yt=Cxt+Dut,where the converted matrices are structured as follows:(16)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1,C=C1−C2A3M−1A2−1A3M−1A1Πr,D=D−C2A3M−1A2−1A3M−1B1,with ΠlTv0=ΠlTv0. The projectors Πl and Πr are composed as follows:(17)Πl=I−M−1A2A3M−1A2−1A3,Πr=I−A2A3M−1A2−1A3M−1.For the general case, we assumeB2≠0. According to Section of reference [45], we consider the velocity vt can be decomposed as follows:(18)vt=v0t+vgt,where v0t satisfies A3v0t=0 and vgt=M−1A2A3M−1A2−1B2ut. Thus, the required converted system can be formed as follows:(19)Ex˙t=Axt+ℬut,yt=Cxt+Dut−C2A3M−1A2−1B2ut,with ΠlTv0=ΠlTv0−vg0 and the converted matrices are as follows:(20)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1−A1M−1A2A3M−1A2−1B2,C=C1−C2A3M−1A2−1A3M−1A1Πr,D=D−C2A3M−1A2−1A3M−1B1.In addition, if bothB2=0 and C2=0, the converted system will be the same as a system (15) but converted matrices will have the following form:(21)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1,C=C1Πr,D=D.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:(22)E˜=WTMV,A˜=WTA1V,B˜=WTB1−WTA1M−1A2A3M−1A2−1B2,C˜=C1V−C2A3M−1A2−1A3M−1A1V,D˜=D−C2A3M−1A2−1A3M−1B1.Here, the projection matricesV and W 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:(23)E˜x˙t=A˜xt+ℬ˜ut,yt=C˜xt+D˜ut.


## 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 atX=Xi and applying the Newton iteration with the condition ΔXi=Xi+1−Xi, we obtain the following:(9)A−BBTXiETXi+1E+ETXi+1A−BBTXiE=−CTC−ETXiBBTXiE.Now, we considerA˜i=A−BBTXiE and Wi=CTETXiB, then formula (9) reduces to a generalized Lyapunov equation such as follows:(10)A˜iTXi+1E+ETXi+1A˜i=−WiWiT.The generalized CALE (10) can be solved forXi+1 by any conventional method, and the corresponding feedback matrix Ki+1=BTXi+1E can be estimated [35]. The summary of the method is given in Algorithm 1.Algorithm 1:Basic Kleinman–Newton method
Input:E,A,B,C and X0 (initial assumption).Output: approximate feedback matrixK(1)
Whilei≤imax do(2)
ComputeA˜i=A−BBTXiE and Wi=CTETXiB.(3)
Solve the Lyapunov (10) for Xi+1.(4)
ComputeKi+1=BTXi+1E.(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 [36–38] and references therein. In this approach, only the Cholesky factor Z of the solution X in needs to store. The summary of LRCF-ADI techniques is given in Algorithm 2.Algorithm 2: LRCF-ADI for generalized systems
Input:E,A,C, τ (tolerance), imax (number of iterations) and shift parameters μjj=1imax.Output: low-rank Cholesky-factorZ such that X≈ZZT.(1)
ConsiderW0=CT,Z0=, and i=1.(2)
whileWi−1Wi−1T≥τ or i≤imax do(3)
SolveVi=A˜T+μiET−1Wi−1.(4)
ifImμi=0 then(5)
UpdateZi=Zi−1−2μiVi,(6)
ComputeWi=Wi−1−2μiETVi.(7)
else(8)
Assumeγi=−2Reμi, δi=Reμi/Imμi,(9)
UpdateZi+1=Zi−1γiReVi+δiImViγiδi2+1ImVi,(10)
ComputeWi+1=Wi−1+2γi2ETReVi+δiImVi.(11)
i=i+1(12)
end if(13)
i=i+1(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 projectorV∈ℝn×m spanned by the m-dimensional rational Krylov subspace needs to construct for finding the an approximated low-rank CARE defined as follows:(11)A^TX^E^+E^TX^A^−E^TX^B^B^TX^E^+C^TC^=0,where X^=VTXV,E^=VTEV,A^=VTAV,B^=VTB and C^=CV. The low-rank CARE (11) can be solved for X^∈ℝm×m by any conventional method or MATLAB care command, where X^ is taken as a low-rank approximation of X. Then, the factored solution Z of the solution X is estimated to find the optimal feedback matrix Ko.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.Algorithm 3: RKSM for generalized systems.
Input:E,A,B,C,imax (number of iterations) and μi (initial shifts).Output: optimal feedback matrixKo(1)
ComputeQ0R0=CT (QR factorization).(2)
ChooseV0=Q0.(3)
while not converged orm≤imax do(4)
Solvevm+1=AT−μm+1ET−1Vm.(5)
Compute shift for the next iteration.(6)
Using the Arnoldi algorithm orthogonalizevm+1 against Vm to obtain v^m+1, such that Vm+1=Vm,v^m+1.(7)
AssumingE^=Vm+1TEVm+1, A^=Vm+1TAVm+1, B^=Vm+1TB and C^=CVm+1, for X^ solve the Riccati (11).(8)
Computeℛmrelative for convergence.(9)
end while(10)
Compute eigenvalue decompositionX^=TΛTT=T1T2Λ100Λ2T1TT2T.(11)
For negligible eigenvalues truncateΛ2 and compute Z=Vm+1T1Λ11/2.(12)
ComputeKo=BTZZTE.


## 2.3. Iterative Rational Krylov Algorithm

Iterative rational Krylov algorithm (IRKA) can be applied to construct anr-dimensional r≪n reduced-order model(12)E^x˙^t=A^x^t+B^u^t,y^t=C^x^t+D^u^t,such that its transfer function G^s=C^sE^−A^−1B^+D^ interpolates the original one, Gs, 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 matricesV and W are to be constructed; hence, by the Petrov–Galerkin condition, we can construct the reduced-order matrices in formula (12) as follows:(13)E^=WTEV,A^=WTAV,B^=WTB,C^=CV,D^=D.Using the reduced-order matrices defined in formula (13), the reduced-order CARE can be attained as formula (11) and can be solved for X^ as previously. Then, the stabilizing feedback matrix K^=B^TX^E^ corresponding to the ROM (12) can be estimated; hence, for stabilizing the full model (7), the approximated optimal feedback matrix Ko can be retrieved by the scheme of reverse projection as follows:(14)Ko=B^TX^E^VTE=K^VTE.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.Algorithm 4: IRKA for generalized systems.
Input:E,A,B,C,DOutput: optimal feedback matrixKo.(1)
Make the initial selection of the interpolation pointsαii=1r and the tangential directions bii=1r and cii=1r.(2)
Construct the projection matricesV=α1E−A−1Bb1,…,αrE−A−1Bbr, W=α1E−A−TCTc1,…,αrE−A−TCTcr.(3)
while (not converged) do(4)
Compute the reduced-order matrices as (13).(5)
ComputeA^zi=λiE^zi and yi∗A^=λiyi∗E^ for αi←−λi, bi∗←−yi∗B^ and ci∗←C^zi∗, for i=1,…,r.(6)
Repeat step 2.(7)
Repeat step 4.(8)
Solve the Riccati (11) for X^.(9)
ComputeK^=B^TX^E^ and hence Ko=K^VTE.(10)
end while.


## 2.4. Conversion of the Index-2 Descriptor System to Generalized System

Let us AssumeM is a nonsymmetric positive definite and A2T≠A3. Initially, we consider B2=0. According to Section 3 of reference [44], the algebraic part of the system (1), pt can be expressed as vt. Then, by proper elimination and substitution, system (1) can be converted to an equivalent form:(15)Ex˙t=Axt+ℬut,yt=Cxt+Dut,where the converted matrices are structured as follows:(16)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1,C=C1−C2A3M−1A2−1A3M−1A1Πr,D=D−C2A3M−1A2−1A3M−1B1,with ΠlTv0=ΠlTv0. The projectors Πl and Πr are composed as follows:(17)Πl=I−M−1A2A3M−1A2−1A3,Πr=I−A2A3M−1A2−1A3M−1.For the general case, we assumeB2≠0. According to Section of reference [45], we consider the velocity vt can be decomposed as follows:(18)vt=v0t+vgt,where v0t satisfies A3v0t=0 and vgt=M−1A2A3M−1A2−1B2ut. Thus, the required converted system can be formed as follows:(19)Ex˙t=Axt+ℬut,yt=Cxt+Dut−C2A3M−1A2−1B2ut,with ΠlTv0=ΠlTv0−vg0 and the converted matrices are as follows:(20)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1−A1M−1A2A3M−1A2−1B2,C=C1−C2A3M−1A2−1A3M−1A1Πr,D=D−C2A3M−1A2−1A3M−1B1.In addition, if bothB2=0 and C2=0, the converted system will be the same as a system (15) but converted matrices will have the following form:(21)x=x1,E=ΠlMΠr,A=ΠlA1Πr,ℬ=ΠlB1,C=C1Πr,D=D.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:(22)E˜=WTMV,A˜=WTA1V,B˜=WTB1−WTA1M−1A2A3M−1A2−1B2,C˜=C1V−C2A3M−1A2−1A3M−1A1V,D˜=D−C2A3M−1A2−1A3M−1B1.Here, the projection matricesV and W 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:(23)E˜x˙t=A˜xt+ℬ˜ut,yt=C˜xt+D˜ut.


## 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 ℋ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:(24)Bn=MA2A30−1B1B2,Cn=MTA3TA2T0−1C1TC2T.AssumingPn=MBnℬ˜T and Qn=−MCnC˜, the required sparse-dense Sylvester equations can be formed as follows:(25)AXE˜+EXA˜+Pn=0,AYE˜+EYA˜+Qn=0.


### 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 X only. To generate the Krylov subspace QZ-decomposition of E˜ and A˜ is to be computed, such that E˜=QSZT and A˜=QTZT, 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:(26)SjjA1+TjjMTSjjA2SjjA30X˜j=−Pn^j,and the successive columns for the solution matrix X can be generated as follows:(27)Pn^j=−Pn˜j−∑i=1j−1FijX˜i,where Pn˜j=PnjZ and Fij=SijA1+TijMTSijA2SijA30. The updated sparse form of the desired techniques are provided in Algorithm 5.Algorithm 5: Solution of generalized sparse-dense Sylvester equation.
Input:M,A1,A2,A3,B1,B2,E˜,A˜,ℬ˜, jmax (number of iterations).Output: approximate solutionX of the Sylvester (25).(1)
Compute the QZ-decompositionE˜=QSZT and A˜=QTZT.(2)
ConstructBn according to (24) and compute Pn=MBnℬ˜T.(3)
ComputePn˜=PnZ.(4)
whilej≤jmax do(5)
ComputePn^j by formula (27).(6)
Solve the linear system (26) for X˜j.(7)
end while(8)
ComputeX=X˜QT.


### 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 X and Y, respectively. Then, the approximated projector matrices V˜ and W˜ can be computed by the QR-decomposition of the matrices X and Y, respectively. Those crude projector matrices need to be refined as V=V˜ and W=W˜V˜TMTW˜−1, respectively. Implementing V and W 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:(28)A˜TX˜E˜+E˜TX˜A˜−E˜TX˜ℬ˜ℬ˜TX˜E˜+C˜TC˜=0,which can be solved forX˜ by MATLAB care command, and then, the reduced-order feedback matrix can be estimated as K˜=ℬ˜TX˜E˜. Finally, the optimal feedback matrix for system (1) can be approximated as follows:(29)Ko=ℬ˜TX˜E˜VTE˜=K˜VTM.The whole process is summed up in Algorithm6.Algorithm 6: TSIA for the optimal feedback matrix of index-2 descriptor systems.
Input:M,A1,A2,A3,B1,B2,C1,C2,E˜,A˜,ℬ˜,C˜Output:E˜,A˜,ℬ˜,C˜ satisfying the Wilson conditions.(1)
Construct the Sylvester equations defined in (25) and solve for X and Y, respectively.(2)
Compute the QR decompositionX=V˜βv and Y=W˜βw.(3)
ComputeV=V˜ and W=W˜V˜TMTW˜−1.(4)
Construct the Wilson conditions satisfying reduced-order matrices defined in (22).(5)
Form the reduced-order Riccati (28) and solve for X˜.(6)
Approximate the optimal feedback matrix from (29).


### 3.4. Stabilization of the Index-2 Descriptor System

Plugging the optimal feedback matrixKo, the unstable incompressible Navier–Stokes flow given by (1) can be optimally stabilized as follows:(30)Mv˙t=A1−B1Kovt+A2pt+B1ut,0=A3vt+B2ut,yt=C1vt+C2pt+Dut.


### 3.5.ℋ2 Norm of the Error System

Let us considerGs and G˜s are the transfer functions of the converted generalized system (15) and (23), respectively. Then, the associated error system can be formulated as follows:(31)Gerr=Gs−G˜s=CerrsEerr−Aerr−1ℬerr,where we have taken the reduced-order matrices of formula (22) into account and the system matrices are formed as follows:(32)Eerr=E00E˜,Aerr=A00A˜,ℬerr=ℬℬ˜,andCerr=C−C⌢.Authors in reference [49] derived the way of estimating ℋ2 norm of the error system (31) as follows:(33)Gerrℋ2=Gsℋ22+G˜sℋ22+2traceℬTQsℬ˜.Here,ℋ2 norm of the full model is defined as Gsℋ2 and needs to be evaluated for once, which is inconvenient for the conventional simulation solvers. Assuming Zq 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. Then, Gsℋ2 can be written as follows:(34)Gsℋ2=traceB1TZqB1TZqT+B2TZqB2TZqT.Again,G˜sℋ2 is the ℋ2 norm of the ROM that can be evaluated by the Gramian Q˜ of the low-rank Lyapunov equation(35)A˜TQ˜E˜+E˜TQ˜A˜+C˜TC˜=0,due to the involvement of the reduced-order matrices, the Lyapunov (35) is solvable by the MATLAB library command lyap.Finally,traceℬTQsℬ˜ can be computed by the Gramian Qs of the sparse-dense Sylvester equation(36)ATQsE˜+ETQsA˜+CTC˜=0,that can be conveniently solved by reforming the techniques provided in Algorithm 5.


## 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 ℋ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:(24)Bn=MA2A30−1B1B2,Cn=MTA3TA2T0−1C1TC2T.AssumingPn=MBnℬ˜T and Qn=−MCnC˜, the required sparse-dense Sylvester equations can be formed as follows:(25)AXE˜+EXA˜+Pn=0,AYE˜+EYA˜+Qn=0.


## 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 X only. To generate the Krylov subspace QZ-decomposition of E˜ and A˜ is to be computed, such that E˜=QSZT and A˜=QTZT, 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:(26)SjjA1+TjjMTSjjA2SjjA30X˜j=−Pn^j,and the successive columns for the solution matrix X can be generated as follows:(27)Pn^j=−Pn˜j−∑i=1j−1FijX˜i,where Pn˜j=PnjZ and Fij=SijA1+TijMTSijA2SijA30. The updated sparse form of the desired techniques are provided in Algorithm 5.Algorithm 5: Solution of generalized sparse-dense Sylvester equation.
Input:M,A1,A2,A3,B1,B2,E˜,A˜,ℬ˜, jmax (number of iterations).Output: approximate solutionX of the Sylvester (25).(1)
Compute the QZ-decompositionE˜=QSZT and A˜=QTZT.(2)
ConstructBn according to (24) and compute Pn=MBnℬ˜T.(3)
ComputePn˜=PnZ.(4)
whilej≤jmax do(5)
ComputePn^j by formula (27).(6)
Solve the linear system (26) for X˜j.(7)
end while(8)
ComputeX=X˜QT.


## 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 X and Y, respectively. Then, the approximated projector matrices V˜ and W˜ can be computed by the QR-decomposition of the matrices X and Y, respectively. Those crude projector matrices need to be refined as V=V˜ and W=W˜V˜TMTW˜−1, respectively. Implementing V and W 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:(28)A˜TX˜E˜+E˜TX˜A˜−E˜TX˜ℬ˜ℬ˜TX˜E˜+C˜TC˜=0,which can be solved forX˜ by MATLAB care command, and then, the reduced-order feedback matrix can be estimated as K˜=ℬ˜TX˜E˜. Finally, the optimal feedback matrix for system (1) can be approximated as follows:(29)Ko=ℬ˜TX˜E˜VTE˜=K˜VTM.The whole process is summed up in Algorithm6.Algorithm 6: TSIA for the optimal feedback matrix of index-2 descriptor systems.
Input:M,A1,A2,A3,B1,B2,C1,C2,E˜,A˜,ℬ˜,C˜Output:E˜,A˜,ℬ˜,C˜ satisfying the Wilson conditions.(1)
Construct the Sylvester equations defined in (25) and solve for X and Y, respectively.(2)
Compute the QR decompositionX=V˜βv and Y=W˜βw.(3)
ComputeV=V˜ and W=W˜V˜TMTW˜−1.(4)
Construct the Wilson conditions satisfying reduced-order matrices defined in (22).(5)
Form the reduced-order Riccati (28) and solve for X˜.(6)
Approximate the optimal feedback matrix from (29).


## 3.4. Stabilization of the Index-2 Descriptor System

Plugging the optimal feedback matrixKo, the unstable incompressible Navier–Stokes flow given by (1) can be optimally stabilized as follows:(30)Mv˙t=A1−B1Kovt+A2pt+B1ut,0=A3vt+B2ut,yt=C1vt+C2pt+Dut.


## 3.5.ℋ2 Norm of the Error System

Let us considerGs and G˜s are the transfer functions of the converted generalized system (15) and (23), respectively. Then, the associated error system can be formulated as follows:(31)Gerr=Gs−G˜s=CerrsEerr−Aerr−1ℬerr,where we have taken the reduced-order matrices of formula (22) into account and the system matrices are formed as follows:(32)Eerr=E00E˜,Aerr=A00A˜,ℬerr=ℬℬ˜,andCerr=C−C⌢.Authors in reference [49] derived the way of estimating ℋ2 norm of the error system (31) as follows:(33)Gerrℋ2=Gsℋ22+G˜sℋ22+2traceℬTQsℬ˜.Here,ℋ2 norm of the full model is defined as Gsℋ2 and needs to be evaluated for once, which is inconvenient for the conventional simulation solvers. Assuming Zq 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. Then, Gsℋ2 can be written as follows:(34)Gsℋ2=traceB1TZqB1TZqT+B2TZqB2TZqT.Again,G˜sℋ2 is the ℋ2 norm of the ROM that can be evaluated by the Gramian Q˜ of the low-rank Lyapunov equation(35)A˜TQ˜E˜+E˜TQ˜A˜+C˜TC˜=0,due to the involvement of the reduced-order matrices, the Lyapunov (35) is solvable by the MATLAB library command lyap.Finally,traceℬTQsℬ˜ can be computed by the Gramian Qs of the sparse-dense Sylvester equation(36)ATQsE˜+ETQsA˜+CTC˜=0,that 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 numberRe=300,400,500, 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 B2 and C2 both are the sparse zero matrices. The full specifications of the target models are given in Table 1.Table 1
Structure of the target Navier–Stokes models.
DimensionReynolds number (Re)StatesAlgebraic variablesInputOutput13003142453400500230085681123400500330019770261527400500430044744578340050053009805412566400500


### 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ℋ2-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 withRe=500.In Figure1, 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.Figure 1
Comparison of full model and ROM of 3-dimensional model forRe=500. (a) Sigma plot. (b) Absolute error. (c) Relative error.
(a)(b)(c)


### 4.4.ℋ2 Norm of the Error System for the ROMs

Here,ℋ2 norm of the error system for the ROMs of the target models will be provided.Table2 conveys the desired ℋ2 error norm of the ROMs. From the measured ℋ2 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.Table 2
ℋ2 error norm of the ROMs of the target models.
Reynolds NumberDimensionℋ2 error norm30012.11×10−0327.52×10−0437.90×10−0448.12×10−0458.49×10−0440011.55×10−0223.25×10−0333.57×10−0343.95×10−0354.15×10−0350012.62×10−0121.22×10−0231.74×10−0241.88×10−0259.46×10−03Figure 2
Unstable eigenvalues of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 3
Stabilized eigenvalues of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(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 3-dimensional models.


### 4.6. Stabilization of the Eigenvalues

The subfigures of Figure2 exhibit the eigenvalues of the original (unstable) 3-dimensional systems with Reynolds number Re=300,400,500, 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 Figures4 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.Figure 4
Unstable step response for 1st input/1st output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 5
Unstable step response for 2nd input/7th output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 6
Stabilized step response for 1st input/1st output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 7
Stabilized step response for 2nd input/7th output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(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 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ℋ2 error norm of the ROMs of selected 3-dimensional models. Table 3 reveals the required information on computation time and ℋ2 error norm.Table 3
Comparison of computation time andℋ2 error norm of the ROMs of 3-dimensional models archived by TSIA and IRKA.
Reynolds NumberTimeℋ2error normTSIAIRKATSIAIRKA3006.81×1039.83×1037.90×10−049.97×10−044006.86×1039.65×1033.57×10−036.30×10−035006.94×1039.68×1031.74×10−025.13×10−02From 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ℋ2 error norm of the ROMs manifests the betterment of TSIA over IRKA.


## 4.1. Model Description

In this work, we have considered some unstable Navier–Stokes models with the Reynolds numberRe=300,400,500, 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 B2 and C2 both are the sparse zero matrices. The full specifications of the target models are given in Table 1.Table 1
Structure of the target Navier–Stokes models.
DimensionReynolds number (Re)StatesAlgebraic variablesInputOutput13003142453400500230085681123400500330019770261527400500430044744578340050053009805412566400500


## 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ℋ2-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 withRe=500.In Figure1, 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.Figure 1
Comparison of full model and ROM of 3-dimensional model forRe=500. (a) Sigma plot. (b) Absolute error. (c) Relative error.
(a)(b)(c)


## 4.4.ℋ2 Norm of the Error System for the ROMs

Here,ℋ2 norm of the error system for the ROMs of the target models will be provided.Table2 conveys the desired ℋ2 error norm of the ROMs. From the measured ℋ2 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.Table 2
ℋ2 error norm of the ROMs of the target models.
Reynolds NumberDimensionℋ2 error norm30012.11×10−0327.52×10−0437.90×10−0448.12×10−0458.49×10−0440011.55×10−0223.25×10−0333.57×10−0343.95×10−0354.15×10−0350012.62×10−0121.22×10−0231.74×10−0241.88×10−0259.46×10−03Figure 2
Unstable eigenvalues of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 3
Stabilized eigenvalues of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(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 3-dimensional models.


## 4.6. Stabilization of the Eigenvalues

The subfigures of Figure2 exhibit the eigenvalues of the original (unstable) 3-dimensional systems with Reynolds number Re=300,400,500, 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 Figures4 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.Figure 4
Unstable step response for 1st input/1st output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 5
Unstable step response for 2nd input/7th output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 6
Stabilized step response for 1st input/1st output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(a)(b)(c)Figure 7
Stabilized step response for 2nd input/7th output of 3-dimensional models. (a) Reynolds number 300. (b) Reynolds number 400. (c) Reynolds number 500.
(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 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ℋ2 error norm of the ROMs of selected 3-dimensional models. Table 3 reveals the required information on computation time and ℋ2 error norm.Table 3
Comparison of computation time andℋ2 error norm of the ROMs of 3-dimensional models archived by TSIA and IRKA.
Reynolds NumberTimeℋ2error normTSIAIRKATSIAIRKA3006.81×1039.83×1037.90×10−049.97×10−044006.86×1039.65×1033.57×10−036.30×10−035006.94×1039.68×1031.74×10−025.13×10−02From 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ℋ2 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ℋ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 ℋ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. 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.



---

*Source: 4435167-2022-11-22.xml*