Attention mechanism-based multi-column neural network effluent total phosphorus concentration multi-step prediction method

Abstract

The invention discloses an attention mechanism-based multi-column neural network effluent total phosphorus concentration multi-step prediction method, relates to the field of artificial intelligence, and is directly applied to the field of sewage treatment. The soft measurement method is used for multi-step prediction of the effluent total phosphorus concentration, and the problems that in the sewage treatment process, the effluent total phosphorus concentration trend is difficult to master, and the prediction cost is high are solved. Aiming at the problem that potential importance difference between a target variable and each auxiliary variable in a prediction task is neglected in traditional multi-step prediction, an input reconstruction mechanism based on an attention mechanism is provided, and the generalization performance of a prediction model is improved; meanwhile, for the problem that only a single error evaluation index is considered in traditional multi-step prediction task division, a loss function combining prediction errors and shape errors is provided for model training, the prediction precision of the model is improved, and management operation of the sewage treatment plant is promoted.

Description

Technical Field

[0001] The present invention is a multi-step prediction method for effluent total phosphorus concentration based on a multi-column neural network with an attention mechanism, which aims to achieve accurate multi-step prediction of effluent total phosphorus concentration in the sewage treatment process. It involves the field of artificial intelligence and is directly applicable to the field of sewage treatment. Background Art

[0002] Total phosphorus is a key water quality indicator. Excessive total phosphorus concentration may lead to eutrophication of water bodies, which in turn may cause excessive growth of algae and affect the health of the ecosystem. In the sewage treatment process, whether the total phosphorus concentration of effluent meets the standard is one of the important indicators for evaluating the effluent quality of sewage treatment plants. This is not only related to the self-repair ability of the water body, but also affects the ecological balance and biodiversity of the surrounding water bodies.

[0003] The actual situation in the sewage treatment process often has complex dynamic changes, and single-step prediction cannot fully consider long-term trends and potential fluctuations. Multi-step prediction can identify the changing trend of the total phosphorus concentration in the effluent in advance and provide managers with early warning information for a longer period of time. This method is of great significance because it allows preventive measures to be taken when facing potential problems in the sewage treatment process, rather than making adjustments after the problem occurs. To solve this problem, the present invention proposes a multi-step prediction method for the total phosphorus concentration in the effluent based on a multi-column neural network with an attention mechanism, which aims to achieve accurate multi-step prediction of the total phosphorus concentration in the effluent during sewage treatment. Summary of the invention

[0004] 1. Technical problems that the present invention needs and can solve:

[0005] This paper proposes a multi-step prediction method for effluent total phosphorus concentration based on a multi-column neural network with an attention mechanism, which realizes accurate multi-step prediction of effluent total phosphorus concentration based on data collected during the actual sewage treatment process. This method effectively solves the problem of difficulty in grasping the trend of total phosphorus concentration changes during sewage treatment, and significantly improves the flexibility of control decisions in urban sewage treatment plants.

[0006] 2. Specific technical solutions of the present invention:

[0007] The present invention proposes a multi-step prediction method for effluent total phosphorus concentration based on a multi-column neural network with an attention mechanism. The algorithm comprises the following steps:

[0008] Step 1: Initialization; set the maximum number of iterations I max and population size Q, prediction step number M, and set the initial iteration number i = 1;

[0009] Step 2: Reconstruct the input sample based on the attention mechanism;

[0010] Collect actual water quality parameter data of the sewage treatment plant, predict variable (total phosphorus in effluent) Y, auxiliary variable X=[X1,…,X n ,…,X N ], where X n represents the nth auxiliary variable, and N represents the number of auxiliary variables;

[0011] Step 2.1: Data normalization;

[0012] Normalize the obtained auxiliary variables X and predictor variables Y to [0,1]:

[0013]

[0014] in, and Respectively represent the normalized nth auxiliary variable and output variable; min(X n ) and max(X n ) represent the minimum and maximum values ​​of the nth auxiliary variable, min(Y) and max(Y) represent the minimum and maximum values ​​of the predictor variable, respectively;

[0015] Step 2.2: Calculate the attention score for each auxiliary variable; for the nth auxiliary variable Calculate its effect on the predictor variable The attention score ξ n for:

[0016]

[0017] Step 2.3: Based on the attention score, assign attention weights to each auxiliary variable; use the softmax function to calculate according to formula (4) The attention weight α n for:

[0018] α n =softmax(ξ n ) (4)

