Current transformer state analyzing and monitoring method based on interval modeling

Abstract

The invention discloses a current transformer state analysis monitoring method based on interval modeling, and the method comprises the steps: collecting the ratio error data of an electronic current transformer, carrying out the preprocessing, and dividing the preprocessed data into a training set, a test set, and a verification set; the method comprises the following steps of: constructing a hierarchical deep learning model BITCN-BiLSTM-MHA; the preprocessed data are input into the constructed hierarchical deep learning model BITCN-BiLSTM-MHA, and model parameters after training are obtained; according to the trained model, measuring errors of the current transformer at different moments in the future are predicted, and deterministic prediction results of different time periods are generated; and according to the prediction result, carrying out probability distribution modeling on the prediction result by using an adaptive window width kernel density function, and generating prediction intervals of different confidence intervals. According to the method, in a dynamic complex signal environment, the variation trend of the measurement ratio error interval of the electronic current transformer in the future can be well predicted.

Description

Technical Field

[0001] The present invention relates to the technical field of evaluation of electronic current transformers, and particularly to a method for analyzing and monitoring the state of current transformers based on interval modeling. Background Art

[0002] Electronic current transformers are basic measurement devices in the power system. By measuring high-precision voltages, they provide support for power metering and relay protection, etc. Their measurement accuracy is crucial for the operation stability and safety of the power system. However, during long-term operation, the complex industrial environment and equipment aging have led to the gradual accumulation of measurement errors in electronic current transformers, affecting their measurement accuracy and posing a potential threat to the safety and reliability of the power grid. Therefore, it is necessary to evaluate the measurement accuracy of electronic current transformers.

[0003] The traditional modeling method is to accurately model the transmission line or system to predict the measurement error of the current transformer. Although the error estimation can be achieved theoretically, with the wide access of new energy, the access of intermittent power sources such as wind power and photovoltaic power has exacerbated the randomness of non-stationary signals in the power system. The traditional modeling and simulation methods are difficult to distinguish the gradual measurement errors caused by the deterioration of the current transformer performance and the fluctuations of the sensing signals, and have poor adaptability to dynamic and multimodal signals, seriously affecting the accuracy and reliability of error evaluation.

[0004] To overcome the above limitations, time series modeling methods represented by recurrent neural networks (RNNs) and long short-term memory networks (LSTMs) have emerged in recent years. By analyzing the non-linear characteristics of data and dealing with the temporal dependence, significant results have been achieved in both the new energy generation side (photovoltaic power generation, wind power) and the load side (power load forecasting). These technologies use the historical data of the equipment to predict the state information at future moments, improving the timeliness of risk warning, maintenance and repair, etc. However, directly applying deep learning technology to EVT error prediction still faces many challenges: First, the existing models have limited ability to model the global characteristics of the non-linear dynamic signals of the current transformer and cannot fully explore the internal laws of complex signals, resulting in poor prediction ability. Second, the outputs of the above models are all deterministic results, that is, the output is a single predicted value. However, due to the inevitable uncertainty of neural network prediction errors, there is always a deviation between the predicted value and the actual value. Summary of the Invention

[0005] To solve the problem of insufficient accuracy caused by prediction errors in traditional deterministic prediction methods, the present invention provides a method for analyzing and monitoring the state of current transformers based on interval modeling. This method can better predict the future change trend of the measurement error (ratio error) interval of electronic current transformers in a dynamic and complex signal environment.

[0006] The technical solution adopted by the present invention is as follows:

[0007] A method for analyzing and monitoring the state of a current transformer based on interval modeling, comprising the following steps:

[0008] Step 1: Collect the ratio error data of the electronic current transformer and perform preprocessing. Divide the preprocessed data into a training set, a test set, and a validation set;

[0009] Step 2: Construct a hierarchical deep learning model BITCN-BiLSTM-MHA;

[0010] Step 3: Input the data preprocessed in Step 1 into the hierarchical deep learning model

[0011] BITCN-BiLSTM-MHA constructed in Step 2 to obtain the trained model parameters;

[0012] Step 4: According to the model trained in Step 3, predict the measurement errors of the current transformer at different future times, and generate deterministic prediction results for different time periods;

[0013] Step 5: According to the prediction results in Step 4, use an adaptive window width kernel density function to model the probability distribution of the prediction results and generate prediction intervals with different confidence intervals.

[0014] In the said Step 1, the preprocessing includes the following steps:

[0015] S1.1. Eliminate abnormal data points and missing values greater than three times the standard deviation in the collected data of the electronic current transformer, and interpolate the vacant values after elimination by the linear interpolation method;

[0016] For the linear interpolation method, the calculation formula is as follows:

[0017]

[0018] In the formula, x i and x i+1 are the values at the corresponding time indices t i and t i+1 of the known points respectively; t is the time index point to be interpolated, satisfying t i ≤t≤t i+1 ; x is the interpolated value.

[0019] S1.2. Use an improved variational mode decomposition model to decompose the data processed in S1.1 to obtain the intrinsic mode components and the decomposition residuals.

[0020] S1.3. Eliminate the decomposition residuals in S1.2. At this time, the processed data consists of M + 3 columns, where M represents the number of intrinsic mode components after decomposition, and the remaining 3 columns are the interpolated specific difference data in S1.1, that is, the data before decomposition, and the corresponding date and specific time; the preprocessing is completed.

[0021] Divide the preprocessed data into a training set, a test set, and a validation set.

[0022] In S1.2, the construction of the improved variational mode decomposition model includes the following steps:

[0023] S1.2.1: Variational mode decomposition (VMD) has significant advantages in overcoming mode aliasing and endpoint effects and shows excellent performance in the field of signal filtering. Its core goal is to decompose the real-valued input signal x(t) into several discrete sub-signals u k (t), and iteratively match the optimal center frequency and bandwidth of each mode in the frequency domain. Its mathematical expression is:

[0024]

[0025] In the formula, u k (t) represents the transformation of the modal signal at time t, k represents the number of components, δ(t) is the Hilbert transform, j represents the analytic signal, represents the cyclic frequency, {u k} and {w k} represent the k modal functions and their center frequencies respectively, represents a complex exponential function used to capture the oscillation characteristics at frequency ω k and time t; x(t) represents the real-valued input signal, that is, the input data, and in the present invention, it refers to the interpolated specific difference data in S1.1.

[0026] S1.2.2: To solve the constrained variational problem, an augmented Lagrangian function is introduced, and the expression is:

[0027]

[0028] In the formula, α represents the quadratic penalty factor, λ represents the Lagrange multiplier, and L({u k}, {w k}, λ) represents converting the constrained problem in Equation (2) into an unconstrained problem; represents the inner product of the Lagrange multiplier and the constraint condition.

[0029] S1.2.3: Use the alternating direction multiplier algorithm and Fourier isometric transformation for u k and w k , and iterate repeatedly until the cut-off condition is met;

[0030]

[0031] where ε represents the cut-off precision, represents the frequency domain representation of the k-th mode at the (n + 1)-th iteration; represents the frequency domain representation of the k-th mode at the n-th iteration; g is a summation index representing the summation range from 1 to k.

[0032] S1.2.4: Select the optimization algorithm PRO as the optimization model for variational mode decomposition (VMD). According to the formulas in S1.2.1 - S1.2.3, in addition to selecting k and α, ε and the fidelity constraint G are additionally introduced tau ;

[0033] Set the above parameters: the decomposition number k, the penalty factor α, the cut-off precision ε, and the fidelity constraint G tau , and transform the decomposition problem of variational mode decomposition (VMD) into a four-dimensional optimization problem. The conditional formula is as follows:

[0034]

