The present disclosure generally relates to electric power systems and synchronized phasor (or synchrophasor) measurement methods for rapid situational awareness of electric power systems.
Synchronized measurements of currents and voltages at nodes throughout electric power systems, also known as power grids, are typically limited by slow responses and low reporting rates of devices used for the measurements. An example of such devices is a Phasor Measurement Unit (PMU). Conventional PMUs employ measurement algorithms that are generally based on discrete Fourier transform (DFT) and that estimate a frequency and/or a phase angle from either current or voltage measurements over a window spanning six line cycles (i.e., 0.1 seconds for a 60-Hz power grid). Moreover, given that the IEEE C37.118 standard for PMUs prescribes a low reporting rate requirement, commercial PMUs are designed and manufactured to have typical data rates ranging from 10 Hz to 60 Hz, which are not fast enough for dynamic response prediction and transient stability control of power grids. Consequently, neither fast dynamic response (e.g., within one line cycle) nor reliable and accurate measurements during transients (due to sudden changes in power generation and/or power-consuming loads in the power grids) may be provided by current commercial PMUs. Therefore, conventional PMUs are not suited for transient stability control.
Low-latency frequency and angle measurements may benefit power grid situational awareness, event analysis, and transient stability control. Furthermore, they may extend power grid visibility during transients, dynamics, oscillation, first swing, frequency instability, voltage instability, etc. High data rate may allow for dynamic response prediction and may increase the accuracy of power system model validation. Therefore, the inventors recognized a need in the art for fast synchronized measurement methods for rapid situational awareness of power grids.
An embodiment of the present disclosure provides a synchrophasor measurement method for a device configured to take synchronized measurements in a power system. The synchrophasor measurement method includes receiving global positioning system (GPS)-synchronized samples of a signal sensed by the device from the power system; determining a level of distortion of the signal; selecting, based on the level of distortion, a computation method, the computation method being one of an improved zero-crossing (IZC) method and an enhanced phase-lock-loop (EPLL) method; performing the selected computation method to determine at least one parameter of the signal at a reporting frequency, which is at least twice a line frequency of the power system; and outputting, at the reporting frequency, the at least one parameter to an operator of the power system to allow the operator to perform at least one of a monitoring and a controlling of at least one element of the power system.
Another embodiment of the present disclosure provides an electronic device for a power system. The electronic device comprise a sensor, a GPS receiver, one or more processors, and storing instructions adapted to be executed by the one or more processors to perform the synchrophasor measurement method.
The synchronized measurement devices 110 are typically sparsely installed to perform the local measurements at different locations across the power grid. The synchronized measurement devices 110 generally include voltage and/or current sensors to make voltage and/or current measurements. The synchronized measurement devices 110 may include one or more processors and memory storing instructions to be executed by the one or more processors to determine the frequency, phase angle, etc. of the voltage and/or current measurements locally. It is to be appreciated that the synchronized measurement devices 110 are not limited to any particular device, and may refer to any sensor that uses synchrophasor measurement technology. Each synchronized measurement devices 110 may include a GPS receiver providing a GPS-based synchronization signal (e.g., pulse per second (PPS), inter-range instrumentation group B (IRIG-B)), which may be used for sampling control in the synchronized measurement device 110.
The Internet 120 may serve as a wide-area communication network (WAN) 122 with a plurality of firewalls/routers 124 to connect the synchronized measurement devices 110 and one or more clients 126 to the IMS 130. The IMS 130 may collect the sampled measured data and/or the computed data from the synchronized measurement devices 110, store the data in databases in data storage devices 132, and provide a platform for analyses of the data either before or after storing the data. The servers 134-137 in the IMS 130 may include a plurality of processors to manipulate and analyze the stored data serially and/or in parallel. The servers 134-137 may be centrally or distributedly located. Data generated from the analyses of the stored data may also be stored in the data storage devices 132. The data storage devices 132 may include secondary or tertiary storage to allow for non-volatile or volatile storage of the measured, computed, and generated data. The IMS 130 may be entirely contained at one location or may also be implemented across a closed or local network, an internet-centric network, or a cloud platform.
The IMS 130 may provide the measured and/or computed data, either before or after storing the data, to the one or more clients 126. Alternatively, the sampled measured data and/or the computed data may be transmitted directly from the synchronized measurement devices 110 to the clients 126 via the Internet 120. The data received at the clients 126 may be displayed on one or more displays and/or further manipulated by one or more processors for analysis, monitoring, and control of the power grid either automatically by computer systems or visually and manually by one or more operators of the power grid. Based on the synchrophasor data, one or more devices/equipment such as generators, transformers, switches, capacitors, transmission lines, power-consuming loads, etc. in the power grid may be controlled automatically by the computer systems or manually by the one or more operators, in order to maintain the stability and safe operation of the power grid.
To allow for fast synchrophasor measurements for rapid situational awareness, event analysis, and transient stability control, it is desirable for the synchronized measurement devices 110 to have sampling rates, which are higher than those of conventional PMUs, for faster response. For example, the sampling rates of the synchronized measurement devices 110 may range from 720 Hz to 100 kHz. The synchronized measurement devices 110 may be configured to adjust their sampling rates, for example, based on the quality of the signal(s) being measured. Furthermore, it is desirable for the synchronized measurement devices 110 to be equipped with measurement methods/algorithms that are capable of continuously measuring/estimating the frequency, phase angle, etc. of the power grid as sample(s) of voltage and/or current signals become available. Such synchronized measurement devices 110 thus do not require a window spanning six line cycles to estimate the frequency, phase angle, etc. and may have reporting rates of 120 Hz or higher, which are much faster than the typical PMU data rates that range from 10 Hz to 60 Hz.
(x(i)−D)×(x(i+1)−D)≤0 (1)
In equation (1), x(i) and x(i+1) represent two consecutive samples of the measured signal x. The parameter D represents a direct-current (DC) offset in the measured signal x and may be determined, for example, by continuously or periodically averaging samples of the signal x over at least one line cycle. A zero crossing location may be identified to be between the ith and (i+1)th samples whenever the inequality in equation (1) is satisfied. From consecutive zero crossing locations, frequency and phase angle of the measured signal may be determined.
However, it is common for power grid signals to contain ˜40-80 dB white noise. As the noise level of a measured signal increases, its waveform may not be a monotonic function and there may be several samples around a zero crossing, as illustrated in
x(n)=α0+α1n (2)
At step 504, given a sample vector x=[x(1), x(2), . . . , x(m)] of m samples around each zero crossing location, the IZC method generates a constant-valued matrix M by converting equation (2) into matrix form as in equation (3).
x=αM (3)
where a coefficient vector α=[α0, α1]. The constant-valued matrix M may thus be expressed as:
Since the matrix M and the sample vector x are known, the IZC method 500 may determine, at step 506, the coefficient vector α by the following pseudo-inverse computation:
α=[MTM]−1MTx (5)
Once the coefficient vector α determined, at step 508, the IZC method 500 interpolates between the m samples to generate a polynomial function and determines a zero crossing location as the point where the polynomial function is zero. Every time the IZC method 500 determines a zero crossing location, at step 510, the IZC method 500 also determines a first time interval T1 between the zero crossing location (i.e., the last zero crossing location) and the second-to-last zero crossing location, and a second time interval T2 between the second-to-last zero crossing location and the third-to-last zero crossing location, as in equation (6).
T
1
=T
s×(sample count of last zero crossing−sample count of second-to-last zero crossing)
T
2
=T
s×(sample count of second-to-last zero crossing−sample count of third-to-last zero crossing) (6)
In equation (6), Ts is the sampling period of the synchronized measurement device.
At step 512, the IZC method 500 computes the frequency f of the measured signal x using equation (7).
After calculating the frequency f the IZC method 500 calculates the phase angle φ of the measured signal x as follows at step 514:
In equation (8), T represents the time between the last zero-crossing location and a reference time, e.g., a PPS signal from a GPS receiver of the synchronized measurement device.
Since the measured signal x may have a DC offset, the minimum window length of the IZC method is one line cycle (i.e., 1/60 seconds for a 60-Hz power grid). It is important for the IZC method 500 to, continuously or periodically, remove the DC offset from the measure signal x prior to determining the zero crossing locations. The DC offset in the measured signal x may be determined, for example, by continuously or periodically averaging samples of the measured signal x over at least one line cycle. Given that the IZC method 500 calculates the frequency f and the angle φ at every zero crossing location, the maximum output rate of the IZC method 500 is 120 Hz, which is faster than conventional DFT output rate of 10 Hz to 60 Hz. Accordingly, the IZC method 500 is suitable for rapid power grid situational awareness, event analysis, and transient stability control.
It is to be appreciated that, for the IZC method 500, the resolution is directly related to the sampling rate. The resolution r of the IZC method can be expressed as:
where f0 is the nominal power grid frequency (e.g., 60 Hz), fs is the sampling frequency (i.e., fs=1/Ts), and “floor( )” is a function that rounds its operand down to the nearest integer. The limitation of the resolution r is due to the integer restraint when counting the samples between two zero-crossing locations and the fact that, typically, a sample is not identically zero. Increasing the sampling rate is the most effective method to improve the frequency measurement accuracy for the IZC method.
It is to be also appreciated that the implementation of the IZC method 500 is not limited to any particular programming language or execution environment, and the method 200 may be applied to any computer programming languages or logic.
An EPLL method estimates the magnitude, phase angle, angular frequency, and DC offset of a measured signal u, which may be represented by equation (9), as one or more samples of the measured signal u become available. The EPLL method generates an estimated signal y represented by equation (10), and corrects the estimated signal y based on an error e between the measured signal u and the estimated signal y such that the error e is minimized.
u=A cos(ωt+ϕ)+D=A cos(θ)+D (9)
y=Â cos({circumflex over (ω)}t+{circumflex over (ϕ)})+{circumflex over (D)}=Â cos({circumflex over (θ)})+{circumflex over (D)} (10)
In equation (9), A is the magnitude, ω is the angular frequency, ϕ is the initial phase angle, θ=ωt+ϕ is the phase angle, and D is the DC offset of the measured signal u. Similarly, in equation (10), Â is the estimated magnitude, {circumflex over (ω)} is the estimated angular frequency, {circumflex over (ϕ)} is the estimated initial phase angle, {circumflex over (θ)}={circumflex over (ω)}t+{circumflex over (ϕ)} is the estimated phase angle, and {circumflex over (D)} is the estimated DC offset of the estimated signal y.
The error e between the measured signal u and the estimated signal y may be denoted as:
e=u−y (11)
An objective function J may be formulated as in equation (12) to minimize the error e.
J=e
2=(u−y)2 (12)
At a minimum value of the objective function J, the following first-order optimality conditions are satisfied:
where P=[A ω θ D]. Using an optimization method, such as the steepest descent method, the differentiation of each variable in P may be obtained as in equations (14)-(17).
{dot over (A)}=2e cos {circumflex over (θ)} (14)
{dot over (ω)}=−2eA sin {circumflex over (θ)} (15)
{dot over (θ)}=ω (16)
{dot over (D)}=2e (17)
Based on equations (14)-(17), in the discrete domain, the estimated magnitude  may be determined first by using equation (18), where e(n) may be determined from equation (11), n being the discrete time index.
Â(n)=k1e(n)cos θ(n−1)Ts+Â(n−1) (18)
The integral of the estimated angular frequency {circumflex over (ω)}, {circumflex over (ω)}int may then be determined as follows:
{circumflex over (ω)}int(n)=−k3e(n)Â(n)sin {circumflex over (θ)}(n−1)Ts+{circumflex over (ω)}int(n−1) (19)
The estimated angular frequency {circumflex over (ω)} may be determined using equation (20).
{circumflex over (ω)}(n)=−k2e(n)Â(n)sin {circumflex over (θ)}(n−1)+{circumflex over (ω)}int(n) (20)
Equation (21) may then be used to determine the estimated phase angle {circumflex over (θ)}.
{circumflex over (θ)}(n)={circumflex over (ω)}(n)Ts+{circumflex over (θ)}(n−1) (21)
The integral of the estimated DC offset {circumflex over (D)}, {circumflex over (D)}int may be determined as in equation (22) and the estimated DC offset {circumflex over (D)} may be determined as in equation (23).
{circumflex over (D)}
int(n)=k5e(n)Ts+{circumflex over (D)}int(n−1) (22)
{circumflex over (D)}(n)=k4e(n)+{circumflex over (D)}int(n) (23)
In the equations (18)-(23), Ts is the sampling period of the synchronized measurement device. The coefficients k1, k2, k3, k4, and k5 determine the speed of convergence and the steady-state stability of the EPLL method. An estimated frequency {circumflex over (f)} of the estimated signal y may also be determined as follows:
{circumflex over (f)}(n)=2π{circumflex over (ω)}(n) (24)
The EPLL method 600a starts at step 602 by initializing the estimated variables described above and also intermediate variables used in generating the estimated variables. The estimated variables are initialized to nominal values. For example, for a power grid with a 120-V 60-Hz nominal voltage, the estimated magnitude  is initialized to 120 (i.e., Â(0)=120); the estimated angular frequency {circumflex over (ω)} is initialized to 120π (i.e., {circumflex over (ω)}(0)=120π); the estimated phase angle {circumflex over (θ)} is initialized to 0 (i.e., {circumflex over (θ)}(0)=0); the estimated DC offset {circumflex over (D)} is initialized to 0 (i.e., {circumflex over (D)}(0)=0). Any intermediate variable is generally initialized to 0. For example, the integral {circumflex over (ω)}int of estimated angular frequency is initialized to 0 (i.e., {circumflex over (ω)}int(0)=0) and the integral {circumflex over (D)}int of the estimated DC offset is initialized to 0 (i.e., {circumflex over (D)}int(0)=0).
Once the variables initialized, the EPLL method 600a receives one or more samples of the measured signal u at step 604. The EPLL method 600a generates the estimated signal y at step 606 according to equation (10) and calculates the error e at step 608 according to equation (11). With the calculated error e, the EPLL method 600a calculates the estimated magnitude  at step 610 using equation (18), and the estimated DC offset {circumflex over (D)} at step 614 using equations (22) and (23). With the calculated error e from step 608 and the calculated estimated magnitude  from step 610, the EPLL method 600a calculates the estimated angular frequency {circumflex over (ω)} at step 612 based on equations (19) and (20). The EPLL method 600a may also calculate the estimated frequency {circumflex over (f)} at step 612 according to equation (24). At step 616, the EPLL method 600a then calculates the estimated phase angle {circumflex over (θ)} according to equation (21).
The EPLL method 600a updates the variables at step 618 and outputs the estimation results (for the synchronized measurement device to output/transmit externally, for example), in particular the estimated phase angle {circumflex over (θ)} and the estimated frequency {circumflex over (f)}, at step 622. The EPLL method 600a then repeats steps 604-622 when subsequent one or more samples of the measured signal u are received.
k
j′(n)=ajkj(n−1)+rje(n)2 (25)
In equation (25), 0<aj<1 and rj>1. The parameters aj and rj may be chosen based on a desired performance of the EPLL algorithm and characteristics of the measured signal u. For example, a1=0.9, r1=4×10−4, a2=0.9, r2=4.8×10−2, a3=0.9, r3=6×10−2a4=0.9, r4=4.8×10−3, as=0.9, and r5=4.8×10−4.
Once the intermediate coefficients kj′ calculated, the EPLL method 600b determines the coefficients kj as follows:
In equation (26), kjmin and kjmax respectively represent lower and upper boundaries for the corresponding coefficient kj, and may be enforced to keep the EPLL algorithm stable and accurate in situations where there may be disturbances or computational errors, for example. An example of a set of boundaries is k1max=2000, k1min=400, k2max=0.7, k2min=0.3, k3max=60, k3min=25, k4max=0.08, k4min=0.02, k5max=0.002, and k5min=0.0005.
Thereafter, the EPLL method 600b outputs the estimation results at step 622, and then repeats steps 604-622 when subsequent one or more samples of the measured signal u are received. Steps 604-618 and 622 of the EPLL method 600b are identical to steps 604-618 and 622 of the EPLL method 600a. Therefore, descriptions of steps 604-618 and 622 of the EPLL method 600b will be omitted.
As illustrated in
The maximum number of iterations M may be chosen based on a desired performance and a desired computation speed of the EPLL algorithm. A larger M may improve the dynamic performance of the EPLL algorithm, but may increase the computation burden. An exemplary value of M is 10.
With an iteration index i, equations (10), (11), and (18)-(24) may be modified as equations (27)-(35), respectively, as follows:
y
i(n)=ÂM(n−1)cos({circumflex over (θ)}M(n−1))+{circumflex over (D)}M(n−1) (27)
e
i(n)=u(n)−yi(n) (28)
Â
i(n)=k1(n)ei(n)cos θ(n−1)Ts+Âi-1(n) (29)
{circumflex over (ω)}inti(n)=−k3(n)ei(n)Âi(n)sin {circumflex over (θ)}i-1(n)Ts+ωint(i-1)(n) (30)
{circumflex over (ω)}i(n)=−k2(n)ei(n)Âi(n)sin {circumflex over (θ)}i-1(n)+{circumflex over (ω)}int(i-1)(n) (31)
{circumflex over (θ)}i(n)={circumflex over (ω)}i(n)Ts+{circumflex over (θ)}M(n−1) (32)
{circumflex over (D)}
inti(n)=k5(n)ei(n)Ts+{circumflex over (D)}int(i-1)(n) (33)
{circumflex over (D)}
i(n)=k4(n)ei(n)+{circumflex over (D)}int(i-1)(n) (34)
{circumflex over (f)}
i(n)=2π{circumflex over (ω)}i(n) (35)
Therefore, the EPLL method 600c invokes equations (27)-(35) accordingly in steps 606-616 for i={1, 2, . . . , M} iterations. The coefficients k1, k2, k3, k4 and k5 may be fixed as in EPLL method 600a, in which case the step 620 is not needed in the EPLL method 600c. However, the coefficients k1, k2, k3, k4 and k5 may also be updated in every iteration i, as described above with respect to the EPLL method 600b, in which case the EPLL method 600c includes the step 620.
It is to be appreciated that the EPLL methods 600a-600c may compute and output the estimated frequency {circumflex over (f)} and phase angle {circumflex over (θ)} whenever a sample of the measured signal u is available. Thus, if the sampling rate is between 720 Hz and 36 kHz for any of the EPLL methods 600a-600c, the corresponding reporting rate will also be between 720 Hz and 36 kHz, which is much faster than the conventional DFT output rate of 10 Hz to 60 Hz. Accordingly, the EPLL methods 600a-600c are suitable for power grid situational awareness, event analysis, and transient stability control. It is to be also appreciated that the implementations of the EPLL methods 600a-600c are not limited to any particular programming language or execution environment, and the method 200 may be applied to any computer programming languages or logic.
Embodiments of the disclosure are specifically illustrated and/or described herein. However, it will be appreciated that modifications and variations of the disclosure are covered by the above teachings and within the purview of the appended claims without departing from the spirit and intended scope of the disclosure. Further variations that are consistent with the principles described above are permissible.
This invention was made with government support under the U.S. Department of Energy Grid Modernization Laboratory Consortium (GMLC) and EEC-1041877 awarded by the U.S. National Science Foundation and the U.S. Department of Energy. The U.S. Government has certain rights in this invention.