[0019] Step 2.4: Reconstruct the input sample; the network input Z is determined by the sum of the weighted auxiliary variables and the predictor variables, and the calculation formula is:

[0020]

[0021] Among them, λ is the weight coefficient, which is used to adjust the influence of the prediction variable and the auxiliary variable, and is set in the range of [0,0.4] according to experience; therefore, for the multi-step prediction of M steps, assuming that the input variable Z at the current moment and the previous D-1 time points is Z = [z t-D+1 ,…,z t-D+d ,…,z t ]T The predictor variable Y used to predict the next M steps is [y t+1 ,…,y t+m ,…,y t+M ] T , then the input dimension is D, and the corresponding output dimension is M;

[0022] Step 3: Task coding;

[0023] For the M-step prediction task, binary particle encoding is used, and for each particle, its position is represented as B = [b1,…,b m ,…,b M ], where each element represents a prediction step number, with a value of 0 or 1; at the initial moment, the positions of all particles are randomly generated in binary form, and their corresponding velocities V are randomly initialized in the range of [-4,4];

[0024] Step 4: Task decoding;

[0025] For a given position B, consecutive elements with the same value are considered to belong to the same subtask; assuming that B generates C subtasks, a vector representing the subtask allocation will be obtained, denoted by S = [s1,...,s c ,...,s C ](c=1,2,...,C), where s c represents the prediction step size of the cth subtask;

[0026] Step 5: Design a multi-column neural network model;

[0027] After task decoding, each subtask is predicted by a sub-network, thus forming a multi-column neural network consisting of C sub-networks, and the output dimensions of each sub-network are s1,...,s c ,...,s C ; Taking the cth sub-network as an example, its network structure is designed as follows;

[0028] Step 5.1: Initialize network parameters;

[0029] The network does not initially set up cross-layer direct connection between the input layer and the output layer, that is, the connection weight matrix V between the input neuron and the output neuron is 0, and the number of neurons in the hidden layer of the network is H = 0; the maximum number of neurons in the hidden layer is set to H max , the maximum number of iterations S max ;

[0030] Step 5.2: Determine the network input and output;

[0031] Taking the cth subnetwork as an example, the input dimension D c Determined by the following formula:

[0032]

[0033] Among them, D min and D max D c The boundary value of the network, the output dimension is s c , from m c-1 +1 to m c ,in For simplicity, use represents the network input, Represents network output;

[0034] Step 5.3: Build the input and output direct connection structure;

[0035] For all input neurons and output neurons, a direct connection structure is constructed, and the connection weight is expressed as Among them, v d represents the connection weights between d input neurons and all output neurons, and is set to 1 accordingly;

[0036] Step 5.4: Network growth process;

[0037] Hidden layer neurons are added in pairs to compensate for the maximum error, represented by neurons H+1 and H+2, respectively, and their centers, widths, and weights are set as:

[0038]

[0039] σ H+1 =σ H+1 =1 (8)

[0040] w H+1 =1, w H+2 =-1 (9)

[0041] in, is the input sample with the largest network error, k max The calculation formula is:

[0042]

[0043] Where E=[E1,…,E k ,…,E K ] T is the network error vector calculated for all samples, and the error of the kth sample is obtained by the following formula:

[0044]

[0045] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, o s(k) and are the expected output and the predicted output respectively; therefore, the number of hidden neurons H increases by 2;

[0046] Step 5.5: Online learning;

[0047] The network parameters are updated using the improved Levenberg-Marquardt (LM) algorithm; let Δ contain all the parameters of the network to be trained (including the center C and width σ of the neuron, the connection weight W from the hidden layer to the output layer, and the connection weight V from the input layer to the output layer):

[0048]

[0049] Compute the Jacobian matrix row vector j s (k):

[0050]

[0051] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, c h is the center of the hth hidden layer neuron, σ h is the width of the hth hidden layer neuron, w h,s is the weight from the hth hidden layer neuron to the sth output neuron, v d,s is the weight from the dth input neuron to the sth output neuron.

[0052] The Jacobian vector components can be obtained according to the chain rule as follows:

[0053]

[0054] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, c h is the center of the hth hidden layer neuron, σ h is the width of the hth hidden layer neuron, w h,s is the weight from the hth hidden layer neuron to the sth output neuron, v d,s is the weight from the dth input neuron to the sth output neuron,

[0055] is the kth input sample, is the d-th dimension input in the k-th input sample, φ h is the radial basis function.

[0056] Calculate the Hessian matrix Q. In order to reduce the space complexity and time complexity of the Jacobi matrix calculation process, the Hessian matrix is ​​converted into multiple Hessian sub-matrices q s(k) and:

[0057]

[0058] Among them, j s (k) is the row vector of the Jacobian matrix. Similarly, the calculation of the gradient vector g is also converted into K gradient sub-vectors η s (k) and:

[0059]

[0060] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, j s (k) is the Jacobian matrix row vector;

[0061] The network structure parameter update rules are as follows:

[0062] Δ i+1 =Δ i -(Q i +μ i I) -1 g i (twenty three)

[0063] Among them, Δ i+1 and Δ i are the network structure parameters for the i+1th and ith trainings, respectively, Q i is a Hessian-like matrix, I is the identity matrix, μ i is the learning rate parameter for the i-th training, which is adaptively adjusted according to the change of training error (when RMSE increases, μ i Zooming in ten times, when RMSE decreases, μ i Reduced tenfold), g i is the gradient vector; the training process is repeated iteratively until i reaches the maximum number of iterations S max (Set according to experience, generally set to 100); calculate the root mean square error RMSE at the end of training. If RMSE converges (less than the set threshold 0.01) or the network hidden layer presets the maximum structure H max (according to experience, it is usually set to 10), then the training stops; otherwise, it returns to the growth type construction of the neural network;

[0064] Step 5.6: Calculate multiple columns of neural network output;

[0065] After the C sub-networks are constructed, the sub-network outputs are integrated to obtain the M-step output of the multi-column neural network.

[0066] Step 6: Calculate the fitness function;

[0067] Calculate the fitness function of all particles; the formula is as follows:

[0068] f=L RMSE +L DTW (twenty four)

[0069] Among them, L RMSE is the individual loss, L DTW is the trend loss; L RMSE The average RMSE of all forecast horizons is used to evaluate the forecast error, and the calculation formula is:

[0070]

[0071] Among them, y t+m (k) and are the expected output and predicted output of the mth step of the kth sample, respectively, K is the number of training samples; L DTW The shape difference between the expected time series and the predicted time series is measured by the dynamic time warping (DTW) algorithm, and the calculation formula is:

[0072]

[0073] Among them, m1=1,2,...,M, m2=1,2,...,M are the step positions of expected output and predicted output, Predict the output of step m1 for the k-th sample, is the expected output of the m2th step of the kth sample, is the cumulative distance of the optimal path from (1,1) to (M,M), and the goal of DTW is to obtain the minimum Euclidean distance between two sequences. value;

[0074] After calculating the fitness values ​​of all particles, the optimal position P of the particle can be determined. best and the optimal position of the population G best ;

[0075] Step 7: Particle parameter update;

[0076] The update rules for the velocity V of all particles are as follows:

[0077] V(i+1)=wV(i)+C1R1(P best (i)+B(i))+C2R2(G best (i)+B(i))(27)

[0078] Among them, R1 and R2 are two random numbers uniformly distributed in the range [0,1], w is the inertia weight, C1 and C2 are two acceleration coefficients, P best is the optimal position of the particle, G bestis the optimal position of the population, B i is the position vector of the current particle; by using the sigmoid transformation, the speed is limited to the interval [0,1]:

[0079]

[0080] Among them, v m Represents the mth element in vector V (m=1,2,…,M), and the calculation formula is:

[0081]

[0082] Here, v max is a constant value, set to 4; then, the position of the mth component of the particle is:

[0083]

[0084] Where ζ is a randomly generated real number between 0 and 1;

[0085] After the positions and velocities of all particles are updated, i=i+1;

[0086] Step 8: Determine the termination condition of multi-step prediction;

[0087] If i = I max Or the fitness function f converges (less than the set threshold 0.01), then the optimal task decomposition G for multi-step prediction is output best , and the M-step predicted value; if not, return to step 4 until the termination condition is reached;

[0088] Step 9: Prediction of total phosphorus in effluent;

[0089] The test sample data is used as the input of the trained multi-column neural network based on the attention mechanism to obtain the output of the network, which is then denormalized to obtain the predicted value of the effluent total phosphorus concentration.

[0090] Compared with the prior art, the present invention has the following obvious advantages and beneficial effects:

[0091] (1) The present invention uses a soft measurement method to perform multi-step prediction of the total phosphorus concentration in the effluent, which solves the problems of complex measurement operation and difficulty in grasping the change trend of the total phosphorus concentration in the effluent during sewage treatment. A multi-column neural network multi-step prediction method for the total phosphorus concentration in the effluent based on the attention mechanism is proposed, which has the advantages of simple model, high measurement accuracy, and strong generalization ability.

[0092] (2) Aiming at the problem that the traditional multi-step prediction ignores the potential importance difference between the target variable and each auxiliary variable in the prediction task, this paper proposes an input reconstruction mechanism based on the attention mechanism. By assigning weighted inputs to different auxiliary variables, the generalization performance of the prediction model is improved.

[0093] (3) In order to address the problem of only considering a single error evaluation index in the traditional multi-step prediction task division, the present invention proposes a loss function that combines the prediction error and the shape error for model training. At the same time, it considers the individual differences and trend differences of the multi-step prediction results, thereby improving the prediction accuracy of the model. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a model structure diagram of the present invention;

[0095] Figure 2 This is a multi-step prediction effect diagram of the effluent total phosphorus concentration of the present invention; DETAILED DESCRIPTION

[0096] The present invention obtains a multi-step prediction method for effluent total phosphorus concentration based on a multi-column neural network with an attention mechanism, and realizes accurate multi-step prediction of effluent total phosphorus concentration according to data collected during the actual sewage treatment process; this method effectively solves the problem that it is difficult to grasp the trend of effluent total phosphorus concentration changes during the sewage treatment process, and improves the flexibility of sewage treatment plant decision-making;

[0097] The experimental data comes from 2000 sets of samples from a sewage treatment plant in Beijing, China from April to November 2021. The collected variables include (1) influent flow, (2) effluent flow, (3) return flow, (4) regulating tank water level, (5) influent pH, (6) influent chemical oxygen demand (COD), (7) effluent COD, (8) effluent total nitrogen, and (9) effluent total phosphorus, among which effluent total phosphorus is the predictor variable and the remaining variables are auxiliary variables; 1000 sets of data are randomly selected as training samples, and the remaining 1000 sets of data are used as test samples;

[0098] Step 1: Initialization; set the maximum number of iterations I max And the population size Q, the number of prediction steps M, and set the initial number of iterations i = 1; I max Usually the value is 50-100, here we take 50; Q usually takes 20-40, here we take 20; M takes 18;

[0099] Step 2: Reconstruct the input sample based on the attention mechanism;

[0100] Collect actual water quality parameter data of sewage treatment plants, predict variable Y, auxiliary variable X=[X1,…,X n ,…,X N ], where X nrepresents the nth auxiliary variable, and N represents the number of auxiliary variables. In this embodiment, the number of auxiliary variables N is 8, including (1) inlet flow rate, (2) outlet flow rate, (3) return flow rate, (4) regulating tank water level, (5) inlet pH, (6) inlet chemical oxygen demand (COD), (7) outlet COD, and (8) outlet total nitrogen. The predicted variable Y is the outlet total phosphorus.

[0101] Step 2.1: Data normalization;

[0102] Normalize the obtained auxiliary variables X and predictor variables Y to [0,1]:

[0103]

[0104] in, and Respectively represent the normalized nth auxiliary variable and output variable; min(X n ) and max(X n ) represent the minimum and maximum values ​​of the nth auxiliary variable, min(Y) and max(Y) represent the minimum and maximum values ​​of the predictor variable, respectively;

[0105] Step 2.2: Calculate the attention score for each auxiliary variable; for the nth auxiliary variable Calculate its effect on the predictor variable The attention score ξ n for:

[0106]

[0107] Step 2.3: Based on the attention score, assign attention weights to each auxiliary variable; use the softmax function to calculate according to formula (4) The attention weight α n for:

[0108]

[0109] Step 2.4: Reconstruct the input sample; the network input Z is determined by the sum of the weighted auxiliary variables and the predictor variables, and the calculation formula is:

[0110]

[0111] Among them, λ is the weight coefficient, which is used to adjust the influence of the prediction variable and the auxiliary variable, and is set in the range of [0,0.4] according to experience; therefore, for the multi-step prediction of M steps, assuming that the input variables Z at the current time and the previous D time points are Z = [z t-D+1 ,…,z t-D+d ,…,z t ] TPredictor variables used to predict the next M steps Then the input dimension is D, and the corresponding output dimension is M. In this embodiment, λ is 0.2, M is the prediction step size, which is 18 here, and the maximum input dimension of D is 20;

[0112] Step 3: Task coding;

[0113] For the M-step prediction task, binary particle encoding is used, and for each particle, its position is represented as B = [b1,…,b m ,…,b M ], where each element represents a prediction step number, with a value of 0 or 1; at the initial moment, the positions of all particles are randomly generated in binary form, and their corresponding velocities V are randomly initialized in the range of [-4,4];

[0114] Step 4: Task decoding;

[0115] For a given position B, consecutive elements with the same value are considered to belong to the same subtask; assuming that B generates C subtasks, a vector representing the subtask allocation will be obtained, denoted by S = [s1,...,s c ,...,s C ](c=1,2,...,C), where s c represents the prediction step of the cth subtask; assuming a particle is B = [10 0 1 1 0 1 1 10 0 1 0 1 1 0 0 0], the prediction task M = 18 is divided into C = 10 subtasks, and the task set is S = [1, 2, 2, 1, 3, 2, 1, 1, 2, 3];

[0116] Step 5: Design a multi-column neural network model;

[0117] After task decoding, each subtask is predicted by a sub-network, thus forming a multi-column neural network consisting of C sub-networks, and the output dimensions of each sub-network are s1,...,s c ,...,s C ; Taking the cth sub-network as an example, its network structure is designed as follows;

[0118] Step 5.1: Initialize network parameters;

[0119] The network does not initially set up cross-layer direct connection between the input layer and the output layer, that is, the connection weight matrix V between the input neuron and the output neuron is 0, and the number of neurons in the hidden layer of the network is H = 0; the maximum number of neurons in the hidden layer is set to H max , the maximum number of iterations S max In this embodiment, H max Take 10, the maximum number of iterations S max Take 100;

[0120] Step 5.2: Determine the network input and output;

[0121] Taking the cth subnetwork as an example, the input dimension D c Determined by the following formula:

[0122]

[0123] Among them, D min and D max D c The input dimension of the sub-network will increase as the prediction range expands, thereby reducing the waste of input computing resources for short-term predictions and incorporating more historical information for long-term predictions; the output dimension of the network is s c , from m c-1 +1 to m c ,in For simplicity, use represents the network input, Represents the network output; in this embodiment, D min Take 3, D max Take 20;

[0124] Step 5.3: Build the input and output direct connection structure;

[0125] For all input neurons and output neurons, a direct connection structure is constructed, and the connection weight is expressed as Among them, v d represents the connection weights between d input neurons and all output neurons, and the corresponding initial value is set to 1;

[0126] Step 5.4: Network growth process;

[0127] Hidden layer neurons are added in pairs to compensate for the maximum error, represented by neurons H+1 and H+2, respectively, and their centers, widths, and weights are set as:

[0128]

[0129] σ H+1 =σ H+1 =1 (8)

[0130] w H+1 =1, w H+2 =-1 (9)

[0131] in, is the input sample with the largest network error, k max The calculation formula is:

[0132]

[0133] Where E=[E1,…,E k ,…,E K ] T is the network error vector calculated for all samples, and the error of the kth sample is obtained by the following formula:

[0134]

[0135] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, o s (k) and are the expected output and the predicted output respectively; therefore, the number of hidden neurons H increases by 2;

[0136] Step 5.5: Online learning;

[0137] The network parameters are updated using the improved Levenberg-Marquardt (LM) algorithm; let Δ contain all the parameters of the network to be trained (including the center C and width σ of the neuron, the connection weight W from the hidden layer to the output layer, and the connection weight V from the input layer to the output layer):

[0138]

[0139] Compute the Jacobian matrix row vector j s (k):

[0140]

[0141] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, c h is the center of the hth hidden layer neuron, σ h is the width of the hth hidden layer neuron, w h,s is the weight from the hth hidden layer neuron to the sth output neuron, v d,s is the weight from the dth input neuron to the sth output neuron.

[0142] The Jacobian vector components can be obtained according to the chain rule as follows:

[0143]

[0144]

[0145] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, c h is the center of the hth hidden layer neuron, σ his the width of the hth hidden layer neuron, w h,s is the weight from the hth hidden layer neuron to the sth output neuron, v d,s is the weight from the dth input neuron to the sth output neuron, is the kth input sample, is the d-th dimension input in the k-th input sample, φ h is the radial basis function.

[0146] Calculate the Hessian matrix Q. In order to reduce the space complexity and time complexity of the Jacobi matrix calculation process, the Hessian matrix is ​​converted into multiple Hessian sub-matrices q s (k) and:

[0147]

[0148] Among them, j s (k) is the row vector of the Jacobian matrix. Similarly, the calculation of the gradient vector g is also converted into K gradient sub-vectors η s (k) and:

[0149]

[0150] Among them, e s (k) is the error of the s-th output neuron for the k-th sample, j s (k) is the Jacobian matrix row vector;

[0151] The network structure parameter update rules are as follows:

[0152] Δ i+1 =Δ i -(Q i +μ i I) -1 g i (twenty three)

[0153] Among them, Δ i+1 and Δ i are the network structure parameters for the i+1th and ith trainings, respectively, Q i is a Hessian-like matrix, I is the identity matrix, μ i is the learning rate parameter for the i-th training, which is adaptively adjusted according to the change of training error (when RMSE increases, μ i Zooming in ten times, when RMSE decreases, μ i Reduced tenfold), g i is the gradient vector; the training process is repeated iteratively until i reaches the maximum number of iterations S max(Set according to experience, generally set to 100); calculate the root mean square error RMSE at the end of training. If RMSE converges (less than the set threshold 0.01) or the network hidden layer presets the maximum structure H max (according to experience, it is usually set to 10), then the training stops; otherwise, it returns to the growth type construction of the neural network;

[0154] Step 5.6: Calculate multiple columns of neural network output;

[0155] After the C sub-networks are constructed, the sub-network outputs are integrated to obtain the M-step output of the multi-column neural network.

[0156] Step 6: Calculate the fitness function;

[0157] Calculate the fitness function of all particles; the formula is as follows:

[0158] f=L RMSE +L DTW (twenty four)

[0159] Among them, L RMSE is the individual loss, L DTW is the trend loss; L RMSE The average RMSE of all forecast horizons is used to evaluate the forecast error, and the calculation formula is:

[0160]

[0161] Among them, y t+m (k) and are the expected output and predicted output of the mth step of the kth sample, respectively, K is the number of training samples; L DTW The shape difference between the expected time series and the predicted time series is measured by the dynamic time warping (DTW) algorithm, and the calculation formula is:

[0162]

[0163] Among them, m1=1,2,...,M, m2=1,2,...,M are the step positions of expected output and predicted output, Predict the output of step m1 for the k-th sample, is the expected output of the m2th step of the kth sample, is the cumulative distance of the optimal path from (1,1) to (M,M), and the goal of DTW is to obtain the minimum Euclidean distance between two sequences. value;

[0164] The fitness value is determined by the prediction accuracy of the sub-network. After calculating the fitness values ​​of all particles, the optimal position P of the particle can be determined. bestand the optimal position of the population G best ;

[0165] Step 7: Particle parameter update;

[0166] The update rules for the velocity V of all particles are as follows:

[0167] V(i+1)=wV(i)+C1R1(P best (i)+B(i))+C2R2(G best (i)+B(i))(27)

[0168] Among them, R1 and R2 are two random numbers uniformly distributed in the range [0,1], w is the inertia weight, which is 0.8 here; C1 and C2 are two acceleration coefficients, which are 1.5 here, P best is the optimal position of the particle, G best is the optimal position of the population, B i is the position vector of the current particle; by using the sigmoid transformation, the speed is limited to the interval [0,1]:

[0169]

[0170] Among them, v m Represents the mth element in vector V (m=1,2,…,M), and the calculation formula is:

[0171]

[0172] Here, v max is a constant value, set to 4; then, the position of the mth component of the particle is:

[0173]

[0174] Where ζ is a randomly generated real number between 0 and 1;

[0175] After the positions and velocities of all particles are updated, i=i+1;

[0176] Step 8: Determine the termination condition of multi-step prediction;

[0177] If i = I max Or the fitness function f converges (less than the set threshold 0.01), then the optimal task decomposition G for multi-step prediction is output best , and the M-step prediction value; if not, return to step 4 until the termination condition is reached; usually, the target f value of the real data set ranges from 0.001 to 0.1, and this experiment takes 0.035;

[0178] Step 9: Prediction of total phosphorus in effluent;

[0179] The test sample data is used as the input of the trained multi-column neural network based on the attention mechanism to obtain the output of the network, which is then denormalized to obtain the predicted value of the effluent total phosphorus concentration.

