Direct current system relay protection method and device based on equivalent filter parameter identification

Abstract

This invention relates to a method and apparatus for DC system relay protection based on equivalent filter parameter identification. The method includes: constructing a lossy Begeron model to extract fault characteristics outside the DC system area; constructing an IIR filter model and identifying its parameters; calculating the confidence interval of the IIR filter model parameter vector and obtaining the protection interval by upper and lower envelopes; calculating the start-up decision value vector and determining whether a fault has occurred based on the start-up decision value vector; using first-order difference to back-calculate the standard unit impulse input sequence and the original unit impulse response sequence based on the discrete-time linear time-invariant system theory; calculating the real-time IIR filter model parameter vector based on the standard unit impulse input sequence and the original unit impulse response sequence; and determining the fault interval of the DC system based on the real-time IIR filter model parameter vector and the protection interval. This invention improves the reliability, speed, sensitivity, and selectivity of the protection.

Description

Technical Field

[0001] This invention relates to the field of power system transmission line relay protection, specifically to a DC system relay protection method and device based on equivalent filter parameter identification. Background Technology

[0002] High-voltage direct current (HVDC) transmission systems offer numerous advantages, including long-distance and high-capacity power transmission, making them highly suitable for my country's resource-rich environment and enabling optimized energy allocation. Utilizing increasingly sophisticated relay protection devices for rapid response in fault conditions, coupled with a reliable relay protection and control system to ensure accurate fault identification and isolation, the systems can fully guarantee the safe and stable operation of HVDC transmission corridors. If line protection fails to respond correctly (failure to operate or maloperation), it can not only cause HVDC system shutdown but also easily trigger large-scale power loss accidents leading to forced outages, severely impacting the safety and stability of the power system.

[0003] Currently, the academic and engineering communities have proposed many fault identification and protection schemes for DC lines, which can extract transient characteristics after a fault occurs and use them for fault location and isolation. However, the following problems still exist: (1) Traditional single-ended or double-ended protection schemes: including differential undervoltage protection, low voltage protection and longitudinal current differential protection, etc. These schemes have certain effects in actual complex operating environments, but they are generally difficult to achieve good selective fault isolation. Moreover, longitudinal current differential protection has to avoid abnormal differential current in fault scenarios outside the zone, so it has to sacrifice speed and lengthen the response speed. It can only be used as backup protection and is difficult to meet the requirements of rapid isolation of line faults.

[0004] (2) Traveling wave protection scheme: Extract the characteristics of the reverse traveling wave propagating from the fault point to both ends of the line, thereby accurately locating the location of the fault and triggering the control system to respond to the fault. It has high speed, but the setting depends on simulation and has limited sensitivity (i.e., poor resistance to transition resistance). Some high-speed and sensitive protection schemes have weak performance in terms of noise interference.

[0005] (3) Data-driven protection scheme: Its basic principle is similar to transient quantity feature extraction. The difference is that it uses physical structure to derive analytical expressions or uses machine learning and deep learning networks to automatically extract and fuse features. It is suitable for fault type and location determination under complex working conditions. This scheme has broad application prospects in improving fault identification accuracy. However, the analytical method that relies on mathematical models will lead to a decrease in identification stability when the power grid parameters change. When the data-driven model increases the amount of data and the system structure changes, its feature extraction efficiency, convergence speed and generalization ability are difficult to guarantee. Summary of the Invention

[0006] This invention provides a DC system relay protection method and system based on equivalent filter parameter identification to solve at least one of the above-mentioned technical problems.

[0007] The technical solution of this invention to solve the above-mentioned technical problems is as follows: A DC system relay protection method based on equivalent filter parameter identification, comprising: Step 1: Construct a lossy Begeron model of the DC system and set multiple sets of different input signals according to the most severe external fault conditions; use the lossy Begeron model to calculate each set of input signals to extract the external fault characteristics of the DC system and obtain the corresponding output signals. Step 2: Treat the DC system as a black box network to be identified, construct the IIR filter model of the black box network, and use the Steiglitz-McBride method to identify the parameters of the IIR filter model using multiple sets of input signals and their corresponding output signals to obtain multiple sets of IIR filter model parameter vectors. Step 3: Calculate the confidence interval of each group of IIR filter model parameter vectors, and based on the confidence intervals of multiple groups of IIR filter model parameter vectors, perform upper and lower envelope operations on the IIR filter model parameter vectors to obtain the guard interval; Step 4: Real-time sampling of the positive and negative voltages and currents at the DC system protection installation point, followed by mode decoupling to construct a state matrix. Discrete energy calculation and first-order differential gradient calculation are then performed on the state matrix to construct a gradient-enhanced energy matrix. Weighted coefficient vectors for line, zero-mode voltage, and current are defined. A start-up decision value vector is calculated based on the gradient-enhanced energy matrix and the weighted coefficient vectors. The start-up decision value vector is used to determine whether a fault has occurred in the DC system, and the time of fault occurrence is recorded when a fault is determined to have occurred in the DC system. Step 5: Set up the fault experiment and obtain the experimental observation signal by inputting the experimental input vector into the DC system; use the normalized least mean square algorithm to transform the experimental observation signal into a fixed parameter system based on the frequency-varying system, and iteratively update the IIR filter model parameter vector; then predict the output of the DC system based on the iteratively updated IIR filter model parameter vector to obtain the final predicted output vector; based on the discrete-time linear time-invariant system theory, use the experimental input vector and the final predicted output vector to back-calculate using the first-order difference to obtain the standard unit impulse input sequence and the original unit impulse response sequence. Step 6: Based on the Steiglitz-McBride method, calculate the parameters of the IIR filter model according to the standard unit impulse input sequence and the original unit impulse response sequence to obtain the real-time IIR filter model parameter vector; based on the real-time IIR filter model parameter vector and the protection interval, determine the interval in which the DC system has a fault, and then execute the corresponding protection.

[0008] Based on the above-mentioned DC system relay protection method based on equivalent filter parameter identification, the present invention also provides a DC system relay protection device based on equivalent filter parameter identification.

[0009] A DC system relay protection device based on equivalent filter parameter identification includes a processor, a memory, and a computer program stored in the memory. When the computer program is executed by the processor, it implements the DC system relay protection method based on equivalent filter parameter identification as described above.

[0010] The beneficial effects of this invention are: the DC system relay protection method and device based on equivalent filter parameter identification improves the reliability, speed, sensitivity, and selectivity of the protection. Specifically: In terms of reliability: First, this invention does not rely on a two-end communication system for data synchronization. For high-voltage DC transmission lines that have experienced faults, the implementation of this invention can complete the start-up decision based solely on local positive and negative voltage and current at the protection installation point. This is protection based on local quantities, and its implementation method is relatively reliable. Second, this invention only uses experimental observation signals during operation for online adaptive learning. Furthermore, before fault identification, it uses a normalized least mean square algorithm to iteratively update the IIR filter model parameters online, eliminating the risk of parameter fluctuations caused by changes in power grid operation mode or topology reconfiguration, and significantly improving the adaptive capability and operational reliability of the protection system under complex operating conditions. In terms of speed: In the fault feature extraction and start-up decision stages, this invention introduces discrete energy calculation and first-order differential gradient calculation to construct a gradient-enhanced energy matrix, which can reliably trigger start-up within the first few milliseconds of a fault occurrence using only a very short time data window; In addition, this invention cleverly transforms the complex frequency-varying nonlinear DC system into a discrete-time linear time-invariant system with fixed parameters, and uses a normalized least mean square algorithm with extremely low computational complexity for iteration; This "dimensionality reduction and linearization" process minimizes the computing power consumption of the relay protection's underlying microprocessor, improving the protection's operating speed from an algorithmic perspective; In terms of sensitivity and selectivity: This invention calculates the out-of-zone fault conditions of DC systems based on the lossy Begeron model and constructs the IIR filter model of the black box to be identified. The upper and lower envelopes of the confidence interval generated by parameter identification are used as protection intervals to measure whether the real-time parameter vector of the system exceeds the limit (whether an in-zone fault occurs). Compared with the traditional analytical method or pure data-driven method that strongly depends on the accuracy of the model, it has a high resistance to transition resistance to cope with in-zone faults and reliably does not operate in the out-of-zone fault scenario. Attached Figure Description

[0011] Figure 1 This is a flowchart of the DC system relay protection method based on equivalent filter parameter identification according to the present invention. Figure 2 This is a schematic diagram of the three-terminal DC system in the embodiment; Figure 3 This is a schematic diagram of the lossy Begeron model constructed in the embodiment; Figure 4 This is a unit impulse response curve under Experiment 1 in the embodiment. Figure 5 This is a unit impulse response curve for Experiment 2 in the embodiment. Figure 6 This is a unit impulse response curve for Experiment 3 in the embodiment; Figure 7 This is a unit impulse response curve for Experiment 4 in the embodiment; Figure 8 This is a unit impulse response curve for Experiment 5 in the embodiment; Figure 9 This is a unit impulse response curve for Experiment 6 in the embodiment; Figure 10 This is a schematic diagram of the 95% confidence interval of the IIR filter model parameter vector in the embodiment; Figure 11 This is a schematic diagram of the protection zone in the embodiment; Figure 12 This is a schematic diagram of the startup decision value and the startup dynamic decision threshold curve in the embodiment; Figure 13 This is a schematic diagram of the fitting curve for a metallic single-pole grounding fault within the zone in the embodiment. Figure 14 This is a schematic diagram illustrating the fault determination process within the area in the embodiment. Figure 15 This is a schematic diagram illustrating the fault determination situation outside the zone in the embodiment. Detailed Implementation

[0012] The principles and features of the present invention are described below with reference to the accompanying drawings. The examples given are only for explaining the present invention and are not intended to limit the scope of the present invention.

[0013] like Figure 1 As shown, the DC system relay protection method based on equivalent filter parameter identification includes: Step 1: Construct a lossy Begeron model of the DC system and set multiple sets of different input signals according to the most severe external fault conditions; use the lossy Begeron model to calculate each set of input signals to extract the external fault characteristics of the DC system and obtain the corresponding output signals. Step 2: Treat the DC system as a black box network to be identified, construct the IIR filter model of the black box network, and use the Steiglitz-McBride method to identify the parameters of the IIR filter model using multiple sets of input signals and their corresponding output signals to obtain multiple sets of IIR filter model parameter vectors. Step 3: Calculate the confidence interval of each group of IIR filter model parameter vectors, and based on the confidence intervals of multiple groups of IIR filter model parameter vectors, perform upper and lower envelope operations on the multiple groups of IIR filter model parameter vectors to obtain the guard interval; Step 4: Real-time sampling of the positive and negative voltages and currents at the DC system protection installation point and modal decoupling to construct a state matrix; discrete energy calculation and first-order differential gradient calculation are performed on the state matrix to construct a gradient-enhanced energy matrix; weighted coefficient vectors for line, zero-mode voltage, and current are defined; and a start-up decision value vector is calculated based on the gradient-enhanced energy matrix and the weighted coefficient vectors. The start-up decision value vector is used to determine whether a fault has occurred in the DC system, and the time of fault occurrence is recorded when a fault is determined to have occurred in the DC system. Step 5: Set up the fault experiment and obtain the experimental observation signal by inputting the experimental input vector into the DC system; use the normalized least mean square algorithm to transform the experimental observation signal into a fixed parameter system based on the frequency-varying system, and iteratively update the IIR filter model parameter vector; then predict the output of the DC system based on the iteratively updated IIR filter model parameter vector to obtain the final predicted output vector; based on the discrete-time linear time-invariant system theory, use the experimental input vector and the final predicted output vector to back-calculate using the first-order difference to obtain the standard unit impulse input sequence and the original unit impulse response sequence. Based on the Steiglitz-McBride method, the parameters of the IIR filter model are calculated according to the standard unit impulse input sequence and the original unit impulse response sequence to obtain the real-time IIR filter model parameter vector; based on the real-time IIR filter model parameter vector and the protection interval, the fault interval of the DC system is determined, and then the corresponding protection is executed.

[0014] Figure 2This is a three-terminal flexible DC system. Sides A and B are MMC (Multi-Center) stations, and side C is an LCC (Limited-Center) station. The two ends of the two DC lines are connected to MMC-A and MMC-B, and MMC-B and LCC, respectively. The current measurement point is located at... , , , Set in Single-pole grounding or inter-pole faults can occur at any time. This embodiment uses... Figure 2 The three-terminal flexible DC system shown is used to illustrate the steps of the method of the present invention.

[0015] In step 1 of this embodiment, a lossy Bergeron model of the DC system at a fixed frequency is first constructed. The entire DC system is considered as a black-box network. Depending on the fault type, transition resistance, and specific frequency, multiple sets (e.g., ...) are configured. Groups of different off-site fault experiments: i.e., setting up A set of input signals (specifically, unit pulse excitation signals) are input to the DC system under different operating conditions with zero initial state. The unit pulse response signal of the DC system (which is the output signal) is calculated using the lossy Begeron model.

[0016] Preferably, in step 1, constructing a lossy Begeron model of the DC system specifically includes: Step 11: Equivalent the inductor in the DC system to a parallel connection of a resistor and a current source to construct the Bergeron model of the inductor.

[0017] Specifically, according to the definition of inductance L: (1) In the formula, , The nodes represent the two sides of the inductor L. , For nodes , exist Voltage to ground at any given moment. for From the node Flow to Node The current.

[0018] Perform the following steps on both sides of equation (1) arrive Integrals: (2) In the formula, for From the node Flow to Node The current.

[0019] when When the value is sufficiently small, the integral term in equation (2) is approximated by a trapezoidal rule: (3) Substituting equation (3) into equation (2), we get: (4) In the formula, we can know from the first term that Depend on and Therefore, it is obtained through superposition. It is the equivalent current source of inductor L. Let L be the equivalent resistance of the inductor.

[0020] From the first term of equation (4), we can obtain: (5) Substituting equation (5) into the second term of equation (4), we can obtain the equivalent current source of inductor L. Rewritten as: (6) After rewriting It can be removed from use. Perform calculations, instead of This makes calculations easier.

[0021] Step 12: Equivalently represent the capacitor in the DC system as a parallel connection of a resistor and a current source to construct a capacitor Bergeron model.

[0022] Specifically, assume that the nodes on both sides of capacitor C are respectively , The capacitance C is defined as: (7) Since capacitors and inductors are a pair of dual components, following the derivation process of the inductor Bergeron model in step 11, the capacitor Bergeron model can be obtained: (8) From this formula, we can see that Depend on and Therefore, it is obtained through superposition. This is the equivalent current source for capacitor C. Let be the equivalent resistance of capacitor C.

[0023] Step 13: In the modal domain, the fault point of the DC system is equivalent to a resistor and voltage source connected in series and grounded to construct the equivalent circuit of the line zero mode fault component.

[0024] Specifically, definition , , , These represent the positive and negative voltages and the positive and negative currents, respectively. Assume... A positive ground fault can occur at any time, and the fault point on the transmission line is... When a fault occurs, the effect of the fault point on the DC system can be equivalent to a voltage source (negative terminal grounded) with a series transition resistor drawn at the fault point. Therefore, at the fault point... The polar boundary conditions at point are: (9) In the formula, Fault point The transition resistance at the point, Fault point The positive terminal rated operating voltage, , Fault point The positive and negative currents at the point. Fault point The positive voltage at that point.

[0025] Since there is coupling between the bipolar poles, a transformation matrix is ​​introduced to decouple them: (10) definition , , , These represent the line-mode and zero-mode voltages and currents, respectively. The transformation relationship between the mode-domain components and the polar-domain physical quantities can be expressed as: (11) in , , , .

[0026] By performing the decoupling transformation of equation (11) on the fault boundary condition equation (9), we can obtain: (12) In the formula, , for Linear and zero-mode currents at the location , for Line-mode and zero-mode voltages at the location.

[0027] From equation (12), the equivalent network for the boundary conditions can be derived as follows: the zero-mode network and the line-mode network are connected through a voltage source. and resistance The series branches are connected in series, including the voltage source. The negative terminal is grounded, and .

[0028] Step 14: Construct a time-domain equivalent calculation model of the Bergeron transmission line considering frequency-dependent losses for the DC system.

[0029] Specifically, select the transmission line in the DC system. Its wave impedance is set to The line length is The transmission line connects the nodes. and nodes ,in , They represent time Side and voltage on the side, , They represent Time flows by Side and The current on the side is defined as having a positive direction from the node towards the interior of the line. According to the laws of mixed wave propagation, Side The mixed waves emitted at any time, after Arriving after time transmission delay Side. Among them For traveling waves in transmission lines One-way transmission time in This is the traveling wave speed. Due to transmission line losses, the wave attenuates during transmission. Definition side Voltage at time and current and side The electrical quantity relationships at any given time are as follows: First, define the transmission line. side Time and side Mixed wave variables at time and : (13) The two are linked through the attenuation coefficient: (14) in, For transmission lines The attenuation coefficient. In the frequency domain, this attenuation coefficient... The expression is: (15) In the formula, , , , Transmission lines The unit resistance, unit inductance, unit conductance, and unit capacitance. Because... , , , All are frequency-varying parameters. It cannot be directly represented using simple analytical mathematics.

[0030] Define amplitude attenuation coefficient and phase attenuation coefficient : (16) In the formula, , For a specific frequency.

[0031] To achieve time-domain calculations, the attenuation coefficient... Approximation using rational functions in the complex frequency domain: (17) In the formula, and These are the amplitude attenuation coefficient and the phase attenuation coefficient, respectively, and both are constants; , Substituting this approximation into the mixed-wave relation and performing a Laplace transform, we obtain the complex frequency domain expression: (18) In the formula, For transmission lines Mixed wave variables in the complex frequency domain, For transmission lines Mixed wave variables in the side complex frequency domain.

[0032] To solve the above equation, an intermediate equivalent computational circuit is constructed to characterize the transfer function relationship. This intermediate equivalent circuit consists of a resistor with a resistance of... The resistor is connected in series with an operational inductor, where the inductance of the operational inductor is [value missing]. The voltage source applied across the series branch is of the magnitude of (corresponding to the complex frequency domain) In this circuit, The voltage drop across the resistor corresponds to the required value. Time-domain discretization and final equivalent model derivation address the operational inductance in the aforementioned intermediate equivalent calculation circuit. Using the Bergeron model of the inductor element from step 11 above, the inductor is equivalent to a parallel circuit, which consists of a resistor with a resistance of... A resistor and a historical current source Composition, in which To calculate the step size, based on this discretized circuit, applying Kirchhoff's laws to the intermediate equivalent circuit, we can derive... and The time-domain recurrence relation: (19) In the formula, For transmission lines side Mixed wave variables at time Effect on The current generated across the resistor. Rearranging the above equations and combining them with the initially defined expression for the mixed wave, we can obtain the transmission line... The final time-domain equivalent computational model equation for the side is: (20) Therefore, the transmission line is in The side port can be equivalent to the Thevenin model: this model consists of a resistance value of The resistor is connected in series with a controlled voltage source, the voltage of which is [value missing]. Similarly, according to the duality principle, the transmission line in The equivalent model of the side port is also given by a resistance value. A resistor and a controlled voltage source The transmission lines are connected in series. This completes the construction of the equivalent time-domain calculation model for the Bergeron transmission line, which includes line loss characteristics.

[0033] Step 15: Based on the inductor Begeron model, capacitor Begeron model, equivalent circuit of line-zero mode fault component, and Begeron transmission line time-domain equivalent calculation model constructed in steps 11 to 14, set all equivalent power supplies to zero, and for solving the Begeron model with pulse signal as input excitation, replace the step signal of the fault component with a unit pulse signal to construct a lossy Begeron model of the DC system in zero initial state; the lossy Begeron model includes line-mode network and zero-mode network.

[0034] Specifically, assuming it occurs Figure 2 middle For a positive ground fault at a certain location, construct the corresponding lossy Bergeron model as follows: Figure 3 As shown (also known as the lossy Begeron equivalent circuit).

[0035] After the lossy Begeron model is constructed, multiple sets are set according to the most severe out-of-zone fault conditions (e.g., Groups) different experiments (sets) Each group of experiments uses a unit pulse excitation signal as its input signal. The response of each group of unit pulse excitation signals is calculated using a lossy Begeron model to extract the fault characteristics outside the DC system, that is, to calculate the current at the protection installation point (unit pulse response signal) when the DC system is faulty, thereby obtaining the corresponding output signal.

[0036] Preferably, the lossy Begeron model is used to calculate the external fault characteristics of the DC system for each group of input signals to obtain the corresponding output signal, specifically including step 16: For the line-mode network, the zero-mode network is disconnected and the fault component is grounded to construct the node admittance matrix and current source vector of the line-mode network; For zero-mode networks, the line-mode network is disconnected and the fault components are grounded to construct the node admittance matrix and current source vector of the zero-mode network; Construct the mutual admittance contribution matrix of the fault point in the online mode network and the zero mode network, and construct the node admittance matrix of the DC system based on the node admittance matrix of the online mode network, the node admittance matrix of the zero mode network and the mutual admittance contribution matrix. Construct the DC system current source vector based on the current source vector of the linear mode network and the current source vector of the zero mode network; Based on the nodal voltage method, a nodal voltage calculation model is constructed to calculate the voltage of each node in the DC system, according to the nodal admittance matrix of the DC system and the current source vector of the DC system. A calculation model for the positive current at the protection installation point in a DC system is constructed to calculate the positive current at the protection installation point at the time of a fault. The input signal is substituted into the node voltage calculation model and the positive current calculation model at the protection installation point to obtain the output signal.

[0037] Specifically, assuming in Figure 2 of A positive ground fault occurs at point A. Based on the circuit, the node voltage method is used to determine the fault. Figure 3 The voltages of each node (node ​​1 to node 20) are solved. For the line-mode network and fault component, the zero-mode network is disconnected and the fault component is grounded. From the lossy Begeron model, it can be seen that the DC system can be roughly divided into three parts, and the node admittance matrices are constructed from left to right: ;(twenty one) In the formula, The admittance is connected in parallel between node 1 and the left grounding point. The admittance is the series connection between node 1 and the right-side grounding point. The admittance is the parallel connection between node 1 and node 2. The admittance is the parallel connection between node 2 and node 3. and The admittance is the parallel connection between node 3 and node 4. and The admittance is the parallel connection between node 4 and the grounding point.

[0038] ;(twenty two) In the formula, The admittance connected in series between node 5 and the left grounding point (from Figure 3 It can be seen that the admittance connected in series between node 1 and the right grounding point, and the admittance connected in series between node 5 and the left grounding point, both adopt... express), The admittance is the series connection between node 5 and the right-side grounding point. The admittance is the parallel connection between nodes 5 and 6. The admittance is the parallel connection between node 6 and node 7. The admittance is the parallel connection between node 7 and the grounding point.

[0039] ;(twenty three) In the formula, The admittance connected in series between node 8 and the left grounding point (from Figure 3 It can be seen that the admittance connected in series between node 5 and the right grounding point, and the admittance connected in series between node 8 and the left grounding point, both adopt... express), The admittance is the parallel connection between nodes 8 and 9. The admittance is the series connection between node 9 and node 19. The admittance is the parallel connection between node 9 and node 10. The admittance is the parallel connection between node 10 and the grounding point.

[0040] Construct the current source vectors from left to right: ;(twenty four) In the formula, This is a current source connected in parallel between node 1 and the left-side grounding point. This is a current source connected in parallel between node 1 and node 2. This is a voltage source connected in series between node 1 and the grounding point on the right. This is a current source connected in parallel between node 2 and node 3. and This is a current source connected in parallel between node 3 and node 4. and This is a current source connected in parallel between node 4 and the grounding point.

[0041] (25) In the formula, This is a voltage source connected in series between node 5 and the left-side grounding point. This is a voltage source connected in series between node 5 and the grounding point on the right. This is a current source connected in parallel between node 5 and node 6. This is a current source connected in parallel between node 6 and node 7. This is a current source connected in parallel between node 7 and the grounding point.

[0042] (26) In the formula, This is a voltage source connected in series between node 8 and the left-side grounding point. This is a current source connected in parallel between node 8 and node 9. and These are the pulse voltage source and resistor connected in series between node 9 and node 19, respectively. This is a current source connected in parallel between node 9 and node 10. This is a current source connected in parallel between node 10 and the grounding point.

[0043] Similarly, for the zero-mode network and the fault component, the linear-mode network is disconnected and the fault component is grounded. The node admittance matrix is ​​also constructed sequentially from left to right. , , and current vector , , .

[0044] Consider the mutual admittance contribution matrix between node 9 and node 19 in the original network. The matrix dimension is 20×20, and the mutual admittance contribution matrix is... The elements in are defined as: (27) The complete nodal admittance matrix of the DC system can now be written. : (28) Current source vector of DC system for: (29) Define the node voltage vector as ( For nodes n voltage, n =1,2,...,20), the superscript T denotes matrix transpose, and the voltages at each node are calculated using the nodal voltage method: (30) Calculation protection installation location , , , Positive current: (31) In the formula, , , , for time , , , The positive current, , , , Before the fault , , , The positive current, , , They represent in The voltages at time nodes 1, 5, and 8 This represents the impedance connected in series between node 15 and the right-side grounding point, and also represents the impedance connected in series between node 18 and the left-side grounding point. This represents the impedance connected in series between node 15 and the left grounding point, and also represents the impedance connected in series between node 11 and the right grounding point.

[0045] Under the condition that the system excitation source is a unit pulse, the unit pulse response of the positive current at the protection installation point can be obtained using equations (30) and (31). Different faults are used to represent the most extreme external fault conditions under different operating conditions. A unit pulse excitation signal (i.e., input signal) is input to the DC system, and the established lossy Begeron model is used to calculate the unit pulse response signal of the current at the protection installation point under the most severe external fault.

[0046] Preferably, the sampling frequency is set. Select data window length Based on the time of the failure Set the time window for data extraction as the time reference point. The total number of sampling points within the time window is Extract the discrete sequence of the output signal within the time window. Define the discrete sequence of the output signal. The maximum value in is The minimum value is The final unit impulse response sequence is constructed as the output sequence. For the output sequence The first in element The calculation formula is as follows: (32) In the formula, Represents the discrete sequence of the output signal within the time window. The first in Each element.

[0047] Similarly, the input signal is discretized using the aforementioned time window to obtain the input sequence. In this embodiment, the input sequence is a unit impulse excitation sequence, and the output sequence is a unit impulse response sequence. The input and output sequences serve as the data for system identification in step 2.

[0048] In this embodiment, the fault location is set as Figure 2 In The sampling frequency is set to [value]. That is, the sampling period Based on the time of the failure Set the time window for data extraction as the time reference point. The total number of sampling points within the time window is According to the Nyquist sampling theorem, the upper limit of its effective observation bandwidth is... Select The data window length, calculated to have an inherent frequency resolution of approximately [missing information], is given. For ease of engineering calculations, this resolution is standardized to... Thus, the final frequency range was determined to be Six sets of experiments were set up, and the specific experimental conditions are shown in Table 1 below: Table 1: 6 sets of experimental conditions The selected specific frequencies are respectively , Step 14 yields the equivalent time-domain calculation model of the Begeron transmission line at the corresponding frequency. Based on the lossy Begeron model, theoretical calculations of the fault current traveling wave are performed on six sets of experiments. Taking a single-pole metallic grounding operation at a specific frequency of 350 Hz as an example, the normalized output signal curve is calculated as follows: Figures 4 to 9 As shown.

[0049] Preferably, step 2 specifically includes: Step 21: Based on the Z-transform of the input sequence and its corresponding output sequence, the discrete-time transfer function of the black-box network is represented in rational fraction form to construct an IIR filter model.

[0050] Specifically, the unit pulse excitation sequence (i.e., the discrete sequence of the input signal) is used as the input sequence. The calculated unit impulse response sequence (i.e., the discrete sequence of the output signal) is used as the output sequence. For the first Group( Given the input and output sequences, the corresponding IIR filter model structure is defined as follows: (33) In the formula: For the discrete-time transfer function of a black-box network, and These are the output sequences (unit impulse response sequences). and the input sequence (unit pulse excitation sequence). of Transformation; and Express the numerator and denominator polynomials of the IIR filter model respectively; and They represent delays respectively. and One sampling period; and Let be the orders of the denominator polynomial and the numerator polynomial, respectively. Based on the general case, assume that... ; ( )and ( ) represents the parameters of the denominator and numerator of the IIR filter model to be identified; Define the parameter vector of the IIR filter model to be solved. for: (34) The corresponding time-domain difference equation is expressed as: (35) In the formula, For the first The output value at each time point For the first The output value at each time point For the first The input value at each moment.

[0051] In this embodiment, a 7th-order IIR filter model is constructed, which in... The specific form under the domain is: (36) Step 22: Construct a regression matrix and an output vector based on the input sequence and its corresponding output sequence, and construct a parametric equation for the IIR filter model parameter vector based on the regression matrix and the output vector. Solve the parametric equation using the linear least squares method to obtain the initial values ​​of the IIR filter model parameter vector.

[0052] Specifically, to achieve faster convergence and more reliable convergence results, the initial values ​​of the IIR filter model parameter vector are first obtained using the linear least squares method before starting the Steiglitz-McBride iteration. .

[0053] First, construct a regression matrix based on the input and output sequences. and output vector : (37) express from The output sequence that is truncated afterward has the following specific form: (38) Regression Matrix The specific form is: (39) Then, the problem of finding the initial values ​​of the parameters is transformed into solving linear matrix equations: (40) In the formula, This is the residual vector.

[0054] Finally, the initial parameter estimates are obtained by solving the normal equations: (41) Step 23: Based on the initial values ​​of the IIR filter model parameter vector, the Steiglitz-McBride method is used to iteratively calculate the IIR filter model parameter vector to obtain the IIR filter model parameter vector.

[0055] Specifically, the Steiglitz-McBride method (also known as the Steiglitz-McBride iterative algorithm, or SM algorithm for short) is used to correct the bias caused by measurement noise or model mismatch through iterative methods.

[0056] In the In the next iteration ( ), perform the following operations: Step 231: Construct a pre-filter and perform signal filtering. In the... In the iteration, using the first The prefilter is constructed from the denominator polynomial coefficients obtained in the next iteration. : (42) For the original input sequence respectively and output sequence Filtering is performed to obtain the pre-filtered signal. and : (43) In the time domain, this is equivalent to solving the following difference equation: (44) Step 232: Update parameter estimates. Use the pre-filtered signal. and Construct a new regression matrix and output vector Its structure is the same as in step 22, except that the original data is replaced with filtered data. Parameters are updated by solving a system of linear equations. : (45) Step 233: Convergence Test. After obtaining the... Parameter estimates for the next iteration Then, first calculate the output error cost function corresponding to this set of parameters. The output error cost function is defined as the measurement output. With model output Mean square error between: (46) in, Indicates using the current iteration parameters For input signal The model prediction output is obtained after filtering. Based on the convergence characteristics of the SM algorithm, when the error is small, the output error sequence satisfies the monotonically non-increasing property, i.e. Therefore, the iteration termination condition is set as the difference in cost function between two adjacent iterations being less than a preset threshold. : (47) If this condition is met, or the maximum number of iterations is reached, then stop iterating and output the final parameters. Otherwise, let Return to step 231 and continue execution.

[0057] Step 24: Traverse each set of input sequences and their corresponding output sequences, and repeat steps 21 to 23 to obtain multiple sets of IIR filter model parameter vectors.

[0058] Specifically, for the product generated in step 1 The input sequence and its corresponding output sequence are processed repeatedly, and steps 21 to 23 are executed to finally identify the input sequence and its corresponding output sequence. The IIR filter model parameter vectors are grouped together to form the IIR filter model parameter vector set. Each set of IIR filter model parameter vectors It uniquely characterizes the frequency domain response features of the DC system black box network under the Bergeron model under different most extreme external fault conditions, providing a data foundation for the subsequent construction of internal and external protection zones.

[0059] In this embodiment, the final identification is obtained The parameter vectors of the group IIR filter model are shown in Tables 2 and 3 below: Table 2: Parent Parameters of Each Component Table 3: Molecular parameters for each group Preferably, step 3 specifically includes: Step 31: Based on the difference between the actual observed output and the model predicted output, calculate the residual vector of the parameter vector of each group of IIR filter models, and calculate the unbiased estimator of the noise variance of the parameter vector of each group of IIR filter models based on the residual vector of each group of IIR filter models, and calculate the covariance matrix of the estimated parameter vector of each group of IIR filter models based on the unbiased estimator of the noise variance of the parameter vector of each group of IIR filter models.

[0060] Specifically, for the first identified using the Steiglitz-McBride method in step 2... Group( IIR filter model parameter vector First, it is necessary to evaluate its estimation error. Assuming the DC system output contains zero-mean Gaussian white noise, the noise variance and parameter covariance matrix are estimated using the identification residuals. Define the... The residual vector of the IIR filter model parameter vector The difference between the actual observed output and the model's predicted output: (48) Wherein: number of sampling points , This represents the total number of sampling points within the time window. Let be the order of the IIR filter model; For the first Group output sequence observation vector; For the first The regression matrix corresponding to the input sequence and its corresponding output sequence contains the lag terms of the input and output signals. Based on the residual vector, the ... Unbiased estimator of the noise variance of the parameter vector of the IIR filter model : (49) in, For the first Group IIR filter model parameter vector The dimension of . Then, calculate the . Covariance matrix of parameter vector estimation for group IIR filter model : (50) The diagonal elements of the covariance matrix represent the variance of each parameter in the parameter vector of each IIR filter model, while the off-diagonal elements represent the correlation between the parameters.

[0061] Step 32: Extract the diagonal elements from the covariance matrix of the parameter vector estimation of each group of IIR filter models to calculate the standard error of each parameter in the parameter vector of each group of IIR filter models.

[0062] Specifically, from the first Covariance matrix of parameter vector estimation for group IIR filter model Extract the diagonal elements and calculate the first... Group IIR filter model parameter vector The Middle Parameters The standard error of . Definition variance vector for: (51) Then the first The first group of IIR filter model parameter vectors Standard error of each parameter The calculation is as follows: (52) Step 33: Based on each parameter and its standard error in the parameter vector of each group of IIR filter models, calculate the 95% confidence interval of each parameter in the parameter vector of each group of IIR filter models.

[0063] Specifically, based on the assumption that the parameter estimates follow an asymptotically normal distribution or a t-distribution under small samples, a significance level is set. (Corresponding to 95% confidence level). Search degrees of freedom are... Two-sided critical value of Student's t-distribution For the first The group identification results are used to calculate the first group. Group IIR filter model parameter vector Each parameter 95% confidence interval : Upper bound calculation: (53) Lower bound calculation: (54) Extending the above calculations to the parameter vectors of the entire IIR filter model, we obtain the first... The upper and lower bounds of the confidence interval for the parameter vectors of the IIR filter model: (55) (56) In this embodiment, the 95% confidence intervals of the constructed IIR filter model parameter vectors are as follows: Figure 10 As shown.

[0064] Step 34: Based on the 95% confidence intervals of each parameter in the parameter vectors of all groups of IIR filter models, construct the global upper and lower envelopes to obtain the guard interval.

[0065] Specifically, in order to cover the uncertainties caused by factors such as measurement noise and slight drift of the system operating point, the above calculations are used. Construct parameter protection intervals by grouping confidence intervals. Define the global upper envelope vector between parameter protection intervals. and global envelope vector For the first group of IIR filter model parameter vectors... The envelope of the parameters is calculated as follows: (57) (58) Ultimately, the obtained protection interval Represented as A hyperrectangle in 3D space is mathematically represented as: (59) This protected area This will serve as the benchmark threshold range for determining whether real-time monitoring parameters have experienced abnormal deviations (i.e., whether the fault is within or outside the zone) in subsequent steps.

[0066] In this embodiment, the constructed protection interval is as follows: Figure 11 As shown.

[0067] Preferably, step 4 specifically includes: Step 41: Real-time sampling of the positive voltage, negative voltage, positive current, and negative current at the DC system protection installation point is performed under the sliding time window to obtain the state vectors of the four phase domains at each sampling time. These vectors are then combined to form the state matrix of the phase domains at each sampling time. The phase mode transformation matrix is ​​then used to perform phase mode transformation on the state matrix of the phase domains at each sampling time to decouple the fault characteristics. This results in the state matrix of the module domains at each sampling time, which is composed of the state vectors of the four module domains. The state vectors of the four module domains include the zero-mode voltage vector, the line-mode voltage vector, the zero-mode current vector, and the line-mode current vector.

[0068] Specifically, set the sampling rate. 20kHz, long time window 500μs, based on the current time Using time as a reference, a sliding time window The data length is The positive and negative voltages, as well as the positive and negative currents, are simultaneously collected at the protection installation point using a high-frequency sampling device. Complete data collection at all times Within the time window, the state vectors of the four phase domains sampled can be represented as a state matrix of one phase domain: (60) In the formula, , , , Time windows The sampled positive voltage sequence, negative voltage sequence, positive current sequence, and negative current sequence are the state vectors of the four phase domains. , , , They are respectively Positive voltage, negative voltage, positive current, and negative current at any given time.

[0069] To decouple fault characteristics and eliminate the effects of inter-electrode coupling, a phase mode transformation matrix is ​​constructed. Phase mode transformation matrix for: (61) The sampling time is obtained through linear transformation. The state matrix of the modulus domain , including zero-mode voltage vector Line-mode voltage vector Zero-mode current vector and line mode current vector : (62) Step 42: Perform discrete operations on the line-mode voltage, zero-mode voltage, line-mode current, and zero-mode current in the state matrix of the mode domain at each sampling time based on the absolute gradient to obtain the gradient-enhanced energy matrix; normalize the gradient-enhanced energy matrix to obtain the normalized gradient-enhanced energy matrix.

[0070] Specifically, a composite feature matrix integrating "frequency abrupt changes" and "wavefront steepness" is constructed to obtain strong sensitivity. Definition They are respectively The line-mode voltage, zero-mode voltage, line-mode current, and zero-mode current at time t, for Discrete TEO (TEO: Gradient-Enhanced) Energy Vector of the Modular Domain State Based on Absolute Gradient at Each Moment Perform the calculation: (63) In the formula, , , , They represent The gradient-enhanced energy vectors of zero-mode voltage, line-mode voltage, zero-mode current, and line-mode current at time t; TEO energy matrix under sliding time window Represented as: (64) In the formula, , , , Let represent the gradient-enhanced energy vectors of the zero-mode voltage, line-mode voltage, zero-mode current, and line-mode current at the sampling time within the sliding time window, respectively. Define the vectors. maximum value ,vector minimum value Similarly, the definition The maximum values ​​are respectively , The minimum values ​​are respectively For the TEO energy matrix Normalization is performed: (65) In the formula, This is the normalized TEO energy matrix.

[0071] Step 43: Define the normalized weighted coefficient vectors for line, zero-mode voltage, and current, and calculate the start-up decision value vector for each sampling time based on the normalized gradient-enhanced energy matrix and the normalized weighted coefficient vectors.

[0072] Specifically, the sum of the elements in the normalized weighted coefficient vector is 1. The normalized weighted coefficient vector is used to aggregate features. for: (66) in, Representing line-mode voltage Zero-mode voltage Linear current and zero-mode current Normalized weights. Generally set as follows: and Since grounding faults are mainly manifested in the zero-mode component, and the zero-mode component is mainly affected by the system structure, grounding method and fault type, and is less affected by the control system, a higher weighting coefficient should be set to improve sensitivity.

[0073] Using normalized weighted coefficient vectors Calculate the initiation decision value vector : (67) In the formula, for The vector of start-up decision values ​​within the time window.

[0074] Step 44: Calculate the mean of the start-up decision value vector at each sampling time to obtain the mean of the start-up decision value at each sampling time; calculate the start-up dynamic decision threshold at each sampling time based on the reliability coefficient and the mean of the start-up decision value at the previous sampling time; determine whether the mean of the start-up decision value at each sampling time is greater than or equal to the start-up dynamic decision threshold at the corresponding sampling time. If not, determine that no fault has occurred in the DC system; if so, determine that a fault has occurred in the DC system and record the time of the fault occurrence.

[0075] Specifically, for the initiation decision value vector Calculate the mean: (68) In the formula, China The average value of the start-up decision at any given time. (Setting) Dynamic decision threshold for activation at any time For reliability coefficient and The product of the average of the initiation decision values ​​at time points can be expressed as: (69) In the formula, Let be the reliability coefficient. The criterion logic is: if the following equation is satisfied: (70) If a fault occurs in the DC system, the first time this condition is met is recorded as the fault occurrence time. Proceed to step 5. If the criterion is not met, repeat step 4.

[0076] In this embodiment, the midpoint of the line is set. A single-pole grounding fault occurs at the location, and a reliability factor is set. Normalized weights: , Fault diagnosis results are as follows: Figure 12 As shown.

[0077] Depend on Figure 12 It can be seen that, in If the criteria are met at any time, a fault is determined to have occurred, the time of the fault is recorded, and the process proceeds to step 5.

[0078] Preferably, step 5 specifically includes: Step 51: Input the falling edge signal as the experimental input signal into the DC system to conduct a fault experiment, and collect the fault current traveling wave at the protection installation point in the DC system as the experimental observation signal with the fault occurrence time as the time reference point. Normalize the experimental input signal and the experimental observation signal respectively to obtain the experimental input vector and the experimental observation vector.

[0079] Specifically, the DC system adopts the same sampling method as in step 1, and the sampling frequency is set to [value missing]. Based on the time of the failure Set the time window for data extraction as the time reference point. The total number of sampling points within the time window is The setup failure occurred Figure 2 In In the experiment, the experimental input signal can be equivalent to a falling edge signal, which is then normalized to obtain the experimental input vector for subsequent system identification. Meanwhile, at the protective installation location Acquire current signals (i.e., fault current traveling waves) as experimental observation vectors for DC systems. Experimental input vector and experimental observation vector It can be represented as: (71) (72) In the formula, exist The sudden amplitude changes at any moment characterize the broadband excitation source of the DC system.

[0080] In this embodiment, the fault current traveling wave detection location is the protection installation point. For the location of a single-pole grounding fault (For external faults), a metallic single-pole grounding fault experiment was set up to obtain the corresponding fault current traveling wave.

[0081] Step 52: Define a regression vector based on the experimental input vector and the time-domain difference equation of the output signal of the lossy Begeron model; regard the experimental observation vector as the desired signal and the output sequence of the lossy Begeron model as the estimated signal; calculate the prior estimation error based on the regression vector and the experimental observation vector; and iteratively update the IIR filter model parameter vector based on the prior estimation error to obtain the converged IIR filter model parameter vector.

[0082] Specifically, the parameter vector of the IIR filter model to be solved is defined according to equation (34) in step 2. The time-domain difference equation shown in equation (35) is used to transform the frequency-varying system into a fixed-parameter system using the NLMS algorithm. During this process, the experimental observation vectors are... Treating it as the desired signal will impair the output sequence of the Bergeron model. Considered an estimated signal. Definition Regression vector at time step (Including historical inputs and historical outputs): (73) The Normalized Least Mean Square (NLMS) algorithm is used to... IIR filter model parameter vector at time 1 Perform iterative updates to minimize the output sequence. middle Output at time With experimental observation vector middle Experimental observations at time The error between them. Calculate. Prior estimation error at time : (74) According to the NLMS criterion, the update iterative formula for the IIR filter model parameter vector is: (75) in: The step size factor has a range of values. This is used to control the convergence speed and steady-state error; for Time regression vector The square of the Euclidean norm is used to normalize the input energy; To prevent tiny positive numbers with a denominator of zero (regularization parameter).

[0083] Step 53: Use the converged IIR filter model parameter vector to predict the experimental input vector over the entire time period to obtain the predicted output vector, and calculate the fitting degree between the predicted output vector and the experimental observation vector.

[0084] Specifically, after completing iterative training of the IIR filter model parameter vector within the time window, a converged IIR filter model parameter vector is obtained. The parameter vector of the converged IIR filter model is used to evaluate the experimental input vector. Perform full-time prediction to obtain the prediction output vector. Calculate the predicted output vector. With experimental observation vector The fit between The normalized goodness-of-fit formula is defined as follows: (76) In the formula, To predict the output vector middle Predicted output at time step Experimental observation vector The mean.

[0085] Step 54: Determine whether the fitting degree is greater than or equal to the preset fitting degree threshold, or determine whether the number of iterations of the IIR filter model parameter vector has reached the preset maximum number of iterations. If not, adjust the step size factor of the IIR filter model and repeat steps 52 to 53. If yes, use the currently obtained predicted output vector as the final predicted output vector and execute steps 55 to 56.

[0086] Specifically, when If the maximum number of iterations is reached, the system identification is considered complete, and the converged IIR filter model parameter vector is determined. Accurately reflects the frequency-varying characteristics of the black-box network; otherwise, adjust the step size factor of the IIR filter. Repeat steps 52 to 53.

[0087] In this embodiment, the final model fitting result is as follows: Figure 13 As shown.

[0088] Step 55: Perform first-order difference calculation on the experimental input at two adjacent time points in the experimental input vector to deduce the standard unit impulse input sequence; perform first-order difference calculation on the predicted output at two adjacent time points in the final predicted output vector to deduce the original unit impulse response sequence.

[0089] Specifically, based on the theory of discrete-time linear time-invariant systems, the final predicted output vector is used. and experimental input vector The standard unit pulse input sequence is derived by using first-order difference. and the original unit impulse response .

[0090] Using the first-order difference to reverse the process Standard unit pulse input of time : (77) Input the standard unit pulse at each moment Combined to form a standard unit pulse input sequence ; Using the first-order difference to reverse the process The raw unit impulse response at time t. : (78) The raw unit impulse response at each time point Combined to form the original unit impulse response sequence .

[0091] Step 56: Normalize the standard unit impulse input sequence and the original unit impulse response sequence to obtain the normalized standard unit impulse input sequence and the normalized original unit impulse response sequence.

[0092] Specifically, to eliminate the influence of the initial amplitude of the fault and to meet the requirements of the subsequent system identification algorithm for the input data range, the standard unit impulse input sequence and the original unit impulse response sequence obtained in step 55 are standardized. The original unit impulse response sequence... For example, define the original unit impulse response sequence. The maximum value in is The minimum value is Construct the normalized original unit impulse response sequence The calculation formula is as follows: (79) In the formula, for The normalized original unit impulse response at time t. The normalized original unit impulse response sequence. The range of values ​​for is strictly limited to between.

[0093] The normalized standard unit pulse input sequence and the normalized original unit impulse response sequence The dataset is combined and used for system identification in step 6.

[0094] Preferably, step 6 specifically includes: Step 61: Construct a regression vector for real-time monitoring data based on the normalized standard unit impulse input sequence and the normalized original unit impulse response sequence.

[0095] Specifically, based on a set of observational data obtained in step 5 after the fault occurred, namely the normalized standard unit pulse input sequence and the normalized original unit impulse response sequence The observed data are organized into a matrix form to construct a data structure suitable for IIR filter model identification. The data length of the time window is... Real-time monitoring data regression vector for: (80) in, and These are the denominator and numerator orders of the IIR filter model, respectively.

[0096] Step 62: Based on the Steiglitz-McBride method, iteratively calculate the regression vector of the real-time monitoring data to obtain the real-time IIR filter parameter vector that reflects the dynamic characteristics of the current DC system.

[0097] Specifically, using the Steiglitz-McBride iterative algorithm described in step 2, the deviation caused by output noise is eliminated through pre-filtering. After iterative convergence, the real-time IIR filter parameter vector reflecting the dynamic characteristics of the current DC system is obtained, denoted as... : (81) in, The total dimension of the feature parameters.

[0098] Step 63: Based on the parameter comparison method, determine whether each parameter in the real-time IIR filter parameter vector falls within the protection interval. If so, determine that an external fault has occurred and control the protection device to remain inactive; otherwise, determine that an internal fault has occurred and control the protection device to operate.

[0099] Specifically, the parameter protection range pre-built in step 3 is invoked. The protected area is in In 3D space, it forms a hyper-rectangle, which is composed of the global envelope vector. and global envelope vector Unique definition: (82) (83) The parameter vector of the real-time IIR filter is determined by parameter-by-parameter comparison. Each parameter in Does it satisfy the following inequality constraints: (84) If the real-time IIR filter parameter vector All parameters satisfy this inequality constraint, that is This indicates that the current system parameter offset is still within the preset range for out-of-zone fault fluctuations. In this case, the judgment is an out-of-zone fault, and the protection device does not operate; conversely, if the real-time IIR filter parameter vector... There exist one or more elements in the inequality that do not satisfy the inequality constraint, i.e. ,lead to This indicates that the transmission characteristics of the system have changed drastically, exceeding the characteristic boundary of an external fault. At this point, the output result is determined to be an internal fault, and the corresponding protection trip logic is initiated.

[0100] In this embodiment, the location where the fault occurs is... (Internal fault) and (External faults) Consider conducting single-pole grounding fault and inter-pole short-circuit fault experiments respectively. For these four operating conditions, use steps 5 and 6 to obtain the results. and The system IIR filter parameters at the time of a fault can be used to determine whether the fault occurs within or outside the zone. Figure 14 and Figure 15 As shown.

[0101] Depend on Figure 14 It can be seen that during faults within the protection zone, both single-pole metallic grounding and inter-pole short circuits have system identification parameters falling outside the protection zone, thus indicating faults within the protection zone. Figure 15It can be seen that when there is a fault outside the protection zone, the system identification parameters for both single-pole metallic grounding and inter-pole short circuit fall within the protection zone, and the fault is determined to be outside the protection zone.

[0102] In summary, in the DC system relay protection method based on equivalent filter parameter identification of this invention, a lossy Beeron model of the DC system is first constructed, and the operating conditions of the most severe external fault zone are set. A series of experiments were conducted to extract fault characteristics outside the protection zone. Each experiment used a unit pulse signal as the input signal and employed a lossy Beigeron model to calculate the current at the protection installation point during a DC system fault, setting the sampling frequency accordingly. Select data window length Based on the time of the failure Based on this, the time window for data extraction is set as follows: The current pulse response is discretized and used as the output signal. Construction The IIR filter model has a total of One parameter, to Group Experiments The input and output signal sequences were analyzed, and the Steiglitz-McBride method was used for system identification to obtain... Group IIR filter parameters. Calculate the upper and lower interval matrices of the 95% confidence interval for each IIR coefficient in each group. Each parameter is used as a boundary, taking the maximum / minimum value of each IIR coefficient in the upper / lower interval matrix as the bound. This boundary is then used to... The parameters of each IIR filter are used to construct upper and lower envelopes, and the envelope region serves as the protection interval for distinguishing faults inside and outside the detection zone. Using a sliding time window, four state variables—positive and negative voltage and current—at the protection installation location are collected and converted into moduli using a sampling device. The enhanced energy (EEG) of these four state variables is calculated using TEO energy and a first-order differential gradient, and then normalized. A decision value is calculated for each of the four state variables using a normalized weight matrix. The average decision value within the sliding time window is calculated. The average decision value from the previous time step is multiplied by a reliability coefficient to serve as the dynamic decision threshold for the next time step. When the decision value at the next time step exceeds this threshold, a fault is determined to have occurred, and the fault time is recorded. A fault current traveling wave was obtained by setting up a fault experiment in a DC system. In the experiment, the input signal sequence can be equivalent to a falling edge signal. Both the fault current traveling wave and the falling edge signal were normalized. Since the experimental data is a dispersive signal, the NLMS method was used to filter the experimental signal. The unit impulse excitation signal and unit impulse response were obtained by reverse engineering the falling edge signal and the filtered experimental signal, respectively. These were used as the input and output signals, and the Steiglitz-McBride method was used for system identification to obtain the real-time IIR filter parameters. The system was considered safe when all real-time parameters fell within the protection range. When any parameter does not fall within the protection range, it is determined to be an external fault; When this occurs, it is determined to be a fault within the affected area.

[0103] Based on the above-mentioned DC system relay protection method based on equivalent filter parameter identification, the present invention also provides a DC system relay protection system based on equivalent filter parameter identification.

[0104] A DC system relay protection system based on equivalent filter parameter identification includes a processor, a memory, and a computer program stored in the memory. When the computer program is executed by the processor, it implements the DC system relay protection method based on equivalent filter parameter identification as described above.

[0105] The DC system relay protection method and device based on equivalent filter parameter identification of this invention improves the reliability, speed, sensitivity, and selectivity of the protection. Specifically: In terms of reliability: First, this invention does not rely on a two-end communication system for data synchronization. For high-voltage DC transmission lines that have experienced faults, the implementation of this invention can complete the start-up decision based solely on local positive and negative voltage and current at the protection installation point. This is protection based on local quantities, and its implementation method is relatively reliable. Second, this invention only uses experimental observation signals during operation for online adaptive learning. Furthermore, before fault identification, it uses a normalized least mean square algorithm to iteratively update the IIR filter model parameters online, eliminating the risk of parameter fluctuations caused by changes in power grid operation mode or topology reconfiguration, and significantly improving the adaptive capability and operational reliability of the protection system under complex operating conditions. In terms of speed: In the fault feature extraction and start-up decision stages, this invention introduces discrete energy calculation and first-order differential gradient calculation to construct a gradient-enhanced energy matrix, which can reliably trigger start-up within the first few milliseconds of a fault occurrence using only a very short time data window; In addition, this invention cleverly transforms the complex frequency-varying nonlinear DC system into a discrete-time linear time-invariant system with fixed parameters, and uses a normalized least mean square algorithm with extremely low computational complexity for iteration; This "dimensionality reduction and linearization" process minimizes the computing power consumption of the relay protection's underlying microprocessor, improving the protection's operating speed from an algorithmic perspective; In terms of sensitivity and selectivity: This invention calculates the out-of-zone fault conditions of DC systems based on the lossy Begeron model and constructs the IIR filter model of the black box to be identified. The upper and lower envelopes of the confidence interval generated by parameter identification are used as protection intervals to measure whether the real-time parameter vector of the system exceeds the limit (whether an in-zone fault occurs). Compared with the traditional analytical method or pure data-driven method that strongly depends on the accuracy of the model, it has a high resistance to transition resistance to cope with in-zone faults and reliably does not operate in the out-of-zone fault scenario.

[0106] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A DC system relay protection method based on equivalent filter parameter identification, characterized in that, include: Step 1: Construct a lossy Beerlong model of the DC system and set multiple different input signals according to the most severe external fault conditions; The lossy Begeron model is used to calculate the external fault characteristics of the DC system for each group of input signals to obtain the corresponding output signal. Step 2: Treat the DC system as a black box network to be identified, construct the IIR filter model of the black box network, and use the Steiglitz-McBride method to identify the parameters of the IIR filter model using multiple sets of input signals and their corresponding output signals to obtain multiple sets of IIR filter model parameter vectors. Step 3: Calculate the confidence interval of each group of IIR filter model parameter vectors, and based on the confidence intervals of multiple groups of IIR filter model parameter vectors, perform upper and lower envelope operations on the multiple groups of IIR filter model parameter vectors to obtain the guard interval; Step 4: Real-time sampling of the positive and negative voltages and currents at the DC system protection installation point and mode decoupling to construct a state matrix; and discrete energy calculation and first-order differential gradient calculation are performed on the state matrix to construct a gradient-enhanced energy matrix. Define a weighted coefficient vector for line, zero-mode voltage, and current, and calculate a start-up decision value vector based on the gradient-enhanced energy matrix and the weighted coefficient vector. Determine whether a fault has occurred in the DC system based on the start-up decision value vector, and record the time of fault occurrence when a fault is determined to have occurred in the DC system. Step 5: Set up the fault experiment and obtain the experimental observation signal by inputting the experimental input vector into the DC system; use the normalized least mean square algorithm to transform the experimental observation signal into a fixed parameter system based on the frequency-varying system, and iteratively update the IIR filter model parameter vector; then predict the output of the DC system based on the iteratively updated IIR filter model parameter vector to obtain the final predicted output vector; based on the discrete-time linear time-invariant system theory, use the experimental input vector and the final predicted output vector to back-calculate using the first-order difference to obtain the standard unit impulse input sequence and the original unit impulse response sequence. Step 6: Based on the Steiglitz-McBride method, calculate the parameters of the IIR filter model according to the standard unit impulse input sequence and the original unit impulse response sequence to obtain the real-time IIR filter model parameter vector; based on the real-time IIR filter model parameter vector and the protection interval, determine the interval in which the DC system has a fault, and then execute the corresponding protection.

2. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, In step 1, the lossy Begeron model of the DC system is constructed, specifically including: Step 11: Equivalent the inductor in the DC system to a parallel connection of a resistor and a current source to construct the inductor Bergeron model; Step 12: Equivalently represent the capacitor in the DC system as a parallel connection of a resistor and a current source to construct a capacitor Begeron model; Step 13: In the modal domain, the fault point of the DC system is equivalent to a resistor and voltage source connected in series and grounded to construct the equivalent circuit of the line zero mode fault component. Step 14: Construct a time-domain equivalent calculation model of the Begeron transmission line considering frequency-dependent losses for the DC system; Step 15: Based on the inductor Begeron model, capacitor Begeron model, equivalent circuit of line-zero mode fault component, and Begeron transmission line time-domain equivalent calculation model constructed in steps 11 to 14, set all equivalent power supplies to zero, and for solving the Begeron model with pulse signal as input excitation, replace the step signal of the fault component with a unit pulse signal to construct a lossy Begeron model of the DC system in zero initial state.

3. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, The lossy Begeron model includes linear model networks and zero-model networks; In step 1, the lossy Begeron model is used to calculate the external fault characteristics of the DC system for each group of input signals to obtain the corresponding output signal, specifically including step 16: Disconnect the zero-mode network and ground the fault component to construct the node admittance matrix and current source vector of the line-mode network; The linear mode network is disconnected and the fault component is grounded to construct the node admittance matrix and current source vector of the zero-mode network; Construct the mutual admittance contribution matrix of the fault point in the online mode network and the zero mode network, and construct the node admittance matrix of the DC system based on the node admittance matrix of the online mode network, the node admittance matrix of the zero mode network and the mutual admittance contribution matrix. Construct the DC system current source vector based on the current source vector of the linear mode network and the current source vector of the zero mode network; Based on the nodal voltage method, a nodal voltage calculation model is constructed to calculate the voltage of each node in the DC system, according to the nodal admittance matrix of the DC system and the current source vector of the DC system. A calculation model for the positive current at the protection installation point in a DC system is constructed to calculate the positive current at the protection installation point at the time of a fault. The input signal is substituted into the node voltage calculation model and the positive current calculation model at the protection installation point to calculate the output signal. Specifically, the input signal is a unit pulse excitation signal, and the output signal is specifically the unit pulse response signal at the protection installation location under the most severe external fault zone.

4. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, Before parameter identification, step 2 also includes: Set sampling frequency Select data window length Based on the time of the failure Using the time reference point, a time window for data extraction is constructed as follows: ; Extract the discrete sequence of the output signal within the time window. And define the discrete sequence of the output signal within the time window. The maximum value in is The minimum value is ; Based on the discrete sequence of the output signal within the time window Given its maximum and minimum values, construct the output sequence. The output sequence The first in element Represented as: ； In the formula, Represents the discrete sequence of the output signal within the time window. The first in One element; Extract the discrete sequence of the input signal within the time window to obtain the input sequence; The input sequence and its corresponding output sequence are used as the data for parameter identification of the IIR filter model in step 2.

5. The DC system relay protection method based on equivalent filter parameter identification according to claim 4, characterized in that, Step 2 specifically includes: Step 21: Based on the Z-transform of the input sequence and its corresponding output sequence, the discrete-time transfer function of the black-box network is represented in rational fraction form to construct an IIR filter model; Step 22: Construct a regression matrix and an output vector based on the input sequence and its corresponding output sequence, and construct a parametric equation for the IIR filter model parameter vector based on the regression matrix and the output vector. Solve the parametric equation using the linear least squares method to obtain the initial values ​​of the IIR filter model parameter vector. Step 23: Based on the initial values ​​of the IIR filter model parameter vector, the Steiglitz-McBride method is used to iteratively calculate the IIR filter model parameter vector to obtain the IIR filter model parameter vector. Step 24: Traverse each set of input sequences and their corresponding output sequences, and repeat steps 21 to 23 to obtain multiple sets of IIR filter model parameter vectors.

6. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, Step 3 specifically includes: Step 31: Based on the difference between the actual observed output and the model predicted output, calculate the residual vector of the parameter vector of each group of IIR filter models, and calculate the unbiased estimator of the noise variance of the parameter vector of each group of IIR filter models based on the residual vector of each group of IIR filter models, and calculate the covariance matrix of the estimated parameter vector of each group of IIR filter models based on the unbiased estimator of the noise variance of the parameter vector of each group of IIR filter models. Step 32: Extract the diagonal elements from the covariance matrix of the parameter vector estimation of each group of IIR filter models to calculate the standard error of each parameter in the parameter vector of each group of IIR filter models. Step 33: Based on each parameter and its standard error in the parameter vector of each group of IIR filter models, calculate the 95% confidence interval of each parameter in the parameter vector of each group of IIR filter models. Step 34: Based on the 95% confidence intervals of each parameter in the parameter vectors of all groups of IIR filter models, construct the global upper and lower envelopes to obtain the guard interval.

7. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, Step 4 specifically includes: Step 41: Real-time sampling of the positive voltage, negative voltage, positive current, and negative current at the DC system protection installation point is performed under a sliding time window to obtain the state vectors of the four phase domains at each sampling time. These vectors are then combined to form the state matrix of the phase domains at each sampling time. The phase mode transformation matrix is ​​then used to perform phase mode transformation on the state matrix of the phase domains at each sampling time to decouple the fault characteristics. This results in the state matrix of the module domains at each sampling time, which is composed of the state vectors of the four module domains. The state vectors of the four module domains include the zero-mode voltage vector, the line-mode voltage vector, the zero-mode current vector, and the line-mode current vector. Step 42: Perform discrete operations on the mode domain state variables based on absolute gradient for the line-mode voltage, zero-mode voltage, line-mode current, and zero-mode current in the mode domain state matrix at each sampling time to obtain the gradient-enhanced energy matrix; normalize the gradient-enhanced energy matrix to obtain the normalized gradient-enhanced energy matrix. Step 43: Define the normalized weighted coefficient vectors for line, zero-mode voltage, and current, and calculate the start-up decision value vector for each sampling time based on the normalized gradient-enhanced energy matrix and the normalized weighted coefficient vectors. Step 44: Calculate the mean of the start-up decision value vector at each sampling time to obtain the mean of the start-up decision value at each sampling time; calculate the start-up dynamic decision threshold at each sampling time based on the reliability coefficient and the mean of the start-up decision value at the previous sampling time; determine whether the mean of the start-up decision value at each sampling time is greater than or equal to the start-up dynamic decision threshold at the corresponding sampling time. If not, determine that no fault has occurred in the DC system; if so, determine that a fault has occurred in the DC system and record the time of the fault occurrence.

8. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, Step 5 specifically includes: Step 51: Input the falling edge signal as the experimental input signal into the DC system to conduct a fault experiment, and collect the fault current traveling wave at the protection installation point in the DC system as the experimental observation signal with the fault occurrence time as the time reference point. Normalize the experimental input signal and the experimental observation signal respectively to obtain the experimental input vector and the experimental observation vector. Step 52: Define a regression vector based on the experimental input vector and the time-domain difference equation of the output signal of the lossy Begeron model; regard the experimental observation vector as the expected signal and the output sequence of the lossy Begeron model as the estimated signal; calculate the prior estimation error based on the regression vector and the experimental observation vector; and iteratively update the IIR filter model parameter vector based on the prior estimation error to obtain the converged IIR filter model parameter vector. Step 53: Use the converged IIR filter model parameter vector to predict the experimental input vector over the entire time period to obtain the predicted output vector, and calculate the fitting degree between the predicted output vector and the experimental observation vector. Step 54: Determine whether the fitting degree is greater than or equal to the preset fitting degree threshold, or determine whether the number of iterations of the IIR filter model parameter vector has reached the preset maximum number of iterations. If not, adjust the step size factor of the IIR filter model and repeat steps 52 to 53. If yes, use the currently obtained predicted output vector as the final predicted output vector and execute steps 55 to 56. Step 55: Perform first-order difference calculation on the experimental input at two adjacent time points in the experimental input vector to deduce the standard unit impulse input sequence; perform first-order difference calculation on the predicted output at two adjacent time points in the final predicted output vector to deduce the original unit impulse response sequence. Step 56: Normalize the standard unit impulse input sequence and the original unit impulse response sequence to obtain the normalized standard unit impulse input sequence and the normalized original unit impulse response sequence.

9. The DC system relay protection method based on equivalent filter parameter identification according to claim 1, characterized in that, Step 6 specifically includes: Step 61: Construct a regression vector for real-time monitoring data based on the normalized standard unit impulse input sequence and the normalized original unit impulse response sequence; Step 62: Based on the Steiglitz-McBride method, iteratively calculate the regression vector of the real-time monitoring data to obtain the real-time IIR filter parameter vector that reflects the dynamic characteristics of the current DC system. Step 63: Based on the parameter comparison method, determine whether each parameter in the real-time IIR filter parameter vector falls within the protection interval. If so, determine that an external fault has occurred and control the protection device to remain inactive; otherwise, determine that an internal fault has occurred and control the protection device to operate.

10. A DC system relay protection device based on equivalent filter parameter identification, characterized in that, The system includes a processor, a memory, and a computer program stored in the memory, wherein the computer program, when executed by the processor, implements the DC system relay protection method based on equivalent filter parameter identification as described in any one of claims 1 to 9.

Images