# Estimation of the Dynamic Moduli of Viscoelastic Asphalt Mixtures Using the Extended Kalman Filter Algorithm.

1. Introduction

The dynamic modulus is the major linear viscoelastic (LVE) parameter of asphalt and has been widely used to study the asphalt/concrete mixtures used to make pavements and rehabilitate roads. Thus, laboratory determination of the dynamic moduli of LVE asphalts is very important. Uniaxial compression is used to this end but requires heavy and expensive testing machines, such as material testing system (MTS) loading frames, equipped with either a 25 or 8.9 kN loading cell depending on the purpose of the test, and a 30 kN Universal Testing Machine (UTM; IPC Global). Normally, dynamic moduli are derived using an MTS or UTM because the values obtained have better global credibility than impact resonance test (IRT) data. However, certain problems are apparent. First, equipment maintenance and operation is expensive. Second, a sinusoidal loading must be directly applied. Thus, it is difficult to both define and adjust the linear range of the applied load during testing; it is essential that the material properties of the asphalt remain unaltered.

Use of the IRT ensures that specimens are not damaged during testing; impact force is delivered by a small steel ball. However, the IRT is regarded less favorably than the MTS and UTM tests because few discrete data points are obtained. Thus, we developed a predictive model of the dynamic modulus using the EKF algorithm and IRT test data. When the EKF algorithm was employed to determine the coefficients and the dynamic moduli of hot-mix asphalt/concrete using IRT data, a reasonably appropriate master curve was created compared to that afforded by UTM measurements. We created predictive models of the dynamic modulus by frequency (e.g., at 0.1, 0.5, 1, 5, 10, and 20 Hz) at the temperatures of 5.3, 20.1, 40.1, and 53.7[degrees]C. We used the IRT test to obtain dynamic moduli at 20, 30, 35, 40, and 50[degrees]C.

The remainder of this paper is organized as follows: Section 2 describes the EKF algorithm and the Taylor series, which are related to the application of an optimization algorithm in this study. Section 3 explains the IRT and conventional, dynamic UTM modulus derivation, which have been used for characterizing asphalt concrete mixtures. Section 4 discusses the test results and comparisons, based on developing a nonlinear master curve and a predictive model for a nonlinear, sigmoidal master curve reflecting the dynamic modulus. Section 5 implements the EKF algorithm according to the master curve; furthermore, the efficiency of the EKF algorithm is discussed in terms of the comparison between the applied EKF and conventional regression algorithms. Finally, we conclude this research, based on discussing the accuracy of the applied EKF algorithm.

2. The Extended Kalman Filter Algorithm

The Kalman filter algorithm is an optimization algorithm analogizing posterior values by adding various weights to stochastic moving averages using a recursive data processing method. The extended Kalman filter (EKF) algorithm, related to the above algorithm, is derived from the linearized Kalman filter algorithm for a linear system and can be applied to nonlinear systems such as those of the present study. The EKF algorithm features two major steps: estimation and prediction [1-3].

At first, initial values, such as white and observational noise and the true state, are used to execute the EKF algorithm. Using these values, a priori estimates and error covariance can be calculated during prediction. The dynamic modulus values are predicted by a Taylor series (the higher order differentials are ignored) that expresses the modulus with respect to the frequency of the linear state. This can be written as follows:

[[??].sub.[bar.k]] = F([[??].sub.k-1], [fre.sub.k-1]), [[bar.E].sub.k] = h([[??].sub.[bar.k]],[V.sub.k]), (1)

where [[??].sub.[bar.k]] and [[bar.E].sub.k] are the estimate and the measurement respectively; fre is the frequency; and [W.sub.k] and [V.sub.k] are the white and observation noise, respectively (and which are normally distributed). We thus obtain

P([W.sub.k]) ~ N(0,Q), P([V.sub.k]) ~ N(0,R). (2)

In practice, Q and R represent the process and measurement noise covariances, respectively. The state Jacobian used to linearize a nonlinear model (such as the linearization Kalman filter algorithm) can be expressed as follows:

[mathematical expression not reproducible]. (3)

This Jacobian, and the a priori error and white noise covariance, is used to derive a posterior error covariance:

[P.sub.[bar.k]] = [AP.sub.k-1][A.sup.T] + Q([W.sub.k]). (4)

