An IES attack detection method based on non-invasive thermal load detection modeling

Abstract

This invention designs an IES attack detection method based on non-invasive heat load detection modeling, belonging to the field of integrated energy technology. Addressing the thermal inertia of the heating network in an integrated electric-thermal energy system, a load redistribution attack model is constructed. Based on this, a non-invasive detection model is used to establish dual models on the load side: a radiator model and a building storage heat dissipation model, to detect amplified attack deviations in the room temperature state matrix. State prediction is performed using a sliding window method for state matrix similarity matching, replacing the attacked and contaminated state matrix with the predicted state matrix. The MILP method is used to enable rapid load recovery after an attack on the integrated electric-thermal energy system. This invention's detection method is fast and effective, enabling online detection, and has a certain theoretical basis and practical engineering significance.

Description

Technical Field

[0001] This invention relates to the field of integrated energy technology, and specifically to an IES attack detection method based on non-intrusive heat load detection modeling. Background Technology

[0002] Integrated Energy Systems (IES) improve energy efficiency and reduce environmental pollution through the complementary and tiered utilization of multiple energy sources such as electricity, heat, and gas. Advanced Information and Communication Technologies (ICT) play a crucial role in ensuring the secure, efficient, clean, and flexible operation of industrial information systems and promoting deep coupling between their network and physical systems. However, the coupling between various energy systems and their supporting information systems increases system complexity, introduces more network vulnerabilities, and poses greater cybersecurity challenges to the secure and efficient operation of IES.

[0003] In recent years, energy system failures caused by cyberattacks have occurred frequently. Therefore, it is urgent to study the impact of cyber threats on IES operations, especially the cascading effects of cyberattacks from one system to another.

[0004] Existing anomaly detection methods can be categorized into bias-based detection and feature-based detection methods based on their identification criteria. Bias-based methods typically select one or more variables strongly correlated with attacks based on the target system. When the values ​​of these variables deviate significantly from their normal range during operation, an attack is considered to have occurred. Feature-based detection methods extract features from the system during normal operation and under attack through physical mechanism analysis or artificial intelligence methods, and determine whether an attack has occurred by comparing these features during detection. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention proposes an IES attack detection method based on non-invasive thermal load detection modeling.

[0006] An IES attack detection method based on non-invasive thermal load detection modeling specifically includes the following steps:

[0007] Step 1: Construct a load redistribution attack model to address the thermal inertia of the integrated electric and thermal energy system's heating network;

[0008] The establishment of the load redistribution attack model is specifically as follows:

[0009] Indoor temperature setpoint attack, also known as an indoor temperature setpoint attack, is a type of HLR attack that targets measured indoor temperatures. It works by tampering with the indoor temperature... indoor setpoint temperature The mismatch between them induces the heating system to continuously supply inappropriate electricity to the heat load:

[0010]

[0011] Where e is a constant, Γ represents indoor temperature, t represents time, and Γ represents [missing information]. a The attack time range λ is the indoor temperature setpoint. a1 Here, a1 represents the attack parameters, and t0 represents the attack start time.

[0012] Step 2: Based on the load redistribution attack model built in Step 1, a non-intrusive detection model is adopted to establish a dual model on the load side: a radiator model and a building storage heat dissipation model, and to detect the amplification attack deviation of the room temperature state matrix.

[0013] Step 2.1: Establish a load-side radiator model;

[0014] The heat exchanger is modeled and analyzed, and the heat fluid energy equation is as follows:

[0015]

[0016]

[0017] Among them, C h For the heat capacity of the heat fluid, T h,i T represents the inlet water temperature. h,e A represents the outlet water temperature. h h is the contact area between the fluid and the radiator. h T is the heat transfer coefficient. h,m T is the average temperature of the thermal fluid. w The effective wall temperature is given by t, where t is time.

[0018] Wall energy equation:

[0019]

[0020] Among them, U f Let A be the total heat dissipation coefficient of the fins. f,b T represents the effective area of ​​the heatsink. f For the effective fin area, C W For the heat capacity of the pipe wall;

[0021] The area of ​​the heat sink is much larger than the area of ​​the pipe wall. Considering only the area of ​​the heat sink, the energy equation for the heat sink is:

[0022]

[0023] Among them, h c Let A be the heat transfer coefficient of the cold fluid. c T represents the effective area of ​​the cold fluid. c,m C is the average temperature of the cold fluid.f For fin heat capacity;

[0024] Assuming the temperatures of the hot and cold fluids within the heat exchanger are linearly distributed, the energy equation for the cold fluid is:

[0025]

[0026] Among them, C c For the heat capacity of the cold fluid, T c,i The initial temperature of the cold fluid;

[0027] Step 2.2: Establish a building storage heat dissipation model;

[0028] A modeling analysis of building heat transfer was performed, and the heat lost from the building was:

[0029] Q WH =S1U1(T c,i -T out (1+x1)(1+x2)

[0030] Q CH =S2U2(T c,i -T out (1+x1)(1+x2)

[0031] Q J =C1m1(T c,i -T out )

[0032] Among them, Q WH S1 is the heat dissipation power of the enclosure, S1 is the enclosure area, U1 is the heat transfer coefficient of the enclosure, and T is the heat dissipation power of the enclosure. out The outdoor temperature is represented by x1, the building orientation correction factor is x2, and the floor level correction factor is Q. CH S2 is the heat dissipation power of the window, S2 is the window area, U2 is the window heat transfer coefficient, and Q is the heat dissipation power of the window. J The heat dissipation due to the exchange of air between indoors and outdoors, C1 is the heat capacity of air, and m1 is the mass of air exchanged between indoors and outdoors per unit time;

[0033] The heat stored in the walls and indoor air is:

[0034]

[0035] Q T =C1m3T c,i

[0036] Among them, Q q The wall stores heat, C2 is the wall's heat capacity, m2 is the wall's mass, and Q is the heat storage capacity. T Air stores heat, C1 is the heat capacity of air, and m3 is the mass of air;

[0037] The amount of heat entering the room is:

[0038] Q R =Cm4(T h,i -T h,e )

[0039] Among them, Q R The heat transferred into the room by the radiator is C, where C is the specific heat capacity of water and m4 is the water mass flow rate.

[0040] The difference between the heat input into the room and the heat loss from the room:

[0041] ΔQ=Q R -Q WH -Q CH -Q J -Δ(Q q +Q T )

[0042] The state matrix data at each moment includes the difference between the heat input into the room and the heat lost from the room, the inlet temperature, the outlet temperature, and the outdoor temperature;

[0043] Step 3: Perform state prediction by matching similarity days of the state matrix using a sliding window, and replace the attacked and corrupted state matrix with the predicted state matrix;

[0044] Step 3.1: Under normal operation of the integrated energy system, status data is generated every moment and stored as a historical database. The most similarity matrix is ​​searched in the historical database through multiple matching.

[0045] Using a multi-state prediction method, the constructed historical state variables are arranged in chronological order to form a q-row, p-column multidimensional time series matrix T. q×p Where q represents the number of node state variables, and p represents the number of state variable acquisition points; T q×p The set of L continuous state quantity acquisition points in the middle is used as the state model of the power grid, which represents the changing trend of the state quantities of each node of the power grid over time; where L represents the length of the time window, and the state quantity x is a q-row, 1-column matrix; a multidimensional time series refers to a set of observation values ​​obtained at various times on the same time axis for a set of indicators of the observed object, and arranged in a time series.

[0046] In the historical database, assuming at time k (X) V ) q×L Let X be the current state matrix, and let (X1) be the matrix adjacent to the current state matrix. q×L ;with (X1) q×L Based on the baseline, with a sliding interval length of w and a time window length of L, (X1) q×LSlide backwards along the time axis to obtain the state matrix X for the j-th time window. j =[x k-2L-w(j-1)+1 ,x k-2L-w(j-1)+2 ,…,x k-L-w(j-1) In the formula Indicates less than or equal to The largest integer that can be obtained;

[0047] By using a sliding time window model, the discrete time series is divided into a set of multiple q-row, L-column historical state matrices, i.e. Generally speaking, the time window length L is greater than the sliding interval length w. Under the sliding time window model, multiple matching state prediction refers to finding the state matrix with the highest similarity to the current state matrix from the historical state matrix X. Its following state matrix has the same trend as the following state matrix of the current state matrix, and is used as the state prediction result.

[0048] Step 3.1.1: Detect the current state variables based on similarity analysis;

[0049] In terms of similarity analysis, based on the actual situation of power grid state variables and the Pearson correlation coefficient, a similarity metric suitable for heating networks is proposed; any two state matrices A′=[a′] are selected from the historical database. q×L And B′=[b′] q×L The similarity metric is defined as follows:

[0050]

[0051] In the formula: A′ and B′ represent the state matrix with q rows and L columns in the two history matrices;

[0052] Step 3.1.2: For the state matrix set obtained from similarity analysis, the density space clustering algorithm DBSCAN is used to obtain the optimal matching result of the state matrix. The DBSCAN method essentially clusters samples into clusters where points within the same cluster are density-connected, and points between different clusters are not connected. Mahalanobis distance is selected as the clustering metric. Based on the results obtained from DBSCAN, K clusters are obtained, with the cluster center C = {C1, C2, ..., C...}. K The condition that must be satisfied is shown in the following formula, which is to ensure that the sum of the Mahalanobis distances between the cluster center and the other state matrices in the cluster is minimized;

[0053]

[0054] In the formula: R K U is the number of state matrix samples in the Kth cluster; t M(C) is the state matrix; K) represents the Mahalanobis distance between the center points of the i-th cluster;

[0055] Exact matching in cluster analysis is achieved by comparing the current state matrix (X). V ) q×L The final matching result is achieved by selecting the state matrix with the smallest difference from the state matrices of each cluster center.

[0056] Step 3.1.3: To measure the overall difference in state variables between two state matrices, a difference index is used to compare and measure the state matrices; the difference is judged by combining the eigenvectors and the distance between the characteristic trends of the matrices.

[0057] For a state matrix X q×L Let the eigenvector F be (f x ,f y Let represent the maximum and minimum differences between the vector at a certain moment in the state matrix and the average value vector of the state matrix; let C be the difference between the vector at a certain moment in the state matrix and the average value vector of the state matrix. 1×L ,but:

[0058]

[0059] Where X i (t j ) indicates at time point t j The value corresponding to the i-th state variable; the feature vector F is represented as:

[0060] F = (f x ,f y )=(max(C) 1×L ,min(C) 1×L )

[0061] For a state matrix C q×L Its characteristic trend distance is its L2 norm D, which is expressed as:

[0062]

[0063] For any state matrix, it is represented by a pair G = (F, D) consisting of the distance between the eigenvector and the eigentrend.

[0064] For any two state matrices X a and X b Difference X ab Represented as:

[0065]

[0066] The difference degree is used to measure the overall difference between two state matrices; the greater the difference degree, the lower the similarity of corresponding elements in the state matrices, and vice versa.

[0067] Step 3.2: Based on the most similar matrix matched in Step 3.1, compare it with the current state matrix to detect whether the current time is under attack;

[0068] An attack is detected when the following equation is satisfied;

[0069] K≤max((X v ) 1,end -(X j ) 1,end ,,(X v ) i,end -(X j ) i,end )

[0070] K is the detection threshold, X v Let X be the matrix at the current time. j The matrix with the most similarity to the historical matrix;

[0071] Step 3.3: Replace the attacked and contaminated state matrix with the predicted state matrix;

[0072] For the detected attack matrix, it is replaced with historical data and stored in the historical database to avoid misscheduling by the scheduling center and to provide data support for subsequent state prediction.

[0073] Step 4: Rapid load recovery of the integrated electrothermal energy system after an attack based on the MILP method;

[0074] Step 4.1: Overall Optimization Objective: The overall optimization objective is to achieve the fastest recovery of the heat load state while keeping costs as low as possible;

[0075] The overall optimal objective function includes the following costs: 1) gas cost g, 2) grid-related costs e, 3) environmental emission costs, and 4) user dissatisfaction. Adding a load temperature difference integral penalty term to the objective function effectively improves the load recovery speed.

[0076]

[0077] in, These represent the hourly electricity price, the hourly gas price, and the prices related to emissions from the power grid and emissions from the gas network, respectively. Let represent the power purchased from the main grid during the t-th time interval, and let represent the power associated with the cogeneration unit and boiler unit during the t-th time interval. t1 and t2 represent the times from the start of the power attack to the system's return to stability. setIndoor temperature setting;

[0078] In this framework, electricity demand is met through electricity purchased from the main grid, electricity generated by CHP units, electricity generated by WT units, and electricity emitted by ES equipment; in addition, demand response (DR) services are considered to provide load reduction when necessary; reducing some electricity demand within a specified time period and increasing electricity demand at other time intervals; the following formula represents the electricity demand of the developed EH:

[0079]

[0080] In the formula P t Electrical Let be the electricity demand in time interval t; Let be the power generation of the wind turbine generator in the t-th time period; and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively. and These represent the power of the energy storage device during the charging and discharging processes in the t-th time interval; and These are the electrical efficiency of the transformer, the gas-electric efficiency of the CHP unit, and the electrical efficiency of the wind turbine generator set, respectively.

[0081] Heat load requirements:

[0082]

[0083] In the formula, P t Thermal The heat demand for the t-th time interval; and The respective gas thermal efficiency of the boiler and the gas thermal efficiency of CHP; finally, This refers to the gas purchased from the gas network at regular intervals. and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively.

[0084] Step 4.2: Quickly restore constraints;

[0085] When the system performs optimized scheduling, there will be some upper and lower capacity limits for basic equipment; the power transmission constraints of transmission lines are:

[0086]

[0087] p l,line To transmit power to transmission lines in real time, This is the lower limit of the transmission power of the transmission line. This is the upper limit of the transmission power of the power transmission line;

[0088] The relevant parameters on the heating network side indicate that the heat loss from the pipeline is as follows:

[0089]

[0090] T end T represents the water outlet temperature of the pipe. start T represents the inlet temperature of the pipe. a Let λ be the ambient temperature, λ be the heat transfer coefficient of the pipe, L be the pipe length, and C be the temperature. p The specific heat capacity of water, For mass flow rate;

[0091]

[0092] m l,line This represents the real-time quality flow rate of the pipeline. This represents the upper limit of the pipeline's mass flow rate.

[0093] From the law of conservation of nodal energy:

[0094]

[0095] These are the mixing temperatures of the water supply pipe nodes and the return water pipe nodes, respectively. These are the mass flow rates of the water supply pipeline node and the return water pipeline node, respectively. These are the temperatures of the supply and return water pipes, respectively; according to the law of conservation of flow at the nodes, we can obtain:

[0096]

[0097] These are the mass flow rates of the heat exchanger and the heat source at different times;

[0098] CHP unit thermo-electric operational domain constraints and ramp-up constraints:

[0099]

[0100] Let the output power of the i-th unit at time t be... This represents the lower limit of the output electrical power of the i-th unit at time t, which is related to thermal power. This represents the upper limit of the output electrical power of the i-th unit at time t, which is related to thermal power.

[0101]

[0102] Let the output thermal power of the i-th unit at time t be... This represents the lower limit of the output thermal power of the i-th unit at time t, which is related to the electrical power. This represents the upper limit of the output thermal power of the i-th unit at time t, which is related to the electrical power.

[0103]

[0104] For the output power of the i-th unit at time t-1, The range of power change for the i-th unit within a time interval.

[0105] The beneficial effects of this invention are:

[0106] This invention designs an online defense process for integrated electric and thermal energy network security. It comprises two parts: establishing a heat load model and mining state variables, as well as processing malicious data. State variable mining includes historical state variable acquisition and malicious data processing, which includes malicious data detection, removal, and correction. When acquiring state variables, the privacy of heat load user data is considered, and a non-intrusive detection method is adopted. The established heat load model effectively protects user privacy while ensuring efficient acquisition of state variables. The heat load model involved in the non-intrusive detection model proposed in this invention can effectively amplify attack bias, thereby improving the sensitivity of attack detection. The detection method of this invention is fast and effective, enabling online detection, and has a certain theoretical basis and practical engineering significance. Attached Figure Description

[0107] Figure 1 The experimental simulation topology of the attack detection method in this embodiment of the invention;

[0108] Figure 2 An attack comparison diagram of an IES attack detection method based on non-invasive thermal load detection modeling according to an embodiment of the present invention;

[0109] Figure 3 Comparison chart of attack detection methods for IES attacks based on non-invasive thermal load detection modeling in this invention; Detailed Implementation

[0110] The present invention will be further described below with reference to the accompanying drawings and embodiments;

[0111] An IES attack detection method based on non-invasive thermal load detection modeling specifically includes the following steps:

[0112] Step 1: Construct a load redistribution attack model to address the thermal inertia of the integrated electric and thermal energy system's heating network;

[0113] The establishment of the load redistribution attack model is specifically as follows:

[0114] Indoor temperature setpoint attack, also known as an indoor temperature setpoint attack, is a type of HLR attack that targets measured indoor temperatures. It works by tampering with the indoor temperature... indoor setpoint temperature The mismatch between them induces the heating system to continuously supply inappropriate electricity to the heat load:

[0115]

[0116] Where e is a constant, Γ represents indoor temperature, t represents time, and Γ represents [missing information]. a The attack time range λ is the indoor temperature setpoint. a1 Here, a1 represents the attack parameters, and t0 represents the attack start time.

[0117] Step 2: Based on the load redistribution attack model built in Step 1, a non-intrusive detection model is adopted to establish a dual model on the load side: a radiator model and a building storage heat dissipation model, and to detect the amplification attack deviation of the room temperature state matrix.

[0118] Step 2.1: Establish a load-side radiator model;

[0119] The heat exchanger is modeled and analyzed, and the heat fluid energy equation is as follows:

[0120]

[0121]

[0122] Among them, C h For the heat capacity of the heat fluid, T h,i T represents the inlet water temperature. h,e A represents the outlet water temperature. h h is the contact area between the fluid and the radiator. h T is the heat transfer coefficient. h,m T is the average temperature of the thermal fluid. w The effective wall temperature is given by t, where t is time.

[0123] Wall energy equation:

[0124]

[0125] Among them, U f Let A be the total heat dissipation coefficient of the fins. f,b T represents the effective area of ​​the heatsink. f For the effective fin area, C W For the heat capacity of the pipe wall;

[0126] The area of ​​the heat sink is much larger than the area of ​​the pipe wall. Considering only the area of ​​the heat sink, the energy equation for the heat sink is:

[0127]

[0128] Among them, h c Let A be the heat transfer coefficient of the cold fluid. c T represents the effective area of ​​the cold fluid. c,m C is the average temperature of the cold fluid. f For fin heat capacity;

[0129] Assuming the temperatures of the hot and cold fluids within the heat exchanger are linearly distributed, the energy equation for the cold fluid is:

[0130]

[0131] Among them, C c For the heat capacity of the cold fluid, T c,i The initial temperature of the cold fluid;

[0132] Step 2.2: Establish a building storage heat dissipation model;

[0133] A modeling analysis of building heat transfer was performed, and the heat lost from the building was:

[0134] Q WH =S1U1(T c,i -T out (1+x1)(1+x2)

[0135] Q CH =S2U2(T c,i -T out (1+x1)(1+x2)

[0136] Q J =C1m1(T c,i -T out )

[0137] Among them, Q WH S1 is the heat dissipation power of the enclosure, S1 is the enclosure area, U1 is the heat transfer coefficient of the enclosure, and T is the heat dissipation power of the enclosure. out The outdoor temperature is represented by x1, the building orientation correction factor is x2, and the floor level correction factor is Q. CH S2 is the heat dissipation power of the window, S2 is the window area, U2 is the window heat transfer coefficient, and Q is the heat dissipation power of the window. J The heat dissipation due to the exchange of air between indoors and outdoors, C1 is the heat capacity of air, and m1 is the mass of air exchanged between indoors and outdoors per unit time;

[0138] The heat stored in the walls and indoor air is:

[0139]

[0140] Q T =C1m3T c,i

[0141] Among them, Q qThe wall stores heat, C2 is the wall's heat capacity, m2 is the wall's mass, and Q is the heat storage capacity. T Air stores heat, C1 is the heat capacity of air, and m3 is the mass of air;

[0142] The amount of heat entering the room is:

[0143] Q R =Cm4(T h,i -T h,e )

[0144] Among them, Q R The heat transferred into the room by the radiator is C, where C is the specific heat capacity of water and m4 is the water mass flow rate.

[0145] The difference between the heat input into the room and the heat loss from the room:

[0146] ΔQ=Q R -Q WH -Q CH -Q J -Δ(Q q +Q T )

[0147] The state matrix data at each moment includes the difference between the heat input into the room and the heat lost from the room, the inlet temperature, the outlet temperature, and the outdoor temperature;

[0148] Step 3: Perform state prediction by matching similarity days of the state matrix using a sliding window, and replace the attacked and corrupted state matrix with the predicted state matrix;

[0149] Step 3.1: Under normal operation of the integrated energy system, status data is generated every moment and stored as a historical database. The most similarity matrix is ​​searched in the historical database through multiple matching.

[0150] Using a multi-state prediction method, the constructed historical state variables are arranged in chronological order to form a q-row, p-column multidimensional time series matrix T. q×p Where q represents the number of node state variables, and p represents the number of state variable acquisition points; T q×p The set of L continuous state quantity acquisition points in the middle is used as the state model of the power grid, which represents the changing trend of the state quantities of each node of the power grid over time; where L represents the length of the time window, and the state quantity x is a q-row, 1-column matrix; a multidimensional time series refers to a set of observation values ​​obtained at various times on the same time axis for a set of indicators of the observed object, and arranged in a time series.

[0151] In the historical database, assuming at time k (X) V ) q×LLet X be the current state matrix, and let (X1) be the matrix adjacent to the current state matrix. q×L ;with (X1) q×L Based on the baseline, with a sliding interval length of w and a time window length of L, (X1) q×L Slide backwards along the time axis to obtain the state matrix X for the j-th time window. j =[x k-2L-w(j-1)+1 ,x k-2L-w(j-1)+2 ,…,x k-L-w(j-1) In the formula Indicates less than or equal to The largest integer that can be obtained;

[0152] By using a sliding time window model, the discrete time series is divided into a set of multiple q-row, L-column historical state matrices, i.e. Generally speaking, the time window length L is greater than the sliding interval length w. Under the sliding time window model, multiple matching state prediction refers to finding the state matrix with the highest similarity to the current state matrix from the historical state matrix X. Its following state matrix has the same trend as the following state matrix of the current state matrix, and is used as the state prediction result.

[0153] Step 3.1.1: Detect the current state variables based on similarity analysis;

[0154] In terms of similarity analysis, based on the actual situation of power grid state variables and the Pearson correlation coefficient, a similarity metric suitable for heating networks is proposed; any two state matrices A′=[a′] are selected from the historical database. q×L And B′=[b′] q×L The similarity metric is defined as follows:

[0155]

[0156] In the formula: A′ and B′ represent the state matrix with q rows and L columns in the two history matrices;

[0157] Step 3.1.2: For the state matrix set obtained from similarity analysis, the density spatial clustering of application with noise (DBSCAN) algorithm is used to obtain the optimal matching result of the state matrix. DBSCAN essentially clusters samples into clusters where points within the same cluster are density-connected, and points between different clusters are not connected. Mahalanobis distance is selected as the clustering metric. Based on the results obtained from DBSCAN, K clusters are obtained, with the cluster center C = {C1, C2, ..., C...}. KThe condition that must be satisfied is shown in the following formula, which is to ensure that the sum of the Mahalanobis distances between the cluster center and the other state matrices in the cluster is minimized;

[0158]

[0159] In the formula: R K U is the number of state matrix samples in the Kth cluster; t M(C) is the state matrix; K ) represents the Mahalanobis distance between the center points of the i-th cluster;

[0160] Exact matching in cluster analysis is achieved by comparing the current state matrix (X). V ) q×L The final matching result is achieved by selecting the state matrix with the smallest difference from the state matrices of each cluster center.

[0161] Step 3.1.3: To measure the overall difference in state variables between two state matrices, a difference index is used to compare and measure the state matrices; the difference is judged by combining the eigenvectors and the distance between the characteristic trends of the matrices.

[0162] For a state matrix X q×L Let the eigenvector F be (f x ,f y Let represent the maximum and minimum differences between the vector at a certain moment in the state matrix and the average value vector of the state matrix; let C be the difference between the vector at a certain moment in the state matrix and the average value vector of the state matrix. 1×L ,but:

[0163]

[0164] Where X i (t j ) indicates at time point t j The value corresponding to the i-th state variable; the feature vector F is represented as:

[0165] F = (f x ,f y )=(max(C) 1×L ,min(C) 1×L )

[0166] For a state matrix C q×L Its characteristic trend distance is its L2 norm D, which is expressed as:

[0167]

[0168] For any state matrix, it is represented by a pair G = (F, D) consisting of the distance between the eigenvector and the eigentrend.

[0169] For any two state matrices X a and X b Difference X ab Represented as:

[0170]

[0171] The difference degree is used to measure the overall difference between two state matrices; the greater the difference degree, the lower the similarity of corresponding elements in the state matrices, and vice versa.

[0172] Step 3.2: Based on the most similar matrix matched in Step 3.1, compare it with the current state matrix to detect whether the current time is under attack;

[0173] An attack is detected when the following equation is satisfied;

[0174] K≤max((X v ) 1,end -(X j ) 1,end ,,(X v ) i,end -(X j ) i,end )

[0175] K is the detection threshold, X v Let X be the matrix at the current time. j The matrix with the most similarity to the historical matrix;

[0176] Step 3.3: Replace the attacked and contaminated state matrix with the predicted state matrix;

[0177] For the detected attack matrix, it is replaced with historical data and stored in the historical database to avoid misscheduling by the scheduling center and to provide data support for subsequent state prediction.

[0178] Step 4: Rapid load recovery of the integrated electrothermal energy system after an attack based on the MILP method;

[0179] Step 4.1: Overall Optimization Objective: The overall optimization objective is to achieve the fastest recovery of the heat load state while keeping costs as low as possible;

[0180] The overall optimal objective function includes the following costs: 1) gas cost g, 2) grid-related costs e, 3) environmental emission costs, and 4) user dissatisfaction. Adding a load temperature difference integral penalty term to the objective function effectively improves the load recovery speed.

[0181]

[0182] in, These represent the hourly electricity price, the hourly gas price, and the prices related to emissions from the power grid and emissions from the gas network, respectively. Let represent the power purchased from the main grid during the t-th time interval, and let represent the power associated with the cogeneration unit and boiler unit during the t-th time interval. t1 and t2 represent the times from the start of the power attack to the system's return to stability. set Indoor temperature setting;

[0183] In this framework, electricity demand is met through electricity purchased from the main grid, electricity generated by CHP units, electricity generated by WT units, and electricity emitted by ES equipment; in addition, demand response (DR) services are considered to provide load reduction when necessary; reducing some electricity demand within a specified time period and increasing electricity demand at other time intervals; the following formula represents the electricity demand of the developed EH:

[0184]

[0185] In the formula P t Electrical Let be the electricity demand in time interval t; Let be the power generation of the wind turbine generator in the t-th time period; and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively. and These represent the power of the energy storage device during the charging and discharging processes in the t-th time interval; and These are the electrical efficiency of the transformer, the gas-electric efficiency of the CHP unit, and the electrical efficiency of the wind turbine generator set, respectively.

[0186] Heat load requirements:

[0187]

[0188] In the formula, P t Thermal The heat demand for the t-th time interval; and The respective gas thermal efficiency of the boiler and the gas thermal efficiency of CHP; finally, This refers to the gas purchased from the gas network at regular intervals. and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively.

[0189] Step 4.2: Quickly restore constraints;

[0190] When the system performs optimized scheduling, there will be some upper and lower capacity limits for basic equipment; the power transmission constraints of transmission lines are:

[0191]

[0192] p l,line To transmit power to transmission lines in real time, This is the lower limit of the transmission power of the transmission line. This is the upper limit of the transmission power of the power transmission line;

[0193] The relevant parameters on the heating network side indicate that the heat loss from the pipeline is as follows:

[0194]

[0195] T end T represents the water outlet temperature of the pipe. start T represents the inlet temperature of the pipe. a Let λ be the ambient temperature, λ be the heat transfer coefficient of the pipe, L be the pipe length, and C be the temperature. p The specific heat capacity of water, For mass flow rate;

[0196]

[0197] m l,line This represents the real-time quality flow rate of the pipeline. This represents the upper limit of the pipeline's mass flow rate.

[0198] From the law of conservation of nodal energy:

[0199]

[0200] These are the mixing temperatures of the water supply pipe nodes and the return water pipe nodes, respectively. These are the mass flow rates of the water supply pipeline node and the return water pipeline node, respectively. These are the temperatures of the supply and return water pipes, respectively; according to the law of conservation of flow at the nodes, we can obtain:

[0201]

[0202] These are the mass flow rates of the heat exchanger and the heat source at different times;

[0203] CHP unit thermo-electric operational domain constraints and ramp-up constraints:

[0204]

[0205] Let the output power of the i-th unit at time t be... This represents the lower limit of the output electrical power of the i-th unit at time t, which is related to thermal power. This represents the upper limit of the output electrical power of the i-th unit at time t, which is related to thermal power.

[0206]

[0207] Let the output thermal power of the i-th unit at time t be... This represents the lower limit of the output thermal power of the i-th unit at time t, which is related to the electrical power. This represents the upper limit of the output thermal power of the i-th unit at time t, which is related to the electrical power.

[0208]

[0209] For the output power of the i-th unit at time t-1, The range of power change for the i-th unit within a time interval.

[0210] A simulation model of an integrated electric and thermal energy network was established. The distribution network topology is shown below. Figure 1 As shown, the topology adopts the Bali Island electric heating system model, which includes two systems, electric and heating, coupled together through combined heat and power (CHP). The network pipeline parameters are shown in Table 1. The network contains two CHP units and an energy hub as the power and heat supply sources, with source 1 serving as the balance node for both power and heat. The supply and return water temperatures of the heating network pipelines at each node at a certain time are shown in Table 2, and the pipeline water mass flow rate at a certain time is shown in Table 3.

[0211] Table 1. Pipeline parameters from the Bali case study.

[0212]

[0213] Table 2 Supply and return water temperatures at nodes in the heating pipeline

[0214]

[0215] Table 3 Mass flow rate of supply and return water in heating pipelines

[0216]

[0217] When launching an attack on the network, the chosen target was the hot-load node 3, where the final value of the attack was doubled. The resulting fluctuations were as follows: Figure 2 As shown. The prediction results of the multiple matching method in this invention are as follows. Figure 3 As shown, this demonstrates the high accuracy of the prediction results obtained by this invention. It can effectively monitor deviations caused by attacks.

Claims

1. An IES attack detection method based on non-invasive thermal load detection modeling, characterized in that, Specifically, the following steps are included: Step 1: To address the thermal inertia of the integrated electric and thermal energy system's heating network, a load redistribution attack model is constructed, specifically as follows: Indoor temperature setpoint attack, also known as an indoor temperature setpoint attack, is a type of HLR attack that targets measured indoor temperatures. It works by tampering with the indoor temperature... With indoor setpoint temperature The mismatch between them induces the heating system to continuously supply inappropriate electricity to the heat load: Where e is a constant, Indoor temperature, For time, The attack time range Set the indoor temperature value. For attack strength parameters, For attack rate parameters, This is the moment the attack begins; Step 2: Based on the load redistribution attack model built in Step 1, a non-intrusive detection model is adopted to establish a dual model on the load side: a radiator model and a building storage heat dissipation model, and to detect the amplification attack deviation of the room temperature state matrix. Step 2.1: Establish a load-side radiator model; The heat exchanger is modeled and analyzed, and the heat fluid energy equation is as follows: in, For heat capacity of the heat fluid, The inlet water temperature, The outlet temperature, The contact area between the fluid and the radiator. The heat transfer coefficient is... The average temperature of the thermal fluid. The effective wall temperature is given by t, where t is time. Wall energy equation: in, The total heat dissipation coefficient of the fins. The effective area of ​​the heat sink For effective fin area, For the heat capacity of the pipe wall; The area of ​​the heat sink is much larger than the area of ​​the pipe wall. Considering only the area of ​​the heat sink, the energy equation for the heat sink is: in, The heat transfer coefficient of the cold fluid, For the effective area of ​​the cold fluid, The average temperature of the cold fluid. For fin heat capacity; Assuming the temperatures of the hot and cold fluids within the heat exchanger are linearly distributed, the energy equation for the cold fluid is: in, For the heat capacity of the cold fluid, The initial temperature of the cold fluid; Step 2.2: Establish a building storage heat dissipation model; A modeling analysis of building heat transfer was performed, and the heat lost from inside the building was: in, To protect the heat dissipation capacity, For the enclosure area, For the heat transfer coefficient of the building envelope, Outdoor temperature This is a correction factor for building orientation. This is the floor correction factor. To improve the heat dissipation capacity of windows, For window area, The heat transfer coefficient of the window. Heat is dissipated through the exchange of air between indoors and outdoors. For the heat capacity of air, The mass of indoor and outdoor air exchanged per unit time; The heat stored in the walls and indoor air is: in, The wall stores heat. For the heat capacity of the wall, For the quality of the wall, Air stores heat. For the heat capacity of air, For air quality; The amount of heat entering the room is: in, To transfer indoor heat to the radiator, The specific heat capacity of water, Water mass flow rate; The difference between the heat input into the room and the heat loss from the room: The state matrix data at each moment includes the difference between the heat input into the room and the heat lost from the room, the inlet temperature, the outlet temperature, and the outdoor temperature; Step 3: Perform state prediction by matching similarity days of the state matrix using a sliding window, and replace the attacked and corrupted state matrix with the predicted state matrix; Step 4: Rapid load recovery of the integrated electrothermal energy system after an attack based on the MILP method.

2. The IES attack detection method based on non-invasive thermal load detection modeling according to claim 1, characterized in that, Step 3 specifically involves: Step 3.1: Under normal operation of the integrated energy system, status data is generated every moment and stored as a historical database. The most similarity matrix is ​​searched in the historical database through multiple matching. Step 3.2: Based on the most similar matrix matched in Step 3.1, compare it with the current state matrix to detect whether the current time is under attack; An attack is detected when the following equation is satisfied; For the detection threshold, The matrix at the current time. The matrix with the most similarity to the historical matrix; Step 3.3: Replace the attacked and contaminated state matrix with the predicted state matrix; For the detected attack matrix, it is replaced with historical data and stored in the historical database to avoid misscheduling by the scheduling center and to provide data support for subsequent state prediction.

3. The IES attack detection method based on non-invasive thermal load detection modeling according to claim 2, characterized in that, Step 3.1: Using the multi-state prediction method, the constructed historical state variables are arranged in chronological order to form a... OK Columnar multidimensional time series matrix ;in Represents the number of node state variables. This represents the number of state variable acquisition points; middle A set of continuous state variable acquisition points serves as the state model of the power grid, characterizing the changing trend of the state variables of each node in the power grid over time; among which... Represents the length of the time window, where the state variables... All A matrix with 1 row and 1 column; a multidimensional time series refers to a set of observations of a set of indicators of an observed object at different times on the same time axis, arranged in a time series. In historical databases, assuming time... Down Let be the current state matrix, and let the matrix adjacent to the current state matrix be... ;by Based on the reference, with a sliding interval length of The premise of time window length Next, Slide backwards along the time axis to obtain the first... The state matrix under each time window is: In the formula ; Indicates less than or equal to The largest integer that can be obtained; By modeling with a sliding time window, the discrete time series is divided into multiple... OK The set of historical state matrices of columns, i.e. ; Generally speaking, the length of the time window Greater than the sliding interval length In the sliding time window model, multiple matching state prediction refers to predicting the state from the historical state matrix. Find the state matrix with the highest similarity to the current state matrix, whose following state matrix has the same trend as the following state matrix, and use it as the state prediction result.

4. The IES attack detection method based on non-invasive thermal load detection modeling according to claim 2, characterized in that, Step 3.1 specifically involves: Step 3.1.1: Detect the current state variables based on similarity analysis; In terms of similarity analysis, based on the actual situation of power grid state variables and the Pearson correlation coefficient, a similarity metric suitable for heating networks is proposed; any two state matrices from the historical database are selected. and The similarity metric is defined as follows: In the formula: , , and Represents the two history matrices OK Column state matrix; Step 3.1.2: For the state matrix set obtained from similarity analysis, the density space clustering algorithm DBSCAN is used to obtain the optimal matching result of the state matrix. The DBSCAN method essentially clusters samples into clusters where points within the same cluster are density-connected, and points in different clusters are not connected. Mahalanobis distance is selected as the clustering metric, and the optimal matching result is obtained based on the DBSCAN results. There are 1 cluster, and the center point of each cluster is... The condition is as shown in the following formula, which ensures that the sum of the Mahalanobis distances between the cluster center and the other state matrices in the cluster is minimized; In the formula: For the first Number of state matrix samples in the cluster; The state matrix; Let the sum of the Mahalanobis distances to the center points of the first type of cluster be ; Exact matching in cluster analysis is achieved by comparing the state matrix at the current time step. The final matching result is achieved by selecting the state matrix with the smallest difference from the state matrices of each cluster center. Step 3.1.3: To measure the overall difference in state variables between two state matrices, a difference index is used to compare and measure the state matrices; the difference is judged by combining the eigenvectors and the distance between the characteristic trends of the matrices. For a state matrix Let the eigenvectors be... for Let represent the maximum and minimum differences between the vector at a certain time point in the heating network state matrix and the average value vector of the state matrix; let the difference between the vector at a certain time point in the state matrix and the average value vector of the state matrix be . ,but: in Indicates the time point , No. The numerical values ​​corresponding to each state variable; feature vector Represented as: For a state matrix Its characteristic trend distance is its 2-norm. The expression is: For any state matrix, it is represented by a pair consisting of the distance between the eigenvector and the eigentrend. express; For any two state matrices and Difference Represented as: The difference degree is used to measure the overall difference between two state matrices; the greater the difference degree, the lower the similarity of corresponding elements in the state matrices, and vice versa.

5. The IES attack detection method based on non-invasive thermal load detection modeling according to claim 1, characterized in that, Step 4 is as follows: Step 4.1: Overall Optimization Objective: The overall optimization objective is to achieve the fastest recovery of the heat load state while keeping costs as low as possible; The overall optimal objective function includes the following costs: 1) gas cost g, 2) grid-related costs e, 3) environmental emission costs, and 4) user dissatisfaction. Adding a load temperature difference integral penalty term to the objective function effectively improves the load recovery speed. in, , , , These represent the hourly electricity price, the hourly gas price, and the prices related to emissions from the power grid and emissions from the gas network, respectively. , , Let represent the power purchased from the main grid in the t-th time interval, and let represent the power associated with the cogeneration unit and boiler unit in the t-th time interval, respectively. , This indicates the time from the start of the attack to the system stabilizing. Indoor temperature setting; In this framework, electricity demand is met through electricity purchased from the main grid, electricity generated by CHP units, electricity generated by WT units, and electricity emitted by ES equipment; in addition, demand response (DR) services are considered to provide load reduction when necessary; reducing some electricity demand within a specified time period and increasing electricity demand at other time intervals; the following formula represents the electricity demand of the developed EH: In the formula Let be the electricity demand in time interval t; Let be the power generation of the wind turbine generator in the t-th time period; and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively. and These represent the power of the energy storage device during the charging and discharging processes in the t-th time interval, respectively. , and These are the electrical efficiency of the transformer, the gas-electric efficiency of the CHP unit, and the electrical efficiency of the wind turbine generator set, respectively. Heat load requirements: In the formula, The heat demand for the t-th time interval; and The respective gas thermal efficiency of the boiler and the gas thermal efficiency of CHP; finally, This refers to the gas purchased from the gas network at regular intervals. and The power is shifted upwards and downwards by the demand response procedure at time interval t, respectively. Step 4.2: Quickly restore constraints; When the system performs optimized scheduling, there will be some upper and lower capacity limits for basic equipment; the power transmission constraints of transmission lines are: To transmit power to transmission lines in real time, This is the lower limit of the transmission power of the transmission line. This is the upper limit of the transmission power of the power transmission line; The relevant parameters on the heating network side indicate that the heat loss in the pipeline is as follows: The water temperature at the pipe outlet. The water inlet temperature of the pipe. For ambient temperature, The heat transfer coefficient of the pipe is... For the length of the pipe, The specific heat capacity of water, For mass flow rate; This represents the real-time quality flow rate of the pipeline. This represents the upper limit of the pipeline's mass flow rate. From the law of conservation of nodal energy: , These are the mixing temperatures of the water supply pipe nodes and the return water pipe nodes, respectively. , These are the mass flow rates of the water supply pipeline node and the return water pipeline node, respectively. , These are the temperatures of the supply and return water pipes, respectively; according to the law of conservation of flow at the nodes, we can obtain: , These are the mass flow rates of the heat exchanger and the heat source at different times; CHP unit thermal-electric operating domain constraints and ramp-up constraints: Let the output power of the i-th unit at time t be... This represents the lower limit of the output electrical power of the i-th unit at time t, which is related to thermal power. This represents the upper limit of the output electrical power of the i-th unit at time t, which is related to thermal power. Let the output thermal power of the i-th unit at time t be... This represents the lower limit of the output thermal power of the i-th unit at time t, which is related to the electrical power. This represents the upper limit of the output thermal power of the i-th unit at time t, which is related to the electrical power. For the output power of the i-th unit at time t-1, The range of power change for the i-th unit within a time interval.

Images