An adaptive interval prediction method for power load distribution drift

Abstract

The application discloses a kind of adaptive interval prediction methods for power load distribution drift, it is related to power load prediction technical field, including: collection history prediction residual set, based on history prediction residual set, construct local moment state vector, and utilize local moment state vector with adjacent distributed node to carry out information interaction, obtain neighborhood state information, and update history prediction residual set according to neighborhood state information, output calibration prediction residual set;Calculate the skewness coefficient of calibration prediction residual set, and based on skewness coefficient, through dynamic asymmetric kernel density estimation, the probability density of calibration prediction residual set is reconstructed, and dynamic quantile is output;Through conformal prediction, linear combination is carried out to dynamic quantile and power load prediction value, and power load prediction interval is generated.The application effectively supports the reliable dispatching and decision of power grid operation.

Description

Technical Field

[0001] This invention relates to the field of power load forecasting technology, and in particular to an adaptive interval forecasting method for power load distribution drift. Background Technology

[0002] With the deepening of smart grid construction and the integration of a high proportion of renewable energy into new power systems, the volatility and randomness of power load are exhibiting more complex non-stationary characteristics. In current industrial practice, the technical approaches for load forecasting mainly focus on point forecasting and probabilistic interval forecasting. Among these, statistical methods, represented by deep learning and ensemble learning, have shown significant advantages in handling high-dimensional time-series features. Traditional interval forecasting techniques largely rely on quantile regression or parameter estimation methods based on the assumption of a normal distribution. By mining the inherent evolutionary patterns of historical load curves, they construct confidence intervals that can cover future load fluctuations. These techniques, under the assumption of stationary stochastic processes, can effectively characterize load uncertainty, providing important auxiliary decision-making support for dispatching departments.

[0003] However, when faced with extreme weather changes, fluctuations in social activities, or distribution network topology reconfiguration, the temporal characteristics of load data will undergo nonlinear shifts, leading to severe skewness and heteroscedasticity in the distribution of prediction residuals. Existing static prediction architectures struggle to capture and quantify the dynamic evolution of this distribution drift in real time. Furthermore, when dealing with spatial multi-node correlations, they often perform isolated single-point predictions, lacking a neighborhood information interaction mechanism based on distributed control logic. This results in low interval coverage or excessively wide prediction intervals in distribution drift scenarios, failing to provide accurate prediction ranges while maintaining confidence levels. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides an adaptive interval prediction method for power load distribution drift to solve the problems of low accuracy and poor reliability of interval prediction.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] This invention provides an adaptive interval prediction method for power load distribution drift, comprising:

[0008] Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data to generate time-series feature vectors;

[0009] The time-series feature vector is input into the pre-trained power load prediction model for forward calculation, and the power load prediction value is output.

[0010] Collect historical prediction residuals, construct local moment state vectors based on historical prediction residuals, and use local moment state vectors to interact with neighboring distributed nodes to obtain neighborhood state information. Update the historical prediction residuals based on the neighborhood state information and output the calibration prediction residuals.

[0011] Calculate the skewness coefficient of the calibration prediction residual set, and reconstruct the probability density of the calibration prediction residual set based on the skewness coefficient through dynamic asymmetric kernel density estimation, and output the dynamic quantile;

[0012] By using conformal prediction, dynamic quantiles and power load forecast values ​​are linearly combined to generate power load forecast intervals.

[0013] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the specific steps for collecting multi-source heterogeneous data and preprocessing the multi-source heterogeneous data to generate time-series feature vectors are as follows:

[0014] Dimensionless processing is performed on multi-source heterogeneous data to output normalized feature data;

[0015] Interpolation and time-scale resampling are performed on the normalized feature data to output spatiotemporally aligned data;

[0016] Sliding window sampling is performed on spatiotemporally aligned data to output temporal feature vectors.

[0017] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the power load value prediction model is constructed based on a Riemannian geometry embedding layer, a manifold tangent space transformation layer, and a curvature adaptive aggregation layer.

[0018] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the specific steps of inputting the time-series feature vector into a pre-trained power load prediction model for forward calculation and outputting the power load prediction value are as follows:

[0019] The time-series feature vector is input into the power load prediction model, and a Riemannian geometric embedding layer is used to perform nonlinear spatial mapping on the time-series feature vector to output the manifold feature matrix.

[0020] By utilizing the manifold tangent space transformation layer, the manifold feature matrix is ​​locally linearized and projected to output a high-dimensional tangent vector sequence;

[0021] By utilizing a curvature adaptive aggregation layer, spatiotemporal correlation modeling and curvature correction are performed on high-dimensional tangent vector sequences to output predicted power load values.

[0022] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the steps of collecting historical prediction residual sets, constructing local moment state vectors based on these residual sets, and using these local moment state vectors to interact with neighboring distributed nodes to obtain neighborhood state information are as follows:

[0023] The lower-order raw moments and higher-order central moments are calculated based on the historical prediction residual set, and the lower-order raw moments and higher-order central moments are constructed into local moment state vectors through dimension concatenation.

[0024] Local dual variables are set based on the local moment state vector, and the local moment state vector and local dual variables are synchronized to the adjacent distributed nodes, and the neighborhood state information is received back.

[0025] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the specific steps for updating the historical prediction residual set based on neighborhood state information and outputting a calibrated prediction residual set are as follows:

[0026] An augmented Lagrangian function is constructed based on neighborhood state information, and the local moment state vector and local dual variable are iteratively updated using the augmented Lagrangian function to obtain the consensus moment vector.

[0027] Using the consensus moment vector as the moment matching constraint for the historical prediction residual set, a maximum entropy optimization problem is constructed.

[0028] The Lagrange multiplier method is used to solve the maximum entropy optimization problem and obtain the optimal probability weights.

[0029] The historical prediction residual set is resampled based on the optimal probability weights to generate a calibration prediction residual set.

[0030] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the specific steps of calculating the skewness coefficient of the calibration prediction residual set and reconstructing the probability density of the calibration prediction residual set based on the skewness coefficient through dynamic asymmetric kernel density estimation to output dynamic quantiles are as follows:

[0031] Using the mean of the calibration prediction residual set as a benchmark, the skewness coefficient of the calibration prediction residual set is calculated, and the AMISE objective function is constructed based on the skewness coefficient.

[0032] By searching for the global minimum, the AMISE objective function is jointly optimized, and the optimal left bandwidth and optimal right bandwidth are output.

[0033] Using the optimal left bandwidth and optimal right bandwidth as scaling parameters, an asymmetric Gaussian kernel function is selected to reconstruct the probability density of the calibration prediction residual set, generating the residual probability density function.

[0034] Numerical integration is performed on the residual probability density function to obtain the cumulative distribution function, and the upper and lower tail quantiles are solved for the cumulative distribution function to obtain the dynamic quantiles; the dynamic quantiles include the upper quantile and the lower quantile.

[0035] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the step of jointly optimizing the AMISE objective function by retrieving the global minimum value and outputting the optimal left bandwidth and optimal right bandwidth is as follows:

[0036] Based on the statistical distribution characteristics of the calibration prediction residual set, the range of optimal left bandwidth and optimal right bandwidth is defined, and the range is discretized at equal intervals to form a two-dimensional search grid.

[0037] The bandwidth combination in the two-dimensional search grid is substituted into the AMISE objective function for numerical calculation to generate the bandwidth loss matrix.

[0038] The global minimum is retrieved within the bandwidth loss matrix, and the bandwidth combination corresponding to the global minimum is extracted from the two-dimensional search grid to obtain the optimal left bandwidth and the optimal right bandwidth.

[0039] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the steps of using the optimal left bandwidth and optimal right bandwidth as scale parameters, and employing an asymmetric Gaussian kernel function to reconstruct the probability density of the calibration prediction residual set to generate a residual probability density function are as follows:

[0040] An asymmetric Gaussian kernel function is constructed with the residual samples in the calibration prediction residual set as the center and the optimal left bandwidth and optimal right bandwidth as the supporting radii.

[0041] The local probability contribution value of each residual sample is calculated using an asymmetric Gaussian kernel function, and the arithmetic mean and normalization of all local probability contribution values ​​are performed to generate the residual probability density function.

[0042] As a preferred embodiment of the adaptive interval prediction method for power load distribution drift described in this invention, the step of generating a power load prediction interval by linearly combining dynamic quantiles and power load prediction values ​​through conformal prediction is as follows:

[0043] The electricity load forecast is linearly summed with the lower quantile and the upper quantile respectively to form the initial forecast interval;

[0044] The lower and upper quantiles are used to construct the boundary of the residual interval, and for each residual sample in the calibration prediction residual set, the maximum non-negative deviation beyond the boundary of the residual interval is calculated as the inconsistency score.

[0045] Arrange all non-consistent scores in ascending order to form an ordered score sequence, and use the non-consistent scores at the quantile rank in the ordered score sequence as the conformal correction amount;

[0046] The conformal correction is linearly compensated to the lower and upper boundaries of the initial prediction interval to generate the power load prediction interval.

[0047] The beneficial effects of this invention are as follows: by deeply mining the nonlinear spatiotemporal characteristics of load data through Riemannian geometric embedding and curvature correction technology to improve the accuracy of point prediction, and by combining moment consensus collaboration and maximum entropy optimization under the distributed control concept to calibrate the prediction residuals in real time, and finally using dynamic asymmetric kernel density estimation and conformal prediction technology to reconstruct the residual distribution, so that even when the load distribution drifts, it can still generate an accurate prediction interval with rigorous statistical coverage guarantee, effectively supporting the reliable scheduling and decision-making of power grid operation. Attached Figure Description

[0048] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. 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.

[0049] Figure 1 This is a flowchart of an adaptive interval prediction method for power load distribution drift.

[0050] Figure 2 A flowchart for generating power load forecast intervals.

[0051] Figure 3 A flowchart for generating the calibration prediction residual set.

[0052] Figure 4 A flowchart for obtaining dynamic quantiles. Detailed Implementation

[0053] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0054] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0055] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0056] Reference Figures 1-4 This is one embodiment of the present invention, which provides an adaptive interval prediction method for power load distribution drift, comprising the following steps:

[0057] S1. Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data to generate time-series feature vectors.

[0058] S1.1 Perform dimensionless processing on multi-source heterogeneous data and output normalized feature data.

[0059] Furthermore, outlier detection is performed on the multi-source heterogeneous data to remove abnormal sampling points that exceed the preset statistical distribution range, thus obtaining pure heterogeneous data. A linear mapping is then performed on the pure heterogeneous data based on the maximum and minimum values ​​to output normalized feature data.

[0060] It should be noted that the statistical distribution interval is set based on the statistical characteristics of multi-source heterogeneous data.

[0061] Multi-source heterogeneous data includes power load data, meteorological data, substation and power plant equipment operation data from different sources, of different types and formats, as well as socio-economic data (power load is affected not only by weather and equipment status, but also by socio-economic factors such as industrial production, commercial operation and holidays).

[0062] S1.2 Perform interpolation and time-scale resampling on the normalized feature data to output spatiotemporally aligned data.

[0063] Furthermore, for the missing numerical locations in the normalized feature data, the weighted average of the neighboring non-missing points is used to fill in the missing values ​​and generate a complete feature sequence; the complete feature sequence is then fitted with cubic spline curves and resampled to output spatiotemporally aligned data.

[0064] S1.3 Perform sliding window sampling on the spatiotemporally aligned data and output the temporal feature vector.

[0065] Furthermore, continuous observation segments are extracted from the spatiotemporally aligned data according to a preset time span, and the observation segments are flattened to output a one-dimensional mapping sequence of multi-dimensional features; the one-dimensional mapping sequence is stacked into tensors to output a temporal feature vector.

[0066] It should be noted that the time span is set based on the sampling frequency of the spatiotemporally aligned data.

[0067] S2. Input the time-series feature vector into the pre-trained power load prediction model for forward calculation and output the power load prediction value.

[0068] Furthermore, the power load prediction model is constructed based on a Riemannian geometry embedding layer, a manifold tangent space transformation layer, and a curvature adaptive aggregation layer.

[0069] The overall architecture of the power load prediction model is as follows:

[0070] The power load prediction model consists of a Riemannian geometry embedding layer, a manifold tangent space transformation layer, and a curvature adaptive aggregation layer cascaded in sequence. The overall data stream is as follows: the input time-series feature vector first enters the Riemannian geometry embedding layer, and its output manifold feature matrix serves as the sole input to the manifold tangent space transformation layer. The high-dimensional tangent vector sequence output by this layer, together with the manifold feature matrix, serves as the input to the curvature adaptive aggregation layer, and finally outputs the predicted power load value.

[0071] The specific architecture of the power load prediction model is as follows:

[0072] The Riemannian geometry embedding layer contains a set of learnable metric kernels. Each metric kernel stores the product of a lower triangular matrix and its transpose as the kernel matrix, along with a scalar offset parameter. After receiving the input vector, the layer performs a quadratic operation on it with each metric kernel matrix. The resulting scalar is added to the corresponding offset and then passed through the hyperbolic tangent function to obtain an element. All the elements generated by the metric kernels are arranged in a row, and this process is repeated according to a preset embedding depth (the embedding depth is set based on the trade-off between the dimension of the temporal feature vector and the target rank of the manifold feature matrix) to generate multiple rows, thus obtaining the manifold feature matrix.

[0073] The manifold tangent space transformation layer maintains a tangent bundle basis dictionary and a connection coefficient tensor. The tangent bundle basis dictionary consists of multiple basis entries, each containing a centroid vector, a set of left tangent basis vectors, and a set of right tangent basis vectors. The connection coefficient tensor is a three-dimensional array, where the first two dimensions correspond to any two basis entries, and the third dimension stores the affine transformation parameters. For each row of the manifold feature matrix, this layer calculates its Euclidean distance to all centroids and selects the nearest basis entry. The local tangent representation is obtained by performing an inner product between the basis vector of this entry and the row vector. Subsequently, the local tangent representation is successively passed and aligned according to the row order using the affine parameters of the nearest adjacent basis entry pairs in the connection coefficient tensor. All the passed and aligned representations are stacked vertically to form a high-dimensional tangent vector sequence.

[0074] The curvature adaptive aggregation layer includes a discrete curvature calculator and a rotation memory aggregation module. The discrete curvature calculator reads the manifold feature matrix, calculates the angle defect value spanned by the difference between the preceding and following vectors for every three consecutive rows, and outputs the curvature sequence. The rotation memory aggregation module holds a learnable complex memory matrix, whose number of columns is equal to the dimension of the high-dimensional tangent vector sequence, and the number of rows is equal to the preset number of memory slots (set based on the time span of the high-dimensional tangent vector sequence and the fluctuation frequency of the curvature sequence), and a complex position encoding matrix is ​​injected row by row. The module first writes the high-dimensional tangent vector sequence into the memory matrix by slot, and then maps each scalar in the curvature sequence to a complex rotation factor through a parameterized rotation gate mechanism. The complex rotation controlled by the factor is applied to each row of the memory matrix slot by slot. All curvature values ​​are sequentially driven to rotate and update. After the memory matrix is ​​updated, the real part operation and column summation are performed, and then the power load prediction value is output through a linear projection layer.

[0075] The training process for the electricity load prediction model is as follows:

[0076] A training dataset containing historical time-series feature vectors and corresponding actual power load values ​​is constructed. The historical time-series feature vectors are used as model inputs, and the predicted power load values ​​are obtained by forward calculation through a Riemannian geometry embedding layer, a manifold tangent space transformation layer, and a curvature adaptive aggregation layer.

[0077] The mean square error between the predicted and actual power load values ​​is used as the loss function. The gradient of the loss with respect to the learnable parameters of each layer is calculated using the backpropagation algorithm. The gradient-based optimizer is then used to update the metric kernel matrix and scalar offset parameters of the Riemann geometry embedding layer, the tangent bundle basis dictionary and connection coefficient tensor of the manifold tangent space transformation layer, and the complex memory matrix and rotation gate parameters of the curvature adaptive aggregation layer. The forward calculation, loss calculation, gradient backpropagation, and parameter update steps are repeated until the loss function converges or the preset number of training rounds is reached.

[0078] The number of training rounds is set based on the size of the training data and the convergence trend of the loss function.

[0079] It should be noted that the tangent bundle basis dictionary refers to a set of multiple basis entries. Each entry contains a centroid vector, a set of left tangent basis vectors attached to the centroid, and a set of right tangent basis vectors. The centroid vector is used for local neighborhood assignment, and the left and right tangent basis vectors span the in-tangent space and out-tangent space at the centroid, respectively.

[0080] The parameterized rotation gate mechanism refers to transforming each scalar in the curvature sequence into the argument of a complex rotation factor through a learnable scaling factor and offset factor, and embedding this argument into a unit complex exponential form as a rotation gate, driving each row of the memory matrix to perform the unitary rotation corresponding to this argument on the complex plane.

[0081] It should also be noted that the learnable parameters of each metric kernel are a lower triangular matrix, with the number of rows equal to the embedding depth and the number of columns equal to the dimension of the temporal feature vector. A symmetric positive definite kernel matrix is ​​obtained by matrix multiplication of this lower triangular matrix and its transpose. The embedding depth is equal to the number of metric kernels, and the embedding depth is equal to the number of rows in the manifold feature matrix.

[0082] The initial values ​​of the centroid vectors are obtained by randomly sampling from the training data or by using k-means clustering. The dimension of each centroid vector is consistent with the dimension of the temporal feature vector. The dimensions of the left and right tangent basis vectors are both set to the embedding depth and generated using orthogonal random initialization. To ensure that the tangent basis vectors maintain orthogonality during training, a Gram-Schmidt reorthogonalization operation is performed after each parameter update.

[0083] The specific rules for propagation alignment are as follows: Given the local cut representation of the current row, the local cut representation of the next row, and the linear transformation matrix and translation vector between the corresponding basis entries, first perform matrix multiplication with the current local cut representation, then perform vector addition with the product result and the translation vector to obtain the transformed representation vector; then perform weighted fusion with the local cut representation of the next row, and the fusion result is used as the new current local cut representation to continue to propagate.

[0084] For the curvature scalar at each time step, it is first mapped to an intermediate value through a learnable linear transformation containing a weight matrix and a bias vector. This intermediate value is then compressed to the interval between zero and one using a sigmoid function and multiplied by twice the value of pi to obtain the rotation angle. Using this rotation angle as the argument, a complex rotation factor with a modulus of one is constructed. This complex rotation factor is then used to perform complex multiplication on each element of the complex elements in each row of the memory matrix.

[0085] For the symmetric positive definite kernel matrix in the metric kernel, after each gradient update, if non-positive eigenvalues ​​appear, eigenvalue decomposition is performed, and all eigenvalues ​​smaller than the smallest positive number are raised to that positive number before the matrix is ​​reconstructed. For the tangent basis vectors in the tangent bundle basis dictionary, an orthogonal regularization term can be added to the loss function, which is the sum of the squares of the inner products between all pairs of left tangent basis vectors, or Gram-Schmidt reorthogonalization can be performed after each parameter update.

[0086] S2.1 Input the time-series feature vector into the power load prediction model, and use the Riemannian geometric embedding layer to perform nonlinear spatial mapping on the time-series feature vector, and output the manifold feature matrix.

[0087] Furthermore, the Riemannian geometric embedding layer maintains a set of learnable metric kernels. Each metric kernel contains a symmetric positive definite kernel matrix obtained by multiplying a lower triangular matrix by its transpose, and a scalar offset parameter. For the input temporal feature vector, a quadratic form is calculated with each metric kernel matrix in sequence. The resulting scalar is added to the corresponding scalar offset parameter and then passed through a hyperbolic tangent activation function to obtain the activation scalar. The activation scalars generated by all metric kernels are concatenated into a row in sequence. This row is repeatedly generated according to the embedding depth and stacked vertically to obtain the manifold feature matrix.

[0088] S2.2 Utilize the manifold tangent space transformation layer to perform local linearization projection on the manifold feature matrix and output a high-dimensional tangent vector sequence.

[0089] Furthermore, within the manifold tangent space transformation layer, the input manifold feature matrix is ​​processed row-wise. For each row vector, its Euclidean distance to all centroid vectors in the tangent bundle basis dictionary is calculated. The basis entry with the smallest distance is selected as the nearest basis. The left tangent basis vector of this basis is then used to perform an inner product with the row vector to obtain the local tangent representation. Subsequently, in row order, the local tangent representations of each row are successively propagated and aligned using the affine transformation parameters indexed by the nearest basis entries in adjacent rows in the connection coefficient tensor. All propagated and aligned tangent representations are stacked vertically row-wise to form a high-dimensional tangent vector sequence.

[0090] S2.3. Using the curvature adaptive aggregation layer, the high-dimensional tangent vector sequence is modeled for spatiotemporal correlation and curvature correction, and the predicted power load value is output.

[0091] Furthermore, in the curvature adaptive aggregation layer, the high-dimensional tangent vector sequence is written into a complex memory matrix by slot. Then, a parameterized rotation gate maps each curvature value in the curvature sequence to a complex rotation factor. Complex rotations are applied to each row of the memory matrix slot by slot, and all curvature values ​​are sequentially updated cyclically. After the update is complete, the real part is taken and summed column by column to obtain a real vector. This real vector is then mapped to a scalar, i.e., the predicted power load value, through a linear projection layer.

[0092] S3. Collect historical prediction residuals, construct local moment state vectors based on historical prediction residuals, and use local moment state vectors to interact with neighboring distributed nodes to obtain neighborhood state information. Update the historical prediction residuals based on the neighborhood state information and output the calibration prediction residuals.

[0093] S3.1 Calculate the low-order raw moments and high-order central moments based on the historical prediction residual set, and construct the local moment state vector by dimensional concatenation of the low-order raw moments and high-order central moments.

[0094] Furthermore, the historical output values ​​of the power load forecasting model are subtracted from the actual measured power load values ​​at the corresponding times, and the results are summarized according to the time series to obtain the historical forecast residual set. The arithmetic mean of all residual samples in the historical forecast residual set is taken, and this result is used as the lower-order raw moment. The square of the difference between each residual sample and the lower-order raw moment is squared and then the arithmetic mean is calculated to obtain the variance. The cube of the difference between each residual sample and the lower-order raw moment is cubed and then the arithmetic mean is calculated, and then compared with the cube of the variance to obtain the skewness. The fourth power of the difference between each residual sample and the lower-order raw moment is cubed and then the arithmetic mean is calculated, and then compared with the square of the variance to obtain the kurtosis. The variance, skewness, and kurtosis are combined to form the higher-order central moment. The lower-order raw moment and the higher-order central moment are concatenated as vectors to obtain the local moment state vector.

[0095] S3.2. Set local dual variables based on the local moment state vector, synchronize the local moment state vector and local dual variables to the adjacent distributed nodes, and receive the returned neighborhood state information.

[0096] Furthermore, a local dual variable with the same dimension as the local moment state vector is established, and the local dual variable is initialized to a zero vector in the first iteration. In subsequent iterations, the local dual variable updated in the previous iteration is used.

[0097] Read the numbers of the adjacent distributed nodes directly connected to the current node in the preset distributed topology, encapsulate the current node number, the current iteration round, the local moment state vector, and the local dual variable into a state synchronization data packet, and send them to each of the adjacent distributed nodes respectively.

[0098] Listen to and receive the state synchronization data packets returned by each adjacent distributed node, verify the node number, iteration round, data dimension and data integrity in them, and remove data packets that fail the verification or have inconsistent iteration rounds. Write the adjacent node moment state vector and adjacent node dual variable in the data packets that pass the verification into the neighborhood cache according to the node number to obtain the neighborhood state information.

[0099] It should be noted that distributed topology is set based on the physical or logical connection relationship of distributed nodes and their communication network structure.

[0100] S3.3 Construct an augmented Lagrangian function based on the neighborhood state information, and use the augmented Lagrangian function to iteratively update the local moment state vector and the local dual variable to obtain the consensus moment vector.

[0101] Furthermore, the local moment state vector of the current node is used as the variable to be optimized, the moment state vectors of each neighboring distributed node in the neighborhood state information are used as the consensus reference, and the augmented Lagrangian function is constructed by combining the local dual variable and the dual variable of the neighboring node.

[0102] In each iteration, the neighborhood state information and the current dual variable are fixed, and the augmented Lagrangian function is minimized with respect to the local moment state vector to obtain the updated local moment state vector. Based on the deviation between the updated local moment state vector and the neighborhood moment state vector, the local dual variable is cumulatively updated. The updated local moment state vector and local dual variable are then synchronized to adjacent distributed nodes, and a new round of neighborhood state information is received. The process of updating the local moment state vector, updating the dual variable, and synchronizing the neighborhood is repeated. When the change in the local moment state vector is lower than a preset convergence threshold in two consecutive iterations, and the consistency deviation between the local moment state vector and the neighborhood moment state vector is lower than a preset consensus threshold, the iteration stops, and the local moment state vector obtained in the last iteration is used as the consensus moment vector.

[0103] The augmented Lagrangian function can be expressed as:

[0104] ;

[0105] in, It is the augmented Lagrangian function; It is a node The local moment state vector; It is a node The neighborhood moment state vector; It is a node For nodes The local dual variable; It is a secondary penalty parameter; Is with the current node A set of neighboring node numbers that have direct communication links; It is a Euclidean norm operator; It is the index of the currently processed node; It is an index of adjacent nodes; express and The inner product;

[0106] It should be noted that the quadratic penalty parameter is a positive real scalar, which can be set based on the trade-off between convergence speed and consensus accuracy during the iteration process.

[0107] The convergence threshold is set based on the change in the local moment state vector of the distributed nodes during the iterative update process. An example value range is... The basis for the range of values ​​in the example is that numerical calculations need to balance computational accuracy and convergence speed. Too small a threshold will cause the number of iterations to surge and cause oscillations, while too large a threshold cannot guarantee that the local moment state vector will reach a steady state. Selecting this range can ensure that the state vector stops iterating quickly while meeting the engineering accuracy requirements.

[0108] The consensus threshold is set based on the consistency deviation term between the local moment state vector and the neighborhood moment state vector among distributed nodes. An example value range is... The basis for the example value range is: since each node needs to achieve global state consistency in a distributed topology, the consensus threshold directly determines the quality of information interaction between nodes. Setting it within this range can take into account both the communication bandwidth limit and the strictness of global parameter synchronization, and avoid the failure of prediction residual calibration due to excessive differences between nodes.

[0109] The augmented Lagrangian function is a constrained optimization objective function with a quadratic penalty term and a dual variable. It is used to gradually achieve numerical consistency of the local moment state vectors of multiple distributed nodes through iteration of dual variables and penalty constraints, thereby obtaining the consensus moment vector.

[0110] S3.4. Using the consensus moment vector as the moment matching constraint for the historical prediction residual set, construct the maximum entropy optimization problem.

[0111] Furthermore, a probability weight to be solved is set for each residual sample in the historical prediction residual set, and all the probability weights to be solved are uniformly recorded as the set of optimization variables;

[0112] Following the order of mean, variance, skewness, and kurtosis in the consensus moment vector, moment constraint entries are generated sequentially by associating residual samples, unsolved probability weights, and corresponding moment types. Each moment constraint entry is used to ensure that the weighted residual statistical moments are consistent with the corresponding moment components in the consensus moment vector. For each unsolved probability weight, a non-negativity constraint is applied item by item, stipulating that the value of each unsolved probability weight is not less than zero. Based on this, a corresponding probability weight constraint entry is generated, which is based on the fact that probability weights, as probability values, must satisfy non-negativity by definition. At the same time, a normalization constraint is applied to all unsolved probability weights as a whole, stipulating that the sum of all unsolved probability weights is strictly equal to the unit total. Based on this, a normalization constraint entry is generated, which is based on the probability distribution requirement that the sum of all weights is one.

[0113] The entropy value of the probability weights to be solved is taken as the objective term, and the objective term is set as the direction of maximization. The set of optimization variables, the objective term, the moment constraint entries, the probability weight constraint entries, and the normalization constraint entries are written into the same optimization task structure to form a maximum entropy optimization problem.

[0114] S3.5. Use the Lagrange multiplier method to solve the maximum entropy optimization problem and obtain the optimal probability weights.

[0115] Furthermore, in the maximum entropy optimization problem, Lagrange multipliers are configured for the normalization constraints and each moment constraint item, and the maximization objective term, normalization constraint items, and moment constraint items are combined into a Lagrange function; stationary point conditions are solved for each probability weight to be solved, and the exponential mapping relationship between the probability weight and the corresponding residual sample, the value of each moment type, and the Lagrange multiplier is obtained; this exponential mapping relationship is substituted into the regression normalization constraints and each moment constraint item to form a nonlinear equation system containing only Lagrange multipliers.

[0116] The nonlinear equation system is numerically solved using Newton iteration, quasi-Newton iteration, or damped gradient iteration. In each iteration, the Lagrange multipliers are updated based on the constraint residuals until the normalized constraint residuals and the moment matching constraint residuals are all less than a preset solution threshold. The converged Lagrange multipliers are substituted back into the exponential mapping relationship to calculate the probability weight corresponding to each residual sample. The calculated probability weights are then checked for non-negativity (i.e., each weight is checked to see if its value is less than zero, and those weights are forced to be adjusted to zero). All probability weights after non-negativity check are summed. If the sum deviates from the unit total, a scaling factor is applied to each probability weight. This scaling factor is the ratio of the unit total to the current sum, ensuring that all probability weights are not less than zero and the total value remains at the unit total. The corrected probability weights of each residual sample are arranged in the original order of the residual samples to obtain the optimal probability weights.

[0117] It should be noted that the threshold value is set based on the normalized constraint residuals and moment matching constraint residuals generated during the solution of the maximum entropy optimization problem using the Lagrange multiplier method. The example value range is... The basis for the range of values ​​in the example is that, since the numerical solution of nonlinear equations requires a very high degree of satisfaction of the constraints, a smaller threshold can ensure that the optimal probability weight accurately reflects the statistical characteristics of the historical prediction residuals, and ensure that the calibration prediction residual set after resampling can strictly follow the constraints of the consensus moment vector, thereby improving the reliability of subsequent probability density reconstruction.

[0118] S3.6. Use the optimal probability weights to resample the historical prediction residual set based on importance, and generate the calibration prediction residual set.

[0119] Furthermore, following the original arrangement order of residual samples in the historical prediction residual set, each residual sample and its corresponding optimal probability weight are read, and each residual sample is bound to the optimal probability weight at the same position. The optimal probability weights are normalized and verified. If the total weight deviates from the unit total, the optimal probability weights are proportionally corrected according to the total weight to obtain a resampling probability distribution that can be used for sampling. The cumulative sampling probability of each residual sample is calculated based on the resampling probability distribution, and the sampling number is set to a preset resampling scale. A random number generator is used to generate sampling markers within the cumulative sampling probability range. Each sampling marker is matched with the cumulative sampling probability in intervals, and the residual sample corresponding to the interval is selected as a resampling result based on the interval matching result. This ensures that residual samples with larger optimal probability weights have a higher probability of being selected, and allows the same residual sample to be selected repeatedly. The sampling marker generation, interval matching, and residual sample selection process are repeated until the number of selected residual samples reaches the preset resampling scale. All selected residual samples are arranged according to the sampling generation order to generate a calibration prediction residual set.

[0120] It should be noted that the resampling size is set based on the sample size required for calibrating the prediction residual set and the limitations of computational resources.

[0121] S4. Calculate the skewness coefficient of the calibration prediction residual set, and based on the skewness coefficient, reconstruct the probability density of the calibration prediction residual set through dynamic asymmetric kernel density estimation, and output the dynamic quantile.

[0122] S4.1 Using the mean of the calibration prediction residual set as a benchmark, calculate the skewness coefficient of the calibration prediction residual set, and construct the AMISE objective function based on the skewness coefficient.

[0123] Furthermore, all residual samples in the calibration prediction residual set are read. First, the mean of the residual samples is calculated as the center position. Then, the deviation of each residual sample from the center position is calculated. Based on the deviation, the dispersion and third-order asymmetry of the residual samples are statistically analyzed. The third-order asymmetry is standardized according to the dispersion to obtain the skewness coefficient of the calibration prediction residual set.

[0124] The skewness coefficient is used to determine the skewness direction and intensity of the distribution of the calibration prediction residual set. When the skewness coefficient is greater than zero, it is determined that the right tail extension is greater than the left. When the skewness coefficient is less than zero, it is determined that the left tail extension is greater than the right. When the skewness coefficient is close to zero, it is determined that the left and right sides are approximately symmetrical.

[0125] The left and right bandwidths are denoted as the bandwidth variables to be optimized. Bandwidth adjustment rules are set for the left and right bandwidths respectively using skewness coefficients, so that the side with higher tail expansion corresponds to a larger smoothing tolerance range, and the side with lower tail expansion corresponds to a smaller smoothing tolerance range. Left smoothing bias and left estimation variance evaluation terms are set for the left bandwidth, and right smoothing bias and right estimation variance evaluation terms are set for the right bandwidth. The bandwidth adjustment rules corresponding to the skewness coefficients are written into the weight configurations of the left smoothing bias, right smoothing bias, left estimation variance, and right estimation variance evaluation terms, respectively. Finally, the left smoothing bias, right smoothing bias, left estimation variance, and right estimation variance evaluation terms are merged into an error evaluation objective, which is then specified as the minimization search object for the left and right bandwidths, resulting in the AMISE objective function, expressed as:

[0126] ;

[0127] ; ;

[0128] in, It is the AMISE objective function; It is the bandwidth variable on the left. It is the bandwidth variable on the right. It is the skewness coefficient for calibrating the predicted residual set; Indicates the weighting factor on the left; Indicates the weighting factor on the right; This represents the deviation coefficient of the left-hand smoothing deviation evaluation item; This represents the deviation coefficient of the right-hand smoothing deviation evaluation item; This represents the variance coefficient of the left-hand side estimate of the variance evaluation term; This represents the variance coefficient of the estimated variance evaluation term on the right side; This indicates the number of residual samples in the calibration prediction residual set; It is the skewness adjustment intensity coefficient; This represents the maximum value selection function, which selects the largest value from a set of values ​​given in parentheses as the output.

[0129] It should be noted that, It is set based on the probability density curvature characteristics of the calibration prediction residual set in the left residual value region and the strength of the left kernel smoothing bias;

[0130] It is set based on the probability density curvature characteristics of the calibration prediction residual set in the right residual value region and the strength of the right kernel smoothing bias;

[0131] It is set based on the sample dispersion of the calibration prediction residual set, the shape of the left kernel function, and the variance strength of the left probability density estimation;

[0132] It is set based on the sample dispersion of the calibration prediction residual set, the shape of the right-hand kernel function, and the strength of the variance of the right-hand probability density estimation;

[0133] The settings are based on the skewness sensitivity requirements of the calibration prediction residual set, the adjustment range of the left and right bandwidths, and the stability of the probability density estimation.

[0134] S4.2 Based on the statistical distribution characteristics of the calibration prediction residual set, the range of optimal left bandwidth and optimal right bandwidth is defined, and the range is discretized at equal intervals to form a two-dimensional search grid.

[0135] Furthermore, all residual samples in the calibration prediction residual set are read, and the number of residual samples, residual value range, residual dispersion, interquartile range, and skewness coefficient are statistically analyzed. The basic scale of the bandwidth is determined based on the residual value range and residual dispersion; the coarseness of the bandwidth search is determined based on the number of residual samples; and the relative expansion direction of the left and right bandwidth intervals is determined based on the skewness coefficient. When the skewness coefficient is greater than zero, the upper limit of the right bandwidth interval is set to be greater than the upper limit of the left bandwidth interval; when the skewness coefficient is less than zero, the upper limit of the left bandwidth interval is set to be greater than the upper limit of the right bandwidth interval; when the skewness coefficient is close to zero, the upper limit of the left bandwidth interval and... The right-side bandwidth range is set to be the same or approximately the same range; positive lower and upper limits are set for the left and right bandwidths respectively, so that the lower limit is used to avoid oscillations in probability density estimation caused by an overly narrow kernel function, and the upper limit is used to avoid excessive smoothing of distribution details caused by an overly wide kernel function; the left-side and right-side bandwidth ranges are discretely divided at equal intervals according to the preset number of divisions to obtain left-side candidate bandwidth sequences and right-side candidate bandwidth sequences; each candidate value in the left-side candidate bandwidth sequence is paired with each candidate value in the right-side candidate bandwidth sequence in turn, and all pairing results are arranged according to the left-side bandwidth dimension and the right-side bandwidth dimension to form a two-dimensional search grid.

[0136] It should be noted that the statistical distribution characteristics of the calibration prediction residual set refer to the number of residual samples, the range of residual values, the degree of residual dispersion, the interquartile range, and the skewness coefficient.

[0137] The number of segments is set based on a balance between the discretization accuracy requirements of the two-dimensional search grid and the consumption of computational resources.

[0138] S4.3 Substitute each bandwidth combination in the two-dimensional search grid into the AMISE objective function for numerical calculation to generate the bandwidth loss matrix.

[0139] Furthermore, all bandwidth combinations in the two-dimensional search grid are read, and the current bandwidth combination is selected one by one according to the grid arrangement order of the left and right bandwidth dimensions. For each current bandwidth combination, the left candidate bandwidth is written into the left bandwidth variable position of the AMISE objective function, and the right candidate bandwidth is written into the right bandwidth variable position of the AMISE objective function. The left smoothing bias evaluation term, right smoothing bias evaluation term, left estimation variance evaluation term, and right estimation variance evaluation term configured in the AMISE objective function are called for numerical calculation. The calculation results of each evaluation term are merged into the AMISE loss value corresponding to the current bandwidth combination, and the AMISE loss value is written into the matrix cell at the same position as the current bandwidth combination in the two-dimensional search grid.

[0140] Repeat the process of bandwidth combination selection, bandwidth variable writing, evaluation item calculation, loss value merging, and matrix cell writing until all bandwidth combinations in the two-dimensional search grid have been numerically calculated; store the AMISE loss values ​​in all matrix cells according to the arrangement relationship between the left and right bandwidth candidate sequences to generate a bandwidth loss matrix.

[0141] It should be noted that the bandwidth loss matrix is ​​a two-dimensional numerical array with left and right candidate bandwidths as its two dimensions and the AMISE objective function value corresponding to each bandwidth combination as its element. It is used to evaluate the estimation error of different bandwidth combinations in the two-dimensional search grid of asymmetric kernel density estimation, and to locate the optimal left and right bandwidth combinations that minimize the overall error.

[0142] S4.4. Search for the global minimum value within the bandwidth loss matrix, and extract the bandwidth combination corresponding to the global minimum value from the two-dimensional search grid to obtain the optimal left bandwidth and the optimal right bandwidth.

[0143] Furthermore, the bandwidth loss matrix and the two-dimensional search grid corresponding to its row and column positions are read, and each AMISE loss value is accessed one by one according to the arrangement order of the matrix units. The first AMISE loss value accessed is set as the current minimum loss value, and its row and column positions are recorded. The subsequent AMISE loss values ​​are compared with the current minimum loss value in turn. If the subsequent AMISE loss value is less than the current minimum loss value, the current minimum loss value is replaced with the AMISE loss value, and the corresponding row and column position records are updated synchronously. If there are multiple identical minimum loss values, one of the row and column positions is selected as the target position according to the preset selection rules.

[0144] After completing the comparison of all matrix units, the row and column positions of the final record are determined as the positions of the global minimum. The bandwidth combination corresponding to the row and column positions is read in the two-dimensional search grid. The left bandwidth in the bandwidth combination is determined as the optimal left bandwidth, and the right bandwidth in the bandwidth combination is determined as the optimal right bandwidth.

[0145] It should be noted that the selection rules are set based on the search preferences or smoothness priorities of the bandwidth loss matrix, such as prioritizing the combination with smaller bandwidth, prioritizing the combination with smaller difference between the left and right bandwidth, or prioritizing the combination that was retrieved first.

[0146] S4.5. Using the residual samples in the calibration prediction residual set as the center and the optimal left bandwidth and optimal right bandwidth as the supporting radii, construct an asymmetric Gaussian kernel function.

[0147] Furthermore, each residual sample in the calibration prediction residual set is read, and the value of the current residual sample is set as the center position of the corresponding kernel function. Using the current residual sample as the dividing point, the space of residual values ​​to be estimated is divided into a left region and a right region. Positions less than or equal to the current residual sample are assigned to the left region, and positions greater than the current residual sample are assigned to the right region. The optimal left bandwidth is configured as the smoothing scale and support radius of the left region, and the optimal right bandwidth is configured as the smoothing scale and support radius of the right region, so that the kernel function controls the decay rate according to the optimal left bandwidth to the left of the center position, and controls the decay rate according to the optimal right bandwidth to the right of the center position. A Gaussian kernel branch centered on the current residual sample and controlling the diffusion range with the optimal left bandwidth is set in the left region, and a Gaussian kernel branch centered on the current residual sample and controlling the diffusion range with the optimal right bandwidth is set in the right region, with the two Gaussian kernel branches continuously connected at the center position of the current residual sample. Scale correction and overall normalization are performed on the two Gaussian kernel branches to ensure that the total integral of the kernel function in the entire space of residual values ​​to be estimated is a unit total.

[0148] Repeatedly perform center position setting, left and right region division, bandwidth configuration, left and right Gaussian kernel branch setting, contiguous connection, and normalization processing on all residual samples in the calibration prediction residual set to obtain the asymmetric Gaussian kernel function centered on each residual sample, with the expression:

[0149] ;

[0150] in, It is an asymmetric Gaussian kernel function; It is a natural exponential function; It is the deviation between the values ​​of the sampling points to be estimated and the residual samples in the calibration prediction residual set; It is the optimal left-side bandwidth; It is the optimal right-side bandwidth; It is pi;

[0151] S4.6 Calculate the local probability contribution value of each residual sample using the asymmetric Gaussian kernel function, and perform arithmetic averaging and normalization on all local probability contribution values ​​to generate the residual probability density function.

[0152] Furthermore, density estimation sampling points are set within a preset residual value range with a fixed step size, and the asymmetric Gaussian kernel function centered on each residual sample is read. For each density estimation sampling point, the function response value of the sampling point on each asymmetric Gaussian kernel function is calculated sequentially, and the function response value is used as the local probability contribution value of the corresponding residual sample at that sampling point. The arithmetic mean of the local probability contribution values ​​generated by all residual samples at the same density estimation sampling point is used to obtain the initial probability density value at that sampling point.

[0153] Repeat the process of selecting sampling points, calculating kernel function response values, obtaining local probability contribution values, and arithmetic averaging until all density estimation sampling points have obtained initial probability density values. Perform non-negativity verification on all initial probability density values, and calculate the total integral of the initial probability density curve based on the fixed step size of the density estimation sampling points. Then, perform uniform proportional correction on the initial probability density values ​​at each sampling point according to the total integral, so that the total integral of the corrected probability density curve is a unit total over the entire residual value range. Arrange the corrected sampling point positions and their corresponding probability density values ​​in ascending order of residual values ​​to generate the residual probability density function.

[0154] It should be noted that the range of residual values ​​is set based on the boundary extreme values ​​of the calibration prediction residual set and the coverage requirements of the probability density estimation.

[0155] S4.7. Perform numerical integration on the residual probability density function to obtain the cumulative distribution function, and solve for the upper and lower tail positions of the cumulative distribution function to obtain the dynamic quantiles.

[0156] Furthermore, the sampling points and their corresponding probability density values ​​in the residual probability density function are read in ascending order of residual value. Starting from the smallest residual sampling point, numerical integration is performed sequentially along the direction of residual value. The probability density surface before each sampling point is accumulated to obtain the cumulative probability value corresponding to that sampling point. Monotonicity verification and endpoint normalization correction are performed on the cumulative probability value to ensure that the cumulative probability remains non-decreasing as the residual value increases, and the cumulative probability corresponding to the largest residual sampling point reaches the unit total, thereby obtaining the cumulative distribution function.

[0157] The lower tail probability target and upper tail probability target are determined based on the preset interval confidence level. The sampling point that first reaches the lower tail probability target is retrieved from the cumulative distribution function in ascending order, and the residual value corresponding to this sampling point is determined as the lower quantile. Similarly, the sampling point that first reaches the upper tail probability target is retrieved from the cumulative distribution function in ascending order, and the residual value corresponding to this sampling point is determined as the upper quantile. If the probability target falls between two adjacent sampling points, linear interpolation is performed based on the residual values ​​of the adjacent sampling points and the cumulative probability value to obtain the quantiles under the corresponding tail probability target. The obtained lower and upper quantiles are used together as the dynamic quantiles.

[0158] It should be noted that the interval confidence level is set based on the target coverage probability of the power load forecast interval, the business risk tolerance, and the stability of the historical residual distribution.

[0159] S5. Through conformal prediction, the dynamic quantile and the power load forecast value are linearly combined to generate the power load forecast interval.

[0160] S5.1. Linearly sum the predicted power load value with the lower quantile and the upper quantile respectively to form the initial prediction interval.

[0161] Furthermore, the predicted power load value is linearly summed with the lower quantile to obtain the lower boundary of the initial prediction interval; the predicted power load value is linearly summed with the upper quantile to obtain the upper boundary of the initial prediction interval; if the lower boundary is greater than the upper boundary, the positions of the lower boundary and the upper boundary are swapped to ensure that the interval boundary order is valid; the lower boundary and the upper boundary with the valid order are combined to form the initial prediction interval.

[0162] S5.2 Construct residual interval boundaries using lower and upper quantiles, and for each residual sample in the calibration prediction residual set, calculate the maximum non-negative deviation beyond the residual interval boundaries as the inconsistency score.

[0163] Furthermore, the lower quantile is set as the lower boundary of the residual interval, and the upper quantile is set as the upper boundary of the residual interval. If the lower quantile is greater than the upper quantile, their positions are swapped to ensure the validity of the residual interval boundary order. The residual samples in the calibration prediction residual set are read one by one, and the current residual sample is compared with the lower and upper boundaries of the residual interval. When the current residual sample is less than the lower boundary of the residual interval, the difference between it and the lower boundary is taken as the lower deviation, and the upper deviation is set to zero. When the current residual sample is greater than the upper boundary of the residual interval, the difference between it and the upper boundary is taken as the upper deviation, and the lower deviation is set to zero. When the current residual sample is between the lower and upper boundaries of the residual interval, both the lower and upper deviations are set to zero. The larger non-negative value between the lower and upper deviations is selected as the inconsistency score corresponding to the current residual sample.

[0164] Repeat the boundary comparison, deviation calculation and selection of larger non-negative values ​​process for all residual samples in the calibration prediction residual set to obtain the non-consistency score corresponding to each residual sample.

[0165] S5.3 Arrange all non-consistent fractions into an ordered fraction sequence in ascending order, and use the non-consistent fractions at the quantile rank in the ordered fraction sequence as conformal corrections.

[0166] Furthermore, the non-consistency scores corresponding to each residual sample in the calibration prediction residual set are read, and the validity of the non-consistency scores is verified. Null values ​​or non-numerical results are removed, and all valid non-consistency scores are retained. The valid non-consistency scores are sorted in ascending order of numerical value to form an ordered score sequence. The quantile rank is determined based on the interval confidence level and the number of valid non-consistency scores, and the quantile rank is limited to the valid position range of the ordered score sequence. When the quantile rank corresponds to an integer position, the non-consistency score at that position in the ordered score sequence is directly read. When the quantile rank falls between two adjacent positions, a sequence position not lower than the quantile rank is selected according to the preset rounding rule, and the non-consistency score at that position is read to ensure that the correction amount covers the interval confidence level requirement. The read non-consistency scores are determined as the conformal correction amount.

[0167] It should be noted that the rounding rule is set based on the discrete properties of ordered fractional sequences.

[0168] S5.4. Linearly compensate the conformal correction amount to the lower and upper boundaries of the initial prediction interval to generate the power load prediction interval.

[0169] Furthermore, the conformal correction is checked for non-negativity; if the conformal correction is less than zero, it is corrected to zero. The conformal correction is used as the interval expansion compensation amount. Negative linear compensation is performed on the lower boundary of the initial prediction interval to ensure that the compensated lower boundary is not higher than the lower boundary of the initial prediction interval. Positive linear compensation is performed on the upper boundary of the initial prediction interval to ensure that the compensated upper boundary is not lower than the upper boundary of the initial prediction interval. The order of the compensated lower boundary and compensated upper boundary is checked. If the compensated lower boundary is greater than the compensated upper boundary, their positions are swapped to ensure that the interval boundary order is valid. The compensated lower boundary and compensated upper boundary with valid order are combined to generate the power load prediction interval.

[0170] In summary, this invention improves point prediction accuracy by deeply mining the nonlinear spatiotemporal characteristics of load data through Riemannian geometric embedding and curvature correction techniques, and performs real-time calibration of prediction residuals by combining moment consensus collaboration and maximum entropy optimization under the distributed control concept. Finally, it reconstructs the residual distribution using dynamic asymmetric kernel density estimation and conformal prediction techniques, thereby generating accurate prediction intervals with rigorous statistical coverage guarantees even when load distribution drifts, effectively supporting reliable scheduling and decision-making in power grid operation.

[0171] 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. An adaptive interval prediction method for power load distribution drift, characterized in that, include: Collect multi-source heterogeneous data and preprocess the multi-source heterogeneous data to generate time-series feature vectors; The time-series feature vector is input into the pre-trained power load prediction model for forward calculation, and the power load prediction value is output. Collect historical prediction residuals, construct local moment state vectors based on historical prediction residuals, and use local moment state vectors to interact with neighboring distributed nodes to obtain neighborhood state information. Update the historical prediction residuals based on the neighborhood state information and output the calibration prediction residuals. Calculate the skewness coefficient of the calibration prediction residual set, and reconstruct the probability density of the calibration prediction residual set based on the skewness coefficient through dynamic asymmetric kernel density estimation, and output the dynamic quantile; By using conformal prediction, dynamic quantiles and power load forecast values ​​are linearly combined to generate power load forecast intervals.

2. The adaptive interval prediction method for power load distribution drift as described in claim 1, characterized in that, The process of collecting multi-source heterogeneous data and preprocessing it to generate time-series feature vectors involves the following steps: Dimensionless processing is performed on multi-source heterogeneous data to output normalized feature data; Interpolation and time-scale resampling are performed on the normalized feature data to output spatiotemporally aligned data; Sliding window sampling is performed on spatiotemporally aligned data to output temporal feature vectors.

3. The adaptive interval prediction method for power load distribution drift as described in claim 2, characterized in that, The power load prediction model is constructed based on a Riemannian geometry embedding layer, a manifold tangent space transformation layer, and a curvature adaptive aggregation layer.

4. The adaptive interval prediction method for power load distribution drift as described in claim 3, characterized in that, The steps for inputting the time-series feature vector into the pre-trained power load prediction model for forward calculation and outputting the power load prediction value are as follows: The time-series feature vector is input into the power load prediction model, and a Riemannian geometric embedding layer is used to perform nonlinear spatial mapping on the time-series feature vector to output the manifold feature matrix. By utilizing the manifold tangent space transformation layer, the manifold feature matrix is ​​locally linearized and projected to output a high-dimensional tangent vector sequence; By utilizing a curvature adaptive aggregation layer, spatiotemporal correlation modeling and curvature correction are performed on high-dimensional tangent vector sequences to output predicted power load values.

5. The adaptive interval prediction method for power load distribution drift as described in claim 1, characterized in that, The process involves collecting historical prediction residuals, constructing local moment state vectors based on these residuals, and using these local moment state vectors to interact with neighboring distributed nodes to obtain neighborhood state information. The specific steps are as follows: The lower-order raw moments and higher-order central moments are calculated based on the historical prediction residual set, and the lower-order raw moments and higher-order central moments are constructed into local moment state vectors through dimension concatenation. Local dual variables are set based on the local moment state vector, and the local moment state vector and local dual variables are synchronized to the adjacent distributed nodes, and the neighborhood state information is received back.

6. The adaptive interval prediction method for power load distribution drift as described in claim 1, characterized in that, The steps for updating the historical prediction residual set based on neighborhood state information and outputting the calibration prediction residual set are as follows: An augmented Lagrangian function is constructed based on neighborhood state information, and the local moment state vector and local dual variable are iteratively updated using the augmented Lagrangian function to obtain the consensus moment vector. Using the consensus moment vector as the moment matching constraint for the historical prediction residual set, a maximum entropy optimization problem is constructed. The Lagrange multiplier method is used to solve the maximum entropy optimization problem and obtain the optimal probability weights. The historical prediction residual set is resampled based on the optimal probability weights to generate a calibration prediction residual set.

7. The adaptive interval prediction method for power load distribution drift as described in claim 6, characterized in that, The calculation of the skewness coefficient of the calibration prediction residual set, and the probability density reconstruction of the calibration prediction residual set based on the skewness coefficient through dynamic asymmetric kernel density estimation, outputting dynamic quantiles, are detailed in the following steps: Using the mean of the calibration prediction residual set as a benchmark, the skewness coefficient of the calibration prediction residual set is calculated, and the AMISE objective function is constructed based on the skewness coefficient. By searching for the global minimum, the AMISE objective function is jointly optimized, and the optimal left bandwidth and optimal right bandwidth are output. Using the optimal left bandwidth and optimal right bandwidth as scaling parameters, an asymmetric Gaussian kernel function is selected to reconstruct the probability density of the calibration prediction residual set, generating the residual probability density function. Numerical integration is performed on the residual probability density function to obtain the cumulative distribution function, and the upper and lower tail quantiles are solved for the cumulative distribution function to obtain the dynamic quantiles; the dynamic quantiles include the upper quantile and the lower quantile.

8. The adaptive interval prediction method for power load distribution drift as described in claim 7, characterized in that, The steps for jointly optimizing the AMISE objective function by retrieving the global minimum and outputting the optimal left and right bandwidths are as follows: Based on the statistical distribution characteristics of the calibration prediction residual set, the range of optimal left bandwidth and optimal right bandwidth is defined, and the range is discretized at equal intervals to form a two-dimensional search grid. The bandwidth combination in the two-dimensional search grid is substituted into the AMISE objective function for numerical calculation to generate the bandwidth loss matrix. The global minimum is retrieved within the bandwidth loss matrix, and the bandwidth combination corresponding to the global minimum is extracted from the two-dimensional search grid to obtain the optimal left bandwidth and the optimal right bandwidth.

9. The adaptive interval prediction method for power load distribution drift as described in claim 7, characterized in that, The steps involve using the optimal left and right bandwidths as scaling parameters, employing an asymmetric Gaussian kernel function to reconstruct the probability density of the calibration prediction residual set, and generating a residual probability density function. The specific steps are as follows: An asymmetric Gaussian kernel function is constructed with the residual samples in the calibration prediction residual set as the center and the optimal left bandwidth and optimal right bandwidth as the supporting radii. The local probability contribution value of each residual sample is calculated using an asymmetric Gaussian kernel function, and the arithmetic mean and normalization of all local probability contribution values ​​are performed to generate the residual probability density function.

10. The adaptive interval prediction method for power load distribution drift as described in claim 1, characterized in that, The process of generating an electricity load forecast interval by linearly combining dynamic quantiles and electricity load forecast values ​​through conformal prediction is as follows: The electricity load forecast is linearly summed with the lower quantile and the upper quantile respectively to form the initial forecast interval; The lower and upper quantiles are used to construct the boundary of the residual interval, and for each residual sample in the calibration prediction residual set, the maximum non-negative deviation beyond the boundary of the residual interval is calculated as the inconsistency score. Arrange all non-consistent scores in ascending order to form an ordered score sequence, and use the non-consistent scores at the quantile rank in the ordered score sequence as the conformal correction amount; The conformal correction is linearly compensated to the lower and upper boundaries of the initial prediction interval to generate the power load prediction interval.

Images