Multi-Delay Identification of Rare Earth Extraction Process Based on Improved Time-Correlation Analysis

The rare earth extraction process has significant time delay characteristics, making it challenging to identify the time delay and establish an accurate mathematical model. This paper proposes a multi-delay identification method based on improved time-correlation analysis. Firstly, the data are preprocessed by grey relational analysis, and the time delay sequence and time-correlation data matrix are constructed. The time-correlation analysis matrix is defined, and the H∞ norm quantifies the correlation degree of the data sequence. Thus the multi-delay identification problem is transformed into an integer optimization problem. Secondly, an improved discrete state transition algorithm is used for optimization to obtain multi-delay. Finally, based on an Neodymium (Nd) component content model constructed by a wavelet neural network, the performance of the proposed method is compared with the unimproved time delay identification method and the model without an identification method. The results show that the proposed algorithm improves optimization accuracy, convergence speed, and stability. The performance of the component content model after time delay identification is significantly improved using the proposed method, which verifies its effectiveness in the time delay identification of the rare earth extraction process.


Introduction
The rare earth extraction process includes dozens or even hundreds of extraction tanks. The mixing speed and time of each group of agitators are different, which affects the reaction and transmission time of materials, leading to multi-delay, so a large amount of data cannot be effectively utilized. Current modeling studies of the rare earth extraction process do not consider time delay or take it as a constant [1][2][3], resulting in a particular gap between the model and the extraction site. Therefore, it is significant to study how to identify multi-delay.
A series of solutions have been proposed for the problem of time delay identification. The step response method was used in the time delay identification of systems [4,5], but this method is susceptible to noise and requires a filter to remove high-frequency noise. Rad et al. [6] estimated the time delay of input and output signals based on the crosscorrelation function, which cannot reflect their specific relationship, so the results are not ideal. Some scholars used a recursive least squares algorithm to identify the system with time delay [7][8][9], which is assumed to be known and may not be established in engineering practice. Neural networks were used to identify time delay, which have the problem of long training times and easily falling into local optima [10][11][12]. Liu et al. [13] developed a compressed sensing recovery algorithm for the multiple input single output finite impulse response systems with unknown time delay, but it was not easy to choose the optimal threshold. Chen et al. [14] proposed an effective identification model based on the Bayesian theorem for systems with unknown time delay. Wang et al. [15] proposed a parameter identification method of a fractional-order time delay system based on the Legendre wavelet, which reduces the effect of noise on the accuracy of parameter identification. Hofmann et al. [16] proposed an offline time-delay identification strategy based on falling film evaporator pilot plant experiments and obtained good results in both validation experiments, with and without evaporation. Prasad et al. [17] used fractional-order modeling technique to identify the parameters of Hammerstein structured nonlinear systems with discontinuous asymmetric (two segment piecewise-linear with a dead-zone) nonlinearity and input time delay. Li et al. [18] designed a discrete-time robust adaptive estimator to identify the time delay and sandwich system parameters. Meanwhile, they reconstructed the observation and augmented data to obtain the explicit expression of the delay parameter. To achieve an effective system control strategy and accurate response prediction, Liu et al. [19] proposed a new method to identify the parameters of linear time-delay differential systems by analyzing the frequency domain response of complex systems. To solve the influence of time delay on HVAC systems, Li et al. [20] introduced transfer entropy and proposed a model-free identification method based on the information theory framework. Ni et al. [21] studied the parameter estimation problem for a class of linear time-delay systems. Based on the frequency responses and harmonic balance methods and by means of the gradient search, a two-stage stochastic gradient and gradient-based iterative algorithm was developed by using the collected data under the sinusoidal excitation. The maximum likelihood method [22], variable structure observer [23], particle swarm optimization [24], and other methods have also been applied.
Industrial processes have become complex with the rapid development of science and technology, resulting in multi-delay, and the abovementioned methods have been unable to meet practical requirements. Xie et al. [25] improved the genetic algorithm, applied it to the identification of multi-delay, and obtained their optimal estimate, but the method requires an accurate system model. Huang et al. [26] proposed an improved cross-correlation function method and realized multi-delay identification of the alumina carbon separation process. Wang et al. [27] proposed a trend similarity analysis method and realized the multi-delay identification of the hydrocracking process. All the abovementioned methods have the problems of high computational redundancy and time consumption. This paper draws on the successful application of the time-correlation analysis method [28] in the alumina carbon separation and evaporation processes. We propose an improved time-correlation analysis method for the rare earth extraction process. Based on field data and the improved method, the multi-delay identification problem is transformed into an unconstrained integer optimization problem without changing the time delay relationship. Since the discrete state transition algorithm [29] can effectively solve the unconstrained integer optimization problem [30], we adopt the improved algorithm to solve it. Experimental comparison and analysis show that the improved time-correlation analysis method has high speed, high accuracy, and good stability, and is suitable for the multi-delay identification of the rare earth extraction process.

Improvement of Time-Correlation Analysis Method
Time-correlation analysis is a time delay identification method based on the relationship between data sequences, which has the advantage of high efficiency. This paper improves its shortcomings in data preprocessing and selection.

Grey Relational Analysis
Grey relational analysis (GRA) is derived from the grey system theory in system science [31]. Its basic idea is to judge the tightness of sequence connections according to the similarity between the geometric shapes of sequence curves. The closer the curves, the more significant the correlation between sequences.
Compared with traditional multi-factor analysis methods (such as canonical correlation analysis and multiple linear regression), this method has lower data requirements and less computational burden. It is suitable for quantitative analysis of the dynamic develop-ment process and hence can be used to analyze the correlation between a large amount of process data obtained from a production site and the content of rare earth elements.
Let the time base of a critical process variable in N work units be d = [d 1 , d 2 , . . ., d i , . . ., d N ], where i = 1, 2, . . ., N, d i = τ i /T. T is the sampling period, the delay sequence is Γ = [τ 1 , τ 2 , . . ., τ i , . . ., τ N ], and τ i is the time delay of the ith work unit. The content of n rare earth elements and m process variables is obtained by k times sampling. The element component content U j (t) = [u 1 (t), u 2 (t), . . ., u j (t), . . ., u n (t)] are used as the reference sequence in the correlation analysis, where 1 ≤ j ≤ n, 1 ≤ t ≤ k, and u j (t) represents the jth rare earth element component content. The process variable data E l (t) = [e 1 (t), e 2 (t), . . ., e l (t), . . ., e m (t)] are used as the comparison sequence, where 1 ≤ l ≤ m, 1 ≤ t ≤ k, and e l (t) represents the lth process variable data.
The original data are normalized to eliminate the influence of different dimensions on the results. Standard processing methods include initialization and averaging. This paper adopts averaging, where U j (t) is the processed reference sequence data, and E l (t) is the processed comparison sequence data. The correlation coefficient is used to express the degree of closeness between the index values of the comparison and reference sequence in grey relational analysis. The higher the value, the greater the degree of proximity, where ξ lj (t) is the correlation coefficient of the lth characteristic variable corresponding to the content of the jth rare earth element component, and ρ ∈ [0, 1] is the resolution coefficient, which we take as 0.5. According to the correlation coefficient, the correlation degree between each process variable and the content of rare earth elements can be obtained as The correlation degree is sorted from large to small. If r 11 < r 21 , then the correlation degree between the comparison sequence e 2 (t) and the content of the first rare earth element component is greater than that of comparison sequence e 1 (t).
where e 0, * and e i, * are the time series of the inlet process variables and the ith unit outlet process variables, respectively. F ≥ ∑ N 1 d i , so that the data in the time-correlation matrix contain information about the entire process cycle, from inlet to outlet.
The matrix E describes the degree of correlation between the data corresponding to the time sequences. Multiple time sequences are described by the time-correlation analysis matrix, where cov(E) is the covariance matrix of E, whose standard deviation of column i is σ i . R is a time-correlation analysis matrix, which reflects the correlation degree of multiple groups of sequences under multi-delay. The H ∞ norm can quantify it and is expressed as A maximum value of the H ∞ norm indicates the maximum correlation between each time sequence in the data matrix. At this time, the corresponding delay sequence Γ consists of the multi-delay to be solved.
According to the extraction production site experience, each unit's time delay during operation will fluctuate within a fixed range. Based on this, the time-base value range of key process variables can be determined as d i ∈ [d imin , d imax ]. Thus, the time delay identification problem of Equation (6) can be transformed into an unconstrained integer optimization problem, In summary, the steps of the improved method of combining grey correlation analysis with time-correlation analysis are as follows.
Step 1: Based on the original data of the rare earth extraction process, the multi-delay sequence Γ is constructed; Step 2: Grey correlation analysis is used to identify key process variables and construct a time-correlation data matrix according to Equation (4); Step 3: According to Equation (5), the time-correlation analysis matrix is defined, and the correlation degree of the data sequence is quantified by the H ∞ norm; Step 4: The time delay identification results are obtained using Equation (7). There are many methods to solve the time delay identification. The discrete state transition algorithm (DSTA) has been successfully applied to typical discrete optimization problems such as Boolean integer programming [32] and staff assignment [30]. We use this method to solve Equation (7).

Adaptive Chaotic Discrete State Transition Algorithm
Discrete state transition algorithm is an individual-based optimization algorithm. Its basic idea is to regard the solution of an optimization problem as a state, and the process of updating the solution is called state transition. The standard form of the discrete state transition algorithm can be described as where x s ∈ Z is a current state; A s , B s are transformation operators; u s ∈ Z is a control variable; ⊕ is an operation; f (·) is the evaluation function, which is used to measure the quality of x s . The four special transformation operators [32] are as follows.
(1) Swap transformation: where A swap s ∈ R n×n is a random 0-1 matrix with swap action, called a swap transformation matrix, and m a is a swap factor that can control the number of swap elements in the solution. Swap transformation is called local exploration and global exploration when m a = 2 and m a ≥ 3, respectively.
(2) Shift transformation: where A shi f t s ∈ R n×n is a random 0-1 matrix with shift action, called a shift transformation matrix, and a shift factor, m b , can control the continuous shift of elements in the solution. If m b = 1, the shift transformation is regarded as local exploitation, and if m b ≥ 2, the shift transformation is regarded as global exploration.
(3) Symmetry transformation: where A sym s ∈ R n×n is a random 0-1 matrix with symmetry action, called a symmetry transformation matrix, and a symmetry factor, m c , can control the continuous symmetry of elements in the solution. Symmetry transformation is intrinsically called global exploration. (4) Substitute transformation: where A sub s , B sub s ∈ R n×n is a substitute transformation matrix, and m d is a constant integer, called a substitute factor, to control the maximum number of positions to be substituted. If m d = 1, the substitute transformation is regarded as local exploitation, and if m d ≥ 2, the substitute transformation is regarded as global exploration.
The initial solution is given randomly in the discrete state transition algorithm, and the solution significantly affects the convergence performance of the algorithm. DSTA easily falls into local optima in the iterative process. Therefore, an adaptive chaotic discrete state transition algorithm (ACDSTA) is proposed by introducing the opposition-based learning strategy, chaotic perturbation strategy, and adaptive recovery strategy.

Initialization Method Based on Opposition-Based Learning Strategy
The initialization method of the discrete state transition algorithm has the problem of uneven distribution, which somewhat affects its optimization efficiency. Therefore, an opposition-based learning strategy (OBL) is introduced to initialize the discrete state transition algorithm. OBL [33] is a machine-learning method whose idea is to generate a reverse solution based on the forward solution, compare their fitness values, and select the optimal solution as the initial solution, thereby improving the optimization speed of the algorithm. Let where L a and L b are the upper and lower bounds, respectively, of the value range of the target vector. Based on the idea of OBL, the best initial solution is selected as where f (·) is the fitness function, y o is the fitness value corresponding to X o , and X o is the initial solution into the subsequent iterative process.

Chaos Perturbation Strategy
Chaos comes from nonlinear dynamic systems. Because of its unique randomness, ergodicity, and complexity, it can effectively prevent the algorithm from falling into local optima, and it is widely used in optimization problems [34]. We propose the chaotic perturbation strategy to improve the problem that the discrete state transition algorithm can easily fall into local optima and is described as follows. During the algorithm iteration, when a value recurs a certain number of times, it indicates that the algorithm has fallen into a local optimum. Chaotic perturbation is applied to obtain a chaotic variable sequence, which is inversely mapped to the original solution space to obtain the perturbation solution, which is substituted in the next iteration so that the calculation exits the local extremum. When the termination condition is satisfied, the algorithm ends the iteration, and the final solution is the global optimum.
There are many rules to generate chaos, among which logistic mapping is common, but it has the problem of uneven frequency distribution. Zhou et al. [35] combined logistic mapping and tent mapping based on the uniform ergodicity of tent mapping to realize logistic-tent mapping, where X n ∈ [0, 1] is a chaotic variable, and the mod1 operation ensures that its output data are in the range of [0, 1], and α ∈ (0, 4] is a chaotic factor, which we take as 3.99. The chaotic variables generated by Equation (15) cannot be directly used for iterative calculation of the algorithm. So, chaotic variables are mapped to the solution space of the objective function, where X new is a new solution generated after chaotic perturbation.

Adaptive Recovery Strategy
In the iterative process of the discrete state transition algorithm, the chaos perturbation strategy is introduced to generate new solutions, which can effectively improve its ability to jump out of local optima. However, the new solution directly enters the next iteration, which will decrease the algorithm's convergence performance. To only use the greedy criterion can no longer meet the convergence requirements. We propose adaptive recovery, adopting the greedy criterion to ensure the general convergence of the algorithm. The current value is restored to the preserved historical best value with adaptive probability to further improve convergence performance.
The discrete state transition algorithm is in the stage of rapid optimization in the early stages of iteration, which greatly decreases fitness. The recovery probability should be small at this time, so as not to affect the searchability of the algorithm in the early stage of iteration. In the middle and later stages, the searchability decreases, and the optimal historical value should be restored with a large probability. We use a nonlinear adaptive adjustment method [36], where P ∈ [0, 1], P a is the maximum value of the recovery probability P, P b is the minimum value of P, iter is the current number of iterations, iter max is the maximum number of iterations, and µ is the adaptive factor, which we take as 2.
The steps of the proposed adaptive chaotic discrete state transition algorithm are as follows.
Step 1: Set relevant parameters such as swap factor m a , shift factor m b , symmetry factor m c , substitution factor m d , chaos factor α, and adaptive factor µ; Step 2: Generate the reverse solution using the initialization method according to Equation (13). The fitness of the forward and reverse solutions is compared by Equation (14) to select the best initial solution; Step 3: Update the solution by swap, shift, symmetry, and substitution transformations (Equations (9)-(12)) in turn; Step 4: According to Equation (17), judge whether the adaptive recovery strategy is satisfied, and if so, assign the optimal historical value to the current solution; Step 5: Determine whether the condition of the chaotic perturbation strategy is satisfied. If so, generate the chaotic variable according to Equation (15), and the perturbation solution X new according to Equation (16) for the next iteration. Otherwise, go to step 6.
Step 6: Determine whether the iteration termination condition is satisfied. If so, terminate the search and output the final optimization result. Otherwise, return to step 3.
The maximum numbers of iterations of EXP1-EXP3 were set to 100, 500, and 1000, respectively. The remaining parameters were set as follows: PSO learning probability c 1 = c 2 = 1.5, initial population 120, inertia weight 0.8. The ACDSTA parameters were set to SE = 30, m a = 2, m b = 1, m c = 0, m d = 1, and m max = 20. The parameter settings of DSTA and DSTAI were the same as in Zhou et al.'s study [32].
The average error, average value, and accuracy of optimization (S/20) were selected as performance evaluation indices of each algorithm. Table 1 compares the results of quadratic integer programming problems EXP1, EXP2, and EXP3, where the optimal values of each function index are in boldface.    Table 1 and Figure 1, we can see the following: (1) For the unconstrained integer optimization problem EXP1 with a low dimension, PSO, DSTAI, and ACDSTA can all find the optimal solution, and all the indices are the best; (2) For the unconstrained integer optimization problems EXP2 and EXP3 with higher dimensions, although PSO and DSTAI can also find the optimal solution, each index of ACDSTA is better. Among them, the optimization accuracies of EXP3 were improved by 65%, 75%, and 40% compared with PSO, DSTA, and DSTAI, respectively, and the respective average values were 20.3%, 133.1%, and 6.7% higher. Moreover, the average error was only 2.19%, which is significantly lower than for the other comparison algorithms. This indicates that the higher the dimensionality of the integer optimization problem and the larger the search space of the solution, the more prominent the advantage of ACDTSA; (3) Whether low-dimensional EXP1 or highdimensional EXP3, ACDSTA was superior to the other three algorithms in terms of initial average fitness and convergence speed and could find the optimal value faster. This shows that ACDSTA can approach the optimal value faster and improve convergence performance through the opposition-based learning strategy and adaptive recovery strategy; (4) For EXP3, the optimization accuracy of ACDSTA was better than that of the other three algorithms. Although there was a tendency to fall into local optima in late iterations, it could effectively jump out of local optima and obtain better optimization accuracy because of the chaos perturbation strategy.
In summary, ACDSTA was superior to the other algorithms in the test of unconstrained integer optimization problems, especially when the dimension and optimization range were extensive, which can better reflect the advantages of ACDSTA. This demonstrates the feasibility, superiority, and applicability of ACDSTA when solving practical engineering problems.

Rare Earth Extraction Process Analysis
Rare earth extraction is the obtaining of a single rare earth product from rare earth liquid, extractant, and detergent through specific equipment. We use the praseodymium/ neodymium (Pr/Nd) extraction and separation process as an example, as shown in Figure 2. The extractant is added from the first stage, and it flows from left to right through the stirrer. The detergent is added from the n + m stage and flows from right to left. The organic phase is added to the feed liquid from the nth stage. The solution in the tank is divided into two layers by stirring and clarifying. The upper layer is the organic phase, and the lower layer is the water phase. The difficult extraction product Pr is obtained in the lower layer of the first stage and the easy extraction product Nd from the upper layer of the n + m stage.
There are many extraction stages, and it is necessary to control the content distribution of Pr/Nd components in the tank to ensure the stability of product quality. It takes several hours or more to cause changes in the content of the export grade Pr/Nd component due to changes in the flow rate of the feed solution, extractant, and detergent. Therefore, it is necessary to obtain the residence time of materials in each unit group for timely control.

Time Delay Identification of Rare Earth Extraction Process
To verify the ability of the proposed method to solve the engineering delay problem, the delay identification of the 30-stage Pr/Nd production process in a rare earth extraction plant was carried out. During the Pr/Nd extraction process, the content of Pr/Nd components in different tanks changes with time, which leads to a change in the color of the solution. Therefore, the characteristic color variable of the solution image can be selected as an auxiliary variable to identify the time delay.
The sampling period was 5 min, and 220 groups of data of Nd component content and color characteristic variables (H, S, I, R, G, B) were selected in the continuous stable production process. The grey correlation method was used to analyze the correlation between Nd component content and color characteristic variables. The results are shown in Table 2, where the B component has the highest correlation degree, and the H component has the lowest correlation degree. Therefore, the B component is regarded as the critical process variable. In the actual extraction site, every five stages of the extraction tank share a set of agitators, i.e., every five stages constitute a unit group. Therefore, the 30-stage extraction tank was constructed as six groups of units for identification. According to the flow direction of the extractant, the first-stage inlet sampling data and each group of outlet sampling data were recorded as e 0 , . . ., e 6 , and the original data matrix E was obtained. According to the operation experience of the extraction separation site, the time delay range of stirring and clarification of each stage extraction tank is [3,8] min. Given the abovementioned construction method, the time delay range of each unit group is [15,40] min. Therefore, the value range of the time-base sequence is [3,8]. According to Equation (4), the time-correlation data matrix E is constructed, and the time delay sequence solution is transformed into the optimization problem by Equations (5) and (6), as shown in Equation (7).
ACDSTA is used to solve the abovementioned optimization problem. A certain number of time-based sequences is generated according to the range of time-based values, and the fitness function is constructed according to Equation (7). A new sequence is generated after the operation of the time-based sequence by the swap, shift, symmetry, and substitution operators. At the same time, the fitness value is calculated, and the current optimal value and the optimal time-based sequence are retained. The abovementioned operation is repeated until the iteration termination condition is satisfied, and the global optimal fitness value and global optimal time-based sequence are obtained. The solution process of ACDSTA is shown in Figure 3. Here, the maximum number of iterations is set to 100, SE = 5, and the remaining parameters are set according to Section 3.4.

Time Delay Identification Results and Method Verification
To verify the accuracy of the improved time delay identification method, we conducted the following experiments. Firstly, the characteristic components H, S, and I of the solution image of the rare earth extraction tank are considered auxiliary variables. A prediction model of Nd component content that meets the requirements of the extraction site is established by the wavelet neural network and used as a verification model. Then the maximum relative error (MAXRE), mean relative error (MEANRE), and mean absolute error (MAE) are determined to measure the performance of the model, where z is the predicted component content value of the wavelet neural network, and z is the actual component content value. Finally, the following comparative experiments are designed. Based on the data processed by the improved time delay identification method (Improved method), time-correlation analysis method (Original method), and Unused method, the performance of the component content model based on wavelet neural network (WNN) is compared and analyzed. WNN is a neural network based on wavelet analysis theory that combines the excellent time-frequency localization property of wavelet function and the powerful self-learning function of the neural network. WNN uses the wavelet basis function as the activation function, which has a strong prediction ability than the back propagation neural network. Lu et al. [38] showed a nonlinear relationship between the color characteristic components H, S, and I of the rare earth extraction solution image and Nd component content. Therefore, we use the WNN to model the component content of the rare earth extraction process, using the Morlet function as the wavelet basis function.
The parameter settings are as follows. The learning probability is 0.01, the momentum factor is 0.001, and the maximum number of iterations is 1000. The data samples consist of 220 groups, and 190 groups are randomly selected for model training. The remaining are used to verify the model, as shown in Figure 4. The maximum relative error of the model prediction is 4.76%, which meets the accuracy requirements of the rare earth extraction site for the prediction model, so it can be used as a verification model for time delay identification. Table 3 shows the model evaluation indices based on the Unused, Original, and Improved methods, where the bold data are the optimal values, and the corresponding relative error curve is shown in Figure 5. From Table 3 and Figure 5, we can see the following: (1) The performance of the component content model based on the Original method is better than that based on the Unused method. However, its maximum relative error is greater than 5%, which does not meet actual requirements. This shows that although the Original method can improve the model's performance to a certain extent, due to the randomness of its data selection and the lack of data preprocessing, the method fails to accurately identify the real-time delay.
(2) The mean relative error of the component content model based on the Improved method is better than that based on the Unused and Original methods, which is reduced by 70.1% and 63.9%, respectively. This shows that the model based on the improved method is more stable. (3) Compared with the Unused and Original methods, the mean absolute error of the model based on the improved method is reduced by 60.6% and 58.6%, respectively, indicating that the prediction accuracy of the model is higher.
In summary, the improved method significantly improved the model's performance. This shows that the improved method based on grey correlation analysis can effectively select the data closest to the real-time delay and improve the accuracy of the delay identification results. At the same time, the maximum relative error of the model based on the improved method is less than 5%, which meets the actual requirements. This shows that the improved time delay identification method proposed in this paper is suitable for the time delay identification of rare earth extraction process.

Conclusions
Rare earth extraction is a typical nonlinear, large-time-delay industrial process. The existence of time delay precludes the effective use of much field data and leads to a gap between the model describing the rare earth extraction process and the actual situation. We did the following work in response to this problem: Based on the standard discrete state transition algorithm, an improved algorithm (ACDSTA) was proposed, using an opposition-based learning strategy to initialize and accelerate the convergence of the algorithm, and an adaptive recovery strategy to improve its convergence performance. A chaotic perturbation strategy can improve the ability of the algorithm to jump out of local optima. An experimental comparison with three unconstrained integer optimization problems showed that ACDSTA can effectively solve such problems, and verified its effectiveness, superiority, and applicability; An improved time delay identification method was proposed to solve the problem that the data are not preprocessed and are randomly selected in the time-correlation analysis method. The method was applied to the time delay identification of a rare earth extraction process. The superiority and effectiveness of the proposed improved time delay identification method were verified by comparing the time-correlation analysis method and the data without the identification method under the same Nd component content. To sum up, the proposed time delay can provide a reference for modeling the rare earth extraction process and has guiding significance for the improvement of the online detection accuracy of component content.