IEGS probabilistic energy flow monitoring method based on sparse arbitrary chaotic polynomial model

Abstract

This invention relates to the field of IEGS probabilistic energy flow technology, and particularly to an IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model, comprising: S1, establishing a deterministic energy flow model for an integrated energy system (IEGS); S2, modeling random input variables and configuring corresponding non-Gaussian probability density functions; S3, generating P multidimensional orthogonal aPC basis functions; S4, generating M0 initial collocation points; S5, transforming the problem of solving the sparse aPC coefficient vector A into an l1-l2 norm minimization optimization model; solving for the sparse aPC coefficient vector A under the current collocation point set; and constructing a sparse aPC surrogate model; S6, verification; outputting the final sparse aPC surrogate model; S7, analytically calculating the statistical moments of key output quantities and reconstructing the probability density functions of key output quantities; S8, evaluating the safety margin of IEGS operation. This method can effectively reduce model dimensionality and improve computational efficiency by combining sparse recovery techniques while retaining the good adaptability of the aPC method to non-Gaussian inputs.

Description

Technical Field

[0001] This invention relates to the field of IEGS probabilistic energy flow technology, and in particular to an IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model. Background Technology

[0002] Energy is a crucial cornerstone of global sustainable development. With continuous socio-economic progress and sustained growth in energy demand, the limitations of traditional single-energy systems are becoming increasingly apparent. Historically, power systems and natural gas pipeline networks have long been operated by independent entities, lacking effective coordination and coupling. This has led to problems such as low energy resource utilization efficiency, insufficient dispatch flexibility, and weak capacity to cope with fluctuating loads. In recent years, as the significant advantages of gas-fired power generation in terms of economics and environmental performance have become increasingly apparent, its penetration rate in power systems has been continuously increasing, driving the development of Integrated Electricity-Gas Systems (IEGS). IEGS, through deep integration of the power grid and natural gas network, achieves multi-energy complementarity and synergistic optimization, demonstrating enormous potential in enhancing the absorption capacity of renewable energy sources (RES), improving system operational flexibility, and increasing overall energy utilization efficiency.

[0003] However, the complexity of IEGS also brings new challenges. Its operation involves numerous uncertainties, primarily including random fluctuations in electricity and natural gas loads, the highly intermittent and unpredictable output of renewable energy sources such as wind and solar power, and the impact of external environmental changes on equipment operation. These uncertainties cause system state variables (such as node voltage, branch power flow, and pipeline pressure) to exhibit complex probability distribution characteristics, making it difficult for traditional deterministic energy flow analysis methods to accurately reflect the system's actual operating conditions. Therefore, effectively quantifying the impact of uncertainties on IEGS operation has become a crucial prerequisite for achieving safe, stable, and efficient system operation.

[0004] Probabilistic Energy Flow (PEF), as an important tool for quantifying uncertainty, is widely used to evaluate the probabilistic characteristics of key outputs (Quantities of Interest, QoIs) of a system. PEF methods are mainly divided into two categories: one is the random sampling method represented by Monte Carlo Simulation (MCS), and the other is the model-driven method based on mathematical analysis. While MCS can provide high-precision probability density functions (PDFs) or cumulative distribution functions (CDFs), its computational cost is extremely high, especially in high-dimensional input spaces, where the required sample size grows exponentially, making it difficult to meet the needs of real-time or online applications. In contrast, analytical methods such as the method of moments (MoM) and point estimation methods require only a small number of samples to estimate the statistical moments of QoIs and can approximate the PDF using Edgeworth or Cornish-Fisher series expansions. However, these methods typically assume that the output follows a Gaussian or near-Gaussian distribution; when the actual distribution exhibits significant skewness, kurtosis anomalies, or multimodality, their approximation accuracy decreases significantly.

[0005] To overcome the aforementioned limitations, Polynomial Chaos Expansion (PCE) has been introduced into PEF analysis. PCE represents the system response as a linear combination of orthogonal polynomial basis functions of the input random variables, thereby achieving efficient modeling from the propagation of input uncertainty to the statistical properties of the output. However, traditional PCE requires all input variables to be converted to standard Gaussian variables first, which has inherent limitations when dealing with non-Gaussian distributed inputs, restricting its applicability in real-world scenarios. Therefore, Arbitrary Polynomial Chaos (aPC) has emerged. It no longer relies on distributional assumptions about variables, but directly constructs non-Gaussian probability density functions based on the statistical properties of historical data or prediction errors, making the model more adaptable to real-world scenarios.

[0006] Although aPC has made breakthroughs in dealing with non-Gaussian uncertainties, it still faces severe challenges in the highly complex multi-source coupled system of IEGS. First, the large number of random variables involved in IEGS, covering multiple dimensions such as power load, natural gas demand, wind speed, light intensity, etc., leads to an exponential growth of the number of polynomial terms in the aPC model with the number of variables, causing the "curse of dimensionality" problem. Second, traditional aPC methods require a large number of collocation points for solution, resulting in a heavy computational burden and difficulty in meeting the requirements of engineering practice for computational efficiency. Although sparse aPC techniques (such as partial least squares, analysis of variance, least angle regression, etc.) can reduce the model complexity by eliminating low-contribution terms, there is still a lack of an effective means that can both maintain high accuracy and significantly reduce the sample requirements.

[0007] Compressive Sensing (CS) theory provides a new idea for solving this problem. CS can recover the original signal with far fewer samples than required by the Nyquist sampling theorem by utilizing the sparsity of the signal. In the aPC framework, if the coefficient vector has a sparse structure, efficient modeling can be achieved by solving an underdetermined linear system. Among them, l0 norm minimization is the most ideal sparse recovery criterion, but its practical application is limited due to its NP-hard nature. Therefore, researchers proposed l1 minimization as a convex relaxation alternative, but its sparse recovery ability is limited; furthermore, non-convex optimization models such as l p (0 < p < 1) and l1-l2 norm minimization were proposed. The latter shows a stronger sparse promotion ability because its contour is closer to the l0 norm and has become a current research hotspot.

[0008] However, although existing technologies have alleviated the problem of uncertainty quantification in IEGS to a certain extent, in the face of the contradiction between non-Gaussian inputs, high-dimensional random variables, and computational efficiency, there is still a lack of a unified solution that combines high accuracy, strong robustness, and low computational cost. Especially in IEGS, due to the strong system coupling, high variable dimensions, and complex distributions, traditional aPC methods are difficult to balance model accuracy and computational feasibility. At the same time, the stability and convergence of sparse recovery algorithms still need to be improved when the measurement matrix conditions do not meet the ideal assumptions.

[0009] Therefore, how to effectively reduce the model dimension and improve the computational efficiency by combining sparse recovery techniques while retaining the good adaptability of aPC methods to non-Gaussian inputs, so as to achieve high-precision, low-cost, and real-time analysis of the probabilistic energy flow in IEGS, has become an urgent problem to be solved. Summary of the Invention

[0010] Aiming at the above deficiencies in the existing technologies, the purpose of the present invention is to provide a method for monitoring the probabilistic energy flow of IEGS based on a sparse arbitrary chaos polynomial model, which can, while retaining the good adaptability of the aPC method to non-Gaussian inputs, effectively reduce the model dimension and improve the calculation efficiency by combining sparse recovery techniques, so as to achieve high-precision, low-cost, and real-time analysis of the probabilistic energy flow of IEGS.

[0011] In order to solve the above technical problems, the technical solution adopted by the present invention is as follows:

[0012] The method for monitoring the probabilistic energy flow of IEGS based on a sparse arbitrary chaos polynomial model includes the following steps:

[0013] S1. Establish a deterministic model of the energy flow of the integrated energy system IEGS for solving the system state variables under given input conditions;

[0014] S2. Based on the statistical characteristics of historical operation data or prediction errors, identify the key uncertainty sources affecting the energy flow in the integrated energy system, model the prediction deviation of each key uncertainty source as a random input variable, so as to form a set containing d mutually independent random input variables, and configure a corresponding non-Gaussian probability density function for each random input variable according to its actual prediction error distribution; where d is the number of key uncertainty sources;

[0015] S3. According to the number d of random input variables and the preset maximum expansion order p, use the non-Gaussian probability density functions of the random input variables, and respectively construct orthogonal univariate polynomial bases for each variable through a numerical orthogonalization method based on data-driven sample moments; then generate P multi-dimensional orthogonal aPC basis functions through tensor product operations on the orthogonal univariate polynomial bases;

[0016] S4. Set an initial number of collocation points M0, M0 < P; adopt a preset sampling method to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions;

[0017] S5. Substitute the random input variables corresponding to the current collocation points into the energy flow deterministic model in S1 for solution to obtain the corresponding system state variables, and form a sample vector Y by arranging all the obtained system state variables in columns as the key output quantity; for the i-th collocation point, calculate the basis function values of all P aPC basis functions in S3 at this point, and use this P-dimensional row vector as the i-th row of the measurement matrix Φ, and gather the basis function value row vectors of all collocation points to form a complete measurement matrix Φ; transform the problem of solving the sparse aPC coefficient vector A into the following l1-l2 norm minimization optimization model:

[0018] ;

[0019] In the formula, λ is the regularization parameter, and λ>0;

[0020] Solve the model to obtain the sparse aPC coefficient vector A under the current collocation set; and based on the expression Construct a sparse aPC proxy model; where, Let d represent a random input vector consisting of d random input variables; Denotes the basis function of the i-th multivariable orthogonal polynomial; Represents the expansion coefficients corresponding to the i-th aPC basis function;

[0021] S6. Based on the sparse aPC surrogate model of S5, calculate the relative squared error on the additional set of verification collocations generated in the joint probability space defined by d random input variables and their non-Gaussian probability density functions; if the relative squared error is greater than the preset threshold γ, increase the number of collocations by the preset increment and return to step S5; if the relative squared error is less than or equal to γ, terminate the iteration and output the final sparse aPC surrogate model.

[0022] S7. Utilizing the orthogonality of the aPC basis functions in the final sparse aPC surrogate model of S6, the statistical moments of the key output quantities are analytically calculated through the sparse aPC coefficient vector A, and the probability density function of the key output quantities is reconstructed based on the statistical moments.

[0023] S8, based on the probability density function of the reconstructed key output quantity of S7, evaluates the safety margin of IEGS operation, and generates a risk warning signal when the safety margin is lower than the preset safety threshold, which is used to trigger real-time control or day-ahead scheduling adjustment.

[0024] Compared with the prior art, the present invention has the following advantages:

[0025] 1. This method achieves effective modeling of non-Gaussian input variables, improving the model's real-world adaptability. It employs a non-Gaussian probability density function based on historical operating data or the statistical characteristics of prediction errors, configuring a corresponding probability distribution for each key uncertainty source and constructing an arbitrary chaotic polynomial (aPC) basis function. Unlike traditional PCE methods that rely on Gaussianization transformations, this technique directly utilizes the probability characteristics driven by actual data, avoiding modeling errors caused by biased distribution assumptions. Especially in IEGS, electricity, gas loads, and renewable energy output often exhibit skewed, heavy-tailed, and other non-Gaussian properties. This method can more realistically reflect the statistical behavior of input uncertainties, thereby significantly improving the accuracy of the surrogate model.

[0026] 2. The sparse aPC framework effectively alleviates the "curse of dimensionality" and reduces model complexity. Addressing the issue of a large number of random variables and exponential growth of polynomial terms in the IEGS system, this method introduces the idea of ​​sparse recovery, transforming the aPC coefficient solution problem into an l1-l2 norm minimization optimization model. Compared to traditional full-order expansion methods, this method retains only a few polynomial terms that significantly contribute to the system response, greatly reducing the effective degrees of freedom. Simultaneously, by combining the construction of the measurement matrix and an iterative increment mechanism, the computational scale is further controlled. Compared to conventional aPC or high-dimensional sparse regression methods, this strategy significantly compresses the model dimensionality while maintaining model accuracy, improving scalability.

[0027] 3. An adaptive collocation generation and validation mechanism enhances the robustness and convergence of the model. This method sets initial collocations and dynamically adjusts the number of samples based on the relative squared error on the validation set, forming a closed-loop iterative optimization process. If the model accuracy is insufficient, collocations are automatically added to improve the estimation quality; otherwise, the iteration terminates. This adaptive strategy avoids the resource waste caused by manually setting too many samples and prevents underfitting due to insufficient samples. Compared to fixed sampling strategies or one-time large-sample methods, this mechanism minimizes unnecessary computational overhead while ensuring model accuracy, improving the algorithm's intelligence and practicality.

[0028] 4. This method achieves an efficient mapping from deterministic energy flow models to probabilistic outputs, supporting real-time risk assessment and early warning. It embeds a sparse aPC surrogate model into the IEGS energy flow analysis process, analytically calculating the statistical moments of key outputs and reconstructing their probability density functions to assess the system's operational safety margin. Because the surrogate model possesses rapid evaluation capabilities, it can complete uncertainty propagation analysis across numerous scenarios in a very short time, making it suitable for real-time monitoring and day-ahead scheduling scenarios. Compared to Monte Carlo simulations, which require thousands of simulations to obtain reliable results, this method significantly reduces computation time, enabling probabilistic energy flow analysis to have online application potential.

[0029] In summary, this method can effectively reduce model dimensionality and improve computational efficiency by combining sparse recovery techniques while retaining the good adaptability of the aPC method to non-Gaussian inputs, thereby achieving high-precision, low-cost, and real-time analysis of IEGS probabilistic energy flow.

[0030] Preferably, in S5, the alternating direction multiplier method is used to solve the l1-l2 norm minimization optimization model, and the process includes:

[0031] Will Rewritten as G(A)-H(A), where G(A) and H(A) are convex:

[0032] ;

[0033] Linearize H(A) as follows, where A (n) ≠0:

[0034] ;

[0035] Will Simplified to:

[0036] ;

[0037] In the formula, z = Φ T Y +λA (n) / ||A (n) ||2;

[0038] Introduce an auxiliary variable B to decouple the non-smooth l1 norm sum. The smooth portion; the constraint A = B is enforced by an augmented Lagrange formula, where B handles the l1 penalty term, while A manages the remaining terms; the augmented Lagrange formula is as follows:

[0039] ;

[0040] right Taking the partial derivatives with respect to A, B, and u, we get:

[0041] ;

[0042] After giving the initial values, set the three partial derivatives to 0, and then iteratively solve A, B and u until the error is less than the preset error threshold; thus obtaining the sparse aPC coefficient vector A under the current collocation set.

[0043] This approach effectively addresses non-smooth optimization problems, improving solution stability and convergence. The original optimization model contains non-smooth l1 and l2 norm terms, making direct solution difficult. This scheme decomposes the objective function into two convex functions, G(A) and H(A), linearizes H(A), and introduces an auxiliary variable B to separate the l1 term, constructing an augmented Lagrangian function. This structured approach avoids the numerical instability caused by directly handling non-smooth terms, making the optimization process smoother and more controllable, significantly improving the algorithm's robustness and convergence performance.

[0044] 2. Achieve efficient iterative solutions and reduce computational complexity. Using the ADMM framework, the complex joint optimization problem is decomposed into three subproblems concerning A, B, and the dual variable u. Each subproblem can be solved analytically or iteratively. For example, updating A involves a regularized quadratic optimization problem, which can be quickly solved analytically or using the conjugate gradient method; updating B only requires a soft thresholding operation. This block-based solution strategy avoids large-scale matrix inversions or global searches, significantly reducing the computational burden of each iteration, making it suitable for real-time applications in high-dimensional IEGS systems.

[0045] Preferably, in the iterative solution of A, B, and u, l1 minimizes the following solution as the initial value A. (1) :

[0046] .

[0047] This setup provides high-quality initial values, accelerates algorithm convergence, and improves solution stability. Since l1 norm minimization is a classic method for sparse recovery, it can effectively approximate the optimal sparse solution under certain conditions (such as restricting isometric properties). This solution is used as the initial value A for the ADMM iteration. (1) This approach allows the optimization process to start from a point close to a real sparse structure, avoiding getting trapped in local poor solutions or prolonged oscillations. Compared to random initialization, this strategy significantly improves iteration efficiency, reduces the number of steps required to reach convergence, and enhances the overall robustness of the algorithm.

[0048] Preferably, in S1, the deterministic energy flow model of IEGS includes the gas flow balance equation in the gas network, the power flow equation in the power grid, the gas flow equation in the gas pipeline, the compressor model, and the gas turbine model; the key outputs include the voltage amplitude of the power grid node, the phase angle, the pressure of the gas network node, and / or the pipeline flow rate.

[0049] This setup, by constructing a high-fidelity deterministic model covering the core components of electro-pneumatic coupling, enhances the integrity and physical realism of the IEGS system modeling, providing a solid and reliable underlying support for subsequent probabilistic energy flow analysis.

[0050] Preferably, the gas flow balance equation in the gas network is:

[0051] ;

[0052] In the formula, This indicates the incoming airflow from the gas source; This indicates the gas flow in natural gas pipeline ij. This indicates the airflow to the compressor between the natural gas pipelines IK; This represents the natural gas load consumed by the gas turbine unit connected to node i; This represents the natural gas load at node i; This indicates the total number of natural gas pipelines in the gas network; Indicates the number of compressors in IEGS;

[0053] The power flow equations are:

[0054] ;

[0055] In the formula, and These represent the active and reactive power generation at bus i, respectively; This indicates the active power consumed by the compressor in the natural gas network; δ represents the voltage magnitude at bus i; ij =δ i -δ j δ i This represents the voltage phase angle at bus i; This represents the voltage amplitude at bus j; This represents the active power load demand at bus i; This represents the total number of buses in the power grid; and Let i and j represent the real and imaginary parts of the nodal admittance matrix between bus i and j, respectively;

[0056] The gas flow equation in a gas pipeline is:

[0057] ;

[0058] In the formula, This represents the natural gas pressure at node i in the natural gas network; This represents the natural gas pressure at node j in the natural gas network; This represents the flow coefficient of pipe ij.

[0059] This setup enables a precise characterization of the core physical processes of the electro-gas coupling system, enhancing the model's engineering applicability. The gas flow balance equations for the gas network consider key aspects such as gas source injection, pipeline transportation, compressor pressurization, and gas consumption by gas turbine units, comprehensively reflecting the supply, demand, and transmission characteristics of natural gas. The power flow equations for the power grid, based on the standard power balance principle, accurately describe the relationship between node voltage and power flow. The gas pipeline flow equations introduce a nonlinear flow model driven by pressure difference, reflecting the physical nature of gas flow in the pipeline. Together, these three constitute a fundamental IEGS model with strong physical consistency, effectively supporting subsequent uncertainty analysis and operational optimization tasks, significantly enhancing the model's applicability and reliability in practical engineering scenarios.

[0060] Preferably, the compressor model is as follows:

[0061] ;

[0062] In the formula, This indicates the active power consumed by the compressor. and These represent the pressures at the compressor inlet and outlet nodes, respectively. K represents the compressor power conversion factor. G Indicates the adiabatic index of natural gas;

[0063] The gas turbine unit model is as follows:

[0064] ;

[0065] In the formula, This indicates the active power generated by the gas turbine unit; , , These are the quadratic, linear, and constant coefficients of the fuel consumption function, obtained by fitting unit operating data, which together characterize the nonlinear relationship between unit power and gas consumption.

[0066] This setup accurately reflects the nonlinear characteristics of the electro-gas energy conversion process, improving the accuracy of system modeling. As a key regulating device in the natural gas pipeline network, the compressor's energy consumption exhibits a strongly nonlinear relationship with its pressure ratio, making traditional linear approximations prone to error accumulation. Furthermore, the fuel consumption of gas turbine units shows a typical quadratic curve characteristic with output variation, exhibiting significant efficiency decline, especially in low-load areas. This solution employs a compressor model and a secondary fuel consumption model based on thermodynamic characteristics, which more realistically reflects the operating behavior of these devices under different operating conditions. This effectively improves the accuracy of the overall energy flow calculation in IEGS, particularly in multi-energy synergistic optimization and uncertainty propagation analysis.

[0067] Preferably, in S3, the process of generating P multidimensional orthogonal aPC basis functions includes:

[0068] (a) Using the observed values ​​of random input variables corresponding to historical operating data or prediction error samples in S2 as driving sample data, calculate the first 2P+1 order statistical moments u of each random input variable. m m = 0, 1, …, 2P; and form a one-dimensional statistical moment matrix. :

[0069] ;

[0070] (b) By Cholesky decomposition This yields the upper triangular matrix. :

[0071] ;

[0072] (c) According to Calculate the recursive parameter b from the elements. k and c k :

[0073] ;

[0074] In the formula, r 0,0 = 1 and r 0,1 = 0; b k and c k This represents the coefficients determined by random input, k = 1, 2, …, P;

[0075] (d) Using parameter b k and c k One-dimensional orthogonal polynomial basis functions are constructed through three recurrence relations. This ensures that the orthogonality condition is satisfied on the domain D:

[0076] ;

[0077] In the formula, For d-dimensional random input vectors The cumulative distribution function;

[0078] (e) Perform steps (a)-(d) above on d random input variables respectively to obtain d sets of univariate orthogonal bases, and then construct P multidimensional orthogonal aPC basis functions through tensor product.

[0079] This setup achieves two key advantages: 1) Data-driven construction of orthogonal basis functions for non-Gaussian inputs, enhancing modeling flexibility and applicability. Traditional aPC methods typically rely on known probability distributions (such as Beta, Gamma, etc.) to predetermine orthogonal polynomial families, making them difficult to apply directly to unknown or complex-distributed input variables. This approach, however, only requires the statistical moments of the input variables (estimated from historical data) to generate orthogonal polynomial bases through Cholesky decomposition and recursion. This allows the method to be applied to any non-Gaussian uncertainties, avoiding strong assumptions about the distribution form and significantly improving its applicability in complex systems such as IEGS.

[0080] 2. The orthogonality of the basis functions is maintained, which is beneficial for subsequent sparse recovery and numerical stability. The one-dimensional orthogonal polynomial constructed through Cholesky decomposition and third-order recurrence relations ensures that the generated basis functions satisfy the inner product orthogonality condition in their domains. This orthogonality guarantees that the measurement matrix Φ has a good condition number, which helps improve the numerical stability of the sparse aPC coefficient solution process and reduces redundant correlations between basis functions, thereby enhancing the accuracy and robustness of the surrogate model.

[0081] Preferably, in S5, the sample vector For: Y = [y 1 ,…,y i …, y M ] T Where M is the current number of points; y i This is the key output value corresponding to the i-th collocation point;

[0082] Measurement Matrix Let M×P be a matrix, and let the i-th row be... ;in, Let d be the d-dimensional random input vector for the i-th collocation point. Let j be the j-th multidimensional aPC basis function generated in S3;

[0083] And based on and Establish a linear relationship Where A = [a0, a1, …, a P-1 ] T Let be the vector of aPC coefficients to be determined.

[0084] This setup, by constructing a linear relationship between the output and the basis functions, enables the algebraic representation of complex IEGS models, providing crucial support for the efficient solution of sparse aPC coefficients and the rapid evaluation of surrogate models.

[0085] Preferably, in S7, the statistical moments of the key output quantities analytically calculated through the sparse aPC coefficient vector A include:

[0086] Using the orthogonality and normalization properties of aPC basis functions, calculate:

[0087] Expected value of key output ;

[0088] Variance of key outputs ;

[0089] Where a0 is the coefficient of the constant term; a1 to a p-1 These are the sparse nonzero coefficients corresponding to higher-order basis functions.

[0090] This setup enables closed-form analytical computation of statistical moments, significantly improving the efficiency of uncertainty quantification. Traditional Monte Carlo methods require extensive sampling and repeated solutions to the energy flow model to estimate the expectation and variance, resulting in extremely high computational costs. In contrast, this scheme, based on the orthogonality of aPC expansion, only requires the use of sparse coefficients a0, a1, …, a P-1The expected value and variance of the output can be obtained directly without additional simulation. This analytical approach transforms complex probabilistic analysis into simple algebraic operations, greatly reducing the computational burden. It is particularly suitable for rapid uncertainty assessment in high-dimensional, strongly coupled scenarios in IEGS systems.

[0091] Preferably, in S4, an improved Latin hypercube sampling method is used to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions; wherein, the improvement of the improved Latin hypercube sampling method includes: on the basis of traditional Latin hypercube sampling, introducing a hierarchical weighting strategy based on the probability density function to improve the sampling representativeness of non-uniformly distributed regions.

[0092] This setup improves sampling accuracy in non-uniformly distributed regions and enhances the generalization ability of the surrogate model. Traditional LHS performs stratified sampling within uniform intervals, making it difficult to effectively capture the characteristics of non-Gaussian input variables in high probability density regions. This scheme, however, introduces a weighting strategy based on the probability density function, allocating more sampling points to regions with higher probability density, thereby increasing the sampling density in these key areas. This enables the subsequently constructed aPC surrogate model to more accurately reflect the system's response behavior under high-probability conditions, significantly improving the accuracy and reliability of model predictions. Attached Figure Description

[0093] To make the objectives, technical solutions, and advantages of the invention clearer, the invention will now be described in further detail with reference to the accompanying drawings, wherein:

[0094] Figure 1 This is a flowchart of the method;

[0095] Figure 2 PDFs of the pressure at node 7 calculated using different methods in Example 2;

[0096] Figure 3 The aPC coefficient represents the pressure at node 7 in Example 2;

[0097] Figure 4 This is a box plot showing the relative pressure error at each node in Example 2. Detailed Implementation

[0098] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0099] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, not all of them. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0100] It should be noted that similar reference numerals and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. In the description of this invention, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the figures, or the orientation or positional relationship commonly used when the product is in use. They are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first," "second," and "third," etc., are only used to distinguish descriptions and should not be construed as indicating or implying relative importance. In addition, the terms "horizontal," "vertical," etc., do not indicate that the component is required to be absolutely horizontal or suspended, but can be slightly tilted. For example, "horizontal" simply means that its direction is more horizontal than "vertical," and does not mean that the structure must be completely horizontal, but can be slightly tilted. In the description of this invention, it should also be noted that, unless otherwise explicitly specified and limited, the terms "set," "install," "connect," and "link" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0101] Example 1

[0102] like Figure 1 As shown, this invention provides an IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model, comprising the following steps:

[0103] S1. Establish a deterministic energy flow model for the Integrated Energy System (IEGS) to solve for the system state variables under given input conditions.

[0104] In practice, the IEGS energy flow deterministic model includes the gas flow balance equation in the gas network, the power flow equation in the power grid, the gas flow equation in the gas pipeline, the compressor model, and the gas turbine model; the key outputs include the voltage amplitude of the power grid node, the phase angle, the pressure of the gas network node, and / or the pipeline flow rate.

[0105] The gas flow balance equation in the gas network is:

[0106] ;

[0107] In the formula, This indicates the incoming airflow from the gas source; This indicates the gas flow in natural gas pipeline ij. This indicates the airflow to the compressor between the natural gas pipelines IK; This represents the natural gas load consumed by the gas turbine unit connected to node i; This represents the natural gas load at node i; This indicates the total number of natural gas pipelines in the gas network; Indicates the number of compressors in IEGS;

[0108] The power flow equations are:

[0109] ;

[0110] In the formula, and These represent the active and reactive power generation at bus i, respectively; This indicates the active power consumed by the compressor in the natural gas network; δ represents the voltage magnitude at bus i; ij =δ i -δ j δ i This represents the voltage phase angle at bus i; This represents the voltage amplitude at bus j; This represents the active power load demand at bus i; This represents the total number of buses in the power grid; and Let i and j represent the real and imaginary parts of the nodal admittance matrix between bus i and j, respectively;

[0111] The gas flow equation in a gas pipeline is:

[0112] ;

[0113] In the formula, This represents the natural gas pressure at node i in the natural gas network; This represents the natural gas pressure at node j in the natural gas network; This represents the flow coefficient of pipe ij.

[0114] The gas flow balance equations for the gas network consider key aspects such as gas source injection, pipeline transportation, compressor pressurization, and gas consumption by gas turbine units, comprehensively reflecting the supply, demand, and transmission characteristics of natural gas. The power flow equations for the power grid, based on the standard power balance principle, accurately describe the relationship between node voltage and power flow. The gas pipeline flow equations introduce a nonlinear flow model driven by pressure difference, reflecting the physical nature of gas flow in pipelines. Together, these three constitute a fundamental IEGS model with strong physical consistency, effectively supporting subsequent uncertainty analysis and operational optimization tasks, significantly enhancing the model's applicability and reliability in practical engineering scenarios.

[0115] The compressor model is as follows:

[0116] ;

[0117] In the formula, This indicates the active power consumed by the compressor. and These represent the pressures at the compressor inlet and outlet nodes, respectively. K represents the compressor power conversion factor. G Indicates the adiabatic index of natural gas;

[0118] The gas turbine unit model is as follows:

[0119] ;

[0120] In the formula, This indicates the active power generated by the gas turbine unit; , , These are the quadratic, linear, and constant coefficients of the fuel consumption function, obtained by fitting unit operating data, which together characterize the nonlinear relationship between unit power and gas consumption.

[0121] This approach accurately reflects the nonlinear characteristics of the electro-gas energy conversion process, improving the accuracy of system modeling. As a key regulating device in natural gas pipeline networks, the compressor's energy consumption exhibits a strongly nonlinear relationship with its pressure ratio, making traditional linear approximations prone to error accumulation. Furthermore, the fuel consumption of gas turbine units varies with output, exhibiting a typical quadratic curve characteristic, with significant efficiency reduction, especially in low-load areas. This solution employs a compressor model and a secondary fuel consumption model based on thermodynamic characteristics, which more realistically reflects the operating behavior of these devices under different operating conditions. This effectively improves the accuracy of IEGS overall energy flow calculations, particularly in multi-energy synergistic optimization and uncertainty propagation analysis.

[0122] S2. Based on the statistical characteristics of historical operating data or prediction errors, identify the key uncertainty sources affecting energy flow in the integrated energy system, model the prediction deviation of each key uncertainty source as a random input variable, thereby forming a set containing d mutually independent random input variables, and configure a corresponding non-Gaussian probability density function for each random input variable according to its actual prediction error distribution; where d is the number of key uncertainty sources.

[0123] S3. Based on the number of random input variables d and the preset maximum expansion order p, using the non-Gaussian probability density function of each random input variable, construct an orthogonal univariate polynomial basis for each variable through a data-driven sample moment numerical orthogonalization method; then generate P multidimensional orthogonal aPC basis functions by performing tensor product operations on each orthogonal univariate polynomial basis.

[0124] In practice, the process of generating P multidimensional orthogonal aPC basis functions includes:

[0125] (a) Using the observed values ​​of random input variables corresponding to historical operating data or prediction error samples in S2 as driving sample data, calculate the first 2P+1 order statistical moments u of each random input variable. m m = 0, 1, …, 2P; and form a one-dimensional statistical moment matrix. :

[0126] ;

[0127] (b) By Cholesky decomposition This yields the upper triangular matrix. :

[0128] ;

[0129] (c) According to Calculate the recursive parameter b from the elements. k and c k :

[0130] ;

[0131] In the formula, r 0,0 = 1 and r 0,1 = 0; b k and c k This represents the coefficients determined by random input, k = 1, 2, …, P;

[0132] (d) Using parameter b k and c k One-dimensional orthogonal polynomial basis functions are constructed through three recurrence relations. , making it satisfy the orthogonality condition on the domain D:

[0133] ;

[0134] where, is the cumulative distribution function of the d-dimensional random input vector ;

[0135] (e) Perform the above steps (a)-(d) on d random input variables respectively to obtain d groups of univariate orthogonal bases, and then construct P multi-dimensional orthogonal aPC basis functions through tensor product:

[0136] .

[0137] Traditional aPC methods usually rely on known probability distributions (such as Beta, Gamma, etc.) to preset the orthogonal polynomial family, and it is difficult to be directly applied to input variables with unknown or complex distributions. However, this scheme only requires the statistical moment information of the input variables (which can be estimated from historical data), and can generate orthogonal polynomial bases that satisfy orthogonality through Cholesky decomposition and recurrence. This makes the method applicable to any non-Gaussian distributed uncertainty source, avoiding strong assumptions about the distribution form, and greatly improving the applicability in complex systems such as IEGS. Moreover, it also maintains the orthogonality of the basis functions, which is beneficial to subsequent sparse recovery and numerical stability. The one-dimensional orthogonal polynomials constructed through Cholesky decomposition and third-order recurrence relations ensure that the generated basis functions satisfy the inner product orthogonality condition on the domain. This orthogonality ensures that the measurement matrix Φ has a good condition number, helps to improve the numerical stability of the sparse aPC coefficient solving process, and reduces the redundant correlation between basis functions, thereby enhancing the accuracy and robustness of the surrogate model.

[0138] S4. Set the initial number of collocation points M0, where M0 < P; adopt a preset sampling method to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions.

[0139] Specifically, when implementing, an improved Latin hypercube sampling method is used to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions; among them, the improvement of the improved Latin hypercube sampling method includes: on the basis of traditional Latin hypercube sampling, a hierarchical weighting strategy based on the probability density function is introduced to improve the sampling representativeness of non-uniform distribution regions.

[0140] Traditional LHS (Local Hierarchical Sampling) performs stratified sampling within uniform intervals, making it difficult to effectively capture the characteristics of non-Gaussian input variables in high probability density regions. This approach, however, introduces a weighting strategy based on the probability density function, allocating more sampling points to regions with higher probability density, thereby increasing the sampling density in these critical areas. This enables the subsequently constructed aPC surrogate model to more accurately reflect the system's response behavior under high-probability conditions, significantly improving the accuracy and reliability of model predictions.

[0141] S5. Substitute the random input variables corresponding to each current collocation point into the energy flow deterministic model of S1 to obtain the corresponding system state variables. Organize all the obtained system state variables into a sample vector Y, which serves as the key output. For the i-th collocation point, calculate the basis function values ​​of all P aPC basis functions in S3 at that location, and use this P-dimensional row vector as the i-th row of the measurement matrix Φ. Collect the basis function value row vectors of all collocation points to form the complete measurement matrix Φ. The problem of solving the sparse aPC coefficient vector A is transformed into the following l1-l2 norm minimization optimization model:

[0142] ;

[0143] In the formula, λ is the regularization parameter, and λ>0;

[0144] Solve the model to obtain the sparse aPC coefficient vector A under the current collocation set; and based on the expression Construct a sparse aPC proxy model; where, Let d represent a random input vector consisting of d random input variables; Denotes the basis function of the i-th multivariable orthogonal polynomial; Let represent the expansion coefficients corresponding to the i-th aPC basis function.

[0145] In practical implementation, sample vector For: Y = [y 1 ,…,y i …, y M ] T Where M is the current number of points; y i This is the key output value corresponding to the i-th collocation point;

[0146] Measurement Matrix Let M×P be a matrix, and let the i-th row be... ;in, Let d be the d-dimensional random input vector for the i-th collocation point. Let j be the j-th multidimensional aPC basis function generated in S3;

[0147] And based on and Establish a linear relationship Where A = [a0, a1, …, a P-1 ] T Let be the vector of aPC coefficients to be determined.

[0148] In this way, by constructing a linear relationship between the output and the basis functions, the algebraic representation of the complex IEGS model is realized, which provides key support for the efficient solution of sparse aPC coefficients and the rapid evaluation of the surrogate model.

[0149] In practice, the alternating direction multiplier method is used to solve the l1-l2 norm minimization optimization model. The process includes:

[0150] Will Rewritten as G(A)-H(A), where G(A) and H(A) are convex:

[0151] ;

[0152] Linearize H(A) as follows, where A (n) ≠0:

[0153] ;

[0154] Will Simplified to:

[0155] ;

[0156] In the formula, z = Φ T Y +λA (n) / ||A (n) ||2;

[0157] Introduce an auxiliary variable B to decouple the non-smooth l1 norm sum. The smooth portion; the constraint A = B is enforced by an augmented Lagrange formula, where B handles the l1 penalty term, while A manages the remaining terms; the augmented Lagrange formula is as follows:

[0158] ;

[0159] right Taking the partial derivatives with respect to A, B, and u, we get:

[0160] ;

[0161] After giving the initial values, set the three partial derivatives to 0, and then iteratively solve A, B and u until the error is less than the preset error threshold; thus obtaining the sparse aPC coefficient vector A under the current collocation set.

[0162] The original optimization model contains non-smooth l1 and l2 norm terms, making direct solution difficult. This scheme decomposes the objective function into two convex functions, G(A) and H(A), linearizes H(A), and introduces an auxiliary variable B to separate the l1 term, constructing an augmented Lagrangian function. This structured approach avoids the numerical instability caused by directly handling non-smooth terms, making the optimization process smoother and more controllable, significantly improving the algorithm's robustness and convergence performance. Furthermore, using the ADMM framework, the complex joint optimization problem is decomposed into three subproblems concerning A, B, and the dual variable u, each of which can be solved analytically or iteratively. For example, updating A involves a regularized quadratic optimization problem, which can be quickly solved analytically or using the conjugate gradient method; updating B only requires a soft threshold operation. This block-based solution strategy avoids large-scale matrix inversion or global search, significantly reducing the computational burden of each iteration, making it suitable for real-time applications in high-dimensional IEGS systems.

[0163] In the iterative solution of A, B, and u, l1 minimizes the following solution as the initial value A. (1) :

[0164] .

[0165] Since l1 norm minimization is a classic method for sparse recovery, it can effectively approximate the optimal sparse solution under certain conditions (such as restricting the isometric property). This solution is used as the initial value A for the ADMM iteration. (1) This approach allows the optimization process to start from a point close to a real sparse structure, avoiding getting trapped in local poor solutions or prolonged oscillations. Compared to random initialization, this strategy significantly improves iteration efficiency, reduces the number of steps required to reach convergence, and enhances the overall robustness of the algorithm.

[0166] To facilitate better understanding, the following explanation is provided.

[0167] The dimension of aPC increases dramatically due to tensor products with many input variables. Choosing M CPs and computing the QoIs of the original system is impractical. A dimensionality reduction method for constructing an aPC model using fewer CPs is described below. Consider the case where the number of CPs is less than the dimension of the aPC, i.e., M < N. Then the equation becomes an ill-conditioned problem with no unique solution. To alleviate this difficulty, a unique solution can be obtained by adding appropriate constraints. In this invention, compressed sensing technology is applied to recover the sparse signal to solve the underdetermined system of equations. This method can be considered a dimensionality reduction method for the aPC model. "Sparseness" refers to the number of non-zero elements in A, defined as follows:

[0168] ;

[0169] To find a sparse solution to A in the above equation, CS is described as minimizing l0:

[0170] ;

[0171] Since the above is an NP-hard problem, we use l1-l2 minimization to obtain an approximate solution.

[0172] ;

[0173] Its unconstrained equivalent modulus is:

[0174] This refers to the l1-l2 norm minimization optimization model transformed in S5 of this method.

[0175] The subsequent details, such as the specific process of S5, will not be repeated here.

[0176] Table 1 shows the pseudocode for minimizing the aPC coefficients using l1-l2.

[0177] Table 1

[0178]

[0179] In the algorithm in Table 1, shrink is defined as:

[0180] .

[0181] S6. Based on the sparse aPC surrogate model of S5, calculate the relative squared error on the additional set of verification collocations generated in the joint probability space defined by d random input variables and their non-Gaussian probability density functions; if the relative squared error is greater than the preset threshold γ, increase the number of collocations by the preset increment and return to step S5; if the relative squared error is less than or equal to γ, terminate the iteration and output the final sparse aPC surrogate model.

[0182] In practice, the sampling method used in S6 is the same as that used in S4.

[0183] S7. Utilizing the orthogonality of the aPC basis functions in the final sparse aPC surrogate model of S6, the statistical moments of the key outputs are analytically calculated through the sparse aPC coefficient vector A, and the probability density function of the key outputs is reconstructed based on the statistical moments.

[0184] In specific implementation, the statistical moments of the key output quantities analytically calculated through the sparse aPC coefficient vector A include:

[0185] Using the orthogonality and normalization properties of aPC basis functions, calculate:

[0186] Expected value of key output ;

[0187] Variance of key outputs ;

[0188] Where a0 is the coefficient of the constant term; a1 to a p-1 These are the sparse nonzero coefficients corresponding to higher-order basis functions.

[0189] Traditional Monte Carlo methods require extensive sampling and repeated solutions to the energy flow model to estimate the expectation and variance, resulting in extremely high computational costs. In contrast, this scheme, based on the orthogonality of aPC expansion, only requires the use of sparse coefficients a0, a1, …, a P-1 The expected value and variance of the output can be obtained directly without additional simulation. This analytical approach transforms complex probabilistic analysis into simple algebraic operations, greatly reducing the computational burden. It is particularly suitable for rapid uncertainty assessment in high-dimensional, strongly coupled scenarios in IEGS systems.

[0190] S8, based on the probability density function of the reconstructed key output quantity of S7, evaluates the safety margin of IEGS operation, and generates a risk warning signal when the safety margin is lower than the preset safety threshold, which is used to trigger real-time control or day-ahead scheduling adjustment.

[0191] This method employs a non-Gaussian probability density function based on historical operating data or the statistical characteristics of prediction errors to configure a corresponding probability distribution for each key uncertainty source, constructing an arbitrary chaotic polynomial (aPC) basis function. Unlike traditional PCE methods that rely on Gaussianization transformations, this technique directly utilizes the probability characteristics driven by actual data, avoiding modeling errors caused by biased distribution assumptions. Especially in IEGS, electricity, gas loads, and renewable energy output often exhibit non-Gaussian properties such as skewness and heavy tails. This method can more realistically reflect the statistical behavior of input uncertainties, thereby significantly improving the accuracy of the surrogate model. In addition, addressing the problem of a large number of random variables and exponential growth of polynomial terms in IEGS systems, this method introduces the idea of ​​sparse recovery, transforming the aPC coefficient solution problem into an l1-l2 norm minimization optimization model. Compared to traditional full-order expansion methods, this method retains only a few polynomial terms that significantly contribute to the system response, greatly reducing the effective degrees of freedom. Simultaneously, by combining the construction of the measurement matrix and an iterative point-addition mechanism, the computational scale is further controlled. Compared to conventional aPC or high-dimensional sparse regression methods, this strategy significantly compresses the model dimensionality while maintaining model accuracy, thus improving scalability.

[0192] This method establishes a closed-loop iterative optimization process by setting initial collocation points and dynamically adjusting the sampling quantity based on the relative squared error on the validation set. If the model accuracy is insufficient, collocation points are automatically added to improve the estimation quality; otherwise, the iteration terminates. This adaptive strategy avoids the resource waste caused by manually setting too many samples and prevents underfitting due to insufficient samples. Compared to fixed sampling strategies or one-time large-sample methods, this mechanism minimizes unnecessary computational overhead while ensuring model accuracy, enhancing the algorithm's intelligence and practicality. Furthermore, this method embeds a sparse aPC surrogate model into the IEGS energy flow analysis process, calculating the statistical moments of key outputs analytically and reconstructing their probability density functions to evaluate the system's operational safety margin. Because the surrogate model has rapid evaluation capabilities, it can complete uncertainty propagation analysis in a large number of scenarios in a very short time, making it suitable for real-time monitoring and day-ahead scheduling scenarios. Compared to Monte Carlo simulations, which require thousands of simulations to obtain reliable results, this method significantly shortens the computation time, enabling probabilistic energy flow analysis to have online application potential.

