Improved Hybrid Model for Predicting Concrete Crack Openings Based on Chaos Theory

Conventional statistical models provide inaccurate predictions of concrete crack openings because they do not consider the nonlinear temperature response and the residual characteristics of concrete. To address this problem, this study introduces a nonlinear temperature factor and develops an improved statistical model of crack openings.*e chaotic characteristics of residual time series of the improved statistical model are analyzed based on chaos theory and phase-space reconstruction theory. *ese theories are integrated with back-propagation (BP) artificial neural networks and genetic algorithms (GAs) to establish a GA-BP neural network model for predicting residuals. Finally, a hybrid model is developed for predicting the concrete crack opening behavior. *e predictions of the conventional statistical model, the statistical model considering nonlinear temperature component, and the hybrid model are compared using the case study on the crack openings of a regulating sluice.*e results show that the proposed hybrid model in this study for predicting concrete crack openings is significantly more accurate than the conventional statistical model and the statistical model considering nonlinear temperature component.


Introduction
Cracks are a common structural defect in hydraulic concrete. e operational behavior of macroscopic cracks in hydraulic concrete structures is directly related to the safety of the structure and is an important indicator of structural stability. To determine the operational behavior of concrete crack openings, the momentum misplacement is often monitored in real time by installing a crack meter group. Given the conditions of practical projects, concrete cracks are affected by various complex factors, including water pressure, temperature, and time. Hence, it is vital to establish an accurate mathematical model to aid in monitoring concrete crack openings and to quantitatively analyze, evaluate, and predict their behavior.
Concrete cracks fall into the category of concrete deformation. At present, mathematical models based on conventional statistical models of concrete deformations are widely used in engineering monitoring and parameter inversion [1]. Such models usually include temperature component, hydraulic component, and aging component. Experiments demonstrate that temperature is crucial to control concrete deformation [2,3] and that variations in temperature are an important factor in determining the opening and closing of concrete cracks [4]. e temperature component is thus considerably sensitive to the temperature variations [1]. Although the hydraulic and aging components of conventional statistical deformation models generally fit well to experimental results, the temperature component is underfitted because conventional statistical models [5,6] generally consider only the linear response of the system to temperature but neglect the complex nonlinear relationship between the actual affected variables and environmental factors in practical engineering. Previous studies show that variations in temperature are related nonlinearly to concrete deformation, and including such nonlinearities in prediction models should provide more accurate results than conventional models [7,8].
To further improve the accuracy of statistical models, chaos theory has recently been harnessed to analyze and predict the residual series of statistical models, effectively establishing a hybrid model that provides more accurate predictions than the statistical models. For instance, Wei et al. [9] established a back-propagation autoregressive integrated moving-average model based on the wavelet analysis that directly predicts chaotic residual time series. e nonlinear method predicts the residual, which is generated by a phase-space reconstruction, with the result being a hybrid model based on a conventional statistical model [10,11].
e results show that hybrid models based on residual prediction are more accurate than conventional models. Bao and Wu [10] improved the mathematical model for monitoring concrete cracks by combining chaos theory with a conventional statistical model to establish a hybrid model of concrete crack openings. e results indicate that the stress field at the crack tip generates singularities when exposed to water pressure and temperature variations, which irreversibly deforms the concrete at the crack tip. Hence, crack openings also involve an irrecoverable part caused by nonlinear factors, suggesting that conventional statistical models cannot predict the deformation behavior as a function of water pressure and temperature. is part of the deformation contains significant information about the crack evolution process, which is usually described in terms of a residual series, that is, a chaotic residual time series. Hence, it is necessary to analyze the residual time series of statistical models of concrete crack openings.
At present, there are some reports on the prediction methods of time series [12]. e studies show that methods based on phase-space reconstruction are commonly used to predict chaotic time series. e essence of this approach is to transform the extension of time series into the interpolation of phase space. Based on the Takens embedding theorem [13], the optimized delay time and minimum embedding dimension are determined according to the actual time series. e one-dimensional (1D) time series can then be reconstructed into a phase space with the same topological meaning as the motive power system. A reasonable selection of the delay time and embedding dimension then allows reconstruction of the phase space with satisfactory adaptability, thereby producing the internal nonlinear mapping. By exploiting the properties of the Lyapunov exponent, this approach can recognize and make short-term predictions of chaotic time series [14,15]. Currently, the methods for determining optimized delay time and minimum embedding dimension depend on the amount of data they are fed, the amount of computation, the anti-noise performance, and the objectivity with which the parameters are selected. e methods available for selecting optimized delay time typically include the autocorrelation coefficient method [16], the mutual information method [17], and the average displacement (AD) method [18], and the methods used to calculate the minimum embedding dimension typically include the geometric invariant method [14], the pseudoneighbor method [19], and the Cao method [20].
Casdagli [21] proposed the basic principle for predicting chaotic time series, which was used to develop prediction methods involving artificial neural networks [9,22,23]. BP neural networks are suitable for nonlinear mapping, making them good candidates for predicting chaotic time series. However, single neural networks converge slowly, resulting in local minima and overfitting. Ding et al. [24] proposed GA-BP neural networks and demonstrated that genetic algorithms (GAs) provide obvious improvements over BP neural networks. Hence, phase-space reconstruction combined with GA-BP neural networks is expected to accurately predict chaotic time series.
In this study, the conventional statistical model is firstly improved by considering the nonlinear temperature factor. Further, the improved statistical model is combined with the GA-BP neural network model using chaos theory and phasespace reconstruction theory, to develop a hybrid model for predicting concrete crack openings. Finally, a detailed case study on the crack openings of a regulating sluice is presented to illustrate the application of the hybrid model.

