Power grid operation delay control method under emergency situation

Abstract

The invention discloses a power grid operation delay control method under an emergency situation, and relates to the technical field of power grid control and protection. Comprising the following steps: S1, carrying out time compensation on a PMU device, inserting an SCADA virtual timestamp, and carrying out synchronous time compensation; s2, performing data processing in different working condition modes according to the working condition modes and synchronization time compensation; s3, determining a total time error needing to be compensated according to the working condition mode, the equipment response time and the compensation coefficient, and determining corrected instruction issuing time; and S4, combining the corrected instruction issuing time with the predicted compensation time to determine the final instruction issuing time. According to the method, the local clock error of the PMU is compensated through GPS time, the SCADA virtual timestamps are inserted to realize time axis alignment of the multi-source data, and meanwhile, the preset number of virtual timestamps are inserted between SCADA sampling points, so that the time resolution difference between the SCADA sampling points and the PMU data is reduced, and the synchronization error is reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of power grid control and protection, and specifically to a method for controlling the operation delay of a power grid under emergency conditions. Background Art

[0002] With the large-scale access of renewable energy (such as photovoltaic and wind power) and the rapid development of DC loads, low-voltage DC distribution systems have gradually become an important part of modern power grids. However, such systems face many challenges during operation, especially voltage stability problems caused by emergencies (such as sudden changes in renewable energy output, load fluctuations, equipment failures, or communication interruptions). Due to the intermittency and uncertainty of distributed energy sources in the system, as well as the dynamic characteristics of DC loads, the system voltage is prone to severe fluctuations, over-limit, or even oscillatory divergence. In severe cases, it may lead to the failure of the system's safe operation, threatening power quality and power supply reliability.

[0003] Traditional control methods (such as droop control and master-slave control) perform well in dealing with normal operating conditions, but under sudden disturbances, their response speed and robustness are often insufficient. Especially in multi-terminal DC systems, the time-delay problem of the communication network further exacerbates the performance degradation of the control system. For example, when there is a delay in the communication link, the controller cannot obtain real-time status information in a timely manner, resulting in the disconnection between control instructions and the actual system state, thereby amplifying voltage fluctuations and even causing system instability. In addition, some existing control strategies (such as feedback control based on PID) are difficult to simultaneously consider multi-objective optimization (such as minimizing control costs and voltage fluctuations) and are difficult to achieve global optimality under complex operating conditions.

[0004] The Chinese invention patent with the publication number CN114172262A discloses an intelligent substation sampling data quality comprehensive evaluation method and system, including: obtaining the merged unit data, measuring and controlling device data, and PMU data in the substation to form multi-source redundant three-phase measurements under a unified time scale in the substation; using the multi-source redundant three-phase measurements in the substation, as well as the substation model and the real-time switch state of the circuit breaker, to perform three-phase linear state estimation of the substation, and using the weighted least squares algorithm to solve the state estimation problem; performing comprehensive comparison based on the state estimation results, using the Chi-square test to determine whether there are bad data in the measured quantities, and applying the maximum normalized residual LNR method to identify the most likely bad measurement quantities, and obtaining the comprehensive analysis result of data quality. This invention realizes the comprehensive evaluation of the sampling data quality analysis of an intelligent substation based on the multi-interval, multi-type, and available sampling data in the existing integrated monitoring system in the station, and improves the operation reliability of the substation.

[0005] However, during the process of comprehensive quality assessment of the existing sampling data in smart substations, since the PMU (Phasor Measurement Unit) captures the details of the power grid electromagnetic transient at a high frequency of thousands of Hz, while the SCADA (Supervisory Control and Data Acquisition) is limited by the low-speed sampling period of seconds, the synchronization accuracy between the heterogeneous data streams of the two is insufficient due to cross-scale timestamps. At the same time, when aligning the data, a fixed interpolation filtering rule for offline calibration is forced to be adopted, and the non-stationary signal characteristics at the moment of failure cannot be dynamically adapted, resulting in a delay of 80 - 150 ms in the real-time data cleaning link. Summary of the Invention

[0006] The purpose of the present invention is to provide a method for controlling the delay of power grid operation under sudden situations to solve the problems raised in the above background technology.

[0007] To achieve the above purpose, the present invention provides the following technical solution: A method for controlling the delay of power grid operation under sudden situations, including:

[0008] S1: Create a unified time baseline: According to the designed PMU channels and SCADA channels, perform time compensation for the PMU device, insert SCADA virtual timestamps, and perform synchronous time compensation;

[0009] S2: Dynamically select an interpolation filter bank: According to the operating condition mode and synchronous time compensation, perform data processing under different said operating condition modes, including:

[0010] S2.1: Determine the operating condition mode: According to the voltage, frequency, power analog quantity, and total harmonic distortion rate, determine the corresponding operating condition mode, and the operating condition mode includes a steady state mode and a transient mode;

[0011] S2.2: Steady state mode processing: When the operating condition mode is the steady state mode, determine the virtual interpolation data in the SCADA channel according to the preset number of virtual time points, and at the same time, according to the complex frequency variable and the cut-off angular frequency of the filter, obtain the transfer function of the low-pass filter;

[0012] S2.3: Transient mode processing: When the operating condition mode is the transient mode, determine the virtual interpolation data in the SCADA channel through the wavelet basis function and the preset number of virtual time points, and at the same time, according to the constructed state variable and observation model, determine the state equation and the observed value, and adaptively adjust the process noise according to the reference noise;

[0013] S3: Determine the control instruction timestamp error: According to the operating condition mode, device response time, and compensation coefficient, determine the total time error that needs to be compensated, and at the same time, according to the total time error that needs to be compensated, determine the corrected instruction issuance time;

[0014] S4: timestamp advance correction: the final instruction issuing time is determined based on the corrected instruction issuing time and the predicted compensation time, specifically:

[0015] t″ cmd =t′ cmd -Δt pre

[0016] Where: t″ cmd The time when the final revised instruction is issued, t′ cmd is the corrected instruction issuing time, Δt pre Compensate time for predictions.

[0017] Furthermore, synchronization time compensation is performed, including:

[0018] S1.1: Design separation channels: Design PMU channels and SCADA channels separately through optical fiber channels and dispatching data networks;

[0019] S1.2: Perform time compensation of the PMU device: determine the mean square error of all the PMU nodes through the local clock and GPS time at each PMU node, and obtain the total compensation amount of the PMU channel according to the correction item of optical fiber transmission, and adjust the local clock;

[0020] S1.3: inserting SCADA virtual timestamp: determining a preset number of virtual time points inserted into the SCADA channel according to the sampling period and target time resolution of the SCADA channel, and obtaining a corrected timestamp;

[0021] S1.4: Perform synchronization time compensation: Obtain the ambient temperature change rate of the clock module through the temperature sensor, and obtain the synchronized timestamp based on the total compensation amount of the PMU channel, the original uncompensated timestamp recorded by the local clock, and the ambient temperature change rate. Specifically:

[0022]

[0023] Where: t syn is the timestamp after synchronization, t raw is the original uncompensated timestamp recorded by the local clock, Δt PMU is the total compensation of the PMU channel, δ is the crystal temperature drift coefficient, Δt obs is the observation time window, is the rate of change of ambient temperature.

[0024] Furthermore, the local clock is adjusted, including:

[0025] S1.2.1: Perform clock error compensation: At each PMU node, obtain the local clock and, based on the local clock and GPS time, obtain the time difference, specifically:

[0026]

[0027] Where: Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, is the timestamp recorded by the local clock at the i-th PMU node, is the GPS timestamp received at the i-th PMU node, and i is the index of the PMU node;

[0028] S1.2.2: Obtain PMU node error: Based on the time difference at each PMU node, obtain the mean square error of all PMU nodes, specifically:

[0029]

[0030] Where: t s is the mean square error of all PMU nodes, N is the total number of PMU nodes, and Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, and i is the index of the PMU node;

[0031] S1.2.3: Determine delay correction: Based on the standard deviation of the fiber channel transmission delay, determine the correction term for fiber transmission, specifically:

[0032]

[0033] Where: t k is the transmission correction term, k is the compensation coefficient, and σ d is the standard deviation of the fiber channel transmission delay;

[0034] S1.2.4: Determine the total compensation amount: Based on the mean square error of all PMU nodes and the transmission correction term, determine the total compensation amount of the PMU channel, and based on the total compensation amount of the PMU channel, distribute the total compensation amount to each PMU node to dynamically adjust the corresponding local clock, specifically:

[0035]

[0036] Where: f new is the adjusted crystal oscillator output frequency, Δθ is the adjustment amount of the phase angle, and Δt PMU is the total compensation amount of the PMU channel, f0 is the initial nominal frequency of the crystal oscillator, T obs is the observation window time length, and T clkis the period of the crystal oscillator output clock.

[0037] Furthermore, the acquisition formula for the total compensation amount of the PMU channel is specifically:

[0038] Δt PMU = t s + t k

[0039] where: Δt PMU is the total compensation amount of the PMU channel, t s is the mean square error of all PMU nodes, t k is the transmission correction term.

[0040] Furthermore, the compensation coefficient is adjusted through the master station residual of the PMU channel, including:

[0041] M1: Determine the master station residual: According to the total compensation amount of the PMU channel and the actual clock deviation at the PMU node, obtain the master station residual, specifically:

[0042] ε = |Δt′ i - Δt PMU |

[0043] where: ε is the master station residual, Δt′ i is the actual clock deviation at the i-th PMU node, Δt PMU is the total compensation amount of the PMU channel, and i is the index of the PMU node;

[0044] M2: Adaptive adjustment: Compare the master station residual with the preset residual. When the master station residual is not greater than the preset residual, the compensation coefficient is not adjusted. Otherwise, the compensation coefficient is adjusted through a PID controller, specifically:

[0045]

[0046] where: k new is the adjusted compensation coefficient, k is the compensation coefficient, ε is the master station residual, K p is the proportional gain coefficient, K i is the integral gain coefficient, K d is the derivative gain coefficient, and t is the sampling time.

[0047] Furthermore, obtaining the corrected timestamp includes:

[0048] S1.3.1: Expand the time axis: According to the sampling period of the SCADA channel and the target time resolution, insert a preset number of virtual time points between two original sampling points on the SCADA channel. The acquisition formula for the preset number is specifically:

[0049]

[0050] Where: N vir is the preset number of virtual time points, T SCA is the sampling period of the SCADA channel, Δt tar is the target time resolution;

[0051] S1.3.2: Determine the interpolated timestamp: According to the magnitude of the preset number, obtain the corrected timestamp, specifically:

[0052] t′ SCA (t′) = t SCA (t′) + 0.1×N vir

[0053] Where: t′ SCA (t′) is the timestamp corresponding to time t′, t SCA (t′) is the corrected timestamp at time t′, N vir is the preset number of virtual time points, and t′ is the current time point of the virtual data.

[0054] Furthermore, according to the judgment formula of the steady state mode, when all the judgment conditions in the judgment formula are satisfied, the working condition mode is the steady state mode. According to the judgment formula of the transient mode, when any one of the judgment conditions in the judgment formula is satisfied, the working condition mode is the transient mode;

[0055] The judgment formula of the steady state mode is specifically:

[0056]

[0057] Where: ΔV is the voltage change, Δf is the frequency change, Δt is the time interval, and THD is the total harmonic distortion rate;

[0058] The judgment formula of the transient mode is specifically:

[0059]

[0060] Where: ΔV is the voltage change, Δt is the time interval, f dom is the dominant oscillation frequency.

[0061] Furthermore, when the working condition mode is the steady state mode, determine the virtual interpolation data and the transfer function, specifically:

[0062]

[0063] Where: x vir (t j) is the SCADA channel at time t j The virtual data value inserted at j-1 ) is the SCADA channel at time t j-1 The measured value at x(t j ) is the SCADA channel at time t j The measured value at T SCA is the sampling period of the SCADA channel, Δt vir is the time interval between the virtual timestamp and the previous actual sampling point, H(s) is the transfer function of the filter, S is the complex frequency variable, ω c is the cutoff frequency of the filter;

[0064] When the operating mode is a transient mode, virtual interpolation data, state equations, observation values ​​and adaptive noise are obtained, specifically:

[0065]

[0066] Where: x tua (t) is the wavelet interpolation in the SCADA channel at time t j The dummy data value inserted at m,o is the coefficient of the scaling function, φ m，o is the scale function, d m,o is the coefficient of the wavelet function, ψ m，o is the wavelet function, M is the maximum number of layers of wavelet decomposition, O is the translation parameter, m is the number of decomposition layers, x d is the state vector at time d, x d-1 is the state vector at time d-1, ν d-1 is the state change rate at time d-1, Δt is the time interval, ω k is the process noise, z d is the observed value at time d, ν d is the state change rate at time d, Q is the process noise covariance matrix, Q0 is the reference process noise, and ΔV is the voltage change.

[0067] Furthermore, the revised instruction issuance time is determined, including:

[0068] S3.1: Determine the total time error compensation: According to the working mode, equipment response time and compensation coefficient, determine the total time error that needs to be compensated, specifically:

[0069] Δ tcom =K mode ·(τ mode +σ delay )

[0070] Where: Δtcom is the total time error to be compensated, K mode is the compensation coefficient in the operating mode, τ mode is the average response time in the operating mode, σ delay is the standard deviation of the communication transmission delay;

[0071] S3.2: Determine the corrected timestamp: According to the total time error to be compensated and the expected execution time of the original control instruction, determine the corrected instruction issuance time, specifically:

[0072] t cmd = t cmd -Δt com

[0073] where: t′ cmd is the corrected instruction issuance time, t cmd is the expected execution time of the original control instruction, Δt com is the total time error to be compensated.

[0074] Furthermore, through the transient feature prediction module, voltage change amount, and time interval, determine the predicted compensation time, specifically:

[0075]

[0076] where: Δt pre is the predicted compensation time, ΔV is the voltage change amount, Δt is the time interval, T hor is the prediction time window, Δt th is the delay threshold.

[0077] Compared with the prior art, the beneficial effects of the present invention are:

[0078] First: The present invention compensates the local clock error of the PMU through GPS time and inserts the SCADA virtual timestamp to achieve the time axis alignment of multi-source data. At the same time, a preset number of virtual timestamps are inserted between SCADA sampling points, thereby narrowing the time resolution gap with PMU data, reducing the synchronization error, and providing an accurate time reference;

[0079] Second: Based on the voltage / frequency change rate and total harmonic distortion rate, the present invention divides the operating mode into steady state and transient modes. In the steady state mode, high-frequency noise is suppressed through linear interpolation and a low-pass filter. In the transient mode, transient features are captured through wavelet basis function interpolation and combined with the state equation to adaptively adjust the process noise;

[0080] Thirdly: The present invention dynamically calculates the compensation amount according to the working condition mode coefficient to correct the command issuance time. At the same time, in combination with the voltage change rate, it predicts the future delay, thereby obtaining the final command time, significantly improving the control real-time performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 is a schematic flow chart of the power grid operation delay control method under emergency conditions of the present invention;

[0082] Figure 2 is a comparison chart of the time compensation effect of PMU nodes of the present invention;

[0083] Figure 3 is a comparison chart of signal interpolation of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0084] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0085] Embodiment 1

[0086] Refer to Figures 1 - 3 , this embodiment provides a power grid operation delay control method under emergency conditions. The power grid operation delay control method under emergency conditions includes the following steps:

[0087] Step S1: Create a unified time baseline. That is, through the optical fiber channel and the dispatching data network, the PMU channel and the SCADA channel are respectively designed. At the same time, according to the designed PMU channel, time compensation of the PMU device is performed. According to the designed SCADA channel, the SCADA virtual time stamp is inserted, and synchronous time compensation is performed according to the time compensation of the PMU device and the inserted SCADA virtual time stamp. Specifically as follows:

[0088] Step S1.1: Design separate channels. That is, through the optical fiber channel, the PMU channel is designed. Through the dispatching data network, the SCADA channel is designed.

[0089] In this embodiment, data is transmitted through single-mode G.652.D optical fiber, and at the same time, a dual-ring redundant topology network architecture is adopted, where the number of nodes in the ring network does not exceed 32 to ensure that the self-healing time is less than 50 ms, the maximum transmission distance is set to 80 km, and the relay station interval does not exceed 40 km.

[0090] Furthermore, data transmission is carried out through a dual-star redundant topology structure, where the central node is set as the scheduling master station, with dual-machine hot standby and a switching time less than 200 ms. At the same time, the access node is set as the substation RTU, and the single-station bandwidth is not less than 10 Mbps.

[0091] Step S1.2: Perform time compensation for the PMU device. That is, at the local clock and GPS time at each PMU node, obtain the time difference at each PMU node, so as to determine the mean square error of all PMU nodes. At the same time, determine the correction term for optical fiber transmission through the standard deviation of the optical fiber channel transmission delay, and then determine the total compensation amount of the PMU channel, and adjust the corresponding local clock according to the determined total compensation amount of the PMU channel. Specifically as follows:

[0092] Step S1.2.1: Perform clock error compensation. That is, at each PMU node, obtain the local clock at a preset time, and at the same time, according to the obtained local clock and GPS time, obtain the time difference between the two, specifically:

[0093]

[0094] where: Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, is the timestamp recorded by the local clock at the i-th PMU node, is the GPS timestamp received at the i-th PMU node, and i is the index of the PMU node.

[0095] Step S1.2.2: Obtain the PMU node error. That is, according to the time difference at each PMU node obtained in Step S1.2.1, calculate and obtain the mean square error of all PMU nodes, specifically:

[0096]

[0097] where: t s is the mean square error of all PMU nodes, N is the total number of PMU nodes, and Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, and i is the index of the PMU node.

[0098] In the process of specific implementation, the total number of PMU nodes is set to 120, and at the same time, the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node is 2.1 μs That is to say, the mean square error of all PMU nodes is 4.41 μs 2 .

[0099] Step S1.2.3: Determine the delay correction. That is, according to the standard deviation of the optical fiber channel transmission delay, determine the correction term for optical fiber transmission, specifically:

[0100]

[0101] Where: t k is the transmission correction term, k is the compensation coefficient, and σ d is the standard deviation of the optical fiber channel transmission delay.

[0102] In the process of specific implementation, the standard deviation of the optical fiber channel transmission delay is 0.2 ms, and the compensation coefficient is set to 1.5. That is to say, the transmission correction term is 0.06 ms 2 .

[0103] Step S1.2.4: Determine the total compensation amount. That is, according to the mean square error of all PMU nodes obtained in step S1.2.2 and the transmission correction term obtained in step S1.2.3, determine the total compensation amount of the PMU channel, specifically:

[0104] Δt PMU = t s + t k

[0105] Where: Δt PMU is the total compensation amount of the PMU channel, t s is the mean square error of all PMU nodes, and t k is the transmission correction term.

[0106] Furthermore, according to the total compensation amount of the PMU channel obtained, distribute the total compensation amount to each PMU node to dynamically adjust the corresponding local clock. Specifically, according to the total compensation amount of the PMU channel obtained, adjust the crystal oscillator control voltage through a digital-to-analog converter, and adjust the phase angle through the MMCM (Mixed-Mode Clock Manager) module of the FPGA, specifically:

[0107]

[0108] Where: f new is the adjusted crystal oscillator output frequency, Δθ is the adjustment amount of the phase angle, Δt PMU is the total compensation amount of the PMU channel, f0 is the initial nominal frequency of the crystal oscillator, T obs is the observation window time length, and T clk is the period of the crystal oscillator output clock.

[0109] Step S1.3: Insert the SCADA virtual timestamp. That is, according to the sampling period of the SCADA channel and the target time resolution, determine the preset number of virtual time points to be inserted in the SCADA channel, and obtain the corrected timestamp. Specifically as follows:

[0110] Step S1.3.1: Expand the time axis. That is, according to the sampling period of the SCADA channel and the target time resolution, insert a preset number of virtual time points between two original sampling points on the SCADA channel. In this embodiment, the preset number of virtual time points is specifically:

[0111]

[0112] Where: N vir is the preset number of virtual time points, T SCA is the sampling period of the SCADA channel, Δt tar is the target time resolution.

[0113] Step S1.3.2: Determine the interpolated timestamp. That is, according to the preset number determined in Step S1.3.1, obtain the corrected timestamp, specifically:

[0114] t′ SCA (t′) = t SCA (t′) + 0.1 × N vir

[0115] Where: t′ SCA (t′) is the timestamp corresponding to the time t′, t SCA (t′) is the corrected timestamp at the time t′, N vir is the preset number of virtual time points, and t′ is the current time point of the virtual data.

[0116] Step S1.4: Perform synchronous time compensation. That is, set the temperature sensor at the crystal oscillator of the clock module to obtain the temperature data at the clock module, so as to determine the ambient temperature change rate. At the same time, according to the total compensation amount of the PMU channel determined in Step S1.2.4, the original uncompensated timestamp recorded by the local clock, and the ambient temperature change rate, obtain the synchronized timestamp, specifically:

[0117]

[0118] Where: t syn is the synchronized timestamp, t raw is the original uncompensated timestamp recorded by the local clock, Δt PMU is the total compensation amount of the PMU channel, δ is the crystal oscillator temperature drift coefficient, Δt obs is the observation time window, is the environmental temperature change rate.

[0119] Step S2: Dynamically select the interpolation filter bank. That is, according to the working condition mode and synchronous time compensation, data processing under different working condition modes is carried out. Specifically as follows:

[0120] Step S2.1: Determine the working condition mode. That is, through the SCADA channel designed in step S1.1, the corresponding voltage, frequency, and power analog quantities are obtained, and according to the PMU channel designed in step S1.1, the obtained voltage, frequency, and power analog quantities are aligned with the PMU time baseline. At the same time, through FFT analysis, the corresponding total harmonic distortion rate is obtained. And according to the voltage, frequency, and power analog quantities and the total harmonic distortion rate, the corresponding working condition mode is determined.

[0121] In this embodiment, the working condition modes include a steady state mode and a transient mode. Further, the judgment formula for the steady state mode is specifically:

[0122]

[0123] Where: ΔV is the voltage change amount, Δf is the frequency change amount, Δt is the time interval, and THD is the total harmonic distortion rate.

[0124] Further, the judgment formula for the transient mode is specifically:

[0125]

[0126] Where: ΔV is the voltage change amount, Δt is the time interval, and f dom is the dominant oscillation frequency.

[0127] That is to say, according to the judgment formula of the steady state mode, when all judgment conditions are met, the current working condition mode is the steady state mode. At the same time, according to the judgment formula of the transient mode, when any judgment condition is met, the current working condition mode is the transient mode.

[0128] Step S2.2: Steady state mode processing. That is, when the current working condition mode is the steady state mode, according to the preset number of virtual time points determined in step S1.3.1, the corresponding number of virtual interpolation data is inserted into the SCADA channel. At the same time, according to the complex frequency variable and the cut-off angular frequency of the filter, the transfer function of the low-pass filter is obtained.

[0129] In this embodiment, the acquisition formulas for the virtual interpolation data and the transfer function are specifically:

[0130]

[0131] Where: x vir (t j) is the virtual data value inserted at time t in the SCADA channel j at this point, x(t j-1 ) is the measured value at time t in the SCADA channel j-1 at this point, x(t j ) is the measured value at time t in the SCADA channel j at this point, T SCA is the sampling period of the SCADA channel, Δt vir is the time interval between the virtual timestamp and the previous actual sampling point, H(s) is the transfer function of the filter, S is the complex frequency variable, ω c is the cut-off angular frequency of the filter.

[0132] In the process of specific implementation, the measured value at 1.0 s in the SCADA channel is 220 kV, and the measured value at 2.0 s in the SCADA channel is 218 kV. Then the virtual data value inserted at 1.5 s in the SCADA channel is: 220+(218 - 220) / 1*0.5 = 219 kV.

[0133] Step S2.3: Transient mode processing. That is, when the current working mode is the transient mode, a corresponding number of virtual interpolation data are inserted into the SCADA channel through the wavelet basis function and the preset number of virtual time points determined in step S1.3.1. At the same time, the corresponding state equation and observation value are determined through the constructed state variable and observation model. And the process noise is adaptively adjusted according to the reference noise.

[0134] In this embodiment, the acquisition formulas for virtual interpolation data, state equation, observation value, and adaptive noise are specifically:

[0135]

[0136] Where: x tua (t) is the virtual data value inserted at time t in the SCADA channel through wavelet interpolation j at this point, c m，o is the coefficient of the scaling function, φ m，o is the scaling function, d m，o is the coefficient of the wavelet function, ψ m，o is the wavelet function, M is the maximum number of layers of wavelet decomposition, O is the translation parameter, m is the decomposition layer, x d is the state vector at time d, x d-1 is the state vector at time d - 1, ν d-1 is the state change rate at time d - 1, Δt is the time interval, ω k is the process noise, z d is the observation value at time d, ν dis the state change rate at time d, Q is the process noise covariance matrix, Q0 is the reference process noise, and ΔV is the voltage change amount.

[0137] Step S3: Determine the control instruction timestamp error. That is, based on the working condition mode, device response time, and compensation coefficient, determine the total time error that needs to be compensated, and at the same time, based on the total time error that needs to be compensated, determine the corrected instruction issuance time. Specifically as follows:

[0138] Step S3.1: Determine the total time error compensation. That is, based on the working condition mode determined in Step S2.1, determine the corresponding transient parameters and steady-state parameters. Further, when the working condition mode changes from the steady-state mode to the transient mode, obtain the transient parameters, that is, the device response time and compensation coefficient in the transient mode. When the working condition mode changes from the transient mode to the steady-state mode, obtain the steady-state parameters, that is, the device response time and compensation coefficient in the steady-state mode. Specifically, the compensation coefficient in the transient mode is set to 1.2, and the compensation coefficient in the steady-state mode is set to 0.95.

[0139] In this embodiment, based on the obtained device response time and compensation coefficient, determine the total time error that needs to be compensated, specifically:

[0140] Δt com =K mode ·(τ mode +σ delay )

[0141] Where: Δt com is the total time error that needs to be compensated, K mode is the compensation coefficient in the working condition mode, τ mode is the average response time in the working condition mode, and σ delay is the standard deviation of the communication transmission delay.

[0142] In the process of specific implementation, the standard deviation of the communication transmission delay is 0.2 ms, and at the same time, the device response time in the transient mode is 2 ms, and the device response time in the steady-state mode is 10 ms. That is to say, the total time error that needs to be compensated in the steady-state mode is: 0.95 * (10 + 0.2) = 9.69 ms, and the total time error that needs to be compensated in the transient mode is: 1.2 * (2 + 0.2) = 2.64 ms.

[0143] Further, in this embodiment, the standard deviation of the optical fiber transmission delay is 0.2 ms, and the operation time of the circuit breaker is 8 ms, if the data processing delay time is 1.5 ms, the total time error is: 0.2 + 8 + 1.5 = 9.7 ms. That is to say, combined with the total time error to be compensated in the steady state mode, the total time error is reduced from 9.7 ms to 9.69 ms, that is, the overall decrease is 0.01 ms. Similarly, combined with the total time error to be compensated in the transient mode, the total time error is reduced from 9.7 ms to 2.64 ms, that is, the overall decrease is 7.06 ms.

[0144] Step S3.2: Determine the corrected timestamp. That is, according to the total time error to be compensated obtained in Step S3.1 and the expected execution time of the original control instruction, determine the corrected instruction issuance time, specifically:

[0145] t cmd = t cmd -Δt com

[0146] where: t′ cmd is the corrected instruction issuance time, t cmd is the expected execution time of the original control instruction, Δt com is the total time error to be compensated.

[0147] In the process of specific implementation, the expected execution time of the original control instruction is 100 ms. That is to say, the corrected timestamp in the steady state mode is: 100 - 9.69 = 90.31 ms, and the corrected timestamp in the transient mode is: 100 - 2.64 = 97.36 ms.

[0148] Step S4: Timestamp advanced correction. That is, according to the corrected timestamp determined in Step S3.2, combine it with the predicted compensation time to determine the final instruction issuance time, specifically:

[0149] t″ cmd = t′ cmd -Δt pre

[0150] where: t″ cmd is the finally corrected instruction issuance time, t′ cmd is the corrected instruction issuance time, Δt pre is the predicted compensation time.

[0151] In this embodiment, the predicted compensation time is determined through the transient feature prediction module, the voltage change amount, and the time interval, specifically:

[0152]

[0153] where: Δt preTo predict the compensation time, ΔV is the voltage change, Δt is the time interval, and T hor is the prediction time window, and Δt th is the delay threshold.

[0154] In the process of specific implementation, the prediction compensation time is set to 10 ms. That is to say, the finally corrected timestamp in the steady state mode is: 90.31 - 10 = 80.31 ms, and the finally corrected timestamp in the transient state mode is: 97.36 - 10 = 87.36 ms .

[0155] Reference Figure 2 , Figure 2 is the comparison chart of the time compensation effect of the PMU node. It can be seen from Figure 2 that the deviations of each PMU node are randomly distributed within the range of ±4 μs (standard deviation 2 μs), conforming to the characteristics of normal distribution.

[0156] Reference Figure 3 , Figure 3 is the interpolation comparison chart of the signal. The upper part is the interpolation comparison chart of the steady state signal, and the lower part is the interpolation comparison chart of the transient signal. It can be seen from Figure 3 that for the steady state signal, the original signal feature is a 50 Hz sine wave, and the SCADA device samples at 1 Hz. While in this application, the sparse sampling data is interpolated and reconstructed through a low-pass filter, making its variance decrease from 0.048 to 0.012, effectively suppressing high-frequency noise. For the transient signal, the original signal feature superimposes a 100 Hz oscillation at 0.5 s, and in this application, high-frequency components are captured through wavelet interpolation, making the variance decrease from 0.217 to 0.012. That is to say, the accuracy is improved by 60.8% in the transient state mode, and the 100 Hz oscillation waveform is successfully restored.

[0157] Embodiment 2

[0158] This embodiment provides a method for controlling the grid operation delay in case of emergencies. The specific implementation method is the same as that of Embodiment 1, and the difference is that during the process of determining the optical fiber transmission correction term, the compensation coefficient is continuously adjusted. The following is an example of the present invention in combination with the specific implementation manner of this embodiment.

[0159] In this embodiment, the compensation coefficient is continuously adjusted. That is to say, the compensation coefficient is adjusted through the master station residual of the PMU channel. Specifically as follows:

[0160] Step M1: Determine the master station residual. That is, according to the total compensation amount of the PMU channel obtained in step S1.2.4, compare it with the actual clock deviation obtained at the PMU node to obtain the master station residual. Specifically:

[0161] ε = |Δt′ i - Δt PMU |

[0162] Where: ε is the master station residual, Δt′ i is the actual clock deviation at the i-th PMU node, Δt PMU is the total compensation amount of the PMU channel, and i is the index of the PMU node.

[0163] Step M2: Adaptive adjustment. That is, compare the master station residual obtained in Step M1 with a preset residual (which can be specifically set according to actual operation, so it is not specifically elaborated in this embodiment). When the master station residual is not greater than the preset residual, the compensation coefficient is not adjusted. Otherwise, the compensation coefficient is adjusted through a PID controller. In this embodiment, the adjusted compensation coefficient is specifically:

[0164]

[0165] Where: k new is the adjusted compensation coefficient, k is the compensation coefficient, ε is the master station residual, K p is the proportional gain coefficient, K i is the integral gain coefficient, K d is the derivative gain coefficient, and t is the sampling time.

[0166] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended embodiments and their equivalents.

Claims

1. A method for controlling the time delay of power grid operation under sudden conditions, characterized in that, It includes: S1: Create a unified time baseline: According to the designed PMU channels and SCADA channels, perform time compensation for the PMU device, insert SCADA virtual timestamps, and perform synchronous time compensation. S2: Dynamically select an interpolation filter bank: According to the operating mode and synchronous time compensation, perform data processing under different operating modes, including: S2.1: Determine the operating mode: According to voltage, frequency, power analog quantities, and total harmonic distortion rate, determine the corresponding operating mode, and the operating mode includes a steady-state mode and a transient mode. S2.2: Steady-state mode processing: When the operating mode is the steady-state mode, determine the virtual interpolation data in the SCADA channel according to the preset number of virtual time points, and at the same time, according to the complex frequency variable and the cut-off angular frequency of the filter, obtain the transfer function of the low-pass filter. S2.3: Transient mode processing: When the operating mode is the transient mode, determine the virtual interpolation data in the SCADA channel through wavelet basis functions and the preset number of virtual time points, and at the same time, according to the constructed state variables and observation model, determine the state equation and observation values, and adaptively adjust the process noise according to the reference noise. S3: Determine the control instruction timestamp error: According to the operating mode, device response time, and compensation coefficient, determine the total time error to be compensated, and at the same time, according to the total time error to be compensated, determine the corrected instruction issuing time. S4: Timestamp lead correction: Combine the corrected instruction issuing time and the predicted compensation time to determine the final instruction issuing time, specifically: t″ cmd = t′ cmd - Δt pre where: t″ cmd is the time when the finally corrected instruction is issued, t′ cmd is the time when the corrected instruction is issued, Δt pre is the predicted compensation time.

2. The method for controlling the grid operation delay in case of emergencies according to claim 1, characterized in that, Perform synchronous time compensation, including: S1.1: Design separate channels: Design PMU channels and SCADA channels respectively through fiber optic channels and dispatching data networks. S1.2: Perform time compensation for the PMU device: Determine the mean square error of all PMU nodes through the local clock and GPS time at each PMU node, and obtain the total compensation amount of the PMU channel according to the correction term of fiber optic transmission, and adjust the local clock. S1.3: Insert SCADA virtual timestamps: According to the sampling period and target time resolution of the SCADA channel, determine the preset number of virtual time points inserted in the SCADA channel, and obtain the corrected timestamp. S1.4: Perform synchronous time compensation: Obtain the ambient temperature change rate of the clock module through a temperature sensor, and at the same time, according to the total compensation amount of the PMU channel, the original uncompensated timestamp recorded by the local clock, and the ambient temperature change rate, obtain the synchronized timestamp, specifically: Where: t syn is the timestamp after synchronization, t raw is the original uncompensated timestamp recorded by the local clock, Δt PMU is the total compensation amount of the PMU channel, δ is the crystal oscillator temperature drift coefficient, Δt obs is the observation time window, is the environmental temperature change rate.

3. A method for controlling the time delay of power grid operation in case of emergencies according to claim 2, characterized in that, Adjust the local clock, including: S1.2.1: Perform clock error compensation: At each PMU node, obtain the local clock, and obtain the time difference according to the local clock and GPS time, specifically: where: Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, is the timestamp recorded by the local clock at the i-th PMU node, is the GPS timestamp received at the i-th PMU node, and i is the index of the PMU node; S1.2.2: Obtain PMU node errors: According to the time differences at each PMU node, obtain the mean square error of all PMU nodes, specifically: where: t s is the mean square error of all PMU nodes, N is the total number of PMU nodes, and Δt i is the instantaneous deviation between the local clock and the GPS reference clock at the i-th PMU node, and i is the index of the PMU node; S1.2.3: Determine delay correction: Determine the correction term of fiber optic transmission according to the standard deviation of the fiber optic channel transmission delay, specifically: where: t k is the transmission correction term, k is the compensation coefficient, and σ d is the standard deviation of the transmission delay of the optical fiber channel; S1.2.4: Determine the total compensation amount: Based on the mean square error of all PMU nodes and the transmission correction term, determine the total compensation amount of the PMU channel, and distribute the total compensation amount to each PMU node according to the total compensation amount of the PMU channel to dynamically adjust the corresponding local clock. Specifically: Where: f new is the adjusted crystal oscillator output frequency, Δθ is the adjustment amount of the phase angle, and Δt PMU is the total compensation amount of the PMU channel, f0 is the initial nominal frequency of the crystal oscillator, and T obs is the observation window time length, and T clk is the period of the crystal oscillator output clock.

4. A method for controlling the delay of power grid operation in case of emergencies according to claim 3, characterized in that The formula for obtaining the total compensation amount of the PMU channel is specifically: Δt PMU = t s + t k where: Δt PMU is the total compensation amount of the PMU channel, t s is the mean square error of all PMU nodes, t k is the transmission correction term.

5. A method for controlling the delay of power grid operation in case of emergencies according to claim 3, characterized in that, Adjust the compensation coefficient through the master station residual of the PMU channel, including: M1: Determine the master station residual: Based on the total compensation amount of the PMU channel and the actual clock deviation at the PMU node, obtain the master station residual. Specifically: ε = |Δt′ i - Δt PMU | Where: ε is the master station residual, Δt′ i is the actual clock deviation at the i-th PMU node, Δt PMU is the total compensation of the PMU channel, and i is the index of the PMU node; M2: Adaptive adjustment: Compare the master station residual with the preset residual. When the master station residual is not greater than the preset residual, the compensation coefficient is not adjusted. Otherwise, adjust the compensation coefficient through a PID controller. Specifically: Where: k new is the adjusted compensation coefficient, k is the compensation coefficient, ε is the master station residual, K p is the proportional gain coefficient, K i is the integral gain coefficient, K d is the derivative gain coefficient, and t is the sampling time.

6. A method for controlling the time delay of power grid operation in case of emergencies according to claim 2, characterized in that, Obtain the corrected timestamp, including: S1.3.1: Expand the time axis: Based on the sampling period of the SCADA channel and the target time resolution, insert a preset number of virtual time points between two original sampling points on the SCADA channel. The formula for obtaining the preset number is: Where: N vir is the preset number of virtual time points, T SCA is the sampling period of the SCADA channel, Δt tar is the target time resolution; S1.3.2: Determine the interpolated timestamp: Based on the size of the preset number, obtain the corrected timestamp. Specifically: t′ SCA (t′) = t SCA (t′) + 0.1×N vir where: t' SCA (t') is the time stamp corresponding to time t', t SCA (t') is the corrected time stamp at time t', N vir is the preset number of virtual time points, and t' is the current time point of the virtual data.

7. A method for controlling the grid operation delay in case of emergencies according to claim 1, characterized in that, According to the judgment formula of the steady state mode, when all judgment conditions in the judgment formula are satisfied, the operating mode is the steady state mode. According to the judgment formula of the transient mode, when any judgment condition in the judgment formula is satisfied, the operating mode is the transient mode; The judgment formula of the steady state mode is specifically: Where: ΔV is the voltage change, Δf is the frequency change, Δt is the time interval, and THD is the total harmonic distortion rate; The judgment formula of the transient mode is specifically: Where: ΔV is the voltage change, Δt is the time interval, and f dom is the dominant oscillation frequency.

8. A method for controlling the time delay of power grid operation in case of emergencies according to claim 1 or 7, characterized in that, When the operating mode is the steady state mode, determine the virtual interpolation data and the transfer function. Specifically: where: x vir (t j ) is the virtual data value inserted at time t j in the SCADA channel, x(t j-1 ) is the measured value at time t j-1 in the SCADA channel, x(t j ) is the measured value at time t j in the SCADA channel, T SCA is the sampling period of the SCADA channel, Δt vir is the time interval between the virtual timestamp and the previous actual sampling point, H(s) is the transfer function of the filter, S is the complex frequency variable, ω c is the cut-off angular frequency of the filter; When the operating mode is the transient mode, obtain the virtual interpolation data, the state equation, the observed value, and the adaptive noise. Specifically: where: x tua (t) is the virtual data value inserted at time t through wavelet interpolation in the SCADA channel j , c m，o is the coefficient of the scaling function, φ m，o is the scaling function, d m，o is the coefficient of the wavelet function, ψ m，o is the wavelet function, M is the maximum number of levels of wavelet decomposition, O is the translation parameter, m is the decomposition level, x d is the state vector at time d, x d-1 is the state vector at time d - 1, v d-1 is the state change rate at time d - 1, Δt is the time interval, ω k is the process noise, z d is the observation value at time d, v d is the state change rate at time d, Q is the process noise covariance matrix, Q0 is the reference process noise, and ΔV is the voltage change amount.

9. A method for controlling the time delay of power grid operation in case of emergencies according to claim 1, characterized in that Determine the corrected command issuance time, including: S3.1: Determine the total time error compensation: Based on the operating mode, the device response time, and the compensation coefficient, determine the total time error to be compensated. Specifically: Δt com = K mode ·(τ mode + σ delay ) where: Δt com is the total time error to be compensated, K mode is the compensation coefficient under the operating mode, τ mode is the average response time under the operating mode, σ delay is the standard deviation of the communication transmission delay; S3.2: Determine the corrected timestamp: Based on the total time error to be compensated and the expected execution time of the original control command, determine the corrected command issuance time. Specifically: t′ cmd = t cmd - Δt com Where: t' cmd is the time when the corrected instruction is issued, t cmd is the expected execution time of the original control instruction, Δt com is the total time error that needs to be compensated.

10. A method for controlling the delay of power grid operation in case of emergencies according to claim 1, characterized in that, Determine the predicted compensation time through the transient feature prediction module, the voltage change amount, and the time interval. Specifically: where: Δt pre is the prediction compensation time, ΔV is the voltage change, Δt is the time interval, T hor is the prediction time window, Δt th is the delay threshold.

Images