Research on time series forecasting method of ion concentration in produced fluid

In the research on strong alkali‐surfactant‐polymer flooding scaling prediction, the variation characteristics of ion concentration in the produced fluid are obviously consistent with those of scale component content. Consequently, studying the variation tendency of ion concentration in produced fluid is helpful for further revealing the scaling tendency. Given the periodicity and chaos characteristics of the ion concentration data of the produced fluid, an echo state network (ESN) is used to realize the relevant time series forecasting. Simultaneously, to cope with the failure of the ESN in selecting suitable reservoir parameters according to different characteristics of time series, a modified discrete particle swarm optimization algorithm based on objective space decomposition is used to optimize the reservoir parameters. The experimental results indicate that the improved ESN presents the lowest error and is the closest to the target value.


INTRODUCTION
The scaling phenomenon caused by the injection of a strong alkali-surfactant-polymer (SASP) flooding system into the formation will lead to the blockage of some pores in the reservoir, as well as frequent failures of the lifting equipment of the production wells, which will seriously threaten the normal operation of crude oil extraction.Therefore, it is of paramount importance to predict the scaling tendency of reservoirs.Relevant studies have found that the variation characteristics of ion concentration in the produced fluid are obviously consistent with those of scale component content.With the continuous injection of SASP and the periodic addition of anti-scaling chemicals, the water ion concentration in the oilfield presents short-term predictability and long-term unpredictability, that is, chaos characteristics.Against this background, the problem to be solved in this research is to predict the tendency of time series data with chaos characteristics and periodicity. 1n the time series forecasting model, given the incompleteness of data, we selected the adaptive prediction model to conduct research. 2The neural network is regarded as an effective tool for the approximation and modeling of nonlinear systems. 3,4Feedforward neural network (FFNN), such as multilayer perceptron network (MLP) [5][6][7] and radial basis function (RBF), 8,9 is the earliest neural network model applied to nonlinear time series forecasting, achieving preferable application effects.Subsequently, the support vector machine (SVM) and kernel method are progressively applied to pattern classification and regression problems. 10It is worth acknowledging that SVM embodies a host of incomparable advantages over FFNN and therefore is widely utilized in the forecasting of nonlinear time series.Nonetheless, in the case of establishing a nonlinear mapping relationship, MLP, RBF, and SVM, without exception, are all processing methods for static feedforward, essentially simplifying the dynamic modeling problem to a problem related to the approximation of static functions.Therefore, the foregoing methods inevitably expose a certain degree of limitation.
The emergence of the recurrent neural network provides an effective tool for solving the time series problem. 11Conceptually, the recurrent neural network is defined as a dynamic neural network with a delayed feedback link, in which the feedback connection enables the output of the network to be affected not only by the current input but also by the historical input, thus effectively depicting the dynamic characteristics.In the recurrent network, however, the learning algorithm of the network has consistently been a difficulty in the current research.
The reservoir network is a new type of recurrent neural network, which, to a certain extent, solves the deficiencies related to the difficulty in training the traditional recurrent neural network.The reservoir network involves three specific implementations, namely, the echo state network (ESN) proposed by Jaeger in 2001, the liquid state machine (LSM) proposed by Mass and Natschlaeger et al. in 2002, and backpropagation decorrelation (BPDC) proposed by Steil 12 in 2004 13,14 .Specifically, LSM focuses on the simulation of biological neural networks; ESN with linear or sigmoid-type analog neurons focuses on engineering applications 12 ; BPDC uses a relatively complex learning algorithm with a single-step output to the state feedback. 15eep learning models have been widely applied in the fields of healthcare, image processing, machine translation, 15,16 and autonomous driving, 17 and have achieved excellent performance.However, a deep neural networks is a complex structure with high nonlinearity, which can be regarded as a "black box" model in essence.Obviously, this opaque modeling technology cannot be trusted in many fields, which severely limits the use of deep learning models in many sensitive or high-risk fields.
In view of the periodic and chaotic characteristics of ion concentration data of produced liquid, ESN has the most advantages.Unlike traditional neural networks, ESN will not collapse due to input data containing noise.Instead, it can adapt to noise and learn from it.ESN can handle dynamic data because it has recursiveness, and it can capture and utilize time correlation in time series data, unlike other static neural networks.Most important of all, it can use known information to predict future states, thereby reducing the amount of computation needed.
Therefore, we use ESN for time series prediction.To solve the problem that echo state network cannot select suitable parameters of reserve pool according to different characteristics of time series, particle swarm optimization algorithm based on object space decomposition is used to optimize the parameters of reserve pool.The specific implementation plan is as follows.

Principle of ESN
The ESN structure mainly includes an input layer, a reservoir, and an output layer, as shown in Figure 1, where the input layer and the reservoir are connected by the input weight matrix, while the internal units of the reservoir are connected by the internal weight matrix to form state feedback.The state reservoir and the output layer are connected by the output weight vector.
ESN introduces a reservoir calculation model to replace the original hidden layer.In particular, the connection state of neurons in the reservoir is random without manual control.Meanwhile, the connection weight as a fixed value not only further reduces the calculation time of training in the application process of the reservoir but also avoids the generation of the local minimum solution of the gradient descent optimization algorithm to a certain extent, thereby endowing the network structure with dynamics, adaptability, and self-organization. 1 Assuming that the obtained discrete time series is s(1), s(2), … , s(t), … , in which the observed value s(t) = [s 1 (t), s 2 (t), … , s l (t)] T at time t is an l-tuple vector.The time series is modeled on the basis that the input at time t is u in (t) = s(t) and the output is y(t) = s(t + ), and then the prediction of different time domains is realized by changing the value of .Additionally, a multi-input and single-output prediction model with y(t) as the scalar is established, with the first-dimension observation value s 1 (t + ) of the sequence s(t + ) being selected.
The state equation of ESN is as follows: where (⋅) is the activation function of the neuron, which usually selects the Sigmod or tanh function, and W in and W x are the input connection matrix and the reservoir internal connection matrix, respectively, which remain unchanged after initialization.In addition, T is the internal state vector of the N-tuple reservoir at time t, which is generated under the excitation of the current input u in (t) and the last-time state x(t − 1).The output equation can thus be obtained as follows: where W out denotes the output weight vector of the reservoir and is also the only parameter to be solved by training.Generally, the spectral radius of W x should satisfy the condition of (w x ) < 1 to ensure the stable operation of the ESN network.
The ESN uses the least squares method to solve the output weights, minimizing the objective function as shown in the following formula: where T , where p represents the length of the training data, and

Setting of reservoir parameters
The performance of the reservoir depends on the scale N of the reservoir, the spectral radius R, the sparsity D, and the input transformation factor S. Admittedly, the key to constructing ESN is to select the optimal solution of the above four parameters regarding their combination.As a decisive parameter affecting the prediction performance of ESN, the scale N of the reservoir refers to the number of neurons in the reservoir.An excessively large value of N tends to cause an overfit, whereas an excessively small value of N tends to cause an underfit.
The spectral radius R refers to the absolute value of the maximum eigenvalue of the internal connection matrix W x of the reservoir, and the value range is usually between [0, 1].Moreover, different application scenarios of ESN will lead to different values of R.
The sparsity D represents the degree of sparseness in connection with the internal neurons of the reservoir, and the neurons are partially connected.Generally, most of the elements in the connection matrix W x of the reservoir are set to zero.The proponent of ESN should set the reservoir's power to be satisfied under the condition of D ∈ [0, 0.1].
The input transformation factor S indicates a coefficient scale that is scaled before the signal is input to the reservoir, which characterizes the range of values for the input connection weights.According to Formula (1), S determines not only the working range of the activation function but also the effect intensity of the input on the reservoir state variable, with a value range of [0, 1].
Given that the ESN varies according to the characteristics of the training data, it is necessary to design a modified discrete particle swarm optimization (MDPSO) algorithm based on objective space decomposition to obtain the optimal solution of the ESN parameters, to improve the prediction performance of the ESN.

Analysis of Multi-objective Optimization Problem
Multi-objective optimization problem 18,19 (MOP) aims to obtain the optimal solution for multiple objectives in a given region, and its mathematical expression is as follows: where ) represents all objective functions; g i (x) ≤ 0 represents the i-th inequality constraint; h j (x) = 0 represents the j-th equality constraint; Y = {F(x)|x ∈ Ω} represents the objective space.The particle swarm optimization algorithm is one of the methods to solve MOP.To put it concretely, the algorithm is initialized as a group of random particles and iterates to find the optimal solution.The velocity of the k-th particle is expressed as v k = (v k1 , v k2 , … , v kn ), with its position expressed as x k = (x k1 , x k2 , … , x kn ).In addition, the current optimal position (pbest) of the particle and the optimal neighborhood particle position (gbest) are used to update the velocity vector of the particle.Further, the update of the velocity and position of the k-th particle is as follows: where j = 1, 2, … , n, t is algebra; c 1 and c 2 are learning factors, provided that c 1 = c 2 = 2; r 1 and r 2 are random numbers between [0, 1]; w is the inertia weight.In addition, the particle is limited by the maximum speed v max when adjusting its position according to the speed, and will be limited to v max when v t kj exceeds v max .In the interest of maintaining the diversity and improving the rate of convergence of the algorithm, the particle swarm optimization algorithm based on objective decomposition is used to solve the MOP.Moreover, space decomposition, solution classification, and update strategy are used to maintain the diversity of the solutions of the objective function.Meanwhile, interlace operation and particle swarm update strategy are used to realize decision space search, with selection strategy utilized to realize local or global search.

Objective space decomposition and solution classification
The objective space Y of the MOP is decomposed into sub-objective spaces based on a set of direction vectors, including Y 1 , Y 2 , … , Y N , and the resulting solutions are classified in the sub-objective spaces such that each sub-objective space corresponds to at least one solution.Therefore, the current solution set POP of the vector set direction vectors is further obtained.The objective space decomposition and solution classification can be expressed by the following formulas: where is the cosine of the angle between the direction vectors  i and F(x) − Z.In case p i (1 ≤ i ≤ N) is empty, the largest solution x i of Δ ( F(x),  j ) in POP is selected to be put into p i , so that each sub-objective space corresponds to at least one solution, thus maintaining the diversity of the obtained solutions.
By Formula (7), the solution set POP is divided into N classes, namely, p 1 , p 2 , … , p N .Likewise, the objective space Y is divided into N sub-objective spaces, namely, Y 1 , Y 2 , … , Y N , by using Formula (8).

Update strategy of population classification
The elitism preservation strategy is typically beneficial to speed up the rate of convergence and improve the quality of solutions.Therefore, an update strategy based on the elitism preservation strategy is used to maintain the diversity of solutions and improve the quality of solutions.Details of the update strategy are as follows: If the objective vector of the current solution p i is not in the sub-objective space Y i , and the objective vector of the new solution is in Y i or the new solution dominates the current solution, the new solution is the current solution of Y i .
If p i is in Y i , this solution is selected when there is only one solution in Y i ; when there are two or more solutions in Y i , the objective vector of the new solution must be in Y i , and the new solution dominates p i .If the new solution and p i do not dominate each other, the angle between the objective vector of the new solution and the direction vector  i is smaller than the angle between the objective vector of p i and  i .
Interlace operation and neighborhood correction are employed to avoid strategies that may result in some better solutions being discarded.The interlace operation of PSO can be used to achieve a global search, while the neighborhood correction can be used to perform a local search.In the interlace operation, the interlace operation of the PSO uses Formulas ( 5) and (6).In contrast, neighborhood correction performs a search through the neighborhood of the current particle, which improves the quality of the solution and accelerates its convergence to a good Pareto front.Also, it can be used to improve the ability of local search.In this work, a particle x generates a new solution according to the following correction operation: where x ′ is the best neighborhood particle position of x, and x ′′ is the current optimal position of x.
According to the method, the MDPSO algorithm is designed as follows: Step 1: initialization.Given N direction vectors (  1 ,  2 , … ,  N ) , t = 0, the initial population POP(t) and initial velocity set V(t) are randomly generated, the size of which is N, and Step 2: generate a new solution.Formulas ( 5), (6), and ( 9) are used to produce generations of the temporary population pop, and then a new speed set V(t + 1) is obtained, with the generations set marked as o.
Step 3: update solutions.For each Specifically, first, the population o is classified by Formula (7).Furthermore, the update strategy is used to select N optimal solutions to be put into population POP(t + 1) and update t = t + 1.
Step 4: terminate the algorithm.If the termination condition is satisfied, the algorithm is terminated; otherwise, Step 2 should continue to be executed.
In this algorithm, the pbest and gbest of each particle do not need to be stored, thus reducing the storage space.Meanwhile, the algorithm ensures that each sub-region corresponds to a solution in each generation and improves the diversity of the population during the search process.

Optimization of ESN reservoir parameters based on MDPSO
The MDPSO algorithm is employed to perform a joint solution on {N, R, D, S}, thereby transforming the optimization problem for the output weights into an optimization problem for multiple objectives, with root mean square error (RMSE), symmetric mean absolute percentage error (SMAPE), and normalized root mean squared error (NRMSE) as optimization functions. 20,21Therefore, the optimization objectives of MDPSO are as follows: where y(t) is the true value of the test data; y ′ (t) is the network prediction values; y(t) is the average of the true values; and n is the test set capacity.The steps of the ESN time series forecasting algorithm based on MDPSO are described as follows: Step 1: Initialize the iterations of PSO.Given N direction vectors (  1 ,  2 , … ,  N ) , t = 0, an initial population POP(t) and an initial velocity set V(t) are randomly generated.The size of which is N, and The initial value of the reservoir parameter {N, R, D, S} of ESN should be set.
Step 2: Establish the ESN model and process the training data of ESN i by using the normalization formula.The objective optimization function shown in Formula ( 10) is used to calculate the fitness values of all particles in the population to obtain pbest and gbest.
Step 4: Repeat Steps 2 and 3 until the iterations are met, thus acquiring the optimal solution that meets the objective function.
Step 5: Output the global optimal solution as the optimal solution of ESN optimization parameters and use ESN I to conduct prediction to obtain prediction results.

Performance comparison experiment of MDPSO algorithm
We choose non-dominated sorting genetic algorithm-II (NSGA-II), wild geese algorithm (WGA), and white shark optimizer (WSO) to compare with MDPSO.The test functions used in the experiment include four ZDT problems 22 and tow DTLZ problems 23 to verify the convergence and diversity of the algorithm.Experimental parameters are set as follows: 1.The population size is 100 individuals.2. The crossover probability is 0.8; The probability of variation is 0.3.The number of iterations is 100.
Table 1 shows the generational distance (GD), inverted generational distance (IGD), and hypervolume of the solutions obtained by 30 independent runs of the four algorithms (HV) statistical results of average and optimal values of indicators.
In Table 1, for metric GD, MDPSO is inferior to WGA and NSGA-II in ZDT1 problem, but superior to WSO; in ZDT3 problem, MDPSO is inferior to NSGA-II, superior to the other two algorithms; for measuring IGD, MDPSO is only disadvantageous on problem ZDT1; for measuring HV, MDPSO is inferior to WGA for ZDT1 and ZDT3, but superior to the other two algorithms.These results show that MDPSO is better than NSGA-II, WSO and WGA in maintaining the convergence and diversity of solutions among the six test functions.

Comparison experiment of prediction performance of MDPSO-ESN model
To verify the validity of the model proposed in this research, the pH value in the ion assay data is selected as the experimental data.The data volume is 500, 70% of the data are experimental data, whereas the remaining 30% are target data (i.e., test data).Due to the large amount of data, Figure 2 takes one data from five points for display, showing a total of 100 points.Meanwhile, echo state network based on improved differential evolution (IDE-ESN), 24 PSO-ESN, and ESN are adopted as comparative models, with an ESN spectral radius of 0.98, a sparsity of 0.02, an input transformation factor of 0.5, and a reservoir scale of 200.The mean square error (MSE) is used as an evaluation index 24 to verify the performance of different algorithms.The specific formula is as follows: where T is the number of samples; y d (t) represents the forecast value; and y(t) represents the prediction value.Figure 3 compares the prediction results of the four models and the target data.As can be seen from Figure 3, the prediction values and tendencies of the model proposed in this article and the IDE-ESN model are closer to the target data.From the comparison of running time shown in Table 2, it can be seen that the running time of the  ESN model without the optimization algorithm is the shortest.In contrast, the running time of the MPSO/D-ESN algorithm proposed in this article is slightly higher than that of ESN and PSO-ESN prediction models.Nevertheless, from the results of the MSE evaluation index and the prediction errors shown in Figure 4, it can be seen that the prediction performance of MPSO/D-ESN proposed in this article is better than that of the other three comparative models.More specifically, the performance of IDE-ESN is better than PSO-ESN, followed by ESN.It is proved that the prediction effect of the model proposed in this article is better than that of the other three comparative models.

CONCLUSION
In order to realize the prediction of periodic and chaotic time series data, a particle swarm optimization algorithm based on object space decomposition is proposed in this article to optimize the parameters of the reserve pool.In order to solve the problems of premature convergence of the particle swarm optimization algorithm and its tendency to fall into local optimality, the target space decomposition and solution classification steps are designed to improve the convergence speed of the algorithm.Finally, the prediction performance of echo state network is improved.The prediction results show that MDPSO-ESN has the lowest error and is closer to the real data than IDE-ESN and PSO-ESN, which verifies the effectiveness of the improved method.

F I G U R E 2
Experimental data of pH valueF I G U R E 3 Prediction results of four models TA B L E 1 GD, IGD, and HV measurements of the four algorithms Note:The bolded values represent the optimal values for the experiment.By bolding the optimal values, we can intuitively see which algorithm works best.
TA B L E 2 MSE prediction and running schedule of four models F I G U R E 4 Prediction errors of four models