CN116011200A - IES attack detection method based on thermal load non-invasive detection modeling - Google Patents
IES attack detection method based on thermal load non-invasive detection modeling Download PDFInfo
- Publication number
- CN116011200A CN116011200A CN202211654470.0A CN202211654470A CN116011200A CN 116011200 A CN116011200 A CN 116011200A CN 202211654470 A CN202211654470 A CN 202211654470A CN 116011200 A CN116011200 A CN 116011200A
- Authority
- CN
- China
- Prior art keywords
- state
- power
- attack
- heat
- matrix
- Prior art date
- Legal status (The legal status 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 status listed.)
- Pending
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 55
- 239000011159 matrix material Substances 0.000 claims abstract description 135
- 238000000034 method Methods 0.000 claims abstract description 29
- 230000017525 heat dissipation Effects 0.000 claims abstract description 19
- 238000005485 electric heating Methods 0.000 claims abstract description 13
- 238000010438 heat treatment Methods 0.000 claims abstract description 9
- 230000003321 amplification Effects 0.000 claims abstract description 4
- 238000003199 nucleic acid amplification method Methods 0.000 claims abstract description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 46
- 239000012530 fluid Substances 0.000 claims description 33
- 230000005540 biological transmission Effects 0.000 claims description 21
- 239000013598 vector Substances 0.000 claims description 21
- 238000012546 transfer Methods 0.000 claims description 18
- 238000004458 analytical method Methods 0.000 claims description 16
- 238000011084 recovery Methods 0.000 claims description 9
- 230000004044 response Effects 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 7
- 230000005611 electricity Effects 0.000 claims description 7
- 238000005457 optimization Methods 0.000 claims description 6
- 230000009194 climbing Effects 0.000 claims description 3
- 238000007621 cluster analysis Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims description 3
- 238000007599 discharging Methods 0.000 claims description 3
- 238000004146 energy storage Methods 0.000 claims description 3
- 230000007613 environmental effect Effects 0.000 claims description 3
- 238000005338 heat storage Methods 0.000 claims description 3
- 238000012423 maintenance Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 230000005855 radiation Effects 0.000 claims description 3
- 230000006641 stabilisation Effects 0.000 claims description 3
- 238000011105 stabilization Methods 0.000 claims description 3
- 230000003211 malignant effect Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 3
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000005065 mining Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 101150034459 Parpbp gene Proteins 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000003912 environmental pollution Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention designs an IES attack detection method based on thermal load non-invasive detection modeling, and relates to the technical field of comprehensive energy sources; constructing a load redistribution attack model aiming at thermal inertia existing in a heating network of an electric heating comprehensive energy system; on the basis, a non-invasive detection model is adopted, a radiator model and a building storage and heat dissipation model double model are built on the load side, and the amplification attack deviation of a room temperature state matrix is detected; carrying out state prediction by a method of carrying out state matrix similarity day matching through a sliding window, and replacing the state matrix polluted by the attack with a predicted state matrix; based on the MILP method, the load is quickly recovered after the electric heating comprehensive energy system is attacked; the detection method is quick and effective, can realize online detection, and has certain theoretical basis and engineering practical significance.
Description
Technical Field
The invention relates to the technical field of comprehensive energy, in particular to an IES attack detection method based on thermal load non-invasive detection modeling.
Background
The comprehensive energy system (IES) improves the energy utilization efficiency and reduces the environmental pollution through multi-energy complementation and cascade utilization of various energy sources such as electricity, heat, gas and the like. Advanced Information Communication Technology (ICT) plays a vital role in realizing safe, efficient, clean and flexible operation of industrial information systems and promoting deep coupling of industrial information system network systems and physical systems. However, the coupling between various energy systems and their associated information systems increases the complexity of the system, introduces more network vulnerability points, and presents more network security challenges for the safe and efficient operation of IES.
In recent years, failure of an energy system due to network attacks has occurred. Thus, there is an urgent need to study the impact of cyber-threats on IES operation, particularly the cascading effect of cyber-attacks from one system to another.
The existing anomaly detection methods can be classified into deviation-based detection and feature-based detection methods according to the identification basis. Among them, bias-based approaches typically select one or more variables strongly related to the attack based on the defending target system, and an attack is considered to occur when it is detected that the values of these variables deviate too much from the normal range in operation. The feature-based detection method extracts the features of the system in normal operation and under attack through physical mechanism analysis or artificial intelligence method, and judges whether the attack occurs or not through comparing the features in the detection.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides an IES attack detection method based on thermal load non-invasive detection modeling.
An IES attack detection method based on thermal load non-invasive detection modeling specifically comprises the following steps:
step 1: constructing a load redistribution attack model aiming at thermal inertia existing in a heating network of an electric heating comprehensive energy system;
the establishment of the load redistribution attack model is specifically as follows:
indoor temperature set point attack, which is an HLR attack against the measured value of indoor temperature by tampering with indoor temperature
Temperature +.>
Mismatch between, induces the heating system to continue to provide inadequate power to the thermal load:
wherein e is a constant value,
is the indoor temperature, t is the time, Γ a For attack time interval, ++>
Lambda is the indoor temperature set value a1 For attack parameters, a 1 To attack parameter, t 0 The attack starting time is the attack starting time;
step 2: on the basis of constructing a load redistribution attack model in the step 1, a non-invasive detection model is adopted, a radiator model and a building storage and radiation model double model are built on the load side, and the amplification attack deviation of a room temperature state matrix is detected;
step 2.1: establishing a load side radiator model;
modeling analysis is carried out on the radiator, and a thermal fluid energy equation is as follows:
wherein ,Ch Is the heat capacity of the heat fluid, T h,i T is the temperature of the water inlet h,e At the water outlet temperature, A h Is the contact area of the fluid and the radiator, h h T is the heat transfer coefficient h,m For the average temperature of the hot fluid, T w The effective wall temperature is obtained, and t is time;
wall energy equation:
wherein ,Uf Is the total heat dissipation coefficient of the fins, A f,b T is the effective area of the radiating fin f C as effective fin area W Is the heat capacity of the pipe wall;
the area of the radiating fins is far larger than the area of the pipe wall, only the area of the radiating fins is considered, and the energy equation of the radiating fins is as follows:
wherein ,hc For the heat transfer coefficient of the cold fluid, A c For cold fluid effective area, T c,m For cold fluid average temperature, C f The heat capacity of the fins is obtained;
for the cold and hot fluid temperatures within the heat exchanger, assuming a linear distribution, the cold fluid energy equation is:
wherein ,Cc For cold fluid heat capacity, T c,i Is the cold fluid initial temperature;
step 2.2: building a building storage heat dissipation model;
modeling analysis is carried out on heat transfer of a building, and the heat dissipated in a room is as follows:
Q WH ＝S 1 U 1 (T c,i -T out )(1+x 1 )(1+x 2 )
Q CH ＝S 2 U 2 (T c,i -T out )(1+x 1 )(1+x 2 )
Q J ＝C 1 m 1 (T c,i -T out )
wherein ,QWH For enclosing heat dissipation power S 1 Is a space enclosing area, U 1 To enclose the heat transfer coefficient T out Is the outdoor temperature x 1 For building orientation correction factor, x 2 For floor correction factor, Q CH For the heat dissipation power of the window S 2 For the window area, U 2 For window heat transfer coefficient, Q J Heat dissipation for indoor and outdoor gas exchange, C 1 Is the heat capacity of air, m 1 Exchanging quality for indoor and outdoor air in unit time;
the heat stored by the wall body and the indoor air is as follows:
Q T ＝C 1 m 3 T c,i
wherein ,Qq Wall body heat storage, C 2 For heat capacity of wall, m 2 Is wall mass, Q T Air stores heat, C 1 Is the heat capacity of air, m 3 Is the air quality;
the heat transferred into the room is as follows:
Q R ＝Cm 4 (T h,i -T h,e )
wherein ,QR For the heat transferred into the room by the radiator, C is the specific heat capacity of water, m 4 Is the water mass flow;
difference between heat input into the room and heat dissipated from the room:
ΔQ＝Q R -Q WH -Q CH -Q J -Δ(Q q +Q T )
the state matrix data at each moment comprises a difference value between the heat input into the room and the heat dissipated from the room, a water inlet temperature, a water outlet temperature and an outdoor temperature;
step 3: carrying out state prediction by a method of carrying out state matrix similarity day matching through a sliding window, and replacing the state matrix polluted by the attack with a predicted state matrix;
step 3.1: state data can be generated at any moment under the normal operation of the comprehensive energy and stored as a historical database, and the most similar matrix is searched in the historical database through multiple matching;
the constructed historical state quantity is arranged in time sequence by a multi-state prediction method to form a q-row p-column multi-dimensional time sequence matrix T q×p The method comprises the steps of carrying out a first treatment on the surface of the Wherein q represents the number of node state quantity, and p represents the number of state quantity acquisition points; t (T) q×p The collection of the L groups of continuous state quantity acquisition points serves as a state model of the power grid, and the change trend of the state quantity of each node of the power grid along with time is represented; wherein L represents the length of a time window, and the state quantity x is q rows and 1 columns of matrix; the multidimensional time sequence is a digital set which is formed by arranging a series of observation values obtained by a group of indexes of an observation object at each moment on the same time axis according to the time sequence;
in the history database, it is assumed that at time k (X V ) q×L For the current state matrix, let the matrix adjacent to the current state matrix be (X 1 ) q×L The method comprises the steps of carrying out a first treatment on the surface of the By (X) 1 ) q×L For reference, on the premise that the sliding interval length is w and the time window length is L, the method comprises the steps of (X 1 ) q×L Sliding reversely along the time axis to obtain the state matrix X in the jth 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 
Representing less than or equal to
The maximum integer that can be obtained;
by sliding time window modeling, a discrete time sequence is divided into a set of a plurality of q rows and L columns of history state matrix, namely
In general, the time window length L is greater than the sliding interval length w; under the sliding time window model, the multiple matching state prediction refers to searching 1 state matrix with highest similarity with the current state matrix from the historical state matrix X, wherein the following state matrix and the following state matrix have the same change trend, so as to be used as a state prediction result;
step 3.1.1: detecting the state quantity at the current moment based on similarity analysis;
in the aspect of similarity analysis, combining with the actual condition of the state quantity of the power grid, based on the pearson correlation coefficient thought, providing a similarity measurement index applicable to the heat supply network; any two state matrixes A ' = [ a ' in the history database are selected '] q×L And B '= [ B ]'] q×L The similarity metric is defined as follows:
in the formula ：
a 'and B' represent q rows and L columns of state matrices in the two history matrices;
step 3.1.2: for a state matrix set obtained by similarity analysis, obtaining an optimal matching result of a state matrix by adopting a density space clustering algorithm DBSCAN; the DBSCAN method essentially clusters samples into density connection points among the same clusters, and different clusters are not connected; selecting a mahalanobis distance as a clustering metric index, and obtaining K clustering clusters based on a DBSCAN obtained result, wherein the center point C= { C of each clustering cluster 1 ,C 2 ,…,C K The condition satisfied by the following formulaEnsuring the Marshall distance sum of the central point of the cluster and other state matrixes in the cluster to be minimum;
in the formula ：RK The number of the state matrix samples in the K-th cluster; u (U) t Is a state matrix; m (C) K ) The sum of the mahalanobis distance of the center point of the first cluster;
the exact match of the cluster analysis is achieved by comparing the current moment state matrix (X V ) q×L The state matrix with the smallest difference degree is selected as a final matching result to realize the difference degree with the state matrix of each cluster center;
step 3.1.3: for measuring the comprehensive difference of state quantity between two state matrixes, comparing and measuring the state matrixes by adopting a difference index; the difference degree is comprehensively judged through the feature vectors among the matrixes and the feature trend distance result;
for a state matrix X q×L Let feature vector F be (F x ,f y ) Representing the maximum and minimum difference between a vector of a certain moment point in the hot network state matrix and a state matrix average value vector; let the difference between the vector of a certain point of the state matrix and the average value of the state matrix be C 1×L Then:
wherein Xi (t j ) Indicated at the time point t j The value corresponding to the i-th state quantity; the feature vector F is expressed as:
F＝(f x ,f y )＝(max(C) 1×L ,min(C) 1×L )
for a state matrix C q×L The characteristic trend distance is two norms D, and the expression is:
for any state matrix, the state matrix is represented by a binary group G= (F, D) formed by the feature vector and the feature trend distance;
for any 2 state matrices X a and Xb Degree of difference X ab Expressed as:
the comprehensive difference between the 2 state matrixes is measured through the difference degree; the larger the difference degree is, the lower the similarity of the corresponding elements in the state matrix is, and the higher the opposite is;
step 3.2: according to the most similar matrix matched in the step 3.1, comparing the most similar matrix with the state matrix at the current moment to detect whether the current moment is attacked or not;
indicating that an attack is detected when the following equation is satisfied;
K≤max((X v ) 1,end -(X j ) 1,end ,,(X v ) i,end -(X j ) i,end )
k is a detection threshold value, X v For the current moment matrix, X j Is the most similar matrix of history;
step 3.3: replacing the polluted state matrix under attack with a predicted state matrix;
for the detected attack matrix, replacing the attack matrix by using historical data and storing the attack matrix in a historical database, so that a dispatching center is prevented from carrying out error dispatching, and data support is provided for the later state prediction;
step 4: based on the MILP method, the load is quickly recovered after the electric heating comprehensive energy system is attacked;
step 4.1: overall optimization objective: the overall optimization goal is to achieve the fastest recovery of the thermal load state with the lowest maintenance cost;
the overall optimal objective function includes the following costs: 1) gas cost g, 2) grid related cost e, 3) environmental emission cost, 4) user dissatisfaction; the load temperature difference integral penalty term is added in the objective function, so that the load recovery speed is effectively improved;
wherein ,
represents the price per hour, the price per hour of gas and the price related to the emissions produced by the grid and the price related to the emissions produced by the gas network, +.>
The power purchased from the main power grid in the t time interval is represented respectively, and the power associated with the cogeneration unit and the boiler unit in the t time interval is represented. t is t 1 ，t 2 Time t representing start of power attack and system recovery stabilization set An indoor temperature set point;
in this framework, the power demand is satisfied by the power purchased from the main grid, the power generated by the CHP unit, the power generated by the WT, and the power discharged by the ES device; in addition, demand Response (DR) services are also considered to provide load curtailment if necessary; reducing part of electricity demand in a specified time, and increasing the electricity demand at other time intervals; the following formula represents the electrical requirements of the developed EH:
in the formula Pt Electrical Power demand for the t-th time interval;
generating energy of the wind generating set in the t-th time period;
and
By demand response programs at time intervals t, respectivelyMoving power up and down;
and
Respectively the power of the energy storage device in the charging and discharging process in the t-th time interval;
and
The electric efficiency of the transformer, the gas-electric efficiency of the CHP unit and the electric efficiency of the wind generating set are respectively;
demand for thermal load:
in the formula ,Pt Thermal Heat demand for the t-th time interval;
and
The gas-heat efficiency of the boiler and the gas-heat efficiency of the CHP are respectively; finally, let(s)>
Means gas purchased from a gas network at intervals, < >>
and
Moving power up and down by the demand response program at time intervals t, respectively;
step 4.2: quickly recovering constraint conditions;
the base equipment has some upper and lower capacity limits when the system performs optimal scheduling; the power transmission constraint of the transmission line is as follows:
p l,line the power is transmitted in real time for the transmission line,
for the transmission power lower limit of the transmission line, +.>
The upper limit of transmission power of the transmission line is set;
the related parameter conditions of the heat supply network side are that the heat dissipation of the pipeline is as follows:
T end t is the temperature of the water outlet of the pipeline start T is the temperature of a water inlet of a pipeline a Lambda is the heat transfer coefficient of the pipeline, L is the length of the pipeline, C p Is the specific heat capacity of water,
is mass flow;
m l,line is the real-time mass flow rate of the pipeline,
is the upper limit of the pipeline mass flow;
from node conservation of energy law:
CHP unit thermo-electric operable domain constraint and climbing constraint:
The beneficial effects of the invention are as follows:
the invention designs a set of online defense flow for electric heating comprehensive energy network security. The method comprises the steps of building a thermal load model, mining state quantity and processing malignant data, wherein the state quantity mining comprises two parts of historical state quantity acquisition and malignant data processing, and the malignant data processing comprises malignant data detection, rejection and correction. When the state quantity is acquired, the privacy of the heat load user data is considered, non-invasive detection is adopted, and a heat load model is established, so that the user privacy is effectively protected, and the state quantity can be effectively acquired. The thermal load model related to the non-invasive detection model can effectively amplify attack deviation so as to improve the sensitivity of attack detection; the detection method is quick and effective, can realize online detection, and has certain theoretical basis and engineering practical significance.
Drawings
FIG. 1 illustrates an experimental simulation topology of an attack detection method according to an embodiment of the present invention;
FIG. 2 is a comparison graph of IES attack detection methods based on thermal load non-invasive detection modeling according to an embodiment of the present invention;
FIG. 3 is a comparison of attack detection by IES attack detection method based on thermal load non-invasive detection modeling in accordance with an embodiment of the present invention;
Detailed Description
The invention is further described below with reference to the drawings and examples;
an IES attack detection method based on thermal load non-invasive detection modeling specifically comprises the following steps:
step 1: constructing a load redistribution attack model aiming at thermal inertia existing in a heating network of an electric heating comprehensive energy system;
the establishment of the load redistribution attack model is specifically as follows:
indoor temperature set point attack, which is an HLR attack against the measured value of indoor temperature by tampering with indoor temperature
Temperature +.>
Mismatch between, induces the heating system to continue to provide inadequate power to the thermal load:
wherein e is a constant value,
is the indoor temperature, t is the time, Γ a For attack time interval, ++>
Lambda is the indoor temperature set value a1 For attack parameters, a 1 To attack parameter, t 0 The attack starting time is the attack starting time;
step 2: on the basis of constructing a load redistribution attack model in the step 1, a non-invasive detection model is adopted, a radiator model and a building storage and radiation model double model are built on the load side, and the amplification attack deviation of a room temperature state matrix is detected;
step 2.1: establishing a load side radiator model;
modeling analysis is carried out on the radiator, and a thermal fluid energy equation is as follows:
wherein ,Ch Is the heat capacity of the heat fluid, T h,i T is the temperature of the water inlet h,e At the water outlet temperature, A h Is the contact area of the fluid and the radiator, h h T is the heat transfer coefficient h,m For the average temperature of the hot fluid, T w The effective wall temperature is obtained, and t is time;
wall energy equation:
wherein ,Uf Is the total heat dissipation coefficient of the fins, A f,b T is the effective area of the radiating fin f C as effective fin area W Is the heat capacity of the pipe wall;
the area of the radiating fins is far larger than the area of the pipe wall, only the area of the radiating fins is considered, and the energy equation of the radiating fins is as follows:
wherein ,hc For the heat transfer coefficient of the cold fluid, A c For cold fluid effective area, T c,m For cold fluid average temperature, C f The heat capacity of the fins is obtained;
for the cold and hot fluid temperatures within the heat exchanger, assuming a linear distribution, the cold fluid energy equation is:
wherein ,Cc For cold fluid heat capacity, T c,i Is the cold fluid initial temperature;
step 2.2: building a building storage heat dissipation model;
modeling analysis is carried out on heat transfer of a building, and the heat dissipated in a room is as follows:
Q WH ＝S 1 U 1 (T c,i -T out )(1+x 1 )(1+x 2 )
Q CH ＝S 2 U 2 (T c,i -T out )(1+x 1 )(1+x 2 )
Q J ＝C 1 m 1 (T c,i -T out )
wherein ,QWH For enclosing heat dissipation power S 1 Is a space enclosing area, U 1 To enclose the heat transfer coefficient T out Is the outdoor temperature x 1 For building orientation correction factor, x 2 For floor correction factor, Q CH For the heat dissipation power of the window S 2 For the window area, U 2 For window heat transfer coefficient, Q J Heat dissipation for indoor and outdoor gas exchange, C 1 Is the heat capacity of air, m 1 Exchanging quality for indoor and outdoor air in unit time;
the heat stored by the wall body and the indoor air is as follows:
Q T ＝C 1 m 3 T c,i
wherein ,Qq Wall body heat storage, C 2 For heat capacity of wall, m 2 Is wall mass, Q T Air stores heat, C 1 Is the heat capacity of air, m 3 Is the air quality;
the heat transferred into the room is as follows:
Q R ＝Cm 4 (T h,i -T h,e )
wherein ,QR For the heat transferred into the room by the radiator, C is the specific heat capacity of water, m 4 Is the water mass flow;
difference between heat input into the room and heat dissipated from the room:
ΔQ＝Q R -Q WH -Q CH -Q J -Δ(Q q +Q T )
the state matrix data at each moment comprises a difference value between the heat input into the room and the heat dissipated from the room, a water inlet temperature, a water outlet temperature and an outdoor temperature;
step 3: carrying out state prediction by a method of carrying out state matrix similarity day matching through a sliding window, and replacing the state matrix polluted by the attack with a predicted state matrix;
step 3.1: state data can be generated at any moment under the normal operation of the comprehensive energy and stored as a historical database, and the most similar matrix is searched in the historical database through multiple matching;
the constructed historical state quantity is arranged in time sequence by a multi-state prediction method to form a q-row p-column multi-dimensional time sequence matrix T q×p The method comprises the steps of carrying out a first treatment on the surface of the Wherein q represents the number of node state quantity, and p represents the number of state quantity acquisition points; t (T) q×p The collection of the L groups of continuous state quantity acquisition points serves as a state model of the power grid, and the change trend of the state quantity of each node of the power grid along with time is represented; wherein L represents the length of a time window, and the state quantity x is q rows and 1 columns of matrix; the multidimensional time sequence is a digital set which is formed by arranging a series of observation values obtained by a group of indexes of an observation object at each moment on the same time axis according to the time sequence;
in the history database, it is assumed that at time k (X V ) q×L For the current state matrix, let the matrix adjacent to the current state matrix be (X 1 ) q×L The method comprises the steps of carrying out a first treatment on the surface of the By (X) 1 ) q×L For reference, on the premise that the sliding interval length is w and the time window length is L, the method comprises the steps of (X 1 ) q×L Sliding reversely along the time axis to obtain the state matrix X in the jth 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 
Representing less than or equal to
The maximum integer that can be obtained;
by sliding time window modeling, a discrete time sequence is divided into a set of a plurality of q rows and L columns of history state matrix, namely
In general, the time window length L is greater than the sliding interval length w; under the sliding time window model, the multiple matching state prediction refers to searching 1 state matrix with highest similarity with the current state matrix from the historical state matrix X, wherein the following state matrix and the following state matrix have the same change trend, so as to be used as a state prediction result;
step 3.1.1: detecting the state quantity at the current moment based on similarity analysis;
in the aspect of similarity analysis, combining with the actual condition of the state quantity of the power grid, based on the pearson correlation coefficient thought, providing a similarity measurement index applicable to the heat supply network; any two state matrixes A ' = [ a ' in the history database are selected '] q×L And B '= [ B ]'] q×L The similarity metric is defined as follows:
in the formula ：
a 'and B' represent q rows and L columns of state matrices in the two history matrices;
step 3.1.2: for the state matrix set obtained by similarity analysis, a density space clustering algorithm (density spatial clustering of application withnoise, DBSCAN) to obtain an optimal matching result of the state matrix; the DBSCAN method essentially clusters samples into density connection points among the same clusters, and different clusters are not connected; selecting a mahalanobis distance as a clustering metric index, and obtaining K clustering clusters based on a DBSCAN obtained result, wherein the center point C= { C of each clustering cluster 1 ,C 2 ,…,C K The satisfied condition is shown as the following formula, namely, the Markov distance sum of the central point of the cluster and other state matrixes in the cluster is guaranteed to be minimum;
in the formula ：RK The number of the state matrix samples in the K-th cluster; u (U) t Is a state matrix; m (C) K ) The sum of the mahalanobis distance of the center point of the first cluster;
the exact match of the cluster analysis is achieved by comparing the current moment state matrix (X V ) q×L The state matrix with the smallest difference degree is selected as a final matching result to realize the difference degree with the state matrix of each cluster center;
step 3.1.3: for measuring the comprehensive difference of state quantity between two state matrixes, comparing and measuring the state matrixes by adopting a difference index; the difference degree is comprehensively judged through the feature vectors among the matrixes and the feature trend distance result;
for a state matrix X q×L Let feature vector F be (F x ,f y ) Representing the maximum and minimum difference between a vector of a certain moment point in the hot network state matrix and a state matrix average value vector; let the difference between the vector of a certain point of the state matrix and the average value of the state matrix be C 1×L Then:
wherein Xi (t j ) Indicated at the time point t j The value corresponding to the i-th state quantity; the feature vector F is expressed as:
F＝(f x ,f y )＝(max(C) 1×L ,min(C) 1×L )
for a state matrix C q×L The characteristic trend distance is two norms D, and the expression is:
for any state matrix, the state matrix is represented by a binary group G= (F, D) formed by the feature vector and the feature trend distance;
for any 2 state matrices X a and Xb Degree of difference X ab Expressed as:
the comprehensive difference between the 2 state matrixes is measured through the difference degree; the larger the difference degree is, the lower the similarity of the corresponding elements in the state matrix is, and the higher the opposite is;
step 3.2: according to the most similar matrix matched in the step 3.1, comparing the most similar matrix with the state matrix at the current moment to detect whether the current moment is attacked or not;
indicating that an attack is detected when the following equation is satisfied;
K≤max((X v ) 1,end -(X j ) 1,end ,,(X v ) i,end -(X j ) i,end )
k is a detection threshold value, X v For the current moment matrix, X j Is the most similar matrix of history;
step 3.3: replacing the polluted state matrix under attack with a predicted state matrix;
for the detected attack matrix, replacing the attack matrix by using historical data and storing the attack matrix in a historical database, so that a dispatching center is prevented from carrying out error dispatching, and data support is provided for the later state prediction;
step 4: based on the MILP method, the load is quickly recovered after the electric heating comprehensive energy system is attacked;
step 4.1: overall optimization objective: the overall optimization goal is to achieve the fastest recovery of the thermal load state with the lowest maintenance cost;
the overall optimal objective function includes the following costs: 1) gas cost g, 2) grid related cost e, 3) environmental emission cost, 4) user dissatisfaction; the load temperature difference integral penalty term is added in the objective function, so that the load recovery speed is effectively improved;
wherein ,
represents the price per hour, the price per hour of gas and the price related to the emissions produced by the grid and the price related to the emissions produced by the gas network, +.>
The power purchased from the main power grid in the t time interval is represented respectively, and the power associated with the cogeneration unit and the boiler unit in the t time interval is represented. t is t 1 ，t 2 Time t representing start of power attack and system recovery stabilization set An indoor temperature set point;
in this framework, the power demand is satisfied by the power purchased from the main grid, the power generated by the CHP unit, the power generated by the WT, and the power discharged by the ES device; in addition, demand Response (DR) services are also considered to provide load curtailment if necessary; reducing part of electricity demand in a specified time, and increasing the electricity demand at other time intervals; the following formula represents the electrical requirements of the developed EH:
in the formula Pt Electrical Power demand for the t-th time interval;
generating energy of the wind generating set in the t-th time period;
and
Moving power up and down by the demand response program at time intervals t, respectively;
and
Respectively the power of the energy storage device in the charging and discharging process in the t-th time interval;
and
The electric efficiency of the transformer, the gas-electric efficiency of the CHP unit and the electric efficiency of the wind generating set are respectively;
demand for thermal load:
in the formula ,Pt Thermal Heat demand for the t-th time interval;
and
The gas-heat efficiency of the boiler and the gas-heat efficiency of the CHP are respectively; finally, let(s)>
Means gas purchased from a gas network at intervals, < >>
and
Moving power up and down by the demand response program at time intervals t, respectively;
step 4.2: quickly recovering constraint conditions;
the base equipment has some upper and lower capacity limits when the system performs optimal scheduling; the power transmission constraint of the transmission line is as follows:
p l,line the power is transmitted in real time for the transmission line,
for the transmission power lower limit of the transmission line, +.>
The upper limit of transmission power of the transmission line is set; />
The related parameter conditions of the heat supply network side are that the heat dissipation of the pipeline is as follows:
T end t is the temperature of the water outlet of the pipeline start T is the temperature of a water inlet of a pipeline a Lambda is the heat transfer coefficient of the pipeline, L is the length of the pipeline, C p Is the specific heat capacity of water,
is mass flow;
m l,line is the real-time mass flow rate of the pipeline,
is the upper limit of the pipeline mass flow;
from node conservation of energy law:
CHP unit thermo-electric operable domain constraint and climbing constraint:
And establishing an electrothermal comprehensive energy network simulation model. The topological diagram of the distribution network is shown in fig. 1, and the topological structure adopts a pari island electric heating system model, wherein the electric heating system model comprises two electric heating systems, and the electric heating systems are coupled together through cogeneration. The network heat supply network pipeline parameter is shown in table 1, and the network comprises two CHP units and energy hinges as power and heat supply sources, and the electric heat uses a source 1 as a balance node. The temperature of the heat supply and return water of each node heat supply network pipeline at a certain moment is shown in table 2, and the mass flow of the pipeline at a certain moment is shown in table 3.
Table 1 parameters of the tubing for the barisland case study
TABLE 2 supply and return water temperature at nodes in heating pipelines
TABLE 3 mass flow of supply and return water in heating pipeline
In the case of an attack on the network, a network attack with a final value of double the detection value is selected at the thermal load node 3, and the fluctuation caused by the attack is shown in fig. 2. The result of the multiple matching prediction method in the invention is shown in fig. 3. The prediction result of the invention can be shown to have higher precision. Deviations caused by attacks can be effectively monitored.
Claims (7)
1. An IES attack detection method based on thermal load non-invasive detection modeling is characterized by comprising the following steps:
step 1: constructing a load redistribution attack model aiming at thermal inertia existing in a heating network of an electric heating comprehensive energy system;
step 2: on the basis of constructing a load redistribution attack model in the step 1, a non-invasive detection model is adopted, a radiator model and a building storage and radiation model double model are built on the load side, and the amplification attack deviation of a room temperature state matrix is detected;
step 3: carrying out state prediction by a method of carrying out state matrix similarity day matching through a sliding window, and replacing the state matrix polluted by the attack with a predicted state matrix;
step 4: and the load is quickly recovered after the electric heating comprehensive energy system is attacked based on the MILP method.
2. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 1, wherein the building of the load redistribution attack model in step 1 specifically includes:
indoor temperature set point attack, which is an HLR attack against the measured value of indoor temperature by tampering with indoor temperature
Temperature +.>
Mismatch between, induces the heating system to continue to provide inadequate power to the thermal load:
wherein e is a constant value,
is the indoor temperature, t is the time, Γ a For attack time interval, ++>
Lambda is the indoor temperature set value a1 For attack parameters, a 1 To attack parameter, t 0 Is the attack start time.
3. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 1, wherein step 2 specifically comprises:
step 2.1: establishing a load side radiator model;
modeling analysis is carried out on the radiator, and a thermal fluid energy equation is as follows:
wherein ,Ch Is the heat capacity of the heat fluid, T h,i T is the temperature of the water inlet h,e At the water outlet temperature, A h Is the contact area of the fluid and the radiator, h h T is the heat transfer coefficient h,m For the average temperature of the hot fluid, T w The effective wall temperature is obtained, and t is time;
wall energy equation:
wherein ,Uf Is the total heat dissipation coefficient of the fins, A f,b T is the effective area of the radiating fin f C as effective fin area W Is the heat capacity of the pipe wall;
the area of the radiating fins is far larger than the area of the pipe wall, only the area of the radiating fins is considered, and the energy equation of the radiating fins is as follows:
wherein ,hc For the heat transfer coefficient of the cold fluid, A c For cold fluid effective area, T c,m For cold fluid average temperature, C f The heat capacity of the fins is obtained;
for the cold and hot fluid temperatures within the heat exchanger, assuming a linear distribution, the cold fluid energy equation is:
wherein ,Cc For cold fluid heat capacity, T c,i Is the cold fluid initial temperature;
step 2.2: building a building storage heat dissipation model;
modeling analysis is carried out on heat transfer of a building, and the heat dissipated in a room is as follows:
Q WH ＝S 1 U 1 (T c,i -T out )(1+x 1 )(1+x 2 )
Q CH ＝S 2 U 2 (T c,i -T out )(1+x 1 )(1+x 2 )
Q J ＝C 1 m 1 (T c,i -T out )
wherein ,QWH For enclosing heat dissipation power S 1 Is a space enclosing area, U 1 To enclose the heat transfer coefficient T out Is the outdoor temperature x 1 For building orientation correction factor, x 2 For floor correction factor, Q CH For the heat dissipation power of the window S 2 For the window area, U 2 For window heat transfer coefficient, Q J Heat dissipation for indoor and outdoor gas exchange, C 1 Is the heat capacity of air, m 1 Exchanging quality for indoor and outdoor air in unit time;
the heat stored by the wall body and the indoor air is as follows:
Q T ＝C 1 m 3 T c,i
wherein ,Qq Wall body heat storage, C 2 For heat capacity of wall, m 2 Is wall mass, Q T Air stores heat, C 1 Is the heat capacity of air, m 3 Is the air quality;
the heat transferred into the room is as follows:
Q R ＝Cm 4 (T h,i -T h,e )
wherein ,QR For the heat transferred into the room by the radiator, C is the specific heat capacity of water, m 4 Is the water mass flow;
difference between heat input into the room and heat dissipated from the room:
ΔQ＝Q R -Q WH -Q CH -Q J -Δ(Q q +Q T )
the state matrix data at each moment comprises the difference value between the heat input into the room and the heat dissipated from the room, the inlet temperature, the outlet temperature and the outdoor temperature.
4. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 1, wherein step 3 specifically comprises:
step 3.1: state data can be generated at any moment under the normal operation of the comprehensive energy and stored as a historical database, and the most similar matrix is searched in the historical database through multiple matching;
step 3.2: according to the most similar matrix matched in the step 3.1, comparing the most similar matrix with the state matrix at the current moment to detect whether the current moment is attacked or not;
indicating that an attack is detected when the following equation is satisfied;
K≤max((X v ) 1,end -(X j ) 1,end ,,(X v ) i,end -(X j ) i,end )
k is a detection threshold value, X v For the current moment matrix, X j Is the most similar matrix of history;
step 3.3: replacing the polluted state matrix under attack with a predicted state matrix;
for the detected attack matrix, the attack matrix is replaced by historical data and stored in a historical database, so that the error scheduling of a scheduling center is avoided, and data support is provided for the subsequent state prediction.
5. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 4, wherein step 3.1 is to arrange the constructed historical state quantities in time sequence by a multiple state prediction method to form a q-row p-column multidimensional time series matrix T q×p The method comprises the steps of carrying out a first treatment on the surface of the Wherein q represents the number of node state quantity, and p represents the number of state quantity acquisition points; t (T) q×p The collection of L groups of continuous state quantity acquisition points in the middle is used as a state model of the power grid to represent the state quantity of each node of the power gridTrend over time; wherein L represents the length of a time window, and the state quantity x is q rows and 1 columns of matrix; the multidimensional time sequence is a digital set which is formed by arranging a series of observation values obtained by a group of indexes of an observation object at each moment on the same time axis according to the time sequence;
in the history database, it is assumed that at time k (X V ) q×L For the current state matrix, let the matrix adjacent to the current state matrix be (X 1 ) q×L The method comprises the steps of carrying out a first treatment on the surface of the By (X) 1 ) q×L For reference, on the premise that the sliding interval length is w and the time window length is L, the method comprises the steps of (X 1 ) q×L Sliding reversely along the time axis to obtain the state matrix X in the jth 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 
Representing less than or equal to
The maximum integer that can be obtained;
by sliding time window modeling, a discrete time sequence is divided into a set of a plurality of q rows and L columns of history state matrix, namely
In general, the time window length L is greater than the sliding interval length w; in the sliding time window model, multiple matching state prediction refers to searching 1 state matrix with highest similarity with the current state matrix from the historical state matrix X, wherein the following state matrix and the following state matrix have the same change trend, so as to be used as a state prediction result.
6. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 4, wherein step 3.1 specifically comprises:
step 3.1.1: detecting the state quantity at the current moment based on similarity analysis;
in the aspect of similarity analysis, combining with the actual condition of the state quantity of the power grid, based on the pearson correlation coefficient thought, providing a similarity measurement index applicable to the heat supply network; any two state matrixes A ' = [ a ' in the history database are selected '] q×L And B '= [ B ]'] q×L The similarity metric is defined as follows:
in the formula ：
a 'and B' represent q rows and L columns of state matrices in the two history matrices;
step 3.1.2: for a state matrix set obtained by similarity analysis, obtaining an optimal matching result of a state matrix by adopting a density space clustering algorithm DBSCAN; the DBSCAN method essentially clusters samples into density connection points among the same clusters, and different clusters are not connected; selecting a mahalanobis distance as a clustering metric index, and obtaining K clustering clusters based on a DBSCAN obtained result, wherein the center point C= { C of each clustering cluster 1 ,C 2 ,…,C K The satisfied condition is shown as the following formula, namely, the Markov distance sum of the central point of the cluster and other state matrixes in the cluster is guaranteed to be minimum;
in the formula ：RK The number of the state matrix samples in the K-th cluster; u (U) t Is a state matrix; m (C) K ) The sum of the mahalanobis distance of the center point of the first cluster;
the exact match of the cluster analysis is achieved by comparing the current moment state matrix (X V ) q×L Differences from the state matrix of each cluster centerThe dissimilarity degree is realized by selecting a state matrix with the smallest dissimilarity degree as a final matching result;
step 3.1.3: for measuring the comprehensive difference of state quantity between two state matrixes, comparing and measuring the state matrixes by adopting a difference index; the difference degree is comprehensively judged through the feature vectors among the matrixes and the feature trend distance result;
for a state matrix X q×L Let feature vector F be (F x ,f y ) Representing the maximum and minimum difference between a vector of a certain moment point in the hot network state matrix and a state matrix average value vector; let the difference between the vector of a certain point of the state matrix and the average value of the state matrix be C 1×L Then:
wherein Xi (t j ) Indicated at the time point t j The value corresponding to the i-th state quantity; the feature vector F is expressed as:
F＝(f x ,f y )＝(max(C) 1×L ,min(C) 1×L )
for a state matrix C q×L The characteristic trend distance is two norms D, and the expression is:
for any state matrix, the state matrix is represented by a binary group G= (F, D) formed by the feature vector and the feature trend distance;
for any 2 state matrices X a and Xb Degree of difference X ab Expressed as:
the comprehensive difference between the 2 state matrixes is measured through the difference degree; the larger the degree of difference, the lower the similarity of the corresponding elements in the state matrix, and conversely, the higher the degree of difference.
7. The IES attack detection method based on thermal load non-invasive detection modeling according to claim 1, wherein step 4 specifically comprises:
step 4.1: overall optimization objective: the overall optimization goal is to achieve the fastest recovery of the thermal load state with the lowest maintenance cost;
the overall optimal objective function includes the following costs: 1) gas cost g, 2) grid related cost e, 3) environmental emission cost, 4) user dissatisfaction; the load temperature difference integral penalty term is added in the objective function, so that the load recovery speed is effectively improved;
wherein ,
represents the price per hour, the price per hour of gas and the price related to the emissions produced by the grid and the price related to the emissions produced by the gas network, +.>
Respectively representing the power purchased from the main power grid in the t time interval, and representing the power associated with the cogeneration unit and the boiler unit in the t time interval; t is t 1 ，t 2 Time t representing start of power attack and system recovery stabilization set An indoor temperature set point;
in this framework, the power demand is satisfied by the power purchased from the main grid, the power generated by the CHP unit, the power generated by the WT, and the power discharged by the ES device; in addition, demand Response (DR) services are also considered to provide load curtailment if necessary; reducing part of electricity demand in a specified time, and increasing the electricity demand at other time intervals; the following formula represents the electrical requirements of the developed EH:
in the formula Pt Electrical Power demand for the t-th time interval;
generating energy of the wind generating set in the t-th time period;
and
Moving power up and down by the demand response program at time intervals t, respectively;
And
respectively the power of the energy storage device in the charging and discharging process in the t-th time interval;
and
The electric efficiency of the transformer, the gas-electric efficiency of the CHP unit and the electric efficiency of the wind generating set are respectively;
demand for thermal load:
in the formula ,Pt Thermal Heat demand for the t-th time interval;
and
The gas-heat efficiency of the boiler and the gas-heat efficiency of the CHP are respectively; finally, let(s)>
Means gas purchased from a gas network at intervals, < >>
and
Moving power up and down by the demand response program at time intervals t, respectively;
step 4.2: quickly recovering constraint conditions;
the base equipment has some upper and lower capacity limits when the system performs optimal scheduling; the power transmission constraint of the transmission line is as follows:
p l,line the power is transmitted in real time for the transmission line,
for the transmission power lower limit of the transmission line, +.>
The upper limit of transmission power of the transmission line is set;
the related parameter conditions of the heat supply network side are that the heat dissipation of the pipeline is as follows:
T end t is the temperature of the water outlet of the pipeline start T is the temperature of a water inlet of a pipeline a Lambda is the heat transfer coefficient of the pipeline, L is the length of the pipeline, C p Is the specific heat capacity of water,
is mass flow;
m l,line is the real-time mass flow rate of the pipeline,
is the upper limit of the pipeline mass flow;
from node conservation of energy law:
CHP unit thermo-electric operable domain constraint and climbing constraint:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211654470.0A CN116011200A (en) | 2022-12-22 | 2022-12-22 | IES attack detection method based on thermal load non-invasive detection modeling |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211654470.0A CN116011200A (en) | 2022-12-22 | 2022-12-22 | IES attack detection method based on thermal load non-invasive detection modeling |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116011200A true CN116011200A (en) | 2023-04-25 |
Family
ID=86026086
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211654470.0A Pending CN116011200A (en) | 2022-12-22 | 2022-12-22 | IES attack detection method based on thermal load non-invasive detection modeling |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116011200A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116781429A (en) * | 2023-08-24 | 2023-09-19 | 国网冀北电力有限公司 | Method, device and equipment for detecting invisible attack of power system |
| CN118112980A (en) * | 2024-04-28 | 2024-05-31 | 南京邮电大学 | Multi-type intelligent building energy system scheduling method considering FDI attack detection |
-
2022
- 2022-12-22 CN CN202211654470.0A patent/CN116011200A/en active Pending
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116781429A (en) * | 2023-08-24 | 2023-09-19 | 国网冀北电力有限公司 | Method, device and equipment for detecting invisible attack of power system |
| CN116781429B (en) * | 2023-08-24 | 2023-10-31 | 国网冀北电力有限公司 | Method, device and equipment for detecting invisible attack of power system |
| CN118112980A (en) * | 2024-04-28 | 2024-05-31 | 南京邮电大学 | Multi-type intelligent building energy system scheduling method considering FDI attack detection |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116011200A (en) | IES attack detection method based on thermal load non-invasive detection modeling | |
| CN109345012B (en) | Park energy Internet operation optimization method based on comprehensive evaluation indexes | |
| CN113822496A (en) | Multi-unit thermal power plant heat supply mode and parameter online optimization method | |
| CN110094802A (en) | A kind of heat pump and heat storage electric boiler united heat load distribution method and device | |
| CN108303898B (en) | Intelligent scheduling method of novel solar-air energy coupling cold and heat cogeneration system | |
| CN117689178B (en) | Method and device for dispatching and optimizing long-period operation of combined type ground source heat pump system | |
| CN116067188B (en) | Waste heat recovery system for lithium hexafluorophosphate preparation | |
| CN108717488A (en) | A kind of multi-objective optimization design of power method of forced air cooling cooling system heat structure | |
| CN108563877A (en) | The Holistic modeling of solar energy BrLi chiller and optimum control integral method | |
| Liu et al. | Optimization of intelligent heating ventilation air conditioning system in urban building based on BIM and artificial intelligence technology | |
| CN117556696A (en) | Design optimization method of medium-deep geothermal coupling phase-change heat storage system with high calculation efficiency | |
| CN111898856A (en) | Extreme learning machine-based physical-data fusion building analysis method | |
| Rao et al. | Design optimization of rotary regenerator using artificial bee colony algorithm | |
| CN117764456A (en) | Park comprehensive energy system optimization method based on co-evolution algorithm | |
| Guo et al. | A thermal response time ahead energy demand prediction strategy for building heating system using machine learning methods | |
| Chen et al. | GCN-and GRU-based intelligent model for temperature prediction of local heating surfaces | |
| Kou et al. | Generation prediction of ultra-short-term wind farm based on quantum genetic algorithm and fuzzy neural network | |
| Mei et al. | Energy efficiency prediction of screw chillers on BP neural network optimized by improved genetic algorithm | |
| Cao et al. | Energy Management using Optimal Fuzzy Logic Control in Wireless Sensor Network. | |
| Fong et al. | A robust evolutionary algorithm for HVAC engineering optimization | |
| CN111402074B (en) | Comprehensive optimization method for mass energy of circulating water system | |
| CN110595008A (en) | Multi-equipment collaborative optimization method and system for ground source heat pump air conditioning system | |
| CN109657885A (en) | A kind of electric-thermal net coordinated operation optimization method based on more searcher algorithms | |
| CN110516867A (en) | A kind of integrated study load forecasting method based on principal component analysis | |
| CN113077125B (en) | Energy efficiency-considered integrated energy system typical scene generation method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |