Industrial process operation index modeling method for cross-time scale sampling

Abstract

The invention discloses a cross-time scale sampling-oriented industrial process operation index modeling method, which comprises the following steps of: acquiring process variables and operation indexes of an industrial process under different sampling frequencies, and performing preliminary preprocessing to obtain process data D *; a maximum information coefficient time lag technology based on disc discretization processing is utilized to identify time lag in the process data D * after preliminary preprocessing and align the time lag, and process data D after overall data cleaning is obtained; dividing the process data D after overall data cleaning into a process variable set X and an operation index set Y, and establishing an industrial process operation index soft measurement model for cross-time scale sampling. According to the method, the consistency of input and output quantities of the model can be ensured when the model faces multi-sampling-period data, meanwhile, geometric structure information of a process data sample without operation indexes is fully utilized to assist model construction, and the risk of model over-fitting is reduced.

Description

Technical Field

[0001] The present invention belongs to the technical field of industrial process multi-sampling period data modeling, and specifically relates to an industrial process operation indicator modeling method for cross-time scale sampling. Background Art

[0002] With the rapid development of industrial technology and the deepening of digital transformation, the data involved in industrial processes have shown unprecedented complexity and multi-periodicity. On the one hand, the industrial process is a continuous process with multiple processes and multiple equipment working in coordination. Along with the processing of equipment and the flow of materials, the industrial process data contains noise from the comprehensive effects of the surrounding environment, the operation of processing equipment, etc. In view of the problem of noise in process data, the median filter is used to reduce the impact of noise by local sorting and taking the median, thereby achieving data smoothing and noise reduction. There is a lag in the time from the production operation of each process to the final product, which leads to a time delay between the collected process data, that is, there is an erroneous time series correspondence between the process data. If not handled in time, the accuracy of process data modeling will be greatly reduced. On the other hand, process variables and operating indicators are collected by sensors and manual sampling and testing respectively, so that the collected data presents multi-sampling period characteristics. In practical problems, multi-sampling period data is manifested as the simultaneous existence of process data with operating indicators and process data without operating indicators in the data. The data input model will cause the model input and output to mismatch. Most of the existing modeling methods only use process data samples with operating indicators for modeling, ignoring the potential information contained in a large number of process data samples without operating indicators, which greatly reduces the utilization rate of sampled data. Summary of the invention

[0003] In view of the problems in the prior art that the model input and output do not match when facing multi-sampling period data and the potential information of a large number of process data samples without operating indicators is ignored, the present invention provides an industrial process operation indicator modeling method for cross-time scale sampling, which can ensure the consistency of the input and output quantities of the model when facing multi-sampling period data, and at the same time make full use of the geometric structure information of the process data samples without operating indicators to assist model construction, thereby reducing the risk of model overfitting.

[0004] The present invention discloses a method for industrial process operation index modeling for cross-time scale sampling, the method comprising the following steps:

[0005] S1, collect process variables and operating indicators of industrial processes at different sampling frequencies, and obtain process data D after preliminary preprocessing * ;

[0006] S2, using the maximum information coefficient time-lag technique based on discretization processing, identifies the process data D after preliminary preprocessing *The time lag in is aligned to obtain the process data D after overall data cleaning;

[0007] S3, the process data D after overall data cleaning is divided into a process variable set X and an operation indicator set Y, and a soft sensor model of industrial process operation indicators for cross-time scale sampling is established.

[0008] Step S1 further comprises:

[0009] S11, obtain N process variables and 1 operating index at different sampling frequencies through sensors and manual sampling, and obtain the de-noised process data D after median filtering. ** ;

[0010] S12, using the formula For the filtered process data D ** Perform normalization processing to obtain normalized process data D * , where μ D , σ D denote mean and variance respectively;

[0011] Step S2 further comprises:

[0012] S21, select the smallest diameter R among n samples min , maximum diameter R max , minimum polar angle θ min and the maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, respectively, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,θ max ];

[0013] S22, in the defined polar range [R min ,R max ] is divided into a equal parts, and [θ min ,θ max ]It is divided into b equal parts, and a*b grids are obtained after division;

[0014] S23, for each set of a and b values, calculate each sample of the i-th process variable under the corresponding grid division and the corresponding running indicator variable samples The mutual information value of The sum of all sample mutual information values ​​of the i-th process variable

[0015] S24, normalize the calculated mutual information value, record the division situation corresponding to the maximum mutual information value, and calculate the MIC value under this division situation, the formula is as follows:

[0016]

[0017] In the formula, log 2 (min{a,b}) represents the discretization area of ​​the disk in the polar coordinate system [R min ,R max ] and [θ min ,θ max ] is divided into a*b grids, and the normalized value range is [0,1];

[0018] S25, run the indicator data y * In the time domain relative to the process variable Gradually translate forward, the translation step length is The sampling period is , and the i-th process variable is calculated once every translation step and operation index y * and record the MIC value each time; select the corresponding translation step number t when the MIC value reaches the maximum value, and determine it as the i-th process variable and operation index y * Optimal delay time t;

[0019] S26, will run indicator y * In the time domain, the ith process variable x is obtained by translating forward t time steps in the time domain. i and running indicator y.

[0020] Step S21 further includes:

[0021] S211, in the plane rectangular coordinate system, it is known that the process variable value and the operating index value corresponding to the jth sample of the i-th process variable are and The corresponding rectangular coordinates are The value range of i is an integer in [1, N], and N is the normalized process data D * The total number of process variables in the process, the value range of j is an integer in [1,n], and n is the total number of process data samples;

[0022] S212, calculate the rectangular coordinates according to the polar coordinate transformation formula The corresponding polar coordinates (R ij ,θ ij ), where R ij and θ ijdenote the polar diameter and polar angle of the jth sample of the ith process variable respectively;

[0023] S213, select the smallest diameter R among n samples min , maximum diameter R max , minimum polar angle θ min and the maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,] max ].

[0024] Step S22 further includes:

[0025] In the delineated polar axis range [R min ,R max ] is divided into a equal parts, and a+1 points on the polar axis are represented by r 0 ,r 1 ,...,r a , according to the a+1 points in the defined polar angle range [θ min ,θ max ] to draw an arc counterclockwise; at the same time, [θ min ,θ max ] is divided into b equal parts, and [θ min ,θ max ], with the point as the center and a length of [R min ,R max ], denoted as ω 0 ,ω 1 ,...,ω b ; Finally, we get a*b grids after division; where the value range of a and b is (0,c a ) and (0,c b ), c a and c b They are the upper limits of the values ​​of a and b, and each time the division is performed, the values ​​of a and b satisfy a*b<n 0.6 .

[0026] Step S3 further comprises:

[0027] S31, prepare multi-sampling period data: according to whether the process variable has an operating indicator, all sample process variable sets of the training set X={x 1 ,x 2 ,...,x N}, x i =[x 1,i ,x 2,i ,...,x d,i ]T ∈R d It is divided into a process data sample set with operation indicators and a process data sample set without operation indicators; the process data sample set with operation indicators is The corresponding operating indicators are The process data sample set without operation index is Where N A is the total number of samples, N U is the number of samples without running indicators, N M is the number of samples with running indicators, satisfying N A =N U +N M ;

[0028] S32, construct cross-scale objective function and supervision mechanism: for multi-sampling period data, add cross-time scale constraints to the incremental network, obtain a new objective function, and solve the output weight β L Update the formula and then build a node update supervision mechanism;

[0029] S33, model parameter initialization: set the maximum number of hidden nodes of the established model to L max , the expected error is ∈ tol , the maximum number of random configurations is T max , the parameter λ candidate set ψ is: {λ min ,...,λ max}, the distribution interval of random parameters is ±λ, and the candidate set φ of parameter r is: {r min ,...,r max};L 2 The regularization coefficient is C, the cross-time scale constraint coefficient is R, and the time step distance constraint parameter is U;

[0030] S34, randomly generate candidate parameters and calculate hidden layer output: for each newly added node L, from [-λ,λ] and [-λ,λ] d Randomly generate T max Group candidate weight w L and bias b L , and calculate the hidden output corresponding to each group as well as where g(·) is the sigmoid function, is the hidden layer output with running indicator samples, H is the hidden layer output of the sample without running indicators, L is the hidden layer output corresponding to all samples;

[0031] S35, select the best candidate parameters: calculate Select the largest one in the candidate hidden layer node pool The corresponding candidate hidden layer node parameters As the optimal candidate hidden layer node parameters, using the selected optimal candidate hidden node parameters Calculate the corresponding hidden Update the hidden layer mapping matrix Add the optimal candidate hidden node to the current network model;

[0032] S36. Calculate the model residual after adding the new node: Calculate the residual modulus value of the optimal candidate hidden layer node: In the formula, represents the residual when the current number of hidden layer nodes is L, represents the residual when the network has L - 1 hidden nodes, represents the output matrix of the current network, β L is the output weight of the network;

[0033] S37. Judge whether the model meets the expectation after adding the new node: When is less than the given expectation tolerance the model establishment process ends; otherwise, return to step S34 and add nodes to the network in the form of dot increment, that is, L = L + 1, until L reaches the maximum number of hidden layer nodes L max up to.

[0034] Furthermore, in step S32, the constructed objective function expression is as follows:

[0035]

[0036] In the formula, the first term is the empirical error term of the process data sample with the operation index, the second term is L 2 regularization, the third term is the cross - time - scale constraint term, f(.) represents the output of the network, ||·|| represents the modulus in the form of the second norm, h i,j is the cross - time - scale constraint term coefficient, obtained from the distance between the process data sample set X M with the operation index and the process data sample set X M without the operation index in the manifold space. The specific calculation formula is as follows:

[0037]

[0038] In the formula, dist(x i , x j ) represents the Euclidean distance between the samples x i and x j ;

[0039] Furthermore, in step S32, the matrix - form expression of the constructed objective function is as follows:

[0040]

[0041] In the formula, Indicates the hidden layer output corresponding to the running indicator process data sample, represents the hidden layer output corresponding to the process data sample without operation index, β L represents the output weight of the Lth node, ||·|| F represents the Frobenius norm modulus, O, P, Q represent the cross-time scale correlation matrix, and the calculation formula is as follows:

[0042]

[0043] Furthermore, the following formula is used to construct the output weight update formula of the Lth hidden node:

[0044]

[0045] Where I is the unit matrix,

[0046] Furthermore, in step S3, combined with the idea of ​​error convergence, it is determined whether the output of each sample hidden layer of the Lth newly added node satisfies the supervision constraint. The specific supervision mechanism expression is as follows:

[0047]

[0048] In the formula, k = 1, 2, ..., m, m represents the dimension corresponding to the operating index in the multi-sampling period data, <·, ·> represents the inner product of the vector, and the superscript "T" represents the transposition operation. represents the supervision constraint corresponding to the kth operating indicator sample when the number of current hidden layer nodes is L, μ L =(1-r) / (L+1), It represents the error corresponding to the kth running index sample when the number of current hidden layer nodes is L-1; if T max The randomly generated weights and biases of the group cannot meet the supervision constraints, and the increase of μ L , r=r+Δr, Δr is the incremental parameter of the parameter r allocation interval; if all values ​​in the parameter r allocation interval are traversed, and the randomly generated weights and biases still cannot meet the supervision constraints, then the candidate parameter λ allocation interval is expanded, λ=λ+Δλ, Δλ is the incremental parameter of the allocation interval.

[0049] The beneficial effects of the present invention are:

[0050] First, to address the problem of time lags between data, the industrial process operation indicator modeling method for cross-time scale sampling of the present invention, based on the maximum information coefficient time lag identification technology of disk discretization processing, identifies and aligns the time lags in industrial process data, and corrects the erroneous time series correspondence between process data.

[0051] Second, in order to solve the problem that multi-sampling periodic data leads to model input-output mismatch and potential information of process data samples without operation indicators is difficult to utilize, the industrial process operation indicator modeling method for cross-time scale sampling of the present invention designs cross-time scale constraints based on the manifold hypothesis, constructs a manifold structure diagram of process data samples, mines the correlation between samples in the manifold structure, and constrains process data samples without operation indicators to have similar or identical outputs with the nearest process data samples with operation indicators. This not only solves the problem that multi-sampling periodic data leads to model input-output mismatch, but also makes full use of the geometric structure information in process data samples without operation indicators, thereby improving the accuracy of modeling of multi-sampling periodic data of industrial processes.

[0052] Third, the industrial process operation indicator modeling method for cross-time scale sampling of the present invention can ensure the consistency of the input and output quantities of the model when facing multi-sampling period data, make full use of the geometric structure information of process data samples without operation indicators to assist model construction, and reduce the risk of model overfitting. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 It is a schematic diagram of the overall process of the present invention;

[0054] Figure 2 It is a schematic diagram of multi-sampling period data in the present invention;

[0055] Figure 3 It is a schematic diagram of interval determination based on disc discretization processing in the present invention;

[0056] Figure 4 This is a flow chart for modeling industrial process operation indicators for sampling across time scales in the present invention. DETAILED DESCRIPTION

[0057] The following examples will enable those skilled in the art to more fully understand the present invention, but are not intended to limit the present invention in any way.

[0058] The overall process diagram of the present invention is as follows Figure 1As shown in the figure, industrial process data are collected by sensors and manual sampling and transmitted to the data storage terminal. Then, the median filtering and normalization techniques are used to preliminarily preprocess the process data. The maximum information coefficient time-delay identification technology based on disk discretization is used to identify and align the time lags in the process data after preliminary preprocessing. Finally, a soft measurement model of industrial process operation indicators for cross-time scale sampling is established to train and test multi-sampling period data of industrial processes.

[0059] The present invention discloses a method for modeling industrial process operation indicators for sampling across time scales, comprising the following steps:

[0060] Step 1: Industrial process multi-sampling cycle data collection and preliminary preprocessing: Use sensors and manual sampling to collect industrial process multi-sampling cycle data, and combine median filtering and normalization to process it to obtain the process data D after preliminary preprocessing. * ;

[0061] Step 2: Time lag identification and alignment of process data after preliminary preprocessing: Using the maximum information coefficient time lag technology based on discretization, the process data D after preliminary preprocessing is identified. * The time lag in is aligned to obtain the process data D after overall data cleaning;

[0062] Step 3: Model establishment: Divide the data set D after overall data cleaning into a process variable set X and an operation indicator set Y, and establish a soft measurement model of industrial process operation indicators for cross-time scale sampling.

[0063] In step 1, industrial process multi-sampling period data collection and preliminary preprocessing include the following steps:

[0064] 1.1 The process variables and operating indicators at different sampling frequencies are obtained through sensors and manual sampling and then processed by median filtering to obtain the noise-reduced process data D ** , which includes N process variables and 1 operating indicator;

[0065] 1.2 Use For the filtered process data D ** Perform normalization processing to obtain normalized process data D * , where μ D , σ D represent the mean and variance respectively.

[0066] In step 2, the maximum information coefficient time lag identification technology based on disc discretization is used to identify and align the time lag in the process data after preliminary preprocessing (only process data samples with operating indicators are considered when calculating the time lag), and then the process data set D after overall data cleaning is obtained. The specific steps of maximum information coefficient time lag identification and alignment based on disc discretization are as follows:

[0067] 2.1 Calculate the discretization range: First, in the plane rectangular coordinate system, the process variable value and operating index value corresponding to the jth sample of the i-th process variable are known to be and The corresponding rectangular coordinates are The value range of i is an integer in [1, N], and N is the normalized process data D * The total number of process variables in the equation, j is an integer in the range [1, n], and n is the total number of process data samples; secondly, according to the polar coordinate transformation formula, the rectangular coordinates are calculated. The corresponding polar coordinates (R ij ,θ ij ), where R ij and θ ij Respectively represent the polar diameter and polar angle of the jth sample of the ith process variable; select the smallest polar diameter R among n samples min , maximum diameter R max , minimum polar angle θ min , maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, respectively, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,θ max ];

[0068] 2.2 Discretization range grid processing: In the defined polar axis range [R min ,R max ] is divided into a equal parts, and a+1 points on the polar axis are represented by r 0 ,r 1 ,...,r a , according to the a+1 points in the defined polar angle range [θ min ,θ max ] to draw an arc counterclockwise. At the same time, [θ min ,θ max ] is divided into b equal parts, and [θ min ,θ max ], with the point as the center and a length of [R min ,R max ], denoted as ω 0 ,ω1 ,...,ω b Finally, we get a*b grids after division. The value range of a and b is (0,c a ) and (0,c b ), c a and c b They are the upper limits of the values ​​of a and b, and each time a and b are divided, they must satisfy a*b<n 0.6 ;

[0069] 2.3 Calculate the mutual information value under different partitions: For each set of a and b values, calculate each sample of the i-th process variable under the grid partition and the corresponding running indicator variable samples The mutual information value of The sum of all sample mutual information values ​​of the i-th process variable

[0070] 2.4 Calculate the maximum information coefficient value: Normalize the calculated mutual information value, record the division situation corresponding to the maximum mutual information value, and calculate the MIC value under this division situation. The formula is as follows:

[0071]

[0072] In the formula, log 2 (min{a,b}) represents the discretization area of ​​the disk in the polar coordinate system [R min ,R max ] and [θ min ,θ max ] is divided into a*b grids. The normalized value range is [0,1].

[0073] 2.5 Lag Estimation: Run the indicator data y * In the time domain relative to the process variable Gradually translate forward, the translation step length is The sampling period is , and the i-th process variable is calculated once every translation step and operation index y * The MIC value between the two values ​​is recorded. The corresponding translation step number t when the MIC value reaches the maximum value is selected and determined as the i-th process variable and operation index y * Optimal delay time t;

[0074] 2.6 Time-delay alignment: Align the running indicator y * In the time domain, the ith process variable x is obtained by translating forward t time steps in the time domain. i and running indicator y.

[0075] In step 3, model building includes the following steps:

[0076] 3.1 Prepare multi-sampling period data: According to whether the process variable has an operating indicator, all sample process variable sets of the training set X = {x 1 ,x 2 ,...,x N}, x i =[x 1,i ,x 2,i ,...,x d,i ] T ∈R d Divided into process data sample sets with operation indicators and process data sample sets without operation indicators. Process data sample sets with operation indicators: The corresponding operating indicators are: Sample set of process data without operating indicators: Where N A is the total number of samples N U , is the number of samples without running indicators, N M is the number of samples with running indicators, satisfying N A =N U +N M .

[0077] 3.2 Construction of cross-scale objective function and supervision mechanism: For multi-sampling period data, add cross-time scale constraints to the incremental network to obtain a new objective function and solve the output weight β L Update the formula and then build a node update supervision mechanism.

[0078] Preferably, the constructed objective function expression is as follows:

[0079]

[0080] In the formula, the first term is the empirical error term of the process data sample with operation indicators, and the second term is L 2 Regularization, the third term is the cross-time scale constraint term, f(.) represents the output of the network, ||·|| represents the modulus in the form of the two norm, h i,j is the coefficient of the cross-time scale constraint term, that is, the similarity matrix based on the manifold space, which can be obtained by the process data sample set X with operation indicators. M And the process data sample set X without operation index M The distance in the manifold space is obtained, and the specific calculation formula is as follows:

[0081]

[0082] In the formula, dist(x i ,x j ) represents the sample xi and x j The Euclidean distance between them is the Gaussian weighted Euclidean distance formula.

[0083] In particular, the matrix form of the cross-time scale constraint terms is constructed, and after converting the remaining terms into matrix form, the objective function matrix is ​​expressed as follows:

[0084]

[0085] In the formula, Indicates the hidden layer output corresponding to the running indicator process data sample, represents the hidden layer output corresponding to the process data sample without operation index, β L represents the output weight of the Lth node, ||·|| F represents the Frobenius norm modulus, O, P, Q represent the cross-time scale correlation matrix, and the calculation formula is as follows:

[0086]

[0087] Preferably, the expression for constructing the output weight update formula of the Lth hidden node is as follows:

[0088]

[0089] Where I is the unit matrix,

[0090] 3.3 Model parameter initialization: Set the maximum number of hidden nodes of the established model to L max , the expected error is ∈ tol , the maximum number of random configurations is T max , the parameter λ candidate set ψ is: {λ min ,...,λ max}, the distribution interval of random parameters is ±λ, and the candidate set φ of parameter r is: {r min ,...,r max};L 2 The regularization coefficient is C, the cross-time scale constraint coefficient is R, and the time step distance constraint parameter is U;

[0091] 3.4 Randomly generate candidate parameters and calculate hidden layer output: For each newly added node L, from [-λ,λ] and [-λ,λ] d Randomly generate T max Group candidate weight w L and bias b L , and calculate the hidden output corresponding to each group as well as where g(·) is the sigmoid function, is the hidden layer output with running indicator samples, H is the hidden layer output of the sample without running indicators, L is the hidden layer output corresponding to all samples;

[0092] 3.5 Select the optimal candidate parameters: Calculate the maximum value among the candidate hidden layer node pool Corresponding candidate hidden layer node parameters As the best candidate hidden layer node parameters, using the selected best candidate hidden node parameters Calculate the corresponding implicit Update the hidden layer mapping matrix Finally, add the best candidate hidden node to the current network model;

[0093] 3.6 Calculate the model residual after adding new nodes: Calculate the residual modulus value of the best candidate hidden layer node: In the formula, Represents the residual when the number of hidden layer nodes is L. represents the residual when the network has L-1 hidden nodes, represents the output matrix of the current network, β L is the output weight of the network;

[0094] 3.7 Determine whether the model expectations are met after adding new nodes: Less than a given desired tolerance ∈ tol When , the model building process ends, otherwise return to step 3.4 and add nodes to the network in the form of point increments, that is, L = L + 1, until L reaches the maximum number of hidden layer nodes L max until.

[0095] Preferably, combined with the idea of ​​error convergence, the supervision constraints that each newly added node needs to satisfy can be obtained. The specific supervision mechanism expression is as follows:

[0096]

[0097] In the formula, k = 1, 2, ..., m, m represents the dimension corresponding to the operating index in the multi-sampling period data, <·, ·> represents the inner product of the vector, and the superscript "T" represents the transposition operation. represents the supervision constraint corresponding to the kth operating indicator sample when the number of current hidden layer nodes is L, μ L =(1-r) / (L+1), It represents the error corresponding to the kth running index sample when the number of current hidden layer nodes is L-1; if T max The randomly generated weights and biases of the group cannot meet the supervision constraints, and the increase of μ L, that is, r = r + Δr, Δr is the incremental parameter of the parameter r allocation interval; if all values ​​in the parameter r allocation interval are traversed, and the randomly generated weights and biases still cannot meet the supervision constraints, then the candidate parameter λ allocation interval is expanded, that is, λ = λ + Δλ, Δλ is the incremental parameter of the allocation interval.

[0098] Examples

[0099] The present invention provides an industrial process operation index modeling method for cross-time scale sampling, and the specific steps are as follows:

[0100] Step 1: Industrial process data collection and preliminary preprocessing: Use sensors and manual sampling to collect industrial process data, including 4 process variables and 1 operating indicator, a total of 12013 sets of data, where the process variable set X = {x 1 ,x 2 ,x 3 ,x 4 The ratio of the acquisition frequency of the variable Y to the operating indicator variable is 20:1. The median filter is used to reduce the noise of the collected process data under the four sampling ratios. For the filtered process data set D ** Normalization is performed in to obtain the normalized process data set D * , where μ D , σ D Respectively represent the mean and variance; the present invention uses multi-input single-output multi-sampling period data such as Figure 2 As shown, within the same time interval, the process variable is a fast sampling period variable under a fast sampling frequency, and the operating indicator variable is a slow sampling period variable under a slow sampling frequency. The difference in the sampling periods of the two results in the existence of both process data samples with operating indicators and process data samples without operating indicators in the collected process data.

[0101] Step 2, identification and alignment of time delays in process data after preliminary preprocessing: using the maximum information coefficient time delay identification technology based on disc discretization processing, identify and align the time delays in the process data set after preliminary preprocessing (only consider process data samples with operating indicators when calculating the time delay), and then obtain the process data set D after the overall data cleaning. The present invention provides a maximum information coefficient time delay identification technology based on disc discretization processing, which is used for time delay identification and alignment in process data. In this technology, the range of data discretization processing is determined based on the disc, as shown in the following example. Figure 3 The specific steps of maximum information coefficient time lag identification and alignment based on disk discretization are as follows:

[0102] 2.1 Calculate the discretization range: First, in the plane rectangular coordinate system, the process variable value and operating index value corresponding to the jth sample of the i-th process variable are known to be and The corresponding rectangular coordinates are The value range of i is an integer in [1,4], and N is the normalized process data D * The total number of process variables in the equation, j is an integer in the range of [1,12013], and 12013 is the total number of process data samples; secondly, according to the polar coordinate transformation formula, the rectangular coordinates are calculated The corresponding polar coordinates (R ij ,θ ij ), where R ij and θ ij Respectively represent the polar diameter and polar angle of the jth sample of the ith process variable; select the smallest polar diameter R among 12013 samples min , maximum diameter R max , minimum polar angle θ min , maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, respectively, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,θ max ];

[0103] 2.2 Discretization range grid processing: Let a=16,b=17 in the defined polar axis range [R min ,R max ] is divided into 16 equal parts, and the 16+1 points on the polar axis are represented by r 0 ,r 1 ,...,r 16 According to these 16+1 points in the defined polar angle range [θ min ,θ max ] to draw an arc counterclockwise. At the same time, [θ min ,θ max ] is divided into 17 equal parts, and [θ min ,θ max ] 17+1, with the dot as the center, and the length is [R min ,R max ], denoted as ω 0 ,ω 1 ,...,ω 17 Finally, we get a 16*17 grid after division. The value range of a and b is (0,c a ) and (0,c b ), c a and cb are the upper limits of the values ​​of a and b respectively;

[0104] 2.3 Calculate the mutual information value under different partitions: Change the values ​​of a and b. For each set of a and b values, calculate each sample of the i-th process variable under the grid partition and the corresponding running indicator variable samples The mutual information value of The sum of all sample mutual information values ​​of the i-th process variable

[0105] 2.4 Calculate the maximum information coefficient value: Normalize the calculated mutual information value, record the division situation corresponding to the maximum mutual information value, and calculate the MIC value under this division situation. The formula is as follows:

[0106]

[0107] In the formula, log 2 (min{a,b}) represents the discretization area of ​​the disk in the polar coordinate system [R min ,R max ] and [θ min ,θ max ] is divided into a*b grids. The normalized value range is [0,1].