[0180] In this embodiment, Figure 1 The sub-graphs are the results of the 1st, 4th, 7th, 10th, 13th, 16th, and 18th steps of the multi-step prediction respectively; X-axis: test sample, unit is piece; Y-axis: predicted value of effluent total phosphorus concentration, unit is mg / L, the solid line is the measured value of effluent total phosphorus concentration, and the dotted line is the predicted value of effluent total phosphorus concentration; it can be seen that the predicted value of effluent total phosphorus concentration and its actual value curve of the present invention are basically fitted, with high accuracy, and accurate multi-step prediction of effluent total phosphorus can be achieved. Table 1 is the effluent total phosphorus concentration predicted by the present invention.

[0181] Table 1. Total phosphorus concentration in effluent predicted by the present invention (mg / L)

[0182]

[0183]

[0184]

[0185]

Claims

1. A multi-step prediction method for effluent total phosphorus concentration based on a multi-column neural network with an attention mechanism, characterized in that: include: Step 1: Initialization; set the maximum number of iterations I max and population size Q, prediction step number M, and set the initial iteration number i = 1; Step 2: Reconstruct the input sample based on the attention mechanism; Collect actual water quality parameter data of the sewage treatment plant, the predicted variable is the effluent total phosphorus Y, and the auxiliary variable X=[X1,…,X n ,…,X N ], where X n represents the nth auxiliary variable, and N represents the number of auxiliary variables; Step 2.1: Data normalization; Normalize the obtained auxiliary variables X and predictor variables Y to [0,1]: in, and Respectively represent the normalized nth auxiliary variable and output variable; min(X n ) and max(X n ) represent the minimum and maximum values ​​of the nth auxiliary variable, min(Y) and max(Y) represent the minimum and maximum values ​​of the predictor variable, respectively; Step 2.2: Calculate the attention score for each auxiliary variable; for the nth auxiliary variable Calculate its effect on the predictor variable The attention score ξ n for: Step 2.3: Based on the attention score, assign attention weights to each auxiliary variable; use the softmax function to calculate according to formula (4) The attention weight α n for: a n =softmax(ξ n (4) Step 2.4: Reconstruct the input sample; the network input Z is determined by the sum of the weighted auxiliary variables and the predictor variables, and the calculation formula is: Among them, λ is the weight coefficient, which is used to adjust the influence of the prediction variable and the auxiliary variable, and is set in the range of [0,0.4] according to experience; therefore, for the multi-step prediction of M steps, assuming that the input variable Z at the current moment and the previous D-1 time points is Z = [z t-D+1 ,…,z t-D+d ,…,z t ] T The predictor variable Y used to predict the next M steps is [y t+1 ,…,y t+m ,…,y t+M ] T , then the input dimension is D, and the corresponding output dimension is M; Step 3: Task coding; For the M-step prediction task, binary particle encoding is used, and for each particle, its position is represented as B = [b1,…,b m ,…,b M ], where each element represents a prediction step number, with a value of 0 or 1; at the initial moment, the positions of all particles are randomly generated in binary form, and their corresponding velocities V are randomly initialized in the range of [-4,4]; Step 4: Task decoding; For a given position B, consecutive elements with the same value are considered to belong to the same subtask; assuming that B generates C subtasks, a vector representing the subtask allocation will be obtained, denoted by S = [s1,...,s c ,...,s C ](c=1,2,...,C), where s c represents the prediction step size of the cth subtask; Step 5: Design a multi-column neural network model; After task decoding, each subtask is predicted by a sub-network, thus forming a multi-column neural network consisting of C sub-networks, and the output dimensions of each sub-network are s1,...,s c ,...,s C ; Taking the cth sub-network as an example, its network structure is designed as follows; Step 5.1: Initialize network parameters; The network does not initially set up cross-layer direct connection between the input layer and the output layer, that is, the connection weight matrix V between the input neuron and the output neuron is 0, and the number of neurons in the hidden layer of the network is H = 0; the maximum number of neurons in the hidden layer is set to H max , the maximum number of iterations S max ; Step 5.2: Determine network input and output; Taking the cth subnetwork as an example, the input dimension D c Determined by the following formula: Among them, D min and D max D c The boundary value of the network output dimension is s c , from m c-1 +1 to m c ,in For simplicity, use represents the network input, Represents network output; Step 5.3: Build the input and output direct connection structure; For all input neurons and output neurons, a direct connection structure is constructed, and the connection weight is expressed as Among them, v d represents the connection weights between d input neurons and all output neurons, and is set to 1 accordingly; Step 5.4: Network growth process; Hidden layer neurons are added in pairs to compensate for the maximum error, represented by neurons H+1 and H+2, respectively, and their centers, widths, and weights are set as: s H+1 =s H+1 =1(8) In H+1 =1,in H+2 =-1(9) in, is the input sample with the largest network error, k max The calculation formula is: Where E=[E1,…,E k ,…,E K ] T is the network error vector calculated for all samples, and the error of the kth sample is obtained by the following formula: Among them, e s (k) is the error of the s-th output neuron for the k-th sample, o s (k) and are the expected output and the predicted output respectively; therefore, the number of hidden neurons H increases by 2; Step 5.5: Online learning; The network parameters are updated using the improved Levenberg-Marquardt (LM) algorithm; let Δ contain all the parameters of the network to be trained (including the center C and width σ of the neuron, the connection weight W from the hidden layer to the output layer, and the connection weight V from the input layer to the output layer): Compute the Jacobian matrix row vector j s (k): Among them, c h is the center of the hth hidden layer neuron, σ h is the width of the hth hidden layer neuron, w h,s is the weight from the hth hidden layer neuron to the sth output neuron, v d,s is the weight from the dth input neuron to the sth output neuron; The Jacobian vector components can be obtained according to the chain rule as follows: in, is the kth input sample, is the d-th dimension input in the k-th input sample; Calculate the Hessian matrix Q. In order to reduce the space complexity and time complexity of the Jacobi matrix calculation process, the Hessian matrix is ​​converted into multiple Hessian sub-matrices q s (k) and: Similarly, the calculation of the gradient vector g is also converted into K gradient sub-vectors η s (k) and: The network structure parameter update rules are as follows: D i+1 =D i -(Q i +m i I) -1 g i (23) Among them, Q i is a Hessian-like matrix, I is the identity matrix, μ is the learning rate parameter, g i is the gradient vector; the training process is repeated iteratively until the maximum number of iterations S is reached max , set to 100; calculate the root mean square error RMSE at the end of training. If RMSE converges or the maximum structure H of the network hidden layer is preset max , then stop the training; otherwise, return to the growth type construction of the neural network; Step 5.6: Calculate multiple columns of neural network output; After the C sub-networks are constructed, the sub-network outputs are integrated to obtain the M-step output of the multi-column neural network. Step 6: Calculate the fitness function; Calculate the fitness function of all particles; the formula is as follows: f=L RMSE +L DTW (24) Among them, L RMSE is the individual loss, L DTW is the trend loss; L RMSE The average RMSE of all forecast horizons is used to evaluate the forecast error, and the calculation formula is: Among them, y t+m (k) and are the expected output and predicted output of the mth step of the kth sample, respectively, K is the number of training samples; L DTW The shape difference between the expected time series and the predicted time series is measured by the dynamic time warping (DTW) algorithm, and the calculation formula is: Among them, m1=1,2,...,M, m2=1,2,...,M are the step positions of expected output and predicted output, Predict the output of step m1 for the k-th sample, is the expected output of the m2th step of the kth sample, is the cumulative distance of the optimal path from (1,1) to (M,M), and the goal of DTW is to obtain the minimum Euclidean distance between two sequences. value; After calculating the fitness values ​​of all particles, the optimal position P of the particle can be determined. best and the optimal position of the population G best ; Step 7: Particle parameter update; The update rules for the velocity V of all particles are as follows: V(i+1)=wV(i)+C1R1(P best (i)+B(i))+C2R2(G best (i)+B(i))(27) Among them, R1 and R2 are two random numbers uniformly distributed in the range [0,1], w is the inertia weight, C1 and C2 are two acceleration coefficients, P best is the optimal position of the particle, G best is the optimal position of the population; by using the sigmoid transformation, the speed is limited to the interval [0,1]: Among them, v m Represents the mth element in vector V (m=1,2,…,M), and the calculation formula is: Here, v max is a constant value, set to 4; then, the position of the mth component of the particle is: Where ζ is a randomly generated real number between 0 and 1; After the positions and velocities of all particles are updated, i=i+1; Step 8: Determine the termination condition of multi-step prediction; If i = I max Or the fitness function f converges, that is, it is less than the set threshold 0.01, then the optimal task decomposition G for multi-step prediction is output best , and the M-step predicted value; if not, return to step 4 until the termination condition is reached; Step 9: Prediction of total phosphorus in effluent; The test sample data is used as the input of the trained multi-column neural network based on the attention mechanism to obtain the output of the network, which is then denormalized to obtain the predicted value of the effluent total phosphorus concentration.

Images