A Pattern Recognition and Optimization Method for Refining Processes Based on Big Data

Abstract

The present invention relates to a method for pattern recognition and optimization of a refining process based on big data, comprising the following steps: preprocessing the historical data collected during the refining process to obtain standardized data; processing the standardized data by using the principal component analysis method, establishing a model by using the score matrix and drawing a confidence ellipse; calculating the score matrix for the newly collected real-time samples, and projecting them onto the confidence ellipse; the samples projected inside the ellipse are normal points, which can be added to the historical data to establish a new model to achieve adaptive update of the model; the samples projected outside the ellipse are abnormal points, and fault tracing can be carried out according to the fault contribution rate; further, according to the original variables corresponding to the points in the confidence ellipse, an improved path optimization algorithm can be adopted to solve the adjustment method and path of the operating variables.

Description

Technical Field

[0001] The present invention belongs to the technical field of refining process monitoring, and particularly relates to a method for pattern recognition and optimization of refining process based on big data. Background Art

[0002] With the continuous improvement of modern information technology, a large amount of data can be collected in the refining process by a data acquisition system through various measuring instruments. The changes of these data are often related to different production modes of the process. Therefore, effectively monitoring these data is of great significance for improving the production efficiency of the refining process and ensuring its production safety.

[0003] Production processes such as catalytic reforming, catalytic cracking, sulfur recovery, residue hydrotreating, atmospheric and vacuum distillation, and hydrocracking in the refining process have characteristics such as multi-variables, strong interference, large lag, and strong coupling. They are very complex large industrial systems. It is difficult to extract from the historical process data containing numerous operating variables and raw material property variables the influence that can fully reflect each parameter variable on the production process, clarify the types of device production modes, and distinguish the excellent and poor operating regions of the device under various production modes. In addition, it is a challenging task to judge the level of the device operating state according to the current device operation data and timely adjust the process parameters to achieve the optimized operation of the production process.

[0004] At present, the process operation method based on the state at a single time point and the single-factor curve within a time period can no longer meet the requirement of real-time monitoring of the changes in multiple modes of the refining process. In order to comprehensively consider the characteristics that the data collected in the process has numerous variables and the importance of each variable for process monitoring changes with time, it is necessary to consider a multi-factor mode control method that can change with time to realize the collective performance of the mode in the time dimension and the space dimension, and can dynamically monitor each key variable, and uniformly process fault diagnosis, fault prediction, operation safety, dynamic optimization, and static optimization under the concept of the mode. Summary of the Invention

[0005] The purpose of the present invention is to propose a method for pattern recognition and optimization of refining process based on big data aiming at the deficiencies existing in the existing methods. By introducing the principal component analysis method, the key features of numerous variables in the process are extracted, and a confidence ellipse is drawn. The online monitoring of the refining process is realized by using the mapping position of the real-time collected data in the ellipse. This method can be used for monitoring the production mode of the refining process, fault detection and traceability. An improved path optimization algorithm can be further adopted. Through the path optimization between the current working condition position and the corresponding position of the optimal benefit in the confidence ellipse, the optimization direction of the key variables can be given.

[0006] Specifically, the first aspect of the present invention provides a method for recognizing the refining process pattern based on big data, and the method includes the following steps:

[0007] (1) Compose the historical data collected during the refining process into a training sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples in the sample set and n is the number of variables in the sample set;

[0008] (2) Preprocess the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a mean of 0 and a variance of 1;

[0009] (3) Apply the principal component analysis method to X to reduce its dimension from n - dimensional to k - dimensional, and obtain a score matrix T ∈ R m×k and a loading matrix P ∈ R n×k ;

[0010] (4) Use the first two columns of the score matrix T to draw a two - dimensional confidence ellipse;

[0011] (5) Collect new online real - time data Y ∈ R N×n , and preprocess Y using the sample mean and sample variance obtained when preprocessing the training samples in step (2) to obtain standardized data Ym ∈ R N×n ;

[0012] (6) Multiply Ym by the first two columns of the loading matrix P obtained in step (3) to obtain the first two - column score matrix scorey ∈ R N×2 of Ym according to the training samples;

[0013] (7) Use the first column of scorey as the data on the x - axis and the second column of scorey as the data on the y - axis, and map scorey to the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the working condition of the refining process at this time is normal; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the refining process at this time.

[0014] In one or more embodiments, in step (2), the preprocessing method adopts the Z - score standardization method, and the calculation formula is:

[0015]

[0016] where Z = [z1, z2,..., z mis the training data matrix, X represents the standardized data matrix, μ is the mean of the training data, σ is the standard deviation of the training data, and the calculation formulas for μ and σ are:

[0017]

[0018]

[0019] In one or more embodiments, in step (3), principal component analysis is used to perform dimensionality reduction on the obtained X after preprocessing, and the specific steps are as follows:

[0020] (3-a) Calculate the covariance matrix of matrix X, and the calculation formula for the covariance matrix is:

[0021]

[0022] X is an m×n matrix, m is the number of training samples, n is the number of features, T represents the transpose, so the obtained covariance matrix C is an n×n-dimensional matrix;

[0023] (3-b) Calculate the eigenvalues λ i and eigenvectors p i , and sort them in descending order of eigenvalues, and the calculation formula is:

[0024] det(C - Iλ) = 0 (5)

[0025] Cp i = λ i p i (6)

[0026] In formula (5), I is the identity matrix, and we get:

[0027] Eigenvalue matrix: where λ1 > λ2 >... > λ n

[0028] Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0029] (3-c) Establish a principal component model, and the calculation formula is:

[0030] Calculate the information ratio included in each principal component:

[0031]

[0032] Calculate the information ratio included in the largest k principal components:

[0033]

[0034] (3-d) Retain the eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,..., p k ∈ R n×k , and the calculation formula for the score matrix is:

[0035] t i = Xp i , i = 1, 2,..., k (9)

