Systematic method for evaluating bearing capacity of regional multi-energy system based on time coupling consideration

Abstract

The invention belongs to the field of regional multi-energy system bearing capacity evaluation, and discloses a time coupling consideration-based systematic method for evaluating the bearing capacity of a regional multi-energy system, which adopts an EH method to model the system and introduces an EH steady-state safe region concept as a means for evaluating the bearing capacity of the system. In order to calculate a safe area, a multi-parameter planning algorithm is provided, and a mathematical model is established by using a standardized matrix modeling method. For a high-dimensional polyhedron generated by time coupling, a dimension reduction strategy that a low-dimensional polyhedron approximates an original polyhedron is adopted in the method. These low-dimensional polygons are only associated with a limited number of time periods, thereby reducing computational complexity. According to the method, better balance can be achieved between calculation precision and efficiency in the application process, and meanwhile the inherent time coupling characteristics of different time period combinations are reserved.

Description

Technical Field

[0001] The present invention belongs to the field of evaluation of the carrying capacity of regional multi - energy systems, and more specifically, relates to a systematic method for evaluating the carrying capacity of regional multi - energy systems considering time coupling. Background Art

[0002] In recent years, multi - energy systems (MES) have become increasingly important because they can integrate various energy sectors such as electricity, fuel, heat, and cooling, and can reduce costs and emissions compared to individual energy systems. By enabling interactions between different energies, MES can unlock the flexibility of conversion across multiple energy sectors, which has become a promising solution to the challenges of limited energy and climate change. Therefore, there is an increasing interest in studying MES and exploring its potential to provide a more sustainable and reliable energy system for the future.

[0003] Evaluating the carrying capacity of MES plays a crucial role in various aspects such as determining the feasibility of a planning scheme, identifying energy supply bottlenecks, and introducing optimization strategies. For various reasons, it is essential, including system optimization, reliability assessment, and future capacity expansion planning. By evaluating the carrying capacity, valuable insights are obtained into the system's ability to accommodate additional loads or power sources while maintaining optimal performance and stability. Therefore, evaluating the power supply capacity of MES has become an essential part of distribution network planning and MES planning and construction.

[0004] Although the definitions of power system security regions and MES are similar, their carrying capacity evaluations differ in several aspects. The calculation of power system security regions is mainly affected by system topology, so a power - flow - based method is usually used to evaluate the carrying capacity of power systems. In contrast, MES involves multiple forms of energy conversion, so the typical security region of MES is expected to exhibit several basic properties: 1) Since multi - energy coupling is not considered in MES, traditional single - network flow analysis is insufficient. Therefore, the conversion and transmission of various energy forms must be considered; 2) Modeling the operating characteristics of different types of energy converters and storage devices is crucial. Energy storage devices facilitate energy transfer in the spatial and temporal dimensions, and their operating characteristics significantly affect the performance of MES; 3) The time - continuity of the operation of energy storage devices should be considered because all forms of energy storage devices exhibit time - coupled operating characteristics.

[0005] Existing methods do not consider the time - variability of energy demand and supply in MES, which can significantly affect the carrying capacity of the system. Therefore, a systematic and effective method is needed to evaluate the carrying capacity of MES, taking into account the dynamics and complexity of the system. Summary of the Invention

[0006] In view of the deficiencies of the prior art and the improvement requirements, the present invention proposes a systematic method for evaluating the carrying capacity of a multi - energy system in an assessment area considering time coupling to improve its planning. First, a mathematical form of the EH is established using a standardized matrix modeling method. Then, a multi - parameter programming (MPP) algorithm is proposed to calculate the safe region as a spatial projection problem. To handle the high - dimensional polyhedron generated by time coupling, some low - dimensional polyhedra are used to approximate the original polyhedron, and each polyhedron is only related to a few periods to reduce the computational burden. Compared with the existing methods, the proposed method provides an improved balance between computational accuracy and efficiency while retaining the inherent time - coupling characteristics between different period combinations.

[0007] To achieve the above object, the present invention discloses a multi - parameter programming algorithm (abbreviation: MPP) to effectively solve the safe region calculation problem.

[0008] The steps to solve the MPP problem are as follows: 1) Optimal partitioning: The first step in solving the MPP problem is to partition the feasible region into subsets, each corresponding to a unique solution. This is achieved by identifying the extreme points of the feasible region. The number of extreme points is equal to the dimension of the problem. For example, a two - dimensional problem has four extreme points, determined by the maximum and minimum values along each axis. All the above - mentioned constraint sets can be divided into active constraints and non - active constraints, which are respectively expressed as: A A x * (θ)+B A θ=C A A I x * (θ)+B I θ<C I 2) Critical region: The critical region represents the set of parameter values θ that share a common "optimal partition", which is the region where the active constraints change. These boundaries are usually determined by solving a system of linear equations that represent the active constraints when transitioning from one polyhedron to another. Each θ corresponds to a specific combination of active and non - active constraints.

[0009] For any θ i ∈Θ, the optimal solution is obtained by affine transformation: Substituting the above formula into the non - active constraints, we get That is: Therefore, the critical region associated with θ i is defined as: It can be expressed as: Where: The present invention also proposes an innovative algorithm for calculating the steady-state safety region in EH, approximating a high-dimensional polytope as a low-dimensional polytope, reducing the computational burden while retaining the original time coupling of various cycle combinations: The space of the output variables is high-dimensional, which means that exploring the high-dimensional safety region with time coupling in the exploration space will cause the computational burden to increase exponentially with the increase in the number of time periods. There are methods to exhaustively list all possible combinations of time periods to explain the original time coupling. However, this will lead to a large number of combination enumerations, expressed as Where represents the number of combinations, j represents the number of selected time periods, n out represents the number of system output streams, and T is the total number of time periods. This indicates that the calculation process involves a large number of calculations of low-dimensional polytopes. Therefore, the present invention discloses a fast approximation method to evaluate the carrying capacity of MES. The fast approximation method is based on the concept of decomposing an exact high-dimensional polytope into the intersection of several low-dimensional polytopes. These low-dimensional polytopes approximate the safety region by reflecting the time coupling between them. The coupling between specific time periods or power, heat, and cooling outputs. Therefore, the method first needs to select the i-th low-dimensional polytope θ i The number of time periods involved. Thus, the high-dimensional polytope θ0 can be decomposed into low-dimensional polytopes θ i . Suppose the method needs to combine k arbitrary time cycles, where k is a natural number not greater than nT. In this case, the total number of required low-dimensional polytopes is given by: Compared with the original method, the overall computational burden is significantly reduced.

[0010] Meanwhile, based on the above algorithm, the present invention also proposes a precision measurement method and precision index based on the difference between the high-dimensional safety region and its approximation to evaluate the precision loss brought by the proposed approximation: Regarding the constraint Θ = {V out ∈R K |V ∈ (1)-(12)} is the projection of the safety region, which is essentially the constraint G = {(V out , V) ∈ R K ×R B |V ∈ (1)-(12)} represented by the constraint space R×R projected onto the output space R, that is, the former is an approximation of the latter. Therefore, using the latter as a benchmark can effectively measure the precision loss of the proposed approximation without the need to calculate the exact safety region, thus avoiding a significant computational burden.

[0011] It should be noted that the most accurate accuracy evaluation method is to compare the deviations of all coordinate points on the boundary plane of the feasible region. However, a limitation of this method is that when approximating the original polytope with low-dimensional polytopes, the resulting high-dimensional polytope composed of multiple low-dimensional polytopes does not necessarily have the same number of vertices as the original polytope. In this case, it becomes challenging to find variables that can accurately evaluate the accuracy. To overcome the above challenges, a hyperplane D T V out can be constructed, where D = [d1, d2,..., d K T ∈ R K , and two optimization problems are designed related to the above two constraints: The vector D represents a specific combination of outputs in the given problem. The objective function g in both problems shows the maximum output that can be achieved in the case of a specific combination between the outputs described by the given vector D. We can quantify the accuracy of the proposed approximate safety region by comparing the vertices obtained from the constraint G = {(V out , V) ∈ R K × R B | V ∈ (1)-(12)} and those obtained from the constraint Θ = {V out ∈ R K | V ∈ (1)-(12)}. This comparison is crucial for ensuring the accuracy of the approximation method.

[0012] The Monte Carlo method is used to randomly generate the elements in D in order to randomly compare the vertices within the polytopes related to the above two constraints. This method ensures that the elements generated in D are random and within the range of [-1, 1], where the norm of the vector in this method is equal to 1. Brief Description of the Drawings

[0013] Figure 1 is a structural illustration diagram of nodes, ports, and branches in the Energy Hub (EH).

[0014] Figure 2 is the equivalent model of the MPP algorithm proposed by the present invention.

[0015] Figure 3 is the implementation of the proposed fast approximation method in a three-output system.

[0016] Figure 4 is the flowchart of the parameter solution algorithm of the proposed method.

[0017] Figure 5 ​and Figure 6 are respectively the structural diagrams of the energy hubs of the example systems 1 and 2.

[0018] Figure 7 and Figure 8 are the comparison diagrams of the steady-state safety regions of System 1 with and without considering time coupling. Figure 7 represents the safety region of System 1 considering time coupling. Figure 8 represents the safety region of System 1 without considering time coupling.

[0019] Figure 9 and Figure 10 are respectively the actual mapping and approximate mapping diagrams of the steady-state power safety regions during the 1st, 2nd, and 3rd time periods. Detailed implementation manner

[0020] In order to make the objectives, technical solutions, and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0021] The present invention provides a systematic method for evaluating the load-carrying capacity of a multi-energy system in an evaluation area considering time coupling, which is used to evaluate the load-carrying capacity of the MES. The method includes formulating an accuracy index to evaluate the calculation accuracy of the proposed method.

[0022] As shown in the appendix Figure 1 EH has multiple inputs and outputs, including four stages: energy input, conversion, storage, and output. Energy, power grids, and gas networks can be used as inputs to EH, while distribution networks, heat networks, and loads are usually outputs of EH. In energy network analysis, the matrix representation of a directed graph can be used to conveniently establish a mathematical model to describe the relationship between energy flows in the network. Elements such as components and energy transmission pipelines in EH directly correspond to the nodes and branches of network graph theory. Therefore, MES can be mathematically expressed using a directed graph, and matrices that can be used to describe the relationship between energy flows and components within the EH system are introduced.

[0023] The input (output) port-branch incidence matrix is: The energy conversion characteristics of node g are captured by the converter characteristic matrix H. Assume that node g has K ports and its characteristics can be described by P equations, then the corresponding dimension of matrix H is P×K: where η represents the efficiency of the energy conversion process s.

[0024] The energy conversion matrix can be obtained by calculating based on the first two types of operation matrices. The port-branch incidence matrix and the converter characteristic matrix represent the input-output relationship of energy devices and the energy conversion relationship in the EH respectively. Since the sub-matrix of the converter characteristic matrix is a diagonal matrix and its dimension is equal to that of the port-branch incidence matrix, their product can represent a set of energy conversion relationships of specific energy devices in the EH during the energy input-output process. Accordingly, the energy conversion matrix of node g can be deduced as follows: By aggregating the energy conversion matrices of all nodes, the energy conversion matrix of the EH can be obtained: Z = [Z1; Z2;....; Z g Based on the predefined matrix, the energy conversion equation can be given as: ZV = 0 where Z is defined in the energy conversion matrix and V is the energy flow vector containing the energy flows of all branches.

[0025] The energy input flow is represented by V, while V represents the energy output flow. The input (output) incidence matrices are as follows: Then the matrix form of the entire energy flow equation can be expressed as: where C in and C out represent the input and output incidence matrices respectively, V is the energy flow vector containing the energy flows of all branches, V in and V out represent the energy input and output flows respectively.

[0026] The energy flow relationship in the EH, as shown in the above matrix form of the energy flow equation, only captures the real-time energy balance and does not consider the energy storage components. Energy storage devices such as batteries, compressed air energy storage, chilled water storage, thermal energy storage, and gas storage are not reflected.

[0027] To comprehensively establish the mathematical model of the network graph theory of the EH, energy storage nodes need to be included. The energy flow entering the energy storage node is regarded as the charging variable, while the energy flow leaving the energy storage node is regarded as the discharging variable. The charge-discharge equation of the energy storage device can be expressed as: where is g th ​The number of (output) ports of the storage component, S(t) is the state of charge (SoC) of the storage device within the time period t, and ΔS is the charge and discharge energy of the g storage component. The node feature matrix of the storage device can be defined as: The port-branch incidence matrix of the storage device is enhanced to: A mathematical expression for the SoC of the storage device. The energy conversion matrix of the storage device can be expressed as: Therefore, for the storage node, the energy balance equation can be derived from the energy charge and discharge equation as: where represents all the elements in V representing the energy flow in V, -I represents the identity matrix of dimension units, and ΔS represents all the elements in ΔS g in.

[0028] The energy flow should obey the branch energy constraint. In addition, the capacity of each converter should comply with the limitations of the input or output branches: where and are the input and output capacities of all nodes respectively, V Bmin and V B，max represent the upper and lower branch energy constraints respectively.

[0029] In addition to the SoC limitations of the energy storage device discussed above, the MES also faces ramp rate limitations imposed on the energy storage device and other devices. These constraints result in strong time coupling within the system. The change in the node output between two consecutive time intervals should comply with the ramp rate constraint, i.e., where and represent the ramp down and ramp up rates within the time period t. Generally speaking, when the storage capacity is large enough or the SoC is moderate, the SoC constraint of the energy storage device can be relaxed. However, the ramp rate limitation is usually non-negotiable and cannot be relaxed.

[0030] The steady-state security region of the power system is related to the power flow equation and power constraints. Similarly, for the EH, the security region is related to the device operation and safety constraints. The output port of the EH is the key to evaluating the load capacity of the MES. The security region represents the set of energy output vectors that satisfy the relevant constraint conditions. The present invention defines the hyperspace R K ×RB , which is constructed by v out (K×1) and V(B×1), where B is the number of branches in EH considering time coupling. All the above constraints form a polyhedron G together in the R K ×R B hyperspace: where (1)-(12) represent the first 12 constraints mentioned above. In addition, in order to obtain the safety region θ’ the polyhedron G in R×R should be projected onto R°. For this, the present invention proposes a multi-parameter programming algorithm (MPP) to effectively solve the safety region calculation problem.

[0031] MPP is an optimization problem where the decision variables are subject to linear constraints and the objective function is quadratic or piecewise linear. MPP has applications in various fields such as process control, robotics, and power systems. According to the constraints introduced above, the coupling variable of the system is the output power, represented by the parameter vector θ = V out denoted. The energy flow vector V is represented as the optimization variable vector x. Therefore, the EH model described by the above constraints becomes a typical MPP problem with the parameter θ = V out . The above optimization problem can be concisely expressed as: min z(x) s.t.Ax + Bθ ≤ C where θ is the parameter vector. The parameter vector θ is selected as the output vector V, the optimization variable vector x is the energy flow vector V, and z(x) is the objective function of the optimization problem. A, B, and C are the coefficient matrices of the above operation constraints. The concept of the proposed equivalent model is as Figure 2 shown.

[0032] The steps to solve the MPP problem are as follows: 1) Optimal partitioning: The first step in solving the MPP problem is to partition the feasible region into subsets, each corresponding to a unique solution. This is achieved by identifying the extreme points of the feasible region. The number of extreme points is equal to the dimension of the problem. For example, a two-dimensional problem has four extreme points determined by the maximum and minimum values along each axis. All the above constraint sets can be divided into active constraints and non-active constraints, respectively represented as: A A x * (θ) + B A θ = C A A I x * (θ) + B I θ < C I 2) Critical Region: The critical region represents the set of parameter values θ that share a common "optimal partition", which is the region where the active constraints change. These boundaries are typically determined by solving a system of linear equations that represent the active constraints when transitioning from one polyhedron to another. Each θ corresponds to a specific combination of active and inactive constraints.

[0033] For any θ i ∈ Θ, the optimal solution is obtained by affine transformation: Substituting the above equation into the inactive constraints, we get That is: Therefore, the critical region associated with θ i is defined as: It can be expressed as: where: As Figure 3 shown, the present invention provides an innovative algorithm for calculating the steady-state safety region in EH, approximating a high-dimensional polyhedron as a low-dimensional polyhedron, reducing the computational burden while retaining the original time coupling of various cycle combinations.

[0034] The space of output variables is high-dimensional, which means that exploring the high-dimensional safety region with time coupling in the exploration space will cause the computational burden to increase exponentially with the increase in the number of time periods. There are methods to exhaustively list all possible combinations of time periods to explain the original time coupling. However, this will result in a large number of combination enumerations, expressed as M = where represents the number of combinations, j represents the number of selected time periods, n out represents the number of system output streams, and T is the total number of time periods. This indicates that the calculation process involves a large number of calculations of low-dimensional polyhedra. For this reason, the present invention discloses a fast approximation method to evaluate the carrying capacity of MES. The fast approximation method is based on the concept of decomposing an exact high-dimensional polyhedron into the intersection of several low-dimensional polyhedra. These low-dimensional polyhedra approximate the safety region by reflecting the time coupling between them. The coupling between specific time periods or between power, heat, and cooling outputs. Therefore, the method first needs to select the i-th low-dimensional polyhedron θ i involving the number of time periods. Thus, the high-dimensional polyhedron θ0 can be decomposed into low-dimensional polyhedra θ iSuppose our method needs to combine k arbitrary time periods, where k is a natural number not greater than nT. In this case, the total number of required low-dimensional polytopes is given by: Compared with the original method, the overall computational burden is significantly reduced.

[0035] Meanwhile, based on the above algorithm, the present invention also proposes a precision measurement method and a precision index based on the difference between a high-dimensional safety region and its approximation to evaluate the precision loss brought by the proposed approximation.

[0036] Regarding the constraint Θ = {V out ∈R K |V ∈ (1)-(12)} which is the projection of the safety region and is essentially the constraint G = {(V out , V) ∈ R K ×R B |V ∈ (1)-(12)} represented by the constraint space R×R projected onto the output space R, that is, the former is an approximation of the latter. Therefore, using the latter as a benchmark can effectively measure the precision loss of the proposed approximation without the need to calculate the exact safety region, thus avoiding a significant computational burden.

[0037] It should be noted that the most accurate precision evaluation method is to compare the deviations of all coordinate points on the boundary plane of the feasible region. However, a limitation of this method is that when approximating the original polytope with low-dimensional polytopes, the resulting high-dimensional polytope composed of multiple low-dimensional polytopes does not necessarily have the same number of vertices as the original polytope. In this case, it becomes challenging to find variables that can accurately evaluate the accuracy.

[0038] To overcome the above challenges, a hyperplane D T V out , where D = [d1, d2,..., d K T ∈R K , can be constructed and two optimization problems related to the above two constraints are designed: The vector D represents a specific combination of outputs in a given problem. The objective function g in both problems shows the maximum output that can be achieved under a specific combination of outputs described by the given vector D. Among them, (15) represents the region determined by Θ = {V out ∈R K |V ∈ (1)-(12)}, and (16) represents the constraint G = {(V out , V) ∈ R K ×R B ​The region determined by {V ∈ (1)-(12)}. We can quantify the accuracy of the proposed approximate safety region by comparing the vertices obtained from the constraint G = {(V out , V) ∈ R K ×R B |V ∈ (1)-(12)} and those obtained from the constraint Θ = {V out ∈ R K |V ∈ (1)-(12)}. This comparison is crucial for ensuring the accuracy of the approximation method.

[0039] The Monte Carlo method is used to randomly generate elements in D for randomly comparing the vertices within the polytopes associated with the above two constraints. This method ensures that the elements generated in D are random and within the range [-1, 1], where the norm of the vectors in this method equals 1.

[0040] To demonstrate the effectiveness of the proposed algorithm, this embodiment conducts a case study on two test systems.

[0041] For Figure 5 shows the MES used in Test System 1, while Figure 6 shows the MES used in Test System 2. Compressed electric refrigerator group is abbreviated as CERG, combined heat and power unit is abbreviated as CHP, absorption electric refrigerator group is abbreviated as WARG, auxiliary boiler is abbreviated as AB, electric heat pump is abbreviated as EHP, cold storage is abbreviated as CS, and heat storage is abbreviated as HS. The basic configuration is shown in Table Ⅰ. Table 1 Test System Parameters The first embodiment considered is Figure 5 Test System 1 in, which shows only two time periods. During the intraday operation, time coupling is considered by taking into account the ramp rate and storage limits. Within these two time periods, the system has a total of 2×3 = 6 possible output states: V(1), V(2), V(1), V(2), V(1), V(2). In this example, A is a 214×90 matrix, B is a 214×6 matrix, and C is a 214×1 matrix. Due to the nature of the system operation constraints described by equations (1)-(14), the two matrices A, B, and C are sparse matrices.

[0042] For Figure 7 and Figure 8 show the safety regions within the rolling scheduling range of 2 cycles in Test System 1, with and without ramp rate and storage limits. Obviously, when time coupling is considered, the area of the safety region is much smaller. Ignoring time coupling may result in considering some infeasible power combinations as feasible ones, which may pose a significant risk to intraday operation safety.

[0043] In the fast approximation method, in a system characterized by a low-dimensional output variable Vand that includes only two time periods, the method uses multiple low-dimensional polyhedra to construct a high-dimensional polyhedron. Precision exponents ω f，max , ω f，ave , ω v，max and ω v，ave are 10.5%, 0.96%, 9.30%, and 0.74% respectively. Through this strategy, it can accurately approximate the complex load-bearing capacity of MES.

[0044] To prove the efficiency of the proposed method in achieving a better balance between computational accuracy and efficiency compared with existing methods, we selected "2" as the number of time periods for our proposed fast approximation method. Table I Computation times of different methods Table II Precision evaluation of different methods 1) Comparison methods: We compared our proposed method with the following four methods: M0: Using the original multi-period constraint; M1: Adopting the fast approximation method provided by the present invention; M2: Adopting an existing method that ignores time coupling; M3: Extending the method used in M2 to calculate the maximum output power combination of the EH output ports within each time period. 2) Computation time: As mentioned before, due to the significant increase in the number of vertices involved, the combination of time coupling and the accurate identification of high-dimensional safety regions may pose significant computational challenges. In the case of the M0 method, a large number of high-dimensional vertices result in a computation time of more than one hour, as shown in Table II.

[0045] The computation time for M2 is much shorter, with the computation times for the two test systems being 6.42 seconds and 14.07 seconds respectively. The reduction in computation time can be attributed to the removal of the coupling between different time periods in the M2 method, resulting in a significant reduction in the number of vertices of each polyhedron compared with M0. However, the M2 method ignores the coupling between different time periods, resulting in an inaccurate computed safety region over multiple time periods. In addition, the computation time for the M3 method is also favorable. For example, it only takes 9.38 seconds to compute the safety region of test system 2. This is because M3 uses a specific output power combination at the output ports. However, this method also leads to significant inaccuracies.

[0046] For the proposed method, compared with M0, its calculation time is shorter. Although the calculation time of the M1 method may be longer than that of the M2 and M3 methods because multiple polyhedra related to different shortening cycle combinations are calculated, M1 exhibits significantly superior precision performance.

[0047] 3) Computational accuracy: Table III presents the accuracy evaluation results of the two test systems. In these two test systems, the proposed M1 method is significantly superior to the M2 method, and the metrics (ω f，max , ω f，ave , ω v，max and ω v，ave ) are significantly reduced. This indicates that the error caused by the M2 method is significantly reduced. Specifically, in test system 2, our method only achieves an exponential ω f，max of 9.2%, while the exponential ω f，max of the M2 method is significantly higher, at 51.11%, demonstrating the excellent accuracy of our method. In test system 2, compared with the M3 method, the exponential values (ω f，max , ω f，ave , ω v，max and ω v，ave ) of the M2 method are also lower. However, the M3 method lacks in precision because it only focuses on the power transfer ability of the output port. For example, in the same test system 2, the exponential ω f，max obtained by the M3 method is 41.00%.

[0048] The M0 method can accurately capture the safe region, including time coupling, and the corresponding metrics ω f，max , ω f，ave , ω v，max and ω v，ave are zero. However, when constructing polyhedra with a large number of high-dimensional vertices, the M0 method faces a huge computational burden, making it unrealistic to produce results within a reasonable time range (i.e., within one hour).

[0049] As described above, the number of selected time periods is k, where k is a natural number not exceeding T. Here, the influence of different time period selections will be compared.

[0050] The present invention proposes a new method for approximating a high-dimensional Θ feasible region using low-dimensional polyhedra. The computational complexity of this approximation is affected by the dimension and the combination of selected time periods. To evaluate the required calculation time, Table IV provides a table showing the relationship between different selected time periods and the calculation time of the system response. It is expected that the computational complexity of the proposed method will increase with the increase in the number of selected time periods. Table I Calculation times for different selected time periods Table II Precision for Different Selection Time Periods Table IV shows that when two time periods are selected, the calculation times of the two test systems are 10.96 seconds and 21.23 seconds respectively. This difference is due to the limitation of the M2 method, which only accommodates 2-dimensional decoupled safety regions and combination numbers. It should be noted that it can be observed and indicate a direct correlation between the calculation time and the number of selected time periods. As shown in Table IV, when the selected time periods are 3 and 4, the calculation times of Test System 1 are 34.72 s and 2018.89 s respectively, with a difference of 58 times. At the same time, in terms of the number of combinations, and show that the difference between them is only 3.75 times. Therefore, the dimension of the low-dimensional polyhedron, that is, the selected time period, has a significant impact on the calculation time.

[0051] Figure 9 and Figure 10 demonstrate the effectiveness of the proposed method in approximating high-dimensional safety regions. In addition, the reduced-dimensional safety region shows a high degree of precision in approximating the original high-dimensional safety region. Figure 9 The reduced-dimensional safety region is illustrated by showing the power safety regions for time periods 1, 2, and 3. In addition, Figure 10 shows the approximate projections of the power safety regions for time periods 1, 2, and 3 when two time periods are selected. These visual representations clearly show that the reduced-dimensional safety region provides an accurate approximation of the high-dimensional safety region. In addition, the accuracy evaluation results given in Table V verify that selecting two time periods can achieve the desired accuracy.

[0052] The present invention aims to solve the problem of evaluating the MES carrying capacity in a systematic manner and demonstrate its potential for improving MES planning. First, the EH analysis method is adopted, which is a successful method for MES operation analysis. The proposed method uses the EH safety region to evaluate the MES carrying capacity by first defining the attributes of the MES, and then establishing a mathematical representation by using a standardized matrix modeling method. To solve the problems of high dimension and time coupling, the present invention provides a multi-parameter programming algorithm for approximating a high-dimensional polyhedron as a low-dimensional polyhedron. Compared with other existing methods, the proposed method can achieve a better balance between calculation accuracy and efficiency, especially for MES with multiple devices and time coupling. It should be emphasized that for small-scale MES or systems with time decoupling, the direct application of the MPP algorithm is sufficient. The method provided by the present invention can identify key components and bottlenecks in the system and contribute to formulating effective MES planning and operation strategies, making contributions to the effective optimization of MES and its successful integration with the existing market.

[0053] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code. The solutions in the embodiments of the present application can be implemented in various computer languages. For example, object-oriented programming languages such as Java and interpreted scripting languages such as JavaScript, etc.

[0054] The present application is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to the embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram, as well as the combination of flows and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0055] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing devices to work in a specific manner, such that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means implement the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0056] These computer program instructions can also be loaded onto a computer or other programmable data processing devices, such that a series of operation steps are executed on the computer or other programmable devices to generate a computer-implemented process. Thus, the instructions executed on the computer or other programmable devices provide steps for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0057] Although the preferred embodiments of the present application have been described, those skilled in the art can make additional changes and modifications once they know the basic creative concepts. Therefore, the appended claims are intended to be construed as including the preferred embodiments and all changes and modifications falling within the scope of the present application.

[0058] Obviously, those skilled in the art can make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalent technologies, this application is also intended to include these modifications and variations.

Claims

1. A systematic method for evaluating the load-carrying capacity of a multi-energy system in an assessment area considering time coupling, which is used to evaluate the load-carrying capacity of the MES, is characterized in that: The EH method is used to model the system, and the EH steady-state security region is introduced as a means to evaluate the system's carrying capacity; when calculating the security region, a multi-parameter programming algorithm is adopted, and its mathematical model is established using the standardized matrix modeling method; in the multi-parameter programming algorithm, an innovative algorithm for calculating the steady-state security region in EH is adopted; and based on the above innovative algorithm, the accuracy loss brought by the proposed approximation value is evaluated.

2. The system method for evaluating the carrying capacity of a multi - energy system in an assessment area considering time coupling according to claim 1, characterized in that: EH has multiple inputs and outputs, including four stages: energy input, conversion, storage, and output; Energy, power grid, and gas network serve as the inputs of EH, while the distribution network, heat network, and load are the outputs of EH. In the analysis of the energy network, the matrix representation of the directed graph is used to establish a mathematical model to describe the relationship between the energy flows within the network. The elements in EH, such as components and energy transmission pipelines, directly correspond to the nodes and branches of the network graph theory. MES uses the directed graph for mathematical expression, and the matrix used to describe the relationship between the energy flow and components within the EH system; The input / output port-branch incidence matrix is: The energy conversion characteristics of node g are captured by the converter characteristic matrix H. Assuming that node g has K ports and its characteristics are described by the P equation, the corresponding dimension of matrix H is P×K: where η represents the efficiency of the energy conversion process s; Based on the calculation of the first two types of operation matrices, the energy conversion matrix is obtained. The port-branch incidence matrix and the converter characteristic matrix represent the input-output relationship of energy devices and the energy conversion relationship in EH respectively. Since the sub-matrix of the converter characteristic matrix is a diagonal matrix and its dimension is equal to the dimension of the port-branch incidence matrix, the product represents a set of energy conversion relationships of a specific energy device in the energy input-output process in EH. Accordingly, the energy conversion matrix of node g is derived as follows: Aggregating the energy conversion matrices of all nodes, the energy conversion matrix of EH is obtained: Z = [Z1; Z2;...; Z g ​ Based on the predefined matrix, the energy conversion equation is given as: ZV = 0 where Z is defined in the energy conversion matrix, V is the energy flow vector containing the energy flows of all branches, the energy input flow is represented by V, and V represents the energy output flow. The input / output incidence matrix is as follows: The matrix form of the entire energy flow equation is expressed as: Among them, C in and C out respectively represent the input and output incidence matrices, V represents the energy flow vector containing the energy flows of all branches, V in and V out respectively represent the energy input and output flows; The energy flow relationship within EH, as shown in the above energy flow equation matrix, only captures the real-time energy balance and does not consider energy storage components. Energy storage devices such as batteries, compressed air energy storage, chilled water storage, thermal energy storage, and gas storage are not reflected; To comprehensively establish the network graph theory mathematical model of EH, energy storage nodes need to be included. The energy flow entering the energy storage node is regarded as the charging variable, while the energy flow leaving the energy storage node is regarded as the discharging variable. The charge-discharge equation of the energy storage device is expressed as: where is g th the number of output ports of the storage component, S(t) is the state of charge (SoC) of the storage device within the time period t, ΔS is the charge and discharge energy of the g storage component, and the node feature matrix of the storage device is defined as: The port-branch incidence matrix of the storage device is enhanced to: The mathematical expression for the SoC of the storage device; The energy conversion matrix of the storage device is represented as: For the storage node, the energy balance equation is derived from the energy charge-discharge equation as: Among them denotes all the elements in that represent the energy flow in V, -I represents the identity matrix of dimension units, and ΔS represents all the elements in ΔS g in; The energy flow obeys the branch energy constraint. In addition, the capacity of each converter conforms to the limit of the input or output branch: Among them and are the input and output capacities of all nodes, respectively, V B，min and V B，max represent the upper and lower branch energy constraints, respectively; In addition to the SoC limitations of the energy storage devices discussed above, the MES also faces ramp rate limitations imposed on the energy storage devices and other devices. These constraints lead to strong time coupling within the system. The change in the node output between two consecutive time intervals should comply with the ramp rate constraint, i.e., Among them and represent the ramp-down and ramp-up rates during the time period t. When the storage capacity is large enough or the SoC is moderate, the SoC constraint of the energy storage device is relaxed, but the ramp rate limit cannot be relaxed; The steady-state security region of a power system is related to power flow equations and power constraints. Similarly, for an EH, the security region is related to device operation and safety constraints. The output port of the EH is the key to evaluating the load capacity of the MES. The security region represents the set of energy output vectors that satisfy the relevant constraint conditions and defines the hyperspace R K ×R B , constructed by v out (K×1) and V(B×1), where B is the number of branches in the EH considering time coupling. All of the above constraint conditions jointly form a polyhedron G in the hyperspace R K ×R B : G = {(V out , V) ∈ R K × R B | V ∈ (1)-(12)} where (1)-(12) represent the above-mentioned first 12 constraint conditions. To obtain the safe region θ, the polyhedron G in R×R is projected onto R.

3. A systematic method for evaluating the carrying capacity of a multi - energy system in an assessment area considering time coupling, as claimed in claim 2, wherein: A multi-parametric programming algorithm is adopted to effectively solve the safe region calculation problem. This algorithm is called MPP, and can also be referred to as: a multi-parametric programming algorithm. The steps to solve the MPP problem are as follows: 1) Optimal partitioning: The first step in solving the MPP problem is to partition the feasible region into subsets, each corresponding to a unique solution. This is achieved by identifying the extreme points of the feasible region. The number of extreme points is equal to the dimension of the problem. All the above constraint sets are divided into active constraints and non-active constraints, which are represented as: A A x * (θ)+B A θ = C A A I x * (θ) + B I θ < C I 2) Critical region: The critical region represents the set of parameter values θ that share a common optimal partition. This is the region where the active constraints change, and the boundaries are determined by solving a system of linear equations that represent the active constraints when transitioning from one polyhedron to another. Each θ corresponds to a specific combination of active and non-active constraints; For any θ i ∈ Θ, the optimal solution is obtained by an affine transformation: Substituting the above equation into the non-active constraints, we get That is: Therefore, the critical region associated with θ i is defined as: Represented as: Where:

4. A systematic method for evaluating the carrying capacity of a multi - energy system in an assessment area considering time coupling, as claimed in claim 3, wherein: In the multi-parameter programming algorithm, an innovative algorithm for calculating the steady-state safety region in the EH is adopted. In the innovative algorithm, a high-dimensional polytope is approximated as a low-dimensional polytope, reducing the computational burden while retaining the original time coupling of various cycle combinations. The space of output variables is high-dimensional, meaning that exploring the high-dimensional safety region with time coupling in the exploration space will cause the computational burden to increase exponentially with the increase in the number of time periods. In the existing methods, all possible time period combinations are enumerated as much as possible to explain the original time coupling. However, this will result in a large number of combination enumerations, expressed as where represents the number of combinations, j represents the number of selected time periods, and n out represents the number of system output streams, and T is the total number of time periods; this indicates that the calculation process involves a large number of calculations of low-dimensional polytopes. Therefore, a fast approximation method is adopted to evaluate the carrying capacity of the MES. The fast approximation method is based on the concept of decomposing an exact high-dimensional polytope into the intersection of several low-dimensional polytopes. The low-dimensional polytopes approximate the safety region by reflecting the time coupling between them, the coupling between specific time periods or between power, heat, and cooling outputs. First, the i-th low-dimensional polytope θ i needs to be selected, which involves the number of time periods. The high-dimensional polytope θ0 is decomposed into low-dimensional polytopes θ i . Assuming that k arbitrary time cycles need to be combined, where k is a natural number not greater than nT, in this case, the total number of required low-dimensional polytopes is given by the following formula: Compared with the original method, the overall computational burden is significantly reduced.

5. A systematic method for evaluating the carrying capacity of a multi - energy system in an assessment area considering time coupling, as claimed in claim 4, wherein: Based on the above innovative algorithm, a precision measurement method and precision index based on the difference between a high-dimensional safety region and its approximation are obtained to evaluate the precision loss brought by the proposed approximation; where: Regarding the constraint Θ = {V out ∈R K |V ∈ (1)-(12)} is the projection of the safety region, and the constraint G = {(V out , V) ∈ R K ×R B |V ∈ (1)-(12) to the output space R, that is, the former is an approximation of the latter; using the latter as a benchmark to measure the precision loss of the proposed approximation without the need to calculate the exact safety region, thus avoiding a significant computational burden; the most accurate precision evaluation method is to compare the deviations of all coordinate points on the boundary plane of the feasible region; however, a limitation of this method is that when approximating the original polytope with low-dimensional polytopes, the resulting high-dimensional polytope composed of multiple low-dimensional polytopes does not necessarily have the same number of vertices as the original polytope. In this case, construct a hyperplane D T V out , where D = [d1, d2, …, d K T ∈R K , and design two optimization problems related to the above two constraints:​ The vector D represents a specific combination of outputs in a given problem, and the objective function g in both problems shows the maximum output achievable in the case of a specific combination between the outputs described by the given vector D; by comparing the vertices obtained from the constraint G = {(V out , V) ∈ R K ×R B | V ∈ (1)-(12)}, the accuracy of the proposed approximate safety region and that obtained from the constraint Θ = {V out ∈ R K | V ∈ (1)-(12)} are quantified; the Monte Carlo method is adopted to randomly generate the elements in D, and the vertices within the polytope related to the above two constraints are randomly compared; it is ensured that the elements generated in D are random and within the range of [-1, 1], where the norm of the vector in this method is equal to 1. Among them, (15) represents the region determined by Θ = {V out ∈ R K | V ∈ (1)-(12)}, and (16) represents the region determined by the constraint G = {(V out , V) ∈ R K ×R B | V ∈ (1)-(12)}.

Images