Journal of Mechanical and Electrical Intelligent System (JMEIS)

Linear system identification of 1-dof vibratory system based on
the Maximum Likelihood Estimation using the analyticalsolution of Fokker-Planck equation
Soichiro Takata1,a,*
1Department of Mechanical Engineering, National Institute of Technology, Tokyo College, 1220-2

Kunugida-machi, Hachioji-shi, 193-0997, Japan*Corresponding author
a<takata@tokyo-ct.ac.jp>

Keywords: random vibration, maximum likelihood estimation, Fokker-Planck equation

Abstract. This paper discusses a new identification method of linear single-degree-of-freedom systemusing Gaussian random vibration response. The propose method is based on the method of Maximum
Likelihood Estimation (MLE). The likelihood function of proposed method is composed from theanalytical solution of Fokker-Planck equation. The estimation formulas of unknown parameter are
obtained by maximization of the original likelihood function. The obtained estimators correspond withpopulation variance estimation of multivariate Gaussian model. Furthermore, the numerical
identifications are conducted using the random vibration response by calculation result of 4th RungeKutta method. In the result, the estimation performance of the propose method is confirmed in termsof dependency of sample number and dependency of damping coefficient. Especially, the proposed
method is implied the application to identification problem of large damping system. Quantificationof large damping characteristic is important problem, because it is very difficult problem inconventional identification method. Moreover, the benchmark tests are conducted with Half-Power
Method (HPM) based on the spectral analysis and Auto-Regressive Method (ARM) based on the timeseries analysis, respectively. The results of benchmark are shown in the accuracy of propose method
is higher than its of HPM and ARM, respectively. Finally, the expansion to recursive estimation
algorithm is conducted using MLE estimator of recurrence form. In addition, operation of the recursivealgorithm is confirmed.
1. Introduction
The system identification method using random vibration response is traditionally used in the fieldof mechanical structural vibration. The method using the auto-regressive time series analysis is themost widely known method [1,2,3,4]. Because it only needs output data and obtaining the input datais generally difficult in field test, Auto-Regressive Method (ARM) is often used in actual field data.
However, because of its dependency on output data, its accuracy in variance of prediction error is low.
In recent years, the ARX model which uses input-output data solves the above problem [5,6].
Operational modal analysis (OMA) is recently developing in the field of stochastic signal
processing technologies [7,8,9,10]. In OMA, the mode shape visualization and the estimation of modalparameters are conducted using frequency response based on modal expansion, which uses white noise
excitation as input. The modal shape in actual structure can then be understood using OMA. It isavailable to data analysis in actual field sensor data. However, because it is based on the frequency
method, OMA cannot be applied to large damping structural and material systems, especially inresonance spectrum [11].
In this paper, we propose an identification method based on the probability density function. The
shape of the probability density function obeys system parameters such as spring constant, damping
constant, and diffusion coefficient of input white noise. The characteristics are unrelated to frequency
spectrum and correlation function. Thus, a method based on the probability density function is
52
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
expected to apply to large damping vibratory system. In linear and weakly nonlinear systems that are
subjected to white random noise excitation, the analytical solution of Fokker-Planck equation can also
be obtained [12,13,14]. In this study, we propose a new function using the analytical solution of
Fokker-Planck equation. Moreover, the unknown parameter estimation formula has been derived using
Maximum Likelihood Estimation (MLE).
2. Derivation of the Identification Algorithm
In this chapter, the proposed identification algorithms are derived based on the method of MLE.
The equation of motion in single-degree-of-freedom (dof) system in case of unit mass is as follows:𝑑 𝑥(𝑡)𝑑𝑥(𝑡)+𝑐+ 𝑘𝑥(𝑡) = 𝑤(𝑡)
𝑑𝑡
𝑑𝑡

(1)

Here, 𝑤(𝑡) denotes the white noise, 𝑁(0, 𝜎 ), whereas, 𝜎 represents the variance of input
white noise. Its probabilistic differential equation is shown in the Eq. (2) with state variables, 𝑥 =𝑥(𝑡) and 𝑥 (𝑡) = 𝑑𝑥 (𝑡)/𝑑𝑡.
𝑑 𝑥
𝑑𝑡 𝑥

𝑥
= −𝑘𝑥 − 𝑐𝑥

+

0
𝑤(𝑡)
1

(2)

