Mixed control of vibrational systems

TL;DR

This paper introduces new $H_2$ norm-based performance measures for vibrational systems, demonstrating their advantages over standard measures through theoretical analysis and numerical experiments to optimize damping strategies.

Contribution

It proposes novel $H_2$ norm-based performance measures for vibrational systems and evaluates their effectiveness in optimizing damping configurations.

Findings

New measures outperform standard ones in certain scenarios

Optimal damping positions vary with new measures

Numerical experiments validate theoretical advantages

Abstract

We consider new performance measures for vibrational systems based on the $H_2$ norm of linear time invariant systems. New measures will be used as an optimization criterion for the optimal damping of vibrational systems. We consider both theoretical and concrete cases in order to show how new measures stack up against the standard measures. The quality and advantages of new measures as well as the behaviour of optimal damping positions and corresponding damping viscosities are illustrated in numerical experiments.

Tables

Table 1. Table 1 : Mass and stiffness parameters of the structure

$i$th index	1	2	3	4	5
$m_i$	4000	3000	2000	1000	800
$k_i$	$3.375\times {10}^6$	$3.75\times {10}^6$	$3.375\times {10}^6$	$3\times {10}^6$	$2.25\times {10}^6$

Equations

$$G=\{\begin{array}{rl}& M\ddot{q}+D\dot{q}+Kq=B_{2}u,\\ & q(0)=q_0,\dot{q}(0)={\dot{q}}_0,\\ & y=\left[\begin{array}{c}{C}_{1}q\\ C_2\dot{q}\end{array}\right].\end{array}$$

$$\Vert G{\Vert }_2={\left(\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{Tr}(\hat{G}(i\omega {)}^{\ast }\hat{G}(i\omega ))\mathrm{d}\omega \right)}^{1/2},$$

$$\mathrm{Tr}\left(\left[\begin{array}{cc}{C}_1^{\ast }{C}_2& 0\\ 0& C_2^{\ast }{C}_2\end{array}\right]X\right),$$

$$\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}D\end{array}\right]X+X{\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}D\end{array}\right]}^{\ast }=-\left[\begin{array}{cc}0& 0\\ 0& M^{-1}{B}_2{B}_2^{\ast }{M}^{-1}\end{array}\right].$$

$$\hat{y}=\hat{G}\hat{u}=\left[\begin{array}{c}{C}_1(s^{2}M+sD+K{)}^{-1}\\ C_{2}s(s^{2}M+sD+K{)}^{-1}\end{array}\right]B_2\hat{u}.$$

$$\hat{y}$$

$$=\left[\begin{array}{cc}{C}_1& 0\\ 0& C_2\end{array}\right]\left[\begin{array}{cc}\ast & (s^2+s{M}^{-1}D+M^{-1}K{)}^{-1}\\ \ast & s(s^2+s{M}^{-1}D+M^{-1}K{)}^{-1}\end{array}\right]\left[\begin{array}{c}0\\ M^{-1}{B}_2\end{array}\right]\hat{u}$$

$$=\left[\begin{array}{cc}{C}_1& 0\\ 0& C_2\end{array}\right]{\left(sI-\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}D\end{array}\right]\right)}^{-1}\left[\begin{array}{c}0\\ M^{-1}{B}_2\end{array}\right]\hat{u}.$$

$$\{\begin{array}{rl}\dot{x}& =\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}D\end{array}\right]x+\left[\begin{array}{c}0\\ M^{-1}{B}_2\end{array}\right]u,\\ y& =\left[\begin{array}{cc}{C}_1& 0\\ 0& C_2\end{array}\right]x.\end{array}$$

$$A=\left[\begin{array}{cc}0& T_1{T}_2^{-1}\\ -T_2{M}^{-1}K{T}_1^{-1}& -T_2{M}^{-1}D{T}_2^{-1}\end{array}\right],C=\left[\begin{array}{cc}{C}_1{T}_1^{-1}& 0\\ 0& C_2{T}_2^{-1}\end{array}\right],B=\left[\begin{array}{c}0\\ T_2{M}^{-1}{B}_2\end{array}\right].$$

$$\{\begin{array}{rl}\dot{x}& =Ax+Bu,\\ y& =Cx.\end{array}$$

$$\sqrt{\mathrm{Tr}(C^{\ast }CX)}\text{ where }AX+X{A}^{\ast }=-B{B}^{\ast }.$$

$$\min \{\Vert G(D){\Vert }_2:D\text{ positive semidefinite and }A(D)\text{ a stable matrix}\}$$

$$\inf \{\Vert G(D){\Vert }_2:D\text{ positive semidefinite and }A(D)\text{ a stable matrix}\}=0.$$

$$\tilde{A}=\left[\begin{array}{cc}0& \Omega \\ -\Omega & -\tilde{D}\end{array}\right],\tilde{C}=\left[\begin{array}{cc}{\tilde{C}}_1& 0\\ 0& {\tilde{C}}_2\end{array}\right],\tilde{B}=\left[\begin{array}{c}0\\ {\tilde{B}}_2\end{array}\right],$$

$$\Omega X_{12}(\alpha {)}^{\ast }+X_{12}(\alpha )\Omega =0$$

$$\Omega X_{22}(\alpha )-X_{11}(\alpha )\Omega -\alpha X_{12}(\alpha )\Omega =0$$

$$-\Omega X_{12}(\alpha )-\alpha \Omega X_{22}(\alpha )-X_{12}^{\ast }(\alpha )\Omega -\alpha X_{22}(\alpha )\Omega =-\tilde{B}{\tilde{B}}^{\ast }.$$

$$-\Omega X_{22}(\alpha )-X_{22}(\alpha )\Omega =-\frac{1}{\alpha }\tilde{B}{\tilde{B}}^{\ast }+\frac{1}{\alpha }\Omega X_{12}(\alpha )+\frac{1}{\alpha }{X}_{12}^{\ast }(\alpha )\Omega \to 0,$$

$$\mathrm{Tr}(X_{11}(\alpha )\Omega )=\mathrm{Tr}(\Omega X_{22}(\alpha ))-\alpha \mathrm{Tr}(X_{12}(\alpha )\Omega )=\mathrm{Tr}(\Omega X_{22}(\alpha ))\to 0\text{ as }\alpha \to \infty .$$

$$-\frac{1}{{\alpha }^{k-1}}\Omega X_{22}(\alpha )-\frac{1}{{\alpha }^{k-1}}{X}_{22}(\alpha )\Omega =-\frac{1}{{\alpha }^k}\tilde{B}{\tilde{B}}^{\ast }+\frac{1}{{\alpha }^k}\Omega X_{12}(\alpha )+\frac{1}{{\alpha }^k}{X}_{12}^{\ast }(\alpha )\Omega \to 0,$$

$$\sum _{k=1}^m\Vert y_i{\Vert }_2^2$$

$$=\sum _{k=1}^m\mathrm{Tr}(e_i^{\ast }{B}^{\ast }XB{e}_i)=\mathrm{Tr}(B^{\ast }XB)=\Vert G{\Vert }_2^2,$$

$$e(t)={\Vert \left[\begin{array}{c}{C}_{1}q(t)\\ C_2\dot{q}(t)\end{array}\right]\Vert }^2={\Vert C\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\end{array}\right]\Vert }^2={\Vert C{\mathrm{e}}^{At}T\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right]\Vert }^2={\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right]}^{\ast }{T}^{\ast }{\mathrm{e}}^{A^{\ast }t}{C}^{\ast }C{\mathrm{e}}^{At}T\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right].$$

$${\int }_0^{\infty }e(t)\mathrm{d}t={\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right]}^{\ast }{T}^{\ast }{\int }_0^{\infty }{\mathrm{e}}^{A^{\ast }t}{C}^{\ast }C{\mathrm{e}}^{At}\mathrm{d}tT\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right]={\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right]}^{\ast }{T}^{\ast }XT\left[\begin{array}{c}{q}_0\\ {\dot{q}}_0\end{array}\right],$$

$$A^{\ast }X+XA=-C^{\ast }C.$$

$$(T^{-1}AT{)}^{\ast }(T^{\ast }XT)+(T^{\ast }XT)(T^{-1}AT)=-\left[\begin{array}{cc}{C}_1^{\ast }{C}_1& 0\\ 0& C_2^{\ast }{C}_2\end{array}\right].$$

$$(T^{-1}AT{)}^{\ast }Y+Y(T^{-1}AT{)}^{\ast }=-\left[\begin{array}{cc}{C}_1^{\ast }{C}_1& 0\\ 0& C_2^{\ast }{C}_2\end{array}\right].$$

$$\Vert G{\Vert }_{2,\mathrm{hom}}^2:={\int }_{\Vert q_0{\Vert }^2+\Vert {\dot{q}}_0{\Vert }^2=1}{\int }_0^{\infty }e(t;q_0,{\dot{q}}_0)\mathrm{d}t\mathrm{d}\sigma ={\int }_{\Vert q_0{\Vert }^2+\Vert {\dot{q}}_0{\Vert }^2=1}{\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]}^{\ast }Y\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]\mathrm{d}\sigma .$$

$$Y\mapsto {\int }_{\Vert q_0{\Vert }^2+\Vert {\dot{q}}_0{\Vert }^2=1}{\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]}^{\ast }Y\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]\mathrm{d}\sigma$$

$${\int }_{\Vert q_0{\Vert }^2+\Vert {\dot{q}}_0{\Vert }^2=1}{\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]}^{\ast }Y\left[\begin{array}{c}{q}_0\\ \dot{q_0}\end{array}\right]\mathrm{d}\sigma =\mathrm{Tr}(Y{\hat{Z}}_{\sigma })=\mathrm{Tr}(T^{\ast }XT{\hat{Z}}_{\sigma })=\mathrm{Tr}(T{\hat{Z}}_{\sigma }{T}^{\ast }X).$$

$$\sqrt{\mathrm{Tr}(C^{\ast }CY)}\text{ where }AY+Y{A}^{\ast }=-Z_{\sigma }.$$

$$\sqrt{\mathrm{Tr}(C^{\ast }CX)},\text{ where }AX+X{A}^{\ast }=-p{Z}_{\sigma }-(1-p)B{B}^{\ast }.$$

Full text

1 Introduction

In this paper we are concerned with the minimization of vibrations of an abstract time invariant vibrational system [16] described by

[TABLE]

Here $M$ denotes the mass matrix, $D$ denotes the damping matrix and $K$ denotes the stiffness matrix. We assume that $M$ is a non-singular matrix and that matrices $M$, $D$ and $K$ are real matrices of order $n$. Vector $q\in\mathbb{R}^{n}$ contains the state variables, while the vector $y\in\mathbb{R}^{r}$ denotes the observed output which together with the output matrices $C_{1},C_{2}\in\mathbb{R}^{p\times n}$ determine the system displacements and velocities of interest. System disturbances are denoted by the vector $u\in\mathbb{R}^{m}$ and the matrix $B_{2}\in\mathbb{R}^{n\times m}$. The first equation in (1) is usually called the state equation.

Systems of the form (1) are used as a linearized model for a large class of vibrational systems. Vibrations are typical and mostly unwanted phenomenon in mechanical systems, since resonance and sustained oscillations can have undesired effects such as energy waste, noise creation and even structural damage. Thus, the minimization of vibrations is a widely studied topic. There is a vast literature in this field of research, particularly in the engineering and applied mathematics. For a brief insight we give just a few references: [2, 23, 24, 13, 33, 9, 15]. All these references contain the topic of the minimization of dangerous vibrations from different aspects and for different problems.

In our setting, all matrices except $D$ are fixed and our main focus is on the damping matrix $D$. The damping matrix $D$ can be modeled in several different ways. Usually it is modeled as a sum of internal and external damping, that is, $D=D_{\mathrm{int}}+D_{\mathrm{ext}}$, where $D_{\mathrm{int}}$ represents internal damping and the external damping part $D_{\mathrm{ext}}$ depends on positive real parameters $v_{i}$, $i=1,\ldots,k$ (called viscosities) and corresponding damping positions. This means that external damping encodes damping positions with corresponding damping viscosities. The internal damping $D_{\mathrm{int}}$ can be modeled in different ways, for example as a small multiple of the critical damping or the proportional damping. Overview of different options for modeling internal damping can be found in [17].

We are interested in choosing a damping matrix $D$ in such a way that the vibrations of the system are as small as possible. To be able to pose this problem as an optimization problem, we need to choose an optimization criterion.

In damping optimization setting different optimization criteria are used, depending on different applications. One criterion for systems described by (1) would be to minimize the $H_{2}$ norm of the system as the penalty function. In particular, damping optimization using the $H_{2}$ norm was considered in [7, 37] as well as in [25, 3], where the authors considered model order reduction approaches in order to determine the optimal damping parameters efficiently. Also, one can consider damping optimization using the $H_{\infty}$ norm. This can be also calculated efficiently using the model order reduction approaches, for more details see, e.g., [1, 6]. Moreover, damping optimization in mechanical systems with external force was also considered in [33], [28] and [18] where the authors used optimization criteria based on average energy amplitude and average displacement amplitude over the considered time. For closed systems, without external forces, there exists a multitude of optimization criteria. Overview of such criteria can be found, e.g., in [33] or [19]. Some of them are based on eigenvalues, such as spectral abscissa criterion (for more details, see, e.g., [12], [27], [21], [34], [14]), while other criteria are based on the total energy of the system, such as the total average energy. Total average energy was considered widely in the last two decades, more details can be found in, e.g., [31], [32], [29], [30], [26], [8] and [20]. The total average energy was considered also in [4] and [5], where the authors employed dimension reduction techniques that allowed efficient calculation of the total average energy.

In this paper we will show that for the system given by (1), the standard $H_{2}$ norm may lead to an optimization problem which is not well posed in general. As the main contribution of the paper we propose an alternative criterion, which can be seen as an introduction of a constraint or alternatively as a use of a mixed norm combining $H_{2}$ norm of the closed, homogeneous system with initial data and $H_{2}$ norm of the open, non-homogeneous systems without initial data. Moreover, using the proposed mixed norm as a criterion for the optimization, we show that the problem of global optimization problem in the case when $B_{2}$, $C_{1}$ and $C_{2}$ are full rank matrices has a unique solution which belongs in the class of modal damping matrices. The uniqueness of the solution and its affiliation to the class of modal matrices suggest that the new criterion should be a viable alternative to the standard approaches. We illustrate this also by some numerical examples.

The paper is organized as follows. In Section 2 we present a standard criterion which is based on $H_{2}$ norm of a vibrational system and we show its drawbacks. In Section 3 we construct an $H_{2}$ norm for the homogeneous system while Section 4 considers new criterion which is based on mixed $H_{2}$ norms of homogeneous and non-homogeneous vibrational systems. In Section 5 we compare considered criteria and illustrate advantages of the new criterion as well as the optimal parameters that may arise in different cases.

2 $H_{2}$ norm of a vibrational system

Since the $H_{2}$ norm of a dynamical system will play the central role in the paper, we start with its definition.

Definition 2.1**.**

Suppose that the linear time invariant system $G$ is defined by (1). Then the $H_{2}$ norm of this system, denoted by $\lVert G\rVert_{2}$, is defined as

[TABLE]

where $\hat{G}$ is the transfer function of the system $G$, i.e. after the application of the Laplace transform to (1) we arrive at the equation $\hat{y}(s)=\hat{G}(s)\hat{u}(s)$.

Note that the $H_{2}$ norm does not depend on the initial data $q_{0}$, $\dot{q}_{0}$, hence in the case of the $H_{2}$ norm we can always assume $q_{0}=0$, $\dot{q}_{0}=0$.

By $\lVert\cdot\rVert_{2}$ we will denote both the $H_{2}$ norm of a system and the $L_{2}$ norm of the function. The precise meaning will always be clear from the context.

To use control–theoretic methods on (1), we can either apply Definition 2.1 or we can linearise (1) and then calculate $H_{2}$ norm of the corresponding linearised system. By linearization here we mean writing the system as a first order matrix ODE. In sequel, the term linearization will always have this meaning. Linearised system is then given in the standard form $\dot{x}=Ax+Bu,y=Cx$. In the following proposition we show that in both cases we obtain the same expression for the $H_{2}$ norm. This is a folklore result, we give a short proof for the reader’s convenience.

Proposition 2.2**.**

Assume that the system (1) has a finite $H_{2}$ norm. Then the square of the $H_{2}$ norm of the system (1) is given by

[TABLE]

where $X$ is the solution of the Lyapunov equation

[TABLE]

Proof.

Applying the Laplace transform to (1) we obtain

[TABLE]

Then

[TABLE]

Hence the system (1) has the same transfer function as the system

[TABLE]

Since the system (4) can be obtained from (1) by the linearization $x_{1}=q$, $x_{2}=\dot{q}$ and the multiplication to the left of the corresponding state equation with the matrix $\operatorname{diag}(0,M^{-1})$, the statement follows from the well–known state–space formula for the $H_{2}$ norm (see, for example, [36, Lemma 4.4]). ∎

Remark 2.3*.*

Note that by the assumption of Proposition 2.2, the system matrix $\left[\begin{smallmatrix}0&I\\ -M^{-1}K&-M^{-1}D\end{smallmatrix}\right]$ is stable and hence the matrix $X$ is symmetric positive semidefinite.

The choice of different linearizations of (1) amounts to differnt state transformations of the system (4). Indeed, for regular matrices $T_{1}$ and $T_{2}$ let $x_{1}=T_{1}q$, $x_{2}=T_{2}\dot{q}$ be a linearization of (1). Then one calculates the corresponding matrices

[TABLE]

Now, for system matrices from (4) we have $A=T\begin{bmatrix}0&I\\ -M^{-1}K&-M^{-1}D\end{bmatrix}T^{-1}$, $B=T\begin{bmatrix}0\\ M^{-1}B_{2}\end{bmatrix}$ and $C=\begin{bmatrix}C_{1}&0\\ 0&C_{2}\end{bmatrix}T^{-1}$, where $T=\mathrm{diag}(T_{1},T_{2})$. Sometimes we will write $A=A(D)$ and $G=G(D)$ to denote that the matrix $A$ and the system $G$ depend on $D$. Thus, in general, our linear time invariant system can be written as

[TABLE]

Hence, using [36, Lemma 4.4], the $H_{2}$ norm of the system (1) can be calculated as

[TABLE]

In the next proposition we show that the optimization problem $\min_{D}\lVert G(D)\rVert_{2}$ does not in general admit a solution. We treat a special case of a vibrational system without gyroscopic forces and assume that the damping is passive, so no additional energy is introduced into the system. These assumptions imply that $M$ and $K$ are positive definite symmetric matrices and that matrices $D$ are positive semidefinite symmetric matrices. We additionally assume that the matrices $A(D)$ are stable, so that the Lyapunov equation from (7) always has a solution.

In the sequel, the relation $X\geq 0$ will be used to denote the fact that the symmetric matrix $X$ is positive semidefinite.

Proposition 2.4**.**

Assume that $M$ and $K$ are symmetric positive definite matrices. Assume also that $C_{2}B_{2}\neq 0$. Then the optimization problem

[TABLE]

does not have a solution and moreover

[TABLE]

Proof.

Since the choice of linearization does not change the $H_{2}$ norm, we will work with a particulary convenient linearization given in (8) below. Let $M=L_{2}L_{2}^{*}$ and $K=L_{1}L_{1}^{*}$ be Cholesky factorizations of $M$ and $K$ and let $L_{2}^{-1}L_{1}=U_{2}\Omega U_{1}^{*}$ be a SVD decomposition of $L_{2}^{-1}L_{1}$. Note that the diagonal elements of $\Omega=\mathrm{diag}(\omega_{1},\ldots,\omega_{n})$ are square roots of the eigen–frequencies of the corresponding undamped system (with $D=0$). We choose $\widetilde{T}=\mathrm{diag}(\widetilde{T}_{1},\widetilde{T}_{2})=\mathrm{diag}(U_{1}^{*}L_{1}^{*},U_{2}^{*}L_{2}^{*})$. Then we obtain for the corresponding system matrices

[TABLE]

where $\widetilde{D}=U_{2}^{*}L_{2}^{-1}DL_{2}^{-*}U_{2}$, $\widetilde{C}_{1}=C_{1}L_{1}^{-*}U_{1}$, $\widetilde{C}_{2}=C_{2}L_{2}^{-*}U_{2}$ and $\widetilde{B}_{2}=U_{2}^{*}L_{2}^{-1}B_{2}$.

We will show that for the damping matrices $\widetilde{D}(\alpha)=\alpha\Omega$ the $H_{2}$ norm of the corresponding system $G(\alpha):=G(\widetilde{D}(\alpha))$ tends to zero when $\alpha\to\infty$. First note that $\widetilde{A}(\widetilde{D}(\alpha))$ is stable since $\widetilde{D}(\alpha)$ is positive definite (see, for example, [33, Corollary 15.7]). Let us partition the solution of the Lyapunov equation $\widetilde{A}(\widetilde{D}(\alpha))X+X\widetilde{A}(\widetilde{D}(\alpha)^{*}=-\widetilde{B}\widetilde{B}^{*}$ as $2\times 2$ block matrix with $n\times n$ entries $X_{ij}(\alpha)$, which depend on $\alpha$. Then we have

[TABLE]

First we assume that $\alpha^{-1}X_{12}(\alpha)\to 0$ as $\alpha\to\infty$. Then (11) reads

[TABLE]

hence $X_{22}(\alpha)\to 0$ as $\alpha\to\infty$. Since $X(\alpha)\geq 0$ it follows $X_{12}(\alpha)\to 0$ as $\alpha\to\infty$. From (9) it follows that $\operatorname{Tr}(X_{12}(\alpha)\Omega)=0$, and from (10)

[TABLE]

But this implies $X_{11}(\alpha)\to 0$, hence $X(\alpha)\to 0$ as $\alpha\to\infty$. Then $\lVert G(\alpha)\rVert_{2}\to 0$ as $\alpha\to\infty$.

In the other case, when $\alpha^{-1}X_{12}(\alpha)\nrightarrow 0$ as $\alpha\to\infty$, let $k\in\mathbb{N}$ be the smallest integer such that $\alpha^{-k}X_{12}(\alpha)\to 0$ as $\alpha\to\infty$. Such $k$ exists since $X(\alpha)$ is a rational function of $\alpha$, which follows from (7). Since $\alpha^{-k+1}X_{12}(\alpha)\nrightarrow 0$ as $\alpha\to\infty$, the non–negativity of $X(\alpha)$ implies that $\alpha^{-k+1}X_{22}(\alpha)\nrightarrow 0$ as $\alpha\to\infty$. Dividing (11) with $\alpha^{k}$ we obtain

[TABLE]

as $\alpha\to\infty$. But this implies $\alpha^{-k+1}X_{22}(\alpha)\to 0$ as $\alpha\to\infty$, a contradiction. To finish the proof, note that $\lVert G(D)\rVert_{2}=0$ if and only if $G(D)=0$ which implies that $CA^{k}B=0$, $k=0,1,\ldots$ (see, for example, [10, Lemma 2.26]). This implies $C_{2}B_{2}=0$, a contradiction to our assumption $C_{2}B_{2}\neq 0$. ∎

Remark 2.5*.*

One way to look at Proposition 2.4 is to note that one can interpret the $H_{2}$ norm of the system as the square root of the sum of the norms of the system responses to the initial data of the form $x(0)=Be_{i}$ with zero external force, where $e_{i}$ denotes the $i$th canonical vector. Indeed, in that case the corresponding system responses are $y_{i}(t)=C\mathrm{e}^{At}Be_{i}$, so

[TABLE]

where $X$ solves the Lyapunov equation $A^{*}X+XA=-C^{*}C$. Hence in our case, the $H_{2}$ norm only measures the responses to initial data only consisting of velocities, initial displacements do not play any role. This seems to be a general issue when the (first order) control system is obtain by a linearization from the higher order systems.

Remark 2.6*.*

The assumption $C_{2}B_{2}\neq 0$ could be relaxed, the system (1) is zero-system under a much stronger condition than $C_{2}B_{2}=0$. Since the precise statement and its proof are complicated, we will not formulate it here.

3 $H_{2}$ norm of a homogeneous system

In this section we will generalize the total energy approach for the measurement of unwanted vibrations of a homogeneous vibrational system which is in a way counterpart to the $H_{2}$ norm of the system given in preceding section.

For the system given by (1) we take $u=0$ but now the initial conditions $q(0)=q_{0}$, $\dot{q}(0)=\dot{q}_{0}$ will play a role. With $e(t;q_{0},\dot{q}_{0})=\lVert y(t;q_{0},\dot{q}_{0})\rVert^{2}$ we denote the energy of the output of the system. We want to average the total energy of the output, given by $\int_{0}^{\infty}e(t)\,\mathrm{d}t$ over all initial data $(q_{0},\dot{q}_{0})$ on the corresponding unit sphere.

We calculate, using a linearization given by $T=\mathrm{diag}(T_{1},T_{2})$,

[TABLE]

Now

[TABLE]

where $X$ is the solution of the Lyapunov equation

[TABLE]

Multiplying this equation on the left with $T^{*}$ and on the right with $T$ we obtain

[TABLE]

Let $Y=T^{*}XT$. Then $Y$ is the solution of the Lyapunov equation

[TABLE]

Let us choose $\sigma$, a surface measure on the unit sphere $\mathbb{R}^{2n}$. The $H_{2}$ norm of the homogeneous system is defined by

[TABLE]

Since

[TABLE]

is a linear functional on the space of symmetric matrices, by the Riesz representation theorem there exists a symmetric positive semidefinite matrix $\hat{Z}_{\sigma}$ ([33, Proposition 21.1]) such that

[TABLE]

Note that if we chose the Lebesgue measure for $\sigma$, then $\hat{Z}_{\sigma}=\frac{1}{2n}I$. Here by the Lebesgue measure on the unit sphere we mean the surface measure obtained by the Minkowski formula taking Lebesgue measure as the ambient measure, for the comprehensive treatment see [11]. A formula for the matrix $\hat{Z}_{\sigma}$ when the measure $\sigma$ is Gaussian is given in [22].

Let us denote $Z_{\sigma}=T\hat{Z}_{\sigma}T^{*}$. Then succinctly the $H_{2}$ norm of the homogeneous system can be calculated as $\lVert G\rVert_{2,\mathrm{hom}}=\sqrt{\operatorname{Tr}(Z_{\sigma}X)}$, where $X$ solves (12). In the sequel it will be useful to calculate the same norm by the use of the dual Lyapunov equation:

[TABLE]

Note that (7) and (13) differ only in the right hand side of the Lyapunov equations, where right hand sides encode the information about dangerous external forces and initial conditions, respectively.

4 Mixed control of vibrational systems

In this section we will combine the two approaches given in Sections 2 and 3 to create a norm which takes into account both external forces and initial data and we will show that using this norm as a criterion for the optimization we can prove that the global minimum exists and is obtained in the class of modal damping matrices.

The issue with (13) is that it does not carry any information about the external forces, and the issue with (7) is that it does not carry all the needed information about the initial data. A natural choice is to try to combine these two norms by taking their (generalized) quadratic mean. Let $0<p<1$. We define the $p$–mixed $H_{2}$ norm of the system $G$ by $\lVert G\rVert_{2,p}^{2}=(1-p)\lVert G\rVert_{2}^{2}+p\lVert G\rVert_{2,\mathrm{hom}}^{2}$. Using (7) and (13), we see that the $p$–mixed $H_{2}$ norm can be calculated as

[TABLE]

Note that (14) does not depend on the choice of the linearization. One can also think of (14) as the standard $H_{2}$ norm with an additional geometric constraint. Indeed, we can think of the function $G\mapsto p\lVert G\rVert_{2,\mathrm{hom}}^{2}$ as a barrier function, constraining the set of feasible systems $G$.

Let us now assume that $M$ and $K$ are positive definite symmetric matrices. When modeling vibrational systems, the matrix $B_{2}$ is usually designed as a band–pass filter where only the dangerous frequencies are passed. In the linearization given in (8), this would mean that $\widetilde{B}_{2}=Z_{1}:=\mathrm{diag}(I_{r},0_{n-r})$, where we assumed that we chose SVD decomposition $L_{2}^{-1}L_{1}=U_{2}\Omega U_{1}^{*}$ in such a way that the dangerous frequencies of the undamped system are exactly $\omega_{1},\ldots,\omega_{r}$. Hence in that case $B_{2}=L_{2}U_{2}Z_{1}$.

The measure $\sigma$ typically is chosen in such a way that it attenuates frequencies which are not dangerous. In particular, if the surface measure is chosen in such a way that it corresponds to Lebesgue measure on the subspace spanned by the vectors $[x_{i},0]^{\operatorname{T}}$ and $[0,x_{i}]^{\operatorname{T}}$, $i=1,\ldots,r$, where $x_{i}$ are the eigenvectors of the first $r$ undamped eigenfrequencies and on the rest of $\mathbb{R}^{2n}$ it corresponds to the Dirac measure concentrated at zero, then, in the linearization given by (8), we have $\widetilde{Z}_{\sigma}=\frac{1}{2n}Z$, where $Z=\mathrm{diag}(Z_{1},Z_{1})$, $Z_{1}=\mathrm{diag}(I_{r},0_{n-r})$.

Since the output $y$ usually corresponds to the energy corresponding to the unwanted vibrations, a usual choice for matrices $C_{1}$ and $C_{2}$ is $C_{1}=\frac{1}{\sqrt{2}}L_{1}^{*}$, $C_{2}=\frac{1}{\sqrt{2}}L_{2}^{*}$. Then $\lVert y(t)\rVert^{2}$ equals the energy of the system at the time $t$. This corresponds to $\widetilde{C}_{1}=\frac{1}{\sqrt{2}}U_{1}$, $\widetilde{C}_{2}=\frac{1}{\sqrt{2}}U_{2}$ and hence $\widetilde{C}^{*}\widetilde{C}=\frac{1}{2}I$. If we are (as is usually the case) interested only in the energy output corresponding to a part of the system, then we have $\widetilde{C}^{*}\widetilde{C}=\frac{1}{2}Z$. But in some cases there would be other appropriate choices for $C_{1}$ and $C_{2}$. One natural setting is when the matrix $C$ is such that $\frac{\mathrm{d}}{\mathrm{d}t}\lVert y(t)\rVert\leq 0$, which is equivalent to the assumption that $C^{*}CA$ is a stable matrix. Note that this is not automatically satisfied as the following example shows.

Example 4.1*.*

Let

[TABLE]

Then $A$ is stable, but $SA$ is not.

From the practical point of view, it is important to note that $C$ usually will have low rank.

One typical situation is that $B_{2}$, $C$ and the measure $\sigma$ are such that

[TABLE]

If we drop the constants which are not relevant for the optimization purposes, we obtain

[TABLE]

The next theorem shows that one can calculate the global minimum in this particular case.

Theorem 4.2**.**

Assume that $M$ and $K$ are symmetric positive definite matrices. Let $Z_{1}=I$ in (16). Let

[TABLE]

Then for all $0<p<1$ there exists a unique global minimum of the following optimization problem:

[TABLE]

The minimum is attained at $\widetilde{D}=\sqrt{\frac{2(1+p)}{p}}\Omega$.

Proof.

To emphasize the dependence of $X$ to the parameter $\widetilde{D}$ we will sometimes write $X(\widetilde{D})$.

Let $\widetilde{D}\in\mathcal{D}_{s}$ be arbitrary. From [33, Corollary 15.7] it follows that $\widetilde{D}\in\mathcal{D}_{s}$ if and only if all diagonal entries of $\widetilde{D}$ are non-zero. Let $Z_{(i)}$ be a diagonal matrix with all diagonal entries zero except the $i$-th which is $1$. Set $Z^{i}=\left(\begin{smallmatrix}pZ_{(i)}&0\\ 0&Z_{(i)}\end{smallmatrix}\right)$. Let $X_{i}$ be the solution of the Lyapunov equation

[TABLE]

Then it is easy to see that the solution of the Lyapunov equation $\widetilde{A}(\widetilde{D})X+X\widetilde{A}(\widetilde{D})^{\ast}=-\left[\begin{smallmatrix}pI&0\\ 0&I\end{smallmatrix}\right]$ is

[TABLE]

Our aim is to show

[TABLE]

Here the right hand side is chosen in such a way that when we sum it for $i=1,\ldots,n$, we obtain the minimal value of the optimization problem. Observe that by simple permutation argument we can assume $i=1$. Let us decompose a matrix $X_{1}\in\mathbb{R}^{2n\times 2n}$ in the following way:

[TABLE]

where $x_{11},x_{33},x_{13}\in\mathbb{R}$, $X_{12},X_{14},X_{34}\in\mathbb{R}^{1\times(n-1)}$, $X_{22},X_{24},X_{44}\in\mathbb{R}^{(n-1)\times(n-1)}$, and $X_{23}\in\mathbb{R}^{(n-1)\times 1}$. Next we partition the Lyapunov equation

[TABLE]

in the same way as we did the matrix $X_{1}$. We obtain

[TABLE]

[TABLE]

where $\omega_{1},\widetilde{d}_{11}\in\mathbb{R}$, $\widetilde{D}_{12}\in\mathbb{R}^{1\times(n-1)}$, and $\widetilde{D}_{22},\Omega_{n-1}\in\mathbb{R}^{(n-1)\times(n-1)}$.

From (1,1) we obtain $x_{13}=-\frac{p}{2\omega_{1}}$. Since $\widetilde{D}\geq 0$, one can easily see that $\widetilde{d}_{11}=0$ implies $\widetilde{D}_{12}=0$, hence (3,3) reads $p=-1$, a contradiction. Hence, $\widetilde{d}_{11}>0$. From (3,3) we now get

[TABLE]

The relation (4,4), together with the facts $X_{44}\geq 0$, $\widetilde{D}_{22}\geq 0$, implies

[TABLE]

and the relation (2,2) implies $\operatorname{Tr}(X_{24}\Omega_{n-1})=0$, hence we obtain

[TABLE]

From the relation (1,3) we obtain

[TABLE]

From relation (2,4) we obtain

[TABLE]

hence

[TABLE]

From the relation (2,2) we obtain

[TABLE]

where $S\in\mathbb{R}^{(n-1)\times(n-1)}$ is a skew–symmetric matrix.

Hence

[TABLE]

From the relation (1,2) follows $X_{23}=-\frac{1}{\omega_{1}}\Omega_{n-1}X_{14}^{\ast}$, hence

[TABLE]

Now (21) and (22) imply

[TABLE]

The last inequality follows from the following observation. Let us define the function $g(x)=\frac{1+p}{x}+\frac{px}{2\omega_{1}^{2}}$. Then the function $g$ attains its unique minimum $\frac{\sqrt{2p(1+p)}}{\omega_{1}}$ in $x=\sqrt{\frac{2(1+p)}{p}}\omega_{1}$.

Hence, we have shown (19). Now (18) implies

[TABLE]

Since

[TABLE]

this is indeed the global minimum.

Assume that $\widetilde{D}\in D_{s}$ is such that $\operatorname{Tr}(X(\widetilde{D}))=\sqrt{2p(1+p)}\sum_{i=1}^{n}\omega_{i}^{-1}$. Then (23) and (18) imply $\operatorname{Tr}X_{1}=\sqrt{2p(1+p)}\omega_{1}^{-1}$. Observe that the matrix $X_{1}$ is decomposed as in (20). Then (23) implies $X_{44}=0$. Since $X_{1}\geq 0$, it follows $X_{14}=X_{24}=X_{34}=0$. From the relation (1,2) it follows $X_{23}=0$, from relation (2,4) it follows $X_{22}=0$, and from relation (2,3) it follows $X_{12}=0$. Finally, from (1,4) now it follows $\widetilde{D}_{12}=0$.

By repeating this procedure for $i=2,\ldots,n$ we obtain that $\widetilde{D}$ is a diagonal matrix. From (23) it follows that $\widetilde{D}=\alpha\Omega^{-1}$ for some $\alpha$ and then it is easy to check that $\alpha=\sqrt{2(1+p)}/\sqrt{p}$. ∎

Remark 4.3*.*

Theorem 4.2 shows that the $p$–mixed $H_{2}$ norm should be well suited for the use in damping optimization problems for vibrational systems, at least in the case when there are no gyroscopic forces. The unique global minimizer is obtained in the class of the modal matrices, which is a very natural result.

5 Numerical experiments

In this section, we consider numerical examples in order to illustrate the quality and advantages of new performance measures as well as the comparison with standard performance measures. In these examples the corresponding Lyapunov equation was solved by Matlab’s function lyap.

Example 5.1*.*

We consider an $n$-mass oscillator or an oscillator ladder given by 1. The oscillator describes the mechanical system of $n$ masses and $n+1$ springs with two grounded dampers. Similar models were considered e.g. in [4], [5], [27] and [33]. For this mechanical system, the mathematical model is given by (1), where the mass and stiffness matrices are defined by

[TABLE]

Coefficients of mass and stiffness matrices are given as

[TABLE]

The damping matrix is given as $D=D_{\mathrm{int}}+D_{\mathrm{ext}},$ where the internal damping $D_{int}$ is small multiple of the critical damping, $D_{\mathrm{int}}=0.04\cdot M^{1/2}\sqrt{M^{-1/2}KM^{-1/2}}M^{1/2}$. Here $\sqrt{X}$ denotes the unique symmetric positive definite matrix such that $(\sqrt{X})^{2}=X$.

We consider two grounded dampers with the positions $i$ and $j$ that have the viscosities $v_{1},v_{2}\geq 0$, respectively. That is, the external damping is defined by $D_{ext}=v_{1}e_{i}e_{i}^{T}+v_{2}e_{j}e_{j}^{T}$, where $1\leq i<j\leq n$. We are interested in the damping of all undamped eigenfrequencies.

In order to compare optimal positions when considering different performance measures, we will calculate the optimal viscosities for all possible damping configurations where feasible interval for each viscosity will be $[\,0,\,5000\,]$. The viscosities were optimized by Matlab’s function fminsearchbnd and they are rounded to two digits.

Regarding the input and the output matrices we consider two cases. The first case considers one particular choice of the input and output matrices, while the second case considers matrices given by (15). In both cases we will calculate new performance measure using the linearization given by (8).

Let us consider the first case where for the penalty function we will use $p$–mixed $H_{2}$ norm of the system defined by (14).

Here the input matrix $B_{2}\in\mathbb{R}^{n\times 5}$ from (1) is such that

[TABLE]

while all the other entries are equal to zero. This means that the input is applied to the first 5 masses on the left-hand side of the considered $n$-mass oscillator, hence the masses on the left-hand side od the oscillator have a larger influence on the input.

Moreover, we are interested in the 10 displacements and velocities in the middle of the $n$-mass oscillator. This means that the matrices $C_{1},C_{2}\in\mathbb{R}^{10\times n}$ from (1) are defined by

[TABLE]

with all other entries being equal to zero.

First, in Figure 2 we illustrate the behaviour of the optimal parameters and the magnitude of the $p$–mixed $H_{2}$ norm of the system defined by (14). In more details, for given damping positions $i,j\in\{1,2,\ldots,n\}$ with $i<j$, we have optimized the viscosities and the optimal $p$–mixed $H_{2}$ norm is plotted. Figure 2 presents four different subplots that correspond to parameters $p=0,\frac{1}{3},\frac{2}{3},1$. Here, for each damping configuration $(i,j)$ where $i<j$, we have optimized the viscosities and the optimal viscosities are denoted by $(v_{1}^{0},v_{2}^{0})$ while the optimal position for given $p$ is denoted by $(i^{0},j^{0})$. In this figure we can also see the magnitudes of optimal damping parameters for considered parameters $p$.

Figure 3 shows similar results, but here on the $z$ axis instead of the $p$–mixed $H_{2}$ norm we present the optimal first viscosity. This means that for each damping configuration $(i,j)$ where $i<j$, we have optimized viscosities $v_{1}$ and $v_{2}$ and the optimal first viscosity ($v_{1}^{0}$) is plotted. Similarly we can plot the optimal second viscosity, but surface plots illustrate a similar behaviour.

From Figures 2 and 3 we can conclude that the parameter $p$ has a very strong influence on optimal damping parameters. First, we can see that optimal damping positions vary significantly, while we change the parameter $p$ and the difference in optimal viscosities is even bigger. Moreover, for $p=0$ from Figure 3 (which corresponds to the standard $H_{2}$ norm), we can see that there exists a whole area where the optimal first viscosity was equal to $5000$ (the right hand side of the feasible segment), which means that the optimal viscosity should be as large as possible.

Let us now consider the second case where we consider the new performance measure given by (16).

Now, we will consider the same example, but thew optimization criterion will be based on matrices given by (15). This means that we will consider the criterion given by (16).

The same optimization process was performed as in the first case. First of all, we should emphasize that we did not obtain a significant difference in optimal damping parameters. In particular, for all parameters $p=0,\frac{1}{3},\frac{2}{3},1$ the optimal configuration of damping positions was equal to $(27,53)$. On the other hand, we observe a slight change in optimal parameters, that is, optimal viscosities are (234.57,222.08), (229.05,217.41), (225.99,214.72) and (224.01,213.06) for $p=0,\frac{1}{3},\frac{2}{3},1$, respectively.

In this case there isn’t a significant difference in the surface plot for this penalty function. Thus, on Figure 4 we present surface subplots obtained by using the parameter $p=\frac{1}{3}$. On the first subplot we present the surface plot for the new performance measure given by (16), while on the second subplot we show magnitudes of the optimal first viscosity.

*In this case our performance measure given by (16) is not significantly influenced by a change of the parameter $p$. Hence in this example the optimal damping parameters do not depend on the particular choice of the performance measure as long as $p>0$. This shows that for some vibrational systems the choice of the parameter $p$ is not important. *

Example 5.2*.*

In the second example we will consider a five-story shear frame structure from [35] which is shown on Figure 5. Similarly to the first example, in this example the structure can be modeled as a lumped mass system with equation (1) where

[TABLE]

Configuration of masses and stiffnesses in this example corresponds to a structure from [35] with mass and stiffness parameters given in Table 1.

The damping matrix is modeled as $D=D_{\mathrm{int}}+D_{\mathrm{ext}},$ where the internal damping $D_{int}$ is again a small multiple of the critical damping as it was in the first example. In this example we consider a damper that connects second and third floor of the five-story shear frame structure, therefore external damping is defined by $D_{ext}=v(e_{2}-e_{3})(e_{2}-e_{3})^{T}$, where the parameter $v$ represents the viscosity parameter.

Similarly to the previous example, we are interested in the damping of all undamped eigenfrequencies.

Here the input matrix from (1) is determined by the matrix B_{2}=\left(\begin{array}[]{ccccc}5000&0&0&0&0\end{array}\right)^{T} which means that the input is applied on the mass that is closest to the ground. This is natural since that this corresponds to disturbance that comes from ground. On the other hand, we will be interested in the stabilization of oscillations at the highest floor and therefore we are interested in the displacement and velocity of the highest mass, hence we take

[TABLE]

In Figure 6 we show results for the new performance measure where the penalty function corresponds to the $p$–mixed $H_{2}$ norm of the system defined by (14). The parameter $p$ varies from [math] to $1$ and the first subplot shows the magnitude of the optimal viscosity parameter, while the second subplot shows the magnitude of the $p$–mixed $H_{2}$ norm at the optimal viscosity. We can see from this figure the behaviour of the optimal value and the optimal function value. In particular, the optimal viscosity varies from $1.09\cdot 10^{5}$ to $1.44\cdot 10^{5}$ which means that the parameter $p$ has significant influence on the magnitude of optimal viscosities.

**Acknowledgments

**Supported in part by the National Science Foundation under the project “Control of Dynamical Systems”, Grant No. IP-2016-06-2468 and under the project “Optimization of parameter dependent mechanical systems” (IP-2014-09-9540), Grant No. 9540.