Power sale quantity risk prediction method, device and system based on causal interpretation and online self-learning

Abstract

The invention provides an electricity sales risk prediction method, device and system based on causal interpretation and online self-learning, and the method comprises the steps: obtaining historical electricity sales records and associated power grid operation observation data, building a sample data set containing a time stamp, carrying out the causal recognition of the sample data set, and carrying out the prediction of the electricity sales risk. Analyzing direct causal association and indirect causal association among variables through a structural causal model, generating a causal path weight set, performing Bayesian network fusion on the causal path weight set and time sequence pattern features in a sample data set, constructing a time sequence representation vector with causal logic constraints, and inputting the time sequence representation vector into a meta-learning framework model to obtain a meta-learning framework model; and generating an initial risk predictor through parameter fine tuning, obtaining a deviation value between the real-time power sale quantity observation value and the probability distribution sequence so as to adjust parameter configuration of the initial risk predictor, updating a causal intensity quantized value in the causal path weight set, and generating a power sale quantity risk assessment result. The accuracy and long-term stability of the risk prediction result can be improved.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a method, apparatus and system for predicting electricity sales risk based on causal explanation and online self-learning. Background Technology

[0002] With the continuous development of the electricity market and the increasing complexity of power grid operation, analyzing historical data and related influencing factors can assess the likelihood of abnormal fluctuations in electricity sales over a future period, providing decision support for power grid dispatch, resource allocation, and risk prevention. Currently, most electricity sales risk prediction methods are based on historical electricity sales data and related influencing factors (such as economic indicators and meteorological data). Predictive relationships are constructed using statistical or machine learning models. These methods typically rely on the statistical correlations between data for modeling, training model parameters using historical data to predict future electricity sales. However, existing methods often face the problem of imprecise variable relationship modeling in practical applications, making them susceptible to spurious correlations that lead to prediction biases. Furthermore, when actual operating conditions change, the prediction models are prone to insufficient adaptability, resulting in decreased risk prediction accuracy and impacting the accuracy of power grid operation decisions. Summary of the Invention

[0003] In view of this, the present invention provides a method, apparatus, and system for predicting electricity sales risk based on causal explanation and online self-learning. The technical solution of the embodiments of the present invention is implemented as follows:

[0004] On one hand, embodiments of the present invention provide a method for predicting electricity sales risk based on causal explanation and online self-learning. The method includes: acquiring historical electricity sales records and associated power grid operation observation data, establishing a sample dataset containing time stamps, wherein the sample dataset has temporal continuity and variable correlation; identifying causal relationships in the sample dataset, analyzing direct and indirect causal relationships between variables through a structural causal model, and generating a causal path weight set containing causal strength quantification values; fusing the causal path weight set with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logic constraints; inputting the time series representation vector into a meta-learning framework model, generating an initial risk predictor through parameter fine-tuning, wherein the initial risk predictor outputs a probability distribution sequence of electricity sales fluctuations; acquiring the deviation value between the real-time electricity sales observation value and the probability distribution sequence, adjusting the parameter configuration of the initial risk predictor based on the deviation value, and updating the causal strength quantification values ​​in the causal path weight set to generate an electricity sales risk assessment result.

[0005] On the other hand, embodiments of the present invention provide a power sales risk prediction device, comprising: a data acquisition module, used to acquire historical power sales records and associated power grid operation observation data, and establish a sample dataset containing time stamps, wherein the sample dataset has temporal continuity and variable correlation; a causal identification module, used to identify causal relationships in the sample dataset, analyze direct and indirect causal relationships between variables through a structural causal model, and generate a causal path weight set containing causal strength quantification values; a vector construction module, used to fuse the causal path weight set with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logic constraints; a parameter fine-tuning module, used to input the time series representation vector into a meta-learning framework model, generate an initial risk predictor through parameter fine-tuning, wherein the initial risk predictor outputs a probability distribution sequence of power sales fluctuations; and a risk assessment module, used to acquire the deviation value between the real-time power sales observation value and the probability distribution sequence, adjust the parameter configuration of the initial risk predictor based on the deviation value, update the causal strength quantification values ​​in the causal path weight set, and generate a power sales risk assessment result.

[0006] In another aspect, embodiments of the present invention provide a computer system including a memory and a processor, wherein the memory stores a computer program that can run on the processor, and the processor executes the program to implement the steps in the methods described above.

[0007] This invention establishes a sample dataset containing temporal continuity and variable correlation, analyzes direct and indirect causal relationships between variables using a structural causal model, and generates a set of causal path weights. This set is then fused with time series pattern features through a Bayesian network to construct a time series representation vector with causal logic constraints. This vector is input into a meta-learning framework model to generate an initial risk predictor. Based on online self-learning, the parameter configuration and causal path weights are dynamically adjusted. This method introduces causal logic into the electricity sales risk prediction model by constructing a set of causal path weights, effectively distinguishing between direct and indirect influences between variables. It avoids the spurious correlation problem caused by traditional data-driven models relying solely on statistical correlations, improving the inherent logical reliability of risk prediction. Simultaneously, online self-learning enables real-time monitoring of prediction deviations and collaborative updates of model parameters and causal path weights, allowing the prediction model to dynamically adapt to new data pattern changes and avoiding the decline in prediction accuracy caused by model solidification. Ultimately, through the organic combination of causal explanation and online self-learning, this method achieves a technological leap in electricity sales risk prediction from simple data fitting to causal-driven and dynamically evolving approaches, improving the accuracy and long-term stability of risk prediction results. Attached Figure Description

[0008] Figure 1This is a schematic diagram illustrating the implementation process of a method for predicting electricity sales risk based on causal explanation and online self-learning, provided in an embodiment of the present invention.

[0009] Figure 2 This is a schematic diagram of the composition of a power sales risk prediction device provided in an embodiment of the present invention.

[0010] Figure 3 This is a schematic diagram of the hardware entity of a computer system provided in an embodiment of the present invention. Detailed Implementation

[0011] This invention provides a method for predicting electricity sales risk based on causal explanation and online self-learning. This method can be executed by a processor of a computer system. The computer system can refer to a device with data processing capabilities, such as a server, laptop, tablet, or desktop computer.

[0012] Figure 1 This is a schematic diagram illustrating the implementation process of a method for predicting electricity sales risk based on causal explanation and online self-learning, as provided in an embodiment of the present invention. Figure 1 As shown, the method includes: Step S100: Obtain historical electricity sales records and associated power grid operation observation data, and establish a sample dataset containing time stamps. The sample dataset has time continuity and variable correlation.

[0013] Historical electricity sales records are data on the amount of electricity sold over a past period, reflecting historical electricity sales and demonstrating trends in electricity sales over different time periods. Grid operation observation data is data obtained by observing various parameters and states during grid operation, covering information such as voltage, current, power, and frequency. This data reflects the grid's operating status and performance. Time stamps are added to the data to identify the data collection time, giving the data a temporal order and a temporal dimension. The sample dataset consists of historical electricity sales records and associated grid operation observation data, including time stamps, and exhibits temporal continuity and variable correlation. Temporal continuity indicates that the data is continuous in time, without significant time intervals or gaps; variable correlation indicates that there is a relationship between different variables in the dataset, such as electricity sales being related to grid load, voltage, and other factors.

[0014] Step S200: Identify causal relationships in the sample dataset, analyze the direct and indirect causal relationships between variables through structural causal modeling, and generate a set of causal path weights containing quantified causal strength values.

[0015] Causal relationship identification involves finding causal relationships between different variables in a sample dataset, determining which variables are causal variables and which are outcome variables. Structural causal models are models used to analyze causal relationships. They represent causal relationships between variables by constructing causal graphs, clearly showing direct and indirect causal associations. A direct causal association is a direct connection between variables without intermediate nodes; a change in one variable directly causes a change in another. An indirect causal association is a transitive relationship between variables through at least one intermediate node; a change in one variable indirectly causes a change in another through the transit of that node. Causal strength quantification is a quantitative representation of the strength of a causal relationship, reflecting the degree of influence of the causal variable on the outcome variable. A causal path weight set is a set of causal paths containing causal strength quantification values. These causal paths are dynamically adjusted over time, reflecting the causal relationships between variables over different time periods.

[0016] In one implementation, step S200 may specifically include the following steps S210 to S260: Step S210: Perform multi-scale time series decomposition on the sample dataset, and separate the intraday fluctuation component, intraweek trend component and intramonth cycle component according to the nested relationship of time period to obtain a multi-dimensional time series set with hierarchical time structure, and each component maintains the time label correspondence of the original data.

[0017] Multi-scale time series decomposition involves breaking down time series data according to different time scales, dividing it into components of varying frequencies to better analyze and understand the characteristics of the time series. Nested time periods refer to the nested relationships between different time periods; for example, intraday volatility is nested within weekly trend components, and weekly trend components are nested within monthly periodic components. Intraday volatility reflects the fluctuations of the time series within a single day, indicating changes in different time periods. Weekly trend reflects the trend of the time series within a week, indicating changes in different days within a week. Monthly periodic components reflect the periodic changes of the time series within a month, indicating periodic characteristics in different time periods within a month. A hierarchical multi-dimensional time series set consists of intraday volatility, weekly trend, and monthly periodic components, forming a multi-dimensional time series structure. Maintaining the original time stamp correspondence for each component ensures that each component retains the original time stamp during the decomposition process, making the components temporally corresponding.

[0018] In one implementation, step S210 may specifically include the following steps S211 to S216: Step S211: Adaptively decompose the historical electricity sales records in the sample dataset using the ensemble empirical mode decomposition algorithm to generate multiple intrinsic mode function components. Each component is arranged from high to low frequency to represent the fluctuation characteristics at different time scales.

[0019] Ensemble Empirical Mode Decomposition (EEMD) is an improved algorithm based on Empirical Mode Decomposition (EMD). By adding white noise to the original signal, it overcomes the mode aliasing problem of EMD and can more accurately decompose a time series into multiple Intrinsic Mode Function (IMF) components. These IMF components are obtained through EEMD, each possessing specific frequency and fluctuation characteristics, representing the fluctuations of the time series at different time scales. Arranging the components by frequency from high to low involves sorting the IMF components according to their frequency, with higher-frequency components appearing first and lower-frequency components appearing last.

[0020] For example, white noise of a certain intensity can be added to the original historical electricity sales record signal to obtain a noisy signal. Then, empirical mode decomposition (EMD) is performed on the noisy signal, decomposing it into multiple intrinsic mode function (IMF) components and a residual component. This process is repeated multiple times, adding different white noise each time, to obtain multiple sets of decomposition results. Finally, the multiple sets of decomposition results are averaged to obtain the final IMF components. For instance, for a time series containing one year of historical electricity sales records, an ensemble EMD algorithm can decompose it into multiple IMF components. The high-frequency IMF components may reflect the fluctuations in electricity sales within a day, while the low-frequency IMF components may reflect the trend changes in electricity sales over a month or a year.

[0021] Step S212: Identify the periodic characteristics of the intrinsic mode function components, calculate the principal period length of each component through the autocorrelation function, and classify the components with principal period lengths in different ranges into intraday fluctuation components, intraweek trend components, and intramonth periodic components.

[0022] For example, the autocorrelation function of each intrinsic mode function component can be calculated first. The autocorrelation function can be a general autocorrelation function (ACF). Then, the peak points of the autocorrelation function are found, and the time delays corresponding to these peak points are the periods of the components. Finally, based on the length of the main period, the components are classified into intraday fluctuation components, weekly trend components, and monthly period components.

[0023] Step S213: Calculate the energy percentage of each periodic component. The energy percentage is the ratio of the component variance to the original sequence variance. Components with an energy percentage exceeding a preset ratio are retained, and noise components with an energy percentage below a preset ratio threshold are removed.

[0024] Energy percentage is the proportion of energy of each periodic component in the total energy of the original sequence, reflecting the contribution of that component to the original sequence. Component variance is the variance of the periodic component, measuring its volatility. Original sequence variance is the variance of the original time series, reflecting the overall volatility of the original sequence. Preset proportion is a pre-defined threshold used to determine whether the periodic component has sufficient energy. Noise components are those with an energy percentage below the preset threshold; these components are usually caused by noise or interference and have a relatively small impact on the original sequence, so they can be removed.

[0025] Step S214: Align the intraday fluctuation component, weekly trend component, and monthly cycle component with time stamps so that each component has a corresponding observation value at the same time node, and construct a multi-dimensional time series set containing three time scales.

[0026] Specifically, a unified time range and time interval are determined as the benchmark for time stamping. Then, for each periodic component, its time stamp is adjusted to match the benchmark time stamp. If a component lacks an observation at a certain time point, interpolation methods such as linear interpolation or spline interpolation can be used to supplement it. Finally, the aligned intraday volatility component, intraweekly trend component, and intramonthly periodic component are combined to construct a multi-dimensional time series set containing three time scales.

[0027] Step S215: Analyze the interaction between different periodic components, calculate the modulation coefficient of intraday fluctuation component on intraweek trend component. The modulation coefficient characterizes the influence of short-term fluctuation on medium-term trend. The modulation coefficient is incorporated as an additional feature into the multi-dimensional time series set.

[0028] The interaction between different periodic components refers to the interaction and influence between intraday volatility components, weekly trend components, and monthly periodic components. For example, changes in intraday volatility components may affect the trend of weekly trend components, and changes in weekly trend components may also affect monthly periodic components. The modulation coefficient is an indicator used to measure the strength of the influence of intraday volatility components on weekly trend components, reflecting the modulating effect of short-term fluctuations on medium-term trends. Incorporating the modulation coefficient as an additional feature into a multi-dimensional time series dataset involves adding the calculated modulation coefficient as a new feature to the dataset, allowing the impact of short-term fluctuations on medium-term trends to be considered in subsequent analysis and forecasting.

[0029] Specifically, the intraday volatility component and the weekly trend component are first normalized to ensure they have the same scale and range. Then, regression analysis is used to establish a regression model between the intraday volatility component and the weekly trend component. For example, a linear regression model y = β0 + β1x + γ can be used, where y is the weekly trend component, x is the intraday volatility component, β0 and β1 are regression coefficients, and γ is the error term. The regression coefficient β1 is estimated using the least squares method; this coefficient is the modulation coefficient of the intraday volatility component on the weekly trend component. Finally, the modulation coefficient is added as an additional feature to the multi-dimensional time series dataset.

[0030] Step S216: Map each periodic component to a unified time axis using a time-scale reconstruction algorithm, so that each component in the multi-dimensional time series set has the same time resolution.

[0031] For example, a uniform time resolution is first determined as the target time resolution. Then, for each periodic component, time scaling is performed using interpolation or sampling methods based on its original time resolution and the target time resolution. If the original time resolution is higher than the target time resolution, sampling can be used to select an observation at regular time intervals. If the original time resolution is lower than the target time resolution, interpolation methods, such as linear interpolation or spline interpolation, can be used to estimate observations at missing time points. Finally, the transformed periodic components are mapped onto a uniform time axis to give them the same time resolution.

[0032] Step S220: Perform time-varying analysis of variable associations on the multi-dimensional time series set, calculate the dynamic association degree between variables within the sliding time window, and generate a time-varying curve of association strength as the time window moves. The time-varying curve of association strength represents the changing trend of variable associations in different time periods.

[0033] Because the relationships between variables can change over time, time-varying analysis is needed to capture these changes. A sliding window is a method for local analysis on a time series, analyzing the data within a fixed-length window by sliding it across the time series. Dynamic correlation strength is the degree of association between variables calculated within the sliding window, reflecting the correlation between variables within the current time window. A time-varying correlation strength curve is formed by calculating dynamic correlation strength within different time windows and connecting these correlation strength values ​​in chronological order; it visually demonstrates the changing trend of the correlation strength between variables over time.

[0034] For example, the length of the sliding time window and the sliding step size can be determined first. The length of the sliding time window determines the data range for each analysis, and the sliding step size determines the speed at which the window moves across the time series. Then, the window is slid across the time series, and for the data within each window, the dynamic correlation between variables is calculated. Methods such as correlation coefficients and mutual information, for example, the Pearson correlation coefficient, can be used to calculate the dynamic correlation. Finally, the dynamic correlations calculated for each window are connected in chronological order to generate a time-varying curve of the correlation strength.

[0035] Step S230: Based on the time-varying curve of correlation strength, perform hierarchical division of causal paths, and distinguish between direct causal paths and indirect causal paths through path length analysis. Direct causal paths are direct connections between variables without intermediate nodes, while indirect causal paths contain transitive connections with at least one intermediate node.

[0036] The hierarchical division of causal paths is based on the causal relationships between variables, dividing causal paths into different levels to more clearly demonstrate the structure of causal relationships. Path length analysis is a method used to distinguish between direct and indirect causal paths, determining the path type by calculating the number of edges contained in the causal path. A direct causal path is a direct connection between variables without intermediate nodes; that is, a change in one variable directly causes a change in another variable. An indirect causal path is a transitive relationship between variables through at least one intermediate node; that is, a change in one variable indirectly causes a change in another variable through the transit of an intermediate node.

[0037] In one implementation, step S230 may specifically include the following steps S231 to S236: Step S231: Construct the adjacency matrix of the variable association network. The matrix elements represent the maximum value of the time-varying curve of the association strength of the corresponding variable pair. The larger the value, the stronger the association between the variables.

[0038] A variable association network is a network used to represent the relationships between variables. It consists of nodes and edges, where nodes represent variables and edges represent the relationships between variables. An adjacency matrix is ​​a matrix used to represent a graph. For a variable association network, the rows and columns of the adjacency matrix correspond to the nodes (i.e., variables) in the network, and the matrix elements represent the association strength between the corresponding nodes. The maximum value of the association strength time-varying curve is the maximum association strength between a pair of variables. This maximum value reflects the strongest association between the pair of variables over the entire time span. A larger value indicates a stronger association between the variables, meaning a greater mutual influence between the two variables.

[0039] For example, first, based on the time-varying curve of association strength, find the maximum value of the association strength for each variable pair. Assuming there are n variables in the variable association network, we need to find the maximum value of the association strength for n×n variable pairs. Then, fill these maximum values ​​into the adjacency matrix, where the element in the i-th row and j-th column corresponds to the maximum value of the association strength between variable i and variable j.

[0040] Step S232: Traverse the adjacency matrix, starting from each variable node to explore other reachable variable nodes, and record all possible variable connection paths and their corresponding path lengths. The path length is the number of edges contained in the path.

[0041] Specifically, we can start from the first row of the adjacency matrix and visit each row sequentially. For each row's corresponding variable node, we use it as the starting point and use a graph search algorithm, such as Depth-First Search (DFS) or Breadth-First Search (BFS), to explore other reachable variable nodes. During the search, we record the paths traversed and their lengths. For example, using the DFS algorithm, starting from variable node A, if the elements of A and B in the adjacency matrix are not 0, it means there is a connection between A and B. We record the path AB with a length of 1. Then, starting from B, we continue the search. If there is a connection between B and C, we record the path ABC with a length of 2. We repeat the above process until all rows of the adjacency matrix have been traversed, recording all possible variable connection paths and their corresponding lengths. For example, for the variable association network containing three variables A, B, and C, starting from A, possible paths include AB (path length 1), AC (path length 1), ABC (path length 2), etc.

[0042] Step S233: Obtain a preset path length threshold, mark the connection paths with a path length equal to the path length threshold as direct causal path candidate set, and mark the path paths with a path length greater than the path length threshold as indirect causal path candidate set, thus initially distinguishing between the two types of paths.

[0043] The preset path length threshold is a pre-defined standard used to distinguish between direct and indirect causal paths. Paths with a path length equal to the threshold are considered candidates for direct causal paths because their paths are relatively short, potentially indicating that the variables are not directly connected through many intermediate nodes. Paths with a path length greater than the threshold are considered candidates for indirect causal paths because their paths are relatively long, potentially indicating that the variables are correlated through the passage of intermediate nodes.

[0044] Step S234: Perform time sequence verification on the candidate set of direct causal paths, determine the causal direction between variables through Granger causality test, and retain the path where the timestamp of the cause variable is earlier than the timestamp of the result variable.

[0045] Time chronology verification checks whether the causal relationship between variables conforms to a chronological order, meaning the causal variable should occur earlier than the outcome variable. The Granger causality test is a statistical test used to determine whether one time series can help predict another. The Granger causality test can determine the causal direction between variables, i.e., which variable is the causal variable and which is the outcome variable. Paths where the causal variable's timestamp is earlier than the outcome variable's timestamp are retained from the direct causal path candidate set; only those paths that conform to the chronological order are kept, while paths that do not conform to the chronological order are excluded.

[0046] When verifying the temporal sequence of a candidate set of direct causal paths and determining the causal direction using the Granger causality test, for each path in the candidate set, the timestamps of the two variables in the path are obtained. Then, the Granger causality test is used to determine the causal direction between the variables. The basic idea of ​​the Granger causality test is that if the lagged value of a variable X significantly improves the predictive ability of another variable Y, then X can be considered a Granger cause of Y. The specific process can be found in existing techniques and will not be elaborated here. In summary, two regression models are established (one including the lagged term of X, and one not), and the F-test is used to compare the goodness of fit of the two models. If the goodness of fit of the model including the lagged value of X is significantly better than that of the model not including the lagged value of X, then X can be considered a Granger cause of Y. Finally, it is checked whether the timestamp of the causal variable is earlier than the timestamp of the outcome variable. If so, the path is retained; otherwise, the path is excluded.

[0047] Step S235: Verify the validity of intermediate nodes in the candidate set of indirect causal paths, calculate the conditional correlation degree of intermediate nodes, and if the correlation strength of the original variable pair decreases by more than a preset proportion after removing an intermediate node, then the intermediate node is determined to be a valid transit node, and the indirect causal path containing the valid transit node is retained.

[0048] Conditional correlation strength refers to the intensity of the correlation between pairs of original variables, considering intermediate nodes. A decrease in the correlation strength of the original variable pairs after removing an intermediate node exceeds a preset percentage. If this condition is met, it indicates that the intermediate node played a significant role in transmitting the correlation between the original variable pairs, and thus is considered a valid transmission node. Retaining indirect causal paths containing valid transmission nodes means that only paths containing valid transmission nodes are retained from the candidate set of indirect causal paths, while paths without valid transmission nodes are excluded.

[0049] Specifically, for each path in the candidate set of indirect causal paths, intermediate nodes are identified. Then, the association strength of the original variable pair with and without the intermediate node is calculated. Association strength can be calculated using correlation coefficients, mutual information, or other methods. Next, the decrease in association strength is calculated as (r1-r2) / r1, where r1 is the association strength with and without the intermediate node, and r2 is the association strength after removing the intermediate node. Finally, the decrease in association strength is compared to a preset percentage. If the decrease exceeds the preset percentage, the intermediate node is considered a valid transit node, and the indirect causal path containing this valid transit node is retained.

[0050] Step S236: Integrate the verified direct and indirect causal paths to construct a hierarchical causal path list that includes path type, path length, and correlation strength.

[0051] Specifically, the verified direct and indirect causal paths can be merged into a set. Then, for each path in the set, its path type, path length, and association strength are recorded. The path type can be represented as "direct" or "indirect," the path length is the number of edges contained in the path, and the association strength can be the maximum value of the time-varying association strength curve used when constructing the adjacency matrix. Finally, this information is arranged in a certain order to form a hierarchical causal path list.

[0052] Step S240: Model the attenuation of the transmission effect of the indirect causal path, calculate the transmission efficiency coefficient of each intermediate node in the path, the transmission efficiency coefficient decreases exponentially with the increase of the path length, and generate the initial value of the indirect causal strength containing the attenuation factor.

[0053] Modeling the attenuation of transmission effects involves modeling the transmission effects of intermediate nodes in an indirect causal path, considering that the causal effect gradually diminishes as the path length increases. The transmission efficiency coefficient is an indicator used to measure the efficiency of intermediate nodes in transmitting causal effects, reflecting their ability to pass on the causal effect to the next node. The transmission efficiency coefficient decays exponentially with increasing path length; that is, the longer the path, the lower the transmission efficiency. The initial value of the indirect causal strength, including the attenuation factor, is the initial strength value of the indirect causal path calculated considering the attenuation of transmission effects. This value takes into account the transmission efficiency coefficients of each intermediate node in the path.

[0054] For example, suppose the relationship between the transmission efficiency coefficient e and the path length L is e = λ LWhere λ is the attenuation factor, 0 < λ < 1. The value of λ can be adjusted according to the actual situation, for example, determined through experiments or experience. Then, for each intermediate node in the indirect causal path, the transmission efficiency coefficient is calculated based on the length of its path. Next, the initial value of the indirect causal strength is calculated based on the transmission efficiency coefficient. Assuming the initial correlation strength of the indirect causal path is r, the initial value of the indirect causal strength after transmission through intermediate nodes is r × e1 × e2 × … × e n Where e1, e2, ..., e n It is the transmission efficiency coefficient of each intermediate node in the path.

[0055] Step S250: Eliminate interfering factors in hierarchical causal paths using the dynamic backdoor criterion, identify and block non-causal paths within different time windows, and retain valid causal paths that meet the time sequence constraints.

[0056] The dynamic backdoor criterion is used to identify and eliminate interfering factors in causal paths, helping to distinguish between causal and non-causal relationships. Identifying and blocking non-causal paths within different time windows involves analyzing hierarchical causal paths across different timeframes to identify and block those paths that are not true causal relationships. Retaining valid causal paths that conform to temporal constraints means only retaining those causal paths where the timestamp of the cause variable is earlier than the timestamp of the result variable, and which conform to causal logic.

[0057] For example, within each time window, a causal graph is constructed based on the hierarchical list of causal paths. A causal graph is a graph used to represent causal relationships between variables, clearly showing the causal structure between them. Then, backdoor paths in the causal graph are identified according to the dynamic backdoor criterion. Backdoor paths are those paths that connect cause and effect variables through non-causal associations. For example, in a causal graph, if there exists a path from a cause variable to an effect variable, and this path contains a non-causal intermediate node, then this path is a backdoor path. Next, the causal graph is adjusted to block backdoor paths, i.e., edges on backdoor paths are removed, so that only true causal paths are retained in the causal graph. Finally, it is checked whether the retained causal paths meet the time sequence constraint, i.e., whether the timestamp of the cause variable is earlier than the timestamp of the effect variable. If not, the causal graph is further adjusted until all retained causal paths meet the time sequence constraint.

[0058] Step S260: Integrate the intensity quantification value of the direct causal path with the attenuated intensity value of the indirect causal path, and perform weighted fusion according to the time window weight to generate a set of causal path weights that are dynamically adjusted over time.

[0059] The strength quantification value of a direct causal path is a quantitative representation of the influence of the causal variable on the outcome variable within that path, reflecting the strength of the direct causal relationship. The attenuated strength value of an indirect causal path is obtained after modeling the attenuation of transmission effects, taking into account the transmission efficiency of intermediate nodes and the path length. Time window weights are assigned based on the importance or reliability of different time windows. Different time windows may contribute differently to the final causal relationship determination, thus requiring weighting. Time window weights can be determined based on various factors, such as data reliability and the length of the time window.

[0060] Step S300: The causal path weight set and the time series pattern features in the sample dataset are fused using a Bayesian network to construct a time series representation vector with causal logic constraints.

[0061] In one implementation, step S300 may specifically include the following steps S310 to S360: Step S310: Perform multimodal feature extraction on the time series pattern features in the sample dataset. Extract voltage stability features, load distribution features, and line loss features from the power grid operation observation data. Combine these features with the time series pattern features of historical electricity sales records to form a multimodal input feature set.

[0062] For example, for the time series pattern characteristics of historical electricity sales records, time series analysis methods, such as the Autoregressive Integral Moving Average (ARIMA) model and seasonal decomposition methods, can be used to extract features such as trends, periodicity, and seasonality. For instance, the seasonal decomposition method can be used to decompose the electricity sales time series into trend components, seasonal components, and residual components, extracting features such as the slope of the trend component and the period and amplitude of the seasonal component. Then, voltage stability features, load distribution features, and line loss features are extracted from the power grid operation observation data. For voltage stability features, statistics such as the standard deviation, maximum value, and minimum value of voltage can be calculated to reflect voltage fluctuations. For load distribution features, indicators such as peak value, valley value, average value, and load factor can be statistically analyzed to describe the load distribution. For line loss features, the power loss of the line can be calculated based on parameters such as line resistance and current. Finally, the extracted time series pattern characteristics, voltage stability features, load distribution features, and line loss features of historical electricity sales records are combined to form a multimodal input feature set.

[0063] Step S320: Based on the causal path weight set, perform causal correlation screening on the multimodal input feature set, retain the features that have a direct causal correlation with the electricity sales risk prediction, and remove redundant features that have no causal correlation.

[0064] For example, the causal relationship and causal strength between each feature and electricity sales risk prediction are determined based on the causal path weight set. This can be achieved by reviewing the causal path list to identify paths related to electricity sales risk prediction and obtaining their weight values. Then, a causal strength threshold is set, retaining features with causal strength greater than the threshold, indicating a direct causal relationship between these features and electricity sales risk prediction. Features with causal strength less than the threshold are discarded, considered redundant features with no causal relationship.

[0065] Step S330: Perform time granular alignment on the filtered multimodal features, and convert features with different sampling frequencies into feature sequences with the same time interval, with the time interval consistent with the time stamp interval of the sample dataset.

[0066] Temporal granularity alignment adjusts features with different sampling frequencies to have the same time interval, ensuring temporal consistency. Different sampling frequencies arise from variations in data acquisition equipment or methods, resulting in different sampling time intervals for each feature. For example, some features might be sampled hourly, while others might be sampled daily. Converting features with different sampling frequencies into a feature sequence with the same time interval is achieved through interpolation or sampling methods. Maintaining consistency between the time interval and the time stamp interval of the sample dataset ensures that the time interval of the converted feature sequence matches the time stamp interval of the sample dataset, guaranteeing a temporal correspondence between the feature sequence and the sample dataset.

[0067] In one implementation, step S330 may specifically include the following steps S331 to S336: Step S331: Perform time sampling characteristic analysis on the selected multimodal features, statistically analyze the original sampling interval and sampling point distribution density of each feature, and generate a feature-sampling characteristic correspondence table containing feature identifier, sampling interval duration and data integrity index. The data integrity index represents the proportion of missing observations of the feature at each time node.

[0068] Specifically, each selected multimodal feature is first iterated through. For each feature, its start and end sampling times are recorded, and the original sampling interval is calculated by determining the time difference between adjacent sampling points. Simultaneously, the number of sampling points throughout the entire sampling period is counted, and the sampling point distribution density is calculated based on the length of the sampling period. For calculating the data integrity index, it is necessary to check whether each time node has an observation, count the number of time nodes with missing observations, and divide this number by the total number of time nodes to obtain the data integrity index. Finally, the feature identifier, sampling interval duration, and data integrity index are recorded in a feature-sampling characteristic correspondence table.

[0069] Step S332: Based on the feature-sampling characteristic correspondence table, the multimodal features are divided into two categories: short-interval sampling features and long-interval sampling features. The sampling interval of short-interval sampling features is less than the target time interval, and the sampling interval of long-interval sampling features is greater than or equal to the target time interval, thereby realizing the classification of the temporal characteristics of the features.

[0070] Short-interval sampling features, due to their sampling interval being shorter than the target time interval, have a relatively large data volume and may contain more detailed information on feature changes. For example, a current feature sampled every 5 minutes is a short-interval sampling feature relative to a target time interval of 1 hour. When unifying the time granularity of these features, compression processing is needed to reduce the data volume and extract key information. Long-interval sampling features have a sampling interval greater than or equal to the target time interval, resulting in relatively sparse data. For example, a device maintenance cost feature sampled monthly is a long-interval sampling feature relative to a target time interval of 1 day. These features require expansion processing to fill in the time intervals between adjacent sampling points, aligning them with the target time interval.

[0071] Specifically, the target time interval is first determined. Then, each feature in the feature-sampling characteristic correspondence table is traversed, and its sampling interval is compared with the target time interval. If the sampling interval is less than the target time interval, the feature is classified as a short-interval sampling feature; if the sampling interval is greater than or equal to the target time interval, the feature is classified as a long-interval sampling feature.

[0072] Step S333: Perform time-granularity compression processing on the short-interval sampling features, and use the sliding window trend extraction algorithm to calculate the feature fluctuation trend parameters within each target time interval. The trend parameters include the slope of the rising segment, the slope of the falling segment, and the fluctuation amplitude. The feature extreme points within the window are retained, and a compressed feature sequence containing trend parameters and extreme points is generated. The time interval of the compressed feature sequence is consistent with the target time interval.

[0073] For example, the window length is set based on the target time interval and the sampling interval of the short-interval sampling features. For instance, if the target time interval is 1 hour and the sampling interval of the short-interval sampling features is 15 minutes, the window length can be set to 4 sampling points. Then, the window slides across the short-interval sampling feature sequence. For the data within each window, the following operations are performed: First, identify the rising and falling segments within the window. The rising segment is the time period where the feature value gradually increases, and the falling segment is the time period where the feature value gradually decreases. The rising and falling segments are obtained by calculating the ratio of the change in feature value to the change in time between the start and end points of the rising and falling segments. The fluctuation amplitude is the difference between the maximum and minimum values ​​of the feature within the window. Simultaneously, the feature extrema, i.e., the maximum and minimum values, are recorded within the window. Finally, the trend parameters and extrema of each window are arranged according to the target time interval to generate a compressed feature sequence.

[0074] Step S334: Perform time granularity expansion processing on the long-interval sampling features. Construct a trend extrapolation model based on the feature value sequence and sampling interval of historical sampling points. Extrapolate the time interval between adjacent historical sampling points using the trend extrapolation model to generate a preliminary estimate of the intermediate time feature under the target time interval. Calculate the fluctuation compensation amount based on the historical data of this long-interval sampling feature. The fluctuation compensation amount is determined by analyzing the fluctuation patterns between historical sampling points. Superimpose the fluctuation compensation amount onto the preliminary estimate generated by the trend extrapolation to generate an extended feature sequence containing trend components and fluctuation components.

[0075] For long-interval sampling features, due to the large sampling interval and relatively sparse data, time granularity expansion processing is required to fill the time intervals between adjacent sampling points. Various methods can be used to construct models based on the feature value sequences and sampling intervals of historical sampling points, such as linear regression, multinomial regression, and exponential smoothing.

[0076] However, relying solely on trend extrapolation may not accurately reflect the true changes in characteristics, as characteristics often fluctuate during actual changes. Therefore, it is necessary to calculate fluctuation compensation. The calculation of fluctuation compensation is based on historical data of the long-interval sampled feature, determined by analyzing fluctuation patterns between historical sampling points. For example, observing the fluctuations in sales over the past few months reveals the patterns and magnitudes of these fluctuations. The differences between adjacent historical sampling points can be calculated, and the distribution of these differences can be analyzed to determine the range and frequency of fluctuations. Based on the fluctuation patterns, the fluctuation compensation amount for each intermediate time point is calculated.

[0077] By superimposing the volatility compensation onto the initial estimate generated by trend extrapolation, we obtain an extended feature sequence that includes both trend and volatility components. This extended feature sequence more accurately reflects the true changes of long-interval sampling features over the target time interval.

[0078] Step S335: Perform time stamp alignment verification on the compressed feature sequence and the expanded feature sequence. Using the time stamp of the sample dataset as a reference, check whether each feature sequence has a valid feature value at the corresponding time node. For missing time nodes, fill them with a weighted average of the feature values ​​of the preceding and following time nodes. The weight is inversely proportional to the time distance.

[0079] Specifically, the time stamps of the compressed and expanded feature sequences are mapped one-to-one with the time stamps of the sample dataset. For each time point, the existence of valid feature values ​​is checked. If missing time points exist, imputation values ​​are calculated based on the feature values ​​and time distances of the preceding and following time points. When calculating the imputation values, it is necessary to ensure that valid feature values ​​exist for the preceding and following time points. If missing values ​​also exist for the preceding and following time points, it may be necessary to further expand the search scope to find more suitable preceding and following time points. Finally, the imputed compressed and expanded feature sequences are perfectly aligned with the time stamps of the sample dataset.

[0080] Step S336: Calculate the time synchronization error of each feature sequence after alignment. The time synchronization error is the sum of squares of the deviations between the actual time stamp and the target time stamp of the feature sequence. If the error exceeds the preset threshold, readjust the trend extraction window size or extrapolation model parameters until all feature sequences meet the time synchronization requirements and generate a set of multimodal feature sequences with uniform time granularity.

[0081] The preset threshold is an upper limit for error set in advance based on actual needs and data characteristics. If the time synchronization error exceeds the preset threshold, it indicates a significant problem with the temporal alignment of the feature sequences, requiring readjustment of the trend extraction window size or extrapolation model parameters. For short-interval sampling features, adjusting the trend extraction window size can change the length of the sliding window, thus affecting the generation of compressed feature sequences. For long-interval sampling features, adjusting the extrapolation model parameters can alter the trend extrapolation results, making them more consistent with reality.

[0082] Step S340: Input the aligned multimodal features and causal path weight set into the evidence layer of the Bayesian network, activate the feature states of each time node in chronological order, and constrain the probability dependency direction between nodes through causal path weights.

[0083] Activating the feature states of each time node sequentially is to simulate the dynamic changes in time series data. At each time node, the multimodal feature value of that time node is input into the evidence layer of the Bayesian network to activate the corresponding node. Constraining the probabilistic dependency direction between nodes through causal path weights determines the probabilistic dependency between nodes in the Bayesian network based on the weight values ​​of each causal path in the causal path weight set. Larger weight values ​​indicate a stronger causal relationship and a higher degree of probabilistic dependency between nodes.

[0084] In practical implementation, the structure of the Bayesian network can be constructed first, determining the relationships between nodes and edges. Nodes represent multimodal features and other relevant variables, while edges represent causal relationships between variables. Then, the aligned multimodal features and the set of causal path weights are input into the evidence layer of the Bayesian network. At each time point, the state of the corresponding node in the evidence layer is set according to the multimodal feature value at that time point. For example, if the value of a feature is high at that time point, the corresponding node state is set to high. Next, the conditional probability table between nodes in the Bayesian network is adjusted according to the set of causal path weights. The conditional probability table describes the probabilistic dependencies between nodes, which can be adjusted by the causal path weights to ensure that the probabilistic dependencies between nodes conform to causal relationships. Finally, the feature states at each time point are processed sequentially in chronological order, completing the activation process of the Bayesian network.

[0085] Step S350: Perform probabilistic inference on the activated Bayesian network using a variational inference algorithm, calculate the posterior probability distribution of the latent variables at each time point, and comprehensively characterize the joint influence of multimodal features and causal paths to generate an intermediate feature vector containing probability distribution parameters.

[0086] Variational inference algorithms are methods for approximating complex probability distributions. In Bayesian networks, they can help compute the posterior probability distribution of latent variables. Latent variables are variables that cannot be directly observed in a Bayesian network and have probabilistic dependencies on observable multimodal features and causal paths. The posterior probability distribution is the probability distribution of the latent variables after considering the observed multimodal features and causal path information; it comprehensively represents the joint influence of multimodal features and causal paths. Generating intermediate feature vectors containing probability distribution parameters is to quantify and represent the information of the posterior probability distribution. These intermediate feature vectors contain parameters of the posterior probability distribution, such as mean and variance, which reflect the characteristics and state of the latent variables.

[0087] In one implementation, step S350 may specifically include the following steps S351 to S356: Step S351: Initialize the variational distribution parameters of the latent variables in the Bayesian network. Determine the initial values ​​of the mean vector and covariance matrix based on the prior probability distribution of the latent variables. The mean vector represents the initial estimate of the central tendency of the latent variables, and the covariance matrix represents the initial correlation between the dimensions of the latent variables.

[0088] For example, first determine the method for calculating the mean and covariance based on the type of the prior probability distribution of the latent variable. If the prior probability distribution is Gaussian, then the mean and covariance of the Gaussian distribution can be directly used as initial values. For example, for a one-dimensional latent variable, its prior probability distribution is a Gaussian distribution N(μ0,σ0). 2 If the initial value of the mean vector is μ0, then the initial value of the covariance matrix is ​​σ0. 2 For multidimensional latent variables, the mean vector and covariance matrix need to be calculated based on the multidimensional form of their prior probability distribution. Then, the calculated mean vector and covariance matrix are used as the initial values ​​for the variational distribution parameters.

[0089] Step S352: At each time point, calculate the likelihood function value based on the aligned multimodal feature sequence and the causal path weight set. The likelihood function value represents the probability density of the corresponding value of the latent variable under the current feature state. Adjust the contribution ratio of different modal features to the likelihood function through the causal path weight.

[0090] First, determine the form of the likelihood function. An appropriate likelihood function, such as a Gaussian likelihood function or a Bernoulli likelihood function, can be selected based on the distribution type of the multimodal features. For each modal feature, determine its contribution to the likelihood function based on the causal path weight between it and the latent variables. For example, if a modal feature has a large causal path weight with the latent variables, it indicates that the feature has a significant impact on the latent variables, and a larger weight should be given to this feature when calculating the likelihood function. At each time point, combine the likelihood function values ​​of each modal feature to obtain the likelihood function value for that time point. The combination method can be chosen according to specific circumstances, such as multiplication or addition. For example, for multiple independent modal features, their likelihood function values ​​can be multiplied to obtain the total likelihood function value.

[0091] Step S353: Construct a variational lower bound function as the optimization objective. The variational lower bound function includes the expectation term of the likelihood function and the entropy term of the variational distribution of the latent variables. Maximize the variational lower bound function to complete the approximate estimation of the posterior probability distribution.

[0092] For example, first, based on the definitions of the likelihood function and the variational distribution, the expectation term of the likelihood function is calculated. Calculating the expectation term requires integrating over the variational distribution, which can typically be approximated using the Monte Carlo method. Then, the entropy term of the latent variable variational distribution is calculated. Entropy is a measure of the uncertainty of a probability distribution; for common distributions such as the Gaussian distribution, there are corresponding entropy formulas. Finally, the expectation term of the likelihood function and the entropy term of the latent variable variational distribution are added together to obtain the variational lower bound function. When maximizing the variational lower bound function, optimization algorithms, such as the stochastic gradient ascent algorithm, can be used. By iteratively updating the parameters of the variational distribution, the value of the variational lower bound function gradually increases until convergence.

[0093] Step S354: Iteratively update the variational distribution parameters using the stochastic gradient ascent algorithm. In each iteration, randomly sample time window samples from the multimodal feature sequence, calculate the gradient of the variational lower bound function with respect to the mean vector and covariance matrix, and adjust the parameter values ​​according to the gradient direction until the function values ​​converge.

[0094] Each iteration randomly samples a time window sample from the multimodal feature sequence, thus reducing computational burden while maintaining data diversity. The gradient of the variational lower bound function with respect to the mean vector and covariance matrix is ​​calculated; the gradient represents the direction of change of the variational lower bound function with respect to the current parameter values. The parameter values ​​are adjusted according to the gradient direction, i.e., increasing the parameter values ​​along the positive direction of the gradient, so that the value of the variational lower bound function gradually increases. First, the number of iterations and the learning rate are set. The learning rate controls the step size of parameter updates in each iteration. In each iteration, a time window sample is randomly selected from the multimodal feature sequence, and the gradient of the variational lower bound function with respect to the mean vector and covariance matrix is ​​calculated based on this sample. Then, the parameter values ​​of the mean vector and covariance matrix are updated according to the gradient and the learning rate. This process is repeated until the value of the variational lower bound function converges, i.e., the change in the function value in consecutive iterations is less than a preset threshold.

[0095] Step S355: After the variational lower bound function converges, extract the mean vector and the diagonal elements of the covariance matrix of the variational distribution. The mean vector elements represent the optimal estimates of each dimension of the latent variables, and the diagonal elements of the covariance matrix represent the degree of uncertainty of the estimates of each dimension.

[0096] Once the variational lower bound function converges, it indicates that the variational distribution has approximated the true posterior probability distribution as closely as possible. At this point, extracting the mean vector and the diagonal elements of the covariance matrix of the variational distribution yields the optimal estimates of the latent variables and the degree of uncertainty in the estimates of each dimension.

[0097] The elements of the mean vector represent the optimal estimates of the latent variables in each dimension. The diagonal elements of the covariance matrix represent the variance of the estimates in each dimension; the larger the variance, the higher the uncertainty of the estimate in that dimension. The specific process of extracting the mean vector and the diagonal elements of the covariance matrix is ​​as follows: After the variational lower bound function converges, record the current mean vector and covariance matrix. Then, extract the diagonal elements from the covariance matrix to obtain the variance of each dimension.

[0098] Step S356: Concatenate the mean vector and covariance diagonal elements of each time node in chronological order to obtain a feature matrix containing the time dimension. The rows of the matrix correspond to the time nodes, and the columns correspond to the distribution parameters of the latent variables. This feature matrix is ​​the set of intermediate feature vectors containing probability distribution parameters.

[0099] Concatenating the mean vector and covariance diagonal elements of each time point in chronological order integrates the information of latent variables at different time points, forming a feature matrix that includes a time dimension. The rows of the matrix correspond to time points, and the columns correspond to the distribution parameters of the latent variables, thus clearly showing the state and characteristics of the latent variables at different time points.

[0100] Step S360: Perform temporal correlation enhancement processing on the intermediate feature vector. Capture the dependency relationship between features at different time nodes through a temporal convolutional network to generate a temporal representation vector that integrates causal logic and temporal correlation. Each dimension of the temporal representation vector corresponds to the comprehensive feature state of the corresponding time node.

[0101] Temporal Convolutional Networks (TCNs) are convolutional neural networks specifically designed for processing time-series data. They can effectively capture the dependencies between different time points in a time series. Temporal correlation enhancement processing of intermediate feature vectors aims to better extract the temporal information and causal logic contained within them.

[0102] Temporal convolutional networks (TCNNs) process intermediate feature vectors along the temporal dimension through convolution operations. The convolutional kernel slides along the temporal dimension, performing convolution operations on features within each time window to extract local temporal features. Through multiple layers of convolution and pooling operations, the receptive field can be gradually expanded to capture dependencies over a longer time span. Techniques such as causal convolution and dilated convolution can also be introduced into temporal convolutional networks to better adapt to the causality and sparsity of time-series data. Causal convolution ensures that when processing features at a given time point, the network only uses features from that time point and earlier, conforming to the causal relationships in time-series data. Dilated convolution can expand the receptive field without increasing the kernel size, improving the network's ability to process long sequences. Generating temporal representation vectors that integrate causal logic and temporal correlation is the ultimate goal of temporal convolutional networks. Each dimension of the temporal representation vector corresponds to the comprehensive feature state of the corresponding time point, containing not only information from intermediate feature vectors but also incorporating dependencies and causal logic between different time points.

[0103] Step S400: Input the time series representation vector into the meta-learning framework model, generate an initial risk predictor through parameter fine-tuning, and output the probability distribution sequence of electricity sales fluctuations.

[0104] In one implementation, step S400 may specifically include the following steps S410 to S460: Step S410: Construct a meta-learning framework model that includes a feature encoding module, a meta-knowledge storage module, a task adaptation module, and a prediction output module. The feature encoding module is used for nonlinear transformation of time-series representation vectors, the meta-knowledge storage module stores empirical parameters of historical prediction tasks, the task adaptation module realizes cross-task knowledge transfer, and the prediction output module generates probability distribution sequences.

[0105] The meta-knowledge storage module is one of the core components of the meta-learning framework, storing empirical parameters from historical prediction tasks. These empirical parameters are learned by the model while processing multiple related tasks, containing common knowledge about different tasks. By storing and utilizing these empirical parameters, the model can quickly adjust to new tasks, reducing learning time and data requirements.

[0106] The task adaptation module implements cross-task knowledge transfer functionality. Based on the characteristics of the current task, it selects appropriate empirical parameters from the meta-knowledge storage module and adjusts the weights of these parameters to enable the model to adapt to new prediction tasks. The attention mechanism plays a crucial role in the task adaptation module, dynamically adjusting parameter weights according to the importance of different historical meta-knowledge points to the current task.

[0107] The prediction output module generates a probability distribution sequence of electricity sales fluctuations based on the output of the task adaptation module. This module can use probability distribution fitting methods, such as Gaussian mixture models or kernel density estimation, to convert the output of the task adaptation module into a probability distribution sequence.

[0108] When constructing a meta-learning framework model, the structure of the feature encoding module is designed first. Neural network structures such as Multilayer Perceptrons (MLPs) or Convolutional Neural Networks (CNNs) can be used to perform nonlinear transformations on the temporal representation vectors. Next, a meta-knowledge storage module is constructed, selecting an appropriate data structure to store empirical parameters from historical tasks, such as dictionaries or matrices. Then, a task adaptation module is designed, introducing an attention mechanism to adjust the weight allocation of meta-parameters. Finally, a prediction output module is constructed, selecting an appropriate probability distribution fitting method to generate the probability distribution sequence. Combining these modules forms a complete meta-learning framework model.

[0109] Step S420: Input the temporal representation vector into the feature encoding module, and perform deep extraction of temporal features through a stacked long short-term memory network. The hidden states of the long short-term memory network contain feature information at different time scales, generating a high-dimensional abstract feature vector.

[0110] Stacked LSTMs consist of multiple LSTM layers stacked sequentially. Each LSTM layer processes the input features, extracting temporal features at different levels. By stacking multiple layers, the receptive field can be gradually expanded, capturing dependencies over a longer time span.

[0111] The hidden states of Long Short-Term Memory (LSTM) networks contain feature information at different time scales. At each time step, the LSTM updates the current hidden state based on the current input and the hidden state of the previous time step. The hidden state not only contains feature information from the current time step but also retains important information from past time steps. By stacking LSTMs, the hidden states of different layers can capture features at different time scales, from short-term local features to long-term global features.

[0112] The high-dimensional abstract feature vector is the output of the stacked LSTM. This vector contains rich temporal feature information inherent in the time series representation vector and has undergone nonlinear transformation, giving it higher expressive power. The high-dimensional abstract feature vector can be used as input to subsequent modules for further analysis and processing.

[0113] Step S430: Input the high-dimensional abstract feature vector into the meta-knowledge storage module and perform dynamic similarity matching with the stored historical task meta-knowledge. The similarity matching is based on the comprehensive calculation of the cosine similarity of the feature vector and the causal structure similarity of the task scenario, and extracts the meta-parameters of the most similar historical task.

[0114] In one implementation, step S430 may specifically include the following steps S431 to S435: Step S431: Read the meta-knowledge records of historical prediction tasks from the meta-knowledge storage module. Each meta-knowledge record contains the high-dimensional abstract feature vector of the task, the set of causal path weights, and the corresponding model parameters.

[0115] The meta-knowledge storage module stores information related to historical prediction tasks. Reading the meta-knowledge records of historical prediction tasks is the first step in dynamic similarity matching. Each meta-knowledge record contains key information about the task; the high-dimensional abstract feature vector is the feature representation of the task, containing its main feature information. The causal path weight set describes the causal relationships and strengths between variables in the task, reflecting the task's causal structure. The corresponding model parameters are those learned by the model when processing this historical task; these parameters contain empirical knowledge about the task.

[0116] Step S432: Calculate the cosine similarity between the high-dimensional abstract feature vector of the current task and the feature vectors of each historical task. The cosine similarity measures the degree of similarity in the feature space. The larger the value, the closer the feature distributions are.

[0117] Step S433: Extract the causal path weight set of the current task and the historical task, and calculate the structural similarity between the two. The structural similarity is calculated by the graph edit distance algorithm to measure the similarity of the causal network topology. The smaller the graph edit distance, the more similar the structure.

[0118] Specifically, the causal path weights of the current and historical tasks are first converted into a graph representation. Then, a graph edit distance algorithm is used to calculate the edit distance between the two graphs. Existing graph edit distance algorithms can be used, such as dynamic programming-based or search-based algorithms. Finally, the graph edit distance is converted into a structural similarity metric, for example, 1 / (1+d), where d is the graph edit distance. The resulting structural similarity metric ranges from [0,1], with values ​​closer to 1 indicating greater structural similarity.

[0119] Step S434: Perform weighted fusion of cosine similarity and structural similarity, and use the fusion result as a comprehensive similarity index. Sort the historical tasks according to the comprehensive similarity index and select the top few historical tasks with the highest similarity.

[0120] Weighted fusion is a method that comprehensively considers cosine similarity and structural similarity. By assigning different weights to each, their contribution to the overall similarity index can be adjusted according to the actual situation. The overall similarity index can more comprehensively measure the similarity between the current task and historical tasks. The specific formula for weighted fusion of cosine similarity and structural similarity is: S = αC + (1-α)Ss Where S is the comprehensive similarity index, C is the cosine similarity, and S s This refers to structural similarity, where α is the weighting coefficient, ranging from [0,1]. The value of α can be adjusted according to the actual situation. If more emphasis is placed on the similarity of the feature space, α should be set larger; if more emphasis is placed on the similarity of the causal network topology, α should be set smaller. Historical tasks are sorted according to the comprehensive similarity index, with the historical tasks with the highest similarity placed first. The top few historical tasks with the highest similarity are selected, and these historical tasks will serve as the basis for subsequent extraction of meta-parameters.

[0121] Step S435: Perform a weighted average of the meta-parameters of the selected historical tasks, with the weights being the normalized result of the comprehensive similarity index of each task, to generate the initial values ​​of the meta-parameters for the current task.

[0122] For each meta-parameter, the corresponding meta-parameter of the selected historical task is multiplied by its normalized weight, and then all results are summed.

[0123] Step S440: Initialize the network parameters of the task adaptation module based on the extracted meta-parameters. The task adaptation module adjusts the weight distribution of the meta-parameters through an attention mechanism. The weight distribution represents the applicability of different historical meta-knowledge to the current prediction task.

[0124] Specifically, the attention mechanism determines the importance of each historical meta-knowledge by calculating the similarity between current task features and historical task features. The higher the similarity, the greater the weight of the corresponding historical meta-knowledge in the current task, and the more significant its impact on the model.

[0125] Step S450: Fine-tune the parameters of the task adaptation module through the inner loop optimization algorithm. The inner loop takes the prediction loss of the current task as the optimization target and updates the network parameters through gradient descent. During the learning process, the parameters of the feature encoding module and the meta-knowledge storage module are fixed.

[0126] Fixing the parameters of the feature encoding module and the meta-knowledge storage module during the learning process ensures that the model can fully utilize previously learned general knowledge while quickly adapting to new tasks. The feature encoding module has already effectively extracted features from the input temporal representation vector, and the meta-knowledge storage module stores empirical parameters from historical tasks. Fixing their parameters avoids destroying this knowledge during fine-tuning, allowing the model to focus on adjusting the parameters of the task adaptation module to suit the characteristics of the current task.

[0127] In one implementation, step S450 specifically includes the following steps S451 to S456: Step S451: Perform preheating training on the task adaptation module based on the initial values ​​of the meta-parameters, calculate the initial prediction loss by inputting some temporal representation vector samples, and generate a task complexity index by combining the feature entropy value. The feature entropy value is calculated by the negative logarithmic expectation of the feature probability distribution, and the task complexity index is the weighted sum of the initial prediction loss and the feature entropy value.

[0128] Warm-up training is to allow the task adaptation module to undergo preliminary adaptation to the current task before formal parameter fine-tuning. The task adaptation module is initialized based on initial meta-parameter values, and then forward propagation is performed using a subset of temporal representation vector samples to obtain initial prediction results. These prediction results are then compared with the true labels to calculate the initial prediction loss. The initial prediction loss reflects the model's prediction error for the current task under the current parameters.

[0129] Feature entropy is calculated using the negative log-expectation of the feature's probability distribution. Feature entropy measures the degree of uncertainty of a feature; a higher entropy value indicates a more dispersed distribution and higher uncertainty. Calculating feature entropy first requires estimating the feature's probability distribution. Methods such as kernel density estimation and histograms can be used to estimate this distribution. Then, the negative log-expectation is calculated based on the probability distribution to obtain the feature entropy value.

[0130] The task complexity metric is a weighted sum of the initial prediction loss and the feature entropy value. This weighted sum comprehensively considers both the model's prediction error and the uncertainty of the features. The weight allocation can be adjusted according to the specific situation. If more attention is paid to the model's prediction error, the weight of the initial prediction loss can be set larger; if more attention is paid to the uncertainty of the features, the weight of the feature entropy value can be set larger.

[0131] Step S452: Compare the task complexity index with a preset threshold. When the task complexity index exceeds the preset threshold, use incremental learning rate scheduling. The learning rate is adjusted according to the reciprocal of the square root of the number of iterations. When the task complexity index does not exceed the preset threshold, use constant learning rate scheduling. Divide the batch sample set according to the ratio of the complexity index to the threshold. When the ratio is greater than 1, refine the granularity of the sample set division to improve the targeting of parameter updates.

[0132] When the task complexity metric exceeds a preset threshold, it indicates that the current task is relatively complex, and the model needs more fine-tuning of its parameters. A progressive learning rate scheduling method is used, where the learning rate is adjusted according to the reciprocal of the square root of the iteration count, i.e., the learning rate η. t=η0 / √t, where η0 is the initial learning rate and t is the number of iterations. This allows the model to converge quickly with a larger learning rate in the early stages of training. As the number of iterations increases, the learning rate gradually decreases, preventing the model from oscillating or failing to converge in the later stages. When the task complexity metric does not exceed the preset threshold, it indicates that the current task is relatively simple, and constant learning rate scheduling can be used. Constant learning rate scheduling means that the learning rate remains unchanged throughout the training process. Simultaneously, the batch sample set is divided according to the ratio of the complexity metric to the threshold. If the ratio is greater than 1, it indicates that although the task complexity does not exceed the threshold, it is relatively high. In this case, the granularity of the sample set partitioning is refined, that is, the batch sample set is divided into smaller subsets for training. This improves the targeting of parameter updates, allowing the model to adapt more accurately to the current task.

[0133] Step S453: Associate and map the network parameters of the task adaptation module with the set of causal path weights. Based on the comparison results of the weight values ​​of the causal paths corresponding to the parameters and the preset weight thresholds, divide the key parameter layer and the non-key parameter layer. The core causal path features have weight values ​​greater than the preset weight thresholds, and the auxiliary causal path features have weight values ​​less than or equal to the preset weight thresholds.

[0134] The preset weight threshold is a pre-defined critical value used to distinguish between core causal paths and auxiliary causal paths. Based on the comparison between the weight values ​​of the causal paths corresponding to parameters and the preset weight threshold, the network parameters are divided into critical parameter layers and non-critical parameter layers. Causal paths corresponding to parameters with weight values ​​greater than the preset weight threshold are considered core causal paths; these paths have a significant impact on electricity sales risk prediction, and their corresponding parameter layers are the critical parameter layers. Causal paths corresponding to parameters with weight values ​​less than or equal to the preset weight threshold are auxiliary causal paths, and their corresponding parameter layers are the non-critical parameter layers.

[0135] This division allows for different optimization strategies to be applied to the critical parameter layer and the non-critical parameter layer, improving the training efficiency and performance of the model. For example, a more refined optimization method can be used for the critical parameter layer to ensure that the model can accurately learn the information of the core causal path; while a more relaxed optimization method can be used for the non-critical parameter layer to reduce computational cost.

[0136] Step S454: Use momentum optimization strategy to update gradients for key parameter layers. Adjust the current gradient update direction by accumulating historical gradient directions to maintain the convergence trend of parameters on the core causal path. Use gradient pruning strategy for non-key parameter layers to limit the gradient update magnitude within the preset gradient norm to maintain the overall stability of the network.

[0137] Employing momentum optimization strategies at key parameter layers can help the model converge better on the core causal path. Specifically, momentum optimization introduces a momentum term, which is a weighted sum of historical gradients. With each parameter update, not only the current gradient but also the cumulative effect of historical gradients is considered. This reduces gradient oscillations, making parameter updates smoother and accelerating convergence. Gradient clipping strategies are used to limit the magnitude of gradient updates. Applying gradient clipping at non-key parameter layers limits the gradient update magnitude to a preset gradient norm. When the gradient norm exceeds a preset threshold, the gradient is scaled to make its norm equal to the threshold. This avoids gradient explosion and maintains the overall stability of the network. During training, gradient explosion can cause drastic changes in model parameters, preventing convergence. Gradient clipping effectively controls the magnitude of gradients, ensuring a more stable training process.

[0138] Step S455: During the fine-tuning process, calculate the Barthel distance between the current task feature distribution and the historical task feature distribution in the meta-knowledge storage module in real time. The Barthel distance is calculated by the square root integral of the feature distribution probability density function. When the Barthel distance exceeds the preset distance threshold, extract the supplementary similar task meta-parameters from the meta-knowledge storage module and incorporate them into the current parameter fine-tuning process.

[0139] The preset distance threshold is a critical value for Bach distance pre-set based on experiments and experience. When the Bach distance exceeds the preset threshold, it indicates that the feature distribution of the current task differs significantly from the feature distribution of historical tasks, and the model may not be able to fully utilize the meta-knowledge of historical tasks. In this case, supplementary meta-parameters of similar tasks are extracted from the meta-knowledge storage module and incorporated into the current parameter fine-tuning process. These supplementary meta-parameters can provide the model with more information relevant to the current task, helping the model to better adapt to new tasks.

[0140] When extracting supplementary meta-parameters for similar tasks, historical tasks can be sorted according to Bach distance, and meta-parameters from historical tasks with smaller Bach distances to the current task can be selected. These meta-parameters are then fused with the parameters of the current model. The fusion method can be a weighted average, where the weights can be adjusted based on the Bach distance; the smaller the Bach distance, the greater the weight of the corresponding historical task meta-parameter.

[0141] Step S456: Monitor the continuous rate of change of the predicted loss of the validation set and the average value of the parameter update magnitude. When the rate of change of the validation set loss is less than the preset ratio and the average value of the parameter update magnitude is lower than the minimum adjustment threshold for multiple consecutive rounds, stop the parameter fine-tuning iteration and save the parameter configuration of the current task adaptation module.

[0142] Monitoring the continuous rate of change of the validation set prediction loss and the mean of the parameter update magnitude is crucial for determining whether the model has converged. The continuous rate of change of the validation set prediction loss reflects the model's performance changes on the validation set. If the rate of change of the validation set loss is less than a preset proportion for multiple consecutive iterations, it indicates that the model's performance has stabilized, and further parameter adjustments may not yield significant performance improvements. The mean of the parameter update magnitude reflects the degree of change of the model parameters in each iteration. When the mean of the parameter update magnitude is lower than the minimum adjustment threshold, it indicates that the model's parameters have stabilized and no further significant adjustments are needed. When both conditions are met—a continuous rate of change of the validation set loss being less than a preset proportion and the mean of the parameter update magnitude being lower than the minimum adjustment threshold—the parameter fine-tuning iteration stops. At this point, the model is considered to have converged, and the parameter configuration of the current task adaptation module is saved.

[0143] Step S460: After fine-tuning, input the output features of the task adaptation module into the prediction output module, and generate a probability distribution sequence containing the probabilities of different fluctuation ranges through probability distribution fitting, thus completing the construction of the initial risk predictor.

[0144] Probability distribution fitting involves finding a suitable probability distribution to describe the distribution patterns of known data or features. Various probability distribution models can be used in this process, such as Gaussian, Poisson, and gamma distributions. Choosing the appropriate probability distribution model requires considering the characteristics of the data and the needs of the practical problem. For example, if electricity sales fluctuations exhibit an approximately symmetrical distribution, a Gaussian distribution might be a suitable choice; if the data is count-based, a Poisson distribution might be more appropriate.

[0145] When fitting a probability distribution, it is necessary to estimate the parameters of the probability distribution. Methods such as maximum likelihood estimation and moment estimation can be used to estimate these parameters. Maximum likelihood estimation is a commonly used parameter estimation method that finds the most probable parameter values ​​by maximizing the likelihood function. The likelihood function is the probability of observing data given the parameters. By solving for the maximum value of the likelihood function, the estimated values ​​of the parameters can be obtained. Through probability distribution fitting, a probability distribution sequence containing probabilities of different fluctuation ranges is generated. This probability distribution sequence can provide decision-makers with information about the uncertainty of electricity sales fluctuations. For example, decision-makers can calculate the probability of different fluctuation ranges based on the probability distribution sequence, assess the risk level of electricity sales fluctuations, and thus formulate corresponding decision-making strategies. After generating the probability distribution sequence, the initial risk predictor is completed, which can be used to predict and assess the risks of future electricity sales fluctuations.

[0146] Step S500: Obtain the deviation value between the real-time electricity sales observation value and the probability distribution sequence, adjust the parameter configuration of the initial risk predictor based on the deviation value, update the causal strength quantification value in the causal path weight set, and generate the electricity sales risk assessment result.

[0147] In one implementation, step S500 may specifically include the following steps S510 to S560: Step S510: Receive continuous electricity sales observations through the real-time data interface, divide the observations into equal-length data segments according to the time stamp, each segment contains observations from multiple consecutive time nodes, and obtain real-time data blocks that match the time granularity of the probability distribution sequence. The time stamps of the data blocks correspond one-to-one with the time nodes of the probability distribution sequence.

[0148] A real-time data interface is a channel used to receive real-time electricity sales observations. It can connect to power system monitoring equipment, data acquisition systems, etc., to obtain electricity sales data in real time. Continuous electricity sales observations are electricity sales data collected over a period of time, and these data have temporal order and continuity.

[0149] Dividing observations into equal-length segments based on time stamps is a way to group consecutive observations to match the time granularity of the probability distribution sequence. Time stamps can be specific dates, timestamps, etc., used to identify the time each observation was collected. The length of each equal-length segment can be determined based on the time granularity of the probability distribution sequence.

[0150] Step S520: Align the real-time data block with the probability distribution sequence on the time axis, calculate the absolute difference between the observation value at each time node and multiple quantile values ​​of the corresponding node probability distribution, and generate a deviation sequence containing multiple deviation components, with the number of deviation components being the same as the number of quantile values.

[0151] Timeline alignment involves matching real-time data blocks and probability distribution sequences in time, ensuring that the observed and predicted values ​​at each time point correspond correctly. During timeline alignment, it's necessary to check if the time stamps of the real-time data blocks and the probability distribution sequences are consistent. If time stamp mismatches exist, interpolation or sampling may be required to ensure timeline consistency.

[0152] In one implementation, step S520 may specifically include the following steps S521 to S526: Step S521: Extract the time stamp set of real-time data blocks and probability distribution sequences, determine common time nodes by time stamp matching, and supplement time nodes with missing observations by interpolation calculation of observations from adjacent valid time nodes. The interpolation weight is dynamically allocated according to the time distance.

[0153] The timestamp set contains the time information for each data point. By comparing these two sets, their common time nodes can be found. Common time nodes are those that exist in both the real-time data block and the probability distribution sequence; these nodes form the basis for subsequent comparisons.

[0154] For time nodes with missing observations, interpolation calculations using observations from adjacent valid time nodes are used to fill the gaps. Interpolation is a method of estimating unknown data points using known data points. In this case, observations from adjacent valid time nodes are used to estimate the observations at the missing time node. The interpolation weights are dynamically allocated based on time distance; that is, the closer a valid time node is to the missing time node, the greater the weight of its observations in the interpolation calculation. For example, if the missing time node is located between two valid time nodes and is closer to the previous valid time node, then the observations from the previous valid time node will have a relatively larger weight in the interpolation calculation.

[0155] Step S522: Read multiple quantile values ​​for each common time node from the probability distribution sequence. The quantile values ​​are determined by the cumulative probability corresponding to the probability distribution function, reflecting the prediction boundary values ​​under different confidence levels.

[0156] Quantile values ​​are important indicators in probability distributions, reflecting the prediction boundary values ​​at different confidence levels. For example, the 50th quantile (median) indicates a 50% probability that the observed value will be less than that value; the 90th quantile indicates a 90% probability that the observed value will be less than that value.

[0157] This process involves reading multiple quantile values ​​for each common time point from a probability distribution sequence. When reading quantile values, the cumulative probability of the probability distribution function must be used to determine them. For common probability distributions, such as the Gaussian distribution, the corresponding cumulative distribution function table or calculation method can be used to determine the quantile values.

[0158] Step S523: Subtract the observed value of each common time node from the multiple quantile values ​​of the corresponding node, take the absolute value of the difference, generate multiple deviation components for each time node, and arrange them in chronological order to form a multidimensional deviation matrix.

[0159] The observations at each common time point are subtracted from the corresponding quantile values, resulting in multiple difference values. The absolute values ​​of these differences are taken to eliminate the sign and focus only on their magnitude. The multiple bias components generated at each time point reflect the degree of difference between the observed value and the probability distribution at different quantiles. These bias components are arranged in chronological order to form a multidimensional bias matrix. Each row of the multidimensional bias matrix corresponds to a time point, and each column corresponds to a quantile value.

[0160] Step S524: Identify and correct the extreme deviation components in the multidimensional deviation matrix. Mark the deviation components that exceed the normal fluctuation range using the deviation degree detection algorithm, and replace them with the statistical reference values ​​of the deviation components in the same dimension. The statistical reference values ​​are determined based on the historical distribution of deviation components.

[0161] Extreme deviation components are those that deviate significantly from the normal fluctuation range in a multidimensional deviation matrix. These extreme deviation components may be caused by data noise, outliers, or model errors. If left unaddressed, these extreme deviation components may mislead subsequent analysis and decision-making.

[0162] Deviation components that exceed the normal fluctuation range are identified using a deviation detection algorithm. This algorithm can determine the normal fluctuation range based on the distribution of historical deviation components. For example, statistical methods can be used, such as calculating the mean and standard deviation of historical deviation components, and marking deviation components that exceed the mean plus or minus a certain multiple of the standard deviation as extreme deviation components.

[0163] Replace the extreme deviation components with statistical reference values ​​for the deviation components in the same dimension. These statistical reference values ​​are determined based on the historical distribution of the deviation components; for example, the median, mode, or other statistical measures of the historical deviation components can be used. This reduces the impact of extreme deviation components on the overall analysis, making the deviation series more stable and reliable.

[0164] Step S525: Perform feature fusion on the corrected multidimensional deviation matrix, and combine multiple deviation components into a comprehensive deviation value by weighted averaging. The weights are dynamically adjusted according to the importance of the quantile values ​​to generate a comprehensive deviation sequence that includes the time dimension.

[0165] Multiple deviation components in the corrected multidimensional deviation matrix are combined using a weighted average method. The weights of the weighted average are dynamically adjusted based on the importance of the quantile values. Different quantile values ​​may have different importance in practical applications; for example, the median may better represent the central trend of the data, while the 90th quantile may focus more on the upper limit of the data. Different weights can be assigned to each quantile value according to specific needs and application scenarios. A comprehensive deviation sequence containing a time dimension is generated. The comprehensive deviation sequence is a one-dimensional sequence, where each element corresponds to the comprehensive deviation value at a specific time point.

[0166] Step S526: Associate the comprehensive deviation sequence with the time markers to obtain a deviation sequence containing time nodes and corresponding deviation values.

[0167] Associating the composite deviation series with time markers clarifies the time node corresponding to each composite deviation value. This imbues the deviation series with temporal information, facilitating time series analysis and trend assessment. The composite deviation series and time markers can be combined into a two-dimensional array or data structure, where one column represents the time node and the other column represents the corresponding composite deviation value.

[0168] Step S530: Perform cumulative trend analysis on the deviation sequence, calculate the cumulative deviation value by weighted summation through a sliding window, and start the parameter adjustment process when the cumulative deviation value exceeds the preset trigger condition; otherwise, maintain the current configuration of the initial risk predictor.

[0169] The preset trigger condition is a critical value for accumulated deviation, pre-set based on extensive experiments and experience. When the accumulated deviation exceeds the preset trigger condition, it indicates that the model's prediction error has reached an unacceptable level, requiring the parameter adjustment process to be initiated. The parameter adjustment process can adjust the parameters of the initial risk predictor to improve the model's prediction accuracy. If the accumulated deviation does not exceed the preset trigger condition, it means that the model's current configuration is still well adapted to the actual situation, and the current configuration of the initial risk predictor can be maintained.

[0170] Step S540: After the parameter adjustment process is started, an optimization objective function is constructed with the deviation sequence as input. The optimization objective function includes the cumulative deviation value and the model parameter regularization term. The network parameters of the initial risk predictor are updated through the gradient descent algorithm, and the parameter update direction is consistent with the deviation direction.

[0171] The optimization objective function includes the cumulative bias and the model parameter regularization term. The cumulative bias reflects the model's prediction error and is the primary objective of optimization. The model parameter regularization term is introduced to prevent overfitting. Regularization terms constrain the magnitude of model parameters, preventing them from becoming too large and improving the model's generalization ability. Examples of regularization terms include L1 and L2 regularization. L1 regularization makes the model parameters sparser, while L2 regularization makes them smoother.

[0172] The network parameters of the initial risk predictor are updated using the gradient descent algorithm. Gradient descent is a commonly used optimization algorithm that calculates the gradient of the objective function with respect to the network parameters and then updates the parameters in the opposite direction of the gradient. The direction of parameter updates is consistent with the direction of bias, i.e., adjusting the parameters in the direction of reducing the accumulated bias. In each iteration, the values ​​of the network parameters are updated based on the magnitude of the gradient and the learning rate. Through continuous iteration, the value of the objective function is gradually reduced, improving the model's prediction accuracy.

[0173] Step S550: Based on the parameter adjustment magnitude, identify the network layer that contributes significantly to the bias, trace the causal path corresponding to the input features of the network layer, and adjust the intensity quantization value of the corresponding causal path according to the parameter adjustment magnitude ratio.

[0174] In one implementation, step S550 may specifically include the following steps S551 to S556: Step S551: Based on the difference between the updated parameters and the parameters before the update, obtain the parameter adjustment range, calculate the comprehensive value of the parameter adjustment range of each network layer, the comprehensive value reflects the overall influence of the layer's parameters on the deviation, and mark the network layers with the highest comprehensive values ​​as deviation-sensitive layers.

[0175] The parameter adjustment magnitude is obtained based on the difference between the updated and unupdated parameters. During parameter adjustment, the parameter values ​​before and after the update for each network layer are recorded, and the difference between them is calculated. This difference reflects the adjustment magnitude of the network layer's parameters in one iteration. A comprehensive value of the adjustment magnitude of each network layer's parameters is calculated. This value can be the sum of the absolute values ​​of the parameter adjustment magnitudes, the sum of squares, or other statistical measures; the specific calculation method can be chosen according to the actual situation. The comprehensive value reflects the overall influence of the network layer's parameters on bias. The larger the comprehensive value, the greater the contribution of the network layer's parameter adjustment to reducing the cumulative bias. Network layers with the highest comprehensive values ​​are marked as bias-sensitive layers. Bias-sensitive layers are network layers that significantly contribute to bias; the parameter adjustments of these network layers have a significant impact on the model's predictive performance. By marking bias-sensitive layers, attention can be focused on these key network layers to further analyze their input features and corresponding causal paths.

[0176] Step S552: Analyze the source of input features of the bias-sensitive layer, determine the causal path identifier corresponding to the feature by reverse tracing the feature propagation path, and establish a mapping relationship table between network layers and causal paths. The mapping relationship table contains the sensitive layer number and the corresponding causal path set.

[0177] Analyze the sources of input features for the bias-sensitive layer to understand from which data sources or intermediate layers these features are passed. The feature propagation path is the process by which features are transferred within the model. By tracing the feature propagation path in reverse, the causal path identifier corresponding to the input features can be determined. The causal path identifier is a unique number or name that identifies each causal path.

[0178] Establish a mapping table between network layers and causal paths. Each row of the mapping table corresponds to a bias-sensitive layer, and each column corresponds to a causal path. The mapping table clearly shows which causal paths each bias-sensitive layer is associated with. The table includes the sensitive layer number and the corresponding set of causal paths, which helps in subsequently adjusting the intensity quantization value of the corresponding causal path according to the parameter adjustment amplitude ratio.

[0179] Step S553: ​​Based on the parameter adjustment range and mapping relationship table, calculate the contribution ratio of each causal path to the deviation. The contribution ratio is the ratio of the parameter adjustment range of the corresponding feature of the path to the total adjustment range. Sort the paths by contribution ratio from largest to smallest to generate a path adjustment priority list.

[0180] Based on the parameter adjustment magnitude and mapping table, the contribution ratio of each causal path to the bias is calculated. For each causal path, the bias-sensitive layer containing its corresponding feature is identified, and the parameter adjustment magnitude of that network layer is obtained. The ratio of the parameter adjustment magnitude of the path's corresponding feature to the total adjustment magnitude is taken as the contribution ratio of that causal path to the bias. The total adjustment magnitude is the sum of the parameter adjustment magnitudes of all bias-sensitive layers. A path adjustment priority list is generated by sorting the paths from largest to smallest contribution ratio. The path adjustment priority list helps determine which causal paths need to have their intensity quantization values ​​adjusted first. Causal paths with larger contribution ratios have a greater impact on the bias and should be adjusted first.

[0181] Step S554: Adjust the intensity quantization value of the causal path in order according to the priority list. The adjustment amount is the product of the original intensity value and the corresponding contribution ratio. The adjustment direction is consistent with the parameter adjustment direction.

[0182] The intensity quantization values ​​of causal paths are adjusted sequentially according to a priority list. Starting with the causal path with the largest contribution, the intensity quantization values ​​are adjusted in turn. The adjustment amount is the product of the original intensity value and the corresponding contribution proportion. The adjustment direction is consistent with the parameter adjustment direction. If the parameter adjustment aims to reduce the cumulative bias, then the adjustment of the causal path intensity quantization values ​​should also aim to make the model more accurately reflect causal relationships. By adjusting the intensity quantization values ​​of causal paths sequentially according to priority, the set of causal path weights can more accurately reflect the causal relationships between variables, improving the model's interpretability and predictive performance.

[0183] Step S555: Perform logical verification on the adjusted causal path weight set, check the correlation and coordination between path strength values, and make fine adjustments if there are logical conflicts so that the weight set as a whole conforms to the transmission characteristics of causal relationships.

[0184] Logical verification of the adjusted causal path weight set is necessary to ensure its rationality and consistency. Causal relationships have a transitive property: if A leads to B, and B leads to C, then A should indirectly lead to C through B. After adjusting the quantified values ​​of causal path strength, it is necessary to check the correlation and consistency between the path strength values ​​to ensure they conform to the transitive property of causality. Checking the correlation and consistency between path strength values ​​can be done in several ways. For example, it can be checked whether the strength values ​​of causal paths are non-negative and whether there are unreasonable negative strength values. It can also be checked whether the relative strengths between causal paths are reasonable; for example, if the strength value of a direct causal path is much smaller than the strength value of its indirect causal path, there may be a logical conflict. If a logical conflict exists, fine-tuning is performed. Fine-tuning can be done by appropriately adjusting the quantified values ​​of the relevant causal paths based on the logic of the causal relationship and the actual situation. For example, the strength values ​​of some causal paths can be appropriately increased or decreased to conform to the transitive property of causality. Through logical verification and fine-tuning, the overall causal path weight set becomes more reasonable and reliable.

[0185] Step S556: Store the verified causal path weight set into the preset parameter library, overwrite the original weight set, and generate a weight update record, which includes the change in intensity value before and after adjustment and the corresponding timestamp.

[0186] Step S560: Input the updated network parameters and causal path weight set into the initial risk predictor to generate a sales volume risk assessment result that includes risk level classification and fluctuation range description. The risk level classification is determined based on the quantile value interval of the probability distribution sequence, and the fluctuation range description reflects the sales volume fluctuation range corresponding to different risk levels.

[0187] The updated network parameters and causal path weights are input into the initial risk predictor, enabling the model to make predictions using the latest parameters and causal relationship information. The updated network parameters improve the model's predictive accuracy, while the updated causal path weights allow the model to better understand the causal relationships between variables, enhancing the model's interpretability.

[0188] Generate a risk assessment result for electricity sales that includes risk level classification and fluctuation range description. Risk level classification is determined based on the quantile intervals of the probability distribution sequence. For example, different quantile intervals of the probability distribution sequence can be divided into different risk levels, such as low risk (e.g., 0-20% quantile interval), medium risk (20%-80% quantile interval), and high risk (80%-100% quantile interval). The fluctuation range description reflects the electricity sales fluctuation range corresponding to different risk levels; for example, the electricity sales fluctuation range corresponding to a low-risk level may be smaller, while the electricity sales fluctuation range corresponding to a high-risk level may be larger.

[0189] Based on the foregoing embodiments, this invention provides a power sales risk prediction device. The various units and modules included in the device can be implemented by a processor in a computer device; of course, they can also be implemented by specific logic circuits. In the implementation process, the processor can be a central processing unit (CPU), a microprocessor unit (MPU), a digital signal processor (DSP), or a field programmable gate array (FPGA), etc.

[0190] Figure 2 This is a schematic diagram of the composition structure of a power sales risk prediction device provided in an embodiment of the present invention, as shown below. Figure 2 As shown, the electricity sales risk prediction device 200 includes: The data acquisition module 210 is used to acquire historical electricity sales records and related power grid operation observation data, and to establish a sample dataset containing time stamps. The sample dataset has time continuity and variable correlation. Causal identification module 220 is used to identify causal relationships in sample datasets, analyze direct and indirect causal relationships between variables through structural causal model, and generate a set of causal path weights containing causal strength quantification values. The vector construction module 230 is used to fuse the causal path weight set with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logic constraints. The parameter fine-tuning module 240 is used to input the time series representation vector into the meta-learning framework model, generate an initial risk predictor through parameter fine-tuning, and output the probability distribution sequence of electricity sales fluctuations. The risk assessment module 250 is used to obtain the deviation value between the real-time electricity sales observation value and the probability distribution sequence, adjust the parameter configuration of the initial risk predictor based on the deviation value, update the causal strength quantification value in the causal path weight set, and generate the electricity sales risk assessment result.

[0191] The descriptions of the apparatus embodiments above are similar to those of the method embodiments above, and have similar beneficial effects. In some embodiments, the functions or modules included in the apparatus provided by the present invention can be used to perform the methods described in the method embodiments above. For technical details not disclosed in the apparatus embodiments of the present invention, please refer to the descriptions of the method embodiments of the present invention for understanding.

[0192] Figure 3A hardware entity diagram of a computer system provided as an embodiment of the present invention, such as... Figure 3 As shown, the hardware entity of the computer system 1000 includes a processor 1001 and a memory 1002, wherein the memory 1002 stores a computer program that can run on the processor 1001, and the processor 1001 executes the program to implement the steps in the method of any of the above embodiments.

Claims

1. A method for predicting electricity sales risk based on causal explanation and online self-learning, characterized in that, The method includes: Historical electricity sales records and associated power grid operation observation data are obtained to establish a sample dataset containing time stamps, wherein the sample dataset has temporal continuity and variable correlation. The sample dataset is used to identify causal relationships. Direct and indirect causal relationships between variables are analyzed using a structural causal model to generate a set of causal path weights containing quantified causal strength values. The causal path weight set is fused with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logic constraints. The time series representation vector is input into the meta-learning framework model, and an initial risk predictor is generated through parameter fine-tuning. The initial risk predictor outputs a probability distribution sequence of electricity sales fluctuations. Obtain the deviation value between the real-time electricity sales observation value and the probability distribution sequence, adjust the parameter configuration of the initial risk predictor based on the deviation value, update the causal strength quantification value in the causal path weight set, and generate the electricity sales risk assessment result.

2. The method according to claim 1, characterized in that, The step involves identifying causal relationships in the sample dataset, analyzing direct and indirect causal associations between variables using a structural causal model, and generating a set of causal path weights containing quantified causal strength values, including: The sample dataset is decomposed into a multi-scale time series, and the intraday fluctuation component, the intraweek trend component, and the intramonth cycle component are separated according to the nested relationship of time period to obtain a multi-dimensional time series set with hierarchical time structure, and each component maintains the time label correspondence of the original data. Time-varying analysis of variable association is performed on the multi-dimensional time series set. The dynamic association degree between variables is calculated within a sliding time window, and a time-varying curve of association strength is generated as the time window moves. The time-varying curve of association strength represents the changing trend of variable association in different time periods. Based on the time-varying curve of the correlation strength, the causal path is hierarchically divided, and the direct causal path and indirect causal path are distinguished by path length analysis. The indirect causal path is modeled for transmission effect attenuation. The transmission efficiency coefficient of each intermediate node in the path is calculated. The transmission efficiency coefficient decreases exponentially with the increase of path length, and an initial value of indirect causal strength containing attenuation factor is generated. The dynamic backdoor criterion is used to remove interference factors from the hierarchical causal path, identify and block non-causal paths in different time windows, and retain the valid causal path that meets the time sequence constraint. The intensity quantification values ​​of direct causal paths and the attenuated intensity values ​​of indirect causal paths are integrated and weighted by time window weights to generate a set of causal path weights that dynamically adjust over time.