Next, the a priori estimate and error covariance derived during prediction are used to in turn derive an estimate with error covariance. As indicated in Figure 1, the Kalman gain is initially derived. The role played by this Kalman gain becomes clear on comparing the next step in the calculation of the dynamic modulus with the first-order low-pass filter outcome for a high-frequency signal. This outcome and calculation of the estimated dynamic modulus are shown below:

[[bar.E].sub.k] = (1 - K)[[bar.E].sub.k-1] + K[E.sub.k], (5a)

[[??].sub.k] = [[??].sub.-k] + [K.sub.k]([[bar.E].sub.k] - h([[??].sub.[bar.k]])). (5b)

The estimated value ([[??].sub.k]) can be derived from the function (h([[??].sub.[bar.k]])), which reflects the relationship between the measured value ([[bar.E].sub.k]) and both Kalman gain and state variable. In equation (5a), K is the weight used to change the value. Thus, comparison of the two equations above shows that the Kalman gain is the weight of the change. Therefore, the Kalman gain can be derived from the estimation error covariance, the linear model, and the noise matrix (e.g., R of Equation (6)) associated with the measured variable. A linear model is obtained from the linearized function reflecting the relationship between the measured value and the state variable:

[K.sub.k] = [P.sub.[bar.k]][H.sup.T]/H[P.sup.[bar.k]][H.sup.T] + R. (6)

When the EKF algorithm is iterated, the error covariance is derived at the last step; this covariance is used to measure the stochastic accuracy of the dynamic modulus estimate at the frequency in question.

The posterior estimate, calculated during estimation, can be used to determine the a priori estimate for the next step (in which the posterior estimate is applied to the previous a priori estimate). This repeated prediction-to-estimation process continues to the final frequency (5,000 Hz). The procedure is shown in Figure 1. Based on the recurrent relationships among moving averages obtained after first-order low-pass filtering, the EKF algorithm changes the weight of the Kalman gain to predict the result of repeated estimation and prediction processes, based on covariance and the standard normal deviation.

3. Impact Resonance and the Typical Compressive Test

Here, we used the impact resonance test (IRT) to determine the dynamic moduli of hot-mix asphalt (HMA) mixtures. When the material properties were defined, a discrete dynamic modulus in the low-frequency domain was determined via the IRT; this was compared to the dynamic modulus obtained by UTM testing. Dynamic moduli were obtained at different frequencies; the IRT was used to gather data at various temperatures.

3.1. Impact Resonance Test. The IRT is conducted in a hygrostat chamber of uniform temperature and humidity (Figure 2). The specimen, suspended from two strings, is impacted by a steel ball. A Briiel and Kj\$r pulse device receives a vibrational analog signal upon impact and digitizes the signal. Then, a frequency response function (FRF) is calculated. The longitudinal vibration of a uniform, LVE specimen can be represented by a one-dimensional wave equation:

[[partial derivative].sup.2](x,t)/[partial derivative][x.sup.2] = [rho]/E* [[partial derivative].sup.2]u(x,t)/[partial derivative][t.sup.2], (7)

where u (x, t) is the longitudinal displacement at position x along the string at time t and [rho]/E* is the reciprocal of the complex wave velocity c* as follows:

c* = [square root of E*/[rho]], (8a)

E* = E' + IE" = E'(1 + i[eta]), (8b)

where E' and E", respectively, are the real and imaginary parts of the complex modulus E* (e.g., the dynamic modulus of the present study); furthermore, i is [(-1).sup.1/2] and [eta] is the loss factor. To solve the above two equations, separation of variables, imposition of initial and boundary conditions, and use of a Fourier series are required. The equation is divided into a local and a time function:

u(x,t) = [[A.sub.n] cos(n[pi]/l [square root of E*/[rho] t)] + [B.sup.n]sin(n[pi]/l[square root of E*/[rho]t])sin(n[pi]/l x), (9)

where l is the effective length of the string wave; n is the wave number; [rho] is the vibrational density after impact; and [A.sub.n] and [B.sub.n] are the real components of mutual complex conjugates. In the above equation, the sum of the time functions is the amplitude of the wave equation.

Many researchers [4-8] performed the IRT and found that when the impact wave penetrated the full width of the specimen, the input mass ([m.sub.1]) of the impacted surface differed from the output mass ([m.sub.2]) of the opposite surface because of a cancellation effect. Use of a general transmissibility function showed that the input mass was the sum of a mass [m.sub.trans] attributable to offsetting and the output mass. If we exclude the effect of accelerometer weight on the peak transmissibility associated with the cancellation effect, the transmissibility functional amplitude is

[absolute value of [[tau].sub.amp]] = 1/[[square root of [cos.sup.2](p)[cosh.sup.2]([epsilon]p)+[sin.sup.2](p)[sinh.sup.2]([epsilon]p)], (10)

where p is the product of the complex wave number (n*) and the length of the complex as follows:

n* = [square root of [pi][[omega].sup.2]/E*] = [omega][square root of [rho]/E'(1+i[eta])] = p(1 + i[epsilon]), (11)

where [epsilon] is related to the loss factor and [omega] is the angular frequency. p can be derived by separating the equation into imaginary and real parts:

[epsilon] = [square root of (D-1)/(D+1)], (12a)

D = [square root of (1+[[eta].sup.2]) = [absolute value of E*]/E' = E/E', (12b)

sin(2p) = 0,

p = ([omega][iota])/D)[square root of (D+1)[rho]/2E']. (12c)

When a perfect (no-loss) elastic string is used, the transmissibility will be maximal and given by

p = (2n-1)[pi]/2, n = 1, 2, 3, ..., m. (13)

The relationship between the dynamic modulus ([E.sub.n]) and the modal frequency ([f.sub.n]) can be derived using p of Equations (8a), (8b) and (12a)-(12c) as follows:

[f.sub.n] = 2n-1/4 [square root of 2D[E.sub.n]/(1+D)[rho]], n = 1, 2, 3 ..., [infinity]. (14)

In IRTs performed when both ends of the specimen are longitudinally fixed, the frequency resonance function (FRF) around the first maximum excitation (thus over a very short time) is derived.

3.2. Typical Dynamic Modulus Test. We used the UTM to determine the dynamic modulus following the AASHTO-62 procedure . The procedure features six loading rates (20, 10, 5, 1, 0.5, and 0.1 Hz) at 5.3, 20.1, 40.1, and 53.7[degrees]C; the specimen was that previously used in the IRTs. In the conventional compressive derivation of the dynamic modulus, linear variable differential transformers (LVDTs) are attached to the sides of the specimen. A sigmoidal function [10-15] can be used to construct a master curve predicting the dynamic modulus of a hot-mix asphalt (HMA) specimen as follows:

log(E*) = [delta] + [alpha]/[1+exp([beta]+[gamma] log([f.sub.R])), (15)

where E* is the dynamic modulus; [delta] is the minimum value of that modulus; [alpha], [beta], and [gamma] are the model parameters; and [f.sub.R] is the reduced frequency (derived by time/temperature superimposition). The reference temperature for shift factor establishment is 20[degrees]C. The reduced frequency ([f.sub.R]) can be obtained using the shift factor as follows:

[mathematical expression not reproducible]. (16)

The shift factor of an HMA specimen can be modeled as follows:

log ([a.sub.T]) = a x [f.sup.2] + b x f + c. (17)

4. Test Results

We drew master curves from low to high frequencies based on linear viscoelastic (LVE) parameters; we used both UTM (universal testing machine) and the impact resonance test (IRT) to this end. However, the IRT yielded only a few discrete points. To compare UTM test and IRT data, master IRT curves were obtained by analyzing IRT discrete point data using the extended Kalman filter (EKF) algorithm; the curves were then compared. The hot-mix asphalt (HMA) specimen contained 10% air void. The asphalt binder of PG (performance grade) 76-22 was used. The gradation used for the mixture is shown in Table 1.

The IRT specimen weighed with 2.57 kg and was 150 mm in height and 100 mm in diameter. We performed dynamic modulus testing under AASHTO-62 conditions . The test machine is shown in Figure 3. The IRT was used to obtain dynamic moduli at the temperatures of 20, 30, 35, 40, and 50[degrees]C; UTM testing was performed at 5.3,20.1,40.1, and 53.7[degrees]C using the IRT specimen. The IRT setup is shown in Figure 4.

4.1. Master Curve Determination Using the UTM Test. A shift factor based on the time-temperature superimposition principle (TTSP) can be used to construct a master curve of the dynamic modulus at any temperature using data from the UTM test. The reference temperature used to construct the regression model for shift factor determination was set to 20[degrees]C, estimated using the least squares error method of the second-order polynomial equation (Equation (17)), as shown in Figure 5.

The table below shows the results of the nonlinear regression analysis (based on the least squares error method) used to determine the parameters of the predictive model, where RMSE is the root-mean-square error; R square is the simple Pearson correlation; MPE is the mean percentage error; and the p value is the stochastic significance level. Using a standard sigmoidal function, the comparison between estimated and measured values can be represented as shown in Table 2 and Figure 6

4.2. Master Curve Determination Using the IRT. In the IRT, both sides of the specimen are fixed when the steel ball impacts; the temperature and humidity are held constant. The ball generates an impulse wave measured by both input and output accelerometers. The frequency response function (FRF) can be calculated using a fast Fourier transform (FFT).

Viscoelastic material can be mathematically modeled in terms of the spring constant and the damping coefficient. However, it is not guaranteed that resonance will occur at the highest FRF. Therefore, a coherence function should be used to evaluate response function reliability:

[C.sub.xy]([omega]) = [P.sub.xy]([omega])/[square root of [P.sub.xx]([omega])[P.sub.yy]([omega])], (18)

where [P.sub.xx]([omega]) and [P.sub.yy] ([omega]) are the power spectrum densities of the input and output signals, respectively, and [P.sub.xy]([omega]) is the cross-spectral density. The dynamic modulus can be calculated using the IRT, based on Equation (14):

[E.sub.n] = 16[rho] x [f.sup.2.sub.n] x [l.sup.2]/[(2n-1).sup.2](1+D/2D). (19)

The dynamic modulus is expressed in terms of the resonance frequencies shown in Table 3. The shift factor for IRT was earlier  shown to be

log ([a.sub.T]) - 0.0006[f.sup.2] - 0.157f + 3.5399. (20)

Sigmoidal function parameters can be estimated using the generalized reduced gradient (GRG) method to in turn determine the master curve. The reduced dynamic modulus and its curve is shown in Table 3 and Figure 7.

4.3. Comparisons. The two master curves derived are shown in Figure 8 and were compared to improve confidence in the IRT data. Next, a combined master curve was created to exploit all available data. The root-mean-square error (RMSE) and R square differences between measured and predicted moduli were lowest and highest, respectively, for the dynamic modulus test, reflecting the lower error and a more appropriate trend. However, the IRT MPE (mean percentage error) value was lower. Also, the master curves of the IRT and dynamic modulus test exhibited higher MPEs than the curves of the other tests. However, the extent of agreement was about 99%.

5. Model Using the Extended Kalman Filter Algorithm to Predict the Dynamic Modulus

The accuracy of predictions afforded by the extended Kalman filter (EKF) algorithm was checked by comparing the regression algorithms of the sigmoidal functions, using data from the UTM test and the IRT. We initially fixed data obtained from the sigmoidal function. The standard normal deviation of the dynamic modulus data was first calculated, and the error covariance and average preexisting values of that deviation were determined as

E([E.sub.N]) = [[summation].sup.N.sub.n=1][E.sub.n]/N, [P.sub.k] = E{[([E.sub.k] - [[??].sub.k])([E.sub.k] - [[??].sub.k]).sup.T]}, (21)

where E([E.sub.N]) is the average dynamic modulus and [P.sub.k] is the error covariance. Figure 9 shows the standard normal deviations of the dynamic moduli derived by the UTM test and the IRT.

An a priori estimate was obtained by applying the Taylor series, which assumes that the relationship between the dynamic modulus to be calculated and the previous dynamic modulus is linear:

E([fre.sub.k]) = E([fre.sub.k-1]) + E'([fre.sub.diff]) x [fre.sub.diff], (22)

where k is the number of matrix elements and diff is the difference in the proportions of linear relationships among frequencies. During the iterative calculations, the estimated distance between frequencies was 0.1 Hz, and this distance was multiplied by 10 to allow integral permutations to the number of matrix elements. The next procedure is shown in Figure 1, and a result was obtained after 57 loops; Figure 10 shows the results.

Optimization based on the EKF algorithm increased the prediction accuracy. As shown in Table 4, the RMSE and MPE values decreased, and the R square values increased, compared to those of the regression approach. Figures 8 and 10 show both results of nonlinear regression method and extended Kalman filter (EKF) optimization algorithm. In comparison, the nonlinear regression method (Figure 8) shows relatively large discrepancy between continuous master curves of UTM-30 and IRT, when compared to the continuous master curves determined from the EKF optimization algorithm (Figure 10). Thus, small improvement of the RMSE, MPE, and R square can result in the nice master curves of two different testing methods because the continuous master curves of dynamic moduli are required to represent the wide range of frequency from 0.0000001 to 100000 Hz.

6. Conclusions

Asphalt concrete materials widely used in pavement construction industries can be complicatedly characterized with respect to viscoelasticity, when compared to other elastic materials such as steel and cement concrete, which are well known in civil and construction engineering. Thus, the frequency-dependent or loading rate-dependent dynamic modulus should be determined in order to represent the linear viscoelastic (LVE) material parameter of asphalt concrete mixture. In this study, two different testing systems were used to determine the LVE material parameters, based on the conventional loading frame system of 30 kN Universal Testing Machine (UTM) and impact resonance test (IRT).

Thus, the LVE material parameters were determined through the extended Kalman filter (EKF) optimization algorithm and typical regression analysis. The EKF algorithm implements two procedures of prediction and estimation; finally, the error between prediction and measurement reduces during the iterations of two procedures. On the contrary, the typical nonlinear regression analysis is based on the least squares error method.

Finally, we used the EKF algorithm to predict the LVE properties of hot asphalt concrete mixtures, resulting in obtaining the dynamic modulus parameters through UTM and IRT. Thus, we compared dynamic moduli derived by regression and application of the EKF algorithm to nonlinear sigmoidal functions, using data obtained from the UTM test and the IRT. The EKF algorithm, with omission of the higher order Taylor series, reduced more errors, affording higher R squared and lower RMSE and MPE values, when compared to the typical nonlinear regression analysis.

https://doi.org/10.1155/2018/3089085

Data Availability

The testing data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the Research Program funded by the SeoulTech (Seoul National University of Science and Technology).

References

 F. Gao and Y. Lu, "A Kalman-filter based time-domain analysis for structural damage diagnosis with noisy signals," Journal of Sound and Vibration, vol. 297, no. 3-5, pp. 916-930, 2006.

 R. Wu, J. W. Choi, and J. T. Harvey, Extended Kalman Filter and its Application in Pavement Engineering, Intelligent and Soft Computing in Infrastructure Systems Engineering, Vol. 259, 2009.

 J. W. Choi, R. Wu, and J. M. Pestana, "New layer-moduli back-calculation method based on the constrained extended Kalman filter," Journal of Transportation Engineering, vol. 136, no. 1, pp. 20-30, 2010.

 S. O. Oyadiji and G. R. Tomlinson, "Determination of the complex moduli of viscoelastic structural elements by resonance and non-resonance methods," Journal of Sound and Vibration, vol. 101, no. 3, pp. 277-298, 1985.

 M. M. Rodriguez, Analysis of Structural Damping, Thesis, Lulea University of Technology, Lulea, Sweden, 2006.

 N. Ryden, "Resonant frequency testing of cylindrical asphalt samples," European Journal of Environmental and Civil Engineering, vol. 15, no. 4, pp. 587-600, 2011.

 A. Gudmarsson, N. Ryden, and B. Birgisson, "Characterizing the low strain complex modulus of asphalt concrete specimens through optimization of frequency response functions," Journal of the Acoustical Society of America, vol. 132, no. 4, pp. 2304-2312, 2012.

 S. Mun, "Determining the dynamic modulus of a viscoelastic asphalt mixture using an Impact Resonance Test with damping effect," Journal of Research in Nondestructive Evaluation, vol. 26, no. 4, pp. 189-207, 2015.

 American Association of State Highway and Transportation Officials (AASHTO), Standard Method of Test for Determining Dynamic Modulus of Hot Mix Asphalt (HMA), Vol. 62, AASHTO TP, Washington, DC, USA, 2007.

 M. W. Witczak, K. Kaloush, T. Pellinen, M. El-Basyouny, and H. V. Quintus, "Simple Performance test for superpave mix design," NCHRP Report No. 465, Transportation Research Board, National Research Council, National Academy Press, Washington, DC, USA, 2002.

 J. S. Daniel and Y. R. Kim, "Development of a simplified fatigue test and analysis procedure using a viscoelastic continuum damage model," Journal of Association of Asphalt Paving Technologists, vol. 71, pp. 619-650, 2002.

 G. R. Ghehab, Y. R. Kim, R. A. Schapery, M. W. Witczak, and R. Bonaquist, "Time-temperature superposition principle for asphalt concrete mixtures with growing damage in tension," Journal of Association of Asphalt Paving Technologists, vol. 71, pp. 559-593, 2002.

 S. Mun, G. R. Chehab, and Y. R. Kim, "Determination of time-domain viscoelastic functions using optimized interconversion techniques," Road Materials and Pavement Design, vol. 8, no. 2, pp. 351-365, 2007.

 S. Mun and Z. W. Geem, "Determination of viscoelastic and damage properties of hot mix asphalt concrete using a harmony search algorithm," Mechanics of Materials, vol. 41, no. 3, pp. 339-353, 2009.

 S. Lee, S. Mun, and Y. R. Kim, "Fatigue and rutting performance of lime-modified hot-mix asphalt mixtures," Construction and Building Materials, vol. 25, no. 11, pp. 4202-4209, 2011.

Yu-Seok Gong, (1) Dowan Kim, (2) and Sungho Mun [ID], (3)

(1) Highway and Transportation Technology Institute, Korea Expressway Corporation, 208-96, Dongbu-daero 922 Beon-Gil, Dongtan-Myeon, Hwaseong-Si 18489, Gyeonggi-Do, Republic of Korea

(2) Department of Road & Airport, Kunhwa Engineering & Consulting Co., Ltd., 321 Teheran-Ro, Gangnam-Gu, Seoul, Republic of Korea

(3) Department of Civil Engineering, Seoul National University of Science & Technology, 232 Gongneung-Ro, Nowon-Gu, Seoul 01811, Republic of Korea

Correspondence should be addressed to Sungho Mun; smun@seoultech.ac.kr

Received 4 July 2018; Revised 25 October 2018; Accepted 4 November 2018; Published 10 December 2018

Caption: Figure 1: Flow chart of the EKF algorithm.

Caption: Figure 2: Equipment setup of the IRT.

Caption: Figure 3: The testing machine of UTM.

Caption: Figure 4: Setup of IRT.

Caption: Figure 5: Shift factor of the tested dynamic modulus.

Caption: Figure 6: Comparison between measured data of UTM and predicted sigmoidal function.

Caption: Figure 7: Comparison between measured data and predicted sigmoidal function from IRT.

Caption: Figure 8: Comparison between UTM and IRT.

Caption: Figure 9: Standard normal distribution of the dynamic modulus of UTM and IRT.

Caption: Figure 10: Master curves using the EKF algorithm.
```Table 1: Gradation used for the mixture.

Sieve size        % passing

1/2" (12.5 mm)       100
3/8" (9.5 mm)       51.1
#4 (4.75 mm)        20.0
#8 (2.36 mm)        15.6
#30 (0.600 mm)      10.0
#50 (0.300 mm)       7.7
#100 (0.150 mm)      5.4
#200 (0.075 mm)      4.3
Pan                  0.0

Table 2: Parameters of the sigmoidal function, e.g., Equation (15))
and statistical results.

Parameter   Value   Statistical result     Value

[delta]     2.36        RMSE (MPa)       491.18 MPa
[alpha]     2.01         R square          0.997
[beta]      1.62           MPE             0.081
[gamma]     0.47         p value          7.83E-26

Table 3: Dynamic modulus and resonance frequencies of IRT.

Temperature ([degrees]C)     20        30        35        40

Resonance                  14.0035   1.65295   0.67565   0.39734
frequency (Hz)
Dynamic modulus            16965.2   14508.4   11511.5   13540.4
(MPa)

Temperature ([degrees]C)     50

Resonance                  0.17150
frequency (Hz)
Dynamic modulus            10720.5
(MPa)

Table 4: Improvement of EKF algorithm.

UTM                         IRT

Regression   EKF algorithm   Improvement   Regression

[delta]       2.36          2.36            --           2.36
[alpha]       2.01          2.01                         2.01
[beta]        1.62          1.71                         1.98
[gamma]       0.47          0.44                         0.48
RMSE         491.18        481.62          9.56         902.29
MPE          0.081          0.067          0.014        0.0524
R square     0.9969        0.9971         0.0002        0.9138

IRT

EKF algorithm   Improvement

[delta]        2.52            --
[alpha]        1.86
[beta]         1.86
[gamma]        0.49
RMSE          901.65          0.64
MPE           0.0470         0.0054
R square      0.9140         0.0002
```