An Iterative Model Reduction Scheme for Quadratic-Bilinear Descriptor Systems with an Application to Navier-Stokes Equations

TL;DR

This paper develops an efficient iterative model reduction method for quadratic-bilinear descriptor systems, extending interpolation-based techniques to these systems and demonstrating effectiveness on Navier-Stokes equations.

Contribution

It introduces a transformation approach to apply QBODE reduction techniques to descriptor systems, avoiding explicit projector computation for improved efficiency.

Findings

Successfully reduces Navier-Stokes based systems

Produces near-optimal reduced models

Enhances model reduction reliability for descriptor systems

Abstract

We discuss model reduction for a particular class of quadratic-bilinear (QB) descriptor systems. The main goal of this article is to extend the recently studied interpolation-based optimal model reduction framework for QBODEs [Benner et al. '16] to a class of descriptor systems in an efficient and reliable way. Recently, it has been shown in the case of linear or bilinear systems that a direct extension of interpolation-based model reduction techniques to descriptor systems, without any modifications, may lead to poor reduced-order systems. Therefore, for the analysis, we aim at transforming the considered QB descriptor system into an equivalent QBODE system by means of projectors for which standard model reduction techniques for QBODEs can be employed, including aforementioned interpolation scheme. Subsequently, we discuss related computational issues, thus resulting in a modified algorithm that allows us to construct \emph{near}--optimal reduced-order systems without explicitly computing the projectors used in the analysis. The efficiency of the proposed algorithm is illustrated by means of a numerical example, obtained via semi-discretization of the Navier-Stokes equations.

Figures

Equations

$$E_{11}\dot{v}(t)$$

$$0$$

$$y(t)$$

$$[left=\Sigma :\text{empheqlbrace}]\dot{x}(t)$$

$$y(t)$$

$$[left=\hat{\Sigma }:\text{empheqlbrace}]\dot{\hat{x}}(t)$$

$$\hat{y}(t)$$

$$W^T\left(V\dot{\hat{x}}(t)-AV\hat{x}(t)+HV\hat{x}(t)\otimes V\hat{x}(t)+\sum _{k=1}^m{\hat{N}}_k\hat{x}(t)u_k(t)+Bu(t)\right)=0.$$

$$\hat{A}$$

$$\hat{H}$$

$$\mathrm{vec}\left(XYZ\right)=(Z^T\otimes X)\mathrm{vec}\left(Y\right),$$

$$j=1+\sum _{m=1,m\ne k}^m(i_m-1)J_m,\text{where}{J}_m:=\prod _{g=1,g\ne m}^{l-1}{I}_g.$$

$$\Vert \Sigma {\Vert }_{{\mathcal{H}}_2}:=\sqrt{\mathrm{tr}\left(\sum _{i=1}^{\infty }{\int }_0^{\infty }\cdots {\int }_0^{\infty }{f}_i(t_1,\dots ,t_i)f_i(t_1,\dots ,t_i{)}^T\right)},$$

$$f_i(t_1,\dots ,t_i)=C{g}_i(t_1,\dots ,t_i),$$

$$g_1(t_1)$$

$$g_2(t_1,t_2)$$

$$g_i(t_1,\dots ,t_i)$$

$$\displaystyle~{}~{}\qquad\begin{bmatrix}N_{1},\ldots,N_{m}\end{bmatrix}\left(I_{m}\otimes g_{i-1}(t_{1},\ldots,t_{i-1})\right)\big{]},\qquad i\geq 3.$$

$$\Vert \Sigma {\Vert }_{{\mathcal{H}}_2}=\sqrt{\mathrm{tr}\left(CP{C}^T\right)}=\sqrt{\mathrm{tr}\left(B^{T}QB\right)},$$

$$AP+P{A}^T+H(P\otimes P)H^T+\sum _{k=1}^m{N}_{k}P{N}_k^T+B{B}^T=0,$$

$$A^{T}Q+QA+H^{(2)}(P\otimes Q){\left(H^{(2)}\right)}^T+\sum _{k=1}^m{N}_k^{T}Q{N}_k+C^{T}C=0.$$

$$\Vert \Sigma {\Vert }_{{\mathcal{H}}_2^{\mathcal{T}}}:=:=\sqrt{\mathrm{tr}\left(\sum _{i=1}^3{\int }_0^{\infty }\cdots {\int }_0^{\infty }{\tilde{f}}_i(t_1,\dots ,t_i){\tilde{f}}_i(t_1,\dots ,t_i{)}^{T}d{t}_1\cdots d{t}_i\right)},$$

$${\tilde{f}}_i(t_1,\dots ,t_i)=C{\tilde{g}}_i(t_1,\dots ,t_i),$$

$${\tilde{g}}_1(t_1)$$

$${\tilde{g}}_2(t_1,t_2)$$

$${\tilde{g}}_3(t_1,t_2,t_3)$$

$$\Vert \Sigma {\Vert }_{{\mathcal{H}}_2^{\mathcal{T}}}=\sqrt{\mathrm{tr}\left(C{P}_{\mathcal{T}}{C}^T\right)}=\sqrt{\mathrm{tr}\left(B^T{Q}_{\mathcal{T}}{Q}^T\right)},$$

$$A{P}_{\mathcal{T}}+P_{\mathcal{T}}{A}^T$$

$$A^T{Q}_{\mathcal{T}}+Q_{\mathcal{T}}A$$

$$A{P}_1+P_1{A}^T$$

$$-\left(\Lambda \otimes E+I_r\otimes A\right)\mathrm{vec}\left(V_1\right)=\mathrm{vec}\left(B{\tilde{B}}^T\right),$$

$$E_{11}\dot{v}(t)$$

$$0$$

$$y(t)$$

$$0$$

$$=A_{21}{E}_{11}^{-11}\left(A_{11}v(t)+A_{12}p(t)+Hv(t)\otimes v(t)+\sum _{k=1}^m{N}_{k}v(t)u_k(t)+B_{1}u(t)\right)$$

$$(\text{using eqrefeq:differentialpart}).$$

$$p(t)=-S^{-1}{A}_{21}{E}_{11}^{-1}\left(A_{11}v(t)+Hv(t)\otimes v(t)+\sum _{k=1}^m{N}_{k}v(t)u_k(t)+B_{1}u(t)\right),$$

Full text

1 Introduction

We investigate model order reduction for quadratic-bilinear (QB) descriptor systemsxx of the form

[TABLE]

where $E_{11}$, $A_{11}$, $N_{k}\in\mathbb{R}^{n_{v}\times n_{v}}$, $H\in\mathbb{R}^{n_{v}\times n_{v}^{2}}$, $A_{12}$, $A_{21}^{T}\in\mathbb{R}^{n_{v}\times n_{p}}$, $B_{1}\in\mathbb{R}^{n_{v}\times m}$, $B_{2}\in\mathbb{R}^{n_{p}\times m}$, $C_{1}\in\mathbb{R}^{p\times n_{v}}$, $C_{2}\in\mathbb{R}^{p\times n_{p}}$; $v(t)\in\mathbb{R}^{n_{v}}$ and $p(t)\in\mathbb{R}^{n_{p}}$ are the state vectors; $u(t)\in\mathbb{R}^{m}$ and $y(t)\in\mathbb{R}^{p}$ are the control input and measured output vectors of the system, respectively, and $u_{k}(t)$ is the $k$th component of the input vector; $v_{0}\in\mathbb{R}^{n_{v}}$ is an initial value for $v(t)$. Furthermore, we assume that $E_{11}$ and $A_{21}E_{11}^{-1}A_{12}$ are invertible. Hence, the linear part of the system (1) ($H=0$, $N_{k}=0$) has an index-2 structure, e.g., see HaiW02 . The special structure for QB descriptor systems as in (1), particularly, appears in the modeling of the incompressible Navier-Strokes equations.

A high-fidelity modeling of dynamical systems is required to have a better understanding of underlying dynamical behaviors of a system. However, numerical simulations of such high-fidelity systems are expensive and often inefficient. Thus, it is not a straightforward task, or sometimes not even possible, to perform engineering design studies using these high-fidelity systems. One approach to circumvent this problem is model order reduction (MOR), aiming at constructing surrogate models (reduced-order models) which are less complicated and replicate the important dynamics of the high-fidelity system.

MOR techniques for linear systems are now very well-established and are widely applied in numerous applications, e.g., see morAnt05 ; morBauBF14 ; morBenMS05 ; morSchVR08 . Several of those techniques have been successfully extended to special classes of nonlinear ODE systems, namely bilinear and QB systems, see, e.g., morBenB12b ; morBenB12 ; morBenB15 ; morBenD11 ; morBTQBgoyal ; flagg2015multipoint ; morGu09 ; morZhaL02 . These techniques can be classified mainly into two categories, namely, trajectory-based methods, and system theoretic approaches such as balanced truncation and moment-matching. Briefly mentioning the primary ideas of these methods, trajectory-based methods rely on a set of snapshots of the state solutions for training inputs, which is then used to determine a Galerkin projection, for more details, see, e.g., morAstWWetal08 ; morChaS10 ; morHinV05 ; morKunV02 ; morKunV08 . Furthermore, the idea of balanced truncation is to find states which are hard to control as well as hard to observe, and truncating such states gives us a reduced-order system. This method for QB systems has been recently studied in morBTQBgoyal . Interpolation-based methods aim at constructing reduced-order systems which approximate the input-output behavior of the system. With this intent, such a problem for QB systems was first considered in morGu09 , where a one-sided projection method to obtain an interpolating reduced-order system is proposed. Later on, a similar problem was addressed in morBenB12a ; morBenB15 for single-input single-output (SISO) QB systems, where a two-sided projection method was proposed, ensuring a higher number of moment to be matched. However, the main challenges for this method are a good selection of interpolation points and the application to multi-input multi-output (MIMO) QB systems. To address these issues, an interpolation-based optimal model reduction problem for QB systems was recently addressed in morBenGG16 , where a reduced-order system is constructed, aiming at minimizing a system norm of the error system, e.g., a truncated $\mathcal{H}_{2}$-norm.

However, there has been little attention paid to descriptor systems (DAEs) which involve algebraic constraints as well along with a differential part, and this still requires further research. In the direction of MOR for nonlinear DAEs, interpolation-based methods for specially structured bilinear DAEs have been investigated, e.g., in morBenG16 ; goyal2016iterative . Furthermore, a moment-matching method for SISO QB systems as in (1) was studied in MPIMD15-18 . However, the first challenge in this method is the choice of interpolation points which plays a crucial role in determining the quality of a reduced-order system, and secondly, it is applicable only to SISO systems which is certainly very restrictive from a real-world application point of view.

In this work, we aim to extend an $\mathcal{H}_{2}$-optimal model reduction framework for QBODEs morBenGG16 to QBDAEs, having the structure as in (1). To that end, we first recall an interpolation-based $\mathcal{H}_{2}$-optimal model reduction technique for QBODEs and the corresponding iterative scheme to construct reduced-order systems, in the subsequent section. In Section 3, we first present the transformation of the system (1) into an equivalent ODE system by means of projectors. We further investigate how the iterative scheme can efficiently be applied to the equivalent ODE systems to obtain reduced-order systems without computing the projectors explicitly. Finally, in Section 4, we illustrate the proposed methodology using a lid driven cavity model, which is obtained by semi-discretized Navier-Stokes equations.

2 Model Reduction for Quadratic-Bilinear ODEs and Related Work

In this section, we briefly discuss an $\mathcal{H}_{2}$-optimal model reduction problem for QBODEs. We begin by introducing the problem setting for QBODEs. For this, we consider a QB system in the state-space form:

[TABLE]

where $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$ and $y(t)\in\mathbb{R}^{p}$ are the state, input, and output vectors at the time instant $t$; all other matrices are real and are of appropriate dimensions. It is assumed that the matrix $A$ is stable. Furthermore, one may also consider a general nonsingular matrix $E$ in front of $\dot{x}(t)$; however, to keep the discussion simple in this section, we consider it to be an identity matrix.

In the context of model reduction, our aim is to replace the system (2) with a simpler and reliable reduced system, having a form:

[TABLE]

where $\hat{x}(t)\in\mathbb{R}^{r}$ with $r\ll n$, all other reduced matrices are of appropriate dimensions, and the reduced output $\hat{y}(t)$ approximates the corresponding original output $y(t)$ for a wide range of admissible $L_{2}$-bounded system inputs.

We focus on constructing a reduced system (3) by means of Petrov-Galerkin projections. For this, we need to identify two appropriate subspaces $V\in\mathbb{R}^{n\times r}$ and $W\in\mathbb{R}^{n\times r}$ such that the state $x(t)$ can be approximated as $x(t)\approx V\hat{x}(t)$ where $\hat{x}(t)$ is a reduced state and satisfies a Petrov-Galerkin condition as follows:

[TABLE]

Assuming $W^{T}V$ is invertible, it yields a reduced system (3) with the realization computed as

[TABLE]

As noted in the introduction, there are several methods proposed in the literature to determine these projection matrices $V$ and $W$ or the corresponding subspaces. However, our primary objective is to extend the $\mathcal{H}_{2}$-optimal model reduction technique morBenGG16 to the QBDAE (1). Before we proceed further, we first recall important properties of the $\mathop{\mathrm{vec}}\left(\cdot\right)$ operator, stacking the columns of a matrix on top of each other:

[TABLE]

where $\otimes$ denotes the Kronecker (tensor) product. Secondly, we note down an important concept from tensor theory, that is, matricization of a tensor.

Definition 1

(e.g., koldatensor09 ) Consider a $K$-th order tensor $\mathcal{L}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{K}}$. Then, the mode-$k$ matricization of the tensor $\mathcal{L}$, denoted by $L^{(k)}$, is determined by mapping the elements $(i_{1},i_{2},\ldots,i_{K})$ of the tensor onto the matrix entries $(i_{k},j)$ as follows:

[TABLE]

Using the above concept, we can construct a 3rd-order tensor $\mathcal{H}^{n\times n\times n}$ such that its mode-1 matricization of $\mathcal{H}$ is the same as the Hessian $H$ in (2), and let $H^{(2)}$ and $H^{(3)}$ denote the mode-2 and mode-3 matricization of the tensor $\mathcal{H}$. Furthermore, without loss of generality, we can assume that $H(w_{1}\otimes w_{2})=H(w_{2}\otimes w_{1})$ for given vectors $w_{1},w_{2}\in\mathbb{R}^{n}$, see morBenB15 .

Having said that, we note down the definition of the $\mathcal{H}_{2}$-norm and its truncated version, the so-called truncated $\mathcal{H}_{2}$-norm for QB systems (2).

Definition 2

morBenGG16 Consider a QB system (2). Assume the $\mathcal{H}_{2}$-norm of the system exists, the $\mathcal{H}_{2}$-norm is defined as

[TABLE]

in which

[TABLE]

where the $g_{i}(t_{1},\ldots,t_{i})$ satisfy

[TABLE]

Remark 1

We would like to mention that the $g_{i}(t_{1},\ldots,t_{i})$ are the kernels of the Volterra series of the QB system (2) that maps the input-to-output of the system morBenGG16 . For more details on Volterra series expressions for QB systems, we refer to morBenGG16 .

Furthermore, the connection between the above $\mathcal{H}_{2}$-norm definition for QB systems and the recently proposed algebraic Gramians for QB systems morBTQBgoyal has been studied in morBenGG16 . Therein, it is shown that the $\mathcal{H}_{2}$-norm can also be computed in terms of the Gramians as follows:

[TABLE]

where $P$ and $Q$ are the controllability and observability Gramians for QB systems which satisfy the following quadratic Lyapunov equations:

[TABLE]

However, investigating an optimal model reduction problem for QB systems using the $\mathcal{H}_{2}$-norm is not a trivial task. Therefore, to ease the optimization problem, a concept of truncated $\mathcal{H}_{2}$-norms for QB systems has also been introduced in morBenGG16 , which mainly relies on the first three terms of the corresponding Volterra series. Precisely, one possible truncated $\mathcal{H}_{2}$-norm $\|\Sigma\|_{\mathcal{H}_{2}^{\mathcal{T}}}$ is defined as follows:

[TABLE]

in which

[TABLE]

where $\tilde{g}_{i}(t_{1},\ldots,t_{i})$ satisfy

[TABLE]

Similar to the $\mathcal{H}_{2}$-norm expression, a connection between a truncated $\mathcal{H}_{2}$-norm (10) and the truncated Gramians for QB systems, introduced in morBTQBgoyal , has also been established in morBenGG16 . Thus, an alternative way to compute the truncated $\mathcal{H}_{2}$-norm is as follows:

[TABLE]

where $P_{\mathcal{T}}$ are $Q_{\mathcal{T}}$ are, respectively, truncated versions of the controllability and observability Gramians, satisfying

[TABLE]

in which $P_{1}$ and $Q_{1}$ are the unique solutions of the conventional Lyapunov equations:

[TABLE]

Using the truncated $\mathcal{H}_{2}$ measure, the aim is to construct a reduced-order system such that the measure of the error system is minimized. In other words, we need to determine a reduced-order system such that $\|\Sigma-\hat{\Sigma}\|_{\mathcal{H}_{2}^{\mathcal{T}}}$ is minimized. This problem has been studied in detail in morBenGG16 , where, first, necessary conditions for optimality are derived. Based on these conditions, an iterative scheme is proposed, which on convergence, yields a reduced-order system that approximately satisfies the derived optimality conditions. The iterative scheme is referred to as truncated quadratic-bilinear rational Krylov algorithm (TQB-IRKA). A brief sketch of TQB-IRKA is given in Algorithm 1 which considers a generalized nonsingular matrix $E$ as well in the system (2) .

Remark 2

To apply Algorithm 1, we need to compute Kronecker products such as $H(V_{1}\otimes V_{1})\tilde{H}^{T}$. In a large-scale setting, a direct computation of such Kronecker products is infeasible. As a remedy, in morBenB15 ; morBenGG16 , alternative methods are proposed to perform these computations efficiently by using some tools from tensor theory or by using the special structure of the Hessian $H$, arising from semi-discretization of PDEs.

Remark 3

For example, Step 1 in Algorithm 1 which solves for $V_{1}$ is presented in a Sylvester equation form. However, using the $\mathop{\mathrm{vec}}\left(\cdot\right)$ property (4), one can write it into a linear system, that is,

[TABLE]

where $I_{r}$ is the identity matrix of size $r\times r$. Similar expressions can also be derived for the matrices $V_{2}$, $W_{1}$ and $W_{2}$ in Algorithm 2. The above $\mathop{\mathrm{vec}}\left(\cdot\right)$ form to compute the projection matrices is very useful for the latter part of the paper.

3 Transformation of a QBDAE into a QBODE and Model Reduction

Our next task is to investigate how TQB-IRKA (Algorithm 1) can be applied to QBDAEs, having the structure as in (1). For this, we first transform the system (1) into an equivalent ODE system by means of projections. Such a transformation is widely done in the literature for Navier-Stokes type equations, see, e.g. MPIMD15-18 ; morHeiSS08 . For completeness, we show the necessary steps that transform the system (1) into an equivalent ODE system. We begin by considering $B_{2}=0$ and the zero initial condition $v(0)=0$ in (1), that is:

[TABLE]

From (13b), we get $A_{21}\tfrac{d}{dt}v(t)=0$, thus, leading to

[TABLE]

As a result, we obtain an explicit expression for the pressure $p(t)$ as

[TABLE]

where $S=A_{21}E_{11}^{-1}A_{12}$. Substituting for $p(t)$ using (14) in (13a) and (13c) yields

[TABLE]

where

[TABLE]

Recall the properties of $\Pi_{l}$ from morHeiSS08 , which are as follows:

[TABLE]

This implies that $\Pi_{l}$ is an oblique projector. Moreover, for later purpose, we also define another oblique projector

[TABLE]

First, we note down a relation between $\Pi_{l}$ and $\Pi_{r}$, that is, $\Pi_{l}E_{11}=E_{11}\Pi_{r}$. Moreover, it can be verified that

[TABLE]

Using the properties of the projectors $\Pi_{l}$ and $\Pi_{r}$, we obtain

[TABLE]

with the zero initial condition, i.e., $v(0)=0$. Moreover, the dynamics of the system (16) lies in the $n_{v}{-}n_{p}$ dimensional subspace, that is nothing but the null space of $\Pi_{l}$. Next, we decompose the projectors $\Pi_{l}$ and $\Pi_{r}$ as follows:

[TABLE]

in which $\theta_{\{l,r\}}$, $\phi_{\{l,r\}}\in\mathbb{R}^{n_{v}\times n_{v}-n_{p}}$ satisfy

[TABLE]

As a result, if one defines $\tilde{v}(t):=\phi_{r}^{T}v(t)$, we consequently obtain an equivalent ODE system of the system (16) as follows:

[TABLE]

where $\tilde{v}(0)=0$. Thus, one can apply Algorithm 1 to obtain projection matrices which give a near-optimal reduced system, having neglected nonlinear terms in the output equation (18b). However, there are two major issues: one is related to computations of $\phi_{l}$, $\theta_{r}$, or $\Pi_{\{r,l\}}$ which are expensive to compute, and in case we are able to determine these projectors and their decompositions efficiently, the matrix coefficients of the system (18), e.g., $\phi_{l}^{T}E_{11}\theta_{r}$, $\phi_{l}^{T}A_{11}\theta_{r}$, might be dense, thus making model reduction techniques numerically inefficient and expensive. Therefore, in the following, we discuss how to employ Algorithm 1 to the system (18) without requiring explicit computations of the projection matrices and their decompositions.

Computational Issues

In order to determine the projection matrices to compute a reduced system, we consider the following associated QB system, which is the system (18) but some neglected nonlinear terms in the output equation:

[TABLE]

where

[TABLE]

Here, we need to neglect the nonlinear terms appearing in the output equation due to the transformation to an ODE system as Algorithm 1 is developed so far only for linear output equations. The full incorporation of the nonlinearities in the output equation requires further work. But note that the nonlinear terms are again included when the projections are applied to obtain the reduced-order model, see Remark 4.

As a first step, we aim at determining the projection matrices $V$ and $W$ such that the original matrices like $E_{11}$, $A_{11}$ can be used to compute the reduced matrices instead of using, e.g., $\bar{E}$, $\bar{A}$, i.e., a reduced matrix $\hat{E}$ can be computed as $W^{T}E_{11}V$, and so on.

For this purpose, let $\bar{V}$ and $\bar{W}$ be the solutions of the Sylvester equations in Algorithm 1 (steps 6–1) using the matrices $\bar{E}$, $\bar{A}$, etc. Furthermore, we define the matrices $V$ and $W$, satisfying:

[TABLE]

Then, it can easily be verified that the reduced matrices computed using the quantities denoted by bar, e.g., $\bar{V}$, $\bar{W}$, $\bar{E}$ are the same as the reduced matrices computed using $V$ and $W$ and the original matrices. In other words, $\bar{W}^{T}\bar{E}\bar{V}=W^{T}E_{11}V$ and so on. Using this formulation, we sketch an algorithm for model reduction of (19) which gives near-optimal reduced systems on convergence by projecting the original system matrices, see Algorithm 2.

However, to compute the projection matrices at each iteration in Algorithm 2, we still require the basis matrices $\phi_{l}$ and $\theta_{r}$. Therefore, our next goal is to determine these projection matrices without involving $\phi_{l}$ and $\theta_{r}$, or more precisely, we aim at computing the matrices $V$ and $W$ using only the original matrices such as $E_{11}$, $A_{11}$, $A_{12}$.

Luckily, a similar problem has been studied in goyal2016iterative , where for a symmetric case, i.e., $\Pi_{l}=\Pi_{r}^{T}$, it is shown how to compute $(I_{r}\otimes\theta_{r})L(I_{r}\otimes\phi_{l}^{T})f$ for given arbitrary vector $f$ and matrix $L$, without explicitly forming $\theta_{r}$ and $\phi_{l}$. However, for the unsymmetric case, i.e., $\Pi_{l}\neq\Pi_{r}^{T}$, the following result can be developed in a similar fashion as in goyal2016iterative , and, therefore, a detailed proof is omitted.

Lemma 1

Consider $\phi_{l}$ and $\theta_{r}$ as in (17), and assume $\mathcal{X}=(I_{r}\otimes\phi_{l}^{T})\mathcal{T}(I_{r}\otimes\theta_{r})$ is invertible for a given $\mathcal{T}$. Furthermore, let $\mathcal{G}$ and $\mathcal{G}_{T}$ be $(I_{r}\otimes\theta_{r})\mathcal{X}^{-1}(I_{r}\otimes\phi_{l}^{T})$ and $(I_{r}\otimes\theta_{l})\mathcal{X}^{-T}(I_{r}\otimes\phi_{r}^{T})$, respectively. Then, the vector

[TABLE]

solves

[TABLE]

Similarly, the vector

[TABLE]

solves

[TABLE]

By making use of the result of Lemma 1, it can be shown that the steps 7-10 in Algorithm 2 can be performed without an explicit requirement of the basis matrices $\phi_{l}$ and $\theta_{r}$. We rather need to solve appropriate saddle point problems to compute $V_{\{1,2\}}$ and $W_{\{1,2\}}$. All these analyses lead to Algorithm 3 that does not require any explicit computation of the basis matrices $\phi_{l}$ and $\theta_{r}$ even for computing projection matrices.

Remark 4

Recall that Algorithm 3 gives on convergence a near-optimal reduced system for the system (18), having neglected the nonlinear terms in the output equation in the model reduction process. Nonetheless, we reduce these nonlinear terms in the output equation using the projection matrix $V$, obtained on convergence, i.e., $\hat{\mathcal{C}}_{H}=\mathcal{C}_{H}(V\otimes V)$ and $\hat{\mathcal{C}}_{N_{k}}=\mathcal{C}_{N_{k}}V$. Furthermore, if the output of the system (13) is given only by linear combinations of the velocity $v(t)$, then all nonlinear terms in the output equation of the system (18) are zero.

Remark 5

Throughout the above discussion, we have assumed that $B_{2}=0$ in (1). However, as discussed, e.g., in morHeiSS08 , the general case $B_{2}\neq 0$ can be converted into the case $B_{2}=0$ by an appropriate change of variables for $v(t)$. For this, one needs to decompose $v(t)$ as follows:

[TABLE]

where $v_{u}(t)=\Omega u(t)$ in which $\Omega:=-E_{11}^{-1}A_{12}(A_{21}E_{11}A_{12})B_{2}$. If one substitutes for $v(t)$ from (22) into (1), then it can be easily seen that $A_{21}v_{m}(t)=0$. Furthermore, we assume an initial condition to be consistent. Performing similar algebraic calculations as done for the problem $B_{2}=0$ leads to the following system:

[TABLE]

where

[TABLE]

in which $S:=A_{21}E_{11}^{-1}A_{12}$. Thus, we obtain a system equivalent to (1) which has a similar form as in (16). However, the system (23) contains some extra terms which are functions of the input $u(t)$, e.g., $u_{p}\cdot u_{q}$, $\{p,q\}\in\{1,\ldots,m\}$ and the derivative of the input $u(t)$. Although they are functions of the input in the forward simulation, we consider them formally as different inputs as far as the model reduction problem is concerned. Hence, Algorithm 3 can readily be applied to the system (23) to determine reduced-order systems.

4 Numerical Experiment

In this section, we test the efficiency of Algorithm 3 using a numerical example of a lid-driven cavity, obtained using semi-discretization of the Navier–Stokes equation. We initialize Algorithm 3 randomly and iterate it until the relative change in the eigenvalues of the reduced matrix $\hat{A}$ is less than $10^{-5}$. All the simulations were done on a board with 4 Intel® Xeon® E7-8837 CPUs with a 2.67- GHz clock speed using MATLAB® 8.0.0.783 (R2012b).