[0108] 2.5 Lag Estimation: Run the indicator data y * In the time domain relative to the process variable Gradually translate forward, the translation step length is The sampling period is , and the i-th process variable is calculated once every translation step and operation index y * The MIC value between the two values ​​is recorded. The corresponding translation step number t when the MIC value reaches the maximum value is selected and determined as the i-th process variable and operation index y * Optimal delay time t;

[0109] 2.6 Time-delay alignment: Align the running indicator y * In the time domain, the ith process variable x is obtained by translating forward t time steps in the time domain. i and operation index y;

[0110] Step 3: Model building: Use the data set D after the overall data cleaning of the industrial process operation indicator soft measurement model for cross-time scale sampling to build a model, see Figure 4 , the specific modeling steps are as follows:

[0111] 3.1 Divide the training set and test set: Divide the data set D after the overall data cleaning into a training set and a test set. The training set is the first 70% of the overall data set, a total of 8409 sets of data, and the test set is the first 3000 sets of data in the last 30% of the data set;

[0112] 3.2 Prepare multi-sampling period data: According to whether the process variable has an operating indicator, all sample process variables in the training set X = {x 1 ,x 2 ,...,x 12013}, x i =[x 1,i ,x 2,i ,x 3,i ,x 4,i ] T ∈R 4 Divided into process data sample sets with operation indicators and process data sample sets without operation indicators. Process data sample sets with operation indicators: The corresponding operating indicators are: Sample set of process data without operating indicators: Total number of samples N A 8409, no running indicator sample number N U is 421, and the number of running index samples is N M is 7988, satisfying N A =N U +N M ; Do the same for the test set;

[0113] 3.3 Calculate the similarity matrix based on manifold space: Calculate the sample set X of the process data of the operating index M And the process data sample set X without operation index M The distance in the manifold space is calculated as follows:

[0114]

[0115] In the above formula, dist(x i ,x j ) represents the sample x i and x j The Euclidean distance between them is the Gaussian weighted Euclidean distance formula.

[0116] 3.4 Calculation of cross-time scale correlation matrix: The cross-time scale correlation matrix O, P, Q is expressed as follows:

[0117]

[0118] 3.5 Model parameter initialization: Set the maximum number of hidden nodes of the established model to 10, the expected error to 0.01, the maximum number of random configurations to 20, the parameter λ candidate set ψ to: {1:5:100}, the random parameter allocation interval to ±λ, and the parameter r candidate set φ to: {0.1:0.1:0.9}; L 2 The regularization coefficient C is 0.9, the cross-scale constraint coefficient R and the time step distance constraint parameter U are in the process variable set X = {x 1 ,x 2 ,x 3 ,x 4} and the operating indicator variable Y when the ratio of collection frequency is 20:1 is 8 and 64 respectively;

[0119] 3.6 Randomly generate candidate parameters and calculate hidden layer output: For each newly added node L, from [-λ,λ] and [-λ,λ] d Randomly generate 20 sets of candidate weights w corresponding to the Lth hidden node L and bias b L , and calculate the hidden output corresponding to each group as well as where g(·) is the sigmoid function, is the hidden layer output with running indicator samples, H is the hidden layer output of the sample without running indicators, L is the hidden layer output corresponding to all samples;

[0120] 3.7 Determine whether the supervision mechanism is satisfied: Determine whether the hidden layer output of each sample of the Lth newly added node satisfies the supervision constraint. The specific supervision mechanism expression is as follows:

[0121]

[0122] Where k = 1, 2, ..., 4, k represents the dimension corresponding to the operating index in the multi-sampling period data, <·, ·> represents the inner product of the vector, and the superscript "T" represents the transposition operation; represents the supervision constraint corresponding to the kth running index sample of each training set when the number of current hidden layer nodes is L, μ L =(1-r) / (L+1), represents the error corresponding to the kth operating indicator sample when the number of current hidden layer nodes is L-1; if the 20 sets of randomly generated weights and biases cannot meet the supervision constraints, expand μ LThe calculation space of r, that is, r = r + Δr, Δr is the incremental parameter of the parameter r allocation interval; if all values ​​in the parameter r allocation interval are traversed, and the randomly generated weights and biases still cannot meet the supervision constraints, then the candidate parameter λ allocation interval is expanded, that is, λ = λ + Δλ, Δλ is the incremental parameter of the allocation interval;

[0123] 3.8 Selecting the best candidate parameters: Calculation Select the largest one in the candidate hidden layer node pool Corresponding candidate hidden layer node parameters As the best candidate hidden layer node parameters, using the selected best candidate hidden node parameters Calculate the corresponding implicit Update the hidden layer mapping matrix Finally, the best candidate hidden node is added to the current network 3.9 to calculate the output weight β L : Calculate the output weight of the Lth hidden node. The expression is as follows:

[0124]

[0125] 3.10 Calculate the model residual after adding new nodes: Calculate the residual modulus value of the best candidate hidden layer node: In the formula, Represents the residual when the number of hidden layer nodes is L. represents the residual when the network has L-1 hidden nodes, Represents the output matrix of the current network;

[0126] 3.11 Determine whether the model expectations are met after adding new nodes: When the expected tolerance is less than 0.01, the model building process ends. Otherwise, it returns to 3.6 and adds nodes to the network in the form of point increments, that is, L = L + 1. This modeling process is repeated until It is greater than the given expected tolerance 0.01 or L reaches the maximum number of hidden layer nodes 10.

[0127] Step 4: Model testing: Divide the data set into a training set and a test set, and use the test set to test the accuracy of the trained cross-time scale neural network operation indicator model.

[0128] In order to fully illustrate the performance of this model in modeling industrial multi-sampling cycle data, the model method is compared with a shallow network model without cross-time scale constraints and a learning model that can only constrain the nearest neighbor samples without operating indicators in terms of root mean square error, mean absolute error, and the best value in 20 running results. The "mean value ± standard deviation" of 20 independent repeated experiments is taken, and the test performance results of different models are shown in the following table. It can be seen from the table that the model method proposed in the present invention is superior to the other two models in terms of test accuracy on industrial process multi-sampling cycle data.

[0129] Table 1 Test results of 3 models when the sampling frequency ratio of process variables and operating indicators is 20:1

[0130]

[0131] The present invention establishes a soft measurement model of industrial process operation indicators for sampling across time scales, which can ensure the consistency of the input and output quantities of the model when facing multi-sampling period data, and at the same time make full use of the geometric structure information of process data samples without operation indicators to assist model construction, making up for the shortcomings of traditional modeling methods in terms of inconsistent input and output when facing multi-sampling period data and unutilized potential information of process data samples without operation indicators. The model is particularly suitable for the field of multi-sampling period data modeling of industrial processes.

[0132] The above are only preferred embodiments of the present invention. The protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.

Claims

1. A method for modeling industrial process operation indicators for sampling across time scales, characterized in that: The method comprises the following steps: S1, collect process variables and operating indicators of industrial processes at different sampling frequencies, and obtain process data D after preliminary preprocessing * ; S2, using the maximum information coefficient time-lag technique based on discretization processing, identifies the process data D after preliminary preprocessing * The time lag in is aligned to obtain the process data D after overall data cleaning; S3, the process data D after overall data cleaning is divided into a process variable set X and an operation indicator set Y, and a soft sensor model of industrial process operation indicators for cross-time scale sampling is established.

2. The industrial process operation index modeling method for cross-time scale sampling according to claim 1 is characterized in that: Step S1 further comprises: S11, obtain N process variables and 1 operating index at different sampling frequencies through sensors and manual sampling, and obtain the de-noised process data D after median filtering. ** ; S12, using the formula For the filtered process data D ** Perform normalization processing to obtain normalized process data D * , where μ D , σ D represent mean and variance respectively.

3. The industrial process operation index modeling method for cross-time scale sampling according to claim 1 is characterized in that: Step S2 further comprises: S21, select the smallest diameter R among n samples min , maximum diameter R max , minimum polar angle θ min and the maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, respectively, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,θ max ]; S22, in the defined polar range [R min ,R max ] is divided into a equal parts, and [θ min ,θ max ]It is divided into b equal parts, and a*b grids are obtained after division; S23, for each set of a and b values, calculate each sample of the i-th process variable under the corresponding grid division and the corresponding running indicator variable samples The mutual information value of The sum of all sample mutual information values ​​of the i-th process variable S24, normalize the calculated mutual information value, record the division situation corresponding to the maximum mutual information value, and calculate the MIC value under this division situation, the formula is as follows: Where log2(min{a,b}) represents the discretization area of ​​the disk in the polar coordinate system [R min ,R max ] and [θ min ,θ max ] is divided into a*b grids, and the normalized value range is [0,1]; S25, run the indicator data y * In the time domain relative to the process variable Gradually translate forward, the translation step length is The sampling period is , and the i-th process variable is calculated once every translation step and operation index y * and record the MIC value each time; select the corresponding translation step number t when the MIC value reaches the maximum value, and determine it as the i-th process variable and operation index y * Optimal delay time t; S26, will run indicator y * In the time domain, the ith process variable x is obtained by translating forward t time steps in the time domain. i and running indicator y.

4. The industrial process operation index modeling method for cross-time scale sampling according to claim 3 is characterized in that: Step S21 further comprises: S211, in the plane rectangular coordinate system, it is known that the process variable value and the operating index value corresponding to the jth sample of the i-th process variable are and The corresponding rectangular coordinates are The value range of i is an integer in [1, N], and N is the normalized process data D * The total number of process variables in the process, the value range of j is an integer in [1,n], and n is the total number of process data samples; S212, calculate the rectangular coordinates according to the polar coordinate transformation formula The corresponding polar coordinates (R ij ,θ ij ), where R ij and θ ij denote the polar diameter and polar angle of the jth sample of the ith process variable respectively; S213, select the smallest diameter R among n samples min , maximum diameter R max , minimum polar angle θ min and the maximum polar angle θ max As the boundary values ​​of the discretization range on the polar radius and polar angle, the range on the polar axis is defined [R min ,R max ] and the range of polar angles [θ min ,θ max ].

5. The industrial process operation index modeling method for cross-time scale sampling according to claim 3 is characterized in that: Step S22 further includes: In the delineated polar axis range [R min ,R max ] is divided into a equal parts, and a+1 points on the polar axis are represented as r0, r1, ..., r a , according to the a+1 points in the defined polar angle range [θ min ,θ max ] to draw an arc counterclockwise; at the same time, [θ min ,θ max ] is divided into b equal parts, and [θ min ,θ max ], with the point as the center and a length of [R min ,R max ], denoted as ω0,ω1,...,ω b ; Finally, we get a*b grids after division; where the value range of a and b is (0,c a ) and (0,c b ), c a and c b They are the upper limits of the values ​​of a and b, and each time the division is performed, the values ​​of a and b satisfy a*b<n 0.6 .

6. The industrial process operation index modeling method for cross-time scale sampling according to claim 1 is characterized in that: Step S3 further comprises: S31, prepare multi-sampling period data: according to whether the process variable has an operating indicator, all sample process variable sets of the training set X={x1,x2,...,x N }, x i =[x 1,i ,x 2,i ,...,x d,i ] T ∈R d It is divided into a process data sample set with operation indicators and a process data sample set without operation indicators; the process data sample set with operation indicators is The corresponding operating indicators are The process data sample set without operation index is Where N A is the total number of samples, N U is the number of samples without running indicators, N M is the number of samples with running indicators, satisfying N A =N U +N M ; S32, construct cross-scale objective function and supervision mechanism: for multi-sampling period data, add cross-time scale constraints to the incremental network, obtain a new objective function, and solve the output weight β L Update the formula and then build a node update supervision mechanism; S33, model parameter initialization: set the maximum number of hidden nodes of the established model to L max The expected error is The maximum number of random configurations is T max 、The candidate set ψ of parameter λ is: {λ min ,...,λ max }, the distribution interval of random parameters is ±λ, and the candidate set φ of parameter r is: {r min ,...,r max }; The L2 regularization coefficient is C, the cross-time scale constraint coefficient is R, and the time step distance constraint parameter is U; S34, randomly generate candidate parameters and calculate hidden layer output: for each newly added node L, from [-λ,,λ] and [-λ,,λ] d Randomly generate T max Group candidate weight w L and bias b L , and calculate the hidden output corresponding to each group as well as where g(·) is the sigmoid function, is the hidden layer output with running indicator samples, is the hidden layer output of the sample without running index, H L is the hidden layer output corresponding to all samples; S35, select the best candidate parameters: calculate Select the largest one in the candidate hidden layer node pool Corresponding candidate hidden layer node parameters As the best candidate hidden layer node parameters, using the selected best candidate hidden node parameters Calculate the corresponding implicit Update the hidden layer mapping matrix H L =[H L-1 ,g L ], Add the best candidate hidden node to the current network model; S36, calculate the model residual after adding new nodes: calculate the residual modulus value of the best candidate hidden layer node: In the formula, Represents the residual when the number of hidden layer nodes is L. represents the residual when the network has L-1 hidden nodes, represents the output matrix of the current network, β L is the output weight of the network; S37, determine whether the model expectations are met after the new node is added: Less than a given desired tolerance ∈ tol When , the model building process ends, otherwise it returns to step S34 and adds nodes to the network in the form of point increments, that is, L = L + 1, until L reaches the maximum number of hidden layer nodes L max until.

7. The industrial process operation index modeling method for cross-time scale sampling according to claim 6 is characterized in that: In step S32, the constructed objective function expression is as follows: In the formula, the first term is the empirical error term of the process data sample with operation indicators, the second term is L2 regularization, and the third term is the cross-time scale constraint term. is the true output value of the i-th process data sample with operating indicators, β j , b j is the output weight, random weight (the "T" in the upper right corner indicates transposition) and random bias of the jth hidden layer node, f(.) indicates the output of the network, ‖·‖ indicates the modulus in the form of the second norm, and h i,j is the cross-time scale constraint coefficient, which is composed of the process data sample set X with operation index M And the process data sample set X without operation index M The distance in the manifold space is obtained, and the specific calculation formula is as follows: In the formula, dist(x i ,x j ) represents the sample x i and x j The Euclidean distance between .

8. The industrial process operation index modeling method for cross-time scale sampling according to claim 7 is characterized in that: In step S32, the matrix form expression of the constructed objective function is as follows: Where Y M Indicates the actual output of the process data sample with running indicators, Indicates the hidden layer output corresponding to the running indicator process data sample, represents the hidden layer output corresponding to the process data sample without operation index, β L represents the output weight of the Lth node. The "T" in the upper right corner represents transposition. F represents the Frobenius norm modulus, tr(·) represents the trace of the matrix, O, P, Q represent the cross-time scale correlation matrix, and the calculation formula is as follows:

9. The industrial process operation index modeling method for cross-time scale sampling according to claim 8 is characterized in that: The following formula is used to construct the output weight update formula of the Lth hidden node: Where I is the unit matrix, 10. The industrial process operation index modeling method for cross-time scale sampling according to claim 6 is characterized in that: In step S3, combined with the idea of ​​error convergence, it is determined whether the output of each sample hidden layer of the Lth newly added node satisfies the supervision constraint. The specific supervision mechanism expression is as follows: In the formula, k = 1, 2, ..., m, m represents the dimension corresponding to the operating index in the multi-sampling period data, 〈·, ·> represents the inner product of the vector, and the superscript "T" represents the transposition operation. represents the supervision constraint corresponding to the kth operating indicator sample when the number of current hidden layer nodes is L, μ L =(1-r) / (L+1), It represents the error corresponding to the kth running index sample when the number of current hidden layer nodes is L-1; if T max The randomly generated weights and biases of the group cannot meet the supervision constraints, and the increase of μ L , r=r+Δr, Δr is the incremental parameter of the parameter r allocation interval; if all values ​​in the parameter r allocation interval are traversed, and the randomly generated weights and biases still cannot meet the supervision constraints, then the candidate parameter λ allocation interval is expanded, λ=λ+Δλ, Δλ is the incremental parameter of the allocation interval.

Images