CN113675951A - Power system dynamic state estimation method based on PMU (phasor measurement Unit) and SCADA (Supervisory control and data acquisition) mixed measurement - Google Patents

Power system dynamic state estimation method based on PMU (phasor measurement Unit) and SCADA (Supervisory control and data acquisition) mixed measurement Download PDF

Info

Publication number
CN113675951A
CN113675951A CN202110975112.9A CN202110975112A CN113675951A CN 113675951 A CN113675951 A CN 113675951A CN 202110975112 A CN202110975112 A CN 202110975112A CN 113675951 A CN113675951 A CN 113675951A
Authority
CN
China
Prior art keywords
measurement
pmu
state
rtu
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
Application number
CN202110975112.9A
Other languages
Chinese (zh)
Other versions
CN113675951B (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.)
Nantong University
Original Assignee
Nantong 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 Nantong University filed Critical Nantong University
Priority to CN202110975112.9A priority Critical patent/CN113675951B/en
Publication of CN113675951A publication Critical patent/CN113675951A/en
Application granted granted Critical
Publication of CN113675951B publication Critical patent/CN113675951B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J13/00—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network
    • H02J13/00001—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network characterised by the display of information or by user interaction, e.g. supervisory control and data acquisition systems [SCADA] or graphical user interfaces [GUI]
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J13/00—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network
    • H02J13/00002—Circuit arrangements for providing remote indication of network conditions, e.g. an instantaneous record of the open or closed condition of each circuitbreaker in the network; Circuit arrangements for providing remote control of switching means in a power distribution network, e.g. switching in and out of current consumers by using a pulse code signal carried by the network characterised by monitoring
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00—Circuit arrangements for ac mains or ac distribution networks
    • H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
    • H02J3/242—Arrangements for preventing or reducing oscillations of power in networks using phasor measuring units [PMU]
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
    • H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
    • 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
    • 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/40—Display of information, e.g. of data or controls

Landscapes

  • Engineering & Computer Science (AREA)
  • Power Engineering (AREA)
  • Human Computer Interaction (AREA)
  • Supply And Distribution Of Alternating Current (AREA)

Abstract

The invention provides a dynamic state estimation method of a power system based on PMU and SCADA mixed measurement, and belongs to the technical field of state monitoring and control of the power system. The technical scheme is as follows: the method comprises the steps of firstly reading real-time power grid data through real-time network parameters and topological structures of the power system and dynamic changes of load output, secondly constructing a power system dynamic model, an RTU (remote terminal Unit) measurement equation, a PMU (phasor measurement Unit) measurement equation and a mixed measurement expression under the normal operation of the power system, then performing multisource asynchronous dynamic state estimation on the power system by combining an extended Kalman filtering basis after drawing the measurement equation and configuring the PMU, and finally judging an estimation cut-off condition. The invention has the beneficial effects that: the invention overcomes the problem of inconsistent acquired information frequency under hybrid measurement, and the accuracy of system state estimation can be improved by the hybrid measurement and the introduction of PMU, so that a real-time and better estimation algorithm is obtained, and the invention has the characteristics of safer, more accurate and real-time power control.

Description

Power system dynamic state estimation method based on PMU (phasor measurement Unit) and SCADA (Supervisory control and data acquisition) mixed measurement
Technical Field
The invention relates to the technical field of power system state monitoring and control, in particular to a power system dynamic state estimation method based on PMU and SCADA mixed measurement.
Background
As the most important key infrastructure of the country and society, the operation safety and reliability of the power system are the primary concerns. With the rapid development of the scientific and technological era, the power grid structure is increasingly complex, the operation is increasingly heavy, the power grid problem is easily caused, the power grid electric energy loss is increased, and the power grid disaster can be caused under severe conditions. Once the problem occurs, the power grid problem not only affects the normal life of people, but also even brings obstruction to the economic development of the society.
The traditional method for estimating the static state of the power system only determines the estimation of the state quantity at a certain moment according to measured data at the moment, cannot meet the requirements on real-time economic dispatching and preventive control of the power system, and the early basic weighted least square estimation algorithm for estimating the state of the power system is common.
In recent years, however, a global positioning system based PMU is gradually applied to a power system, and the PMU can measure a node voltage and provide a phase angle measurement of a related branch current, and the measurement is linear with a system state. Compared with the traditional nonlinear state estimation method, the linear state estimation method based on PMU measurement has higher precision and faster calculation speed. Therefore, for the complex power system, it has important practical significance and application value to develop the power system dynamic state estimation research under the hybrid measurement.
With the development of the dynamic state estimation technology of the power system, a power system dynamic state estimation method based on distributed extended kalman filtering under PMU/SCADA hybrid measurement is proposed to solve the above problems.
Disclosure of Invention
The invention aims to provide a dynamic state estimation method of a power system based on PMU and SCADA mixed measurement, which can accurately measure the state of the power system, thereby providing necessary basis for the stable operation of the power system and better monitoring the operation state (voltage amplitude and phase angle) of each node of the system.
The invention is realized by the following measures: the method for estimating the dynamic state of the power system based on PMU and SCADA mixed measurement comprises the following steps:
a. describing a dynamic state equation of the power system through a quasi-steady state model;
b. the RTU measurement equation expressed by the sampling rate is constructed by combining RTUs with different sampling rates, PMU measurement devices, active power among all node lines, node active power injection and voltage amplitude values, and then according to the configured PMU, the PMU measurement equation expressed by the sampling rate is constructed and a mixed measurement equation of the RTU and the PMU is constructed; (ii) a
c. And carrying out dynamic state estimation on the power system on the basis of extended Kalman filtering by using a multi-source asynchronous mixed measurement equation.
Further, in the step a, a quasi-steady state model of the power system is described under a standard S node system study (the voltage phase angle of the default system node 1 is a reference phase angle, and thus is not taken as a state variable):
xk+1＝Fkxk+uk+wk
wherein ,xkIs (2S-1) x 1 dimensional state quantity at the moment k, namely the voltage phase angle and the voltage amplitude of each node of the system; fk(. is an adaptive state transfer function; u. ofkIs (2S-1) x 1 dimension input parameter item; w is akIs (2S-1). times.1 dimensional process noise and follows a zero mean normal distribution, i.e., w to N (0, Q), with Q being a ((2S-1). times.1). times. ((2S-1). times.1) dimensional system error covariance matrix.
Further, the quantifying measurement equation of the RTU in the step b mainly includes constructing a measurement vector, and quantifying the measurement expression by combining the sampling frequency of the RTU:
constructing a measurement vector: the active power between lines and the active power of the nodes are expressed as follows:
Pij(θij,Uij)＝Ui 2(gsi+gij)-UiUj(gijcosθij+bijsinθij)
Figure BDA0003227379120000021
wherein ,θijIs the voltage phase angle between bus i and j; u shapei、UjThe voltage amplitudes of the buses i and j, respectively; gij、bijRespectively, the conductance and susceptance of line ij; gsiIs the conductance of the shunt leg ij; gij、BijRespectively the real and imaginary parts of the elements ij in the admittance matrix.
Combined node voltage amplitude Ui(Ui) Constructing a measurement vector of the RTU:
Figure BDA0003227379120000022
linearizing the research object, omitting terms of second order and above by using Taylor expansion, and obtaining a linearized RTU measurement vector:
Figure BDA0003227379120000023
and quantifying the measurement expression by combining the RTU sampling frequency: considering the existence of measurement errors, the RTU measurement equation is:
zR(k)＝HRxR(k)+vR(k)
wherein ,zR(k) Is (m)RX 1) measurement vector of dimension RTU; hRIs (m)RState quantity coefficient of dimension x ((2S-1) × 1))A matrix; x is the number ofR(k) The state quantity is (2S-1) multiplied by 1 dimension when only the RTU measuring device is installed; v. ofR(k) Is (m)RX 1) dimensional measurement of error vector and obeying a zero-mean normal distribution, i.e. vR～N(0,RR)，RRIs (m)R×mR) Measuring an error variance matrix in a dimension mode; m isRMeasuring the number of RTUs during independent installation; x is the number of0Is a linearization point.
Suppose the whole sampling period is h time sections, zR(k) Has a sampling period of TR＝LRh；zP(k) Has a sampling period of TP＝LPh, let M be LR and LPLeast common multiple of ql＝M/LlAnd l is equal to R, and P is the sampling number of the RTU and the PMU in a common period respectively, and the RTU measurement equation is described as follows:
ZR(ks)＝HRxR(ks)+vR(ks)
wherein ,ks＝LRk represents a slow rate time index.
Furthermore, quantifying the measurement equation of the PMU in step b mainly includes constructing a measurement vector, and quantifying the measurement expression by combining the sampling frequency of the PMU:
constructing a measurement vector: under the research of standard S node system HpFrom the element aijFormed and at a node
Figure BDA0003227379120000031
The PMU is installed, and the measurement vector is constructed as follows:
Figure BDA0003227379120000032
Figure BDA0003227379120000033
in the formula ,zP(k) Is a measurement vector of (2 lx 1) -dimensional PMU, HPA state quantity coefficient matrix of (2 lx((2S-1) × 1)) dimension, xP(k) The state quantity (including voltage amplitude and phase angle) v is the (2S-1) × 1 dimension when PMU measurement device is installed independentlyP(k) Is a (2 lx 1) -dimensional measurement error vector and follows a normal distribution of zero mean, i.e., vP～N(0,RP)，RPIs a (2l × 2l) dimensional measurement error variance matrix. Note: i is 2l and by default the voltage phase angle of node 1 is taken as the reference phase angle, where the voltage phase angle is 0.
And quantifying the measurement expression by combining the sampling frequency of the PMU: z is a radical ofP(k) Showing the PMU measurements at time k. Assuming that l PMUs are installed, in
Figure BDA0003227379120000034
And the node is installed, and is mainly used for measuring the amplitude and the phase angle of the bus voltage and the amplitude and the phase angle of branch current at a system node where the measuring device is installed. Taking into account the presence of measurement noise vPThe PMU measurement equation is easily obtained:
zP(ks')＝HPxP(ks')+vP(ks')
wherein ,ks'＝LPk。
Further, in the step b, a mixing measurement is constructed: at the moment of hybrid state estimation, the SCADA data and the PMU data are updated, the hybrid measurement is represented by the measurement equations of the RTU and the PMU together, and the PMU is used for every LPSampling once every moment, RTU every LRSampling is carried out once at a moment, the PMU and the SCADA/RTU can acquire measured values at every M moments, and the expression of the mixed measurement is described as follows:
Figure BDA0003227379120000035
in the formula ,
Figure BDA0003227379120000036
Zmix(k) is ((m)RA +2l) x 1) -dimensional mixture measurement vector, x (k) is a state quantity under (n x 1) -dimensional mixture measurement, vmix(k) Is ((m)RMixed measurement error vector of +2l) x 1) dimension and obeying zeroNormal distribution of the mean, i.e. vmix～N(0,Rm)，RmIs ((m)R+2l)×(mR+2l)) dimension measure error variance matrix.
In addition, when only a fast sampling rate measuring device PMU is used for data updating and only a slow sampling rate measuring device RTU is used for data updating, the multisource asynchronous measurement is easily expressed as follows:
Figure BDA0003227379120000037
further, in the step c, on the premise of obtaining the RTU, PMU measurement equation and hybrid measurement equation, the dynamic state estimation of the power system is performed by using the multisource asynchronous hybrid measurement equation in combination with the extended kalman filter, including parameter identification, state prediction and state filtering:
parameter identification: the state transition F of the system in step akAnd an entry ukIt is difficult to quantify the specific mathematical form, and a Holt's two-parameter exponential smoothing model is proposed to identify the state transition quantity FkAnd input parameter item uk。
And (3) state prediction: state of system under standard S node system research
Figure BDA0003227379120000041
Sum covariance matrix Pk+1|kThe prediction result is as follows:
Figure BDA0003227379120000042
Pk+1|k＝FkPk|kFk T+Qk
wherein ,
Figure BDA0003227379120000043
the predicted value of the state at the time k +1 by (2S-1) × 1 dimension,
Figure BDA0003227379120000044
updating the correction value for the state of time k versus time k, Pk+1|kThe matrix is predicted for the error covariance at time k versus time k + 1.
And (3) state filtering: using the measurements at the multi-rates in step c, real-time measurements Z for a group of power systems have been obtainedkFiltering to obtain new state estimation vector by calculating Kalman gain and corresponding innovation value
Figure BDA0003227379120000045
The state filtered value is:
Kk+1＝Pk+1|kHθ T[HθPk+1|kHθ T+Rθ]-1
Figure BDA0003227379120000046
Figure BDA0003227379120000047
Pk+1|k+1＝[In-Kk+1Hθ]Pk+1|k
where θ ∈ { R, P, mix }.
Compared with the prior art, the invention has the beneficial effects that:
1) the PMU is introduced into the measurement device, so that real-time and accurate voltage amplitude value and phase angle equivalent measurement values can be provided for the system, the measurement redundancy of the system is increased by adopting a measurement algorithm of mixing the RTU and the PMU, the PMU sampling interval is short, and the state estimation precision of the power system is improved.
2) The invention constructs a step-by-step measurement expression by a time-interval analysis method aiming at two measurement devices with different sampling rates, namely an RTU (remote terminal Unit) and a PMU (Power management Unit), and easily realizes the mixed measurement sampling condition of a multi-rate sampling device by using a simple mathematical expression. The method can accurately measure the state of the power system, thereby further promoting the development of power companies.
3) In the prior art, node voltage amplitude and phase angle data accurately measured by a PMU are directly involved in state estimation, and the problem of asynchronous sampling frequency of the PMU measuring device and the RTU measuring device is not considered.
4) The invention establishes the dynamic state estimation model of the power system from the perspective of solving the asynchronous problem of RTU and PMU mixed measurement sampling, and better conforms to the actual operation condition of reading data by the power measurement device.
Drawings
The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the principles of the invention and not to limit the invention.
FIG. 1 is an overall flow chart of the present invention.
FIG. 2 is a diagram of the pi-type equivalent circuit measurement calculation employed in the present invention.
FIG. 3 is a diagram illustrating a fusion process of RTU and PMU mixed measurement data according to the present invention.
FIG. 4 is an IEEE9 test chart in an embodiment of the invention.
FIG. 5 is a graph of the voltage amplitude estimation at node 7 of the present invention.
FIG. 6 is a graph of average error percentage of the voltage amplitude of the 7-node according to the present invention.
FIG. 7 is an overall performance index plot for the algorithm of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. Of course, the specific embodiments described herein are merely illustrative of the invention and are not intended to be limiting.
Example 1
Referring to fig. 1 to 7, the present invention provides a technical solution that the present invention provides a method for estimating a dynamic state of a power system based on PMU and SCADA hybrid measurement, and a flow is shown in fig. 1, and includes the following steps:
(1) reading grid data
In this step, real-time grid data is read through real-time network parameters and topology of the power system and dynamic changes of load output.
(2) The system comprises a power system dynamic model, an RTU measurement equation, a PMU measurement equation and a mixed measurement expression, wherein under the condition of normal operation of the power system, the state quantity of the power system changes along with the change of a load, and the state quantity x mainly comprises a voltage amplitude and a phase angle. In this step, the power system is described as follows using a quasi-steady state model (assuming that the reference phase angle is the voltage phase angle of node 1 and is not used as a state variable):
xk+1＝Fkxk+uk+wk
wherein ,xkIs the (n x 1) dimensional state quantity at the time k, namely the voltage phase angle and the voltage amplitude of each node of the system; fk(. h) is an (n × 1) -dimensional state transfer function; u. ofkIs an (n × 1) dimensional input parameter item; w is akIs (N × 1) -dimensional process noise and follows a zero-mean normal distribution, i.e., w to N (0, Q), Q being a (N × N) -dimensional system error covariance matrix.
Preferably, the measurement vector of the RTU consists of the active power between the lines of the nodes, the active power of the nodes and the node voltage amplitude, as illustrated below by the measurement function used for the correlation measurement of the typical pi-type equivalent circuit of fig. 2, gij、bijRespectively, the conductance and susceptance of line ij; gsiIs the conductance of the shunt leg ij; gij、BijRespectively the real and imaginary parts of the elements ij in the admittance matrix. The active power between lines and the active power equation of the node are expressed as:
Pij(θij,Uij)＝Ui 2(gsi+gij)-UiUj(gijcosθij+bijsinθij)
Figure BDA0003227379120000051
combined node voltage amplitude Ui(Ui) Constructing a measurement vector of the RTU:
Figure BDA0003227379120000052
linearizing the research object, omitting terms of second order and above by using Taylor expansion to obtain a linearized RTU measurement vector, considering the existence of measurement errors, and the RTU measurement equation is as follows:
Figure BDA0003227379120000061
zR(k)＝HRxR(k)+vR(k)
wherein ,zR(k) Is (m)RX 1) measurement vector of dimension RTU; hRIs (m)RA state quantity coefficient matrix of dimension x ((2S-1) × 1)); x is the number ofR(k) The state quantity is (2S-1) multiplied by 1 dimension when only the RTU measuring device is installed; v. ofR(k) Is (m)RX 1) dimensional measurement of error vector and obeying a zero-mean normal distribution, i.e. vR～N(0,RR)，RRIs (m)R×mR) Measuring an error variance matrix in a dimension mode; m isRMeasuring the number of RTUs during independent installation; x is the number of0Is a linearization point.
Preferably, the general RTU measurement equation is written as follows, considering the different sampling frequencies of the different measurement devices:
ZR(ks)＝HRxR(ks)+vR(ks)
wherein ,ks＝LRk represents a slow rate time index, the whole sampling period is h time sections, zR(k) Has a sampling period of TR＝LRh；zP(k) Has a sampling period of TP＝LPh, let M be LR and LPLeast common multiple of ql＝M/LlAnd l is R, and P is the sampling number of the RTU and the PMU in a common period respectively.
Preferably, buildingPMU measurement equation under standard S node system research HpFrom the element aijFormed and at a node
Figure BDA0003227379120000065
The PMU is installed, and the measurement vector is constructed as follows:
Figure BDA0003227379120000062
Figure BDA0003227379120000063
in the formula ,zP(k) Is a measurement vector of (2 lx 1) -dimensional PMU, HPA state quantity coefficient matrix of (2 lx((2S-1) × 1)) dimension, xP(k) The state quantity (including voltage amplitude and phase angle) v is the (2S-1) × 1 dimension when PMU measurement device is installed independentlyP(k) Is a (2 lx 1) -dimensional measurement error vector and follows a normal distribution of zero mean, i.e., vP～N(0,RP)，RPIs a (2l × 2l) dimensional measurement error variance matrix. Note: i is 2l and by default the voltage phase angle of node 1 is taken as the reference phase angle, where the voltage phase angle is 0.
And quantifying the measurement expression by combining the sampling frequency of the PMU: z is a radical ofP(k) Showing the PMU measurements at time k. Assuming that l PMUs are installed, in
Figure BDA0003227379120000066
And the node is installed, and is mainly used for measuring the amplitude and the phase angle of the bus voltage and the amplitude and the phase angle of branch current at a system node where the measuring device is installed. Taking into account the presence of measurement noise vPThe PMU measurement equation is easily obtained:
zP(ks')＝HPxP(ks')+vP(ks')
wherein ,ks'＝LPk。
Preferably, a mixing measurement expression is constructed, and SCADA data and PMU number are calculated at the moment of estimating the mixing stateAccording to the updating, the mixed measurement is represented by the measurement equations of RTU and PMU, and PMU is used for every LPSampling once every moment, RTU every LRSampling is carried out once at a moment, the PMU and the SCADA/RTU can acquire measured values at every M moments, and the expression of the mixed measurement is described as follows:
Figure BDA0003227379120000064
in the formula ,
Figure BDA0003227379120000071
Zmix(k) is ((m)RA +2l) x 1) -dimensional mixture measurement vector, x (k) is a state quantity under (n x 1) -dimensional mixture measurement, vmix(k) Is ((m)R+2l) x 1) dimensional mixed measurement error vector and obeys a zero-mean normal distribution, i.e., vmix～N(0,Rm)，RmIs ((m)R+2l)×(mR+2l)) dimension measure error variance matrix.
Besides, when only data updating of a fast sampling rate measuring device PMU and only data updating of a slow sampling rate measuring device RTU are considered, it is easy to show that multi-source asynchronous measurement is as follows:
Figure BDA0003227379120000072
(3) power system dynamic state estimation
After a measurement equation is described and PMU is configured, multisource asynchronous dynamic state estimation is carried out on the power system on the basis of extended Kalman filtering.
Amount of state transition F of systemkAnd an entry ukIt is difficult to quantify the specific mathematical form, and a Holt's two-parameter exponential smoothing model is proposed to identify the state transition quantity FkAnd input parameter item uk. State of system under standard S node system research
Figure BDA0003227379120000073
Hexie (harmony partner)Difference matrix Pk+1|kThe prediction result is as follows:
Figure BDA0003227379120000074
Pk+1k＝FkPkkFk T+Qk
wherein ,
Figure BDA00032273791200000710
the predicted value of the state at the time k +1 by (2S-1) × 1 dimension,
Figure BDA0003227379120000076
updating the correction value for the state of time k versus time k, Pk+1|kThe matrix is predicted for the error covariance at time k versus time k + 1.
Using measurements at multiple rates, real-time measurements Z have been obtained for a set of power systemskFiltering to obtain new state estimation vector by calculating Kalman gain and corresponding innovation value
Figure BDA0003227379120000077
The state filtered value is:
Kk+1＝Pk+1|kHθ T[HθPk+1|kHθ T+Rθ]-1
Figure BDA0003227379120000078
Figure BDA0003227379120000079
Pk+1|k+1＝[In-Kk+1Hθ]Pk+1|k
and updating different measurement values corresponding to different sampling moments, wherein theta is belonged to { R, P, mix }.
(4) Estimation cut-off determination
And if the state of the power system is estimated in N time sections, when k is larger than N, finishing simulation, and outputting the voltage amplitude and the phase angle of each node, otherwise, turning to the step 3.
In the MatlabR2018a environment, the method designed by the invention is verified by taking the standard 9 node testing system shown in fig. 4 as an example. The node load output curve is composed of linear trend and random disturbance, the disturbance is zero-mean normal distribution Gaussian white noise, the variance is 0.01, the true value of the node state dynamically changes along with the real-time change of the node load output value, and the true value of the state is obtained by combining Matpower in Matlab. The measurement vector consists of 27 observed values, which are the active power P between branches1-4,P4-5,P5-6,P3-6,P6-7,P7-8,P8-2,P8-9,P9-4(ii) a Active power P between nodes1,P2,P3,P4,P5,P6,P7,P8,P9(ii) a Amplitude of voltage V1,V2,V3,V4,V5,V6,V7,V8,V9. The nodes 4 and 8 are configured with PMUs, assuming that standard deviations of voltage and current phasors of the PMUs are both 0.001, a standard deviation of measurement data of the SCADA system is 0.01, and a standard deviation of the measurement data is 0.01, and in addition, active injection measurement of each node of the system, active power injection between branches of each node, and voltage amplitude measurement can be calculated by using a related measurement function.
First, Fk、ukThe model can be updated in real time according to a Holt's two-parameter exponential smoothing model, namely:
Figure BDA0003227379120000081
wherein the related parameter is assigned with an initial value a1＝[0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8；0.8]，α＝0.65,β＝0.1，b1＝[0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2；0.2]. Next, assume that the sampling frequency of RTU is 1L sampling of 5 time slicesRThe sampling frequency of PMU is 1 sampling times L for 1 time section (5)P1, M is LR and LPThe least common multiple M of (5). Then, with the 1 st time section as the initial sampling time, as shown in fig. 3, at this time, both PMU and RTU update the measurement value, the RTU between 2 nd to 5 th time sections does not have the measurement value, PMU updates the measurement value in real time, PMU and RTU update the measurement value at the 6 th time section, and so on. The mixing measurements are described as:
Figure BDA0003227379120000082
wherein ,
Figure BDA0003227379120000083
n＝0,1,2,...，
Figure BDA0003227379120000084
Figure BDA0003227379120000085
Pij(θR,UR)、Pi(θR,UR)、URcan be calculated from the correlation measurement function.
And then, carrying out dynamic state estimation on the power system by combining with the extended Kalman filtering: when in use
Figure BDA0003227379120000086
Time, state prediction estimate
Figure BDA0003227379120000087
Wherein, it is required to
Figure BDA0003227379120000088
and P11And assigning an initial value.
When in use
Figure BDA0003227379120000089
And meanwhile, the SCADA updates data and PMU data are synchronously updated, and the state prediction estimation value is as follows:
Figure BDA0003227379120000091
then, by constructing a dynamic model of the power system, an RTU measurement equation, a PMU measurement equation and a mixed measurement expression, parameter identification and state prediction filtering, the dynamic state estimator can estimate the running state of the system according to the known measurement. And finally, if the state of the power system is estimated in 30 time sections, when k is larger than 30, the simulation is finished, the voltage amplitude and the phase angle of the node 7 are output, and otherwise, the dynamic state estimation is continued.
Furthermore, an index is defined: average relative error value epsilon of voltage amplitude estimated valuevi_gjAnd overall system performance index JkThe description is as follows:
Figure BDA0003227379120000092
Figure BDA0003227379120000093
wherein ,Vi fil(k) Is the filtered estimate of the ith voltage amplitude at time k, Vi T(k) Is the true value of the ith voltage amplitude at time k.
The results show that:
the test results are shown in fig. 5, fig. 6 and fig. 7, and fig. 5 shows that the voltage amplitude of the node 7 of the standard 9-node test system is used as the state monitoring quantity to carry out dynamic state estimation on the voltage amplitude, so that the difference between the state estimation value obtained by the method and the actual state value is very small, which indicates that the estimation precision is high, the real-time predictability is strong, and the running state of the system can be tracked well; FIG. 6 shows a graph of the average error percentage of the voltage amplitude of the node 7 in a certain sampling period, wherein it is seen that the error is kept within 0.03, and the smaller the error is, the more accurate the estimation precision is; fig. 7 shows that the variation trend of the overall performance index of the system is calculated by taking the measurement value of the node 7 as an example, the variation trend is actually the variance value between the state estimation value and the true value, the smaller the variance value is, the more accurate the estimation precision is, and the variance tends to 0 in the figure under a certain sampling period, which shows that the accuracy of the dynamic state estimation of the system is improved, and a better estimation algorithm is obtained. In addition, the method has short simulation time and high estimation speed in the whole state estimation process.
In summary, on the basis of the dynamic state estimation research of the power system based on the extended Kalman filtering, the method and the system utilize the measurement data from the SCADA and the PMU to quickly predict the operation state of each node of the power system in real time, overcome the problem that different measurement devices have different sampling frequencies in the process of estimating the state of the power grid after RTU and PMU are measured in a mixed mode, so that the sampling frequency is asynchronous when multi-source data are fused, ensure the consistency of measurement information, have high estimation precision and strong real-time predictability, and provide data support for improving the measurement precision of the power distribution network and better predicting the future operation state information of the system, and the characteristics have important significance for the construction of a future intelligent power grid.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.

Claims (4)

1. The method for estimating the dynamic state of the power system based on PMU and SCADA mixed measurement is characterized by comprising the following steps of:
a. describing a dynamic state equation of the power system through a quasi-steady state model;
b. the RTU measurement equation expressed by the sampling rate is constructed by combining RTUs with different sampling rates, PMU measurement devices, active power among all node lines, node active power injection and voltage amplitude values, and then the PMU measurement equation expressed by the sampling rate is constructed according to the configured PMU, so that a mixed measurement equation of the RTU and the PMU is constructed;
c. and carrying out dynamic state estimation on the power system on the basis of extended Kalman filtering by using a multi-source asynchronous mixed measurement equation.
2. The method for estimating the dynamic state of the power system based on PMU and SCADA hybrid measurement according to claim 1, wherein the quasi-steady state model of the power system under the standard S-node system research in the step a is constructed as follows:
xk+1＝Fkxk+uk+wk
wherein ,xkIs (2S-1) x 1 dimensional state quantity at the moment k, namely the voltage phase angle and the voltage amplitude of each node of the system; fk() Is an adaptive state transfer function; u. ofkIs (2S-1) x 1 dimension input parameter item; w is akIs (2S-1). times.1 dimensional process noise and follows a zero mean normal distribution, i.e., w to N (0, Q), with Q being a ((2S-1). times.1). times. ((2S-1). times.1) dimensional system error covariance matrix.
3. The method for estimating the dynamic state of a power system based on PMU and SCADA hybrid measurement according to claim 1, wherein the RTU measurement equation, the PMU measurement equation and the hybrid measurement expression are constructed in step b as follows:
the RTU measurement vector is composed of active power among lines of each node, active power of each node and node voltage amplitude, and the active power among the lines and the active power equation of the node are expressed as follows:
Pij(θij,Uij)＝Ui 2(gsi+gij)-UiUj(gijcosθij+bijsinθij)
Figure FDA0003227379110000011
wherein ,gij、bijRespectively, the conductance and susceptance of line ij; gsiIs the conductance of the shunt leg ij; gij、BijRespectively the real part and the imaginary part of the element ij in the admittance matrix;
combined node voltage amplitude Ui(Ui) Constructing a measurement vector of the RTU:
Figure FDA0003227379110000012
linearizing the research object, omitting terms of second order and above by using Taylor expansion to obtain a linearized RTU measurement vector, considering the existence of measurement errors, and the RTU measurement equation is as follows:
Figure FDA0003227379110000013
zR(k)＝HRxR(k)+vR(k)
wherein ,zR(k) Is (m)RX 1) measurement vector of dimension RTU; hRIs (m)RA state quantity coefficient matrix of dimension x ((2S-1) × 1)); x is the number ofR(k) The state quantity is (2S-1) multiplied by 1 dimension when only the RTU measuring device is installed; v. ofR(k) Is (m)RX 1) dimensional measurement of error vector and obeying a zero-mean normal distribution, i.e. vR～N(0,RR)，RRIs (m)R×mR) Dimensional measurement error variance matrix, x0Is a linearization point;
considering the different sampling frequencies of different measuring devices, the general RTU measurement equation is written as follows:
ZR(ks)＝HRxR(ks)+vR(ks)
wherein ,ks＝LRk represents a slow rate time index, the whole sampling period is h time sections, zR(k) Has a sampling period of TR＝LRh；zP(k) Has a sampling period of TP＝LPh，Let M be LR and LPLeast common multiple of ql＝M/LlL is R, and P is the sampling number of RTU and PMU in a common period respectively;
constructing PMU measurement equation, and performing standard S node system researchpFrom the element aijFormed and at node l1,l2,...,llThe PMU is installed, and the measurement vector is constructed as follows:
Figure FDA0003227379110000021
Figure FDA0003227379110000022
in the formula ,zP(k) Is a measurement vector of (2 lx 1) -dimensional PMU, HPA state quantity coefficient matrix of (2 lx((2S-1) × 1)) dimension, xP(k) Is the state quantity of (2S-1) × 1D separately installed PMU measuring device, vP(k) Is a (2 lx 1) -dimensional measurement error vector and follows a normal distribution of zero mean, i.e., vP～N(0,RP)，RPMeasuring an error variance matrix in (2l multiplied by 2l) dimensions;
and quantifying the measurement expression by combining the sampling frequency of the PMU: assuming that l PMUs are installed, at l1,l2,...,llThe node is installed, and the amplitude and the phase angle of the bus voltage and the amplitude and the phase angle of the branch current at the system node where the measuring device is installed are mainly measured; taking into account the presence of measurement noise vPThe PMU measurement equation is easily obtained:
zP(ks')＝HPxP(ks')+vP(ks')
wherein ,ks'＝LPk；
Constructing a mixed measurement expression, updating SCADA data and PMU data at the moment of mixed state estimation, wherein the mixed measurement is represented by measurement equations of RTU and PMU together, and PMU is used for every LPSampling once every moment, RTU every LRTime of day samplingOnce, the PMU and the SCADA/RTU can acquire measured values at every M moments, and at the moment, the mixed measurement expression is described as follows:
Figure FDA0003227379110000023
in the formula ,
Figure FDA0003227379110000024
Zmix(k) is ((m)RA +2l) x 1) -dimensional mixture measurement vector, x (k) is a state quantity under (n x 1) -dimensional mixture measurement, vmix(k) Is ((m)R+2l) x 1) dimensional mixed measurement error vector and obeys a zero-mean normal distribution, i.e., vmix～N(0,Rm)，RmIs ((m)R+2l)×(mR+2l)) dimension measurement error variance matrix;
besides, when only data updating of a fast sampling rate measuring device PMU and only data updating of a slow sampling rate measuring device RTU are considered, it is easy to show that multi-source asynchronous measurement is as follows:
Figure FDA0003227379110000025
4. the method for estimating the dynamic state of the electric power system based on PMU and SCADA measurement according to claim 1, wherein in the step c, the dynamic state estimation is performed on the electric power system based on the extended Kalman filtering by using the multisource asynchronous mixed measurement equation, and the state transition quantity F is identified by using a Holt's two-parameter index smoothing modelkAnd input parameter item ukState of the system under standard S-node system research
Figure FDA0003227379110000031
Sum covariance matrix Pk+1|kThe prediction result is as follows:
Figure FDA0003227379110000032
Pk+1|k＝FkPk|kFk T+Qk
wherein ,
Figure FDA0003227379110000038
the predicted value of the state at the time k +1 by (2S-1) × 1 dimension,
Figure FDA0003227379110000034
updating the correction value for the state of time k versus time k, Pk+1|kAn error covariance prediction matrix for the time k to the time k + 1;
using measurements at multiple rates, real-time measurements Z have been obtained for a set of power systemskFiltering to obtain new state estimation vector by calculating Kalman gain and corresponding innovation value
Figure FDA0003227379110000035
The state filtered value is:
Kk+1＝Pk+1|kHθ T[HθPk+1|kHθ T+Rθ]-1
Figure FDA0003227379110000036
Figure FDA0003227379110000037
Pk+1|k+1＝[In-Kk+1Hθ]Pk+1|k
if the state of the power system is estimated in N time sections, when k is larger than N, the simulation is finished, the voltage amplitude and the phase angle of each node are output, and otherwise, the state is continuously updated and filtered.
CN202110975112.9A 2021-08-24 2021-08-24 Power system dynamic state estimation method based on PMU and SCADA hybrid measurement Active CN113675951B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110975112.9A CN113675951B (en) 2021-08-24 2021-08-24 Power system dynamic state estimation method based on PMU and SCADA hybrid measurement

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110975112.9A CN113675951B (en) 2021-08-24 2021-08-24 Power system dynamic state estimation method based on PMU and SCADA hybrid measurement

Publications (2)

Publication Number Publication Date
CN113675951A true CN113675951A (en) 2021-11-19
CN113675951B CN113675951B (en) 2023-05-12

Family

ID=78545630

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110975112.9A Active CN113675951B (en) 2021-08-24 2021-08-24 Power system dynamic state estimation method based on PMU and SCADA hybrid measurement

Country Status (1)

Country Link
CN (1) CN113675951B (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106707061A (en) * 2016-12-16 2017-05-24 湖南大学 Hybrid measurement based power distribution network dynamic state estimation method
CN113255959A (en) * 2021-04-09 2021-08-13 广东电网有限责任公司电力调度控制中心 Power system dynamic state estimation method and system

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106707061A (en) * 2016-12-16 2017-05-24 湖南大学 Hybrid measurement based power distribution network dynamic state estimation method
CN113255959A (en) * 2021-04-09 2021-08-13 广东电网有限责任公司电力调度控制中心 Power system dynamic state estimation method and system

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
卓亮;陈利跃;何星;: "基于混合量测的动态状态估计算法研究" *
李波;何星;方兴其;: "混合量测下序贯抗差状态估计算法研究" *
葛泉波;汪国安;汤天浩;文成林;: "基于有理数倍采样的异步数据融合算法研究" *
邱爱兵;文成林;姜斌;: "基于异步多传感器采样量测的最优状态融合估计" *

Also Published As

Publication number Publication date
CN113675951B (en) 2023-05-12

Similar Documents

Publication Publication Date Title
CN110471024B (en) Intelligent electric meter online remote calibration method based on measurement data analysis
CN108448568B (en) Power distribution network hybrid state estimation method based on multiple time period measurement data
CN106921156A (en) A kind of active distribution network method for estimating state based on many sampling period hybrid measurements
CN106707061A (en) Hybrid measurement based power distribution network dynamic state estimation method
CN108155648A (en) Method for estimating state based on the infinite Extended Kalman filter of adaptive H
CN107565553A (en) A kind of power distribution network robust dynamic state estimator method based on UKF
CN107658881A (en) Voltage stability critical point determination methods based on Thevenin's equivalence method
CN100554976C (en) Area voltage stability monitoring method based on synchronous phasor measurement
CN113595058A (en) Processing method of power distribution network state estimation measurement data
CN110707693B (en) Aggregate Kalman filtering dynamic state estimation method based on AMI full measurement point partition
CN111581768A (en) Power distribution network distributed state estimation method based on hybrid measurement
CN109560550A (en) The mains by harmonics method for estimating state measured based on optimization
CN115081299A (en) UPF-based robust auxiliary prediction state estimation method for power system
CN115453193A (en) Power distribution network harmonic state estimation method based on cooperation of PQM, TTU and SM measurement data
CN106372440B (en) A kind of adaptive robust state estimation method of the power distribution network of parallel computation and device
CN104537233A (en) Distribution network pseudo measurement generating method based on kernel density estimation
CN113962053A (en) Power distribution network state evaluation method based on multi-section intelligent instrument data
CN110417009A (en) Power distribution network hybrid robust state estimation method based on data of different sampling periods
CN109586289B (en) Power distribution network multi-time scale recursive dynamic state estimation method and system
CN113675951A (en) Power system dynamic state estimation method based on PMU (phasor measurement Unit) and SCADA (Supervisory control and data acquisition) mixed measurement
Yang et al. Dynamic variable-weight least squares for state estimation of distribution network based on data fusion
Zeraati et al. Meter placement algorithms to enhance distribution systems state estimation: Review, challenges and future research directions
CN110021931A (en) It is a kind of meter and model uncertainty electric system assist predicted state estimation method
CN113595078B (en) Power distribution network state estimation method and device based on multi-source mixed data fusion
CN115693763A (en) State prediction-based alternating current-direct current power distribution network linear interval state estimation method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant