A block tangential Lanczos method for model reduction of large-scale first and second order dynamical systems

TL;DR

This paper introduces an adaptive block tangential Lanczos method for reducing large-scale first and second order dynamical systems, improving computational efficiency while maintaining accurate system behavior approximation.

Contribution

The paper proposes a novel adaptive block tangential Lanczos algorithm for model reduction of large-scale MIMO systems, with new algebraic properties and demonstrated effectiveness.

Findings

The ABTL algorithm effectively reduces system complexity.

Numerical experiments show high approximation accuracy.

The method is applicable to large-scale systems with multiple inputs and outputs.

Abstract

In this paper, we present a new approach for model reduction of large scale first and second order dynamical systems with multiple inputs and multiple outputs (MIMO). This approach is based on the projection of the initial problem onto tangential Krylov subspaces to produce a simpler reduced-order model that approximates well the behavior of the original model. We present an algorithm named: Adaptive Block Tangential Lanczos-type (ABTL) algorithm. We give some algebraic properties and present some numerical experiences to show the effectiveness of the proposed algorithms.

Tables

Table 1. Table 1: Matrix Tests

Model										$n$		$p$		$s$
FDM10000										n = 40 000		p = 6		s = 3
FDM90000										n = 90 000		p = 6		s = 3
1DBeam-LF1000										n = 19 998		p = 4		s = 2
1DBeam-LF5000										n = 19 994		p = 4		s = 2
RAIL79841										n = 17 361		p = 12		s = 2

Table 2. Table 2: The calculation time and the ℋ ∞ subscript ℋ {\mathcal{H_{\infty}}} error-norm

Model		ABTL	IRKA	TRKSM
		Time  Err-${\mathscr{H}}_{\infty }$	Time  Err-${\mathscr{H}}_{\infty }$	Time  Err-${\mathscr{H}}_{\infty }$
FDM40.000	m=20	$13.29s5.39\times {10}^{-4}$	$126.28s1.06\times {10}^{-4}$	$34.89s7.94\times {10}^{-4}$
	m=30	$20.79s7.87\times {10}^{-5}$	$188.91s2.24\times {10}^{-5}$	$36.82s6.42\times {10}^{-5}$
	m=40	$27.21s1.48\times {10}^{-5}$	$269.36s1.56\times {10}^{-5}$	$38.12s1.93\times {10}^{-5}$
FDM90.000	m=20	$43.291.27\times {10}^{-2}$	$>2000s$  – –	$126.97s2.61\times {10}^{-1}$
	m=30	$64.111.36\times {10}^{-3}$	$>2000s$  – –	$179.55s1.43\times {10}^{-1}$
	m=40	$87.196.25\times {10}^{-4}$	$>2000s$  – –	$186.38s4.61\times {10}^{-2}$

The file name		Degrees of freedom	Num. nodes	Dimension $n$
1DBeam-LF100		flexural (${\Theta }_z$ and $y$)	10000	n = 19998
1DBeam-LF5000		( ${\Theta }_z$ and $y$ ), ( ${\Theta }_x$ and $x$ )	50000	n = 19994

Equations

$$\Sigma:=\left\{\begin{array}[]{lllll}\dot{x}(t)&=&Ax(t)+Bu(t)\\ y(t)&=&Cx(t),\\ \end{array}\right.$$

$$H(\omega ):=C(\omega I_n-A{)}^{-1}B\in {\mathbb{R}}^{p\times p}.$$

$$\Sigma_{m}:\left\{\begin{array}[]{lllll}\dot{x}_{m}(t)&=&A_{m}x_{m}(t)+B_{m}u(t)\\ y_{m}(t)&=&C_{m}x_{m}(t),\\ \end{array}\right.$$

$$H_m(\omega ):=C_m(\omega I_m-A_m{)}^{-1}{B}_m\in {\mathbb{R}}^{p\times p},$$

$${\tilde{\mathcal{K}}}_m(A,B)=Range\{({\sigma }_1{I}_n-A{)}^{-1}B{R}_1,...,({\sigma }_m{I}_n-A{)}^{-1}B{R}_m\},$$

$${\tilde{\mathcal{K}}}_m(A^T,C^T)=Range\{({\mu }_1{I}_n-A{)}^{-T}{C}^T{L}_1,...,({\mu }_m{I}_n-A{)}^{-T}{C}^T{L}_m\},$$

$$H(\omega )=\sum _{i=0}^{\infty }{\eta }_i^{(\sigma )}\frac{(\omega -\sigma {)}^i}{i!},$$

$${\eta }_i^{(\sigma )}=C(\sigma I_n-A{)}^{-(i+1)}B=(-1{)}^i\frac{d^i}{d{\omega }^i}H(\omega ){|}_{\omega =\sigma },i=0,1,...$$

$${\eta }_i=C{A}^{i}B.$$

$$\left\{\begin{array}[]{lllll}H_{m}(\sigma_{i})R_{i}&=&H(\sigma_{i})R_{i}\\ &&&for\ i=1,...,m.\\ L_{i}^{T}H_{m}(\mu_{i})&=&L_{i}^{T}H_{(}\mu_{i})\\ \end{array}\right.$$

$$(\omega I_n-A)X=B,\mathrm{and}(\omega I_n-A{)}^{T}Y=C^T.$$

$$X_m^i\in Range\{V_1,...,V_m\}and{Y}_m^i\in Range\{W_1,...,W_m\}$$

$$R_B^i(\omega )\perp Range\left\{W_1,...,W_m\right\},i=1,...,p$$

$$R_C^i(\omega )\perp Range\left\{V_1,...,V_m\right\},i=1,...,p$$

$$X_m={\mathbb{V}}_m(\omega I_{ms}-A_m{)}^{-1}{\mathbb{W}}_m^{T}B,$$

$$Y_m={\mathbb{W}}_m(\omega I_{ms}-A_m{)}^{-T}{\mathbb{V}}_m^T{C}^T,$$

$$H_m(\omega )=C_m(\omega I_{ms}-A_m{)}^{-1}{B}_m,$$

$$R_B(\omega )=B-(\omega I_n-A){\mathbb{V}}_m(\omega I_{ms}-A_m{)}^{-1}{\mathbb{W}}_m^{T}B.$$

$$R_C(\omega )=C^T-(\omega I_n-A{)}^T{\mathbb{W}}_m(\omega I_{ms}-A_m{)}^{-T}{\mathbb{V}}_m^T{C}^T.$$

$${\mathcal{K}}_m(A,B)=Range\{({\sigma }_1{I}_n-A{)}^{-1}B{R}_1,...,({\sigma }_m{I}_n-A{)}^{-1}B{R}_m\},$$

$${\tilde{\mathcal{K}}}_m(A^T,C^T)=Range\{({\mu }_1{I}_n-A{)}^{-T}{C}^T{L}_1,...,({\mu }_m{I}_n-A{)}^{-T}{C}^T{L}_m\},$$

$$\left\{\begin{array}[]{llll}W_{i}^{T}V_{j}=I,&i=j,\\ W_{i}^{T}V_{j}=0,&i\neq j.\\ \end{array}\right.$$

$$\widetilde{\mathbb{H}}^{(j)}=\left[\begin{array}[]{cccc}H_{1,j}\\ \vdots\\ H_{j,j}\\ H_{j+1,j}\\ \textbf{0}\end{array}\right]\,{\rm and}\ \ \widetilde{\mathbb{F}}^{(j)}=\left[\begin{array}[]{cccc}F_{1,j}\\ \vdots\\ F_{j,j}\\ F_{j+1,j}\\ \textbf{0}\end{array}\right],\quad{\rm for}\ j=0,...,m.$$

$${\tilde{D}}_{m+1}^{(1)}=D_{m+1}^{(1)}\otimes I_s,{\tilde{D}}_{m+1}^{(2)}=D_{m+1}^{(2)}\otimes I_s$$

$$A{\mathbb{V}}_{m+1}=\left[{\mathbb{V}}_{m+1}{\mathbb{G}}_{m+1}{\tilde{D}}_{m+1}^{(1)}-B{\mathbb{R}}_{m+1}\right]{\mathbb{G}}_{m+1}^{-1},$$

$$A^T{\mathbb{W}}_{m+1}=\left[{\mathbb{W}}_{m+1}{\mathbb{Q}}_{m+1}{\tilde{D}}_{m+1}^{(2)}-C^T{\mathbb{L}}_{m+1}\right]{\mathbb{Q}}_{m+1}^{-1}.$$

$${\mathbb{T}}_{m+1}=\left[({\sigma }_{1}I-A{)}^{-1}B{R}_1,...,({\sigma }_{m+1}I-A{)}^{-1}B{R}_{m+1}\right]\mathrm{and}$$

$${\mathbb{Y}}_{m+1}=\left[({\mu }_{1}I-A{)}^{-T}{C}^T{L}_1,...,({\mu }_{m+1}I-A{)}^{-T}{C}^T{L}_{m+1}\right],$$

$${\mathbb{T}}_{m+1}={\mathbb{V}}_{m+1}{\mathbb{G}}_{m+1}\mathrm{and}{\mathbb{Y}}_{m+1}={\mathbb{W}}_{m+1}{\mathbb{Q}}_{m+1},$$

$$V_{j+1}{H}_{j+1,j}=({\sigma }_{j+1}{I}_n-A{)}^{-1}B{R}_{j+1}-\sum _{i=1}^j{V}_i{H}_{i,j}j=1,...,m,$$

Full text

1 Introduction

Consider a linear time-invariant multi-input, multi-output linear time independent (LTI) dynamical system described by the state-space equations

[TABLE]

where $x(t)\in\mathbb{R}^{n}$ denotes the state vector, $u(t)$ and $y(t)$ are the input and the output signal vectors, respectively. The matrix $A\in\mathbb{R}^{n\times n}$ is assumed to be large and sparse and $B,\ C^{T}\in\mathbb{R}^{n\times p}$. The transfer function associated to the system in (1) is given as

[TABLE]

The goal of our model reduction approach consists in defining two orthogonal matrices $V_{m}$ and $W_{m}\in\mathbb{R}^{n\times m}$ (with $m\ll n$) to produce a much smaller order system $\Sigma_{m}$ with the state-space form

[TABLE]

and its transfer function is defined by

[TABLE]

where $A_{m}=W_{m}^{T}AV_{m}\in\mathbb{R}^{m\times m}$, $B_{m}=W_{m}^{T}B\in\mathbb{R}^{m\times p}$ and $C_{m}=CV_{m}\in\mathbb{R}^{p\times m}$, such that the reduced system $\Sigma_{m}$ will have an output $y_{m}(t)$ as close as possible to the one of the original system to any given input $u(t)$, which means that for some chosen norm, $\|y-y_{m}\|$ is small.

Various model reduction techniques, such as Padé approximation [17, 27], balanced truncation [28], optimal Hankel norm [16] and Krylov subspace methods, [8, 9, 14, 21] have been used for large multi-input multi-output (MIMO) dynamical systems, see [25, 16, 20]. Balanced Truncation Model Reduction (BTMR) method is a very popular method; [1, 15]; the method preserves the stability and provides a bound for the approximation error. In the case of small to medium systems, (BTMR) can be implemented efficiently. However, for large-scale settings, the method is quite expensive to implement, because it requires the computation of two Lyapunov equations, and results in a computational complexity of $\mathcal{O}(n^{3})$ and a storage requirement of $\mathcal{O}(n^{2})$, see [1, 5, 18]. In this paper, we project the system (1) onto the following block tangential Krylov subspaces defined as,

[TABLE]

[TABLE]

in order to obtain a small scale dynamical system. The $\{\sigma_{i}\}_{i=1}^{m}$ and $\{\mu_{i}\}_{i=1}^{m}$ are the right and left interpolation points, the $\{R_{i}\}_{i=1}^{m}$ and $\{L_{i}\}_{i=1}^{m}$ are the right and left blocks tangent directions with $R_{i},\ L_{i}\in\mathbb{R}^{p\times s}$ with $s\leq p$. Later, we will show how to choose these tangent interpolation points and directions.

The paper is organized as follows: In section 2 we give some definition used later and we introduce the tangential interpolation. In Section 3, we present the tangential block Lanczos-type method and the corresponding algorithm. Section 4 is devoted to the selection of the interpolation points and the tangential directions that are used in the construction of block tangential Krylov subspaces, and we present briefly the adaptive tangential block Lanczos-type algorithm. Section 5 we treat the model reduction of second-order systems. The last section is devoted to numerical tests and comparisons with some well known model order reduction methods.

2 Moments and interpolation

We first give the following definition.

Definition 1**.**

Given the system $\Sigma$, its associated transfer function $H(s)=C(\omega I_{n}-A)^{-1}B$ can be decomposed through a Laurent series expansion around a given $\sigma\in\mathbb{R}$ (shift point), as follows

[TABLE]

where $\eta_{i}^{(\sigma)}\in\mathbb{R}^{p\times p}$ is called the $i$-th moments at $\sigma$ associated to the system and defined as follows

[TABLE]

In the case where $\sigma=\infty$ the moments are called Markov parameters and are given by

[TABLE]

Problem: Given a full-order model (1) and assume that the following parameters are given:

• •

Left interpolation points $\{\mu_{i}\}^{m}_{i=1}\subset\mathbb{C}$ & block tangent directions $\{L_{i}\}^{m}_{i=1}\subset\mathbb{C}^{p\times s}$.

• •

Right interpolation points $\{\sigma_{i}\}^{m}_{i=1}\subset\mathbb{C}$ & block tangent directions $\{R_{i}\}^{m}_{i=1}\subset\mathbb{C}^{p\times s}$.

The main problem is to find a reduced-order model (3) such that the associated transfer function, $H_{m}$ in (4) is a tangential interpolant to $H$ in (2), i.e.

[TABLE]

The interpolation points and tangent directions are selected to realize the model reduction goals described later.

3 The block tangential Lanczos-type method

Let the original transfer function $H(\omega)=C(\omega I-A)^{-1}B$ be expressed as $H(\omega)=CX=Y^{T}B$ where $X$ and $Y$ are such that,

[TABLE]

Given a system of matrices $\{V_{1},\ldots,V_{m}\}$ and $\{W_{1},\ldots,W_{m}\}$ where $V_{i},\ W_{i}\in\mathbb{R}^{n\times s}$, the approximate solution $X_{m}$ and $Y_{m}$ of $X$ and $Y$ are computed such that

[TABLE]

and

[TABLE]

[TABLE]

where $X_{m}^{i}$, $Y_{m}^{i}$, $R_{B}^{i}$ and $R_{C}^{i}$ are the $i$-th columns of $X_{m}$, $Y_{m}$, $R_{B}=B-(\omega I_{n}-A)X_{m}$ and $R_{C}=C^{T}-(\omega I_{n}-A)^{T}Y_{m}$, respectively. If we set $\mathbb{V}_{m}=[V_{1},\ldots,V_{m}]$ and $\mathbb{W}_{m}=[W_{1},\ldots,W_{m}]$, then from (9), (10) and (11) , we obtain

[TABLE]

[TABLE]

which gives the following approximate transfer function

[TABLE]

where $A_{m}=\mathbb{W}_{m}^{T}A\mathbb{V}_{m}$, $B_{m}=\mathbb{W}_{m}^{T}B$ and $C_{m}=C\mathbb{V}_{m}$. The matrices $\mathbb{V}_{m}=[V_{1},...,V_{m}]$ and $\mathbb{W}_{m}=[W_{1},...,W_{m}]$ are bi-orthonormal, where the $V_{i},\ W_{i}\in\mathbb{R}^{n\times s}$. Notice that the residuals can be expressed as

[TABLE]

[TABLE]

3.1 Tangential block Lanczos-type algorithm

This algorithm consists in constructing two bi-orthonormal bases, spanned by the columns of $\{V_{1},V_{2},\ldots,V_{m}\}$ and $\{W_{1},W_{2},\ldots,W_{m}\}$, of the following block tangential Krylov subspaces

[TABLE]

and

[TABLE]

where $\{\sigma_{i}\}_{i=1}^{m}$ and $\{\mu_{i}\}_{i=1}^{m}$ are the right and left interpolation points respectively and $\{R_{i}\}_{i=1}^{m}$, $\{L_{i}\}_{i=1}^{m}$ are the right and left tangent directions with $R_{i},\ L_{i}\in\mathbb{R}^{p\times s}$. We present next the Block Tangential Lanczos (BTL) algorithm that allows us to construct such bases. It is summarized in the following steps.

Let $\mathbb{V}_{m}=[V_{1},V_{2},\ldots,V_{m}]$ and $\mathbb{W}_{m}=[W_{1},W_{2},\ldots,W_{m}]$. Then, we should have the bi-orthogonality conditions for $i,j=1,\ldots,m$:

[TABLE]

Here we suppose that we already have the set of interpolation points $\sigma=\{\sigma_{i}\}_{i=1}^{m+1}$, $\mu=\{\mu_{i}\}_{i=1}^{m+1}$ and the tangential matrix directions $R=\{R_{i}\}_{i=1}^{m+1}$ and $L=\{L_{i}\}_{i=1}^{m+1}$. The upper block upper Hessenberg matrices $\widetilde{\mathbb{H}}_{m}=\left[\widetilde{\mathbb{H}}^{(1)},...,\widetilde{\mathbb{H}}^{(m)}\right]$ and $\widetilde{\mathbb{F}}_{m}=\left[\widetilde{\mathbb{F}}^{(1)},...,\widetilde{\mathbb{F}}^{(m)}\right]\in\mathbb{R}^{(m+1)s\times ms}$ are obtained from the BTL algorithm, with

[TABLE]

The matrices $H_{i,j}$ and $F_{i,j}$ constructed in Step 3 of Algorithm 1 are of size $s\times s$ and 0 is the zero matrix of size $(m-j)s\times s$. We define also the following matrices,

[TABLE]

where $D_{m+1}^{(1)}=Diag\{\sigma_{1},...,\sigma_{m+1}\}$ and $D_{m+1}^{(2)}=Diag\{\mu_{1},...,\mu_{m+1}\}$. With all those notations, we have the following theorem.

Theorem 2**.**

Let $\mathbb{V}_{m+1}$ and $\mathbb{W}_{m+1}$ be the bi-orthonormal matrices of $\mathbb{R}^{n\times(m+1)s}$ constructed by Algorithm 1. Then we have the following relations

[TABLE]

and

[TABLE]

Let $\mathbb{T}_{m+1}$ and $\mathbb{Y}_{m+1}$ be the matrices,

[TABLE]

[TABLE]

then we have

[TABLE]

where $\mathbb{R}_{m+1}=\left[R_{1},...,R_{m+1}\right]$ and $\mathbb{L}_{m+1}=\left[L_{1},...,L_{m+1}\right]$. The matrices $\mathbb{G}_{m+1}=\left[\widetilde{\mathbb{H}}^{(0)}\ \ \widetilde{\mathbb{H}}_{m}\right]$ and $\mathbb{Q}_{m+1}=\left[\widetilde{\mathbb{F}}^{(0)}\ \ \widetilde{\mathbb{F}}_{m}\right]$ are block upper triangular matrices of sizes $(m+1)s\times(m+1)s$ and are assumed to be non-singular.

Proof.

From Algorithm 1, we have

[TABLE]

multiplying (20) on the left by $(\sigma_{j+1}I_{n}-A)$ and re-arranging terms, we get

[TABLE]

which gives

[TABLE]

that written as

[TABLE]

where 0 is the zero matrix of size $(m-j)s\times s$. Then for $j=1,\ldots,m$, we have

[TABLE]

now, since $V_{1}H_{1,0}=(\sigma_{1}I_{n}-A)^{-1}BR_{1}$, we can deduce from (22), the following expression

[TABLE]

which ends the proof of (17). The same proof can be done for the relation (18).

For the proof of (19), we first use (20) to obtain

[TABLE]

which gives

[TABLE]

it follows that

[TABLE]

which ends the proof of the first relation of (19). In the same manner, we can prove the second relation. ∎

Theorem 3**.**

Let $\sigma,\ \mu\in\mathbb{C}$ be such that $(\omega I-A)$ is invertible for $\omega=\sigma,\ \mu$. Let $\mathbb{V}_{m}=\left[V_{1},...,V_{m}\right]$ and $\mathbb{W}_{m}=\left[W_{1},...,W_{m}\right]$ have full-rank, where the $V_{i},\ W_{i}\in\mathbb{R}^{n\times s}$. Let $R=[r_{1},...,r_{s}],\ L=[l_{1},...,l_{s}]\in\mathbb{R}^{p\times s}$ be a chosen tangential matrix directions. Then,

If $(\sigma I-A)^{-1}Br_{i}\in Range\{V_{1},...,V_{m}\}$ for $i=1,...,s$, then

[TABLE] 2. 2.

If $(\mu I-A)^{-T}C^{T}l_{i}\in Range\{W_{1},...,W_{m}\}$ for $i=1,...,s$, then

[TABLE] 3. 3.

If both (1) and (2) hold and in addition we have $\sigma=\mu$, then,

[TABLE]

Proof.

1. We follow the same techniques as those given in [2] for the non-block case. Define

[TABLE]

and

[TABLE]

It is easy to verify that $\mathcal{P}_{m}(\omega)$ and $\mathcal{Q}_{m}(\omega)$ are projectors. Moreover, for all $\omega$ in a neighborhood of $\sigma$ we have

[TABLE]

and

[TABLE]

Observe that

[TABLE]

Evaluating the expression (23) at $\omega=\sigma$ and multiplying by $r_{i}$ from the right, yields the first assertion, and evaluating the same expression at $\omega=\mu$ and multiplying by $l_{i}^{T}$ from the left, yields the second assertion

2. Now if both (1) and (2) hold and $\sigma=\mu$, notice that

[TABLE]

and

[TABLE]

Therefore, evaluating (23) at $s=\sigma+\varepsilon$, multiplying by $l_{j}^{T}$ and $r_{i}$, from the left and the right respectively, for $i,j=1,...,s$, we get

[TABLE]

Now notice that since $l_{j}^{T}H(\sigma)r_{i}=l_{j}^{T}H_{m}(\sigma)r_{i}$, we have

[TABLE]

which proves the third assertion. ∎

In the following theorem, we give the exact expression of the residual norms in a simplified and economical computational form.

4 An adaptive choice of the interpolation points and tangent directions

In the section, we will see how to chose the interpolation points $\{\sigma_{i}\}_{i=1}^{m}$, $\{\mu_{i}\}_{i=1}^{m}$ and tangential directions $\{R_{i}\}_{i=1}^{m}$, $\{L_{i}\}_{i=1}^{m}$, where $R_{i},\ L_{i}\in\mathbb{R}^{p\times s}$. In this paper we adopted the adaptive approach, inspired by the work in [10]. For this approach, we extend our subspaces $\mathcal{K}_{m}(A,B)$ and $\widetilde{\mathcal{K}}_{m}(A^{T},C^{T})$ by adding new blocks $\widetilde{V}_{m+1}$ and $\widetilde{W}_{m+1}$ defined as follows

[TABLE]

where the new interpolation point $\sigma_{m+1}$, $\mu_{m+1}$ and the new tangent direction $R_{m+1}$, $L_{m+1}$ are computed as follows

[TABLE]

[TABLE]

Here $S_{m}\subset\mathbb{C}^{+}$ is defined as the convex hull of $\{-\lambda_{1},...,-\lambda_{m}\}$ where $\{\lambda_{i}\}_{i=1}^{m}$ are the eigenvalues of the matrix $A_{m}$. For solving the problem (25), we proceed as follows. First we compute the next interpolation point, by computing the norm of $R_{B}(\omega)$ for each $\omega$ in $S_{m}$, i.e we solve the following problem,

[TABLE]

Then the tangent direction $R_{m+1}$ is computed by evaluating (25) at $\omega=\sigma_{m+1}$,

[TABLE]

We can easily prove that the tangent matrix direction $R_{m+1}$ is given as

[TABLE]

where the $r_{i}^{(m+1)}$’s are the right singular vectors corresponding to the $s$ largest singular values of the matrix $R_{B}(\sigma_{m+1})$. This approach of maximizing the residual norm, works efficiently for small to medium matrices, but cannot be used for large scale systems. To overcome this problem, we give the following proposition.

Proposition 4**.**

Let $R_{B}(\omega)=B-(\omega I_{n}-A)\mathbb{V}_{m}U_{m}^{B}(\omega)$ and $R_{C}(\omega)=C^{T}-(\omega I_{n}-A)^{T}\mathbb{W}_{m}U_{m}^{C}(\omega)$ be the residuals given in (12) and (23), where $U_{m}^{B}(\omega)=(\omega I-A_{m})^{-1}\mathbb{W}_{m}^{T}B$ and $U_{m}^{C}(\omega)=(\omega I-A_{m})^{-T}\mathbb{V}_{m}^{T}C^{T}$. Then we have the following new expressions

[TABLE]

and

[TABLE]

Proof.

The residual $R_{B}(\omega)$ can be written as

[TABLE]

using Equation (17), we get

[TABLE]

which gives,

[TABLE]

which proves (29). In the same way we can prove (30). ∎

The expression of $R_{B}(\omega)$ given in (29) allows us to significantly reduce the computational cost while seeking the next pole and direction. In fact, applying the skinny QR decomposition

[TABLE]

we get the simplified residual norm

[TABLE]

This means that, solving the problem (25) requires only the computation of matrices of size $ms\times ms$ for each value of $\omega$.

The next algorithm, summarizes all the steps of the adaptive choice of tangent interpolation points and tangent directions.

5 Model Reduction of Second-Order Systems

Linear PDEs modeling structures in many areas of engineering (plates, shells, beams …) are often second order in time see for example [25, 26, 29]. The spatial semi-discretization of its models by a method of finite elements leads to systems that write in the form:

[TABLE]

where $M\in\mathbb{R}^{n\times n}$ is the mass matrix, $D\in\mathbb{R}^{n\times n}$ is the damping matrix and $K\in\mathbb{R}^{n\times n}$ the stiffness matrix. When the source term $Bu(t)$ is null, the system is said to be free, otherwise, it is said forced. If $D=0$, the system is said to be undamped. We assume that the mass matrix $M$ is invertible, then the system (32) can be written as

[TABLE]

where $D_{M}=M^{-1}D$, $K_{M}=M^{-1}K$ and $B_{M}=M^{-1}B$, for simplicity we still denote $K,\ D,$ and $B$ instead of $D_{M}$, $K_{M}$ and $B_{M}$. The transfer function associated with the system (33) is given by using the Laplace transform as:

[TABLE]

Usually, it’s difficult to have the efficient solution of various control or simulation tasks because the original system is too large to allow it. In order to solve this problem, methods that produce a reduced system of size $m\ll n$ that preserves the essential properties of the full order model have been developed. The reduced model have the following form:

[TABLE]

where $D_{m},\ K_{m}\in\mathbb{R}^{m\times m}$, $B_{m},\ C_{m}^{T}\in\mathbb{R}^{m\times p}$ and $q_{m}(t)\in\mathbb{R}^{m}$. The transfer function associated to the system (35) is given by:

[TABLE]

Second-order systems (33) can be written as a first order linear systems. In fact,

[TABLE]

which is equivalent to

[TABLE]

with x(t)=\left[\begin{array}[]{ccc}q(t)\\ \dot{q}(t)\\ \end{array}\right], \mathcal{A}=\left[\begin{array}[]{ccc}0&I_{n}\\ -K&-D\\ \end{array}\right], \mathcal{B}=\left[\begin{array}[]{ccc}0\\ B\\ \end{array}\right] and \mathcal{C}=\left[\begin{array}[]{ccc}C&0\\ \end{array}\right].

Thus, the corresponding transfer function is defined as,

[TABLE]

We note that $\mathcal{F}(\omega)=F(\omega)$. In fact, setting

[TABLE]

wich gives $\mathcal{F}(\omega)=\mathcal{C}X$, where $X$ verifies $(\omega I_{2n}-\mathcal{A})X=\mathcal{B}$. Using the expressions of the matrices $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$, we get,

[TABLE]

Hence

[TABLE]

We can reduce the second-order system (33) by applying linear model reduction technique presented in the previous section, to $(A,\ B,\ C)$ to yield a small linear system $(A_{m},\ B_{m},\ C_{m})$. Unfortunately, there is no guarantee that the matrices defining the reduced system have the necessary structure to preserve the second-order form of the original system. For that we follow the model reduction techniques of second-order structure-preserving, presented in [4, 6, 7].

5.1 The structure-preserving of the second-order reduced model

Using the Krylov subspace-based methods discussed in the previous section do not guaranty the second-order structure when applied to the linear system (38). the authors in [4, 7] proposed a result, that gives a simple sufficient condition to satisfy the interpolation condition and produce a second order reduced system.

Lemma 5.1**.**

Let $(\mathcal{A},\ \mathcal{B},\ \mathcal{C})$ be the state space realization defined in (38). If we project the state space realization with $2n\times 2ms$ bloc diagonal matrices

[TABLE]

where $\mathcal{V}_{m}^{1},\ \mathcal{V}_{m}^{2},\ \mathcal{W}_{m}^{1}\text{and}\ \mathcal{W}_{m}^{2}\in\mathbb{R}^{n\times ms}$, then the reduced transfer function

[TABLE]

is a second order transfer function, on condition that the matrix $(\mathcal{W}_{m}^{2})^{T}\mathcal{V}_{m}^{1}$ is invertible.

Theorem 5**.**

Let $\mathcal{F}(\omega):=\mathcal{C}(\omega I_{2n}-\mathcal{A})^{-1}\mathcal{B}=C(\omega^{2}I_{n}+\omega D+K)^{-1}B$, with

[TABLE]

be a second order transfer function. Let $\mathbb{V}_{m},\ \mathbb{W}_{m}\in\mathbb{R}^{2n\times ms}$ be defined as:

[TABLE]

where $\mathbb{V}_{m}^{1},\ \mathbb{V}_{m}^{2},\ \mathbb{W}_{m}^{1}\ \text{and}\ \mathbb{W}_{m}^{2}\in\mathbb{R}^{n\times ms}$, with $(\mathbb{W}_{m}^{1})^{T}\mathbb{V}_{m}^{1}=(\mathbb{W}_{m}^{2})^{T}\mathbb{V}_{m}^{2}=I_{ms}$. Let us construct the $2n\times 2ms$ projecting matrices as

[TABLE]

Define the second order transfer function of order $m$ by

[TABLE]

if we have

[TABLE]

and

[TABLE]

where $\sigma_{i},\ \mu_{i}$, are the interpolation points, and $R_{i},\ L_{i}\in\mathbb{R}^{p\times s},$ for $i=1,...m$ are the tangential directions. Then the reduced order transfer function $\mathcal{F}_{m}(\omega)$ interpolates the values of the original transfer function $\mathcal{F}(\omega)$ and preserves the structures of the second-order model provided that the matrix $(\mathbb{W}_{m}^{2})^{T}\mathbb{V}_{m}^{1}$ is non-singular.

6 Numerical experiments

In this section, we present some numerical examples to show the effectiveness of the Adaptive Block Tangential Lanczos (ABTL) algorithm. All the experiments were carried out using the CALCULCO computing platform, supported by SCoSI/ULCO (Service Commun du Système d’Information de l’Unive-rsité du Littoral Côte d’Opale). The algorithms were coded in Matlab R2018a. We used the following functions from LYAPACK [24]:

• •

lp$\_$lgfrq: Generates a set of logarithmically distributed frequency sampling points.

• •

lp$\_$gnorm: Computes $\|H(j\omega)-H_{m}(j\omega)\|_{2}$.

We used various matrices from LYAPACK and from the Oberwolfach collection111Oberwolfach model reduction benchmark collection 2003. http://www.imtek.de/simulation/benchmark. These matrix tests are reported in Table 1 with different values of $p$ and the used values of $s$.

6.1 Example 1: FDM model

The finite differences semi-discretized heat equation will serve as the most basic test example here. Its corresponding matrix $A$, is obtained from the centered finite difference discretization of the operator,

[TABLE]

on the unit square $[0,1]\times[0,1]$ with homogeneous Dirichlet boundary conditions with

[TABLE]

The matrices $B$ and $C$ were random matrices with entries uniformly distributed in $[0,1]$. The dimension of $A$ is $n=n_{0}^{2}$, where $n_{0}^{2}$ is the number of inner grid points in each direction. The advantages of this model are:

• •

It’s easy to understand.

• •

The discretization using the finite difference method (FDM) is easy to implement.

• •

It allows for simple generation of almost arbitrary size test problems.

In Table 2, we compared the execution times and the ${\mathcal{H_{\infty}}}$ norm $\parallel H-H_{m}\parallel_{{\mathcal{H}}_{\infty}}$ of the ABTL algorithm with the Iterative Rational Krylov Algorithm (IRKA [30]) and the adaptive tangential method represented by Druskin and Simonsini (TRKSM) see for more details [12], with different values of $m$. We notice that the obtained timing didn’t contain the execution times used to obtain the errors. As can be seen from the results in Table 2, the cost of IRKA method is much higher than the cost required with the adaptive block tangential Lanczos method.

6.2 Example 2: Linear 1D Beam

Moving structures are an essential part for many micro-system devices, among them fluidic components like pumps and electrically controllable valves, sensing cantilevers, and optical structures. While the single component can easily be simulated on a usual desktop computer, the calculation of a system of many coupled devices still presents a challenge. This challenge is raised by the fact that many of these devices show a nonlinear behavior. This model describes a slender beam with four degrees of freedom per node: ”$x$ the axial displacement”, ”$\Theta_{x}$ the axial rotation”, ”$y$ the flexural displacement” and ”$\Theta_{z}$ the flexural rotation”. The model is from the Oberwolfach collection. The matrices are obtained by using the finite element discretization presented in [30]. We used two examples of linear 1D Beam model:

The Figures 2 above represent the norm of the original transfer function $\|H(j\omega)\|_{2}$ and the norm of the reduced transfer function $\|H_{m}(j\omega)\|_{2}$ versus the frequencies $\omega\in[10^{-3},\;10^{3}]$ of the 1Dbeam-LF50000 model and it is a second-order model of dimension $2\times n=39988$ with one input and one output. The Figure 2 represents the exact error $\|H(j\omega)-H_{m}(j\omega)\|_{2}$ versus the frequencies.

The plots in Figures 4 and 4, represent the original transfer function $\|H(j\omega)\|_{2}$ and the norm of the reduced transfer function $\|H_{m}(j\omega)\|_{2}$ where we modified the matrices B and C (random matrices) to get a MIMO system with four inputs and four outputs.

The plots in Figures 6 and 6, represent the exact error $\|H(j\omega)-H_{m}(j\omega)\|_{2}$ versus the frequencies $\omega\in[10^{-6},\;10^{6}]$ of the 1Dbeam-LF10000 model $2\times n=39996$ with one input and one output.

6.3 Example 3: Butterfly Gyroscope

The structural model of the gyroscope has been done in ANSYS (the global leader in engineering simulation) using quadratic tetrahedral elements. The model used here is a simplified one with a coarse mesh as it is designed to test the model reduction approaches. It includes the pure structural mechanics problem only. The load vector is composed from time-varying nodal forces applied at the centers of the excitation electrodes. The Dirichlet boundary conditions have been applied to all degree of freedom of the nodes belonging to the top and bottom surfaces of the frame. This benchmark is also part of the Oberwolfach Collection. It is a second-order model of dimension $n=17361$, (then the matrix $\mathcal{A}$ is of size $2\times n=34722$) with the matrix $B=C^{T}$ to get a MIMO system with 12 inputs and 12 outputs.

The plots in Figure 8 represent $\|H(j\omega)\|_{2}$ and the norm of the reduced transfer function $\|H_{m}(j\omega)\|_{2}$. Figure 8 represent the exact error $\|H(j\omega)-H_{m}(j\omega)\|_{2}$ versus the frequencies. The dimension of the reduced model is $m=40$. The execution time was 41.83 seconds with $\mathcal{H}_{\infty}$-err norms equal to $1.73\times 10^{-3}$.

The plots in Figure 10 represent $\|H(j\omega)\|_{2}$ and the norm of the reduced transfer function $\|H_{m}(j\omega)\|_{2}$. Figure 10 represent the exact error $\|H(j\omega)-H_{m}(j\omega)\|_{2}$ versus the frequencies with $m=80$. The execution time was 96.26 seconds with $\mathcal{H}_{\infty}$-err norms equal to $1.52\times 10^{-2}$.

we notice that the tree methods coincide, with an execution time almost the same (TRKSM: 9.26 seconds, ABTL: 9.76 seconds, ABTA: 10.45 seconds)

7 Conclusion

In the present paper, we proposed a new approach based on block tangential Krylov subspaces to compute low rank approximation of large-scale first and second order dynamical systems with multiple inputs and multiple outputs (MIMO). The method constructs sequences of orthogonal blocks from matrix tangential Krylov subspaces using the block Lanczos-type approach. The interpolation shifts and the tangential directions are selected in an adaptive way by maximizing the residual norms. We gave some new algebraic properties and compared our algorithms with well knowing methods to show the effectiveness of this latter.