[0035] where f VMD-fitness represents performing VMD decomposition on the input data; k min represents the minimum number of intrinsic mode components; k max represents the maximum number of intrinsic mode components; α min represents the minimum value of the penalty factor; α max represents the maximum value of the penalty factor; represents the minimum value of the fidelity constraint; represents the maximum value of the fidelity constraint; ε min represents the minimum value of the cut-off precision represented by ε, max represents the maximum value of the cut-off precision. After setting the value ranges of the above parameters, use the optimization algorithm PRO to solve in order to obtain the optimal values of the above 4 parameters.

[0036] S1.2.5: In the optimal values obtained during the PRO solution of the VMD decomposition process in S1.2.4, the objective function needs to be set. In order to ensure the independence of each component, the maximum information coefficient (MIC) is selected as the basic objective function. In order to ensure the independence of information between each mode function and their relevance to the original information (data before decomposition), while minimizing the MIC between each component, we maximize the MIC between each component and the original sequence x(t) at the same time. The improved calculation formula is as follows:

[0037]

[0038] where fMIC (x(t), S) represents the improved objective function; mode(i) represents the i-th intrinsic mode sub-component after decomposition; mode(i + 1) represents the (i + 1)-th intrinsic mode sub-component after decomposition; mode(j) represents the j-th intrinsic mode sub-component after decomposition; MIC(mode(i), mode(i + 1)) represents the MIC value between the i-th intrinsic mode sub-component and its adjacent (i + 1)-th intrinsic mode sub-component; MIC(mode(j), x(t)) represents the MIC value between the j-th intrinsic mode sub-component and the original sequence x(t); S is a parameter of the decomposition algorithm and can be a matrix. This fitness function is used to evaluate the quality of the decomposition effect and guide parameter optimization.

[0039] S1.2.6: Based on the above S1.2.1 - S1.2.5, the solution problem of PRO-VMD can be transformed into the following formula:

[0040]

[0041] In the formula, (k * , α * , ε * , G * tau ) is the optimal parameter combination under the minimum f MIC , where k * represents the optimal decomposition number after solution; α * represents the value of the optimal penalty factor after solution; ε * represents the value of the optimal cut-off precision after solution; G * tau represents the value of the optimal fidelity constraint after solution.

[0042] In step 2, the hierarchical deep learning model BITCN-BiLSTM-MHA includes:

[0043] The BiTCN module is used for feature extraction;

[0044] The BiLSTM module is used to capture the time-dependent relationship of the input sequence;

[0045] The MHA module is used to reduce information forgetting during the processing of the BiLSTM module;

[0046] After the MHA module, there is a fully connected layer, which is used to output a definite prediction result.

[0047] In step 2, the construction of the hierarchical deep learning model BITCN-BiLSTM-MHA includes the following steps:

[0048] First, as the first layer of the model, the BiTCN module extracts deep temporal features from the improved variational mode decomposition data, and its calculation formula is as follows:

[0049]

[0050] In the formula, x(t T ) represents the input time series; y(t) represents the output features after convolution; w i represents the convolution kernel weight; r represents the dilation rate, which is used to control the expansion of the receptive field; k er represents the size of the convolution kernel; b represents the bias term. When multiple BiTCN layers are stacked, the output can be recursively defined as:

[0051] z (l+1) = f(z (l) ; θ (l) ) + z (l)

[0052] In the formula, z represents the output after residual connection; z (l) represents the output of the l-th layer; θ (l) represents the parameter set of the l-th layer; f represents the combined operation of causal convolution and dilation convolution.

[0053] Secondly, the multi-dimensional tensor output by the BiTCN module is flattened as follows:

[0054] Convert the output high-dimensional tensor into a two-dimensional tensor, and its formula is as follows:

[0055] T flattened = reshape(T);

[0056] In the formula, represents the input tensor, with dimensions of time step n TCN , number of channels d TCN , and feature dimension k TCN .

[0057] After flattening the tensor

[0058] The high-dimensional data is converted into a two-dimensional tensor and then passed to the BiLSTM module in the second layer for further processing. The calculation formula of the BiLSTM module is as follows:

[0059] h t = f(W f ·x t + U f ·h t-1 + b f )

[0060] ht = f(W b ·x t + U b ·h t-1 + b b )

[0061] where h t represents the forward hidden state; h t represents the backward hidden state; W f , U f , b f represent the forward LSTM parameters; W b , U b , b b represent the backward LSTM parameters; x t represents the input sequence; f represents the activation function. The output of BiLSTM concatenates the forward and backward states:

[0062]

[0063] where represents the concatenation operation.

[0064] Again, based on the output of the BiLSTM module, MHA can focus on key feature points from a global perspective through parallel computing of multiple attention heads.

[0065]

[0066] head i = Attention(QW i Q , KW i K , VW i V )

[0067] where Q, K, and V are the query matrix, key matrix, and value matrix respectively; h Att is the number of attention heads; W i Q , W i K , W i V , W o represent the projection matrices for the head and output; Concat refers to the concatenation operation; Attention refers to the operation of calculating attention scores.

[0068] Finally, connect the fully connected layer as the final output layer of the model. The fully connected layer receives the processing results from the MHA module and generates the final prediction value of the model through mapping. Its calculation formula is as follows:

[0069]

[0070] In the formula, h MHA represents the output of the MHA; W ol represents the weight matrix; b o represents the bias; σ represents the activation function.

[0071] In step 3, the training set in S1.4 is input into the BiTCN module. After passing through the BiTCN module, the BiLSTM module, the MHA module, and the fully connected layer, the predicted ratio difference is output. The trained hierarchical deep learning model BiTCN - BiLSTM - MHA is verified using the test set and the validation set, and the model under the optimal index training parameters is obtained by means of multiple iterations.

[0072] In step 4, the test set is aggregated into a 1 - hour period. To achieve the period conversion, the original acquisition data with a 10 - minute period is converted into 1 - hour period data by calculating the mean value. The calculation formula is as follows:

[0073]

[0074] In the formula, x(t i ) is the 1 - hour period data after conversion; x(t i ) is the 10 - minute period data before conversion; where i ranges from 0 to 5, representing 6 data points within the hour. t i is the start time of the current 1 - hour period. The mean value processing can effectively smooth the data fluctuations and more clearly reflect the change of the overall trend of the ratio difference.

[0075] Step 5 includes the following steps:

[0076] S5.1. Calculate the probability density distribution of the prediction error, and use adaptive window - width kernel density estimation to estimate the obtained probability density distribution respectively.

[0077] S5.2. Set different confidence interval scores and calculate the prediction intervals under different confidence intervals.

[0078] S5.3. Superimpose the prediction intervals under different confidence intervals to obtain the final ratio difference prediction interval.

[0079] Step 5 includes the following steps:

[0080] 1). Calculate the probability density distribution of the prediction error. First, define the PDF as f(x) and the cumulative distribution function as F(x). Then, for a given sample ξ, then:

[0081]

[0082] Wherein, f(ξ) represents the probability density value of the sample point ξ, which is used to measure the possibility of the sample appearing at ξ; h is the window width parameter of the kernel function, which determines the degree of smoothing; F(ξ + h) represents the value of the cumulative distribution function at the point ξ + h, indicating the probability that the random variable takes a value less than or equal to ξ + h; F(ξ - h) represents the value of the cumulative distribution function at the point ξ - h, indicating the probability that the random variable takes a value less than or equal to ξ - h.

[0083] 2), Assume that the data set to be estimated by the mutual inductor is X, X = {ξ 1, ξ2, ξ3, …, ξ n}, ξ 1, ξ2, ξ3, …, ξ n respectively represent the independent and identically distributed samples of the mutual inductor data; by introducing the indicator function I(·), the probability density function can be obtained

[0084]

[0085] Wherein, ξ i is an independent and identically distributed sample point; n is the sample size; K(·) is a kernel function that satisfies continuity, non-negativity, and integration equal to 1. Here, the Cauchy kernel is selected.

[0086] 3), First, minimize the mean integrated square error as the optimization objective, and solve the optimal window width h opt within the domain to provide a preliminary reference for density estimation. The solution formula is as follows:

[0087]

[0088] Wherein, represents the mean integrated square of the error between the estimated density and the true density f, reflecting the accuracy of the estimation. represents the expected value of the square integral of the density estimation error, measuring the global difference between the estimated density and the true density.

[0089] By transforming the MISE optimization problem into finding its extreme point, the cross-validation method is used to solve the estimation error to obtain the global optimal window width h opt , and the simplified expression CV(h) of the cross-validation method is as follows:

[0090]

[0091] Select the initial bandwidth h min and h max , and define the step size Δh. For each candidate window width h k ∈[h min , h max , calculate CV(h k) Select CV(h k ) The smallest h k As the optimal window width h opt :

[0092]

[0093] 4), Based on the obtained global optimal window width, calculate the initial density estimate value of each data point

[0094]

[0095] 5), According to the initial density estimate value, calculate the adaptive window width of each data point:

[0096]

[0097] In the formula, λ h Is an adjustment parameter used to adjust the reference window width. In the present invention, the adjustment parameter is determined by an improved golden section search. h(ξ i ) Is the adaptive window width at ξ i ; τ is a sensitivity control parameter, taking a positive value, used to adjust the sensitivity of the window width to the change of regional density, so as to reflect local characteristics.

[0098] 6), Finally, based on the mathematical expression constructed by adaptive kernel density estimation Is:

[0099]

[0100] The different confidence interval scores mentioned refer to 97.5%, 95%, 90% confidence intervals.

[0101] In order to avoid the golden section method falling into a local optimal solution, the golden section method is improved by adding perturbations at each division point calculation to jump out of the local trap:

[0102]

[0103] In the formula, Represents the trial point close to the left end point of the interval corresponding to the golden section ratio in the nth iteration; Represents the trial point close to the right end point of the interval corresponding to the golden section ratio in the nth iteration; a (n) Represents the left end point of the interval in the nth iteration; b (n) The right end point of the interval in the nth iteration; And Are random perturbations, following a normal distribution, that is As the number of iterations increases, gradually reduce the perturbation amplitude:

[0104] σ n = σ0·exp(-λ δ n)(23);

[0105] where σ n represents the perturbation amplitude in the nth iteration; σ0 is the initial perturbation amplitude, where σ0 ~ 0.01·(b - a); a represents the left endpoint of the initial interval; b represents the right endpoint of the initial interval; λ δ is the decay rate, where λ δ ∈ [0.05, 0.15]; n represents the current iteration number.

[0106] The method of adding local interpolation is used to further improve the accuracy, and the steps are as follows:

[0107] 1), Perform local smooth sampling on the obtained new interval [a (n+1) , b (n+1) :

[0108]

[0109] where h i represents the window width of the ith local sampling point; m h is the total number of sampling points; i is the index of the sampling point, i = 1, 2,, m h ; a (n+1) represents the left endpoint of the interval obtained in the (n + 1)th iteration; b (n+1) represents the right endpoint of the interval obtained in the (n + 1)th iteration.

[0110] 2), After calculating the objective function f(h i ) at the sampling points, perform moving average smoothing:

[0111]

[0112] where is the estimated value of the objective function after moving average at the ith sampling point, used to smooth the noise; w' j is the smoothing weight, and a uniform weight is used here; k h is the smoothing window size; f(h i+j ) is the value of the objective function at the (i + j)th sampling point, reflecting the function value at h i+j .

[0113] 3), Use the smoothed sampling point values for interpolation optimization. In the first step, use cubic spline interpolation to fit the smoothed objective function values:

[0114]

[0115] 4), Calculate the minimum point: Let S′(h) = 0 to obtain the minimum point h min . If h min ∈[a (n+1) , b (n+1) , and S(h min ) < min{f(a (n+1) ), f(b (n+1) ), f(x1 (n+1) ), f(x2 (n+1) )}, then accept h min as the final optimization result; otherwise, the interpolation is invalid.

[0116] In the formula, f(a (n+1) ) represents the value of the objective function at the left endpoint a (n+1) of the interval; f(b (n+1) ) represents the value of the objective function at the right endpoint b (n+1) of the interval; f(x1 (n+1) ) and f(x2 (n+1) ) respectively represent the values of the objective function at the current golden section points.

[0117] 5), In the case of invalid interpolation, if |b (n+1) - a (n+1) | > δ / 5, continue the golden section iteration to further narrow the interval. If |b (n+1) - a (n+1) | ≤ δ / 5, then return the midpoint (a (n+1) + b (n+1) ) / 2 of the interval as the final output. δ is the stopping condition of the golden section method.

[0118] The technical effects of a current transformer status analysis and monitoring method based on interval modeling according to the present invention are as follows:

[0119] 1) A current transformer status analysis and monitoring method based on interval modeling according to the present invention aims to solve the problem of insufficient prediction accuracy caused by prediction errors in traditional deterministic prediction methods. This method analyzes the distribution characteristics of prediction results by using adaptive kernel density estimation, dynamically adjusts the bandwidth of kernel density estimation, and constructs a prediction interval to reflect the future change trend of the measurement error of electronic current transformers. This interval prediction improves the reliability and accuracy of prediction results by showing the fluctuation range of the measurement error of electronic current transformers, and provides a more accurate basis for judging the early faults of electronic current transformers.

[0120] 2) The present invention improves variational mode decomposition, transforming the cut-off accuracy, penalty factor, number of components, and fidelity constraint that affect the decomposition effect into an optimization problem. Through refined decomposition, the non-linear measurement ratio error under the influence of dynamic complex signals is decomposed into simple modal components, while the decomposition residuals are eliminated, thereby reducing the difficulty of the prediction model in learning non-stationary ratio error signals.

[0121] 3) The present invention enhances the modeling ability of the input signal through a hierarchical BiTCN-BiLSTM-MHA model. This model accurately captures the non-linear characteristics and time dependence of the signal, thereby improving the accuracy of error prediction.

[0122] 4) The present invention generates error prediction intervals with different confidence intervals through an adaptive kernel density estimation method, enabling the prediction results to not only include a single deterministic value but also reflect the probability distribution characteristics of measurement errors, effectively quantifying the uncertainty of the prediction.

[0123] 5) The present invention improves the traditional golden section method by introducing random perturbations to avoid local optimal traps. At the same time, by combining local smooth sampling and cubic spline interpolation, the accuracy of bandwidth selection is further improved, thereby optimizing the smoothing effect of kernel density estimation.

[0124] 6) The present invention provides high-confidence interval error assessment results, providing strong data support for the operation and maintenance management, fault prediction, and health assessment of power grid equipment, and supporting the measurement and control requirements of future intelligent power systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0125] The present invention will be further described below in conjunction with the drawings and examples;

[0126] Figure 1 It is a flowchart of a current transformer status analysis and monitoring method based on interval modeling.

[0127] FIG. 2(a) is a prediction result diagram of BiTCN-BiLSTM-MHA under a 10-minute sampling period (deterministic prediction result);

[0128] FIG. 2(b) is a prediction result diagram of BiTCN-BiLSTM-MHA under a 1-hour sampling period (deterministic prediction result).

[0129] FIG. 3(a) is a distribution diagram of the uncertainty of the prediction result over time;

[0130] FIG. 3(b) is a probability density function distribution diagram of the uncertainty of the prediction result;

[0131] The sampling periods corresponding to FIGS. 3(a) and 3(b) are 10 minutes.

[0132] Figure 4(a) shows the prediction interval of the adaptive kernel function estimation under the improved golden section method Figure 1 ;

[0133] Figure 4(b) shows the second figure of the prediction interval of the adaptive kernel function estimation under the improved golden section method;

[0134] In Figure 4(a) and Figure 4(b), the 97.5%, 95% and 90% confidence intervals are superimposed, and the corresponding sampling periods are 10 minutes and 1 hour.

[0135] Figure 5 It is the structural schematic diagram of BiTCN-BiLSTM-MHA.

[0136] Figure 6 It is the result diagram of the prediction interval from the perspective of evaluation and monitoring. Specific implementation manners

[0137] For the current transformer status analysis and monitoring method based on interval modeling, first, the ratio error data of the electronic current transformer is collected and preliminarily preprocessed. Then, the improved variational mode decomposition method is used to decompose the data into different intrinsic mode components and remove the residuals, thereby reducing the modeling difficulty of the prediction model for non-stationary signals. Next, a deep learning model BiTCN-BiLSTM-MHA based on a hierarchical structure is constructed to enhance the modeling ability for the time dependence and non-linear characteristics of the input signal. Finally, to further quantify the uncertainty of the prediction results, the present invention adopts the method of adaptive kernel density estimation, and generates prediction intervals under different confidence intervals by dynamically adjusting the window width of the kernel function to characterize the probability distribution characteristics of the prediction error of the deep learning model. In addition, the golden section optimization algorithm combines random perturbation and local interpolation techniques, significantly improving the smoothing effect of the kernel density estimation. The present invention provides reliable data support for the operation and maintenance management, fault prediction and health assessment of power grid equipment.

[0138] The current transformer status analysis and monitoring method based on interval modeling includes the following steps:

[0139] Step 1: Collect the ratio error data of the electronic current transformer, perform preprocessing, and divide the preprocessed data into a training set, a test set, and a validation set;

[0140] Step 2: Construct the hierarchical deep learning model BITCN-BiLSTM-MHA;

[0141] Step 3: Input the data preprocessed in Step 1 into the hierarchical deep learning model BITCN-BiLSTM-MHA constructed in Step 2 to obtain the trained model parameters;

[0142] Step 4: According to the model trained in Step 3, predict the measurement errors of the current transformer at different future times, and generate deterministic prediction results for different time periods;

[0143] Step 5: According to the prediction results in Step 4, use the adaptive window width kernel density function to model the probability distribution of the prediction results and generate prediction intervals with different confidence intervals.

[0144] In Step 1, since noise may interfere with the true ratio error data of the electronic current transformer during the acquisition process, preprocessing is required. The preprocessing includes the following steps:

[0145] S1.1: Eliminate the abnormal data points and missing values greater than three standard deviations in the collected data of the electronic current transformer, and interpolate the vacant values after elimination by the linear interpolation method;

[0146] For the linear interpolation method, the calculation formula is as follows:

[0147]

[0148] In the formula, x i and x i+1 are the values at the corresponding time indices t i and t i+1 of the known points respectively; t is the time index point to be interpolated, satisfying t i ≤t≤t i+1 ; x is the interpolated value.

[0149] S1.2: Use the improved variational mode decomposition model to decompose the data processed in S1.1 to obtain the intrinsic mode components and the decomposition residuals.

[0150] S1.3: Eliminate the decomposition residuals in S1.2. At this time, the processed data consists of M + 3 columns, where M represents the number of intrinsic mode components after decomposition, and the remaining 3 columns are the ratio error data after interpolation in S1.1, that is, the data before decomposition, as well as the corresponding date and the corresponding specific time; the preprocessing ends.

[0151] Divide the preprocessed data into a training set, a test set, and a validation set.

[0152] In S1.2, the construction of the improved variational mode decomposition model includes the following steps:

[0153] S1.2.1: Variational mode decomposition (VMD) has significant advantages in overcoming mode aliasing and endpoint effects and shows excellent performance in the field of signal filtering. Its core goal is to decompose the real-valued input signal x(t) into several discrete sub-signals u k(t), and the optimal center frequency and bandwidth of each mode are matched through frequency-domain iteration. Its mathematical expression is:

[0154]

[0155] In the formula, u k (t) represents the transformation of the modal signal at time t, k represents the number of components, δ(t) is the Hilbert transform, j represents the analytic signal, represents the cyclic frequency, {u k} and {w k} represent the k modal functions and their center frequencies respectively, represents a complex exponential function used to capture the oscillation characteristics at frequency ω k and time t; x(t) represents the real-valued input signal, that is, the input data, and in the present invention, it refers to the interpolated ratio difference data in S1.1.

[0156] S1.2.2: To solve the constrained variational problem, an augmented Lagrangian function is introduced, and the expression is:

[0157]

[0158] In the formula, α represents the quadratic penalty factor, λ represents the Lagrange multiplier, L({u k}, {w k}, λ) represents converting the constrained problem in Equation (2) into an unconstrained problem; represents the inner product of the Lagrange multiplier and the constraint condition.

[0159] S1.2.3: For u k and w k , the alternating direction multiplier algorithm and Fourier equidistant transformation are adopted, and iterative calculations are performed repeatedly until the cut-off condition is satisfied;

[0160]

[0161] In the formula, ε represents the cut-off accuracy, represents the frequency-domain representation of the kth mode at the (n + 1)th iteration; represents the frequency-domain representation of the kth mode at the nth iteration; g is a summation index representing the summation range from 1 to k.

[0162] S1.2.4: Select the optimization algorithm PRO as the optimization model for variational mode decomposition VMD. According to the formulas in S1.2.1 to S1.2.3, in addition to selecting k and α, ε and the fidelity constraint G tau are additionally introduced;

[0163] Set the above parameters: the decomposition number k, the penalty factor α, the cut-off accuracy ε, and the fidelity constraint Gtau , the decomposition problem of variational mode decomposition (VMD) is transformed into a four-dimensional optimization problem, and the conditional formula is as follows:

[0164]

[0165] In the formula, f VMD-fitness represents performing VMD decomposition on the input data; k min represents the minimum number of intrinsic mode components; k max represents the maximum number of intrinsic mode components; α min represents the minimum value of the penalty factor; α max represents the maximum value of the penalty factor; represents the minimum value of the fidelity constraint; represents the maximum value of the fidelity constraint; ε min represents the minimum value of the cut-off precision, ε max represents the maximum value of the cut-off precision. After setting the value ranges of the above parameters, use the optimization algorithm PRO to solve, in order to obtain the optimal values of the above 4 parameters.

[0166] S1.2.5: In the optimal values during the PRO solution of the VMD decomposition process in S1.2.4, the objective function needs to be set. To ensure the independence of each component, the maximum information coefficient (MIC) is selected as the basic objective function. To ensure the independence of the information between each modal function and their relevance to the original information (data before decomposition), on the basis of minimizing the MIC between each component, we simultaneously maximize the MIC between each component and the original sequence x(t). The improved calculation formula is as follows:

[0167]

[0168] In the formula, f MIC (x(t), S) represents the improved objective function; mode(i) represents the i-th intrinsic mode sub-component after decomposition; mode(i + 1) represents the (i + 1)-th intrinsic mode sub-component after decomposition; mode(j) represents the j-th intrinsic mode sub-component after decomposition; MIC(mode(i), mode(i + 1)) represents the MIC value between the i-th intrinsic mode sub-component and its adjacent (i + 1)-th intrinsic mode sub-component; MIC(mode(j), x(t)) represents the MIC value between the j-th intrinsic mode sub-component and the original sequence x(t); S is a parameter of the decomposition algorithm and can be a matrix. This fitness function is used to evaluate the quality of the decomposition effect and guide parameter optimization.

[0169] S1.2.6: Based on the above S1.2.1 - S1.2.5, the solution problem of PRO-VMD can be transformed into the following formula:

[0170]

[0171] In the formula, (k * , α * , ε * , G * tau ) is the optimal parameter combination under the minimum f MIC , where k * represents the optimal number of decompositions after solution; α * represents the value of the optimal penalty factor after solution; ε * represents the value of the optimal cut-off accuracy after solution; G * tau represents the value of the optimal fidelity constraint after solution.

[0172] The PRO algorithm described in the above features is an evolutionary optimization algorithm called Partial Reinforcement Optimizer proposed in 2023. This algorithm updates the behavior priority by simulating the partial reinforcement process to obtain a more effective response, thereby improving the solution efficiency of global optimization problems.

[0173] The maximum information coefficient (MIC) described in the above features is a calculation method based on the concept of mutual information, used to capture the non-linear relationship and non-functional dependence between data. Its discrete form is as follows:

[0174]

[0175] In the formula, A and B are random variables, p(A,B) is the joint probability distribution of A and B, and p(A) and p(B) are the marginal probabilities of A and B respectively. The MIC algorithm divides the values of A and B into d A and d B two intervals and combines them into a grid G. Calculate the mutual information on each grid, then the maximum mutual information of the fixed grid G can be defined as:

[0176] I * (A,B,d A ,d B ) = maxI(A,B|G) (9);

[0177] In the formula, A,B|G represents the sample space divided by the grid G. Normalize the maximum mutual information of each grid and record it in the feature matrix:

[0178]

[0179] Therefore, MIC is defined as

[0180]

[0181] Where N is the number of samples, and B(N) represents the upper limit of the grid size.

[0182] In the said step 2, the hierarchical deep learning model BITCN-BiLSTM-MHA includes:

[0183] The BiTCN module is used for feature extraction;

[0184] The BiLSTM module is used to capture the temporal dependencies of the input sequence;

[0185] The MHA module is used to reduce information forgetting during the processing of the BiLSTM module;

[0186] A fully connected layer is connected after the MHA module, which is used to output the determined prediction result.

[0187] The structural schematic diagram of the hierarchical deep learning model BITCN-BiLSTM-MHA is as Figure 5 shown.

[0188] In the said step 2, the construction of the hierarchical deep learning model BITCN-BiLSTM-MHA includes the following steps:

[0189] First, as the first layer of the model, the BiTCN module extracts the deep temporal features of the improved variational mode decomposition data, and its calculation formula is as follows:

[0190]

[0191] Where x(t T ) represents the input time series; y(t) represents the output feature after convolution; w i represents the convolution kernel weight; r represents the dilation rate, which is used to control the expansion of the receptive field; k er represents the size of the convolution kernel; b represents the bias term. When multiple layers of BiTCN are stacked, the output can be recursively defined as:

[0192] z (l+1) = f(z (l) ; θ (l) ) + z (l)

[0193] Where z represents the output after residual connection; z (l) represents the output of the l-th layer; θ (l) represents the parameter set of the l-th layer; f represents the combined operation of causal convolution and dilated convolution.

[0194] The BiTCN module significantly expands the perception range through causal convolution and dilated convolution, enabling the model to capture key features over a longer time span. The introduction of a multi-layer residual structure alleviates the problem of vanishing gradients that may occur during the training of deep networks and effectively avoids the adverse effects of over-smoothing on the feature expression ability.

[0195] Secondly, the multi-dimensional tensor output by the BiTCN module is flattened as follows:

[0196] The high-dimensional tensor of the output is converted into a two-dimensional tensor, and its formula is expressed as follows:

[0197] T flattened = reshape(T);

[0198] In the formula, represents the input tensor with dimensions of time step n TCN , number of channels d TCN , and feature dimension k TCN .

[0199] The flattened tensor

[0200] After converting the high-dimensional data into a two-dimensional tensor, it is passed to the BiLSTM module in the second layer for further processing. The calculation formula of the BiLSTM module is as follows:

[0201] h t = f(W f ·x t + U f ·h t-1 + b f )

[0202] h t = f(W b ·x t + U b ·h t-1 + b b )

[0203] In the formula, h t represents the forward hidden state; h t represents the backward hidden state; W f , U f , b f represent the forward LSTM parameters; W b , U b , b b represent the backward LSTM parameters; x t represents the input sequence; f represents the activation function. The output of the BiLSTM concatenates the forward and backward states:

[0204]

[0205] In the formula, represents the splicing operation.

[0206] The BiLSTM module makes full use of the forward and backward information of the time series with its bidirectional cyclic characteristics, enhancing the ability to capture time dependencies. By processing bidirectional time series information in parallel, the BiLSTM not only optimizes the model's information flow control ability but also expands the model's information reception range, and provides high-quality feature inputs for subsequent modules.

[0207] Again, based on the output of the BiLSTM module, the MHA can focus on key feature points from a global perspective through parallel calculations of multiple attention heads. During the process of the BiLSTM module processing time series information, some important information may be weakened or lost due to long-term dependencies or information redundancy, and the introduction of the MHA enhances the model's ability to capture this information. The calculation formula of the MHA is as follows:

[0208]

[0209] head i = Attention(QW i Q , KW i K , VW i V )

[0210] In the formula, Q, K, and V are the query matrix, key matrix, and value matrix respectively; h Att is the number of attention heads; W i Q , W i K , W i V , W o represent the projection matrices of the head and the output; Concat refers to the splicing operation; Attention refers to the operation of calculating the attention score.

[0211] Finally, connect the fully connected layer as the final output layer of the model. The fully connected layer receives the processing results from the MHA module and generates the final prediction value of the model through mapping. The calculation formula is as follows:

[0212]

[0213] In the formula, h MHA represents the output of the MHA; W ol represents the weight matrix; b o represents the bias; σ represents the activation function.

[0214] In step 3, the training set in S1.4 is input into the BiTCN module. After passing through the BiTCN module, the BiLSTM module, the MHA module, and the fully connected layer, the predicted ratio difference is output. The trained hierarchical deep learning model BiTCN - BiLSTM - MHA is verified using the test set and the validation set, and the model under the optimal metric training parameters is obtained through multiple iterations.

[0215] In step 4, the test set is aggregated into a 1 - hour period. To achieve the period conversion, the original acquisition data with a 10 - minute period is converted into 1 - hour period data by calculating the mean value. The calculation formula is as follows:

[0216]

[0217] In the formula, x(t i ) is the converted 1 - hour period data; x(t i ) is the 10 - minute period data before conversion; where i ranges from 0 to 5, representing 6 data points within the hour. t i is the start time of the current 1 - hour period. Mean processing can effectively smooth data fluctuations and more clearly reflect the change in the overall trend of the ratio difference.

[0218] Step 5 includes the following steps:

[0219] S5.1. Calculate the probability density distribution of the prediction error, and use adaptive window - width kernel density estimation to estimate the obtained probability density distribution respectively.

[0220] S5.2. Set different confidence interval scores and calculate the prediction intervals under different confidence intervals.

[0221] S5.3. Superimpose the prediction intervals under different confidence intervals to obtain the final ratio - difference prediction interval.

[0222] The adaptive window - width kernel density estimation is a non - parametric method for estimating the probability density function (PDF). The core lies in revealing the probability distribution law based on the characteristics of the data itself, avoiding prior assumptions about the distribution form.

[0223] Step 5 includes the following steps:

[0224] 1). Calculate the probability density distribution of the prediction error. First, define the PDF as f(x) and the cumulative distribution function as F(x). Then, for the given sample ξ, then:

[0225]

[0226] In the formula, f(ξ) represents the probability density value of the sample point ξ, which is used to measure the possibility of the sample appearing at ξ; h is the window width parameter of the kernel function, which determines the degree of smoothing; F(ξ + h) is the value of the cumulative distribution function at the point ξ + h, representing the probability that the random variable takes a value less than or equal to ξ + h; F(ξ - h) represents the value of the cumulative distribution function at the point ξ - h, representing the probability that the random variable takes a value less than or equal to ξ - h.

[0227] 2), Assume that the data set to be estimated by the mutual inductor is X, X = {ξ 1, ξ2, ξ3, …, ξ n}, ξ 1, ξ2, ξ3, …, ξ n respectively represent independent and identically distributed samples of the mutual inductor data; by introducing the indicator function I(·), the probability density function can be obtained

[0228]

[0229] In the formula, ξ i is an independent and identically distributed sample point; n is the sample size; K(·) is a kernel function that satisfies continuity, non-negativity, and integral equal to 1. Here, the Cauchy kernel is selected.

[0230] 3), Adaptive window width kernel density estimation can dynamically adjust the window width h according to the local density characteristics of the data, avoiding possible underestimation of density in sparse regions and possible over-smoothing in high-density regions. First, the minimum mean square integral error is used as the optimization objective to solve the optimal window width h within the domain of definition opt , to provide a preliminary reference for density estimation. The solution formula is as follows:

[0231]

[0232] In the formula, represents the mean square integral of the error between the estimated density and the true density f, reflecting the accuracy of the estimation. represents the expected value of the square integral of the density estimation error, measuring the global difference between the estimated density and the true density.

[0233] By transforming the MISE optimization problem into finding its extreme point, the cross-validation method is used to solve the estimation error to obtain the global optimal window width h opt , and the simplified expression of the cross-validation method CV(h) is as follows:

[0234]

[0235] Select the initial bandwidth h min and h max , and define the step size Δh. For each candidate window width hk ∈ [h min , h max , calculate CV(h k ). Select the h for which CV(h k ) is the smallest k as the optimal window width h opt :

[0236]

[0237] 4), Based on the obtained global optimal window width, calculate the initial density estimate value of each data point

[0238]

[0239] 5), According to the initial density estimate value, calculate the adaptive window width of each data point:

[0240]

[0241] In the formula, λ h is an adjustment parameter used to adjust the reference window width. In the present invention, the adjustment parameter is determined by an improved golden section search. h(ξ i ) is the adaptive window width at ξ i ; τ is a sensitivity control parameter, taking a positive value, used to adjust the sensitivity of the window width to the change of regional density, so as to reflect local characteristics.

[0242] 6), Finally, based on the mathematical expression constructed by adaptive kernel density estimation is:

[0243]

[0244] The different confidence interval scores mentioned refer to the 97.5%, 95%, and 90% confidence intervals.

[0245] The so-called golden section search aims to divide the interval using the golden section ratio within the search interval, compare the function values at the division points, and gradually iterate to narrow the range to approach the minimum value of the objective function. Let the initial interval of the local bandwidth be [a, b]. For the nth iteration, the calculation formula is as follows:

[0246]

[0247] If f(x1 (n) ) < f(x2 (n) ), update the interval to If f(x1 (n) ) ≥ f(x2 (n) ), then update the interval to Set the termination condition to ò. When the updated new interval |b (n+1) -a (n +1) | < ò, terminate the iteration, and the obtained optimal solution is:

[0248] h optimal =(a (n+1) +b (n+1) ) / 2(20);

[0249] To avoid the golden section method falling into a local optimal solution, improve the golden section method by adding perturbations at each calculation of the division point to jump out of the local trap:

[0250]

[0251] In the formula, represents the trial point close to the left endpoint of the interval corresponding to the golden section ratio in the nth iteration; represents the trial point close to the right endpoint of the interval corresponding to the golden section ratio in the nth iteration; a (n) represents the left endpoint of the interval in the nth iteration; b (n) the right endpoint of the interval in the nth iteration; and are random perturbations, following a normal distribution, that is As the number of iterations increases, gradually reduce the perturbation amplitude:

[0252] σ n =σ0·exp(-λ δ n)(23);

[0253] In the formula, σ n represents the perturbation amplitude in the nth iteration; σ0 is the initial perturbation amplitude, where σ0 ~ 0.01·(b - a); a represents the left endpoint of the initial interval; b represents the right endpoint of the initial interval; λ δ is the attenuation rate, where λ δ ∈[0.05, 0.15]; n represents the current number of iterations.

[0254] Further improve the accuracy by adding the method of local interpolation. The steps are as follows:

[0255] 1), Perform local smooth sampling on the obtained new interval [a (n+1) , b (n+1) :

[0256]

[0257] In the formula, h i represents the window width of the i-th local sampling point; m his the total number of sampling points; i is the index of the sampling point, i = 1, 2,..., m h ; a (n+1) represents the left endpoint of the interval obtained in the (n + 1)-th iteration; b (n+1) represents the right endpoint of the interval obtained in the (n + 1)-th iteration.

[0258] 2), Calculate the objective function f(h i ) at the sampling points and then perform moving average smoothing:

[0259]

[0260] In the formula, is the estimated value of the objective function after moving average at the i-th sampling point, used to smooth the noise; w' j is the smoothing weight, and uniform weight is used here; k h is the smoothing window size; f(h i+j ) is the value of the objective function at the (i + j)-th sampling point, reflecting the function value at h i+j .

[0261] 3), Use the smoothed sampling point values for interpolation optimization. In the first step, use cubic spline interpolation to fit the smoothed objective function values:

[0262]

[0263] 4), Calculate the minimum point: Let S′(h) = 0 to obtain the minimum point h min . If h min ∈[a (n+1) , b (n+1) , and S(h min ) < min{f(a (n+1) ), f(b (n+1) ), f(x1 (n+1) ), f(x2 (n+1) )}, then accept h min as the final optimization result, otherwise the interpolation is invalid.

[0264] In the formula, f(a (n+1) ) represents the value of the objective function at the left endpoint a (n+1) of the interval; f(b (n+1) ) represents the value of the objective function at the right endpoint b (n+1) of the interval; f(x1 (n+1) ) and f(x2 (n+1) ) respectively represent the values of the objective function at the current golden section points.

[0265] 5), In the case of invalid interpolation, if |b (n+1) - a(n+1) | > ε / 5, continue the golden section iteration to further narrow the interval. If |b (n+1) - a (n+1) | ≤ ε / 5, then return the midpoint of the interval (a (n+1) + b (n+1) ) / 2 as the final output. ε is the stopping condition of the golden section method.

[0266] Figure 2(a) is the prediction result graph of BiTCN - BiLSTM - MHA under a 10 - minute sampling period (deterministic prediction result);

[0267] Figure 2(b) is the prediction result graph of BiTCN - BiLSTM - MHA under a 1 - hour sampling period (deterministic prediction result).

[0268] It can be seen from Figure 2(a) and Figure 2(b) that the prediction using the improved VMD combined with BiTCN - BiLSTM - MHA is higher, and compared with directly using BiTCN - BiLSTM - MHA, the prediction result is closer to the true value.

[0269] Figure 3(a) shows the distribution graph of the difference between the predicted value and the true value of the prediction model over time, and Figure 3(b) shows the corresponding probability density distribution.

[0270] It can be seen from Figure 4(a) and Figure 4(b) that compared with the deterministic prediction results in Figure 2(a) and Figure 2(b), the present invention further constructs a probability interval of the model prediction error based on the deterministic prediction results. The interval almost contains all the prediction errors, providing more interpretability and scientificity when making decisions.

[0271] Figure 6 It is the prediction interval result graph (1 - hour week) under the monitoring perspective. According to the national standard regulations, the measurement deviation of a 0.2 - class current transformer shall not exceed 0.2%. Figure 6 The predicted value is the red line. By the method of the present invention, prediction intervals under different confidence intervals are constructed. This interval band is the deviation distribution between the predicted value and the true value. Therefore, it can provide more comprehensive and reliable results.

Claims

1. A method for analyzing and monitoring the state of a current transformer based on interval modeling, characterized in that It includes the following steps: Step 1: Collect the ratio error data of the electronic current transformer, and perform preprocessing. Divide the preprocessed data into a training set, a test set, and a validation set; Step 2: Construct the hierarchical deep learning model BITCN - BiLSTM - MHA; Step 3: Input the data preprocessed in Step 1 into the hierarchical deep learning model BITCN - BiLSTM - MHA constructed in Step 2 to obtain the trained model parameters; Step 4: According to the model trained in Step 3, predict the measurement errors of the current transformer at different future moments, and generate deterministic prediction results for different time periods; Step 5: According to the prediction results in Step 4, use the adaptive window width kernel density function to perform probability distribution modeling on the prediction results, and generate prediction intervals with different confidence intervals.

2. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 1, wherein: In the said Step 1, the preprocessing includes the following steps: S1.1: Eliminate the abnormal data points and missing values greater than three times the standard deviation in the collected data of the electronic current transformer, and interpolate the vacant values after elimination by the linear interpolation method; For the linear interpolation method, the calculation formula is as follows: where x i and x i+1 are the values of the known points at the corresponding time indices t i and t i+1 respectively; t is the time index point to be interpolated, satisfying t i ≤ t ≤ t i+1 ; x is the interpolated value; S1.2: Use the improved variational mode decomposition model to decompose the data processed in S1.1 to obtain the intrinsic mode components and the decomposition residual; S1.3: Eliminate the decomposition residual in S1.

2. At this time, the processed data consists of M + 3 columns, where M represents the number of the intrinsic mode components after decomposition, and the other 3 columns are the ratio error data after interpolation in S1.1, that is, the data before decomposition, as well as the corresponding date and the corresponding specific time; The preprocessing ends. Divide the preprocessed data into a training set, a test set, and a validation set.

3. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 2, wherein: In the said S1.2, the construction of the improved variational mode decomposition model includes the following steps: S1.2.1: Variational Mode Decomposition (VMD) is used to decompose the real-valued input signal x(t) into several discrete sub-signals u k (t), and iteratively match the optimal center frequency and bandwidth of each mode in the frequency domain; its mathematical expression is as follows: where u k (t) represents the transformation of the modal signal at time t, k represents the number of components, δ(t) is the Hilbert transform, and j represents the analytic signal, represents the cyclic frequency, {u k} and {w k} represent the k modal functions and their central frequencies respectively, represents a complex exponential function used to capture the oscillation characteristics at frequency ω k and time t; x(t) represents the real-valued input signal, i.e., the input data; S1.2.2: In order to solve the constrained variational problem, an augmented Lagrangian function is introduced, and the expression is: where α represents the quadratic penalty factor, λ represents the Lagrange multiplier, and L({u k},{w k},λ) represents transforming the constrained problem in Equation (2) into an unconstrained problem; represents the inner product of the Lagrange multiplier and the constraint condition; S1.2.3: For u k and w k Use the alternating direction multiplier algorithm and the Fourier equidistant transformation, and iterate repeatedly until the cut-off condition is met; where ε represents the cut-off accuracy, represents the frequency-domain representation of the k-th mode at the (n + 1)-th iteration; represents the frequency-domain representation of the k-th mode at the n-th iteration; g is a summation index representing the summation range from 1 to k; S1.2.4: Select the optimization algorithm PRO as the optimization model of variational mode decomposition (VMD). According to the formulas in S1.2.1 - S1.2.3, in addition to selecting k and α, ε and the fidelity constraint G are additionally introduced. tau ; Set the above-mentioned parameter decomposition number \(k\), penalty factor \(\alpha\), cut-off accuracy \(\varepsilon\) and fidelity constraint \(G\). tau Then, transform the decomposition problem of variational mode decomposition (VMD) into a four-dimensional optimization problem. The conditional formula is as follows: Where, f VMD-fitness represents performing VMD decomposition on the input data; k min represents the minimum number of intrinsic mode components; k max represents the maximum number of intrinsic mode components; α min represents the minimum value of the penalty factor; α max represents the maximum value of the penalty factor; represents the minimum value of the fidelity constraint; represents the maximum value of the fidelity constraint; ε min represents the minimum value of the cut-off accuracy, ε max represents the maximum value of the cut-off accuracy; S1.2.5: In order to ensure the independence of each component, the maximum information coefficient (MIC) is selected as the basic objective function; In order to ensure the independence of the information between each modal function and their relevance to the data before decomposition, on the basis of minimizing the MIC between each component, we maximize the MIC between each component and the original sequence x(t) at the same time. The improved calculation formula is as follows: where f MIC (x(t), S) represents the improved objective function; mode(i) represents the i-th intrinsic mode sub-component after decomposition; mode(i + 1) represents the (i + 1)-th intrinsic mode sub-component after decomposition; mode(j) represents the j-th intrinsic mode sub-component after decomposition; MIC(mode(i), mode(i + 1)) represents the MIC value between the i-th intrinsic mode sub-component and its adjacent (i + 1)-th intrinsic mode sub-component; MIC(mode(j), x(t)) represents the MIC value between the j-th intrinsic mode sub-component and the original sequence x(t); S is the parameter of the decomposition algorithm; S1.2.6: Based on the above S1.2.1 - S1.2.5, the solution problem of PRO - VMD can be transformed into the following formula: In the formula, (k * , α * , ε * , G * tau ) is the optimal parameter combination under the minimum f MIC , where k * represents the optimal number of decompositions after solution; α * represents the value of the optimal penalty factor after solution; ε * represents the value of the optimal cut-off accuracy after solution; G * tau represents the value of the optimal fidelity constraint after solution.

4. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 1, wherein: In the said Step 2, the construction of the hierarchical deep learning model BITCN - BiLSTM - MHA includes the following steps: First, the BiTCN module, as the first layer of the model, extracts the deep temporal features of the improved variational mode decomposition data, and its calculation formula is as follows: where \(x(t T )\) represents the input time series; \(y(t)\) represents the output feature after convolution; \(w i \) represents the convolution kernel weight; \(r\) represents the dilation rate, which is used to control the expansion of the receptive field; \(k er \) represents the size of the convolution kernel; \(b\) represents the bias term; when multiple layers of BiTCN are stacked, the output can be recursively defined as: z (l+1) = f(z (l) ; θ (l) ) + z (l) where z represents the output after residual connection; z (l) represents the output of the l-th layer; θ (l) represents the set of parameters of the l-th layer; f represents the combined operation of causal convolution and dilated convolution; Secondly, the multi - dimensional tensor output by the BiTCN module is flattened, specifically as follows: Convert the output high-dimensional tensor into a two-dimensional tensor, and its formula is expressed as follows: In the formula, represents the input tensor with dimensions of time step n TCN , number of channels d TCN , and feature dimension k TCN ; the flattened tensor After converting the high - dimensional data into a two - dimensional tensor, it is passed to the BiLSTM module in the second layer for further processing; The calculation formula of the BiLSTM module is as follows: In the formula, represents the forward hidden state; represents the backward hidden state; W f , U f , b f represent the forward LSTM parameters; W b , U b , b b represent the backward LSTM parameters; x t represents the input sequence; f represents the activation function; the output of the BiLSTM concatenates the forward and backward states: In the formula, represents a splicing operation; Thirdly, based on the output of the BiLSTM module, MHA can focus on the key feature points from a global perspective through the parallel calculation of multiple attention heads; The calculation formula of MHA is as follows: Where Q, K, and V are the query matrix, key matrix, and value matrix respectively; h Att is the number of attention heads; W o represents the projection matrix for the head and the output; Concat refers to the concatenation operation; Attention refers to the operation of calculating attention scores; Finally, connect the fully connected layer as the final output layer of the model; the fully connected layer receives the processing results from the MHA module and generates the final predicted value of the model through mapping. The calculation formula is as follows: Where h MHA represents the output of the MHA; W ol represents the weight matrix; b o represents the bias; σ represents the activation function.

5. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 4, wherein: In step 3, the training set is input into the BiTCN module, and after passing through the BiTCN module, the BiLSTM module, and the MHA module, the predicted ratio difference is output by the fully connected layer. The trained BiTCN-BiLSTM-MHA is verified using the test set and the validation set, and the model under the optimal metric training parameters is obtained through multiple iterations.

6. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 1, characterized in that: In step 4, to achieve cycle conversion, the original acquisition data with a 10-minute cycle is converted into 1-hour cycle data by calculating the mean value. The calculation formula is as follows: In the formula, is the converted 1-hour cycle data; x(t i ) is the 10-minute cycle data before conversion; where i ranges from 0 to 5, representing 6 data points within that hour; t i is the start time of the current 1-hour period.

7. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 1, wherein: Step 5 includes the following steps: S5.

1. Calculate the probability density distribution of the prediction error, and use adaptive window width kernel density estimation to estimate the obtained probability density distribution respectively; S5.

2. Set different confidence interval scores and calculate the prediction intervals under different confidence intervals; S5.

3. Superimpose the prediction intervals under different confidence intervals to obtain the final ratio difference prediction interval.

8. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 1, wherein: Step 5 includes the following steps: 1). Calculate the probability density distribution of the prediction error. First, define the PDF as f(x) and the cumulative distribution function as F(x). Then, for a given sample ξ, then: In the formula, f(ξ) represents the probability density value of the sample point ξ, which is used to measure the possibility of the sample appearing at ξ; h is the window width parameter of the kernel function, which determines the smoothness; F(ξ + h) is the value of the cumulative distribution function at the point ξ + h, representing the probability that the random variable takes a value less than or equal to ξ + h; F(ξ - h) represents the value of the cumulative distribution function at the point ξ - h, representing the probability that the random variable takes a value less than or equal to ξ - h; 2), Let the data set to be estimated by the mutual inductor be X, X = {ξ 1, ξ2, ξ3, …, ξ n}, ξ 1, ξ2, ξ3, …, ξ n respectively represent the independent and identically distributed samples of the mutual inductor data; introducing the indicator function I(·), the probability density function can be obtained where ξ i is an independent and identically distributed sample point; n is the sample size; K(·) is a kernel function that satisfies continuity, non-negativity, and an integral of 1. 3) The adaptive window width kernel density estimation can dynamically adjust the window width h according to the local density characteristics of the data. First, the minimization of the mean square integral error is used as the optimization objective, and the optimal window width h within the domain of definition is solved as follows opt to provide a preliminary reference for density estimation. The solution formula is as follows: In the formula, represents the estimated density The mean square integral of the error between and the true density f reflects the accuracy of the estimation; represents the expected value of the squared integral of the density estimation error, which measures the global difference between the estimated density and the true density; By transforming the MISE optimization problem into finding its extreme points and using the cross-validation method to solve for the estimation error, the globally optimal window width h is obtained. opt , the simplified expression of the cross-validation method formula CV(h) is as follows: Select the initial bandwidth h min and h max , and define the step size Δh; for each candidate window width h k ∈[h min , h max , calculate CV(h k ); select the h k for which CV(h k ) is the smallest as the optimal window width h opt : 4) Calculate the initial density estimate value of each data point based on the obtained global optimal window width 5). Calculate the adaptive window width of each data point according to the initial density estimation value: where λ h is an adjustment parameter for adjusting the reference window width, and the adjustment parameter is determined by an improved golden section search; h(ξ i ) is the adaptive window width at ξ i ; τ is the sensitivity control parameter, taking a positive value, which is used to adjust the sensitivity of the window width to the regional density change, so as to reflect the local characteristics; 6), Finally, the mathematical expression constructed based on adaptive kernel density estimation is as follows:

9. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 8, wherein: To avoid the golden section method falling into the local optimal solution, improve the golden section method by adding perturbations at each division point calculation to jump out of the local trap: wherein, represents the trial point close to the left endpoint of the interval corresponding to the golden ratio in the n-th iteration; represents the trial point close to the right endpoint of the interval corresponding to the golden ratio in the n-th iteration; a (n) represents the left endpoint of the interval in the n-th iteration; b (n) the right endpoint of the interval in the n-th iteration; is a random perturbation, following a normal distribution, that is as the number of iterations increases, the perturbation amplitude gradually decreases: σ n = σ0·exp(-λ δ n) (23); where σ n represents the perturbation amplitude in the nth iteration; σ0 is the initial perturbation amplitude; a represents the left endpoint of the initial interval; b represents the right endpoint of the initial interval; λ δ is the attenuation rate; n represents the current iteration number.

10. The method for analyzing and monitoring the state of a current transformer based on interval modeling according to claim 8, characterized in that: The method of adding local interpolation is further used to improve the accuracy; it includes the following steps: 1) Perform local smooth sampling on the obtained new interval [a (n+1) , b (n+1) : where h i represents the window width of the i-th local sampling point; m h is the total number of sampling points; i is the index of the sampling point, i = 1, 2, …, m h ; a (n+1) represents the left endpoint of the interval obtained in the (n + 1)-th iteration; b (n+1) represents the right endpoint of the interval obtained in the (n + 1)-th iteration; 2), calculate the objective function f(h i ) and then perform moving average smoothing processing: In the formula, The estimated value of the objective function after moving average at the i-th sampling point, which is used to smooth the noise; w' j is the smoothing weight, and uniform weight is used here; k h is the smoothing window size; f(h i+j ) is the value of the objective function at the (i + j)-th sampling point, which reflects the function value at h i+j location; 3), Use the smoothed sampled point values Perform interpolation optimization. First step, use cubic spline interpolation to fit the smoothed objective function values: 4), Calculate the minimum point: Let S′(h) = 0 to obtain the minimum point h min ; If h min ∈[a (n+1) , b (n+1) , and S(h min ) < min{f(a (n+1) ), f(b (n+1) ), f(x1 (n+1) ), f(x2 (n+1) )}, then accept h min as the final optimization result; otherwise, the interpolation is invalid; Wherein, f(a (n+1) ) represents the value of the objective function at the left endpoint a of the interval (n+1) ; f(b (n+1) ) represents the value of the objective function at the right endpoint b of the interval (n+1) ; f(x1 (n+1) ) and f(x2 (n+1) ) respectively represent the values of the objective function at the current golden section points; 5) In the case where interpolation is invalid, if |b (n+1) - a (n+1) | > ε / 5, continue the golden section iteration to further narrow the interval; if |b (n+1) - a (n+1) | ≤ ε / 5, then return the midpoint of the interval (a (n+1) + b (n+1) ) / 2 as the final output; ε is the stopping condition of the golden section method.

Images