[0193] This method can effectively reduce model dimensionality and improve computational efficiency by combining sparse recovery techniques while retaining the good adaptability of the aPC method to non-Gaussian inputs, thereby achieving high-precision, low-cost, and real-time analysis of IEGS probabilistic energy flow.

[0194] Example 2

[0195] To better illustrate the effectiveness of this method, the following economic benefit calculations and simulation experiments are conducted to demonstrate its effectiveness.

[0196] Test System I consists of an IEEE 39-bus power system and a Belgian 20-bus gas system. Forecasts for gas and electricity loads are matched to their original loads in System I (both active loads), while other loads remain constant. Renewable energy and electricity load forecast error data are derived from open-access day-ahead forecasts for Belgian loads and onshore wind power from October 31 to November 30, 2024, with a resolution of 15 minutes. Gas load forecast data also ensures a consistent 15-minute resolution.

[0197] It is assumed that all data are independently and identically distributed across all buses / nodes and generators. This assumption is supported by the Pearson correlation coefficients among the prediction error samples: RES is correlated with electrical load at -0.1185, with natural gas load at -0.0574, and with natural gas load at -0.0336, indicating that linear dependence is negligible. In this invention, since load and RES are considered fixed inputs, their correlation is not considered. If the prediction error exhibits significant dependence, methods such as Copula-based modeling, whitening transformation, or Nataf transformation can be applied.

[0198] All case studies were conducted on a laptop equipped with an Intel™ Core Ultra 5 125H CPU 1.20 GHz and 32GB of RAM.

[0199] The computational efficiency and dimensionality reduction performance of the aPC-CS method on System I were evaluated. The baseline method was a Monte Carlo simulation using 10,000 samples. Comparisons were also made with the following models: an aPC model based on the stochastic response surface methodology and an analysis of variance model with ϵ=0.99. CPs were sampled from a probability box model obtained through a two-level sampling method. All PDFs were fitted using kernel density estimation. The maximum order of aPC was set to 2, and M0 and ΔM were set to 100 and 25, respectively. λ, σ, and γ were set to 1×10⁻⁶. -8 1×10 -2 and 1×10 -3 .

[0200] System I's power grid has 39 buses, 10 generators, and 46 branches, while the natural gas network has 20 nodes, 4 gas wells, 2 gas storage facilities, and 19 pipelines. There are 3 compressors and 3 gas turbine units between the subsystems. The RES generators and their predicted values ​​are as follows: 120MW at bus 6, 80MW at bus 14, and 100MW at bus 19. The pressure (p7) at node 7 is selected as QoI. Figure 2 The PDF shows the pressure at node 7 calculated using different aPC class methods. From Figure 2 As can be seen, the PDFs of all aPC-class methods are almost identical to those of MCS. This indicates that aPC-class methods can provide the same high accuracy as MCS.

[0201] Table 2. Relative errors of pressure at node 7 calculated by different methods.

[0202]

[0203] Table 2 shows the relative errors for dimension, computation time, mean, variance, skewness, and kurtosis. MCS was used as the baseline. To quantify the effectiveness of precision and sparsity, the relative errors and variance retention rates were calculated:

[0204] ;

[0205] ;

[0206] A relative error close to 0 indicates high consistency with MCS, while a variance retention rate close to 1 suggests that the sparse model retains sufficient information, thus supporting the sparsity assumption. Table 2 shows that aPC-CS significantly reduces dimensionality while maintaining a small relative error. Lower dimensionality means fewer computational steps are needed, greatly reducing computation time. Although aPC-CS requires more computation time than aPC-ANOVA due to the iterative optimization algorithm used, the trade-off is that aPC-CS has smaller relative errors in mean, variance, skewness, and kurtosis. Furthermore, the variance retention rates of both aPC-ANOVA (0.9839) and aPC-CS (0.9905) are close to 1, indicating that both sparse models retain sufficient information.

[0207] Figure 3 The aPC coefficients for the pressure at node 7 are shown for three aPC methods. The left column represents the first-order aPC coefficients, and the right column represents the second-order aPC coefficients. Specifically, (a), (c), and (e): first-order coefficients for aPC, aPC-ANOVA, and aPC-CS; (b), (d), and (f): second-order coefficients for aPC, aPC-ANOVA, and aPC-CS. Figure 3 It can be seen that the first-order coefficients are much larger than the second-order coefficients. This indicates that higher-order basis functions have a weaker influence on the approximate QoI. The main difference in coefficients among these methods lies in the second-order coefficients. Many coefficients in aPC-ANOVA and aPC-CS are small or even zero, reflecting the inherent sparsity of the expansion. Since the proposed method uses fewer CPs to compute the coefficients, the aPC model has lower dimensionality and sparsity. According to the definition of sparsity, Figure 3 This reflects that the proposed method has more zero-valued second-order coefficients.

[0208] Figure 4 The relative error distribution of gas pressure at the gas load nodes in System I was compared. It can be seen that the overall error is less than 0.3. MCS verification shows that the error distribution is highly concentrated near the mean, validating the accuracy of this method. However, some outliers still appeared at nodes 6, 19, and 20, indicating significant deviations in the distribution of these nodes. These outliers contribute to a relatively large skewness in the RE% index. Figure 4This explains why the skewness error is significantly larger than the error in lower-order moments. First, the true skewness in a nearly symmetric input distribution is very small, which amplifies the relative error. Second, while the sparse aPC model retains the mean and variance information during dimensionality reduction, it loses most of the skewness information, leading to skewness bias.

[0209] Table 2 also evaluates the dimensionality reduction performance of aPC-CS and aPC-ANOVA on system I. As shown in Table 2, aPC-ANOVA achieved a dimensionality reduction of approximately 66% for system I, from 595 to 199. On the other hand, aPC-CS reduced the dimensionality of system I from 595 to 107, a dimensionality reduction of 80%. For f... 33 aPC-CS reduces the dimensionality of the aPC model from 7,750 to only 3, while aPC-ANOVA reduces it to 127. These results demonstrate that aPC-CS is more efficient than aPC-ANOVA in dimensionality reduction.

[0210] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit the technical solutions. Those skilled in the art should understand that any modifications or equivalent substitutions to the technical solutions of the present invention without departing from the spirit and scope of the present invention should be covered within the scope of the claims of the present invention.

Claims

1. A probabilistic energy flow monitoring method for IEGS based on a sparse arbitrary chaotic polynomial model, characterized in that, It includes the following steps: S1. Establish a deterministic model of the energy flow of the integrated energy system (IEGS) to solve the system state variables under given input conditions. S2. Based on the statistical characteristics of historical operation data or prediction errors, identify the key uncertainty sources affecting the energy flow in the integrated energy system. Model the prediction deviation of each key uncertainty source as a random input variable, thus forming a set containing d mutually independent random input variables, and configure the corresponding non-Gaussian probability density function for each random input variable according to its actual prediction error distribution. Here, d is the number of key uncertainty sources. S3. According to the number d of random input variables and the preset maximum expansion order p, use the non-Gaussian probability density functions of each random input variable to construct orthogonal univariate polynomial bases for each variable respectively through a numerical orthogonalization method based on data-driven sample moments. Then, generate P multi-dimensional orthogonal aPC basis functions through tensor product operations on the orthogonal univariate polynomial bases. S4. Set an initial number of collocation points M0, where M0 < P. Use a preset sampling method to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions. S5. Substitute the random input variables corresponding to the current collocation points into the energy flow deterministic model in S1 for solution to obtain the corresponding system state variables. Form a sample vector Y by arranging all the obtained system state variables in columns as the key output quantity. For the i-th collocation point, calculate the basis function values of all P aPC basis functions in S3 at this point, and use this P-dimensional row vector as the i-th row of the measurement matrix Φ. Pool the basis function value row vectors of all collocation points to form a complete measurement matrix Φ. Transform the problem of solving the sparse aPC coefficient vector A into the following l1-l2 norm minimization optimization model: ； In the formula, λ is the regularization parameter, and λ > 0. Solve the model to obtain the sparse aPC coefficient vector A under the current collocation set; and based on the expression Construct a sparse aPC proxy model; where, Let d represent a random input vector consisting of d random input variables; Denotes the basis function of the i-th multivariable orthogonal polynomial; Represents the expansion coefficients corresponding to the i-th aPC basis function; S6. Based on the sparse aPC surrogate model in S5, calculate the relative squared error on an additional generated verification collocation point set in the joint probability space defined by d random input variables and their non-Gaussian probability density functions. If the relative squared error is greater than the preset threshold γ, increase the number of collocation points by a preset increment and return to step S5. If the relative squared error is less than or equal to γ, terminate the iteration and output the final sparse aPC surrogate model. S7. Utilize the orthogonality of the aPC basis functions in the final sparse aPC surrogate model in S6 to analytically calculate the statistical moments of the key output quantity through the sparse aPC coefficient vector A, and reconstruct the probability density function of the key output quantity based on the statistical moments. S8. Based on the reconstructed probability density function of the key output quantity in S7, evaluate the operation safety margin of the IEGS, and generate a risk warning signal when the safety margin is lower than the preset safety threshold to trigger real-time control or day-ahead scheduling adjustment.

2. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 1, characterized in that: In S5, use the alternating direction multiplier method to solve the l1-l2 norm minimization optimization model, and the process includes: Will Rewritten as G(A)-H(A), where G(A) and H(A) are convex: ； Linearize H(A) as follows, where A (n) ≠0: ； Will Simplified to: ； In the formula, z = Φ T Y +λA (n) / ||A (n) ||2; Introduce an auxiliary variable B to decouple the non-smooth l1 norm sum. The smooth portion; the constraint A = B is enforced by an augmented Lagrange formula, where B handles the l1 penalty term, while A manages the remaining terms; the augmented Lagrange formula is as follows: ； right Taking the partial derivatives with respect to A, B, and u, we get: ； After giving the initial values, set the three partial derivatives to 0, and then iteratively solve A, B and u until the error is less than the preset error threshold; thus obtaining the sparse aPC coefficient vector A under the current collocation set.

3. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 2, characterized in that: In the iterative solution of A, B, and u, l1 minimizes the following solution as the initial value A. (1) : 。 4. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 1, characterized in that: In S1, the deterministic energy flow model of IEGS includes the gas flow balance equation in the gas network, the power flow equation in the power grid, the gas flow equation in the gas pipeline, the compressor model, and the gas turbine model; the key outputs include the voltage amplitude of the power grid node, the phase angle, the pressure of the gas network node, and / or the pipeline flow.

5. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 4, characterized in that: The gas flow balance equation in the gas network is: ； In the formula, This indicates the incoming airflow from the gas source; This indicates the gas flow in natural gas pipeline ij. This indicates the airflow to the compressor between the natural gas pipelines IK; This represents the natural gas load consumed by the gas turbine unit connected to node i; This represents the natural gas load at node i; This indicates the total number of natural gas pipelines in the gas network; Indicates the number of compressors in IEGS; The power flow equations are: ； In the formula, and These represent the active and reactive power generation at bus i, respectively; This indicates the active power consumed by the compressor in the natural gas network; δ represents the voltage magnitude at bus i; ij =δ i -δ j δ i This represents the voltage phase angle at bus i; This represents the voltage amplitude at bus j; This represents the active power load demand at bus i; This represents the total number of buses in the power grid; and Let i and j represent the real and imaginary parts of the nodal admittance matrix between bus i and j, respectively; The gas flow equation in a gas pipeline is: ； In the formula, This represents the natural gas pressure at node i in the natural gas network; This represents the natural gas pressure at node j in the natural gas network; This represents the flow coefficient of pipe ij.

6. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 5, characterized in that: The compressor model is as follows: ； In the formula, This indicates the active power consumed by the compressor. and These represent the pressures at the compressor inlet and outlet nodes, respectively. K represents the compressor power conversion factor. G Indicates the adiabatic index of natural gas; The gas turbine unit model is as follows: ； In the formula, This indicates the active power generated by the gas turbine unit; , , These are the quadratic, linear, and constant coefficients of the fuel consumption function, obtained by fitting unit operating data, which together characterize the nonlinear relationship between unit power and gas consumption.

7. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 1, characterized in that: In S3, the process of generating P multidimensional orthogonal aPC basis functions includes: (a) Using the observed values ​​of random input variables corresponding to historical operating data or prediction error samples in S2 as driving sample data, calculate the first 2P+1 order statistical moments u of each random input variable. m m = 0, 1, …, 2P; and form a one-dimensional statistical moment matrix. : ； (b) By Cholesky decomposition This yields the upper triangular matrix. : ； (c) According to Calculate the recursive parameter b from the elements. k and c k : ； In the formula, r 0,0 = 1 and r 0,1 = 0; b k and c k This represents the coefficients determined by random input, k = 1, 2, …, P; (d) Using parameter b k and c k One-dimensional orthogonal polynomial basis functions are constructed through three recurrence relations. This ensures that the orthogonality condition is satisfied on the domain D: ； In the formula, For d-dimensional random input vectors The cumulative distribution function; (e) Perform steps (a)-(d) above on d random input variables respectively to obtain d sets of univariate orthogonal bases, and then construct P multidimensional orthogonal aPC basis functions through tensor product.

8. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 7, characterized in that: In S5, the sample vector For: Y = [y 1 ,…,y i …, y M ] T Where M is the current number of points; y i This is the key output value corresponding to the i-th collocation point; Measurement Matrix Let M×P be a matrix, and let the i-th row be... ;in, Let d be the d-dimensional random input vector for the i-th collocation point. Let j be the j-th multidimensional aPC basis function generated in S3; And based on and Establish a linear relationship Where A = [a0, a1, …, a P-1 ] T Let be the vector of aPC coefficients to be determined.

9. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 8, characterized in that: In S7, the statistical moments for the key output quantities are analytically calculated using the sparse aPC coefficient vector A, including: Using the orthogonality and normalization properties of aPC basis functions, calculate: Expected value of key output ; Variance of key outputs ; Where a0 is the coefficient of the constant term; a1 to a p-1 These are the sparse nonzero coefficients corresponding to higher-order basis functions.

10. The IEGS probabilistic energy flow monitoring method based on a sparse arbitrary chaotic polynomial model as described in claim 1, characterized in that: In S4, an improved Latin hypercube sampling method is used to generate M0 initial collocation points in the joint probability space defined by d random input variables and their non-Gaussian probability density functions. The improvement of the improved Latin hypercube sampling method includes: on the basis of traditional Latin hypercube sampling, introducing a hierarchical weighting strategy based on the probability density function to improve the sampling representativeness of non-uniformly distributed regions.

Images