Grey-box nonlinear state-space modelling
for mechanical vibrations identification

arXiv:1610.09145v1 [] 28 Oct 2016

J.P. Noël ∗,∗∗ J. Schoukens ∗ G. Kerschen ∗∗
∗
Department ELEC
Vrije Universiteit Brussel, Brussels, Belgium
e-mail: jp.noel@ulg.ac.be, johan.schoukens@vub.ac.be
∗∗
Space Systems and Structures Laboratory (S3L)
Aerospace and Mechanical Engineering Department
University of Liège, Liège, Belgium
e-mail: g.kerschen@ulg.ac.be

Abstract: In the present paper, a flexible and parsimonious model of the vibrations of nonlinear
mechanical systems is introduced in the form of state-space equations. It is shown that the
nonlinear model terms can be formed using a limited number of output measurements. A twostep identification procedure is derived for this grey-box model, integrating nonlinear subspace
initialisation and maximum likelihood optimisation. The complete procedure is demonstrated on
the Silverbox benchmark, which is an electrical mimicry of a single-degree-of-freedom mechanical
system with one displacement-dependent nonlinearity.
Keywords: Nonlinear system identification; mechanical systems; state-space modelling;
nonlinear subspace initialisation; maximum likelihood optimisation; Silverbox benchmark.
1. INTRODUCTION
Nonlinear system identification constantly faces the compromise between the flexibility of the fitted model and its
parsimony. Flexibility refers to the ability of the model
to capture complex nonlinearities, while parsimony is its
quality to possess a low number of parameters. In this
regard, a nonlinear state-space representation

ẋ(t) = A x(t) + B u(t) + E g(x, u)
(1)
y(t) = C x(t) + D u(t) + F h(x, u)
can be classified as very flexible but little parsimonious,
two features typically shared by black-box models. In
Eqs. (1), A ∈ R ns ×ns , B ∈ R ns ×m , C ∈ R l×ns and
D ∈ R l×m are the linear state, input, output and direct
feedthrough matrices, respectively; x(t) ∈ R ns is the state
vector; y(t) ∈ R l and u(t) ∈ R m are the output and
input vectors, respectively. The linear-in-the-parameters
expressions E g(x, u) ∈ R ns and F h(x, u) ∈ R l are
the nonlinear model terms coupling the state and input
variables. The order of the model, i.e. the dimension of
the state space, is noted ns .
In the present paper, it is shown that, in the case of
mechanical systems where nonlinearities are physically
localised, the model structure in Eqs. (1) can be drastically
simplified. More specifically, Section 2 demonstrates that
the nonlinear terms in Eqs. (1) can be constructed using a
limited number of output measurements, and so without
resorting to the state and input vectors. This makes the
resulting grey-box state-space model a parsimonious representation of nonlinear mechanical systems. An efficient,
two-step identification procedure for this model is derived
in Section 3, integrating nonlinear subspace initialisation

and maximum likelihood optimisation. Finally, the complete procedure is applied in Section 4 to an electrical circuit mimicking the behaviour of a single-degree-of-freedom
mechanical system with one displacement-dependent nonlinearity.
2. GREY-BOX NONLINEAR STATE-SPACE
MODELLING OF MECHANICAL VIBRATIONS
In the analysis of mechanical vibrations, one very often
distinguishes localised nonlinearities, which are physically
confined to a small area, from nonlinearities distributed
throughout (some large region of) the entire structure.
Localised elements are arguably the most common in
engineering practice, since structural nonlinearities typically arise from the complex dynamics of joints interfacing
subcomponents. Many meaningful examples of this reality
are to be found in the aerospace sector. For instance,
nonlinearities resulting from the appearance of gaps in
the truss supports of the Huygens probe were attested
during the modal survey of the Cassini spacecraft [Carney
et al. (1997)]. Nonlinearities were also reported during
ground vibration testing of the Airbus A400M, and were
attributed to the elastomeric mounts supporting the four
turboprop engines of the aircraft [Ahlquist et al. (2010)].
Moreover, the analysis of in-orbit data of the International
Space Station highlighted that the opening of a pin connection in the assembly of its solar arrays led to severe
nonlinearity [Laible et al. (2013)].
Assuming localised nonlinearities, the vibrations of a np degree-of-freedom mechanical system obey Newton’s second law written in the form

M q̈(t)+Cv q̇(t)+Kq(t)+

s
X

ca ga (qnl (t), q̇nl (t)) = p(t),

a=1

(2)
where M, Cv , K ∈ R np ×np are the mass, linear viscous
damping and linear stiffness matrices, respectively; q(t)
and p(t) ∈ R np are the generalised displacement and
external force vectors, respectively; the nonlinear restoring force term is written using s basis function vectors
ga (t) ∈ R np associated with coefficients ca . The subset
of generalised displacements and velocities involved in
the construction of the basis functions are denoted qnl (t)
and q̇nl (t), respectively. Physically, they correspond to
mechanical degrees of freedom located on both sides of
the localised nonlinearities in the system.

where B = [B E] and D = [D F]; the extended input

T
vector u(t) is similarly defined as u(t)T g(t)T , where
T is the transpose operation.
3. IDENTIFICATION PROCEDURE

The dynamics governed by Eq. (2) is conveniently interpreted by moving the nonlinear restoring force term to the
right-hand side, i.e.
s
X
M q̈(t)+Cv q̇(t)+Kq(t) = p(t)−
ca ga (qnl (t), q̇nl (t)),

The simplified, linear-like structure of the grey-box statespace model given in Eqs. (5) lends itself to an efficient,
two-step identification procedure. First, initial estimates
of the A, B, C and D matrices are obtained using
a nonlinear subspace identification method. Second, the
subspace estimates are optimised in maximum likelihood
sense by applying a nonlinear minimisation routine. The
complete procedure is carried out in the frequency domain, which opens the possibility to apply user-defined
weighting functions in specific frequency bands. For the
sake of comparison, the identification of the black-box
polynomial nonlinear state-space model in Ref. [Paduart
et al. (2010)] requires a four-step process, including linear
subspace parameter estimation and two nonlinear optimisation searches.

which leads to the block diagram in Fig. 1.

3.1 Initialisation using a nonlinear subspace method

a=1

(3)

p(t)
+

Underlying linear

q(t), q̇(t)

system: M, Cv , K

Nonlinear feedback:
ca , ga (qnl (t), q̇nl (t))

Fig. 1. Feedback interpretation of Newton’s law in Eq. (2).
The feedback structure of this diagram suggests that localised nonlinearities in mechanical systems act as additional inputs applied to the underlying linear system. This,
in turn, reveals that black-box nonlinear terms in a statespace model, such as E g(x, u) and F h(x, u) in Eqs. (1),
are overly complex to address mechanical vibrations. A
more parsimonious description of nonlinearities is achieved
by translating Eq. (3) in state space, which provides the
grey-box model

ẋ(t) = A x(t) + B u(t) + E g(ynl (t), ẏnl (t))
(4)
y(t) = C x(t) + D u(t) + F g(ynl (t), ẏnl (t)),
where g(t) ∈ R s is a vector concatenating the nonzero elements in the basis function vectors ga (t), and E ∈ R np ×s
and F ∈ R l×s are the associated coefficient matrices;
ynl (t) and ẏnl (t) are the subsets of the measured displacements and velocities located close to nonlinearities,
respectively.
For the sake of conciseness, one adopts the concatenated
equations

ẋ(t) = A x(t) + B u(t)
(5)
y(t) = C x(t) + D u(t),

Subspace algorithms are well known for solving linear
system identification problems [Van Overschee and De
Moor (1996); McKelvey et al. (1996)]. In recent years,
nonlinear generalisations of subspace identification have
also emerged, following the original idea of Lacy and
Bernstein [Lacy and Bernstein (2005)]. The present study
exploits the frequency-domain nonlinear subspace method
proposed in Ref. [Noël and Kerschen (2013)], termed
FNSI method,
to calculate an initial set of matrices

A, B, C, D . A consistent estimate of the set is obtained
if both the input u(t) and the nonlinear basis functions
in g(ynl (t), ẏnl (t)) are noiseless, as inferred from the results in Refs. [McKelvey et al. (1996); Pintelon (2002)].
The use of FNSI offers two important advantages. Firstly,
the method constructs a fully nonlinear model of the
system under test from the beginning, in contrast with
the linearised model which serves as a starting point in
Ref. [Paduart et al. (2010)]. Secondly, the method embeds
an intuitive model order selection capability via stabilisation diagrams [Noël et al. (2014)].
3.2 Nonlinear optimisation of the subspace model in the
maximum likelihood framework
In the case of output measurements with a low signal-tonoise ratio, consistency of the initial state-space matrices
A, B, C, D is lost. Unbiased parameter estimates can be
recovered by optimising the subspace model in maximum
likelihood (ML) sense. The ML framework also guarantees
the lowest possible uncertainty on the model parameters,
i.e. the efficiency of the estimates [Pintelon and Schoukens
(2001)]. To formulate the ML cost function, the state-space
Eqs. (5) are recast in the frequency domain as

ξk X(k) = A X(k) + B U(k)
(6)
Y(k) = C X(k) + D U(k),
where k is the frequency line, ξk the Laplace and ztransform variable depending on whether a continuous- or
discrete-time formulation is selected, and Y(k), X(k) and


∂y(t)
∂
=
C x(t) + D u(t)
∂Ai j
∂Ai j
∂x(t)
∂u(t)
= C
+D
∂Ai j
∂Ai j
∂u(t) ∂y(t)
∂x(t)
+D
.
= C
∂Ai j
∂y(t) ∂Ai j

U(k) the discrete Fourier transforms (DFTs) of y(t), x(t)
and u(t), respectively.
Assumption on the noise model
The extended input spectrum is assumed to be noiseless, i.e. observed without errors and independent of the
output noise. The output disturbing noise term NY (k)
is Gaussian distributed, has zero mean E (NY (k)) = 0,
where E is the expectation operator, and has a covariance matrix with only nonzero diagonal elements equal to
2
σ 2Y (k) = E |NY (k)| , as described in Ref. [Schoukens
and Renneboog (1986)].
Maximum likelihood cost function
The Gaussianity of the output noise implies that the ML
cost function simplifies to a weighted least-squares estimator in the frequency domain [Pintelon and Schoukens
(2001)]. Introducing the vector of model parameters θ as



θ = vec (A) ; vec B ; vec (C) ; vec D ,
(7)
where the operation denoted vec stacks the columns of a
matrix on top of each other, the cost function to minimise
hence writes
F
X
V(θ) =
ǫH (k, θ) W2 (k) ǫ(k, θ),
(8)
k=1

where F is the number of processed lines, H the Hermitian
transpose, and W(k) a weighting function chosen equal to
σ Y−1 (k). The model error vector ǫ ∈ R l is defined as the
complex-valued difference
ǫ(k, θ) = Ym (k, θ) − Y(k),

The first term in the right-hand side of Eq. (12) is obtained
by taking the derivative of the state relation in Eqs. (5)
with respect to Ai j , that is

∂ ẋ(t)
∂
=
A x(t) + B u(t)
∂Ai j
∂Ai j
∂x(t)
=A
+ Iinjs ×ns x(t) + B
∂Ai j
∂x(t)
+ Iinjs ×ns x(t) + B
=A
∂Ai j

Analytical calculation of the Jacobian matrix
In practice, the nonlinear least-squares cost function in
Eq. (8) is minimised using a Levenberg-Marquardt optimisation algorithm, which combines the large convergence
region of the gradient descent method with the fast convergence of the Gauss-Newton method [Levenberg (1944);
Marquardt (1963)]. This algorithm requires the calculation
of the Jacobian matrix J(k, θ) associated with the cost
function or, equivalently, with the error function in Eq. (9),
i.e.
∂ǫ(k, θ)
∂Ym (k, θ)
J(k, θ) =
=
.
(10)
∂θ
∂θ
Given the nonlinear relationship which exists between
Y(k) and U(k), it may not be practical to compute the
elements of J(k, θ) directly in the frequency domain. An
alternative approach consists in carrying out the computation of the Jacobian matrix in the time domain, and
then in applying the DFT. One first focuses on the determination of the element JAij (t) ∈ R l of the time-domain
Jacobian defined as
∂y(t)
.
(11)
JAij (t) =
∂Ai j
The derivative of the output relation in Eqs. (5) with
respect to Ai j is given by

∂u(t)
(13)
∂Ai j
∂u(t) ∂y(t)
,
∂y(t) ∂Ai j

where Iinjs ×ns is a zero matrix with a single element equal
to one at entry (i, j).
The element JAij (t) is therefore given by the solution of
the two equations

∂ ẋ(t)
∂u(t) ∂y(t)
∂x(t)


= A
+ Iinjs ×ns x(t) + B

∂Ai j
∂Ai j
∂y(t) ∂Ai j
∂u(t) ∂y(t)
∂x(t)
∂y(t)


= C
+D
.

∂Ai j
∂Ai j
∂y(t) ∂Ai j
(14)
Introducing the notations
x∗ (t) =

(9)

where Ym (k, θ) and Y(k) are the DFTs of the modelled
and measured outputs, respectively.

(12)

u∗ (t) =



∂x(t)
∂y(t)
; y∗ (t) =
;
∂Ai j
∂Ai j
T

x(t)



∂u(t) ∂y(t)
∂y(t) ∂Ai j

T T

(15)

and

∗
A∗ = A ; B = Iinjs ×ns B ;

∗
C∗ = C ; D = 0 l×ns D ,

Eqs. (14) can be recast in the form
 ∗
∗
ẋ (t) = A∗ x∗ (t) + B u∗ (t)
∗
y∗ (t) = C∗ x∗ (t) + D u∗ (t).

(16)

(17)

Eqs. (17) reveal that the elements of the Jacobian matrix
associated with the parameters in A are solutions of an
auxiliary
state-space
model defined by the four matrices


∗
∗
∗
∗
A , B , C , D . The first term in the auxiliary extended
input u∗ (t) in Eq. (15) is the state vector x(t). It is obtained by simulating in time the original model in Eqs. (5)
with the estimated parameters of the previous LevenbergMarquardt iteration. The second term in u∗ (t) depends on
∂u(t)/∂y(t), which is formed using the derivatives of the
nonlinear basis functions g(ynl (t), ẏnl (t)) with respect to
y(t).
The determination of the element JB ij (t) ∈ R l is conducted similarly to JAij (t). The result is given in Eqs. (18),
where JB ij (t) is seen to be the solution of another auxiliary
state-space model,


∂ ẋ(t)
∂u(t) ∂y(t)
∂x(t)
n ×(m+sl)


u(t) + B
= A
+ Ii js


∂y(t) ∂B i j
 ∂B i j
∂B i j

The computation of JCij (t) ∈ R and JDij (t) ∈ R l is
easier because they do not involve time integration, as
shown in Eq. (19) and Eq. (20), respectively,
∂y(t)
∂u(t) ∂y(t)
s
x(t) + D
= Iil×n
;
j
∂Ci j
∂y(t) ∂Ci j

(19)

∂u(t) ∂y(t)
∂y(t)
l×(m+sl)
u(t) + D
= Ii j
.
∂y(t) ∂Di j
∂Di j

(20)

4. EXPERIMENTAL DEMONSTRATION
ON THE SILVERBOX BENCHMARK
The identification procedure described in Section 3 is
demonstrated herein using experimental measurements
acquired on the Silverbox circuit mimicking the behaviour
of a single-degree-of-freedom nonlinear mechanical system.
Ideally, this system should exhibit the dynamics of a
Duffing oscillator with cubic nonlinearity, as prescribed by
the equation
M q̈(t) + Cv q̇(t) + K q(t) + c1 q 3 (t) = p(t).
(21)
In practice, it is also known to be characterised by, at least,
an additional quadratic stiffness term c2 q 2 (t). The system was excited using random phase multisines [Pintelon
and Schoukens (2001)] considering equivalent root-meansquared (RMS) amplitudes of 5 and 150 mN . The input
frequency spectrum was limited to 0 – 300 Hz, excluding
the DC component, with a sampling frequency of 2441 Hz.
Experiments were conducted over 25 periods of 8192 samples, removing the first 5 periods to achieve steady-state
conditions. Table 1 reports the underlying linear modal
properties of the benchmark estimated using a subspace
analysis at 5 mN RMS. Fig. 2 depicts the comparison
between frequency response functions (FRFs) measured
at 5 and 150 mN RMS. The two curves reveal that the
Silverbox vibrates in a strongly nonlinear regime of motion
at high level, as a shift of the resonance frequency of more
than 13 Hz is noticed at 150 mN RMS together with
severe noisy-like distortions.
Natural frequency (Hz)
68.57

Amplitude (dB)

(18)

l

10

−5

−20
0

50
100
Frequency (Hz)

150

Fig. 2. Comparison of FRFs measured at 5 (dashed line)
and 150 (solid line) mN RMS.
and linear terms only in the output equation. This led to
a nonlinear model with 37 parameters.
The time- and frequency-domain errors associated with
the two steps of the identification procedure are plotted in
Figs. 3 and 4, respectively. In the time domain in Fig. 3,
the RMS error of the initial subspace model is equal to 6.24
10−7 m, compared to the signal RMS value of 10.54 10−7
m. The error is decreased down to 0.40 10−7 m after ML
optimisation. The analysis of the reconstructed spectra in
Fig. 4 is also interesting. The error of the final model is
generally 30 dB below the measured output spectrum and
20 dB lower than the subspace model. However, it does
not reach the noise level, which is most probably due to an
imperfect representation of the nonlinearity in the system.
Errors in the final model are particularly visible in Fig. 4
at the resonance location around 83 Hz, and close to third
harmonics of the system around 250 Hz, proving that they
are related to the modelling of the nonlinearity.
5

Displacement (m)



∂y(t)
∂x(t)
∂u(t) ∂y(t)


= C
+D
.

∂y(t) ∂B i j
∂B i j
∂B i j

25

x 10

−6

0

Damping ratio (%)
4.68

Table 1. Natural frequency and damping ratio
of the Silverbox estimated at 5 mN RMS.
A state-space model in the grey-box form of Eqs. (5) is
constructed using as nonlinear terms quadratic and cubic
functions of the measured output displacement ynl (t). A
model order equal to 2 is obviously selected and the resulting vector θ in Eq. (7) thus consists of 15 parameters.
Comparatively, in Ref. [Paduart et al. (2010)], a blackbox, second-order, state-space model, as in Eqs. (1), was
adopted considering a third-degree multivariate polynomial in the state equation with all cross products included

−5
0

1

2
Time (s)

3

Fig. 3. Time-domain errors: measured response (black);
initial model error (orange); final model error (blue).
The final estimates of the state-space parameters can
be converted into estimates of the coefficients c1 and c2
associated with the cubic and quadratic basis functions
in the model, respectively. This is achieved using the
conversion scheme proposed in Ref. [Marchesiello and
Garibaldi (2008)], which yields frequency-dependent and

20

Amplitude (dB)

−20

−60

−100
0

100

200

300

Frequency (Hz)

Fig. 4. Frequency-domain errors over 0 – 300 Hz: measured spectrum (black cross); initial subspace model error (orange
plus); final ML model error (blue circle); noise level (grey square).
0.007

3

Imaginary part (N/m )

Real part (N/m 3)

4.4

4

3.6
0

100
200
Frequency (Hz)

300

−0.007

−0.021

−0.035
0

1

x 10

100
200
Frequency (Hz)

300

−3

Imaginary part (N/m 2)

2

Real part (N/m )

−0.2

−0.25

−0.3
0

100
200
Frequency (Hz)

300

−1

−3
0

100
200
Frequency (Hz)

300

Fig. 5. Complex and frequency-dependent nonlinear coefficients c1 (top left – right) and c2 (bottom left – right) obtained
by converting the initial (in orange) and final (in blue) state-space model parameters.

complex-valued coefficients. As a result of the consistency
of the identification procedure of Section 3, and in the
absence of modelling errors, one expects the real parts of
the coefficients to converge asymptotically to their exact
values with no frequency dependence, and the imaginary
parts to converge accordingly to zero. This makes the
significance of the frequency variations and imaginary
parts of the coefficients a particularly convenient means
to assess the quality of the identification results.
The estimated nonlinear coefficients c1 and c2 are displayed versus frequency in Fig. 5 (a – d). Table 2 lists
the frequency averages of their real parts and the ratios
between real and imaginary parts in logarithmic scaling.
The real parts of the coefficients given by the initial and
final models are almost equal. They are also found to be
satisfactorily stable versus frequency, and remain more
than 2 orders of magnitude greater than the corresponding
imaginary parts. It should however be noted that the
substantial decrease of the error from the initial to the
final model in Figs. 3 and 4 does not translate into noticeably improved estimates of the nonlinear coefficients. This
deserves more investigation to precisely understand the
relation between the quality of physical-space parameters
and the overall quality of the grey-box state-space model.
Real part c1 (N/m3 )
Log10 (real/imag.)
Real part c2 (N/m2 )
Log10 (real/imag.)

Initial subspace model
3.95
2.71
-0.25
3.38

Final model
3.93
2.41
-0.25
2.54

Table 2. Estimates of the nonlinear coefficients
c1 and c2 obtained by converting the initial and
final state-space model parameters.
5. CONCLUSION
The objective of the present paper was to propose a
grey-box state-space modelling framework to support the
identification of nonlinear mechanical vibrations. It was
shown that this framework paves the way for an important
decrease of the number of model parameters with respect
to classical black-box state-space modelling. The Silverbox
benchmark was considered as an experimental case study
demonstrating the derived identification procedure, which
combines nonlinear subspace parameter initialisation and
maximum likelihood optimisation.
Additional work should focus on studying the sensitivity
of the final parameter estimates to the amplitude of excitation and to the quality of the initial subspace model. The
convergence of the parameters throughout the maximum
likelihood iterations should also be analysed in more details. More advanced research prospects include the introduction of a spline-based representation of nonlinearities
in the grey-box framework and the calculation of reliable
confidence bounds on the model parameters.
ACKNOWLEDGEMENTS
The author J.P. Noël is a Postdoctoral Researcher of the
Fonds de la Recherche Scientifique – FNRS which is gratefully acknowledged. This work was also supported in part
by the Fund for Scientific Research (FWO-Vlaanderen),

by the Flemish Government (Methusalem), by the Belgian
Government through the Inter university Poles of Attraction (IAP VII) Program, and by the ERC advanced grant
SNLSID, under contract 320378.
REFERENCES
Ahlquist, J., Carreño, J., Climent, H., de Diego, R., and
de Alba, J. (2010). Assessment of nonlinear structural
response in A400M GVT. In Proceedings of the 28th
International Modal Analysis Conference (IMAC). Jacksonville, FL.
Carney, K., Yunis, I., Smith, K., and Peng, C. (1997).
Nonlinear dynamic behavior in the Cassini spacecraft
modal survey. In Proceedings of the 15th International
Modal Analysis Conference (IMAC). Orlando, FL.
Lacy, S. and Bernstein, D. (2005). Subspace identification for non-linear systems with measured-input nonlinearities. International Journal of Control, 78, 906–
926.
Laible, M., Fitzpatrick, K., and Grygier, M. (2013). International Space Station 2A array modal analysis. In
Proceedings of the 31st International Modal Analysis
Conference (IMAC). Garden Grove, CA.
Levenberg, K. (1944). A method for the solution of
certain problems in least squares. Quarterly of Applied
Mathematics, 2, 164–168.
Marchesiello, S. and Garibaldi, L. (2008). A time domain
approach for identifying nonlinear vibrating structures
by subspace methods. Mechanical Systems and Signal
Processing, 22, 81–101.
Marquardt, D. (1963). An algorithm for least-squares
estimation of nonlinear parameters. Journal of the
Society for Industrial and Applied Mathematics, 11(2),
431–441.
McKelvey, T., Akçay, H., and Ljung, L. (1996). Subspacebased multivariable system identification from frequency response data. IEEE Transactions on Automatic
Control, 41(7), 960–979.
Noël, J. and Kerschen, G. (2013). Frequency-domain subspace identification for nonlinear mechanical systems.
Mechanical Systems and Signal Processing, 40, 701–717.
Noël, J., Marchesiello, S., and Kerschen, G. (2014).
Subspace-based identification of a nonlinear spacecraft
in the time and frequency domains. Mechanical Systems
and Signal Processing, 43, 217–236.
Paduart, J., Lauwers, L., Swevers, J., Smolders, K.,
Schoukens, J., and Pintelon, R. (2010). Identification
of nonlinear systems using Polynomial Nonlinear State
Space models. Automatica, 46, 647–656.
Pintelon, R. (2002). Frequency-domain subspace system
identification using non-parametric noise models. Automatica, 38, 1295–1311.
Pintelon, R. and Schoukens, J. (2001). System Identification: A Frequency Domain Approach. IEEE Press,
Piscataway, NJ.
Schoukens, J. and Renneboog, J. (1986). Modeling the
noise influence on the Fourier coefficients after a discrete
Fourier transform. IEEE Transactions on Instrumentation and Measurement, 35(3), 279–286.
Van Overschee, P. and De Moor, B. (1996). Subspace
Identification for Linear Systems: Theory, Implementation and Applications. Kluwer Academic Publishers,
Dordrecht, The Netherlands.

