The response of grandstands driven by filtered Gaussian white noise
processes
Ondřej Rokoša,∗, Jiřı́ Mácaa

arXiv:1509.04250v1 [] 14 Sep 2015

a

Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague,
Thákurova 7, 166 29 Prague 6, Czech Republic.

Abstract
This paper presents a semi-analytical estimate of the response of a grandstand occupied
by an active crowd and by a passive crowd. Filtered Gaussian white noise processes are
used to approximate the loading terms representing an active crowd. Lumped biodynamic
models with a single degree of freedom are included to reflect passive spectators occupying
the structure. The response is described in terms of the first two moments, employing the Itô
formula and the state augmentation method for the stationary time domain solution. The
quality of the approximation is compared on the basis of three examples of varying complexity
using Monte Carlo simulation based on a synthetic generator available in the literature. For
comparative purposes, there is also a brief review of frequency domain estimates.
Keywords: grandstand response, random vibration, active crowd, white noise process,
filtration
1. Introduction
The response of a grandstand can be resolved quite easily by linear dynamics methods,
if we neglect all the randomness of the system. However, for a more accurate description, at
least the most significant uncertainties need to be taken into account. The main uncertainties
include
• forcing terms resulting from active crowd movements, especially synchronized jumping,
• the uncertainties of the parameters in discrete biodynamic models—randomness of
stiffness, mass and damping matrices,
• the size and spatial distribution of an active crowd and a passive crowd.
Further generalizations can take into account various kinds of nonlinearities, e.g. geometrical
and material non-linearities and non-linearities of biodynamic models (Huang and Griffin,
2008). However, the following restrictions will be assumed for the purposes of this paper:
∗

Corresponding author.
Email addresses: ondrej.rokos@fsv.cvut.cz (Ondřej Rokoš), maca@fsv.cvut.cz (Jiřı́ Máca)

Preprint submitted to arxiv

March 17, 2018

the material parameters of the structure are treated as deterministic, since their influence is
negligible in comparison with the sources listed above and the scope of the overall response;
the spatial distribution of the crowd is fixed; biodynamic models of the passive crowd are
treated as deterministic.
The results of measurements carried out on simple structures indicate that an active
spectator might be represented as a time dependent force, cf (Ellis and Ji, 1994). In other
words, under certain circumstances, an active spectator does not influence the properties of
the structure. Several generators of normalized artificial load processes have been developed
in the literature, e.g. (Sim, 2006) and (Racic and Pavic, 2010), which can be supplemented
by human body weights, e.g. (Hermanussen et al., 2001). However, the passive part of the
crowd is assumed to be stationary in space and in permanent contact with the structure,
excluding accelerations that exceed the gravity of the earth. Thus a passive spectator can be
modeled as a biodynamic system. A survey of models of this type can be found in (Sachse
et al., 2003). These remarks allow us to model the structure considered here in terms of
random vibrations of linear systems.
Monte Carlo simulation (MC) is one of the ways used in advanced grandstand design
procedures for reflecting all kinds of randomness. MC is a general method, but it has considerable disadvantages in its original form, e.g. rather slow convergence when estimating low
probabilities, and high computational demands. MC will be employed in this paper mainly
for controlling the accuracy and the performance of semi-analytical methods introduced in
terms of stochastic differential equations. These simplified methods can be applicable and
useful in the preliminary stages of design procedures, when quick and approximate solutions
are sufficient.
2. Stochastic differential equations
Employing the Finite Element Method, and taking into account the above-mentioned
assumptions, a mathematical model of the structural system considered here can be written
as a set of hyperbolic differential equations
M Z̈(t) + C Ż(t) + KZ(t) = GY (t),

t≥0

(1)

0

where Z and Y are Rd/2 and Rd -valued stochastic processes, and M , C, K and G are
(d/2, d/2) and (d/2, d0 ) matrices of mass, damping, stiffness and input distribution. Overdot denotes a derivative with respect to time ġ(t) = d/dt g(t). To simplify subsequent
expressions, let us apply expectation operator E in equation (1) and substract the result
from (1). We arrive at the set of two equations
M µ̈Z (t) + C µ̇Z (t) + KµZ (t) = GµY (t), t ≥ 0
¨ + C Z̃(t)
˙
M Z̃(t)
+ K Z̃(t) = GỸ (t), t ≥ 0

(2)
(3)

for the mean value µZ (t) = EZ(t) and centered process Z̃(t) = Z(t) − µZ (t), µY (t) =
EY (t) and Ỹ (t) = Y (t) − µY (t). As will become apparent in section 3, under certain
2

conditions process Z̃(t) is approximately normal, and since differential equation (3) is linear
with deterministic coefficients, it is reasonable to accept Gaussian approximation also for
Ỹ (t). This consideration leads us to stochastic differential equations and to the Itô calculus.
Under Gaussian assumptions, the response will be completely specified by its mean µZ (t)
and covariance cZ (t, s) = E[Z̃(t)Z̃(s)T ]. The required quantities can be obtained from the
time domain or from the frequency domain.
2.1. Solution in the time domain
Equation (2) can be solved by direct integration or, more conveniently, by a Fourier series
(or Fourier transform) assuming periodic mean µY (t), cf section 2.2. Let us rewrite (3) as
"
# 
# 
"

Z̃(t)
d Z̃(t)
0
I
0
=
+
Ỹ (t)
(4)
˙
˙
−M −1 K −M −1 C
M −1 G
dt Z̃(t)
Z̃(t)
or, in a more compact form,
˙
X̃(t)
= aX̃(t) + bỸ (t),

t≥0

(5)

where X̃ is an Rd -valued state-space vector stochastic process with zero mean, and a and b
are (d, d) and (d, d0 )-matrices. The solution of this differential equation is given in the form
Z t
θ(t − s)bỸ (s) ds
(6)
X̃(t) = θ(t)X̃(0) +
0

where θ(t − s) denotes a Green function or the unit impulse response satisfying
∂θ(t − s)
= aθ(t − s),
∂t

t ≥ s ≥ 0,

(7)

θ(0) = I the identity and θ(t − s) = exp[a(t − s)] can be expressed as a matrix exponential,
cf (Soong and Grigoriu, 1993). Initial conditions X̃(0) will be set to zero for simplicity.
Forcing term Ỹ (t) can also satisfy its own stochastic differential equation driven by Gaussian
white noise W (t) = dB(t)/dt. For example, let Ŷ1 (t) be a continuous-time Gaussian autoregression scalar process of order p, denoted as AR(p), cf (Brockwell et al., 2007). Then
Ŷ1 (t) = eT1 S 1 (t) where the state vector S 1 (t) = [S1,1 (t), . . . , S1,p (t)]T satisfies the Itô equation





A1 = 



0
0
..
.

1
0
..
.

0
0
−ap −ap−1

dS 1 (t) = A1 S 1 (t)dt + b1 dB(t),

 
0
... 0
1
 0 
1
... 0 

 
..
..  , e =  .. 
...
 . 
1
.
. 

 

 0 
0
... 1
−ap−2 . . . −a1
0

(8)


0
0
..
.




and b1 = 

 0
a0





.



Processes of this kind are also called filtered white noise processes or colored processes, and
they have a specific frequency content. Let us assume that Ỹi (t) of Ỹ (t) = [Ỹ1 (t), . . . , Ỹd0 (t)]T
3

are mutually independent AR(pi ) processes. Then we
obtain one coupled system

 
Z̃(t)
0
I
0
 ˙
  −M −1 K −M −1 C M −1 Gd1 eT
1
 Z̃(t)  

 
0
0
A
1
d  S 1 (t)  = 

 
.

 
..
0
0
0
0
0
0
S d0 (t)


0 ... 0
 0 ... 0 




+  b1 . . . 0  dB(t)


 0 ... 0 
0 . . . bd0

can merge equations (4) and (8) to

...
0
−1
. . . M Gdd0 eTd0
...
0
..
.
0
...

 


Z̃(t)
 ˙

 Z̃(t) 


 S 1 (t)  dt +


..


.
S d0 (t)








Ad0

(9)

0

where di are column vectors with the unit in i-th position, and B(t) is an Rd -valued Brownian motion. This approach P
is called a state augmentation method (Grigoriu, 2002). An
extension to the case Ỹi (t) = nk=1 Ŷk (t), where Ŷk (t) are mutually independent AR(p) processes, is carried out in an obvious manner. This methodology will be employed in section 3
for AR(2) processes. Equation (9) can again be rewritten in compact form
dX(t) = aX(t)dt + bdB(t),

t ≥ 0,

(10)

and employing the Itô formula for semimartingales we arrive at the system of evolutionary
equations for the response mean µX (t) and covariance cX (t, s)
µ̇X (t) = aµX (t), t ≥ 0,
ċX (t, t) = acX (t, t) + cX (t, t)aT + bbT ,
∂cX (t, s)
= acX (t, s), t > s ≥ 0.
∂t

t ≥ 0,

(11)
(12)
(13)

Since the driving forces dB(t) are Gaussian white noise and the coefficients are constant in
time, the solution is an Ornstein-Uhlenbeck process with an existing stationary solution. In
our case, stationary mean µX = 0 and covariance ċX (t, t) = ċX (t − t) = ċX = 0 which leads
to the so-called continuous Lyapunov equation
0 = acX + cX aT + bbT .

(14)

For details and further developments, see (Grigoriu, 2002). Since the stationary matrix cX
contains only response displacements and velocities, the variances of the acceleration have to
be computed through the following formulas which are valid for weakly stationary processes
cẊ = −

d2 cX (t)
dt2
4

= −a2 cX ,
t=0

(15)

where a2 denotes matrix power and cẊ denotes the stationary covariance matrix of velocities
and accelerations. Equation (15) is evaluated employing (13), which in our special case is
simplified to
cX (t) = θ(t)cX = exp[at]cX .
(16)
2.2. Solution in the frequency domain
Taking the Fourier transform of equation (2) leads to
µ̂Z (ω) = H(ω)Gµ̂Y (ω)

(17)

where the FRF (Frequency Response Function) H(ω) is
H(ω) = [−ω 2 M + iωC + K]−1 ,

(18)

i denotes a complex unit,

ω denotes angular frequency and ĝ(ω)
√ the Fourier transform
P denotes
harm
of function g(t). Assuming µY (t) in periodic form µ̂Y (ω) = nk=1
2π/2 r k [exp{iϕk }δ(ω −
ω̄k ) + exp{−iϕk }δ(ω + ω̄k )] where r k , ϕk and ω̄k are amplitude, phase shift and angular
frequency of the k-th harmonic, we obtain the mean response as a solution of the complex
linear system of equations.
Employing spectral decomposition of stationary random processes, the response variances
are acquired in terms of spectral density matrices S Ỹ Ỹ and S Z̃ Z̃ , (SỸ Ỹ (ω))ii = fˆỸ (ω), where
fˆỸ (ω) is a spectral density estimate of the centered forcing term
Z ∞
ˆ
b(x − ω)IT (ω) dx,
(19)
fỸ (ω) = E
−∞

b(x) is some weight function, cf (Anděl, 1976), and IT (ω) denotes the corresponding periodogram
Z T
2
1
−itω
IT (ω) =
Ỹ (t)e
dt ,
−∞ < ω < ∞.
(20)
2πT 0
The diagonal form of S Ỹ Ỹ suggests that we treat all input processes as independent. Knowing
the spectral density matrix of the input vector stochastic process, we obtain the spectral
density matrix of the output process according to Soong and Grigoriu (1993)
S Z̃ Z̃ (ω) = H(ω)GS Ỹ Ỹ (ω)GT H † (ω)

(21)

where H † (ω) denotes a Hermitian transpose to H(ω). The variance of the stationary scalar
process Z̃(t) with two-sided spectral density fZ̃ (ω) or with one-sided spectral density gZ̃ (ω)
is evaluated as
Z ∞
Z ∞
2
σZ̃ =
fZ̃ (ω) dω =
gZ̃ (ω) dω
(22)
−∞

0

˙
and the variance of time derivative Z̃(t)
Z ∞
Z
2
2
2
σ ˙ = σ̇Z̃ =
ω fZ̃ (ω) dω =
Z̃

−∞

∞

ω 2 gZ̃ (ω) dω.

0

By analogy for higher time derivatives, applying higher powers of angular frequency ω.
5

(23)

2.3. Transformation to modal coordinates
It is desirable to reduce the system unknowns in equations (1), (2), (14) and (21) for
MC, mean value, time domain and frequency domain solutions. The usual modal transform technique can be employed directly when a structure with proportional damping is
loaded by active spectators only. Some problems arise when the passive part of the crowd
is introduced. All matrices are extended according to the added degrees of freedom (dofs).
Biodynamic models possess non-proportional damping, so complex eigenvectors should be
used. Simultaneously, all these models have the first eigenfrequency close to 5 Hz (and the
second eigenfrequency close to 8 Hz if using two-degrees-of-freedom models). This implies
that the number of eigenvectors employed in the transform is considerably enlarged. Another
possibility is partial modal transformation. The system matrices can be decomposed into
four parts according to


ASS ASH
A=
(24)
AHS AHH
where A stands either for mass, stiffness or damping matrix. Sub-matrix ASS corresponds
to the structure, AHH to the passive crowd, and ASH or AHS represents the mutual interactions. Let us note that sub-matrix AHH is a diagonal or a band, and ASS is sparse.
Computing the eigenvectors corresponding to an empty structure, i.e. corresponding to ASS ,
arranged in Ṽ and writing


Ṽ 0
V =
,
(25)
0 I
all equations can be partially transformed by V , as in the standard procedure. Then the dofs
inherent to a passive crowd are unchanged, but the structure response is described through
several modal coordinates.
Concerning the Lyapunov equation (14) partially transformed into modal coordinates,
we can further reduce the computational effort by some prior information. Let us assume
that the system response corresponds to the equation (9) with pi = 2. Splitting all matrices
in (14) leads to
 



 
 T
0 0
a11 a12
c11 c12
c11 c12
a11 aT21
+
= 0,
(26)
+
aT12 aT22
0 Q22
a21 a22
c21 c22
c21 c22
where


a11 =

0
I
−[V T M V ]−1 [V T KV ] −[V T M V ]−1 [V T CV ]


,

(27)

a21 = 0 and remaining sub-matrices have obvious structure. Dropped subscript at covariance
matrix c emphasizes partial transformation to modal coordinates. Since bbT resp. Q22 are
symmetric, in fact diagonal, the solution will be also symmetric, c12 = cT21 . In the case
of AR(2) processes, sub-matrix c22 can be computed explicitly; single AR(2) process has
uncorrelated state variables S1 , S2 and varS1 = a20 /(2a1 a2 ), varS2 = a20 /(2a1 ), based on
moment equations, thus c22 is a diagonal matrix. Introduced considerations reduce the

6

system of four equations (26) in expanded form, to the set of two coupled equations
a11 c12 + c12 aT22 + a12 c22 = 0
a11 c11 + c11 aT11 + c12 aT12 + a12 cT12 = 0

(28)
(29)

for unknowns c12 and c11 . The set resembles Sylvester and Lyapunov equations respectively with reduced size. Backward transformation c̄X = V c̄V T gives covariance matrix for
displacement or velocity vector, here c̄ denotes an appropriate sub-matrix of c storing the
modal displacements or velocities, c̄X then contains the nodal displacements or velocities.
2.4. Crossings of Gaussian processes
One of the measures of system performance is level crossing. Under some circumstances,
in a stationary case, it can be shown that up-crossing is directly connected with the reliability
of the system. The x-up-crossing rate of Gaussian process X(t) with non-stationary mean
value µ(t) and stationary variance σ 2 is estimated as (Soong and Grigoriu, 1993)

 

 

µ̇(t)
µ̇(t)
µ̇(t)
x − µ(t)
σ̇
+
Φ
φ
+
φ
,
(30)
νx (t) =
σ
σ̇
σ̇
σ̇
σ
√
where
νx+ (t) is the x-up-crossing rate of level x at time t, φ(α) = 1/ 2π exp −α2 /2, Φ(u) =
Ru
φ(α) dα, σ 2 = varX(t), σ̇ 2 = varẊ(t). The total mean number of upcrossings in time
−∞
interval [0, T ] is computed according to
Z T
+
νx+ (t) dt.
(31)
nx (T ) =
0

The relations can be generalized to D-out-crossings of a d-valued stochastic process, where
D is some set in Rd .
Another measure of system performance from the point of view of serviceability is the
root mean square value estimated as
s
s
Z T
Z
1
1 T
2
RM S =
X(t) dt =
µX (t)2 dt + σ 2 ,
(32)
T 0
T 0
where σ 2 denotes the stationary variance of centered process X̃(t), and µX (t) its mean value.
Analogous formulas are valid for velocity and acceleration.
3. Applications to the response of grandstands
As was noted in section 1, an active spectator can be treated as a time-dependent process.
Figure 1 shows a single realization and the spectral density of a unit process, i.e. of the
process with GH = 1, where GH denotes the weight of a spectator, jumping frequency
f¯ = 2.67 Hz. Spectral densities were computed from (19) with Parzen weight. The realization
was generated according to Sim (2006). Since this function is highly periodic, we will search
7

power [dB]

unit process

50

4

2

0

0

1

2

3

4

0
−50
−100
−150

5

0

t [s]

(a)

10

20
30
frequency [Hz]

40

50

(b)

Figure 1: Single time history (a) and power spectral density of 10 000 realizations (b) of
forcing term Y (t) generated according to Sim (2006).
the mean value in the form
µY (t) = α0 +

p
X

αk cos(k · 2π f¯t) + βk sin(k · 2π f¯t).

(33)

k=1

Then vector α̂ of the estimated parameters α̂0 , α̂1 , . . . , α̂p , β̂1 , . . . , β̂p can be found by the
linear Least Squares Method as
α̂ = (ΦT Φ)−1 ΦT µ̄Y ,

(34)

where



Φ=


1 cos(2π f¯t1 ) sin(2π f¯t1 ) . . . cos(p · 2π f¯t1 ) sin(p · 2π f¯t1 )
1 cos(2π f¯t2 ) sin(2π f¯t2 ) . . . cos(p · 2π f¯t2 ) sin(p · 2π f¯t2 )
..
..
..
..
..
...
.
.
.
.
.
1 cos(2π f¯tn ) sin(2π f¯tn ) . . . cos(p · 2π f¯tn ) sin(p · 2π f¯tn )




,


t1 , . . . , tn is a fine enough and equidistant partition of the time interval, µ̄Y = [µ̄Y (t1 ), . . . , µ̄Y (tn )]T
with
N
1 X
Yk (ti ),
µ̄Y (ti ) =
N k=1
are means over N realizations Yk (ti ) in time instants ti . Generating 10 000 trajectories
provides the coefficients given in tab. 1. Realization centered with the mean value according
to equation (33) and the coefficients from table 1 is depicted in figure 2, together with
the spectral density and a normalized histogram,
i.e. a √
histogram of the process Ȳ (t) =
p
[Y (t) − µ(t)]/σY = Ỹ (t)/σY where σY = varY (t) = 0.7627 is a stationary standard
deviation. The centered process resembles non-Gaussian colored noise. To be more precise,
for Ȳ (t) we have µȲ = 0 for the mean, varȲ = 1 for the variance, γ3,Ȳ = 0.424 for the
coefficient of skewness, and γ4,Ȳ = 4.076 for the coefficient of kurtosis. For comparison, the
standard Gaussian process has coefficients 0, 1, 0 and 3.
Let us briefly analyze the response of a harmonic oscillator with unit mass forced by
jumping process Y (t) with f¯ = 2.67 Hz, employing MC to justify the normality assumptions.
8

Table 1: Coefficients α̂ for an approximation of the mean value in equation (33) for f¯ = 2.67
Hz.
α̂0
α̂1
α̂2
α̂3
α̂4

0.9958
0.2939 β̂1 1.1170
-0.2471 β̂2 0.0984
-0.0037 β̂3 -0.0153
-0.0008 β̂4 -0.0001

power [dB]

intensity

2

0

−2

0

10

20

30
time [s]

40

50

0
−50
−100
−150

60

0

10

20
30
frequency [Hz]

(a)
1.5
histogram
std normal density

1

Density

Density

50

(b)

1.5

0.5
0

40

−2

0

2
Data

4

0.5
0

6

(c)

histogram
std normal density

1

−2

0

2
Data

4

6

(d)

Figure 2: Single centered time history of Ỹ (t) (a), the corresponding power spectral density
of 10 000 realizations (b), a normalized histogram with standard normal density (c) and a
normalized histogram of a non-unit process scaled with GH [kN] ∼ N (0.7709, 0.0167) (d).

9

displacement

velocity
0.2

γ3

0.4
0.2

0

0
−0.2

−0.2
5

10

15

20

15

0
15

10

10

5

5

5

10

15

20

5

10
f1 [Hz]

15

20

γ4

0

0

0

5

10
f1 [Hz]

15

20

0

0

Figure 3: Coefficient of skewness γ3 and coefficient of kurtosis γ4 for the normalized displacement and velocity of the harmonic oscillator forced by the jumping process as functions of
eigenfrequency f1 . Solid line - viscous damping ζ = 0.001; dash-dot line - ζ = 0.07; dashed
line - corresponding values for the Gaussian random variable.
The coefficients of skewness and kurtosis of the state vector X(t) = [Z(t), Ż(t)]T as functions
of the oscillator eigenfrequency f1 for two different values of viscous damping ζ are depicted
in figure 3. Note that for frequency range 0.5 − 7 Hz the response is approximately Gaussian.
As was expected, worse convergence is achieved for higher damping values, cf the Rosenblatt
theorem (Grigoriu, 1995). Normalized histograms of displacement for eigenfrequencies f1 = 4
and 12 Hz and both damping values are depicted in figure 4. Other techniques can be
applied for approximations of the response outside this frequency range, e.g. memoryless
transformations of Brownian colored noise, but this lies beyond the scope of our paper. Based
on heuristic arguments and the Central Limit Theorem, we can assume that the higher the
number of active spectators, and the more complex the grandstand geometry is, the more
Gaussian the response will be.
A spectral density approximation of the forcing process for the frequency domain solution,
figure 2 (b), cannot be further simplified. This is because the FRF of the structure has sharp
peaks, and thus exact function values are needed. Any approximation employing indicator
functions in the vicinities of significant harmonics preserving variance would be inaccurate.
However, we can employ filtered white noise processes AR(2), which arise as a solution
of the second order Itô equation
¨
˙
c2,i Ŷi (t) + c3,i Ŷi (t) + c1,i Ŷi (t) = W (t)

(35)

with spectral density
si (ω) =

1
[c1,i − c2,i

ω 2 ]2

+ (ωc3,i )2

(36)

where (c1,i , c2,i , c3,i ) correspond to the stiffness, massp
and damping of a harmonic oscillator.
This function has a sharp peak positioned at f1 = c1 /c2 /2π when we neglect shifts due
P
to damping effects. For a closer approximation, we assume Ỹ (t) ≈ ni=1 Ŷi (t), where Ŷi (t)
10

Density

Density

0.4
0.2
0
−5

0
Data

0.4
0.2
0
−5

5

0
Data

5

(a) f1 = 4 Hz, ζ = 0.001 (b) f1 = 12 Hz, ζ = 0.001

Density

Density

1
0.4
0.2
0

−5

0
Data

0.5
0
−5

5

(c) f1 = 4 Hz, ζ = 0.07

0
Data

5

(d) f1 = 12 Hz, ζ = 0.07

Figure 4: Histograms of normalized displacement Z̄(t) with standard normal density based
on 1000 MC realizations.
are mutually independent AR(2) processes. Identification leads to a nonlinear optimization
problem: find such ck,i , k = 1, 2, 3 and i = 1, . . . , n, that minimize L2 norm || • ||L2 of the
difference
e(ω, c) = fˆỸ (ω) −

n
X
i=1

min

c={c1,1 ,...,c3,n }

1
[c1,i − c2,i ω 2 ]2 + (ωc3,i )2

||e(ω, c)||L2 , ω ∈ A ⊂ R+ compact,

(37)
(38)

where fˆỸ denotes a spectral density estimate of the centered force term Ỹ (t). The problem
can be solved by the Nonlinear Least Squares method, by Simulated Annealing etc, with
easily estimated initial vector c0 . Optimized coefficients for n = 6 of the centered process
Ỹ (t) are presented in table 2, and the corresponding spectral density and spectral distribution
function are presented in figure 5. Note that the spectral density is two-sided, and only
one half was integrated in the spectral distribution function, thus the variance indicated is
0.7627/2 = 0.3814.
The performance of the structure was briefly quantified in section 2.4, where the two
expressions (30) and (32) depended on the response variance. Thus an alternative approach
is to optimize the response variance directly across some eigenfrequency range of a harmonic
oscillator. Such coefficients are summarized in table 2, with indices var, spectral density and
distribution function are depicted in figure 5, frequency range 0.5 − 10 Hz.
Let us briefly note the situation when the forcing process Y (t) is not a unit process,
i.e. GH 6= 1. Since we are limiting our considerations to Gaussian approximation, only the
first two moments of GH will apply. The deterministic weight corresponds to a singular
case varGH = 0. Then the forcing term has the form YG (t) = GH Y (t), mean response
2
µZ,G (t) = E[GH ]µZ (t) and stationary response variance σZ,G
= E[G2H ]σZ2 , where µZ (t) and
11

Table 2: Coefficients ck,i , k = 1, 2, 3 and i = 1, . . . , 6 of the six independent AR(2) members
used for approximation of the centered forcing term Ỹ (t) in frequency range 0.5 − 10 Hz.
p
i
c1,i
c2,i
c3,i
c1,i /c2,i /2π
1
90.8657 0.3227 0.0148
2.67
2
35.9464 0.1276 0.1167
2.67
3 283.6701 0.2520 0.0118
5.34
4
74.1544 0.0661 0.0737
5.33
5 913.9890 1.1120 21.5186
4.56
6 228.5270 0.0907 0.1576
7.99
1var 91.0909 0.3237 0.0076
2.67
2var 40.3066 0.1430 0.0804
2.67
3var 281.0462 0.2490 0.0209
5.35
4var 83.7399 0.0746 0.0573
5.33
5var 914.0015 0.9376 21.6921
4.97
6var 228.7642 0.0908 0.1552
7.99

power [dB]

0.4
spect. density
spect. approx
var approx

0

−100
0

5

10
frequency [Hz]

15

spect. dist. funct.
spect. approx
var approx

0.2

−50

20

(a) spectral density

0

0

5

10
15
frequency [Hz]

20

(b) spectral distribution function

Figure 5: Spectral P
density and spectral distribution function of centered forcing Ỹ (t) and of
its approximation 6i=1 Ŷi (t), where Ŷi (t) are independent AR(2) processes with the coefficients in table 2 based on spectral and variance optimization.