4.1 Lid Driven Cavity

Here, we consider a lid-driven cavity as shown in Figure 1, which is modeled using the incompressible Navier-Stokes equations in the velocity $\breve{v}$ and the pressure $\breve{p}$ on the unit square $\Omega=(0,1)^{2}$. The governing equations and boundary conditions are:

[TABLE]

and boundary and initial conditions are as follows:

[TABLE]

where $\mathop{\mathrm{Re}}$ is the Reynolds number, $g$ represents the Dirichlet boundary conditions, and we set $\breve{v}=\begin{bmatrix}1\\ 0\end{bmatrix}$ for the upper boundary and zero at the other boundary, and $\breve{v}_{0}$ is an initial condition. Furthermore, the system is subject to a control input $u(t)$ in the domain $\Omega_{c}$. As a result, we ensure that the $x$ and $y$ components of the velocity in the domain $\Omega_{c}$ are $u(t)$.

A similar problem has been considered in the model reduction framework in MPIMD15-18 for a SISO system and with a different control setting. Having applied a finite element scheme using the Taylor-Hood scheme on a uniform mesh, we obtain a discretized system in velocity $v$ and pressure $p$ as follows:

[TABLE]

where $E_{11},A_{11}\in\mathbb{R}^{n_{v}\times n_{v}}$, $A_{12}\in\mathbb{R}^{n_{v}\times n_{p}}$, $H\in\mathbb{R}^{n_{v}\times n_{v}^{2}}$, $f\in\mathbb{R}^{n_{v}}$ and $v_{s}$ is the steady-state solution of the system when the control input $u$ is set to zero. Furthermore, we observe the velocity in the domain $\Omega_{o}$, leading to an output equation as follow:

[TABLE]

Since the initial condition is considered to be the steady-state solution of the system, we perform a change of variables as $v_{\delta}=v+v_{s}$ and $p_{\delta}=p+p_{s}$ to ensure the zero initial condition of the transformed system. This results in the following system:

[TABLE]

where $X:=H(v_{s}\otimes I+I\otimes v_{s})$. We choose the degrees of freedom for the velocity and the pressure $3042$ and $440$, respectively, i.e., $n_{v}=3042$ and $n_{p}=440$. We set the Reynolds number $\mathop{\mathrm{Re}}=100$. Furthermore, we select four nodes in the domain $\Omega_{o}$ on which we measure the $x$ and $y$ components of the velocity; thus, we have $8$ outputs. Next, we employ Algorithm 3 to obtain a reduced-order system of order $r=140$. To check the accuracy of the reduced system, we perform time-domain simulations for the original and reduced systems for an arbitrary control input $u(t)=2t^{2}\exp(-t/2)\sin(2\pi t/5)$ and observe the outputs which are plotted in Figure 2. This shows that the reduced system captures the dynamics of the system faithfully. Furthermore, we also plot the velocity at the full grid in Figure 3, although the considered model reduction aims at capturing only the input-to-output behavior of the system. The figure indicates that the reduced-order system can replicate the visual dynamics of the original system at the full grid as well.

5 Conclusions

In this work, we have studied a model reduction problem for a particular class of quadratic-bilinear descriptor systems, especially arising from semi-discretization of the Navier-Stokes equations. In particular, we have investigated how one can employ the iterative model reduction scheme proposed in morBenGG16 for quadratic-bilinear ODEs to the latter class of descriptor systems. To that end, we have first transformed the quadratic-bilinear descriptor system into an equivalent ODE system by means of projectors. Furthermore, for practical computations, we have proposed an efficient iterative scheme for the considered quadratic-bilinear descriptor systems to avoid the explicit computation of the projectors that are used to transform the system into an ODE system. Using as numerical example a semi-discretized Navier-Stokes equations, we have shown the efficiency of the proposed method.

As a further research topic, an extension of the balanced truncation method for quadratic-bilinear ODE systems morBTQBgoyal to descriptor systems will be very useful. Furthermore, it will be a significant contribution to extended these methods to other classes of descriptor systems, particularly, appearing in mechanical systems.