Improved Statistical Model for Concrete Crack Openings
Hydraulic concrete structures with macroscopic cracks are affected by water pressure, temperature, and time. It is illustrated in Figure 1. Consider a hydraulic concrete structure with K measuring points near a macroscopic crack. A conventional statistical model describing concrete crack openings (denoted as model I) is as follows: where δ(t) is the measured crack opening, δ T (t) is the temperature component, δ H (t) is the hydraulic component, δ θ (t) is the aging component, δ ε (t) is the residual, T i is the measured temperature of the ith measuring point, K is the number of temperature-measuring points, H is the water depth, θ � t/100, t is time, and a 0 , a i , b i , c 1 , and c 2 are all regression coefficients. Wang et al. [5] show that the temperature component in equation (1) does not consider a nonlinear temperature response. Since the concrete openings as a function of temperature is a complex nonlinear process, ignoring this nonlinearity will lead to underfitting. In contrast, including the environmental variables leading to nonlinearity alleviates this underfitting. Hence, we propose an improved statistical model of concrete crack openings, which is called statistical model considering nonlinear temperature component (denoted as model II) that accounts for the nonlinear temperature response [6][7][8]: where L represents the order of the temperature polynomial reflecting the influence of the nonlinear temperature response (L is generally a positive integer between 1 and 4 and is determined by fitting to measured values); a ip is the regression coefficient corresponding to the pth power of the measured temperature at the ith measurement point, and the other symbols have the same meaning as in equation (1). e regression coefficients of the statistical model are determined using regression analysis or optimization algorithms. e residual time series δ ε (t) is obtained by separating the measured data of the concrete crack openings.

Prediction of Residual Based on Chaos Theory
Under the action of water pressure and temperature variation, the stress field at the crack tip manifests singularity, which makes the concrete at the crack tip produce irreversible deformation. at is, the crack opening contains an irreversible part. Since the aging component cannot reflect the irreversible deformation accurately because of the limitations of its simple functions, the irreversible deformation is retained in the residual. In particular, when the crack opening suddenly changes due to overloading or drastic temperature variation, the hydraulic component, temperature component, and aging component of statistical model can hardly reflect this mutation of crack opening. In other words, the residual sequence contains a lot of information about the evolution of cracks [10] and some studies showed [9,14] that this residual sequence contains chaotic characteristics. erefore, chaos theory is used to analyze and predict the residual.

Phase-Space Reconstruction of Residual Sequence.
Previous studies [10,13] have shown that a 1D residual time series x(t) � {x(t i )|i � 1, 2, 3, . . ., N} can be reconstructed into the following form in m-dimensional phase space R m : where M � N − (m − 1)τ is the number of phase points, m > 0 is the integer embedding dimension, and τ is the delay time, which is usually a positive integer multiple of the sampling interval of the residual time series x(t). e first-phase point in equation (3) can be expressed as follows:

(4)
According to equations (3) and (4), the key of phasespace reconstruction lies in the determination of the delay time τ and the embedding dimension m.

Determination of Delay Time.
Selecting the optimized delay time t 0 usually involves the autocorrelation coefficient method [16], the mutual information method [17], or the AD method [18]. e autocorrelation coefficient method is not suitable for nonlinear systems [25] and imposes a subjective selection of the method to calculate the descent coefficient. Although the mutual information method overcomes the shortcomings of the autocorrelation coefficient method, its calculation is relatively complicated. Hence, this study uses the AD method [18] to select the optimized delay time; it is an efficient calculation and avoids information redundancy and the complete irrelevance of delay coordinates caused by improper delay time. Under the 2-norm, the average displacement S 2 (m, τ) is defined as follows:

Mathematical Problems in Engineering
When given different embedding dimensions (≥2), S 2 (m, τ) tends to be stable with increasing τ. Rosenstein et al. [18] suggested that optimized delay time can be determined according to the S 2 (m, τ) vs. τ graph under different embedding dimension m; that is, when the slope of S 2 (m, τ) vs. τ decays to 40% of the initial slope, the corresponding τ is the optimized delay time. Research shows that this method is reliable for small, noisy datasets and that the computation time is short.

Determination of Embedding Dimension.
Methods to calculate the minimum embedding dimension m e include the geometric invariant method, the pseudo-neighbor method, and the Cao method. However, the geometric invariant method requires a noise-free time series, which is barely possible in practical applications. Cao [20] proposed the pseudo-neighbor point method, but this method is very subjective when selecting the decision threshold. ey also proposed a practical method (called the "Cao method") for determining the minimum embedding dimension.
is method not only is robust against noisy data but also reduces the interference of subjective factors and adverse effects caused by slight reductions in data size. For this reason, this study adopts the Cao method to calculate the minimum embedding dimension.
When an optimized time delay is given, the following parameters are defined in m-dimensional phase space: a(i,m) is defined as follows: When the embedding dimension is incremented from m to m + 1, E1(m) is as follows: Cao et al. [20] report that, given a fixed-point attractor in the residual time series, E1(m) would stop changing once the embedding dimension m exceeds a certain value m 0 . In this case, m 0 + 1 is the minimum embedding dimension, which may be determined by a graphical analysis of E1(m) vs. m.
In actual calculations, it is difficult to determine whether E1(m) increases slowly or remains constant as m increases. For a group of random residual series, E1(m) increases with increasing m and converges for a deterministic residual series E1(m). Hence, we must introduce E2(m) to distinguish random residual series from deterministic residual series. E2(m) is given as follows: where For random residual series, E2(m) � 1 for any m; for deterministic residual series, E2(m) is a nonconstant function of m; that is, some m must exist so that E2(m) ≠ 1.
erefore, Cao et al. [20] suggested that E1(m) and E2(m) be calculated at the same time to determine the minimum embedding dimension and distinguish between random residual series and deterministic residual series. For notational convenience, E1(m) and E2(m) are denoted as E1 and E2, respectively.

Identification of Chaotic Characteristics.
e maximum Lyapunov exponent λ 1 is of important significance for determining the chaotic characteristics contained in residual time series. For phase space, λ 1 > 0 indicates that the system is chaotic. Rosenstein et al. [26] analyzed common methods of calculating the maximum Lyapunov exponent. For problems of unreliable calculation, intensive computation, and difficult-to-implement existing methods for small datasets, a practical method (called the "Rosenstein method" and also known as the "small-dataset method") was proposed for small datasets. e advantage of this method is that it is robust against the reconstruction dimension, time delay, and noise, thereby producing a more accurate maximum Lyapunov exponent. According to the theory of chaotic dynamics, the reciprocal λ − 1 1 of the maximum Lyapunov exponent is the maximum prediction time and can be used as a reliability indicator of a short-term prediction [27].
When adopting the Rosenstein method to calculate the maximum Lyapunov exponent, dj(k) is used in the phase space reconstructed by the optimized delay time and minimum embedding dimension to express the Euclidean distance between phase point j, X(t j ), and its initial nearest point X(t j ) (the meaning is the same as the nearest point in (5)) at moment kΔt: where Δt is the sampling interval and is usually taken as Δt � 1 for calculational convenience. e calculation result for d j (k) is expressed in matrix form as follows: d j (k) is expressed as [26].
Taking the logarithm of both sides of equation (13) gives the following: where C j is the initial bifurcation and its value determined for j is constant. y(k) is defined as follows: where 〈ln d j (k)〉 is the average over all j, that is, the average of the natural logarithm of the nonzero elements in row k of the matrix in equation (11).
In this case, we take the part of y(k) that is close to linear with respect to k according to equations (13) and (14) and use the least-squares method to obtain a linear fit. e slope of the line is λ 1 , that is, the maximum Lyapunov exponent.

Predicted Value of Residual Based on GA-BP Neural
Network. In the phase-space reconstruction of the residual time series, the "trajectory" in the reconstructed phase space is "dynamically equivalent" to the original system due to the concept of differential homeomorphism. erefore, there is a unique mapping relation: where P t is the prediction time, and the mapping relation F can be calculated by approximating all phase points in phase space. In this study, the nonlinear mapping relation is constructed using artificial neural networks with the whole domain method. e GA is a heuristic random evolutionary search algorithm, which has global search ability and can solve the error function without gradient information.
Previous studies demonstrated that the GA-BP neural networks established by GAs provide better fits to nonlinear functions and more accurate short-term predictions of typical chaotic time series, even when compared with the predictions of BP neural networks [23,28]. erefore, this study uses GAs to optimize BP neural networks and thereby improves simple BP neural networks to attain the global maximum.
Setting the fitness function is quite important with GAs. In this study, we use the fitness function [29]as follows: which shows that a large R 2 and small root-mean-square error (RMSE) indicate a better fit. e optimal fitness is attained when the fitness becomes stable as a function of the evolutionary generations. e equations for R 2 and the RMSE are as follows: where N t is the number of samples in the test set, δ t is the measured value, δ ∧ is the predicted value, and δ t is the average measured value.
In GA-BP neural networks, the input vector is as follows: and the output vector is as follows: e separated residual time series is trained to optimize the neural network model, which gives the weight and threshold of the GA-BP neural network model. e residual time series is then predicted.

Hybrid Model for Predicting Concrete Crack Opening.
Model II and the GA-BP neural network (15) are integrated to establish the improved hybrid model for predicting concrete crack openings (denoted as "hybrid model"). Mathematically, the model is expressed as follows: Mathematical Problems in Engineering where δ T (t + P t ), δ H (t + P t ), and δ θ (t + P t ) predict the temperature, hydraulic, and aging components of P t , respectively; x(t + P t ) is the predicted value of the residual based on GA-BP neural network. Figure 2 presents a flow chart of the procedure used by the hybrid model for predicting concrete crack openings.

Crack Distribution.
Hubei Province has undertaken a large water control project designed for power generation and shipping, irrigation, aquaculture, and tourism. e project consists of a regulating sluice located in the main river with an earth-rockfill dam in the main riverbed, an earth-rockfill dam in the Gucheng Section, a powerhouse in the old river, a ship lock, a concrete gravity dam, and an earth-rockfill dam on both sides of the embankment.
On-site inspection has revealed three cracks at the top of the upstream traffic bridge of the gate storehouse on the left bank of the sluice (see Figure 3). Two of these cracks run upstream and downstream and stop at the lowest part of the breast wall. Meanwhile, a long-term calcium precipitation trace appears on the vertical surface of the retaining wall upstream of the gate storehouse and extends to the floor of the upstream wall of the gate storehouse, which is located between the two penetrating cracks at the top of the sluice and threatens the normal operation of gate storehouse.

Monitoring by Crack Meter.
To diagnose the evolution of the concrete crack in the regulating sluice, nine sets of bidirectional vibrating string crack meter groups were installed to monitor the openings and misplacement momentum (see Figure 4). Mkj-1 and Mkj-2 were installed on the upstream side of the sluice, Mkj-3-Mkj-5 were installed near the top of the sluice, Mkj-6 and Mkj-7 were installed near the upstream surface of the sluice, and Mkj-8 and Mkj-9 were installed on the top of the sluice near the downstream section. From June 12 to 19, 2020, nine sets of bidirectional crack meters were installed and debugged. Figure 5 shows the measured crack openings and temperature lines of crack meter Mkj-3.
An obvious negative correlation appears between the crack opening at the sluice and the temperature. When the temperature increases, the crack closes and vice versa. Next, we use the measured crack opening in the hybrid model for concrete crack openings.

Determination of L in Model II.
e various measuring points may give different results due to the nonlinear temperature response. To determine L in model II according to the principle mentioned above, we use nine temperature measurements by the crack meter. Using equation (2), we establish model II for concrete crack openings for L � 1, 2, 3, 4, 5.
e model regression coefficient corresponding to different results for L is obtained by regression analysis. Rsquare and RMSE of different improved models were calculated. e results are shown in Figure 6.
When L increases from 1 to 3, both R 2 and RMSE improve significantly; when L ≥ 3, the two parameters remain stable. From the point of view of statistical models, this shows that temperature exerts a nonlinear effect on the concrete crack opening. erefore, when model II with nonlinear temperature effects is included, R 2 � 0.891 and RMSE � 0.04571. Table 1 lists the regression coefficients obtained through regression analysis. Using these, model II for concrete crack openings is as follows. Model I considers only the linear response to temperature, as highlighted by the ellipses in Figure 7(a). Note that the fit near the extreme values is not ideal. In Figure 7(b), the nonlinear response to temperature is considered in model II, which makes the fit closer to the measured trajectory. is is seen in the decrease in the residual error of model II (Figure 7(c)).

Comparison with Other Models.
Notably, limitations of the statistical models shown in Figures 7(b) and 7(c) still result in periods where the fit is poor, which will affect the accuracy of the models' predictions to some extent. Hence, the prediction methods based on the statistical model may still be improved. For this reason, we apply chaos theory to further analyze the residual, which leads to the more accurate model.

Determination of Embedding Dimension and Delay
Time.
e residual time series from model II for crack openings is obtained by separating the measured values. From the residual time series, 1825 samples from 22:00 on June 23, 2020, to 7:00 on February 11, 2021, are selected as the training dataset, and 24 samples from 10:00 on February 11, 2021, to 7:00 on February 14, 2021, are selected as a posteriori prediction dataset. Figure 8 shows S 2 (m, τ) as a function of delay time τ constructed using the AD plotting method, which is based on phase-space reconstruction theory. Figure 9 plots E1 and E2 (determined by the Cao method) as a function of embedding dimension m.
When the embedding dimension is constant, the slope of S 2 (m, τ) vs. τ stabilizes as τ increases. e delay time is optimal when the slope of S 2 (m, τ) vs. τ decreases to 40% of its initial value, which gives τ � 2 as the optimal delay time. For a given optimal delay time, E1 stabilizes with increasing embedding dimension, and E2 ≠ 0. erefore, when E1 stabilizes, the corresponding embedding dimension (i.e., the minimum embedding dimension) is eight.

Phase-Space Reconstruction of Residual Time Series and Determination of Prediction
Step. Given the embedding dimension and delay time, phase-space reconstruction is applied to the 1D residual time series via equation (3), as per the principle mentioned above. e maximum Lyapunov exponent is calculated via the Rosenstein method, indicating that the system contains chaotic components. e reciprocal of the system is then calculated, giving a maximum prediction step of 40. erefore, the posterior step is 24 (i.e., 24 samples at 7:00 on February 11, 2021, and 7:00 on February 14, 2021), and the duration is three days.

Hybrid Model
Considering nonlinear temperature response Observation date of hydraulic structure macro crack opening Accuracy satisfied? Mathematical Problems in Engineering 2020, to 7:00 on February 11, 2021, to train the GA-BP neural networks. e structure of the BP neural network is 8-12-8, the GA's maximum number of generations is 30, the population is 40, the crossover probability is 0.8, and the variation probability is 0.1. Figure 10 shows the resulting fit.
When the fitness stabilizes (i.e., for >5 generations) and satisfies the precision requirements of neural networks, the optimized BP neural network is created using the weights and threshold parameters output at this time to train and predict the reconstructed residual time series.

Hybrid Model.
Based on the length of the posteriori predicted dataset of residual time series, the prediction is made starting at 10:00 on February 11, 2021, and the prediction step is 24 (three days). e residual time series predicts that x(t + 1), x(t + 2), . . ., x(t + 24) are obtained, and the hybrid model is as follows: where δ(t + P t ) is the predicted concrete crack opening at time P t , where P t � 1, 2, . . ., 24.

Prediction Accuracy.
Mean square error (MSE), mean absolute error (MAE), and normalized mean square error (NRMSE) are introduced to evaluate the prediction accuracy of the hybrid model. ese are given by: where σ is the standard deviation of the test set, N t is the number of samples in the test set, δ t is the measured concrete crack opening, and δ ∧ is the predicted concrete crack opening. e closer the three evaluation indexes are to zero, the more accurate is the prediction.
According to the posteriori prediction of the residual time-series dataset, equation (23) gives the prediction of the hybrid model for concrete crack openings starting from 10: 00 on February 11, 2021, with a prediction step of 24 (three days). Figure 11 compares this prediction with those of model I and model II. Table 2 gives the evaluation indicators for different models.
Unlike the values predicted by the conventional statistical model, those predicted by the statistical model considering nonlinear temperature component, which considers the nonlinear response to temperature, are closer to the measured values, although the improvement is minimal near the extreme values. However, the difference between the predicted values and the measured values near the extrema is significantly reduced when using the hybrid model. Compared with model I, the MSE, MAE, and NRMSE of the hybrid model decreased by 87.0%, 68.5%, and 63.9%, respectively. Compared with model II, the MSE, MAE, and NRMSE of the hybrid model decreased by 73.4%, 51.4%, and 48.4%, respectively.
In summary, a positive Lyapunov exponent indicates that concrete crack openings have chaotic components, so the hybrid model makes significantly improved predictions of concrete crack openings.

Dynamic Prediction Based on Hybrid Model of Concrete
Crack Openings. Early warning is a dynamic process in real engineering because concrete crack openings evolve over time. Given new monitoring data, the parameters of the model must be updated and the dynamic prediction must be redone. To verify the reliability of the hybrid model, ten consecutive dynamic predictions were made, each with a time step of 24 (three days), and MSE, MAE, and NRMSE of model I, model II, and the hybrid model were statistically analyzed. Based on these three metrics, Table 2 gives the      calculated improvement for the hybrid model relative to model I and model II, and Figure 12 shows the box plots for a comparative analysis.
For the hybrid model, the three evaluation metrics are significantly improved compared with the two statistical models. e improvement for the hybrid model relative to model I is 45.9%-70%, and the maximum improvement in the MSE exceeds 80%. e average improvement for the hybrid model relative to model II is 23.4%-40.1%, and the maximum improvement in the MSE exceeds 70%. is indicates that the hybrid model provides more accurate predictions of concrete crack openings than does model II.

Conclusions
In this study, the conventional statistical model for concrete crack opening is improved by considering the nonlinear temperature factor. Using chaos theory and phase-space reconstruction theory, a hybrid model that combines the improved statistical model with a GA-BP neural network residual model is finally developed to predict concrete crack openings in hydraulic structures.
It is found that concrete crack openings are affected by the nonlinear temperature response. However, the conventional statistical model cannot accurately reflect the temperature component of the system. e improved statistical model makes a significant improvement in this area by introducing the second-and third-order temperature components. e residual series of the statistical model considering nonlinear temperature component contains useful information of crack evolution; the presence of a chaotic component is demonstrated by calculating the maximum Lyapunov exponent. e phase-space reconstruction of the residual series is combined with a GA-BP neural network and the statistical model considering nonlinear temperature component to establish a hybrid model to predict the residuals. e predictions of the hybrid model are more accurate than those of the statistical model. ese results of the case study show that the hybrid model produces more reliable short-term dynamic predictions of concrete crack openings, making it of practical use for generating early warnings of concrete crack openings.
Data Availability e data in this article were calculated and analyzed by the authors.  Figure 12: Dynamic prediction of hybrid model of concrete crack openings compared with that of statistical models. (a) Improvement due to hybrid model relative to model I. (b) Improvement due to hybrid model relative to model II.

Supplementary Materials
Supplementary 1. Concise description for Figure 5. Figure 5: concrete crack opening (blue) and temperature (red). Column A is time. Column B is concrete crack opening. Column C is measured temperature. Supplementary 2. Concise description for Figure 6.     Figure 12: dynamic prediction of hybrid model of concrete crack openings compared with that of statistical models. Column A is times of dynamic prediction. Column B is improvement of MSE due to hybrid model relative to model I. Column C is improvement of MAE due to hybrid model relative to model I. Column D is improvement of NRMSE due to hybrid model relative to model I. Column E is improvement of MSE due to hybrid model relative to model II. Column F is improvement of MAE due to hybrid model relative to model II. Column G is improvement of NRMSE due to hybrid model relative to model II. (Supplementary Materials)