CEEMDAN algorithm
The CEEMDAN algorithm24 evolved from the Empirical Mode Decomposition (EMD), Ensemble EMD (EEMD), and Complementary Ensemble EMD (CEEMD) algorithms. It effectively suppressed the mode-mixing phenomenon in EMD, generated reconstructed signals with lower residual noise than EEMD, and resolved the misalignment or errors caused by inconsistent decomposition results across subsequence groups in CEEMD. The theoretical steps of the CEEMDAN algorithm are as follows:
(1) Add Gaussian white noise of the same length to the original signal sequence \(f
(1)
Where, \({\varepsilon _0}\) is the weight coefficient of Gaussian white noise, and \({\delta _i}
(2)
(3) Add Gaussian white noise to the residual component obtained in the k-th stage after decomposition, and continue to apply EMD for further decomposition. This process can be represented as:
(4) Repeat step (3) until the residual component becomes a monotonic signal and can no longer be decomposed, at which point the iteration ends. Ultimately, the original signal sequence is decomposed into N modal components and a residual component.
VMD algorithm
VMD is an adaptive, completely non-recursive modal variation signal processing method25. Its fundamental principle is to decompose a signal into multiple components with fixed bandwidths, where each component corresponds to a specific frequency and amplitude within the signal. By optimizing a variational regularization function, VMD can adaptively match the optimal center frequency and limited bandwidth for each mode, thereby achieving effective separation of intrinsic mode functions (IMF), frequency domain partitioning of the signal, and obtaining effective decomposition components of the given signal26.
In addition to decomposing one-dimensional signals, VMD can also be extended to the decomposition of multi-dimensional signals. The research on multi-dimensional VMD aims to extract spatial and temporal features from multi-dimensional signals and apply them to fields such as image processing27,28 and video processing13providing new methods and insights for multi-dimensional data analysis.
The implementation formula is as follows:
$$\mathop {\hbox{min} }\limits_{{\{ {u_k},{\omega _k}\} }} \{ \sum\nolimits_{k} {\left\| {{\partial _t}[\delta
(7)
$$s.t.\sum\nolimits_{k} {{u_k}} =f
(8)
Where, \({u_k}\) is the mode variable of the kth decomposition, \({\omega _k}\) is the center frequency of \({u_k}\), \(\delta
(9)
Where, the penalty factor and Lagrange operator are intended to maintain the strictness of the constraints and ensure the accuracy of signal reconstruction. By continuously iterating using the alternate direction method of multipliers, the optimal solution of formula (9) is obtained. Therefore, the variables \({u_k}\) and \({\omega _k}\) can be updated according to formulas (10) and (11):
$$\omega _{n}^{{n+1}}=\frac{{\int_{0}^{\infty } {\omega {{\left| {{u_k}(\omega )} \right|}^2}d\omega } }}{{\int_{0}^{\infty } {{{\left| {{u_k}(\omega )} \right|}^2}d\omega } }}$$
(11)
Where, n represents the iteration number.
TTAO algorithm
The TTAO algorithm29 is based on the principle of similar triangles. During the iterative process, it searches using the three vertices and one interior point of the triangular topology unit. Optimization is achieved through aggregation within the triangular topology unit and between different triangular topology units. The optimization process consists of three main stages: the construction of triangular topology units, global aggregation, and local aggregation.
(1) Construction of triangular topology units: The number of individuals N can be divided into \([{N \mathord{\left/ {\vphantom {N 3}} \right. \kern-0pt} 3}]\) triangular topology units, where \([\cdot ]\) represents the floor value. The process is as follows:
$$\overrightarrow {{X_{i,1}}} ={r_0} \times (\overrightarrow {UB} – \overrightarrow {LB} )+\overrightarrow {LB}$$
(12)
$$\overrightarrow {{X_{i,2}}} =\overrightarrow {{X_{i,1}}} +l \times \overrightarrow {f(\overrightarrow \theta )}$$
(13)
$$\overrightarrow {{X_{i,3}}} =\overrightarrow {{X_{i,1}}} +l \times \overrightarrow {f(\overrightarrow {\theta +\frac{\pi }{3}} )}$$
(14)
$$\overrightarrow {{X_{i,4}}} ={r_1} \times \overrightarrow {{X_{i,1}}} +{r_2} \times \overrightarrow {{X_{i,2}}} +{r_3} \times \overrightarrow {{X_{i,3}}}$$
(15)
Where, \(\overrightarrow {{X_{i,1}}}\), \(\overrightarrow {{X_{i,2}}}\) and \(\overrightarrow {{X_{i,3}}}\) represent the three vertices of the ith triangular topology unit, and \(\overrightarrow {{X_{i,4}}}\) indicates a randomly chosen interior vertex within the ith triangular topology unit. \(\overrightarrow {UB}\) and \(\overrightarrow {LB}\) are the upper and lower bounds of the search space. l represents the size of the triangular topology unit. Additionally, \(\overrightarrow {f\left( {\overrightarrow \theta } \right)}\) and \(\overrightarrow {f\left( {\overrightarrow {\theta +\frac{\pi }{3}} } \right)}\) denote the direction vectors of the other two edges guided by the first vertex.
(2) Global aggregation: This process represents the exploration phase of the algorithm. By collecting information from excellent individuals within different triangular topology units, new feasible solutions are generated. Specifically, at kth iteration, the optimal vertex \(\overrightarrow {X_{{i,best}}^{k}}\) of each triangular topology unit interacts with the optimal vertex \(\overrightarrow {X_{{rand,best}}^{k}}\) of a randomly selected topology unit to aggregate and form a new feasible solution \(\overrightarrow {X_{{i,new1}}^{{k+1}}}\), as shown in the following formula:
$$\overrightarrow {X_{{i,new1}}^{{k+1}}} ={r_4} \times \overrightarrow {X_{{i,best}}^{k}} +(1 – {r_4}) \times \overrightarrow {X_{{rand,best}}^{k}}$$
(16)
Next, the fitness of the new solution is compared with the fitness of the optimal or suboptimal solution at kth iteration, in order to update the optimal solution, that is:
$$\begin{array}{*{20}{c}} {\overrightarrow {X_{{i,best}}^{{k+1}}} =\overrightarrow {X_{{i,new1}}^{{k+1}}} }&{{f_{\overrightarrow {X_{{i,new1}}^{{k+1}}} }}<{f_{\overrightarrow {X_{{i,best}}^{k}} }}} \\ {\overrightarrow {X_{{i,sbest}}^{{k+1}}} =\overrightarrow {X_{{i,new1}}^{{k+1}}} }&{{f_{\overrightarrow {X_{{i,new1}}^{{k+1}}} }}<{f_{\overrightarrow {X_{{i,sbest}}^{k}} }}} \end{array}$$
(17)
(3) Local aggregation: This process is the development phase of the algorithm. Within each triangular topology unit, perturbations are applied to the optimal solution based on the differences between the optimal and suboptimal solutions, in order to prevent the optimal individuals from getting trapped in local optima, as shown in Eq. (18).
$$\overrightarrow {X_{{i,new2}}^{{k+1}}} =\overrightarrow {X_{{i,best}}^{{k+1}}} +\alpha \times (\overrightarrow {X_{{i,best}}^{{k+1}}} – \overrightarrow {X_{{i,sbest}}^{{k+1}}} )$$
(18)
Where, the \(\alpha\) is continuously reduced to allow the algorithm to gradually approach the optimal solution. To improve convergence, the fitness values of the two vertices before and after the local development are compared to determine the update position. If the new individual is better than the original individual, the position is updated; otherwise, no update is performed, that is:
$$\overrightarrow {X_{{i,best}}^{{k+1}}} =\left\{ {\begin{array}{*{20}{c}} {\overrightarrow {X_{{i,new2}}^{{k+1}}} }&{{f_{\overrightarrow {X_{{i,new2}}^{{k+1}}} }}<{f_{\overrightarrow {X_{{i,best}}^{{k+1}}} }}} \\ {\overrightarrow {X_{{i,best}}^{{k+1}}} }&{{\text{oth}}erwise} \end{array}} \right.$$
(19)
TVMD algorithm
The VMD algorithm improves the predictive accuracy of data by decomposing it to reduce non-stationarity and non-linearity. However, if the number of modes for decomposition is too few, it may leave behind a residual with high complexity, which cannot guarantee prediction accuracy. Conversely, if the number of modes is too high, it can lead to over-decomposition of the data. Therefore, in the VMD algorithm, the number of decomposed modes k and the penalty factor \(\alpha\) are interdependent and jointly affect the decomposition results. Thus, obtaining the optimal parameter combination \([k,\alpha ]\) based on the characteristics of the input signal is key to achieving adaptive decomposition through VMD.
Therefore, this paper employed the TTAO optimization algorithm to optimize the optimal parameter combination of the VMD decomposition algorithm, referred to as the TVMD decomposition algorithm. This algorithm uses the minimum envelope entropy \(\hbox{min} {E_e}\) as the fitness function. After the signal undergoes VMD decomposition, the more noise contained in the subsequences, the greater the envelope entropy value will be. The calculation formula for envelope entropy is as follows:
$$\left\{ {\begin{array}{*{20}{c}} {{E_e}= – \sum\limits_{{j=1}}^{N} {{b_j}\lg {b_j}} } \\ {{b_j}=a(j)/\sum\limits_{{j=1}}^{N} {a(j)} } \end{array}} \right.$$
(20)
Where, \(a(j)\) represents the result obtained by applying Hilbert transform to the IMF component of the original data after VMD decomposition. \({b_j}\) is the normalized form of \(a(j)\).
The steps of the proposed TVMD decomposition algorithm are as follows:
(1) Initialize the TTAO parameters, set the population size, number of iterations, and the value ranges for K and \(\alpha\).
(2) Obtain the fitness value using Eq. (20), and update the positions of the triangular topology units according to Eqs. (17) and (19) until the maximum number of iterations is reached, resulting in the corresponding optimal parameter combination \([k,\alpha ]\).
(3) Perform VMD decomposition on the original precipitation data based on the optimal parameter combination, outputting K intrinsic mode function components and the residual component.
PO algorithm
The Parrot Optimization Algorithm30 is a meta heuristic optimization algorithm proposed in 2024, which solves the optimal parameters by simulating the four key behavioral characteristics of parrots. The solution process is as follows:
(1) Population Initialization: The algorithm initializes a set of candidate solutions as the parrot population, with each parrot representing a potential solution. Assume the population size is N, the maximum number of iterations is \(Ma{x_{iter}}\) and the lower and upper bounds of the search space are lb and ub. The initial positions of the parrots are:
$$X_{w}^{0}=lb+rand(0,1)\cdot (ub – lb)$$
(21)
Where, \(rand(0,1)\) represents a random number in the range \([0,1]\), and \(X_{w}^{0}\) denotes the initial position of the wth parrot.
(2) Foraging behavior: Parrots estimate the approximate location of food by observing its position or the position of its owner, and then fly towards their respective positions, where the position movement follows the following equation:
$$X_{w}^{{d+1}}=(X_{w}^{d} – {X_{best}})\cdot Levy(dim)+rand(0,1)\cdot {(1 – \frac{t}{{Ma{x_{iter}}}})^{\frac{{2d}}{{Ma{x_{iter}}}}}}\cdot X_{{mean}}^{d}$$
(22)
Where, \(X_{w}^{d}\) represents the current position, \(X_{w}^{{d+1}}\) denotes the updated position, \({X_{best}}\) indicates the best position found so far and the master’s current position, \(Levy(dim)\) denotes the Levy distribution, which describes the parrot’s flight, d indicates the current iteration number, and \(X_{{mean}}^{d}\) represents the average position of the current population, i.e.:
$$X_{{mean}}^{d}=\frac{1}{N}\sum\limits_{{k=1}}^{N} {X_{k}^{d}}$$
(23)
(3) Staying Behavior: The parrot suddenly flies to any part of the owner’s body and remains still for a period of time. This process can be represented as:
$$X_{w}^{{d+1}}=X_{w}^{d}+{X_{best}}\cdot Levy(dim)+rand(0,1)\cdot ones(1,dim)$$
(24)
Where, \(ones(1,dim)\) represents an all-ones vector of dimension dim.
(4) Communication Behavior: This behavior involves close interaction within the flock, including both flying towards and not flying towards the group. It is simulated by calculating the average position of the population and adjusting candidate solutions accordingly to promote information sharing and collaboration. This process can be represented as:
$$X_{w}^{{d+1}}=\left\{ {\begin{array}{*{20}{c}} {0.2rand(0,1)\cdot (1 – \frac{d}{{Ma{x_{iter}}}})\cdot (X_{w}^{d} – X_{{mean}}^{d})}&{P \leqslant 0.5} \\ {0.2rand(0,1)\cdot \exp ( – \frac{d}{{rand(0,1)Ma{x_{iter}}}})}&{P>0.5} \end{array}} \right.$$
(25)
Where, when \(P \leqslant 0.5\) represents an individual joining a parrot group for communication, when \(P>0.5\) represents the process of the individual immediately flying out after communication.
(5) Fear behavior: Individuals usually keep a distance from unfamiliar individuals and seek a safe environment together with their owner. This behavior avoids excessive concentration of candidate solutions through a repulsion mechanism, maintaining population diversity. This process can be represented as:
$$\begin{gathered} X_{w}^{{d+1}}=X_{w}^{d}+rand(0,1)\cdot \cos (0.5\pi \cdot \frac{d}{{Ma{x_{iter}}}})\cdot ({X_{best}} – X_{w}^{d}) \hfill \\ – \cos (rand(0,1)\cdot \pi )\cdot {(\frac{d}{{Ma{x_{iter}}}})^{\frac{2}{{Ma{x_{iter}}}}}}\cdot (X_{w}^{d} – {X_{best}}) \hfill \\ \end{gathered}$$
(26)
IPO algorithm
(1) Population initialization with chaotic reverse learning strategy.
Traditional population initialization methods rely heavily on random number generation, which can lead to uneven distribution of individuals in the search population. To maintain population diversity and ensure that the initial population is as evenly distributed as possible, this paper introduced an initialization strategy based on chaotic reverse learning31. This approach helps accelerate the convergence speed of the algorithm. The steps of this strategy are: first, use the Cat chaotic sequence to generate N initial solutions \({X_i}\). For each initial solution, generate the corresponding reverse solution using the following method:
$$O{P_i}=rand(0,1)\cdot (X_{{\hbox{min} }}^{d}+X_{{\hbox{max} }}^{d}) – {X_i}$$
(27)
Where, \(X_{{\hbox{min} }}^{d}\) and \(X_{{\hbox{max} }}^{d}\) represent the minimum and maximum values of the dth dimension vector among all initial solutions.
Finally, the initial solutions with the reverse solutions was combine to sort in ascending order based on their fitness values. The top N solutions with the best fitness values was selected to form the initial population.
(2) Improving Nonlinear Convergence Factor Strategy.
When \(P \leqslant 0.5\) represents an individual joining a parrot group for communication, when \(P>0.5\) represents the process of the individual immediately flying out after communication. Therefore, the value of P is closely related to the communication behavior in the parrot optimization algorithm. However, in the traditional parrot algorithm, the value of P is random, which does not reflect the changes in the optimization algorithm during the iterative process. Therefore, this paper proposed an improved nonlinear formula:
$$P=rand(0,1)\cdot \frac{{Ma{x_{iter}} – d}}{{Ma{x_{iter}}}}$$
(28)
(3) Cauchy-Gaussian variation.
The variation strategy can prevent the algorithm from getting trapped in local optima and also maintain the diversity of the population. To reduce the probability of the parrot algorithm falling into local optima, this paper introduced the Cauchy-Gaussian variation operator32. This mutation operator combines the Cauchy variation operator33 and the Gaussian variation operator34. It allows global search during the early stages of population optimization and local search during the later stages of iteration, thereby significantly enhancing its optimization capability. The expression for this operator is:
$$X_{{new}}^{d}=X_{{best}}^{d} \times [1+{\beta _1}Cauchy(0,1)+{\beta _2}Causs(0,1)]$$
(29)
Where, \(X_{{best}}^{d}\) represents the optimal position of the parrot population in the dth iteration. \(X_{{new}}^{d}\) is the new position generated from the optimal position in the dth iteration using the Cauchy-Gaussian variation strategy. \(Cauchy(0,1)\) and \(Causs(0,1)\)are the random number following a Cauchy distribution and a Gaussian distribution. \({\beta _1}=1 – {\raise0.7ex\hbox{$d$} \!\mathord{\left/ {\vphantom {d {Ma{x_{iter}}}}}\right.\kern-0pt}\!\lower0.7ex\hbox{${Ma{x_{iter}}}$}}\) and \({\beta _2}={\raise0.7ex\hbox{$d$} \!\mathord{\left/ {\vphantom {d {Ma{x_{iter}}}}}\right.\kern-0pt}\!\lower0.7ex\hbox{${Ma{x_{iter}}}$}}\).
BiLSTM network
The LSTM network model35 is proposed based on recurrent neural network36 and is renowned for its excellent ability to process sequential data. The LSTM model introduces gate control units (forget gate, input gate, and output gate) and memory cell states, which address the issues of gradient vanishing and gradient explosion that often occur during long sequence training37.
LSTM network can only encode time series data in a forward direction and cannot learn patterns from both forward and backward information in the sequence. In contrast, BiLSTM neural networks consist of a combination of forward and backward LSTM, allowing them to consider the impact of forward and backward time series data on the current state. Therefore, BiLSTM typically provides better prediction accuracy than LSTM38.