12

10
8
6

Z

4
2
0
−2
−4
−6
−5

0

0

5
10

5
10

15
X

0.2
0
−5

0
Data

5

0.4

0.4

Density

0.4

Density

Density

Y

0.2
0
−5

0
Data

5

0.2
0
−5

0
Data

5

Figure 6: Tested structures, geometries with labeled points of interest and normed histograms
of the displacement, structures occupied by an active crowd only.
Table 3: First eigenvalues corresponding to vertical bending modes [Hz], cf figure 6.
Structure
f1
f2
f3
f4
f5
f6
f7
f8
beam
7.5 22.5 40.3 —
—
— — —
cantilever 5.4 7.0 8.2 25.1 28.4 — — —
grandstand 2.5 2.6 3.8 4.5 4.8 4.9 5.4 6.0
σZ2 are the response mean and the variance of the structure loaded by the unit forcing term
Y (t). However we should be aware that even when process Y (t) was a Gaussian, process
YG (t) as a product of a random variable with a stochastic process, is not Gaussian. To
quantify the influence of such scaling, compare the histograms in figure 2 (c) and 2 (d),
where GH has normal distribution with mean value 0.7709 and variance 0.0167.
4. Numerical examples and comparison
In this section, the quality of the approximation in the time or frequency domain will
be compared with MC simulation. A total of four mechanical systems will be tested: a
harmonic oscillator, a simply supported beam, a simple cantilever grandstand, and a realistic
grandstand with 1, 4, 72 and 630 positions for spectators, respectively. Geometries with
centered normed response histograms for an active crowd only, based on MC simulation,
are depicted in figure 6, and several lowest eigenfrequencies corresponding to the vertical
bending modes are presented in table 3. Note also that assumptions on convergence to
normal distribution are approximately fulfilled. All examples are artificial, not realistic, so
the measured responses would be unacceptable.

13

−3

displ. [m]

4

x 10

MC
Fourier series

2
0
−2

0

2

4

6

8

10

t [s]

Figure 7: Response mean displacement in comparison with an approximation based on the
first four harmonics for a harmonic oscillator, f1 = 5 Hz ζ = 0.07.
4.1. Harmonic oscillator
Let us assume a harmonic oscillator with unit mass, two values of viscous damping ζ
and variable stiffness. The mean value response is presented in figure 7. The response was
acquired by direct integration of (2), starting at t = 0.5 s. The approximation utilizes
equation (33) and the coefficients in table 1. A comparison of the total mean up-crossings
n+
x (T ) in the time interval [0, T ], T = 160 s, as functions of oscillator eigenfrequency f1 for
two values of viscous damping ζ = 0.001 and ζ = 0.07 and for two fixed levels x = 0.002
and x = 0.005 m are depicted in figure 8, employing formulas (30) and (31). The size of
time interval T is based on heuristic considerations about the average length of the musical
compositions. The time domain solution is based on the sum of six independent AR(2)
processes with the coefficients in table 2 for var optimization. The stationary response
variances are computed according to formulas (28), (29) and (15). The frequency domain
approximation employs equations (21), (22) and (23). The spectral density estimate fˆỸ (ω) of
the centered input process has 307 values over the frequency range 0 − 10 Hz, using variable
division. Figure 8 shows that the results are roughly in agreement with MC in the frequency
range 0.2 − 10 Hz. Note that the total number of mean zero-up-crossings for a deterministic
periodic function with frequency 2.67 Hz is 160 · 2.67 ≈ 427, cf figure 8 (c) and (d), where
distinct plateaux are found. The results for MC are based on 1000 realizations 160 s in
length.
4.2. A simply supported beam
The next example is a simply supported beam with Rayleigh damping ζ1 = 0.05 and
ζ2 = 0.08 for the first two vertical modes and total mass 1700 kg. Poor approximations are
anticipated, since the structure is quite stiff with high eigenfrequencies, see table 3 and the
results for the harmonic oscillator. Two cases are studied: a structure occupied by an active
crowd only; and a structure occupied by a mixed crowd. In the second case, the two left
hand side positions are loaded by forces and the two right hand side positions are occupied
by passive spectators. Deterministic biodynamic models according to Coermann are used.
These are simple-degree-of-freedom oscillators with mass 86.2 kg, stiffness 85.25 kN/m and
viscous damping 1.72 kNs/m, eigenfrequency 5 Hz, cf (Sachse et al., 2003). Results are
presented only for the labeled point in figure 6. The total mean up-crossings for the first case
n+
x (160) as a function of level x are depicted in figure 9 (a). The RM S values for acceleration
are 1.947 m/s2 for the MC solution, 2.135 m/s2 for the frequency domain solution, and
1.914 m/s2 for the time domain solution. All input processes are treated as independent, so
14

MC
freq.
time

1000

n+0.005(160)

n+0.002(160)

1500

500
0

0

2

4

6
f1 [Hz]

8

10

MC
freq.
time

1000
500
0

12

0

2

(a) ζ = 0.001

4

6
f1 [Hz]

8

10

12

(b) ζ = 0.001

MC
freq.
time

500

0

0

2

4

6
f1 [Hz]

8

10

n+0.005(160)

n+0.002(160)

1000
MC
freq.
time

400
200
0

12

0

2

(c) ζ = 0.07

6
f1 [Hz]

8

10

12

(d) ζ = 0.07

100

8

MC
freq.
time

50

RMS

RMS

4

MC
freq.
time

6
4
2

0

0

2

4

6
f1[Hz]

8

10

0

12

(e) ζ = 0.001

0

2

4

6
f1[Hz]

8

10

12

(f) ζ = 0.07

Figure 8: Total mean up-crossings n+
x (T ) and acceleration RM S values as functions of f1 for
a harmonic oscillator, distinct levels x and viscous damping ζ, T = 160 s.

15

800
n+x (160)

400

n+x (160)

MC
freq.
time

600

200
0

−5

0
x−level

5

400
200
0

10

2

4

x−level

−3

x 10

(b) 2 active, 2 passive spectators
−3

−3

x 10
displ. mean [m]

6
displ. [m]

0

x 10

x 10
4
2
0
0

−2

−3

(a) active crowd only

−2

MC
freq.
time

600

10

20

30
t [s]

40

50

MC
Fourier series

4
2
0
0

60

2

4

6

8

10

t [s]

(c) single realization

(d) response mean

Figure 9: Total mean up-crossings n+
x (T ) of a simply supported beam as functions of x,
T = 160 s for an active crowd (a) and for a mixed crowd (b), single realization for an active
crowd (c) and mean response for an active crowd (d).
matrix S Ỹ Ỹ in equation (21) has nonzero only diagonal entries. The mean value response
is depicted in figure 9 (d) with a single realization (c). The total up-crossings for a mixed
crowd are depicted in figure 9 (b), and the RM S values are 1.048 m/s2 for the MC solution,
1.094 m/s2 for the frequency domain solution, and 0.990 m/s2 for the time domain solution.
Let us also recall our assumption of fixed spatial distribution of the crowd, mass coefficient
γ = mH /mS = 0.1, where mH denotes the total mass of passive spectators, and mS denotes
the total mass of the structure. Results for MC based on 2000 realizations. The approximate
shape of the total up-crossings, especially for a mixed crowd, differ from the MC simulation,
because of the non-Gaussian response due to high structure eigenfrequencies.
4.3. Cantilever grandstand
This system has total mass 18.2 t, geometry according to figure 6, and is loaded with 72
active spectators in the first case, Rayleigh damping with ζ1 = 0.05 and ζ2 = 0.08 is used for
the first two vertical modes. The total up-crossings of the response displacement are depicted
in figure 10 (a), and single realization and the mean response are depicted in sub-figures (c)
and (d), RM S accelerations 5.583 m/s2 for the MC solution, 5.682 m/s2 for the frequency
domain solution, and 5.616 m/s2 for the time domain solution. For 36 spectators chosen to
be passive according to Coermann with uniformly random but fixed positions, the resulting
up-crossings are depicted in sub-figure (b), mass coefficient γ = 0.17. The acceleration RM S
values in this case appear to be 1.327 m/s2 for the MC solution, 1.347 m/s2 for the frequency
domain solution, and 1.309 m/s2 for the time domain solution. The results for MC are again
based on 2000 realizations 160 s in length.

16

600
MC
freq.
time

500

n+x (160)

n+x (160)

1000

−0.01

0
x−level

0.01

displ. mean [m]

0

10

20

30
t [s]

40

−2

0

2
x−level

4

6

8
−3

x 10

(b) 36 active, 36 passive spectators

0.01

0

−4

0.02

(a) active crowd only

displ. [m]

MC
freq.
time

200
0

0

−0.01

400

50

0
−0.01

60

MC
Fourier series

0.01

0

2

4

6

8

10

t [s]

(c) single realization

(d) response mean

Figure 10: Total mean up-crossings of a cantilever grandstand, active crowd (a) and mixed
crowd (b), single realization (c) and mean response (d) for an active crowd.
4.4. Realistic grandstand
In this concluding example, let us briefly examine the results acquired for a complex
structure. The geometry is sketched in figure 6 with the point of interest labeled, total mass
148.6 t, Rayleigh damping with ζ1 = 0.01 for the first and ζ2 = 0.02 for the sixth vertical
mode used here. The results for the structure loaded by an active crowd only are depicted
in figure 11 (a) (c) (d). The acceleration RM S values appear to be 5.546 m/s2 for the MC
solution, 5.565 m/s2 for the frequency domain solution, and 5.597 m/s2 for the time domain
solution. In the second case, 315 positions are occupied by passive spectators according
to Coermann, mass coefficient γ = 0.18. The results for this case are presented in figure
11 (b), RM S accelerations 1.905 m/s2 for the MC solution, 1.904 m/s2 for the frequency
domain solution, and 1.919 m/s2 for the time domain solution. MC simulation based on
2000 realizations.
4.5. Comparison and performance
The time consumption for the different solution techniques is summarized in table 4,
where the demands of MC simulation are presented for 50 realizations only. This value
is based on the convergence tests of total up-crossings presented in figures 12 (a) and (b)
for a cantilever grandstand. Obviously these tests are highly sensitive to level x, and the
value used can be considered as a lower bound. The size of the time integration step is
chosen to be h = 0.01 s, Newmark integration scheme used. The number of dofs of each
system is also mentioned, together with the size of the Lyapunov equation (14) nL for
a mixed crowd. Obviously nL = 2ndof + 2np + 12na where ndof denotes the number of
dofs of the empty structure, np is the number of dofs of the passive spectators, and na is
the number of active spectators. In the case of a cantilever grandstand we therefore have
17

400
MC
freq.
time

200

0
−0.05

n+x (160)

n+x (160)

400

0
x−level

0

0.05

displ. [m]

0.02
0
−0.02
10

20

30
t [s]

40

−0.01

0
x−level

0.01

0.02

(b) 315 active, 315 passive spectators
displ. mean [m]

(a) active crowd only

0

MC
freq.
time

200

50

MC
Fourier series

0.02
0
−0.02

60

0

(c) single realization

5

10
t [s]

15

20

(d) response mean

Figure 11: Total mean up-crossings of a realistic grandstand, active crowd (a), mixed crowd
(b), single realization (c) and mean response (d) for an active crowd.

n+x0.006(160)

n+x0.005(160)

300
250
200

60
50
40
30

0

500

1000
1500
MC realizations

2000

0

(a) x = 0.005

500

1000
1500
MC realizations

2000

(b) x = 0.006

Figure 12: MC convergence tests of n+
x (160) for a cantilever grandstand, cf figure 10 (b).
2 · 504 + 2 · 36 + 12 · 36 = 1512, since passive spectators are modeled as systems with a single
degree of freedom. For the purposes of comparison, the time consumption for the partial
modal transform are also summarized together with the number of eigenvectors used neig
and according highest frequency feig . All simulations were performed on core i7 with a 16 GB
RAM computer, Matlab R parallel implementation. The results in the table suggest that
when 50 realizations are employed, only structures of moderate size can be effectively solved
by semi-analytical methods.
5. Conclusions
This paper has presented a study of the vibration of grandstands loaded by an active
crowd, using Gaussian approximation of the response. The main results can be summarized
as follows:
1. A mathematical description of the response of a mechanical system employing spectral
18

Table 4: Comparison of computational demands.
Method/Struct. Harm. osc.

Beam Cant. grand. Real. grand.
full system
ndof
1
29
504
4 068
nL
14
86
1512
12 546
MC 50
0.601 s 20.078 s
380.451 s
5 484 s
Freq. domain
0.106 s 2.672 s
74.639 s
4 178 s
Time domain
0.112 s 2.392 s
38.915 s
6 507 s
partial modal transform, cf equations (27), (28) and (29)
neig
–
7
10
50
feig
–
37 Hz
30 Hz
23 Hz
nL
–
42
524
4 510
MC 50
– 13.240 s
36.214 s
307 s
Freq. domain
– 1.764 s
4.452 s
734 s
Time domain
– 1.588 s
2.582 s
67 s

2.

3.

4.

5.

and time domain solutions for weakly stationary Gaussian excitations has been recalled.
Partial modal transformation due to a passive crowd has been briefly discussed.
A motivating example of a harmonic oscillator has shown that the normalized displacement and velocity have approximately normal distribution under the conditions
on eigenfrequencies and damping.
Taking this fact into account, the mean value of the force has been approximated as
a truncated Fourier series, and the spectral density of the centered process has been
estimated employing the Parzen window for the frequency domain solution. For the
time domain solution, we have employed a linear combination of independent autoregression processes of the second order with coefficients optimized to achieve the least
error in the response variance of a harmonic oscillator.
Three different examples of varying complexity have shown the quality of the response
approximation in terms of total displacement up-crossings and acceleration RM S in
comparison with Monte Carlo simulation. Limitations following from a simple oscillator
have been confirmed on multi-degree-of-freedom systems.
The computational demands have been measured and summarized in terms of the time
needed for solution, and the applicability of the techniques has been proved.

Finally, let us note that these methods are approximations and can be further refined. In
particular, the solution can be reformulated for the non-Gaussian processes for which higher
moments can be derived. These apply in improving the mean up-crossing rate estimates for
stiff structures.
Acknowledgement
Financial support for this work from the Czech Technical University in Prague under
project No. SGS12/027/OHK1/1T/11 and from the Czech Science Foundation under project
19

No. GAP105/11/1529 is gratefully acknowledged.
References
Anděl, J., 1976. The Statistical Analysis of Time Series (in Czech). SNTL.
Brockwell, P. J., Davis, R. A., Yang, Y., 2007. Continuous-time gaussian autoregression.
Statistica Sinica 17, 63–80.
Ellis, B. R., Ji, T., 1994. Floor vibration induced by dance-type loads: Verification. The
Structural Engineer 72 (3), 45–50.
Grigoriu, M., 1995. Applied Non-Gaussian Processes: Examples, Theory, Simulation, Linear
Random Vibration, and MATLAB Solutions. Prentice Hall.
Grigoriu, M., 2002. Stochastic Calculus:
Birkhäuser.

Applications in Science and Engineering.

Hermanussen, M., Danker-Hopfe, H., Weber, G. W., 2001. Body weight and the shape of the
natural distribution of weight, in very large samples of german, austrian and norwegian
conscripts. International Journal of Obesity 25, 1550–1553.
Huang, Y., Griffin, M. J., 2008. Nonlinear dual-axis biodynamic response of the semi-supine
human body during longitudinal horizontal whole-body vibration. Journal of Sound and
Vibration 312, 273–295.
Racic, V., Pavic, A., 2010. Mathematical model to generate near-periodic human jumping
force signals. Mechanical Systems and Signal Processing 24, 138–152.
Sachse, R., Pavic, A., Reynolds, P., 2003. Human-structure dynamic interaction in civil
engineering dynamics: A literature review. The Shock and Vibration Digest 35, 3–18.
Sim, J. H., 2006. Human structure interaction in cantilever grandstands.
Soong, T. T., Grigoriu, M., 1993. Random Vibration of Mechanical and Structural Systems.
Prentice Hall, New Jersey.

20

