CN109657613A - Large scale electric network abnormal load recognition methods based on power method and parallel computing - Google Patents

Large scale electric network abnormal load recognition methods based on power method and parallel computing Download PDF

Info

Publication number
CN109657613A
CN109657613A CN201811556503.1A CN201811556503A CN109657613A CN 109657613 A CN109657613 A CN 109657613A CN 201811556503 A CN201811556503 A CN 201811556503A CN 109657613 A CN109657613 A CN 109657613A
Authority
CN
China
Prior art keywords
subregion
matrix
formula
power
signal
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
Application number
CN201811556503.1A
Other languages
Chinese (zh)
Other versions
CN109657613B (en
Inventor
韩松
李洪乾
周忠强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guizhou University
Original Assignee
Guizhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Guizhou University filed Critical Guizhou University
Priority to CN201811556503.1A priority Critical patent/CN109657613B/en
Publication of CN109657613A publication Critical patent/CN109657613A/en
Application granted granted Critical
Publication of CN109657613B publication Critical patent/CN109657613B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Electricity, gas or water supply
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • 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 large scale electric network abnormal load recognition methods based on power method and parallel computing that the invention discloses a kind of, comprising the following steps: step 1: synchronously constructing each partition data source matrixX s,z ；Step 2: determining each zone time window widthT, and set sampling start timet 1 ；Step 3: synchronously obtaining each subregion sliding window matrixX z ；Step 4: by each subregion sliding window matrixX z It is synchronously carried out standardization, obtains the non-Hermitian matrix of each scoping rulesX n,z ；Step 5: synchronously obtaining each subregion sample covariance matrixS z ；Step 6: each subregion sample covariance matrix maximum eigenvalue is quickly estimated using power methodλ max,z ；Step 7: each subregion estimates current time signal-to-noise ratio respectivelyρ z, to obtain the dynamic threshold of the sample covariance matrix maximum eigenvalue of respective partitionγ z ；Step 8: the extremely out-of-limit differentiation of electric network state.The present invention has the characteristics of capable of being obviously improved computational efficiency, the applicability that large scale electric network is applied in enhancing.

Description

Large scale electric network abnormal load recognition methods based on power method and parallel computing
Technical field
The invention belongs to power grid abnormality detection technical field, specifically design a kind of big based on power method and parallel computing Scale power grid abnormal load recognition methods.
Background technique
With reaching its maturity for Wide Area Measurement System (Wide Area Measurement System, WAMS), Yi Jizhi The data volume of the continuous evolution of energy power grid, explosive growth proposes challenge to electric network data processing and knowledge.It will count greatly Power System Analysis is introduced according to thinking, Extracting Knowledge and its value application from electric power data, by High Performance Computing reality The big data thinking analysis and evaluation of existing operation of power networks state, to ensuring that power system security of new generation has important theoretical anticipate Justice.
Random Matrices Theory (Random Matrix Theory, RMT) is a kind of method with universality, without detailed Physical model can recognize the behavioural characteristic of complication system from higher-dimension angle.On the one hand, from the Power System Analysis based on RMT Theory and methods are in progress from the point of view of angle, have document to propose using spectral radius (Mean Spectral average in RMT Radius, MSR) it is used as correlation metric, the inner link carried out between various influence factors and distribution Running State is dug Pick research.There is document proposition to comprehensively consider power grid historical data and real time data, it is quiet to propose a kind of power grid based on MSR index The method that state stablizes Situation Assessment.There is document to propose further to propose using MSR index analysis electric network state a kind of power grid to disturb The method of dynamic positioning.There is document to propose to propose one kind according to RMT principle based on sample covariance matrix maximum eigenvalue (Maximum Eigenvalue of Sample Covariance Matrix's, MESCM) is suitable for low signal-to-noise ratio scene Power grid abnormality recognition methods.However, analysis of cases mostly uses greatly the minisystem of dozens of bus to calculate in above-mentioned document Example, not yet to such method large-scale electrical power system adaptability, MSR such as based on All Eigenvalues solution technique or The problems such as MESCM index computational efficiency is not high, and universe solution efficiency is relatively low conducts a research.
Summary of the invention
It is an object of the invention to overcome disadvantages mentioned above, a kind of efficiency that can improve load abnormality recognition method is proposed, and Promote the large scale electric network abnormal load based on power method and parallel computing of its applicability in large-scale electrical power system Recognition methods.The thinking that calculating synchronous with sub-area division is specifically sought from the partial feature value of big dimension matrix, from calculation Law theory angle accelerates estimation MESCM index by power method, realizes respectively from technology angle by multi-core parallel concurrent computing technique The plesiochronous quick obtaining of subregion index.
A kind of large scale electric network abnormal load recognition methods based on power method and parallel computing of the invention, including with Lower step:
Step 1: synchronously constructing each partition data source matrix Xs,z.There is the extensive power train of w subregion for one System can be acquired by each subregion PMU and be uploaded to the voltage magnitude or phase angle time series data of affiliated subregion dispatching control center main website, Respective partition data source matrix X is constructed respectivelys,z, wherein partition number z=1,2 ..., w.Xs,zExpression such as formula (1) shown in.
Xs,zThere is only the interference of random noise, also receive influence caused by load random fluctuation, abnormality detection mould Type can be expressed as shown in formula (2),
Xs,z=(1+ ψz)Xp,z+mz×ηz (2)
X in formulap,zFor not by the signal matrix of noise pollution, ΨzFor load Stochastic Volatility, fluctuation range is set as ± 1%.ηzFor noise matrix, mzFor noise amplitude.
Step 2: determining each zone time window width T, and set sampling start time t1；
Step 3: synchronously obtaining each subregion sliding window matrix Xz.Parallel computing is relied on, by the calculating of each subregion Task distributes to the completion of its region computing terminal.Using sliding window technique, each zone time window width T obtained by step 2 With sampling start time t1, from the data source matrix X of each subregions,zSynchronously obtain the N of respective partitionz× T ties up sliding window square Battle array Xz, wherein NzIt is the dimension of sampled data, unit: a；T is sliding window width, unit: a；Matrix ranks ratio cz, such as formula (3) shown in,
And cz∈ (0, ∞) is ratio, no unit.
Step 4: each subregion sliding window matrix XzStandardization.In each subregion computing terminal, to sliding window matrix Xz It is normalized, obtains the non-Hermite Matrix of standard (Non-Hermitian Matrix) XN, z, as shown in formula (4).
X in formulazi=(xzi,1,xzi,2,…,xzi,T)；μ(xzi)、σ(xzi) it is respectively xziMean value and standard deviation；μ (xn,zi)、σ(xn,zi) it is respectively non-Hermitian row matrix vector xn,ziMean value and standard deviation.
Step 5: synchronously calculating each subregion sample covariance matrix Sz.Seek the non-Hermite Matrix X of standardn,zSample association Variance matrix Sz, as shown in formula (5),
Subscript H indicates complex conjugate transposition in formula.
Step 6: each subregion sample covariance matrix maximum eigenvalue λmax,zEstimation.In each subregion computing terminal, make With each subregion sample covariance matrix S of power method iterative estimationzMaximum eigenvalue λmax,z, as shown in formula (6),
U in formulak,zIndicate maximum eigenvalue character pair vector in iterative process, vk,zIndicate " normalization " after iteration to Amount, v0,z={ 1,1 ..., 1 }Nz×1.When k is sufficiently big or | | max (uk,z)-max(uk-1,z) | | when < ε, max (uk,z)≈λmax,z。
Step 7: seeking each subregion MESCM dynamic threshold γz。
Step 7.1: the noise score using the classical spectrum estimate method after the correction of Kaiser window function to each subregion Do not estimated, by taking any one subregion as an example.
(1) assume there is N >=1 PMU, wherein the ac voltage signal of i-th (i=1 ... ..., N) a PMU is xi, signal number It is T according to length.Original signal x is filtered by formula (4)iIn DC component obtain xfi。
(2) Kaiser window can customize one group of adjustable window function, at that time domain representation such as formula (8).
In formula, I0(β) is first kind deformation zero Bessel function, the accuracy rate highest wherein estimated when β=38.Simultaneously Calculate window values wi(n) corresponding mean μ (wi), root-mean-square value wRmsi。
(3) by time-domain signal xfi(Discrete Time Fourier is converted by discrete time Fourier Transform, DTFT) it is transformed to frequency-region signal xi(ejω), power spectral density P is calculated in conjunction with Kaiser window valuesi(ejω), Such as formula (9).
Pi(ejω)=(1/T)/| xi(ejω)|2 (9)
(4) when calculating actual signal power and noise power, segment processing need to be carried out to frequency spectrum, calculates power spectral density pair Answer bandwidth such as formula (10).
bwi=wRmsi/(T×μ(wi)2) (10)
(5) every section of bandwidth corresponds to power density average amplitude such as formula (11).
(6) actual signal power P is obtained by Periodical Maximum Methodri, noise power PniSuch as formula (12).
(7) signal-to-noise ratio that current PMU signal can be obtained using the ratio of actual signal power and noise power, such as formula (13)。
(8) when noise amplitude is without acute variation, global signal-to-noise ratio is indicated with ρ, such as formula (14).
Step 7.2, using the threshold gamma setting method based on Spiked model, be equally with wherein any one subregion Example, the setting of threshold gamma such as formula (15),
In formula, ρ is the global signal-to-noise ratio (snr) estimation value of the current time subregion.C is that the dimension of the subregion holds ratio.α(0≤α≤1) For proportionality coefficient, it can be adjusted according to sliding window width T, generally take α=0.5.
Step 8: the extremely out-of-limit differentiation of electric network state.Judge whether maximum eigenvalue λmax,zGreater than threshold gammazIf λmax,z ≥γzIt sets up, then determines that power grid generating state is abnormal, give a warning.Otherwise, current state without exception, enables i=i+1, returns to step Rapid 3 continue to execute abnormal state testing process.
Wherein, the selection of sampled data type described in step 1 can be only bus voltage amplitude V.
Compared with prior art, the present invention having the advantage that
1, compared to the average spectral radius MSR analytic approach of traditional calculations All Eigenvalues, the present invention is based on power method estimating portion The large scale electric network abnormal load method for quickly identifying of dtex value indicative is obviously improved in computational efficiency.
2. the method for the present invention obtains ideal large-scale electrical power system using zone Divided Parallel Calculation technology Acceleration effect.
Detailed description of the invention
Fig. 1 is flow chart of the invention；
Fig. 2 is the main interconnection of 420 machine of Poland, 2736 bus-bar system and subregion schematic diagram in embodiment；
Fig. 3 is the first subregion MESCM of the case one in embodiment and its result of dynamic threshold；
Fig. 4 is other subregions MESCM of the case one in embodiment and its result of dynamic threshold；
Fig. 5 is the MESCM deviation ratio result of the case two in embodiment；
Fig. 6 is 54 machine of IEEE, a 118 bus-bar system wiring schematic diagram of the case three in embodiment.
Specific embodiment
It is described the specific embodiments of the present invention in detail below in conjunction with drawings and examples, but the present invention is not by described specific Embodiment is limited.
As shown in Figure 1, of the invention is a kind of based on the identification of the large scale electric network abnormal load of power method and parallel computing Method, comprising the following steps:
Step 1: synchronously constructing each partition data source matrix Xs,z.There is the extensive power train of w subregion for one System can be acquired by each subregion PMU and be uploaded to the voltage magnitude or phase angle time series data of affiliated subregion dispatching control center main website, Respective partition data source matrix X is constructed respectivelys,z, wherein partition number z=1,2 ..., w.Xs,zExpression such as formula (1) shown in.
Xs,zThere is only the interference of random noise, also receive influence caused by load random fluctuation, abnormality detection mould Type can be expressed as shown in formula (2),
Xs,z=(1+ ψz)Xp,z+mz×ηz (2)
X in formulap,zFor not by the signal matrix of noise pollution, ΨzFor load Stochastic Volatility, fluctuation range is set as ± 1%.ηzFor noise matrix, mzFor noise amplitude.
Step 2: determining each zone time window width T, and set sampling start time t1；
Step 3: synchronously obtaining each subregion sliding window matrix Xz.Parallel computing is relied on, by the calculating of each subregion Task distributes to the completion of its region computing terminal.Using sliding window technique, each zone time window width T obtained by step 2 With sampling start time t1, from the data source matrix X of each subregions,zSynchronously obtain the N of respective partitionz× T ties up sliding window square Battle array Xz, wherein NzIt is the dimension of sampled data, unit: a；T is sliding window width, unit: a；Matrix ranks ratio cz, such as formula (3) shown in,
And cz∈ (0, ∞) is ratio, no unit.
Step 4: each subregion sliding window matrix XzStandardization.In each subregion computing terminal, to sliding window matrix Xz It is normalized, obtains the non-Hermite Matrix of standard (Non-Hermitian Matrix) Xn,z, as shown in formula (4).
X in formulazi=(xzi,1,xzi,2,…,xzi,T)；μ(xzi)、σ(xzi) it is respectively xziMean value and standard deviation；μ (xn,zi)、σ(xn,zi) it is respectively non-Hermitian row matrix vector xn,ziMean value and standard deviation.
Step 5: synchronously calculating each subregion sample covariance matrix Sz.Seek the non-Hermite Matrix X of standardn,zSample association Variance matrix Sz, as shown in formula (5),
Subscript H indicates complex conjugate transposition in formula.
It is worth noting that, sample covariance matrix SzExperience spectral distribution function (Empirical Spectral Distribution, ESD) M-P rule is obeyed, such as formula (6) and (7):
Wherein
As c > 1, f in formulamp(x)Distribution function have quality 1-1/c in origin.A and b are respectively indicated in spectral density function Minimal eigenvalue and maximum eigenvalue, that is to say, that SzFeature Distribution value have a certain range, be bounded.When there is exception The randomness of system is destroyed when event occurs, and will be unsatisfactory for statistical law, and maximum eigenvalue can be more than its bounds, this hair Bright method is exactly the detection that this feature is utilized and carries out abnormality.
Step 6: each subregion sample covariance matrix maximum eigenvalue λmax,zEstimation.In each subregion computing terminal, make With each subregion sample covariance matrix S of power method iterative estimationzMaximum eigenvalue λmax,z, as shown in formula (8),
U in formulak,zIndicate maximum eigenvalue character pair vector in iterative process, vk,zIndicate " normalization " after iteration to Amount, v0,z={ 1,1 ..., 1 }Nz×1.When k is sufficiently big or | | max (uk,z)-max(uk-1,z) | | when < ε, max (uk,z)≈λmax,z。
Step 7: seeking each subregion MESCM dynamic threshold γz。
Step 7.1: the noise score using the classical spectrum estimate method after the correction of Kaiser window function to each subregion Do not estimated, by taking any one subregion in Fig. 2 as an example.
(1) assume there is N >=1 PMU, wherein the ac voltage signal of i-th (i=1 ... ..., N) a PMU is xi, signal number It is T according to length.Original signal x is filtered by formula (4)iIn DC component obtain xfi。
(2) Kaiser window can customize one group of adjustable window function, at that time domain representation such as formula (10).
In formula, I0(β) is first kind deformation zero Bessel function, the accuracy rate highest wherein estimated when β=38.Simultaneously Calculate window values wi(n) corresponding mean μ (wi), root-mean-square value wRmsi。
(3) by time-domain signal xfi(Discrete Time Fourier is converted by discrete time Fourier Transform, DTFT) it is transformed to frequency-region signal xi(ejω), power spectral density P is calculated in conjunction with Kaiser window valuesi(ejω), Such as formula (11).
Pi(ejω)=(1/T)/| xi(ejω)|2 (11)
(4) when calculating actual signal power and noise power, segment processing need to be carried out to frequency spectrum, calculates power spectral density pair Answer bandwidth such as formula (12).
bwi=wRmsi/(T×μ(wi)2) (12)
(5) every section of bandwidth corresponds to power density average amplitude such as formula (13).
(6) actual signal power P is obtained by Periodical Maximum Methodri, noise power PniSuch as formula (14).
(7) signal-to-noise ratio that current PMU signal can be obtained using the ratio of actual signal power and noise power, such as formula (15)。
(8) when noise amplitude is without acute variation, global signal-to-noise ratio is indicated with ρ, such as formula (16).
Step 7.2: applying the threshold gamma setting method based on Spiked model, be equally with any one subregion in Fig. 2 Example, the setting of threshold gamma such as formula (17),
In formula, ρ is the global signal-to-noise ratio (snr) estimation value of the current time subregion.C is that the dimension of the subregion holds ratio.α(0≤α≤1) For proportionality coefficient, it can be adjusted according to sliding window width T, generally take α=0.5.
Step 8: the extremely out-of-limit differentiation of electric network state.Judge whether maximum eigenvalue λmax,zGreater than threshold gammazIf λmax,z ≥γzIt sets up, then determines that power grid generating state is abnormal, give a warning.Otherwise, current state without exception, enables i=i+1, returns to step Rapid 3 continue to execute abnormal state testing process.
Wherein, the selection of sampled data type described in step 1 can be only bus voltage amplitude V.
Embodiment is as follows:
To verify validity and computational efficiency of the proposed method of the present invention under large scale electric network scene, by Matlab 3.0 tool software of R2014a and Power System Toolbox (PST) Version, with 420 machine of Poland, 2736 bus System example carries out dependence test as analog data source.Wherein, which contains 6 subregions, as shown in Figure 2.By Length limits, and the mother of each voltage class in the bus number and each subregion of the section 400kV interconnection and two sides is only provided in figure Line gauge mould.
Case test one: large scale electric network abnormal load identifies validity test
The load that the 110kV bus of number 2733 in the 1st subregion is arranged is abnormal variation, is specifically shown in the following table 1.
For noise jamming and random load fluctuation in simulation on-the-spot test, Gaussian noise is introduced in each area variable or signal Source, signal-to-noise ratio ρ=(30 ± 0.3) dB.In this way, the number of a simulation field measurement PMU signal can be constructed by formula (2) and formula (1) According to source matrix.
Since wider sliding window can make recognition result delay longer, and relatively narrow sliding window can make calculated result Inaccuracy, therefore comprehensively consider the area Hou Ge sliding window and take an identical value.In this way, each area's sliding window can be set by above-mentioned steps 2 Width T is 220, and proportionality coefficient α is 0.5, i.e. includes 220 groups of sampled datas in window, wherein 1 group for current data and 219 groups are historical data.Therefore the MESCM index change curve of the present invention mentioned method acquisition is from sampling instant t220Start.
1 load anomalous variation of table
Further, by the step 3,4 and 5, can obtain the sliding window matrix of each subregion respectively, standard it is non- Hermitian matrix and sample covariance matrix.Then, by step 6, the maximum that each subregion is quickly estimated using power method is special Value indicative, i.e. MESCM index, it is as shown in Figure 3 and Figure 4 in the change curve of the 220th to the 1000th sampled point.Meanwhile by described Step 7, the dynamic threshold that can get the MESCM index of each subregion of the whole network, in Fig. 3 and Fig. 4.It should be noted that due to upper It states under measurement condition, which only has 1st area and 4th area and crosses threshold value, therefore only provide the MESCM index of the 1st subregion and the 4th subregion Threshold value, as shown in Figure 3 and Figure 4.
By Fig. 3 and Fig. 4 as it can be seen that in t300Before sampling instant, even if having the fluctuation of random load and the interference of noise, But due to no exceptions, therefore each area MESCM index has no significant change.And in t301Sampling instant is numbered by the 1st subregion 2733 bus burdens with power are mutated by 100MW to the influence of 200MW, and whole system randomness is broken, so that in Fig. 3 The MESCM index of 1 subregion increases to 31.9 by 23.9 approximated steps, hence it is evident that has crossed the threshold value of its respective partition.And in t600 After sampling instant, influenced by the setting of aforementioned climbing type load, Fig. 3 is in t635Sampling instant MESCM index crosses threshold value and continuous Increase.These illustrate that this method also can not only effectively identify abnormal step change type load to abnormal climbing type load.
Fig. 4 is observed, it can be found that the MESCM index of the 4th subregion is in t in other regions in addition to the 1st subregion804Sampling Moment is more than threshold value, and only continue for 34 sampling instants, and other subregions have no significant change.On the one hand come from internal factor analysis It sees, this is mainly smaller related with the equivalent electrical distance in region 4 with disturbance point in region 1, or can be shown that subregion 1 and subregion 4 Electrical link is more closely compared to other subregions.On the other hand from the point of view of from outer as application, the threshold value presented by Fig. 3 and 4 is out-of-limit Degree difference can tentatively judge that load anomaly occurs in the 1st subregion.
Case tests the analysis of two: MESCM calculating error
Relatively traditional All Eigenvalues acquiring method, the present invention estimate thinking from partial feature value, utilize power method The quick estimation of MESCM index is carried out.Worst error is allowed to be 0.01 in iteration, the situation that maximum number of iterations is 1000 times Under, each subregion MESCM deviation in 220 to 1000 sampled point time windows under above-mentioned 420 machine, 2736 bus Poland electric system example Rate calculated result is as shown in Figure 5.
As it can be seen that the resulting maximum eigenvalue of the method for the present invention is asked compared to All Eigenvalues in MATLAB R2014a software Maximum eigenvalue obtained by function is solved, each subregion deviation ratio is respectively less than 5%, and practical the number of iterations is at 150 times or less.This shows Implement the quick estimation of MESCM using power method for large-scale electrical power system, accuracy is acceptable.
Case test three: time consuming analysis is calculated
Computer using 8 core CPU, 16GB a DRAM memory of outfit AMD-2700 3.2GHz has carried out one The side MSR of 54 machine of IEEE, 118 bus-bar system and 420 machine, 2736 bus Poland electric system under serial and parallel computation The calculating time consuming analysis of method and MESCM method.Wherein, 54 machine of IEEE, 118 bus-bar system wiring schematic diagram is as shown in Figure 6.Parallel Calculation is synchronously completed by the computing terminal of each one subregion of calculating core simulation of central processing unit, i.e. 3 subregions With 3 calculating cores, 6 subregions calculate the parallel computation of core simulation multi partition computing terminal with 6.In this way, for above-mentioned The system of two kinds of scales, MSR method of the tradition based on All Eigenvalues solution technique and this hair under serial and concurrent design conditions The bright calculating time-consuming for proposing the MESCM method based on power method is listed in table 2.
The calculating of tradition MSR method and the mentioned MESCM method of the present invention is time-consuming under the serial and concurrent design conditions of table 2
On the one hand, from the MSR method solved based on All Eigenvalues and the MESCM method based on partial feature value solution From the point of view of comparing angle, by the serial computing result of MSR and MESCM in table 2 it is found that the system for 118 bus scales, sampled point When number is 2300, the MESCM method based on power method calculates time-consuming 33.70s, and the MSR method of traditional calculations All Eigenvalues consumes When 63.19s, the former computational efficiency is compared with 2 times of the latter Gao Liaoyue.For the system of 2736 bus bus scales, sampling number is Time-consuming compared to the 4.18h of the latter when 1000, the former is only 124.45s, and computational efficiency improves nearly 121 times.As it can be seen that due to MSR method needs to calculate All Eigenvalues, and MESCM method only needs to calculate maximum eigenvalue, so, for extensive power train System, the MESCM method based on power method can significantly improve computational efficiency.
On the other hand, from serially compared with parallel efficiency calculation from the point of view of angle, further by table 2MSR and MESCM's and For row calculated result it is found that for 118 bus-bar systems, MESCM method parallel computation time-consuming is 7.58s, compared to serial computing, Speed-up ratio is about 4.4, has reached " superlinearity speed-up ratio ".And for 2736 bus-bar systems, MESCM method parallel computation time-consuming is 65.24s, since the 4th district system scale is total scale half, therefore compared to serial computing, speed-up ratio is about 2, has reached reason The acceleration effect of opinion.
In conclusion being proposed based on MESCM a kind of different suitable for large scale electric network using power method and parallel computing Normal load knows method for distinguishing.By Matlab software, pass through 54 machine of IEEE, 118 bus-bar system and 420 machine of Poland 2736 bus-bar system examples demonstrate the validity and rapidity of this method.Demonstrate the MESCM method based on power method compared to The MSR method of traditional calculations All Eigenvalues is obviously improved in computational efficiency.For large-scale electrical power system, using subregion Parallel computing can get ideal acceleration effect.
The above described is only a preferred embodiment of the present invention, being not intended to limit the present invention in any form, appoint What is to the above embodiments according to the technical essence of the invention any simply to repair without departing from technical solution of the present invention content Change, equivalent variations and modification, all of which are still within the scope of the technical scheme of the invention.

Claims (8)

1. a kind of large scale electric network abnormal load recognition methods based on power method and parallel computing, comprising the following steps:
Step 1: the data source square of each subregion is constructed using the whole network PMU metric data according to the region division of target power system Battle arrayX s,z ；On-site signal is for example simulated, noise and random fluctuation load can be introduced；
Step 2: using sliding window technique, determine each zone time window widthT, and set sampling start timet 1 ；
Step 3: relying on parallel computing, the distribution of computation tasks of each subregion is completed to its region computing terminal；According to institute If window widthT, by the data source matrix of each subregionX s,z Synchronously obtain respective partitionN z ×TTie up sliding window matrixX z ；
Step 4: in each subregion computing terminal, to each subregionX z Row vector be synchronously carried out standardization, it is non-to obtain standard Hermite Matrix (Non-Hermitian Matrix) X n,z ；
Step 5: in each subregion computing terminal, synchronously calculating each subregionX n,z The sample covariance matrix of matrixS z ；
Step 6: in each subregion computing terminal, synchronizing to each subregionS z Power method iterative calculation is carried out, estimates the sample of each subregion Covariance matrix maximum eigenvalue (MESCM), i.e.,S z Maximum eigenvalueλ max,z ；
Step 7: in these region computing terminals, each subregion estimates current time signal-to-noise ratio respectivelyρ z , obtain respective partition MESCM dynamic thresholdγ z ；
Step 8: the extremely out-of-limit differentiation of electric network state: judging whether the maximum eigenvalueλ max,z Greater than threshold valueγ z Ifλ max,z≥ γ z It sets up, then determines that power grid generating state is abnormal, give a warning；Otherwise, current state without exception, return step 3 continue to execute shape State abnormality detection process.
2. the large scale electric network abnormal load recognition methods based on power method and parallel computing as described in claim 1, It is characterized in that: the construction data source matrix in the step 1X s,z : have for onewThe large-scale electrical power system of a subregion, can The voltage magnitude or phase angle time series data of affiliated subregion dispatching control center main website are acquired by each subregion PMU and are uploaded to, respectively Construct respective partition data source matrixX s,z , wherein partition numberz =1,2,…,w；X s,z Expression it is as the formula (1),
(1)
It is described for example to simulate on-site signal, whereinX s,z There is only the interference of random noise, also receive load random fluctuation and make At influence, abnormality detection model can express as the formula,
(2)
In formulaX p,z For not by the signal matrix of noise pollution,Ψ z For load Stochastic Volatility, fluctuation range is set as ± 1%,η z For noise matrix,m z For noise amplitude.
3. the large scale electric network abnormal load recognition methods based on power method and parallel computing as described in claim 1, It is characterized in that: relying on parallel computing in the step 3, obtain respective partition sliding window matrixX z : use sliding window Technology, each zone time window width obtained by step 2TAnd sampling start timet 1, from the data source matrix of each subregionX s,z Together Step ground obtains respective partitionN z ×TTie up sliding window matrixX z , whereinN z It is the dimension of sampled data, unit: a；T is sliding Window width, unit: a；Matrix ranks ratio c z , as the formula (3),
(3)
c z ∈ (0, ∞) is ratio, no unit.
4. the large scale electric network abnormal load recognition methods based on power method and parallel computing as described in claim 1, It is characterized in that: in the step 4 in each subregion computing terminal, by each subregion sliding window matrixX z Standardization: to sliding window Mouth matrixX z It is standardized, obtains the non-Hermite Matrix of standard (Non- Hermitian Matrix) X n,z , such as formula (4) shown in,
(4)
In formulax zi=( x zi,1 ,x zi,2 ,…,x zi,T)；μ(x zi )、σ(x zi ) be respectivelyx zi Mean value and standard deviation；μ(x n,zi )、σ (x n,zi ) it is respectively non-Hermitian matrix row vectorx n,zi Mean value and standard deviation.
5. the large scale electric network abnormal load recognition methods based on power method and parallel computing as described in claim 1, It is characterized in that: in the step 5 in each subregion computing terminal, seeking the non-Hermite Matrix of each scoping rulesX n,z Sample association side Poor matrix S z , as the formula (5),
(5)
Subscript H indicates complex conjugate transposition in formula.
6. the large scale electric network abnormal load recognition methods according to claim 1 based on power method and parallel computing, It is characterized by: in the step 6 in each subregion computing terminal, using each subregion sample covariance matrix of power method iterative estimation Maximum eigenvalueλ max,z , as the formula (6),
(6)
In formulau k,z Indicate maximum eigenvalue character pair vector in iterative process,v k,z Indicate iterative vectorized after " normalization ",v 0,z ={ 1,1 ..., 1 } Nz×1, whenkIt is sufficiently big or | | max (u k,z )-max(u k-1,z ) | | when < ε, max (u k,z )≈λ max,z 。
7. the large scale electric network abnormal load recognition methods according to claim 1 based on power method and parallel computing, It is characterized by: seeking each subregion MESCM dynamic threshold described in step 7γ z : each subregion estimates current time respectively first Then global signal-to-noise ratio goes out the MESCM dynamic threshold of corresponding subregion according to the subregion signal-to-noise ratio computationγ z , specific steps include:
Step 7.1: using the classical spectrum estimate method after the correction of Kaiser window function to the signal-to-noise ratio of each subregion respectively into Row estimation, by taking any one subregion as an example:
(1) assume haveN >=1 PMU, wherein thei(i=1, ……,N) ac voltage signal of a PMU isx i , signal data is long Degree isT；Original signal is filtered by formula (4)x i In DC component obtainx fi ；
(7)
(2) Kaiser window can customize one group of adjustable window function, at that time domain representation such as formula (8)；
(8)
In formula,I 0(β) it is first kind deformation zero Bessel function, whereinβ=The accuracy rate highest estimated when 38；Window is calculated simultaneously Valuew i (n) corresponding mean valueμ(w i ), root-mean-square valuew Rmsi ；
(3) by time-domain signalx fi By discrete time Fourier transformation (Discrete Time Fourier Transform, DTFT) it is transformed to frequency-region signalx i (e jω ), power spectral density is calculated in conjunction with Kaiser window valuesP i (e jω ), such as formula (9),
(9)
(4) when calculating actual signal power and noise power, segment processing need to be carried out to frequency spectrum, calculates power spectral density and corresponds to band Width such as formula (10),
(10)
(5) every section of bandwidth corresponds to power density average amplitude such as formula (11),
(11)
(6) actual signal power is obtained by Periodical Maximum MethodP ri , noise power PniSuch as formula (12),
(12)
(7) signal-to-noise ratio that current PMU signal can be obtained using the ratio of actual signal power and noise power, such as formula (13),
(13)
(8) it when noise amplitude is without acute variation, usesρIndicate global signal-to-noise ratio, such as formula (14),
(14)
Step 7.2: applying the threshold value based on Spiked modelγSetting method sets MESCM dynamic threshold, equally with Wherein for any one subregion, threshold valueγSetting such as formula (15),
(15)
In formula,ρFor current time, global signal-to-noise ratio (snr) estimation value,cHold ratio for dimension,α(0≤αIt≤1) is proportionality coefficient, it can be according to cunning Dynamic window widthTIt is adjusted, generally takesα=0.5。
8. the large scale electric network abnormal load based on power method and parallel computing as described in any one of claims 1 to 7 Recognition methods, it is characterised in that: the type of sampled data can be only bus voltage amplitude in the step 1V。
CN201811556503.1A 2018-12-19 2018-12-19 Large-scale power grid abnormal load identification method based on power method and parallel computing technology Active CN109657613B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811556503.1A CN109657613B (en) 2018-12-19 2018-12-19 Large-scale power grid abnormal load identification method based on power method and parallel computing technology

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811556503.1A CN109657613B (en) 2018-12-19 2018-12-19 Large-scale power grid abnormal load identification method based on power method and parallel computing technology

Publications (2)

Publication Number Publication Date
CN109657613A true CN109657613A (en) 2019-04-19
CN109657613B CN109657613B (en) 2023-05-02

Family

ID=66115805

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811556503.1A Active CN109657613B (en) 2018-12-19 2018-12-19 Large-scale power grid abnormal load identification method based on power method and parallel computing technology

Country Status (1)

Country Link
CN (1) CN109657613B (en)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110137944A (en) * 2019-04-24 2019-08-16 国网山东省电力公司莱芜供电公司 A kind of voltage stability disturbance source locating method based on Random Matrices Theory
CN112230199A (en) * 2019-07-15 2021-01-15 天津大学 Laser radar echo blind denoising method based on high-dimensional characteristic value analysis
CN112307003A (en) * 2020-11-02 2021-02-02 合肥优尔电子科技有限公司 Power grid data multidimensional auxiliary analysis method, system, terminal and readable storage medium
CN112904311A (en) * 2021-01-27 2021-06-04 安徽东风机电科技股份有限公司 Low-cost laser detector echo signal peak value keeps and acquisition circuit
CN113537844A (en) * 2021-09-15 2021-10-22 山东大学 Method and system for analyzing load behaviors of regional energy Internet based on random matrix
CN114062850A (en) * 2021-11-17 2022-02-18 江南大学 Double-threshold power grid early fault detection method
WO2023241326A1 (en) * 2022-06-14 2023-12-21 无锡隆玛科技股份有限公司 Power grid anomaly detection method based on maximum eigenvalue rate of sample covariance matrix
WO2023241327A1 (en) * 2022-06-14 2023-12-21 无锡隆玛科技股份有限公司 Power grid anomaly locating method based on maximum eigenvector
CN117332215A (en) * 2023-12-01 2024-01-02 深圳市大易电气实业有限公司 High-low voltage power distribution cabinet abnormal fault information remote monitoring system

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101685479A (en) * 2008-09-27 2010-03-31 国家电力调度通信中心 Power grid online comprehensive pre-warning method and system based on massively parallel processing
CN104200059A (en) * 2014-07-24 2014-12-10 温州大学 Oil-water well behavior analysis and predication device and method
CN104316768A (en) * 2014-10-28 2015-01-28 国家电网公司 Negative sequence impedance parameter estimation method for locating three-phase unbalanced disturbance source
US20150127819A1 (en) * 2013-11-01 2015-05-07 The Nielsen Company (Us), Llc Methods and apparatus to credit background applications
JP2016009406A (en) * 2014-06-25 2016-01-18 一般財団法人電力中央研究所 Instantaneous value analysis system of electric power system using sparse tableau approach
CN105354656A (en) * 2015-10-09 2016-02-24 珠海许继芝电网自动化有限公司 Partition decoupling based distributed parallel computing method and system for distribution network state estimation
CN108196165A (en) * 2018-01-09 2018-06-22 贵州大学 Power grid abnormal state detection method based on sample covariance matrix maximum eigenvalue

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101685479A (en) * 2008-09-27 2010-03-31 国家电力调度通信中心 Power grid online comprehensive pre-warning method and system based on massively parallel processing
US20150127819A1 (en) * 2013-11-01 2015-05-07 The Nielsen Company (Us), Llc Methods and apparatus to credit background applications
JP2016009406A (en) * 2014-06-25 2016-01-18 一般財団法人電力中央研究所 Instantaneous value analysis system of electric power system using sparse tableau approach
CN104200059A (en) * 2014-07-24 2014-12-10 温州大学 Oil-water well behavior analysis and predication device and method
CN104316768A (en) * 2014-10-28 2015-01-28 国家电网公司 Negative sequence impedance parameter estimation method for locating three-phase unbalanced disturbance source
CN105354656A (en) * 2015-10-09 2016-02-24 珠海许继芝电网自动化有限公司 Partition decoupling based distributed parallel computing method and system for distribution network state estimation
CN108196165A (en) * 2018-01-09 2018-06-22 贵州大学 Power grid abnormal state detection method based on sample covariance matrix maximum eigenvalue

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
刘俊，郝旭东，王旭等: "" 电力系统并行计算综述"", 《智慧电力》 *
周忠强，韩松，李洪乾: ""基于Spiked模型的低信噪比环境电网异常状态检测"", 《电测与仪表》 *
裴健，陈晨: ""电力系统等值阻抗的在线估计法"", 《华北电力大学学报（自然科学版）》 *
贾金伟，吴旭鹏，李启本: ""基于并行计算的大数据挖掘在电网中的应用"", 《电力与能源》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110137944A (en) * 2019-04-24 2019-08-16 国网山东省电力公司莱芜供电公司 A kind of voltage stability disturbance source locating method based on Random Matrices Theory
CN112230199A (en) * 2019-07-15 2021-01-15 天津大学 Laser radar echo blind denoising method based on high-dimensional characteristic value analysis
CN112230199B (en) * 2019-07-15 2022-10-25 天津大学 Laser radar echo blind denoising method based on high-dimensional characteristic value analysis
CN112307003A (en) * 2020-11-02 2021-02-02 合肥优尔电子科技有限公司 Power grid data multidimensional auxiliary analysis method, system, terminal and readable storage medium
CN112904311A (en) * 2021-01-27 2021-06-04 安徽东风机电科技股份有限公司 Low-cost laser detector echo signal peak value keeps and acquisition circuit
CN113537844A (en) * 2021-09-15 2021-10-22 山东大学 Method and system for analyzing load behaviors of regional energy Internet based on random matrix
CN114062850A (en) * 2021-11-17 2022-02-18 江南大学 Double-threshold power grid early fault detection method
CN114062850B (en) * 2021-11-17 2022-09-23 江南大学 Double-threshold power grid early fault detection method
WO2023241326A1 (en) * 2022-06-14 2023-12-21 无锡隆玛科技股份有限公司 Power grid anomaly detection method based on maximum eigenvalue rate of sample covariance matrix
WO2023241327A1 (en) * 2022-06-14 2023-12-21 无锡隆玛科技股份有限公司 Power grid anomaly locating method based on maximum eigenvector
CN117332215A (en) * 2023-12-01 2024-01-02 深圳市大易电气实业有限公司 High-low voltage power distribution cabinet abnormal fault information remote monitoring system
CN117332215B (en) * 2023-12-01 2024-03-15 深圳市大易电气实业有限公司 High-low voltage power distribution cabinet abnormal fault information remote monitoring system

Also Published As

Publication number Publication date
CN109657613B (en) 2023-05-02

Similar Documents

Publication Publication Date Title
CN109657613A (en) Large scale electric network abnormal load recognition methods based on power method and parallel computing
CN106055918B (en) Method for identifying and correcting load data of power system
CN110417011B (en) Online dynamic security assessment method based on mutual information and iterative random forest
CN110162871B (en) Power system dynamic estimation method based on unscented Kalman particle filter
CN103678941A (en) Prediction method for electrode air gap breakdown voltage
CN110059845B (en) Metering device clock error trend prediction method based on time sequence evolution gene model
CN107104442A (en) The computational methods of Probabilistic Load containing wind power plant of meter and parameter fuzzy
CN109659957B (en) APIT-MEMD-based power system low-frequency oscillation mode identification method
CN112307677A (en) Power grid oscillation mode evaluation and safety active early warning method based on deep learning
CN110632455A (en) Fault detection and positioning method based on distribution network synchronous measurement big data
CN110765703A (en) Wind power plant aggregation characteristic modeling method
Xiyun et al. Wind power probability interval prediction based on bootstrap quantile regression method
CN104112062A (en) Method for obtaining wind resource distribution based on interpolation method
CN112688324B (en) Power system low-frequency oscillation mode identification method based on FastICA and TLS-ESPRIT
Yuan et al. A network partition approach for distributed three-phase state estimation of active distribution networks
CN113659564A (en) Low-voltage distribution network topology identification method and system based on voltage fluctuation feature clustering
CN107958120A (en) A kind of system Thevenin&#39;s equivalence calculation method of parameters based on power series expansion
CN111175608A (en) Power distribution network harmonic responsibility quantitative division method based on accelerated independent component analysis
CN105975710B (en) The detection of bad data collection and recognition methods for the identification of synchronous generator on-line parameter
CN114840591A (en) Method and device for determining sectional switch power data
CN114970698A (en) Metering equipment operation performance prediction method based on improved LWPLS
CN111368392B (en) Single-sample non-stationary wind speed simulation method based on MEMD and SRM
Liu et al. A modified PSO algorithm for parameters identification of the double-dispersion Cole Model
CN103258144B (en) Online static load modeling method based on data of fault recorder
CN109390957A (en) A kind of wind power fluctuation induces the detection method of systems force oscillation

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