Low-dimensional approximations and error estimates for systems of linear random ODEs

J. vom Scheidt  H.-J. Starkloff  R. Wunderlich

Abstract

The paper considers large-scale systems of linear ODEs with stationary random excitation terms. Approximations of the correlation and spectral density functions of the stationary solutions are found by means of model reduction techniques. To assess the approximation accuracy suitable error criterions are introduced and error estimates are given.

1 Introduction

We consider systems of linear ordinary differential equations (ODEs) with random parameters for a vector function z = (z_{1}, …, z_{n})^{τ}, n ∈ N,

ż = Az + Bh.  (1)

The random parameter of this equation is the inhomogeneous term containing the excitation function h which is assumed to be a stationary random function with continuous paths defined on R × Ω with values in C^{r}, r ∈ N. Here, R and C denote the sets of real and complex numbers, respectively. (Ω, A, P) is a probability space, where A denotes the σ-algebra of subsets of Ω on which is defined a probability measure P. Further, A is a complex n × n-matrix and B a complex n × r-matrix. The main tasks are the determination of so-called stationary solutions z and the computation of their distribution characteristics, in particular second-order moment functions. These questions arise in the investigation of the long-time behaviour of the response of systems described by Equation (1) to permanently acting random excitations.

The results can be directly applied to discrete vibration systems with damping and an external excitation which are usually described by systems of n real second-order ODEs which can be transformed into systems of 2n first-order ODEs (see Preumont [9], Newland [8], Wirsching, Paez, Ortiz [16], Soong, Grigoriu [14])

A random function z(t, ω) satisfying (1) is called stationary solution of (1) if z and h are stationarily related. Assuming a stable matrix A, i.e. all eigenvalues of A possess negative real parts, then the linear first-order system (1) possesses the unique stationary solution (see Bunke [3], p. 45ff., Arnold, Wihstutz [1])

z(t, ω) = ∫_{−∞}^{t} e^{A(t−s)}Bh(s, ω)ds = ∫_{0}^{∞}e^{As}Bh(t − s, ω)ds.  (2)

The representation (2) can be used for the probabilistic analysis of the stationary solution z(t, ω) and the computation of its moment functions. For the first-order moments it follows m_{z} := E {z(t)} = 0 since h is assumed to be centered. Here, E {.} denotes the expectation with respect to the probability measure P.

For the computation of second-order moment functions as correlation and spectral density function several procedures are known in the literature (see Newland [8], Preumont [9], Soong, Grigoriu [14], vom Scheidt, Fellenberg, Wöhrl [11]). A first group is based on time domain analysis. Here, averages of the representation (2) are used. For the correlation function of z it follows

R_{zz}(τ) := E {(z(t) − m_{z})(z(t + τ) − m_{z})^{∗}}
= ∫_{0}^{∞}∫_{0}^{∞}e^{As_{1}}BE {h(t − s_{1})h^{∗}(t + τ − s_{2})} B^{∗}e^{A^{∗}s_{2}}ds_{1}ds_{2}

= ∫_{0}^{∞}∫_{0}^{∞}e^{As_{1}}BR_{hh}(τ + s_{1} − s_{2})B^{∗}e^{A^{∗}s_{2}}ds_{1}ds_{2}.

Other procedures are based on frequency domain analysis. Here, the spectral density of z is related to the spectral density of h. Assume that the excitation function h possesses a spectral density denoted by S_{hh}(α), α ∈ R. Then, for the spectral density of z it holds (see Soong, Grigoriu [14], p. 194 ff., vom Scheidt, Fellenberg, Wöhrl [11], p. 102 ff.)

S_{zz}(α) := \frac{1}{2π} ∫_{−∞}^{∞} e^{−iατ} R_{zz}(τ) dτ = H(α)S_{hh}(α)H^{∗}(α), α ∈ R,
where H(α) = −(iαI + A)^{−1}B.

The matrix function H(α) is known as the frequency response function or transition function.

The aim of this paper is the approximation of moment functions of stationary solutions of large-scale systems of differential equations by means of model (or dimension) reduction techniques. If the dimension n is very large, then procedures for the computation of moment functions may become intractable because of enormous computational problems. This is the typical case if partial differential equations (PDEs) are approximated by systems of ODEs using semi-discretization techniques (see [17], p. 93 ff.).

For formulating the task of model reduction the first-order system (1) is rewritten in a slightly modified form known from control theory (see e.g. Knobloch, Kwakernaak [6], Zhou, Doyle [18])

ẇ = Aw + Bh  (3)
z = Tw

with T = I_{n}. Here, the first equation is called state equation while the second equation is the output equation. Further, h, w and z are called input, state and output, respectively. An output z is called stationary output if z and h are stationarily related. The system (3) will be referred to as the original system. To find an approximation of the stationary output z one tries to find a system possessing a stable state equation with a dimension m < n, i.e.

ẇ_{m} = A_{m} w_{m} + B_{m} h  (4)
z̃_{m} = T_{m} w_{m}

with an m × m matrix A_{m}, an m × r matrix B_{m}, an n × m matrix T_{m} and an m-dimensional vector function w_{m} : R × Ω → C^{m}. The above system is called reduced-order system to the original system (3) and its stationary output z̃_{m} is called the m-dimensional approximation of z.

There exists a great number of papers concerning the model reduction problem in the context of control theory. Here, h is considered as input or control function and z is the output function which has to be controlled. In some papers the input function is allowed to be random - in most cases multi-dimensional white noise - see e.g. Hyland, Bernstein [5], Moore [7], Wilson [15].

In the present paper the construction of a suitable reduced-order system is based on the assumption that the random dynamics of the stationary output z can be expressed up to a certain degree of accuracy by projections of its values onto a m-dimensional subspace of C^{n}. A similar approach is used in Sirovich, Knight, Rodriguez [13].

To quantify the error of the moment function approximation of the n-dimensional stationary output z by the stationary output z̃_{m} of a reduced-order system in Section 2 error criterions are defined. They are based on L_{1}-distances of the correlation and spectral density functions. For the general case of two stationarily related random vector functions ζ, ζ̃ : R × Ω → C^{n} where ζ̃ is considered as an approximation of ζ, estimates of the error criterions in terms of the second-order moments of ζ̃ and ζ − ζ̃ are given. A central role plays the mean squared Euclidean distance of ζ and ζ̃ defined by E {| ; |ζ − ζ̃| ; |^{2}}. The case of an approximation of ζ by projections of its values onto a m-dimensional subspace of C^{n} is considered in Subsection 2.3. In Subsection 2.4 the problem of an optimal approximation of ζ by those projections is investigated.

In Section 3 the introduced error criterions and the derived estimates are applied to the present case of model reduction approximation. To this end, in Subsection 3.1 starting with the n-dimensional original system (3) differential equations for the m coordinate functions describing the projection of z are derived. It turns out, that these m ODEs also contain the n − m coordinate functions defining the complementary projection of z. Neglecting these terms results in the desired reduced-order system (4) and an additional approximation error, which has to be taken into account.

The goodness of approximation of the correlation and spectral density function of z by those of z̃_{m} is assessed using the error criterions introduced in Section 2. Upper bounds of these criterions in terms of second-order moments of z̃_{m} and z − z̃_{m} are derived in Subsection 3.2. The upper bounds are used to find practicable guidelines for a suitable selection of a reduced-order system. Subsections 3.3 and 3.4 present two methods for deriving proper reduced-order systems which admit small values of the error criterions. The first one is based on mean-square optimal projections while the second method uses eigenvectors of the system matrix A to find suitable subspaces onto which the values of z are projected.

2 Approximation of stationary vector functions

Before low-dimensional approximations of the output z of the system (3) are constructed, the assessment of the approximation goodness of second-order moment functions as R_{zz}(τ) and S_{zz}(α) will be investigated. To this end, error criterions which measure the corresponding approximation error have to be defined. These criterions can be used to quantify the degree of accuracy of a given approximation and to find suitable m-dimensional approximations which guarantee a given level of accuracy. From a practical point of view, the criterions or at least estimates of them should be easy and fast computable.

The starting point is the consideration of two stationarily related random vector functions ζ, ζ̃ : R × Ω → C^{n} where ζ̃ is assumed to be an approximation of ζ. Afterwards, ζ is replaced by the stationary output z of (3) and ζ̃ by the low-dimensional approximation z̃_{m}.

For the considerations in this section the functions ζ and ζ̃ need not to be centered. Therefore, the second-order moment function

Kζζ(τ) := E {ζ(t)ζ^{∗}(t + τ)} = Rζζ(τ) + m_{ζ}m_{ζ}^{∗}

(which is known as covariance function) and the notation of the variance of a scalar stationary random function u

σ_{u}^{2} := E {|u(t) − m_{u}|^{2}} = E {|u(t)|^{2}} − |m_{u}|^{2}

are introduced.

2.1 Error criterions

Definition 1 Let ζ, ζ̃ : R × Ω → C^{n} be stationarily related random vector functions where ζ̃ is considered as an approximation of ζ. The function ζ̂ = ζ̂(t, ω) := ζ(t, ω) − ζ̃(t, ω) is called pathwise approximation error and the quantity

Δ_{0}(ζ, ζ̃) := E {| ; |ζ − ζ̃| ; |^{2}} = E {| ; |ζ̂| ; |^{2}} = ∑_{i=1}^{n} E {|ζ̂_{i}|^{2}}  (5)

is called mean squared Euclidean distance of ζ and ζ̃.

It can be seen that Δ_{0} is related to the second-order moments of ζ̂_{i} = ζ_{i} − ζ̃_{i}, i = 1, …, n, i.e. to the values of the functions K_{ζ̂_{i}ζ̂_{i}}(τ) for τ = 0, only.

In the present case of model reduction approximations for systems of random differential equations not only these moments of the stationary solution but also (cross-) correlation functions R_{ζ_{i}ζ_{j}}(τ) for all τ ∈ R and (cross-) spectral densities S_{ζ_{i}ζ_{j}}(α) for all α ∈ R have to be approximated for i, j = 1, …, n. Therefore, error criterions with respect to the approximation of these moment functions are introduced. They are based on L_{1}-distances of R_{ζ_{i}ζ_{j}}(.) and R_{ζ̃_{i}ζ̃_{j}}(.) as well as of S_{ζ_{i}ζ_{j}}(.) and S_{ζ̃_{i}ζ̃_{j}}(.) which are summed up for i, j = 1, …, n. Subsequently, estimates of these criterions by upper bounds in terms of Δ_{0} are found.

Definition 2

Let W = {Ψ = {Ψ_{ij}}_{i,j=1...n} : Ψ_{ij} ∈ R_{+}, i, j = 1, …, n, ∑_{i,j=1}^{n} Ψ_{ij} = 1} be the set of all

n × n-matrices with non-negative entries whose sum is equal to one. The matrices Ψ ∈ W are called weighting matrices.

Further, let Θ = {θ(.) : R → R_{+}, ∫_{R} θ(τ)dτ = 1} be a set of non-negative scalar func-

tions. The functions θ ∈ Θ are called weighting functions.

Definition 3 Let ζ, ζ̃ : R × Ω → C^{n} be stationarily related random vector functions. Further, let Ψ ∈ W and θ ∈ Θ be a weighting matrix and a weighting function, respectively.

The quantities

Δ^{R}(ζ, ζ̃) := ∫_{R} ∑_{i,j=1}^{n} |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| Ψ_{ij} θ(τ)dτ

and

Δ^{S}(ζ, ζ̃) := ∫_{R} ∑_{i,j=1}^{n} |S_{ζ_{i}ζ_{j}}(α) − S_{ζ̃_{i}ζ̃_{j}}(α)| Ψ_{ij}dα

are called error criterions with respect to the correlation function and spectral density for an approximation of ζ by ζ̃, respectively.

Remark 4 The error criterion Δ^{R}(ζ, ζ̃) is defined using the weighting function θ to allow the subsequent estimates in terms of the variances of ζ̃ and ζ̂. These estimates are based on the Cauchy-Schwarz inequality

|R_{uv}(τ)| ≤ \sqrt{R_{uu}(0)R_{vv}(0)} = σ_{u}σ_{v}

for the cross-correlation function of two stationarily related random functions u and v. With the help of that inequality the terms |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| involved in the integrand of Δ^{R}(ζ, ζ̃) are estimated by positive constants.

2.2 Estimates of the error criterions

In this subsection the introduced error criterions are estimated in terms of secondorder moments of ζ̃ and ζ̂. These estimates allow an easy computation of upper bounds for the error criterions.

Lemma 5 Let ζ, ζ̃ : R × Ω → C^{n} be stationarily related vector functions, ζ̂ := ζ − ζ̃ and θ ∈ Θ. Then for all i, j = 1, …, n it holds

(i)  ∫_{R} |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| θ(τ)dτ ≤ σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j}

(ii)  ∫_{R} |S_{ζ_{i}ζ_{j}}(α) − S_{ζ̃_{i}ζ̃_{j}}(α)| dα ≤ σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j}.

Proof:

(i) Using ζ_{i} = ζ̃_{i} + ζ̂_{i} it follows

R_{ζ_{i}ζ_{j}}(τ) = R_{ζ̃_{i}ζ̃_{j}}(τ) + R_{ζ̃_{i}ζ̂_{j}}(τ) + R_{ζ̂_{i}ζ̃_{j}}(τ) + R_{ζ̂_{i}ζ̂_{j}}(τ)

for i, j = 1, …, n. Cauchy-Schwarz inequality implies

∫_{R} |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| θ(τ)dτ ≤ ∫_{R} (|R_{ζ̃_{i}ζ̂_{j}}(τ)| + |R_{ζ̂_{i}ζ̃_{j}}(τ)| + |R_{ζ̂_{i}ζ̂_{j}}(τ)|) θ(τ)dτ
≤ (σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j} ) ∫_{R} θ(τ)dτ
= σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j}.

(ii) Analogously, in case of the spectral density it follows for i, j = 1, …, n

∫_{R} |S_{ζ_{i}ζ_{j}}(α) − S_{ζ̃_{i}ζ̃_{j}}(α)| dα ≤ ∫_{R} (|S_{ζ̃_{i}ζ̂_{j}}(α)| + |S_{ζ̂_{i}ζ̃_{j}}(α)| + |S_{ζ̂_{i}ζ̂_{j}}(α)|) dα.

Applying the inequality |S_{uv}(α)| ≤ \sqrt{S_{uu}(α)S_{vv}(α)} which holds for cross-- spectral densities of two scalar stationarily related functions u and v (see [2], p. 92) it yields

∫_{R} |S_{ζ_{i}ζ_{j}}(α) − S_{ζ̃_{i}ζ̃_{j}}(α)| dα

≤ ∫_{R} (\sqrt{S_{ζ̃_{i}ζ̃_{i}}(α)S_{ζ̂_{j} ζ̂_{j}}(α)} + \sqrt{S_{ζ̂_{i}ζ̂_{i}}(α)S_{ζ̃_{j} ζ̃_{j}}(α)} + \sqrt{S_{ζ̂_{i}ζ̂_{i}}(α)S_{ζ̂_{j} ζ̂_{j}}(α)}) dα

≤ \sqrt{∫_{R} S_{ζ̃_{i}ζ̃_{i}}(α)dα ∫_{R} S_{ζ̂_{j} ζ̂_{j}}(α)dα} + \sqrt{∫_{R} S_{ζ̂_{i}ζ̂_{i}}(α)dα ∫_{R} S_{ζ̃_{j} ζ̃_{j}}(α)dα}

+\sqrt{∫_{R} S_{ζ̂_{i}ζ̂_{i}}(α)dα ∫_{R} S_{ζ̂_{j} ζ̂_{j}}(α)dα}

= σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j}.

Here, Cauchy-Schwarz inequality and the identity ∫_{R} S_{uu}(α)dα = σ_{u}^{2} for a sta-

tionary scalar function u has been used.

With the help of the assertions of the preceding lemma concerning the approximation of the entries of the correlation function and spectral density matrices the following estimate for the error criterions Δ^{R} and Δ^{S} can be derived. Since both criterions admit identical estimates for the sake of shorter representation the notation Δ^{R,S} is used for ”Δ^{R}, Δ^{S}”.

Theorem 6 Under the assumptions of Lemma 5 it holds for all weighting matrices Ψ ∈ W and all weighting functions θ ∈ Θ

Δ^{R,S}(ζ, ζ̃) ≤ \overline{Δ}

with

\overline{Δ} := \sqrt{∑_{i=1}^{n} σ2 ; ζ̂_{i}} (2\sqrt{∑_{i=1}^{n} σ2 ; ζ̃_{i}} + \sqrt{∑_{i=1}^{n} σ2 ; ζ̂_{i}} ) .

Proof: From the assertion of Lemma 5 and the definition of Δ^{R} and Δ^{S} it follows

Δ^{R}(ζ, ζ̃) = ∑_{i,j=1}^{n} ∫_{R} |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| θ(τ)dτΨ_{ij}

≤ ∑_{i,j=1}^{n} (σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j} ) Ψ_{ij}

(6)

and Δ^{S}(ζ, ζ̃) = ∑_{i,j=1}^{n} ∫_{R} |S_{ζ_{i}ζ_{j}}(α) − S_{ζ̃_{i}ζ̃_{j}}(α)| dαΨ_{ij}

≤ ∑_{i,j=1}^{n} (σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j} ) Ψ_{ij}.

Applying Cauchy-Schwarz inequality and the relation

∑_{i,j=1}^{n} Ψ_{ij}^{2} ≤ ∑_{i,j=1}^{n} Ψ_{ij} = 1

which follows from 0 ≤ Ψ_{ij} ≤ 1, i, j = 1, …, N, it yields

Δ^{R,S}(ζ, ζ̃) ≤ (\sqrt{∑_{i,j=1}^{n} σ2 ; ζ̃_{i}σ2 ; ζ̂_{j}} + \sqrt{∑_{i,j=1}^{n} σ2 ; ζ̂_{i}σ2 ; ζ̃_{j}} + \sqrt{∑_{i,j=1}^{n} σ2 ; ζ̂_{i}σ2 ; ζ̂_{j}} ) \sqrt{∑_{i,j=1}^{n} Ψ_{ij}^{2}}

≤ \sqrt{(∑_{i=1}^{n} σ2 ; ζ̃_{i}) (∑_{i=1}^{n} σ2 ; ζ̂_{i})} + \sqrt{(∑_{i=1}^{n} σ2 ; ζ̂_{i}) (∑_{i=1}^{n} σ2 ; ζ̃_{i})} + \sqrt{(∑_{i=1}^{n} σ2 ; ζ̂_{i}) (∑_{i=1}^{n} σ2 ; ζ̂_{i})}

= \sqrt{∑_{i=1}^{n} σ2 ; ζ̂_{i}} (2\sqrt{∑_{i=1}^{n} σ2 ; ζ̃_{i}} + \sqrt{∑_{i=1}^{n} σ2 ; ζ̂_{i}} )

= \overline{Δ}.

Remark 7 The derived estimates also hold for a modified error criterion with respect to the correlation function defined by

Δ̃^{R}(ζ, ζ̃) := ∑_{i,j=1}^{n} sup_{τ} |R_{ζ_{i}ζ_{j}}(τ) − R_{ζ̃_{i}ζ̃_{j}}(τ)| Ψ_{ij}.

The upper bound \overline{Δ} for the error criterions Δ^{R} and Δ^{S} can be related to the mean squared Euclidean distance Δ_{0}(ζ, ζ̃) given in (5).

Corollary 8 For the upper bound \overline{Δ} given in Theorem 6 it holds

\overline{Δ} ≤ \overline{Δ}_{0} := \sqrt{Δ_{0}(ζ, ζ̃)} (2\sqrt{E \{| ; |ζ̃| ; |^{2}\}} + \sqrt{Δ_{0}(ζ, ζ̃)}) .

Proof: The assertion follows from the representation of \overline{Δ} in Theorem 6 by using the relations

Δ_{0}(ζ, ζ̃) = ∑_{i=1}^{n} E {|ζ̂_{i}|^{2}} ≥ ∑_{i=1}^{n} σ2 ; ζ̂_{i}
and E {| ; |ζ̃| ; |^{2}} = ∑_{i=1}^{n} E {|ζ̃_{i}|^{2}} ≥ ∑_{i=1}^{n} σ2 ; ζ̃_{i}.

For special selections of the weighting matrices Ψ the estimation Δ^{R,S}(ζ, ζ̃) ≤ \overline{Δ} can be improved. The next lemma contains a corresponding assertion for the matrix Ψ = Ψ_{1} = \{\frac{1}{n}δ_{ij}\}_{i,j=1}^{n} which can be used if the auto-correlation functions and autospectral densities have to be approximated, only. Moreover, the case of uniform weights, i.e. Ψ = Ψ_{2} = \{\frac{1}{n^{2}} \}_{i,j=1}^{n} is considered.

Lemma 9 Under the assumptions of Lemma 5 it holds for the weighting matrices

Ψ_{1} = \{\frac{1}{n}δ_{ij}\}_{i,j=1}^{n} and Ψ_{2} = \{\frac{1}{n^{2}} \}_{i,j=1}^{n}

the estimation

Δ^{R,S}(ζ, ζ̃) ≤ \frac{1}{n}\overline{Δ}.

Proof: Using the inequalities (6) found in the proof of Theorem 6 it follows for Ψ = Ψ_{1}

Δ^{R,S}(ζ, ζ̃) ≤ ∑_{i,j=1}^{n} (σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j} ) Ψ_{ij}

= (2 ∑_{i=1}^{n} σζ̂_{i}σζ̃_{i} + ∑_{i=1}^{n} σ2 ; ζ̂_{i}) \frac{1}{n}

≤ (2\sqrt{∑_{i=1}^{n} σ2 ; ζ̂_{i}} \sqrt{∑_{i=1}^{n} σ2 ; ζ̃_{i}} + ∑_{i=1}^{n} σ2 ; ζ̂_{i} ) \frac{1}{n}

= \frac{1}{n}\overline{Δ}

where Cauchy-Schwarz inequality has been applied. Analogously, for Ψ = Ψ_{2} it yields

Δ^{R,S}(ζ, ζ̃) ≤ ∑_{i,j=1}^{n} (σζ̃_{i}σζ̂_{j} + σζ̂_{i}σζ̃_{j} + σζ̂_{i}σζ̂_{j} ) \frac{1}{n^{2}}

= (2 (∑_{i=1}^{n} σζ̂_{i}) (∑_{i=1}^{n} σζ̃_{i}) + (∑_{i=1}^{n} σζ̂_{i})^{2}) \frac{1}{n^{2}}

≤ (2\sqrt{n ∑_{i=1}^{n} σ2 ; ζ̂_{i}} \sqrt{n ∑_{i=1}^{n} σ2 ; ζ̃_{i}} + n ∑_{i=1}^{n} σ2 ; ζ̂_{i} ) \frac{1}{n^{2}}

= \frac{1}{n}\overline{Δ}.

2.3 Projections

The next two subsections deal with the approximation of a random vector function ζ(t, ω) by a function ζ̃(t, ω) which is found by projecting the values of ζ onto a subspace of C^{n} for all t ∈ R and ω ∈ Ω. These projections onto a subspace are defined as follows. Consider the canonical basis (e_{k})_{k=1,...,n} of C^{n} and the representation of ζ(t, ω) with respect to this basis

ζ(t, ω) = ∑_{k=1}^{n} ζ_{k}(t, ω)e_{k}.

Now, let (v_{k})_{k=1,...,n} be a set of linear independent elements of C^{n} with | ; |v_{k}| ; | = 1. Denote by ξ_{1}, …, ξ_{n} the new coordinates of ζ ∈ C^{n} with respect to the basis (v_{k})_{k=1,...,n}, i.e.

ζ(t, ω) = ∑_{k=1}^{n} ξ_{k}(t, ω)v_{k}.  (7)

Defining the C^{n}-valued vector function ξ(t, ω) = (ξ_{1}(t, ω), …, ξ_{n}(t, ω))^{τ} and the matrix V = (v_{1}, …, v_{n}) which is non-singular because of the linear independence of v_{1}, …, v_{n}, the transformation (7) and its inverse transformation can be written as

ζ(t, ω) = Vξ(t, ω)  and  ξ(t, ω) = V^{−1}ζ(t, ω).

Next, sequences (S_{m})_{m=1,...,n} of linear subspaces of C^{n} are considered, where S_{m} is defined by

S_{m} := span(v_{1}, …, v_{m}), m = 1, …, n.

To each sequence (S_{m})_{m=1,...,n} corresponds a sequence (P_{m})_{m=1,...,n} of projectors, where P_{m} : C^{n} → S_{m} is defined by

P_{m}ζ := ∑_{k=1}^{m} ξ_{k}v_{k}, m = 1, …, n,  (8)

as well as a sequence (Q_{m})_{m=1,...,n} of complementary projectors where Q_{m} : C^{n} → span(v_{m+1}, …, v_{n}) is defined by

Q_{m}ζ := ζ − P_{m}ζ, m = 1, …, n.  (9)

Defining for m = 1, …, n − 1 the following subvectors of ξ

ξ_{m} := (ξ_{1}, …, ξ_{m})^{τ} and ξ_{m}^{c} := (ξ_{m+1}, …, ξ_{n})^{τ}

and the submatrices of V

V_{m} := (v_{1}, …, v_{m}) and V_{m}^{c} := (v_{m+1}, …, v_{n})

the relations (8) and (9) for m = 1, …, n − 1 can be written as

P_{m}ζ(t, ω) = V_{m} ξ_{m}(t, ω) and Q_{m}ζ(t, ω) = V_{m}^{c} ξ_{m}^{c}(t, ω).

2.4 Mean-square optimal projections

After defining projections of random vector functions the problem of optimal projections is investigated. Here, the term ”optimal” will be used in the sense of minimization of an error criterion subject to the subspace S_{m} with its defining basis elements v_{1}, …, v_{m}. Natural selections of error criterions which are to be minimized are Δ^{R}(ζ, P_{m}ζ) and/or Δ^{S}(ζ, P_{m}ζ). Because of the unknown moment functions of ζ(t, ω) a general solution of this optimization problem cannot be expected. Here, the objective will be the mean squared Euclidean distance

Δ_{0}(ζ, P_{m}ζ) = E {| ; |ζ − P_{m}ζ| ; |^{2}} = E {| ; |Q_{m}ζ| ; |^{2}}

which only depends on the second-order moments of Q_{m}ζ and is invariant with respect to t ∈ R since ζ is stationary. Similar problems are considered in Gnedin [4], Richter [10].

From Theorem 6 and Corollary 8 it is known, that Δ^{R} and Δ^{S} can be estimated by

Δ^{R,S}(ζ, P_{m}ζ) ≤ \overline{Δ}_{0} = \sqrt{Δ_{0}(ζ, P_{m}ζ)} (2\sqrt{E \{| ; |P_{m}ζ| ; |^{2}\}} + \sqrt{Δ_{0}(ζ, P_{m}ζ)})
= \sqrt{E \{| ; |Q_{m}ζ| ; |^{2}\}} (2\sqrt{E \{| ; |ζ − Q_{m}ζ| ; |^{2}\}} + \sqrt{E \{| ; |Q_{m}ζ| ; |^{2}\}}) .

It can be seen, that for | ; |Q_{m}ζ| ; | ≪ 1 the magnitude of the upper bound \overline{Δ}_{0} is essentially determined by the magnitude of Δ_{0}(ζ, P_{m}ζ) = E {| ; |Q_{m}ζ| ; |^{2}}. Consequently, the minimization of Δ_{0}(ζ, P_{m}ζ) will lead to small upper bounds for the error criterions Δ^{R,S}(ζ, P_{m}ζ).

A similar problem for stationary vector functions which are approximated by projections onto linear submanifolds (instead of subspaces)

S̃_{m} = v_{0} + span(v_{1}, …, v_{m}), m = 1, …, n,

where a translation element v_{0} ∈ C^{n} is involved, is considered in [12].

Since every set of basis elements defines a sequence of subspaces (S_{m})_{m=1,...,n} and a corresponding sequence of projectors (P_{m})_{m=1,...,n} the search for optimal projections is not restricted to a single projector P_{m} for a fixed m but will be considered for sequences of projectors (P_{m})_{m=1,...,n−1}. The case m = n can be excluded since it holds P_{n}ζ = ζ.

Denote by M the set of all sequences (P_{m})_{m=1,...,n−1} of projectors defined in (8). Each sequence can be identified by the n linear independent basis elements v_{1}, …, v_{n} ∈ C^{n},

| ; |v_{k}| ; | = 1, k = 1, …, n.

Definition 10 Let ζ : R × Ω → C^{n} be a stationary random vector function. A sequence \binom{}{P̆_{m}}_{m=1,...,n−1} ∈ M̃ ⊆ M of projectors defined by the basis elements v̆_{1}, …, v̆_{n} ∈ C^{n}

is said to be mean-square optimal with respect to the set M̃ if

Δ_{0}(ζ, P̆_{m}ζ) ≤ Δ_{0}(ζ, P_{m}ζ), m = 1, …, n − 1,

holds for all sequences (P_{m})_{m=1,...,n−1} ∈ M̃.

Next, the case of orthogonal projectors is considered. Denote by N the set of all sequences (P_{m})_{m=1,...,n−1} of projectors with respect to orthonormal basis elements v_{1}, …, v_{n} ∈ C^{n}, i.e. (v_{k}, v_{l}) = δ_{kl}, for 1 ≤ k, l ≤ n. Then N is a subset of M. Subsequently, we consider mean-square optimal sequences \binom{}{P̆_{m}}_{m=1,...,n−1} with respect to N and compare these sequences with the corresponding mean-square optimal sequences for the larger set M.

Lemma 11 [17] p.73 ff., [12] If a sequence \binom{}{P̆_{m}}_{1,...,n−1} ∈ N is mean-square optimal with respect to the set N ⊂ M then it is also mean-square optimal with respect to the set M.

By means of Lemma 11 the search of mean-square optimal projections can be restricted to the case of orthonormal elements v_{1}, …, v_{n}. In this case the transformation matrix V = (v_{1}, …, v_{n}) is an unitary matrix, i.e. V^{−1} = V^{∗}.

The next theorem shows that the mean-square optimal projections are related to the eigenvalues and eigenvectors of the covariance matrix of ζ

Kζζ(0) = E {ζ(t)ζ^{∗}(t)} .

The variance matrix is a Hermitian matrix, i.e. it holds Kζζ(0) = Kζζ^{∗}(0), and it is nonnegative definite. As a standard result of algebraic eigenvalue theory it is known that Kζζ(0) possesses n real and nonnegative eigenvalues σ_{1}^{2} ≥ … ≥ σ_{n}^{2} ≥ 0 and that there exist n pairwise orthogonal eigenvectors u_{1}, …, u_{n} ∈ C^{n}.

Theorem 12 Let ζ : R × Ω → C^{n} be a stationary random function with a covariance matrix Kζζ(0) possessing an eigenvalue decomposition

Kζζ(0) = UΣU^{∗}

with eigenvalues σ_{1}^{2} ≥ … ≥ σ_{n}^{2} forming a diagonal matrix Σ = diag (σ_{1}^{2}, …, σ_{n}^{2}) and corresponding pairwise orthonormal eigenvectors u_{1}, …, u_{n} forming the matrix U = (u_{1}, …, u_{n}).

The sequence \binom{}{P̆_{m}}_{m=1,...,n−1} of projectors defined by the basis elements v̆_{k} = u_{k}, k = 1, …, n is mean-square optimal with respect to M and it holds

Δ_{0}(ζ, P̆_{m}ζ) = ∑_{k=m+1}^{n} σ_{k}^{2}, m = 1, …, n − 1.

Proof: The proof is analogous to the proof of Theorem 10.11 in [17] p.75 ff. where the special case of a centered function ζ is considered.

3 Approximation by model reduction

In the preceding section in a general setup the approximation of stationary random vector functions has been investigated. This section deals with the approximation of the stationary output z of the original system (3) by its m-dimensional approximation z_{m} which is found as the stationary output of a reduced-order system (4).

Subsection 3.1 describes the construction of a reduced-order system and of an mdimensional approximation z̃_{m} which is based on the consideration of projections P_{m}z of z onto m-dimensional subspaces of C^{n}. In Subsection 3.2 the error criterions and its estimates which have been found for the general case of stationary vector functions ζ and ζ̃ are evaluated for the case ζ = z and ζ̃ = z̃_{m}. Based on the derived estimates of Δ^{R}(z, z̃_{m}) and Δ^{S}(z, z̃_{m}) the proper selection of basis elements v_{1}, …, v_{n} defining the projectors P_{m}, m = 1, …, n is investigated. It results in two procedures for the construction of reduced-order systems which are introduced in Subsection 3.3 and 3.4.

3.1 Reduced-order system

A reduced-order system can be derived if an m-dimensional system for the m coordinate functions defining the projection P_{m}z (see Subsection 2.3) is considered. To this end, the original system (3) is written in terms of the transformed coordinates using the transformation w = Vx. Here, x = (x_{1}, …, x_{n})^{τ} contains the coordinates of w with respect to the basis (v_{k})_{k=1,...,n}. Then (3) reads as

ẋ = Cx + Dh  (10)
z = Vx

with C := V^{−1}AV and D := V^{−1}B and is referred to as the transformed original system. Defining

C_{m} := \{C_{ij}\}_{i,j=1,...,m} , C_{m}^{c} := \{C_{ij}\} \begin{matrix}i=1,...,m \\ j=m+1,...,n \end{matrix} , D_{m} := \{D_{ij}\} \begin{matrix}i=1,...,m \\ j=1,...,r \end{matrix} , (11)

i.e.

C = ( \begin{matrix}C_{m} & C_{m}^{c} \\ … & … \end{matrix} ) and D = ( \begin{matrix}D_{m} \\ … \end{matrix} ) ,

for m = 1, …, n − 1, the first m state equations of (10) read as

ẋ_{m} = C_{m}x_{m} + C_{m}^{c}x_{m}^{c} + D_{m}h.  (12)

Assuming a stable matrix C_{m} the system (12) is a stable first-order system for x_{m} with a random excitation term C_{m}^{c}x_{m}^{c} + D_{m}h.

Since z and h and consequently x and h are centered and stationarily related it follows that the excitation term C_{m}^{c}x_{m}^{c} + D_{m}h is stationary and centered. Then, Equation (12) possesses a centered stationary solution which can be represented as

x_{m}(t, ω) = ∫_{0}^{∞}e^{C_{m}s} (D_{m}h(t − s, ω) + C_{m}^{c}x_{m}^{c}(t − s, ω)) ds
= w_{m}(t, ω) + ŵ_{m}(t, ω)

with w_{m}(t, ω) := ∫_{0}^{∞}e^{C_{m}s} D_{m}h(t − s, ω)ds  (13)

ŵ_{m}(t, ω) := ∫_{0}^{∞}e^{C_{m}s} C_{m}^{c}x_{m}^{c}(t − s, ω)ds.  (14)

The m-dimensional equation (12) possesses not the desired form (4) of the state equation of a reduced-order system since its excitation term contains C_{m}^{c}x_{m}^{c}. That means, the stationary solution of the equation (12) and its moment functions cannot be determined exactly without the knowledge of x_{m+1}, …, x_{n} which requires the knowledge of the stationary solution of the large-scale original system.

An exception of this situation is the case of decoupled state equations of the transformed original system (10). This case results if v_{1}, …, v_{n} are the eigenvectors of the diagonalizable matrix A. Then, the matrix C := V^{−1}AV is diagonal which implies C_{m}^{c} = 0_{m,n−m} and C_{m}^{c}x_{m}^{c} = 0_{m}. This exceptional case is described in detail in Subsection 3.4.

To overcome the mentioned drawback of the non-vanishing term C_{m}^{c}x_{m}^{c} a second approximation step is necessary. It consists in the approximation of the variables x_{m+1}, …, x_{n} by their zero means and generates an additional error. Then,

z̃_{m} := V_{m}w_{m}  (15)

instead of the projection P_{m}z = V_{m}x_{m} is used as approximation of z. The matrix V_{m} has been defined in Subsection 2.3. Since the m-dimensional vector function w_{m} given in (13) is the unique stationary solution of

ẇ_{m} = C_{m}w_{m} + D_{m}h

if C_{m} is a stable matrix it follows, that z̃_{m} is the stationary output of the system

ẇ_{m} = C_{m}w_{m} + D_{m}h  (16)
z̃_{m} = V_{m}w_{m}.

The above system may serve as the desired reduced-order system (4) if the settings A_{m} = C_{m}, B_{m} = D_{m} and T_{m} = V_{m} are used. The system (16) possesses the same qualitative properties as the original system. If C_{m} is stable then both systems are stable and possess stationary excitations. That means, procedures for the computation of moment functions of the original system can be applied to the reduced-order system (16), too.

The relation z̃_{m} = V_{m}w_{m} implies that the correlation function R_{z̃_{m}z̃_{m}}(τ) and the spectral density S_{z̃_{m}z̃_{m}}(α) which are used as approximations of the corresponding moment functions of z can be determined from the moment functions of w_{m} using

R_{z̃_{m}z̃_{m}}(τ) = V_{m}R_{w_{m}w_{m}}(τ)V_{m}^{∗}
S_{z̃_{m}z̃_{m}}(α) = V_{m}S_{w_{m}w_{m}}(α)V_{m}^{∗}.

(17)

Now, the pathwise approximation error for the approximation of the stationary solution z(t, ω) of the large-scale original system (3) by z̃_{m}(t, ω) given in (15) which is defined by

ẑ_{m} := ẑ_{m}(t, ω) = z(t, ω) − z̃_{m}(t, ω)  (18)

is considered. It holds

ẑ_{m} = z − z̃_{m} = (P_{m}z + Q_{m}z) − z̃_{m}
= (V_{m}x_{m} + V_{m}^{c}x_{m}^{c}) − V_{m}w_{m}
= V_{m}^{c}x_{m}^{c} + V_{m}(x_{m} − w_{m})
= V_{m}^{c}x_{m}^{c} + V_{m}ŵ_{m} = ^{1}ẑ_{m} +^{2}ẑ_{m}  (19)

with

1ẑ_{m} := V_{m}^{c}x_{m}^{c} = Q_{m}z and ^{2}ẑ_{m} := V_{m}ŵ_{m}.

Here, ^{1}ẑ_{m} results from the first approximation step - the projection of z onto S_{m} - and is called projection error. The second approximation step - the replacement of x_{m+1}, …, x_{n} by their zero means - leads to the error term ^{2}ẑ_{m} which is called defect.

The random function ŵ_{m}(t, ω) given in (14) and used in the definition of the defect ẑ_{m} is the stationary solution of the m-dimensional first-order system

˙ŵ_{m} = C_{m}ŵ_{m} + C_{m}^{c}x_{m}^{c}.

That means, ^{2}ẑ_{m} is the stationary output of the following system

˙ŵ_{m} = C_{m}ŵ_{m} + C_{m}^{c}x_{m}^{c}  (20)
2ẑ_{m} = V_{m}ŵ_{m}

The above system is called defect system.

From the relation ẑ_{m} = z − z̃_{m} = z − V_{m} w_{m} it can be deduced that the pathwise approximation error is the output of the following augmented system with an n + mdimensional state equation

\frac{d}{dt} ( \begin{matrix}z \\ w_{m} \end{matrix} ) = ( \begin{matrix}A & O_{n,m} \\ O_{m,n} & C_{m} \end{matrix} ) ( \begin{matrix}z \\ w_{m} \end{matrix} ) + ( \begin{matrix}B \\ D_{m} \end{matrix} ) h.  (21)
ẑ_{m} = (I_{n}, −V_{m}) ( \begin{matrix}z \\ w_{m} \end{matrix} ) .

The system matrix of the state equation of the above system is stable since A and C_{m} are assumed to be stable. Consequently, it exists a centered stationary output of (21) and the moment functions of the output ẑ_{m} can be determined using the same procedures as in case of the n-dimensional original system (3).

Often, the practical computation of moment functions of the state variables uses the eigenvalue decomposition of the system matrix. In case of the above augmented system the eigenvalue decomposition of the matrix ( \begin{matrix}A & O_{n,m} \\ O_{m,n} & C_{m} \end{matrix} ) is needed. It can

be derived from the corresponding eigenvalue decompositions of the system matrices of the original and the reduced-order sytem. Let these decompositions denoted by

A = VΛV^{−1} and C_{m} = WMW^{−1}

then it yields

( \begin{matrix}A & O_{n,m} \\ O_{m,n} & C_{m} \end{matrix} ) = ( \begin{matrix}V & O_{n,m} \\ O_{m,n} & W \end{matrix} ) ( \begin{matrix}Λ & O_{n,m} \\ O_{m,n} & M \end{matrix} ) ( \begin{matrix}V^{−1} & O_{n,m} \\ O_{m,n} & W^{−1} \end{matrix} ) .

3.2 Estimates of the error criterions

In this subsection the results of Subsection 2.2 for the general problem of moment function approximation for arbitrary stationarily related vector functions ζ and ζ̃ are applied to the present case of model reduction approximation which is obtained with the settings

ζ = z and ζ̃ = z̃_{m} = V_{m}w_{m}.

Now, ζ as well as its approximation ζ̃ are centered functions. The pathwise approximation error can be represented as (see (19))

ζ̂ = ẑ_{m} = ^{1}ẑ_{m} +^{2}ẑ_{m} = V_{m}^{c}x_{m}^{c} + V_{m}ŵ_{m} = V ( ŵ_{m} ; x_{m}^{c} )

or as the stationary output of the augmented system (21). To evaluate the upper bound \overline{Δ}_{0} for the error criterions Δ^{R,S}(z, z̃_{m}) the quantities Δ_{0}(z, z̃_{m}) = ∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\}

and ∑_{i=1}^{n} E {|z̃_{m,i}|^{2}} have to determined. The latter quantity is a by-product of the com-

putation of moment functions of the stationary output of the reduced-order system (16) while for Δ_{0}(z, z̃_{m}) an additional effort is necessary. It consists in the computation of second-order moments of the stationary output ẑ_{m} of the augmented system (21).

The next lemma provides some representations of ∑_{i=1}^{n} E {|z̃_{m,i}|^{2}} and ∑_{i=1}^{n} E {|ẑ_{m,i}|^{2}} in

terms of variances of the stationary solutions of the state equations of the reducedorder system (16) and the augmented system (21). These representations can be useful in the practical computation of the upper bound \overline{Δ}_{0}

Lemma 13 For the m-dimensional approximation z̃_{m} given in (15) and the pathwise approximation error ẑ_{m} given in (18) it holds

∑_{i=1}^{n} E \{|z̃_{m,i}|^{2}\} = ∑_{k,l=1}^{m} E \{w_{m,k}\overline{w_{m,l}}\} (v_{k}, v_{l}) ,  (22)

∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\} = E \{| ; |z| ; |^{2}\} + ∑_{k,l=1}^{m} E \{w_{m,k}\overline{w_{m,l}}\} (v_{k}, v_{l})

−2 Re [∑_{i=1}^{n} ∑_{k=1}^{m} E \{\overline{z_{i}} w_{m,k}\} V_{ik}].  (23)

where (., .) denotes the Euclidean scalar product of C^{n}.

Proof: The relations z̃_{m} = V_{m}w_{m} and ẑ_{m} = z − z̃_{m} imply

z̃_{m,i} = ∑_{k=1}^{m} V_{ik} w_{m,k} and ẑ_{m,i} = z_{i} − ∑_{k=1}^{m} V_{ik} w_{m,k}.

For the second-order moments it follows

E \{|z̃_{m,i}|^{2}\} = ∑_{k,l=1}^{m} E \{w_{m,k}\overline{w_{m,l}}\} V_{ik}\overline{V}_{il}  (24)
E \{|ẑ_{m,i}|^{2}\} = E \{|z_{i}|^{2}\} + ∑_{k,l=1}^{m} E \{w_{m,k}\overline{w_{m,l}}\} V_{ik}\overline{V}_{il}

−2 Re [∑_{k=1}^{m} E \{\overline{z_{i}} w_{m,k}\} V_{ik}].  (25)

Summation for i = 1, …, n and the relation ∑_{i=1}^{n} V_{ik}\overline{V}_{il} = (v_{k}, v_{l}) proves the assertion.

If the model reduction is based on orthogonal elements v_{1}, …, v_{n} which imply orthogonal projectors P_{m}, m = 1, …, n − 1, the property (v_{k}, v_{l}) = δ_{kl}, k, l = 1, …, n, admits simplified representations of the assertions of Lemma 13.

Corollary 14 Under the assumption of Lemma 13 it holds for orthogonal elements v_{1}, …, v_{n}

∑_{i=1}^{n} E {|z̃_{m,i}|^{2}} = ∑_{k=1}^{m} E {|w_{m,k}|^{2}} = E {| ; |w_{m}| ; |^{2}}
∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\} = E \{| ; |z| ; |^{2}\} + E \{| ; |w_{m}| ; |^{2}\} − 2 Re [∑_{i=1}^{n} ∑_{k=1}^{m} E \{\overline{z_{i}} w_{m,k}\} V_{ik}].

The derived representations for the quantities E \{|\left|{ \atop z̃_{m}}\right||^{2}\} = ∑_{i=1}^{n} E \{|z̃_{m,i}|^{2}\} and

Δ_{0} (z, z̃_{m}) = ∑_{i=1}^{n} E {|ẑ_{m,i}|^{2}} given in Lemma 13 and Corollary 14 admit the com-

putation of the upper bound

\overline{Δ}_{0} = \sqrt{Δ_{0} (z, z̃_{m})} (2\sqrt{E \{|\left|{ \atop z̃_{m}}\right||^{2}\}} + \sqrt{Δ_{0} (z, z̃_{m})})

for error criterions Δ^{R} (z, z̃_{m}) and Δ^{S} (z, z̃_{m}) and the assessment of the accuracy of the approximation of moment functions of z. It can be seen that in case of orthogonal elements v_{1}, …, v_{n} only the diagonal entries of the variance matrix of the augmented state \binom{z}{w_{m}} of the system (21) and the submatrix E \{w_{m}z^{∗}\} are needed for the computation of \overline{Δ}_{0}.

In Section 11.2 of [17] additional representations of Δ_{0} (z, z̃_{m}) are derived using an augmented system with the state \binom{x}{ŵ_{m}} instead of system (21) with the state \binom{z}{w_{m}}.

Next, the problem of a proper selection of the reduced-order system is investigated. This leads to the problem of a proper selection of the basis elements v_{1}, …, v_{n} defining the subspaces S_{m} and projectors P_{m}, m = 1, …, n, which are used for the construction of reduced-order systems of dimension m.

Since the goodness of the approximation is measured by the error criterions Δ^{R} (z, z̃_{m}) and Δ^{S} (z, z̃_{m}) those basis elements v_{1}, …, v_{n} are to prefer, which lead to small values of the error criterions or at least of the upper bounds \overline{Δ}_{0}. In Subsection 2.4 an analogous problem has been studied in context with the optimal approximation of z by projections P_{m}z onto a m-dimensional subspace of C^{n}. Again the mean squared Euclidean distance

Δ_{0} (z, z̃_{m}) = E \{||z − z̃_{m}||^{2}\} = E \{|\left|{ \atop ẑ_{m}}\right||^{2}\} = ∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\}

can be used as objective which has to be minimized with respect to v_{1}, …, v_{n}. For Δ_{0} (z, z̃_{m}) ≪ 1 this quantity determines essentially the magnitude of the upper bound \overline{Δ}_{0}.

From the representation of the pathwise approximation error ẑ_{m} as the output of the augmented system (21) it can be seen that the basis elements v_{1}, …, v_{n} are involved in the output equation as well as in the state equation since C_{m} and D_{m} are submatrices of C = V^{−1}AV and D = V^{−1}B, respectively. This implies a rather complicated dependence of the objective Δ_{0} (z, z̃_{m}) = E \{|\left|{ \atop ẑ_{m}}\right||^{2}\} on v_{1}, …, v_{n} and a general solution of the mentioned minimization problem seems to be not available.

To find applicable guidelines for the selection of v_{1}, …, v_{n} the decomposition (see Equation (19))

ẑ_{m} = ^{1}ẑ_{m} +^{2}ẑ_{m}

into the projection error ^{1}ẑ_{m} = V_{m}^{c}x_{m}^{c} = Q_{m}z and the defect ^{2}ẑ_{m} = V_{m}ŵ_{m} is considered. Now, the minimization of E {||^{j}ẑ_{m}||^{2}}, j = 1, 2, subject to v_{1}, …, v_{n} is considered separately. The results are presented in the next two subsections.

3.3 Model reduction based on mean-square optimal projections

From Subsection 2.4 Theorem 12 it is known that

E {| ; |^{1}ẑ_{m}| ; |^{2}} = E {| ; |Q_{m}z| ; |^{2}} = Δ_{0} (z, P_{m}z)

can be minimized if the basis elements v_{1}, …, v_{n} are chosen as the eigenvectors of the variance matrix R_{zz}(0) of z which belong to the eigenvalues σ_{1}^{2} ≥ … ≥ σ_{n}^{2} of R_{zz}(0). The minimal value of the objective is ∑_{k=m+1}^{n} σ_{k}^{2}. It is noted that in the present

case of linear systems with a centered stationary output the variance matrix R_{zz}(0) coincides with the matrix K_{zz}(0) used in Theorem 12 .

The selection of the eigenvectors of R_{zz}(0) as basis elements is not optimal with respect to the objective E \{|\left|{ \atop ẑ_{m}}\right||^{2}\} = E \{||^{1}ẑ_{m} +^{2}ẑ_{m}||^{2}\} but it seems to be a reasonable

selection since it also minimizes E {||x_{m}^{c}||^{2}} and x_{m}^{c} is involved in the input term of the defect system (20) with the output ^{2}ẑ_{m}.

With the selected basis elements the reduced-order system can be derived by executing the following algorithm.

Algorithm 15

Step 1: Find the variance matrix R_{zz}(0) of the output of the original system (3).

Step 2: Determine the eigenvalues σ_{1}^{2} ≥ … ≥ σ_{n}^{2} and corresponding orthonormal eigenvectors v̆_{1}, …, v̆_{n} of R_{zz}(0).

Step 3: Choose m ∈ {1, …, n − 1}.

Step 4: Set

V_{m} = (v̆_{1}, …, v̆_{m})
C_{m} = V_{m}^{∗}AV_{m}
D_{m} = V_{m}^{∗}B.

Step 5: If C_{m} is stable then Step 6 else Step 3.

Step 6: The reduced-order system reads as

ẇ_{m} = C_{m}w_{m} + D_{m}h
z̃_{m} = V_{m}w_{m}.

The moment functions of z̃_{m} = V_{m}w_{m} can be derived from the corresponding moment functions of the m-dimensional stationary solution w_{m} of the stable equation of the reduced-order system in Step 6 applying the transformation (17). For the error criterions with respect to the approximation of the moment functions the following estimate can be obtained.

Theorem 16 For the approximation of the moment functions of the stationary output z of the original system (3) by the moment functions of the stationary output z̃_{m} of the reduced-order system derived in Algorithm 15 it holds

Δ^{R,S}(z, z̃_{m}) ≤ \overline{Δ}_{m} := \sqrt{Δ_{0}(z, z̃_{m})} (2\sqrt{E \{|\left|{ \atop z̃_{m}}\right||^{2}\}} + \sqrt{Δ_{0}(z, z̃_{m})})

where

E \{|\left|{ \atop z̃_{m}}\right||^{2}\} = ∑_{k=1}^{m} E \{|w_{m,k}|^{2}\} = E \{| ; |w_{m}| ; |^{2}\}
Δ_{0}(z, z̃_{m}) = ∑_{k=1}^{n} σ_{k}^{2} + E \{|\left|{ \atop z̃_{m}}\right||^{2}\} − 2 Re [∑_{i=1}^{n} ∑_{k=1}^{m} E \{\overline{z_{i}} w_{m,k}\} V_{ik}].

Proof: The representation of \overline{Δ}_{m} follows from Corollary 8 for the settings ζ = z and ζ̃ = z̃_{m}. Further, the representations of

Δ_{0}(z, z̃_{m}) = ∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\} and E \{|\left|{ \atop z̃_{m}}\right||^{2}\} = ∑_{i=1}^{n} E \{|z̃_{m,i}|^{2}\}

follow immediately from Corollary 14 using E {| ; |z| ; |^{2}} = tr [R_{zz}(0)] = ∑_{k=1}^{n} σ_{k}^{2}.

Theorem 16 shows that the evaluation of the upper bounds \overline{Δ}_{m} requires the knowledge of the second-order moments E {|w_{m,k}|^{2}}, k = 1, …, m, of the eigenvalues σ_{k}^{2}, k = m + 1, …, n of the variance matrix R_{zz}(0) and of the covariance matrix E {w_{m}z^{∗}}. While σ_{k}^{2} and E {|w_{m,k}|^{2}} are by-products of the construction of the reduced-order system and the computation of moment functions for that system an additional effort is necessary for the computation of E {w_{m}z^{∗}}. As it has been pointed out in Subsection 3.1 the practical computation of E {w_{m}z^{∗}} can be performed using the n + m-dimensional augmented system (21).

3.4 Model reduction based on modal decomposition

As it has been announced at the end of Subsection 3.2 in this subsection basis elements v_{1}, …, v_{n} are searched such that to given m the mean squared norm E {||^{2}ẑ_{m}||^{2}} of the defect^{2}ẑ_{m} is minimized. To this end, the eigenvalues λ_{1}, …, λ_{n} and the corresponding eigenvectors v^{∘}_{1}, …, v^{∘}_{n} with | ; |v^{∘}_{k}| ; | = 1, k = 1, …, n, of the matrix A which has been assumed to be stable and diagonalizable are considered. The stability of A implies Re[λ_{i}] < 0, i = 1, …, n. In general, the eigenvectors v^{∘}_{1}, …, v^{∘}_{n} are not orthogonal since A needs not to be a Hermitian matrix.

For basis elements v_{k} = v^{∘}_{k}, k = 1, …, n, the matrices involved in the transformed original system (10)

ẋ = Cx + Dh
z = Vx

are defined by

V = (v^{∘}_{1}, …, v^{∘}_{n})
C = V^{−1}AV = diag (λ_{1}, …, λ_{n}) =: Λ
and D = V^{−1}B.

The transformation z = Vx is referred to as modal decomposition since the state equation of the above system consists of n decoupled equations

ẋ_{i} = λ_{i}x_{i} + ∑_{j=1}^{r} D_{ij}h_{j}, i = 1, …, n.

The matrix C_{m} = {C_{ij}}_{i,j=1,...,m} which is the system matrix of the reduced-order as well as of the defect system is a stable and diagonal matrix, it holds

C_{m} = diag (λ_{1}, …, λ_{m}) .

It follows, that for all m = 1, …, n the state equations of the reduced-order system (16) are stable and coincide with the first m equations of the transformed original system. The latter property is a consequence of the relation

C_{m}^{c} = \{C_{ij}\} \begin{matrix}i=1,...,m \\ j=m+1,...,n \end{matrix} = O_{m,n−m}.

Another consequence of this relation is that for the vector functions ŵ_{m} defined in (14) it holds ŵ_{m} = ŵ_{m}(t, ω) ≡ 0_{m}, i.e. x_{m} = w_{m} + ŵ_{m} = w_{m} and P_{m}z = z̃_{m} = V_{m}x_{m}. Further, for the defect ^{2}ẑ_{m} it follows

2ẑ_{m}(t, ω) = V_{m}ŵ_{m}(t, ω) ≡ 0_{2n}.

It results, that for basis elements v_{k} = v^{∘}_{k}, k = 1, …, n, the value of the objective E {||^{2}ẑ_{m}||^{2}} is equal to zero and therefore minimal.

Remark 17 Since the eigenvectors v^{∘}_{1}, …, v^{∘}_{n} of A need not to be orthogonal the projections P_{m}z defined by those basis elements are oblique projections, in general, and cannot be mean-square optimal (see Lemma 11). It turns out that the mean squared norms of the two terms ^{1}ẑ_{m} and ^{2}ẑ_{m} forming the pathwise approximation error ẑ_{m} cannot be minimized simultaneously.

The property E {||^{2}ẑ_{m}||^{2}} = 0 holds for every permutation of indices of the sequence

of basis elements v^{∘}_{1}, …, v^{∘}_{n} while the value of the mean-squared Euclidean distance

Δ_{0}(z, z̃_{m}) = E \{|\left|{ \atop ẑ_{m}}\right||^{2}\} = E \{| ; |^{1}ẑ_{m} +^{2}ẑ_{m}| ; |^{2}\} = E \{| ; |^{1}ẑ_{m}| ; |^{2}\}

depends on the arrangement of the basis elements. It remains to find an arrangement which leads to a minimal or at least to a small value of Δ_{0}(z, z̃_{m}). A natural suggestion is to arrange v^{∘}_{1}, …, v^{∘}_{n} such that for the variances of the stationary solutions x_{k} of the decoupled state equations of the transformed original system it holds

E {|x_{1}|^{2}} ≥ … ≥ E {|x_{n}|^{2}} .

This arrangement ensures that an m-dimensional reduced-order system contains those m state variables of the decoupled system

ẋ = Λx + V^{−1}Bh

which possess the largest variances.

The above suggestion can be supported by the following estimate of Δ_{0}(z, z̃_{m}) which results from Equation (23). Using ŵ_{m} = 0_{m} and applying Cauchy-Schwarz inequality it follows

Δ_{0}(z, z̃_{m}) = ∑_{i=1}^{n} E \{|ẑ_{m,i}|^{2}\} = ∑_{k,l=m+1}^{n} E \{x_{k}\overline{x_{l}}\} (v^{∘}_{k}, v^{∘}_{l})

≤ ∑_{k,l=m+1}^{n} \sqrt{E \{|x_{k}|^{2}\} E \{|x_{l}|^{2}\}} = ( ∑_{k=m+1}^{n} \sqrt{E \{|x_{k}|^{2}\}})^{2} .

The following algorithm summarizes the above derivation of a reduced-order system.

Algorithm 18

Step 1: Find the eigenvalues λ_{1}, …, λ_{n} and the corresponding eigenvectors v^{∘}_{1}, …, v^{∘}_{n} of the matrix A.

Step 2: Set Λ = diag(λ_{1}, …, λ_{n}) and V = (v^{∘}_{1}, …, v^{∘}_{n}).

Step 3: Determine the variances E {|x_{i}|^{2}}, i = 1, …, n, of the stationary solution of

the system

ẋ = Λx + V^{−1}Bh.

Step 4: Rearrange the equations of the system in Step 3 and the corresponding eigenvalues λ_{1}, …, λ_{n} and eigenvectors v^{∘}_{1}, …, v^{∘}_{n} such that it holds

E {|x_{1}|^{2}} ≥ … ≥ E {|x_{n}|^{2}} .

Step 5: Choose m ∈ {1, …, n − 1}.

Step 6: Set

C_{m} = diag (λ_{1}, …, λ_{m})
V_{m} = (v^{∘}_{1}, …, v^{∘}_{m}) , V = (v^{∘}_{1}, …, v^{∘}_{n})
D_{m} = \{D_{ij}\} \begin{matrix}i=1,...,m \\ j=1,...,r \end{matrix} where D = V^{−1}B.

Step 7: The reduced-order system reads as

ẇ_{m} = C_{m}w_{m} + D_{m}h
z̃_{m} = V_{m}w_{m}.

For the error criterions with respect to the approximation of the moment functions of z the following estimates can be derived.

Theorem 19 For the approximation of the moment functions of the stationary output z of the original system (3) by the moment functions of the stationary output z̃_{m} of the reduced-order system derived in Algorithm 18 it holds

Δ^{R,S}(z, z̃_{m}) ≤ \overline{Δ}_{m} := \sqrt{Δ_{0}(z, z̃_{m})} (2\sqrt{E \{|\left|{ \atop z̃_{m}}\right||^{2}\}} + \sqrt{Δ_{0}(z, z̃_{m})})

where

E \{|\left|{ \atop z̃_{m}}\right||^{2}\} = ∑_{i=1}^{n} E \{|z̃_{m,i}|^{2}\} = ∑_{k,l=1}^{m} E \{x_{k}\overline{x_{l}}\} (v^{∘}_{k}, v^{∘}_{l})
Δ_{0}(z, z̃_{m}) = ∑_{i=m+1}^{n} E \{|ẑ_{m,i}|^{2}\} = ∑_{k,l=m+1}^{n} E \{x_{k}\overline{x_{l}}\} (v^{∘}_{k}, v^{∘}_{l}) ,

where x = (x_{1}, …, x_{n})^{τ} denotes the stationary solution of the decoupled system derived in Step 3 of Algorithm 18.

Proof: The relations ŵ_{m}(t, ω) = 0_{m}, x_{m}(t, ω) = w_{m}(t, ω) and P_{m}z(t, ω) = z̃_{m}(t, ω) = V_{m}w_{m}(t, ω) imply

E \{|\left|{ \atop z̃_{m}}\right||^{2}\} = E \{| ; |P_{m}z| ; |^{2}\} = E \{| ; |∑_{k=1}^{m} x_{k}v^{∘}_{k}| ; |^{2}\} = ∑_{k,l=1}^{m} E \{x_{k} \overline{x_{l}}\} (v^{∘}_{k}, v^{∘}_{l}) ,

and Δ_{0}(z, z̃_{m}) = Δ_{0}(z, P_{m}z) = E {| ; |Q_{m}z| ; |^{2}} = E {| ; | ∑_{k=m+1}^{n} x_{k}v^{∘}_{k}| ; |^{2}}

= ∑_{k,l=m+1}^{n} E \{x_{k} \overline{x_{l}}\} (v^{∘}_{k}, v^{∘}_{l})

Acknowledgement

The support by the Deutsche Forschungsgemeinschaft (DFG) is greatfully appreciated.

References

[1] L. Arnold and V. Wihstutz. Stationary solutions of linear systems with additive and multiplicative noise. Stochastics, 7:133–155, 1982.

[2] A.V. Balakrishnan. Introduction to random processes in engineering. Wiley, New York, 1995.

[3] H. Bunke. Gewöhnliche Differentialgleichungen mit zufälligen Parametern. Akademie-Verlag, Berlin, 1972.

[4] A.V. Gnedin. On mean-square epsilon dimension. Math. Jap., 37(4):623–627, 1992.

[5] D.C. Hyland and D.S. Bernstein. The optimal projection equations for model reductions and the relationships among the methods of Wilson, Skelton, and Moore. IEEE Trans. Autom. Control, AC-30(12):1201–1211, 1985.

[6] H.W. Knobloch and H. Kwakernaak. Lineare Kontrolltheorie. Akademie-Verlag, Berlin, 1986.

[7] B.C. Moore. Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Trans. Autom. Control, AC-26:17–32, 1981.

[8] D.E. Newland. An introduction to random vibrations, spectral analysis & wavelet analysis. Longman Scientific & Technical, Harlow, Essex, 3rd edition, 1993.

[9] A. Preumont. Random vibration and spectral analysis. Solid Mechanics and Its Applications. 33. Kluwer Academic Publishers, Dordrecht, 1994.

[10] M. Richter. Approximation of Gaussian random elements and statistics. Teubner, Leipzig, 1992.

[11] J. vom Scheidt, B. Fellenberg, and U. Wöhrl. Analyse und Simulation stochastischer Schwingungssysteme. Leitfäden der angewandten Mathematik und Mechanik, Band 71. B. G. Teubner, Stuttgart, 1994.

[12] J. vom Scheidt, H.-J. Starkloff, and R. Wunderlich. Optimal low-dimensional approximations of random vector functions. Preprint 98-31, TU Chemnitz, Faculty of Mathematics, http://www.mathematik.tu-chemnitz.de/preprint/, 1998.

[13] L. Sirovich, B.W. Knight, and J.D. Rodriguez. Optimal low-dimensional dynamical approximations. Q. Appl. Math., 48(3):535–548, 1990.

[14] T.T. Soong and M. Grigoriu. Random vibration of mechanical and structural systems. PTR Prentice Hall, Englewood Cliffs NJ, 1993.

[15] D.A. Wilson. Model reduction for multivariable systems. Int. J. Control, 20(1):57–64, 1974.

[16] P.H. Wirsching, T.L. Paez, and K. Ortiz. Random vibrations - theory and praxis. Wiley, New York, 1995.

[17] R. Wunderlich. Systems of random differential equations and model reduction. Habilitation Thesis, TU Chemnitz, 1999.

[18] K. Zhou and J.C. Doyle. Essentials of robust control. Prentice Hall, Upper Saddle River, New Jersey, 1998.