3. The method according to claim 2, characterized in that, The sample dataset is decomposed into multi-scale time series components, separating intraday fluctuation components, weekly trend components, and monthly cycle components according to the nested relationship of time periods, resulting in a multi-dimensional time series set with hierarchical time structure, including: The historical electricity sales records in the sample dataset are adaptively decomposed to generate multiple intrinsic mode function components. Each component is arranged from high to low frequency to characterize the fluctuation characteristics at different time scales. The inherent mode function components are identified by periodicity characteristics. The principal period length of each component is calculated by autocorrelation function. Components with principal period lengths in different ranges are classified as intraday fluctuation components, intraweek trend components, and intramonth periodic components, respectively. Calculate the energy percentage of each periodic component. The energy percentage is the ratio of the component variance to the original sequence variance. Components with an energy percentage exceeding a preset ratio are retained, while noise components with an energy percentage below a preset threshold are removed. The intraday fluctuation component, weekly trend component, and monthly cycle component after classification are aligned with time markers so that each component has a corresponding observation value at the same time node, thus constructing a multi-dimensional time series set containing three time scales. The interaction between different periodic components is analyzed, the modulation coefficient of intraday fluctuation component on intraweek trend component is calculated, and the modulation coefficient is incorporated as an additional feature into a multi-dimensional time series set. Map each periodic component to a unified time axis to give each component in a multi-dimensional time series set the same time resolution. The hierarchical division of causal paths based on the time-varying curve of the correlation strength, and the differentiation of direct causal paths from indirect causal paths through path length analysis, includes: Construct an adjacency matrix for the variable association network. The matrix elements represent the maximum values ​​of the time-varying curves of the association strength between corresponding variable pairs. The larger the value, the stronger the association between the variables. Traverse the adjacency matrix, starting from each variable node, explore other reachable variable nodes, and record all possible variable connection paths and their corresponding path lengths. The path length is the number of edges contained in the path. Obtain a preset path length threshold, mark connection paths with a path length equal to the path length threshold as direct causal path candidate set, and mark path paths with a path length greater than the path length threshold as indirect causal path candidate set, thus initially distinguishing between the two types of paths; Perform time sequence verification on the candidate set of direct causal paths, determine the causal direction between variables, and retain paths where the timestamp of the cause variable is earlier than the timestamp of the result variable; The validity of intermediate nodes in the candidate set of indirect causal paths is verified, and the conditional correlation degree of intermediate nodes is calculated. If the correlation strength of the original variable pair decreases by more than a preset proportion after removing an intermediate node, the intermediate node is determined to be a valid transit node, and the indirect causal path containing the valid transit node is retained. Integrate and verify direct and indirect causal paths to construct a hierarchical causal path list that includes path type, path length, and correlation strength.

4. The method according to claim 1, characterized in that, The step of fusing the causal path weight set with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logical constraints includes: Multimodal feature extraction is performed on the time series pattern features in the sample dataset. Voltage stability features, load distribution features, and line loss features are extracted from the power grid operation observation data. Together with the time series pattern features of historical electricity sales records, they constitute a multimodal input feature set. Based on the causal path weight set, the multimodal input feature set is subjected to causal correlation screening, retaining features that have a direct causal correlation with electricity sales risk prediction and eliminating redundant features that have no causal correlation. The selected multimodal features are aligned at the time granularity, and features with different sampling frequencies are uniformly converted into feature sequences with the same time interval, and the time interval is consistent with the time stamp interval of the sample dataset. The aligned multimodal features and causal path weights are input into the evidence layer of the Bayesian network, and the feature states of each time node are activated sequentially in chronological order. The causal path weights constrain the direction of probability dependence between nodes. Probabilistic inference is performed on the activated Bayesian network, and the posterior probability distribution of the latent variables is calculated at each time point. The posterior probability distribution comprehensively represents the joint influence of multimodal features and causal paths, and generates an intermediate feature vector containing probability distribution parameters. The intermediate feature vector is subjected to temporal correlation enhancement processing. The dependency relationship between features at different time nodes is captured by a temporal convolutional network to generate a temporal representation vector that integrates causal logic and temporal correlation. Each dimension of the temporal representation vector corresponds to the comprehensive feature state of the corresponding time node.

5. The method according to claim 4, characterized in that, The step of aligning the selected multimodal features at the time granularity, converting features with different sampling frequencies into feature sequences with the same time interval, and ensuring that the time interval is consistent with the time stamp interval of the sample dataset, includes: Temporal sampling characteristics analysis was performed on the selected multimodal features. The original sampling interval and sampling point distribution density of each feature were statistically analyzed. A feature-sampling characteristic correspondence table containing feature identifier, sampling interval duration and data integrity index was generated. The data integrity index represents the proportion of missing observations of the feature at each time node. Based on the feature-sampling characteristic correspondence table, multimodal features are divided into two categories: short-interval sampling features and long-interval sampling features. The sampling interval of short-interval sampling features is less than the target time interval, while the sampling interval of long-interval sampling features is greater than or equal to the target time interval, thus realizing the classification of the temporal characteristics of features. The short-interval sampling features are compressed at the time granularity. The sliding window trend extraction algorithm is used to calculate the feature fluctuation trend parameters within each target time interval. The trend parameters include the slope of the rising segment, the slope of the falling segment, and the fluctuation amplitude. The feature extreme points within the window are retained to generate a compressed feature sequence containing trend parameters and extreme points. The time interval of the compressed feature sequence is consistent with the target time interval. The long-interval sampling features are extended at the time granularity. A trend extrapolation model is constructed based on the feature value sequence and sampling interval of historical sampling points. The trend extrapolation model is used to extrapolate the time interval between adjacent historical sampling points to generate a preliminary estimate of the intermediate time feature under the target time interval. The fluctuation compensation amount is calculated based on the historical data of the long-interval sampling features. The fluctuation compensation amount is determined by analyzing the fluctuation pattern between historical sampling points. The fluctuation compensation amount is superimposed on the preliminary estimate generated by the trend extrapolation to generate an extended feature sequence containing trend components and fluctuation components. Time stamp alignment verification is performed on compressed and expanded feature sequences. Based on the time stamp of the sample dataset, it is checked whether each feature sequence has a valid feature value at the corresponding time node. Missing time nodes are filled by a weighted average of the feature values ​​of the preceding and following time nodes. The weight is inversely proportional to the time distance. The time synchronization error of each feature sequence after alignment is calculated. The time synchronization error is the sum of squares of the deviations between the actual time stamp and the target time stamp of the feature sequence. If the error exceeds the preset threshold, the trend extraction window size or extrapolation model parameters are readjusted until all feature sequences meet the time synchronization requirements, generating a set of multimodal feature sequences with uniform time granularity.

6. The method according to claim 4, characterized in that, The process involves performing probabilistic inference on the activated Bayesian network using a variational inference algorithm, calculating the posterior probability distribution of latent variables at each time point. This posterior probability distribution comprehensively represents the joint influence of multimodal features and causal paths, generating an intermediate feature vector containing probability distribution parameters, including: Initialize the variational distribution parameters of the latent variables in the Bayesian network, and determine the initial values ​​of the mean vector and covariance matrix based on the prior probability distribution of the latent variables. The mean vector represents the initial estimate of the central tendency of the latent variables, and the covariance matrix represents the initial correlation between the dimensions of the latent variables. At each time point, the likelihood function value is calculated based on the aligned multimodal feature sequence and the causal path weight set. The likelihood function value represents the probability density of the corresponding value of the latent variable under the current feature state. The contribution ratio of different modal features to the likelihood function is adjusted by the causal path weight. A variational lower bound function is constructed as the optimization objective. The variational lower bound function contains the expectation term of the likelihood function and the entropy term of the variational distribution of the latent variables. The approximate estimation of the posterior probability distribution is completed by maximizing the variational lower bound function. The variational distribution parameters are updated iteratively. In each iteration, time window samples are randomly sampled from the multimodal feature sequence. The gradient of the variational lower bound function with respect to the mean vector and covariance matrix is ​​calculated. The parameter values ​​are adjusted according to the gradient direction until the function value converges. Once the variational lower bound function converges, the diagonal elements of the mean vector and covariance matrix of the variational distribution are extracted. The elements of the mean vector represent the optimal estimates of each dimension of the latent variables, and the diagonal elements of the covariance matrix represent the degree of uncertainty of the estimates of each dimension. By concatenating the mean vector and covariance diagonal elements of each time point in chronological order, a feature matrix containing the time dimension is obtained. The rows of the matrix correspond to the time points, and the columns correspond to the distribution parameters of the latent variables.

7. The method according to claim 1, characterized in that, The step involves inputting the time-series representation vector into the meta-learning framework model, generating an initial risk predictor through parameter fine-tuning, and the initial risk predictor outputting a probability distribution sequence of electricity sales fluctuations, including: A meta-learning framework model is constructed, which includes a feature encoding module, a meta-knowledge storage module, a task adaptation module, and a prediction output module. The feature encoding module is used for nonlinear transformation of temporal representation vectors, the meta-knowledge storage module stores empirical parameters of historical prediction tasks, the task adaptation module realizes cross-task knowledge transfer, and the prediction output module generates probability distribution sequences. The temporal representation vector is input into the feature encoding module, and temporal features are extracted deeply through a stacked long short-term memory network. The hidden state of the long short-term memory network contains feature information at different time scales, generating a high-dimensional abstract feature vector. The high-dimensional abstract feature vector is input into the meta-knowledge storage module and dynamically similar to the stored historical task meta-knowledge to extract the meta-parameters of the most similar historical task. The network parameters of the task adaptation module are initialized based on the extracted meta-parameters. The task adaptation module adjusts the weight distribution of the meta-parameters through an attention mechanism. The weight distribution represents the applicability of different historical meta-knowledge to the current prediction task. The parameters of the task adaptation module are fine-tuned through an inner loop optimization algorithm. The inner loop takes the prediction loss of the current task as the optimization target and updates the network parameters through gradient descent. The parameters of the feature encoding module and the meta-knowledge storage module are fixed during the learning process. After fine-tuning, the output features of the task adaptation module are input into the prediction output module. The probability distribution sequence containing the probability of different fluctuation ranges is generated by probability distribution fitting, thus completing the construction of the initial risk predictor.

8. The method according to claim 1, characterized in that, The process of obtaining the deviation between the real-time electricity sales observation and the probability distribution sequence, adjusting the parameter configuration of the initial risk predictor based on the deviation, updating the causal strength quantification value in the causal path weight set, and generating an electricity sales risk assessment result includes: Continuous electricity sales observations are received through a real-time data interface. The observations are divided into equal-length data segments according to time stamps. Each segment contains observations from multiple consecutive time nodes, resulting in real-time data blocks that match the time granularity of the probability distribution sequence. The time stamps of the data blocks correspond one-to-one with the time nodes of the probability distribution sequence. Align the real-time data blocks with the probability distribution sequence on the time axis, calculate the absolute difference between the observation value at each time node and multiple quantile values ​​of the probability distribution at the corresponding node, and generate a deviation sequence containing multiple deviation components, with the number of deviation components being the same as the number of quantile values. Perform cumulative trend analysis on the deviation sequence, calculate the cumulative deviation value by weighted summation through a sliding window, and start the parameter adjustment process when the cumulative deviation value exceeds the preset trigger condition; otherwise, maintain the current configuration of the initial risk predictor. After the parameter adjustment process is started, the optimization objective function is constructed with the deviation sequence as input. The optimization objective function includes the cumulative deviation value and the model parameter regularization term. The network parameters of the initial risk predictor are updated through the gradient descent algorithm, and the parameter update direction is consistent with the deviation direction. Based on the parameter adjustment magnitude, identify the network layers that significantly contribute to the bias, trace the causal path corresponding to the input features of the network layer, and adjust the intensity quantization value of the corresponding causal path according to the parameter adjustment magnitude ratio. The updated network parameters and causal path weights are input into the initial risk predictor to generate a sales volume risk assessment result that includes risk level classification and fluctuation range description. The risk level classification is determined based on the quantile value interval of the probability distribution sequence, and the fluctuation range description reflects the sales volume fluctuation range corresponding to different risk levels.

9. A power sales risk prediction device, characterized in that, include: The data acquisition module is used to acquire historical electricity sales records and associated power grid operation observation data, and to establish a sample dataset containing time stamps, wherein the sample dataset has time continuity and variable correlation. The causal identification module is used to identify causal relationships in the sample dataset, analyze the direct and indirect causal relationships between variables through a structural causal model, and generate a set of causal path weights containing quantified causal strength values. The vector construction module is used to fuse the causal path weight set with the time series pattern features in the sample dataset using a Bayesian network to construct a time series representation vector with causal logic constraints. The parameter fine-tuning module is used to input the time series representation vector into the meta-learning framework model, generate an initial risk predictor through parameter fine-tuning, and the initial risk predictor outputs a probability distribution sequence of electricity sales fluctuations. The risk assessment module is used to obtain the deviation value between the real-time electricity sales observation value and the probability distribution sequence, adjust the parameter configuration of the initial risk predictor based on the deviation value, update the causal strength quantification value in the causal path weight set, and generate the electricity sales risk assessment result.

10. A computer system comprising a memory and a processor, the memory storing a computer program executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method according to any one of claims 1 to 8.

Images