An abnormality monitoring and early warning method based on low-voltage transformer area power consumption behavior data

Abstract

The application discloses an abnormality monitoring and early warning method based on low-voltage area power consumption behavior data, and belongs to the technical field of power system monitoring. Through the collection of multi-dimensional time sequence power consumption data of a low-voltage area, the data is decomposed into multi-scale components by using variational mode decomposition, an input feature matrix is constructed in combination with time coding, holidays and meteorological characteristics, and a typical power consumption behavior benchmark model is generated by training a CNN-GRU hybrid deep learning model. Real-time residuals are calculated based on the benchmark model, an adaptive detection threshold is dynamically generated by using a peak value exceeding threshold method based on extreme value theory, and effective abnormal events are confirmed in combination with multi-parameter physical correlation cross-validation rules. Through a cloud-edge collaborative architecture, abnormal characteristics are uploaded to the cloud, fault type diagnosis and severity assessment are performed based on a pre-built fault characteristic knowledge base, a three-level grading early warning signal is generated, and protection actions are linked and executed, so that intelligent and accurate monitoring and early warning of low-voltage area electrical faults are realized.

Description

Technical Field

[0001] This invention relates to the field of power system monitoring and control technology, specifically to an abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data, which is particularly suitable for intelligent online monitoring and early fault warning of low-voltage distribution areas at the end of the distribution network. Background Technology

[0002] Low-voltage distribution areas (areas supplied by 380V / 220V distribution transformers) are the "last mile" of the power system facing end users, and their operational status directly determines the power quality and safety for hundreds of millions of users. However, low-voltage distribution areas generally suffer from the following technical challenges:

[0003] 1. Complex Fault Mechanisms and Weak Symptoms: Electrical faults in low-voltage distribution areas exhibit diverse characteristics, including overload, three-phase imbalance, single-phase grounding, harmonic pollution, and intermittent partial discharge caused by equipment insulation aging. In their early stages, these faults often manifest as weak voltage fluctuations, current distortions, or a slow temperature rise. Traditional relay protection devices based on fixed thresholds (such as overcurrent protection and overvoltage protection) cannot trigger operation in the early stages of fault development. By the time the parameters exceed the set limits, equipment damage or safety accidents have often already occurred.

[0004] 2. Load behavior is dynamic and variable, making threshold determination difficult: The load in a distribution area is affected by factors such as user behavior habits, seasonal changes, holidays, and weather, exhibiting strong nonlinear and non-stationary characteristics. The normal current fluctuation range of the same distribution area at different times may differ by several times. Traditional fixed threshold methods either miss early faults due to excessively high threshold settings or frequently generate false alarms during normal load fluctuations due to excessively low threshold settings, seriously affecting operation and maintenance efficiency.

[0005] 3. The coexistence of data silos and massive data: On the one hand, most distribution areas lack effective monitoring methods, and data collection is incomplete and discontinuous, forming data silos; on the other hand, with the popularization of smart meters and IoT terminals, some distribution areas have accumulated massive amounts of time-series data, but due to the lack of efficient intelligent analysis algorithms and collaborative computing architecture, the value of this data has not been fully explored.

[0006] 4. Inaccurate alarms and lack of a tiered alarm mechanism: Most existing alarm systems are simple "exceed limits and then alarm" modes, which cannot distinguish between temporary disturbances, recoverable over-limits, and true early signs of failure. They lack a tiered alarm mechanism based on the severity and urgency of the fault, leading to "alarm fatigue" for maintenance personnel in the face of a massive number of invalid alarms.

[0007] To address the aforementioned issues, existing technologies have significant shortcomings. For example, a Chinese patent application discloses a power system fault monitoring method, but it only uses a single statistical feature for judgment, resulting in poor adaptability to changes in operating conditions. Another patent application uses neural networks for fault diagnosis, but it fails to address the issues of adaptive model updates and computational resource optimization, leading to problems such as poor real-time performance and high false alarm rates in practical deployments.

[0008] Therefore, there is an urgent need for an intelligent monitoring method that can deeply understand the electricity consumption behavior patterns of transformer substations, adaptively capture abnormal signs, and achieve accurate diagnosis and graded early warning. Summary of the Invention

[0009] The purpose of this section is to outline some aspects of embodiments of the present invention and to briefly describe some preferred embodiments. Simplifications or omissions may be made in this section, as well as in the abstract and title of this application, to avoid obscuring the purpose of these documents; however, such simplifications or omissions should not be construed as limiting the scope of the invention.

[0010] In view of the above-mentioned problems, the present invention is proposed.

[0011] Specifically, the present invention provides the following technical solution: an abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data, comprising the following steps: S1: collecting multi-dimensional time-series electricity consumption data of low-voltage distribution areas within a preset historical period, wherein the multi-dimensional time-series electricity consumption data includes three-phase voltage, three-phase current, active power, reactive power, power factor, frequency, zero-sequence current, and ambient temperature; S2: preprocessing the multi-dimensional time-series electricity consumption data to construct a sample dataset, and based on the sample dataset, using variational mode decomposition to decompose the time-series data of each dimension into several intrinsic mode function components, and combining the components with time coding features, holiday identifiers, and meteorological features to form an input feature matrix, and using a hybrid deep learning model of convolutional neural network and gated recurrent unit to train the input feature matrix, generating... S3: Collect the current time-series electricity consumption data of the low-voltage distribution area in real time, perform the same VMD decomposition and feature matrix construction as in S2, and input it into the typical electricity consumption behavior benchmark model. Calculate the difference between the expected normal value output by the model and the actual observed value to generate a multidimensional residual sequence. S4: For the multidimensional residual sequence, use the peak value exceeding threshold method based on extreme value theory to dynamically calculate the detection threshold of each dimension. When the residual of a certain dimension exceeds the dynamic threshold, it is marked as an initial detection anomaly. The initial detection anomaly is verified by multi-parameter cross-validation rules to confirm the valid anomaly event. S5: Upload the feature data of the confirmed valid anomaly event to the cloud. The cloud performs pattern matching and fault type diagnosis based on the pre-built fault feature knowledge base, assesses the severity of the fault, and generates a graded early warning signal.

[0012] Preferably, in step S2, when using VMD for multi-scale time-frequency decomposition, the number of decomposition modes K=8 and the quadratic penalty factor α=2000 are set, and each IMF component and its center frequency are solved iteratively by the alternating direction method of multipliers.

[0013] Preferably, in step S2, the CNN-GRU hybrid model consists of two one-dimensional convolutional layers and two max pooling layers in its CNN module. The first convolutional layer has 64 kernels and a kernel size of 12, and the second convolutional layer has 128 kernels and a kernel size of 8. The GRU module consists of one bidirectional GRU layer and one unidirectional GRU layer. The bidirectional GRU has 128 hidden units, and the unidirectional GRU has 64 hidden units.

[0014] Preferably, the dynamic threshold generation method based on POT in step S4 includes: setting a sliding window length L = 1440 time steps, extracting peak sequences where the absolute value of the residuals within the window exceeds the initial threshold u0 of the 90th percentile, fitting the shape parameter ξ and scale parameter β of the generalized Pareto distribution using maximum likelihood estimation, and then calculating the threshold based on a preset false alarm rate p. FA =1 / 1440 Reverse-push state threshold T dynamic .

[0015] Preferably, the multi-parameter cross-validation rule in step S4 includes: calculating the correlation confidence C when the normalized residuals of at least two electrical parameters simultaneously exceed their respective dynamic thresholds, and these parameters are marked as correlated in a preset physical correlation rule table. ij If C ij Greater than threshold C th If the value is 1.0, it is confirmed as a valid abnormal event; if only a single parameter is abnormal, it is determined to be noise interference.

[0016] Preferably, the fault feature knowledge base in step S5 contains fingerprint templates for at least eight types of typical faults. Each fault template consists of a feature mean vector, a feature standard deviation vector, and a matching similarity threshold. The feature vector includes 36 dimensions of features, including residual sign pattern, relative amplitude, time-series waveform features, and post-fault steady-state changes.

[0017] Preferably, the graded early warning signals in step S5 include three levels: Level 1 general attention warning corresponds to events with Severity < 0.3, which is logged and displayed on the monitoring screen; Level 2 important attention warning corresponds to events with Severity ≤ 0.3 < 0.7, which is sent to maintenance personnel via SMS or APP push notification; Level 3 emergency alarm corresponds to events with Severity ≥ 0.7, which triggers audible and visual alarms, automatically dials a phone number to notify, and links the intelligent circuit breaker to perform trip protection.

[0018] Preferably, step S5 further includes fault elimination confirmation and closed-loop feedback: after tripping, the residual sequence is continuously monitored, and if the residual recovers to below the threshold, it is confirmed that the fault has been isolated; the complete alarm and processing data is recorded to the cloud for periodic updates to the fault feature knowledge base and retraining of the benchmark model.

[0019] This invention provides an anomaly monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data, which has the following beneficial effects:

[0020] 1. By using VMD to adaptively decompose non-stationary power consumption signals into multiple narrowband components, the inherent laws of load fluctuations are revealed at multiple scales (from minute-level to quarter-level). Through the CNN-GRU hybrid model, CNN extracts local short-term patterns, while GRU captures long-term temporal dependencies. The two complement each other to achieve accurate modeling of complex power consumption behavior, and can capture early weak anomalies that traditional methods cannot identify, thus significantly reducing the false negative rate.

[0021] 2. A dynamic threshold generation method based on POT extreme value theory can automatically adjust the detection threshold according to the distribution characteristics of residual data over the past 24 hours. When the load fluctuates significantly, the threshold automatically increases to avoid false alarms; when the load is stable, the threshold automatically decreases to maintain sensitivity, fundamentally solving the technical problem of frequent false alarms or missed alarms under different operating conditions with fixed thresholds. A multi-parameter cross-validation mechanism further filters out false alarms caused by sensor noise and normal load fluctuations through physical association rules, ensuring high reliability of alarms. The overall false alarm rate can be reduced to less than 1 / 10 of that of traditional methods.

[0022] 3. By deploying computationally intensive but real-time-critical VMD decomposition, deep learning inference, and POT threshold calculation at the terminal edge layer, millisecond-level real-time response is achieved, enabling autonomous anomaly detection completely independent of network dependencies. Interaction with the cloud is only required when uploading anomaly events, with normal communication volume approaching zero, significantly saving 4G network bandwidth and data costs (estimated monthly data usage per unit <50MB). The cloud handles computationally intensive but less real-time-critical fault diagnosis and model update tasks, achieving optimal allocation of computing resources.

[0023] 4. The cloud-based fault characteristic knowledge base covers 36-dimensional feature fingerprints of eight common fault types. Using Mahalanobis distance for pattern matching, it can automatically diagnose fault types and assess severity, providing maintenance personnel with clear fault causes and handling suggestions. A three-tiered early warning mechanism adopts differentiated response measures based on fault severity, ranging from logging for minor faults to automatic tripping for severe ones, avoiding overreaction while ensuring rapid handling in emergencies.

[0024] 5. Through a complete closed loop of "monitoring → alarm → diagnosis → processing → feedback → optimization", the system can continuously update the knowledge base and optimize model parameters using new data from each fault handling. The detection performance continuously improves over time and can adapt to long-term changes in the power consumption behavior of the distribution area (such as the increase or decrease of users, equipment aging, etc.), avoiding the maintenance burden of manual parameter adjustment required by traditional methods. Attached Figure Description

[0025] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0026] Figure 1 The flowchart of the abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data provided by the present invention is shown. Detailed Implementation

[0027] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the protection scope of the present invention.

[0028] To address the problems existing in the prior art, this invention provides an anomaly monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data, comprising the following five steps:

[0029] S1: Multi-dimensional time-series electricity consumption data acquisition and preprocessing

[0030] S1.1 Data Acquisition Scheme

[0031] Distributed intelligent monitoring terminals are installed on the low-voltage outgoing lines and key branch lines of the low-voltage distribution transformers in the low-voltage distribution area. These terminals integrate the following sensor modules:

[0032] ① Voltage acquisition module

[0033] A resistor divider network combined with a precision differential amplifier is used to convert the three-phase voltage (U) a U b U cThis function converts 380V / 220V levels to 0-3.3V low-level signals that can be sampled by the ADC. The sampling chip uses a 24-bit high-precision Sigma-Delta ADC (such as Analog Devices' ADE9000), with a sampling frequency set to 6.4kHz.

[0034] Basis for sampling frequency setting: According to the design of this invention, it is necessary to analyze up to the 31st harmonic. The frequency of the 31st harmonic is 31 × 50 Hz = 1550 Hz. According to the Nyquist sampling theorem, to avoid spectral aliasing, the sampling frequency must be no less than twice the highest frequency of the signal, i.e., f0. s ≥ 2 × 1550Hz = 3100Hz. This invention uses 6.4kHz, which is about twice the theoretical minimum sampling rate, leaving sufficient margin to ensure the accuracy of higher harmonic analysis.

[0035] ② Current acquisition module

[0036] Using Rogowski coil current transformers or open-type current transformers, the three-phase current (I a I b I c The current transformer converts the zero-sequence current (I0) into a low-voltage signal proportional to the primary current. The range of the current transformer is selected based on the rated capacity of the transformer in the distribution area, generally taking 1.2 times the rated current as the full-scale range. The measurement accuracy class is 0.5, meaning that the measurement error at the rated current does not exceed ±0.5%. After passing through the signal conditioning circuit (including anti-aliasing filter and amplification circuit), it is input to the ADC.

[0037] ③ Power and electrical energy calculation

[0038] Dedicated metering chips such as the ADE9000 perform digital integration on instantaneous sampled values ​​of voltage and current to calculate active power P, reactive power Q, apparent power S, power factor cosφ, and frequency f in real time. These parameters are read by the MCU via SPI or UART interfaces.

[0039] The calculation principles for each power parameter are as follows:

[0040] Active power P: , where u(n) and i(n) are the nth instantaneous sampled values ​​of voltage and current, respectively, and N is the number of sampling points within one power frequency cycle (20ms);

[0041] Apparent power S: S=U rms ·I rms U rms and I rms These are the effective values ​​of the voltage and current;

[0042] Reactive power Q: Q = (S 2-P 2 ) 0.5 (Undistorted waveform) or calculated by multiplying the current and voltage after a 90° phase shift (distorted waveform).

[0043] Power factor cosφ: cosφ = P / S;

[0044] Frequency f: The period T is obtained by measuring the time difference between two voltage zero crossings, then f = 1 / T.

[0045] ④ Temperature monitoring module

[0046] The DS18B20 digital temperature sensor or the PT100 platinum resistance temperature sensor is used, installed on key heat-generating components such as the surface of the distribution transformer tank, outgoing cable joints, and switch contacts, to monitor the equipment's operating temperature T. The DS18B20 has a temperature measurement range of -55℃ to +125℃ and an accuracy of ±0.5℃, communicating with the MCU via a single-bus protocol. The PT100 has a wider temperature measurement range (-200℃ to +850℃) and higher accuracy (±0.1℃), but requires additional analog signal conditioning circuitry and an ADC.

[0047] ⑤ Supplementary external environment data

[0048] The station terminal uses a built-in 4G communication module to periodically (every 30 minutes) obtain real-time temperature, humidity, and weather conditions (sunny / rainy / snowy, etc.) from the meteorological service platform API.

[0049] ⑥ The data acquisition frequency is set as shown in Table 1 below:

[0050] Table 1: Data Acquisition Frequency Setting Table

[0051]

[0052] S1.2 Data Preprocessing (Data Cleaning)

[0053] The collected raw data has three types of quality problems, requiring targeted cleaning:

[0054] ① Handling missing values

[0055] The causes of missing values ​​include: brief interruption of the communication link, failure to sample during terminal device restart, and occasional sensor failure.

[0056] For short-lived missing values ​​(continuous missing time ≤ 15 minutes), segmented cubic Hermite interpolation is used for filling.

[0057] Suppose there are two known data points (t0, y0) and (t1, y1) before and after the missing point, and the derivative values ​​y'0 and y'1 of these two points are known (obtained by approximation through the difference between adjacent points: y0' = (y0 - y1) / (t0 - y1) / (t1 ... ﹣1 ) / (t0-t ﹣1 ), where (t) ﹣1 y ﹣1 (This refers to the data point preceding t0).

[0058] Then the interpolation y(t) at any time t within the missing interval is:

[0059] y(t)=h 00 (s)y0+h 10 (s)(t1-t0)y0′+h 01 (s)y1+h 11 (s)(t1-t0)y1′;

[0060] Where s = (t - t0) / (t1 - t0) is the normalization parameter, representing the relative position of t in the interval [t0, t1], with a value range of 0 to 1.

[0061] Hermite basis functions h 00 (s), h 10 (s), h 01 (s), h 11 The expressions for (s) are as follows:

[0062] h 00 (s)=2s 3 -3s 2 +1; h 10 (s)=s 3 -2s 2 +s;h 01 (s)=﹣2s 3 +3s 2 h 11 (s)=s 3 -s 2 ;

[0063] The meanings of each parameter are as follows:

[0064] t: Any time interval within the missing interval to be interpolated (unit: minutes);

[0065] t0: The moment when the last known data point is before the start of the missing interval;

[0066] t1: The time when the first known data is obtained after the end of the missing interval;

[0067] y0: Electrical parameter value (such as phase A voltage value) corresponding to time t0;

[0068] y1: The electrical parameter value corresponding to time t1;

[0069] s: The relative position of t in the interval [t0, t1], s=0 corresponds to time t0, s=1 corresponds to time t1;

[0070] y0′ and y1′ are the first derivatives at times t0 and t1, respectively, and are approximated by the difference between adjacent points.

[0071] This invention chooses Hermite interpolation instead of simple linear interpolation because Hermite interpolation not only ensures that the interpolation curve passes through the known data points, but also ensures that the derivative is continuous at the known points, making the interpolation curve smoother and closer to the physical characteristics of the continuous change of actual electricity consumption data.

[0072] For consecutive missing segments exceeding 15 minutes, the interpolation error will increase significantly due to the lack of sufficient local trend information. In this case, the missing segment should be marked as "data unavailable" and skipped during model training and real-time detection.

[0073] ② Outlier Correction

[0074] Sources of outliers: transient interference from sensors, jitter during ADC sampling, and communication errors.

[0075] An improved Hampel filter is used to identify and correct outliers. For each data point x... i Take k points before and after the first point (k is usually 5 to 10) to form a sliding window of length 2k+1, and calculate the median M of the data within the window. i And the median absolute deviation (MAD):

[0076] M i =median(x i-k x i-k+1 , ..., x i , ..., x i+k-1 x i+k );

[0077] MAD i =median(∣x i-k -M i |,|x i-k+1 -M i |,...,|x i+k -M i |);

[0078] in:

[0079] The median is calculated by taking the median value. Compared to the mean, the median is less sensitive to outliers—even if the window contains a very large outlier, the median is almost unaffected, making MAD a robust measure of deviation.

[0080] Anomaly detection rule: If the current point x i A value is considered an outlier if it meets the following conditions:

[0081] |x i -M i |>n0×MAD i ×1.4826;

[0082] in:

[0083] 1.4826 is a scaling factor. When the data follows a normal distribution, the expected value of MAD is equal to 1 / 1.4826 times the standard deviation, i.e., E[MAD] = σ / 1.4826. Therefore, multiplying by 1.4826 makes MAD equivalent to the standard deviation σ under the normality assumption, which can be compared to the "3σ principle". Specifically, for a standard normal distribution N(0,1), the theoretical value of MAD is equal to φ. ﹣1 (0.75) ≈ 0.6745 (because the absolute deviation of the median in the standard normal distribution is the distance of the 75th percentile from the median), then the scaling factor = 1 / 0.6745 ≈ 1.4826.

[0084] n0 is a sensitivity coefficient set by the user, typically ranging from 3 to 5. The larger n0 is, the higher the tolerance for outliers; only points with significant deviations will be marked. In this invention, n0 = 3.

[0085] Correction method: For the point x marked as an outlier i Use the median M within the window i Replace the original value x i .

[0086] ③ Data normalization

[0087] Because different electrical parameters have vastly different dimensions and value ranges (e.g., voltage 220V, current 100A, temperature 40℃), directly inputting them into the model would cause large-scale features to dominate the training process, while small-scale features would be ignored. Therefore, Z-score standardization is performed on all dimensions of the time series data.

[0088] x norm =(x-μ train ) / σ train ;

[0089] in:

[0090] μtrain The mean of historical training data for this dimension (such as the mean of historical phase A voltage) is calculated from the training set data;

[0091] σ train The standard deviation of the historical training data for this dimension (such as the standard deviation of the historical phase A voltage) is calculated from the training set data.

[0092] These two statistics are calculated only from the training set and are used consistently in subsequent validation, test, and real-time inference phases, without being recalculated. This is done to ensure consistent data distribution during training and inference, preventing shifts in the model's input distribution due to changes in data distribution.

[0093] S2: Construction of a benchmark model for typical electricity consumption behavior based on VMD-CNN-GRU

[0094] The goal of this step is to train a model that can accurately predict the expected normal values ​​of various electrical parameters at the next moment based on current contextual information (time, season, historical electricity consumption trends, etc.). The model building process consists of three stages: multi-scale time-frequency decomposition, feature matrix construction, and hybrid deep learning model training.

[0095] S2.1 Multi-scale time-frequency decomposition based on variational mode decomposition (VMD)

[0096] ① The electricity consumption data of the distribution area is a typical nonlinear and non-stationary time-series signal, which contains components of multiple different time scales:

[0097] High-frequency components (cycle to minute level): instantaneous load fluctuations, spike noise caused by appliance start-up and shutdown, random disturbances;

[0098] Mid-frequency components (period to hourly): daily load curve changes, user daily routines (such as morning, noon and evening peak electricity consumption);

[0099] Low-frequency components (period ~ day / week / seasonal level): differences between weekdays and rest days, seasonal changes in electricity consumption (cooling in summer, heating in winter).

[0100] If the original signal is modeled directly, the model needs to learn features at all scales simultaneously, which leads to complex model structure, difficult training, and poor generalization performance. VMD adaptively decomposes complex signals into several amplitude-modulated and frequency-modulated signal components with specific center frequencies, thereby separating features at different scales into different components, making it easier for subsequent models to learn the patterns at each scale.

[0101] ② Mathematical model of VMD algorithm

[0102] VMD decomposes the original signal f(t) into K discrete intrinsic mode function (IMF) components u.k (t), k=1,2,...,K. Each component u k (t) can be expressed as:

[0103] u k (t)=A k (t)cos(φ k (t));

[0104] in:

[0105] A k (t) is the instantaneous amplitude of the k-th component, which changes slowly over time and represents the energy of the signal at this scale;

[0106] φ k (t) represents the instantaneous phase of the k-th component, and its derivative is the instantaneous angular frequency ω. k (t)=φ k ′(t), and ω k (t)≥0 and changes slowly over time (i.e., the component is a narrowband signal).

[0107] The optimization objective of VMD is to minimize the sum of the bandwidths of each component, given the constraint that the sum of all components equals the original signal. By introducing Lagrange multipliers and a quadratic penalty term, it is transformed into an unconstrained variational problem, with its augmented Lagrangian function being:

[0108] ;

[0109] The meanings of each parameter are shown in Table 2 below:

[0110] Table 2: Parameter Description Table

[0111]

[0112] Decomposition of the physical meaning of the augmented Lagrange function:

[0113] First item : Used to measure the bandwidth of each component.

[0114] Where (δ(t)+j / (πt))∗u k (t) is u k The analytic signal (i.e., Hilbert transform) of (t) will transform the real signal u k (t) is expanded into a complex signal;

[0115] The instantaneous frequency of the analytic signal is ω. k (t) = d[arg(analytic signal)] / dt;

[0116] Multiply by e ﹣jωktThen the spectrum is shifted to baseband, and the bandwidth of the analyzed signal is the bandwidth of the component.

[0117] Take the derivative of the baseband signal and the L2 norm, then measure its rate of change. The smaller the rate of change, the narrower the component band.

[0118] α controls the strength of the bandwidth penalty.

[0119] Second item The L2 norm measures the decomposition error, which is the deviation of the sum of all components from the original signal. The smaller the L2 norm, the higher the decomposition accuracy. This term encourages the sum of all components to restore the original signal as much as possible.

[0120] Third item : Lagrange constraint term, ensuring the constraint condition ∑u k =f is strictly satisfied. λ is dynamically adjusted during the iteration process; when the decomposition error is large, λ is increased to strengthen the constraint.

[0121] ③ VMD solution method – Alternating Multiplier Direction Method (ADMM)

[0122] VMD uses the ADMM algorithm for iterative solution, and the specific steps are as follows:

[0123] Initialization: Setting , , Initialize the values ​​(usually set to a zero vector or random values); set the number of iterations n=0; set the maximum number of iterations N. max =500, convergence criterion ε=1e-6.

[0124] (Note: Indicate u k The Fourier transform of (t) is more conveniently updated in the frequency domain by VMD because the bandwidth constraint is expressed in a simple algebraic form in the frequency domain.

[0125] Iterative update (n increments from 1):

[0126] Step 1: Update each modal component (Updated in the frequency domain)

[0127] ;

[0128] in:

[0129] This is the Fourier transform of the original signal f(t);

[0130] It is the sum of the spectra of all components except the k-th component in the nth iteration;

[0131] For the Fourier transform of the Lagrange multiplier λ(t) in the nth iteration;

[0132] The center frequency ω of the k-th component is represented by k The square of the nearby frequency offset, this term serves as the kernel of the Wiener filter, making... Energy is concentrated in ω k nearby.

[0133] The model states that the updated value of the current component equals the original signal minus the remainder of all other components (i.e., the input of the Wiener filter), and then passes through a filter with a step size of ω. k The output is obtained by a low-pass filter centered at ω (denominator term). The denominator is 1 + 2α(ω - ω). k ) 2 In ω=ω k When the minimum value of 1 is taken (no signal attenuation), it is far away from ω. k The value is very large (the signal is attenuated), therefore the updated value is... Keep only ω k Nearby frequency components.

[0134] Step 2: Update the center frequency ω k

[0135] ;

[0136] The model: It is the spectral centroid of the k-th component. The molecule is the frequency ω versus the spectral energy. The weighted integral, with the total spectral energy in the denominator, yields the ratio, which represents the "centroid" frequency of the spectrum. This update rule shifts the center frequency to the frequency where the energy of that component is most concentrated, making subsequent Wiener filtering more accurate.

[0137] Step 3: Update the Lagrange multipliers

[0138] ;

[0139] Where τ is the update step size, typically a value between 0 and 1, used to control the speed of constraint recovery. The meaning of this expression: the reconstruction error of all current components (…). The larger the value of τ, the greater the adjustment to λ, thus strengthening the constraint satisfaction in the next iteration. In this invention, τ = 0.3 is chosen.

[0140] Convergence criterion: Calculate the normalized error after each iteration:

[0141] ;

[0142] If error < ε (ε = 1e-6), the iteration is considered converged and the calculation stops; otherwise, the iteration continues until the maximum number of iterations N is reached. max .

[0143] ④ Decomposition Results

[0144] Finally, we obtain K=8 IMF components u1, u2, ..., u K Taking phase A current as an example, the typical physical meanings of each component after decomposition are shown in Table 3 below:

[0145] Table 3: Explanation of Physical Meaning

[0146]

[0147] Each component retains the complete time-domain information of the original signal in a specific frequency band, that is, it retains the amplitude of the signal at each time point in that frequency band.

[0148] S2.2 Construction of Input Feature Matrix

[0149] After VMD decomposition, the original signal in each dimension is expanded into K IMF components. Simultaneously, to incorporate environmental perception mechanisms, contextual features are added. The final input feature matrix comprises the following components, as shown in Table 4:

[0150] Table 4: Meaning of Input Feature Matrix

[0151]

[0152] Construction of sliding window samples:

[0153] Given the current time t, take the historical window length H, that is, use the data of the past H time steps to predict the next time step (t+1 minutes).

[0154] In this invention, H=1440 is chosen, corresponding to the past 1440 minutes, or the past 24 hours. The basis for this choice is that 24 hours can cover a complete daily load cycle (including the three peak electricity consumption periods in the morning, noon, and evening, and the nighttime off-peak), allowing the model to learn the complete daily periodic pattern.

[0155] Therefore, the input feature matrix for each sample is H × 101, which is a two-dimensional matrix with 1440 rows and 101 columns. The rows represent time steps (the past 1440 minutes, from t-1439 to t), and the columns represent 101 feature dimensions.

[0156] The labels of the samples (the target output of model training) are the original values ​​of 12 electrical parameters at the next time step (t+1 minutes), which is a 12-dimensional vector.

[0157] S2.3 Hybrid Deep Learning Model Structure Design

[0158] The model consists of two cascaded parts: a CNN module and a GRU module. The CNN is responsible for extracting local short-term patterns (such as abrupt changes in load shape and harmonic waveform characteristics), while the GRU is responsible for capturing long-term time-series dependencies (such as daily cycle patterns and long-term trends). The two complement each other to achieve accurate modeling of complex electricity consumption behavior.

[0159] ① CNN module (Convolutional Neural Network)

[0160] Mathematical definition of convolution operation:

[0161] The input is a two-dimensional matrix X∈R of size H×101. H×d (d=101 is the feature dimension). The l-th convolutional kernel W l ∈R k×d The convolution output (where k is the length of the convolution kernel on the time axis, i.e., the receptive field size) is:

[0162] ;

[0163] in:

[0164] Z l (t): The output value of the l-th convolutional kernel at time position t, which is a scalar;

[0165] X(t+i,j): The value of the j-th feature dimension at time (t+i) in the input matrix;

[0166] W l (i, j): The weight parameters of the l-th convolutional kernel at position (i, j) are trainable parameters learned by the model through backpropagation during training.

[0167] k: The time window length of the convolution kernel. In this invention, k=12 (i.e., a 12-minute time window). This value enables the convolution kernel to perceive local change patterns within approximately 12 minutes.

[0168] d: Number of feature dimensions: 101;

[0169] b l : The bias term of the l-th convolutional kernel, a trainable parameter;

[0170] f(·): Nonlinear activation function. This invention uses ReLU (Modified Linear Unit), which is defined as f(x) = max(0, x), to introduce nonlinear expressive power, enabling the model to learn complex nonlinear mapping relationships.

[0171] The hierarchical structure of the CNN module is shown in Table 5 below:

[0172] Table 5: Hierarchical Structure of CNN Modules

[0173]

[0174] Note: When H=1440, after two 2x pooling operations, the sequence length is 1440 / 4=360.

[0175] ② GRU module (Gated Loop Unit)

[0176] GRU uses two gating mechanisms, the reset gate and the update gate, to control the flow of information between time steps, selectively "remembering" important long-term dependent information and "forgetting" unimportant short-term noise.

[0177] The forward propagation formula for a GRU cell at time t:

[0178] Given the input x at the current time step t ∈R din (d) in The dimensions of the Flatten layer output (H / 4×128) and the hidden state h from the previous time step. t-1 ∈R dh (d) h (To determine the number of hidden units), the GRU calculation process consists of the following four steps:

[0179] Step 1: Calculate the reset gate r t

[0180] Resetting the door determines the previous hidden state. t-1 How much information should be "reset" or "ignore":

[0181] r t =σ(W ir x t +b ir +W hr h t-1 +b hr );

[0182] in:

[0183] r t ∈[0,1] dh : Reset gate vector, each element takes a value between 0 and 1, independently determining the degree of reset for the corresponding hidden unit. A value of 0 indicates that the information of the hidden unit from the previous time step is completely ignored, and a value of 1 indicates that it is completely preserved;

[0184] σ(·): Sigmoid function, σ(z) = 1 / (1+e ﹣z This function compresses any real number into the (0, 1) interval, making it suitable as a gating signal.

[0185] W ir ∈R dh×din The weight matrix input to the reset gate determines the current input x. t The degree of influence on the reset door is obtained through model training.

[0186] W hr ∈R dh×dh The weight matrix from the hidden state to the reset gate determines the previous hidden state h. t-1 The degree of influence on the reset door is obtained through model training.

[0187] b ir ∈R dh and b hr ∈R dh : Bias vector, learned through model training.

[0188] Step 2: Calculate the candidate hidden states

[0189] The candidate hidden state represents the content that should be "newly generated" at the current moment, based on the current input and the reset historical information:

[0190] ;

[0191] in:

[0192] : Candidate hidden state vector, each element takes a value between -1 and 1, representing the "newly generated" content based on the current input and the reset historical information;

[0193] tanh(·): hyperbolic tangent activation function, tanh(z) = (e^(-tanh(·))) z -e ﹣z ) / (e z +e ﹣z The output is compressed to the interval (-1, 1). Compared with the Sigmoid function, the output of tanh can be both positive and negative, making it more suitable as a candidate value for state updates.

[0194] W in ∈R dh×din The weight matrix input to the candidate state determines the current input x. t The contribution to the candidate state is learned through model training;

[0195] W hn ∈R dh×dh The weight matrix from the hidden state to the candidate state determines the previous hidden state h. t-1 The contribution to the candidate state is learned through model training;

[0196] b in ∈R dh and b hn ∈R dh Bias vector, obtained through model training;

[0197] ⊙: Hadamard product (element-by-element multiplication), i.e. (r t ⊙(W hn h t-1 )) j =(r t ) j ·(W hn h t-1 ) j .

[0198] Key logic explanation: r t ⊙(W hn h t-1 This indicates historical information after the gate has been reset and filtered. When r t When an element approaches 0, the previous history information of that hidden unit is "reset" (ignored), and the candidate hidden state is then... The corresponding element is completely determined by the current input x. t Decision; when r t When an element approaches 1, the hidden unit retains all historical information. This mechanism allows GRU to flexibly choose whether to generate new states based on historical information.

[0199] Step 3: Calculate the update gate z t

[0200] The update gate determines how much old historical information should be retained and how much new candidate information should be introduced into the final hidden state:

[0201] z t =σ(W iz x t +b iz +W hz h t-1 +b hz );

[0202] in:

[0203] z t ∈[0,1] dh Update the gate vector, where each element independently determines the old state h in the corresponding hidden unit. t-1 and new candidate status The mixing ratio;

[0204] W iz ∈R dh×dinThe weight matrix input to the update gate is learned through model training.

[0205] W hz ∈R dh×dh The weight matrix from the hidden state to the update gate is learned through model training;

[0206] b iz ∈R dh and b hz ∈R dh : Bias vector, learned through model training.

[0207] Step 4: Calculate the final hidden state h at the current time. t

[0208] Weighted fusion between the old state and the candidate state is achieved through linear interpolation:

[0209]

[0210] in:

[0211] h t ∈(﹣1,1) dh The final hidden state at the current moment, which contains a compressed representation of all input information up to the current moment;

[0212] 1-z t This represents the proportion of the old state that has been "forgotten".

[0213] An intuitive understanding of gating mechanisms:

[0214] When z t When h approaches 0 t ≈h t-1 In other words, the network chooses to retain most of the old information, allowing long-distance temporal dependencies to be propagated across multiple time steps. This is crucial for learning diurnal periodic patterns—for example, the load pattern at 8 AM may depend on the load pattern from 24 hours ago (i.e., yesterday morning at 8 AM). The GRU can address this by using z-axis data... t The information was kept close to 0 for 1440 time steps to "remember" it.

[0215] When z t When it approaches 1, This means the network chooses to discard old information and fully utilize the new information generated by the current input. This is important for dealing with sudden load changes (such as sudden power outages or load drops), because the old historical information is no longer applicable.

[0216] This flexible gating mechanism enables GRU to effectively capture time-series dependencies ranging from short-term (minutes) to long-term (days or even weeks), making it highly suitable for modeling electricity consumption data with multiple time scales.

[0217] The hierarchical structure of the GRU module is shown in Table 6 below:

[0218] Table 6: Hierarchical Structure of GRU Modules

[0219]

[0220] Note: The sequence length is H / 4 = 360 (when H = 1440), d in =128 (Number of output channels of the Flatten layer).

[0221] ③ Output layer (fully connected layer)

[0222] Map the 64-dimensional hidden state vector output by the GRU module to the 12-dimensional prediction target: ;

[0223] in:

[0224] h final ∈R 64 The final hidden state of the GRU module, i.e. the output of the GRU2 layer, is the deep feature representation of the entire input sequence.

[0225] W out ∈R 12×64 The output weight matrix is ​​obtained through model training, with each row corresponding to a linear combination coefficient of an output dimension.

[0226] b out ∈R 12 Output bias vector, obtained through model training;

[0227] ∈R 12 : Predicted output vectors, corresponding to the 12 electrical parameter values ​​at the next time step, in the following order: [U a U b U c I a I b I c ,P,Q,S,cosφ,f,T].

[0228] S2.4 Model Training and Baseline Model Generation

[0229] ① Loss Function

[0230] The Huber loss function is adopted, which combines the advantages of mean squared error (MSE) and mean absolute error (MAE)—using squared loss when the error is small (providing smooth gradients and fast convergence), and using linear loss when the error is large (insensitive to large errors, robust, and avoiding outliers dominating the gradient): ;

[0231] in:

[0232] y represents the actual observed value (training label). These are model predictions, all of which are normalized values.

[0233] δ is the switching threshold, controlling the switching point between squared loss and linear loss. In this invention, δ = 1.0 (at the normalized data scale).

[0234] The selection criteria for the loss function are as follows: Under normal operating conditions, the deviation between the model's predicted value and the actual value is usually small (|y-ŷ|≤δ). In this case, squared loss is used to provide a second convergence speed, enabling the model to fit finely on normal samples. When extreme outliers occur (such as data labeling errors or rare extreme operating conditions), the deviation may be large (|y-ŷ|>δ). In this case, linear loss is switched to avoid the excessive influence of extreme outliers on the gradient and improve training stability.

[0235] ② Optimizer

[0236] The Adam optimizer (Adaptive Moment Estimation) is employed, characterized by its ability to adaptively adjust the learning rate of each parameter without requiring manual tuning. The parameters are shown in Table 7 below.

[0237] Table 7: Parameter Overview

[0238]

[0239] Adam's update rules (taking the t-th iteration as an example):

[0240] First, calculate the first moment estimate m of the gradient. t (Momentum term) and second moment estimate v t (Adaptive learning rate term):

[0241] m t =β1·m t-1 +(1-β1)·g t ;

[0242] v t =β2·v t-1 +(1-β2)·g t 2 ;

[0243] in:

[0244] g t =∇ θ L(θ t-1 ) is the gradient vector of the loss function with respect to the model parameters θ, calculated using the backpropagation algorithm;

[0245] m t-1 and v t-1 The moment estimate for the previous iteration step is given by initial values ​​m0=0 and v0=0.

[0246] g t 2 g represents the element-wise square. t ⊙g t .

[0247] Because m t and v t The initial value is biased towards 0 (because the initial value is set to 0), so bias correction is needed:

[0248] in and Let β1 and β2 be powers of t, respectively, where t is the current iteration number. As t increases, the denominator approaches 1, and the effect of bias correction gradually decreases.

[0249] The final parameters are updated as follows: ;

[0250] The learning rate for each parameter is adaptively adjusted to... In other words, parameters with drastic gradient changes (large second moments) receive smaller learning rates, while parameters with gentle gradient changes receive larger learning rates, thus achieving more efficient optimization.

[0251] ③ The training parameters are set as shown in Table 8 below:

[0252] Table 8: Training Parameter Settings Table

[0253]

[0254] ④ Dataset partitioning

[0255] The historical dataset is divided into three subsets, as shown in Table 9 below:

[0256] Table 9: List of Subsets

[0257]

[0258] ⑤ Baseline model generation

[0259] After training, the model weight parameters that achieve the best performance on the validation set (i.e., minimize the Huber loss on the validation set) are saved. This model is the benchmark model for typical electricity consumption behavior.

[0260] The core capability of the baseline model is to accurately predict the expected normal values ​​of 12 electrical parameters at the next moment, given a 101-dimensional feature matrix (containing VMD decomposition components, time coding, holiday identifiers, and meteorological features) for the past H=1440 time steps (24 hours).

[0261] The actual physical meaning of the baseline model: This model learns the "pattern" of electricity consumption behavior in a distribution area under normal operating conditions—when the load should be high, when the load should be low, the variation patterns affected by seasons and weather, and the differences between weekdays and holidays, etc. When the system is operating normally, the model's predicted values ​​should be close to the actual values; when the system malfunctions, the actual values ​​will deviate from the model's predicted values, and this deviation is the basis for anomaly detection.

[0262] S3: Real-time residual calculation and multidimensional deviation assessment

[0263] This step is executed in real time on the smart terminal in the distribution area, and the current operating status is continuously monitored using the deployed benchmark model.

[0264] S3.1 Real-time Data Acquisition and Feature Engineering

[0265] The smart terminal in the distribution area performs the following operations every minute (i.e., at each sampling time t):

[0266] Step 1: Collect real-time data at the current time t

[0267] Collect the raw values ​​of the 12-dimensional electrical parameters at the current time t:

[0268] Y real (t)=[Ua real Ub real Uc real Ia real Ib real ,Ic real P real Q real S real cosφ real f real T real ];

[0269] Step 2: Maintain the historical data sliding window cache

[0270] The terminal maintains a circular buffer of length H=1440 in memory to store all data from the past 1440 minutes. Every minute, newly collected data is written to the buffer, while automatically overwriting the oldest data from 1440 minutes ago.

[0271] Step 3: Extract the feature matrix

[0272] Extract 12-dimensional raw electrical parameter history data from the buffer over the past H=1440 minutes and perform the following feature engineering operations:

[0273] VMD decomposition: Perform VMD decomposition on 1440 historical data points in each dimension (parameters are the same as in S2.1: K=8, α=2000), resulting in 12×8=96 IMF component sequences;

[0274] Time encoding: Calculate the sine and cosine codes of the hour based on the timestamp of the current time t;

[0275] Holiday indicators: Query the calendar based on the current date to obtain holiday indicators;

[0276] Meteorological characteristics: Retrieve temperature and weather code from the most recent cached meteorological API response.

[0277] Organize all the above features into a feature matrix X of size H×101. input (t), whose structure is exactly the same as the input feature matrix during training.

[0278] S3.2 Expected Normal Value Prediction

[0279] The characteristic matrix X input (t) Input the baseline model deployed on the terminal for forward inference computation. The forward propagation process of the model includes:

[0280] Forward propagation of the CNN module: Perform convolution, batch normalization, ReLU activation, pooling, and flattening operations layer by layer to convert the H×101 input into a one-dimensional feature vector sequence of length H / 4×128=360×128;

[0281] Forward propagation of the GRU module: Performs GRU gating computation (reset gate, update gate, candidate state, and final state computation) step by step over time, compressing the feature vector sequence into a 64-dimensional final hidden state vector;

[0282] Forward propagation of the output layer: Perform matrix multiplication of the fully connected layer to map the 64-dimensional hidden state to a 12-dimensional output vector.

[0283] Finally, the model outputs a 12-dimensional vector of expected normal values ​​at the current time t:

[0284]

[0285] Where Model(·) represents the complete forward propagation process of the CNN-GRU hybrid model.

[0286] S3.3 Residual Sequence Calculation

[0287] The deviation between the actual observed value and the model prediction value of each electrical parameter is calculated to obtain a 12-dimensional residual vector:

[0288]

[0289] Expanding element by element:

[0290]

[0291] The physical meaning of residuals:

[0292] When the system is running in normal steady state, the benchmark model can accurately predict the expected values ​​of electrical parameters. The residual ΔY(t) should be approximately zero-mean white noise, that is, randomly fluctuating around the value of 0 with a small amplitude (usually within ±0.1 to 0.2 standard deviations).

[0293] When abnormal signs appear (such as sudden load increase, voltage drop, harmonic distortion, equipment overheating, etc.), the actual value will deviate significantly from the "expected normal value" predicted by the model based on historical patterns. The residual ΔY(t) will show significant positive or negative deviation, and its absolute value may reach several times or even tens of times the normal standard deviation.

[0294] The sign of the residual also has important physical indicative meaning. For example:

[0295] ΔI a (t)>0: The actual current is greater than the expected value, which may mean that the load is out of normal range or a short circuit fault has occurred.

[0296] ΔU a (t)<0: The actual voltage is lower than the expected value, which may mean a voltage sag, increased line impedance or a ground fault.

[0297] ΔT(t)>0: The actual temperature is higher than the expected value, which may mean that the equipment is overloaded or poor contact is causing the heat to rise.

[0298] S3.4 Residual Normalization

[0299] To facilitate subsequent cross-dimensional unified threshold processing, the residuals of each dimension are normalized and transformed into standard residuals (dimensionless):

[0300] ΔY norm，i (t)=ΔY i (t) / σres，i ;

[0301] in:

[0302] ΔY i (t) represents the residual of the i-th electrical parameter at time t (i.e., the actual value minus the predicted value);

[0303] σ res，i Let $\begin{pmatrix}$ be the standard deviation of the residuals in the $i$-th dimension of the training set. This standard deviation is pre-calculated and fixed during the model training phase. Specifically, the model inference is run once on the training set to obtain the residual sequence of all training samples, and then the standard deviation of this sequence is calculated.

[0304] The meaning of normalization: ΔY norm，i (t) represents the multiple of the current residual relative to the normal fluctuation range, and is dimensionless. For example, ΔI norm，a The value of (t) = 3.5 indicates that the actual value of the phase A current deviates from the model prediction by 3.5 times the normal standard deviation of the fluctuation. This provides a unified measure of the severity of anomalies among different electrical parameters.

[0305] S3.5 Comprehensive Deviation Index

[0306] To provide a simple and intuitive overall deviation measure, a comprehensive deviation index D(t) is defined:

[0307]

[0308] in:

[0309] w i The weighting coefficient for the i-th electrical parameter reflects the importance of this parameter in fault detection and satisfies... ;

[0310] Weight determination methods: The weights can be automatically determined based on historical data using the entropy weight method, or manually set based on expert experience.

[0311] The steps for determining weights using the entropy weight method are as follows:

[0312] For each dimension i, calculate the information entropy E of its normalized residual on the training set. i : , where p ij Entropy represents the percentage of the absolute residual value of the j-th sample in the i-th dimension among all samples (normalization is required), where N is the total number of samples. A smaller entropy value indicates a more regular variation in the residuals of that dimension (i.e., a stronger indication of anomalies), and should therefore be assigned a larger weight.

[0313] Calculate the weights for each dimension:

[0314] The exemplary weight allocation used in this invention is shown in Table 10 below:

[0315] Table 10: Exemplary Weight Allocation Table

[0316]

[0317] The larger the value of the overall deviation D(t), the greater the deviation of the current system state from the normal baseline, and the higher the probability of an anomaly. Although this indicator is not directly used for anomaly detection (because it is a composite value, losing the independent information of each parameter), it can serve as an auxiliary reference indicator for subsequent graded early warning.

[0318] S4: Anomaly Detection Based on POT Adaptive Threshold and Multi-parameter Cross-validation

[0319] Traditional methods use fixed thresholds (such as setting an alarm if the current fluctuation exceeds ±10%), but this method cannot adapt to the changes in the normal fluctuation range caused by dynamic load changes. When the load is stable at night, a 10% fluctuation may indeed be abnormal; but when the load is at its peak in the morning, a 10% fluctuation may just be a normal phenomenon.

[0320] This invention proposes an adaptive dynamic threshold method based on Extreme Value Theory (EVT), which can automatically adjust the detection threshold according to the recent data distribution characteristics. At the same time, it combines multi-parameter cross-validation to filter out noise interference, thereby achieving high-precision and low-false-alarm anomaly detection.

[0321] S4.1 Sliding Window and Residual Statistics

[0322] Normalized residual sequence ΔY for each dimension norm，i (t) Set a sliding time window with a length of L time steps (L=1440, i.e., 24 hours). The window contains the normalized residual values ​​of the current time and the past L-1 time steps:

[0323] W i (t)={ΔY norm，i (t-L+1), ΔY norm，i (t-L+2), ..., ΔY norm，i (t)};

[0324] The selection of window length L=1440 is based on the fact that a 24-hour window can cover a complete daily load cycle. This means that the window will contain residual values ​​under various typical operating conditions (low point, flat period, peak period), providing sufficient and representative statistical samples for extreme value modeling. If L is too small (e.g., 1 hour), the data within the window may all come from the same operating condition (e.g., all low point periods), and its statistical characteristics cannot represent the normal fluctuation range of the whole day.

[0325] S4.2 Dynamic Threshold Generation Based on POT (Peak Exceeds Threshold) Method

[0326] The POT method states that for any distribution that satisfies certain regularization conditions (not limited to the normal distribution), the distribution of extreme values ​​exceeding a sufficiently high threshold asymptotically follows a generalized Pareto distribution (GPD). This property allows us to reliably infer the probability and magnitude of extreme events by modeling only the extreme values, even without fully understanding the complete picture of the data (i.e., without knowing what distribution the data follows).

[0327] The application steps of the POT method in this invention are as follows:

[0328] Step 1: Extreme Value Sequence Extraction

[0329] For the sliding window W i The absolute values ​​of L=1440 normalized residuals within (t) {|ΔY norm，i (t-L+1)∣,...,∣ΔY norm，i (t)∣}, set an initial threshold u0, extract all peak values ​​exceeding u0 to form an extreme value sequence.

[0330] The initial threshold u0 is set as follows: Take the 90th percentile (i.e., P) of the absolute value sequence of all residuals within the window. 90 90% of the data points are less than this value. This means that only the largest 10% of the residual values ​​are taken as extreme value candidates. The basis for this selection is that the 90th percentile is usually higher than the main range of normal fluctuations, which can screen out the truly noteworthy extreme values, while retaining a sufficient number of samples for GPD fitting (1440 points × 10% = 144 extreme points, which is statistically sufficient).

[0331] The extracted extreme value sequence is:

[0332] {z1, z2, ..., z N}={∣ΔY norm，i (t j )∣-u0>0,j=1,2,...,N};

[0333] in:

[0334] t j This is the moment when the absolute value of the residual exceeds u0;

[0335] z j The excess amount, that is, the magnitude by which the absolute value of the residual exceeds the initial threshold u0, is always a positive number;

[0336] N is the number of extreme points within the window that exceed u0.

[0337] Step 2: Fitting the Generalized Pareto Distribution (GPD)

[0338] For the extreme value sequence {z1, z2, ..., z} N The parameters of the GPD distribution are fitted using the maximum likelihood estimation (MLE) method.

[0339] The cumulative distribution function (CDF) of GPD:

[0340]

[0341] in:

[0342] z is the excess amount exceeding the initial threshold u0, z = |ΔY norm |-u0>0;

[0343] ξ is a shape parameter, which takes the value of a real number and determines the thickness of the tail of the distribution.

[0344] ξ>0: Heavy tail distribution, with a high probability of extreme values ​​occurring, and extreme values ​​can be very large (such as residuals when the load fluctuates drastically).

[0345] ξ=0: Exponential tail distribution (i.e., Gumbel type), extreme values ​​decay exponentially;

[0346] ξ<0: Bounded distribution (Weibull type), extreme values ​​have an upper bound.

[0347] β>0 is the scale parameter, which determines the degree of dispersion of the distribution. The larger β is, the wider the dispersion of extreme values ​​and the greater the range of fluctuation.

[0348] The probability density function (PDF) of GPD (used for subsequent MLE calculations):

[0349]

[0350] Maximum likelihood estimation (MLE):

[0351] For the observed extreme value samples {z1, z2, ..., z...} N Assuming they are independent and identically distributed, their joint probability density (i.e., likelihood function) is:

[0352] L(ξ,β)=∏ j=1 N g(z j ;ξ,β);

[0353] For ease of calculation, the logarithm is usually taken to obtain the log-likelihood function:

[0354] ℓ(ξ,β)=∑ j=1 N logg(z j ;ξ,β);

[0355] Substitute the PDF expression of GPD (taking ξ≠0 as an example):

[0356] ℓ(ξ,β)=﹣Nlogβ﹣(1+1 / ξ)∑ j=1 N log(1+ξz j / β);

[0357] The goal of MLE is to find parameter values ​​that maximize ℓ(ξ, β). This is an unconstrained optimization problem that can be solved using numerical optimization algorithms (such as the L-BFGS-B quasi-Newton method). This algorithm iteratively calculates the gradient (or approximate gradient) to gradually approximate the optimal parameters.

[0358] Step 3: Dynamic threshold back-calculation

[0359] Based on the well-fitted GPD model, the following key question can be answered: "Given an acceptable false alarm rate p..." FA Given this premise, what detection threshold should be set?

[0360] False alarm rate p FA Definition: The probability that a system will erroneously issue an alarm during normal operation. This invention sets p... FA =1 / 1440, meaning that a maximum of one false alarm is allowed per day (1440 minutes).

[0361] Mathematical derivation of threshold inversion:

[0362] First, the empirical probability that any point within the window exceeds the initial threshold u0 is N / L (N is the number of points exceeding u0, and L=1440 is the total number of points).

[0363] We hope to find a new threshold T dynamic >u0, such that any point exceeds T dynamic The probability is equal to the set false alarm rate p. FA :

[0364] P(∣ΔY norm |>T dynamic )=p FA ;

[0365] According to the conditional probability formula:

[0366] P(∣ΔY norm |>T dynamic )=P(∣ΔYnorm |>u0)·P(|ΔY norm |>T dynamic ||ΔY norm |>u0);

[0367] in:

[0368] P(∣ΔY norm |>u0)=N / L, which is the empirical probability of exceeding the initial threshold;

[0369] P(∣ΔY norm |>T dynamic ||ΔY norm |>u0)=P(z>T dynamic -u0)=1−G ξ，β (T dynamic -u0), which is the survival function of GPD (the probability of exceeding a certain value).

[0370] therefore:

[0371] p FA =N / L·[1-G ξ，β (T dynamic -u0)];

[0372] Summarized as follows:

[0373] 1-G ξ，β (T dynamic -u0)=p FA ·L / N;

[0374] Substituting into the survival function of GPD (when ξ≠0):

[0375] (1+ξ(T dynamic -u0) / β) ﹣1 / ξ =p FA ·L / N;

[0376] Take -1 / ξ to the power of both sides:

[0377] 1+ξ(T dynamic -u0) / β=(N / p FA ·L) ξ ;

[0378] Solve for T dynamic :

[0379] When ξ≠0:

[0380] T dynamic =u0+β / ξ·[(N / p FA ·L) ξ -1];

[0381] When ξ=0 (simplified form obtained through limit derivation):

[0382] T dynamic =u0+β·ln(N / p FA ·L);

[0383] The meanings and sources of each parameter are summarized in Table 11 below:

[0384] Table 11: Summary of Meanings

[0385]

[0386] T dynamic The optimal detection threshold is adaptively calculated based on the actual data characteristics of the most recent 24 hours, under a given false alarm rate constraint.

[0387] When the load fluctuates significantly (e.g., during morning and evening peak hours, when current fluctuations are severe): the tail of the extreme value sequence is thicker (ξ is larger) and the dispersion is greater (β is larger), resulting in a higher calculated T. dynamic It will automatically increase, thus avoiding false alarms caused by normal fluctuations;

[0388] When the load is stable (e.g., late at night, when current fluctuations are small): the tail of the extreme value sequence is thinner (ξ is smaller), T dynamic It will automatically decrease, thus enabling it to sensitively detect even the slightest abnormal deviations.

[0389] This adaptive characteristic enables the present invention to maintain optimal detection performance during periods of drastic load fluctuations and periods of stable load, which is something that traditional fixed threshold methods cannot achieve.

[0390] Step 4: Adaptive Update Mechanism

[0391] The sliding window moves forward over time. At each time step t, the window is updated with the latest residual value ΔY. norm，i (t), while discarding the oldest residual value ΔY. norm，i (tL).

[0392] To reduce computational burden, it is unnecessary to re-execute the complete POT calculation at every time step. This invention sets up a complete POT recalculation every 15 time steps (i.e., every 15 minutes). Between two recalculations, the most recently calculated threshold T is used. dynamic Anomaly detection is performed. The selection of this update frequency is based on the fact that the statistical characteristics of the load usually do not change drastically within 15 minutes. Therefore, a 15-minute update interval is sufficient to ensure the timeliness of the threshold, while significantly reducing computational overhead.

[0393] S4.3 Initial Inspection Anomaly Marking

[0394] For the current time t, for each dimension i, compare the absolute value of its normalized residual with the magnitude of the dynamic threshold to generate an initial anomaly marker:

[0395]

[0396] Where T dynamic，i (t) is the dynamic threshold of the i-th dimension at time t (if the current time is between two POT recalculations, the threshold calculated most recently is taken).

[0397] Flag initial，i (t)=1 indicates that the deviation of the i-th electrical parameter at time t exceeds the dynamic threshold and is marked as an initial detection anomaly.

[0398] S4.4 Multi-parameter cross-validation

[0399] Initial anomalies in single-parameter tests may be caused by the following non-fault factors, which need to be filtered out through cross-validation:

[0400] Natural load fluctuations: For example, the normal start-up of a high-power motor will cause a momentary increase in current, but at the same time the voltage will drop accordingly (because the starting current is large and the line voltage drop is large), and the power factor will also change. There is a definite physical relationship among the three: P=UIcosφ. If these three factors change synchronously and conform to physical laws, it is a normal fluctuation rather than a fault.

[0401] Transient interference in sensors: Electromagnetic interference or ADC sampling jitter causes a jump in the reading of a single sensor, while other related parameters do not change synchronously;

[0402] Communication errors: Single point data errors, usually manifested as isolated large jumps, lacking physical consistency.

[0403] This invention designs a cross-validation logic based on physical association rules. Its core principle is that in an electrical system, a real physical fault (such as a short circuit, grounding, or overload) will simultaneously produce observable abnormal responses on multiple physically associated parameters. If only one parameter is abnormal, it is likely noise or interference; if multiple parameters are abnormal simultaneously, and the abnormal pattern conforms to physical laws, then it is confirmed as a real fault.

[0404] The physical association rules are shown in Table 12 below:

[0405] Table 12: Physical Association Rules Table

[0406]

[0407] Cross-validation decision rules:

[0408] Let the set of all parameters marked as initial anomalies at time t be:

[0409] A initial (t)={i∣Flag initial，i (t)=1,i=1,2,...,12};

[0410] Rule 1: Single-parameter veto principle

[0411] If |A initial If (t)∣<2 (meaning only one parameter initially detects an anomaly), it is directly judged as noise interference or sensor false alarm, and no valid anomaly event is triggered:

[0412] Flag valid (t)=0, when |A initial (t)∣<2;

[0413] Rule 2: Multi-parameter correlation verification

[0414] If |A initial If (t)∣≥2 (at least two parameters simultaneously detect anomalies), then calculate the association confidence C between each pair of parameters that detect anomalies. ij :

[0415]

[0416] in:

[0417] |ΔY norm，i (t)∣ / T dynamic，i (t): This is called the deviation ratio of the i-th parameter, reflecting the multiple of the actual deviation of the parameter relative to the dynamic threshold. A ratio > 1 indicates that the threshold is exceeded, and the larger the ratio, the more significant the anomaly.

[0418] C ij (t): The joint outlier confidence of parameters i and j is the product of the two deviation ratios and the physical correlation coefficient;

[0419] s ij For physical association coefficients: According to the physical association rule table above, if parameters i and j are marked as associated in the table, then s ij =1, otherwise s ij =0.

[0420] Final valid anomaly determination rules:

[0421]

[0422] Where C th To determine the confidence threshold, this invention uses C. th =1.0.

[0423] C th=1.0 means and originates from the fact that when the product of the deviation ratios of two parameters is greater than 1, it indicates that at least one parameter's deviation ratio is greater than 1 (i.e., exceeds the threshold), and the other parameter also has at least some deviation. And s ij =1 guarantees that the two parameters are physically correlated. A product greater than 1 means that the two correlated parameters show significant anomalies at the same time, which is almost impossible to be caused by random noise (because the probability of two independent noises deviating significantly at the same time is extremely low). Therefore, it can be confirmed as a real anomaly event with high confidence.

[0424] Rule 3: Handling Special Circumstances

[0425] For certain special but important fault types, such as a sudden increase in zero-sequence current (indicating a single-phase ground fault, seriously endangering personal safety), even if only the zero-sequence current parameter triggers the initial detection anomaly (because the voltage change may still be very small), it is directly confirmed as a valid abnormal event. This is the "fuse tripping" rule set up to ensure safety.

[0426] Flag valid (t)=1, if i=I0 and |ΔY norm，I0 (t)∣>2·T dynamic，I0 (t);

[0427] That is, when the normalized residual of the zero-sequence current exceeds twice the dynamic threshold, an alarm is triggered directly without cross-validation.

[0428] S4.5 Valid Abnormal Event Confirmation and Feature Packaging

[0429] When Flag valid When (t)=1, a valid abnormal event is confirmed to have occurred at time t. The terminal immediately packages the complete characteristic information of the event to form an abnormal characteristic data packet, the structure of which is shown in Table 13 below:

[0430] Table 13: Structural Overview

[0431]

[0432] S4.6 Continuous Anomaly Tracking and Confirmation Mechanism

[0433] To prevent false alarms caused by single noise interference, and to ensure that genuine persistent anomalies are not missed, this invention designs an anomaly persistence tracking mechanism:

[0434] Initial confirmation: When Flag valid When (t)=1 is first true, record that moment as t. onset (At the moment of the abnormality's onset), but instead of immediately uploading an alarm to the cloud, it enters a confirmation waiting period.

[0435] Confirmation Waiting Period: Continuously monitor the Flag over the next M consecutive time steps (M=5, i.e., 5 minutes). valid Status:

[0436] If there are at least 3 time steps with a Flag in M ​​time steps. valid =1 (majority confirmation principle), then the abnormal event is confirmed as a continuously valid abnormal event, and the S5 cloud upload and alarm process is executed immediately;

[0437] If less than 3 time steps are confirmed, it is determined to be a transient disturbance event, and no alarm is triggered; only the log is recorded for future reference.

[0438] Continuous monitoring: For confirmed persistent anomalies, continue monitoring until a flag is set. valid If (t) is 0 for 10 consecutive minutes, the abnormal event is considered to have ended, and the event end time t is recorded. end And total duration.

[0439] This mechanism effectively avoids false alarms caused by momentary disturbances (such as the instantaneous impact of starting a large motor, lightning strikes, etc.), while ensuring that real, persistent faults (such as overload, three-phase imbalance) are not missed.

[0440] S5: Hierarchical early warning and fault diagnosis based on cloud-edge collaboration

[0441] S5.1 Cloud-Edge Collaborative Architecture Design

[0442] This invention adopts a two-level collaborative architecture of terminal edge computing layer and cloud big data analysis layer:

[0443] ① Terminal edge computing layer (deployed on smart terminals in the distribution area)

[0444] Responsible for performing the following real-time computing tasks:

[0445] Real-time data acquisition and preprocessing (S1);

[0446] Real-time VMD decomposition and feature matrix construction (S2);

[0447] Real-time forward inference and residual calculation of the baseline model (S3);

[0448] Adaptive dynamic threshold generation and initial anomaly detection (S4.1~S4.3);

[0449] Multi-parameter cross-validation and valid anomaly confirmation (S4.4~S4.6);

[0450] The hardware platform of the edge computing layer adopts ARM Cortex-M7 or Cortex-A series embedded processors (such as STM32H7 series or NXP i.MX RT series), with a main frequency of 400MHz to 1GHz, equipped with 512KB to 2MB SRAM and 16MB to 64MB Flash, which is sufficient for running lightweight model inference and real-time computing tasks.

[0451] The edge layer only uploads abnormal feature data packets to the cloud after confirming a continuously valid abnormal event. During normal operation, the terminal operates autonomously without data exchange with the cloud. This design greatly reduces the dependence on 4G / 5G network bandwidth and stability, and effectively controls communication costs.

[0452] ② Cloud big data analysis layer (deployed on cloud server)

[0453] It is responsible for performing the following compute-intensive tasks:

[0454] Receiving abnormal feature data packets uploaded by multiple distribution station area terminals;

[0455] Performing deep pattern matching and fault type diagnosis based on the fault feature knowledge base (S5.3);

[0456] Evaluating fault severity and generating hierarchical early warnings (S5.4);

[0457] Storing and statistically analyzing historical abnormal data;

[0458] Periodically retraining, updating and optimizing the model;

[0459] Performing global linkage analysis and early warning across distribution station areas;

[0460] The cloud platform adopts a microservice architecture and is deployed on mainstream cloud platforms (such as Alibaba Cloud, Huawei Cloud, AWS). The bottom layer uses a time-series database (such as InfluxDB) to store massive time-series data, and a relational database (such as PostgreSQL) to store knowledge bases and business data.

[0461] S5.2 Construction of Fault Feature Knowledge Base

[0462] The knowledge base is the "reference dictionary" for fault diagnosis in the cloud, which pre-stores "fingerprint" features of various typical faults.

[0463] ① Data source

[0464] The data sources of the knowledge base include two categories:

[0465] Historical real-world fault records (60%): These are real-world fault cases recorded in the past, extracted from the power grid operation and maintenance management system. They include time-series data of various electrical parameters before and after the fault, fault type labels (manually confirmed), fault causes, and corrective measures. These are the most reliable sources of knowledge.

[0466] Simulated fault data (40%): Using power system simulation software (such as PSCAD / EMTDC, MATLAB / Simulink's Simscape Electrical toolbox), detailed models of the low-voltage distribution area electrical system are built to simulate various fault conditions (short circuit, grounding, overload, harmonics, etc.) and generate simulated fault data. Simulated data can cover fault types and conditions that may be missing from real data, supplementing the completeness of the knowledge base.

[0467] ② Fault type library

[0468] The knowledge base covers the following eight typical fault types, as shown in Table 14 below:

[0469] Table 14: List of Typical Faults

[0470]

[0471] ③ Fault Feature Extraction and Coding

[0472] For each failure case (including real-world and simulation cases), the following feature vectors are extracted:

[0473] Feature group 1: Sign pattern vectors of residuals for each parameter (7 dimensions)

[0474] At the fault occurrence time t0, extract the residual signs of the three-phase voltage, three-phase current, and zero-sequence current:

[0475] S k =[sign(ΔU a ), sign(ΔU b ), sign(ΔU c ), sign(ΔI a ), sign(ΔI b ), sign(ΔI c ), sign(ΔI0)]

[0476] Where sign(x) is the sign function, which takes three values:

[0477] ;

[0478] ε is a small threshold (set to 0.1) used to ignore small random fluctuations. For example, the typical symbol pattern for a single-phase ground fault (assuming phase A is grounded) of FT-01 is:

[0479] S FT﹣01 =[-1, 0, 0, +1, -1, -1, +1];

[0480] Physical meaning: Phase A voltage ↓ (grounding causes voltage drop), Phase A current ↑ (fault current flows into the grounding point), Phase B and C currents slightly ↓ (voltage drop causes non-fault phase load current to slightly decrease), zero-sequence current ↑ (zero-sequence component appears).

[0481] Feature group 2: Relative magnitude vector of residuals for each parameter (12-dimensional)

[0482] Calculate the deviation ratio of each parameter at the time of the fault occurrence:

[0483] R i =∣ΔY norm，i (t0)∣ / T dynamic，i (t0), i=1, 2, ..., 12;

[0484] This feature reflects the relative severity of anomalies in each parameter. R for different fault types i Different distribution patterns, for example:

[0485] Overload fault FT-03: The R value of the three-phase current is 35 times, while the R value of the voltage is only 11.5 times.

[0486] Single-phase ground fault FT-01 (A-phase grounding): A-phase current R value is extremely large (up to 10 times or more), A-phase voltage R value is relatively large, and zero-sequence current R value is extremely large.

[0487] Feature group 3: Time-series waveform features (5 dimensions)

[0488] The following features were extracted from the residual time series data 5 minutes before and after the fault, as shown in Table 15 below:

[0489] Table 15: Feature Extraction Overview

[0490]

[0491] Feature group 4: Post-fault steady-state characteristics (12 dimensions)

[0492] Record the changes in the values ​​of each electrical parameter relative to their pre-fault state under the new steady-state condition after the fault (usually 30 seconds after the fault occurs, when all parameters tend to stabilize):

[0493] ΔY steady，i =Y after，i -Y before，i ;

[0494] For example, after a single-phase ground fault, the steady-state voltage of the faulted phase drops to about 0.5 to 0.7 pu, while the voltage of the non-faulted phase rises to about 1.0 to 1.1 pu.

[0495] ④ Storage structure of the knowledge base

[0496] Each fault type FT﹣j corresponds to a fault fingerprint template, stored in the following structure:

[0497]

[0498] in:

[0499] ID j Fault type number, such as "FT-01";

[0500] Desc j A textual description of the fault, such as "A-phase metallic single-phase ground fault";

[0501] The mean vector calculated from the feature vectors of all j-th type of fault cases has a dimension of M=7+12+5+12=36, representing the typical characteristics of this type of fault.

[0502] : The standard deviation vector calculated from the feature vectors of all j-th type of fault cases, with a dimension of 36, representing the dispersion of the fault features of this type;

[0503] Sim th，j The similarity threshold for this type of fault is set; cases below this threshold are considered mismatches. The calculation method is as follows: on the training set, calculate the inner product similarity between each case of this type of fault and the template, and take the 5th percentile value as the threshold to ensure that at least 95% of cases of this type of fault can be correctly identified.

[0504] Count j The number of samples for this type of failure reflects the statistical reliability of the template.

[0505] S5.3 Cloud-based Fault Diagnosis Process

[0506] Once the cloud platform receives the abnormal data packet uploaded by the terminal, it executes the following diagnostic process:

[0507] Step 1: Feature Vectorization

[0508] The received abnormal feature data packets are converted into 36-dimensional feature vectors identical to those in the knowledge base. ;

[0509] in:

[0510] S testA 7-dimensional symbol pattern vector (extracted from the residual symbols of the outlier parameters);

[0511] R test It is a 12-dimensional relative magnitude vector (extracted from the deviation ratio of outlier parameters);

[0512] W test A 5-dimensional time-series waveform feature vector (extracted from the waveform_segment field);

[0513] ΔY steady，test It is a 12-dimensional steady-state change vector (extracted from the changes in each parameter before and after the anomaly occurred).

[0514] Step 2: Mahalanobis distance similarity calculation

[0515] calculate The Mahalanobis distance between the model and the templates for each fault type in the knowledge base is then converted into similarity. ;

[0516] in:

[0517] M=36 represents the total number of feature dimensions;

[0518] F test，m : The m-th feature value of the abnormal event to be diagnosed;

[0519] F mean，j，m The mean value of the j-th type of fault on the m-th feature in the knowledge base;

[0520] F std，j，m The standard deviation of the j-th type of fault in the knowledge base on the m-th feature;

[0521] ε=1e-8: a very small constant to prevent division by zero.

[0522] The characteristics of Mahalanobis distance are: it takes into account the variance differences of each feature—features with large variance contribute less to the distance, while features with small variance contribute more, which is equivalent to normalizing each feature according to its variance. This allows features with large differences in numerical range (such as sign features with values ​​of {-1, 0, 1} and amplitude features with values ​​of 0 to 10) to be reasonably weighted.

[0523] Convert Mahalanobis distance to similarity (values ​​from 0 to 1, with higher values ​​indicating greater similarity):

[0524] Sim j =1 / (1+D Mahal，j );

[0525] Step 3: Fault Type Determination

[0526] Find all that satisfy Sim j ≥Simth，j The fault type is selected, and the type with the highest similarity is taken as the diagnostic result:

[0527] j diagnose =argmax j Sim j ;

[0528] Simultaneously output the diagnostic confidence level:

[0529] Conf=Sim jdiagnose ;

[0530] If all Sim j All are less than the corresponding Sim th，j If the fault is identified as an unknown fault type, it will be marked as "awaiting manual assessment".

[0531] Step 4: Severity Assessment

[0532] The severity level is determined based on a comprehensive assessment of the following four factors:

[0533] Severity = α1·min(1, maxR) i / R critical )+α2·min(1, duration / D critical )+α3·D equip +α4·D type ;

[0534] in:

[0535] maxR i The maximum deviation ratio among all parameters (i.e., |ΔY) norm | / T dynamic The maximum value reflects the severity of the anomaly;

[0536] R critical =5.0: Critical deviation multiple; exceeding 5 times is considered extremely serious.

[0537] duration: number of minutes the anomaly lasted;

[0538] D critical =60 (minutes): Critical duration; exceeding 60 minutes is considered extremely serious.

[0539] D equip ∈{0.3, 0.6, 1.0}: Equipment importance coefficient, 1.0 for transformer outgoing line side, 0.6 for main branch, and 0.3 for terminal branch;

[0540] D type∈{0.3, 0.6, 1.0}: Fault type hazard factor, grounding / short circuit is 1.0, overload / unbalance / phase loss is 0.6, and general concern is 0.3;

[0541] α1=0.4, α2=0.2, α3=0.2, α4=0.2: weighting coefficients, summing to 1.

[0542] The value of Severity is in the range of [0, 1].

[0543] S5.4 Three-level Early Warning Mechanism

[0544] Based on the severity rating, the alerts are divided into three levels, as shown in Table 16 below:

[0545] Table 16: Overview of Early Warning Levels

[0546]

[0547] S5.5 Linkage Control and Closed-Loop Feedback

[0548] For certain faults in the Level 3 emergency alarm, the system will automatically perform the following protective actions upon triggering the alarm:

[0549] ① Fault isolation

[0550] The cloud sends a trip command to the protection control module of the distribution transformer terminal, triggering the tripping of the smart circuit breaker in the corresponding branch. The tripping logic is divided into two levels:

[0551] Selective tripping: If there is a branch-level circuit breaker, only the faulty branch is disconnected, while the main branch and other normal branches are kept powered (the faulty branch is locked by analyzing the current residual of each branch).

[0552] Main circuit trip: If the faulty branch cannot be identified (such as a fault on the transformer outgoing side), disconnect the main circuit breaker.

[0553] ② Backup power switching

[0554] If the transformer substation is equipped with a backup power line (such as dual power supply) or a tie switch (interconnected with adjacent substations), it will automatically switch to the backup power supply to ensure continuous power supply to non-faulty areas. The switching time is controlled within 500ms (achieved through an automatic transfer switch (ATS)) to ensure that the impact on sensitive loads is minimized.

[0555] ③ Fault elimination confirmation

[0556] After tripping, continuously monitor the residual sequence. If the residual recovers to below the threshold within 10 seconds after tripping, the fault is confirmed to have been successfully isolated, and "protection action successful" is recorded; if the residual still exceeds the standard after tripping, it is determined that the circuit breaker fails to operate or the fault point is located on the power supply side of the circuit breaker (i.e., the fault has not been isolated), and the alarm level is immediately upgraded and manual emergency handling is notified.

[0557] ④ Closed-loop feedback and model optimization

[0558] The complete process data of each alarm and its handling is recorded in a cloud database, including:

[0559] Abnormal feature data packets

[0560] Cloud-based diagnostic results and confidence levels

[0561] Final fault type and cause confirmed manually

[0562] Processing time and personnel

[0563] Treatment measures and effects

[0564] These data are used periodically (once a month):

[0565] Knowledge base update: Newly confirmed fault cases and their actual fault type labels are added to the knowledge base, and the mean vector and standard deviation vector of each fault type are recalculated to make the knowledge base continuously rich and accurate.

[0566] Baseline model optimization: Add newly accumulated normal operation data to the training set and retrain the CNN-GRU model to adapt it to changes in the electricity consumption behavior of the distribution area (such as user capacity increase or decrease, seasonal changes, equipment replacement, etc.).

[0567] Threshold parameter optimization: Statistical analysis of false positives and false negatives, adjusting p in the POT method. FA The (false alarm rate) parameter achieves a better balance between detection sensitivity and false alarm rate.

[0568] Ultimately, a complete closed-loop mechanism of "monitoring → alarm → diagnosis → processing → feedback → optimization" is formed, enabling the system to continuously self-optimize and evolve over time, and continuously improve detection performance.

[0569] System Architecture Overview

[0570] The anomaly monitoring and early warning system involved in this invention consists of the following five physical and logical modules:

[0571] 1. Data Acquisition Module

[0572] Deployment location: Outgoing line side of low-voltage distribution transformer and key branch lines

[0573] Hardware components: High-precision voltage / current sensor (ADE9000 metering chip + Rogowski coil CT), DS18B20 temperature sensor, 4G communication module, ARM Cortex-M7 main control MCU

[0574] Functions: Acquires instantaneous three-phase voltage / current values ​​at a sampling rate of 6.4kHz, records effective values, power parameters, temperature, etc. at 1-minute intervals, and periodically acquires meteorological data via 4G network.

[0575] 2. Data Preprocessing and Model Building Module

[0576] Deployment location: Cloud server (training phase) / Local terminal (inference phase)

[0577] Functions: Perform VMD decomposition, feature matrix construction, and training and inference of CNN-GRU models.

[0578] 3. Edge computing module

[0579] Deployment location: Smart terminal in the distribution area

[0580] Functions: Real-time data preprocessing → VMD decomposition → CNN-GRU forward inference → residual calculation → POT adaptive threshold generation → multi-parameter cross-validation → anomaly detection and feature packaging

[0581] 4. Cloud-based analytics platform

[0582] Deployment location: Cloud server

[0583] Function: Receives abnormal feature data packets uploaded by the terminal, calls the fault feature knowledge base to perform Mahalanobis distance matching diagnosis, assesses the severity, generates graded early warning instructions, and sends them to the terminal to execute alarm actions.

[0584] 5. Early warning output module

[0585] Deployment locations: Substation side (audible and visual alarms, circuit breakers) and cloud side (SMS gateway, APP push, voice calls).

[0586] Functions: Executes corresponding alarm actions based on the warning level, including audible and visual alarms, SMS push notifications, telephone notifications, and automatic tripping. Example 1: Construction of a baseline model based on VMD-CNN-GRU

[0587] This embodiment uses historical operating data from a real low-voltage distribution area (capacity 500kVA, number of power supply users 328) as an example to explain in detail the construction process of the benchmark model.

[0588] Step A: Data Acquisition and Preprocessing

[0589] Operational data from the past 12 months was extracted from the intelligent monitoring terminal of the transformer substation, including 12 electrical parameters such as three-phase voltage, three-phase current, active power, reactive power, power factor, frequency, and temperature. The data sampling interval was 1 minute, resulting in a total of 12 × 30 × 24 × 60 ≈ 518,400 records.

[0590] Missing value imputation was performed: consecutive missing values ​​due to communication interruption for a maximum of 7 minutes (3 times in total) were detected and filled using segmented cubic Hermite interpolation; 127 isolated outliers were detected (0.024%) and corrected using Hampel filter (k=5, n0=3).

[0591] Perform Z-score standardization and statistically obtain the mean and standard deviation of each parameter.

[0592] Step B: VMD Decomposition

[0593] With K=8 and α=2000, VMD decomposition is performed on the time series data for each dimension. Taking phase A current as an example, the decomposed data consists of 8 IMF components:

[0594] IMF1 (center frequency approximately 0.35Hz, corresponding to a period of approximately 2.9 minutes): mainly contains high-frequency random fluctuations in the load;

[0595] IMF4 (center frequency approximately 0.008 Hz, corresponding to a period of approximately 2.08 hours): mainly includes the morning and evening peak variations of daily load;

[0596] IMF8 (center frequency approximately 0.00005Hz, corresponding to a period of approximately 5.8 days): mainly includes weekday / restday cycle patterns.

[0597] The study verified that VMD decomposition can effectively separate electricity consumption behavior characteristics at different time scales.

[0598] Step C: Model Training

[0599] A sliding window of samples was constructed (H=1440, i.e., 24 hours), generating approximately 517,000 samples in total. The samples were divided into a training set (first 8.4 months, approximately 362,000 samples), a validation set (middle 1.8 months, approximately 77,500 samples), and a test set (last 1.8 months, approximately 77,500 samples).

[0600] The training was performed using a CNN (two Conv1D layers, 64 and 128 filters, 12 and 8 kernel sizes) + bidirectional GRU (128 hidden units) + GRU (64 hidden units) + fully connected layer (12 outputs).

[0601] Training process: After approximately 45 epochs, the Huber loss on the validation set stopped decreasing (from approximately 0.85 initially to approximately 0.12), triggering early stopping. The model has approximately 3.2 million parameters, and a single forward inference operation on the local terminal (STM32H743, 480MHz) takes approximately 18ms, meeting the timing requirement of 1 minute of real-time detection.

[0602] Step D: Model Evaluation

[0603] Evaluate model performance on the test set:

[0604] Mean Absolute Prediction Error (MAE): Voltage approximately 0.3V (0.14% error rate relative to 220V), Current approximately 0.5A (1.0% error rate relative to average current of 50A), Active Power approximately 0.8kW (1.0% error rate relative to average power of 80kW).

[0605] Compared with ARIMA: The MAE of the CNN-GRU model is about 35% lower than that of ARIMA, indicating that deep learning models have a significant advantage in capturing complex nonlinear electricity consumption behavior.

[0606] Example 2: Anomaly Detection Based on POT Adaptive Threshold

[0607] This embodiment simulates anomaly detection on the test set based on the benchmark model trained in Embodiment 1.

[0608] Step A: Residual Calculation

[0609] For each sample in the test set, a normalized residual sequence of 12 electrical parameters was calculated. The test set contains approximately 77,500 time points (approximately 54 days of data).

[0610] Step B: Dynamic Threshold Generation

[0611] For the normalized residual sequence of phase A current, setting L=1440 (24 hours), p FA =1 / 1440. Tracking dynamic threshold changes over 54 days:

[0612] Between 2:00 AM and 5:00 AM (the off-peak period, with an average current of approximately 20A and relatively small fluctuations), the tail of the extreme value sequence is thinner (ξ≈0.05), and the dynamic threshold T... dynamic It stabilizes around 2.1σ;

[0613] Between 18:00 and 21:00 (peak load period, average current approximately 80A, with significant fluctuations), the tail of the extreme value sequence thickens considerably (ξ≈0.25), and the dynamic threshold T... dynamic It automatically increases to approximately 3.5σ.

[0614] Verification: If a fixed threshold (set to 3.0σ) is used, a large number of normal fluctuations will be falsely reported as abnormal during peak load periods (false alarm rate of about 5%); while the adaptive threshold always keeps the actual false alarm rate close to the preset 1 / 1440≈0.07, verifying the effectiveness of the POT method.

[0615] Step C: Cross-validation

[0616] In the 54 days of data in the test set, a total of 13 valid abnormal events were detected, all of which were genuine abnormal events (manually verified): 8 were overloads, 3 were three-phase imbalances, 1 was a voltage dip, and 1 was a sensor failure. Meanwhile, 27 initial anomalies were rejected by cross-validation rules (22 of which were single-parameter anomalies, and 5 were two-parameter anomalies that were physically unrelated), thus avoiding false alarms.

[0617] Example 3: Cloud-based Fault Diagnosis Case Study

[0618] The characteristic data packet content of a certain valid abnormal event (occurrence time: June 15, 2024, 14:23):

[0619] List of abnormal parameters: [I a I b U a U b [I0] (The current and voltage of both phases A and B are abnormal, and the zero-sequence current is abnormal)

[0620] Symbol pattern S: [-1, -1, 0, +1, +1, 0, +1]

[0621] Following Example 3: Cloud-based Fault Diagnosis Case

[0622] Relative magnitude vector R (partial): I a 4.8, I b : 4.5, I0: 12.3, U a 2.5, U b 2.3

[0623] Rise time: 0.3 seconds (extremely short, indicating that the fault is developing very rapidly)

[0624] Duration: In progress (uploaded 2 minutes ago, not yet resumed)

[0625] 50Hz fundamental frequency component amplitude: extremely large (significant power frequency component of short-circuit current)

[0626] Diagnostic process:

[0627] The feature vectors are calculated using Mahalanobis distance with the fault type templates in the knowledge base in the cloud. The results are shown in Table 17 below.

[0628] Table 17: Results Statistics Table

[0629]

[0630] The maximum similarity is 0.556 (FT-02 phase-to-phase short circuit fault), which exceeds the threshold of 0.35. Therefore, it is diagnosed as "B and C phase-to-phase short circuit fault" (the currents of phases A and B are abnormal, but according to phase analysis, it is actually a short circuit fault of phases B and C, and the feature vector gives this result).

[0631] Severity assessment:

[0632] max R i = 12.3 (Zero-sequence current deviation ratio; in reality, this should be the short-circuit phase current, taken as 4.8)

[0633] However, the fault type is phase-to-phase short circuit (D) type =1.0), the equipment is located on the main branch (D equip =0.6)

[0634] Duration: 2 minutes (duration=2)

[0635] Severity=0.4×min(1,4.8 / 5.0)+0.2×min(1,2 / 60)+0.2×0.6+0.2×1.0=0.4×0.96+0.2×0.033+0.2×0.6+0.2×1.0=0.384+0.0067+0.12+0.2=0.71;

[0636] Severity = 0.71 ≥ 0.7, therefore it is determined to be a Level 3 emergency alarm.

[0637] Linked actions:

[0638] The command is sent to the terminal in the distribution area to activate the audible and visual alarm.

[0639] The system automatically dials the duty manager's number and announces: "A phase-to-phase short circuit fault has occurred in the XX transformer area. Please handle it immediately."

[0640] A trip command was sent to the corresponding branch circuit breaker (trip time ≤ 200ms), successfully isolating the fault.

[0641] If a backup power source is available, the system will automatically switch to power supply to the non-faulty area.

[0642] From the occurrence of the anomaly to the completion of fault isolation, the entire process takes about 3.5 seconds (including a 1-minute waiting period for confirmation, terminal upload, cloud diagnosis, command issuance, and circuit breaker action), which is much faster than the reaction time of manual inspection and traditional protection devices.

[0643] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for abnormal monitoring and early warning based on electricity consumption behavior data of low-voltage distribution areas, characterized in that, Includes the following steps: S1: Collect multi-dimensional time-series power consumption data of low-voltage distribution areas within a preset historical period. The multi-dimensional time-series power consumption data includes three-phase voltage, three-phase current, active power, reactive power, power factor, frequency, zero-sequence current, and ambient temperature. S2: Preprocess the multi-dimensional time-series electricity consumption data to construct a sample dataset. Based on the sample dataset, use variational mode decomposition to decompose the time-series data of each dimension into several intrinsic mode function components. Combine the components with time coding features, holiday identifiers, and meteorological features to form an input feature matrix. Use a hybrid deep learning model of convolutional neural network and gated recurrent unit to train the input feature matrix and generate a typical electricity consumption behavior benchmark model. S3: Collect the current time-series electricity consumption data of the low-voltage distribution area in real time, perform the same VMD decomposition and feature matrix construction as in S2, and input it into the typical electricity consumption behavior benchmark model. Calculate the difference between the expected normal value output by the model and the actual observed value, and generate a multidimensional residual sequence. S4: For the multidimensional residual sequence, the peak value exceeding the threshold method based on extreme value theory is used to dynamically calculate the detection threshold of each dimension. When the residual of a certain dimension exceeds the dynamic threshold, it is marked as an initial detection anomaly. The initial detection anomaly is verified by the multi-parameter cross-validation rule to confirm the valid anomaly event. S5: Upload the feature data of confirmed valid abnormal events to the cloud. The cloud performs pattern matching and fault type diagnosis based on the pre-built fault feature knowledge base, assesses the severity of the fault, and generates graded early warning signals.

2. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 1, characterized in that: In step S2, when using VMD for multi-scale time-frequency decomposition, the number of decomposition modes K=8 and the quadratic penalty factor α=2000 are set. Each IMF component and its center frequency are solved iteratively by the alternating direction method of multipliers.

3. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 2, characterized in that: The CNN-GRU hybrid model in step S2 consists of a CNN module composed of two one-dimensional convolutional layers and two max pooling layers. The first convolutional layer has 64 kernels and a kernel size of 12, and the second convolutional layer has 128 kernels and a kernel size of 8. The GRU module consists of a bidirectional GRU layer and a unidirectional GRU layer. The bidirectional GRU has 128 hidden units, and the unidirectional GRU has 64 hidden units.

4. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 3, characterized in that, The dynamic threshold generation method based on POT in step S4 includes: setting a sliding window length L = 1440 time steps, extracting peak sequences where the absolute value of the residuals within the window exceeds the initial threshold u0 of the 90th percentile, fitting the shape parameter ξ and scale parameter β of the generalized Pareto distribution using maximum likelihood estimation, and then calculating the false alarm rate p. FA =1 / 1440 Reverse-push state threshold T dynamic .

5. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 4, characterized in that, The multi-parameter cross-validation rules in step S4 include: calculating the correlation confidence C when the normalized residuals of at least two electrical parameters simultaneously exceed their respective dynamic thresholds, and these parameters are marked as correlated in a preset physical correlation rule table. ij If C ij Greater than threshold C th If the value is 1.0, it is confirmed as a valid abnormal event; if only a single parameter is abnormal, it is determined to be noise interference.

6. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 5, characterized in that: The fault feature knowledge base in step S5 contains fingerprint templates for at least eight types of typical faults. Each fault template consists of a feature mean vector, a feature standard deviation vector, and a matching similarity threshold. The feature vector includes 36 dimensions of features, including residual sign pattern, relative amplitude, time-series waveform features, and steady-state changes after the fault.

7. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 6, characterized in that, The graded early warning signals in step S5 include three levels: Level 1 general attention warning corresponds to events with Severity < 0.3, which is logged and displayed on the monitoring screen; Level 2 important attention warning corresponds to events with Severity ≤ 0.3 < 0.7, which is sent to maintenance personnel via SMS or APP push notification; Level 3 emergency alarm corresponds to events with Severity ≥ 0.7, which triggers audible and visual alarms, automatically dials a phone number to notify, and links the intelligent circuit breaker to perform trip protection.

8. The abnormal monitoring and early warning method based on low-voltage distribution area electricity consumption behavior data according to claim 7, characterized in that, Step S5 also includes fault elimination confirmation and closed-loop feedback: after tripping, the residual sequence is continuously monitored. If the residual recovers to below the threshold, the fault is confirmed to have been isolated. The complete alarm and processing data is recorded in the cloud for periodic updates to the fault feature knowledge base and retraining of the benchmark model.

Images