CN109638862B - CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method - Google Patents
CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method Download PDFInfo
- Publication number
- CN109638862B CN109638862B CN201811647265.5A CN201811647265A CN109638862B CN 109638862 B CN109638862 B CN 109638862B CN 201811647265 A CN201811647265 A CN 201811647265A CN 109638862 B CN109638862 B CN 109638862B
- Authority
- CN
- China
- Prior art keywords
- oscillation
- mode
- frequency
- energy
- imf
- 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.)
- Active
Links
- 230000010355 oscillation Effects 0.000 title claims abstract description 155
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 53
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000013016 damping Methods 0.000 claims abstract description 50
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 45
- 238000005259 measurement Methods 0.000 claims abstract description 38
- 238000004364 calculation method Methods 0.000 claims abstract description 21
- 230000003044 adaptive effect Effects 0.000 claims abstract description 8
- 238000005070 sampling Methods 0.000 claims description 15
- 238000012935 Averaging Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000003909 pattern recognition Methods 0.000 claims 1
- 230000009466 transformation Effects 0.000 abstract description 5
- 238000004870 electrical engineering Methods 0.000 abstract description 2
- 238000004458 analytical method Methods 0.000 description 12
- 238000010586 diagram Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 239000001257 hydrogen Substances 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 230000001052 transient effect Effects 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- -1 hydrogen Chemical class 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000011217 control strategy Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003534 oscillatory effect Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
-
- 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]
Abstract
The invention relates to the field of electrical engineering, in particular to a CEEMDAN algorithm-based power system low-frequency oscillation mode identification method. And acquiring actual measurement data of the power system by using the wide area measurement system, wherein the actual measurement data comprises a rotor angle signal of the generator. The Complete set Empirical Mode Decomposition with Adaptive Noise of each group of original low-frequency oscillation signals is decomposed into a plurality of Intrinsic Mode functions and the sum of IMFs by a CEEMDAN algorithm, wherein each IMF component represents an oscillation Mode. And calculating the energy value and the energy weight of each IMF component, identifying the oscillation frequency and the damping ratio of the dominant oscillation mode by using Hilbert-Huang transformation, and comparing the calculation result with a characteristic value method to ensure the safe and stable operation of the power system.
Description
Technical Field
The invention relates to the field of electrical engineering, in particular to a CEEMDAN algorithm-based power system low-frequency oscillation mode identification method.
Background
Along with the continuous increase of the scale of an interconnection system, the large-capacity and ultra-long-distance alternating-current and direct-current power transmission is continuously increased, the running mode is complicated and changeable, and the power grid accident risk and the oscillation control difficulty caused by low-frequency oscillation are continuously increased. The low frequency oscillation poses a great threat to power supply equipment, and even cascading failure can be induced to cause large-area power failure. Therefore, the method for identifying the low-frequency oscillation mode of the power system under the background of continuously expanding the scale of the interconnected power grid has very important practical significance and engineering practical value.
PMU is phase Measurement Unit, and the Phasor Measurement Unit is synchronized. With the large-scale configuration of PMUs, a wide area measurement system (wide area measurement system, WAMS) is widely applied to an interconnected power grid, the WAMS can record real-time operation data of each site of a large-scale interconnected power system at the same reference time, and the data enable low-frequency oscillation analysis and control of the power system to be possible. According to different types of measured signals in the power system, the power system dynamic stability analysis method based on the wide-area measurement information comprises the following steps: a test signal based dynamic stability analysis method, a fault signal based dynamic stability analysis method, and a noise signal based dynamic stability analysis method. The dynamic stability analysis method based on the fault signal comprises the following steps: the method mainly estimates a transfer function matrix, a state matrix or multi-oscillation mode signal fitting of the power system from fault signals of the power system to identify a dominant oscillation mode of the power system, and commonly used algorithms include classical modal decomposition, Hilbert Spectral Analysis (HAS), minimum feature implementation (ERA), Kalman filtering, continuous wavelet transformation, random subspace identification and TLS-ESPRIT identification.
Classical modal Decomposition (EMD) is now recognized as a powerful multi-scale analysis tool that deals with non-stationary and non-linear signals. The local nature of the EMD algorithm may produce oscillations with different scales in one mode or oscillations with similar scales in a different mode. When this phenomenon is not preferable, it is called "mode mixing". In order to overcome the situation, an algorithm which is provided on the basis of an EMD algorithm, namely a Complete set Empirical Mode Decomposition with Adaptive Noise, is provided, and a CEEMDAN algorithm is provided, wherein the CEEMDAN algorithm can separate an accurate reconstruction and a clearer Mode of an original signal. The CEEMDAN algorithm, in combination with Hilbert-Huang transform (HHT), analyzes the application in low frequency oscillations in power systems. The method is helpful for analyzing dynamic phenomena in a desired time interval range, and is a method for researching transient characteristics of non-stationary data, such as local characteristics of transient frequency, damping and the like, so as to characterize the time characteristics of nonlinear non-stationary oscillation. The invention introduces CEEMDAN algorithm to process the non-stationary oscillation signal, and simultaneously introduces Teager energy operator as the identification basis to eliminate the false noise component generated in the CEEMDAN decomposition process, thereby effectively improving the accuracy and efficiency of identifying the dominant low-frequency oscillation mode of the power system, and identifying the oscillation frequency and the damping ratio of the oscillation mode by using the Hilbert-Huang transform algorithm.
Disclosure of Invention
The invention provides a method for identifying a low-frequency oscillation mode of a power system based on a CEEMDAN algorithm, in particular to a method for identifying an application related control strategy in the oscillation mode of the power system based on the CEEMDAN algorithm, wherein a dominant oscillation mode is identified by a method of combining the CEEMDAN algorithm with a Teager energy operator, and the frequency and damping ratio of the dominant oscillation mode are identified by a Hilbert-Huang transform algorithm, which is described in detail as follows:
the invention discloses a CEEMDAN algorithm-based power system low-frequency oscillation mode identification method, which comprises the following steps:
1) acquiring measured data of a power system by using a data acquisition function of a wide area measurement system, decomposing a Complete integrated Empirical Mode Decomposition with an Adaptive Noise of each group of original low-frequency oscillation signals into a plurality of Intrinsic Mode functions according to a CEEMDAN algorithm, and the sum of IMFs, wherein each IMF component represents an oscillation Mode;
2) preprocessing each IMF component, calculating the energy size and the energy weight of each IMF component by using a Teager energy operator, wherein in the generator rotor angular oscillation signal, the relative energy of a dominant oscillation mode is larger, and one or more modes with larger relative energy are respectively identified.
3) And identifying the oscillation frequency and the damping ratio of the dominant oscillation mode by using a Hilbert-Huang transform algorithm, comparing the calculation result with a characteristic value method, and verifying the effectiveness of the calculation result.
The step 1 of acquiring the measured data of the power system by using the data acquisition function of the wide area measurement system includes:
acquiring actual measurement data of the power system by using a data acquisition method of a Phasor Measurement Unit (PMU) and a Wide Area Measurement System (WAMS), wherein the actual measurement data comprises a rotor angle signal of each generator or active power on a connecting line as a signal to be identified;
further comprising: and acquiring state measurement information of the power system from the wide area measurement system, wherein the state measurement information comprises rotor angle signals and active power on a connecting line, and standardizing the state measurement information.
The CEEMDAN algorithm in the step 1 comprises the following steps:
1) adding different types of random white Gaussian noise to the original signal x (t), x (t) + epsilon0ωi(t)＝xi(t), where I ═ 1, … I number of times to achieve white gaussian noise, ∈0Is the standard deviation of Gaussian white noise, omegai(t) is a different type of white gaussian noise. For new signal xi(t) performing empirical mode decomposition, averaging I empirical mode decomposition results, and extracting a first intrinsic mode function component IMF1；
2) Decomposing different types of random Gaussian white noise empirical modes to obtain epsilon1Eiωi(t), wherein: epsilon1Is a standard deviation of, EiIs the first order IMF component, ωi(t) is white gaussian noise of different types;
3) subtracting IMF from the original signal1Then adding epsilon1E1ωi(t) forming a new signal, performing empirical mode decomposition on the new signal, averaging the results of the empirical mode decomposition, and extracting a second intrinsic mode function component IMF2；
4) And analogizing until an iteration stop criterion is met, wherein each IMF component corresponds to an oscillation mode.
In the step 2, the Teager energy operator is used for calculating the energy size and the energy weight of each IMF component, the Teager energy operator is used for calculating the energy of all IMF components, and the IMF component energy is obtained, wherein the method comprises the following steps: calculating the energy size and the energy weight of each IMF component by using a Teager energy operator;
for continuous signal xnThe formula for calculating the energy is as follows:
calculating the energy size and energy weight of each IMF component in the step 2: refers to summing the energy values of all discrete samples of the IMF component as the energy of this IMF component.
And calculating the energy weights of all IMF components, keeping the IMF components with high energy weights as a dominant oscillation mode, and removing the IMF components with low energy weights as a non-dominant oscillation mode.
In the step 3, the oscillation frequency and the damping ratio of the dominant oscillation mode are identified by using a Hilbert-Huang transform algorithm, and Hilbert transform is performed on the original signal x (t);
comparing the calculation result with a characteristic value method, and verifying the validity of the calculation result;
1) the original signal x (t), i.e. an eigenmode function IMF component, is Hilbert transformed:
wherein P represents the cauchy principal value integral,
2) combining x (t) and y (t) into a pair of complex numbers, the following signal z (t) is formed:
z(t)＝x(t)+j y(t)＝A(t)ejφ(t)
where A (t) is the amplitude signal, phi (t) is the phase function, x (t) is the original signal, y (t) is the Hilbert transformed signal, z (t) is the complex number formed by x (t) and y (t), whose modulus and phase angle represent the amplitude and phase of the signal.
Identifying the oscillation frequency of the dominant oscillation mode in the step 3;
by definition, instantaneous angular frequency ω (t) is the derivative of phase Φ (t) with respect to time, and instantaneous frequency f (t) ═ ω (t)/2 π, the expression for instantaneous frequency can be obtained:
in the above formula: phi (t) is the phase function and f (t) is the instantaneous frequency.
The final result is represented by x (t) and y (t) and their derivatives, and the final expression for the instantaneous frequency f (t) can be obtained as:
in the above formula: x (t) is the original signal, y (t) is the Hilbert transformed signal, x & (t) is the derivative of x (t), y & (t) is the derivative of y (t), and f (t) is the instantaneous frequency.
The instantaneous frequency is averaged to obtain the oscillation frequency of the mode.
Identifying the damping ratio of the dominant oscillation mode in the step 3:
in the above formula: ζ is the damping ratio, dA (t)/dt is the derivative of the amplitude function,
is the derivative of the phase function.
And averaging the instantaneous damping ratio to obtain the damping ratio in the mode.
The characteristic value method in the step 3 is to calculate the oscillation frequency and the damping ratio of the same system by using the characteristic value method; comparing the calculation result with a characteristic value method, and verifying the validity of the calculation result;
solving a linear equation of the system to obtain a state matrix A, and solving a characteristic value lambda of the matrix Ai＝σi+jωiEach pair of conjugate complex roots corresponds to an oscillation mode, and the oscillation frequency and the damping ratio of the mode are calculated as follows:
in the above formula: lambda [ alpha ]iIs a conjugate eigenvalue; sigmaiRepresenting the damping characteristic of the oscillation as the real part of the characteristic value; omegaiThe characteristic value is an imaginary part of the characteristic value and is used for representing the frequency characteristic of the oscillation, xi is a damping ratio, and f is the oscillation frequency.
In the step 3, the oscillation frequency of the dominant oscillation mode is identified by using a Hilbert-Huang transform algorithm and compared with a characteristic value method;
theoretically, if N generators exist in the system, N-1 low-frequency oscillation modes should exist; however, the oscillation parameters of the dominant oscillation modes decomposed by the rotor angle signals of different generators are the same or similar, and the dominant oscillation modes belong to the same oscillation mode; classifying modes having the same oscillation parameters as one oscillation mode; and comparing with a characteristic value method to verify the effectiveness of the algorithm.
The technical scheme provided by the invention has the beneficial effects that:
the invention provides a dynamic stability collaborative identification method of a power system based on wide-area measurement information by means of CEEMDAN decomposition, so as to realize identification of the dominant oscillation mode frequency and the damping ratio of the power system. And the power grid operation state information is provided for the power system operation scheduling personnel so as to improve the safe and stable operation of the power system.
Drawings
FIG. 1 is a diagram of a 16-machine 68 node test system;
FIG. 2 is a graph of rotor angle swing for a generator with branch 46-49 faults;
FIG. 3 is a raw signal of generator G2;
FIG. 4 is a waveform diagram of IMF1 of generator G2;
FIG. 5 is a waveform diagram of IMF2 of generator G2;
FIG. 6 is a waveform diagram of IMF3 of generator G2;
FIG. 7 is a waveform diagram of IMF4 of generator G2;
FIG. 8 is a waveform diagram of IMF5 of generator G2;
fig. 9 is a block flow diagram of the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention are described in further detail below.
In order to solve the defects of insufficient research on the dominant oscillation mode based on wide-area measurement information and the defect of an empirical mode decomposition method in the background art, the embodiment of the invention provides a dynamic stability collaborative identification method for an electric power system by means of the application of a CEEMDAN algorithm and a Hilbert-Huang transform algorithm in electric power system oscillation mode identification, so as to realize collaborative identification of the dominant oscillation mode of the electric power system and provide power grid operation state information for an electric power system operation dispatcher so as to improve the safe and stable operation of the electric power system.
Example 1:
the invention discloses a CEEMDAN algorithm-based power system low-frequency oscillation mode identification method, as shown in FIG. 9, FIG. 9 is a flow chart of the invention, including the following steps:
an application of CEEMDAN and Hilbert-Huang transform algorithm based on electric power system oscillation mode identification is disclosed, the calculation method comprises the following steps:
101, acquiring measured data of an electric power system by using a data acquisition function of a wide area measurement system, decomposing each group of original low-frequency oscillation signals into the sum of a plurality of Intrinsic Mode Functions (IMFs) by using a Complete integrated Empirical Mode with Adaptive Noise Decomposition (CEEMDAN) algorithm, wherein each IMF component represents an oscillation Mode;
102, preprocessing each IMF component, calculating the energy size and the energy weight of each IMF component by using a Teager energy operator, wherein in the generator rotor angular oscillation signal, the relative energy of the dominant oscillation mode is larger, and one or more modes with larger relative energy are respectively identified.
103, identifying the oscillation frequency and the damping ratio of the dominant oscillation mode by using a Hilbert-Huang transform algorithm, comparing a calculation result with a characteristic value method, and verifying the validity of the calculation result;
in summary, the embodiment of the present invention perfects and expands the dynamic evaluation direction and content of the power system based on the wide area measurement information through the steps 101 to 103, and realizes the dynamic stable evaluation of the power system based on data completely.
Example 2:
the scheme in embodiment 1 is further described below with reference to specific calculation formulas and examples, and is described in detail in the following description:
201: acquiring state measurement information of the power system from the wide area measurement system, and standardizing the state measurement information;
the normalization is a process of calculating a standard score, which is a process of dividing the difference between a number and a mean by a standard deviation. The standardized data is simple and convenient to compare, the relation between the data standard deviations can be fully displayed, and the original information of the data is reserved.
202: the standardized generator rotor angle signal is used as an input signal to be identified, and the oscillation signal is decomposed by using the complete set of adaptive noise and empirical mode decomposition, so that each oscillation signal is decomposed into a plurality of components which can represent different oscillation modes and are called eigenmode function components, namely IMF components. The CEEMDAN algorithm for the complete set empirical mode decomposition of the adaptive noise comprises the following specific steps:
the specific implementation process of the scheme is as follows:
1) adding a Gaussian white noise pair into an original signal x (t) to form a new signal:
x(t)+ε0ωi(t) (1)
ωi(t) (where i ═ 1, …, NR) is the white gaussian noise achieved differently, NR is the number of times the white gaussian noise is achieved, ε0Is the standard deviation of the noise. And respectively carrying out NR times of classical modal decomposition on the new signals to obtain NR number ofFirst order IMF component:
x(t)+ε0ωi(t)＝imf1i(t)+r1i(t) (2)
x(t)+ε0ωi(t) the i-th white Gaussian noise added to the original signal x (t), imf1i(t) is a number represented by the formula0ωi(t) the first IMF component after empirical mode decomposition, r1i(t) is the first residual component, for NR term imf1i(t) averaging to obtain final first-order IMF components
In the above formula:
is the mean of the first order IMF components; NR is the number of times of realizing Gaussian white noise and is also subjected to NR empirical mode decomposition; imf1i(t) is a number represented by the formula0ωi(t) the first IMF component after empirical mode decomposition, x (t) being the original signal,
is the mean of the first remaining component.
As can be seen from the equations (2) and (3), the average value is calculated by white Gaussian noise ε0ωi(t) interaction, resulting in
The noise is greatly weakened, and meanwhile, the problem of mode confusion is effectively restrained.
2) Then, the final first-order residual component is found by using the following equations (2) and (3):
in the above formula:
is the first order residual component;
is the mean of the first order IMF components.
Decomposing random white Gaussian noise by using a classical modal decomposition algorithm to obtain a first-order IMF component epsilon1E1ωi(t) is added to r1(t), constituting a new signal, namely:
r1(t)+ε1E1ωi(t)＝imf2i(t)+r2i(t) (5)
in the above formula: r is1(t) is the first order residual component; epsilon1The standard deviation of the first-order IMF component noise is obtained by decomposing random white Gaussian noise by using a classical modal decomposition algorithm; e1Being the amplitude of the noise, omegai(t) (where i ═ 1, …, NR) is white gaussian noise of the first order IMF component.
For NR imf2i(t) averaging to obtain final second-order IMF component
In the above formula:
is a second order IMF component; NR is the number of times of realizing Gaussian white noise and is also subjected to NR empirical mode decomposition; imf2i(t) a second-order IMF component obtained by the ith empirical mode decomposition; r is1(t) is the first order residual component; r is2i(t) is a second-order residual component obtained by the ith empirical mode decomposition;
3) according to the above, the k-th residual component can be:
wherein k is 2,3, L rk(t) is the kth order residual component; r isk-1(t) is the k-1 th order residual component; NR is the number of times of realizing Gaussian white noise and is also subjected to NR empirical mode decomposition; r is a radical of hydrogenki(t) is the kth order residual component obtained by the ith empirical mode decomposition;
decomposing random white Gaussian noise by using a classical modal decomposition algorithm, and obtaining a k-order IMF component epsilonkEkωi(t) is added to rk(t), constituting a new signal, namely:
rk(t)+εkEkωi(t) (8)
wherein k is 2,3, L, rk(t) is the kth order residual component; epsilonkThe standard deviation of the k-th order IMF component noise is obtained by decomposing random white Gaussian noise by using a classical modal decomposition algorithm; ekBeing the amplitude of the noise, omegai(t) (where i ═ 1, …, NR) is white gaussian noise of the k-th order IMF component.
And respectively carrying out NR times of multivariate classical modal decomposition on the new signals to obtain NR (k + 1) -th order IMF component IMF(k+1)i(t)：
rk(t)+εkEkωi(t)＝imf(k+1)i(t)+r(k+1)i(t) (9)
In the above formula: r isk(t) is the kth order residual component; epsilonkThe standard deviation of the k-th order IMF component noise is obtained by decomposing random white Gaussian noise by using a classical modal decomposition algorithm; ekBeing the amplitude of the noise, omegai(t) (where i ═ 1, …, NR) is white gaussian noise of the k-th order IMF component; IMF component IMF of order k +1(k+1)i(t)；r(k+1)i(t) is the k +1 th order residual component.
Wherein E iskωiAnd (t) carrying out classical modal decomposition on the random white Gaussian noise signal to obtain a k-order IMF component. For NR term imf(k+1)i(t) averaging to obtain the final (k + 1) th order IMF component
In the above formula:
is IMF component of k +1 order; NR is the number of times of realizing Gaussian white noise and is also subjected to NR empirical mode decomposition; imf(k+1)i(t) is the (k + 1) th order IMF component obtained by the ith empirical mode decomposition; r isk(t) is the kth order residual component; r is(k+1)iAnd (t) is the k +1 th order residual component obtained by the ith empirical mode decomposition.
The remaining components are solved again:
in the above formula: r is a radical of hydrogenk+1(t) is the k +1 th order residual component; r isk(t) is the kth order residual component;
is IMF component of k +1 order; NR is the number of times of realizing Gaussian white noise and is also subjected to NR empirical mode decomposition; r is a radical of hydrogen(k+1)iAnd (t) is the k +1 th order residual component obtained by the ith empirical mode decomposition.
4) Repeating the operation step (3) until the residual component rk(t) the number of extreme points is less than 2, P IMF components are finally obtained, and the CEEMDAN finally obtains a residual error r (t), wherein P is the total number of IMFs, and r (t) is:
the original signal x (t) can be expressed as:
in the above formula: x (t) is the original signal; imfk(t) is the k-th order IMF component; r (t) is the residual component.
203, the relative rotor angle contains oscillation information, the relative rotor angle information can be decomposed into P IMF components through a CEEMDAN algorithm, each IMF component represents an oscillation mode, only the dominant oscillation mode in all the modes has an influence on the stability of the system, and the relative energy of each IMF is identified by using a Teager energy operator single channel as a basis for identifying the dominant oscillation mode, which is specifically as follows:
for an original signal x (t) with amplitude A and frequency f, when the sampling frequency is fsWhen this is done, a discrete form of the following formula can be obtained:
xn＝Acos(nΩ+φ) (14)
in the above formula, Ω is 2 π f/fsPhi is the initial phase angle. A, Ω, φ is obtained by constructing a system of equations, as follows:
according to the above formula, obtain
In the above formula: x is the number ofnIs a discrete form expression of the original signal x (t) and is also the value of the original signal x (t) of the nth sampling point; a is the amplitude of x (t); n is the nth sampling point of the IMF component; 2 pi f/fsF is the frequency of the original signal x (t); f. ofsIs the sampling frequency; phi is the initial phase angle.
When Ω → 0, sin (Ω) ≈ Ω according to the infinitesimal equivalence, and when the sampling frequency fsWhen large enough, the above formula can be written as
In the above formula: 2 pi f/fsF is the frequency of the original signal x (t); f. ofsIs the sampling frequency; a is the amplitude of x (t); x is the number ofnIs the value of the original signal x (t) of the nth sampling point; x is the number ofn-1Is the value of the original signal x (t) of the (n-1) th sampling point; x is the number ofn+1Is the value of the original signal x (t) at the (n + 1) th sampling point.
Then define signal xnThe Teager Energy Operator (TEO) of
ψ＝A2sin2(Ω)≈A2Ω2 (18)
In the above formula: psi is the IMF component xnThe energy of (a); omega-2 pi f/fsF is the frequency of the original signal x (t); f. ofsIs the sampling frequency; a is the amplitude of x (t).
Wherein psi is TEO. And as can be seen from the above formula, the TEO has the following characteristics:
(1) symmetry. I.e. whether n or-n, xnOr x-nThe psi value is unchanged;
(2) the psi value is independent of the initial phase angle phi;
(3) the calculation needs fewer initial values, and only needs three connected sampling points in the signal.
The value ψ of each sample point is obtained from equation (18) by using the IMF component decomposed from the rotor angle signal. I.e. the psi values of the sample points of each IMF component are added, so that the energy sum of each IMF component can be identified. And if the sum of the IMF energy is obviously larger than the sum of other IMF energy, the required IMF component is obtained, and the rest IMF components are removed.
Calculating the energy weight of each IMF component, keeping the IMF component with high energy weight as a dominant oscillation mode, and removing the IMF component with low energy weight as a non-dominant oscillation mode;
where e (i) represents the energy of the ith IMF component, m represents the number of IMF components per channel, and m (i) represents the energy weight of the ith IMF component in the channel. The larger the M value is, the mode corresponding to the ith IMF component is the dominant oscillation mode, the dominant oscillation modes of all channels are calculated, and oscillation parameters contained in all the dominant oscillation modes are calculated.
203 identifying the instantaneous oscillation frequency and the instantaneous damping ratio of the dominant oscillation mode through a Hilbert-Huang transform algorithm, and summing discrete points to obtain the oscillation frequency and the damping ratio, wherein the method specifically comprises the following steps:
for each time-varying signal x (t), i.e. an IMF signal containing a dominant oscillation mode, it can be Hilbert transformed by:
in the above formula: (ii) a Tau is an integral variable; x (τ) is the same value as x (t), except for the integration variable.
Forming a complex conjugate pair of x (t) and y (t) to construct an analytic signal z (t):
z(t)＝x(t)+j y(t)＝a(t)ejφ(t) (21)
in the above formula: z (t) is a structure analytic signal; x (t) is the original signal; y (t) is the Hilbert transformed signal.
Where a (t) is the amplitude signal and φ (t) is a phase function:
by definition, the instantaneous angular frequency ω (t) is the derivative of the phase function φ (t) with respect to time:
in the above formula: ω (t) is the instantaneous angular frequency; phi (t) is a phase function.
Thus, the instantaneous frequency f (t) can be defined as:
in the above formula: f (t) is the instantaneous frequency; phi (t) is a phase function.
From equations (24) and (25), the final expression for the instantaneous frequency f (t) can be obtained as:
in the above formula: x (t) is the original signal, y (t) is the Hilbert transformed signal, x & (t) is the derivative of x (t), y & (t) is the derivative of y (t), and f (t) is the instantaneous frequency.
Knowledge of the instantaneous amplitude and instantaneous frequency of the signal enables further calculation of the instantaneous damping of the signal. Damping characteristics are another useful alternative to analyzing the local behavior of oscillations.
The time response equation z (t) is expressed in terms of eigenvalues according to equation ():
z(t)＝a(t)ejφ(t)＝Λeθ(t)+jφ(t) (27)
z(t)＝Λeλit＝Λe(σ±jω)t (28)
in the above formula: z (t) is a structure analysis signal; a (t) is an amplitude signal; phi (t) is a phase function; Λ is an amplitude constant; lambda [ alpha ]i(t)＝θ(t)+jφ(t)＝(σ±jω)t，λi(t) is a characteristic value that characterizes an oscillation mode; σ describes the damping characteristic of the power system oscillation as the real part of the eigenvalue, and w describes the frequency characteristic of the power system oscillation as the imaginary part of the eigenvalue.
Wherein a pair of conjugate eigenvalues represents an oscillation mode, and the damping ratio of the oscillation mode is calculated by the eigenvalues, wherein the formula for calculating the damping ratio is as follows:
in the above formula: sigma is a real part of the characteristic value to describe the damping characteristic of the oscillation of the power system, and w is an imaginary part of the characteristic value to describe the frequency characteristic of the oscillation of the power system; a (t) is an amplitude signal; θ (t) is a time characteristic of the amplitude signal a (t); ω (t) is the instantaneous angular frequency; phi (t) is a phase function; ζ is the damping ratio.
The above equation is a generalization of the damping ratio modal concept analysis of nonlinear non-stationary signals, and therefore, to ensure that the calculations satisfy the transient modal analysis. And to be able to obtain an estimate of the damping ratio for a series of signals, this can be achieved by an estimate of the average of the instantaneous damping ratios.
In the HHT algorithm, a hilbert transform is applied to each IMF to estimate its instantaneous frequency, instantaneous amplitude, and instantaneous damping ratio. Since the instantaneous frequency of a single frequency signal, a signal containing only one major frequency, is best defined, it is expected that each IMF is single frequency and meaningful. However, as indicated above, the IMF may contain modal aliasing, so adding white gaussian noise on top of the EMD algorithm improves the EMD algorithm, i.e. the full set of adaptive noise empirical mode decomposition cemsdan suppresses mode aliasing while the dominant oscillatory modes are screened out using the Teager energy operator as described above.
The method is combined with CEEMDAN decomposition, the original low-frequency oscillation signal is decomposed into a plurality of IMF components through CEEMDAN decomposition, then the energy of each sampling point of the IMF components is obtained through a Teager energy operator, the energy of each IMF component is summed, the IMF components with obviously-small energy are directly removed to avoid influencing the accuracy of the result, finally the instantaneous frequency and the instantaneous damping ratio of the real IMF components are obtained through Hilbert transformation, and the result is analyzed.
In summary, in the present invention, the dominant oscillation mode of the power system and the frequency and damping ratio of the dominant oscillation mode are identified through the steps 201 to 203, so that the dynamic mode identification and implementation of the power system based on the wide area measurement information are improved and expanded.
Example 3:
the following experimental data are combined to perform feasibility verification for the schemes of examples 1 and 2, as described in detail below:
in this example, a 16-machine 68 node system is taken as an example to perform simulation analysis, and a 16-machine 68 node test system is shown in fig. 1.
In the time domain simulation process, three-phase permanent short-circuit faults are set on the branches 46-49, the faults are set to be 0.1s, the near-end faults are cut off in 0.2s, and the far-end faults are cut off in 0.22 s. Selecting a No. 1 unit as a reference unit, taking a rotor angle signal of each generator unit relative to G1 as an input signal, generating 15 groups of rotor angle signals by 16 generators, wherein the sampling frequency is 0.01s, the rotor angle rocking curve of the generator after the power system is disturbed is shown in figure 2, and performing CEENDAN decomposition on the rotor angle information of each measuring channel respectively to obtain all intrinsic oscillation modes IMF.
Selecting a dominant oscillation mode: firstly, calculating instantaneous energy of IMF sampling points by using a Teager energy algorithm, accumulating and summing the instantaneous energy to obtain the energy of each IMF, and sequencing the IMFs from high to low according to the energy.
Identifying the relative energy of each IMF component by using a Teager energy operator, a series of discrete energies can be obtained for each channel and then summed, preserving the IMF components with relatively high energy values, as shown in FIG. 3, with the generator rotor angle signal of G2 relative to G1 as input, the waveforms of IMF1 components to IMF5 components are shown in FIGS. 4, 5, 6, 7 and 8, and identifying the IMF components with relatively high energy values by the Teager energy operator, as shown below:
wherein the energy values of IMF1 and IMF2 are significantly higher than those of other modes, and IMF1 and IMF2 are identified as dominant oscillation modes.
The oscillation frequency and the damping ratio of the dominant oscillation mode are identified by using the hilbert-yellow transformation, and the dominant oscillation mode is respectively subjected to the hilbert-yellow transformation identification by taking the rotor angle signal of G2 relative to G1 as an example: the oscillation frequency of the IMF1 is 1.986Hz, and the damping ratio is 2.8 percent; the oscillation frequency of IMF2 was 0.947Hz, and the damping ratio was 0.043%.
Taking the rotor angle signal of G3 relative to G1 as an example, utilizing a Teager energy operator to identify the dominant oscillation mode, and respectively carrying out Hilbert-Huang transform identification on the dominant oscillation mode: the oscillation frequency of IMF1 is 2.086Hz, and the damping ratio is 2.1%; the oscillation frequency of IMF2 was 1.006Hz, and the damping ratio was 4.8%. By analogy, the oscillation frequency and the damping ratio of all the generators are estimated, theoretically, the identification result of the dominant oscillation mode of the system is unique, but in an actual system, the identification result of each measurement channel is different due to the observability of the oscillation mode, a common processing method is to identify the measurement information of a plurality of channels by adopting a single-channel identification method, then the average value of the identification result is taken to represent the dynamic stability of the system, and the estimation result is shown in table 2.
To verify the correctness of the identification result in the embodiment of the present invention, the dominant oscillation mode of the system identified by the characteristic value analysis method is adopted in table 3, and the dominant oscillation mode identified by the method in the comparison table can be known as follows: the embodiment of the invention can accurately identify the dominant oscillation mode of the system, and verify the correctness of identifying the dominant oscillation mode of the system.
In conclusion, the CEEMDAN is expanded on the basis of EMD, the Teager energy operator is used for identifying the dominant oscillation mode, and then the oscillation frequency and the damping ratio are calculated, so that the stability mechanism is better analyzed, and the quantitative calculation of a complex system is realized.
TABLE 1 energies of all oscillation modes
| E1 | E2 | E3 | E4 | E5 |
| 2.8e-05 | 3.5e-06 | 1.2e-07 | 2.6e-09 | 1.5e-10 |
TABLE 2 estimation of dominant oscillation mode frequency to damping ratio
TABLE 3 comparison of the recognition results of different methods
Claims (3)
1. A method for identifying a low-frequency oscillation mode of a power system based on a CEEMDAN algorithm is characterized by comprising the following steps: 1) acquiring measured data of a power system by using a data acquisition function of a wide area measurement system, decomposing a Complete integrated Empirical Mode Decomposition with an Adaptive Noise of each group of original low-frequency oscillation signals into a plurality of Intrinsic Mode functions according to a CEEMDAN algorithm, and the sum of IMFs, wherein each IMF component represents an oscillation Mode;
the CEEMDAN algorithm comprises the steps of:
adding different types of random white Gaussian noise into original signal x (t), x (t) + epsilon0ωi(t)＝xi(t), where I ═ 1, … I number of times to achieve white gaussian noise, ∈0Is the standard deviation of Gaussian white noise, ωi(t) is white gaussian noise of different types; for new signal xi(t) performing empirical mode decomposition, averaging I empirical mode decomposition results, and extracting a first intrinsic mode function component IMF1；
Secondly, decomposing different types of random Gaussian white noise empirical modes to obtain epsilon1Eiωi(t), wherein: epsilon1Is a standard deviation of, EiIs the first order IMF component, ωi(t) is white gaussian noise of different types;
subtracting IMF from original signal1Then adding epsilon1E1ωi(t) forming a new signal, performing empirical mode decomposition on the new signal, averaging the results of the empirical mode decomposition, and extracting a second intrinsic mode function component IMF2；
Fourthly, analogizing in sequence until an iteration stop criterion is met, wherein each IMF component corresponds to an oscillation mode;
2) preprocessing each IMF component, calculating the energy size and the energy weight of each IMF component by using a Teager energy operator, wherein in the generator rotor angular oscillation signal, the relative energy of a dominant oscillation mode is larger, and one or more modes with larger relative energy are respectively identified;
3) identifying the oscillation frequency and the damping ratio of the dominant oscillation mode by using a Hilbert-Huang transform algorithm, comparing a calculation result with a characteristic value method, and verifying the effectiveness of the calculation result;
the step 1 of acquiring the measured data of the power system by using the data acquisition function of the wide area measurement system includes:
acquiring actual measurement data of the power system by using a data acquisition method of a Phasor Measurement Unit (PMU) and a Wide Area Measurement System (WAMS), wherein the actual measurement data comprises a rotor angle signal of each generator or active power on a connecting line as a signal to be identified;
further comprising: acquiring state measurement information of the power system from a wide area measurement system, wherein the state measurement information comprises rotor angle signals and active power on a connecting line, and standardizing the state measurement information;
in the step 2, the energy size and the energy weight of each IMF component are calculated by using a Teager energy operator, namely the energy of all IMF components is calculated by using the Teager energy operator to obtain the IMF component energy, and the method comprises the following steps: calculating the energy size and the energy weight of each IMF component by using a Teager energy operator,
for continuous signal xnThe formula for calculating the energy is as follows:
summing the energy values of all discrete sampling points of the IMF component to serve as the energy of the IMF component;
in the step 3, the oscillation frequency and the damping ratio of the dominant oscillation mode are identified by using a Hilbert-Huang transform algorithm: subjecting the original signal x (t) to a Hilbert transform;
comparing the calculation result with a characteristic value method, and verifying the validity of the calculation result;
1) for a continuous signal x (t), an intrinsic mode function IMF component, it is Hilbert transformed:
wherein P represents the Cauchy prime value integral
2) Combining x (t) and y (t) into a pair of complex numbers, the following signal z (t) is formed:
z(t)＝x(t)+jy(t)＝A(t)ejφ(t)
wherein A (t) is a magnitude signal, phi (t) is a phase function, x (t) is an original signal, y (t) is a Hilbert transformed signal, z (t) is a complex number formed by x (t) and y (t), and a modulus and a phase angle thereof represent the magnitude and the phase of the signal;
identifying the oscillation frequency of the dominant oscillation mode in the step 3;
by definition, the instantaneous angular frequency ω (t) is the derivative of the phase Φ (t) with respect to time, and the instantaneous frequency f (t) ═ ω (t)/2 π, the expression for the instantaneous frequency is obtained:
in the above formula: phi (t) is the phase function, f (t) is the instantaneous frequency;
the final result is represented by x (t) and y (t) and their derivatives, and the final expression for the instantaneous frequency f (t) can be obtained as:
in the above formula: x (t) is the original signal, y (t) is the Hilbert transformed signal, x & (t) is the derivative of x (t), y & (t) is the derivative of y (t), and f (t) is the instantaneous frequency; averaging the instantaneous frequency to obtain the oscillation frequency of the mode;
identifying the damping ratio of the dominant oscillation mode in the step 3:
in the above formula: zeta is damping ratio, dA (t)/dt is amplitude functionThe derivative of (a) of (b),
is the derivative of the phase function; averaging the instantaneous damping ratio to obtain the damping ratio in the mode;
the characteristic value method in the step 3 is to calculate the oscillation frequency and the damping ratio of the same system by using the characteristic value method;
identifying the oscillation frequency and the damping ratio of the dominant oscillation mode by using a Hilbert-Huang transform algorithm, comparing a calculation result with a characteristic value method, and verifying the effectiveness of the calculation result;
solving a linear equation of the system to obtain a state matrix A, decomposing the characteristic of the matrix A, and solving a characteristic value lambdai＝σi+jωiEach pair of conjugate complex roots corresponds to an oscillation mode, and the oscillation frequency and the damping ratio of the mode are calculated as follows:
in the above formula: lambda [ alpha ]iIs a conjugate eigenvalue; sigmaiRepresenting the damping characteristic of the oscillation as the real part of the characteristic value; omegaiThe characteristic value is an imaginary part of the characteristic value and is used for representing the frequency characteristic of the oscillation, xi is a damping ratio, and f is the oscillation frequency.
2. The method for oscillation pattern recognition based on CEEMDAN algorithm as claimed in claim 1, wherein the energy magnitude and energy weight of each IMF component are calculated in step 1:
the method comprises the steps of calculating the energy weight of all IMF components, keeping the IMF components with high energy weight as a dominant oscillation mode, and removing the IMF components with low energy weight as a non-dominant oscillation mode.
3. The method for identifying the low-frequency oscillation mode of the power system based on the CEEMDAN algorithm as claimed in claim 1, wherein the step 3 of identifying the oscillation frequency of the dominant oscillation mode by using the Hilbert-Huang transform algorithm is compared with a characteristic value method; theoretically, if N generators exist in the system, N-1 low-frequency oscillation modes should exist; however, the oscillation parameters of the dominant oscillation modes decomposed by the rotor angle signals of different generators are the same or similar, and the dominant oscillation modes belong to the same oscillation mode;
classifying modes having the same oscillation parameters as one oscillation mode; and comparing with a characteristic value method to verify the effectiveness of the algorithm.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811647265.5A CN109638862B (en) | 2018-12-30 | 2018-12-30 | CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811647265.5A CN109638862B (en) | 2018-12-30 | 2018-12-30 | CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109638862A CN109638862A (en) | 2019-04-16 |
| CN109638862B true CN109638862B (en) | 2022-07-05 |
Family
ID=66055189
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811647265.5A Active CN109638862B (en) | 2018-12-30 | 2018-12-30 | CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109638862B (en) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110096998B (en) * | 2019-04-29 | 2023-05-12 | 天津大学 | Decomposition method for extracting machined surface morphology features |
| CN110112757B (en) * | 2019-05-17 | 2022-04-12 | 福州大学 | Low-frequency oscillation analysis method based on SURE wavelet denoising and improved HHT |
| CN110311392B (en) * | 2019-07-21 | 2022-05-13 | 东北电力大学 | Electric power system oscillation mode and mode identification method based on SPDMD |
| CN112329535B (en) * | 2020-09-29 | 2023-03-24 | 国网四川省电力公司经济技术研究院 | CNN-based quick identification method for low-frequency oscillation modal characteristics of power system |
| CN112395947A (en) * | 2020-10-20 | 2021-02-23 | 太原理工大学 | Non-contact type transformer local short circuit detection method and detection device |
| CN112629862A (en) * | 2020-11-04 | 2021-04-09 | 南京工业大学 | Rolling bearing weak fault diagnosis method based on bistable stochastic resonance and CEEMDAN-TEO |
| CN112505477B (en) * | 2020-11-16 | 2023-12-08 | 广东电网有限责任公司广州供电局 | Disturbance initial judgment method based on synchronous phasor data of power distribution network |
| CN113010844B (en) * | 2021-03-09 | 2022-11-11 | 东北电力大学 | Participation factor calculation method based on subspace dynamic mode decomposition |
| CN113654750B (en) * | 2021-08-13 | 2022-05-03 | 华北电力大学（保定） | Vibration feature extraction method for circuit breaker operating mechanism |
| CN114050568A (en) * | 2021-11-09 | 2022-02-15 | 东北电力大学 | Forced oscillation source positioning method and device based on multivariate empirical mode decomposition |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5983162A (en) * | 1996-08-12 | 1999-11-09 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computer implemented empirical mode decomposition method, apparatus and article of manufacture |
| CN103795073A (en) * | 2014-02-27 | 2014-05-14 | 武汉大学 | Working condition On-line identification and early warning method of low frequency oscillation leading model of electric power system |
| CN104753075A (en) * | 2015-03-19 | 2015-07-01 | 中国农业大学 | Identifying method and device of leading oscillating mode of interconnected electric power system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009048964A1 (en) * | 2007-10-09 | 2009-04-16 | Schweitzer Engineering Laboratories, Inc. | Real-time power system oscillation detection using modal analysis |
-
2018
- 2018-12-30 CN CN201811647265.5A patent/CN109638862B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5983162A (en) * | 1996-08-12 | 1999-11-09 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Computer implemented empirical mode decomposition method, apparatus and article of manufacture |
| CN103795073A (en) * | 2014-02-27 | 2014-05-14 | 武汉大学 | Working condition On-line identification and early warning method of low frequency oscillation leading model of electric power system |
| CN104753075A (en) * | 2015-03-19 | 2015-07-01 | 中国农业大学 | Identifying method and device of leading oscillating mode of interconnected electric power system |
Non-Patent Citations (1)
| Title |
|---|
| HHT在电力系统低频振荡模态参数提取中的应用;李天云等;《中国电机工程学报》;20071005(第28期);第79-83页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109638862A (en) | 2019-04-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109638862B (en) | CEEMDAN algorithm-based electric power system low-frequency oscillation mode identification method | |
| CN106845010B (en) | Low-frequency oscillation dominant mode identification method based on improved SVD noise reduction and Prony | |
| Sanchez-Gasca et al. | Computation of power system low-order models from time domain simulations using a Hankel matrix | |
| CN111289796B (en) | Method for detecting subsynchronous oscillation of high-proportion renewable energy power system | |
| CN109787250B (en) | Power system low-frequency oscillation mode identification method based on multivariate empirical mode decomposition | |
| CN109659957B (en) | APIT-MEMD-based power system low-frequency oscillation mode identification method | |
| Tashman et al. | Modal energy trending for ringdown analysis in power systems using synchrophasors | |
| Swiatek et al. | Power system harmonic estimation using neural networks | |
| Bank et al. | Generator trip identification using wide-area measurements and historical data analysis | |
| Yang et al. | Extracting inter-area oscillation modes using local measurements and data-driven stochastic subspace technique | |
| Papadopoulos et al. | Online parameter identification and generic modeling derivation of a dynamic load model in distribution grids | |
| Betancourt et al. | A spatio-temporal processing Padé approach for visualizing harmonic distortion propagation on electrical networks | |
| Zarei et al. | Broken rotor bars detection via Park's vector approach based on ANFIS | |
| Shetty et al. | Low frequency oscillation detection in smart power system using refined prony analysis for optimal allocation of supplementary modulation controller | |
| Kumar et al. | S-Transform based detection of multiple and multistage power quality disturbances | |
| Shetty et al. | Performance analysis of data-driven techniques for detection and identification of low frequency oscillations in multimachine power system | |
| Zhao et al. | The condition monitoring of wind turbine gearbox based on cointegration | |
| CN107345983A (en) | Multi-harmonic Sources system harmonicses transmitting appraisal procedure based on subharmonic source correlation | |
| Jakpattanajit et al. | On-line estimation of power system low frequency oscillatory modes in large power systems | |
| Kuwałek et al. | Multi-point method using effective demodulation and decomposition techniques allowing identification of disturbing loads in power grids | |
| Zhang et al. | Fault discrimination using synchronized sequence measurements under strong white Gaussian noise background | |
| WO2019244422A1 (en) | Power system operation assisting device and method, and oscillation suppression system | |
| Avdakovic et al. | Identification of low frequency oscillations in power system | |
| Jankowski et al. | The study of chaotic behaviour in the marine power system containing nonlinear varying load | |
| Khodaparast et al. | Implementation of the neural network for tracing of spot welders |
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 |