Thermal power generating unit working condition division method based on dynamic principal component time sequence clustering

Abstract

The invention relates to the technical field of thermal power generating unit monitoring and control, in particular to a thermal power generating unit working condition division method based on dynamic principal component time series clustering, which comprises the following steps of: S1, acquiring time series data of key operation parameters of a thermal power generating unit, removing abnormal values in the data, filling missing values and normalizing the data; s2, setting a sliding time window parameter; s3, dynamic principal component extraction and steady-state discrimination index calculation are carried out; s4, carrying out transition state re-division and boundary determination; and S5, working condition similarity calculation and clustering. Compared with the prior art, the method has the advantages of high division precision, high dynamic adaptability, high calculation efficiency and high universality.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of monitoring and control of thermal power generating units, in particular to a method for dividing the working conditions of thermal power generating units based on dynamic principal component time series clustering. BACKGROUND

[0002] Under the background of power marketization peak shaving, deep flexibility reconstruction and heat supply coupling operation, the load curve of the unit presents the characteristics of "high frequency, deep amplitude and multiple steps". In order to ensure the safety and economy of the unit operation, it is crucial to accurately grasp the operating characteristics of the unit under different working conditions. During the operation of the thermal power generating unit, the working condition is complex and changeable due to the influence of factors such as load fluctuation, coal quality change and change of external environmental conditions. According to the change of the operating parameters of the unit, the operating conditions are roughly divided into two categories: steady state and transient state. The steady state is the main state of the unit operation, and the operating parameters of the unit fluctuate little and the operation is stable. The transient state is the process of switching between different steady states, and the operating parameters change significantly, which is the high incidence stage of unit failure.

[0003] At present, the existing thermal power generating unit working condition division technology has the following shortcomings:

[0004] (1) Low division precision: the traditional method is mostly based on single or a small number of operating parameter thresholds for working condition judgment, ignoring the correlation between parameters, and it is difficult to accurately distinguish similar steady states and different types of transient states, and misjudgment and omission are prone to occur.

[0005] (2) Poor dynamic adaptability: unable to effectively respond to dynamic changes in the unit operation, for parameter slow drift, sudden disturbance and other conditions, the working condition division lags significantly, and it is difficult to reflect the real operating state of the unit in real time.

[0006] (3) Poor generalization ability of calculation method: the working condition division methods in different researches or application scenarios are quite different, resulting in that the division results are not comparable and universal, and it is difficult to be widely applied in the industry.

[0007] Therefore, there is an urgent need for a thermal power generating unit working condition division method that can balance precision, real-time performance and universality to solve the problems existing in the prior art and provide a high-confidence data foundation for unit energy-saving optimization, flexibility evaluation and life management. SUMMARY

[0008] The purpose of the present application is to provide a method for dividing the working conditions of thermal power generating units based on dynamic principal component time series clustering, to solve the problems of low division precision, poor dynamic adaptability and poor generalization ability of calculation method existing in the prior art working condition division technology of thermal power generating units.

[0009] To achieve the above object, the present application provides the following technical scheme: a dynamic principal component time series clustering thermal power unit working condition division method, comprising the following steps: S1: collecting time series data of key operation parameters of a thermal power unit, eliminating abnormal values in the data, filling in missing values and normalizing the data;

[0010] S2: setting a sliding time window parameter;

[0011] S3: dynamic principal component extraction and steady state discrimination index calculation;

[0012] S4: transition state re-division and boundary determination;

[0013] S5: working condition similarity calculation and clustering.

[0014] Preferably, the S1 comprises: collecting time series data of key operation parameters of a thermal power unit, including but not limited to generator power, steam turbine speed, steam temperature, steam pressure, condenser temperature, feed water flow, etc.; adopting 3σ criterion or box plot method, identifying and eliminating abnormal data points caused by sensor failure, data transmission error, etc.; if the data missing rate is <5%, linear interpolation method or adjacent time mean method is used for filling; if the missing rate is ≥5%, the data in this time period is discarded, ensuring that each dimension parameter is complete and the time axis is continuous without vacancy, avoiding affecting the subsequent working condition discrimination accuracy; the cleaned operation parameter data is normalized to eliminate the influence of different parameter order of magnitude differences.

[0015] Preferably, the S2 comprises: calculating the duration of different steady states, taking 1.2 to 1.5 times of the shortest steady state length as the sliding time window length, and setting the moving step length as 1 / 2 or 1 / 3 of the sliding window length.

[0016] Preferably, the S3 comprises: performing principal component analysis (PCA) on the normalized operation data in each sliding time window to extract dynamic principal components that can reflect the main variation information of the data; selecting a principal component center point in the window, calculating the Euclidean distance between all principal component points in a time window and the center point; setting a distance threshold θ and a point number threshold δ.

[0017] Preferably, the S4 comprises: dividing the potential transition state period determined in S3 into r small segments with a step length r, each small segment corresponding to a sub-time window; using the dynamic principal component extraction and steady state discrimination method in S3 to perform steady state discrimination on each sub-time window; determining the starting and ending boundary points of the transition state through the sub-window discrimination results.

[0018] Preferably, the S5 comprises: calculating the mean value of each operation parameter of each determined steady state to obtain a steady state mode vector ( The mean value of the i-th parameter) ; select the operating parameter with the largest mean value difference as the difference characteristic parameter, and define the difference degree of the two steady-state modes M1 and M2 as:

[0019]

[0020] The steady-state mode similarity calculation formula is:

[0021] S M (M1, M2) = 1 - D (M1, M2)

[0022] Set the similarity threshold ω M (determined by the historical working condition classification result, usually 0.8-0.9), judge by S M (M1, M2) > ω M .

[0023] Preferably, the S5 comprises: considering the fluctuation characteristics of the transition state operating parameters, constructing a delay matrix to reflect the dynamic changes of the parameters, setting the time series of a certain operating parameter in the transition state time period as x(t), (t = 1, 2,..., m), and the delay order as l (usually 5-10), then the delay matrix is:

[0024]

[0025] Perform principal component analysis on the delay matrix to extract principal component vectors U1, U2,..., Um. k The similarity of the two transition states is determined by calculating the cosine value of the principal component vectors, and the formula is as follows:

[0026]

[0027] Judge by setting the transition state similarity threshold ω PCA (usually 0.7-0.8).

[0028] Compared with the prior art, the beneficial effects of the present application are: by mining the dynamic correlation characteristics between operating parameters through dynamic principal component analysis, combining multi-index steady state discrimination and transition state boundary optimization, effectively distinguishing similar steady states and different types of transition states, reducing the misjudgment and omission rate, using sliding time window technology, real-time tracking of unit operating parameter changes, rapid response to parameter drift, sudden disturbance and other dynamic situations, condition recognition lag time ≤ 30 seconds, meeting the real-time monitoring demand, reducing the amount of calculation data and repeated calculation through data dimension reduction, reasonable setting of window and step parameters, real-time processing can be realized on ordinary industrial computers, single unit condition division calculation time ≤ 1 second / time, establishing a standardized condition division process, suitable for different types (such as extraction condensing type, back pressure type, low pressure cylinder zero output, etc.) thermal power units, and the parameter threshold and discrimination standard can be flexibly adjusted according to the characteristics of the unit, with wide industry promotion value. Compared with the prior art, the division accuracy is high, the dynamic adaptability is strong, the calculation efficiency is high, and the versatility is strong. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1 A novel method flowchart of the present application;

[0030] Figure 2 A time series sliding window diagram;

[0031] Figure 3 A unit condition division example diagram. DETAILED DESCRIPTION

[0032] The technical solutions in the embodiments of the present application will be described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0033] An embodiment provided by the present application: a dynamic principal component time series clustering thermal power unit condition division method, comprising the following steps: S1: collecting time series data of key operating parameters of a thermal power unit, eliminating outliers in the data, filling in missing values and normalizing the data;

[0034] S2: setting sliding time window parameters;

[0035] S3: dynamic principal component extraction and steady state discrimination index calculation;

[0036] S4: transition state redivision and boundary determination;

[0037] S5: condition similarity calculation and clustering.

[0038] Further, S1 comprises collecting time series data of key operating parameters of the thermal power generating unit, including but not limited to generator power, turbine speed, steam temperature, steam pressure, condenser temperature, feed water flow, etc. The data collection frequency is set according to the operating characteristics of the unit, usually 1-10 seconds / time, to ensure that the details of the unit operating condition changes can be captured. The 3σ criterion or the box plot method is used to identify and eliminate abnormal data points caused by sensor failure, data transmission error, etc. If the data missing rate is <5%, linear interpolation or adjacent time average method is used for filling. If the missing rate is ≥5%, the data in this time period is discarded to ensure that each dimension parameter is complete and the time axis is continuous without vacancy, avoiding affecting the subsequent operating condition discrimination accuracy. The cleaned operating parameter data is normalized to eliminate the influence of different parameter magnitude differences. The Min-Max normalization method is used to map the data to the [0, 1] interval, and the calculation formula is as follows:

[0039]

[0040] Further, S2 comprises: counting the duration of different steady-state operating conditions, taking 1.2 to 1.5 times the shortest steady-state length as the sliding time window length, and setting the moving step length as 1 / 2 or 1 / 3 of the sliding window length, to ensure that each window can completely cover the shortest steady-state period, avoiding misjudging the steady state as a transient state due to too short window. If the window length is T w , the step length is or which ensures that there is a certain overlap area between adjacent windows to avoid missing operating condition change information, and reduces unnecessary repeated calculation to improve processing speed.

[0041] Further, S3 comprises: performing principal component analysis (PCA) on the normalized operating data in each sliding time window to extract dynamic principal components that can reflect the main variation information of the data. The data matrix in the window is X m×n (m is the number of data points, n is the number of operating parameters), the covariance matrix is calculated , the eigenvalues and eigenvectors of the covariance matrix are solved, the first k eigenvectors with cumulative contribution rate ≥85% are selected as the principal component directions, and the corresponding principal component score matrix is Y m×k , realizing data dimension reduction and highlighting the dynamic correlation characteristics between parameters; selecting the center point of the principal components in the window, the Euclidean distance between all principal component points and the center point in a time window is calculated, the formula is as follows: D(y i ,y c )=|y i -y c |, wherein y i is the i-th principal component point in the window, y cFor a certain center point; set the distance threshold θ and the point threshold δ, the distance threshold θ is determined by the historical steady-state data statistics, usually take the maximum value of the distance from each point to the center point in the historical steady-state data, and the point threshold δ takes 80% to 90% of the number of data points in the window). If there are at least δ principal component points in the window that satisfy the distance D(y i ,y c )<θ from any one center point, it is preliminarily determined that the time period corresponding to the window is a potential steady state, otherwise, it is determined to be a potential transition state.

[0042] Further, S4 includes: uniformly segmenting the potential transition state time period determined in S3 with a step size r, dividing the transition state into r small segments, each of which corresponds to a sub-time window; using the dynamic principal component extraction and steady state discrimination method in S3, steady state discrimination is performed on each sub-time window, if a certain sub-window is determined to be a steady state, the time period corresponding to the sub-window is classified as an adjacent steady state; if it is still determined to be a transition state, the sub-window is determined to be a transition state core area; through the sub-window discrimination result, the starting and ending boundary points of the transition state are determined, the area within the boundary points is the final transition state, and the area outside the boundary points is merged into the adjacent steady state, avoiding the error of working condition division caused by the fuzzy boundary between the transition state and the steady state.

[0043] Further, S5 includes: for each determined steady state, calculating the mean value of each operating parameter to obtain a steady state mode vector ( is the mean value of the i-th parameter); selecting the operating parameter with the largest mean value difference as the difference characteristic parameter, and defining the difference degree of two steady state modes M1 and M2 as:

[0044]

[0045] The steady state mode similarity calculation formula is:

[0046] S M (M1,M2)=1-D(M1,M2)

[0047] Set the similarity threshold ω M (determined by the historical working condition classification result, usually 0.8 to 0.9), if the similarity S M (M1,M2) of two steady state modes is greater than ω M , it is determined to be a similar steady state and is classified into the same cluster category; otherwise, it is classified into different categories, and finally the steady state mode clustering result is obtained.

[0048] Further, S5 comprises: considering the fluctuation characteristics of the transition state operation parameters, constructing a delay matrix to reflect the dynamic changes of the parameters, setting the time series of an operation parameter in the transition state time period as x(t), (t=1, 2,..., m), and the delay order as l (usually 5-10), then the delay matrix is:

[0049]

[0050] Perform principal component analysis on the delay matrix, extract principal component vectors U1, U2,..., U k The similarity of two transition states is determined by calculating the cosine of the included angle of the principal component vectors thereof, and the formula is as follows:

[0051]

[0052] Set the transition state similarity threshold ω PCA (usually 0.7-0.8), if the similarity is greater than the threshold, it is determined as a similar transition state and is classified into the same cluster category; otherwise, it is classified into different categories.

[0053] Working principle: Taking the heat supply period working condition of a 300MW extraction condensing type thermal power unit as an example, the detailed implementation process of the present application is described. The operation data of the unit during the heat supply period for one month (a total of 720 hours) are collected, including the generator power, the main steam flow, the steam generator main feed water temperature, the steam turbine speed, the condenser A side cooling water inlet temperature, the data collection frequency is 5 seconds / time, and a total of 518400 data records are obtained. An industrial computer (CPU: Intel Core i7-10700K, memory: 32GB) is used, the software environment is Python 3.8, the Pandas library is used for data processing, the Scikit-learn library is used to realize principal component analysis and clustering algorithm, and the Matplotlib library is used for result visualization.

[0054] Step 1: Data preprocessing

[0055] Outlier rejection: The data of the five measuring points are respectively subjected to 3σ criterion for outlier detection, a total of 128 abnormal data points caused by instantaneous sensor failure are removed, and the proportion of outliers is about 0.025%.

[0056] Missing value filling: It is found that the missing rate of the measuring points is about 0.03% after checking the data, and the linear interpolation method is used for filling to ensure that all measuring points have no missing values and the time axis is continuous.

[0057] Data normalization: The Min-Max normalization method is used to map the data of the five measuring points to the [0, 1] interval, and the mean and standard deviation of the normalized data are controlled in the range of 0.5±0.2, so as to eliminate the influence of the order of magnitude difference.

[0058] Step 2: Sliding time window parameter setting

[0059] Window length (T w ) determination: Analyzing the historical operation data of the unit, the shortest steady-state duration is obtained as 300 seconds (corresponding to 60 data points), and the window length is taken as 1.2 times the shortest steady-state duration, i.e., T w = 72 data points (corresponding to 360 seconds).

[0060] Moving step (L) determination: The step is taken as 1 / 2 of the window length, i.e., L = 36 data points (corresponding to 180 seconds), balancing the calculation efficiency and data coverage.

[0061] Step 3: Steady-state discrimination

[0062] Dynamic principal component extraction: PCA analysis is performed on the 5-dimensional normalized data in each sliding window, and the eigenvalues and eigenvectors of the covariance matrix are calculated. The results show that the cumulative contribution rate of the first three principal components reaches 89.2%, so the first three principal components are extracted as dynamic principal components to realize data dimension reduction.

[0063] Steady-state discrimination: 5 principal component points are randomly selected as center points in each window, and the Euclidean distance of all points to the center points is calculated. Through historical steady-state data statistics, the distance threshold θ = 0.3 and the point threshold δ = 58

[0064] (80.6% of the number of data points in the window) are set. If at least 58 points in the window satisfy the distance < 0.3, it is determined as a potential steady state; otherwise, it is a potential transition state. After discrimination, 286 potential steady-state windows and 152 potential transition-state windows are obtained.

[0065] Step 4: Transition state re-division and boundary determination

[0066] Transition state segmentation: The number of transition state re-division r = 4, each potential transition state window (72 data points) is divided into 4 small segments, each segment has 18 data points, corresponding to 90 seconds.

[0067] Sub-window steady-state discrimination: PCA and steady-state discrimination are performed on each sub-window, and 32 sub-windows are found to be steady-state, which are classified as adjacent potential steady-state; the remaining sub-windows are still transition state, which are determined as transition state core area.

[0068] Boundary determination: According to the sub-window discrimination results, the transition state starting and ending boundaries are adjusted, and finally 120 transition state time periods and 218 steady-state time periods are determined, avoiding the division error caused by ambiguous boundaries.

[0069] Step 5: Working condition similarity calculation and clustering

[0070] Steady state mode clustering: Calculate the 5-point mean of each steady state to obtain the steady state mode vector. Select the generator power (point 1) as the difference characteristic parameter (its mean difference is the largest), calculate the steady state mode difference and similarity, and set the similarity threshold ω M = 0.85. After clustering, a total of 11 steady state mode categories are obtained, and the time period corresponding to each category accounts for 0.006-0.297, among which category 1 has the largest number (37) and the time period accounts for 0.297, which is the main steady state working condition of the heating period of the unit.

[0071] Transition state clustering: For each transition state, take the delay order l = 8, construct the delay matrix and perform PCA, extract the first three principal component vectors, and calculate the principal component angle cosine value as the initial similarity. Combined with the consistency of the front and rear steady states, the similarity threshold ω PCA = 0.75. After clustering, a total of 53 transition state categories are obtained, each of which corresponds to a specific transition process (such as load from 30% to 50%, load from 70% to 40%, etc.).

