CN109818349A - A kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching - Google Patents
A kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching Download PDFInfo
- Publication number
- CN109818349A CN109818349A CN201910188703.4A CN201910188703A CN109818349A CN 109818349 A CN109818349 A CN 109818349A CN 201910188703 A CN201910188703 A CN 201910188703A CN 109818349 A CN109818349 A CN 109818349A
- Authority
- CN
- China
- Prior art keywords
- state
- matrix
- measurement
- current
- node
- 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.)
- Granted
Links
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E40/00—Technologies for an efficient electrical power generation, transmission or distribution
- Y02E40/70—Smart grids as climate change mitigation technology in the energy generation sector
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y04—INFORMATION OR COMMUNICATION TECHNOLOGIES HAVING AN IMPACT ON OTHER TECHNOLOGY AREAS
- Y04S—SYSTEMS INTEGRATING TECHNOLOGIES RELATED TO POWER NETWORK OPERATION, COMMUNICATION OR INFORMATION TECHNOLOGIES FOR IMPROVING THE ELECTRICAL POWER GENERATION, TRANSMISSION, DISTRIBUTION, MANAGEMENT OR USAGE, i.e. SMART GRIDS
- Y04S10/00—Systems supporting electrical power generation, transmission or distribution
- Y04S10/22—Flexible AC transmission systems [FACTS] or power factor or reactive power compensating or correcting units
Abstract
The present invention proposes a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching, comprising: seeks node admittance matrix；SCADA instrument measurement and PMU device measurement are obtained, the state duration set estimated based on admixture is constructed；Form historic state amount database；Smart grid node future state is predicted using multi-dimensional state matrix shiding matching method, obtains prediction result；It whether there is abnormal data using improving in power balance equation detection method detection current system, it is abnormal if it exists, judge whether the exception is mutated by measurement using residual distribution irrelevance detection method to cause, if, then it is modified, bad data is eliminated, provides reliable data supporting for status predication from now on.The present invention has better prediction robustness, can efficiently use the information in historical data base, improve the security reliability of smart grid and resist the ability of pernicious data and bad data.
Description
Technical field
The invention belongs to power system security technologies, and in particular to a kind of power grid based on multi-dimensional state matrix shiding matching
Robust state prediction technique.
Background technique
It promotes with intelligent equipment and constantly rises with Automation of Electric Systems level, conventional electric power system is gradually to intelligence
Power grid carries out transition.Smart grid is to provide Networked Virtual as power grid a new generation critical infrastructures, outstanding feature
Source is in conjunction with physics circle entity device deep layer, by integrated computation, communication network and control technology, in power grid primary system and two
Massive information communication construction and information integration channel are established between subsystem, to realize that information system is merged with the height of physical system
Effective way is provided.
The too busy to get away monitoring control of existing smart grid and data acquisition (SCADA) system and phasor measurement unit (PMU)
The auxiliary of information processing, the analysis and decision of data and order transfer function and regulation central energy management system (EMS).For
Guarantee the safe and reliable of above-mentioned function, need to prevent the generation of bad data or network data attack, therefore regulates and controls center and pass through
The mode of status predication obtains forecast data, to carry out accurate to total system and comprehensively monitor.
However there are a little deficiencies for standing state prediction and bad data detection: 1) for the bad of larger difference degree
Data have good detection effect, and the bad data omission factor small for diversity factor is high；2) lack the evil strong to concealment
Property Data Injection Attacks detection mechanism；3) appearance of measurement catastrophe can make prediction incorrectness, and prediction robustness is poor.
Summary of the invention
Based on the above technical deficiency, the present invention is based on SCADA and PMU hybrid measurement scenes in existing smart grid, propose
A kind of smart grid robust state prediction technique based on multi-dimensional state matrix shiding matching, the proposition of this method are that detection is disliked
Property Data Injection Attacks, smaller difference degree bad data and improve status predication robustness and provide a kind of new solution and think
Road.
A kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching, the specific steps are as follows:
Step 1, the specific location for obtaining SCADA instrument and PMU device, determine current smart grid network parameter and topology
Structure seeks node admittance matrix by self-admittance and the definition of transadmittance or branch and node incidence matrix；
The smart grid network parameter includes: road resistance, branch reactance, transformer voltage ratio；
Step 2 is measured based on the different sample frequency of SCADA and PMU, acquisition SCADA instrument measurement and PMU device
Amount constructs the state duration set estimated based on admixture using power grid sliding time window；
Wherein, SCADA instrument measurement includes that node voltage amplitude measures, node injects active measurement, node injects nothing
Function measurement, the active measurement of Branch Power Flow and the idle measurement of Branch Power Flow；PMU device measurement include node voltage amplitude measure,
Node voltage phase angle measurements, branch current magnitudes measure and branch current phase angle measurements；
Step 2.1, at the k moment, PMU and SCADA instrument is sampled, PMU and SCADA admixture estimation side is passed through
Method obtains k moment lower node quantity of state；
Step 2.2 seeks corresponding measurement h (x) to k moment lower node quantity of state x, inscribes system when as k+ Δ k
In measurement, wherein Δ k be static state cycle estimator；
PMU measurement is inscribed when step 2.3, acquisition k+ Δ k, the measurement that step 2.2 is obtained is as the k+ Δ k moment
SCADA instrument measurement is carried out static state estimation with k+ Δ k moment PMU measurement, is obtained using admixture estimation method
K+ Δ k moment lower node quantity of state xk+Δk。
Step 2.4 is sought inscribing the measurement in system when k+ Δ k*2, and method is with step 2.2, finally in SCADA instrument
In sampling interval, a series of node states are obtained, i.e.,Wherein,It is two SCADA instrument of arbitrary continuation
(such as k and k+1 moment) PMU makees the maximum times of static state estimation in the table sampling time；
Step 2.5, at the time of PMU and SCADA instrument is sampled next time, i.e. k+1 moment repeats step 2.1
The method of~step 2.4 obtains the state duration set based on admixture estimation in a period of time；
Step 3, the state duration set that will acquire based on admixture estimation in a period of time are uploaded to regulation center, structure
It builds according to chronological order, the state duration set estimated based on admixture is formed into historic state amount database；
Step 4, regulation center utilize the sliding of multi-dimensional state matrix based on the information stored in historic state amount database
Method of completing the square predicts smart grid node future state, obtains prediction result；
The multi-dimensional state matrix shiding matching method, the specific steps are as follows:
Step 4.1, building multi-dimensional state matrix: history node state amount in database is sequentially arranged, is constituted
One q row p column multidimensional time-series matrix Tq×p, wherein q represents node state amount number, and p represents quantity of state collection point number,
Tq×pThe set of middle L group continuous node state amount collection point is as a state matrix, and wherein L represents time window length；
Step 4.2, sliding time window modeling, obtain current state matrix and historic state set of matrices: assuming that the k moment
Under (Xv)q×LFor current state matrix, wherein Xv=[Xk-L+1,Xk-L+2,…,Xk], enable the matrix adjacent with current state matrix
For (X1)q×L, i.e. X1=[Xk-2L+1,Xk-2L+2,…,Xk-L].With (X1)q×LOn the basis of, it is w, time window in sliding gap length
Under the premise of length L, by (X1)q×LIt is inversely slided along time shaft, state matrix under the ζ time window can be obtained are as follows:
Xζ=[xk-2L-w(ζ-1)+1,xk-2L-w(ζ-1)+2,…,xk-L-w(ζ-1)]
It is modeled by sliding time window, discrete time series is divided into the historic state matrix of multiple q row L column
Set, wherein time window length L is greater than sliding gap length w.
Step 4.3 passes through state matrix similarity measure index, obtains weight matrix optimal value and optimal time window is long
Degree: current state matrix (X is chosenv)q×L, any two state matrix A in the historic state amount database that selecting step 3 obtains
=[a]q×pWith B=[b]q×p, similarity measure index definition is as follows:
Wherein,According to similarity measure index, for given similarity measure threshold value When, claim state matrix A and B to couple similar matrix each other.State matrix A is chosen as current state
Matrix (Xv)q×L, state matrix B is as historic state matrix.Λ is quantity of state weight matrix, Λj(j=1,2 ..., p }) it is q
Ranks vector, λ indicate each element in Λ, Λ=[Λ1 Λ2…Λp]=[λ]q×p, in the case where no measurement is mutated λ=
1.Weight matrix optimal valueWith optimal time length of windowSelection, meet following condition:
Step 4.4, for meetState matrix collection, using density space clustering algorithm ergodic state square
All matrixes that battle array is concentrated, obtain K clustering cluster.Each clustering cluster central point C={ C1,C2,…,CKMeet condition it is as follows:
Wherein RKFor state matrix sample number in K class cluster, UtFor state matrix, M (CK) it is K class cluster central point geneva
Distance and.
Step 4.5 compares current state matrix (Xv)q×LWith the diversity factor of each cluster centre state matrix, diversity factor is chosen
For a smallest state matrix as matching final result, the following state matrix of the state matrix matched is current shape
State prediction result.
The specific solution procedure of diversity factor is as follows: for a state matrix (X)q×LIf feature vector F is (fx,
fy), indicate the maximum, minimal difference of certain moment point vector and state matrix average value vector in electric network state matrix.If state
The difference of matrix moment point vector and state matrix average value is (C)1×L, then:
Wherein, Xi(tj) indicate in moment point tj, i-th of quantity of state correspond to numerical value.Feature vector F is indicated are as follows:
F=(fx,fy)=(Max (C)1×L,Min(C)1×L)
For a state matrix (X)q×L, characteristic trend distance is secondly norm D, expression formula are as follows:
For any one state matrix, it can pass through the binary group G=of feature vector and characteristic trend distance composition
(F, D) is indicated.For any two state matrix XaAnd Xb, diversity factor δabIt may be expressed as:
Step 5, the hybrid measurement quantity of state obtained based on current state prediction result and predicted time point, using improvement function
Rate equilibrium equation detection method, which detects, whether there is abnormal data in current system, if prediction result meets detection criteria, currently
Quantity of state that admixture is estimated exist it is abnormal, if being unsatisfactory for detection criteria, shape that current admixture is estimated
State amount is normal, specifically:
The power grid for being n for number of nodes, the power balance equation that node injecting power meets are as follows:
S=diag (U) (YU)*
WhereinFor node injecting power vector；For complex voltage vector,
Corresponding node quantity of state；Y is the node admittance matrix of corresponding power grid；* complex conjugate is indicated；Diag (U) is using U as diagonal element
Diagonal matrix.By reference mode, it is set as node 1, corresponding element is separated from S, U, Y, then it can obtain:
Y1=[Y11 Y21…Yn1]T
Power balance equation is convertible into following form:
Un-1=(ATA)-1ATB
Wherein,
Based on above-mentioned formula, the admixture estimated result that current predictive time point obtains acquires U, the node that instrument measures
Injecting power acquires S, substitutes into U for U and S as input quantityn-1=(ATA)-1ATB acquires result and is set asIf passing through step 4
Obtained current state predicts that corresponding result isThen detection criteria is as follows:
Wherein, τ1It is related with rate of failing to report for detection threshold value, and the regulation desired rate of false alarm in center.Meet above formula situation, then
There is exception in the quantity of state that current admixture is estimated, if being unsatisfactory for above formula, normally.
If the quantity of state that step 6, current admixture are estimated is normal, current state prediction result is uploaded to and is gone through
History state quantity data library；If the quantity of state that current admixture is estimated has exception, detected using residual distribution irrelevance
Method judges whether the exception is caused by measurement mutation, if being caused by measurement mutation, then slides to multi-dimensional state matrix
Matching process is modified, and revised result is uploaded to historic state amount database, and historic state amount database is by more
Newly；If not being caused by measurement mutation, then current state prediction result is uploaded into historic state amount database, historic state amount
Database is updated；
The residual distribution irrelevance detection method is specific as follows:
The detection of residual distribution irrelevance predominantly detects vectorThe distribution of middle each element is compared with historical data is distributed
With the presence or absence of larger difference.Residual distribution irrelevance detection formula is shown below:
In formula, Γ is residual between the current instrument measurement that current time anomaly data detection calculates and amendment measurement
Difference vector；For the extent of deviation of element in current time vector Γ；E[(x-μΓ)3] indicate current time vector Γ three ranks
Center away from；E[(x-μΓ)2]3/2Indicate the cube of the standard deviation of current time vector Γ；It indicates in historical dataIt is equal
Value；It indicates in historical dataVariance.τ2Selection with detection rate of false alarm and rate of failing to report it is related.If result is greater than τ2, then different
Regular data is since there are SCADA instrument, there are caused by bad data or pernicious data；As a result it is less than τ2, then abnormal data be
As caused by measurement catastrophe.
It is specific as follows for weight update method in state matrix under measurement catastrophe:
The weight column vector Λ at current time will be corresponded in quantity of state weight matrix ΛkInterior each element becomesIts
Middle i is node state amount number,For predicted state amount,For hybrid measurement estimated state amount；
Step 7, this status predication terminate, and using updated historic state amount database, it is pre- to execute next next state
It surveys.
Advantageous effects:
The present invention provides a kind of smart grid robust state prediction technique based on multi-dimensional state matrix shiding matching, the party
Method for standing state prediction technique prediction robustness difference and raw data detection mechanism lack the perfect pernicious data of detection and compared with
The status of the bad data of small diversity factor is initially set up in conjunction with the background of PMU and SCADA hybrid measurement in current smart grid
SCADA and PMU admixture estimates model, makes to analyze result closer in actual electric network, and SCADA and PMU amount is not present
Measure the influence of clock synchronization；Secondly, this method requires to build using multi-dimensional state matrix shiding matching prediction technique proposed by the present invention
Vertical quantity of state historical data base, and the state estimation result obtained every time is stored to more new database, it is power grid in threadiness
State prediction provides a kind of new means；Later, using improve power balance equation detection method detection current time system in whether
There are abnormal data, for the abnormal data detected, using residual distribution irrelevance detection method judge abnormal data whether by
Measurement mutation causes, if so, the weight matrix in amendment status predication, conversely, bad data is rejected；It finally will amendment
State quantity data afterwards is uploaded to historic state amount database, provides reliable data supporting for status predication from now on.It compares
In previous status predication and raw data detection mechanism, the present invention has better prediction robustness, and can reject smaller
While the bad data of diversity factor and perfect pernicious data, the information in historical data base can be efficiently used, is from now on
Status predication improves more reliable data support, improves the security reliability of smart grid and resists pernicious data and umber of defectives
According to ability, have a good application prospect.
Detailed description of the invention
Fig. 1 is the overview flow chart of the embodiment of the present invention；
Fig. 2 is the power grid sliding time window modeling principle figure of the embodiment of the present invention；
Fig. 3 is the power grid multi-dimensional state matrix shiding matching status predication schematic diagram of the embodiment of the present invention；
Fig. 4 is the multi-dimensional state matrix shiding matching status predication algorithm flow chart of the embodiment of the present invention；
Fig. 5 is a predicted voltage amplitude cluster analysis result figure of the embodiment of the present invention；
Fig. 6 is that a predicted voltage amplitude diversity factor of the embodiment of the present invention analyzes result figure；
Fig. 7 is a predicted voltage phase angle cluster analysis result figure of the embodiment of the present invention；
Fig. 8 is that a predicted voltage phase angle diversity factor of the embodiment of the present invention analyzes result figure；
Fig. 9 is that the quantity of state estimated value of the embodiment of the present invention once predicted and predicted value analyze result figure；
Figure 10 is the difference result figure of state quantity prediction value and estimated value in continuous 3 hours of the embodiment of the present invention, wherein scheming
10 (a) be prediction result of the present invention, and Figure 10 (b) is autoregressive model prediction result.
Specific embodiment
Invention is described further with specific implementation example with reference to the accompanying drawing, one kind is slided based on multi-dimensional state matrix
Matched smart grid robust state prediction technique, as shown in Figure 1, the specific steps are as follows:
Step 1, the specific location for obtaining SCADA instrument and PMU device, determine current smart grid network parameter and topology
Structure seeks node admittance matrix by self-admittance and the definition of transadmittance or branch and node incidence matrix；
The smart grid network parameter includes: road resistance, branch reactance, transformer voltage ratio；
Step 2 is measured based on the different sample frequency of SCADA and PMU, acquisition SCADA instrument measurement and PMU device
Amount constructs the state duration set estimated based on admixture using power grid sliding time window, as shown in Figure 2:
Wherein, SCADA instrument measurement includes that node voltage amplitude measures, node injects active measurement, node injects nothing
Function measurement, the active measurement of Branch Power Flow and the idle measurement of Branch Power Flow；PMU device measurement include node voltage amplitude measure,
Node voltage phase angle measurements, branch current magnitudes measure and branch current phase angle measurements；
Step 2.1, at the k moment, PMU and SCADA instrument is sampled, PMU and SCADA admixture estimation side is passed through
Method obtains k moment lower node quantity of state；
Step 2.2 seeks corresponding measurement h (x) to k moment lower node quantity of state x, inscribes system when as k+ Δ k
In measurement, wherein Δ k is static state cycle estimator, generally less than SCADA instrument sampling period.
PMU measurement is inscribed when step 2.3, acquisition k+ Δ k, the measurement that step 2.2 is obtained is as the k+ Δ k moment
SCADA instrument measurement is carried out static state estimation with k+ Δ k moment PMU measurement, is obtained using admixture estimation method
K+ Δ k moment lower node quantity of state xk+Δk。
Step 2.4 is sought inscribing the measurement in system when k+ Δ k*2, and method is with step 2.2, finally in SCADA instrument
In sampling interval, a series of node states are obtained, i.e.,Wherein,It is two SCADA instrument of arbitrary continuation
(such as k and k+1 moment) PMU makees the maximum times of static state estimation in the table sampling time.
Step 2.5, at the time of PMU and SCADA instrument is sampled next time, i.e. k+1 moment repeats step 2.1
The method of~step 2.4 obtains the state duration set based on admixture estimation in a period of time；
Step 3, the state duration set that will acquire based on admixture estimation in a period of time are uploaded to regulation center, structure
It builds according to chronological order, the state duration set estimated based on admixture is formed into historic state amount database；
Step 4, regulation center utilize the sliding of multi-dimensional state matrix based on the information stored in historic state amount database
Method of completing the square predicts smart grid node future state, obtains prediction result, as shown in Figure 3 and Figure 4:
Step 4.1, building multi-dimensional state matrix: history node state amount in database is sequentially arranged, is constituted
One q row p column multidimensional time-series matrix Tq×p, wherein q represents node state amount number, and p represents quantity of state collection point number,
Tq×pThe set of middle L group continuous node state amount collection point is as a state matrix, and wherein L represents time window length；
Step 4.2, sliding time window modeling, obtain current state matrix and historic state set of matrices: assuming that the k moment
Under (Xv)q×LFor current state matrix, wherein Xv=[Xk-L+1,Xk-L+2,…,Xk], enable the matrix adjacent with current state matrix
For (X1)q×L, i.e. X1=[Xk-2L+1,Xk-2L+2,…,Xk-L].With (X1)q×LOn the basis of, it is w, time window in sliding gap length
Under the premise of length L, by (X1)q×LIt is inversely slided along time shaft, state matrix under the ζ time window can be obtained are as follows:
Xζ=[xk-2L-w(ζ-1)+1,xk-2L-w(ζ-1)+2,…,xk-L-w(ζ-1)]
It is modeled by sliding time window, discrete time series is divided into the historic state matrix of multiple q row L column
Set, wherein time window length L is greater than sliding gap length w.
Step 4.3 passes through state matrix similarity measure index, obtains weight matrix optimal value and optimal time window is long
Degree: current state matrix (X is chosenv)q×L, it is based on Pearson correlation coefficients thought, in conjunction with power grid actual conditions, selecting step 3 is obtained
To historic state amount database in any two state matrix A=[a]q×pWith B=[b]q×p, similarity measure index definition is such as
Under:
Wherein,According to similarity measure index, for given similarity measure threshold value When, claim state matrix A and B to couple similar matrix each other.State matrix A is chosen as current state
Matrix (Xv)q×L, state matrix B is as historic state matrix.Λ is quantity of state weight matrix, Λj(j=1,2 ..., p }) it is q
Ranks vector, Λ=[Λ1Λ2…Λp]=[λ]q×p, λ indicate Λ in each element, in the case where no measurement is mutated λ=
1.Weight matrix optimal valueWith optimal time length of windowSelection, meet following condition:
Step 4.4, for meetState matrix collection, using density space clustering algorithm ergodic state square
All matrixes that battle array is concentrated, obtain K clustering cluster.Each clustering cluster central point C={ C1,C2,…,CKMeet condition it is as follows:
Wherein RKFor state matrix sample number in K class cluster, UtFor state matrix, M (CK) it is K class cluster central point geneva
Distance and.
Step 4.5 compares current state matrix (Xv)q×LWith the diversity factor of each cluster centre state matrix, diversity factor is chosen
For a smallest state matrix as matching final result, the following state matrix of the state matrix matched is current shape
State prediction result.
The specific solution procedure of diversity factor is as follows: for a state matrix (X)q×LIf feature vector F is (fx,
fy), indicate the maximum, minimal difference of certain moment point vector and state matrix average value vector in electric network state matrix.If state
The difference of matrix moment point vector and state matrix average value is (C)1×L, then:
Wherein, Xi(tj) indicate in moment point tj, i-th of quantity of state correspond to numerical value.Feature vector F is indicated are as follows:
F=(fx,fy)=(Max (C)1×L,Min(C)1×L)
For a state matrix (X)q×L, characteristic trend distance is secondly norm D, expression formula are as follows:
For any one state matrix, it can pass through the binary group G=of feature vector and characteristic trend distance composition
(F, D) is indicated.For any two state matrix XaAnd Xb, diversity factor δabIt may be expressed as:
The synthesis difference between two state matrixes is measured by diversity factor.
Step 5, the hybrid measurement quantity of state obtained based on current state prediction result and predicted time point, using improvement function
Rate equilibrium equation detection method, which detects, whether there is abnormal data in current system, if prediction result meets detection criteria, currently
Quantity of state that admixture is estimated exist it is abnormal, if being unsatisfactory for detection criteria, shape that current admixture is estimated
State amount is normal, specifically:
The power grid for being n for number of nodes, the power balance equation that node injecting power meets are as follows:
S=diag (U) (YU)*
WhereinFor node injecting power vector；For complex voltage vector,
Corresponding node quantity of state；Y is the node admittance matrix of corresponding power grid；* complex conjugate is indicated；Diag (U) is using U as diagonal element
Diagonal matrix.By reference mode, it is set as node 1, corresponding element is separated from S, U, Y, then it can obtain:
Y1=[Y11 Y21…Yn1]T
Power balance equation is convertible into following form:
Un-1=(ATA)-1ATB
Wherein,
Based on above-mentioned formula, the admixture estimated result that current predictive time point obtains acquires U, the node that instrument measures
Injecting power acquires S, substitutes into U for U and S as input quantityn-1=(ATA)-1ATB acquires result and is set asIf passing through step 4
Obtained current state predicts that corresponding result isThen detection criteria is as follows:
Wherein, τ1It is related with rate of failing to report for detection threshold value, and the regulation desired rate of false alarm in center.Meet above formula situation, then
There is exception in the quantity of state that current admixture is estimated, if being unsatisfactory for above formula, normally.
If the quantity of state that step 6, current admixture are estimated is normal, current state prediction result is uploaded to and is gone through
History state quantity data library；If the quantity of state that current admixture is estimated has exception, detected using residual distribution irrelevance
Method judges whether the exception is caused by measurement mutation, if being caused by measurement mutation, then slides to multi-dimensional state matrix
Matching process is modified, and revised result is uploaded to historic state amount database, and historic state amount database is by more
Newly；If not being caused by measurement mutation, then current state prediction result is uploaded into historic state amount database, historic state amount
Database is updated；
The residual distribution irrelevance detection method is specific as follows:
The detection of residual distribution irrelevance predominantly detects vectorThe distribution of middle each element is compared with historical data is distributed
With the presence or absence of larger difference.Residual distribution irrelevance detection formula is shown below:
In formula, Γ is residual between the current instrument measurement that current time anomaly data detection calculates and amendment measurement
Difference vector；For the extent of deviation of element in current time vector Γ；E[(x-μΓ)3] indicate current time vector Γ three ranks
Center away from；E[(x-μΓ)2]3/2Indicate the cube of the standard deviation of current time vector Γ；It indicates in historical dataIt is equal
Value；It indicates in historical dataVariance.τ2Selection with detection rate of false alarm and rate of failing to report it is related.If result is greater than τ2, then different
Regular data is since there are SCADA instrument, there are caused by bad data or pernicious data；As a result it is less than τ2, then abnormal data be
As caused by measurement catastrophe.
It is specific as follows for weight update method in state matrix under measurement catastrophe:
The weight column vector Λ at current time will be corresponded in quantity of state weight matrix ΛkInterior each element becomesIts
InFor predicted state amount,For hybrid measurement estimated state amount；
Step 7, this status predication terminate, and using updated historic state amount database, it is pre- to execute next next state
It surveys.
Simulation result explanation:
By taking IEEE-14 node system as an example, electric power that the load data of selection is announced from USA New York electrical management office
Load data.Parameter needed for status predication of the present invention has: similarity measure threshold valueCore radius r=used in clustering
The sliding time window number that returns of 0.0011, minimal amount MinPts=2 used in clustering, each status predication are set as 1000, sliding
Dynamic unit w=1, predicts duration L=4.The every 30s of this experiment acquires one group of data, and 360 collection points are amounted in 3 hours.IEEE-
PMU is deployed in node 2,5,6,7 and 9 in 14 node systems.
As shown in Fig. 5, Fig. 6, Fig. 7, Fig. 8, result is analyzed for primary prediction phase angle clustering and diversity factor.Fig. 5 and Fig. 7
In, dot indicates the first cluster class, and plus sige point indicates the second cluster class, and asterisk point indicates third cluster class.Fig. 9 is the state once predicted
It measures estimated value and predicted value analyzes result.It can be obtained by map analysis, predicted value and estimated value diversity factor are little.
To make result with more generality, continuous 3 hours states are carried out in the case where IEEE-14 node system is without attack context
Prediction, and compared with conventional autoregressive model prediction technique, two methods prediction result is as shown in Figure 10, and Figure 10 (a) is this hair
Bright prediction result, Figure 10 (b) are autoregressive model prediction result.Choose 14 phase angle theta of node14Status predication result and estimation tie
Fruit difference is as ordinate, dark dotted line and light color dotted line respectively indicates status predication result and estimated result difference is equal in Figure 10
Value μ and standard deviation sigma.In addition, each node phase angle absolute error mean value and absolute error standard deviation are as shown in table 1.In conjunction with Figure 10 and
Table 1, it may be verified that prediction technique of the present invention is more preferable compared to conventional autoregressive model predicted method prediction effect.
By above analysis and emulation it can be seen that the invention proposes one kind to be based on multi-dimensional state matrix shiding matching
Smart grid robust state prediction technique.Compared to previous research, the present invention has more preferable in the case where can guarantee forecasting accuracy
Prediction robustness, and the bad data of smaller difference degree can be rejected and while perfect pernicious data, can be efficiently used
Information in historical data base provides data supporting for succeeding state prediction.It is reliable pre- that the present invention facilitates the acquisition of regulation center
Count off evidence, and smart grid can be monitored in real time, it can ensure other applications trouble-free operation in EMS.
IEEE-14 node system voltage phase angle absolute difference analysis of the table 1 based on different prediction techniques
Claims (6)
1. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching, which is characterized in that specific steps
It is as follows:
Step 1, the specific location for obtaining SCADA instrument and PMU device determine current smart grid network parameter and topology knot
Structure seeks node admittance matrix by self-admittance and the definition of transadmittance or branch and node incidence matrix；
Step 2, based on the different sample frequency of SCADA and PMU, obtain SCADA instrument measurement and PMU device measurement, adopt
The state duration set estimated based on admixture is constructed with power grid sliding time window；
Step 3, the state duration set that will acquire based on admixture estimation in a period of time are uploaded to regulation center, and building is pressed
According to chronological order, the state duration set estimated based on admixture is formed into historic state amount database；
Step 4, regulation center utilize multi-dimensional state matrix shiding matching side based on the information stored in historic state amount database
Method predicts smart grid node future state, obtains prediction result；
Step 5, the hybrid measurement quantity of state obtained based on current state prediction result and predicted time point, it is flat using power is improved
It whether there is abnormal data in the equation detection method that weighs detection current system, it is current to mix if prediction result meets detection criteria
Quantity of state that state estimation obtains exist it is abnormal, if being unsatisfactory for detection criteria, quantity of state that current admixture is estimated
Normally；
If the quantity of state that step 6, current admixture are estimated is normal, current state prediction result is uploaded into history shape
State amount database；If there is exception in the quantity of state that current admixture is estimated, using residual distribution irrelevance detection method
Judge whether the exception is caused by measurement mutation, if being caused by measurement mutation, then to multi-dimensional state matrix shiding matching
Method is modified, and revised result is uploaded to historic state amount database, and historic state amount database is updated；If
It is not to be caused by measurement mutation, then current state prediction result is uploaded into historic state amount database, historic state amount data
Library is updated；
Step 7, this status predication terminate, and using updated historic state amount database, execute status predication next time.
2. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching according to claim 1,
It is characterized in that, in step 2, the power grid sliding time window constructs the state duration set estimated based on admixture, specific mistake
Journey is as follows:
Step 2.1, at the k moment, PMU and SCADA instrument is sampled, and is obtained by PMU and SCADA admixture estimation method
To k moment lower node quantity of state；
Step 2.2 seeks corresponding measurement h (x) to k moment lower node quantity of state x, inscribes in system when as k+ Δ k
Measurement, wherein Δ k is static state cycle estimator；
PMU measurement is inscribed when step 2.3, acquisition k+ Δ k, the measurement that step 2.2 is obtained is as k+ Δ k moment SCADA
Instrument measurement carries out static state estimation with k+ Δ k moment PMU measurement, obtains k+ Δ k using admixture estimation method
Moment lower node quantity of state xk+Δk；
Step 2.4 is sought inscribing the measurement in system when k+ Δ k*2, and method is finally sampled in SCADA instrument with step 2.2
In interval, a series of node states are obtained, i.e.,Wherein,It is that two SCADA instrument of arbitrary continuation are adopted
PMU makees the maximum times of static state estimation in the sample time；
Step 2.5, at the time of PMU and SCADA instrument is sampled next time, i.e. k+1 moment repeats step 2.1~step
Rapid 2.4 method obtains the state duration set based on admixture estimation in a period of time.
3. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching according to claim 1,
It is characterized in that, in step 4, the multi-dimensional state matrix shiding matching method, the specific steps are as follows:
Step 4.1, building multi-dimensional state matrix: history node state amount in database is sequentially arranged, a q is constituted
Row p column multidimensional time-series matrix Tq×p, wherein q represents node state amount number, and p represents quantity of state collection point number, Tq×p
The set of middle L group continuous node state amount collection point is as a state matrix, and wherein L represents time window length；
Step 4.2, sliding time window modeling, obtain current state matrix and historic state set of matrices: assuming that inscribing when k
(Xv)q×LFor current state matrix, wherein Xv=[Xk-L+1,Xk-L+2,…,Xk], the order matrix adjacent with current state matrix is
(X1)q×L, i.e. X1=[Xk-2L+1,Xk-2L+2,…,Xk-L], with (X1)q×LOn the basis of, it is w in sliding gap length, time window is long
Under the premise of spending L, by (X1)q×LIt is inversely slided along time shaft, state matrix under the ζ time window can be obtained are as follows:
Xζ=[xk-2L-w(ζ-1)+1,xk-2L-w(ζ-1)+2,…,xk-L-w(ζ-1)]
It is modeled by sliding time window, discrete time series is divided into the collection of the historic state matrix of multiple q row L column
It closes, wherein time window length L is greater than sliding gap length w；
Step 4.3 passes through state matrix similarity measure index, obtains weight matrix optimal value and optimal time length of window: choosing
Take current state matrix (Xv)q×L, any two state matrix A=in the historic state amount database that selecting step 3 obtains
[a]q×pWith B=[b]q×p, similarity measure index definition is as follows:
Wherein,According to similarity measure index, for given similarity measure threshold value When, claim state matrix A and B to couple similar matrix each other, chooses state matrix A as current state matrix
(Xv)q×L, for state matrix B as historic state matrix, Λ is quantity of state weight matrix, Λj(j=1,2 ..., p }) it is q ranks
Vector, λ indicate each element in Λ, Λ=[Λ1 Λ2 … Λp]=[λ]q×p, weight matrix optimal valueAnd optimal time
Length of windowSelection, meet following condition:
Step 4.4, for meetState matrix collection, using density space clustering algorithm ergodic state matrix stack
In all matrixes, obtain K clustering cluster, each clustering cluster central point C={ C1,C2,…,CKMeet condition it is as follows:
Wherein RKFor state matrix sample number in K class cluster, UtFor state matrix, M (CK) it is K class cluster central point mahalanobis distance
With；
Step 4.5 compares current state matrix (Xv)q×LWith the diversity factor of each cluster centre state matrix, it is minimum to choose diversity factor
A state matrix as matching final result, the following state matrix of the state matrix matched is that current state is pre-
Survey result；
The specific solution procedure of diversity factor is as follows: for a state matrix (X)q×LIf feature vector F is (fx,fy), it indicates
The maximum, minimal difference of certain moment point vector and state matrix average value vector in electric network state matrix, if when state matrix
The difference of punctum vector and state matrix average value is (C)1×L, then:
Wherein, Xi(tj) indicate in moment point tj, i-th of quantity of state correspond to numerical value, and feature vector F is indicated are as follows:
F=(fx,fy)=(Max (C)1×L,Min(C)1×L)
For a state matrix (X)q×L, characteristic trend distance is secondly norm D, expression formula are as follows:
For any one state matrix, it can pass through the binary group G=(F, D) of feature vector and characteristic trend distance composition
It indicates, for any two state matrix XaAnd Xb, diversity factor δabIt may be expressed as:
4. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching according to claim 1,
It is characterized in that, in step 5, the improvement power balance equation detection method specifically:
The power grid for being n for number of nodes, the power balance equation that node injecting power meets are as follows:
S=diag (U) (YU)*
WhereinFor node injecting power vector；It is corresponding for complex voltage vector
Node state amount；Y is the node admittance matrix of corresponding power grid；* complex conjugate is indicated；Diag (U) is using U as the diagonal of diagonal element
Reference mode is set as node 1 by matrix, and corresponding element is separated from S, U, Y, then can be obtained:
Y1=[Y11 Y21 … Yn1]T
Power balance equation is convertible into following form:
Un-1=(ATA)-1ATB
Wherein,
Based on above-mentioned formula, the admixture estimated result that current predictive time point obtains acquires U, the node injection that instrument measures
Power acquires S, substitutes into U for U and S as input quantityn-1=(ATA)-1ATB acquires result and is set asIf being obtained by step 4
Current state predict that corresponding result isThen detection criteria is as follows:
Wherein, τ1It is related with rate of failing to report for detection threshold value, and the regulation desired rate of false alarm in center, meet above formula situation, then it is current mixed
There is exception in the quantity of state that conjunction state is estimated, if being unsatisfactory for above formula, normally.
5. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching according to claim 1,
It is characterized in that, in step 6, the residual distribution irrelevance detection method is specific as follows:
The detection of residual distribution irrelevance predominantly detects vectorThe distribution of middle each element with historical data be distributed compared with whether
There are larger difference, residual distribution irrelevance detection formula is shown below:
In formula, Γ be current time anomaly data detection calculate current instrument measurement and amendment measurement between residual error to
Amount；For the extent of deviation of element in current time vector Γ；E[(x-μΓ)3] indicate current time vector Γ three rank centers
Away from；E[(x-μΓ)2]3/2Indicate the cube of the standard deviation of current time vector Γ；It indicates in historical dataMean value；
It indicates in historical dataVariance, τ2Selection with detection rate of false alarm and rate of failing to report it is related, if result be greater than τ2, then abnormal number
According to being since there are SCADA instrument, there are caused by bad data or pernicious data；As a result it is less than τ2, then abnormal data be due to
Caused by measurement catastrophe；
The multi-dimensional state matrix shiding matching method, specific as follows:
It is specific as follows for weight update method in state matrix under measurement catastrophe:
The weight column vector Λ at current time will be corresponded in quantity of state weight matrix ΛkInterior each element becomesWherein i is
Node state amount number,For predicted state amount,For hybrid measurement estimated state amount.
6. a kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching according to claim 1,
It is characterized in that, the SCADA instrument measurement includes that node voltage amplitude measures, node injects active measurement, node injects nothing
Function measurement, the active measurement of Branch Power Flow and the idle measurement of Branch Power Flow；The PMU device measurement includes node voltage amplitude amount
Survey, node voltage phase angle measurements, branch current magnitudes measure and branch current phase angle measurements.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910188703.4A CN109818349B (en) | 2019-03-13 | 2019-03-13 | Power grid robust state prediction method based on multidimensional state matrix sliding matching |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910188703.4A CN109818349B (en) | 2019-03-13 | 2019-03-13 | Power grid robust state prediction method based on multidimensional state matrix sliding matching |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109818349A true CN109818349A (en) | 2019-05-28 |
| CN109818349B CN109818349B (en) | 2022-04-22 |
Family
ID=66608808
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910188703.4A Active CN109818349B (en) | 2019-03-13 | 2019-03-13 | Power grid robust state prediction method based on multidimensional state matrix sliding matching |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109818349B (en) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110620399A (en) * | 2019-08-30 | 2019-12-27 | 北方工业大学 | Inverter parallel control method and system based on robust residual generator |
| CN111126846A (en) * | 2019-12-24 | 2020-05-08 | 广东电网有限责任公司 | Method for evaluating differentiation state of overhead transmission line |
| EP3770785A1 (en) * | 2019-07-25 | 2021-01-27 | Siemens Aktiengesellschaft | Method for detecting bad data injections attacks in an industial control system |
| WO2023065905A1 (en) * | 2021-10-19 | 2023-04-27 | International Business Machines Corporation | Pattern detection and prediction using time series data |
| CN116226239A (en) * | 2023-05-06 | 2023-06-06 | 成都瑞雪丰泰精密电子股份有限公司 | Data-driven-based state monitoring method for spindle system of machining center |
| CN116319377A (en) * | 2023-05-15 | 2023-06-23 | 南京邮电大学 | Distributed dynamic state estimation method for power distribution network for resisting network attack |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103324847A (en) * | 2013-06-17 | 2013-09-25 | 西南交通大学 | Method for detecting and identifying dynamic bad data of electric power system |
| CN103326358A (en) * | 2013-06-17 | 2013-09-25 | 西南交通大学 | Electric power system dynamic state estimation method based on synchronous phase-angle measuring device |
| CN103974311A (en) * | 2014-05-21 | 2014-08-06 | 哈尔滨工业大学 | Condition monitoring data stream anomaly detection method based on improved gaussian process regression model |
| US20170146585A1 (en) * | 2015-11-25 | 2017-05-25 | Hitachi, Ltd. | Estimating the locations of power system events using pmu measurements |
| CN107370150A (en) * | 2017-09-06 | 2017-11-21 | 清华大学 | The Power system state estimation Bad data processing method measured based on synchronized phasor |
| CN107453484A (en) * | 2017-08-24 | 2017-12-08 | 国网辽宁省电力有限公司 | A kind of SCADA data calibration method based on WAMS information |
| CN109146336A (en) * | 2018-10-11 | 2019-01-04 | 厦门大学 | A kind of electric system robust exponentially stabilization method based on t distribution |
| CN109274097A (en) * | 2018-11-16 | 2019-01-25 | 四川大学 | A kind of electric power system transient stability method for situation assessment based on Random Matrices Theory |
| US20190056436A1 (en) * | 2016-05-13 | 2019-02-21 | Hitachi, Ltd. | Similarity detection of abnormal waveforms using pmu measurement |
-
2019
- 2019-03-13 CN CN201910188703.4A patent/CN109818349B/en active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103324847A (en) * | 2013-06-17 | 2013-09-25 | 西南交通大学 | Method for detecting and identifying dynamic bad data of electric power system |
| CN103326358A (en) * | 2013-06-17 | 2013-09-25 | 西南交通大学 | Electric power system dynamic state estimation method based on synchronous phase-angle measuring device |
| CN103974311A (en) * | 2014-05-21 | 2014-08-06 | 哈尔滨工业大学 | Condition monitoring data stream anomaly detection method based on improved gaussian process regression model |
| US20170146585A1 (en) * | 2015-11-25 | 2017-05-25 | Hitachi, Ltd. | Estimating the locations of power system events using pmu measurements |
| US20190056436A1 (en) * | 2016-05-13 | 2019-02-21 | Hitachi, Ltd. | Similarity detection of abnormal waveforms using pmu measurement |
| CN107453484A (en) * | 2017-08-24 | 2017-12-08 | 国网辽宁省电力有限公司 | A kind of SCADA data calibration method based on WAMS information |
| CN107370150A (en) * | 2017-09-06 | 2017-11-21 | 清华大学 | The Power system state estimation Bad data processing method measured based on synchronized phasor |
| CN109146336A (en) * | 2018-10-11 | 2019-01-04 | 厦门大学 | A kind of electric system robust exponentially stabilization method based on t distribution |
| CN109274097A (en) * | 2018-11-16 | 2019-01-25 | 四川大学 | A kind of electric power system transient stability method for situation assessment based on Random Matrices Theory |
Non-Patent Citations (3)
| Title |
|---|
| 周建平: "电力系统广域量测不良数据处理研究", 《科技资讯》 * |
| 李强等: "基于混合量测的电力系统状态估计混合算法", 《电力系统自动化》 * |
| 薛辉等: "基于PMU量测数据和SCADA数据融合的电力系统状态估计方法", 《电网技术》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3770785A1 (en) * | 2019-07-25 | 2021-01-27 | Siemens Aktiengesellschaft | Method for detecting bad data injections attacks in an industial control system |
| US11586190B2 (en) | 2019-07-25 | 2023-02-21 | Siemens Aktiengesellschaft | Method for operating a technical or non-technical system, and facility for such systems |
| CN110620399A (en) * | 2019-08-30 | 2019-12-27 | 北方工业大学 | Inverter parallel control method and system based on robust residual generator |
| CN111126846A (en) * | 2019-12-24 | 2020-05-08 | 广东电网有限责任公司 | Method for evaluating differentiation state of overhead transmission line |
| WO2023065905A1 (en) * | 2021-10-19 | 2023-04-27 | International Business Machines Corporation | Pattern detection and prediction using time series data |
| CN116226239A (en) * | 2023-05-06 | 2023-06-06 | 成都瑞雪丰泰精密电子股份有限公司 | Data-driven-based state monitoring method for spindle system of machining center |
| CN116319377A (en) * | 2023-05-15 | 2023-06-23 | 南京邮电大学 | Distributed dynamic state estimation method for power distribution network for resisting network attack |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109818349B (en) | 2022-04-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109818349A (en) | A kind of power grid robust state prediction technique based on multi-dimensional state matrix shiding matching | |
| CN113156917A (en) | Power grid equipment fault diagnosis method and system based on artificial intelligence | |
| CN109921415A (en) | A kind of pernicious online defence method of Data Injection Attacks of power grid towards hybrid measurement | |
| CN114219147A (en) | Power distribution station fault prediction method based on federal learning | |
| Bastos et al. | Machine learning-based prediction of distribution network voltage and sensor allocation | |
| CN110837915A (en) | Low-voltage load point prediction and probability prediction method for power system based on hybrid integrated deep learning | |
| CN109214565A (en) | A kind of subregion system loading prediction technique suitable for the scheduling of bulk power grid subregion | |
| CN110363334A (en) | Grid-connected grid line loss prediction technique based on Grey Neural Network Model | |
| CN111861013A (en) | Power load prediction method and device | |
| Eikeland et al. | Detecting and interpreting faults in vulnerable power grids with machine learning | |
| Kaplan et al. | Fault diagnosis of smart grids based on deep learning approach | |
| Yan et al. | Short-term load forecasting of smart grid based on load spatial-temporal distribution | |
| Rodrigues et al. | A system for analysis and prediction of electricity-load streams | |
| CN116956203B (en) | Method and system for measuring action characteristics of tapping switch of transformer | |
| Reno et al. | Machine learning for rapid qsts simulations using neural networks | |
| Biyun et al. | A Reliability Forecasting Method for Distribution Network Based on Data Mining | |
| CN116167465A (en) | Solar irradiance prediction method based on multivariate time series ensemble learning | |
| CN115994605A (en) | Multi-data fusion photovoltaic power prediction algorithm for comprehensive meteorological factor data | |
| Li et al. | Applications of entropy principles in power systems: A survey | |
| Wang et al. | Data augmentation for intelligent manufacturing with generative adversarial framework | |
| CN114118279A (en) | Method for identifying abnormal data of residential electricity load | |
| Su et al. | On Machine Learning Apporaches towards Dissolved Gases Analyses of Power Transformer Oil Chromatography | |
| McDonald et al. | Investigating Machine Learning for Anomaly Detection in Phasor Measurement Unit Data | |
| Liu et al. | Ultra-short-term wind power forecasting based on stacking model | |
| Dong et al. | Short-Term Wind Power Forecasting with Combined Prediction Based on Chaotic Analysis |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |