Electric heating comprehensive energy system unit combination method based on interior point bender decomposition

Abstract

The invention relates to the technical field of optimal scheduling of an integrated energy system, and discloses an electric heating integrated energy system unit combination method based on interior point Benders decomposition. The method comprises the following steps: firstly, constructing a mixed integer linear programming model, determining decision variables such as power unit start-stop / scheduling, thermodynamic system output and heat supply network temperature, aiming at minimizing the total operation cost of the system, integrating multiple constraint conditions, analyzing a pipeline temperature conduction partial differential equation by adopting a characteristic line method, and improving the temperature dynamic description accuracy; and the model is solved through an improved Benders decomposition algorithm combined with an interior point method, and parallel solving of the main problem and the sub-problems is achieved. Through IEEE standard system simulation verification, the method can realize electric-thermal system collaborative optimization scheduling, adapts to wind power uncertainty, guarantees system stability, reduces operation cost, and improves energy utilization efficiency.

Description

Technical Field

[0001] This invention relates to the problem of optimized scheduling technology for integrated energy systems, specifically to a method for combining units in an integrated electrothermal energy system based on interior-point Bendes decomposition. Background Technology

[0002] Against this backdrop, integrated energy systems have attracted widespread attention due to their significant advantages in improving energy efficiency and reducing carbon emissions. By integrating multiple energy forms such as electricity, heat, and natural gas, integrated energy systems achieve energy complementarity and optimized allocation, thereby improving the flexibility and reliability of the energy system. As an important component of integrated energy systems, the electric-thermal integrated energy system, through coupling the power system and the heating system, not only achieves efficient energy production and distribution but also effectively utilizes renewable energy sources such as solar and wind power, which is of great significance for promoting the rapid economic development of the energy system.

[0003] Unit combination is a core issue in energy system operation, and its research in integrated power-heat (EPC) energy systems is of crucial value because it directly relates to the system's cost-effectiveness and environmental impact. By optimizing unit combination, EPC energy systems can achieve efficient energy distribution, ensuring that energy demand is met while minimizing environmental impact. This optimization involves the production and distribution of electricity and heat, and also covers the effective use of renewable energy sources such as solar and wind power. The intermittency and unpredictability of these energy sources pose challenges to system scheduling. Through rational planning and scheduling, EPC energy systems can flexibly adjust energy supply according to changes in energy demand in different seasons and time periods, achieving efficient energy utilization through continuous technological innovation and system optimization. However, the optimized operation of EPC energy systems faces many challenges, among which unit combination optimization is one of the core issues. Unit combination optimization directly relates to the system's cost-effectiveness and environmental impact, and its optimization is crucial for improving the system's economy and reliability. Developing effective algorithms to optimize unit combination is of great significance for improving the operating efficiency and environmental friendliness of EPC energy systems. Therefore, optimizing it is key to improving economic efficiency and reliability. Summary of the Invention

[0004] This application provides a unit combination method for an integrated electric-thermal energy system based on interior-point Bendez decomposition, which is used to solve the technical problems of unit combination in existing integrated electric-thermal energy systems regarding system cost-effectiveness and environmental impact. Its optimization is crucial to improving the system's economy and reliability. Developing effective algorithms to optimize unit combination is of great significance for improving the operating efficiency and environmental friendliness of integrated electric-thermal energy systems.

[0005] This invention addresses the unit combination problem that considers complete temperature dynamics by constructing a mixed-integer linear programming (MILP) model.

[0006] Therefore, the present invention provides a method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition, characterized in that it includes: Step 101: Construct a mixed-integer linear programming model to clarify the decision variables such as power system unit start-up / shutdown / dispatch, thermal system output, and heating network temperature. With the goal of minimizing the total system operating cost, integrate multiple constraints such as power balance, spinning reserve, ramp rate, and thermal network transmission. In the dynamic modeling of the thermal network, the method of characteristics is used to analyze the partial differential equation of temperature conduction in the pipeline, avoiding discretization errors and improving the accuracy and computational efficiency of the dynamic temperature description. Step 102: Use the improved Pendes decomposition algorithm combined with the interior point method to solve the mixed integer linear programming model and obtain the final result.

[0007] The aforementioned method for combining generator units in an integrated electrothermal energy system based on interior-point Bendez decomposition, wherein the power decision variables of the mixed-integer linear programming model include the combination state of the generator units. and scheduling ;in, Indicates the unit In time The running status, and They are respectively in Startup and shutdown status, It is a generator set In time Power output, It is a wind turbine In time Power output, and They are the generator sets In time Upward and downward rotational reserve capacity; thermal decision variables encompass the thermal output of the combined cogeneration unit and the supply / return water temperature of the district heating network; thermal output represents the heat power provided by the unit while meeting electricity demand; the supply / return water temperature of the district heating network directly affects the efficiency of heat transmission and distribution; Unit Commitment is abbreviated as UC.

[0008] The method for combining units in an integrated electric and thermal energy system based on interior-point Bends decomposition, wherein the objective function aims to minimize the total operating cost of the system, encompassing the operating costs of the power system and the thermal system, as well as the additional costs due to the uncertainty of wind power; Minimize total operating costs: (1) In the cost structure, the first item represents the operating cost of non-CHP units, the second item represents the cost of wind curtailment penalty, and the third item represents the operating cost of CHP units. The English for combined heat and power unit is: Combined Heat and Power, CHP; the English for district heating network is: District Heating Network, DHN. The operating cost of a non-combined thermal power unit is defined as: (2) in, This represents the piecewise linear power generation cost, which approximately describes the nonlinear relationship between fuel cost and unit output; This represents the startup cost, when the unit is in time. A fee is required upon startup, the amount of which depends on the unit type and startup mode. This represents the shutdown cost, when the unit is in time. When shutting down, additional costs need to be considered, including equipment cooling and maintenance. Indicates no-load cost; The formula for calculating the cost of wind curtailment penalty is: (3) in, Indicates wind turbine In time The actual output, Indicates wind turbine In time The planned output, Indicates the wind curtailment penalty coefficient; Let T represent the set of wind turbine units, and let T represent the time set.

[0009] The operating costs of CHP units are similar to those of non-CHP units, including the cost of power generation. Defined as: (4) Among them, piecewise linear power generation cost Describes the fuel consumption of the CHP unit at different output levels; startup cost item. This reflects the additional costs incurred when the CHP unit starts up at time t; the cost item is closed. This covers the costs incurred when the unit is shut down at time t; no-load cost item. This represents the energy consumption and maintenance costs of the CHP unit to maintain basic operation under no-load conditions.

[0010] The unit combination method for the integrated electric and thermal energy system based on interior-point Bends decomposition, wherein the mixed-integer linear programming model is subject to constraints imposed by the operation of the power and heat network, the constraints being: Power balance constraint: Total power generation and total load must remain balanced within each time period. (5) This constraint ensures that, in each time period, the total power generation of all generating units, including conventional units, CHP units, and wind turbines, is balanced with the total power load of the system. in, This represents the actual power output of a conventional or CHP unit at time t. This is the operating status variable of the generator set; it is set to 1 when the generator set is running and 0 otherwise. The rated power of the unit; for wind turbines, the power generation is determined by... Indicates; the right side This is the total power load demand of the system at time t; this constraint is a basic requirement for the operation of the power system, ensuring the real-time balance of power supply and demand. Spinning reserve constraint: Non-CHP units provide the required spinning reserve capacity. (6) (7) (8) (9) (10) Spinning reserve constraints ensure that the system has sufficient reserve capacity to maintain stable operation in the face of load fluctuations or sudden events; among which, Indicates the unit In time The upward rotation of the spare capacity, Indicates the unit In time The downward rotation of the spare capacity, Indicates the unit The maximum starting ramp rate, Indicates the unit The maximum rate of ... It is a generator set Minimum boot time, This indicates the total backup demand for the entire system rotating upwards. This indicates the total backup demand for the entire system rotating downwards. This represents a set of thermal power units; these constraints, by limiting the relationship between the generating output and reserve capacity of the units, ensure that the system can meet the requirements of spinning reserve under different operating conditions. Specifically, the first three inequalities limit the sum of the unit's maximum power output and reserve capacity based on the unit's start-up and shutdown status and operating mode, preventing the unit from exceeding its operating limits; the last two inequalities, on the other hand, specify the minimum total amount of upward and downward rotational reserves in each time period at the system level, ensuring that the system has the ability to cope with load changes. Ramp-up constraint: Within a single time period, the incremental power generation of a thermal power unit is limited by its ramp-up capability. (11) The ramping constraint limits the rate of change in power generation output of thermal power units within adjacent time periods, i.e., the ramping capability of the units; among which, Indicates the unit The upper limit of the upward climbing rate, Indicates the unit The upper limit of the downward climbing rate, This represents the set of CHP units; this constraint reflects the limitation of the physical characteristics of thermal power units on the speed of output adjustment, avoiding equipment damage or operational instability caused by rapid output adjustment of units, and also helps the economical scheduling of the system, preventing increased operating costs due to frequent adjustments of unit output. Wind power output constraints: The power generation of a wind farm is limited by the available wind power. (12) This constraint specifies the power output range of the wind turbine in each time period; the lower limit is 0, indicating that the wind turbine can stop generating electricity; the upper limit... This is determined by the wind conditions at the time, reflecting the renewable energy characteristics and uncertainties of wind power; this constraint ensures that the power generation of wind turbines does not exceed their maximum power generation capacity in a specific period, which is in line with the actual operation of wind power.

[0011] Unit state constraints (13) Generator set state constraints describe the dynamic relationships between the unit's operating state, startup state, and shutdown state; among them, Indicates the unit In time The running status, Indicates the unit In time The startup state, Indicates the unit In time The equation shows that the change in the unit's operating state during the current period is equal to the difference between its start-up and shutdown states, thus ensuring the continuity and logical consistency of the unit's state and providing a basis for accurately modeling the unit's start-up and shutdown process. Minimum start-stop time constraints (14) (15) The minimum start-up and shutdown time constraint ensures that once the unit starts up or shuts down, it will remain in that state for at least a certain period of time; the first inequality stipulates that, within time... Before During the time period, the unit The number of times a unit can be started cannot exceed its operating state; that is, if the unit is currently in an operating state, then the number of times it has been started in the past... The first inequality ensures that the unit is started at least once within a certain time period; the second inequality ensures that if the unit is currently shut down, then it was started at least once in the past. The unit can be shut down no more than once within a given time period; these constraints reflect the physical requirements for hot and cold starts of the unit, avoiding damage to the equipment caused by frequent start-stop cycles, while also helping to reduce operating costs and improve system reliability.

[0012] Variable bounds constraints (16) (17) The upper and lower bound constraints of the variables clearly define the range of values ​​for each decision variable. This applies to the operating status of the unit. It can only take the value 0 or 1, indicating whether the unit is running at time t; startup and shutdown status variables. and The value of is between 0 and 1, but in actual operation it is usually 0 or 1 to ensure the integer programming characteristics of the model; for the rotated spare variable and The values ​​are limited to 0 to the unit's maximum upward and downward ramp rates, ensuring the rationality and feasibility of the reserve capacity.

[0013] Network constraints: The power flow of transmission lines should not exceed the transmission capacity. Equation (18) represents the nodal power equation, Equation (19) represents the branch power flow equation, and Equation (20) defines the phase angle of the reference bus; where, It is with nodes The set of all directly connected nodes It is a node The upper connected conventional generator set, It is a node The collection of combined heat and power units connected to the top, It is a node The upper-connected wind turbine assembly, For nodes The voltage phase angle at time t, For nodes The voltage phase angle at time t, This represents the voltage phase angle of the reference node in the network at time t. Let (i,j) be the reactance of the line. For the line ( i , j Maximum transmission capacity, This is the set of all routes; (18) (19) (20) Network constraints ensure that the power system meets the physical limitations and safety requirements of power transmission during operation; the node power equations link the node voltage phase angle difference with the power transmitted through the lines, reflecting the basic laws of power flow in the power system; the branch power flow equations limit the power transmission capacity of each transmission line, prevent line overload, and ensure the safe and stable operation of the system; the setting of the reference bus phase angle provides a benchmark for the voltage phase angle of the entire system, facilitating calculation and analysis; these constraints together constitute the operating framework of the power network, ensuring the efficiency and reliability of power transmission.

[0014] The unit combination method for the integrated electrothermal energy system based on interior-point Bends decomposition, wherein the dynamic modeling method of the thermal network in step 101 includes the DHN model and its complete analysis method: DHNs are two-layer networks that integrate supply and return water networks. Heat energy is generated from the heat source, transported through the supply pipeline, distributed at the heat exchange station, and then flows to the heat users. After heat exchange at the user end, the water flows back to the heat source through the return water network. A typical DHN model includes both hydraulic and thermal components. The mass regulation mode of DHNs under constant mass flow rate meets the heat load through temperature regulation. Temperature is defined as relative temperature, that is, the absolute temperature difference between the water flow and the ambient temperature. The mass model formula for DHN under mass conditioning is as follows: (twenty one) (twenty two) (twenty three) (twenty four) in, This represents the longitudinal distance of a point in the pipeline from the pipeline inlet, i.e., the spatial coordinates along the pipeline. Indicates time, Indicates the first j The temperature of the pipe, Indicates the first j The water flow velocity in the pipe This indicates the specific heat capacity of water. Indicates the first j Mass flow rate of the pipeline Indicates the first j The heat loss coefficient of the pipeline, Represents a node thermal power, Represents a node Mass flow rate Represents a node Water supply pipe temperature Represents a node Return water pipe temperature Describes the set of pipes flowing into node j. Describes the set of pipes that flow out of node j. This represents the temperature at the outlet of pipe b, where L is the length of the pipe. This represents the mass flow rate of pipe b. This indicates the temperature at the inlet of pipe b. This represents the mixing temperature at node k. Represents a set of pipes. Represents a set of nodes. Let represent the set of pipelines that start at node k and end at node i.

[0015] The temperature conduction equation represents the temperature distribution along a pipe, and this equation is constructed based on the following assumptions: 1) The heat loss of the pipe cross-section is considered constant; 2) Ignore heat conduction between adjacent water flows; The equation models the nodes in the DHN as heat exchangers, thus allowing a linear correlation between mass flow rate and temperature; the nodes follow the law of energy conservation, and the mixing temperature of a node is a linear weighted sum of the outlet temperatures of the pipes injected into that node; finally, the node mixing temperature of the pipes leaving the node is the same as the inlet temperature.

[0016] The aforementioned method for combining units in an integrated electrothermal energy system based on interior-point Pendes decomposition involves a heat conduction equation in the pipes that is a typical first-order linear and quasi-linear partial differential equation. Quasi-linear partial differential equations (PDEs) are commonly used to describe vibration and wave phenomena and their dynamic processes; their representative form is as follows: (25) Where a, b, and c are the coefficients of the first-order linear hyperbolic partial differential equation. Specifically, a corresponds to the rate of change of temperature in space, b corresponds to the rate of change of temperature in time, and c corresponds to the temperature distribution itself. Let x be the state variables of a first-order linear hyperbolic partial differential equation, where x is the spatial coordinate and t is the time variable. This represents the temperature distribution function, i.e., the temperature at pipe location x and time t. These are the initial conditions for hyperbolic partial differential equations. This represents the initial temperature distribution at location x.

[0017] The solution plane of this equation consists of two-dimensional vectors in space and time, and is also known as the Cauchy problem; The core idea of ​​the Characteristic Line Method (CLM) is to find characteristic lines. Along these characteristic lines, the original PDE can be transformed into an ordinary differential equation (ODE). From this perspective, the solution plane can be considered as being formed by the initial condition curves on... The surface formed by translating along the feature line.

[0018] The aforementioned method for combining electrothermal integrated energy system units based on interior-point Bends decomposition includes the following steps in traditional CLM: The differential of u with respect to t is expressed as: (26) Substituting the above expression into the pipe temperature conduction equation, we observe that when... When this happens, equation (25) can be converted to ODE; therefore, The solution is defined as the characteristic line of equation (25), expressed as: (27) Where c is a constant determined by the initial conditions; Replacing the partial derivative with the total differential, the reconstructed ODE is expressed as: (28) The analytical solution to the above equation can be obtained from the coefficients. Different expressions and initial conditions get; To solve the Cauchy problem using CLM, a set of initial conditions is needed as the boundary of the solution plane; in other words, the PDEs that can be solved by traditional CLM must be unilaterally bounded along the feature line direction. Considering a practical scenario, where... and Both are bounded, and the bilaterally bounded hyperbolic PDE is extended as follows: (29) in, Describes the boundary conditions for hyperbolic partial differential equations. This indicates the temperature setpoint at the heat source end of the pipeline.

[0019] Unlike traditional CLM, its solution is determined by both initial and boundary conditions, expressed as: (30) Determined by initial conditions, this indicates the temperature distribution caused by the initial heat storage in the pipeline; Determined by the boundary conditions, it represents the temperature propagation distribution caused by real-time heating from the heat source; the solution plane can be considered as the initial condition plane along... The result of the intersection of the surface formed by the movement and the surface formed by the movement of the boundary conditions; (31) (32) in, Right now , This represents the temperature distribution caused by the initial heat storage in the pipe at the initial time t=0. Right now , This indicates the temperature propagation distribution at the beginning of the pipeline, x=0, where heat is supplied in real time by a heat source; UC: short for Unit Commitment.

[0020] The unit commitment method for the integrated electric and thermal energy system based on interior-point Bends decomposition, wherein a deterministic optimization method is used to construct the integrated electric and thermal energy system (UC-CEHN) model, as specifically described below: , , This model optimizes the system operation using the MILP method, comprehensively considering the start-up and shutdown plans of generator units, power output scheduling, and coordination of heat supply. The constraints in the model cover the power system's power balance, spinning reserve, ramping limits, wind turbine output range, generator unit state transitions, and minimum start-up and shutdown times. At the same time, relevant constraints of the heat network are introduced to ensure the stable operation of the integrated electric-heat energy system. To address the complex problem of optimizing an integrated electricity-heat energy system, a Pendesic decomposition technique was introduced to construct a distributed decision-making framework between the grid operator and the district heating operator. Within this framework, the PGO (Power Grid Operator) is primarily responsible for solving the Unit Commitment Master Problem (UCMP) and checking the feasibility of grid network constraints, while the DHO (District Heating Operator) focuses on checking the feasibility of DHN (District Heating Network) constraints. This decomposition-coordination model not only improves computational efficiency but also protects the information privacy of both networks to some extent. The abbreviation for grid operator is PGO, and the abbreviation for district heating operator is DHO. 1. Main UC problem: Provide a combination and scheduling plan for power generation equipment, while minimizing the total operating cost in (1) while satisfying the constraints in (2)-(20), (21)-(24) and the feasible cutting plane provided by the feasibility check; As a MILP, the main UC problem can be solved by the Lagrange relaxation method or the MILP solver. 2. Network security checks: Solving these sub-problems is to check the plans of a given crew. Feasibility of applying network constraints; since the network constraints at different times are independent, the subproblems representing these constraints can be solved in parallel. Each subproblem is described as follows: (33) (34) (35) in, The optimal value of the network security inspection subproblem represents the total number of violations of power flow constraints under a given unit schedule in time period t. This represents the planned active power output vector of the generating units during time period t; This represents the planned start-up and shutdown status vector of the generating units during time period t; This is the first relaxation variable, representing a virtual positive power injection used to compensate for the power deficit at this node. The second slack variable represents the virtual negative power absorption, used to absorb the power excess at the node. The two slack variables are used to ensure that the node power balance equation is always solvable. if For a given unit plan Network constraints are feasible; if Then, add the following Pendesic cut plane to UCMP: (36) in, This represents the vector of unit output variables for time period t. This represents the vector of unit start-up and shutdown status variables for time period t.

[0021] 3. DHN Feasibility Check Sub-problem: The purpose of solving this sub-problem is to check the plan for a given unit. The feasibility of the following DHN operating constraints; its calculation formula is as follows: (37) (38) Constraints (2)-(20), (21)-(24); in, The optimal value of the DHN feasibility check subproblem is represented by , and the total number of violations of thermal network constraints is represented by , under a given unit plan. Represents the thermal output planning vector. This represents the planned start-up and shutdown state vector of the thermal power unit. It is the third relaxation variable, representing a virtual positive heat power injection used to compensate for the heat power deficit at this node. The fourth relaxation variable represents the virtual negative heat power absorption, used to absorb the excess heat power at this node. The two relaxation variables are used to ensure that the heat power balance equation of thermal node j is always solvable at time t. Represents the set of all nodes in a heat network. This represents the set of all heat sources supplying heat to node j at time t. The minimum thermal output coefficient of heat source g, This represents the mass flow rate of node j at time t. This indicates the specific heat capacity of water. This represents the water supply temperature setpoint for node j at time t. This represents the setpoint for the return water temperature at time t for node j.

[0022] if For a given unit plan DHN constraints are feasible; if Then the following Pendesic cut plane will be fed back into the main problem: (39) in, Represents the vector of thermal output variables. This represents the start-up and shutdown state variable vector of the thermal power unit.

[0023] Through the detailed explanation of the above model and sub-problems, it is known that the optimization model of the integrated electric-thermal energy system not only needs to consider the operating characteristics and cost factors of the power system, but also needs to meet the physical constraints and operating requirements of the thermal network.

[0024] The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition, wherein the specific problem of combining units based on improved Bends decomposition is as follows: 1. Traditional Pendesic decomposition algorithm: The Pendes decomposition algorithm decomposes the original problem into a master problem and subproblems, and then iteratively approaches the optimal solution of the original problem. This method has significant advantages in dealing with large-scale optimization problems because it can effectively reduce the size and complexity of the problem, and is easy to implement and applicable to scenarios with a large number of integer and continuous variables. The master problem is called the Master Problem (MP), and the subproblems are called subproblems (SP). Suppose the original problem is a mixed-integer linear programming problem, which can be expressed in general form as follows: (40) in, It is an integer variable. It is a continuous variable. It is the coefficient vector corresponding to the integer variables in the objective function. express transpose, It is the coefficient vector corresponding to the continuous variables in the objective function. express transpose; It is a constraint matrix corresponding to integer variables. It is the constraint matrix corresponding to continuous variables. It is a constraint vector; The core idea of ​​Pendes decomposition is to decompose the original problem into MP and SP; MP contains only integer variables. SP, on the other hand, deals with continuous variables. By iteratively solving MP and SP, cutting planes (Cuts) are gradually generated, thus approximating the optimal solution to the original problem; the objective of MP is to minimize integer variables. The objective function value, while also considering the optimal objective value of SP, is in the form of: (41) in, This represents the optimal objective value of SP; SP is a continuous variable. The linear programming problem is of the form: (42) The goal of SP is based on what MP gives. Solve for the optimal And calculate the objective function value; In summary, the Pendes decomposition algorithm is a classic optimization technique that solves mixed integer programming problems by decomposing the original problem into MP and SP components. The algorithm's implementation involves five key steps: First, starting from an initial integer variable... First, solve for SP; second, for a given... Solve for SP to obtain the optimal objective value. Next, a cutting plane is generated. If SP has no feasible solution, a feasible cut is added to MP. This cut is generated based on the dual information of SP to ensure that the MP has a feasible solution. This is feasible in SP. If SP has a feasible solution, an optimality cut is added to MP. This cut, also based on the dual information of SP, is used to constrain the numerical value of the objective function of MP. Then, MP is updated, the generated cut plane is added to MP, and MP is solved again to obtain a new solution. Finally, the iterative process repeats the above steps until the convergence condition is met, such as the difference between the target values ​​of MP and SP being less than a certain threshold. Through this series of steps, the Pendes decomposition algorithm can effectively approximate the optimal solution of the original problem. 2. Improved Bends decomposition algorithm using interior point method To address the shortcomings of the traditional Pendes decomposition algorithm in terms of convergence speed and the requirement for high-quality cutting planes, the interior point method and the Pendes algorithm are combined. The interior point method, as an algorithm for solving convex optimization problems, explores the optimal solution path within the feasible region to approach the optimal solution. The Pendes IPM algorithm is implemented through the following steps: First, select an initial interior point. The interior point is usually taken as the middle value of the variable to ensure that the initial point is located inside the feasible region; this choice helps the algorithm to start the search from the central region of the feasible region, which may improve the convergence speed. Next, based on the current interior point Solve the dual SP and add constraints generated by the dual SP to MP; these constraints help guide the search direction of MP, making it closer to the optimal solution; then, after adding the new constraints, solve the relaxed MP, i.e., ignore the integer constraints, to obtain the extreme points. The purpose of this step is to find a continuous variable solution that minimizes the objective function value; then, based on the extreme points... Update Inner Points Use the formula: (43) in is an interior-point update coefficient, representing a weighted average between the current solution and the extreme points to achieve stable convergence; k represents the number of iterations. Let the interior point be represented in the k-th iteration. This represents the interior point at the (k+1)th iteration, i.e., the updated interior point; The extreme point obtained by solving the relaxed MP represents the optimal solution obtained by solving the relaxed master problem after adding the new constraints generated by the dual subproblem. Finally, the above process is repeated until the convergence condition is met.

[0025] The unit combination method for the integrated electrothermal energy system based on interior-point Bends decomposition, wherein the Bends IPM algorithm is used to solve the UC-CEHN model proposed above; MP is responsible for solving the main UC problem and verifying the feasibility of the grid network constraints, and SP is responsible for verifying the feasibility of the DHN constraints; a)MP: The main problem part updates the interior points according to the given interior point update coefficients, checks the feasibility in (21)-(24), and minimizes the total operating cost of the generator set commitment and scheduling plan in (1) under the condition that the feasibility cut provided by the constraints in (2)-(20) is satisfied. b) SP: Check the feasibility of network constraints by solving these subproblems (33)-(39); since network constraints at different times are independent, it means that their subproblems can be solved in parallel.

[0026] UC: Full name Unit Commitment UCMP: Full name for Unit Commitment Master Problem UC-CEHN: Full name is Unit Commitment with Combined Electricity and District Heating Networks.

[0027] As can be seen from the above technical solutions, the embodiments of this application have the following advantages: This invention focuses on the unit combination problem of a comprehensive power and thermal energy system considering complete temperature dynamics, conducts in-depth analysis, establishes a mixed-integer linear programming model, covering power decision variables, thermal decision variables, objective functions, and constraints, and clarifies the optimization direction and limiting factors. Regarding the dynamic model of the thermal network, it explores the basic district heating network (DHN) model and its complete analysis method, verifies the correctness of the application of the characteristic line method in solving the pipe temperature conduction equation, and analyzes its advantages over the discretization method; it utilizes an improved Bends decomposition-based approach to solve the unit combination problem of a comprehensive power and thermal energy system considering complete temperature dynamics. The Bends IPM algorithm of this invention solves the problems of slow convergence speed and unstable cutting plane quality that occur when the traditional Bends decomposition algorithm solves the unit combination problem. The key innovation of this algorithm is the integration of the efficient optimization strategy of the interior point method with the Bends decomposition framework, improving the convergence speed and numerical stability of the algorithm through a carefully designed initial interior point selection and dynamic interior point update mechanism. Regarding the initial interior point selection, intermediate values ​​of variables are used to ensure that the algorithm starts its search from the center of the feasible region. This avoids slow convergence or iteration failure due to inappropriate initial point selection. In each iteration, the algorithm generates constraints by solving the dual subproblem, adds them to the main problem, then solves the relaxation main problem to obtain the extreme points, and dynamically updates the interior points based on these extreme points. This dynamic update mechanism achieves stable convergence by using a weighted average of the interior point update coefficients between the current solution and the extreme points, while also improving the algorithm's efficiency. Attached Figure Description

[0028] Figure 1 This is a flowchart of the method steps of the present invention.

[0029] Figure 2 This is a schematic diagram of DHN.

[0030] Figure 3 This is a schematic diagram of the bilateral feature line method.

[0031] Figure 4 This is a flowchart of the solution of UC-CEHN based on the Pendes decomposition method of this invention.

[0032] Figure 5 This is a schematic diagram illustrating the selection of the initial interior point in this invention.

[0033] Figure 6 This is a schematic diagram of the interior point update method of the present invention.

[0034] Figure 7 This is a flowchart of the process for solving UC-CEHN based on the Pendes IPM algorithm of this invention.

[0035] Figure 8 This is a topology diagram of a 118-node power grid model.

[0036] Figure 9 This is a topology diagram of a 118-node heating network model.

[0037] Figure 10(a) shows the electrical-thermal load profile of the 118 nodes.

[0038] Figure 10(b) shows the wind power output data for the 118 nodes.

[0039] Figure 11 This is a topology diagram of a 30-node power grid model.

[0040] Figure 12 This is a topology diagram of a 30-node thermal network model.

[0041] Figure 13(a) shows the electrical-thermal load profile of 30 nodes.

[0042] Figure 13(b) shows the wind power output data for 30 nodes. Detailed Implementation

[0043] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present application, and not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present application.

[0044] like Figure 1 As shown, this invention discloses a method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition, characterized by comprising: Step 101: Construct a mixed-integer linear programming model to clarify the decision variables such as power system unit start-up / shutdown / dispatch, thermal system output, and heating network temperature. With the goal of minimizing the total system operating cost, integrate multiple constraints such as power balance, spinning reserve, ramp rate, and thermal network transmission. In the dynamic modeling of the thermal network, the method of characteristics is used to analyze the partial differential equation of temperature conduction in the pipeline, avoiding discretization errors and improving the accuracy and computational efficiency of the dynamic temperature description. Step 102: Use the improved Pendes decomposition algorithm combined with the interior point method to solve the mixed integer linear programming model and obtain the final result.

[0045] The invention addresses the unit combination problem that considers complete temperature dynamics by constructing a mixed-integer linear programming (MILP) model.

[0046] I. Unit Combination Model of Integrated Power and Heat System 1. Decision variables The unit configuration of a power system primarily focuses on minimizing operating costs while meeting electricity demand. This involves optimizing the scheduling, output, and coordination of generator units with the thermal system. In combined heat and power (CHP) units, this problem is more complex because it requires simultaneously balancing the supply and demand of electricity and heat, as well as their interrelationships.

[0047] Power decision variables include the combined state of generator units. and scheduling ;in, Indicates the unit In time The running status, and They are respectively in Startup and shutdown status, It is a generator set In time Power output, It is a wind turbine In time Power output, and They are the generator sets In time The upward and downward rotation of reserve capacity. Thermal decision variables encompass the thermal output of the CHP unit and the supply / return water temperatures of the District Heating Network (DHN). Thermal output represents the heat power provided by the unit while meeting electricity demand. The supply / return water temperatures of the DHN directly affect the efficiency of heat transfer and distribution.

[0048] 2. Objective function The objective function aims to minimize the total operating cost of the system, encompassing the operating costs of the power and heating systems, as well as additional costs due to the uncertainties of wind power. This necessitates a comprehensive consideration of the operating characteristics of both systems and the uncertainties of renewable energy sources like wind power in the model.

[0049] Minimize total operating costs: (1) In the cost structure, the first item represents the operating cost of non-CHP units, the second item represents the cost of wind curtailment penalty, and the third item represents the operating cost of CHP units.

[0050] The operating cost of non-CHP units is defined as follows: (2) in, This represents the piecewise linear power generation cost, which approximately describes the nonlinear relationship between fuel cost and unit output; This represents the startup cost, when the unit is in time. A fee is required upon startup, the amount of which depends on the unit type and startup mode. This represents the shutdown cost, when the unit is in time. When shutting down, additional costs need to be considered, including equipment cooling and maintenance. This represents the cost of no-load operation.

[0051] The formula for calculating the cost of wind curtailment penalty is: (3) in, Indicates wind turbine In time The actual output, Indicates wind turbine In time The planned output, Indicates the wind curtailment penalty coefficient; Let T represent the set of wind turbine units, and let T represent the time set.

[0052] The operating costs of CHP units are similar to those of non-CHP units, including the cost of power generation. Defined as: (4) Among them, piecewise linear power generation cost This describes the fuel consumption of the CHP unit at different output levels. Start-up cost item. This reflects the additional costs incurred when the CHP unit starts up at time t. (Cost item closed) This covers the costs incurred when the unit is shut down at time t; no-load cost item. This represents the energy consumption and maintenance costs of the CHP unit to maintain basic operation under no-load conditions.

[0053] 3. Constraints: The model proposed in this paper is constrained by the operational limitations of the power and heat networks.

[0054] (1) Power balance constraint: The total power generation and total load must be kept in balance during each time period. (5) This constraint ensures that, within each time period, the total power generation of all generating units, including conventional units, CHP units, and wind turbines, is balanced with the total power load of the system. Among these, This represents the actual power output of a conventional or CHP unit at time t. This is the operating status variable of the generator set; it is set to 1 when the generator set is running and 0 otherwise. This refers to the rated power of the unit. For wind turbine units, their power generation is determined by... Indicates. The right side. This represents the total power load demand of the system at time t. This constraint is a fundamental requirement for the operation of the power system, ensuring the real-time balance between power supply and demand.

[0055] (2) Spinning Reserve Constraint: Non-CHP units provide the required spinning reserve capacity. (6) (7) (8) (9) (10) Spinning reserve constraints ensure that the system has sufficient reserve capacity to maintain stable operation in the face of load fluctuations or sudden events. Among these, Indicates the unit In time The upward rotation of the spare capacity, Indicates the unit In time The downward rotation of the spare capacity, Indicates the unit The maximum starting ramp rate, Indicates the unit The maximum rate of ... It is a generator set Minimum boot time, This indicates the total backup demand for the entire system rotating upwards. This indicates the total backup demand for the entire system rotating downwards. This represents a set of thermal power units; these constraints, by limiting the relationship between the generating output and reserve capacity of the units, ensure that the system can meet the requirements of spinning reserve under different operating conditions. Specifically, the first three inequalities limit the sum of the unit's maximum power output and reserve capacity based on the unit's start-up and shutdown status and operating mode, preventing the unit from exceeding its operating limits; the last two inequalities, at the system level, specify the minimum total amount of upward and downward rotational reserves in each time period, ensuring that the system has the ability to cope with load changes.

[0056] (3) Ramp-up rate constraint: Within a single time period, the incremental power generation of a thermal power unit is limited by its ramp-up capability. (11) The ramping constraint limits the rate of change in power generation output of a thermal power unit within adjacent time periods, i.e., the unit's ramping capability. Indicates the unit The upper limit of the upward climbing rate, Indicates the unit The upper limit of the downward climbing rate, This represents the set of CHP units; this constraint reflects the limitation of the physical characteristics of thermal power units on the speed of output adjustment, avoiding equipment damage or operational instability caused by rapid output adjustment of units, and also helps the economical scheduling of the system, preventing increased operating costs due to frequent adjustments of unit output. (4) Wind power output constraints: The power generation of wind farms is limited by the available wind power. (12) This constraint specifies the power output range of the wind turbine in each time period; the lower limit is 0, indicating that the wind turbine can stop generating electricity; the upper limit... This is determined by the wind conditions at the time, reflecting the renewable energy characteristics and uncertainties of wind power; this constraint ensures that the power generation of wind turbines does not exceed their maximum power generation capacity within a specific period, which is consistent with the actual operation of wind power. (5) Unit status constraints (13) Generator set state constraints describe the dynamic relationships between the unit's operating state, startup state, and shutdown state. Among them, Indicates the unit In time The running status, Indicates the unit In time The startup state, Indicates the unit In time The equation shows that the change in the unit's operating state during the current period is equal to the difference between its start-up and shutdown states, thus ensuring the continuity and logical consistency of the unit's state and providing a basis for accurately modeling the unit's start-up and shutdown process. (6) Minimum start-stop time constraint (14) (15) The minimum start-up and shutdown time constraint ensures that once the unit starts up or shuts down, it will remain in that state for at least a certain period of time. The first inequality stipulates that, in time... Before During the time period, the unit The number of times a unit can be started cannot exceed its operating state; that is, if the unit is currently in an operating state, then the number of times it has been started in the past... The first inequality ensures that the unit is started at least once within a certain time period; the second inequality ensures that if the unit is currently shut down, then it was started at least once in the past. The unit can be shut down at most once within a given time period; these constraints reflect the physical requirements for hot and cold starts of the unit, avoiding damage to the equipment caused by frequent start-stops, while also helping to reduce operating costs and improve system reliability. (7) Variable boundary constraints (16) (17) The upper and lower bound constraints of variables clearly define the range of values ​​for each decision variable; for the operating status of the unit... It can only take the value 0 or 1, indicating whether the unit is running at time t; startup and shutdown status variables. and The value of is between 0 and 1, but in actual operation it is usually 0 or 1 to ensure the integer programming characteristics of the model; for the rotated spare variable and The values ​​are limited to 0 to the unit's maximum upward and downward ramp rates, ensuring the rationality and feasibility of the reserve capacity.

[0057] (8) Network constraints: The power flow of transmission lines should not exceed the transmission capacity. Equation (18) represents the nodal power equation, Equation (19) represents the branch power flow equation, and Equation (20) defines the phase angle of the reference bus. Wherein, It is with nodes The set of all directly connected nodes It is a node The upper connected conventional generator set, It is a node The collection of combined heat and power units connected to the top, It is a node The upper-connected wind turbine assembly, For nodes The voltage phase angle at time t, For nodes The voltage phase angle at time t, This represents the voltage phase angle of the reference node in the network at time t. Let (i,j) be the reactance of the line. For the line ( i , j Maximum transmission capacity, This is the set of all routes.

[0058] (18) (19) (20) Network constraints ensure that the power system meets the physical limitations and safety requirements of power transmission during operation. Node power equations link the node voltage phase angle difference with the power transmitted through lines, reflecting the fundamental laws of power flow in the power system. Branch power flow equations limit the power transmission capacity of each transmission line, preventing line overload and ensuring the safe and stable operation of the system. The setting of the reference bus phase angle provides a benchmark for the voltage phase angle of the entire system, facilitating calculation and analysis. These constraints collectively constitute the operational framework of the power network, ensuring the efficiency and reliability of power transmission.

[0059] The constraints described above indicate that the unit combination model of a power system must consider the operating characteristics and cost factors of the units, while also meeting the physical limitations and security requirements of the power network. These constraints collectively constitute a complex optimization problem, providing a solid theoretical foundation and model framework for the efficient operation and cost optimization of integrated power and thermal energy systems.

[0060] II. Dynamic Model of Thermal Networks In integrated power-thermal systems, the dynamic characteristics of the thermal network have a significant impact on the system's stable operation and optimal regulation. This section will explore in depth the dynamic modeling methods for thermal networks, including the DHN model and its complete analytical approach.

[0061] DHNs are two-layer networks integrating supply and return water networks. For example... Figure 2 As shown, heat energy is generated from a heat source, transported through water supply pipelines, distributed at a heat exchange station, and then flows to heat users. After heat exchange at the user end, the water flows back to the heat source through the return water network. Therefore, as... Figure 2 As shown, a typical DHN model includes both hydraulic and thermal components. In this paper, the mass regulation mode of DHNs under constant mass flow rate is investigated, using temperature regulation to meet the heat load.

[0062] For simplicity, the temperature mentioned in this paper is defined as relative temperature, that is, the absolute temperature difference between the water flow and the ambient temperature. The mass model formula for DHN under mass regulation is as follows: (twenty one) (twenty two) (twenty three) (twenty four) in, This represents the longitudinal distance of a point in the pipeline from the pipeline inlet, i.e., the spatial coordinates along the pipeline. Indicates time, Indicates the first j The temperature of the pipe, Indicates the first j The water flow velocity in the pipe This indicates the specific heat capacity of water. Indicates the first j Mass flow rate of the pipeline Indicates the first j The heat loss coefficient of the pipeline, Represents a node thermal power, Represents a node Mass flow rate Represents a node Water supply pipe temperature Represents a node Return water pipe temperature Describes the set of pipes flowing into node j. Describes the set of pipes that flow out of node j. This represents the temperature at the outlet of pipe b, where L is the length of the pipe. This represents the mass flow rate of pipe b. This indicates the temperature at the inlet of pipe b. This represents the mixing temperature at node k. Represents a set of pipes. Represents a set of nodes. Let represent the set of pipelines that start at node k and end at node i.

[0063] The temperature conduction equation represents the temperature distribution along a pipe, and this equation is constructed based on the following assumptions: 1) The heat loss of the pipe cross-section is considered constant; 2) Ignore heat conduction between adjacent water flows.

[0064] This equation models the nodes in the DHN as heat exchangers, thus allowing a linear relationship between mass flow rate and temperature. The nodes obey the law of energy conservation, and the mixing temperature at a node is a linear weighted sum of the outlet temperatures of the pipes injected into that node. Finally, the node mixing temperature at the pipes leaving the node is the same as the inlet temperature.

[0065] The heat conduction equation in a pipe is a typical first-order linear and quasi-linear partial differential equation (PDE), commonly used to describe vibration and wave phenomena and their dynamic processes. Its representative form is: (25) Where a, b, and c are the coefficients of the first-order linear hyperbolic partial differential equation. Specifically, a corresponds to the rate of change of temperature in space, b corresponds to the rate of change of temperature in time, and c corresponds to the temperature distribution itself. Let x be the state variables of a first-order linear hyperbolic partial differential equation, where x is the spatial coordinate and t is the time variable. This represents the temperature distribution function, i.e., the temperature at pipe location x and time t. These are the initial conditions for hyperbolic partial differential equations. This represents the initial temperature distribution at location x.

[0066] The solution plane of this equation consists of two-dimensional vectors in space and time, and is also known as the Cauchy problem.

[0067] The core idea of ​​the Characteristic Line Method (CLM) is to find characteristic lines. Along these characteristic lines, the original PDE can be transformed into an ordinary differential equation (ODE); from this perspective, the solution plane can be considered as being formed by the initial condition curves on... The surface formed by translating along the feature line; such as Figure 3 As shown.

[0068] The steps of traditional CLM are as follows: The differential of u with respect to t is expressed as: (26) Substituting the above expression into the pipe temperature conduction equation, we observe that when... When this happens, equation (25) can be converted to ODE; therefore, The solution is defined as the characteristic line of equation (25), expressed as: (27) Where c is a constant determined by the initial conditions.

[0069] Replacing the partial derivative with the total differential, the reconstructed ODE is expressed as: (28) The analytical solution to the above equation can be obtained from the coefficients. Different expressions and initial conditions get.

[0070] To solve the Cauchy problem using CLM, a set of initial conditions is needed as the boundary of the solution plane. In other words, the PDEs that can be solved by traditional CLM must be unilaterally bounded along the direction of the feature line. Consider a practical scenario where... and Both are bounded, and the bilaterally bounded hyperbolic PDE is extended as follows: (29) in, Describes the boundary conditions for hyperbolic partial differential equations. This indicates the temperature setpoint at the heat source end of the pipeline.

[0071] Unlike traditional CLM, its solution is determined by both initial and boundary conditions, expressed as: (30) Determined by initial conditions, this indicates the temperature distribution caused by the initial heat storage in the pipeline; Determined by boundary conditions, it represents the temperature propagation distribution caused by real-time heating from a heat source; such as Figure 3 As shown, the solution plane can be considered as the initial condition plane along... The result of the intersection of the surface formed by the movement and the surface formed by the movement of the boundary conditions; (31) (32) in, Right now , This represents the temperature distribution caused by the initial heat storage in the pipe at the initial time t=0. Right now , This represents the temperature propagation distribution at x=0 at the beginning of the pipeline, where heat is supplied in real time by a heat source.

[0072] Based on the above, the fully analytical method based on CLM has the following advantages over discretization methods: (1) Higher accuracy Due to discretization methods, approximation errors caused by Taylor or Fourier expansions are unavoidable. For the heat conduction equation, the analytical solution is continuous, while discretization methods typically produce discrete solutions. Therefore, analytical solutions are usually more accurate than discretization methods. Methods based on the Expectation-Maximization Algorithm (EM) have errors proportional to the spatial step size and are more sensitive to discrete boundary conditions, which can lead to numerical oscillations. For methods based on the Nelder-Mead Method (NM), the error is proportional to the time step size and is more likely to cause numerical diffusion. For Fourier methods, the approximation error is proportional to the number of sinusoidal components and the influence of initial conditions, which is also a drawback of Laplace-based methods.

[0073] For the Partitioning Around Medoids Algorithm (PAM), the thermodynamic PDE is solved analytically without discretization. Approximation errors and numerical oscillations can be completely avoided.

[0074] (2) Higher computational efficiency From a mathematical perspective, the computational burden of EM-based methods is proportional to the number of differential steps, while the computational burden of NM-based methods is proportional to the length of the historical temperature series. Unlike numerical methods, the pipe outlet temperature can be directly obtained analytically through a function calculation. Calculations regarding pipe cross-sections and temperature series are no longer required. Therefore, the fully analytical method is mathematically the most efficient.

[0075] III. Optimization Model of Integrated Electric and Thermal Energy System Without considering uncertainties, this section employs a deterministic optimization method to construct a model of the Unit Commitment with Combined Electricity and District Heating Networks (UC-CEHN), as detailed below: , ,

[0076] This model optimizes system operation using the MILP method, comprehensively considering generator start-up and shutdown plans, power output scheduling, and coordination of heat supply. The constraints in the model cover power system power balance, spinning reserve, ramp-up limits, wind turbine output range, generator state transitions, and minimum start-up and shutdown times. It also incorporates relevant constraints from the heat network to ensure the stable operation of the integrated electricity-heat energy system.

[0077] To address the complex problem of optimizing integrated electricity-heat energy systems, this paper introduces the Pendesic decomposition technique to construct a distributed decision-making framework between the power grid operator (PGO) and the district heating operator (DHO). Within this framework, the PGO is primarily responsible for solving the main UC problem and checking the feasibility of power grid constraints, while the DHO focuses on checking the feasibility of DHN constraints. This decomposition-coordination model not only improves computational efficiency but also protects the information privacy of both networks to some extent.

[0078] 1. Main UC Problem: Provide a combination and scheduling plan for power generation equipment to minimize the total operating cost in (1) while satisfying the constraints in (2)-(20), (21)-(24) and the feasible cutting plane provided by the feasibility check. As a MILP, the main UC problem can be solved by the Lagrange relaxation method or the MILP solver.

[0079] 2. Network security checks: Solving these sub-problems is to check the plans of a given crew. Feasibility of applying network constraints. Since the network constraints at different times are independent, the subproblems representing these constraints can be solved in parallel. Each subproblem is expressed as follows: (33) (34) (35) in, The optimal value of the network security inspection subproblem represents the total number of violations of power flow constraints under a given unit schedule in time period t. This represents the planned active power output vector of the generating units during time period t; This represents the planned start-up and shutdown status vector of the generating units during time period t; This is the first relaxation variable, representing a virtual positive power injection used to compensate for the power deficit at this node. The second slack variable represents the virtual negative power absorption, used to absorb the power excess at the node. The two slack variables are used to ensure that the node power balance equation is always solvable.

[0080] if For a given unit plan Network constraints are feasible. If Then, add the following Pendesic cut plane to UCMP: (36) in, This represents the vector of unit output variables for time period t. This represents the vector of unit start-up and shutdown status variables for time period t.

[0081] 3. DHN Feasibility Check Sub-problem: The purpose of solving this sub-problem is to check the plan for a given unit. The feasibility of the following DHN operating constraints; its calculation formula is as follows: (37) (38) The constraints are (2)-(20) and (21)-(24).

[0082] in, The optimal value of the DHN feasibility check subproblem is represented by , and the total number of violations of thermal network constraints is represented by , under a given unit plan. Represents the thermal output planning vector. This represents the planned start-up and shutdown state vector of the thermal power unit. It is the third relaxation variable, representing a virtual positive heat power injection used to compensate for the heat power deficit at this node. The fourth relaxation variable represents the virtual negative heat power absorption, used to absorb the excess heat power at this node. The two relaxation variables are used to ensure that the heat power balance equation of thermal node j is always solvable at time t. Represents the set of all nodes in a heat network. This represents the set of all heat sources supplying heat to node j at time t. The minimum thermal output coefficient of heat source g, This represents the mass flow rate of node j at time t. This indicates the specific heat capacity of water. This represents the water supply temperature setpoint for node j at time t. This represents the setpoint for the return water temperature at time t for node j.

[0083] if For a given unit plan DHN constraints are feasible. If Then the following Pendesic cut plane will be fed back into the main problem: (39) in, Represents the vector of thermal output variables. This represents the start-up and shutdown state variable vector of the thermal power unit.

[0084] Through the detailed explanation of the above model and sub-problems, we can see that the optimization model of the integrated electric-thermal energy system not only needs to consider the operating characteristics and cost factors of the power system, but also needs to meet the physical constraints and operational requirements of the thermal network. The above model constitutes a complex optimization problem, providing a solid theoretical foundation and model framework for achieving the safe and economical operation of the integrated electric-thermal energy system.

[0085] IV. Computer Group Assembling Problem Based on Improved Pendesic Decomposition 1. Traditional Pendes decomposition algorithm The Pendes decomposition algorithm breaks down the original problem into a master problem (MP) and subproblems (SP), and then iteratively approximates the optimal solution to the original problem. This method shows significant advantages when dealing with large-scale optimization problems because it can effectively reduce the size and complexity of the problem, while being easy to implement and applicable to scenarios with a large number of integer and continuous variables.

[0086] Suppose the original problem is a mixed-integer linear programming problem, its general form can be expressed as: (40) in, It is an integer variable. It is a continuous variable. It is the coefficient vector corresponding to the integer variables in the objective function. express transpose, It is the coefficient vector corresponding to the continuous variables in the objective function. express transpose; It is a constraint matrix corresponding to integer variables. It is the constraint matrix corresponding to continuous variables. It is a constraint vector; The core idea of ​​Pendes decomposition is to break down the original problem into MP and SP. MP contains only integer variables. SP, on the other hand, deals with continuous variables. By iteratively solving MP and SP, cutting planes (Cuts) are gradually generated, thus approximating the optimal solution to the original problem. The goal of MP is to minimize integer variables. The objective function value, while also considering the optimal objective value of SP, is in the form of: (41) in, This represents the optimal objective value of SP. SP is a continuous variable. The linear programming problem is of the form: (42) The goal of SP is based on what MP gives. Solve for the optimal And calculate the objective function value.

[0087] In summary, the Pendes decomposition algorithm is a classic optimization technique that solves mixed integer programming problems by decomposing the original problem into MP and SP components. The algorithm's implementation involves five key steps: First, starting from an initial integer variable... First, solve for SP; second, for a given... Solve for SP to obtain the optimal objective value. Next, a cutting plane is generated. If the SP has no feasible solution, a feasibility cut is added to the MP. This cut is generated based on the dual information of the SP and is used to ensure that the MP has a feasible solution. This is feasible in SP. If SP has a feasible solution, an optimality cut is added to MP. This cut is also based on the dual information of SP and is used to constrain the numerical value of the objective function of MP. Then, MP is updated, the generated cut plane is added to MP, and MP is solved again to obtain a new solution. Finally, the iterative process repeats the above steps until a convergence condition is met, such as the difference between the target values ​​of MP and SP being less than a certain threshold. Through this series of steps, the Pendes decomposition algorithm can effectively approximate the optimal solution to the original problem.

[0088] UC: Full name Unit Commitment UCMP: Full name for Unit Commitment Master Problem UC-CEHN: Full name is Unit Commitment with Combined Electricity and District Heating Networks.

[0089] The proposed UC-CEHN model can be solved using the traditional Pendesic decomposition algorithm. Figure 4 In the process shown, MP is responsible for solving the main UC problem and verifying the feasibility of the power grid network constraints, while SP is responsible for verifying the feasibility of the DHN constraints.

[0090] a)MP: The main problem provides the commitment and scheduling plan of the generator set to minimize the total operating cost in (1) while satisfying the constraints in (2)-(20), (21)-(24) and the feasibility check subproblems.

[0091] b) SP: Check the feasibility of network constraints by solving subproblems (33)-(39). Since network constraints at different times are independent, their subproblems can be solved in parallel.

[0092] The Pendes decomposition method transforms a complex optimization problem into two relatively independent sub-modules. The main problem handles the optimization process involving integer variables, while the sub-problems aim to solve for continuous variables. This decomposition strategy improves the efficiency of the algorithm in handling large-scale problems, and its advantages are particularly pronounced in cases containing a large number of integer and continuous variables. From an implementation perspective, this method offers high flexibility; the main problem and sub-problems can be handled by different solvers. This modular design greatly simplifies the algorithm's implementation process. Empirical studies demonstrate that this method is suitable for optimization problems with large numbers of integer and continuous variables, effectively improving solution efficiency and optimizing computational resource consumption. However, this method has some limitations. When there is a strong coupling between the main problem and the subproblems, the convergence speed of the algorithm may decrease. This is because the solution to the main problem needs to be iteratively adjusted based on the feedback from the subproblems, and the tight coupling makes this adjustment process more complex. The quality of the cutting plane plays a decisive role in the algorithm's performance. If the constraints generated by the solutions to the subproblems cannot effectively limit the solution space of the main problem, the algorithm will require more iterations to converge. Improving the quality of the cutting plane and accelerating the convergence process have become the core research topics for optimizing the Bends decomposition algorithm. It is evident that although the Pendes decomposition algorithm has theoretical advantages and can address the shortcomings of currently widely used unit combination solution algorithms, it has some limitations in practical applications. For example, the algorithm's convergence speed may be slow when the coupling between MP and SP is strong, and the quality of the cutting plane directly affects the algorithm's efficiency. If the cutting plane cannot effectively constrain the solution space of MP, it may lead to slow convergence. To overcome these limitations, an improved Pendes decomposition algorithm based on the interior point method is proposed.

[0093] 2. Improved Bends decomposition algorithm using interior point method To address the shortcomings of traditional Pendes decomposition algorithms in terms of convergence speed and the requirement for high-quality cutting planes, this study aims to combine the interior point method with the Pendes algorithm. The interior point method, as an algorithm for solving convex optimization problems, explores the optimal solution path within the feasible region to approach the optimal solution. The advantage of this method is that it can handle optimization problems with complex constraints and generally has good numerical stability and convergence speed.

[0094] The Pendes IPM algorithm is implemented through the following steps: First, select an initial interior point. ,like Figure 5 As shown, the interior point is usually taken as the middle value of the variable to ensure that the initial point is located inside the feasible region; this choice helps the algorithm to start the search from the central region of the feasible region, which may improve the convergence speed.

[0095] Next, based on the current interior point Solve the dual SP and add constraints generated by the dual SP to MP. These constraints help guide the search direction of MP, bringing it closer to the optimal solution. Then, after adding the new constraints, solve the relaxed MP (i.e., ignore integer constraints) to obtain the extreme points. The goal of this step is to find a continuous variable solution that minimizes the objective function value. Then, based on the extreme points... Update Inner Points Use the formula: (43) in is an interior-point update coefficient, representing a weighted average between the current solution and the extreme points to achieve stable convergence; k represents the number of iterations. Let the interior point be represented in the k-th iteration. This represents the interior point at the (k+1)th iteration, i.e., the updated interior point; These are the extreme points obtained by solving the relaxed MP problem, representing the optimal solution obtained by solving the relaxed master problem after adding new constraints generated by the dual subproblem. Finally, as... Figure 6 As shown, repeat the above process until the convergence condition is met, for example, the difference between the target values ​​of MP and SP is less than a preset threshold, or the maximum number of iterations is reached.

[0096] In the interior-point method, the interior-point update coefficients The choice of interior point update coefficients can significantly impact algorithm performance; therefore, it is crucial for the algorithm's convergence and stability. An excessively large step size may cause the algorithm to skip the optimal solution during the search process, while an excessively small step size may result in slow convergence. Therefore, selecting suitable interior point update coefficients is key to achieving efficient algorithm operation.

[0097] The Pendes IPM algorithm can be used to solve the UC-CEHN model proposed above. Figure 7 In the process shown, MP is responsible for solving the Unit Commitment Master Problem (UCMP) and verifying the feasibility of the power grid network constraints, while SP is responsible for verifying the feasibility of the DHN constraints.

[0098] a)MP: The main problem part updates the interior points according to the given interior point update coefficients, checks the feasibility in (21)-(24), and minimizes the total operating cost of the generator set commitment and scheduling plan in (1) under the condition that the feasibility cut provided by the constraints in (2)-(20) is satisfied.

[0099] b) SP: Check the feasibility of network constraints by solving subproblems (33)-(39). Since network constraints at different times are independent, their subproblems can be solved in parallel.

[0100] The Bends IPM algorithm addresses the slow convergence speed and unstable cutting plane quality issues inherent in traditional Bends decomposition algorithms for solving unit combination problems. Its key innovation lies in integrating the efficient optimization strategy of the interior-point method with the Bends decomposition framework. Through carefully designed initial interior-point selection and a dynamic interior-point update mechanism, the algorithm improves both convergence speed and numerical stability. Regarding initial interior-point selection, intermediate values ​​of variables are used to ensure the algorithm starts its search from the center of the feasible region, avoiding slow convergence or iteration failures due to inappropriate initial point selection. In each iteration, the algorithm generates constraints by solving the dual subproblem, adds them to the main problem, then solves the relaxation main problem to obtain extreme points, and dynamically updates the interior points based on these extreme points. This dynamic update mechanism achieves stable convergence by using a weighted average of the interior-point update coefficients between the current solution and the extreme points, while simultaneously improving the algorithm's efficiency.

[0101] 3.3 Application Examples of the Invention To verify the performance of this invention, simulations were performed using the IEEE standard system. All case studies were conducted on a personal computer equipped with an AMD Ryzen 9 9950X processor (4.30GHz, 16 cores) and 64GB of RAM. Program development was performed using Matlab R2023b, and all mixed-integer linear programming (MILP) and linear programming (LP) problems were solved using Gurobi 11.0.3. The relative error of the MILP solver was set to within 0.01%.

[0102] The topology diagram of the power grid and heating network of the 118-node test system is as follows: Figure 8 , Figure 9 As shown. The heat load of the district heating network (DHN) is provided by two extraction-condensing CHP units of the thermal power plant. Wind turbines (W1, W2, and W3) are connected to grids Bs77, Bs78, and Bs79, respectively, while CHP units (CHP1 and CHP2) are connected to grids Bs1 and Bs9, respectively. System configuration details are shown in Table 1, and detailed data can be found in Figures 10(a) and 10(b). Both uplink and downlink reserve requirements are 50 MW. The penalty price for wind power reduction is set at the maximum incremental cost of conventional unit generation.

[0103] In Figure 10(a), the trends of electrical load (blue curve) and heat load (red curve) over time show obvious diurnal fluctuations. Electrical load is higher during the daytime, while heat load is higher at night and in the early morning, which is consistent with users' electricity and heat consumption habits. Figure 10(b) shows the changes in wind power output of the three wind turbines over time. It can be seen that the changes in wind power output are quite significant in different time periods, and the output of each turbine has a certain degree of fluctuation.

[0104] Table 1. 118-Node Model Test Model Information Table

[0105] The topology diagram of the power grid and heating network of the 30-node test system is as follows: Figure 11 , Figure 12 As shown. The heat load of the district heating network (DHN) is provided by a single extraction-condensing CHP unit from the thermal power plant. Wind turbines (W1) are connected to grid Bs4, while the CHP unit (CHP) is connected to grid Bs1. System configuration details are shown in Table 2, with further data in Figures 13(a) and 13(b). Both uplink and downlink reserve requirements are 50 MW. The penalty price for wind power reduction is set at the maximum incremental cost of conventional unit generation.

[0106] Table 2. 30-Node Model Test Model Information Table

[0107] 1. 118-node model A centralized approach was used to solve the problem, completing the solution in 0.57 seconds with a total operating cost of $1,120,511.15. The centralized approach employed in this case is a strategy that directly integrates all decision variables into a large-scale linear or mixed-integer linear programming model, using an efficient solver like Gurobi to obtain the solution. This method fully leverages the capabilities of modern solvers to quickly find the global optimum.

[0108] The Pendescendant decomposition algorithm was used to solve the problem. This method successfully converged after 706 iterations, taking a total time of 2286.13 seconds, with a final total running cost of $1,120,561.39. The Pendescendant decomposition algorithm is a classic optimization algorithm that decomposes the original problem into MP and SP components, gradually approaching the optimal solution by alternating between them. While this method performs well in some cases, its number of iterations and computation time can be relatively high.

[0109] The Bendes IPM algorithm, an improved version based on the interior-point method, was used to solve the problem. The interior-point convergence coefficient was tentatively set to 0.3. This method successfully converged after 457 iterations, with a total time of 699.57 seconds, and the final total running cost was $1,120,561.39. The Bendes IPM method introduces interior-point update coefficients within the feasible region, avoiding iteration points touching the feasible region boundary and reducing the number and amplitude of numerical oscillations during iteration. This method improves solution efficiency and ensures solution accuracy. Table 3 below records the solutions of the three algorithms. Comparing the solution results of the three methods, the maximum difference in the total running cost is only 0.045‰, indicating that all three methods can obtain accurate optimal solutions.

[0110] Table 3. Solutions, solution times, and relative errors for each algorithm with 118 nodes.

[0111] Compared to the Bendes decomposition algorithm, the Bendes IPM algorithm shows significant advantages in solving 118-node cases. Specifically, the Bendes IPM algorithm reduces the number of iterations by 35.32% and the computation time by 69.42%.

[0112] 2. 30-node model A centralized algorithm was used to solve the problem. This method completed the solution in 0.09 seconds, with a final total operating cost of $154,611.39.

[0113] The Pendesic decomposition algorithm was used to solve the problem. This method converged successfully after 441 iterations, with a total execution time of 77.42 seconds. The final total running cost was $154,611.39.

[0114] The Bends IPM algorithm, an improved version based on the interior-point method, was used for the solution, with the interior-point convergence coefficient tentatively set to 0.7. The method successfully converged after 299 iterations, with a total execution time of 44.82 seconds. The final total running cost was $154,611.39.

[0115] Table 4 below records the solutions of the three algorithms. By comparing the solution results of the three methods, we find that the total running cost of the three methods is the same. This result shows that all three methods can find the accurate optimal solution.

[0116] Table 4. Solutions, solution times, and relative errors for each of the 30-node algorithms.

[0117] Compared to the Pendes decomposition algorithm, the Pendes IPM algorithm shows advantages in solving 30-node cases. Specifically, the Pendes IPM algorithm reduces the number of iterations by 27.32% and the computation time by 42.16%. These results indicate that the Pendes IPM algorithm can save computer computing resources, obtain the optimal solution in a shorter time, and maintain the accuracy of the solution.

[0118] Therefore, it can be seen that the Bendes IPM algorithm demonstrates higher efficiency and accuracy in solving specific optimization problems, and is a recommended optimization algorithm. Future research can explore the application potential of the Bendes IPM algorithm in other types of optimization problems.

Claims

1. A method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition, characterized in that, include: Step 101: Construct a mixed-integer linear programming model to clarify the decision variables such as power system unit start-up / shutdown / dispatch, thermal system output, and heating network temperature. With the goal of minimizing the total system operating cost, integrate multiple constraints such as power balance, spinning reserve, ramp rate, and thermal network transmission. In the dynamic modeling of the thermal network, the method of characteristics is used to analyze the partial differential equation of temperature conduction in the pipeline, avoiding discretization errors and improving the accuracy and computational efficiency of the dynamic temperature description. Step 102: Use the improved Pendes decomposition algorithm combined with the interior point method to solve the mixed integer linear programming model and obtain the final result.

2. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 1, characterized in that, The power decision variables of the mixed-integer linear programming model include the combined state of the generator units. and scheduling ;in, Indicates the unit In time The running status, and They are respectively in Startup and shutdown status, It is a generator set In time Power output, It is a wind turbine In time Power output, and They are the generator sets In time Upward and downward rotational reserve capacity; thermal decision variables encompass the thermal output of the combined cogeneration unit and the supply / return water temperature of the district heating network; thermal output represents the heat power provided by the unit while meeting electricity demand; the supply / return water temperature of the district heating network directly affects the efficiency of heat transmission and distribution; Unit Commitment is abbreviated as UC.

3. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 1, characterized in that, The objective function aims to minimize the total operating cost of the system, encompassing the operating costs of the power system and the heating system, as well as the additional costs due to the uncertainty of wind power. Minimize total operating costs: （1） In the cost structure, the first item represents the operating cost of non-CHP units, the second item represents the cost of wind curtailment penalty, and the third item represents the operating cost of CHP units. The English for combined heat and power unit is: Combined Heat and Power, CHP; the English for district heating network is: District Heating Network, DHN. The operating cost of a non-combined thermal power unit is defined as: （2） in, This represents the piecewise linear power generation cost, which approximately describes the nonlinear relationship between fuel cost and unit output; This represents the startup cost, when the unit is in time. A fee is required upon startup, the amount of which depends on the unit type and startup mode. This represents the shutdown cost, when the unit is in time. When shutting down, additional costs need to be considered, including equipment cooling and maintenance. Indicates no-load cost; The formula for calculating the cost of wind curtailment penalty is: （3） in, Indicates wind turbine In time The actual output, Indicates wind turbine In time The planned output, Indicates the wind curtailment penalty coefficient; Let T represent the set of wind turbine units, and T represent the time set. The operating costs of CHP units are similar to those of non-CHP units, including the cost of power generation. Defined as: （4） Among them, piecewise linear power generation cost Describes the fuel consumption of the CHP unit at different output levels; startup cost item. This reflects the additional costs incurred when the CHP unit starts up at time t; the cost item is closed. This covers the costs incurred when the unit is shut down at time t; no-load cost item. This represents the energy consumption and maintenance costs of the CHP unit to maintain basic operation under no-load conditions.

4. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 1, characterized in that, The mixed-integer linear programming model is subject to constraints imposed by the operation of the power and heat networks. The constraints are as follows: Power balance constraint: Total power generation and total load must remain balanced within each time period. （5） This constraint ensures that, in each time period, the total power generation of all generating units, including conventional units, CHP units, and wind turbines, is balanced with the total power load of the system. in, This represents the actual power output of a conventional or CHP unit at time t. This is the operating status variable of the generator set; it is set to 1 when the generator set is running and 0 otherwise. The rated power of the unit; for wind turbines, the power generation is determined by... Indicates; the right side This is the total power load demand of the system at time t; this constraint is a basic requirement for the operation of the power system, ensuring the real-time balance of power supply and demand. Spinning reserve constraint: Non-CHP units provide the required spinning reserve capacity. （6） (7) (8) (9) (10) Spinning reserve constraints ensure that the system has sufficient reserve capacity to maintain stable operation in the face of load fluctuations or sudden events; among which, Indicates the unit In time The upward rotation of the spare capacity, Indicates the unit In time The downward rotation of the spare capacity, Indicates the unit The maximum starting ramp rate, Indicates the unit The maximum rate of ... It is a generator set Minimum boot time, This indicates the total backup demand for the entire system rotating upwards. This indicates the total backup demand for the entire system rotating downwards. This represents a set of thermal power units; these constraints, by limiting the relationship between the generating output and reserve capacity of the units, ensure that the system can meet the requirements of spinning reserve under different operating conditions. Specifically, the first three inequalities limit the sum of the unit's maximum power output and reserve capacity based on the unit's start-up and shutdown status and operating mode, preventing the unit from exceeding its operating limits; the last two inequalities, on the other hand, specify the minimum total amount of upward and downward rotational reserves in each time period at the system level, ensuring that the system has the ability to cope with load changes. Ramp-up constraint: Within a single time period, the incremental power generation of a thermal power unit is limited by its ramp-up capability. (11) The ramping constraint limits the rate of change in power generation output of thermal power units within adjacent time periods, i.e., the ramping capability of the units; among which, Indicates the unit The upper limit of the upward climbing rate, Indicates the unit The upper limit of the downward climbing rate, This represents the set of CHP units; this constraint reflects the limitation of the physical characteristics of thermal power units on the speed of output adjustment, avoiding equipment damage or operational instability caused by rapid output adjustment of units, and also helps the economical scheduling of the system, preventing increased operating costs due to frequent adjustments of unit output. Wind power output constraints: The power generation of a wind farm is limited by the available wind power. (12) This constraint specifies the power output range of the wind turbine in each time period; the lower limit is 0, indicating that the wind turbine can stop generating electricity; the upper limit... This is determined by the wind conditions at the time, reflecting the renewable energy characteristics and uncertainties of wind power; this constraint ensures that the power generation of wind turbines does not exceed their maximum power generation capacity in a specific period, which is in line with the actual operation of wind power. Unit state constraints (13) Generator set state constraints describe the dynamic relationships between the unit's operating state, startup state, and shutdown state; among them, Indicates the unit In time The running status, Indicates the unit In time The startup state, Indicates the unit In time The equation shows that the change in the unit's operating state during the current period is equal to the difference between its start-up and shutdown states, thus ensuring the continuity and logical consistency of the unit's state and providing a basis for accurately modeling the unit's start-up and shutdown process. Minimum start-stop time constraints (14) (15) The minimum start-up and shutdown time constraint ensures that once the unit starts up or shuts down, it will remain in that state for at least a certain period of time; the first inequality stipulates that, within time... Before During the time period, the unit The number of times a unit can be started cannot exceed its operating state; that is, if the unit is currently in an operating state, then the number of times it has been started in the past... The first inequality ensures that the unit is started at least once within a certain time period; the second inequality ensures that if the unit is currently shut down, then it was started at least once in the past. The unit can be shut down at most once within a given time period; these constraints reflect the physical requirements for hot and cold starts of the unit, avoiding damage to the equipment caused by frequent start-stops, while also helping to reduce operating costs and improve system reliability. Variable bounds constraints (16) (17) The upper and lower bound constraints of variables clearly define the range of values ​​for each decision variable; for the operating status of the unit... It can only take the value 0 or 1, indicating whether the unit is running at time t; startup and shutdown status variables. and The value of is between 0 and 1, but in actual operation it is usually 0 or 1 to ensure the integer programming characteristics of the model; for the rotated spare variable and The range of values ​​is limited to 0 to the maximum upward and downward ramp rates of the unit, ensuring the rationality and feasibility of the reserve capacity. Network constraints: The power flow of transmission lines should not exceed the transmission capacity. Equation (18) represents the nodal power equation, Equation (19) represents the branch power flow equation, and Equation (20) defines the phase angle of the reference bus; where, It is with nodes The set of all directly connected nodes It is a node The upper connected conventional generator set, It is a node The collection of combined heat and power units connected to the top, It is a node The upper-connected wind turbine assembly, For nodes The voltage phase angle at time t, For nodes The voltage phase angle at time t, This represents the voltage phase angle of the reference node in the network at time t. Let (i,j) be the reactance of the line. For the line ( i , j Maximum transmission capacity, This is the set of all routes; （18) (19) (20) Network constraints ensure that the power system meets the physical limitations and safety requirements of power transmission during operation; the node power equations link the node voltage phase angle difference with the power transmitted through the lines, reflecting the basic laws of power flow in the power system; the branch power flow equations limit the power transmission capacity of each transmission line, prevent line overload, and ensure the safe and stable operation of the system; the setting of the reference bus phase angle provides a benchmark for the voltage phase angle of the entire system, facilitating calculation and analysis; these constraints together constitute the operating framework of the power network, ensuring the efficiency and reliability of power transmission.

5. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 1, characterized in that, The dynamic modeling method for the thermal network in step 101 includes the DHN model and its complete analysis method: DHNs are two-layer networks that integrate supply and return water networks. Heat energy is generated from the heat source, transported through the supply pipeline, distributed at the heat exchange station, and then flows to the heat users. After heat exchange at the user end, the water flows back to the heat source through the return water network. A typical DHN model includes both hydraulic and thermal components. The mass regulation mode of DHNs under constant mass flow rate meets the heat load through temperature regulation. Temperature is defined as relative temperature, that is, the absolute temperature difference between the water flow and the ambient temperature. The mass model formula for DHN under mass conditioning is as follows: (21) (22) (23) (24) in, This represents the longitudinal distance of a point in the pipeline from the pipeline inlet, i.e., the spatial coordinates along the pipeline. Indicates time, Indicates the first j The temperature of the pipe, Indicates the first j The water flow velocity in the pipe This indicates the specific heat capacity of water. Indicates the first j Mass flow rate of the pipeline Indicates the first j The heat loss coefficient of the pipeline, Represents a node thermal power, Represents a node Mass flow rate Represents a node Water supply pipe temperature Represents a node Return water pipe temperature Describes the set of pipes flowing into node j. Describes the set of pipes that flow out of node j. This represents the temperature at the outlet of pipe b, where L is the length of the pipe. This represents the mass flow rate of pipe b. This indicates the temperature at the inlet of pipe b. This represents the mixing temperature at node k. Represents a set of pipes. Represents a set of nodes. This represents the set of pipelines that start at node k and end at node i. The temperature conduction equation represents the temperature distribution along a pipe, and this equation is constructed based on the following assumptions: 1) The heat loss of the pipe cross-section is considered constant; 2) Ignore heat conduction between adjacent water flows; This equation models the nodes in the DHN as heat exchangers, thereby linearly relating mass flow rate and temperature; The nodes follow the law of conservation of energy, and the mixing temperature of a node is a linear weighted sum of the outlet temperatures of the pipes injected into that node; finally, the node mixing temperature of the pipes leaving the node is the same as the inlet temperature.

6. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 5, characterized in that, The heat conduction equation in a pipe is a typical first-order linear and quasi-linear partial differential equation. Quasi-linear partial differential equations (PDEs) are commonly used to describe vibration and wave phenomena and their dynamic processes; their representative form is: (25) Where a, b, and c are the coefficients of the first-order linear hyperbolic partial differential equation. Specifically, a corresponds to the rate of change of temperature in space, b corresponds to the rate of change of temperature in time, and c corresponds to the temperature distribution itself. Let x be the state variables of a first-order linear hyperbolic partial differential equation, where x is the spatial coordinate and t is the time variable. This represents the temperature distribution function, i.e., the temperature at pipe location x and time t. These are the initial conditions for hyperbolic partial differential equations. This represents the initial temperature distribution at position x; The solution plane of this equation consists of two-dimensional vectors in space and time, and is also known as the Cauchy problem; The core idea of ​​the two-sided feature line method is to find the feature lines. Along these characteristic lines, the original PDE can be transformed into an ordinary differential equation (ODE). From this perspective, the solution plane can be considered as being formed by the initial condition curves on... The surface formed by translating along the feature line; the abbreviation for the bilateral feature line method is: CLM.

7. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 6, characterized in that, The steps of traditional CLM are as follows: The differential of u with respect to t is expressed as: (26) Substituting the above expression into the pipe temperature conduction equation, we observe that when... When this happens, equation (25) can be converted to ODE; therefore, The solution is defined as the characteristic line of equation (25), expressed as: (27) Where c is a constant determined by the initial conditions; Replacing the partial derivative with the total differential, the reconstructed ODE is expressed as: (28) The analytical solution to the above equation can be obtained from the coefficients. Different expressions and initial conditions get; To solve the Cauchy problem using CLM, a set of initial conditions is needed as the boundary of the solution plane; in other words, the PDEs that can be solved by traditional CLM must be unilaterally bounded along the feature line direction. Considering a practical scenario, where... and Both are bounded, and the bilaterally bounded hyperbolic PDE is extended as follows: (29) in, Describes the boundary conditions for hyperbolic partial differential equations. This indicates the temperature setpoint at the heat source end of the pipeline; Unlike traditional CLM, its solution is determined by both initial and boundary conditions, expressed as: (30) Determined by initial conditions, this indicates the temperature distribution caused by the initial heat storage in the pipeline; Determined by the boundary conditions, it represents the temperature propagation distribution caused by real-time heating from the heat source; the solution plane can be considered as the initial condition plane along... The result of the intersection of the surface formed by the movement and the surface formed by the movement of the boundary conditions; (31) (32) in, Right now , This represents the temperature distribution caused by the initial heat storage in the pipe at the initial time t=0. Right now , This represents the temperature propagation distribution at x=0 at the beginning of the pipeline, where heat is supplied in real time by a heat source.

8. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 5, characterized in that, A deterministic optimization method is used to construct a model of an integrated electric-thermal energy system, as detailed below: , , The English term for the integrated electric-thermal energy system is Unit Commitment with Combined Electricity and District Heating Networks (UC-CEHN). To solve the optimization model of the integrated electric-thermal energy system, the Pendesic decomposition technique is introduced to construct a distributed decision-making framework between the grid operator and the district heating operator. Under this framework, the Power Grid Operator (PGO) is mainly responsible for solving the main UC problem and checking the feasibility of grid network constraints, while the District Heating Operator (DHO) focuses on checking the feasibility of DHN constraints. This decomposition-coordination mode not only improves computational efficiency but also protects the information privacy of the two networks to a certain extent. The English abbreviation for grid operator is PGO, and the English abbreviation for district heating operator is DHO. The English term for the main UC problem is Unit Commitment Master Problem (UCMP).

1. Main UC problem: Provide a combination and scheduling plan for power generation equipment to minimize the total operating cost in (1) while satisfying the constraints in (2)-(20), (21)-(24) and the feasible cutting plane provided by the feasibility check; As a MILP, the main UC problem can be solved using the Lagrange relaxation method or the MILP solver; 2. Network security checks: Solving these sub-problems is to check the plans of a given crew. Feasibility of applying network constraints; since the network constraints at different times are independent, the subproblems representing these constraints can be solved in parallel. Each subproblem is described as follows: (33) (34) (35) in, The optimal value of the network security inspection subproblem represents the total number of violations of power flow constraints under a given unit schedule in time period t. This represents the planned active power output vector of the generating units during time period t; This represents the planned start-up and shutdown status vector of the generating units during time period t; This is the first relaxation variable, representing a virtual positive power injection used to compensate for the power deficit at this node. The second slack variable represents the virtual negative power absorption, used to absorb the power excess at the node. The two slack variables are used to ensure that the node power balance equation is always solvable. if For a given unit plan Network constraints are feasible; if Then, add the following Pendesic cut plane to UCMP: (36) in, This represents the vector of unit output variables for time period t. This represents the vector of unit start-up and shutdown status variables for time period t; DHN Feasibility Check Subproblem: The purpose of solving this subproblem is to check the plan for a given unit. The feasibility of the following DHN operating constraints; its calculation formula is as follows: (37) （38） Constraints (2)-(20), (21)-(24); in, The optimal value of the DHN feasibility check subproblem is represented by , and the total number of violations of thermal network constraints is represented by , under a given unit plan. Represents the thermal output planning vector. This represents the planned start-up and shutdown state vector of the thermal power unit. It is the third relaxation variable, representing a virtual positive heat power injection used to compensate for the heat power deficit at this node. The fourth relaxation variable represents the virtual negative heat power absorption, used to absorb the excess heat power at this node. The two relaxation variables are used to ensure that the heat power balance equation of thermal node j is always solvable at time t. Represents the set of all nodes in a heat network. This represents the set of all heat sources supplying heat to node j at time t. The minimum thermal output coefficient of heat source g, This represents the mass flow rate of node j at time t. This indicates the specific heat capacity of water. This represents the water supply temperature setpoint for node j at time t. This represents the setpoint for the return water temperature at node j at time t; if For a given unit plan DHN constraints are feasible; if Then the following Pendesic cut plane will be fed back into the main problem: （39） in, Represents the vector of thermal output variables. This represents the vector of start-up and shutdown state variables for a thermal power unit. Through the detailed explanation of the above model and sub-problems, it is known that the optimization model of the integrated electric-thermal energy system not only needs to consider the operating characteristics and cost factors of the power system, but also needs to meet the physical constraints and operating requirements of the thermal network.

9. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 5, characterized in that, The computer group combinatorial problem based on the improved Pendesic decomposition is as follows:

1. Traditional Pendes decomposition algorithm: The Pendes decomposition algorithm decomposes the original problem into a master problem and subproblems, and then uses an iterative approach to gradually approximate the optimal solution of the original problem. The master problem is called the Master Problem (MP), and the subproblems are called subproblems (SP). Suppose the original problem is a mixed-integer linear programming problem, which is expressed as: （40） in, It is an integer variable. It is a continuous variable. It is the coefficient vector corresponding to the integer variables in the objective function. express transpose, It is the coefficient vector corresponding to the continuous variables in the objective function. express transpose; It is a constraint matrix corresponding to integer variables. It is the constraint matrix corresponding to continuous variables. It is a constraint vector; The core idea of ​​Pendes decomposition is to decompose the original problem into MP and SP; MP contains only integer variables. SP, on the other hand, deals with continuous variables. By iteratively solving MP and SP, cutting planes (Cuts) are gradually generated, thus approximating the optimal solution to the original problem; the objective of MP is to minimize integer variables. The objective function value, while also considering the optimal objective value of SP, is in the form of: （41) in, This represents the optimal objective value of SP; SP is a continuous variable. The linear programming problem is of the form: (42) The goal of SP is based on what MP gives. Solve for the optimal And calculate the objective function value; The Pendes decomposition algorithm is a classic optimization technique that solves mixed integer programming problems by decomposing the original problem into MP and SP components. The algorithm's implementation involves five key steps: First, starting from an initial integer variable... First, solve for SP; second, for a given... Solve for SP to obtain the optimal objective value. Next, a cutting plane is generated. If SP has no feasible solution, a feasible cut is added to MP. This cut is generated based on the dual information of SP to ensure that the MP has a feasible solution. This is feasible in SP. If SP has a feasible solution, an optimality cut is added to MP. This cut, also based on the dual information of SP, is used to constrain the numerical value of the objective function of MP. Then, MP is updated, the generated cut plane is added to MP, and MP is solved again to obtain a new solution. Finally, the iterative process repeats the above steps until the convergence condition is met, such as the difference between the target values ​​of MP and SP being less than a certain threshold. Through this series of steps, the Pendes decomposition algorithm can effectively approximate the optimal solution of the original problem.

2. Improved Bends decomposition algorithm using interior point method To address the shortcomings of the traditional Pendes decomposition algorithm in terms of convergence speed and the requirement for high-quality cutting planes, the interior point method and the Pendes algorithm are combined. The interior point method, as an algorithm for solving convex optimization problems, explores the optimal solution path within the feasible region to approach the optimal solution. The Pendes IPM algorithm is implemented through the following steps: First, select an initial interior point. The interior point is usually taken as the middle value of the variable to ensure that the initial point is located inside the feasible region; this choice helps the algorithm to start the search from the central region of the feasible region, which may improve the convergence speed. Next, based on the current interior point Solve the dual SP and add constraints generated by the dual SP to MP; these constraints help guide the search direction of MP, making it closer to the optimal solution; then, after adding the new constraints, solve the relaxed MP, i.e., ignore the integer constraints, to obtain the extreme points. The purpose of this step is to find a continuous variable solution that minimizes the objective function value; then, based on the extreme points... Update Inner Points Use the formula: (43) in is an interior-point update coefficient, representing a weighted average between the current solution and the extreme points to achieve stable convergence; k represents the number of iterations. Let the interior point be represented in the k-th iteration. This represents the interior point at the (k+1)th iteration, i.e., the updated interior point; The extreme point obtained by solving the relaxed MP represents the optimal solution obtained by solving the relaxed master problem after adding the new constraints generated by the dual subproblem. Finally, the above process is repeated until the convergence condition is met.

10. The method for combining units in an integrated electrothermal energy system based on interior-point Bends decomposition according to claim 5, characterized in that, The UC-CEHN model is solved using the Pendes IPM algorithm; MP is responsible for solving the main UC problem and verifying the feasibility of the power grid constraints, while SP is responsible for verifying the feasibility of the DHN constraints. a)MP: The main problem part updates the interior points according to the given interior point update coefficients, checks the feasibility in (21)-(24), and minimizes the total operating cost of the generator set commitment and scheduling plan in (1) under the condition that the feasibility cut provided by the constraints in (2)-(20) is satisfied. b) SP: Check the feasibility of network constraints by solving these subproblems (33)-(39); since network constraints at different times are independent, it means that their subproblems can be solved in parallel.

Images