[0072] Comparing the division results with the unit operation log (such as load adjustment record, equipment start-stop record), the verification shows that the steady state working condition recognition accuracy is 96.3%, the transition state working condition recognition accuracy is 94.7%, and the misjudgment and omission mainly concentrate on the slight transition process with small parameter fluctuation, which meets the engineering application requirements as a whole.

[0073] It is obvious to those skilled in the art that the present application is not limited to the details of the above exemplary embodiments, but can be implemented in other specific forms without departing from the spirit or essential characteristics of the present application. Therefore, the embodiments should be regarded as exemplary and non-limiting, and the scope of the present application is defined by the appended claims rather than the above description, and therefore all changes falling within the meaning and scope of the equivalent elements of the claims are intended to be included in the present application. Any reference signs in the claims should not be regarded as limiting the claims involved.

Claims

1. A method for dividing the working condition of a thermal power unit by dynamic principal component time series clustering, characterized in that, The method comprises the following steps: S1: collecting time series data of key operation parameters of a thermal power unit, eliminating abnormal values in the data, filling in missing values and normalizing the data; S2: setting a sliding time window parameter; S3: extracting dynamic principal components and calculating a steady state discrimination index; S4: re-dividing a transition state and determining a boundary; S5: calculating a working condition similarity and clustering.

2. The dynamic principal component time series clustering method for dividing the working conditions of a thermal power generating unit according to claim 1, characterized in that, The S1 includes: collecting time series data of key operation parameters of the thermal power generating unit, including but not limited to generator power, steam turbine speed, steam temperature, steam pressure, condenser temperature, feed water flow, etc.; using 3σ criterion or box plot method, identifying and eliminating abnormal data points caused by sensor failure, data transmission error, etc.; if the data missing rate is <5%, using linear interpolation method or adjacent time mean method for filling; if the missing rate is ≥5%, discarding the data of this time period, ensuring that each dimension parameter is complete and the time axis is continuous without vacancy, avoiding affecting the subsequent working condition discrimination accuracy; normalizing the cleaned operation parameter data, eliminating the influence of different parameter magnitude differences 3. The method according to claim 1, wherein, The S2 comprises: counting the duration of different steady state working conditions, taking 1.2 to 1.5 times the shortest steady state length as the length of the sliding time window, and setting the moving step length as 1 / 2 or 1 / 3 of the length of the sliding window.

4. The method of claim 1, wherein, The S3 comprises: performing principal component analysis (PCA) on the normalized operation data in each sliding time window to extract dynamic principal components that can reflect the main variation information of the data; selecting a principal component center point in the window, calculating the Euclidean distance between all principal component points in the window and the center point; and setting a distance threshold θ and a point number threshold δ.

5. The dynamic principal component time series clustering method for thermal power unit operating condition classification of claim 1, wherein, The S4 comprises: dividing the potential transition state time period determined in S3 into r small segments by a uniform step length r, each small segment corresponding to a sub-time window; using the dynamic principal component extraction and steady state discrimination method in S3 to discriminate the steady state of each sub-time window; and determining the starting and ending boundary points of the transition state through the sub-window discrimination results.

6. The dynamic principal component time series clustering method for thermal power unit operating condition classification of claim 1, wherein, S5. For each determined steady state, calculate the mean value of each operating parameter to obtain a steady state mode vector For the mean value of the i-th parameter); select the operating parameter with the largest mean value difference as the difference characteristic parameter, and define the difference degree of the two steady state modes M1 and M2 as:​ The steady state mode similarity calculation formula is: S M (M1,M2) = 1 - D(M1,M2) Setting a similarity threshold ω M (M1, M2) > ω M (M1, M2) > ω M Judgment is made.

7. The method of claim 1, wherein the method further comprises: The S5 comprises: considering the fluctuation characteristics of the transition state operation parameters, constructing a delay matrix to reflect the dynamic changes of the parameters, taking the time series of an operation parameter in a transition state time period as x(t) (t=1, 2,..., m), and taking the delay order as l (usually 5-10), then the delay matrix is: The principal component analysis is performed on the delay matrix to extract principal component vectors U1, U2,..., U k The similarity of the two transition states is determined by calculating the cosine value of the included angle of their principal component vectors, and the formula is as follows: By setting a transition state similarity threshold ω PCA (typically taken 0.7-0.8) is made.

Images