WO2023131257A1 - 一种基于大数据的炼油过程模式识别及优化方法 - Google Patents

一种基于大数据的炼油过程模式识别及优化方法 Download PDF

Info

Publication number
WO2023131257A1
WO2023131257A1 PCT/CN2023/070795 CN2023070795W WO2023131257A1 WO 2023131257 A1 WO2023131257 A1 WO 2023131257A1 CN 2023070795 W CN2023070795 W CN 2023070795W WO 2023131257 A1 WO2023131257 A1 WO 2023131257A1
Authority
WO
WIPO (PCT)
Prior art keywords
data
matrix
ellipse
optimization
confidence
Prior art date
Application number
PCT/CN2023/070795
Other languages
English (en)
French (fr)
Inventor
钟伟民
钱锋
杜文莉
李智
杨明磊
隆建
范琛
Original Assignee
华东理工大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 华东理工大学 filed Critical 华东理工大学
Publication of WO2023131257A1 publication Critical patent/WO2023131257A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/10—Analysis or design of chemical reactions, syntheses or processes

Definitions

  • the invention belongs to the technical field of oil refining process monitoring, and in particular relates to a big data-based oil refining process pattern recognition and optimization method.
  • the data acquisition system can collect a large amount of data in the oil refining process through various measuring instruments.
  • the changes of these data are often related to different production modes of the process, so effective monitoring of these data is of great significance to improve the production efficiency of the refining process and ensure its production safety.
  • the production processes such as catalytic reforming, catalytic cracking, sulfur recovery, residue hydrogenation, atmospheric and vacuum, and hydrocracking in the refining process have the characteristics of multi-variable, strong interference, large hysteresis, and strong coupling, and are very complex large-scale industries.
  • System how to dig out from the historical process data containing numerous operating variables and raw material property variables, which can fully reflect the impact of each parameter variable on the production process, clarify the type of production mode of the device, and identify the pros and cons of the device under various production modes Running the zone is a difficult thing.
  • it is a challenging task to judge the level of the device's operating status, adjust the process parameters in time, and realize the optimized operation of the production process.
  • the process operation mode based on the state at a single time point and the single factor curve within a time period can no longer meet the real-time monitoring of the changes of multiple modes in the refining process.
  • the multi-factor modal control method that can change with time, and realize the mode in time dimension and space. It is a collection of dimensions, and can dynamically monitor each key variable, and handle fault diagnosis, fault forecast, operation safety, dynamic optimization and static optimization under the concept of model.
  • the purpose of the present invention is to address the deficiencies in the existing methods, and propose a method for pattern recognition and optimization of the refining process based on big data.
  • the principal component analysis method By introducing the principal component analysis method, the key features of many variables in the process are extracted, and confidence ellipses are drawn.
  • On-line monitoring of the refining process is realized by using the mapping position of the data collected in real time in the ellipse.
  • the method can be used for monitoring, fault detection and traceability of the production mode in the refining process.
  • the improved path optimization algorithm can be further used, and the optimization direction of key variables can be given through the path optimization between the current working condition position and the corresponding position of the optimal benefit in the confidence ellipse.
  • the first aspect of the present invention provides a method for pattern recognition of oil refining processes based on big data, said method comprising the following steps:
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • step (4) Use the first column of scorey as the data of the x-axis, and the second column of scorey as the data of the y-axis, and map the scorey to the confidence ellipse drawn in step (4); if the sample points are mapped to the ellipse, It means that the working condition of the refining process at this time is normal; if the sample points are mapped outside the ellipse, it means that there is an abnormality in the refining process at this time.
  • step (2) the described preprocessing method adopts the Z-score normalization method, and the calculation formula is:
  • Z [z 1 ,z 2 ,...,z m ] is the training data matrix
  • X is the standardized data matrix
  • ⁇ is the mean value of the training data
  • ⁇ is the standard deviation of the training data
  • ⁇ and ⁇ The calculation formula is:
  • step (3) the preprocessed X is subjected to dimensionality reduction using the principal component analysis method, specifically as follows:
  • X is an m ⁇ n matrix, m is the number of samples, n is the number of variables, and T represents transposition, so the obtained covariance matrix C is an n ⁇ n-dimensional matrix;
  • PC i is the information ratio included in the i-th pivot
  • PC k is the information ratio of the largest k pivots.
  • step (4) also includes: according to the performance index of the refining process, setting different labels for the historical data of different performance levels, and projecting the first two columns of the score matrix of the historical data onto the confidence ellipse , divide the confidence ellipse by the labels of the data points in different regions.
  • step (8) the original variable X ori is calculated using formula (15):
  • g(x) is a cost function, which represents the actual cost required to reach the current node x from the starting node;
  • h(x) is a 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.
  • the oil refining process is a catalytic reforming process, a catalytic cracking process, a sulfur recovery process, a residual oil hydrogenation process, an atmospheric and vacuum process or a hydrocracking process, especially a catalytic reforming process.
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • the preprocessed X is subjected to dimensionality reduction using the principal component analysis method, specifically according to the following steps, and the specific operation steps are as described in the first step of the present invention: Described in step (3) of this aspect.
  • the calculation formula of scorey is as described in the step (6) of the first aspect of the present invention.
  • the step (7) further includes adding the collected data under normal operating conditions to historical data for remodeling.
  • the step (7) further includes calculating the SPE contribution rate for the abnormal data to trace the source of the fault, and the specific operation method is as described in the step (7) of the first aspect of the present invention.
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.
  • step (4) according to the set key characterization indicators of the production mode of the sulfur recovery unit, different labels are set for the historical data, when the scoring matrix of the historical data is projected
  • the confidence ellipse is divided into optimized area and non-optimized area by the labels of data points in different areas.
  • the calculation formula of scorey is as described in the step (6) of the first aspect of the present invention.
  • each set of data in scorey is substituted into the 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 less than or equal to 1, it means that the set of data is mapped in the ellipse.
  • step (7) for the sample mapped inside the two-dimensional confidence ellipse, judge whether it is in the optimized area of the defined production mode according to its position; if the mapping point is within If the mapping point is not within the optimization area, it indicates that the device is in the optimized operating state in the defined production mode at this time; if the mapping point is not in the optimized area, it indicates that the device is in the non-optimized operating state in the defined production mode at this time.
  • the method further includes step (8): For online real-time data in a non-optimized operating state, by observing the projected position of the data in a two-dimensional confidence ellipse in a defined production mode, From the components of each principal component, find 1 to 3 variables with the largest coefficients corresponding to each principal component, as the core variables for adjusting the operation status of the sulfur recovery unit, and determine the adjustment according to the correlation between the core variables and the position changes of the projection points. Direction, change the value of the core variables, to achieve the optimization of the production process.
  • 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 the variable, the production mode of the sulfur recovery unit Make adjustments and watch the data move in the desired area with real-time monitoring.
  • the data-driven method for identifying and optimizing the operating state of a sulfur recovery unit of the present invention includes: collecting historical process data including process operation, raw material evaluation, hydrogen composition analysis, and product distribution of the unit during the sulfur recovery unit, and screening Several groups of production data including device operation variables and dependent variables were obtained; principal component analysis was used to reduce the dimension of the above data, and a two-dimensional confidence ellipse was drawn; according to indicators such as sulfur recovery rate and H 2 S content in purified tail gas, set There are multiple production modes of the device, and according to the indicators and historical data, the optimal operation state area of the device in each mode is screened out; the real-time production data of the device operation is collected, and mapped to the two-dimensional confidence interval through the above method to identify whether the current operation state of the device is in operation optimization area; if the device is in a non-optimized operating state, by calculating the SPE contribution rate of each variable, find out the variable with a larger value for regulation, so that the production state of the device can gradually and timely
  • steps (1)-(8) of the third aspect of the invention are as described in any embodiment of the first aspect of the invention.
  • a fourth aspect of the present invention provides a data-driven method for identifying and optimizing the operating state of a residual oil hydrotreating unit, the method comprising the following steps:
  • the training data sample set Z [z 1 ,z 2 ,..., zi ,...,z n ] ⁇ R m ⁇ n , where m is the number of samples of training data, and n is the number of variables of training data;
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • step (4) Use the first column of scorey as the data of the x-axis, and the second column of scorey as the data of the y-axis, and map the scorey to the confidence ellipse drawn in step (4); if the sample points are mapped to the ellipse, It means that the residual oil hydrotreating unit is in normal working condition at this time; if the sample point is mapped outside the ellipse, it means that the operating state of the residual oil hydrotreating unit may be abnormal at this time.
  • the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.
  • the principal component analysis method is used to carry out dimensionality reduction processing on the X obtained through preprocessing, and the specific operation steps are as in the step (3) of the first aspect of the present invention ) mentioned.
  • the operation steps are as described in step (4) of the first aspect of the present invention.
  • step (4) according to the set key characterization indicators of the production mode of the residual oil hydrotreating unit, different labels are set for the historical data, and in the score matrix of the historical data During projection, the confidence ellipse is divided into optimized and non-optimized regions by the labels of data points in different regions.
  • the calculation formula of scorey is as described in the step (6) of the first aspect of the present invention.
  • each set of data in scorey is substituted into the 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 less than or equal to 1, it means that the set of data is mapped in the ellipse.
  • the step (7) further includes calculating the SPE contribution rate for the abnormal data to trace the source of the fault, and the specific operation method is as described in the step (7) of the first aspect of the present invention.
  • the step (7) for the sample mapped inside the two-dimensional confidence ellipse, it is judged according to its position whether it is in the optimized area of the set production mode; if the mapping point If the mapping point is not in the optimal area, it indicates that the device is in a non-optimized operating state under the set production mode at this time.
  • the method further includes step (8): For online real-time data in a non-optimized operating state, by observing the projected position of the data in a two-dimensional confidence ellipse in a defined production mode, From the components of each principal component, find 1 to 3 variables with the largest coefficients corresponding to each principal component, as the core variables for adjusting the operation status of the residual oil hydrotreating unit, according to the correlation between the core variables and the position changes of the projection points , determine the adjustment direction, change the value of the core variable, and realize the optimization of the production process.
  • the key variable is the variable with the largest absolute value in the first column of the load matrix P, and by adjusting the value corresponding to the variable, the production of the residual oil hydrotreating unit Adjust the mode and observe whether the data is moving to the desired area through real-time monitoring.
  • the data-driven method for identifying and optimizing the operating state of a residual oil hydrotreating unit of the present invention includes: collecting process operations, raw material and product property analysis, and unit product Distribution and other historical process data, screen out several groups of production data including device operating variables and dependent variables; according to device operating characteristics and data, define the production optimization operation mode of residual oil hydrotreating device, such as the maximum benefit of the device and the minimum energy consumption; Use the principal component analysis method to process data under different production modes, use the score matrix to build a model and draw a confidence ellipse; project the new sample calculation score matrix obtained in real time into the model ellipse, and judge the current situation according to the position of the sample point in the ellipse.
  • the production process of the device is in an optimized state under the production mode; add the samples projected into the ellipse to the historical data to build a new model, and realize the adaptive update of the model; if the sample point is projected outside the two-dimensional confidence ellipse, it indicates that the device In the non-optimized operation state, it is necessary to analyze the input variables of the model according to the fault contribution rate to realize fault traceability; screen out the variables with a large contribution rate to the fault, and adjust them according to the principal component and score results, so that the production status of the device can be gradually and timely returned to Go to the optimal operation area under the corresponding production mode to realize the optimization of the device.
  • steps (1)-(8) of the fourth aspect of the invention are as described in any embodiment of the first aspect of the invention.
  • a fifth aspect of the present invention provides a method for monitoring and optimizing a large data-based atmospheric and vacuum process mode, the method comprising the following steps:
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • step (4) Use the first column of scorey as the data of the x-axis, and the second column of scorey as the data of the y-axis, and map the scorey to the confidence ellipse drawn in step (4); if the sample points are mapped to the ellipse, It means that the working condition of the atmospheric and decompression process at this time is a normal working condition; if the sample point is mapped outside the ellipse, it means that there is an abnormality in the atmospheric and decompression process at this time.
  • the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.
  • the principal component analysis method is used to carry out dimensionality reduction processing on the X obtained through preprocessing, and the specific operation steps are as in the step (3) of the first aspect of the present invention ) mentioned.
  • the operation steps are as described in step (4) of the first aspect of the present invention.
  • the calculation formula of scorey is as described in the step (6) of the first aspect of the present invention.
  • each set of data in scorey is substituted into the 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 less than or equal to 1, it means that the set of data is mapped in the ellipse.
  • the step (7) further includes adding the collected data under normal operating conditions to historical data for remodeling.
  • the step (7) further includes calculating the SPE contribution rate for the abnormal data to trace the source of the fault, and the specific operation method is as described in the step (7) of the first aspect of the present invention.
  • the step (4) also includes: combining historical data with on-site process knowledge to perform flagging processing on the drawn confidence ellipse to achieve regional division of data of different performance levels, and the regional division is preferably Including: first calculate the KPI value (benefit value) of various performance indicators (including economic benefits, production energy consumption, product yield, etc.) Find the distribution of the historical data corresponding to each level in the ellipse and mark it with different colors (for example, color marking). Benefit value can be found in historical data or calculated based on historical data.
  • step (4) according to the performance index of the atmospheric and decompression process, different labels are set for the historical data, so when the score matrix of the historical data is projected, the data points in different regions can be used
  • the labels of the bands divide the confidence ellipses.
  • the method further includes step (8): for the samples mapped inside the ellipse, judging the current performance according to its location combined with the area division of the ellipse in step (4) grade.
  • the big data-based atmospheric and decompression process mode monitoring and optimization method of the present invention includes: classifying a large amount of historical data collected during the atmospheric and decompression process according to their different physical meanings, and Defined as different production modes; use the principal component analysis method to process data under different production modes, use the score matrix to build a model and draw a confidence ellipse; project the new sample calculation score matrix obtained in real time into the model, and according to the sample points in the ellipse The position in the ellipse is used to determine the current production mode of the process; the samples projected into the ellipse are added to the historical data to build a new model, and the model is adaptively updated; Carry out fault tracing; during modeling, variables that play a major role in process state changes under different production modes can be obtained, and these variables can be adjusted to achieve mode optimization.
  • steps (1)-(8) of the fifth aspect of the invention are as described in any embodiment of the first aspect of the invention.
  • a sixth aspect of the present invention provides a method for pattern recognition and optimization of a hydrocracking process based on big data, the method comprising the following steps:
  • Ym is multiplied by the first two columns of the load matrix P of step (3) gained, obtains Ym according to the first two column scoring matrix scorey ⁇ R N * 2 that training sample gains;
  • step (4) Use the first column of scorey as the data of the x-axis, and the second column of scorey as the data of the y-axis, and map the scorey to the confidence ellipse drawn in step (4); if the sample points are mapped to the ellipse, It indicates that the working condition of the hydrocracking process at this time is a normal working condition; if the sample point is mapped outside the ellipse, it indicates that there is an abnormality in the hydrocracking process at this time.
  • the pretreatment method in step (2) is as described in step (2) of the first aspect of the present invention.
  • the principal component analysis method is used to carry out dimensionality reduction processing on the X obtained through preprocessing, and the specific operation steps are as in the step (3) of the first aspect of the present invention ) mentioned.
  • the operation steps are as described in step (4) of the first aspect of the present invention.
  • the historical data combined with field process knowledge can be used to perform flagging processing on the drawn confidence ellipse to achieve regional division of data of different performance levels.
  • the regional division includes: first calculate the KPI value (benefit value) of each performance index (including economic benefit, production energy consumption, product yield, etc.) in the historical data according to the process knowledge, and then calculate the KPI value based on the actual The process is divided into different levels, and then the distribution of historical data corresponding to each level in the ellipse is found and marked with different marks (such as color marks). Benefit value can be found in historical data or calculated based on historical data.
  • step (4) according to the performance index of the hydrocracking process, different labels are set for the historical data, so when the score matrix of the historical data is projected, different regional data can be used
  • the labels attached to the points divide the confidence ellipse.
  • the calculation formula of scorey is as described in the step (6) of the first aspect of the present invention.
  • each set of data in scorey is substituted into the 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 less than or equal to 1, it means that the set of data is mapped in the ellipse.
  • the step (7) further includes adding the collected data under normal operating conditions to historical data for remodeling.
  • the step (7) further includes calculating the SPE contribution rate for the abnormal data to trace the source of the fault, and the specific operation method is as described in the step (7) of the first aspect of the present invention.
  • the step (8) also includes: for the samples at the non-optimal performance level, converting them to the expected better performance level by adjusting key variables; preferably, the main The variable corresponding to the largest coefficient in the element is set as the key variable to be adjusted during optimization, and then the adjustment direction is determined according to the correlation between the key variable and the position change of the projection point of the online data in the confidence ellipse, so as to realize the optimization of the production process.
  • 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 the variable, the production mode of the hydrocracking process is adjusted, and whether the data is observed through real-time monitoring Move to the desired area.
  • steps (1)-(8) of the sixth aspect of the invention are as described in any embodiment of the first aspect of the invention.
  • Fig. 13 is a projection diagram of the optimized area data corresponding to a certain production mode in the second embodiment on the confidence ellipse.
  • Fig. 16 is an overall flow chart of the method for identifying and optimizing the operating state of a sulfur recovery device according to the present invention.
  • Figure 21 is the optimization trajectory of the sulfur recovery unit in Example 3.
  • Fig. 24 is the distribution of different production modes of the residual oil hydrotreating unit in Example 4.
  • Fig. 25 is a projection diagram of the optimization area data corresponding to the production sample focusing on the yield of hydrogenated heavy oil in Example 4 on the confidence ellipse.
  • Fig. 27 is the optimization trajectory of the residual oil hydrotreating unit in Example 4.
  • Fig. 28 is a schematic flow chart of the atmospheric and vacuum process.
  • Fig. 30 is a confidence ellipse drawn according to normal data in Example 5.
  • Figure 36 is a process flow diagram of hydrocracking.
  • the refining process pattern recognition method based on big data of the present invention comprises the following steps:
  • the refinery process mode includes the working condition and production mode of the refinery process.
  • Working conditions can be divided into normal working conditions (working conditions without faults) and abnormal working conditions (working conditions with faults). Normal working conditions can be divided into different production modes according to different processes and performance evaluation indicators, such as economic benefit mode, device energy consumption mode, product yield mode, etc. According to the different evaluation indicators corresponding to various production modes, each production mode can be divided into optimal operation state and non-optimal operation state. For example, economic benefit mode can be divided into high economic benefit mode and low economic benefit mode, and device energy consumption mode can be divided into Divided into high device energy consumption mode and low device energy consumption mode, product yield mode can be divided into high product yield mode and low product yield mode, etc.
  • step (2) the training samples are standardized by calculating the mean value and standard deviation of the training sample set to obtain standardized data.
  • the method of normalizing the data based on the mean value and standard deviation is conventional in the art.
  • the standardized data presents a normal distribution with mean 0 and variance 1.
  • the preprocessing method adopts the Z-score normalization method, and the calculation formula is:
  • Z [z 1 ,z 2 ,...,z m ] is the training data matrix
  • X is the standardized data matrix
  • ⁇ is the mean value of the training data
  • ⁇ is the standard deviation of the training data
  • ⁇ and ⁇ The calculation formula is:
  • PC i is the information ratio included in the i-th pivot
  • PC k is the information ratio of the largest k pivots.
  • the steps of drawing a confidence ellipse using the matrix xdat are as follows:
  • step (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:
  • a is the major axis and b is the minor axis;
  • step (4-g) Taking xm in step (4-b) as the center point of the ellipse, a confidence ellipse can be drawn according to the center point and the major and minor axes of the ellipse.
  • the formula for the ellipse is:
  • xm1 is the average value of the first column of xdat
  • xm2 is the average value of the second column of xdat
  • the projections of the data in the confidence ellipse under different production modes will be distributed in different areas.
  • different labels can be set for the historical data according to the benefit classification corresponding to the historical working conditions.
  • Benefit values can be found in historical statistical data.
  • step (4) use the first two columns of the scoring matrix T to draw a two-dimensional confidence ellipse.
  • the confidence level is preferably set to 95%.
  • the variable setting and data sorting of key characterization indicators related to the definition of production mode can be completed, such as economic benefits, device energy consumption, product yield, etc.
  • Classification is carried out, the data in the optimized state is screened out, and the samples corresponding to such data are projected on the two-dimensional confidence ellipse to obtain the optimized area in the defined production mode, which is marked with different colors.
  • step (4) further includes dividing the confidence ellipse into regions according to performance levels.
  • the regional division includes: combining historical data with on-site process knowledge to perform flagging processing on the drawn confidence ellipse to achieve regional division of data of different performance levels.
  • the regional division preferably includes: first, based on process knowledge and historical Benefit statistical data, according to the time series to find the benefits corresponding to the working conditions, divide them into different levels according to the level of benefits, and then find the distribution of historical data corresponding to each level in the ellipse and mark them with different colors.
  • regional division includes: according to the performance indicators of the refining process, setting different labels for historical data of different performance levels, projecting the first two columns of the score matrix of historical data onto a confidence ellipse, and passing data points in different regions The labels carried divide the confidence ellipse.
  • step (5) Y is preprocessed with reference to formula (1) using the mean and variance obtained when preprocessing the training samples in step (2).
  • step (6) the calculation formula of scorey is:
  • step (7) by judging whether the sample point is mapped in the ellipse, consider whether to consider this sample point as the normal working condition of the device and include it in the historical data for remodeling to improve the adaptability of the model; if the sample mapping point is in If it is within the ellipse, it will be included in the database to update the model, otherwise it indicates that there may be abnormalities in the production process of the device at this time, and the SPE contribution rate can be calculated for the abnormal data to trace the source of the fault.
  • each set of data in scorey can be substituted into the aforementioned formula (12) for comparison 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 Group data is mapped within the ellipse.
  • step (7) also includes adding the collected data under normal operating conditions to historical data for remodeling; preferably, while adding the collected data under normal operating conditions, the same amount of early historical data.
  • step (7) further includes calculating the SPE contribution rate for the abnormal data to trace the source of the fault.
  • the method known in the art or the method provided by the present invention can be used to calculate the SPE contribution rate.
  • the formula for calculating the SPE contribution rate is:
  • 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
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • step (7) further includes: for the samples mapped inside the ellipse, judging the current performance level according to its location combined with the area division of the ellipse in step (4), if it is in a non- The optimal performance level is converted to the desired better performance level by adjusting the key variable; preferably, the variable with the largest corresponding coefficient in the principal element is set as the key variable to be adjusted during optimization, and then According to the correlation between the key variable and the position change of the projected point of the online data in the confidence ellipse, determine the adjustment direction to realize the optimization of the production process; preferably, the key variable is the one with the largest absolute value in the first column of the load matrix P Variable, by adjusting the value corresponding to the variable, adjust the production mode of the refining process, and observe whether the data moves to the desired area through real-time monitoring.
  • the method of the present invention may enable optimization of the refinery process model.
  • the present invention also includes a big data-based refinery process mode optimization method, which includes the refinery process mode recognition method as described in any embodiment of the present invention and any implementation of the present invention. Step (8) as described in the protocol.
  • each point can obtain the corresponding original variable through the standard deviation and mean value obtained in the principal component analysis process, so as to obtain the benefit value corresponding to the point.
  • the distribution of high and low benefit values can be distinguished in the confidence ellipse.
  • std(X m ) is the variance of the basic sample used in the principal component analysis
  • mean(X m ) is the mean value of the basic sample used in the principal component analysis
  • m is the number of samples
  • Statistical information in the production process includes feed quantity, product output, etc., and simplified input-output benefit indicators can be obtained based on price information.
  • the input-output benefit index is calculated according to formula (16):
  • the points in the confidence ellipse can be transformed by formula (15) to obtain the original working condition and a corresponding benefit value can be obtained according to formula (16), and the point corresponding to the optimal benefit value can be found at the same time.
  • the benefit value the high-benefit area and the low-benefit area can be divided in the confidence ellipse.
  • the path optimization may use 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 finds the optimal path among all reachable paths from the starting point to the target point after the starting point is determined.
  • the conventional A* algorithm determines the search direction and the next arrived node by evaluating the function f(x).
  • the general form of the evaluation function is as follows:
  • g(x) is a cost function, which represents the actual cost required to reach the current node x from the starting node;
  • h(x) is a heuristic function, representing the estimated cost required to reach the target node from the current node x.
  • the improved A* algorithm of the present invention is used to optimize the path, and the improved A* algorithm improves the evaluation function f(x):
  • profit(x) represents the economic benefit corresponding to the selected node
  • g(x) is a cost function, representing the actual cost required to reach the current node x from the starting node
  • h(x) is a heuristic function, representing The estimated cost for node x to reach the target node.
  • the optimization operation can be performed according to the original value of the operating variable obtained through transformation using formula (15).
  • Refinery production contains a large amount of process data, such as temperature, pressure, liquid level, flow rate and material properties, etc. These variables will have a certain impact on the production process. However, in the mass data, some data have a more significant impact on the production process, while others have less 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.
  • the refining process under normal conditions can be divided into different production modes according to its process conditions and performance indicators, such as high economic benefit mode or low economic benefit mode.
  • the data of different production modes will be mapped in different areas within the confidence ellipse,
  • the ellipse is divided by mapping the historical data under these modes in the ellipse.
  • the production mode of the current refining process can be judged by the position of the real-time collected data mapped in the ellipse, and the path optimization and The operating variable adjustment method obtained by the inverse transformation of the mode point guides the optimization of the production operation of the refinery unit.
  • the present invention adopts the A* algorithm as the optimization algorithm, which can find an optimal path in the good and bad intervals, and at the same time can deduce the optimal operating conditions through the points on the path.
  • the big data-based refinery process mode optimization method of the present invention is applied to a catalytic reforming (continuous catalytic reforming, CCR) process.
  • Figure 1 shows the flow chart of the catalytic reforming process.
  • the catalytic reforming process consists of a pre-hydrogenation unit, a reforming unit, and a catalyst regeneration system. For the purpose of producing aromatics, it also includes aromatics extraction and rectification units.
  • the pretreated raw material enters the reforming section, is mixed with circulating hydrogen and heated, and then enters the reactor. There are 3 to 4 reactors connected in series, and a heating furnace is installed between them to compensate for the heat absorbed by the reaction.
  • the material leaving the reactor enters the separator to separate the hydrogen cycle gas (the excess part is discharged), and the obtained liquid is used as reformed gasoline after removing light components from the stabilization tower, which is a high-octane gasoline component, or sent to aromatics extraction
  • the device produces aromatics.
  • the CCR model includes 51 input variables and 16 output variables.
  • the CCR models of different operating states can be obtained by controlling the input variables.
  • the input variables that have a greater impact on the output results include reaction temperature, pressure, feed flow, etc.
  • the performance indicators of the product are mainly obtained by observing the output variables of hydrogen, pure Hydrogen, dry gas, liquefied petroleum gas, carbon 5, carbon 6, carbon 7+, benzene mass flow and the amount of aromatics are used to judge.
  • Table 1 lists the operating variables involved in this embodiment.
  • variable describe unit 1 a reaction temperature °C 2 Second stage reaction temperature °C 3 Three-stage reaction temperature °C 4 Four stage reaction temperature °C 5 Circulating hydrogen STDm 3 /h 6 Total pre-hydrogenation feed tonne/h 7 compressor pressure MPa 8 T201 tray temperature °C 9 T201 return flow tonne/h 10 T201 tray temperature °C 11 T601 tower bottom temperature °C 12 extract xylene tonne/h 13 Reforming Feed Load tonne/h
  • ⁇ is the mean value taken from the training data:
  • ⁇ is the standard deviation taken from the training data:
  • X is an m ⁇ n matrix, m is the number of samples, and n is the number of variables, so the obtained covariance matrix C is an n ⁇ n-dimensional matrix;
  • Each pivot includes information ratios:
  • a is the major axis and b is the minor axis;
  • std(X m ) is the variance of the basic sample used in the principal component analysis
  • mean(X m ) is the mean value of the basic sample used in the principal component analysis
  • m is the number of samples
  • Each pivot includes information ratios:
  • a is the major axis and b is the minor axis;
  • xm1 is the average value of the first column of xdat
  • xm2 is the average value of the second column of xdat
  • the drawn ellipse is shown in Figure 11.
  • the projection of data in different production modes in the confidence ellipse will be distributed in different areas.
  • different labels are set for the historical data, and the optimization target variable—the economic benefit polyline is drawn, and the optimal operation area of the device is selected from it.
  • Project the first two columns of the score matrix of historical data divide the confidence ellipse by the labels of the data points in different regions, and distinguish them with different colors.
  • Fig. 13 shows the division result of the confidence ellipse area in this embodiment, and the area where the light-colored points are located is the area corresponding to the production mode with the best economic benefit.
  • the process After operating the catalytic cracking unit under normal working conditions for a period of time, the process deviates from the normal state by adjusting the density of the lower part of the burnt tank in the operating variable, and collects data to obtain Y ⁇ R N ⁇ n , and uses training on the collected data The average value and standard deviation obtained from data calculation are normalized to obtain Ym ⁇ R N ⁇ n .
  • 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.
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • the SPE contribution rate is shown in Figure 14. According to Figure 14, it can be determined that the cause of the failure is mainly due to the tenth variable (ie, the density of the lower part of the charred tank), which is consistent with the actual operation of step (6).
  • a method for identifying and optimizing the operating state of a sulfur recovery device based on data-driven according to the present invention is used to identify and optimize the operating state of the sulfur recovery device, as shown in FIG. 16 , including the following steps:
  • ⁇ is the mean value of the training data, and its calculation formula is:
  • ⁇ is the standard deviation of the training data, and its calculation formula is:
  • X is an m ⁇ n matrix
  • m is the number of samples
  • n is the number of variables
  • the covariance matrix C is an n ⁇ n-dimensional matrix
  • I is unit matrix
  • Each pivot includes information ratios:
  • a is the major axis and b is the minor axis;
  • xm1 is the average value of the first column of xdat
  • xm2 is the average value of the second column of xdat
  • the drawn ellipse is shown in Figure 17.
  • the projection of data in different production modes in the confidence ellipse will be distributed in different areas.
  • the corresponding optimization target variable line is drawn, and the optimal operating area of the device is selected (sulfur recovery rate> 99.5%).
  • Project the first two columns of the score matrix of historical data divide the confidence ellipse by the labels of the data points in different regions, and distinguish them with different colors.
  • Fig. 19 shows the division result of the confidence ellipse area in this embodiment, and the area where the light-colored points are located is the area corresponding to the production mode that pursues the highest sulfur recovery rate.
  • 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.
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • the SPE contribution rate is shown in Figure 20. According to Figure 20, it can be judged that the cause of the failure is mainly caused by the eleventh variable (namely, the temperature of the sulfur furnace), which is consistent with the actual operation.
  • the 1.7 million tons/year residual oil hydrotreating unit was used as the identification and optimization object, and historical production data from May 1, 2018 to May 1, 2019 were collected, and abnormal fluctuation data were eliminated. After screening, the data were sorted into 3925 groups with 85 variables. Table 6 lists the variables involved in this embodiment.
  • ⁇ is the mean value of the training data, and its calculation formula is:
  • ⁇ is the standard deviation of the training data, and its calculation formula is:
  • Each pivot includes information ratios:
  • 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.
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • the SPE contribution rate is shown in Figure 26. According to Figure 26, it can be judged that the cause of the failure is mainly caused by the 20th variable (namely, the flow rate of cold slag in the tank area), which is consistent with the actual operation situation.
  • the atmospheric and decompression process mode monitoring and optimization method based on big data of the present invention is used for mode monitoring and optimization of the atmospheric and decompression process, as shown in Figure 29, comprising the following steps:
  • ⁇ is the mean value of the collected data, and its calculation formula is:
  • ⁇ is the standard deviation of the collected data, and its calculation formula is:
  • X is an m ⁇ n matrix
  • m is the number of samples
  • n is the number of variables
  • T represents transposition
  • the covariance matrix C is an n ⁇ n-dimensional matrix
  • I is unit matrix
  • V [p 1 ,p 2 ,p 3 ,...,p n ] ⁇ R n ⁇ n ;
  • Each pivot includes information ratios:
  • a is the major axis and b is the minor axis;
  • xm1 is the average value of the first column of xdat
  • xm2 is the average value of the second column of xdat
  • the drawn ellipse is shown in Figure 30.
  • 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.
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • the SPE contribution rate is shown in Figure 35. According to Figure 35, it can be judged that the cause of the failure is mainly caused by the third variable (ie, the pressure at the top of the initial distillation column), which is consistent with the actual operation (variable 1 is the raw material property variable, which is difficult to directly adjust in actual operation) .
  • FIG. 36 shows the process flow chart of hydrocracking.
  • the hydrocracking process consists of hydrorefining, hydrocracking and fractionation.
  • the raw oil and circulating hydrogen respectively meet the conditions of hydroprocessing and pressurization, and then they are mixed into the hydrofinishing reactor, and desulfurization, denitrogenation, deoxygenation and partial dearomatization are carried out under the action of the hydrofinishing catalyst.
  • the refined reaction product is adjusted to the required temperature for the cracking reaction by injection of cold hydrogen, and then enters the hydrocracking reactor for cracking reaction to convert heavy distillate oil into light distillate oil.
  • the reaction product undergoes three-phase separation of gas, oil and water in the cold high-pressure separator, and is sent to the low-pressure separator after being desulfurized by the circulating hydrogen desulfurization tower, where low-component gas is flashed out.
  • the low-fraction oil enters the debutanizer through heat exchange, and the sulfur-containing gas and light hydrocarbons are separated at the top of the tower, and the bottom oil enters the product separation tower for product separation, and products such as light and heavy naphtha, aviation fuel, diesel oil and tail oil are separated .
  • the HCR model includes 32 input variables and 24 output variables. HCR models in different operating states can be obtained by controlling the input variables. 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.
  • the performance indicators of the product are mainly obtained through Observe the output variables of dry gas, liquefied petroleum gas, light end, light naphtha, heavy naphtha, jet fuel, diesel and tail oil out of the device to judge. Table 8 shows the operating variables of the HCR model of this embodiment.
  • variable describe unit variable describe unit 1 Tank area light wax oil flow rate tonne/h 17 Raw oil sulfur content % 2 Atmospheric and vacuum wax oil total flow tonne/h 18 Raw Oil Nitrogen Content ppmwt 3 Cumulative flow of catalytic diesel tonne/h 19 Raw Oil Basic Nitrogen ppmwt 4 Mass flow rate of circulating tail oil tonne/h 20 raw oil initial boiling point C 5 new hydrogen flow STD_m3/h twenty one Raw oil 10% recovery temperature C 6 Hydrogen flow rate for purge cycle STD_m3/h twenty two Raw oil 50% recovery temperature C 7 Hydrogen oil ratio the twenty three Raw oil 90% recovery temperature C 8 R101-bed inlet temperature C twenty four raw oil final boiling point C 9 R101 second bed inlet temperature C 25 Light oil density (20°C) the 10 R101 three-bed inlet temperature C 26 Sulfur content of light diesel oil % 11 R102-bed inlet temperature C 27 Initial boiling point of light diesel oil C 12 R102 second bed inlet temperature C 28 10% recovery temperature of light diesel oil C 13 R102 three-bed
  • pattern recognition and optimization method of the present invention to carry out pattern recognition and optimization of the hydrocracking process, as shown in Figure 37, includes the following steps:
  • ⁇ is the mean value of the collected data, and its calculation formula is:
  • ⁇ is the standard deviation of the collected data, and its calculation formula is:
  • X is an m ⁇ n matrix
  • m is the number of samples
  • n is the number of variables
  • T represents transposition
  • the covariance matrix C is an n ⁇ n-dimensional matrix
  • I is unit matrix
  • V [p 1 ,p 2 ,p 3 ,...,p n ] ⁇ R n ⁇ n ;
  • Each pivot includes information ratios:
  • a is the major axis and b is the minor axis;
  • xm1 is the average value of the first column of xdat
  • xm2 is the average value of the second column of xdat
  • the drawn ellipse is shown in Figure 38.
  • the projection of the data in different production modes in the confidence ellipse will be distributed in different areas.
  • different labels are set for the historical data, and the first two columns of the score matrix of the historical data are projected.
  • the confidence ellipse is divided and framed by the labels attached to the data points in different regions.
  • Fig. 39, Fig. 40 and Fig. 41 respectively show the models corresponding to the optimization target of total liquid recovery and the target of intermediate oil yield.
  • the big box is the top 10% data points selected according to the corresponding target, and the small box is the top 5% data points.
  • 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.
  • the variable with the largest SPE contribution rate is the variable that has failed.
  • the SPE contribution rate is shown in Figure 43. According to Figure 43, it can be judged that the cause of the fault is mainly caused by the first variable (ie hydrogen-oil ratio), which is consistent with the actual operating conditions.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Artificial Intelligence (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Analytical Chemistry (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种基于大数据的炼油过程模式识别及优化方法，属于炼油过程监控技术领域，包括以下步骤：对炼油过程中所采集到的历史数据进行预处理，得到标准化数据；运用主成分分析法处理标准化数据，利用得分矩阵建立模型并绘制置信椭圆；对实时采集得到的新样本计算得分矩阵，将其投影到置信椭圆上；投影到椭圆内的样本为正常点，可将其添加到历史数据中建立新模型，实现模型自适应更新；投影到椭圆外的样本则为异常点，可根据故障贡献率进行故障溯源；进一步可根据置信椭圆中的点对应的原始变量，采用改进的路径优化算法，求解操作变量调整方式和路径。

Description

一种基于大数据的炼油过程模式识别及优化方法
本申请要求于2022年01月10日提交中国专利局、申请号为202210022589.X、发明名称为“一种基于大数据的炼油过程模式识别及优化方法”的中国专利申请的优先权，其全部内容通过引用结合在本申请中。
技术领域
本发明属于炼油过程监控技术领域，具体涉及一种基于大数据的炼油过程模式识别及优化方法。
背景技术
随着现代信息技术的不断完善，数据采集系统可以通过各种测量仪表在炼油过程中收集到大量的数据。这些数据的变化往往与过程的不同生产模态相关，因此对这些数据进行有效地监控对提高炼油过程的生产效益以及保障其生产安全有着重要的意义。
炼油过程中的催化重整、催化裂化、硫磺回收、渣油加氢、常减压、加氢裂化等生产过程具有多变量、强干扰、大滞后、强耦合等特点，是非常复杂的大工业系统，如何从包含众多的操作变量和原料性质变量的历史过程数据中，挖掘出可以充分反映各参数变量对生产过程的影响，明晰装置生产模式类型，甄别在各类生产模式下装置的优劣运行区域，是具有难度的事情。另外，根据当前装置运行数据，判断装置运行状态所处级别，及时调整过程参数，实现生产过程的优化运行，是具有挑战的工作。
目前，基于单个时间点的状态以及时间段内的单因素曲线的过程操作方式已经不能满足对炼油过程多个模态的变化进行实时地监控。为了综合考量过程中所采集的数据具有众多变量以及各个变量对过程的监控重要性随时间而变的特点，需要考虑可以随着时间变化的多因素模态控制方式，实现模式在时间维度和空间维度的集合表现，并且可以动态地监控各个关键变量，将故障诊断、故障预报、操作安全、动态优化和静态优化等在模式的概念下统一处理。
发明内容
本发明的目的在于针对现有方法中存在的不足，提出一种基于大数据的炼油过程模式识别及优化方法，通过引入主元分析方法，提取过程中众多变量的关键特征，并绘制置信椭圆，利用实时采集到的数据在椭圆中的映射位置，实现对炼油过程的在线监控。该方法可用于对炼油过程的生产模式进行监控、故障检测和溯源。可进一步采用改进的路径优化算法，通过置信椭圆中当前工况位置和最优效益对应位置间的路径优化，可以给出关键变量的优化方向。
具体而言，本发明的第一个方面提供一种基于大数据的炼油过程模式识别方法，所述方法包括以下步骤：
(1)将炼油过程中所采集到的历史数据组成建模用的训练样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为样本集的样本个数，n为样本集的变量个数；
(2)对训练样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)对X运用主元分析法将其从n维降到k维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；
(5)采集新的在线实时数据Y∈R N×n，采用步骤(2)对训练样本集进行预处理时所得到的样本平均值和样本方差对Y进行预处理得到标准化数据Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则说明此时炼油过程所处的工况为正常工况；若样本点映射到椭圆外，则说明此时炼油过程存在异常。
在一个或多个实施方案中，步骤(2)中，所述的预处理方法采用Z-score标准化方法，计算公式为：
Figure PCTCN2023070795-appb-000001
式中，Z＝[z 1,z 2,...,z m]为训练数据矩阵，X表示标准化后的数据矩阵，μ为训练数据的均值，σ为训练数据的标准差，μ和σ计算公式为：
Figure PCTCN2023070795-appb-000002
Figure PCTCN2023070795-appb-000003
在一个或多个实施方案中，步骤(3)中，对经过预处理所得到的X运用主元分析法进行降维处理，具体按如下步骤进行：
(3-a)计算矩阵X的协方差矩阵，协方差矩阵的计算公式为：
Figure PCTCN2023070795-appb-000004
X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，所以得到的协方差矩阵C是一个n×n维的矩阵；
(3-b)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序，计算公式为：
det(C-Iλ)＝0                        (5)
Cp i＝λ ip i                          (6)
式(5)中，λ为特征值集合，即λ＝λ 1,λ 2,...,λ n，I为单位矩阵，得：
特征值矩阵：
Figure PCTCN2023070795-appb-000005
其中λ 1＞λ 2＞...＞λ n
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
(3-c)建立主元模型，计算公式为：
计算每个主元包括的信息比：
Figure PCTCN2023070795-appb-000006
计算最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000007
其中，PC i为第i个主元包括的信息比，PC k为最大的k个主元的信息比。
(3-d)保留对应k个最大特征值的特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                       (9)
t i的本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
在一个或多个实施方案中，步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆的步骤如下：
(4-a)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
(4-b)计算xdat每一列的平均值得xm∈R 1×2，并将xdat每一列的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
(4-c)计算公式xd×s·×xd，并对所得矩阵每行进行求和，得到rd∈R m×1；
(4-d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度Confi，对rd进行升序排序后，得到rd第Confi个位置所对应的值r∈R，优选的Confi为95％；即得到rd前95％的数据位置所对应的值；
(4-e)计算步骤(4-a)所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000008
特征向量矩阵为
Figure PCTCN2023070795-appb-000009
(4-f)利用步骤(4-d)中的r以及步骤(4-e)中的D，即可得到置信椭圆的长轴以及短轴，计算公式为：
Figure PCTCN2023070795-appb-000010
Figure PCTCN2023070795-appb-000011
其中，a为长轴，b为短轴；
(4-g)以步骤(4-b)中的xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，椭圆的公式为：
Figure PCTCN2023070795-appb-000012
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值。
在一个或多个实施方案中，步骤(6)中，scorey的计算公式为：
scorey＝Ym×[p 1,p 2]                  (13)
式中，[p 1,p 2]∈R n×2是负载矩阵P的前两列。
在一个或多个实施方案中，步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模。
在一个或多个实施方案中，步骤(7)还包括将收集到的正常工况下的数据加入到历史数据并剔除相同数量的早期历史数据重新建模。
在一个或多个实施方案中，步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源；
优选地，采用如下方法对异常数据计算SPE贡献率来进行故障溯源：
假设异常数据为x∈R 1×n，其SPE贡献率的计算公式为：
Figure PCTCN2023070795-appb-000013
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000014
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量。
在一个或多个实施方案中，步骤(4)还包括：将历史数据结合现场工艺知识对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分，所述区域划分优选包括：首先根据工艺知识和历史效益统计数据，根据时间序列找到工况对应的效益，根据效益高低进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并加以不同的标记，例如用不同的颜色加以标记。
在一个或多个实施方案中，步骤(4)还包括：根据炼油过程的性能指标，对不同性能等级的历史数据设置不同的标签，将历史数据的得分矩阵的前两列投影到置信椭圆上，通过不同区域数据点所带的标签对置信椭圆进行划分。
本发明还提供一种基于大数据的炼油过程模式优化方法，所述基于大数据的炼油过程模式优化方法包括如本发明任一实施方案所述的基于大数据的炼油过程模式识别方法，且所述基于大数据的炼油过程模式识别方法中的步骤(7)还包括：对于映射在椭圆内部的样本，根据其所处位置结合步骤(4)中对椭圆进行的区域划分判断当前所处的性能等级，若其处于非最优区域，则通过调整关键变量使其转换到所期望的更优的性能等级；优选地，将主元中所对应的系数最大的变量设定为优化时所要调整的关键变量，然后根据关键变量与在线数据在置信椭圆中的投影点位置变化之间的相关性，确定调整方向，实现对生产过程的优化；优选地，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对炼油过程的生产模式进行调整，并通过实时监控观察数据是否往所期望的区域移动。
本发明还提供另一种基于大数据的炼油过程模式优化方法，所述基于大数据的炼油过程模式优化方法包括如本发明任一实施方案所述的基于大数据的炼油过程模式识别方法，和步骤(8)：利用主成分分析过程中得到的标准差和均值求得置信椭圆中的点对应的原始变量，从而得到该点对应的效益值；在置信椭圆中区分效益值高低的分布；当工况运行模式位于置信椭圆中的某一点时，采用路径优化算法，获得当前位置向最优位置最快移动轨迹，将轨迹对应的点进行反变换，获得操作条件的改变方式，从而指导生产装置操作优化。
在一个或多个实施方案中，步骤(8)中，利用公式(15)计算原始变量X ori：
X ori＝(defen×PC T)·×std(X m)+mean(X m)          (15)
其中，defen＝(x,y)∈R 1×2为置信椭圆中的点的横纵坐标，PC＝[p 1,p 2]∈R J×2为 主成分分析建模时负载矩阵的前两列，std(X m)为主成分分析采用的基础样本的方差，mean(X m)为主成分分析采用的基础样本的均值，m为样本个数，X ori为椭圆中defen＝(x,y)∈R 1×2位置对应的反变换后的原始变量；
根据公式(16)计算投入产出效益指标：
Figure PCTCN2023070795-appb-000015
其中，
Figure PCTCN2023070795-appb-000016
为第i种产品的产量，
Figure PCTCN2023070795-appb-000017
为第i种产品的价格，
Figure PCTCN2023070795-appb-000018
为第j种原料的投料量，
Figure PCTCN2023070795-appb-000019
为第j种原料的价格，profit为效益指标。
在一个或多个实施方案中，步骤(8)中，所述路径优化算法采用改进的A*算法，通过如式(17)所示的改进的评估函数f(x)来确定搜索方向和下一个到达的节点：
Figure PCTCN2023070795-appb-000020
其中，g(x)是代价函数，代表从起始节点到达当前节点x的实际所需代价；h(x)是启发式函数，代表从当前节点x到达目标节点估计所需的代价；profit(x)表示所选节点对应的经济效益。
在一个或多个实施方案中，所述炼油过程为催化重整过程、催化裂化过程、硫磺回收过程、渣油加氢过程、常减压过程或加氢裂化过程，特别是催化重整过程。
本发明的第二个方面提供一种基于数据驱动的催化裂化装置运行状态识别及优化方法，所述方法包括以下步骤：
(1)将包含操作数据和原料数据的正常工况状态下的催化裂化装置生产过程历史数据组成建模用的训练数据样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为训练数据的样本个数，n为训练数据的变量个数；
(2)对训练数据样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)对X运用主元分析法将其从n维降到k维，并计算得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；定义生产模式，设定定义的生产模式下的评价指标，根据需求对定义的生产模式下的评价指标进行等级划分，筛选出处于优化运行状态的数据，将处于优化运行状态的数据对应的样本投影到二维置信椭圆上，得到定义的生产模式下的优化区域；
(5)采集在线实时数据Y∈R N×n，采用步骤(2)对训练数据进行预处理时所得到的平均值和方差对Y进行预处理得到Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则表示此时催化裂化装置处于正常工况状态；若样本点映射到椭圆外，则表示此时催化裂化装置运行状态可能存在异常；进一步地，对于映射在二维置信椭圆内部的样本，根据其所处位置判断其是否在定义的生产模式的优化区域内；如果映射点在优化区域内，则表明此时装置处于定义的生产模式下的优化运行状态；如果映射点不在优化区域内，则表明此时装置处于定义的生产模式下的非优化运行状态。
在一个或多个实施方案中，所述步骤(2)中的预处理方法如本发明的第一个方面的步骤(2)所述。
在一个或多个实施方案中，所述步骤(3)中，对经过预处理所得到的X运用主元分析法进行降维处理，具体按如下步骤进行，具体操作步骤如本发明的第一个方面的步骤(3)所述。
在一个或多个实施方案中，所述步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆，具体操作步骤如本发明的第一个方面的步骤(4)所述。
在一个或多个实施方案中，所述步骤(4)中，定义生产模式，设定定义的生产模式下的评价指标，根据定义的生产模式的评价指标，对历史数据设置不同的标签，在对历史数据的得分矩阵进行投影时通过不同区域数据点所带的标签对置信椭圆进行区域划分；或根据需求对定义的生产模式下的评价指标进行等级划分，筛选出处于优化运行状态的数据，将处于优化运行状态的数据对应的样本投影到二维置信椭圆上，得到定义的生产模式下的优化区域。
在一个或多个实施方案中，所述步骤(6)中，scorey的计算公式如本发明的第一个方面的步骤(6)所述。
在一个或多个实施方案中，所述步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，所述步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模。
在一个或多个实施方案中，所述步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源，具体操作方法如本发明的第一个方面的步骤(7)所述。
在一个或多个实施方案中，所述步骤(7)中，对于映射在二维置信椭圆内部的样本，根据其所处位置判断其是否在定义的生产模式的优化区域内；如果映射点在优化区域内，则表明此时装置处于定义的生产模式下的优化运行状态；如果映射点不在优化区域内，则表明此时装置处于定义的生产模式下的非优化运行状态。
在一个或多个实施方案中，所述方法还包括步骤(8)：对于处于非优化运行状态的在线实时数据，通过观察该数据在定义的生产模式下在二维置信椭圆中的投影位置，从由装置操作数据和原料数据降维得到的各主元成分中，找到各主元所对应的系数最大的变量1～3个，作为调整催化裂化装置运行状态的核心变量，根据核心变量与投影点位置变化之间的相关性，确定调整方向，改变核心变量的数值，实现对生产过程的优化。
在一个或多个实施方案中，所述步骤(8)中，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对催化裂化装置的生产模式进行调整，并通过实时监控观察数据是否往期望的区域移动。
在一个或多个实施方案中，本发明的第二个方面的步骤(1)-(8)如本发明的第一个方面的任一实施方案所述。
本发明的第三个方面提供一种基于数据驱动的硫磺回收装置运行状态识别及优化方法，所述方法包括以下步骤：
(1)将正常工况状态下的硫磺回收装置生产过程历史数据组成建模用的训练数据样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为训练数据的样本个数，n为训练数据的变量个数；
(2)对训练数据样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)运用主元分析法将X从n维降到k维，计算得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；设定硫磺回收装置的生产模式，根据所设定的生产模式下的关键指标变量和生产需求，对所设定的生产模式进行等级划分，筛选出处于优化运行状态的数据，将处于优化运行状态的数据对应的样本投影到二维置信椭圆上，得到所设定的生产模式下的优化区域；
(5)采集在线实时数据Y∈R N×n，采用步骤(2)对训练数据进行预处理时所得到的平均值和方差对Y进行预处理得到Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则表示此时硫磺回收装置处于正常工况状态；若样本点映射到椭圆外，则表示此时硫磺回收装置运行状态可能存在异常。
在一个或多个实施方案中，所述步骤(2)中的预处理方法如本发明的第一个方面的步骤(2)所述。
在一个或多个实施方案中，所述步骤(3)中，运用主元分析法对经过预处理所得到的X进行降维处理，具体操作步骤如本发明的第一个方面的步骤(3)所述。
在一个或多个实施方案中，所述步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆，具体操作步骤如本发明的第一个方面的步骤(4)所述。
在一个或多个实施方案中，所述步骤(4)中，根据所设定的硫磺回收装置生产模式的关键表征指标，对历史数据设置不同的标签，在对历史数据的得分矩阵进行投影时通过不同区域数据点所带的标签将置信椭圆划分为优化区域和非优化区域。
在一个或多个实施方案中，所述步骤(6)中，scorey的计算公式如本发明的第一个方面的步骤(6)所述。
在一个或多个实施方案中，所述步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，所述步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源，具体操作方法如本发明的第一个方面的步骤(7)所述。
在一个或多个实施方案中，所述步骤(7)中，对于映射在二维置信椭圆内部的样本，根据其所处位置判断其是否在定义的生产模式的优化区域内；如果映射点在优化区域内，则表明此时装置处于定义的生产模式下的优化运行状态；如果映射点不在优化区域内，则表明此时装置处于定义的生产模式下的非优化运行状态。
在一个或多个实施方案中，所述方法还包括步骤(8)：对于处于非优化运行状态的在线实时数据，通过观察该数据在定义的生产模式下在二维置信椭圆中的投影位置，从各主元成分中，找到各主元所对应的系数最大的变量1～3 个，作为调整硫磺回收装置运行状态的核心变量，根据核心变量与投影点位置变化之间的相关性，确定调整方向，改变核心变量的数值，实现对生产过程的优化。
在一个或多个实施方案中，所述步骤(8)中，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对硫磺回收装置的生产模式进行调整，并通过实时监控观察数据是否往期望的区域移动。
在一些实施方案中，本发明的基于数据驱动的硫磺回收装置运行状态识别及优化方法包括：采集硫磺回收装置过程中包含过程操作、原料评价、氢气组成分析以及装置产品分布等历史过程数据，筛选出包含装置操作变量和因变量的生产数据若干组；采取主成分分析法对上述数据进行降维，绘制出二维置信椭圆；根据硫磺回收率、净化尾气中H 2S含量等指标，设定装置多种生产模式，并根据指标和历史数据，筛选出各模式下装置优化运行状态区域；采集装置运行实时生产数据，通过上述方法映射至二维置信区间，识别装置当前运行状态是否处于运行优化区域；若装置处于非优化运行状态，则通过计算各变量的SPE贡献率，找出其数值较大的变量进行调控，使装置生产状态逐步、及时回到对应生产模式下最优运行区域。
在一个或多个实施方案中，本发明的第三个方面的步骤(1)-(8)如本发明的第一个方面的任一实施方案所述。
本发明的第四个方面提供一种基于数据驱动的渣油加氢处理装置运行状态识别及优化方法，所述方法包括以下步骤：
(1)将渣油加氢处理装置生产过程历史数据组成建模用的训练数据样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为训练数据的样本个数，n为训练数据的变量个数；
(2)对训练数据样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)运用主元分析法将X从n维降到k维，计算得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；设定渣油加氢处理装置的生产模式，根据所设定的生产模式下的指标变量，对数据进行等级划分，筛选出处于优化运行状态的数据，将处于优化运行状态的数据对应的样本投影到二维置信椭圆上，得到所设定的生产模式下的优化区域；
(5)采集在线实时数据Y∈R N×n，采用步骤(2)对训练数据进行预处理时所得到的平均值和方差对Y进行预处理得到Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则表示此时渣油加氢处理装置处于正常工况状态；若样本点映射到椭圆外，则表示此时渣油加氢处理装置运行状态可能存在异常。
在一个或多个实施方案中，所述步骤(2)中的预处理方法如本发明的第一个方面的步骤(2)所述。
在一个或多个实施方案中，所述步骤(3)中，运用主元分析法对经过预处理所得到的X进行降维处理，具体操作步骤如本发明的第一个方面的步骤(3)所述。
在一个或多个实施方案中，所述步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆，具体操作步骤如本发明的第一个方面的步骤(4)所述。
在一个或多个实施方案中，所述步骤(4)中，根据所设定的渣油加氢处理装置生产模式的关键表征指标，对历史数据设置不同的标签，在对历史数据的得分矩阵进行投影时通过不同区域数据点所带的标签将置信椭圆划分为优化区域和非优化区域。
在一个或多个实施方案中，所述步骤(6)中，scorey的计算公式如本发明的第一个方面的步骤(6)所述。
在一个或多个实施方案中，所述步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，所述步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源，具体操作方法如本发明的第一个方面的步骤(7)所述。
在一个或多个实施方案中，所述步骤(7)中，对于映射在二维置信椭圆内部的样本，根据其所处位置判断其是否在设定的生产模式的优化区域内；如果映射点在优化区域内，则表明此时装置处于设定的生产模式下的优化运行状态；如果映射点不在优化区域内，则表明此时装置处于设定的生产模式下的非优化运行状态。
在一个或多个实施方案中，所述方法还包括步骤(8)：对于处于非优化运行状态的在线实时数据，通过观察该数据在定义的生产模式下在二维置信椭圆中的投影位置，从各主元成分中，找到各主元所对应的系数最大的变量1～3 个，作为调整渣油加氢处理装置运行状态的核心变量，根据核心变量与投影点位置变化之间的相关性，确定调整方向，改变核心变量的数值，实现对生产过程的优化。
在一个或多个实施方案中，步骤(8)中，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对渣油加氢处理装置的生产模式进行调整，并通过实时监控观察数据是否往期望的区域移动。
在一个或多个实施方案中，本发明的基于数据驱动的渣油加氢装置运行状态识别及优化方法包括：采集渣油加氢处理装置过程中包含过程操作、原料和产品性质分析以及装置产品分布等历史过程数据，筛选出包含装置操作变量和因变量的生产数据若干组；根据装置运行特性和数据，定义渣油加氢处理装置生产优化运行模式，如装置效益最大、能耗最小等；运用主成分分析法处理不同生产模式下的数据，利用得分矩阵建立模型并绘制置信椭圆；将实时采集得到的新样本计算得分矩阵投影到模型椭圆中，根据样本点在椭圆中的位置，判断目前装置生产过程所处在该生产模式下是否处于优化状态；将投影到椭圆内的样本添加到历史数据中建立新模型，实现模型自适应更新；若样本点投影到二维置信椭圆外则表明装置处于非优化运行状态，此时需根据故障贡献率对模型输入变量分析，实现故障溯源；筛选出对故障贡献率大的变量，按照主元和得分结果进行调控，使装置生产状态逐步、及时回到对应生产模式下最优运行区域，实现装置的优化。
在一个或多个实施方案中，本发明的第四个方面的步骤(1)-(8)如本发明的第一个方面的任一实施方案所述。
本发明的第五个方面提供一种基于大数据的常减压过程模式监控和优化方法，所述方法包括以下步骤：
(1)将常减压过程的历史数据组成建模用的训练样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为样本集的样本个数，n为样本集的变量个数；
(2)对训练数据样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)对X运用主元分析法将其从n维降到k维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；
(5)采集新的在线实时数据Y∈R N×n，采用步骤(2)对训练样本进行预处理时所得到的样本平均值和样本方差对Y进行预处理得到标准化数据Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则说明此时常减压过程所处的工况为正常工况；若样本点映射到椭圆外，则说明此时常减压过程存在异常。
在一个或多个实施方案中，所述步骤(2)中的预处理方法如本发明的第一个方面的步骤(2)所述。
在一个或多个实施方案中，所述步骤(3)中，运用主元分析法对经过预处理所得到的X进行降维处理，具体操作步骤如本发明的第一个方面的步骤(3)所述。
在一个或多个实施方案中，所述步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆，具体操作步骤如本发明的第一个方面的步骤(4)所述。
在一个或多个实施方案中，所述步骤(6)中，scorey的计算公式如本发明的第一个方面的步骤(6)所述。
在一个或多个实施方案中，所述步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，所述步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模。
在一个或多个实施方案中，所述步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源，具体操作方法如本发明的第一个方面的步骤(7)所述。
在一个或多个实施方案中，所述步骤(4)还包括：将历史数据结合现场工艺知识对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分，所述区域划分优选包括：首先根据工艺知识计算历史数据中各项性能指标(包括经济效益、生产能耗、产品产率等)的KPI值(效益值)，然后对KPI值依据实际工艺进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并用不同的颜色加以标记(例如颜色标记)。效益值在历史数据中可以查到或根据历史数据计算得到。
在一个或多个实施方案中，步骤(4)中，根据常减压过程的性能指标，对历史数据设置不同的标签，因此在对历史数据的得分矩阵进行投影时可以通过不同区域数据点所带的标签对置信椭圆进行划分。
在一个或多个实施方案中，所述方法还包括步骤(8)：对于映射在椭圆内部的样本，根据其所处位置结合步骤(4)中对椭圆进行的区域划分判断当前所处的性能等级。
在一个或多个实施方案中，所述步骤(8)还包括：对于处于非最优性能等级的样本，通过调整关键变量使其转换到所期望的更优的性能等级；优选地，将主元中所对应的系数最大的变量设定为优化时所要调整的关键变量，然后根据关键变量与在线数据在置信椭圆中的投影点位置变化之间的相关性，确定调整方向，实现对生产过程的优化；优选地，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对常减压过程的生产模式进行调整，并通过实时监控观察数据是否往所期望的区域移动。
在一个或多个实施方案中，本发明的基于大数据的常减压过程模式监控和优化方法包括：将常减压过程中所采集到的大量历史数据，根据其不同物理意义进行分类，并定义为不同的生产模式；运用主成分分析法处理不同生产模式下的数据，利用得分矩阵建立模型并绘制置信椭圆；将实时采集得到的新样本计算得分矩阵投影到模型中，根据样本点在椭圆中的位置判断目前过程所处的生产模式；将投影到椭圆内的样本添加到历史数据中建立新模型，实现模型自适应更新；投影到椭圆外的样本则为异常点，并根据故障贡献率进行故障溯源；建模时可以得到不同生产模式下对过程状态变化起主要作用的变量，对这些变量进行调整实现模式优化。
在一个或多个实施方案中，本发明的第五个方面的步骤(1)-(8)如本发明的第一个方面的任一实施方案所述。
本发明的第六个方面提供一种基于大数据的加氢裂化过程模式识别和优化方法，所述方法包括以下步骤：
(1)利用加氢裂化过程的历史数据构建建模用的训练数据样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为样本集的样本个数，n为样本集的变量个数；
(2)对训练数据样本集进行预处理，得到方差为1、均值为0的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)运用主元分析法将X从n维降到k维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；
(5)采集在线实时数据Y∈R N×n，采用步骤(2)对训练样本进行预处理时所得到的样本平均值和样本方差对Y进行预处理得到标准化数据Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则说明此时加氢裂化过程所处的工况为正常工况；若样本点映射到椭圆外，则说明此时加氢裂化过程存在异常。
在一个或多个实施方案中，所述步骤(2)中的预处理方法如本发明的第一个方面的步骤(2)所述。
在一个或多个实施方案中，所述步骤(3)中，运用主元分析法对经过预处理所得到的X进行降维处理，具体操作步骤如本发明的第一个方面的步骤(3)所述。
在一个或多个实施方案中，所述步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆，具体操作步骤如本发明的第一个方面的步骤(4)所述。
在一个或多个实施方案中，所述步骤(4)中，将历史数据结合现场工艺知识可以对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分。在一些实施方案中，区域划分包括：首先根据工艺知识计算历史数据中各项性能指标(包括经济效益、生产能耗、产品产率等)的KPI值(效益值)，然后对KPI值依据实际工艺进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并加以不同的标记(例如颜色标记)。效益值在历史数据中可以查到或根据历史数据计算得到。
在一个或多个实施方案中，所述步骤(4)中，根据加氢裂化过程的性能指标，对历史数据设置不同的标签，因此在对历史数据的得分矩阵进行投影时可以通过不同区域数据点所带的标签对置信椭圆进行划分。
在一个或多个实施方案中，所述步骤(6)中，scorey的计算公式如本发明的第一个方面的步骤(6)所述。
在一个或多个实施方案中，所述步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一个或多个实施方案中，所述步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模。
在一个或多个实施方案中，所述步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源，具体操作方法如本发明的第一个方面的步骤(7)所述。
在一个或多个实施方案中，所述步骤(4)还包括：将历史数据结合现场工艺知识对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分，所述区域划分优选包括：首先根据工艺知识计算历史数据中性能指标的KPI值，然后对KPI值依据实际工艺进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并加以不同的标记。
在一个或多个实施方案中，所述方法还包括步骤(8)：对于映射在椭圆内部的样本，根据其所处位置结合步骤(4)中对椭圆进行的区域划分判断当前所处的性能等级。
在一个或多个实施方案中，所述步骤(8)还包括：对于处于非最优性能等级的样本，通过调整关键变量使其转换到所期望的更优的性能等级；优选地，将主元中所对应的系数最大的变量设定为优化时所要调整的关键变量，然后根据关键变量与在线数据在置信椭圆中的投影点位置变化之间的相关性，确定调整方向，实现对生产过程的优化；优选地，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对加氢裂化过程的生产模式进行调整，并通过实时监控观察数据是否往所期望的区域移动。
在一个或多个实施方案中，本发明的基于大数据的加氢裂化过程模式识别和优化方法包括：将加氢裂化过程中所采集到的大量历史数据，根据其不同物理意义进行分类，并定义为不同的生产模式；运用主成分分析法处理不同生产模式下的数据，利用得分矩阵建立模型并绘制置信椭圆；将实时采集得到的新样本计算得分矩阵投影到模型中，根据样本点在椭圆中的位置判断目前过程所处的生产模式；将投影到椭圆内的样本添加到历史数据中建立新模型，实现模型自适应更新；投影到椭圆外的样本则为异常点，并根据故障贡献率进行故障溯源；建模时可以得到不同生产模式下对过程状态变化起主要作用的变量，对这些变量进行调整实现模式优化。
在一个或多个实施方案中，本发明的第六个方面的步骤(1)-(8)如本发明的第一个方面的任一实施方案所述。
根据本发明提供的具体实施例，本发明公开了以下技术效果：采用大数据降维方法对炼油过程海量数据进行降维，获取可以表征炼油过程状态的关键变量，有效提高监控和优化方法的效率，降低计算代价。通过绘制置信椭圆对炼 油过程进行实时监测，采用这种可视化的方式，方便操作人员更加直观地了解过程所处的运行状态并进行相关操作。
说明书附图
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是催化重整工艺流程示意图。
图2是本发明的炼油过程模式识别及优化方法的总体流程框图。
图3是实施例1中正常数据所绘制的置信椭圆。
图4是实施例1中置信椭圆根据不同模式数据的分区图。
图5是实施例1中故障点的SPE贡献率图。
图6是实施例1中历史数据对应效益曲线。
图7是实施例1中历史数据在置信椭圆中的投影和高效益点对应的投影区间。
图8是实施例1中历史数据在置信椭圆中的投影和低效益点对应的投影区间。
图9是实施例1中操作变量优化移动轨迹。
图10是本发明的催化裂化装置运行状态识别及优化方法的总体流程框图。
图11是实施例2中根据催化裂化装置的生产数据所绘制的置信椭圆。
图12是实施例2中根据优化目标变量——经济效益折线选出的装置优化运行区域(折线中带圈的点)。
图13是实施例2中某种生产模式对应的优化区域数据在置信椭圆上的投影图。
图14是实施例2中变量对生产模式评价目标的贡献率图。
图15是实施例2中催化裂化装置调优轨迹。
图16是本发明的硫磺回收装置运行状态识别及优化方法的总体流程框图。
图17是实施例3中硫磺回收装置生产数据所绘制的置信椭圆。
图18是实施例3中根据优化目标变量折线(硫回收率)选出的装置优化 运行区域。
图19是实施例3中某种生产模式对应的优化区域数据在置信椭圆上的投影图。
图20是实施例3中变量对生产模式评价目标的贡献率图。
图21是实施例3中硫磺回收装置调优轨迹。
图22是本发明的基于数据驱动的渣油加氢处理装置运行状态识别及优化方法的总体流程框图。
图23是实施例4中根据渣油加氢处理装置生产数据所绘制的置信椭圆。
图24是实施例4中渣油加氢处理装置不同生产模式的分布。
图25是实施例4中注重加氢重油收率的生产样本对应的优化区域数据在置信椭圆上的投影图。
图26是实施例4中变量对生产模式评价目标的贡献率图。
图27是实施例4中渣油加氢处理装置调优轨迹。
图28是常减压工艺流程示意图。
图29是本发明的基于大数据的常减压过程模式监控和优化方法的总体流程框图。
图30是实施例5中根据正常数据所绘制的置信椭圆。
图31是实施例5中置信椭圆根据以能耗为优化目标的模式数据的分区图。
图32是实施例5中置信椭圆根据以轻收调优为优化目标的模式数据的分区图。
图33是实施例5中置信椭圆根据以总拔为优化目标的模式数据的分区图。
图34是实施例5中实时采集的数据在椭圆上的投影。
图35是实施例5中故障点的SPE贡献率图。
图36是加氢裂化的工艺过程流程图。
图37是本发明的基于大数据的加氢裂化过程模式识别和优化方法的总体流程框图。
图38是实施例6中根据正常数据所绘制的置信椭圆。
图39是实施例6中置信椭圆根据以总液收为优化目标的模式数据的分区图。
图40是实施例6中置信椭圆根据以中间油收率为优化目标的模式数据的 分区图。
图41是实施例6中置信椭圆根据以价值增量为优化目标的模式数据的分区图。
图42是实施例6中实时采集的数据在椭圆上的投影。
图43是实施例6中故障点的SPE贡献率图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明的基于大数据的炼油过程模式识别方法包括以下步骤：
(1)将炼油过程中所采集到的历史数据组成建模用的训练样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为样本集的样本个数，n为样本集的变量个数；
(2)对训练数据样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
(3)对X运用主元分析法将其从n维降到k维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k；
(4)利用得分矩阵T的前两列，绘制二维置信椭圆；
(5)采集新的在线实时数据Y∈R N×n，采用步骤(2)对训练样本进行预处理时所得到的样本平均值和样本方差对Y进行预处理得到标准化数据Ym∈R N×n；
(6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
(7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则说明此时炼油过程所处的工况为正常工况；若样本点映射到椭圆外，则说明此时炼油过程存在异常。
本发明中，炼油过程具有本领域常规的含义，通常是指石油炼制过程中原油及其中间产物经历的一系列工艺加工过程，包括但不限于常压蒸馏、减压蒸馏、催化裂化、催化重整、加氢裂化、硫磺回收、渣油加氢、常减压、延迟焦 化、炼厂气加工和烷基化等。在一些实施方案中，炼油过程为催化重整过程、催化裂化过程、硫磺回收过程、渣油加氢过程、常减压过程或加氢裂化过程。
本发明中，炼油过程模式包括炼油过程的工况和生产模式。工况可以分为正常工况(无故障发生的工况)和异常工况(有故障发生的工况)。正常工况可以根据不同的工艺和性能评价指标划分成的不同生产模式，例如经济效益模式、装置能耗模式、产品产率模式等。根据各种生产模式所对应的不同评价指标，可将各生产模式划分为优化运行状态和非优化运行状态，例如经济效益模式可划分为高经济效益模式和低经济效益模式，装置能耗模式可划分为高装置能耗模式和低装置能耗模式，产品产率模式可划分为高产品产率模式和低产品产率模式等。
步骤(1)中，采用炼油过程中所采集到的历史数据作为训练样本集。本发明中，一个样本通常是指在某个时间点对所选取的各个变量进行采集得到的一组数据，包括多个变量。当样本集以矩阵的形式呈现时，矩阵中的一行数据(行向量)通常表示一个样本，矩阵中的各列数据(列向量)通常分别对应于不同的变量，变量的个数(列数)即为数据的维度。变量包括自变量(本发明又称操作变量)和因变量。训练样本集可以来自工厂实时数据库中的历史数据。通常，训练样本为正常运行状态下的历史数据。可根据装置数据记录点位和性质分析项目，批量采集工厂生产过程历史数据，根据先验知识粗略筛选出数据正常、对装置生产过程可能有影响的操作点(操作数据)和原料性质项(原料数据)，将包含这类操作数据和原料数据的催化裂化装置生产过程历史数据作为建模用的训练数据。
步骤(2)中，通过求取训练样本集的平均值和标准差来对训练样本进行标准化处理，得到标准化的数据。本发明中，基于平均值和标准差对数据进行标准化处理的方法是本领域常规的。标准化后的数据呈现均值为0、方差为1的正态分布。在一些实施方案中，预处理方法采用Z-score标准化方法，计算公式为：
Figure PCTCN2023070795-appb-000021
式中，Z＝[z 1,z 2,...,z m]为训练数据矩阵，X表示标准化后的数据矩阵，μ为训练数据的均值，σ为训练数据的标准差，μ和σ计算公式为：
Figure PCTCN2023070795-appb-000022
Figure PCTCN2023070795-appb-000023
步骤(3)中，利用主元分析法对X进行降维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k可以采用本领域已知的方法或本发明提供的方法。在一些实施方案中，对经过预处理所得到的X运用主元分析法进行降维处理按如下步骤进行：
(3-a)计算矩阵X的协方差矩阵，公式为：
Figure PCTCN2023070795-appb-000024
X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，所以得到的协方差矩阵C是一个n×n维的矩阵；
(3-b)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序，计算公式为：
det(C-Iλ)＝0                      (5)
Cp i＝λ ip i                        (6)
式(5)中，λ为特征值集合，即λ＝λ 1,λ 2,...,λ n，I为单位矩阵，得：
特征值矩阵：
Figure PCTCN2023070795-appb-000025
其中λ 1＞λ 2＞...＞λ n
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
(3-c)建立主元模型，选取主元个数原则：将最大的方差归纳到模型空间，将最小的方差留给噪声空间，计算公式为：
每个主元包括的信息比：
Figure PCTCN2023070795-appb-000026
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000027
其中，PC i为第i个主元包括的信息比，PC k为最大的k个主元的信息比。
(3-d)保留对应k个最大特征值的特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                     (9)
其本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
步骤(4)中，利用得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2绘制置信椭圆。在一些实施方案中，利用矩阵xdat绘制置信椭圆的步骤如下：
(4-a)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；xdat的协方差矩阵可参照前述公式(4)计算；
(4-b)计算xdat每一列的平均值得xm∈R 1×2，并将xdat每一列的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
(4-c)计算公式xd×s·×xd，并对所得矩阵每行进行求和，得到rd∈R m×1；
(4-d)绘制xdat的曲线图，得分矩阵的前两列呈现经验分布，根据其经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度Confi，对rd进行升序排序后，以rd前95％的数据位置所对应的值作为r∈R，优选的Confi为95％；
(4-e)计算步骤(4-a)所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000028
特征向量矩阵为
Figure PCTCN2023070795-appb-000029
矩阵S的特征值和对应的特征向量可参照前述公式(5)和公式(6)计算；
(4-f)利用步骤(4-d)中的r以及步骤(4-e)中的D，即可得到置信椭圆的长轴以及短轴，计算公式为：
Figure PCTCN2023070795-appb-000030
Figure PCTCN2023070795-appb-000031
其中，a为长轴，b为短轴；
(4-g)以步骤(4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴即可绘制置信椭圆，椭圆的公式为：
Figure PCTCN2023070795-appb-000032
其中xm1为xdat第一列的平均值，xm2为xdat第二列的平均值。
不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域，根据炼油过程的性能指标，可以根据历史工况对应的效益分级对历史数据设置不同的标签。效益值在历史统计数据中可以查到。将历史数据的得分矩阵的前两列进行投影，可以通过不同区域数据点所带的标签对置信椭圆进行划分。
步骤(4)中，利用得分矩阵T的前两列，绘制二维置信椭圆。置信度优选设置为95％。可以根据历史数据和生产需求，完成与生产模式定义有关的关键表征指标变量设定和数据整理，如经济效益、装置能耗、产品产率等，采用 作图法，根据需求对生产模式指标的等级进行划分，筛选出处于优化状态数据，将此类数据对应的样本在二维置信椭圆上完成投影得到定义生产模式下的优化区域，并以不同的颜色加以标记。
在一些实施方案中，步骤(4)还包括对置信椭圆按照性能等级进行区域划分。在一些实施方案中，区域划分包括：将历史数据结合现场工艺知识对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分，所述区域划分优选包括：首先根据工艺知识和历史效益统计数据，根据时间序列找到工况对应的效益，根据效益高低进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并用不同的颜色加以标记。在一些实施方案中，区域划分包括：根据炼油过程的性能指标，对不同性能等级的历史数据设置不同的标签，将历史数据的得分矩阵的前两列投影到置信椭圆上，通过不同区域数据点所带的标签对置信椭圆进行划分。
步骤(5)中，参照公式(1)采用步骤(2)中对训练样本进行预处理时所得到的均值和方差对Y进行预处理。
步骤(6)中，scorey的计算公式为：
scorey＝Ym×[p 1,p 2]                 (13)
式中，[p 1,p 2]∈R n×2是负载矩阵P的前两列。
步骤(7)中，可以通过判断样本点是否映射在椭圆内，考虑是否将此样本点视为装置正常工况纳入到历史数据中重新建模，提高模型的适应性；若样本映射点在此椭圆内，则纳入数据库中更新模型，否则表明此时装置生产过程可能存在异常，可对异常数据计算SPE贡献率来进行故障溯源。
步骤(7)中，可将scorey中的每一组数据代入前述公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
在一些实施方案中，步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模；优选地，在加入收集到的正常工况下的数据的同时，剔除相同数量的早期历史数据。
在一些实施方案中，步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源。计算SPE贡献率可采用本领域已知的方法或本发明提供的方法。在一些实施方案中，假设异常数据为x∈R 1×n，其SPE贡献率的计算公式为：
Figure PCTCN2023070795-appb-000033
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵 的第i列，T表示转置，
Figure PCTCN2023070795-appb-000034
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量。
在一些实施方案中，步骤(7)还包括：对于映射在椭圆内部的样本，根据其所处位置结合步骤(4)中对椭圆进行的区域划分判断当前所处的性能等级，若其处于非最优的性能等级，则通过调整关键变量使其转换到所期望的更优的性能等级；优选地，将主元中所对应的系数最大的变量设定为优化时所要调整的关键变量，然后根据关键变量与在线数据在置信椭圆中的投影点位置变化之间的相关性，确定调整方向，实现对生产过程的优化；优选地，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对炼油过程的生产模式进行调整，并通过实时监控观察数据是否往所期望的区域移动。在这类实施方案中，本发明的方法可以实现对炼油过程模式的优化。
在本发明的基于大数据的炼油过程模式识别方法的基础上引入如下的步骤(8)也可以实现对炼油过程模式的优化：
(8)通过主成分分析过程中得到的标准差和均值求得置信椭圆中的点对应的原始变量，从而得到该点对应的效益值；在置信椭圆中区分效益值高低的分布；当工况运行模式位于置信椭圆中的某一点时，采用路径优化算法，获得当前位置向最优位置最快移动轨迹，将轨迹对应的点进行反变换，获得操作条件的改变方式，从而指导生产装置操作优化。
因此，本发明还包括基于大数据的炼油过程模式优化方法，所述基于大数据的炼油过程模式优化方法包括如本发明任一实施方案所述的炼油过程模式识别方法和如本发明任一实施方案所述的步骤(8)。
置信椭圆中，每一个点都可以通过主成分分析过程中得到的标准差和均值求得对应的原始变量，从而可以得到该点对应的效益值。效益值高低的分布即可在置信椭圆中进行区分。当工况运行模式位于置信椭圆中的某一点时，采用路径优化算法，获得当前位置向最优位置最快移动轨迹，将轨迹对应的点进行反变换，获得操作条件的改变方式，即可指导生产装置操作优化。
在一些实施方案中，步骤(8)中，针对置信椭圆中的当前工况点defen＝(x,y)∈R 1×2，利用公式(15)计算原始变量X ori：
X ori＝(defen×PC T)·×std(X m)+mean(X m)           (15)
其中，defen＝(x,y)∈R 1×2为置信椭圆中的点的横纵坐标，PC＝[p 1,p 2]∈R J×2为主成分分析建模时负载矩阵的前两列，std(X m)为主成分分析采用的基础样本的方差，mean(X m)为主成分分析采用的基础样本的均值，m为样本个数，X ori为椭圆中defen＝(x,y)∈R 1×2位置对应的反变换后的原始变量。
生产过程中的统计信息包括进料量、产品产量等，根据价格信息即可得到简化的投入产出效益指标。在一些实施方案中，步骤(8)中，根据公式(16)计算投入产出效益指标：
Figure PCTCN2023070795-appb-000035
其中，
Figure PCTCN2023070795-appb-000036
为第i种产品的产量，
Figure PCTCN2023070795-appb-000037
为第i种产品的价格，
Figure PCTCN2023070795-appb-000038
为第j种原料的投料量，
Figure PCTCN2023070795-appb-000039
为第j种原料的价格，profit为效益指标。
置信椭圆中的点均可通过式(15)的变换后得到原始工况并根据式(16)得到一个对应的效益值，同时可以找到效益最优值对应的点。根据效益值可以在置信椭圆中划分高效益区和低效益区。
步骤(8)中，路径优化可以采用本领域已知的算法(例如常规的A*算法)或采用本发明提供的改进的A*算法。A*算法是一种确定起始点后，找到从起始点到目标点所有可达路径中的最优路径的优化算法。常规的A*算法通过评估函数f(x)来确定搜索方向和下一个到达的节点，评估函数的一般形式如下所示：
f(x)＝g(x)+h(x)                   (17-1)
其中，g(x)是代价函数，代表从起始节点到达当前节点x的实际所需代价；h(x)是启发式函数，代表从当前节点x到达目标节点估计所需的代价。
在优选的实施方案中，采用本发明的改进的A*算法进行路径优化，所述改进的A*算法对评估函数f(x)进行了改进：
Figure PCTCN2023070795-appb-000040
其中，profit(x)表示所选节点对应的经济效益，g(x)是代价函数，代表从起始节点到达当前节点x的实际所需代价；h(x)是启发式函数，代表从当前节点x到达目标节点估计所需的代价。采用本发明的改进的A*算法进行路径优化可以使优化过程中保持较高的经济效益。采用本发明的改进的A*算法可以获得带优化的工况点到高效益点的最优路径。
获得优化路径后，可根据采用式(15)进行变换而得到的操作变量的原始值进行优化操作。
本发明的有益效果如下：
1.炼油生产中包含大量过程数据，例如温度、压力、液位、流量以及物料性质等，这些变量都会对生产过程造成一定的影响。但是在海量的数据中，有些数据对生产过程的影响较为显著，有些则影响较小。采用大数据降维方法对炼油过程海量数据进行降维，获取可以表征炼油过程状态的关键变量，可以 有效提高监控和优化方法的效率，降低计算代价。
2.在建模的过程中，不断地将实时采集到的正常数据加入到建模数据中，并剔除早期的历史数据使建模数据总量保持不变，提高了模型适应性，使模型随着过程的运行进行调整，所监控的关键变量也随之改变。
3.通过绘制置信椭圆对炼油过程进行实时监测，采用这种可视化的方式，方便操作人员更加直观地了解过程所处的运行状态并进行相关操作。
4.正常状态下的炼油过程根据其工艺工况和性能指标可以分为不同的生产模式，例如高经济效益模式或低经济效益模式，不同生产模式的数据会映射在置信椭圆内不同的区域，通过将这些模式下的历史数据映射在椭圆中对椭圆进行划分，在实时监控时，通过实时采集的数据映射在椭圆中的位置来判断当前炼油过程所处的生产模式，并且可以通过路径优化和模式点反变换得到的操作变量调整方式指导炼油装置生产操作优化。
5.本发明采用A*算法作为优化算法，可以在优劣区间找到一条最优路径，同时可以通过路径上的点反推出最优操作条件。
下面通过实施例对本发明进行具体描述。有必要在此指出的是，以下实施例只用于对本发明作进一步说明，不能理解为对本发明保护范围的限制，本领域技术人员根据本发明的内容作出的一些非本质的改进和调整，仍属于本发明的保护范围。
实施例1
本实施例将本发明的基于大数据的炼油过程模式优化方法应用到催化重整(continuous catalytic reforming，CCR)工艺过程中。图1给出了催化重整的工艺过程流程图。催化重整工艺由预加氢单元、重整单元以及催化剂再生系统组成，在以生产芳烃为目的时，还包括芳烃抽提和精馏装置。经过预处理后的原料进入重整工段，与循环氢混合并加热之后进入反应器，反应器由3至4个串联，其间设有加热炉，以补偿反应所吸收的热量。离开反应器的物料进入分离器分离出氢循环气(多余部分排出)，所得液体由稳定塔脱去轻组分后作为重整汽油，是高辛烷值汽油组分，或送往芳烃抽提装置生产芳烃。
CCR模型包括51个输入变量和16个输出变量。通过对输入变量进行控制可以得到不同运行状态的CCR模型，其中对输出结果影响较大的输入变量有反应温度、压力、进料流量等，产物的性能指标主要通过观察输出变量中的氢气、纯氢、干气、液化气、碳5、碳6、碳7+、苯质量流量以及芳烃的量来进行判断。表1列出了本实施例所涉及的操作变量。
表1：操作变量
变量 描述 单位
1 一段反应温度 ℃
2 二段反应温度 ℃
3 三段反应温度 ℃
4 四段反应温度 ℃
5 循环氢量 STDm 3/h
6 预加氢进料总量 tonne/h
7 压缩机压力 Mpa
8 T201塔板温度 ℃
9 T201回流量 tonne/h
10 T201塔板温度 ℃
11 T601塔底温度 ℃
12 抽出二甲苯 tonne/h
13 重整进料负荷 tonne/h
采用本发明的基于大数据的炼油过程模式优化方法对CCR过程进行监控和优化，如图2所示，包括以下步骤：
1.采集一段时间CCR过程正常运行状态下的样本数据Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本。
2.对采集的数据进行预处理，得到均值为0、方差为1的标准化数据集X＝[x 1,x 2,...,x n]∈R m×n，计算公式为：
Figure PCTCN2023070795-appb-000041
式中，μ为取自训练数据的均值：
Figure PCTCN2023070795-appb-000042
σ为取自训练数据的标准差：
Figure PCTCN2023070795-appb-000043
3.对经过标准化处理后所得到的X运用主元分析法进行降维处理，具体 按如下步骤进行：
a)计算矩阵X的协方差矩阵，公式为：
Figure PCTCN2023070795-appb-000044
X是m×n的矩阵，m是样本个数，n是变量个数，所以得到的协方差矩阵C是一个n×n维的矩阵；
b)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序，计算公式为：
det(C-Iλ)＝0                       (5)
Cp i＝λ ip i                         (6)
式(5)中，I为单位矩阵，得：
特征值矩阵：
Figure PCTCN2023070795-appb-000045
其中λ 1＞λ 2＞...＞λ n
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)建立主元模型，选取主元个数原则：将最大的方差归纳到模型空间，将最小的方差留给噪声空间，计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000046
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000047
d)保留对应k个最大特征值的特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
其本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用该矩阵具体绘制置信度为95％的置信椭圆的步骤如下：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，因为所要绘制的置信椭圆的置信度为95％，对rd进行升序排序后，将rd第95％个位置所对应的值作为r∈R；
e)计算4-a)中所得矩阵s的特征值和对应的特征向量，计算公式与式(5)和式(6)类似，得特征值矩阵为
Figure PCTCN2023070795-appb-000048
特征向量矩阵为
Figure PCTCN2023070795-appb-000049
f)利用4-d)中的r以及4-e)中的D，可求得置信椭圆的长轴以及短轴，计算公式为：
Figure PCTCN2023070795-appb-000050
Figure PCTCN2023070795-appb-000051
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴即可绘制置信椭圆，椭圆的公式为：
Figure PCTCN2023070795-appb-000052
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图3所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域，根据CCR过程的性能指标，根据历史工况对应的效益分级对历史数据设置不同的标签，效益值在历史统计数据中可以查到，将历史数据的得分矩阵的前两列进行投影，可以通过不同区域数据点所带的标签对置信椭圆进行划分，并以不同颜色进行区分如图4所示。
6.将CCR过程在正常工况下运行一段时间后通过调整操作变量中压缩机压力的值使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用训练数据计算所得的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.将Ym乘以3-d)所得负载矩阵P的前两列，计算公式为：
scorey＝Ym×[p 1,p 2]                  (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
9.假设故障数据点为x∈R 1×n，其SPE贡献率的计算公式为：
Figure PCTCN2023070795-appb-000053
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000054
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡献率如图5所示，根据图5可以判定故障产生的原因主要是由第七个变量即压缩机压力起主要作用，与实际操作情况相符。
10.催化重整生产过程的经济效益分别与液化气流量(x 1)、氢气流量(x 2)、干气流量(x 3)、碳5(x 4)、碳6(x 5)、芳烃(x 6)以及重整进料负荷(x 7)有关，这些物质的价格如表2所示：
表2：催化重整经济效指标各物质价格
物质/m 3·s -1 进料负荷 液化气 氢气 干气 碳5 碳6 芳烃
价格/元 3092 4153 8489.57 2505 3690 4298 4573
因此催化重整生产过程经济效益的计算公式可以表示如下：
profit＝(4153x 1+8489.57x 2+2505x 3+3690x 4+4298x 5+4573x 6)-3092x 7  (15)
根据历史样本即可绘制如图6所示效益值，从而可以在置信椭圆中划分高效益区和低效益区，分别如图7和图8所示。图7中，浅色点所在位置代表高效益区。图8中，浅色点所在位置代表低效益区。
针对置信椭圆中的当前工况点defen＝(x,y)∈R 1×2，可以采用下式获得原始变量：
X ori＝(defen×PC T)·×std(X m)+mean(X m)           (16)
其中，defen＝(x,y)∈R 1×2为置信椭圆中的点的横纵坐标，PC＝[p 1,p 2]∈R J×2为主成分分析建模时负载矩阵的前两列，std(X m)为主成分分析采用的基础样本的方差，mean(X m)为主成分分析采用的基础样本的均值，m为样本个数，X ori为椭圆中defen＝(x,y)∈R 1×2位置对应的反变换后的原始变量。
为了获得该变量到高效益点的最优路径，采用改进的A*算法实现路径优化。该改进方法的评估函数如下所示：
Figure PCTCN2023070795-appb-000055
其中，g(x)是代价函数，它代表从起始节点到达当前节点x的实际所需代价；h(x)是启发式函数，代表从当前节点x到达目标节点估计所需的代价，中profit(x)表示所选节点对应的经济效益，可以使优化过程中保持较高的经济效益。路径优化后得到的模式移动轨迹如图9所示，采用式(15)进行变换即可得到操作变量的原始值。如表3所示。采用上述模式优化方法得到的操作条件进行优化操作，可以将效益从257584元/小时提高到259272元/小时。
表3：优化过程中对应操作变量的值和效益结果
Figure PCTCN2023070795-appb-000056
实施例2
本实施例根据某180万吨/年的工业催化裂化生产历史数据，如表4所示，选择88个装置自变量数据点作为模型输入变量和20个装置因变量数据点作为模型对应的输出变量，对生产历史数据进行采集和预处理。
本实施例采用本发明的一种基于数据驱动的催化裂化装置运行状态识别及优化方法对催化裂化装置运行状态进行识别及优化，如图10所示，包括以下步骤：
1.采集装置实际生产历史数据，相关变量名称见表4。利用催化裂化装置正常运行状态下大量数据组合构建样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本。
表4：催化装置模型变量名称
Figure PCTCN2023070795-appb-000057
Figure PCTCN2023070795-appb-000058
2.对样本集数据进行预处理，得到方差为1、均值为0的标准数据集X＝[x 1,x 2,...,x n]∈R m×n，标准数据集的计算公式为：
Figure PCTCN2023070795-appb-000059
式中，μ为训练数据的均值，其计算公式为：
Figure PCTCN2023070795-appb-000060
σ为训练数据的标准差，其计算公式为：
Figure PCTCN2023070795-appb-000061
3.按如下步骤进行对标准数据集X运用主元分析法进行降维：
a)根据公式(4)计算矩阵X的协方差矩阵：
Figure PCTCN2023070795-appb-000062
其中，X是m×n的矩阵，m是样本个数，n是变量个数，因此协方差矩阵C 是一个n×n维的矩阵；
b)根据公式(5)和公式(6)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序：
det(C-Iλ)＝0                        (5)
Cp i＝λ ip i                          (6)
式(5)中，I为单位矩阵；
特征值矩阵：
Figure PCTCN2023070795-appb-000063
其中λ 1＞λ 2＞...＞λ n；
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)建立主元模型，选取主元个数原则：将最大的方差归纳到模型空间，将最小的方差留给噪声空间，其计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000064
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000065
d)保留与k个最大特征值相对应的k个特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，其中特征向量按照特征值降序排序，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
t i的本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，按照以下步骤利用该矩阵具体绘制置信度为95％的置信椭圆：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据得分矩阵的前两列呈现经验分布的特点，计算 矩阵rd的百分位数，由于所要绘制的置信椭圆的置信度为95％，对rd求第5％个位置所对应的值，得到r∈R；
e)参照式(5)和式(6)计算4-a)中所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000066
特征向量矩阵为
Figure PCTCN2023070795-appb-000067
f)利用4-d)中r以及4-e)中D，按照公式(10)和公式(11)计算置信椭圆的长轴以及短轴：
Figure PCTCN2023070795-appb-000068
Figure PCTCN2023070795-appb-000069
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，置信椭圆的公式为：
Figure PCTCN2023070795-appb-000070
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图11所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域。如图12所示，根据催化裂化过程的性能指标，对历史数据设置不同的标签，进行优化目标变量——经济效益折线绘制，从中选出装置的优化运行区域。将历史数据的得分矩阵的前两列进行投影，通过不同区域数据点所带的标签对置信椭圆进行划分，并以不同颜色进行区分。图13给出了本实施例的置信椭圆区域划分结果，浅色点所在区域为经济效益最优的生产模式对应的区域。
6.将催化裂化装置在正常工况下运行一段时间后，通过调整操作变量中烧焦罐下部密度使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用训练数据计算所得的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.按照公式(13)将Ym乘以3-d)所得负载矩阵P的前两列，得到scorey：
scorey＝Ym×[p 1,p 2]                  (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内；
9.设故障数据点为x∈R 1×n，根据公式(14)计算故障数据点的SPE贡献率：
Figure PCTCN2023070795-appb-000071
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000072
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡 献率如图14所示。根据图14可以判定故障产生的原因主要是由第10个变量(即烧焦罐下部密度)起主要作用，与步骤(6)的实际操作情况相符。
10.根据负载矩阵P的第一列中绝对值最大的烧焦罐下部密度变量，通过调整该变量所对应的值，对催化裂化过程的生产状态进行调整，并通过实时监控可以观察数据是否往期望的区域移动，从而使装置运行状态回到优化运行状态，调整轨迹如图15所示。图15中，浅色点所在区域对应于优化运行状态。
实施例3
本实施例根据某60万吨/年的工业硫磺回收生产历史数据，如表5所示，选择75个装置数据点，其中55个为模型输入变量，20个作为输出变量，对生产历史数据进行采集和预处理。
本实施例采用本发明的一种基于数据驱动的硫磺回收装置运行状态识别及优化方法对硫磺回收装置运行状态进行识别和优化，如图16所示，包括以下步骤：
1.采集硫磺回收装置实际生产历史数据，部分自变量名称见表5。利用硫磺回收装置正常运行状态下大量数据组合构建样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本。
表5：硫磺装置变量名称
Figure PCTCN2023070795-appb-000073
Figure PCTCN2023070795-appb-000074
2.对样本集进行标准化处理，得到方差为1、均值为0的标准化处理后的数据集X＝[x 1,x 2,...,x n]∈R m×n，标准化处理的计算公式为：
Figure PCTCN2023070795-appb-000075
式中，μ为训练数据的均值，其计算公式为：
Figure PCTCN2023070795-appb-000076
σ为训练数据的标准差，其计算公式为：
Figure PCTCN2023070795-appb-000077
3.按如下步骤进行对标准化处理后的数据集X运用主元分析法进行降维：
a)根据公式(4)计算矩阵X的协方差矩阵：
Figure PCTCN2023070795-appb-000078
其中，X是m×n的矩阵，m是样本个数，n是变量个数，因此协方差矩阵C是一个n×n维的矩阵；
b)根据公式(5)和公式(6)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序：
det(C-Iλ)＝0                       (5)
Cp i＝λ ip i                         (6)
式(5)中，I为单位矩阵；
特征值矩阵：
Figure PCTCN2023070795-appb-000079
其中λ 1＞λ 2＞...＞λ n；
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)按照将最大的方差归纳到模型空间、将最小的方差留给噪声空间的原则选取主元个数k，其计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000080
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000081
d)保留与k个最大特征值相对应的k个特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
t i的本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，按照以下步骤利用该矩阵 具体绘制置信度为95％的置信椭圆：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度为95％，对rd进行升序排序后，将rd第95％个位置所对应的值作为r∈R；
e)参照式(5)和式(6)计算4-a)中所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000082
特征向量矩阵为
Figure PCTCN2023070795-appb-000083
f)利用4-d)中r以及4-e)中D，按照公式(10)和公式(11)计算置信椭圆的长轴以及短轴：
Figure PCTCN2023070795-appb-000084
Figure PCTCN2023070795-appb-000085
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，置信椭圆的公式为：
Figure PCTCN2023070795-appb-000086
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图17所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域。如图18所示，根据硫磺回收过程的装置特性，对历史数据设置硫回收率最大的优化目标进行标签后，进行相应优化目标变量折线绘制，从中选出装置的优化运行区域(硫回收率>99.5％)。将历史数据的得分矩阵的前两列进行投影，通过不同区域数据点所带的标签对置信椭圆进行划分，并以不同颜色进行区分。图19给出了本实施例的置信椭圆区域划分结果，浅色点所在区域为追求硫回收率最高的生产模式对应的区域。
6.将硫磺回收装置在正常工况下运行一段时间后，通过调整操作变量中制硫炉温度使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用步骤(2)所得的训练数据的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.按照公式(13)将Ym乘以3-d)所得负载矩阵P的前两列，得到scorey：
scorey＝Ym×[p 1,p 2]                  (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内；
9.设故障数据点为x∈R 1×n，根据公式(14)计算故障数据点的SPE贡献率：
Figure PCTCN2023070795-appb-000087
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000088
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡献率如图20所示。根据图20可以判定故障产生的原因主要是由第11个变量(即制硫炉温度)起主要作用，与实际操作情况相符。
10.根据负载矩阵P的第一列中绝对值最大的制硫炉温度变量，通过调整该变量所对应的值，对硫磺回收过程的生产状态进行调整，并通过实时监控可以观察数据是否往期望的区域移动，从而使硫磺回收装置运行状态回到优化运行状态，调整轨迹如图21所示。图21中，浅色点所在区域对应于优化运行状态。
实施例4
本实施例以170万吨/年的渣油加氢处理装置为识别和优化对象，采集2018年5月1日至2019年5月1日的历史生产数据，并剔除异常波动数据。经筛选处理后，将数据整理为85变量3925组。表6列出了本实施例涉及的变量。
表6：实施例4涉及的变量
Figure PCTCN2023070795-appb-000089
Figure PCTCN2023070795-appb-000090
采用本发明的基于数据驱动的渣油加氢处理装置运行状态识别及优化方法对该渣油加氢处理装置的运行状态进行识别和优化，如图22所示，包括以下步骤：
1.采集渣油加氢处理装置实际生产历史数据。利用渣油加氢处理装置正常运行状态下大量数据组合构建样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本。
2.按照公式(1)对样本集数据进行标准化处理，得到均值为0、方差为1的标准化处理后的数据集X＝[x 1,x 2,...,x n]∈R m×n：
Figure PCTCN2023070795-appb-000091
式中，μ为训练数据的均值，其计算公式为：
Figure PCTCN2023070795-appb-000092
σ为训练数据的标准差，其计算公式为：
Figure PCTCN2023070795-appb-000093
3.按如下步骤进行对经过标准化处理后的数据集X运用主元分析法进行降维：
a)根据公式(4)计算矩阵X的协方差矩阵：
Figure PCTCN2023070795-appb-000094
其中，X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，因此协方差矩阵C是一个n×n维的矩阵；
b)根据公式(5)和公式(6)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序：
det(C-Iλ)＝0                      (5)
Cp i＝λ ip i                        (6)
式(5)中，I为单位矩阵；
特征值矩阵：
Figure PCTCN2023070795-appb-000095
其中λ 1＞λ 2＞...＞λ n；
特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)按照将最大的方差归纳到模型空间、将最小的方差留给噪声空间的原则 选取主元个数k，其计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000096
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000097
d)保留与k个最大特征值相对应的k个特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
其本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用该矩阵具体绘制置信度为95％的置信椭圆的步骤如下：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度为95％，对rd进行升序排序后，将rd第95％个位置所对应的值作为r∈R；
e)参照式(5)和式(6)计算4-a)中所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000098
特征向量矩阵为
Figure PCTCN2023070795-appb-000099
f)利用4-d)中的r以及4-e)中的D，按照公式(10)和公式(11)计算置信椭圆的长轴以及短轴：
Figure PCTCN2023070795-appb-000100
Figure PCTCN2023070795-appb-000101
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，置信椭圆的公式为：
Figure PCTCN2023070795-appb-000102
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图23所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域。如图24所示，根据渣油加氢处理过程的不同性能指标，将历史数据设置以加氢重油收率最大为目标的模式、以装置效益最大为目标的模式和以自装置外来新氢体积流量最小为目标分为三种模式。通过将历史数据的得分矩阵的前两列进行投影，通过不同区域数据点所带的标签对置信椭圆进行划分。图25中浅色点代表注重加氢重油收率的生产样本对应的优化区域数据在置信椭圆上的投影位置，这些位置即为以加氢重油收率最大为目标的模式所对应的优化区域。
6.待渣油加氢处理装置在正常工况下运行一段时间后，通过调整操作变量中罐区冷渣流量使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用根据训练数据计算所得的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.按照公式(13)将Ym乘以3-d)所得负载矩阵P的前两列，得到scorey：
scorey＝Ym×[p 1,p 2]                 (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内；
9.设故障数据点为x∈R 1×n，根据公式(14)计算故障数据点的SPE贡献率：
Figure PCTCN2023070795-appb-000103
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000104
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡献率如图26所示。根据图26可以判定故障产生的原因主要是由第20个变量(即罐区冷渣流量)起主要作用，与实际操作情况相符。
10.根据负载矩阵P的第一列中绝对值最大的罐区冷渣流量变量，通过调整该变量所对应的值，可以对渣油加氢处理过程的生产状态进行调整，并通过实时监控可以观察数据是否往期望的区域移动，从而使装置运行状态回到优化运行状态，调整轨迹见图27中的浅色点。
实施例5
本实施例中，基于大数据的常减压过程模式监控和优化方法被应用于某炼油厂的常减压实际工艺过程中。图28给出了常减压的工艺过程流程简图。常减压工艺过程分为五个部分，分别是电脱盐部分、初馏部分、常压部分、减压 部分、轻烃回收部分，其中最重要的是初馏、常压和减压这三个部分。原油初馏主要为了拔出部分轻馏分，为后续常减压装置分担一定的装置压力，能够降低整套常减压装置的能耗。常压蒸馏装置的主要作用则是拔出沸点较低的石脑油以及煤油和柴油等馏分。从常压塔底馏出的较为重质的常压渣油，会被送至减压塔进行减压蒸馏，分离出润滑油、蜡油、沥青等二次加工的原料。减压塔底馏出的减压渣油会被送往催化裂化、延迟焦化等装置进行深加工。
常减压模型包括60个输入变量和20个输出变量。其中输入变量包括52组装置操作变量和8组原料性质变量。通过对输入变量进行控制可以得到不同运行状态的CDU模型，其中对输出结果影响较大的输入变量有反应温度、压力、进料性质等，产物的性能指标主要通过观察输出变量来进行判断。如在以总拔收率为目标的生产模式中，总拔收率主要由石脑油、常一线、柴油、蜡油和减四线的收率共同组合，而这些收率的计算可以通过输出变量中的气体质量流量、初顶油出装置流量、常顶油出装置流量、常一线抽出流量、常二线抽出流量、常三线抽出流量、减一线抽出流量、冷蜡油抽出流量、去I催蜡油流量、去Ⅱ催蜡油流量、减四线抽出流量、渣油抽出流量来进行计算。表7列出了本实施例涉及的所有操作变量和原料性质变量。
表7：实施例5涉及的操作变量
变量 描述 单位 变量 描述 单位
1 脱后原油一路流量 tone/h 31 密度 g/cm 3
2 脱后原油二路流量 tone/h 32 IBP-80收率 ％
3 脱后原油进初馏塔温度 C 33 80-120收率 ％
4 原油进装置压力 MPa 34 120-180收率 ％
5 常压塔进口温度 C 35 180-240收率 ％
6 减压炉总出口温度 C 36 240-300收率 ％
7 减压塔闪蒸段温度 C 37 300-350收率 ％
8 初馏塔顶压力 MPa 38 350-400收率 ％
9 初馏塔塔底温度 C 39 400-450收率 ％
10 初顶回流罐温度 C 40 450-500收率 ％
11 常顶压力 MPa 41 500-540收率 ％
12 常压塔底温度 C 42 540-FBP收率 ％
13 常顶回流罐温度 C 43 酸值 ％
14 减压塔顶压力 MPa 44 碳含量 ％
15 减压塔顶温度 C 45 氢含量 ％
16 减压塔底温度 C 46 硫含量 ppm
17 初侧油出装置流量 tone/h 47 氮含量 ppm
18 初顶冷凝器压力 MPa 48 残炭 ％
19 初顶循抽出温度 C 49 凝点 C
20 初顶循返回温度 C 50 胶质 ％
21 初顶循抽出流量 tone/h 51 沥青质 ％
22 常顶循抽出流量 tone/h 52 H2O ％
23 常一中抽出流量 tone/h 53 初顶瓦斯比重 No
24 常二中抽出流量 tone/h 54 常顶瓦斯比重 No
25 常顶循抽出温度 C 55 减顶瓦斯比重 No
26 常顶循返回温度 C 56 初顶终馏点 C
27 常一中抽出温度 C 57 常顶终馏点 C
28 常一中返回温度 C 58 常一线终馏点 C
29 常二中抽出温度 C 59 常二线终馏点 C
30 常二中返回温度 C 60 常三线95％回收温度 C
采用本发明的基于大数据的常减压过程模式监控和优化方法对该常减压过程进行模式监控和优化，如图29所示，包括以下步骤：
1.采集一段时间CDU过程正常运行状态下的样本数据Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本；
2.按照公式(1)对采集的数据进行预处理，得到均值为0、方差为1的标准数据集X＝[x 1,x 2,...,x n]∈R m×n：
Figure PCTCN2023070795-appb-000105
式中，μ为采集的数据的均值，其计算公式为：
Figure PCTCN2023070795-appb-000106
σ为采集的数据的标准差，其计算公式为：
Figure PCTCN2023070795-appb-000107
3.按如下步骤进行对经过标准数据集X运用主元分析法进行降维：
a)根据公式(4)计算矩阵X的协方差矩阵：
Figure PCTCN2023070795-appb-000108
其中，X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，因此协方差矩阵C是一个n×n维的矩阵；
b)根据公式(5)和公式(6)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序：
det(C-Iλ)＝0                        (5)
Cp i＝λ ip i                          (6)
式(5)中，I为单位矩阵；
得到特征值矩阵：
Figure PCTCN2023070795-appb-000109
其中λ 1＞λ 2＞...＞λ n；
得到特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)选取主元个数原则：将最大的方差归纳到模型空间，将最小的方差留给噪声空间，计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000110
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000111
d)保留与k个最大特征值相对应的k个特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
其本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用该矩阵具体绘制置信度为95％的置信椭圆的步骤如下：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度为95％，对rd进行升序排序后，将rd第95％个位置所对应的值作为r∈R；
e)参照式(5)和式(6)计算4-a)中所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000112
特征向量矩阵为
Figure PCTCN2023070795-appb-000113
f)利用4-d)中的r以及4-e)中的D，按照公式(10)和公式(11)计算置信椭圆的长轴以及短轴：
Figure PCTCN2023070795-appb-000114
Figure PCTCN2023070795-appb-000115
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，置信椭圆的公式为：
Figure PCTCN2023070795-appb-000116
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图30所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域，根据CDU过程的性能指标，对历史数据设置不同的标签，将历史数据的得分矩阵的前两列进行投影，可以通过不同区域数据点所带的标签对置信椭圆进行划分，并进行框定，图31、图32和图33分别示出了对应以能耗为优化目标的模式、以轻收调优为目标的模式和以总拔为优化目标的模式的划分区间。图中框选出的区域分别对应能耗最低的50％数据点，轻收最高的50％数据点和总拔最高的50％数据点。
6.将CDU过程在正常工况下运行一段时间后，通过调整操作变量初馏塔顶压力的值使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用步骤(2)计算所得的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.按照公式(13)将Ym乘以3-d)所得负载矩阵P的前两列，得到scorey：
scorey＝Ym×[p 1,p 2]                  (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内，如图34所示。
9.设故障数据点为x∈R 1×n，根据公式(14)计算故障数据点的SPE贡献率：
Figure PCTCN2023070795-appb-000117
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000118
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡献率如图35所示。根据图35可以判定故障产生的原因主要是由第3个变量(即初馏塔顶压力)起主要作用，与实际操作情况相符(其中变量1为原料性质变量，在实际操作中难以直接调整)。
实施例6
本实施例中，基于大数据的加氢裂化过程模式识别和优化方法被应用于国 内某炼油企业加氢裂化(HCR)工艺过程中。图36给出了加氢裂化的工艺过程流程图。加氢裂化工艺过程由加氢精制、加氢裂化和分馏部分组成。在单段串联工艺中，原料油和循环氢分别达到加氢加压反应条件后混合进入加氢精制反应器，在加氢精制催化剂的作用下进行脱硫、脱氮、脱氧和部分脱芳烃反应。精制反应产物经注冷氢调整至裂化反应所需温度后进入加氢裂化反应器进行裂化反应，将重质馏分油转化为轻质馏分油。反应产物在冷高压分离器中进行气、油、水三相分离，经过循环氢脱硫塔脱硫后送至低压分离器，闪蒸出低分气体。低分油经换热进入脱丁烷塔，塔顶分离出含硫气体和轻烃，塔底油进入产品分离塔进行产品分离，分离出轻重石脑油、航煤、柴油和尾油等产品。
HCR模型包括32个输入变量和24个输出变量。通过对输入变量进行控制可以得到不同运行状态的HCR模型，其中对输出结果影响较大的输入变量有氢油比、反应温度、压力、进料流量和进料性质等，产物的性能指标主要通过观察输出变量中的干气、液化气、轻端、轻石脑油、重石脑油、航煤、柴油和尾油出装置的量来进行判断。表8给出了本实施例的HCR模型的操作变量。
表8：HCR模型的操作变量
变量 描述 单位 变量 描述 单位
1 罐区轻蜡油流量 tonne/h 17 原料油硫含量 ％
2 常减压蜡油总流量 tonne/h 18 原料油氮含量 ppmwt
3 催化柴油累计流量 tonne/h 19 原料油碱性氮 ppmwt
4 循环尾油质量流量 tonne/h 20 原料油初馏点 C
5 新氢流量 STD_m3/h 21 原料油10％回收温度 C
6 供吹扫循环氢流量 STD_m3/h 22 原料油50％回收温度 C
7 氢油比   23 原料油90％回收温度 C
8 R101一床层入口温度 C 24 原料油终馏点 C
9 R101二床层入口温度 C 25 轻柴油密度(20℃)  
10 R101三床层入口温度 C 26 轻柴油硫含量 ％
11 R102一床层入口温度 C 27 轻柴油初馏点 C
12 R102二床层入口温度 C 28 轻柴油10％回收温度 C
13 R102三床层入口温度 C 29 轻柴油50％回收温度 C
14 R102四床层入口温度 C 30 轻柴油90％回收温度 C
15 R101塔顶压力 MPa 31 轻柴油95％回收温度 C
16 原料油密度(20℃)   32 新氢氢气 ％
采用本发明的基于大数据的加氢裂化过程模式识别和优化方法对该加氢 裂化过程进行模式识别和优化，如图37所示，包括以下步骤：
1.采集一段时间HCR过程正常运行状态下的样本数据Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，z i＝[z 1i,z 2i,...,z mi] T代表第i个测量变量的m个样本；
2.对采集的数据进行预处理，得到均值为0、方差为1的标准数据集X＝[x 1,x 2,...,x n]∈R m×n，计算公式为：
Figure PCTCN2023070795-appb-000119
式中，μ为采集的数据的均值，其计算公式为：
Figure PCTCN2023070795-appb-000120
σ为采集的数据的标准差，其计算公式为：
Figure PCTCN2023070795-appb-000121
3.对标准数据集X运用主元分析法进行降维，按如下步骤进行：
a)根据公式(4)计算矩阵X的协方差矩阵：
Figure PCTCN2023070795-appb-000122
其中，X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，因此协方差矩阵C是一个n×n维的矩阵；
b)根据公式(5)和公式(6)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序：
det(C-Iλ)＝0                       (5)
Cp i＝λ ip i                         (6)
式(5)中，I为单位矩阵；
得到特征值矩阵：
Figure PCTCN2023070795-appb-000123
其中λ 1＞λ 2＞...＞λ n；
得到特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
c)选取主元个数原则：将最大的方差归纳到模型空间，将最小的方差留给噪声空间，计算公式为：
每个主元包括信息比：
Figure PCTCN2023070795-appb-000124
最大的k个主元包括的信息比：
Figure PCTCN2023070795-appb-000125
d)保留与k个最大特征值相对应的k个特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
t i＝Xp i,i＝1,2,...,k                    (9)
其本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
4.令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用该矩阵具体绘制置信度为95％的置信椭圆的步骤如下：
a)参照公式(4)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
b)计算xdat每一列的平均值得xm∈R 1×2，并将每一行的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
c)计算公式xd×s·×xd，并将所得矩阵每行进行求和，得到rd∈R m×1；
d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度为95％，对rd进行升序排序后，将rd第95％个位置所对应的值作为r∈R；
e)参照式(5)和式(6)计算4-a)中所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
Figure PCTCN2023070795-appb-000126
特征向量矩阵为
Figure PCTCN2023070795-appb-000127
f)利用4-d)中的r以及4-e)中的D，按照公式(10)和公式(11)计算置信椭圆的长轴以及短轴：
Figure PCTCN2023070795-appb-000128
Figure PCTCN2023070795-appb-000129
其中，a为长轴，b为短轴；
g)以4-b)中xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，置信椭圆的公式为：
Figure PCTCN2023070795-appb-000130
其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值，绘制的椭圆如图38所示。
5.不同生产模式下的数据在置信椭圆中的投影会分布在不同的区域，根据HCR过程的性能指标，对历史数据设置不同的标签，将历史数据的得分矩阵的前两列进行投影，可以通过不同区域数据点所带的标签对置信椭圆进行划分，并进行框定，图39、图40和图41分别示出了对应以总液收为优化目标的模式、以中间油收率为目标的模式和以价值增量为优化目标的模式的划分区间。其中，大框为按对应目标选取的前10％数据点，小框为前5％数据点。
6.将HCR过程在正常工况下运行一段时间后，通过调整操作变量氢油比的值使过程偏离正常状态，并收集数据得到Y∈R N×n，对所采集的数据运用步骤(2)计算所得的平均值和标准差进行归一化处理得到Ym∈R N×n。
7.将Ym乘以3-d)所得负载矩阵P的前两列得到scorey，计算公式为：
scorey＝Ym×[p 1,p 2]                  (13)
8.将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内，如图42所示。
9.设故障数据点为x∈R 1×n，根据公式(14)计算故障数据点的SPE贡献率：
Figure PCTCN2023070795-appb-000131
其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
Figure PCTCN2023070795-appb-000132
I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量，SPE贡献率如图43所示。根据图43可以判定故障产生的原因主要是由第1个变量(即氢油比)起主要作用，与实际操作情况相符。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (12)

  1. 一种基于大数据的炼油过程模式识别及优化方法，其特征在于，所述方法包括以下步骤：
    (1)将炼油过程中所采集到的历史数据组成建模用的训练样本集Z＝[z 1,z 2,...,z i,...,z n]∈R m×n，其中，m为样本集的样本个数，n为样本集的变量个数；
    (2)对训练样本集进行预处理，得到均值为0、方差为1的标准化数据X＝[x 1,x 2,...,x n]∈R m×n；
    (3)对X运用主元分析法将其从n维降到k维，并得到得分矩阵T∈R m×k和负载矩阵P∈R n×k；
    (4)利用得分矩阵T的前两列，绘制二维置信椭圆；
    (5)采集新的在线实时数据Y∈R N×n，采用步骤(2)对训练样本集进行预处理时所得到的样本平均值和样本方差对Y进行预处理得到标准化数据Ym∈R N×n；
    (6)将Ym乘以步骤(3)所得的负载矩阵P的前两列，得到Ym根据训练样本所得的前两列得分矩阵scorey∈R N×2；
    (7)以scorey的第一列作为x轴的数据，以scorey的第二列作为y轴的数据，将scorey映射到步骤(4)所绘制的置信椭圆上；若样本点映射到椭圆内，则说明此时炼油过程所处的工况为正常工况；若样本点映射到椭圆外，则说明此时炼油过程存在异常；
    (8)利用主成分分析过程中得到的标准差和均值求得置信椭圆中的点对应的原始变量，从而得到该点对应的效益值；在置信椭圆中区分效益值高低的分布；当工况运行模式位于置信椭圆中的某一点时，采用路径优化算法，获得当前位置向最优位置最快移动轨迹，将轨迹对应的点进行反变换，获得操作条件的改变方式，从而指导生产装置操作优化；其中，所述路径优化算法采用改进的A*算法，通过如式(17)所示的改进的评估函数f(x)来确定搜索方向和下一个到达的节点：
    Figure PCTCN2023070795-appb-100001
    其中，g(x)是代价函数，代表从起始节点到达当前节点x的实际所需代价；h(x)是启发式函数，代表从当前节点x到达目标节点估计所需的代价；profit(x)表示所选节点对应的经济效益。
  2. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(2)中，所述的预处理方法采用Z-score标准化方法，计算公式为：
    Figure PCTCN2023070795-appb-100002
    式中，Z＝[z 1,z 2,...,z m]为训练数据矩阵，X表示标准化后的数据矩阵，μ为训练数据的均值，σ为训练数据的标准差，μ和σ计算公式为：
    Figure PCTCN2023070795-appb-100003
    Figure PCTCN2023070795-appb-100004
  3. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(3)中，对经过预处理所得到的X运用主元分析法进行降维处理，具体按如下步骤进行：
    (3-a)计算矩阵X的协方差矩阵，协方差矩阵的计算公式为：
    Figure PCTCN2023070795-appb-100005
    X是m×n的矩阵，m是样本个数，n是变量个数，T表示转置，所以得到的协方差矩阵C是一个n×n维的矩阵；
    (3-b)计算协方差矩阵C的特征值λ i和特征向量p i，并按特征值降序排序，计算公式为：
    det(C-Iλ)＝0   (5)
    Cp i＝λ ip i     (6)
    式(5)中，I为单位矩阵，得：
    特征值矩阵：
    Figure PCTCN2023070795-appb-100006
    其中λ 1＞λ 2＞...＞λ n
    特征向量：V＝[p 1,p 2,p 3,...,p n]∈R n×n；
    (3-c)建立主元模型，计算公式为：
    计算每个主元包括的信息比：
    Figure PCTCN2023070795-appb-100007
    计算最大的k个主元包括的信息比：
    Figure PCTCN2023070795-appb-100008
    (3-d)保留对应k个最大特征值的特征向量，得到负载矩阵P＝[p 1,p 2,p 3,...,p k]∈R n×k，得分矩阵的计算公式为：
    t i＝Xp i,i＝1,2,...,k   (9)
    t i的本质是X向量在p i方向上的投影，得分矩阵T＝[t 1,t 2,...,t k]∈R m×k。
  4. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(4)中，令得分矩阵T的前两列为xdat＝[t 1,t 2]∈R m×2，利用矩阵xdat具体绘制置信椭圆的步骤如下：
    (4-a)计算xdat的协方差矩阵，并对协方差矩阵求逆，得s∈R 2×2；
    (4-b)计算xdat每一列的平均值得xm∈R 1×2，并将xdat每一列的数值减去对应的平均值进行中心化处理，得到xd∈R m×2；
    (4-c)计算公式xd×s·×xd，并对所得矩阵每行进行求和，得到rd∈R m×1；
    (4-d)绘制xdat的曲线图，根据xdat呈经验分布的特点，计算矩阵rd的百分位数，根据所要绘制的置信椭圆的置信度Confi，对rd进行升序排序后，得到rd第Confi个位置所对应的值r∈R，优选的Confi为95％；
    (4-e)计算步骤(4-a)所得矩阵s的特征值和对应的特征向量，得特征值矩阵为
    Figure PCTCN2023070795-appb-100009
    特征向量矩阵为
    Figure PCTCN2023070795-appb-100010
    (4-f)利用步骤(4-d)中的r以及步骤(4-e)中的D，得到置信椭圆的长轴以及短轴，计算公式为：
    Figure PCTCN2023070795-appb-100011
    Figure PCTCN2023070795-appb-100012
    其中，a为长轴，b为短轴；
    (4-g)以步骤(4-b)中的xm为椭圆的中心点，根据中心点以及椭圆的长轴和短轴绘制置信椭圆，椭圆的公式为：
    Figure PCTCN2023070795-appb-100013
    其中，xm1为xdat第一列的平均值，xm2为xdat第二列的平均值。
  5. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(6)中，scorey的计算公式为：
    scorey＝Ym×[p 1,p 2]    (13)
    式中，[p 1,p 2]∈R n×2是负载矩阵P的前两列。
  6. 如权利要求4所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(7)中，将scorey中的每一组数据代入公式(12)中与数值1进行比较，若大于1，则表示该组数据映射在椭圆外，若小于等于1，则表示该组数据映射在椭圆内。
  7. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(7)还包括将收集到的正常工况下的数据加入到历史数据重新建模。
  8. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(7)还包括对异常数据计算SPE贡献率来进行故障溯源；
    采用如下方法对异常数据计算SPE贡献率来进行故障溯源：
    假设异常数据为x∈R 1×n，其SPE贡献率的计算公式为：
    Figure PCTCN2023070795-appb-100014
    其中，contspe(i)表示第i个变量所对应的SPE贡献，ξ i表示n维单位矩阵的第i列，T表示转置，
    Figure PCTCN2023070795-appb-100015
    I为单位矩阵，P为训练样本所得负载矩阵，n为变量数，所求的SPE贡献率最大的变量即为发生故障的变量。
  9. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(4)还包括：将历史数据结合现场工艺知识对所绘制的置信椭圆进行插旗处理实现对不同性能等级数据的区域划分，所述区域划分优选包括：首先根据工艺知识和历史效益统计数据，根据时间序列找到工况对应的效益，根据效益高低进行不同等级的划分，再找到各个等级所对应的历史数据在椭圆中的分布并加以不同的标记；或者，步骤(4)还包括：根据炼油过程的性能指标，对不同性能等级的历史数据设置不同的标签，将历史数据的得分矩阵的前两列投影到置信椭圆上，通过不同区域数据点所带的标签对置信椭圆进行划分。
  10. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(7)还包括：对于映射在椭圆内部的样本，根据其所处位置结合步骤(4)中对椭圆进行的区域划分判断当前所处的性能等级，若其处于非最优的性能等级，则通过调整关键变量使其转换到所期望的更优的性能等级；
    具体地，将主元中所对应的系数最大的变量设定为优化时所要调整的关键变量，然后根据关键变量与在线数据在置信椭圆中的投影点位置变化之间的相关性，确定调整方向，实现对生产过程的优化；具体地，关键变量为负载矩阵P的第一列中绝对值最大的变量，通过调整该变量所对应的值，对炼油过程的生产模式进行调整，并通过实时监控观察数据是否往所期望的区域移动。
  11. 如权利要求1所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，步骤(8)中，利用公式(15)计算原始变量X ori：
    X ori＝(defen×PC T)·×std(X m)+mean(X m)   (15)
    其中，defen＝(x,y)∈R 1×2为置信椭圆中的点的横纵坐标，PC＝[p 1,p 2]∈R J×2为主成分分析建模时负载矩阵的前两列，std(X m)为主成分分析采用的基础样本的 方差，mean(X m)为主成分分析采用的基础样本的均值，m为样本个数，X ori为椭圆中defen＝(x,y)∈R 1×2位置对应的反变换后的原始变量；
    根据公式(16)计算投入产出效益指标：
    Figure PCTCN2023070795-appb-100016
    其中，
    Figure PCTCN2023070795-appb-100017
    为第i种产品的产量，
    Figure PCTCN2023070795-appb-100018
    为第i种产品的价格，
    Figure PCTCN2023070795-appb-100019
    为第j种原料的投料量，
    Figure PCTCN2023070795-appb-100020
    为第j种原料的价格，profit为效益指标。
  12. 如权利要求11所述的基于大数据的炼油过程模式识别及优化方法，其特征在于，所述炼油过程为催化重整过程、催化裂化过程、硫磺回收过程、渣油加氢过程、常减压过程或加氢裂化过程。
PCT/CN2023/070795 2022-01-10 2023-01-06 一种基于大数据的炼油过程模式识别及优化方法 WO2023131257A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202210022589.X 2022-01-10
CN202210022589.XA CN114239321A (zh) 2022-01-10 2022-01-10 一种基于大数据的炼油过程模式识别及优化方法

Publications (1)

Publication Number Publication Date
WO2023131257A1 true WO2023131257A1 (zh) 2023-07-13

Family

ID=80746302

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2023/070795 WO2023131257A1 (zh) 2022-01-10 2023-01-06 一种基于大数据的炼油过程模式识别及优化方法

Country Status (2)

Country Link
CN (1) CN114239321A (zh)
WO (1) WO2023131257A1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116664018A (zh) * 2023-07-28 2023-08-29 华能济南黄台发电有限公司 一种发电厂设备运行状态评价平台
CN117831659A (zh) * 2024-03-04 2024-04-05 山东钢铁股份有限公司 宽厚板质量在线检测的方法、装置、电子设备及存储介质

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239321A (zh) * 2022-01-10 2022-03-25 华东理工大学 一种基于大数据的炼油过程模式识别及优化方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109062196A (zh) * 2018-10-31 2018-12-21 东北大学 一种集成pca-ica的高炉过程监测及故障诊断方法
US10210862B1 (en) * 2016-03-21 2019-02-19 Amazon Technologies, Inc. Lattice decoding and result confirmation using recurrent neural networks
WO2020237729A1 (zh) * 2019-05-31 2020-12-03 东北大学 一种基于模式转移的虚拟机混合备用动态可靠性评估方法
CN113311803A (zh) * 2021-05-17 2021-08-27 北京航空航天大学 一种基于核主元分析的在轨航天器飞轮故障检测方法
WO2021253550A1 (zh) * 2020-06-16 2021-12-23 北京工业大学 一种基于双核t分布随机近邻嵌入的过程监测可视化方法
CN114239321A (zh) * 2022-01-10 2022-03-25 华东理工大学 一种基于大数据的炼油过程模式识别及优化方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9892238B2 (en) * 2013-06-07 2018-02-13 Scientific Design Company, Inc. System and method for monitoring a process
CN105700517B (zh) * 2016-03-09 2018-10-12 中国石油大学(北京) 一种炼化过程自适应数据驱动的早期故障监测方法及装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10210862B1 (en) * 2016-03-21 2019-02-19 Amazon Technologies, Inc. Lattice decoding and result confirmation using recurrent neural networks
CN109062196A (zh) * 2018-10-31 2018-12-21 东北大学 一种集成pca-ica的高炉过程监测及故障诊断方法
WO2020237729A1 (zh) * 2019-05-31 2020-12-03 东北大学 一种基于模式转移的虚拟机混合备用动态可靠性评估方法
WO2021253550A1 (zh) * 2020-06-16 2021-12-23 北京工业大学 一种基于双核t分布随机近邻嵌入的过程监测可视化方法
CN113311803A (zh) * 2021-05-17 2021-08-27 北京航空航天大学 一种基于核主元分析的在轨航天器飞轮故障检测方法
CN114239321A (zh) * 2022-01-10 2022-03-25 华东理工大学 一种基于大数据的炼油过程模式识别及优化方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
LI ZHI, YING YUHUI, YANG MINGLEI, ZHAO LIANG, ZHAO LING, DU WENLI: "Monitoring and path optimization of catalytic reformer in a refinery: Principal component analysis and A* algorithm application", EXPERT SYSTEMS WITH APPLICATIONS, ELSEVIER, AMSTERDAM, NL, vol. 209, 1 December 2022 (2022-12-01), AMSTERDAM, NL, pages 118358, XP093076603, ISSN: 0957-4174, DOI: 10.1016/j.eswa.2022.118358 *
LI, ZHI ET AL.: "Monitoring and Path Optimization of Catalytic Reformer in a Refinery:Principal Component Analysis and Big Data Application", 第32届中国过程控制会议(CPCC2021)论文集 (PROCEEDINGS OF CPCC2021), 31 August 2021 (2021-08-31), pages 1 - 2 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116664018A (zh) * 2023-07-28 2023-08-29 华能济南黄台发电有限公司 一种发电厂设备运行状态评价平台
CN116664018B (zh) * 2023-07-28 2023-10-31 华能济南黄台发电有限公司 一种发电厂设备运行状态评价平台
CN117831659A (zh) * 2024-03-04 2024-04-05 山东钢铁股份有限公司 宽厚板质量在线检测的方法、装置、电子设备及存储介质
CN117831659B (zh) * 2024-03-04 2024-05-03 山东钢铁股份有限公司 宽厚板质量在线检测的方法、装置、电子设备及存储介质

Also Published As

Publication number Publication date
CN114239321A (zh) 2022-03-25

Similar Documents

Publication Publication Date Title
WO2023131257A1 (zh) 一种基于大数据的炼油过程模式识别及优化方法
CN104765346B (zh) 一种炼油过程全流程建模方法
CN103927412B (zh) 基于高斯混合模型的即时学习脱丁烷塔软测量建模方法
CN109814513B (zh) 一种基于数据模型的催化裂化装置优化方法
CN106845796B (zh) 一种加氢裂化流程产品质量在线预测方法
CN112527788A (zh) 变压器监测数据异常值检测与清洗的方法及装置
CN110456756B (zh) 一种适用于连续生产过程全局运行状况在线评估的方法
CN102213949B (zh) 一种乙烯装置价值优化方法
US12049592B2 (en) Predictive control systems and methods with hydrocracker conversion optimization
CN106709654A (zh) 一种加氢裂化流程全局运行状况评估与质量追溯方法
CN104804761B (zh) 一种加氢裂化装置的收率实时预测方法
EP1810162A2 (en) Application of abnormal event detection technology to hydrocracking units
CN104965967A (zh) 一种常减压蒸馏装置的收率实时预测方法
CN106444672A (zh) 针对炼油和石化装置的分子水平的实时优化方法
CN104392098A (zh) 一种预测催化裂化汽油产率的方法
CN104463327A (zh) 一种预测催化裂化焦炭产率的方法
CN103524284A (zh) 一种乙烯裂解原料配置的预测和优化方法
CN110197222A (zh) 一种基于多分类支持向量机变压器故障诊断的方法
CN101169387B (zh) 一种在线确定常压塔顶石脑油质量指标的软测量方法
CN105740960B (zh) 一种工业加氢裂化反应条件的优化方法
CN111475957A (zh) 一种基于装置机理的炼油过程生产计划优化方法
CN112250538A (zh) 一种异丙苯精制工艺流程
CN117826606A (zh) 一种吸收稳定系统的变量优化方法及控制设备
CN111914888A (zh) 一种多工况识别与故障检测一体化的化工过程监测方法
CN111598305B (zh) 一种轻烃分离装置运行状态优化及预测方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23737119

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE