Heat exchange station thermal load modeling method integrated with echo state network

Abstract

The invention discloses a heat exchange station thermal load modeling method integrated with an echo state network, and belongs to the technical field of industrial process intelligent detection. According to the method, dynamic time sequence characteristics and a nonlinear mapping relation of industrial process data are extracted by constructing a multi-layer echo state network (ESN) deep architecture; performing adaptive optimization on the key hyper-parameters of each layer of ESN by adopting a particle swarm optimization algorithm so as to obtain optimal reservoir dynamic characteristics and output weight; through a layer-by-layer progressive deep learning strategy, cascading the storage pool state of the front-layer ESN with an original input variable to serve as input of a subsequent layer, and realizing layer-by-layer abstraction and representation learning of features; and finally, adopting an average integration strategy to fuse the prediction output of the multi-layer ESN, and improving the prediction precision of the model. The deep ESN ensemble learning framework provided by the invention not only can fully mine complex dynamic laws in industrial time series data, but also can effectively inhibit an overfitting phenomenon of a single model, and achieves the purpose of improving the robustness and reliability of a soft measurement model.

Description

Technical Field

[0001] This invention belongs to the field of intelligent detection technology for industrial processes, specifically relating to a method for modeling the heat load of a heat exchange station by integrating an echo state network. Background Technology

[0002] With the increasing demands for modern industrial process control and energy conservation and emission reduction, district heating systems, as a crucial component of urban infrastructure, directly impact residents' quality of life and the achievement of environmental protection goals through their operational efficiency and energy consumption levels. In district heating systems, heat load forecasting is a key technology for achieving intelligent regulation and optimized operation. However, heat load measurement is often affected by various factors, such as measuring equipment failure, communication interruptions, and sensor drift, leading to inaccurate measurements or data gaps, thus affecting the normal operation of the system. Soft measurement technology, as an effective means to solve this problem, can achieve real-time estimation of key parameters by establishing a mathematical model between easily measurable process variables and difficult-to-measure target variables.

[0003] Traditional soft sensor models mainly include mechanism-based models and data-driven models. Mechanism-based models require a deep understanding of the physicochemical principles of the process, making modeling complex and limiting their generalization ability. On the other hand, data-driven models, such as Support Vector Machines (SVM) and Artificial Neural Networks (ANN), can make predictions by learning patterns from historical data, but they face problems such as training difficulties and overfitting.

[0004] Echo State Networks (ESNs), a special type of recurrent neural network, were first proposed by Jaeger in 2001. With their simple structure, efficient training, and strong dynamic memory capabilities, they have shown great potential in the field of time series prediction. ESNs mainly consist of an input layer, a reservoir, and an output layer. The reservoir contains a large number of randomly connected neurons with "echo state" properties, enabling it to capture the temporal dynamics of the input sequence. Unlike traditional neural networks, ESNs only require training the output weights, significantly reducing computational complexity and making them suitable for real-time prediction applications.

[0005] However, single ESN models still have certain limitations when dealing with complex nonlinear systems: the reservoir topology significantly affects model performance, and traditional random connection structures struggle to capture multi-scale features in complex time series; furthermore, the random initialization characteristic of ESNs introduces uncertainty into model performance. To address this issue, Liu et al. (2020) proposed an ESN optimization method based on the binary gray wolf algorithm, which improved model performance by optimizing reservoir parameters; Gong et al. (2020) designed an improved small-world topology to optimize the network connection pattern of ESNs.

[0006] Ensemble ESNs have also been extensively studied as another strategy to improve model performance. Hu et al. (2020) designed an ensemble echo network with a deep architecture, which effectively reduced prediction variance and improved model stability by combining multiple ESNs with different parameters. Liu et al. (2023) proposed a memory-enhanced echo state network, which enhanced the model's ability to utilize historical information by introducing a memory module, showing advantages in multiple time series forecasting tasks. Trierweiler Ribeiro et al. (2020) applied Bayesian optimization methods to ESNs, developing an optimized ESN framework for short-term load forecasting, which improved prediction accuracy.

[0007] Despite these shortcomings, existing ESN methods still suffer from the following problems: a single ESN model struggles to simultaneously account for both short-term fluctuations and long-term trends, leading to inconsistent prediction performance across different time scales; traditional ESN hyperparameter tuning methods rely on expert experience or grid search, making it difficult to find the global optimum; the inter-layer information transfer mechanism is inefficient, failing to fully utilize the advantages of the multi-layer structure; and the model lacks robustness, being highly sensitive to input noise and system disturbances. Therefore, a novel modeling method that can comprehensively leverage the advantages of ESN while overcoming these drawbacks is urgently needed. Summary of the Invention

[0008] In view of the above-mentioned technical problems in the prior art, the present invention proposes a heat load modeling method for heat exchange stations with integrated echo state network. The method is reasonably designed, overcomes the shortcomings of the prior art, and has good results.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] A method for modeling heat load in a heat exchange station using an integrated echo state network includes the following steps:

[0011] Step 1: Collect measured values ​​of variables in the heating process and operating parameters of the heat exchange station, and perform smoothing and normalization processing on the input data;

[0012] Step 2: Construct the first layer of Echo State Network, set a fixed random seed and initialize the input weight matrix and the reservoir connection matrix. Use the particle swarm optimization algorithm to optimize the key hyperparameters of the input, including scaling factor, spectral radius and ridge regression regularization coefficient. Train the output weights through ridge regression to obtain the first layer ESN model.

[0013] Step 3: Concatenate the reservoir state of the first-layer ESN with the original input variables as the input variables of the second-layer ESN, reinitialize the inter-layer connection weight matrix of the second layer, and use the particle swarm optimization algorithm to optimize the input scaling factor, spectral radius and regularization coefficient of the second-layer ESN to train the second-layer ESN model.

[0014] Step 4: Repeat the modeling process in Step 3, and build a deep ESN integrated architecture layer by layer. Each layer will concatenate the state of the previous layer's ESN reservoir with the original input variables as the input of the current layer, and re-optimize the hyperparameters and train the model until a multi-layer ESN network with the preset number of layers is built.

[0015] Step 5: After each layer of ESN training is completed, immediately calculate the independent prediction performance of that layer, and at the same time use the average ensemble strategy to fuse all prediction outputs from the first layer to the current layer with equal weights to achieve multi-model collaborative prediction.

[0016] Step 6: Collect test data and normalize it. Use the trained multi-layer ESN ensemble model to make predictions. Output the single-layer prediction results and the ensemble prediction results for performance comparison and analysis.

[0017] Preferably, the specific process of step 1 is as follows:

[0018] Step 1.1: Let Y be the initial set of collected return water temperatures. org =[y org,1 ,y org,2 ,…,y org,i ,…,y org,n ] T ∈R n Let U be the initial set of process variable measurements collected. org =[u org,1 ,u org,2 ,…,u org,i ,…,u org,n ] T ∈R n×d ;

[0019] Among them, y org,i Let u be the i-th measured value of the return water temperature in the primary network, and n be the total number of samples. org,i Let d be the vector of the i-th process variable, and d be the dimension of the input variables;

[0020] Step 1.2: Perform square root transformation preprocessing on the return water temperature data of the primary network according to formula (1);

[0021] Y trans =(Y org ) 1 / 2 (1);

[0022] In the formula, Y trans The dataset after square root transformation; Y org This is the original dataset;

[0023] Step 1.3: Smooth the input process variables using the sliding window method to obtain u. smooth,i =smooth(u org,i ,n window ), where n window U is the length of the smooth window; where u smooth,i This represents the smoothed process variable vector obtained after smoothing the i-th process variable vector;

[0024] Step 1.4: Process the smoothed input data U smooth and output data Y org Normalization is performed separately, and the normalization formula is as follows:

[0025]

[0026] In the formula, u min u max These are the minimum and maximum values ​​of each input variable in the entire sample, y. min y max The minimum and maximum values ​​of the output variable are respectively, and the input data U is obtained after normalization. nor and output data Y nor ;

[0027] Step 1.5: Divide the normalized data into a training set and a test set, with N training samples. train The number of test samples is N test =nN train The length of the washout section is set to eliminate the influence of the initial state of the ESN.

[0028] Preferably, the specific process of step 2 is as follows:

[0029] Step 2.1: Set the number of layers in the multi-layer ESN network to NL, and the number of neurons in each layer's reserve pool to N. r The first layer input weight matrix is ​​randomly initialized within the range of [-1,1], as shown in formula (4):

[0030]

[0031] In the formula, W in N is the input weight matrix for the first layer of the ESN; r N represents the number of neurons in the reserve pool. u The dimension of the input variables; +1 indicates the bias term;

[0032] Step 2.2: Randomly initialize the internal connection weight matrix of each layer, as shown in formula (5):

[0033]

[0034] In the formula, W (l) Let N be the internal connection weight matrix of the l-th layer ESN; l = 1, 2, ..., NL is the layer index; N r This represents the number of neurons in the reserve pool.

[0035] Step 2.3: Randomly initialize the inter-layer connection weight matrix to connect the (l-1)th layer and the lth layer; as shown in formula (6):

[0036]

[0037] In the formula, N is the inter-layer connection weight matrix of layer l; r N represents the number of neurons in the reserve pool. u The dimension of the input variables; +1 indicates the bias term;

[0038] Step 2.4: Set a random seed to ensure that the experimental results are reproducible, and initialize each weight matrix according to the uniform distribution U(-1,1).

[0039] Preferably, the specific process of step 3 is as follows:

[0040] Step 3.1: For the first layer ESN, set the hyperparameter vector to be optimized as θ1=[α1,ρ1,β1], where α1 is the input scaling factor, ρ1 is the spectral radius, and β1 is the regularization coefficient;

[0041] Step 3.2: Find the optimal combination of hyperparameters using the particle swarm optimization algorithm. The particle velocity update is shown in formula (7):

[0042]

[0043] In the formula, Let be the velocity of the i-th particle in the (k+1)th iteration; The velocity is the velocity of the k-th iteration; w is the inertia weight; c1 and c2 are learning factors; r1 and r2 are random numbers in the range [0,1]; p i g represents the historical best position of the i-th particle; g represents the global best position. Let be the position of the i-th particle in the k-th iteration;

[0044] The particle position is updated as shown in Equation (8):

[0045]

[0046] In the formula, Let be the position of the i-th particle in the (k+1)-th iteration; This represents the position of the k-th iteration. The speed of the (k+1)th iteration;

[0047] Step 3.3: Set the fitness function to the root mean square error of the validation set, as shown in formula (9):

[0048]

[0049] In the formula, RMSE is the root mean square error; N val To verify the sample size; y i The actual value; This is a predicted value;

[0050] Step 3.4: For the l-th layer ESN (l = 2, 3, ..., NL), repeat steps 3.1-3.3 to obtain the optimal hyperparameter θ for each layer. l =[α l ,ρ l ,β l ].

[0051] Preferably, the parameters of the particle swarm optimization algorithm are set as follows: the particle swarm size is set to 20-50, the maximum number of iterations is set to 50-200, the inertia weight w∈[0.4,0.9], the learning factor c1=c2=2.0, the hyperparameter search range is the input scaling factor α∈[0.1,2.0], the spectral radius ρ∈[0.1,1.0], the regularization index β∈[10,30], and the regularization coefficient is calculated as shown in formula (10):

[0052] β reg =2 -round(β) (10);

[0053] In the formula, β reg The regularization coefficient is used.

[0054] Preferably, the specific process of step 4 is as follows:

[0055] Step 4.1: Normalize the first-layer input weight matrix using the obtained optimal input scaling factor, as shown in formula (11):

[0056]

[0057] In the formula, W in Let be the normalized input weight matrix; α1 is the optimal input scaling factor; ||·|| denotes the Frobenius norm of the matrix;

[0058] Step 4.2: Normalize the first layer internal connection weight matrix using the obtained optimal spectral radius, as shown in formula (12):

[0059]

[0060] In the formula, W (1) ρ1 is the normalized weight matrix of the first layer internal connection; eig(·) represents the matrix eigenvalue operation; max(|·|) represents taking the maximum absolute value.

[0061] Step 4.3: Calculate the state matrix of the first-layer ESN. For time t, the state update equation is shown in formula (13):

[0062] X (1) (t)=tanh(W in [u(t);1]+W (1) X (1) (t-1)1) (13);

[0063] In the formula, X (1) u(t) represents the state of the first-layer ESN at time t; u(t) represents the input vector at time t; tanh is the hyperbolic tangent activation function; [u(t); 1] represents the concatenation of the input vector and the bias term; W in W is the input weight matrix. (1) The internal connection weight matrix; initial state X (1) (0) = 0;

[0064] Step 4.4: Augment the first-layer state matrix with the original input and bias terms, as shown in Equation (14):

[0065]

[0066] In the formula, X is the augmented state matrix of the first layer; (v) Here is the first-layer ESN state matrix; u is the original input variable; 1 indicates the bias term;

[0067] The first layer output weights are calculated using ridge regression, as shown in formula (15):

[0068]

[0069] In the formula, This is the output weight matrix for the first layer; Y is the augmented state matrix of the first layer; β1 is the optimal regularization coefficient; I is the identity matrix; Y train Labels for the training set;

[0070] Step 4.5: For the l-th layer ESN (l = 2, 3, ..., NL), use the augmented vector of the previous layer state and the original input as input, as shown in Equation (16):

[0071] Input l =[X (l-1) ,u] (16);

[0072] In the formula, Input l X is the input of the l-th layer; (l-1) Here is the state matrix of the (l-1)th layer; u is the original input variable;

[0073] The inter-layer connection weights are normalized using the obtained optimal inter-layer scaling factor, as shown in Equation (17):

[0074]

[0075] In the formula, This is the normalized inter-layer connection weight matrix; α l Let be the optimal input scaling factor for the l-th layer; ||·|| denotes the Frobenius norm of the matrix;

[0076] The internal connection weights are normalized using the optimal spectral radius, as shown in Equation (18):

[0077]

[0078] In the formula, W (l) ρ is the normalized internal connection weight matrix of the l-th layer; l Let be the optimal spectral radius of the l-th layer; eig(·) represents the matrix eigenvalue operation; max(|·|) represents taking the maximum absolute value;

[0079] Step 4.6: Calculate the state matrix of the l-th layer ESN and the state update equation, as shown in formula (19):

[0080]

[0081] In the formula, X (l) (t) represents the state of the l-th layer ESN at time t; X (l-1) u(t) represents the state of layer l-1 at time t; u(t) represents the input vector at time t; tanh is the hyperbolic tangent activation function; W is the inter-layer connection weight matrix. (l) The internal connection weight matrix; initial state X (l)(0) = 0;

[0082] The output weights of the l-th layer are calculated as shown in formula (20):

[0083]

[0084] In the formula, This is the output weight matrix for the l-th layer; Let β be the augmented state matrix of the l-th layer; l Y is the optimal regularization coefficient for the l-th layer; I is the identity matrix; train Labels for the training set; This indicates that the augmented state matrix of the l-th layer consists of the l-th layer state matrix, the original input, and the bias term.

[0085] Preferably, the specific process of step 5 is as follows:

[0086] Step 5.1: Calculate the prediction output of each layer of ESN. The prediction result of the l-th layer is shown in formula (21):

[0087]

[0088] In the formula, This represents the prediction result for the L-th layer ESN; Let l be the augmented state matrix of the l-th layer; This is the output weight matrix for the l-th layer;

[0089] Step 5.2: The prediction results of each layer are fused using an average ensemble strategy, and the prediction results of each layer are directly averaged, as shown in formula (22):

[0090]

[0091] In the formula, This represents the average ensemble prediction result; NL represents the total number of layers in the ESN. This represents the prediction result for the L-th layer ESN;

[0092] Step 5.3: Average the squared prediction results of each layer, as shown in formula (23):

[0093]

[0094] In the formula, The result is the ensemble prediction average after squared calculation; NL is the total number of layers in the ESN. This represents the prediction result for the L-th layer ESN;

[0095] Step 5.4: Select the better-performing integration method as the final soft measurement prediction model.

[0096] Preferably, the specific process of step 6 is as follows:

[0097] Step 6.1: Collect test data and perform the same normalization and smoothing processes as the training data to obtain the normalized test set U. test ;

[0098] Step 6.2: Calculate the state matrix of each layer of the test set using the trained multilayer ESN model. and augmented state matrix

[0099] Step 6.3: Calculate the test set prediction results using the trained output weights of each layer, as shown in formula (24):

[0100]

[0101] In the formula, The prediction results of the layer l ESN on the test set; Let l be the augmented state matrix of the l-th layer on the test set; The output weight matrix trained for the l-th layer;

[0102] Step 6.4: Integrate the test prediction results of each layer to obtain the final prediction output, as shown in formula (25):

[0103]

[0104] In the formula, The results are the ensemble predictions on the test set; NL represents the total number of ESN layers. The prediction results of the layer l ESN on the test set;

[0105] Step 6.5: Perform an inverse transformation on the prediction results to obtain the actual predicted value of the primary network return water temperature, as shown in formula (26):

[0106]

[0107] In the formula, Y pred This is the actual predicted value; The results are the integrated predictions in the normalized space.

[0108] Preferably, the root mean square error is as shown in formula (27):

[0109]

[0110] In the formula, MSE is the root mean square error; N is the sample size; y i This represents the true value of the i-th sample. Let be the predicted value for the i-th sample;

[0111] The coefficient of determination is shown in formula (28):

[0112]

[0113] In the formula, R 2 The coefficient of determination; SS res Sum of squared residuals; SS tot The total sum of squares;

[0114] in

[0115]

[0116] In the formula, SS res y is the sum of squared residuals; i This represents the true value of the i-th sample. Let be the predicted value of the i-th sample; N is the number of samples;

[0117]

[0118] In the formula, SS tot For the total sum of squares; y i This represents the true value of the i-th sample. The mean of the true values ​​is N; the sample size is N.

[0119] The mean absolute error is shown in formula (31):

[0120]

[0121] In the formula, MAE is the mean absolute error; N is the sample size; y i This represents the true value of the i-th sample. Let be the predicted value for the i-th sample.

[0122] Network structure parameters are set as follows: the number of layers (NL) in the multi-layer ESN network is set to 8-15 layers, and the size of the reservoir in each layer (N) is... r Set the values ​​to 50-400, the training / validation ratio gamma to 0.6-0.8, the washout length to 50-200, the data sampling interval ND to 1, and the smoothing window length to 3.

[0123] The beneficial technical effects of this invention are as follows:

[0124] This invention extracts dynamic temporal features and nonlinear mapping relationships from industrial process data by constructing a multi-layered echo state network (ESN) deep architecture. It employs a particle swarm optimization algorithm to adaptively optimize the key hyperparameters of each ESN layer to obtain the optimal reservoir dynamics and output weights. Through a layer-by-layer progressive deep learning strategy, the reservoir states of the previous ESN layer are concatenated with the original input variables as input to subsequent layers, achieving layer-by-layer feature abstraction and representation learning. Finally, an average ensemble strategy is used to fuse the prediction outputs of the multi-layered ESNs, improving the model's prediction accuracy. The deep ESN ensemble learning framework provided by this invention not only fully explores the complex dynamic patterns in industrial time-series data but also effectively suppresses overfitting of single models, thereby improving the robustness and reliability of soft measurement models. Attached Figure Description

[0125] Figure 1 This is a flowchart of the method of the present invention.

[0126] Figure 2 This is a comparison chart of the performance indicators of the multi-layer ESN integrated model and each single-layer ESN model in the embodiments of the present invention;

[0127] in, Figure 2 In the diagram, (a) represents the R of each layer of the ESN. 2 Figure (a) is a curve comparing the values; Figure (b) is a comparison of the prediction performance of the optimal model on the test set.

[0128] Figure 3 This is a graph showing the predictive performance evaluation of the ESN model on the training data in an embodiment of the present invention.

[0129] in, Figure 3 Figure (a) is a comparison of time series predictions on the training set between a single-layer and a multi-layer ESN; Figure (b) is a histogram of the prediction error distribution on the training set.

[0130] Figure 4 This is a graph showing the predictive performance evaluation of the ESN model on test data in an embodiment of the present invention;

[0131] in, Figure 4 Figure (a) is a comparison of time series predictions on the test set between single-layer and multi-layer ESNs; Figure (b) is a histogram of prediction error distribution on the test set. Detailed Implementation

[0132] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0133] This embodiment uses actual operating data of a centralized heating system in a city in northern my country as an example for detailed explanation.

[0134] Central heating systems are a crucial component of urban infrastructure. Their primary function is to transport hot water generated by a heat source through a primary network to various heat exchange stations, where the heat is then transferred to the secondary network, ultimately providing stable heating services to residents and businesses. Heat exchange stations, acting as heat exchange hubs between the primary and secondary networks, are typically equipped with key devices such as plate heat exchangers, circulating pumps, regulating valves, and temperature and pressure sensors. A host computer system enables real-time monitoring and adjustment of parameters such as water temperature and flow rate. In actual operation, this system uses an automated system to collect key process data such as temperature, pressure, valve opening, and pump frequency every 10 minutes. Due to physical transmission and thermal inertia involved in the flow of hot water from the heat exchange station to the secondary network users and the completion of heat exchange, heat transfer and system response exhibit significant time lags. This means that the primary network return water temperature and the secondary network supply water temperature cannot be directly obtained from parameters such as the current primary network supply water temperature, valve opening, and pump frequency, but are instead influenced by the combined operating conditions of previous moments. Furthermore, heating systems exhibit time-varying characteristics (e.g., climate and weather changes affect heat demand), thermal inertia (long heat transfer times in pipes and equipment), time lag (changes in the primary network cannot be immediately reflected in the secondary network), and strong coupling between equipment, making regulation difficult. Therefore, to estimate the heat load of heat exchange stations in real time and provide a scientific basis for subsequent water temperature and flow rate regulation, it is necessary to establish a soft-sensor model based on time-series data to estimate the system heat load in real time. Additionally, industrial field data is often affected by factors such as sensor accuracy and network fluctuations, resulting in noise and anomalous jumps. To ensure data quality, the raw data needs to be smoothed and outliers removed to improve the stability and accuracy of subsequent modeling.

[0135] The specific steps of this invention will be described in conjunction with the operation process of a centralized heating system, such as... Figure 1 As shown:

[0136] Step S1, Offline Modeling Stage: Use the collected data as the training dataset and perform preprocessing.

[0137] All input variables are smoothed using a sliding window with a width of 3 to eliminate noise and anomalous jumps. Let the normalized process variable measurements be Unor = [u1, u2, ..., u...]. n ] T ∈R 4675×4 The output variable is Ynor=

[0138] [y1,y2,…,y n ] T ∈R 4675 All variables were normalized to the [0,1] interval. The washout length of the training data was set to washout=100, and the sampling interval was ND=1.

[0139] Step S2: Set a fixed random seed and initialize the input weight matrix Win1 and the reservoir connection matrix W_1 of the first-layer ESN. Use Particle Swarm Optimization (PSO) to optimize key hyperparameters such as the input scaling factor α1, spectral radius ρ1, and ridge regression regularization coefficient λ1. The PSO parameters are set as follows: particle swarm size 200, maximum number of iterations 50, learning factors c1 = 1, c2 = 1.3, inertia weight β = 0.4, initial input scaling factor α = 0.05, and regularization parameter λ = 1. Normalize the weight matrix using the optimized hyperparameters. Obtain the reservoir state X1 of the first-layer ESN through an external input-driven network. After augmenting it with the original input, train the output weight Wout1 using ridge regression to obtain the first-layer ESN model.

[0140] Step S3: Concatenate the state matrix X1 of the first-layer ESN with the original input variable u as the input of the second-layer ESN. Reinitialize the inter-layer connection weight matrix Wil1 of the second layer, and optimize the input scaling factor α2, spectral radius ρ2, and ridge regression regularization coefficient β2 of the second layer using the Particle Swarm Optimization (PSO) algorithm. The PSO parameters are set as follows: particle swarm size 200, maximum number of iterations 50, learning factors c1 = 1, c2 = 1.3, and inertia weight β = 0.4. Normalize the weight matrix using the optimized parameters, drive the network to obtain the state matrix X2 of the second-layer ESN, and train the output weight Wout2 using the ridge regression method to complete the training of the second-layer ESN model.

[0141] Step S4: Repeat the modeling process of the second layer to build the deep ESN integration architecture layer by layer. For each layer l (l = 3, ..., 11), the state matrix X of the previous layer ESN is... l-1 The input variable u is concatenated with the original input variable u and used as the input for the current layer. This process is repeated layer by layer until the construction of a multi-layer ESN ensemble network with a preset number of layers L=11 is completed.

[0142] Step S5: After each ESN layer is trained, immediately calculate the independent prediction performance (e.g., R2, RMSE) of that layer on the test set. Using an average ensemble or variance ensemble strategy, fuse all prediction outputs from layer 1 to the current layer with equal weights and calculate the ensemble prediction performance.

[0143] Step S6, Online Testing Phase: Collect test data and obtain the normalized test set Unor after normalization processing. test The model has dimensions of 935×4. A pre-trained multi-layer ESN ensemble model is used to predict the test set, outputting the prediction results of each layer and the ensemble prediction results. The prediction accuracy of each layer and the ensemble model is compared and analyzed, and curves comparing the true and predicted values, as well as the R² change curve, are plotted to verify the effectiveness and superiority of the model.

[0144] This invention uses root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²) to evaluate the estimation performance of soft sensor models. The test set has 935 samples. The formulas for calculating each performance indicator are as follows:

[0145]

[0146] Where, n t =935, y test,i Let y be the true value of the i-th test sample. pred,i These are the model's predicted values. This represents the mean of the true values ​​in the test set.

[0147] This invention compares the performance of three prediction methods: single-layer ESN (without integration), multi-layer ESN average integration, and multi-layer ESN variance integration. The relevant results are shown in Table 1.

[0148] Table 1. Comparison of prediction performance between single-layer ESN and the two integration methods of this invention.

[0149]

[0150] As shown in Table 1, regardless of whether average ensemble or variance ensemble is used, the prediction accuracy (R2) of the model is higher than that of the single-layer ESN (non-ensemble) model, indicating that the multi-layer ensemble strategy can effectively improve the prediction performance and generalization ability of the model.

[0151] Figure 2 The R-values ​​of the ensemble model and each individual ESN layer during the training process are given. 2 In comparison, the ensemble method has better stability and effectively reduces performance fluctuations caused by the randomness of the ESN model.

[0152] Figure 3-4 The error distributions on the training and test sets are presented. Most of the predicted values ​​deviate little from the true values, and the number of samples with large deviations is limited and the range of deviations is narrow. This indicates that the model's estimates are basically consistent with the true values, further verifying the accuracy of this method.

[0153] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for modeling heat load in a heat exchange station using an integrated echo state network, characterized in that, Includes the following steps: Step 1: Collect measured values ​​of variables in the heating process and operating parameters of the heat exchange station, and perform smoothing and normalization processing on the input data; Step 2: Construct the first layer of Echo State Network, set a fixed random seed and initialize the input weight matrix and the reservoir connection matrix. Use the particle swarm optimization algorithm to optimize the key hyperparameters of the input, including scaling factor, spectral radius and ridge regression regularization coefficient. Train the output weights through ridge regression to obtain the first layer ESN model. Step 3: Concatenate the reservoir state of the first-layer ESN with the original input variables as the input variables of the second-layer ESN, reinitialize the inter-layer connection weight matrix of the second layer, and use the particle swarm optimization algorithm to optimize the input scaling factor, spectral radius and regularization coefficient of the second-layer ESN to train the second-layer ESN model. Step 4: Repeat the modeling process in Step 3, and build a deep ESN integrated architecture layer by layer. Each layer will concatenate the state of the previous layer's ESN reservoir with the original input variables as the input of the current layer, and re-optimize the hyperparameters and train the model until a multi-layer ESN network with the preset number of layers is built. Step 5: After each layer of ESN training is completed, immediately calculate the independent prediction performance of that layer, and at the same time use the average ensemble strategy to fuse all prediction outputs from the first layer to the current layer with equal weights to achieve multi-model collaborative prediction. Step 6: Collect test data and normalize it. Use the trained multi-layer ESN ensemble model to make predictions. Output the single-layer prediction results and the ensemble prediction results for performance comparison and analysis.

2. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 1 is as follows: Step 1.1: Let Y be the initial set of collected return water temperatures. org =[y org,1 ,y org,2 ,…,y org,i ,…,y org,n ] T ∈R n Let U be the initial set of process variable measurements collected. org =[u org,1 ,u org,2 ,…,u org,i ,…,u org,n ] T ∈T n×d ; Among them, y org,l Let u be the i-th measured value of the return water temperature in the primary network, and n be the total number of samples. org,i Let d be the vector of the i-th process variable, and d be the dimension of the input variables; Step 1.2: Perform square root transformation preprocessing on the return water temperature data of the primary network according to formula (1); AND trans =(And org ) 1 / 2 (1); In the formula, Y trans The dataset after square root transformation; Y org This is the original dataset; Step 1.3: Smooth the input process variables using the sliding window method to obtain u. smooth,i =smooth(u org,i ,n window ), where n window U is the length of the smooth window; where u smooth,i This represents the smoothed process variable vector obtained after smoothing the i-th process variable vector; Step 1.4: Process the smoothed input data U smooth and output data Y org Normalization is performed separately, and the normalization formula is as follows: In the formula, u min u max These are the minimum and maximum values ​​of each input variable in the entire sample, y. min y max The minimum and maximum values ​​of the output variable are respectively, and the input data U is obtained after normalization. nor and output data Y nor ; Step 1.5: Divide the normalized data into a training set and a test set, with N training samples. train The number of test samples is N test =nN train The length of the washout section is set to eliminate the influence of the initial state of the ESN.

3. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 2 is as follows: Step 2.1: Set the number of layers in the multi-layer ESN network to NL, and the number of neurons in each layer's reserve pool to N. r The first layer input weight matrix is ​​randomly initialized within the range of [-1,1], as shown in formula (4): In the formula, W in N is the input weight matrix for the first layer of the ESN; r N represents the number of neurons in the reserve pool. u The dimension of the input variables; +1 indicates the bias term; Step 2.2: Randomly initialize the internal connection weight matrix of each layer, as shown in formula (5): In the formula, W (l) Let N be the internal connection weight matrix of the l-th layer ESN; l = 1, 2, ..., NL is the layer index; N r This represents the number of neurons in the reserve pool. Step 2.3: Randomly initialize the inter-layer connection weight matrix to connect the (l-1)th layer and the lth layer; as shown in formula (6): In the formula, N is the inter-layer connection weight matrix of layer l; r N represents the number of neurons in the reserve pool. u Dimensions of the input variables; Step 2.4: Set a random seed to ensure that the experimental results are reproducible, and initialize each weight matrix according to the uniform distribution U(-1,1).

4. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 3 is as follows: Step 3.1: For the first layer ESN, set the hyperparameter vector to be optimized as θ1=[α1,ρ1,β1], where α1 is the input scaling factor, ρ1 is the spectral radius, and β1 is the regularization coefficient; Step 3.2: Find the optimal combination of hyperparameters using the particle swarm optimization algorithm. The particle velocity update is shown in formula (7): In the formula, Let be the velocity of the i-th particle in the (k+1)th iteration; The velocity is the velocity of the k-th iteration; w is the inertia weight; c1 and c2 are learning factors; r1 and r2 are random numbers in the range [0,1]; p i g represents the historical best position of the i-th particle; g represents the global best position. Let be the position of the i-th particle in the k-th iteration; The particle position is updated as shown in Equation (8): In the formula, Let be the position of the i-th particle in the (k+1)-th iteration; This represents the position of the k-th iteration. The speed of the (k+1)th iteration; Step 3.3: Set the fitness function to the root mean square error of the validation set, as shown in formula (9): In the formula, RMSE is the root mean square error; N val To verify the sample size; y i The actual value; This is a predicted value; Step 3.4: For the l-th layer ESN, repeat steps 3.1-3.3 to obtain the optimal hyperparameters θ for each layer. l =[α l ,ρ l ,β l ].

5. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 4, characterized in that, The parameters of the particle swarm optimization algorithm are set as follows: the particle swarm size is set to 20-50, the maximum number of iterations is set to 50-200, the inertia weight w∈[0.4,0.9], the learning factor c1=c2=2.0, the hyperparameter search range is the input scaling factor α∈[0.1,2.0], the spectral radius ρ∈[0.1,1.0], the regularization index β∈[10,30], and the regularization coefficient is calculated as shown in formula (10): β reg ＝2 -round(β) (10)； In the formula, β reg is the regularization coefficient.

6. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 4 is as follows: Step 4.1: Normalize the first-layer input weight matrix using the obtained optimal input scaling factor, as shown in formula (11): In the formula, W in Let be the normalized input weight matrix; α1 is the optimal input scaling factor; ||·|| denotes the Frobenius norm of the matrix; Step 4.2: Normalize the first layer internal connection weight matrix using the obtained optimal spectral radius, as shown in formula (12): In the formula, W (1) ρ1 is the normalized weight matrix of the first layer internal connection; eig(·) represents the matrix eigenvalue operation; max(|·|) represents taking the maximum absolute value. Step 4.3: Calculate the state matrix of the first-layer ESN. For time t, the state update equation is shown in formula (13): X (1) (t)=tanh(W in [u(t);1]+W (1) X (1) (t-1)) (13); In the formula, X (1) u(t) represents the state of the first-layer ESN at time t; u(t) represents the input vector at time t; tanh is the hyperbolic tangent activation function; [u(t); 1] represents the concatenation of the input vector and the bias term; W in W is the input weight matrix. (1) The internal connection weight matrix; initial state X (1) (0) = 0; Step 4.4: Augment the first-layer state matrix with the original input and bias terms, as shown in Equation (14): In the formula, X is the augmented state matrix of the first layer; (1) Here is the first-layer ESN state matrix; u is the original input variable; 1 indicates the bias term; The first layer output weights are calculated using ridge regression, as shown in formula (15): In the formula, This is the output weight matrix for the first layer; Y is the augmented state matrix of the first layer; β1 is the optimal regularization coefficient; I is the identity matrix; Y train Labels for the training set; Step 4.5: For the l-th layer ESN, use the augmented vector of the previous layer state and the original input as input, as shown in Equation (16): Input1=[X (l-1) ,u] (16); In the formula, Input l X is the input of the l-th layer; (l-1) Here is the state matrix of the (l-1)th layer; u is the original input variable; The inter-layer connection weights are normalized using the obtained optimal inter-layer scaling factor, as shown in Equation (17): In the formula, This is the normalized inter-layer connection weight matrix; α l Let be the optimal input scaling factor for the l-th layer; ||·|| denotes the Frobenius norm of the matrix; The internal connection weights are normalized using the optimal spectral radius, as shown in Equation (18): In the formula, W (l) ρ is the normalized internal connection weight matrix of the l-th layer; l Let be the optimal spectral radius of the l-th layer; eig(·) represents the matrix eigenvalue operation; max(|·|) represents taking the maximum absolute value; Step 4.6: Calculate the state matrix of the l-th layer ESN and the state update equation, as shown in formula (19): In the formula, X (l) (t) represents the state of the l-th layer ESN at time t; X (l-1) u(t) represents the state of layer l-1 at time t; u(t) represents the input vector at time t; tanh is the hyperbolic tangent activation function; W is the inter-layer connection weight matrix. (l) The internal connection weight matrix; initial state X (l) (0) = 0; The output weights of the l-th layer are calculated as shown in formula (20): In the formula, This is the output weight matrix for the l-th layer; Let β be the augmented state matrix of the l-th layer; l Y is the optimal regularization coefficient for the l-th layer; I is the identity matrix; train Labels for the training set; This indicates that the augmented state matrix of the l-th layer consists of the l-th layer state matrix, the original input, and the bias term.

7. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 5 is as follows: Step 5.1: Calculate the prediction output of each layer of ESN. The prediction result of the l-th layer is shown in formula (21): In the formula, This represents the prediction result for the L-th layer ESN; Let l be the augmented state matrix of the l-th layer; This is the output weight matrix for the l-th layer; Step 5.2: The prediction results of each layer are fused using an average ensemble strategy, and the prediction results of each layer are directly averaged, as shown in formula (22): In the formula, This represents the average ensemble prediction result; NL represents the total number of layers in the ESN. This represents the prediction result for the L-th layer ESN; Step 5.3: Average the squared prediction results of each layer, as shown in formula (23): In the formula, The result is the ensemble prediction average after squared calculation; NL is the total number of layers in the ESN. This represents the prediction result for the L-th layer ESN; Step 5.4: Select the better-performing integration method as the final soft measurement prediction model.

8. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 1, characterized in that, The specific process of step 6 is as follows: Step 6.1: Collect test data and perform the same normalization and smoothing processes as the training data to obtain the normalized test set U. test ; Step 6.2: Calculate the state matrix of each layer of the test set using the trained multilayer ESN model. and augmented state matrix Step 6.3: Calculate the test set prediction results using the trained output weights of each layer, as shown in formula (24): In the formula, The prediction results of the layer l ESN on the test set; Let l be the augmented state matrix of the l-th layer on the test set; The output weight matrix trained for the l-th layer; Step 6.4: Integrate the test prediction results of each layer to obtain the final prediction output, as shown in formula (25): In the formula, The results are the ensemble predictions on the test set; NL represents the total number of ESN layers. The prediction results of the layer l ESN on the test set; Step 6.5: Perform an inverse transformation on the prediction results to obtain the actual predicted value of the primary network return water temperature, as shown in formula (26): In the formula, Y pred This is the actual predicted value; This represents the integrated prediction results in the normalized space.

9. The heat load modeling method for heat exchange stations using an integrated echo state network according to claim 4, characterized in that, The root mean square error is shown in formula (27): In the formula, MSE is the root mean square error; N is the sample size; y i This represents the true value of the i-th sample. Let be the predicted value for the i-th sample; The coefficient of determination is shown in formula (28): In the formula, R 2 The coefficient of determination; SS res Sum of squared residuals; SS tot The total sum of squares; in In the formula, SS res y is the sum of squared residuals; i This represents the true value of the i-th sample. Let be the predicted value of the i-th sample; N is the number of samples; In the formula, SS tot For the total sum of squares; y i This represents the true value of the i-th sample. The mean of the true values ​​is N; the sample size is N. The mean absolute error is shown in formula (31): In the formula, MAE is the mean absolute error; N is the sample size; y i This represents the true value of the i-th sample. Let be the predicted value for the i-th sample.

Images