[0036] t i essentially represents the projection of the X vector in the direction of p i , and the score matrix T = [t1, t2,..., t k ∈ R m×k .

[0037] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and the steps for specifically drawing the confidence ellipse using the matrix xdat are as follows:

[0038] (4-a) Calculate the covariance matrix of xdat and invert the covariance matrix to obtain s ∈ R 2×2 ;

[0039] (4-b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and subtract the corresponding average value from each column of xdat for centering processing to obtain xd ∈ R m×2 ;

[0040] (4-c) Calculate the formula xd × s · × xd, and sum each row of the obtained matrix to obtain rd ∈ R m×1 ;

[0041] (4-d) Draw a curve graph of xdat. According to the empirical distribution characteristics of xdat, calculate the percentile of the matrix rd. According to the confidence level C i of the confidence ellipse to be drawn, after sorting rd in ascending order, obtain the value r ∈ R corresponding to the C i th position of rd. Preferably, C i is 95%;

[0042] (4-e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in step (4-a) to obtain the eigenvalue matrix as and the eigenvector matrix as

[0043] (4-f) Using r in step (4-d) and D in step (4-e), the major axis and minor axis of the confidence ellipse can be obtained, and the calculation formula is:

[0044]

[0045]

[0046] Among them, a is the major axis and b is the minor axis;

[0047] (4-g) Taking xm in step (4-b) as the center point of the ellipse, draw a confidence ellipse according to the center point, the major axis and the minor axis of the ellipse. The formula of the ellipse is:

[0048]

[0049] Among them, xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat.

[0050] In one or more embodiments, in step (6), the calculation formula of scorey is:

[0051] scorey = Ym × [p1, p2] (13)

[0052] In the formula, [p1, p2] ∈ R n×2 are the first two columns of the load matrix P.

[0053] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the set of data is mapped inside the ellipse.

[0054] In one or more embodiments, step (7) further includes adding the data collected under normal conditions to the historical data for re-modeling.

[0055] In one or more embodiments, step (7) further includes adding the data collected under normal conditions to the historical data and removing the same number of early historical data for re-modeling.

[0056] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data for fault tracing;

[0057] Preferably, the following method is used to calculate the SPE contribution rate for abnormal data for fault tracing:

[0058] Assume that the abnormal data is x ∈ R 1×n , and the calculation formula of its SPE contribution rate is:

[0059]

[0060] Among them, contspe(i) represents the SPE contribution corresponding to the i-th variable, and ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the faulty variable.

[0061] In one or more embodiments, step (4) further includes: combining historical data with on-site process knowledge to flag the drawn confidence ellipse to achieve regional division of data with different performance levels. The regional division preferably includes: first, according to process knowledge and historical benefit statistical data, find the benefit corresponding to the working condition according to the time series, perform different levels of division according to the benefit level, then find the distribution of the historical data corresponding to each level in the ellipse and mark them with different marks, such as marking them with different colors.

[0062] In one or more embodiments, step (4) further includes: setting different labels for historical data with different performance levels according to the performance indicators of the refining process, projecting the first two columns of the score matrix of the historical data onto the confidence ellipse, and dividing the confidence ellipse through the labels carried by the data points in different regions.

[0063] The present invention also provides a method for optimizing the refining process mode based on big data. The method for optimizing the refining process mode based on big data includes the method for identifying the refining process mode based on big data as described in any one of the embodiments herein. And step (7) in the method for identifying the refining process mode based on big data further includes: for the samples mapped inside the ellipse, judge the current performance level according to their positions in combination with the regional division of the ellipse in step (4). If it is in a non-optimal region, adjust the key variables to make it switch to the desired better performance level; preferably, set the variable corresponding to the largest coefficient in the principal component as the key variable to be adjusted during optimization, and then determine the adjustment direction according to the correlation between the key variable and the change in the position of the projection point of the online data in the confidence ellipse to achieve the optimization of the production process; preferably, the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the refining process is adjusted, and it is observed whether the data moves to the desired region through real-time monitoring.

[0064] The present invention also provides another method for optimizing the refining process mode based on big data. The method for optimizing the refining process mode based on big data includes the method for identifying the refining process mode based on big data as described in any one of the embodiments herein, and step (8): obtaining the original variables corresponding to the points in the confidence ellipse by using the standard deviation and mean obtained in the principal component analysis process, so as to obtain the benefit value corresponding to this point; distinguishing the distribution of high and low benefit values in the confidence ellipse; when the operating mode of the working condition is at a certain point in the confidence ellipse, adopting a path optimization algorithm to obtain the fastest moving trajectory from the current position to the optimal position, and performing an inverse transformation on the points corresponding to the trajectory to obtain the change mode of the operating conditions, thereby guiding the optimization of the operation of the production device.

[0065] In one or more embodiments, in step (8), the original variable X is calculated by formula (15) ori :

[0066] X ori =(defen×PC T )·×std(X m )+mean(X m ) (15)

[0067] where defen=(x,y)∈R 1×2 is the abscissa and ordinate of the point in the confidence ellipse, PC=[p1,p2]∈R J×2 is the first two columns of the load matrix when performing principal component analysis modeling, std(X m )is the variance of the basic sample used in principal component analysis, mean(X m )is the mean of the basic sample used in principal component analysis, m is the number of samples, and X ori is the original variable after inverse transformation corresponding to the position of defen=(x,y)∈R 1×2 in the ellipse;

[0068] Calculate the input-output benefit index according to formula (16):

[0069]

[0070] where is the output of the i-th product, is the price of the i-th product, is the feeding amount of the i-th raw material, is the price of the i-th raw material, and profit is the benefit index.

[0071] In one or more embodiments, in step (8), the path optimization algorithm adopts an improved A* algorithm, and determines the search direction and the next node to reach through the improved evaluation function f(x) shown in formula (17):

[0072]

[0073] Among them, g(x) is the cost function, representing the actual cost required to reach the current node x from the starting node; h(x) is the heuristic function, representing the estimated cost required to reach the target node from the current node x; profit(x) represents the economic benefit corresponding to the selected node.

[0074] In one or more embodiments, the refining process is a catalytic reforming process, a catalytic cracking process, a sulfur recovery process, a residue hydrotreating process, an atmospheric and vacuum distillation process, or a hydrocracking process, particularly a catalytic reforming process.

[0075] The second aspect of the present invention provides a method for identifying and optimizing the operating state of a catalytic cracking unit based on data-driven, and the method includes the following steps:

[0076] (1) Compose the historical production process data of the catalytic cracking unit under normal operating conditions including operation data and feedstock data into a training data sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples of the training data, and n is the number of variables of the training data;

[0077] (2) Preprocess the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a mean of 0 and a variance of 1;

[0078] (3) Apply the principal component analysis method to X to reduce its dimension from n-dimensional to k-dimensional, and calculate the score matrix T ∈ R m×k and the loading matrix P ∈ R n×k ;

[0079] (4) Use the first two columns of the score matrix T to draw a two-dimensional confidence ellipse; define the production mode, set the evaluation index under the defined production mode, divide the evaluation index under the defined production mode into grades according to requirements, screen out the data in the optimized operating state, and project the samples corresponding to the data in the optimized operating state onto the two-dimensional confidence ellipse to obtain the optimized region under the defined production mode;

[0080] (5) Collect the online real-time data Y ∈ R N×n , and preprocess Y to obtain Ym ∈ R N×n using the mean and variance obtained when preprocessing the training data in step (2);

[0081] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two column score matrix scorey of Ym based on the training samples, where scorey ∈ R N×2 ;

[0082] (7) Use the first column of scorey as the data for the x-axis and the second column of scorey as the data for the y-axis, and map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the catalytic cracking unit is in a normal operating condition at this time; if the sample point is mapped outside the ellipse, it indicates that there may be an abnormality in the operating state of the catalytic cracking unit at this time; further, for the samples mapped inside the two-dimensional confidence ellipse, determine whether they are within the optimized region of the defined production mode according to their positions; if the mapping point is within the optimized region, it indicates that the unit is in an optimized operating state under the defined production mode at this time; if the mapping point is not within the optimized region, it indicates that the unit is in a non-optimized operating state under the defined production mode.

[0083] In one or more embodiments, the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.

[0084] In one or more embodiments, in step (3), principal component analysis is used to perform dimensionality reduction on X obtained after pretreatment, and the specific steps are as follows, and the specific operation steps are as described in step (3) of the first aspect of the present invention.

[0085] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and use the matrix xdat to specifically draw the confidence ellipse, and the specific operation steps are as described in step (4) of the first aspect of the present invention.

[0086] In one or more embodiments, in step (4), define the production mode, set the evaluation index under the defined production mode, set different labels for the historical data according to the evaluation index of the defined production mode, and divide the confidence ellipse into regions by the labels carried by the data points in different regions when projecting the score matrix of the historical data; or divide the evaluation index under the defined production mode into grades according to requirements, screen out the data in the optimized operating state, and project the samples corresponding to the data in the optimized operating state onto the two-dimensional confidence ellipse to obtain the optimized region under the defined production mode.

[0087] In one or more embodiments, in step (6), the calculation formula of scorey is as described in step (6) of the first aspect of the present invention.

[0088] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that this set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this set of data is mapped inside the ellipse.

[0089] In one or more embodiments, step (7) further includes adding the collected data under normal operating conditions to the historical data for re-modeling.

[0090] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data to trace the fault, and the specific operation method is as described in step (7) of the first aspect of the present invention.

[0091] In one or more embodiments, in step (7), for the samples mapped inside the two-dimensional confidence ellipse, it is judged whether they are within the optimized region of the defined production mode according to their positions; if the mapping points are within the optimized region, it indicates that the device is in the optimized operating state under the defined production mode at this time; if the mapping points are not within the optimized region, it indicates that the device is in the non-optimized operating state under the defined production mode at this time.

[0092] In one or more embodiments, the method further includes step (8): for the online real-time data in the non-optimized operating state, by observing the projection position of the data in the two-dimensional confidence ellipse under the defined production mode, among the principal component components obtained by dimensionality reduction of the device operation data and raw material data, find 1-3 variables with the largest coefficients corresponding to each principal component as the core variables for adjusting the operating state of the fluid catalytic cracking unit. According to the correlation between the core variables and the change of the projection point position, determine the adjustment direction and change the values of the core variables to achieve the optimization of the production process.

[0093] In one or more embodiments, in step (8), the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the sulfur recovery unit is adjusted, and it is observed whether the data moves to the desired region through real-time monitoring.

[0094] In one or more embodiments, steps (1)-(8) of the second aspect of the present invention are as described in any one of the embodiments of the first aspect of the present invention.

[0095] The third aspect of the present invention provides a method for identifying and optimizing the operating state of a sulfur recovery unit based on data driving, and the method includes the following steps:

[0096] (1) Composing the historical data of the production process of the sulfur recovery unit under normal operating conditions into a training data sample set Z = [z1, z2,..., zi ,..., z n ∈ R m×n , where m is the number of samples of the training data, and n is the number of variables of the training data;

[0097] (2) Preprocess the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n ;

[0098] (3) Use the principal component analysis method to reduce the dimension of X from n dimensions to k dimensions, and calculate the score matrix T ∈ R m×k and the load matrix P ∈ R n×k ;

[0099] (4) Use the first two columns of the score matrix T to draw a two-dimensional confidence ellipse; set the production mode of the sulfur recovery unit, and according to the key index variables and production requirements under the set production mode, classify the set production mode, screen out the data in the optimized operation state, and project the samples corresponding to the data in the optimized operation state onto the two-dimensional confidence ellipse to obtain the optimized area under the set production mode;

[0100] (5) Collect online real-time data Y ∈ R N×n , and preprocess Y using the average value and variance obtained when preprocessing the training data in step (2) to obtain Ym ∈ R N×n ;

[0101] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two column score matrix scorey ∈ R of Ym according to the training samples N×2 ;

[0102] (7) Use the first column of scorey as the data on the x-axis and the second column of scorey as the data on the y-axis, and map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it means that the sulfur recovery unit is in the normal working condition at this time; if the sample point is mapped outside the ellipse, it means that there may be an abnormality in the operation state of the sulfur recovery unit at this time.

[0103] In one or more embodiments, the preprocessing method in step (2) is as described in step (2) of the first aspect of the present invention.

[0104] In one or more embodiments, in step (3), the principal component analysis method is used to reduce the dimension of X obtained after preprocessing, and the specific operation steps are as described in step (3) of the first aspect of the present invention.

[0105] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and specifically draw a confidence ellipse using the matrix xdat. The specific operation steps are as described in step (4) of the first aspect of the present invention.

[0106] In one or more embodiments, in step (4), according to the key characterization indexes of the set production mode of the sulfur recovery unit, different labels are set for the historical data. When projecting the score matrix of the historical data, the confidence ellipse is divided into an optimization region and a non-optimization region by the labels carried by the data points in different regions.

[0107] In one or more embodiments, in step (6), the calculation formula of scorey is as described in step (6) of the first aspect of the present invention.

[0108] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the set of data is mapped inside the ellipse.

[0109] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data to trace the fault. The specific operation method is as described in step (7) of the first aspect of the present invention.

[0110] In one or more embodiments, in step (7), for the samples mapped inside the two-dimensional confidence ellipse, it is judged whether they are in the defined optimization region of the production mode according to their positions; if the mapping points are in the optimization region, it indicates that the device is in the optimized operation state under the defined production mode at this time; if the mapping points are not in the optimization region, it indicates that the device is in the non-optimized operation state under the defined production mode at this time.

[0111] In one or more embodiments, the method further includes step (8): for the online real-time data in the non-optimized operation state, by observing the projection position of the data in the two-dimensional confidence ellipse under the defined production mode, find 1 to 3 variables with the largest corresponding coefficients for each principal component from each principal component as the core variables for adjusting the operation state of the sulfur recovery unit. According to the correlation between the core variables and the change of the projection point position, determine the adjustment direction and change the values of the core variables to optimize the production process.

[0112] In one or more embodiments, in step (8), the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the sulfur recovery unit is adjusted, and it is observed whether the data moves towards the desired region through real-time monitoring.

[0113] In some embodiments, the method for identifying and optimizing the operating state of a sulfur recovery unit based on data driving of the present invention includes: collecting historical process data during the process of the sulfur recovery unit, including process operations, raw material evaluation, hydrogen composition analysis, and product distribution of the unit, and screening out several groups of production data including unit operating variables and dependent variables; adopting the principal component analysis method to reduce the dimension of the above data and draw a two-dimensional confidence ellipse; setting multiple production modes of the unit according to indicators such as sulfur recovery rate and H2S content in the purified tail gas, and screening out the optimal operating state regions of the unit under each mode according to the indicators and historical data; collecting real-time production data of the unit operation, mapping it to the two-dimensional confidence interval through the above method, and identifying whether the current operating state of the unit is in the operating optimization region; if the unit is in a non-optimized operating state, then by calculating the SPE contribution rate of each variable, find the variables with larger values for regulation and control, so that the production state of the unit gradually and timely returns to the optimal operating region under the corresponding production mode.

[0114] In one or more embodiments, steps (1)-(8) of the third aspect of the present invention are as described in any one of the embodiments of the first aspect of the present invention.

[0115] The fourth aspect of the present invention provides a method for identifying and optimizing the operating state of a residue hydrotreating unit based on data driving, and the method includes the following steps:

[0116] (1) Composing the historical data of the production process of the residue hydrotreating unit into a training data sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples of the training data, and n is the number of variables of the training data;

[0117] (2) Preprocessing the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a mean of 0 and a variance of 1;

[0118] (3) Using the principal component analysis method to reduce the dimension of X from n dimensions to k dimensions, and calculating the score matrix T ∈ R m×k and the load matrix P ∈ R n×k ;

[0119] (4) Using the first two columns of the score matrix T, draw a two-dimensional confidence ellipse; set the production mode of the residue hydrotreating unit, classify the data according to the index variables under the set production mode, screen out the data in the optimal operating state, and project the samples corresponding to the data in the optimal operating state onto the two-dimensional confidence ellipse to obtain the optimal region under the set production mode;

[0120] (5) Collect online real-time data Y ∈ R N×n , and preprocess Y using the mean and variance obtained when preprocessing the training data in step (2) to obtain Ym ∈ R N×n ;

[0121] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two columns of the score matrix scorey ∈ R of Ym according to the training samples N×2 ;

[0122] (7) Use the first column of scorey as the data on the x-axis and the second column of scorey as the data on the y-axis, and map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it means that the residue hydrotreating unit is in the normal operating condition at this time; if the sample point is mapped outside the ellipse, it means that there may be an abnormality in the operating state of the residue hydrotreating unit at this time.

[0123] In one or more embodiments, the preprocessing method in step (2) is as described in step (2) of the first aspect of the present invention.

[0124] In one or more embodiments, in step (3), principal component analysis is used to perform dimensionality reduction on X obtained after preprocessing, and the specific operation steps are as described in step (3) of the first aspect of the present invention.

[0125] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and use the matrix xdat to specifically draw the confidence ellipse, and the specific operation steps are as described in step (4) of the first aspect of the present invention.

[0126] In one or more embodiments, in step (4), according to the key characterization indexes of the set production mode of the residue hydrotreating unit, different labels are set for the historical data, and the confidence ellipse is divided into an optimal region and a non-optimal region by the labels carried by the data points in different regions when projecting the score matrix of the historical data.

[0127] In one or more embodiments, in step (6), the calculation formula of scorey is as described in step (6) of the first aspect of the present invention.

[0128] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that this set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this set of data is mapped inside the ellipse.

[0129] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data to trace the fault, and the specific operation method is as described in step (7) of the first aspect of the present invention.

[0130] In one or more embodiments, in step (7), for the samples mapped inside the two-dimensional confidence ellipse, it is judged whether they are within the optimized area of the set production mode according to their positions; if the mapping points are within the optimized area, it indicates that the device is in the optimized operation state under the set production mode at this time; if the mapping points are not within the optimized area, it indicates that the device is in the non-optimized operation state under the set production mode at this time.

[0131] In one or more embodiments, the method further includes step (8): for the on-line real-time data in the non-optimized operation state, by observing the projection position of the data in the two-dimensional confidence ellipse under the defined production mode, 1 to 3 variables with the largest coefficients corresponding to each principal component are found from each principal component as the core variables for adjusting the operation state of the residue hydrotreating unit. According to the correlation between the core variables and the change of the projection point position, the adjustment direction is determined, and the values of the core variables are changed to optimize the production process.

[0132] In one or more embodiments, in step (8), the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the residue hydrotreating unit is adjusted, and it is observed whether the data moves to the desired area through real-time monitoring.

[0133] In one or more embodiments, the data-driven method for identifying and optimizing the operating state of a residue hydrotreating unit of the present invention includes: collecting historical process data during the residue hydrotreating process, including process operations, raw material and product property analysis, and unit product distribution, and screening out several groups of production data including unit operating variables and dependent variables; defining an optimized operating mode for the residue hydrotreating unit according to the unit operating characteristics and data, such as maximum unit benefit, minimum energy consumption, etc.; using the principal component analysis method to process the data under different production modes, establishing a model using the score matrix and drawing a confidence ellipse; projecting the score matrix calculated from the newly collected real-time samples into the model ellipse, and judging whether the current unit production process is in an optimized state under this production mode according to the position of the sample point in the ellipse; adding the samples projected into the ellipse to the historical data to establish a new model and realizing the adaptive update of the model; if the sample point is projected outside the two-dimensional confidence ellipse, it indicates that the unit is in a non-optimized operating state. At this time, the model input variables need to be analyzed according to the fault contribution rate to realize fault tracing; screening out the variables with large fault contribution rates and adjusting them according to the principal components and score results, so that the unit production state gradually and timely returns to the optimal operating area under the corresponding production mode, realizing the optimization of the unit.

[0134] In one or more embodiments, steps (1)-(8) of the fourth aspect of the present invention are as described in any embodiment of the first aspect of the present invention.

[0135] The fifth aspect of the present invention provides a method for monitoring and optimizing the atmospheric and vacuum distillation process based on big data. The method includes the following steps:

[0136] (1) Composing the historical data of the atmospheric and vacuum distillation process into a training sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples in the sample set and n is the number of variables in the sample set;

[0137] (2) Preprocessing the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a mean of 0 and a variance of 1;

[0138] (3) Applying the principal component analysis method to X to reduce its dimension from n to k, and obtaining a score matrix T ∈ R m×k and a load matrix P ∈ R n×k ;

[0139] (4) Using the first two columns of the score matrix T to draw a two-dimensional confidence ellipse;

[0140] (5) Collecting new online real-time data Y ∈ RN×n When preprocessing Y using the sample mean and sample variance obtained in step (2) during the preprocessing of the training samples, standardized data Ym ∈ R is obtained. N×n ;

[0141] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the score matrix scorey ∈ R of the first two columns of Ym based on the training samples. N×2 ;

[0142] (7) Using the first column of scorey as the data on the x-axis and the second column of scorey as the data on the y-axis, map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the operating condition of the atmospheric and vacuum distillation process at this time is normal; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the atmospheric and vacuum distillation process at this time.

[0143] In one or more embodiments, the preprocessing method in step (2) is as described in step (2) of the first aspect of the present invention.

[0144] In one or more embodiments, in step (3), principal component analysis is used to perform dimensionality reduction on X obtained after preprocessing, and the specific operation steps are as described in step (3) of the first aspect of the present invention.

[0145] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and use the matrix xdat to specifically draw the confidence ellipse, and the specific operation steps are as described in step (4) of the first aspect of the present invention.

[0146] In one or more embodiments, in step (6), the calculation formula of scorey is as described in step (6) of the first aspect of the present invention.

[0147] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that this set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this set of data is mapped inside the ellipse.

[0148] In one or more embodiments, step (7) further includes adding the data collected under normal operating conditions to the historical data for re-modeling.

[0149] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data for fault tracing, and the specific operation method is as described in step (7) of the first aspect of the present invention.

[0150] In one or more embodiments, step (4) further includes: flagging the drawn confidence ellipse by combining historical data with on-site process knowledge to achieve regional division of data with different performance levels. The regional division preferably includes: first calculating the KPI values (benefit values) of various performance indicators (including economic benefits, production energy consumption, product yield, etc.) in the historical data according to process knowledge, then dividing the KPI values into different levels according to the actual process, and then finding the distribution of the historical data corresponding to each level in the ellipse and marking it with different colors (such as color marking). The benefit value can be found in the historical data or calculated based on the historical data.

[0151] In one or more embodiments, in step (4), different labels are set for the historical data according to the performance indicators of the atmospheric and vacuum distillation process. Therefore, when projecting the score matrix of the historical data, the confidence ellipse can be divided by the labels carried by the data points in different regions.

[0152] In one or more embodiments, the method further includes step (8): for the samples mapped inside the ellipse, judge the current performance level according to their positions in combination with the regional division of the ellipse in step (4).

[0153] In one or more embodiments, step (8) further includes: for the samples at non-optimal performance levels, adjust the key variables to make them switch to the desired better performance level; preferably, set the variable corresponding to the largest coefficient in the principal component as the key variable to be adjusted during optimization, and then determine the adjustment direction according to the correlation between the key variable and the position change of the projection point of the online data in the confidence ellipse to achieve optimization of the production process; preferably, the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the atmospheric and vacuum distillation process is adjusted, and it is observed whether the data moves to the desired region through real-time monitoring.

[0154] In one or more embodiments, the big data-based atmospheric and vacuum distillation process pattern monitoring and optimization method of the present invention includes: classifying a large amount of historical data collected in the atmospheric and vacuum distillation process according to their different physical meanings and defining them as different production modes; processing the data under different production modes by using the principal component analysis method, establishing a model by using the score matrix and drawing a confidence ellipse; projecting the score matrix calculated from newly collected real-time samples into the model, and judging the production mode in which the current process is located according to the position of the sample point in the ellipse; adding the samples projected into the ellipse to the historical data to establish a new model to realize the adaptive update of the model; the samples projected outside the ellipse are abnormal points, and fault tracing is carried out according to the fault contribution rate; when modeling, variables that play a major role in the change of the process state under different production modes can be obtained, and these variables are adjusted to realize mode optimization.

[0155] In one or more embodiments, steps (1)-(8) of the fifth aspect of the present invention are as described in any embodiment of the first aspect of the present invention.

[0156] The sixth aspect of the present invention provides a big data-based hydrocracking process pattern recognition and optimization method, and the method includes the following steps:

[0157] (1) Using the historical data of the hydrocracking process to construct a training data sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples in the sample set and n is the number of variables in the sample set;

[0158] (2) Preprocessing the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a variance of 1 and a mean of 0;

[0159] (3) Using the principal component analysis method to reduce X from n dimensions to k dimensions, and obtaining a score matrix T ∈ R m×k and a load matrix P ∈ R n ×k ;

[0160] (4) Using the first two columns of the score matrix T to draw a two-dimensional confidence ellipse;

[0161] (5) Collecting online real-time data Y ∈ R N×n , and preprocessing Y with the sample mean and sample variance obtained when preprocessing the training samples in step (2) to obtain standardized data Ym ∈ R N×n ;

[0162] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two columns of the score matrix scorey ∈ R of Ym according to the training samples. N×2 ;

[0163] (7) Take the first column of scorey as the data on the x-axis and the second column of scorey as the data on the y-axis, and map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the working condition of the hydrocracking process at this time is normal; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the hydrocracking process at this time.

[0164] In one or more embodiments, the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.

[0165] In one or more embodiments, in step (3), the principal component analysis method is used to perform dimensionality reduction on X obtained after pretreatment, and the specific operation steps are as described in step (3) of the first aspect of the present invention.

[0166] In one or more embodiments, in step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and use the matrix xdat to specifically draw the confidence ellipse, and the specific operation steps are as described in step (4) of the first aspect of the present invention.

[0167] In one or more embodiments, in step (4), combining historical data with on-site process knowledge can perform flagging processing on the drawn confidence ellipse to achieve regional division of data with different performance levels. In some embodiments, the regional division includes: first calculating the KPI values (benefit values) of various performance indicators (including economic benefits, production energy consumption, product yield, etc.) in the historical data according to process knowledge, then performing different-level divisions on the KPI values according to the actual process, and then finding the distribution of the historical data corresponding to each level in the ellipse and adding different marks (such as color marks). The benefit values can be found in the historical data or calculated based on the historical data.

[0168] In one or more embodiments, in step (4), different labels are set for the historical data according to the performance indicators of the hydrocracking process, so that when projecting the score matrix of the historical data, the confidence ellipse can be divided by the labels carried by the data points in different regions.

[0169] In one or more embodiments, in step (6), the calculation formula of scorey is as described in step (6) of the first aspect of the present invention.

[0170] In one or more embodiments, in step (7), each set of data in scorey is substituted into formula (12) and compared with the value 1. If it is greater than 1, it means that this set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this set of data is mapped inside the ellipse.

[0171] In one or more embodiments, step (7) further includes adding the data collected under normal conditions to the historical data for re-modeling.

[0172] In one or more embodiments, step (7) further includes calculating the SPE contribution rate for abnormal data to trace the fault, and the specific operation method is as described in step (7) of the first aspect of the present invention.

[0173] In one or more embodiments, step (4) further includes: flagging the drawn confidence ellipse by combining historical data with on-site process knowledge to achieve regional division of data with different performance levels. The regional division preferably includes: first calculating the KPI value of the performance index in the historical data according to process knowledge, then dividing the KPI value into different levels according to the actual process, and then finding the distribution of the historical data corresponding to each level in the ellipse and marking them differently.

[0174] In one or more embodiments, the method further includes step (8): For the samples mapped inside the ellipse, determine the current performance level according to their positions in combination with the regional division of the ellipse in step (4).

[0175] In one or more embodiments, step (8) further includes: For the samples at non-optimal performance levels, adjust the key variables to make them convert to the desired better performance levels; preferably, set the variable corresponding to the largest coefficient in the principal component as the key variable to be adjusted during optimization, and then determine the adjustment direction according to the correlation between the key variable and the change in the position of the projection point of the on-line data in the confidence ellipse to achieve optimization of the production process; preferably, the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, adjust the production mode of the hydrocracking process, and observe whether the data moves to the desired area through real-time monitoring.

[0176] In one or more embodiments, the big data-based hydrocracking process pattern recognition and optimization method of the present invention includes: classifying a large amount of historical data collected during the hydrocracking process according to their different physical meanings and defining them as different production patterns; processing the data under different production patterns using the principal component analysis method, establishing a model using the score matrix and drawing a confidence ellipse; projecting the score matrix calculated from newly collected real-time samples into the model, and judging the production pattern in which the current process is located according to the position of the sample point in the ellipse; adding the samples projected into the ellipse to the historical data to establish a new model to achieve adaptive update of the model; the samples projected outside the ellipse are outliers, and fault tracing is performed according to the fault contribution rate; variables that play a major role in the change of the process state under different production patterns can be obtained during modeling, and these variables are adjusted to achieve pattern optimization.

[0177] In one or more embodiments, steps (1)-(8) of the sixth aspect of the present invention are as described in any embodiment of the first aspect of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0178] Figure 1 is a schematic process flow diagram of catalytic reforming.

[0179] Figure 2 is the overall process flow block diagram of the refinery process pattern recognition and optimization method of the present invention.

[0180] Figure 3 is the confidence ellipse drawn from the normal data in Example 1.

[0181] Figure 4 is the partition diagram of the confidence ellipse according to different mode data in Example 1.

[0182] Figure 5 is the SPE contribution rate diagram of the fault point in Example 1.

[0183] Figure 6 is the benefit curve corresponding to the historical data in Example 1.

[0184] Figure 7 is the projection interval corresponding to the projection of the historical data in the confidence ellipse and the high-benefit points in Example 1.

[0185] Figure 8 is the projection interval corresponding to the projection of the historical data in the confidence ellipse and the low-benefit points in Example 1.

[0186] Figure 9 is the optimized movement trajectory of the operating variables in Example 1.

[0187] Figure 10 is the overall process flow block diagram of the catalytic cracking unit operation state recognition and optimization method of the present invention.

[0188] Figure 11 It is the confidence ellipse drawn based on the production data of the fluid catalytic cracking unit in Example 2.

[0189] Figure 12 It is the optimized operation area (the circled points in the broken line) of the unit selected according to the optimized target variable - economic benefit broken line in Example 2.

[0190] Figure 13 It is the projection diagram of the optimized area data corresponding to a certain production mode in Example 2 on the confidence ellipse.

[0191] Figure 14 It is the contribution rate diagram of variables to the production mode evaluation target in Example 2.

[0192] Figure 15 It is the optimization trajectory of the fluid catalytic cracking unit in Example 2.

[0193] Figure 16 It is the overall process flow block diagram of the method for identifying and optimizing the operating state of the sulfur recovery unit of the present invention.

[0194] Figure 17 It is the confidence ellipse drawn based on the production data of the sulfur recovery unit in Example 3.

[0195] Figure 18 It is the optimized operation area of the unit selected according to the optimized target variable broken line (sulfur recovery rate) in Example 3.

[0196] Figure 19 It is the projection diagram of the optimized area data corresponding to a certain production mode in Example 3 on the confidence ellipse.

[0197] Figure 20 It is the contribution rate diagram of variables to the production mode evaluation target in Example 3.

[0198] Figure 21 It is the optimization trajectory of the sulfur recovery unit in Example 3.

[0199] Figure 22 It is the overall process flow block diagram of the method for identifying and optimizing the operating state of the residue hydrotreating unit based on data driving of the present invention.

[0200] Figure 23 It is the confidence ellipse drawn based on the production data of the residue hydrotreating unit in Example 4.

[0201] Figure 24 It is the distribution of different production modes of the residue hydrotreating unit in Example 4.

[0202] Figure 25It is the projection diagram of the optimization region data corresponding to the production sample with emphasis on the heavy oil hydrogenation yield in Example 4 on the confidence ellipse.

[0203] Figure 26 It is the contribution rate diagram of variables to the production mode evaluation target in Example 4.

[0204] Figure 27 It is the optimization trajectory of the residue hydrotreating unit in Example 4.

[0205] Figure 28 It is the schematic process flow diagram of atmospheric and vacuum distillation.

[0206] Figure 29 It is the overall process flow block diagram of the atmospheric and vacuum distillation process mode monitoring and optimization method based on big data of the present invention.

[0207] Figure 30 It is the confidence ellipse drawn according to the normal data in Example 5.

[0208] Figure 31 It is the partition diagram of the confidence ellipse according to the mode data with energy consumption as the optimization target in Example 5.

[0209] Figure 32 It is the partition diagram of the confidence ellipse according to the mode data with light ends yield optimization as the optimization target in Example 5.

[0210] Figure 33 It is the partition diagram of the confidence ellipse according to the mode data with total draw as the optimization target in Example 5.

[0211] Figure 34 It is the projection of the real-time collected data on the ellipse in Example 5.

[0212] Figure 35 It is the SPE contribution rate diagram of the fault point in Example 5.

[0213] Figure 36 It is the process flow diagram of hydrocracking.

[0214] Figure 37 It is the overall process flow block diagram of the hydrocracking process mode recognition and optimization method based on big data of the present invention.

[0215] Figure 38 It is the confidence ellipse drawn according to the normal data in Example 6.

[0216] Figure 39 It is the partition diagram of the confidence ellipse according to the mode data with total liquid yield as the optimization target in Example 6.

[0217] Figure 40 It is the partition diagram of the confidence ellipse according to the mode data with middle oil yield as the optimization target in Example 6.

[0218] Figure 41 It is the partition diagram of the confidence ellipse in Embodiment 6 according to the pattern data with the value increment as the optimization target.

[0219] Figure 42 It is the projection of the data collected in real time in Embodiment 6 on the ellipse.

[0220] Figure 43 It is the SPE contribution rate diagram of the fault point in Embodiment 6. Detailed implementation manners

[0221] To enable those skilled in the art to understand the features and effects of the present invention, the following provides a general description and definition of the terms and expressions mentioned in the specification and claims. Unless otherwise specified, all technical and scientific terms used herein shall have the ordinary meaning understood by those skilled in the art with respect to the present invention. In case of conflict, the definition in this specification shall prevail.

[0222] In this article, for the sake of brevity of description, not all possible combinations of the technical features in each embodiment or example are described. Therefore, as long as there is no contradiction in the combination of these technical features, the technical features in each embodiment or example can be combined arbitrarily, and all possible combinations should be considered as the scope described in this specification.

[0223] The method for recognizing the pattern of the oil refining process based on big data of the present invention includes the following steps:

[0224] (1) Compose the historical data collected in the oil refining process into a training sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples in the sample set and n is the number of variables in the sample set;

[0225] (2) Preprocess the training data sample set to obtain standardized data X = [x1, x2,..., x n ∈ R m×n with a mean of 0 and a variance of 1;

[0226] (3) Apply the principal component analysis method to X to reduce its dimension from n to k, and obtain a score matrix T ∈ R m×k and a loading matrix P ∈ R n×k ;

[0227] (4) Use the first two columns of the score matrix T to draw a two-dimensional confidence ellipse;

[0228] (5) Collect new online real-time data Y ∈ R N×nWhen preprocessing Y using the sample mean and sample variance obtained during the preprocessing of the training samples in step (2), standardized data Ym ∈ R is obtained. N×n ;

[0229] (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two column score matrix scorey ∈ R of Ym based on the training samples. N×2 ;

[0230] (7) Use the first column of scorey as the data for the x-axis and the second column of scorey as the data for the y-axis, and map scorey onto the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the operating condition of the refinery process at this time is a normal operating condition; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the refinery process at this time.

[0231] In this article, the refinery process has the conventional meaning in the art, and generally refers to a series of process processing processes experienced by crude oil and its intermediate products during the petroleum refining process, including but not limited to atmospheric distillation, vacuum distillation, catalytic cracking, catalytic reforming, hydrocracking, sulfur recovery, residue hydrotreating, atmospheric and vacuum distillation, delayed coking, refinery gas processing, and alkylation, etc. In some embodiments, the refinery process is a catalytic reforming process, a catalytic cracking process, a sulfur recovery process, a residue hydrotreating process, an atmospheric and vacuum distillation process, or a hydrocracking process.

[0232] In this article, the refinery process mode includes the operating condition and production mode of the refinery process. The operating condition can be divided into a normal operating condition (an operating condition without faults) and an abnormal operating condition (an operating condition with faults). The normal operating condition can be divided into different production modes according to different process and performance evaluation indicators, such as an economic benefit mode, a unit energy consumption mode, a product yield mode, etc. According to the different evaluation indicators corresponding to various production modes, each production mode can be divided into an optimized operating state and a non-optimized operating state. For example, the economic benefit mode can be divided into a high economic benefit mode and a low economic benefit mode, the unit energy consumption mode can be divided into a high unit energy consumption mode and a low unit energy consumption mode, the product yield mode can be divided into a high product yield mode and a low product yield mode, etc.

[0233] In step (1), historical data collected during the oil refining process is used as a training sample set. In this article, a sample generally refers to a set of data obtained by collecting each selected variable at a certain time point, including multiple variables. When the sample set is presented in the form of a matrix, a row of data (row vector) in the matrix generally represents a sample, and each column of data (column vector) in the matrix generally corresponds to a different variable, and the number of variables (number of columns) is the dimension of the data. Variables include independent variables (also referred to as operating variables in this article) and dependent variables. The training sample set can come from historical data in the factory real-time database. Usually, the training samples are historical data under normal operating conditions. According to the device data record points and property analysis items, the factory production process historical data can be collected in batches, and the operating points (operating data) and raw material property items (raw material data) with normal data and possible impact on the device production process can be roughly screened out based on prior knowledge, and the catalytic cracking unit production process historical data containing such operating data and raw material data is used as training data for modeling.

[0234] In step (2), the training samples are standardized by obtaining the mean value and standard deviation of the training sample set to obtain standardized data. In the present invention, the method of standardizing data based on the mean value and standard deviation is conventional in the art. The standardized data presents a normal distribution with a mean of 0 and a variance of 1. In some embodiments, the preprocessing method adopts the Z-score standardization method, and the calculation formula is:

[0235]

[0236] Where Z = [z1,z2,...,z m ] is the training data matrix, X represents the standardized data matrix, μ is the mean of the training data, σ is the standard deviation of the training data, and the calculation formulas of μ and σ are:

[0237]

[0238]

[0239] In step (3), the principal component analysis method is used to reduce the dimension of X and obtain the score matrix T∈R m×k and the load matrix P∈R n×k A method known in the art or a method provided by the present invention may be used. In some embodiments, the dimension reduction process of X obtained after preprocessing is performed by principal component analysis according to the following steps:

[0240] (3-a) Calculate the covariance matrix of matrix X, the formula is:

[0241]

[0242] X is an m×n matrix, where m is the number of training samples and n is the number of features. T represents the transpose. Therefore, the obtained covariance matrix C is an n×n-dimensional matrix;

[0243] (3-b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C, and sort them in descending order of eigenvalues. The calculation formula is:

[0244] det(C - Iλ) = 0 (5)

[0245] Cp i = λ i p i (6)

[0246] In formula (5), I is the identity matrix, and we get:

[0247] Eigenvalue matrix: where λ1 > λ2 >... > λ n

[0248] Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0249] (3-c) Establish a principal component model. The principle for selecting the number of principal components is: Incorporate the largest variance into the model space and leave the smallest variance to the noise space. The calculation formula is:

[0250] The information ratio included in each principal component:

[0251]

[0252] The information ratio included in the largest k principal components:

[0253]

[0254] (3-d) Retain the eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,..., p k ∈ R n×k , and the calculation formula for the score matrix is:

[0255] t i = Xp i , i = 1, 2,..., k (9)

[0256] Its essence is the projection of the X vector in the p i direction. The score matrix T = [t1, t2,..., t k ∈ R m×k .

[0257] In step (4), use the first two columns of the score matrix T to obtain xdat = [t1, t2] ∈ R m×2 Draw a confidence ellipse. In some embodiments, the steps of drawing a confidence ellipse using the matrix xdat are as follows:

[0258] (4-a) Calculate the covariance matrix of xdat, and invert the covariance matrix to obtain s ∈ R 2×2 ; The covariance matrix of xdat can be calculated with reference to the aforementioned formula (4);

[0259] (4-b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and subtract the corresponding average value from each column of xdat for centering processing to obtain xd ∈ R m×2 ;

[0260] (4-c) Calculate the formula xd × s · × xd, and sum each row of the obtained matrix to obtain rd ∈ R m×1 ;

[0261] (4-d) Draw a curve graph of xdat. The first two columns of the score matrix present an empirical distribution. According to the characteristics of its empirical distribution, calculate the percentile of the matrix rd. According to the confidence level C i of the confidence ellipse to be drawn, after sorting rd in ascending order, use the value corresponding to the C i -th position of rd as r ∈ R. Preferably, C i is 95%;

[0262] (4-e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in step (4-a) to obtain the eigenvalue matrix as The eigenvector matrix is The eigenvalues and corresponding eigenvectors of the matrix S can be calculated with reference to the aforementioned formulas (5) and (6);

[0263] (4-f) Using r in step (4-d) and D in step (4-e), the major axis and minor axis of the confidence ellipse can be obtained. The calculation formula is:

[0264]

[0265]

[0266] where a is the major axis and b is the minor axis;

[0267] (4-g) Using xm in step (4-b) as the center point of the ellipse, the confidence ellipse can be drawn according to the center point and the major axis and minor axis of the ellipse. The formula of the ellipse is:

[0268]

[0269] Where xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat.

[0270] The projections of data under different production modes in the confidence ellipse will be distributed in different regions. According to the performance indicators of the refining process, different labels can be set for historical data according to the benefit grading corresponding to historical operating conditions. The benefit values can be found in the historical statistical data. By projecting the first two columns of the score matrix of historical data, the confidence ellipse can be divided according to the labels carried by the data points in different regions.

[0271] In step (4), using the first two columns of the score matrix T, a two-dimensional confidence ellipse is drawn. The confidence level is preferably set at 95%. According to historical data and production requirements, key characterization index variables related to production mode definition and data collation can be completed, such as economic benefits, unit energy consumption, product yield, etc. Using the graphical method, the grades of production mode indicators can be divided according to requirements, and the data in the optimized state can be screened out. The samples corresponding to such data are projected on the two-dimensional confidence ellipse to obtain the optimized region under the defined production mode and marked with different colors.

[0272] In some embodiments, step (4) further includes dividing the region of the confidence ellipse according to the performance level. In some embodiments, the region division includes: flagging the drawn confidence ellipse by combining historical data with on-site process knowledge to achieve region division of data with different performance levels. The region division preferably includes: first, according to process knowledge and historical benefit statistical data, finding the benefit corresponding to the operating condition according to the time series, dividing into different levels according to the benefit level, and then finding the distribution of the historical data corresponding to each level in the ellipse and marking it with different colors. In some embodiments, the region division includes: setting different labels for historical data with different performance levels according to the performance indicators of the refining process, projecting the first two columns of the score matrix of historical data onto the confidence ellipse, and dividing the confidence ellipse according to the labels carried by the data points in different regions.

[0273] In step (5), Y is preprocessed with reference to formula (1) using the mean and variance obtained when preprocessing the training samples in step (2).

[0274] In step (6), the calculation formula of scorey is:

[0275] scorey = Ym × [p1, p2] (13)

[0276] In the formula, [p1, p2] ∈ R n×2 are the first two columns of the load matrix P.

[0277] In step (7), by judging whether the sample point is mapped within the ellipse, it is considered whether to regard this sample point as the normal operating condition of the device and include it in the historical data for re-modeling to improve the adaptability of the model; if the sample mapping point is within this ellipse, it is included in the database to update the model; otherwise, it indicates that there may be an abnormality in the production process of the device at this time, and the SPE contribution rate of the abnormal data can be calculated to trace the fault.

[0278] In step (7), each set of data in scorey can be substituted into the aforementioned formula (12) and compared with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse. If it is less than or equal to 1, it means that the set of data is mapped inside the ellipse.

[0279] In some embodiments, step (7) further comprises adding the collected data under normal operating conditions to the historical data for remodeling; preferably, while adding the collected data under normal operating conditions, the same amount of early historical data is eliminated.

[0280] In some embodiments, step (7) further includes calculating the SPE contribution rate of the abnormal data to perform fault tracing. The SPE contribution rate can be calculated using methods known in the art or the method provided by the present invention. In some embodiments, assuming that the abnormal data is x∈R 1×n , the calculation formula of its SPE contribution rate is:

[0281]

[0282] Among them, contspe(i) represents the SPE contribution corresponding to the i-th variable, ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the unit matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the variable with fault.

[0283] In some embodiments, step (7) further includes: for the samples mapped inside the ellipse, determine the current performance level according to their positions in combination with the regional division of the ellipse in step (4). If it is in a non-optimal performance level, adjust the key variables to make it transform to the desired better performance level; preferably, set the variable corresponding to the largest coefficient in the principal components as the key variable to be adjusted during optimization, and then determine the adjustment direction according to the correlation between the key variable and the position change of the projection point of the online data in the confidence ellipse to achieve the optimization of the production process; preferably, the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the oil refining process is adjusted, and it is observed whether the data moves to the desired area through real-time monitoring. In such embodiments, the method of the present invention can achieve the optimization of the oil refining process mode.

[0284] The optimization of the oil refining process mode can also be achieved by introducing the following step (8) on the basis of the oil refining process mode recognition method based on big data of the present invention:

[0285] (8) Obtain the original variables corresponding to the points in the confidence ellipse through the standard deviation and mean obtained in the principal component analysis process, so as to obtain the benefit value corresponding to this point; distinguish the distribution of high and low benefit values in the confidence ellipse; when the operating mode of the working condition is at a certain point in the confidence ellipse, adopt a path optimization algorithm to obtain the fastest moving trajectory from the current position to the optimal position, perform inverse transformation on the points corresponding to the trajectory, obtain the change mode of the operating conditions, and thus guide the optimization of the operation of the production device.

[0286] Therefore, the present invention also includes an oil refining process mode optimization method based on big data. The oil refining process mode optimization method based on big data includes the oil refining process mode recognition method described in any embodiment herein and step (8) described in any embodiment herein.

[0287] In the confidence ellipse, each point can obtain the corresponding original variables through the standard deviation and mean obtained in the principal component analysis process, and thus the benefit value corresponding to this point can be obtained. The distribution of high and low benefit values can be distinguished in the confidence ellipse. When the operating mode of the working condition is at a certain point in the confidence ellipse, adopt a path optimization algorithm to obtain the fastest moving trajectory from the current position to the optimal position, perform inverse transformation on the points corresponding to the trajectory, obtain the change mode of the operating conditions, and then guide the optimization of the operation of the production device.

[0288] In some embodiments, in step (8), for the current working condition point defen=(x, y)∈R in the confidence ellipse 1 ×2 , calculate the original variable X using formula (15) ori :

[0289] X ori = (defen × PC T ) · × std(X m ) + mean(X m ) (15)

[0290] where defen = (x, y) ∈ R 1×2 is the abscissa and ordinate of the points in the confidence ellipse, PC = [p1, p2] ∈ R J×2 are the first two columns of the load matrix during principal component analysis modeling, std(X m ) is the variance of the basic samples used in principal component analysis, mean(X m ) is the mean of the basic samples used in principal component analysis, m is the number of samples, and X ori is the original variable after inverse transformation corresponding to the position where defen = (x, y) ∈ R 1×2 in the ellipse.

[0291] The statistical information in the production process includes the feed rate, product output, etc. The simplified input-output benefit index can be obtained based on the price information. In some embodiments, in step (8), the input-output benefit index is calculated according to formula (16):

[0292]

[0293] where is the output of the i-th product, is the price of the i-th product, is the feed amount of the i-th raw material, is the price of the i-th raw material, and profit is the benefit index.

[0294] All the points in the confidence ellipse can obtain the original working conditions through the transformation of formula (15) and obtain a corresponding benefit value according to formula (16). At the same time, the point corresponding to the optimal benefit value can be found. The high-benefit area and low-benefit area can be divided in the confidence ellipse according to the benefit value.

[0295] In step (8), path optimization can adopt the algorithms known in the art (such as the conventional A* algorithm) or the improved A* algorithm provided by the present invention. The A* algorithm is an optimization algorithm that, after determining the starting point, finds the optimal path among all reachable paths from the starting point to the target point. The conventional A* algorithm determines the search direction and the next node to reach through the evaluation function f(x). The general form of the evaluation function is as follows:

[0296] f(x) = g(x) + h(x) (17-1)

[0297] Among them, g(x) is the cost function, which represents the actual cost required to reach the current node x from the starting node; h(x) is the heuristic function, which represents the estimated cost required to reach the target node from the current node x.

[0298] In a preferred embodiment, the improved A* algorithm of the present invention is used for path optimization, and the improved A* algorithm improves the evaluation function f(x):

[0299]

[0300] Where profit(x) represents the economic benefit corresponding to the selected node, g(x) is the cost function, which represents the actual cost required to reach the current node x from the starting node; h(x) is the heuristic function, which represents the estimated cost required to reach the target node from the current node x. The improved A* algorithm of the present invention is used for path optimization to maintain a high economic benefit during the optimization process. The improved A* algorithm of the present invention can obtain the optimal path from the optimized operating point to the high benefit point.

[0301] After obtaining the optimization path, the optimization operation can be performed based on the original value of the operating variable obtained by transformation using formula (15).

[0302] The beneficial effects of the present invention are as follows:

[0303] 1. Refining production contains a large amount of process data, such as temperature, pressure, liquid level, flow rate, and material properties. These variables will have a certain impact on the process. However, among the massive amount of data, some data have a more significant impact on the process, while others have a smaller impact. Using big data dimensionality reduction methods to reduce the dimensionality of massive data in the refining process and obtain key variables that can characterize the state of the refining process can effectively improve the efficiency of monitoring and optimization methods and reduce computational costs.

[0304] 2. During the modeling process, normal data collected in real time is continuously added to the modeling data, and early historical data is eliminated to keep the total amount of modeling data unchanged, thereby improving the adaptability of the model and allowing the model to adjust as the process runs, and the key variables monitored also change accordingly.

[0305] 3. The refining process is monitored in real time by drawing confidence ellipses. This visual approach allows operators to more intuitively understand the operating status of the process and perform related operations.

[0306] 4. The oil refining process under normal conditions can be divided into different production modes according to its process conditions and performance indicators, such as high economic efficiency mode or low economic efficiency mode. The data of different production modes will be mapped to different regions within the confidence ellipse. By mapping the historical data under these modes into the ellipse to divide the ellipse, during real-time monitoring, the production mode of the current oil refining process can be judged based on the position where the real-time collected data is mapped in the ellipse, and the production operation optimization of the oil refining device can be guided through the operation variable adjustment method obtained by path optimization and mode point inverse transformation.

[0307] 5. The present invention adopts the A* algorithm as the optimization algorithm, which can find an optimal path in the superior and inferior intervals, and at the same time, the optimal operating conditions can be inversely deduced from the points on the path.

[0308] The present invention will be specifically described below through embodiments. It is necessary to point out here that the following embodiments are only used to further illustrate the present invention and cannot be understood as limiting the protection scope of the present invention. Some non-essential improvements and adjustments made by those skilled in the art based on the content of the present invention still fall within the protection scope of the present invention.

[0309] Embodiment 1

[0310] In this embodiment, the oil refining process mode optimization method based on big data of the present invention is applied to the catalytic reforming (CCR) process. Figure 1 The process flow diagram of catalytic reforming is given. The catalytic reforming process consists of a pre-hydrogenation unit, a reforming unit, and a catalyst regeneration system. When the purpose is to produce aromatics, it also includes an aromatics extraction and distillation unit. The raw material after pretreatment enters the reforming section, is mixed with recycle hydrogen and heated, and then enters the reactor. The reactor consists of 3 to 4 in series, with a heating furnace arranged between them to compensate for the heat absorbed by the reaction. The material leaving the reactor enters the separator to separate the hydrogen recycle gas (the excess part is discharged), and the obtained liquid is stripped of light components by the stabilizer tower and used as reformed gasoline, which is a high-octane gasoline component, or sent to the aromatics extraction unit to produce aromatics.

[0311] The CCR model includes 51 input variables and 16 output variables. Different operating states of the CCR model can be obtained by controlling the input variables. Among them, the input variables that have a greater impact on the output results include reaction temperature, pressure, feed flow rate, etc. The performance indicators of the product are mainly judged by observing the mass flow rates of hydrogen, pure hydrogen, dry gas, liquefied gas, C5, C6, C7+, benzene, and the amount of aromatics in the output variables. Table 1 lists the operating variables involved in this embodiment.

[0312] Table 1: Operating Variables

[0313] Variable Description Unit 1 Reaction temperature of the first stage ℃ 2 Reaction temperature of the second stage ℃ 3 Reaction temperature of the third stage ℃ 4 Reaction temperature of the fourth stage ℃ 5 Circulating hydrogen volume <![CDATA[STD_m 3 / h]]> 6 Total feed of pre-hydrogenation tonne / h 7 Compressor pressure Mpa 8 Tray temperature of T201 ℃ 9 Reflux flow of T201 tonne / h 10 Tray temperature of T201 ℃ 11 Bottom temperature of T601 ℃ 12 Extracted xylene tonne / h 13 Reforming feed load tonne / h

[0314] The CCR process is monitored and optimized by using the big data-based refinery process model optimization method of the present invention, as Figure 2 shown, which includes the following steps:

[0315] 1. Collect the sample data Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where z i = [z 1i , z 2i ,..., z mi T represents the m samples of the i-th measurement variable.

[0316] 2. Preprocess the collected data to obtain a standardized data set X = [x1, x2,..., x n ∈ R m×n , and the calculation formula is:

[0317]

[0318] In the formula, μ is the mean value taken from the training data:

[0319]

[0320] σ is the standard deviation taken from the training data:

[0321]

[0322] 3. Perform dimensionality reduction on X obtained after standardization processing by using the principal component analysis method, which is specifically carried out according to the following steps:

[0323] a) Calculate the covariance matrix of matrix X, and the formula is:

[0324]

[0325] X is an m×n matrix, m is the number of training samples, and n is the number of features, so the obtained covariance matrix C is an n×n-dimensional matrix;

[0326] b) Calculate the eigenvalues λ i and eigenvectors p i , and sort them in descending order of eigenvalues. The calculation formula is:

[0327] det(C - Iλ) = 0 (5)

[0328] Cp i = λ i p i ​(6)

[0329] In formula (5), I is the identity matrix, and we get:

[0330] Eigenvalue matrix: where λ1 > λ2 >... > λ n

[0331] Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0332] c) Establish the principal component model. The principle for selecting the number of principal components is: Induce the largest variance into the model space and leave the smallest variance to the noise space. The calculation formula is:

[0333] The information ratio included in each principal component:

[0334]

[0335] The information ratio included in the largest k principal components:

[0336]

[0337] d) Retain the eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,..., p k ∈ R n ×k , and the calculation formula for the score matrix is:

[0338] t i = Xp i , i = 1, 2,..., k (9)

[0339] Its essence is the projection of the X vector in the p i direction, and the score matrix T = [t1, t2,..., t k ∈ R m×k .

[0340] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and the steps to specifically draw the confidence ellipse with a confidence level of 95% using this matrix are as follows:

[0341] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2 ;

[0342] b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and subtract the corresponding average value from each row's numerical value for centering processing to obtain xd ∈ R m×2 ;

[0343] c) Calculate the product \(x_d\times s\times x_d\), and sum each row of the resulting matrix to obtain \(r_d\in R\). m×1 ;

[0344] d) Plot the curve of \(x_{dat}\). According to the empirical distribution characteristics of \(x_{dat}\), calculate the percentile of the matrix \(r_d\). Since the confidence level of the confidence ellipse to be plotted is 95%, after sorting \(r_d\) in ascending order, take the value corresponding to the 95% position of \(r_d\) as \(r\in R\).

[0345] e) Calculate the eigenvalues and corresponding eigenvectors of the matrix \(s\) obtained in 4 - a). The calculation formula is similar to those in Equation (5) and Equation (6), and the eigenvalue matrix is The eigenvector matrix is

[0346] f) Using \(r\) in 4 - d) and \(D\) in 4 - e), the major axis and minor axis of the confidence ellipse can be obtained. The calculation formula is:

[0347]

[0348]

[0349] where \(a\) is the major axis and \(b\) is the minor axis.

[0350] g) Taking \(x_m\) in 4 - b) as the center point of the ellipse, the confidence ellipse can be plotted according to the center point, major axis and minor axis of the ellipse. The formula of the ellipse is:

[0351]

[0352] where \(x_{m1}\) is the average value of the first column of \(x_{dat}\), and \(x_{m2}\) is the average value of the second column of \(x_{dat}\). The plotted ellipse is as Figure 3 shown.

[0353] 5. The projections of data under different production modes in the confidence ellipse will be distributed in different regions. According to the performance indicators of the CCR process, different labels are set for the historical data according to the benefit grading corresponding to the historical working conditions. The benefit value can be found in the historical statistical data. Project the first two columns of the score matrix of the historical data. The confidence ellipse can be divided by the labels carried by the data points in different regions and distinguished by different colors as Figure 4 shown.

[0354] 6. After the CCR process runs under normal working conditions for a period of time, deviate the process from the normal state by adjusting the value of the compressor pressure in the operating variables, and collect data to obtain \(Y\in R\). N×n , and normalize the collected data using the mean and standard deviation calculated from the training data to obtain \(Y_m\in R\).N×n 。

[0355] 7. Multiply the first two columns of the load matrix P obtained by multiplying Ym by 3 - d), and the calculation formula is:

[0356] scorey = Ym × [p1, p2] (13)

[0357] 8. Substitute each set of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the set of data is mapped inside the ellipse.

[0358] 9. Assume that the fault data point is x ∈ R 1×n , and the calculation formula for its SPE contribution rate is:

[0359]

[0360] where contspe(i) represents the SPE contribution corresponding to the i-th variable, ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate sought is the variable with a fault. The SPE contribution rate is as Figure 5 shown. According to Figure 5 it can be determined that the main cause of the fault is mainly due to the seventh variable, i.e., the compressor pressure, which is consistent with the actual operation situation.

[0361] 10. The economic benefits of the catalytic reforming production process are respectively related to the liquefied gas flow rate (x1), hydrogen flow rate (x2), dry gas flow rate (x3), C5 (x4), C6 (x5), aromatics (x6), and reforming feed load (x7). The prices of these substances are shown in Table 2:

[0362] Table 2: Prices of various substances in the economic efficiency indicators of catalytic reforming

[0363] <![CDATA[substance / m 3 ·s -1 > Feed load Liquefied gas Hydrogen Dry gas C5 C6 Aromatics Price / yuan 3092 4153 8489.57 2505 3690 4298 4573

[0364] Therefore, the calculation formula for the economic benefits of the catalytic reforming production process can be expressed as follows:

[0365] profit = (4153x1 + 8489.57x2 + 2505x3 + 3690x4 + 4298x5 + 4573x6) - 3092x7 (16 - 1)

[0366] According to the historical samples, the benefit values as shown in Figure 6 can be plotted, and thus the high-benefit area and low-benefit area can be divided in the confidence ellipse, as shown in Figure 7 andFigure 8 as shown Figure 7 In it, the position where the green dots (light dots) are located represents the high-benefit area. Figure 8 In it, the position where the green dots (light dots) are located represents the low-benefit area.

[0367] For the current operating condition point defen = (x, y) ∈ R in the confidence ellipse 1×2 , the original variables can be obtained using the following formula:

[0368] X ori = (defen × PC T ) · × std(X m ) + mean(X m ) (15)

[0369] where defen = (x, y) ∈ R 1×2 is the abscissa and ordinate of the point in the confidence ellipse, PC = [p1, p2] ∈ R J×2 are the first two columns of the load matrix during principal component analysis modeling, std(X m ) is the variance of the basic sample used in principal component analysis, mean(X m ) is the mean of the basic sample used in principal component analysis, m is the number of samples, and X ori is the original variable after inverse transformation corresponding to the position of defen = (x, y) ∈ R in the ellipse 1×2 .

[0370] To obtain the optimal path from this variable to the high-benefit point, an improved A* algorithm is used to achieve path optimization. The evaluation function of this improved method is as follows:

[0371]

[0372] where g(x) is the cost function, which represents the actual cost required to reach the current node x from the starting node; h(x) is the heuristic function, representing the estimated cost required to reach the target node from the current node x, and profit(x) represents the economic benefit corresponding to the selected node, which can keep a high economic benefit during the optimization process. The pattern movement trajectory obtained after path optimization is as Figure 9 shown. The original value of the operating variable can be obtained by transformation using Equation (15). As shown in Table 3. By performing optimization operations using the operating conditions obtained by the above pattern optimization method, the benefit can be increased from 257,584 yuan / hour to 259,272 yuan / hour.

[0373] Table 3: Values of corresponding operating variables and benefit results during optimization

[0374]

[0375]

[0376] Example 2

[0377] In this example, according to the historical production data of a 1.8 million tons / year industrial catalytic cracking unit, as shown in Table 4, 88 device independent variable data points are selected as the model input variables and 20 device dependent variable data points are selected as the corresponding output variables of the model, and the historical production data is collected and preprocessed.

[0378] In this example, a data-driven method for identifying and optimizing the operating state of a catalytic cracking unit of the present invention is used to identify and optimize the operating state of the catalytic cracking unit, as Figure 10 shown, including the following steps:

[0379] 1. Collect the actual production historical data of the device, and the relevant variable names are shown in Table 4. Use a large number of data combinations under the normal operating state of the catalytic cracking unit to construct a sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where z i = [z 1i , z 2i ,..., z mi T represents m samples of the i-th measurement variable.

[0380] Table 4: Variable names of the catalytic device model

[0381]

[0382]

[0383]

[0384] 2. Preprocess the data of the sample set to obtain a standard data set X = [x1, x2,..., x n ∈ R m×n , and the calculation formula of the standard data set is:

[0385]

[0386] In the formula, μ is the mean of the training data, and its calculation formula is:

[0387]

[0388] σ is the standard deviation of the training data, and its calculation formula is:

[0389] ​

[0390] 3. The dimensionality reduction of the standard dataset X is performed using the principal component analysis method according to the following steps:

[0391] a) Calculate the covariance matrix of matrix X according to formula (4):

[0392]

[0393] where X is an m×n matrix, m is the number of training samples, and n is the number of features. Therefore, the covariance matrix C is an n×n-dimensional matrix;

[0394] b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C according to formula (5) and formula (6), and sort them in descending order of eigenvalues:

[0395] det(C - Iλ) = 0 (5)

[0396] Cp i = λ i p i (6)

[0397] In formula (5), I is the identity matrix;

[0398] Eigenvalue matrix: where λ1 > λ2 >... > λ n ;

[0399] Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0400] c) Establish a principal component model. The principle for selecting the number of principal components is: Induce the largest variance into the model space and leave the smallest variance to the noise space. The calculation formula is:

[0401] The information ratio included in each principal component:

[0402]

[0403] The information ratio included in the largest k principal components:

[0404]

[0405] d) Retain the k eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,... p, k ∈ R n×k , where the eigenvectors are sorted in descending order of eigenvalues. The calculation formula for the score matrix is:

[0406] t i= Xp i , i = 1, 2, ..., k (9)

[0407] t i essentially is the projection of the X vector in the p i direction, and the score matrix T = [t1, t2, ..., t k ∈ R m×k .

[0408] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and specifically draw a confidence ellipse with a confidence level of 95% using this matrix according to the following steps:

[0409] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2

[0410] b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and perform centering processing by subtracting the corresponding average value from each row's value to obtain xd ∈ R m×2 ;

[0411] c) Calculate xd × s · × xd, and sum each row of the resulting matrix to obtain rd ∈ R m×1 ;

[0412] d) Draw a curve graph of xdat, calculate the percentile of the matrix rd according to the characteristics of the empirical distribution presented by the first two columns of the score matrix. Since the confidence level of the confidence ellipse to be drawn is 95%, find the value corresponding to the 5% position of rd to obtain r ∈ R;

[0413] e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in 4 - a) with reference to equations (5) and (6) to obtain the eigenvalue matrix as and the eigenvector matrix as

[0414] f) Use r in 4 - d) and D in 4 - e) to calculate the major axis and minor axis of the confidence ellipse according to formulas (10) and (11):

[0415]

[0416]

[0417] where a is the major axis and b is the minor axis;

[0418] g) Use xm in 4 - b) as the center point of the ellipse, and draw the confidence ellipse according to the center point, major axis, and minor axis of the ellipse. The formula for the confidence ellipse is:

[0419]

[0420] Among them, xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat. The drawn ellipse is as Figure 11 shown.

[0421] 5. The projections of data under different production modes in the confidence ellipse will be distributed in different regions. As Figure 12 shown, according to the performance indicators of the fluid catalytic cracking process, different labels are set for historical data, and the optimization target variable - the economic benefit broken line is drawn, and the optimized operation area of the device is selected from it. The first two columns of the score matrix of historical data are projected, and the confidence ellipse is divided according to the labels carried by data points in different regions and distinguished by different colors. Figure 13 The confidence ellipse region division result of this embodiment is given. The region where the green (light color) points are located is the region corresponding to the production mode with the best economic benefit.

[0422] 6. After the fluid catalytic cracking unit operates under normal conditions for a period of time, the process is deviated from the normal state by adjusting the density at the lower part of the charring pot in the operating variables, and data Y∈R is collected N×n , and the collected data is normalized using the average value and standard deviation calculated from the training data to obtain Ym∈R N×n .

[0423] 7. Multiply the first two columns of the load matrix P obtained by multiplying Ym by (3 - d) according to formula (13) to obtain scorey:

[0424] scorey = Ym × [p1, p2] (13)

[0425] 8. Substitute each group of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that the group of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the group of data is mapped inside the ellipse;

[0426] 9. Let the fault data point be x∈R 1×n , and calculate the SPE contribution rate of the fault data point according to formula (14):

[0427]

[0428] Among them, contspe(i) represents the SPE contribution corresponding to the i-th variable, ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate sought is the faulty variable. The SPE contribution rate is asFigure 14 As shown. According to Figure 14 It can be determined that the main cause of the fault is mainly due to the 10th variable (i.e., the density at the lower part of the charring pot), which is consistent with the actual operation in step (6).

[0429] 10. According to the variable of the density at the lower part of the charring pot with the largest absolute value in the first column of the load matrix P, by adjusting the value corresponding to this variable, the production state of the fluid catalytic cracking process is adjusted, and it can be observed whether the data moves towards the desired area through real-time monitoring, so that the operating state of the device returns to the optimized operating state. The adjustment trajectory is as Figure 15 shown. Figure 15 Among them, the area where the green (light color) points are located corresponds to the optimized operating state.

[0430] Example 3

[0431] In this example, according to the production historical data of a 600,000-ton / year industrial sulfur recovery, as shown in Table 5, 75 device data points are selected, among which 55 are model input variables and 20 are output variables, and the production historical data is collected and preprocessed.

[0432] In this example, a data-driven method for identifying and optimizing the operating state of a sulfur recovery device of the present invention is used to identify and optimize the operating state of the sulfur recovery device, as Figure 16 shown, including the following steps:

[0433] 1. Collect the actual production historical data of the sulfur recovery device. Some independent variable names are shown in Table 5. Use a large number of data combinations in the normal operating state of the sulfur recovery device to construct a sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where z i = [z 1i , z 2i ,..., z mi T represents m samples of the i-th measurement variable.

[0434] Table 5: Variable names of sulfur device

[0435]

[0436]

[0437] 2. Standardize the sample set to obtain a standardized data set X = [x1, x2,..., x n ∈ R m×n , and the calculation formula for standardization is: ​

[0438]

[0439] Wherein, μ is the mean of the training data, and its calculation formula is:

[0440]

[0441] σ is the standard deviation of the training data, and its calculation formula is:

[0442]

[0443] 3. Perform dimensionality reduction on the standardized dataset X using the principal component analysis method according to the following steps:

[0444] a) Calculate the covariance matrix of matrix X according to formula (4):

[0445]

[0446] Where X is an m×n matrix, m is the number of training samples, and n is the number of features. Therefore, the covariance matrix C is an n×n-dimensional matrix;

[0447] b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C according to formula (5) and formula (6), and sort them in descending order of eigenvalues:

[0448] det(C - Iλ) = 0 (5)

[0449] Cp i = λ i p i (6)

[0450] In formula (5), I is the identity matrix;

[0451] Eigenvalue matrix: Where λ1 > λ2 >... > λ n ;

[0452] Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0453] c) Select the number of principal components k according to the principle of attributing the largest variance to the model space and leaving the smallest variance to the noise space. Its calculation formula is:

[0454] Information ratio included in each principal component:

[0455]

[0456] Information ratio included in the largest k principal components:

[0457]

[0458] d) Retain the k eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,... p, k ∈ R n×k , and the calculation formula for the score matrix is:

[0459] t i = Xp i , i = 1, 2,..., k (9)

[0460] t i essentially represents the projection of the X vector in the direction of p i , and the score matrix T = [t1, t2,..., t k ∈ R m×k .

[0461] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and use this matrix to specifically draw a confidence ellipse with a confidence level of 95% according to the following steps:

[0462] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2

[0463] b) Calculate the mean value of each column of xdat to obtain xm ∈ R 1×2 , and perform centering processing by subtracting the corresponding mean value from each row of the numerical values to obtain xd ∈ R m×2 ;

[0464] c) Calculate xd × s · × xd, and sum each row of the obtained matrix to obtain rd ∈ R m×1 ;

[0465] d) Draw a curve graph of xdat. According to the characteristics of the empirical distribution of xdat, calculate the percentile of the matrix rd. According to the confidence level of the confidence ellipse to be drawn being 95%, after sorting rd in ascending order, use the value corresponding to the 95% position of rd as r ∈ R;

[0466] e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in 4-a) with reference to equations (5) and (6) to obtain the eigenvalue matrix as the eigenvector matrix is

[0467] f) Use r in 4-d) and D in 4-e) to calculate the major axis and minor axis of the confidence ellipse according to formulas (10) and (11):

[0468]

[0469]

[0470] Among them, a is the major axis and b is the minor axis;

[0471] g) Taking xm in 4 - b) as the center point of the ellipse, draw a confidence ellipse according to the center point, the major axis and the minor axis of the ellipse. The formula of the confidence ellipse is:

[0472]

[0473] Among them, xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat. The drawn ellipse is as Figure 17 shown.

[0474] 5. The projections of data under different production modes in the confidence ellipse will be distributed in different regions. As Figure 18 shown, according to the device characteristics of the sulfur recovery process, after setting the optimization target of the maximum sulfur recovery rate for historical data and labeling it, draw a broken line of the corresponding optimization target variable, and select the optimized operation region of the device (sulfur recovery rate > 99.5%). Project the first two columns of the score matrix of historical data, divide the confidence ellipse according to the labels carried by data points in different regions, and distinguish them with different colors. Figure 19 The confidence ellipse region division result of this embodiment is given. The region where the green (light - colored) points are located is the region corresponding to the production mode pursuing the highest sulfur recovery rate.

[0475] 6. After the sulfur recovery device operates under normal conditions for a period of time, deviate the process from the normal state by adjusting the temperature of the sulfur - making furnace in the operating variables, and collect data to obtain Y ∈ R N×n , and normalize the collected data using the average value and standard deviation of the training data obtained in step (2) to obtain Ym ∈ R N×n .

[0476] 7. Multiply Ym by the first two columns of the load matrix P obtained in 3 - d) according to formula (13) to obtain scorey:

[0477] scorey = Ym × [p1, p2] (13)

[0478] 8. Substitute each group of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that the group of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the group of data is mapped inside the ellipse;

[0479] 9. Let the fault data point be x ∈ R 1×n, calculate the SPE contribution rate of the fault data point according to formula (14):

[0480]

[0481] Among them, contspe(i) represents the SPE contribution corresponding to the i-th variable, and ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the faulty variable. The SPE contribution rate is as Figure 20 shown. According to Figure 20 it can be determined that the main cause of the fault is mainly due to the 11th variable (i.e., the temperature of the sulfur-making furnace), which is consistent with the actual operation situation.

[0482] 10. According to the sulfur-making furnace temperature variable with the largest absolute value in the first column of the load matrix P, adjust the production state of the sulfur recovery process by adjusting the value corresponding to this variable, and through real-time monitoring, observe whether the data moves to the desired area, so that the operating state of the sulfur recovery device returns to the optimized operating state. The adjustment trajectory is as Figure 21 shown. Figure 21 Among them, the area where the green (light color) points are located corresponds to the optimized operating state.

[0483] Example 4

[0484] In this example, a residue hydrotreating unit with a capacity of 1.7 million tons / year is taken as the object for identification and optimization. Historical production data from May 1, 2018 to May 1, 2019 are collected, and abnormal fluctuation data are excluded. After screening and processing, the data are sorted into 3925 groups of 85 variables. Table 6 lists the variables involved in this example.

[0485] Table 6: Variables involved in Example 4

[0486]

[0487]

[0488] Use the data-driven method for identifying and optimizing the operating state of the residue hydrotreating unit of the present invention to identify and optimize the operating state of the residue hydrotreating unit, as Figure 22 shown, including the following steps:

[0489] 1. Collect the actual production historical data of the residue hydrotreating unit. Use a large number of data combinations in the normal operating state of the residue hydrotreating unit to construct a sample set Z = [z1, z2,..., z i ,..., z n ∈ Rm×n , where z i = [z 1i , z 2i ,..., z mi T represents m samples of the i-th measurement variable.

[0490] 2. Standardize the data of the sample set according to formula (1) to obtain a standardized data set X = [x1, x2,..., x n ∈ R m×n :

[0491]

[0492] In the formula, μ is the mean of the training data, and its calculation formula is:

[0493]

[0494] σ is the standard deviation of the training data, and its calculation formula is:

[0495]

[0496] 3. Dimensionality reduction of the standardized data set X is performed using the principal component analysis method according to the following steps:

[0497] a) Calculate the covariance matrix of matrix X according to formula (4):

[0498]

[0499] Among them, X is an m×n matrix, m is the number of training samples, n is the number of features, T represents transpose, so the covariance matrix C is an n×n-dimensional matrix;

[0500] b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C according to formula (5) and formula (6), and sort them in descending order of eigenvalues:

[0501] det(C - Iλ) = 0 (5)

[0502] Cp i = λ i p i (6)

[0503] In formula (5), I is the identity matrix;

[0504] Eigenvalue matrix: where λ1 > λ2 >... > λ n ; ​

[0505] Eigenvector: V = [p1, p2, p3, ..., p n ∈ R n×n ;

[0506] c) Select the number of principal components k according to the principle of inducing the largest variance into the model space and leaving the smallest variance to the noise space. Its calculation formula is:

[0507] The information ratio included in each principal component:

[0508]

[0509] The information ratio included in the k largest principal components:

[0510]

[0511] d) Retain the k eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3, ... p, k ∈ R n×k , and the calculation formula of the score matrix is:

[0512] t i = Xp i , i = 1, 2, ..., k (9)

[0513] Its essence is the projection of the X vector in the p i direction, and the score matrix T = [t1, t2, ..., t k ∈ R m×k .

[0514] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and the steps to specifically draw a confidence ellipse with a confidence level of 95% using this matrix are as follows:

[0515] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2 ;

[0516] b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and perform centering processing by subtracting the corresponding average value from each row's value to obtain xd ∈ R m×2 ;

[0517] c) Calculate xd × s · × xd, and sum each row of the obtained matrix to obtain rd ∈ R m×1 ;

[0518] d) Plot the curve of xdat, calculate the percentile of matrix rd according to the empirical distribution characteristics of xdat. Based on the 95% confidence level of the confidence ellipse to be plotted, sort rd in ascending order, and take the value corresponding to the 95% position of rd as r ∈ R;

[0519] e) Refer to Equation (5) and Equation (6) to calculate the eigenvalues and corresponding eigenvectors of matrix s obtained in 4-a), and obtain the eigenvalue matrix as The eigenvector matrix is

[0520] f) Use r in 4-d) and D in 4-e) to calculate the major axis and minor axis of the confidence ellipse according to Formula (10) and Formula (11):

[0521]

[0522]

[0523] where a is the major axis and b is the minor axis;

[0524] g) Take xm in 4-b) as the center point of the ellipse, and draw the confidence ellipse according to the center point, the major axis and the minor axis of the ellipse. The formula of the confidence ellipse is:

[0525]

[0526] where xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat. The drawn ellipse is as Figure 23 shown.

[0527] 5. The projections of data under different production modes in the confidence ellipse will be distributed in different regions. As Figure 24 shown, according to different performance indicators of the residue hydrotreating process, the historical data is set into three modes with the maximum heavy oil hydrogenation yield as the target, the maximum plant benefit as the target, and the minimum volume flow of fresh hydrogen from outside the plant as the target. By projecting the first two columns of the score matrix of the historical data, the confidence ellipse is divided by the labels carried by the data points in different regions. Figure 25 The green (light color) points in [Figure] represent the projection positions of the optimized region data corresponding to the production samples focusing on the heavy oil hydrogenation yield on the confidence ellipse, and these positions are the optimized regions corresponding to the mode with the maximum heavy oil hydrogenation yield as the target.

[0528] 6. After the residue hydrotreating unit operates under normal conditions for a period of time, adjust the operating variable, the cold slag flow rate in the tank farm, to make the process deviate from the normal state, and collect data to obtain Y ∈ R N×n, the collected data is normalized using the mean and standard deviation calculated from the training data to obtain Ym ∈ R N×n .

[0529] 7. Multiply the first two columns of the load matrix P obtained by multiplying Ym by 3 - d according to formula (13) to get scorey:

[0530] scorey = Ym × [p1, p2] (13)

[0531] 8. Substitute each set of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the set of data is mapped inside the ellipse;

[0532] 9. Let the fault data point be x ∈ R 1×n , and calculate the SPE contribution rate of the fault data point according to formula (14):

[0533]

[0534] where contspe(i) represents the SPE contribution corresponding to the i - th variable, ξ i represents the i - th column of the n - dimensional identity matrix, T represents transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the variable where the fault occurs. The SPE contribution rate is as Figure 26 shown. According to Figure 26 it can be determined that the main cause of the fault is mainly due to the 20th variable (i.e., the cold slag flow rate in the tank farm), which is consistent with the actual operation situation.

[0535] 10. According to the cold slag flow rate variable in the tank farm with the largest absolute value in the first column of the load matrix P, by adjusting the value corresponding to this variable, the production state of the residue hydrotreating process can be adjusted, and through real - time monitoring, it can be observed whether the data moves to the desired area, so that the device operation state returns to the optimized operation state. The adjustment trajectory is shown in the light - colored points in Figure 27 .

[0536] Example 5

[0537] In this example, the big - data - based atmospheric and vacuum process pattern monitoring and optimization method is applied to the actual atmospheric and vacuum process of a refinery. Figure 28The process flow diagram of the atmospheric and vacuum distillation is given. The atmospheric and vacuum distillation process is divided into five parts, namely the electro - desalting part, the pre - distillation part, the atmospheric distillation part, the vacuum distillation part, and the light hydrocarbon recovery part. Among them, the most important are the three parts of pre - distillation, atmospheric distillation, and vacuum distillation. The main purpose of crude oil pre - distillation is to extract some light fractions, share a certain device pressure for the subsequent atmospheric and vacuum distillation unit, and reduce the energy consumption of the whole atmospheric and vacuum distillation unit. The main function of the atmospheric distillation unit is to extract fractions such as naphtha, kerosene, and diesel with lower boiling points. The relatively heavy atmospheric residue distilled from the bottom of the atmospheric column will be sent to the vacuum column for vacuum distillation to separate raw materials for secondary processing such as lubricating oil, wax oil, and asphalt. The vacuum residue distilled from the bottom of the vacuum column will be sent to units such as catalytic cracking and delayed coking for further processing.

[0538] The atmospheric and vacuum model includes 60 input variables and 20 output variables. Among them, the input variables include 52 sets of device operation variables and 8 sets of raw material property variables. By controlling the input variables, CDU models in different operating states can be obtained. Among them, the input variables that have a greater impact on the output results are reaction temperature, pressure, feed properties, etc. The performance indicators of the products are mainly judged by observing the output variables. For example, in the production mode with the total draw ratio as the target, the total draw ratio is mainly composed of the combined yields of naphtha, the first atmospheric side - stream, diesel, wax oil, and the fourth vacuum side - stream. The calculation of these yields can be carried out through the gas mass flow rate, the flow rate of the first overhead oil out of the device, the flow rate of the atmospheric overhead oil out of the device, the flow rate of the first atmospheric side - stream extraction, the flow rate of the second atmospheric side - stream extraction, the flow rate of the third atmospheric side - stream extraction, the flow rate of the first vacuum side - stream extraction, the flow rate of cold wax oil extraction, the flow rate of wax oil sent to the first catalytic cracking unit, the flow rate of wax oil sent to the second catalytic cracking unit, the flow rate of the fourth vacuum side - stream extraction, and the flow rate of residue extraction in the output variables. Table 7 lists all the operation variables and raw material property variables involved in this embodiment.

[0539] Table 7: Operation variables involved in Example 5

[0540]

[0541]

[0542] Using the method for monitoring and optimizing the atmospheric and vacuum process mode based on big data of the present invention to monitor and optimize the atmospheric and vacuum process, as Figure 29 shown, it includes the following steps:

[0543] 1. Collect sample data Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where z i = [z 1i , z 2i ,..., z mi ​T m samples representing the i-th measurement variable;

[0544] 2. Preprocess the collected data according to formula (1) to obtain a standard data set X = [x1, x2,..., x n ∈ R m×n :

[0545]

[0546] where μ is the mean of the collected data, and its calculation formula is:

[0547]

[0548] σ is the standard deviation of the collected data, and its calculation formula is:

[0549]

[0550] 3. Perform dimensionality reduction on the standardized data set X using the principal component analysis method according to the following steps:

[0551] a) Calculate the covariance matrix of matrix X according to formula (4):

[0552]

[0553] where X is an m×n matrix, m is the number of training samples, n is the number of features, T represents transpose, so the covariance matrix C is an n×n-dimensional matrix;

[0554] b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C according to formula (5) and formula (6), and sort them in descending order of eigenvalues:

[0555] det(C - Iλ) = 0 (5)

[0556] Cp i = λ i p i (6)

[0557] In formula (5), I is the identity matrix;

[0558] Obtain the eigenvalue matrix: where λ1 > λ2 >... > λ n ;

[0559] Obtain the eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0560] c) Principle for selecting the number of principal components: Incorporate the largest variances into the model space and leave the smallest variances to the noise space. The calculation formula is as follows:

[0561] Each principal component includes an information ratio:

[0562]

[0563] Information ratio included in the largest k principal components:

[0564]

[0565] d) Retain the k eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,... p, k ∈ R n×k , and the calculation formula for the score matrix is:

[0566] t i = Xp i , i = 1, 2,..., k (9)

[0567] Its essence is the projection of the X vector in the p i direction. The score matrix T = [t1, t2,..., t k ∈ R m×k .

[0568] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 . The steps to specifically draw a confidence ellipse with a confidence level of 95% using this matrix are as follows:

[0569] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2 ;

[0570] b) Calculate the average value of each column of xdat to obtain xm ∈ R 1×2 , and perform centering processing by subtracting the corresponding average value from each row's numerical value to obtain xd ∈ R m×2 ;

[0571] c) Calculate xd × s · × xd, and sum each row of the resulting matrix to obtain rd ∈ R m×1 ;

[0572] d) Draw a curve graph of xdat. According to the characteristics of the empirical distribution of xdat, calculate the percentile of the matrix rd. Based on the confidence level of 95% for the confidence ellipse to be drawn, after sorting rd in ascending order, take the value corresponding to the 95% position of rd as r ∈ R;

[0573] e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in 4-a) with reference to Equation (5) and Equation (6), and the eigenvalue matrix is The eigenvector matrix is

[0574] f) Using r in 4-d) and D in 4-e), calculate the major axis and minor axis of the confidence ellipse according to Equation (10) and Equation (11):

[0575]

[0576]

[0577] where a is the major axis and b is the minor axis;

[0578] g) Taking xm in 4-b) as the center point of the ellipse, draw a confidence ellipse according to the center point, major axis and minor axis of the ellipse. The formula of the confidence ellipse is:

[0579]

[0580] where xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat. The drawn ellipse is as Figure 30 shown.

[0581] 5. The projections of data under different production modes in the confidence ellipse will be distributed in different regions. According to the performance indicators of the CDU process, different labels are set for the historical data. By projecting the first two columns of the score matrix of the historical data, the confidence ellipse can be divided by the labels carried by the data points in different regions and framed in different colors. Figure 31 、 Figure 32 and Figure 33 show the division intervals corresponding to the mode with energy consumption as the optimization target, the mode with light ends yield optimization as the target, and the mode with total draw as the optimization target respectively. The boxed areas in the figure correspond to the 50% data points with the lowest energy consumption, the 50% data points with the highest light ends yield, and the 50% data points with the highest total draw respectively.

[0582] 6. After the CDU process operates under normal conditions for a period of time, deviate the process from the normal state by adjusting the value of the operating variable, the initial distillation tower top pressure, and collect data to obtain Y ∈ R N×n , and normalize the collected data using the average value and standard deviation calculated in step (2) to obtain Ym ∈ R N×n .

[0583] 7. Multiply Ym by the first two columns of the load matrix P obtained in 3-d) according to Equation (13) to obtain scorey:

[0584] scorey = Ym × [p1, p2] (13)

[0585] 8. Substitute each set of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that this set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this set of data is mapped inside the ellipse, as Figure 34 shown.

[0586] 9. Let the fault data point be x ∈ R 1×n , and calculate the SPE contribution rate of the fault data point according to formula (14):

[0587]

[0588] where contspe(i) represents the SPE contribution corresponding to the i-th variable, ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the variable with a fault. The SPE contribution rate is as Figure 35 shown. According to Figure 35 it can be determined that the main cause of the fault is mainly due to the 3rd variable (i.e., the pressure at the top of the pre-fractionator), which is consistent with the actual operation situation (where variable 1 is the raw material property variable and is difficult to directly adjust in actual operation).

[0589] Example 6

[0590] In this example, the big data-based process pattern recognition and optimization method for hydrocracking is applied to the hydrocracking (HCR) process of a domestic refinery. Figure 36 The process flow chart of the hydrocracking process is given. The hydrocracking process consists of hydrotreating, hydrocracking, and fractionation sections. In the single-stage series process, the feedstock oil and recycle hydrogen reach the hydrocracking reaction conditions respectively and then are mixed and enter the hydrotreating reactor. Under the action of the hydrotreating catalyst, desulfurization, denitrification, deoxidation, and partial dearomatization reactions are carried out. The refined reaction product is adjusted to the temperature required for the cracking reaction by injecting cold hydrogen and then enters the hydrocracking reactor for cracking reaction to convert heavy distillate oil into light distillate oil. The reaction product is separated into gas, oil, and water phases in the cold high-pressure separator, and after desulfurization in the recycle hydrogen desulfurization tower, it is sent to the low-pressure separator to flash off the low-boiling gas. The low-boiling oil enters the debutanizer after heat exchange. The overhead separates out the sulfur-containing gas and light hydrocarbons, and the bottom oil of the tower enters the product separation tower for product separation to separate products such as light and heavy naphtha, jet fuel, diesel, and tail oil.

[0591] The HCR model includes 32 input variables and 24 output variables. By controlling the input variables, HCR models in different operating states can be obtained. Among them, the input variables that have a greater impact on the output results include hydrogen-oil ratio, reaction temperature, pressure, feed flow rate, and feed properties, etc. The performance indicators of the products are mainly judged by observing the amounts of dry gas, liquefied gas, light ends, light naphtha, heavy naphtha, aviation kerosene, diesel, and tail oil leaving the unit in the output variables. Table 8 gives the operating variables of the HCR model in this embodiment.

[0592] Table 8: Operating Variables of the HCR Model

[0593] Variable Description Unit Variable Description Unit 1 Flow rate of light wax oil in the tank farm tonne / h 17 Sulfur content of the feedstock oil ％ 2 Total flow rate of atmospheric and vacuum wax oil tonne / h 18 Nitrogen content of the feedstock oil ppmwt 3 Accumulative flow rate of catalytic diesel tonne / h 19 Basic nitrogen of the feedstock oil ppmwt 4 Mass flow rate of recycled tail oil tonne / h 20 Initial boiling point of the feedstock oil C 5 Flow rate of fresh hydrogen STD_m3 / h 21 10% recovery temperature of the feedstock oil C 6 Flow rate of circulating hydrogen for purging STD_m3 / h 22 50% recovery temperature of the feedstock oil C 7 Hydrogen-oil ratio 23 90% recovery temperature of the feedstock oil C 8 Inlet temperature of the first bed of R101 C 24 Final boiling point of the feedstock oil C 9 Inlet temperature of the second bed of R101 C 25 Density of light diesel oil (20℃) 10 Inlet temperature of the third bed of R101 C 26 Sulfur content of light diesel oil ％ 11 Inlet temperature of the first bed of R102 C 27 Initial boiling point of light diesel oil C 12 Inlet temperature of the second bed of R102 C 28 10% recovery temperature of light diesel oil C 13 Inlet temperature of the third bed of R102 C 29 50% recovery temperature of light diesel oil C 14 Inlet temperature of the fourth bed of R102 C 30 90% recovery temperature of light diesel oil C 15 Top pressure of R101 MPa 31 95% recovery temperature of light diesel oil C 16 Density of the feedstock oil (20℃) 32 Fresh hydrogen ％

[0594] Using the method for pattern recognition and optimization of the hydrocracking process based on big data of the present invention to perform pattern recognition and optimization on the hydrocracking process, as Figure 37 shown, it includes the following steps:

[0595] 1. Collect sample data Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where z i = [z 1i , z 2i ,..., z mi T represents m samples of the i-th measurement variable;

[0596] 2. Preprocess the collected data to obtain a standard data set X = [x1, x2,..., x n ∈ R m×n , and the calculation formula is:

[0597]

[0598] In the formula, μ is the mean of the collected data, and its calculation formula is:

[0599]

[0600] σ is the standard deviation of the collected data, and its calculation formula is:

[0601]

[0602] 3. Use the principal component analysis method to reduce the dimension of the standard data set X, and proceed as follows:

[0603] a) Calculate the covariance matrix of matrix X according to formula (4):

[0604] ​

[0605] Among them, X is an m×n matrix, where m is the number of training samples and n is the number of features. T represents the transpose. Therefore, the covariance matrix C is an n×n-dimensional matrix;

[0606] b) Calculate the eigenvalues λ i and eigenvectors p i of the covariance matrix C according to formulas (5) and (6), and sort them in descending order of eigenvalues:

[0607] det(C - Iλ) = 0 (5)

[0608] Cp i = λ i p i (6)

[0609] In formula (5), I is the identity matrix;

[0610] Obtain the eigenvalue matrix: where λ1 > λ2 >... > λ n ;

[0611] Obtain the eigenvectors: V = [p1, p2, p3,..., p n ∈ R n×n ;

[0612] c) Principle for selecting the number of principal components: Incorporate the largest variances into the model space and leave the smallest variances to the noise space. The calculation formula is:

[0613] The information ratio included in each principal component:

[0614]

[0615] The information ratio included in the largest k principal components:

[0616]

[0617] d) Retain the k eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,... p, k ∈ R n×k , and the calculation formula for the score matrix is:

[0618] t i = Xp i , i = 1, 2,..., k (9)

[0619] Its essence is the projection of the X vector in the p i direction, and the score matrix T = [t1, t2,..., t k ∈ Rm×k 。

[0620] 4. Let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and the steps to specifically draw a confidence ellipse with a confidence level of 95% using this matrix are as follows:

[0621] a) Calculate the covariance matrix of xdat with reference to formula (4), and invert the covariance matrix to obtain s ∈ R 2×2 ;

[0622] b) Calculate the mean value of each column of xdat to obtain xm ∈ R 1×2 , and centralize the values of each row by subtracting the corresponding mean value to obtain xd ∈ R m×2 ;

[0623] c) Calculate xd × s · × xd, and sum the values of each row of the resulting matrix to obtain rd ∈ R m×1 ;

[0624] d) Draw a curve graph of xdat. According to the characteristics of the empirical distribution of xdat, calculate the percentile of the matrix rd. According to the confidence level of 95% of the confidence ellipse to be drawn, after sorting rd in ascending order, take the value corresponding to the 95% position of rd as r ∈ R;

[0625] e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in 4 - a) with reference to formulas (5) and (6) to obtain the eigenvalue matrix as The eigenvector matrix is

[0626] f) Use r in 4 - d) and D in 4 - e) to calculate the major axis and minor axis of the confidence ellipse according to formulas (10) and (11):

[0627]

[0628]

[0629] where a is the major axis and b is the minor axis;

[0630] g) Take xm in 4 - b) as the center point of the ellipse, and draw a confidence ellipse according to the center point, major axis, and minor axis of the ellipse. The formula of the confidence ellipse is:

[0631]

[0632] where xm1 is the mean value of the first column of xdat, xm2 is the mean value of the second column of xdat, and the drawn ellipse is as Figure 38 shown.

[0633] 5. Data under different production modes will be projected onto different regions within the confidence ellipse. According to the performance metrics of the HCR process, different labels are set for historical data. By projecting the first two columns of the score matrix of historical data, the confidence ellipse can be divided based on the labels carried by data points in different regions and framed in different colors. Figure 39 and Figure 40 and Figure 41 respectively show the division intervals corresponding to the mode with the total liquid recovery as the optimization target, the mode with the middle oil yield as the target, and the mode with the value increment as the optimization target. Among them, the large red box represents the top 10% of the data points selected according to the corresponding target, and the small red box represents the top 5% of the data points.

[0634] 6. After the HCR process operates under normal conditions for a period of time, the process is deviated from the normal state by adjusting the value of the operating variable hydrogen-oil ratio, and data is collected to obtain Y ∈ R N×n , and the collected data is normalized using the mean and standard deviation calculated in step (2) to obtain Ym ∈ R N×n .

[0635] 7. Multiply the first two columns of the load matrix P obtained in 3 - d) by Ym to get scorey, and the calculation formula is:

[0636] scorey = Ym × [p1, p2] (13)

[0637] 8. Substitute each set of data in scorey into formula (12) and compare it with the value 1. If it is greater than 1, it means that the set of data is mapped outside the ellipse; if it is less than or equal to 1, it means that the set of data is mapped inside the ellipse, as Figure 42 shown.

[0638] 9. Let the faulty data point be x ∈ R 1×n , and calculate the SPE contribution rate of the faulty data point according to formula (14):

[0639]

[0640] where contspe(i) represents the SPE contribution corresponding to the i-th variable, ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from training samples, n is the number of variables, and the variable with the largest SPE contribution rate sought is the faulty variable. The SPE contribution rate is as Figure 43 shown. According to Figure 43 it can be determined that the main cause of the fault is mainly due to the first variable (i.e., the hydrogen-oil ratio), which is consistent with the actual operation situation.

Claims

1. A method for pattern recognition and optimization of the oil refining process based on big data, characterized in that, The method includes the following steps: (1) Compose the historical data collected during the oil refining process into a training sample set Z = [z1, z2,..., z i ,..., z n ∈ R m×n , where m is the number of samples in the sample set and n is the number of variables in the sample set; (2) Preprocess the training sample set to obtain standardized data \(X = [x_1, x_2, \cdots, x\) with a mean of 0 and a variance of 1. n \(\in R\) m×n ; (3) Apply the principal component analysis method to X to reduce its dimension from n to k, and obtain the score matrix T ∈ R m×k and the load matrix P ∈ R n ×k ; (4) Using the first two columns of the score matrix T, draw a two-dimensional confidence ellipse; (5) Collect new online real-time data Y ∈ R N×n , and preprocess Y using the sample mean and sample variance obtained when preprocessing the training sample set in step (2) to obtain standardized data Ym ∈ R N×n ; (6) Multiply Ym by the first two columns of the load matrix P obtained in step (3) to obtain the first two column score matrix scorey ∈ R of Ym according to the training samples N×2 ; (7) Using the first column of scorey as the data for the x-axis and the second column of scorey as the data for the y-axis, map scorey to the confidence ellipse drawn in step (4); if the sample point is mapped inside the ellipse, it indicates that the operating condition of the refinery process at this time is normal; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the refinery process at this time; (8) Using the standard deviation and mean obtained during the principal component analysis process to obtain the original variables corresponding to the points in the confidence ellipse, thereby obtaining the benefit value corresponding to this point; distinguish the distribution of high and low benefit values in the confidence ellipse; when the operating condition mode is at a certain point in the confidence ellipse, use the path optimization algorithm to obtain the fastest moving trajectory from the current position to the optimal position, and perform inverse transformation on the points corresponding to the trajectory to obtain the change method of the operating conditions, thereby guiding the optimization of the production device operation; wherein, the path optimization algorithm uses an improved A* algorithm, and determines the search direction and the next node to reach through the improved evaluation function shown in formula (17): where g(x) is the cost function, representing the actual cost required to reach the current node x from the starting node; h(x) is the heuristic function, representing the estimated cost required to reach the target node from the current node x; profit(x) represents the economic benefit corresponding to the selected node.

2. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, wherein, In step (2), the preprocessing method uses the Z-score normalization method, and the calculation formula is: where \(Z = [z_1, z_2, \cdots, z m \) is the training data matrix, \(X\) represents the standardized data matrix, \(\mu\) is the mean of the training data, \(\sigma\) is the standard deviation of the training data, and the calculation formulas for \(\mu\) and \(\sigma\) are:

3. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that In step (3), perform dimensionality reduction processing on X obtained through preprocessing using the principal component analysis method, and the specific steps are as follows: (3-a) Calculate the covariance matrix of matrix X, and the calculation formula of the covariance matrix is: X is an m×n matrix, m is the number of training samples, n is the number of features, and T represents transpose, so the obtained covariance matrix C is an n×n-dimensional matrix; (3-b) Calculate the eigenvalues λ of the covariance matrix C i and the eigenvectors p i , and sort them in descending order of eigenvalues. The calculation formula is: det(C - Iλ) = 0 (5) Cp i = λ i p i (6) In formula (5), I is the identity matrix, and we get: Eigenvalue matrix: where λ1 > λ2 >... > λ n Eigenvector: V = [p1, p2, p3,..., p n ∈ R n×n ; (3-c) Establish a principal component model, and the calculation formula is: Calculate the information ratio included in each principal component: Calculate the information ratio included in the largest k principal components: (3-d) Retain the eigenvectors corresponding to the k largest eigenvalues to obtain the load matrix P = [p1, p2, p3,..., p k ∈ R n ×k , and the calculation formula for the score matrix is: t i = Xp i , i = 1, 2, ..., k (9) t i is the projection of the X vector in the p i direction, and the score matrix T = [t1, t2,..., t k ∈ R m×k .

4. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, wherein In step (4), let the first two columns of the score matrix T be xdat = [t1, t2] ∈ R m×2 , and the steps of specifically drawing the confidence ellipse using the matrix xdat are as follows: (4-a) Calculate the covariance matrix of xdat and invert the covariance matrix to obtain s ∈ R 2×2 ; (4-b) Calculate the mean value of each column of xdat to obtain xm ∈ R 1×2 , and centralize the values of each column of xdat by subtracting the corresponding mean value to obtain xd ∈ R m×2 ; (4-c) Calculate the formula \(x_d\times s\times x_d\), and sum each row of the resulting matrix to obtain \(r_d\in R\). m×1 ; (4-d) Plot the curve of xdat. According to the characteristic that xdat presents an empirical distribution, calculate the percentile of the matrix rd. According to the confidence level C of the confidence ellipse to be plotted i , after sorting rd in ascending order, obtain the value r ∈ R corresponding to the C i -th position of rd. Preferably, C i is 95%; (4 - e) Calculate the eigenvalues and corresponding eigenvectors of the matrix s obtained in step (4 - a), and obtain the eigenvalue matrix as The eigenvector matrix is (4-f) Using r in step (4-d) and D in step (4-e), the major axis and minor axis of the confidence ellipse can be obtained, and the calculation formula is: where a is the major axis and b is the minor axis; (4-g) Using xm in step (4-b) as the center point of the ellipse, draw a confidence ellipse according to the center point, the major axis and the minor axis of the ellipse, and the formula of the ellipse is: where xm1 is the average value of the first column of xdat, and xm2 is the average value of the second column of xdat.

5. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that, In step (6), the calculation formula of scorey is: scorey = Ym×[p1, p2] (13) where [p1, p2] ∈ R n×2 are the first two columns of the load matrix P.

6. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 4, characterized in that In step (7), substitute each group of data in scorey into the following formula and compare it with the value 1. If it is greater than 1, it means that this group of data is mapped outside the ellipse; if it is less than or equal to 1, it means that this group of data is mapped inside the ellipse; 7. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, wherein Step (7) also includes adding the data collected under normal operating conditions to the historical data for re-modeling.

8. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that, Step (7) further includes calculating the SPE contribution rate for abnormal data to trace the fault source; Preferably, the following method is adopted to calculate the SPE contribution rate for abnormal data to trace the fault source: Suppose the abnormal data is \(x\in R\) 1×n , and the calculation formula for its SPE contribution rate is as follows: Among them, contspe(i) represents the SPE contribution corresponding to the i-th variable, and ξ i represents the i-th column of the n-dimensional identity matrix, T represents the transpose, I is the identity matrix, P is the load matrix obtained from the training samples, n is the number of variables, and the variable with the largest SPE contribution rate is the faulty variable.

9. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that Step (4) further includes: combining historical data with on-site process knowledge to flag the drawn confidence ellipse to achieve regional division of data with different performance levels. The regional division preferably includes: first, according to process knowledge and historical benefit statistics data, finding the corresponding benefit of the working condition according to the time series, dividing into different levels according to the benefit level, and then finding the distribution of the historical data corresponding to each level in the ellipse and marking them differently; or, step (4) further includes: setting different labels for historical data with different performance levels according to the performance indicators of the refining process, projecting the first two columns of the score matrix of historical data onto the confidence ellipse, and dividing the confidence ellipse through the labels carried by the data points in different regions.

10. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that, Step (7) further includes: for the samples mapped inside the ellipse, judging the current performance level according to their positions in combination with the regional division of the ellipse in step (4). If it is in a non-optimal performance level, adjust the key variables to make it switch to the desired better performance level; preferably, set the variable corresponding to the largest coefficient in the principal component as the key variable to be adjusted during optimization, and then determine the adjustment direction according to the correlation between the key variable and the change in the position of the projection point of the on-line data in the confidence ellipse to achieve the optimization of the production process; preferably, the key variable is the variable with the largest absolute value in the first column of the load matrix P. By adjusting the value corresponding to this variable, the production mode of the refining process is adjusted, and it is observed whether the data moves to the desired area through real-time monitoring.

11. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 1, characterized in that In step (8), the original variable X is calculated using formula (15). ori : X ori = (defen × PC T ) · × std(X m ) + mean(X m ) (15) where defen = (x, y) ∈ R 1×2 is the abscissa and ordinate of the points in the confidence ellipse, PC = [p1, p2] ∈ R J×2 are the first two columns of the load matrix during principal component analysis modeling, std(X m ) is the variance of the basic samples used in principal component analysis, mean(X m ) is the mean of the basic samples used in principal component analysis, m is the number of samples, X ori is the original variable after the inverse transformation corresponding to the position where defen = (x, y) ∈ R 1×2 in the ellipse; Calculate the input-output benefit index according to formula (16): Among them, is the output of the i-th product, is the price of the i-th product, is the feeding amount of the i-th raw material, is the price of the i-th raw material, and profit is the benefit index.

12. The method for pattern recognition and optimization of the oil refining process based on big data according to claim 11, wherein, The refining process is a catalytic reforming process, a catalytic cracking process, a sulfur recovery process, a residue hydrotreating process, an atmospheric and vacuum distillation process or a hydrocracking process.

Images