Using the above equations, the Fokker-Planck equation is as follows:𝜕𝑓(𝑥 , 𝑥 , 𝑡)𝜕𝑓(𝑥 , 𝑥 , 𝑡)𝜕𝑓(𝑥 , 𝑥 , 𝑡)
𝜕 𝑓(𝑥 , 𝑥 , 𝑡)
= −𝑥
+ (𝑘𝑥 + 𝑐𝑥 )+ 𝑐𝑓(𝑥 , 𝑥 , 𝑡) + 𝐷
(3)
𝜕𝑡
𝜕𝑥
𝜕𝑥
𝜕𝑥
Here, 𝑓(𝑥 , 𝑥 , 𝑡) and 𝐷 represent the probability density function of the stochastic response andthe diffusion coefficient, respectively.
The analytical solution of stationary Fokker-Planck equation is as follows:𝑓 (𝑥 , 𝑥 ) =

𝑐√𝑘
𝑐
exp −
(𝑘𝑥 + 𝑥 )
2𝜋𝐷
2𝐷

(4)

Eq. (4) contains the parameters 𝑘 and 𝐷/𝑐 and the stationary distribution is represented by𝑓 (𝑥 , 𝑥 ).
In this paper, we will propose an estimation method of unknown parameters, 𝑘 and 𝐷/𝑐, basedon MLE. Here, the likelihood function is defined by the analytical solution of the Fokker-Planck
defined in Eq. (4) and is observed as time series data. The observed time series data contain thedisplacement,𝑥,andvelocity,
𝑥
;
data
sets
are
𝐃=
𝑥 , , 𝑥 , , 𝑥 , , 𝑥 , , … , 𝑥 , , 𝑥 , … , 𝑥 , , 𝑥 , . Therefore, the likelihood function is defined
as follows:
𝐿=

𝑐√𝑘
𝑐
exp −
(𝑘𝑥
2𝜋𝐷
2𝐷

+𝑥 ) =

𝑐√𝑘
2𝜋𝐷

exp −

𝑐
2𝐷

(𝑘𝑥

+𝑥 )

(5)

Using MLE, the negative log-likelihood function is assumed as follows:

53
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
𝐽 = − log 𝐿 = − log

𝑐√𝑘
2𝜋𝐷

+

𝑐
2𝐷

(𝑘𝑥

+𝑥 )

(6)

In the process of MLE, the negative log-likelihood function is partially differentiated using the
unknown parameters, 𝑘 and 𝐷/𝑐. Moreover, because the case of differentiated likelihood is 0, themaximum likelihood estimators are obtained as follows:
𝑘

=

∑
∑

𝑥
𝑥

(7)

𝐷
𝑐

=

1
𝑁

𝑥

(8)

Eq. (7) and Eq. (8) are equal to the variances of analytical solution with respect to displacement,
𝑥 , and velocity, 𝑥 , in the stationary Fokker-Planck equation. Therefore, our proposed method isessentially equal to the variance estimation problem.
3. Numerical Simulation
In this chapter, we will verify operation of the proposed identification method by using thenumerical consideration. In addition, the proposed method is compared with the Half Power Method(HPM) that is based on the frequency response method, and ARM that is based on the time series
analysis.
3.1 Conditions
To identify the unknown system parameters, the time series is generated using the method of 4-thorder Runge-Kutta. In HPM, the half spectrum is calculated by Fast Fourier Transform (FFT, there is
no window function when the number of samples is 8192), associated with auto-correlation functionof displacement response of 1-dof subjected to white noise excitation. Furthermore, the spring constantis estimated by the resonance frequency and the ratio between diffusion coefficient; while the damping
constant is estimated by the resonance amplitude and the width of half power point in the half spectrum.
In ARM, AR coefficient and prediction error are calculated by the Yule-Walker method using thedisplacement response of 1-dof subjected to white noise excitation. The spring constant and damping
constant are estimated by eigen frequency and damping ratio using AR coefficients. Furthermore, thediffusion coefficient is estimated using the variance of prediction error signal.
3.2 Result and Discussion
3.2.1 Dependency on the Number of Samples
The dependencies of the number of samples are shown in Fig. 1: the spring constant is 𝑘 = 1.2;
the damping constant is 𝑐 = 0.1; variance of input white noise excitation is 𝜎 = 1; initial conditions
are 𝑥(0) = 0 and 𝑣(0) = 0; and the sampling period is Δ𝑡 = 0.1. The left graph of (a) represents
the estimation result of spring constant; the right graph of (b) represents the estimation result of theratio between diffusion coefficient and damping constant. In addition, the blue, red, green, and black
dotted lines show the result of the proposed method, ARM, HPM, and the true values, respectively.
Focus on the estimation result of spring constant. The estimation values converge to the true values
when the samples are at least 2000 in all estimation methods. On the other hand, the estimation values
are not in agreement with the true values when there are less than 100 samples in all methods. The
54
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
proposed method has the least error values which is less than 10% when the number of samples ranges
from 100 to 2000.
In the case of HPM, the error between estimation and true values is approximately 20%. HPM iscalculated using the FFT of the auto-correlation function. The auto-correlation function consists of thesum of the noise signal and the periodic signal of the forced vibration. Thus, the FFT of the periodic
signal is independent of the noise signal. As the results showed, the decline of accuracy does not occur
when the samples are low. Here, the accuracy of HPM is less than the proposed method.
In the case of ARM, high estimation errors of more than 90% have been observed. Generally, ARM
has higher accuracy of estimation performance than the FFT method in low sample conditions.
However, the estimation accuracy decreased when there were less than 500 samples. In addition, theaccuracy of ARM is less than the HPM and the proposed method. Therefore, the result shows theapplicable limits of ARM.
Furthermore, focus on the estimation result of the ratio between the diffusion coefficient and thedamping constant. The proposed method converges to the true value with increasing number ofsamples. The estimation errors between the estimation and true values are at the smallest in our
proposed method. The proposed method has an estimation error of 30% in the range of over 1000
samples. Between 100 to 1000 samples, the estimation error is approximately 40%. Moreover, withsamples less than 100, the estimation error is over 100%.
In the case of HPM, the fluctuation of estimation values increased with increasing number ofsamples; especially, in the case of samples less than 100, the estimation value is not obtained in orderto estimate the resonance frequency is 0. On the other hand, in the case of ARM, errors between
estimation and true values occurred in the 103 number of samples.
5.0
MLE&FPeq
AR2
HalfPower
True

3.0

D/c [Nm]k [N/m]

4.0

2.0

10

4

10

2

MLE&FPeq
AR2
HalfPower
True

100
10

-2

10-4

1.0
0.0
1
10

10

2

10

3

10

4

Number of sample(a) the case of spring constant

10-6
1
10

10

2

10

3

10

4

Number of sample(b) the case of ratio between diffusion
coefficient and damping coefficientFig. 1. Benchmark test results in terms of sample number dependency.
3.2.2 Dependency on the Damping Coefficients
The dependencies of the damping constant are shown in Fig. 2: the spring constant is 𝑘 = 1.2;
variance of input white noise excitation is 𝜎 = 1; initial conditions are 𝑥(0) = 0 and 𝑣(0) = 0;
sampling period is Δ𝑡 = 0.1.; and the damping constants are in the range of 𝑐 = 0.01 to 𝑐 = 10.
The left graph of (a) represents the estimation result of spring constant while the right one of (b)represents the estimation result of the ratio between diffusion coefficient and damping constant. In
55
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
addition, the blue, red, green, and black dotted lines show the result of proposed method, ARM, HPM,
and the true values, respectively.
Here, the focus is on the estimation result of spring constant (in case of (a)). Our proposed method
is in good agreement with the true values throughout the range of values of the damping constant. On
the other hand, HPM and ARM are not in agreement with the true values when the damping constantis between 0.5 and 1.8. Our proposed method is based on the variance ratio of probability density
function between displacement and velocity. The probability density functions are obtained both inlow signal and high noise levels. As a result, the variance ratio of probability density function between
displacement and velocity can be precisely calculated. The proposed method estimates using the raw
time series data and is independent of the auto-correlation function and frequency spectrum data.
On the other hand, HPM and ARM are based on frequency and auto-correlation function domains.
In these methods, the noise component is more dominant than the signal. As a result, the spring
constant evaluated is less than true values.
Furthermore, focus on the estimation result of the ratio between the diffusion coefficient and thedamping constant (in case of (b)). The proposed method is in good agreement with the true values
throughout the entire range of ratio between diffusion coefficient and damping constant. HPM is notin agreement with the true values when 𝑐 > 0.8. ARM is not in agreement with the true values
throughout the entire range of the D/c parameter.
2.5

MLE&FPeq
AR2
HalfPower
True

1.5

D/c [Nm]k [N/m]

2.0

1.0
0.5
0.0
0.0

5.0

10.0

c [Ns/m]

10

4

10

2

MLE&FPeq
AR2
HalfPower
True

100
10

-2

10

-4

10-6
10-2

10-1

100

101

c [Ns/m](a) the case of spring constant(b) the case of ratio between diffusion
coefficient and damping coefficientFig. 2. Benchmark test result in terms of damping coefficient dependency.
4. Extension to the Recursive Algorithm
This chapter discusses the expansion of the above proposed method to the recursive algorithm tosave the capacity of the memory. The recurrence formula is derived by the relationship of the root
mean square as follows:𝑘𝑓[𝑗] =

[𝑗] =

𝑓[𝑗]
𝑔[𝑗]

(9)

1
𝑗−1
(10)
𝑓[𝑗 − 1] + 𝑥 [𝑗]𝑗𝑗56
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
𝑔[𝑗] =𝑗−1
1
𝑔[𝑗 − 1] + 𝑥 [𝑗]
𝑗
𝑗

(11)

Here, the initial values of 𝑓[0] and 𝑔[0] are zeros. In addition, variable 𝑗 represents an integer
to distinguish it from the continuous time variable, 𝑡. The ratio between the diffusion coefficient andthe damping constant is as follows:
𝐷
𝑐

[𝑗] =

𝑗 − 1𝐷𝑗 𝑐

1[𝑗 − 1] + 𝑥 [𝑗]𝑗(12)

The initial value is obtained as a result of the above batch processing algorithm. The estimation
results using Eq. (9) ~ (12) are shown in Fig. 3. Here, the horizontal axes represent the time; the verticalaxes represent the spring constant (a) and the ratio between the diffusion coefficient and the damping
constant (b). The blue solid and black dotted lines show the estimation values of the unknown
parameter and true values, respectively. Here, the spring constant is 𝑘 = 1.2; the damping constant is
𝑐 = 0.5; variance of input white noise excitation is 𝜎 = 1; initial conditions are 𝑥(0) = 0 and𝑣(0) = 0; and the sampling period is Δ𝑡 = 0.1.
The estimation results of both the spring constant and the ratio between the diffusion coefficientand damping constant converge to the true values. Therefore, the operation of recursive algorithm isverified. With the proposed online estimation algorithm, t should be greater than 500 to converge withtrue values, with the number of samples being 5000. On the other hand, 1000 samples are necessary
for the above batch processing. Therefore, the number of samples for the proposed recursive algorithm
is approximately the same as number of samples of the above batch algorithm.

(a) the case of spring constant(b) the case of ratio between diffusion
coefficient and damping coefficientFig. 3. Estimation result using the proposed online estimation algorithm.
5. Damping Coefficients Dependency on Recursive Algorithm
In this chapter, the damping coefficient dependency of the recursive algorithm is considered. Here,the spring constant is 𝑘 = 1.2; variance of input white noise excitation is 𝜎 = 1; initial conditions
are 𝑥(0) = 0 and 𝑣(0) = 0; sampling period is Δ𝑡 = 0.1; and the damping constants are between𝑐 = 0.1 to 𝑐 = 10. The recursive algorithm is performed till 2000 s; then, the ratio between the trueand estimation value is calculated by the converged estimation value using the proposed algorithm.
The estimation results are shown in Fig. 4 and 5. In Fig. 4, the horizontal axis represents the damping
coefficient while the vertical axis represents the ratio between true value and estimation value. The
57
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)

Ratio between true and estimate of kcenters of the ratio is near the value of 1 and fluctuates in the 0.8 ~ 1.2 range. The operation of the
proposed algorithm to estimate the spring constant is confirmed in the large damping system.
Furthermore, focus on the estimation result of ratio between the diffusion coefficient and the damping
constant in Fig. 5. The ratio between the diffusion coefficient and the damping constant decreases as
the damping constant increases and slightly fluctuates between 0.9 and 1.05. The operation of the
proposed algorithm to estimate the ratio between the diffusion coefficient and the damping constant isconfirmed in the large damping system.

Fig. 4. Damping coefficient dependency on estimation accuracy in case of spring constant.
2.0

1.5

1.0

0.5

0.0
0.0

2.0

4.0

6.0

8.0

10.0

c [Ns/m]Fig. 5. Damping coefficient dependency on estimation accuracy in case of ratio between
diffusion coefficient and damping coefficient.

58
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)
6. Conclusions
In this paper, we discussed the identification method based on the MLE using analytical solution of
Fokker-Planck equation. As the result, the following results are obtained.
(1) Estimation formulas of the spring constant and the ratio between the diffusion coefficient andthe damping constant were derived using MLE method using Analytical solution of Fokker-Planck
equation.
(2) The benchmark tests in terms of estimation accuracy were conducted in cases of HPM and ARM,
respectively. As the result, our proposed method showed the most accuracy in cases of the number ofsamples dependence and damping constant dependence.
(3) The expansion to the recursive algorithm were conducted beyond the batch processing algorithm.
In addition, the estimation experiments were conducted using the numerical simulation.
(4) The verification of proposed method in terms of estimation accuracy were conducted. As theresult, our proposed method showed the applicability to large damping system.
References[1] K. Suzuki and A. Nakajima, “Practical method to estimate damping ratio based on thecharacteristics of AR spectral analysis technique”, Transactions of the Japan Society of
Mechanical Engineers, Series C, Vol. 48, No. 433, pp. 1389 – 1397, 1982 (in Japanese).
[2] K. Suzuki and Y. Harada, “Analysis of structural nonstationary behavior and dynamic system
identification based on the auto-regressive model fitting”, Transactions of the Japan Society of
Mechanical Engineers, Series C, Vol. 50, No. 458, pp. 2013 – 2020, 1984 (in Japanese).
[3] H. Sohn and C. R. Farrar, “Damage diagnosis using time series analysis of vibration signals”,
Smart Materials and Structures, Vol. 10, No. 3, pp. 446 – 451, 2001.
[4] N. K. Krishnan and A. S. Kiremidjian, “Time series based structural damage detection algorithm
using Gausian mixtures modeling”, Transactions of the ASME, Journal of Dynamic system,measurement and control, Vol.129, No. 3, pp. 285 – 293, 2007.
[5] S. Adachi, Fundamentals of system identification, Tokyo Denki University Press, Tokyo/Japan,
2009, pp. 98 – 99 (in Japanese).
[6] Y. Lu and F. Gao, “A novel time-domain auto-regressive model for structural damage diagnosis”,
Journal of Sound and Vibration, Vol. 283, No. 3-5, pp. 1031 – 1049, 2005.
[7] Y. C. Zhu and S. K. Au, “Bayesian operational modal analysis with asynchronous data, part I:
Most probable value”, Mechanical Systems and Signal Processing, Vol. 98, pp. 652-666, 2018.
[8] Y. C. Zhu and S. K. Au, “Bayesian operational modal analysis with asynchronous data, part II:
Posterior uncertainty”, Mechanical Systems and Signal Processing, Vol. 98, pp. 920-935, 2018.
[9] M. H. Masjedian and M. Kershmiri, “A review on operational modal analysis researches:classification of methods and applications”, Proc. IOMAC 2009 the 3rd International Operational
Modal Analysis Conference (Portonovo, Italy) May 2009.
[10] R. Brincker, P. Anderson and L. Zhang, ‘Modal identification from ambient responses using
frequency domain decomposition’, Proc. IMAC 2000 the 18th International Modal Analysis
Conference (San Antonio, Texas) February 2000.
[11] Y. C. Zhu, S. K. Au and J. M. W. Brownjohn, “Bayesian operational modal analysis with buriedmodes”, Mechanical Systems and Signal Processing, Vol.121, pp. 246 – 263, 2019.
59
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020

Journal of Mechanical and Electrical Intelligent System (JMEIS)[12] T. K. Caughey, “Derivation and application of the Fokker-Planck equation to discrete nonlinear
dynamic systems subjected to white random excitation”, The Journal of the Acoustical Societyof America, Vol. 35, No. 11, pp. 1683 – 1692, 1963.
[13] P. Kumar, S. Narayanan, S. Adhikari and M. I. Friswell, “Fokker-Planck equation analysis ofrandomly excited nonlinear energy harvester”, Journal of Sound and Vibration, Vol. 333, No. 7,
pp. 2040 – 2053, 2014.
[14] P. Kumar and S. Narayanan, “Solution of Fokker-Planck equation by finite element and finite
difference methods for nonlinear systems”, Sadhana, Vol. 31, pp. 445 – 461, 2006.

60
J. Mech. Elect. Intel. Syst., Vol.3, No.2, 2020