Balanced-based model reduction of uncertain systems with interval parameters

Abstract

Model reduction is a significant issue in dynamic system simulation and control, as a consequence of the unmanageable levels of storage and computational requirements for large-scale systems. In this paper, the concept of a balanced truncation approximation method is extended to large-scale systems with interval uncertainties to get the reduced-order model with uncertainties. In order to get the balanced system, the balancing transformation matrix is introduced by using the nominal system, and the reduced-order model with uncertainties is obtained by using balanced truncation. A major characteristic of this model reduction method is that the reduced-order model obtained in this way is also as uncertain as the original model. The closeness of the reduced-order model to the original model relies on the upper bounds of the ignored Hankel singular values. To compare the original model and the reduced-order model, a perturbation method is proposed to give the interval bounds of the responses of the original model and the reduced-order model. As applications of the proposed method, three numerical examples are given.

Keywords

Model reduction balanced truncation uncertainty interval analysis interval perturbation

1. Introduction

The state-space method based on the finite element method (FEM) is popular in the field of dynamics and control, for example structural dynamics analysis, fluid-solid coupling vibration analysis, vibro-acoustic analysis and active control of vibration and noise. However, for a complete, large-scale finite element model, it is time-consuming and inefficient to obtain the structural responses to dynamic loads. Therefore, an approximate method, which can significantly improve the computational efficiency, is necessary (Antoulas, 2005). Model reduction is one such model approximation method, looking for a low-order model whose dynamic characteristics are close to the original model (Besselink et al., 2013).

A lot of research about model reduction techniques has been done in the dynamics and control field where the original complex models are completely deterministic (Dolgin and Zeheb, 2005; Ghommem et al., 2012). The first model reduction method is Guyan reduction, which is exact for a static model (Guyan, 1965). An improved reduced system (IRS) was proposed by O’Callaghan (1989) to improve the static reduction method. Friswell et al. (1995, 1998) proposed the iterative scheme to improve the accuracy of the IRS and proved its convergence. Xia and Lin (2004a, 2004b) revised the iterative IRS and proved that it converges more rapidly. However, in systems and control field, the most popular model reduction approach is the balanced truncation method, which was first introduced by Mullis and Roberts (1976) and developed by Moore (1981). Essentially, the balanced truncation method is selecting the states that have the largest contribution to the system response and neglecting the states that have the smallest contribution to the system response. The most notable feature of this method is that it can keep the stability of the original complex model and give an error bound for quality assessment of the reduced-order model. The closeness of the reduced-order model to the original model can be measured by this error bound. Another two related model reduction methods are optimal Hankel norm approximation and balanced singular perturbation approximation, which are built on singular value decomposition (Kokotovic et al., 1976; Glover, 1984; Liu and Anderson, 1989). Like the balanced truncation method, proper orthogonal decomposition is widely used to deal with model reduction problem in dynamics analysis, which focuses on nonlinear partial differential equations (Holmes et al., 1998; Bellizzi and Sampaio, 2009; Khattak et al., 2010; Ritto et al., 2011). In addition to the balanced truncation method, the modal truncation method is also a widespread reduction method. Modal truncation and modal superposition methods are used together in structural dynamic analysis (Rayleigh, 1896; Craig and Kurdila, 2006). In recent years, many methods have been proposed to improve the accuracy and efficiency of the modal truncation method, such as mode acceleration (Wilson et al., 1982). The state residualization method is another model reduction method based on singular perturbation approximation (Benner et al., 2000). The model reduction methods described above are all based on singular value decomposition. There are additional model reduction methods based on Krylov, which are different from the methods based on singular perturbation approximation (Antoulas et al., 2001). More and more scholars (Gugercin et al., 2003; Gallivan et al., 2004; Antoulas, 2005; Gugercin, 2008; Komarasamy et al., 2012) are researching how to combine the Krylov-based methods and the singular value decomposition-based method to get a novel model reduction method. However, due to manufacturing errors, modeling errors, measurement errors as well as other factors, the parameters of the state-space model are generally inexact. Therefore, the dynamic behaviors of the systems are also uncertain. A model reduction method that is able to consider the influences of those uncertainties, should be studied further. That is, the reduced-order model obtained by a model reduction method considering the uncertainties should not only characterize the nominal value of the dynamic behaviors of the original model, but also represent the uncertainties of the dynamic behaviors.

Few papers on model reduction technology considering the uncertainties of the high-order model have been studied. Qiu et al. (1996) proposed an approximate method for solving the eigenvalue problem of interval matrix, which can be used in the modal superposition method, which is a popular method for dealing with structural dynamic problems. In this method, only the first few modes that make much greater contributions to the responses were taken into account. Thus modal superposition method is also a kind of model reduction technology (Wilson et al., 1982; Craig and Kurdila, 2006). Based on a linear fractional transformation, a model reduction method for a multidimensional uncertain system was proposed by Beck et al. (1996). A framework for Routh-Pade approximation of an interval system was proposed by Bandyopadhyay et al. (1994) in the frequency domain based on interval arithmetic. Dolgin and Zeheb (2005) proposed an approximate method based on Chebyshev for a complex model with interval parameters. In Wu and Jaramillo (2002), model reduction problems with multidimensional uncertainties were solved.

In this paper, based on the balanced truncation method, a model reduction strategy considering the uncertainties of the original system is proposed. In order to evaluate the quality of the reduced-order model, the interval perturbation method for solving the eigenvalue problem of the interval matrix is presented. The organization of the rest of the paper is as follows. In Section 2, based on balanced truncation methodology, a model reduction method for a system with interval uncertainties is proposed. In Section 3, the calculation method for the dynamic responses of the interval systems is proposed to compare the reduced-order model and the original model. Three numerical simulations are given to demonstrate the validity of the presented model reduction technique in Section 4 and the paper is concluded in Section5.

2. Balanced-based model reduction of uncertain systems with interval parameters

A generalization of the balanced-based model reduction strategy is proposed for complex systems with interval uncertainties. The closeness of the original model and the reduced-order model is quantified by the upper bounds of the $H \infty$ norm. Firstly, the decomposition of the system matrix with interval parameters is introduced.

2.1. The decomposition of the system matrix with interval parameters

For the stability analysis in engineering, the matrices of a system can be expressed as a linear function of the system parameters $b = (b i), i = 1, 2, \dots m$ , as follows:

S (b) = S 0 + \sum_{i = 1}^{m} b i S i

(1)

where

S i

are positive semi-definite (Qiu et al., 2005).

Due to the complexity of the actual engineering systems, the uncertainties are often involved in the parameters of the systems. The traditional approaches to dealing with uncertain problems are the probabilistic method and fuzzy-sets-based technique. Both methods require sufficient information to identify the probability density function or membership function. However, in practical engineering, it is expensive and time-consuming to achieve this information. To overcome the disadvantages of the probabilistic method and the fuzzy set method, interval analysis was developed based on interval mathematics. The advantages of the interval method are that only the bounds of the uncertain variables are required, instead of the probability density function or the membership function.

In this paper, considering that the actual situation where available information of the system parameters $b = (b i), i = 1, 2, \dots m$ is not sufficient to give the probability density function or membership function, we assume that each element of the system parameter vector b belongs to an interval number. Thus, the elements $s ij$ of the system matrix S are determined by b. Based on the interval natural extension principle (Moore, 1979), the following expression can be obtained:

S I = [\underline{S}, S ¯] = S 0 + \sum_{i = 1}^{m} b_{i}^{I} S i

(2)

where

b_{i}^{I} = [\underline{b} i, b ¯ i], i = 1, 2, \dots m

are the interval uncertain parameters of the system. This matrix decomposition is often used in engineering. For example, consider the following vibration equation with n degrees of freedom:

M x ·· (t) + P x \cdot (t) + Kx (t) = F (t)

(3)

where

M = (m ij), P = (c ij)

K = (k ij)

F (t) = (f i (t))

are the mass matrix, damping matrix, stiffness matrix and external load vector, respectively.

x (t) = (x i (t)), x \cdot (t) = (x \cdot i (t))

and

x ·· (t) = (x ·· i (t))

are the displacement vector, velocity vector, and acceleration vector, respectively. The vibration equation (3) can be transformed into a state-space equation as follows:

X \cdot = A \cdot X + B \cdot U Y = C \cdot X + D \cdot U

(4)

where

X = [x, x \cdot] T

U = [F (t)]

, and Y are the state vector, input and output vector, respectively, and

A = [0 I - M - 1 K - M - 1 P], B = [0 M - 1], C = [\underset{i - 1}{\underset{︸}{0 \dots 0}} 1 \underset{2 n - i}{\underset{︸}{0 \dots 0}}], D = [0]

(5)

The system matrices A and B are determined by the mass matrix, damping matrix and stiffness matrix of the structures, which can be decomposed into the following form, with respect to the structural parameter vector $b = (b i)$ :

M (b) = \sum_{i = 1}^{m} b i M i, P (b) = \sum_{i = 1}^{m} b i P i, K (b) = \sum_{i = 1}^{m} b i K i

(6)

The system matrices A and B as the functions of the $M, P$ and K can also be written as equation (1). For a survey of the non-negative decomposition of the matrices of structures, the reader may consult Qiu and Wang (2005); Zhang et al. (2007); Schmitendorf, (1988); and Petersen and Hollot (1986). Schmitendorf (1988) and Petersen and Hollot (1986) also report that an arbitrary $n \times n$ matrix can always be decomposed as the sum of rank-one matrices.

2.2. Balanced truncation of uncertain system

Consider the following state-space model:

X \cdot (t) = A \cdot X (t) + B \cdot U (t) Y \cdot (t) = C \cdot X (t) + D \cdot U (t)

(7)

where

A \in ℝ 2 n \times 2 n, B \in ℝ 2 n \times t, C \in ℝ 1 \times 2 n, D \in ℝ 1 \times t

. As mentioned in Section 1, because of various errors, uncertainties in the complex system are inevitable. Therefore, the system matrices are determined by the uncertain parameters, such as

b = (b i), i = 1, 2, \dots m

, and can be expressed as follows:

A = A (b) = (a ij (b)), B = B (b) = (b ij (b)) C = C (b) = (c ij (b)), D = D (b) = (d ij (b))

(8)

where

b = (b i)

is an m-dimensional vector. Thus, equation (7) can be rewritten as

X \cdot (t) = A (b) X (b, t) + B (b) U (t) Y (t) = C (b) X (b, t) + D (b) U (t)

(9)

Because there is not enough data to determine the probability density function or membership function, we assume that the uncertain parameter $b = (b i)$ belongs to an interval vector as follows:

b \in b I = [\underline{b}, b ¯] = (b_{i}^{I}), b i \in b_{i}^{I} = [\underline{b} i, b ¯ i], i = 1, 2, \dots m

(10)

where

\underline{b} = (\underline{b} i)

and

b ¯ = (b ¯ i)

are the lower and upper bounds of uncertain parameters

b = (b i)

, respectively. Based on interval arithmetic, the nominal value vector of the uncertain parameters vector can be defined as (Alefeld and Herzberger, 1984)

b c = (b_{i}^{c}) = m (b I) = \frac{(\underline{b} + b ¯)}{2}, b_{i}^{c} = m (b_{i}^{I}) = \frac{(\underline{b} i + b ¯ i)}{2}, i = 1, 2, \dots, m

(11)

and the uncertain radius vector of the uncertain parameters can be defined as

Δ b = (Δ b i) = rad (b I) = \frac{(b ¯ - \underline{b})}{2}, Δ b i = rad (b_{i}^{I}) = \frac{(b ¯ i - \underline{b} i)}{2}, i = 1, 2, \dots, m

(12)

Expanding the state matrix $A (b)$ around the deterministic part yields

A (b) = A (b c) + \sum_{j = 1}^{m} (\frac{\partial A (b c)}{\partial b j}) (b j - b_{j}^{c}) = A c + \sum_{j = 1}^{m} (\frac{\partial A (b c)}{\partial b j}) (b j - b_{j}^{c})

(13)

With the help of first-order Taylor series expansion, the system matrix can be approximately expressed as equation (13). When the uncertainties of the parameters are small, the above conclusion is appropriate. However, when the uncertainties are large, this method will no longer be suitable. In this case, the second-order terms in equation (13) should be considered. That is, the method proposed in this paper is limited to the case where the uncertainties of the parameters are small. Based on the interval natural extension principle, equation (13) can be written as follows:

A I = A c + Δ A \cdot e

(14)

where

Δ A = \sum_{j = 1}^{m} | (\frac{\partial A (b c)}{\partial b j}) | Δ b j

and

e = [- 1, 1]

Similar expressions exist for the other system matrices

B I = B c + Δ B \cdot e = B c + Δ B I C I = C c + Δ C \cdot e = C c + Δ C I D I = D c + Δ D \cdot e = D c + Δ D I

(15)

where

Δ B = \sum_{j = 1}^{m} | (\frac{\partial B (b c)}{\partial b j}) | \cdot Δ b j, Δ C = \sum_{j = 1}^{m} | (\frac{\partial C (b c)}{\partial b j}) | \cdot Δ b j, Δ D = \sum_{j = 1}^{m} | (\frac{\partial D (b c)}{\partial b j}) | \cdot Δ b j

Assume $Σ c = (A c B c C c D c)$ represents the deterministic part of the system or model, which is the nominal value of the uncertain model. $δ A$ , $δ B$ , $δ C$ and $δ D$ represent the uncertain part of the model. They are uncertain, but the bounds of the uncertainties are deterministic as follows:

δ A \in Δ A I, δ B \in Δ B I, δ C \in Δ C I, δ D \in Δ D I

(16)

Then, the system with uncertainties can be rewritten as follows:

Σ = (A B C D) = (A c B c C c D c) + (δ A δ B δ C δ D)

(17)

Let $L c$ and $L o$ denote the Gramians of the system shown as equation (17). The infinite reachability and infinite observability Gramians are

L c = \int_{0}^{\infty} e A τ BB T e A T τ d τ L o = \int_{0}^{\infty} e A T τ C T C e A τ d τ

(18)

which satisfy the following Lyapunov equations:

AL c + L c A T + BB T = 0 A T L o + L o A + C T C = 0

(19)

Readers are encouraged to consult a basic text on the concept of infinite Gramians such as Approximation of Large-scale Dynamical Systems (Antoulas, 2005) for a more detailed understanding. Because $δ A$ , $δ B$ , $δ C$ and $δ D$ are small, the Gramians matrices of the uncertain system can be written in the following form based on the principle of perturbation:

L c = L_{c}^{c} + δ L c, L o = L_{o}^{c} + δ L o

(20)

Substituting equation (20) into equation (19) yields

(A c + δ A) (L_{c}^{c} + δ L c) + (L_{c}^{c} + δ L c) (A c + δ A) T + (B c + δ B) (B c + δ B) T = 0 (A c + δ A) T (L_{o}^{c} + δ L o) + (L_{o}^{c} + δ L o) (A c + δ A) + (C c + δ C) T (C c + δ C) = 0

(21)

Expanding equation (21), and neglecting the high-order terms yields

A c δ L c + δ {AL}_{c}^{c} + L_{c}^{c} δ A T + δ L c (A c) T + B c δ B T + δ B (B c) T = 0 (A c) T δ L o + δ A T L_{o}^{c} + L_{o}^{c} δ A + δ L o A c + (C c) T δ C + δ C T C c = 0

(22)

and

A c L_{c}^{c} + L_{c}^{c} (A c) T + B c (B c) T = 0 (A c) T L_{o}^{c} + L_{o}^{c} A c + (C c) T C c = 0

(23)

Rearrangement of equation (23) yields the equations:

A c δ L c + δ L c (A c) T + (B c δ B T + δ B (B c) T + δ {AL}_{c}^{c} + L_{c}^{c} δ A T) = 0 (A c) T δ L o + δ L o A c + ((C c) T δ C + δ C T C c + δ A T L_{o}^{c} + L_{o}^{c} δ A) = 0

(24)

which are also Lyapunov equations.

Noticeably, the nominal Gramians $L_{c}^{c}$ and $L_{o}^{c}$ of the system can be obtained by solving equation (23) directly. Simultaneously, solving the Lyapunov equations equation (24) yields the $δ L c$ and $δ L o$ which are the perturbations of the Gramians. Based on the interval natural extension principle, the interval Gramians $L_{c}^{I}$ and $L_{o}^{I}$ can be obtained as

L c \in L_{c}^{I} = L_{c}^{c} + Δ L_{c}^{I} = [L_{c}^{c} - Δ L c, L_{c}^{c} + Δ L c] L o \in L_{o}^{I} = L_{o}^{c} + Δ L_{o}^{I} = [L_{o}^{c} - Δ L o, L_{o}^{c} + Δ L o]

(25)

Finding a transformation matrix T is the key to make the system balanced. Based on the Cholesky decomposition and the eigenvalue decomposition theorem of matrix, the nominal Gramians can be factored into

L_{c}^{c} = VV T, V T L_{o}^{c} V = G (Σ c) 2 G T

(26)

where

Σ c = diag (σ_{1}^{c}, σ_{2}^{c}, \dots, σ_{2 n}^{c})

and

σ_{k}^{c} (Σ c) = \sqrt{λ k (L_{c}^{c} L_{o}^{c})}, k = 1, 2, \dots, 2 n

λ k (\cdot)

denotes the kth eigenvalue of

(\cdot)

. Then, the transformation matrix T can be achieved as follows:

T = (Σ c) 1 / 2 G T V - 1

(27)

After transformation, the system matrices of the balanced system with interval uncertainties can be obtained as follows:

A ~ I = TA I T - 1, B ~ I = TB I, C ~ I = C I T - 1, D ~ I = D I

(28)

Thus, the balanced system $Σ ~ I$ with interval uncertainties is

Σ ~ I = (A ~ B ~ C ~ D ~) \in (A ~ I B ~ I C ~ I D ~ I) = ([A ~ c - Δ A ~, A ~ c + Δ A ~] [B ~ c - Δ B ~, B ~ c + Δ B ~] [C ~ c - Δ C ~, C ~ c + Δ C ~] [D ~ c - Δ D ~, D ~ c + Δ D ~])

(29)

where

Δ A ~ = rad (A ~ I), Δ B ~ = rad (B ~ I), Δ C ~ = rad (C ~ I), Δ D ~ = rad (D ~ I)

In order to achieve the objective of reducing the order, some state variables need to be eliminated. Considering the characteristics of the balanced model, it is known that the state variables are sorted according to the influence on the input-output behavior. Thus, portioning the state vector $X ~$ of the balanced model with interval uncertainties as $X ~ 1 = [X ~ 1, X ~ 2, \dots, X ~ k] T \in ℝ k$ containing the state components with the largest influence on the input-output behaviors and $X ~ 2 = [X ~ k + 1, X ~ k + 2, \dots, X ~ 2 n] T \in ℝ 2 n - k$ containing the state components with the smallest influence on the input-output behaviors yields

Σ ~ I = [{Σ ~}_{1}^{I} 00 {Σ ~}_{2}^{I}], A ~ I = [{A ~}_{11}^{I} {A ~}_{12}^{I} {A ~}_{21}^{I} {A ~}_{22}^{I}], B ~ I = [{B ~}_{1}^{I} {B ~}_{2}^{I}], C ~ I = [{C ~}_{1}^{I}, {C ~}_{2}^{I}]

(30)

Then, a reduced-order model with interval uncertainties of order k can be obtained as follows

Σ \land I = [{A ~}_{11}^{I} {B ~}_{1}^{I} {C ~}_{1}^{I} D ~ I] = [[A ~ 11 - Δ A ~ 11, A ~ 11 + Δ A ~ 11] [B ~ 1 - Δ B ~ 1, B ~ 1 + Δ B ~ 1] [C ~ 1 - Δ C ~ 1, C ~ 1 + Δ C ~ 1] [D ~ c - Δ D ~, D ~ c + Δ D ~]]

(31)

As can be seen from the above equation, the system matrices of the reduced-order model are interval matrices rather than deterministic matrices. The reduced-order model with interval uncertainties can be used to analyze the robustness of the original system. However, the reduced-order model obtained by conventional model reduction technologies cannot meet this requirement.

2.3. The error bound of reduced-order model with interval uncertainties

Antoulas (2005) has demonstrated that the quality of the reduced-order model obtained by balanced truncation can be determined by the $H \infty$ error bound, which depends on the eliminated Hankel singular values

‖ Σ - Σ \land ‖ H \infty \leq 2 (σ k + 1 + \dots + σ 2 n)

(32)

Because the model reduction method proposed in the previous section for a system with interval uncertainties is based on the balanced truncation, equation (32) can also be used to measure the quality of the reduced-order model, but the neglected Hankel singular values are interval numbers. Thus, in order to assess the quality of the reduced-order model, we should first calculate the intervals of the discarded Hankel singular values. In the previous section, the Gramians $L_{c}^{I}$ and $L_{o}^{I}$ of uncertain system with interval parameters have been obtained. For the sake of convenience, the product of $L_{c}^{I}$ and $L_{o}^{I}$ is denoted by $Λ I$

Λ I = L_{c}^{I} \cdot L_{o}^{I} = Λ c + Δ Λ I = [Λ c - Δ Λ, Λ c + Δ Λ]

(33)

Using the interval algorithm, $Λ I$ can be obtained. To get the Hankel singular values of the system $Σ I$ with interval uncertainties, the following interval perturbation method for eigenvalue calculation of interval matrix is introduced.

Assume $λ i, i = 1, 2, \dots, 2 n$ are the eigenvalues of the nominal matrix $Λ c$ of an interval matrix $Λ c + Δ Λ I$ and $v_{i}^{c}$ are the corresponding eigenvectors. The eigenvalues of the matrix with interval uncertainties $Λ c + Δ Λ I$ are (Qiu et al., 2001)

λ_{i}^{I} = [\underline{λ} i, λ ¯ i], i = 1, 2, \dots, 2 n

(34)

where

\underline{λ} i = λ_{i}^{c} - | v_{i}^{c} | T \cdot Δ Λ \cdot | v_{i}^{c} | λ ¯ i = λ_{i}^{c} + | v_{i}^{c} | T \cdot Δ Λ \cdot | v_{i}^{c} |

(35)

v_{i}^{c}

and

λ_{i}^{c}

satisfy the following equations:

Λ c v_{i}^{c} = λ_{i}^{c} v_{i}^{c}, i = 1, 2, \dots, 2 n

(36)

Qiu et al. had proved that the interval perturbation algorithm requires minimal computational effort compared with Deif’s solution (Deif, 1991). Then, the interval bounds of the Hankel singular values will be

σ_{i}^{I} (Σ I) = \sqrt{λ_{i}^{I}} = [\sqrt{\underline{λ} i}, \sqrt{λ ¯ i}]

(37)

where

\underline{λ} i

and

λ ¯ i

can be obtained by equation (35). To facilitate our problem, we can utilize the interval perturbation algorithm to get the interval Hankel singular values of the complex model with interval parameters. Then, a novel method to get the

H \infty

error bound for model reduction with interval parameters is proposed. Considering equations (32) and (37), we know that the neglected Hankel singular values

σ i

i = k + 1, \dots, 2 n

are in

2 n - k

intervals, which means

σ i \in [\underline{σ} i, σ ¯ i]

. Then, using the interval algorithm equation (37), we can obtain the following equation:

σ_{k + 1}^{I} + \dots + σ_{2 n}^{I} = [\underline{σ} k + 1, σ ¯ k + 1] + \dots + [\underline{σ} 2 n, σ ¯ 2 n] = [\underline{σ} k + 1 + \dots + \underline{σ} 2 n, σ ¯ k + 1 + \dots + σ ¯ 2 n]

(38)

which means that

(σ k + 1 + \dots + σ 2 n) \in [\underline{σ} k + 1 + \dots + \underline{σ} 2 n, σ ¯ k + 1 + \dots + σ ¯ 2 n]

Therefore, the error bound of the reduced-order model with interval uncertainties is

‖ Σ I - Σ \land I ‖ H \infty \leq 2 (σ k + 1 + \dots + σ 2 n) \leq 2 σ ¯ k + 1 + \dots + 2 σ ¯ 2 n

(39)

The error bound confirms that a good reduction will be taken when the upper bounds of the Hankel singular values which are neglected $σ ¯ k + 1, \dots, σ ¯ n$ are small. The balanced-based model reduction method of uncertain systems with interval parameters is different from deterministic reduction methods where the error bound depends on the Hankel singular values that are neglected.

Based on balanced truncation, a model reduction method for complex systems with interval uncertainties is proposed. The reduced-order model obtained by this method can hold the stability of the original model and the error bound of the reduced-order model can be obtained by using the upper bounds of the neglected Hankel singular values. In order to compare the original model and the reduced-order model, dynamic responses of the reduced-order model and the original model are selected as an intuitive criterion. By comparing the upper and lower bounds of the dynamic responses of the original and reduced-order model, we can evaluate the quality of the proposed method in this paper. Thus, an algorithm to get the responses of the interval system will be introduced in the next section.

3. Dynamic response of system with interval uncertainties

Consider the interval state-space model as follows (Chen, 2004):

X \cdot = A (b) X + B (b) U

(40)

where

b \in b I

. Due to the uncertainty of the parameter b, the dynamic responses of the system determined by uncertain parameter b are also uncertain. Because b is bounded in an interval vector, there must be a response set containing all possible responses

Γ = {X (b, t) : X \cdot = A (b) X + B (b) U, b \in b I}

(41)

In general, finding out the complex response set is very difficult and not necessary. In order to reduce the amount of computation, we just need to find an interval vector containing the response set. That is, we need to find out the upper and lower bounds of the dynamic response.

X I (b, t) = [\underline{X} (b, t), X ¯ (b, t)] = (X_{i}^{I} (b, t))

(42)

where

\underline{X} (b, t) = (\underline{X} i (b, t))

and

X ¯ (b, t) = (X ¯ i (b, t))

, and

\underline{X} (b, t) = min {X (b, t) : X \cdot = A (b) X + B (b) U, b \in b I} X ¯ (b, t) = max {X (b, t) : X \cdot = A (b) X + B (b) U, b \in b I}

(43)

Assuming that $δ A$ and $δ B$ represent small perturbations from $A c$ and $B c$ , the state-space equation can be rewritten as

X \cdot = (A c + δ A) X + (B c + δ B) U

(44)

Assume that the solution X is of the form

X = X c + δ X

(45)

Substituting equation (45) into equation (44) yields

X \cdot c + δ X \cdot = (A c + δ A) (X c + δ X) + (B c + δ B) U

(46)

Expanding equation (46) and neglecting the higher-order terms yields

X \cdot c = A c X c + B c U δ X \cdot = A c δ X + δ AX c + δ BU

(47)

where

δ X

can be seen as a function of the elements in

δ A

and

δ B

. Based on the interval extension principle, the second equation in equation (47) can be written as

Δ X \cdot = A c \cdot Δ X + Δ A \cdot X c + Δ BU

(48)

where

Δ A

and

Δ B

are the radius of interval matrices

A I

and

B I

, respectively. This method is limited to the case where the uncertainties of the parameters are small, because only the first-order perturbation is used in equation (47). If the uncertainties of the parameters of systems are large, the second-order perturbation should be considered. Solving the first equation in equation (47) and (48) yields

X I = X c + Δ X I

(49)

where

Δ X I = Δ X \cdot [- 1, 1]

Then, the upper bounds and the lower bounds of the solution to equation (40) can be obtained:

X ¯ = X c + Δ X \underline{X} = X c - Δ X

(50)

Thus, the upper and lower bounds of dynamic responses of the reduced-order model and the original model are obtained. By comparing the bounds of the responses, we can quantitatively analyze the closeness of the reduced-order model and the original model.

4. Numerical examples

In this section, the proposed method is applied to three numerical simulations and the results are discussed later.

Example 1

Consider a fourth-order single-input-single-output linear time-invariant system with a state-space realization $(A, B, C, D)$ , where $A \in A I, B \in B I, C \in C I, D \in D I$ . $Δ A, Δ B, Δ C, Δ D$ represent the bounds of the system uncertainties, namely, what is unknown or poorly understood about the system.

A c = [- 10 - 4.375 - 3.125 - 1.5 800002000010], B c = [2000] C c = [0.5 0.4375 0.75 0.75], D c = 0 Δ A = | A c | \times 0.05 Δ B = | B c | \times 0.05 Δ C = | C c | \times 0.05

The reachability Gramian $L_{c}^{I}$ and the observability Gramian $L_{o}^{I}$ of this uncertain system with interval parameters can be calculated as follows:

L_{c}^{c} = [0.2397 0 - 0.1270 00 0.5079 0 - 0.2032 - 0.127 0 0.4063 00 - 0.2032 0 0.5079], Δ L c = [0.0407 0 0.0343 00 0.1371 0 0.0548 0.0343 0 0.096 00 0.0548 0 0.1286]

L_{o}^{c} = [0.1107 0.1228 0.1982 0.1875 0.1228 0.1371 0.2207 0.2083 0.1982 0.2207 0.3572 0.3382 0.1875 0.2083 0.3382 0.3208], Δ L o = [0.0254 0.0172 0.0267 0.0253 0.0141 0.0027 0.0036 0.0036 0.0241 0.0063 0.008 0.0074 0.0328 0.0178 0.0254 0.0233]

Using equation (28) and interval mathematics, the balanced system with interval uncertainties can be obtained

A ~ c = [- 0.82 0.3146 - 0.7301 0.0766 - 0.3146 - 0.4480 3.7879 - 0.2365 - 0.7301 - 3.7879 - 7.1090 1.3934 0.0766 0.2365 1.3934 - 1.6231], B ~ c = [0.9216 0.1663 0.4201 - 0.0431] Δ A ~ = [0.1196 1.0001 2.3658 5.3110 0.0187 0.1572 0.5654 0.8926 0.0377 0.2800 0.6178 1.5320 0.0050 0.0426 0.1024 0.2221], Δ B ~ = [0.0461 0.0083 0.0210 0.0022] C ~ c = [0.9216 - 0.1663 0.4201 - 0.0431] Δ C ~ = [0.0461 0.2065 0.2718 1.3348]

Seen from the Table 1, the first two eigenvalues are much larger than the other eigenvalues. Thus, the third and fourth state variables can be deleted according to the balanced truncation approximation method, and the reduced-order model (with

X \land \in ℝ 2

an approximation of

X ~ 1 \in ℝ 2

) is given by the state-space realization:

{̂ \cdot X = {A ~}_{11}^{I} X ̂ + {B ~}_{1}^{I} U Y ̂ = {C ~}_{1}^{I} X ̂ + {D ~}_{1}^{I} U

where

{A ~}_{11}^{I} = [[- 0.9396, - 0.7003] [- 0.6855, 1.3148] [- 0.3334, - 0.2959] [- 0.6052, - 0.2907]], D ~ I = 0 {B ~}_{1}^{I} = [[0.8754, 0.9677] [0.1579, 0.1746]], {C ~}_{1}^{I} = [[0.8754, 0.9677] [- 0.3728, 0.0402]]

Table 1.

Interval eigenvalues of the original system with interval parameters.

$λ_{i}^{I}$	$λ c$	$Δ λ i$
$λ_{1}^{I}$	0.2682	0.1520
$λ_{2}^{I}$	0.95e-3	0.62e-4
$λ_{3}^{I}$	0.15e-3	0.30e-3
$λ_{4}^{I}$	0.33e-6	0.34e-3

In order to compare the complex model and the reduced-order model, the input $u = 2 \sin (2 π t)$ is given. Based on the method of calculating the dynamic response for an interval system presented in the previous section, the interval dynamic response of the complex model and the reduced-order model are obtained.

As can be seen from Figure 1, the response of the reduction-order model (order 2) is a good reproduction of the original system (order 4) response. Ignoring the uncertainties, we can also get the above result using the conventional balanced truncation method. This shows that the proposed method can be degraded to the traditional method of balanced truncation.

Figure 1.

The nominal output of the original system and reduced-order system (reduction from order 4 to order 2).

Figure 2 gives the details of responses of the complex model and the reduced model. Qualitatively, we can conclude that the reduced-order model closely tracks the pattern of the original model. As shown in Figure 2, the proposed method in this paper can get both the nominal response of the original model and the ranges of response of the complex model.

Figure 2.

The bounds of responses of the original system and reduced-order system (reduction from order 4 to order 2).

Example 2

Considering the following mass-spring-damper vibration system shown in Figure 3, we assume that the mass, stiffness and damping are all uncertain quantities. Assume that there is one sinusoidal excitation $F = 20 N$ , $ω = 2 π f$ , $f = 10 Hz$ acting on the node m₃ in the horizontal direction, and the initial conditions are $x = 0, x \cdot = 0$ . Due to errors, the mass, stiffness and damping exhibit some uncertainties, and are considered to be interval parameters. Their intervals are taken as $m_{1}^{I} = [1.15, 1.25] kg$ , $m_{2}^{I} = [0.95, 1.05] kg$ , $m_{3}^{I} = [1.5, 1.7] kg$ , $k_{1}^{I} = [349.7, 350.3] Nm - 1$ , $k_{2}^{I} = [299.7, 300.3] Nm - 1$ , $k_{3}^{I} = [319.7, 320.3] Nm - 1$ , $c_{1}^{I} = [3.96, 4.04] Nsm - 1$ , $c_{2}^{I} = [2.97, 3.03] Nsm - 1$ , $c_{3}^{I} = [3.17, 3.23] Nsm - 1$ . The output is the horizontal displacement of m₃.

Figure 3.

The three-degrees-of-freedom system (mass-spring-damper vibration system).

Applying the Lagrange equation, the mathematical model of the system in Figure 3 is

[m 1 000 m 2 000 m 3] [x ·· 1 x ·· 2 x ·· 3] + [c 1 + c 2 - c 2 0 - c 2 c 2 + c 3 - c 3 0 - c 3 c 3] [x \cdot 1 x \cdot 2 x \cdot 3] + [k 1 + k 2 - k 2 0 - k 2 k 2 + k 3 - k 3 0 - k 3 k 3] [x 1 x 2 x 3] = [f 1 f 2 f 3]

Using equations (5) and (6), we obtain the means and radius of the system matrices. Then, the proposed model reduction method is used to get the reduced-order model. Table 2 shows that the first two eigenvalues are much larger than the other four eigenvalues. Thus, the third ∼ sixth state variables can be deleted, and the reduced-order model will be derived based on the balanced truncation.

Table 2.

Interval eigenvalues of the original system with interval parameters.

$λ_{i}^{I}$	$λ c$	$Δ λ i$
$λ_{1}^{I}$	3.91e-3	1.14e-1
$λ_{2}^{I}$	3.38e-3	8.63e-2
$λ_{3}^{I}$	9.66e-7	2.05e-4
$λ_{4}^{I}$	6.11e-7	1.68e-4
$λ_{5}^{I}$	4.07e-9	4.48e-5
$λ_{6}^{I}$	3.15e-9	4.55e-6

The nominal horizontal displacement of m₃ of the original model and the reduced-order model are shown in Figure 4. As can be seen, even without considering the uncertainties, the reduced-order model can give a high-quality estimate of the original model. A comparison between the complex model and the reduced-order model with interval uncertainties is given in Figure 5. It shows that the reduced model can not only track the nominal response of the complex model, but also the response ranges of the complex model with a small deal of computational efforts. The reduced-order model with uncertainties will be appropriate when the uncertainties of the original model need to be considered, such as stability analysis, robustness analysis or reliability analysis.

Figure 4.

The nominal displacement of m₃ of the original system and reduced system (reduction from order 6 to order 2).

Figure 5.

The bounds of the displacement of the original system and reduced system (reduction from order 6 to order 2).

Example 3

In order to demonstrate the applicability of the model reduction method proposed in this paper for the engineering system with uncertainties, we apply this model reduction method to a two-dimensional truss structure shown in Figure 6. The material properties are as follows: Young’s modulus $E = 210 GPa$ and the density of the material is $ρ = 7300 kg / m 3$ . Geometry and boundary conditions of the truss are shown in Figure 6. An impulse excitation act on the truss vertically. The sectional areas of the rod members are uncertain variables $A I = [A c - η A c, A c + η A c]$ and all trusses have the same cross section $A c = 1 \times 10 - 4 m 2$ and $η = 0.01$ . Proportional damping $C = α M + β K$ is adopted with $α = 0.001$ and $β = 0.003$ . The output is the vertical displacement of Node 2.

Figure 6.

A seven-bar two-dimensional truss.

Table 3 shows that the first two eigenvalues are much larger than the other 12 eigenvalues. The reduced-order model is derived by deleting the third ∼14th state variables. The nominal values of the dynamic responses of the original model and the reduced-order model of the seven-bar two-dimensional truss are shown in Figure 7. The bounds of the displacement of Node 2 are shown in Figure 8. Qualitatively, we conclude that the reduced-order model reproduces the original model well, and the proposed model reduction method in this paper is also feasible for practical engineering problems.

Figure 7.

The nominal displacement of Node 2 of the original system and reduced system (reduction from order 14 to order 2).

Figure 8.

The bounds of the displacement of Node 2 (reduction from order 14 to order 2).

Table 3.

Interval eigenvalues of the original system with interval parameters.

$λ_{i}^{I}$	$λ c$	$Δ λ i$
$λ_{1}^{I}$	1.69e-10	2.84e-12
$λ_{2}^{I}$	5.40e-11	1.02e-12
$λ_{3}^{I}$	3.82e-14	1.61e-14
$λ_{4}^{I}$	9.99e-15	6.58e-16
$λ_{5}^{I}$	8.93e-16	1.40e-14
$λ_{6}^{I}$	5.70e-16	2.93e-13
$λ_{7}^{I}$	3.56e-17	8.76e-15
$λ_{8}^{I}$	2.88e-18	1.15e-15
$λ_{9}^{I}$	4.87e-22	6.81e-15
$λ_{10}^{I}$	1.83e-22	2.08e-15
$λ_{11}^{I}$	1.00e-30	1.01e-13
$λ_{12}^{I}$	1.00e-31	3.27e-14
$λ_{13}^{I}$	1.00e-32	3.88e-14
$λ_{14}^{I}$	1.00e-33	2.56e-14

5. Conclusion

Based on interval arithmetic, a model reduction method for complex models with interval uncertainties was proposed. By using the interval perturbation method, the upper bounds of the Hankel singular values of the uncertain system, which can be used to assess the quality of the reduced-order model, were derived. In order to give an intuitive evaluation of the quality of the reduced-order model, the matrix perturbation method was adopted to calculate the bounds of the dynamic response of the original model and the reduced-order model. The difference between the model reduction method proposed in this paper and conventional methods is that the coefficients of the obtained reduced-order model are interval numbers instead of fixed numbers. This method is useful when the uncertainties of original systems must be considered, such as robustness analysis, reliability analysis and stability analysis. The results of numerical examples showed that the proposed model reduction method could lead to a significant reduction in the computation time and guarantee precision. As a future research direction, other model reduction methods would be extended to complex models with uncertainties.

Footnotes

Acknowledgements

The authors would like to thank their funding bodies.

Funding

This work was supported by the National Nature Science Foundation of China (grant number 11372025), Defense Industrial Technology Development Program (grant numbers B2120110011 and JCKY2013601B), and the Aeronautical Science Foundation of China (grant number A2012ZA51010).

References

Alefeld

Herzberger

(1984) Introduction to Interval Computation, New York: Academic Press.

Antoulas

(2005) Approximation of Large-scale Dynamical Systems, Philadelphia: SIAM.

Antoulas

Sorensen

Gugercin

(2001) A survey of model reduction methods for large-scale systems. Contemporary Mathematics 280: 193–220.

Bandyopadhyay

Ismail

Gorez

(1994) Routh-Padé approximation for interval systems. IEEE Transactions on Automatic Control 39(12): 2454–2456.

Beck

Doyle

Glover

(1996) Model reduction of multidimensional and uncertain systems. IEEE Transactions on Automatic Control 41(10): 1466–1477.

Bellizzi

Sampaio

(2009) Smooth Karhunen–Loève decomposition to analyze randomly vibrating systems. Journal of Sound and Vibration 325(3): 491–498.

Benner P, Quintana-Orti E and Quintana-Ortí G (2000) Singular perturbation approximation of large, dense linear systems. IEEE International Symposium on Computer-Aided Control System Design, 2000. CACSD 2000, 25–27 September 2000, Anchorage, Alaska, pp. 255–260.

Besselink

Tabak

Lutowska

(2013) A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control. Journal of Sound and Vibration 332(19): 4403–4422.

Chen

Guo

Chen

(2004) A method for estimating upper and lower bounds of eigenvalues of closed-loop systems with uncertain parameters. Journal of Sound and Vibration 276(3): 527–539.

10.

Craig

Kurdila

(2006) Fundamentals of Structural Dynamics, New York: Wiley.

11.

Deif

(1991) The interval eigenvalue problem. Journal of Applied Mathematics and Mechanics 71(1): 61–64.

12.

Dolgin

Zeheb

(2005) Model reduction of uncertain systems retaining the uncertainty structure. Systems & Control Letters 54(8): 771–779.

13.

Friswell

Garvey

Penny

JET

(1995) Model reduction using dynamic and iterated IRS techniques. Journal of Sound and Vibration 186(2): 311–323.

14.

Friswell

Garvey

Penny

JET

(1998) The convergence of the iterated IRS method. Journal of Sound and Vibration 211(1): 123–132.

15.

Gallivan

Vandendorpe

Van Dooren

(2004) Sylvester equations and projection-based model reduction. Journal of Computational and Applied Mathematics 162(1): 213–229.

16.

Ghommem

Nayfeh

Choura

(2012) Model reduction and analysis of a vibrating beam microgyroscope. Journal of Vibration and Control 19((8): 1240–1249.

17.

Glover

(1984) All optimal Hankel-norm approximations of linear multivariable systems and their L,∞-error bounds. International Journal of Control 39(6): 1115–1193.

18.

Gugercin

(2008) An iterative SVD-Krylov based method for model reduction of large-scale dynamical systems. Linear Algebra and its Applications 428: 1964–1986.

19.

Gugercin

Sorensen

Antoulas

(2003) A modified low-rank Smith method for large-scale Lyapunov equations. Numerical Algorithms 32(1): 27–55.

20.

Guyan

(1965) Reduction of stiffness and mass matrices. AIAA Journal 3(2): 380.

21.

Holmes

Lumley

Berkooz

(1998) Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge: Cambridge University Press.

22.

Khattak

Garvey

Popov

(2010) Proper orthogonal decomposition of the dynamics in bolted joints. Journal of Sound and Vibration 329(9): 1480–1498.

23.

Kokotovic

O’Malley

Sannuti

(1976) Singular perturbations and order reduction in control theory – an overview. Automatica 12(2): 123–132.

24.

Komarasamy

Albhonso

Gurusamy

(2012) Order reduction of linear systems with an improved pole clustering. Journal of Vibration and Control 18(12): 1876–1885.

25.

Liu

Anderson

(1989) Singular perturbation approximation of balanced systems. International Journal of Control 50(4): 1379–1405.

26.

Moore

(1981) Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Transactions on Automatic Control 26(1): 17–32.

27.

Moore

(1979) Methods and Applications of Interval Analysis, Philadelphia: SIAM.

28.

Mullis

Roberts

(1976) Synthesis of minimum roundoff noise fixed point digital filters, circuits and systems. IEEE Transactions on Automatic Control 23(9): 551–562.

29.

O’Callahan J (1989) A procedure for an improved reduced system (IRS) model. In: Proceedings of the 7th International Modal Analysis Conference, Las Vegas, Nevada, February, pp. 17–21.

30.

Petersen

Hollot

(1986) A Riccati equation approach to the stabilization of uncertain linear systems. Automatica 22(4): 397–411.

31.

Qiu

Wang

(2005) Parameter perturbation method for dynamic responses of structures with uncertain-but-bounded parameters based on interval analysis. International Journal of Solids and Structures 42(18): 4958–4970.

32.

Qiu

Chen

Elishakoff

(1996) Bounds of eigenvalues for structures with an interval description of uncertain-but-non-random parameters. Chaos, Solitons & Fractals 7(3): 425–434.

33.

Qiu

Müller

Frommer

(2001) An approximate method for the standard interval eigenvalue problem of real non-symmetric interval matrices. Communications in Numerical Methods in Engineering 17(4): 239–251.

34.

Qiu

Wang

Friswell

(2005) Eigenvalue bounds of structures with uncertain-but-bounded parameters. Journal of Sound and Vibration 282(1): 297–312.

35.

Rayleigh

BJWS

(1896) The Theory of Sound, London: Macmillan.

36.

Ritto

Buezas

Sampaio

(2011) A new measure of efficiency for model reduction: application to a vibroimpact system. Journal of Sound and Vibration 330(9): 1977–1984.

37.

Schmitendorf

(1988) Designing stabilizing controllers for uncertain systems using the Riccati equation approach. IEEE Transactions on Automatic Control 33(4): 376–379.

38.

Wilson

Yuan

Dickens

(1982) Dynamic analysis by direct superposition of Ritz vectors. Earthquake Engineering & Structural Dynamics 10(6): 813–821.

39.

Wu F and Jaramillo J (2002) Numerical algorithm for model reduction of polytopic uncertain linear systems. In: Proceedings of the 2002, IEEE, American Control Conference, 8–10 May 2002, Anchorage, Alaska, pp. 1748–1752.

40.

Xia

Lin

(2004a) A new iterative order reduction (IOR) method for eigensolutions of large structures. International Journal for Numerical Methods in Engineering 59(1): 153–172.

41.

Xia

Lin

(2004b) Improvement on the iterated IRS method for structural eigensolutions. Journal of Sound and Vibration 270(4): 713–727.

42.

Zhang

Ding

Chen

(2007) Interval finite element method for dynamic response of closed-loop system with uncertain parameters. International Journal for Numerical Methods in Engineering 70(5): 543–562.