Linear approximation natural gas pipeline network scheduling optimization method considering hydraulic constraint

By linearly approximating the nonlinear terms in hydraulic constraints, large-scale natural gas transportation problems are transformed into linear planning problems, solving problems that are too long in the existing technology, and quickly solving high-quality natural gas transportation plans are achieved.

CN120181272APending Publication Date: 2025-06-20PETROCHINA CO LTD
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202311741510.X
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2023-12-18
Publication Date
2025-06-20

AI Technical Summary

Technical Problem

In large-scale natural gas transportation pipeline networks, it is difficult for the existing technology to quickly solve natural gas transportation plans, resulting in too long solution time and reducing the real-time and flexibility of the system.

Method used

By performing a linear approximation of the first-order Taylor expansion of the nonlinear term in the hydraulic constraint near the basis point, the nonlinear mixed integer programming problem is degenerated into a linear mixed integer programming problem, so as to quickly find the local optimal solution of the original problem.

Benefits of technology

It has achieved high-quality natural gas transportation plans in large-scale natural gas pipelines within hourly times, which has significantly improved the real-time and flexibility of the system.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN120181272A_ABST
    Figure CN120181272A_ABST
Patent Text Reader

Abstract

The invention provides a linear approximation natural gas pipeline network scheduling optimization method considering hydraulic constraint, and the method comprises the steps: obtaining parameters needed for solving a nonlinear mixed integer programming problem of natural gas pipeline transportation, and substituting the parameters into the nonlinear mixed integer programming problem of natural gas pipeline transportation for solving; the nonlinear mixed integer programming problem comprises discontinuous factors and steady-state hydraulic constraints; the discontinuous factors at least comprise a steady-state hydraulic constraint, a stepped supply price function and a pipeline direction; the steady-state hydraulic constraint at least comprises a non-linear term about pressure variables at the two ends of the pipeline and a non-linear term about flow variables; the step of solving the nonlinear mixed integer programming problem comprises the substeps that linear approximation is conducted on nonlinearity in hydraulic constraint near a base point, and the nonlinear mixed integer programming problem is degraded into a linear mixed integer programming problem. According to the method provided by the invention, the optimal large-scale natural gas pipeline network scheduling method can be solved within the hour level.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This document relates to, but is not limited to, the field of natural gas scheduling, and in particular, but not limited to, a linear approximation method for optimizing the scheduling of natural gas pipeline networks considering hydraulic constraints. Background Art

[0002] Affected by the macro - economy, temperature, and the rising prices of alternative energy sources, as well as the strong promotion of environmental protection policies such as "coal - to - gas conversion" and "clean heating", natural gas, as a clean energy source, accounts for an increasing proportion in China's energy consumption structure, and the absolute consumption of natural gas is also increasing. During the pipeline transportation of natural gas, the gas throughput of the pipeline and the pressures at both ends of the pipeline are the main decision variables, with complex physical relationships. This physical relationship is generally referred to as hydraulic constraints, which are relatively complex non - linear relationships. As the proportion of natural gas in energy use increases, the scale of natural gas pipeline networks is also constantly expanding. Secondly, since the supply price of natural gas is generally a step - function related to the purchase volume, the formulation of natural gas transportation plans also includes discontinuous variables such as pipeline directions and compressor station switches. Therefore, there are a large number of combinatorial factors in natural gas pipeline transportation problems. The huge scale combined with complex physical relationships makes it very difficult to solve the transportation plan of natural gas pipeline networks.

[0003] Most of the optimization problems of natural gas pipeline transportation plans are represented by non - linear mixed - integer programming problems (NLMIP). In existing applications or research, there has been no successful case of solving the natural gas transportation plan for the entire network in such a large - scale natural gas transportation pipeline network. Most of the application objects or research objects are small and medium - sized national or regional natural gas transportation pipeline networks. In these studies, a network flow model is established for the natural gas transportation pipeline network, and the hydraulic constraints under steady state are added, and a commercial solver for quadratic mixed - integer programming is used to solve the problem.

[0004] For the above - mentioned solution, if applied to China's large - scale natural gas transportation pipeline network, when solved with the current most advanced quadratic mixed - integer programming solver, the solving time needs to reach more than ten hours. This greatly reduces the real - time performance and flexibility of the solution. Moreover, during the process of users formulating natural gas transportation plans using the natural gas pipeline network transportation plan optimization system, they often need to try to adjust the parameters to obtain an ideal transportation plan. The overly long solving time greatly reduces the usability of the system. Summary of the Invention

[0005] The following is an overview of the topics described in detail in this document. This overview is not intended to limit the scope of protection of the claims.

[0006] The present disclosure formulates a non - linear mixed - integer programming problem for natural gas pipeline transportation issues, which includes steady - state hydraulic constraints. These constraints contain non - linear terms related to pressure variables at both ends of the pipeline and non - linear terms related to flow variables. Secondly, in this problem, discontinuous factors such as step - supply price functions and pipeline directions are also included. When solving this problem, the non - linearity in the hydraulic constraints is linearly approximated by a first - order Taylor expansion near the base point. In this way, the non - linear mixed - integer programming problem is reduced to a linear mixed - integer programming problem, and the local optimal solution of the original non - linear problem can be quickly found, so as to achieve the purpose of quickly (hour - level) obtaining a high - quality natural gas transportation plan for large - scale natural gas pipeline networks.

[0007] The present disclosure provides a linear - approximation natural gas pipeline network scheduling optimization method considering hydraulic constraints. The method includes the following steps:

[0008] Obtain the parameters required for solving the non - linear mixed - integer programming problem of natural gas pipeline transportation, and substitute them into the non - linear mixed - integer programming problem of natural gas pipeline transportation for solution;

[0009] The non - linear mixed - integer programming problem includes discontinuous factors and steady - state hydraulic constraints;

[0010] The discontinuous factors at least include: steady - state hydraulic constraints, step - supply price functions, and pipeline directions.

[0011] The steady - state hydraulic constraints at least include: non - linear terms related to pressure variables at both ends of the pipeline and non - linear terms related to flow variables;

[0012] Solving the non - linear mixed - integer programming problem includes: linearly approximating the non - linearity in the hydraulic constraints by a first - order Taylor expansion near the base point, and reducing the non - linear mixed - integer programming problem to a linear mixed - integer programming problem.

[0013] In some embodiments provided by the present application, for the linear - approximation natural gas pipeline network scheduling optimization method considering hydraulic constraints, the method includes the following steps:

[0014] Set the original problem as maximizing the sales profit = total sales amount at natural gas demand points - total purchase amount at natural gas supply points - transportation cost, that is, formula (1);

[0015] 1) Read natural gas pipeline network data and supply - demand price data;

[0016] 2) Take all flow values as 0 as the initial base point, set the expansion radius ε as 1, and substitute it into the objective function of the first stage:

[0017] The constraint formulas that the objective function in the first stage obeys are as follows: Formula (2), Formula (3), Formula (4), Formula (5), Formula (8), Formula (9), Formula (6), Formula (13), Formula (14), Formula (15), Formula (16), Formula (17), Formula (18), Formula (19), Formula (20), Formula (22.1), Formula (23), Formula (24), Formula (25);

[0018] In the objective function of the first stage, represents the relaxed positive term, represents the relaxed negative term (here, this value will be comprehensively obtained through continuous iterative solution and approximated to 0 as much as possible);

[0019] Use a linear mixed-integer planner to solve the problem and record the optimal value of the problem;

[0020] 3) Determine whether the optimal value of the objective function of the first stage solved in step 2) is close enough to the optimal value in the previous iteration when solving the objective function of the first stage (for example, when the gap between two iterations is small enough (generally taking 0.0001), it is considered relatively stable). If it is close enough, reduce the expansion radius ε to 0.5 times that in the previous iteration. Determine whether the expansion radius is small enough. If it is not small enough, update the base point to the pipeline flow value in the current optimal solution and substitute it into step 3) to solve again;

[0021] 4) If the expansion radius described in step 3) is small enough, determine whether it is in the first stage (during the calculation process, an identifier will be set for each stage, and whether it is in the first stage can be judged through the identifier; it is identified by setting variables in the corresponding stage during the algorithm solution process. Each time a new stage is entered, the current stage is set to the new stage). If it is not in the first stage, determine whether it is in the second stage. If it is not in the second stage, then complete the solution of the original problem, and the parameters at this time can maximize the sales profit;

[0022] 5) If the conclusion in step 4) that determines whether it is in the first stage is yes, then determine whether the optimal value is 0. If not, then the original problem has no solution. If it is, substitute it into the objective function of the second stage: In this case, update the base point to the pipeline flow value in the current optimal solution (using the optimal solution when the first-stage problem terminates as the starting base point);

[0023] The constraint formulas that the objective function in the second stage obeys are as follows: Formula (2), Formula (3), Formula (4), Formula (5.1), Formula (8), Formula (9), Formula (6), Formula (13), Formula (14), Formula (15), Formula (16), Formula (17), Formula (18), Formula (19), Formula (20), Formula (22), Formula (23), Formula (24), Formula (25);

[0024] Update the base point to the pipeline flow value in the current optimal solution, and substitute it into Steps 2) to 4) for solution;

[0025] 6) If the conclusion in Step 4) that it is in the second stage is yes, then judge whether the optimal value is 0. If not, then the original problem has no solution; if so, then substitute it into the objective function of the third stage: Formula (1). Update the base point to the pipeline flow value in the current optimal solution (using the optimal solution when the second-stage problem terminates as the starting base point), and substitute it into Step 1) for solution;

[0026] The constraint formulas that the said Formula (1) obeys are as follows: Formula (2), Formula (3), Formula (4), Formula (5), Formula (8), Formula (9), Formula (6), Formula (13), Formula (14), Formula (15), Formula (16), Formula (17), Formula (18), Formula (19), Formula (20), Formula (22), Formula (23), Formula (24), Formula (25);

[0027] Formula (1) is as follows:

[0028]

[0029] In Formula (1), I is the set of demand nodes; O is the set of gas source nodes, PI is the set of ordinary pipe segments (excluding compressor stations and regulating valves); CP is the set of pipe segments containing compressor stations; R ijm is the total sales amount in the previous m - 1 stages when the usage reaches the mth ladder (parameter), with the unit of ten thousand yuan; pijm is the unit price (parameter) of the jth demand of node i at the mth ladder, with the unit of yuan per ten thousand cubic meters; is the supply unit price (parameter) of the jth gas source of node i, with the unit of yuan per ten thousand cubic meters; is the transportation unit price (parameter) from node i to node j, with the unit of yuan per ten thousand cubic meters;

[0030] Formula (2) is as follows:

[0031]

[0032] In Formula (2), the said v ijm represents the satisfaction amount d of the jth demand of node i ijWhether it is in the m-th ladder. If it is in the m-th ladder, take the value of 1; if it is not in the m-th ladder, take the value of 0, dimensionless;

[0033] Formula (3) is as follows:

[0034]

[0035] In formula (3), q ijm is the quantity (parameter) of the j-th demand of node i at the m-th ladder, with the unit of 10,000 cubic meters;

[0036] Formula (4) is as follows:

[0037]

[0038] In formula (4), Q ijm is the total sales volume (parameter) in the previous m - 1 stages when the usage reaches the m-th ladder, with the unit of 10,000 cubic meters; d ij is the satisfied quantity of the j-th demand of node i, with the unit of 10,000 cubic meters;

[0039] Formula (5) is as follows. The d in formula (4) ij satisfies formula (5):

[0040]

[0041] In formula (5), is the lower limit (parameter) of the j-th demand of node i, with the unit of 10,000 cubic meters; is the upper limit (parameter) of the j-th demand of node i, with the unit of 10,000 cubic meters;

[0042] Formula (5.1) is as follows. The d in formula (4) ij satisfies formula (5.1):

[0043]

[0044] In formula (5.1), is the lower limit (parameter) of the j-th demand of node i, with the unit of 10,000 cubic meters; is the upper limit (parameter) of the j-th demand of node i, with the unit of 10,000 cubic meters; ξ ij is the slack term of the demand, with the unit of 10,000 cubic meters;

[0045] Formula (6) is as follows. In formula (9), o i is the total procurement volume of node i, with the unit of 10,000 cubic meters, and is calculated by formula (6):

[0046]

[0047] Formula (7) is as follows. In formula (9), di is the total demand of node i, in 10,000 cubic meters, calculated by formula (7):

[0048]

[0049] Formula (8) is as follows:

[0050]

[0051] In formula (8), o ij is the procurement volume of the j-th gas source of node i, in 10,000 cubic meters; is the lower supply limit (parameter) of the j-th gas source of node i, in 10,000 cubic meters; is the upper supply limit (parameter) of the j-th gas source of node i, in 10,000 cubic meters;

[0052] Formula (9) is as follows:

[0053]

[0054] In formula (9), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, in 10,000 cubic meters (f in formula (1)); ij is a scalar without considering direction. In formula (9), f ij and f ji represent the flow rates in both directions of the pipe section; N is the set of all network nodes; VP is the set of pipe sections containing regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable;

[0055] Formula (13) is as follows: P in formula (22), formula (22.1), formula (23) and formula (24); i and P j satisfy formula (13):

[0056]

[0057] In formula (13), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe sections containing compressor stations; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable;

[0058] Formula (14) is as follows: P in formula (22), formula (22.1), formula (23) and formula (24); i and P j also satisfy formula (14):

[0059]

[0060] In formula (14), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe segments containing compressor stations; x ij takes a value of 0 or 1. If the gas flow is from i to j, the value is 1; otherwise, it is 0. It is the gas flow direction variable;

[0061] Formula (15) is as follows: P in formula (22), formula (22.1), formula (23), and formula (24) i and P j also satisfy formula (15):

[0062]

[0063] In formula (15), P i is the pressure value of node i; P j is the pressure value of node j; P i ub is the upper limit of the pressure value of node i (parameter); y ij takes a value of 0 or 1. If the compressor station bypasses, the value is 1; otherwise, it is 0; x ij takes a value of 0 or 1. If the gas flow is from i to j, the value is 1; otherwise, it is 0. It is the gas flow direction variable;

[0064] Formula (16) is as follows: P in formula (22), formula (22.1), formula (23), and formula (24) i and P j also satisfy formula (16):

[0065]

[0066] In formula (16), P i is the pressure value of node i; P j is the pressure value of node j; is the lower limit of the compression ratio of compressor station i (parameter); y ij takes a value of 0 or 1. If the compressor station bypasses, the value is 1; otherwise, it is 0;

[0067] Formula (17) is as follows: P in formula (22), formula (22.1), formula (23), and formula (24) i satisfies formula (17):

[0068]

[0069] In formula (17), N is the set of all network nodes; PI is the set of ordinary pipe segments (excluding compressor stations and regulating valves); P i is the pressure value of node i; P i lb is the lower limit (parameter) of the pressure value of node i; P i ub is the upper limit (parameter) of the pressure value of node i;

[0070] Formula (18) is as follows:

[0071]

[0072] In formula (18), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 cubic meters; VP is the set of pipe segments containing regulating valves; PI is the set of ordinary pipe segments (excluding compressor stations and regulating valves); is the upper limit (parameter) of the transportation volume of the pipe segment from node i to node j;

[0073] Formula (19) is as follows:

[0074]

[0075] In formula (19), the v ijm represents whether the satisfaction amount d ij of the jth demand of node i is in the mth ladder. If it is in the mth ladder, it takes a value of 1; if it is not in the mth ladder, it takes a value of 0, dimensionless;

[0076] Formula (20) is as follows:

[0077]

[0078] In formula (20), the u ijm is a non - negative real - valued variable representing the amount of d ij falling into the mth ladder, with the unit of 10,000 cubic meters;

[0079] Formula (22) is as follows:

[0080]

[0081] In formula (22), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 cubic meters; x ij takes a value of 0 or 1. If the air flow is from i to j, it takes a value of 1; otherwise, it takes a value of 0, which is the air flow direction variable; f ij,0 is the flow value at the base point, in 10,000 cubic meters; P i is the pressure value of node i, with the unit of Pascal;j The pressure value of node j, with the unit of Pascal; where, λ is the friction coefficient between the gas in the pipeline and the inner wall of the pipe, with the unit of dimensionless; Z is the gas compressibility factor, with the unit of dimensionless; Δ * is the relative density of the gas, with the unit of dimensionless; T is the average temperature of the pipeline, with the unit of Kelvin; L is the length of the pipeline, with the unit of meter; C0 takes the constant 0.03848; D is the inner diameter of the pipeline, with the unit of meter; where, the parameter where, g is the acceleration of gravity; Z is the gas compressibility factor, dimensionless; T is the average temperature of the pipeline, with the unit of Kelvin; R is the gas constant; h i is the elevation of the i-th section of the pipeline, with the unit of meter; L i is the length of the i-th section of the pipeline, with the unit of meter; β = a·Δh, where, a is the same as a in calculating θ, and Δh is the difference in elevation between the starting point and the ending point of the pipe section, with the unit of meter;

[0082] Formula (22.1) is as follows:

[0083]

[0084] The parameter meanings in Formula (22.1) are exactly the same as those in Formula (22);

[0085] Formula (23) is as follows:

[0086]

[0087] In Formula (23), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of the direction, with the unit of 10,000 cubic meters; VP is the set of pipe sections containing regulating valves; x ij takes the value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable; f ij,0 is the flow value at the reference point, with the unit of 10,000 cubic meters; P i is the pressure value of node i, with the unit of Pascal; P j is the pressure value of node j, with the unit of Pascal; γ, θ, and β in Formula (23) are the same as those in Formula (22);

[0088] Formula (24) is as follows:

[0089]

[0090] In Formula (24), VP is the set of pipe sections containing regulating valves; x ij takes the value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable; f ij,0is the flow value at the base point, with the unit of 10,000 m³; P i is the pressure value of node i, with the unit of Pa; P j is the pressure value of node j, with the unit of Pa; γ, θ, and β in formula (24) are the same as those in formula (22);

[0091] Formula (25) is as follows:

[0092]

[0093] In formula (25), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 m³; VP is the set of pipe segments with regulating valves; f ij,0 is the flow value at the base point, with the unit of 10,000 m³; ε is set artificially according to the convergence situation during the calculation, representing the difference between the unsolved solution and the initial solution; f ij,0 is the flow value at the base point, with the unit of 10,000 m³.

[0094] In some embodiments provided by the present disclosure, in step 3), the "close" means that the solution result is less than 10 -4 to the optimal value in the previous iteration.

[0095] In some embodiments provided by the present disclosure, in step 3), the "sufficiently small" means that the solution result is less than 10 -4 to the optimal value in the previous iteration.

[0096] In some embodiments provided by the present disclosure, the value range of ε is from 0.005 to 0.0005, and it can also be appropriately adjusted according to the conventional calculation accuracy requirements in the art.

[0097] Other features and advantages of the present disclosure will be described in the subsequent description, and some of them will become obvious from the description, or be understood by implementing the present disclosure. Other advantages of the present disclosure can be achieved and obtained through the solutions described in the description. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] The drawings are used to provide an understanding of the technical solutions of the present disclosure, and constitute a part of the description. Together with the embodiments of the present disclosure, they are used to explain the technical solutions of the present disclosure, and do not constitute a limitation to the technical solutions of the present disclosure.

[0099] Figure 1 is the flowchart of the embodiment of the present disclosure. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0100] To make the objectives, technical solutions, and advantages of the present disclosure more clear and understandable, the embodiments of the present disclosure are described in detail below. It should be noted that, without conflict, the embodiments in the present disclosure and the features in the embodiments may be arbitrarily combined with each other.

[0101] In an exemplary embodiment of the present disclosure, a linear approximation natural gas pipeline network scheduling optimization method considering hydraulic constraints is provided. The method includes the following steps:

[0102] Obtain the parameters required for solving the non-linear mixed integer programming problem of natural gas pipeline transportation, and substitute them into the non-linear mixed integer programming problem of natural gas pipeline transportation for solution;

[0103] The non-linear mixed integer programming problem includes discontinuous factors and steady-state hydraulic constraints;

[0104] The discontinuous factors at least include: steady-state hydraulic constraints, step supply price function, pipeline direction;

[0105] The steady-state hydraulic constraints at least include: non-linear terms regarding the pressure variables at both ends of the pipeline and non-linear terms of the flow rate variables;

[0106] Solving the non-linear mixed integer programming problem includes: making a linear approximation of the non-linearity in the hydraulic constraints by a first-order Taylor expansion near the base point, and degrading the non-linear mixed integer programming problem into a linear mixed integer programming problem.

[0107] The original problem can be set as maximizing the sales profit = total sales amount at the natural gas demand points - total procurement amount at the natural gas supply points - transportation cost, that is, formula (1), which is a linear approximation problem of the natural gas pipeline network transportation scheduling problem;

[0108] As Figure 1 shown, the solution process of the original problem can be divided into three stages. In the first stage, the hydraulic constraint formula (22.1) of the pipeline is relaxed. In the second stage, the demand lower bound constraint formula (5.1) is relaxed.

[0109] In the first stage, taking the point where all flow values are zero as the base point, a feasible solution for the hydraulic constraints is found by minimizing the relaxation variables of the hydraulic constraints. The following is the relaxation problem in the first stage (hereinafter simply referred to as the first-stage problem):

[0110] Stage 1:

[0111] By iteratively obtaining the optimal solution within the neighborhood of the base point and updating the base point until the optimal value of the objective function reaches 0, a base point that meets the constraint formula (22.1) is obtained. Otherwise, it indicates that the original non-linear problem has no solution. Enter the second stage. In the second stage, taking the optimal solution at the end of the first-stage problem as the base point, a base point that meets the demand lower bound constraint is obtained.

[0112] The following is the relaxation problem in Phase II (hereinafter referred to as the Phase II problem):

[0113] By iteratively obtaining the optimal solution within the neighborhood of the base point and updating the base point, when the objective value does not change significantly compared to the optimal objective value in the previous iteration, the neighborhood radius is shrunk. When it shrinks to a small enough value, the solution of the Phase II problem is terminated. If the optimal solution of the model at termination is 0, a base point that satisfies Constraint Formula 5 can be obtained; otherwise, it indicates that the original problem has no solution. Enter Phase III, and use the optimal solution at the termination of the Phase II problem as the starting base point to solve the approximate problem. After reaching the termination conditions described above, the local optimal solution of the original non-linear problem can be obtained after completing Phase III.

[0114] In this method, by obtaining the optimal solution of the problem within a certain range near the base point, then updating the base point, and obtaining the optimal solution near the new base point again. When the optimal solution at the current base point is the same as the optimal solution at the previous base point, the neighborhood radius ε is shrunk. When the neighborhood radius ε is less than the preset value, it indicates that the original non-linear problem has reached the local optimal solution.

[0115] The specific steps may include:

[0116] 1) Read the natural gas pipeline network data and supply-demand price data;

[0117] 2) Use all flow values being 0 as the initial base point, set the expansion radius ε to 1, and substitute it into the objective function of Phase I: It is a method for calculating the initial value that satisfies the hydraulic constraints in the linear approximation problem of the natural gas pipeline network transportation scheduling problem;

[0118] The constraint formulas that the objective function of Phase I obeys are as follows: Formula (2), Formula (3), Formula (4), Formula (5.1), Formula (8), Formula (9), Formula (6), Formula (13), Formula (14), Formula (15), Formula (16), Formula (17), Formula (18), Formula (19), Formula (20), Formula (22.1), Formula (23), Formula (24), Formula (25);

[0119] In the objective function of Phase I, represents the positive term of relaxation, represents the negative term of relaxation (here this value will be integrated through continuous iterative solution and approximated to 0 as much as possible);

[0120] Use a linear mixed-integer programming solver to solve the problem and record the optimal value of this problem;

[0121] 3) Determine whether the optimal value of the objective function in Phase 1 obtained in Step 2 is close enough to the optimal value in the previous iteration when solving the objective function in Phase 1 (for example, when the gap between two iterations is small enough (usually taken as 0.0001), it is considered relatively stable). If it is close enough, reduce the expansion radius ε to 0.5 times that in the previous iteration. Determine whether the expansion radius is small enough. If it is not small enough, update the base point to the pipeline flow value in the current optimal solution, and substitute it into Step 3) to solve again;

[0122] 4) If the expansion radius described in Step 3) is small enough, determine whether it is in Phase 1 (during the calculation process, an identifier is set for each phase, and whether it is in Phase 1 can be determined through the identifier. During the algorithm solving process, variables corresponding to the respective phases are set for identification. When entering a new phase, set the current phase as the new phase). If it is not in Phase 1, determine whether it is currently in Phase 2. If it is not in Phase 2, then complete the solution of the original problem, and the parameters at this time can maximize the sales profit;

[0123] 5) If the conclusion in Step 4) that determines whether it is in Phase 1 is yes, then determine whether the optimal value is 0. If not, then the original problem has no solution. If so, substitute it into the objective function of Phase 2: In it, update the base point to the pipeline flow value in the current optimal solution (using the optimal solution when the Phase 1 problem terminates as the starting base point), which is a method for calculating the initial value of the lower bound constraint of the operator demand in the linear approximation problem of the natural gas pipeline network transportation scheduling problem;

[0124] The constraint formulas that the objective function of Phase 2 obeys are as follows: Formula (2), Formula (3), Formula (4), Formula (5.1), Formula (8), Formula (9), Formula (6), Formula (13), Formula (14), Formula (15), Formula (16), Formula (17), Formula (18), Formula (19), Formula (20), Formula (22), Formula (23), Formula (24), Formula (25);

[0125] Update the base point to the pipeline flow value in the current optimal solution, and substitute it into Steps 2) to 4) to solve;

[0126] 6) If the conclusion in Step 4) that determines whether it is in Phase 2 is yes, then determine whether the optimal value is 0. If not, then the original problem has no solution; if so, substitute it into the objective function of Phase 3: In Formula (1), update the base point to the pipeline flow value in the current optimal solution (using the optimal solution when the Phase 2 problem terminates as the starting base point), and substitute it into Step 1) to solve;

[0127] The constraint formulas that the formula (1) obeys are as follows: formula (2), formula (3), formula (4), formula (5), formula (8), formula (9), formula (6), formula (13), formula (14), formula (15), formula (16), formula (17), formula (18), formula (19), formula (20), formula (22), formula (23), formula (24), formula (25);

[0128] Stage 1 and stage 2 can be regarded as the process of solving the initial feasible solution of stage 3. Their association is through the system of constraint equations. In stage 1, the pressure-flow constraint is minimized, and in stage 2, the customer slack is minimized. After these two conditions are basically satisfied, the optimal solution of the original problem (profit problem) is then obtained on this basis.

[0129] The formula (1) is as follows:

[0130]

[0131] In formula (1), I is the set of demand nodes; O is the set of gas source nodes, PI is the set of ordinary pipe segments (excluding compressor stations and regulating valves); CP is the set of pipe segments including compressor stations; R ijm is the total sales amount in the previous m - 1 stages when the usage reaches the mth ladder (parameter), with the unit of ten thousand yuan; pijm is the unit price of the jth demand of node i at the mth ladder (parameter), with the unit of yuan per ten thousand cubic meters; is the supply unit price of the jth gas source of node i (parameter), with the unit of yuan per ten thousand cubic meters; is the transportation unit price from node i to node j (parameter), with the unit of yuan per ten thousand cubic meters;

[0132] The formula (2) is as follows:

[0133]

[0134] In formula (2), the v ijm represents whether the satisfaction amount d ij of the jth demand of node i is in the mth ladder. If it is in the mth ladder, it takes the value of 1; if it is not in the mth ladder, it takes the value of 0, dimensionless;

[0135] Formula (2) ensures that the demand of each natural gas demand point must fall within and only within one sales ladder;

[0136] The formula (3) is as follows:

[0137]

[0138] In formula (3), q ijm is the quantity of the jth demand of node i at the mth ladder (parameter), with the unit of ten thousand cubic meters;

[0139] q ijm is a conventional parameter in this field. For example, the demand quantity signed between the Shanghai terminal of the node and Shanghai Gas may include two steps. For example, the quantity less than or equal to the signed contract quantity is the quantity of step 1, and the additional demand quantity exceeding the contract agreement is the quantity of step 2. q can be determined by the specific number of steps and the quantity of each step respectively ijm .

[0140] Formula (3) ensures u ijm can only take a value when the sales volume d ij falls within a certain sales step, otherwise it is 0. Specifically, when the sales volume falls within step m, v ijm takes 1, and at this time u ijm ≤q ijm , which limits the sales volume u ijm in stage m to be less than the upper limit of the sales volume within this step; when the sales volume does not fall within step m, v ijm takes 0, and at this time u ijm ≤0. Also, due to the constraint u ijm ≥0, so u ijm = 0;

[0141] Formula (4) is as follows:

[0142]

[0143] In formula (4), Q ijm is the total sales volume (parameter) in the previous m - 1 stages when the usage reaches the mth step, with the unit of 10,000 cubic meters; d ij is the satisfied quantity of the jth demand at node i, with the unit of 10,000 cubic meters;

[0144] Formula (4) ensures that if the sales volume d ij falls within step m, then the sales volume d ij is equal to the total sales volume Q ijm in the previous m - 1 stages plus the remaining sales volume u ijm within step m;

[0145] Formula (5) is as follows. d ij in formula (4) satisfies formula (5):

[0146]

[0147] In formula (5), is the lower limit (parameter) of the jth demand at node i, with the unit of 10,000 cubic meters; is the upper limit (parameter) of the jth demand at node i, with the unit of 10,000 cubic meters;

[0148] Formula (5) ensures that the sales volume d ijWithin the required upper and lower limits;

[0149] Equation (5.1) is as follows. The d in Equation (4) ij Satisfies Equation (5.1):

[0150]

[0151] In Equation (5.1), Is the lower limit (parameter) of the j-th demand of node i, in units of 10,000 cubic meters; Is the upper limit (parameter) of the j-th demand of node i, in units of 10,000 cubic meters; ξ ij Is the slack term of the demand, in units of 10,000 cubic meters;

[0152] Equation (6) is as follows. In Equation (9), o i Is the total procurement volume of node i, in units of 10,000 cubic meters, calculated from Equation (6):

[0153]

[0154] Equation (7) is as follows. In Equation (9), d i Is the total demand volume of node i, in units of 10,000 cubic meters, calculated from Equation (7):

[0155]

[0156] Equation (8) is as follows:

[0157]

[0158] In Equation (8), o ij Is the procurement volume of the j-th gas source of node i, in units of 10,000 cubic meters; Is the supply lower limit (parameter) of the j-th gas source of node i, in units of 10,000 cubic meters; Is the supply upper limit (parameter) of the j-th gas source of node i, in units of 10,000 cubic meters;

[0159] The purpose of Equation (8) is to ensure that the procurement volume o of each gas source ij Is within the required upper and lower limits;

[0160] Equation (9) is as follows:

[0161]

[0162] In Equation (9), f ij Is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, in units of 10,000 cubic meters (the f in Equation (1) ij Is a scalar without considering direction. In Equation (9), f ij And f jiFlow rates in two directions of the represented pipe segments); N is the set of all network nodes; VP is the set of pipe segments with regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0. It is the air flow direction variable;

[0163] The purpose of formula (9) is to represent the flow balance constraints of each network node, including demand points, gas source points, pressure gas station nodes, etc. Specifically, the inflow gas flow rate of a node plus its total purchase volume o i should be equal to the outflow gas flow rate of the node plus the total demand volume d i ;

[0164] Formula (13) is as follows: P in formula (22), formula (22.1), formula (23), and formula (24) i and P j satisfy formula (13):

[0165]

[0166] In formula (13), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe segments of the pressure gas station; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0. It is the air flow direction variable;

[0167] Formula (14) is as follows: P in formula (22), formula (22.1), formula (23), and formula (24) i and P j also satisfy formula (14):

[0168]

[0169] In formula (14), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe segments of the pressure gas station; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0. It is the air flow direction variable;

[0170] Both formula (13) and (14) ensure that the air flow in the pipe segments split from the pressure gas station satisfies the relationship between the transmission direction and the upstream and downstream air pressures. When the transmission direction is from i to j, the upstream air pressure P i should be less than the downstream air pressure P j , and at this time x ij is 1; when the transmission direction is from j to i, the downstream air pressure P i should be less than the upstream air pressure P j , and at this time xij is 0;

[0171] Equation (15) is as follows: P in Equation (22), Equation (22.1), Equation (23) and Equation (24) i and P j also satisfy Equation (15):

[0172]

[0173] In Equation (15), P i is the pressure value of node i; P j is the pressure value of node j; P i ub is the upper limit of the pressure ratio of node i (parameter); y ij takes a value of 0 or 1. The value is 1 when the compressor station bypasses, otherwise it is 0; x ij takes a value of 0 or 1. If the gas flow is from i to j, the value is 1, otherwise it is 0, which is the gas flow direction variable;

[0174] Equation (15) represents the upper limit restriction of the pressure ratio of the compressor station. When bypassing, y ij = 1. When the transmission direction is from i to j, x ij is 1, and the inequality is: P i ≤P j When the transmission direction is from j to i, x ij is 0, and the inequality is: P j ≤P i When not bypassing, y ij = 0. When the transmission direction is from i to j, x ij is 1, and the inequality is: When the transmission direction is from j to i, x ij is 0, and the inequality is:

[0175] Equation (16) is as follows: P in Equation (22), Equation (22.1), Equation (23) and Equation (24) i and P j also satisfy Equation (16):

[0176]

[0177] In Equation (16), P i is the pressure value of node i; P j is the pressure value of node j; is the lower limit of the pressure ratio of compressor station i (parameter); y ij takes a value of 0 or 1. The value is 1 when the compressor station bypasses, otherwise it is 0;

[0178] Equation (16) represents the lower limit restriction of the compression ratio of the compressor station;

[0179] Equation (17) is as follows: P in Equation (22), Equation (22.1), Equation (23), and Equation (24) i Satisfies Equation (17):

[0180]

[0181] In Equation (17), N is the set of all network nodes; PI is the set of ordinary pipeline segments (excluding compressor stations and regulating valves); P i Is the pressure value of node i; P i lb Is the lower limit (parameter) of the pressure value of node i; P i ub Is the upper limit (parameter) of the pressure value of node i;

[0182] Equation (17) ensures that the pressures of all nodes satisfy the upper and lower limit restrictions;

[0183] Equation (18) is as follows:

[0184]

[0185] In Equation (18), f ij Is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 cubic meters; VP is the set of pipeline segments containing regulating valves; PI is the set of ordinary pipeline segments (excluding compressor stations and regulating valves); Is the upper limit (parameter) of the transportation volume of the pipeline segment from node i to node j;

[0186] Equation (18) ensures that the pipeline transportation capacity of each pipeline segment does not exceed the upper limit;

[0187] Equation (19) is as follows:

[0188]

[0189] In Equation (19), the v ijm Represents the satisfaction amount d of the jth demand of node i ij Whether it is in the mth ladder. If it is in the mth ladder, take the value of 1; if it is not in the mth ladder, take the value of 0, dimensionless;

[0190] Equation (20) is as follows:

[0191]

[0192] In Equation (20), the u ijm Is a non - negative real - valued variable, representing d ijThe quantity falling onto the m-th step, in units of 10,000 cubic meters;

[0193] Equation (22) is as follows:

[0194]

[0195] In Equation (22), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, independent of direction, in units of 10,000 cubic meters; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable; f ij,0 is the flow value at the base point, in units of 10,000 cubic meters; P i is the pressure value of node i, in units of Pascal; P j The pressure value of node j, in units of Pascal; Among them, λ is the friction coefficient between the gas and the inner wall of the pipe, dimensionless; Z is the gas compressibility factor, dimensionless; Δ * is the relative density of the gas, dimensionless; T is the average temperature of the pipe, in Kelvin; L is the length of the pipe, in meters; C0 takes the constant 0.03848; D is the inner diameter of the pipe, in meters; Among them, the parameter Among them, g is the acceleration due to gravity; Z is the gas compressibility factor, dimensionless; T is the average temperature of the pipe, in Kelvin; R is the gas constant; h i is the altitude of the i-th section of the pipe, in meters; L i is the length of the i-th section of the pipe, in meters; β = a·Δh, where a is the same as a in calculating θ, and Δh is the difference in altitude between the starting point and the ending point of the pipe section, in meters;

[0196] Equation (22) describes the physical relationship between air pressure and air flow on each ordinary transmission pipeline (pipeline without compressor stations and regulating valves). x ij is a 0 / 1 variable indicating whether the transmission direction is from i to j. When its value is 0 (i.e., the transmission direction is from j to i), the value of (2x ij -1) on the left side of the constraint is -1, which serves to reverse the physical relationship formula.

[0197] Equation (22.1) is as follows:

[0198]

[0199] The parameter meanings in Equation (22.1) are exactly the same as those in Equation (22);

[0200] Equation (23) is as follows:

[0201]

[0202] In formula (23), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, independent of direction, with the unit of 10,000 cubic meters; VP is the set of pipe segments with regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0, which is the air flow direction variable; f ij,0 is the flow value at the reference point, with the unit of 10,000 cubic meters; P i is the pressure value of node i, with the unit of Pascal; P j is the pressure value of node j, with the unit of Pascal; γ, θ, and β in formula (23) are the same as those in formula (22).

[0203] Formula (23) represents the physical relationship between air pressure and air flow in the transmission pipeline with regulating valves;

[0204] Formula (24) is as follows:

[0205]

[0206] In formula (24), VP is the set of pipe segments with regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0, which is the air flow direction variable; f ij,0 is the flow value at the reference point, with the unit of 10,000 cubic meters; P i is the pressure value of node i, with the unit of Pascal; P j is the pressure value of node j, with the unit of Pascal; γ, θ, and β in formula (24) are the same as those in formula (22);

[0207] Formula (24) represents the physical relationship between air pressure and air flow in the transmission pipeline with regulating valves;

[0208] Formula (25) is as follows:

[0209]

[0210] In formula (25), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, independent of direction, with the unit of 10,000 cubic meters; VP is the set of pipe segments with regulating valves; f ij,0 is the flow value at the reference point, with the unit of 10,000 cubic meters; ε is set artificially according to the convergence situation during the calculation, representing the difference between the unsolved solution and the initial solution; f ij,0 is the flow value at the reference point, with the unit of 10,000 cubic meters;

[0211] Equation (25) states that the error between the final solution and the initial solution should be less than a certain range, which is a limitation on the solution result, ensuring that the error of the linear approximation is within a reasonable range and controlling that the flow value does not exceed a certain range near the base point;

[0212] Exemplarily, in step 3), the "close" means that the solution result is less than 10 compared to the optimal value in the previous iteration -4 .

[0213] Exemplarily, in step 3), the "sufficiently small" means that the solution result is less than 10 -4 .

[0214] Exemplarily, the value range of ε is from 0.005 to 0.0005, and it can also be appropriately adjusted according to the conventional calculation accuracy requirements in the art.

[0215] Exemplarily, the nodes in the present disclosure are the stations in the actual production process.

[0216] Exemplarily, the base point in the present disclosure represents the initial starting point when the model is solved. The domain radius in the present disclosure represents the amplitude by which the variable can move near the base point. The non - linear mixing in the semi - disclosure is a non - linear mixed - integer programming problem (NLMIP), which means that a mathematical formula with a quadratic or higher - order term is added to the objective function or constraints of a general linear programming model. The natural gas demand points in the present disclosure indicate whether the set of natural gas demand points is the set of demand nodes I in the formula. The natural gas supply points in the present disclosure represent the set O of gas source nodes. The sales ladder in the present disclosure, for example, the m - th ladder is the m - th sales ladder, and so on.

[0217] Exemplarily, in step 1), reading the natural gas pipeline network data and supply - demand price data may include reading:

[0218] 1. Supply category:

[0219] 1.1. Domestic and foreign purchased gas: Upper limit of supply volume (10,000 m³ / day), lower limit of supply volume (10,000 m³ / day), purchase price (yuan / 10,000 m³).

[0220] 1.2. Domestically produced gas: Upper limit of supply volume (10,000 m³ / day), lower limit of supply volume (10,000 m³ / day), purchase price (yuan / 10,000 m³).

[0221] 1.3. Imported gas: Upper limit of supply volume (10,000 m³ / day), lower limit of supply volume (10,000 m³ / day), purchase price (yuan / 10,000 m³). (Imported gas is purchased from abroad, and domestic and foreign purchased gas is purchased from other domestic producers)

[0222] 2. Pipeline network category:

[0223] 2.1. Station: Station name, compressor information, upper pressure limit, lower pressure limit, altitude.

[0224] 2.2, Pipe section: Revised outer diameter, absolute roughness, gas transmission efficiency, ambient temperature, overall heat transfer coefficient, wall thickness, resistance coefficient, mileage, pipe transmission capacity (10,000 m³ / day), maximum allowable pressure, unit operating cost (yuan / (10,000 m³·km)), freight rate (yuan / (10,000 m³·km)).

[0225] 2.3, Compressor: Compression ratio, maximum downstream pressure, minimum upstream pressure, minimum downstream temperature, maximum downstream temperature, maximum power.

[0226] 3. Requirements:

[0227] 3.1, Customer: Minimum sales volume (10,000 m³ / day), maximum sales volume (10,000 m³ / day), sales price (yuan / 10,000 m³).

[0228] 4. Others:

[0229] 4.1, Gas storage reservoir: Connected station code, connected station service code, minimum gas production volume (10,000 m³ / day), maximum gas production volume (10,000 m³ / day), minimum gas injection volume (10,000 m³ / day), maximum gas injection volume (10,000 m³ / day), unit gas injection cost (yuan / 10,000 m³), unit gas production cost (yuan / 10,000 m³), capacity cost (yuan / month).

[0230] All the parameters used in the method provided by the present disclosure can be obtained based on the above natural gas pipeline network data and supply-demand price data.

[0231] In the mathematical model involved in the present disclosure, the non-linear constraints of pressure and flow in the pipeline and integer variables such as the direction compressor switch are simultaneously considered, and for integer-type variables, they are solved together with continuous variables. Therefore, theoretically, compared with some practices of distributed solution, such as: first ignoring the non-linear physical relationship of the pipeline, only considering the supply-demand constraints to solve the integer-type variables and fixing them, and then adding non-linear constraints to solve; a better objective value can be obtained.

[0232] Secondly, the linear approximation method proposed in the present disclosure can obtain a high solution accuracy. As the expansion radius gradually decreases, the moving range of the variable near the base point becomes smaller, and the error between the linear approximation and the non-linear function at the base point is zero. Therefore, as the expansion radius approaches 0, the error between the linear approximation and the non-linear function also continuously approaches 0.

[0233] Taking the production and operation plan calculation of a certain month as an example, in the scenario where the sales of the long-distance natural gas pipeline are 293 million cubic meters per day under the same calculation, when using the linear optimization model, since only the pipeline transportation capacity is considered and factors such as pipeline transportation pressure are not considered, the calculated pipeline transportation fee is 4.38 billion yuan. After using the non-linear optimization model, due to considering the relationship between pipeline transportation and flow rate, the capacity of some pipeline sections cannot reach the corresponding flow rate due to pressure factors, and it will consider transporting through the corresponding longer pipelines, and the pipeline transportation fee is 4.708 billion yuan. Moreover, it can calculate the pressure of the corresponding stations.

Claims

1. A linear approximation method for optimizing the scheduling of natural gas pipelines considering hydraulic constraints, characterized in that, The method includes the following steps: Obtain the parameters required for solving the non - linear mixed - integer programming problem of natural gas pipeline transportation, and substitute them into the non - linear mixed - integer programming problem of natural gas pipeline transportation for solution; The non - linear mixed - integer programming problem contains discontinuous factors and steady - state hydraulic constraints; The discontinuous factors at least include: steady - state hydraulic constraints, step - supply price function, pipeline direction; The steady - state hydraulic constraints at least include: non - linear terms regarding the pressure variables at both ends of the pipeline and non - linear terms of the flow variables; Solving the non - linear mixed - integer programming problem includes: making a linear approximation of the non - linearity in the hydraulic constraints near the base point, and degrading the non - linear mixed - integer programming problem into a linear mixed - integer programming problem.

2. The method according to claim 1, characterized in that, The method includes the following steps: Set the original problem as maximizing the sales profit = total sales amount at natural gas demand points - total purchase amount at natural gas supply points - transportation cost, that is, formula (1); 1) Read the natural gas pipeline network data and supply - demand price data; 2) Take all flow values as 0 as the initial base point, set the expansion radius ε as 1, and substitute it into the objective function of stage one: The constraint formulas that the objective function of stage one obeys are as follows: formula (2), formula (3), formula (4), formula (5), formula (8), formula (9), formula (6), formula (13), formula (14), formula (15), formula (16), formula (17), formula (18), formula (19), formula (20), formula (22.1), formula (23), formula (24), formula (25); In the objective function of Phase 1, represents the relaxed positive terms, represents the relaxed negative terms; solve the problem using a linear mixed-integer planner and record the optimal value of the problem; 3) Judge whether the optimal value of the objective function of stage one solved in step 2) is close enough to the optimal value in the previous iteration when solving the objective function of stage one. If it is close enough, reduce the expansion radius ε to 0.5 times that in the previous iteration. Judge whether the expansion radius is small enough. If it is not small enough, update the base point to the pipeline flow value in the current optimal solution, and substitute it into step 3) for re - solution; 4) If the expansion radius in step 3) is small enough, judge whether it is in stage one. If it is not in stage one, judge whether it is in stage two currently. If it is not in stage two, then complete the solution of the original problem, and the parameters at this time can maximize the sales profit; 5) If the conclusion of determining whether it is in Phase 1 in step 4) is yes, then determine whether the optimal value is 0. If not, the original problem has no solution. If so, substitute it into the objective function of Phase 2: Update the base point to the pipeline flow value in the current optimal solution; The constraint formulas that the objective function of stage two obeys are as follows: formula (2), formula (3), formula (4), formula (5.1), formula (8), formula (9), formula (6), formula (13), formula (14), formula (15), formula (16), formula (17), formula (18), formula (19), formula (20), formula (22), formula (23), formula (24), formula (25); Update the base point to the pipeline flow value in the current optimal solution, and substitute it into steps 2) to 4) for solution; 6) If the conclusion in step 4) that it is in stage two is yes, then determine whether the optimal value is 0. If not, then the original problem has no solution; if so, substitute it into the objective function of stage three: in formula (1), update the base point to the pipeline flow value in the current optimal solution, and substitute it into step 1) for solution; The constraint formulas that the formula (1) obeys are as follows: formula (2), formula (3), formula (4), formula (5), formula (8), formula (9), formula (6), formula (13), formula (14), formula (15), formula (16), formula (17), formula (18), formula (19), formula (20), formula (22), formula (23), formula (24), formula (25); Stage one and stage two are the solution processes of the initial feasible solution of stage three. Their association is through the system of constraint equations. In stage one, the pressure-flow constraint is minimized, and in stage two, the customer relaxation is minimized. After these two conditions are basically met, the optimal solution of the original problem is then obtained on this basis; Formula (1) is as follows: In formula (1), I is the set of demand nodes; O is the set of gas source nodes, PI is the set of ordinary pipeline segments; CP is the set of pipeline segments with compressor stations; R ijm is the total sales amount in the first m-1 stages when the usage reaches the mth ladder, in ten thousand yuan; p ijm is the unit price of the jth demand of node i at the mth ladder, in yuan per ten thousand cubic meters; is the supply unit price of the jth gas source of node i, in yuan per ten thousand cubic meters; is the transportation unit price from node i to node j, in yuan per ten thousand cubic meters; Formula (2) is as follows: In formula (2), the v ijm represents the satisfaction amount d of the j-th requirement of node i ij whether it is in the m-th ladder. If it is in the m-th ladder, take a value of 1; if it is not in the m-th ladder, take a value of 0, dimensionless Formula (3) is as follows: In formula (3), q ijm is the quantity of the m-th ladder of the j-th demand of node i, with the unit of 10,000 cubic meters; Formula (4) is as follows: In formula (4), Q ijm is the total sales volume in the first m - 1 stages when the usage reaches the mth step, with the unit of 10,000 cubic meters; d ij is the satisfied quantity of the jth demand at node i, with the unit of 10,000 cubic meters; Formula (5) is as follows: In formula (5), is the lower limit of the j-th demand of node i, with the unit of 10,000 cubic meters; is the upper limit of the j-th demand of node i, with the unit of 10,000 cubic meters; Formula (5.1) is as follows: In formula (5.1), is the lower limit of the j-th demand of node i, in ten thousand cubic meters; is the upper limit of the j-th demand of node i, in ten thousand cubic meters; ξ ij is the slack term of the demand, in ten thousand cubic meters; Formula (6) is as follows: In formula (6), o i is the total procurement volume of node i, with the unit of 10,000 cubic meters; Formula (7) is as follows: In formula (7), d i is the total demand of node i, with the unit of 10,000 cubic meters; Formula (8) is as follows: In formula (8), o ij is the procurement volume of the j-th gas source at node i, with the unit of 10,000 m³; is the lower supply limit of the j-th gas source at node i, with the unit of 10,000 m³; is the upper supply limit of the j-th gas source at node i, with the unit of 10,000 m³; Formula (9) is as follows: In formula (9), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of the direction, with the unit of 10,000 m³; N is the set of all network nodes; VP is the set of pipe segments with regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable; Formula (13) is as follows: P in Formula (22), Formula (22.1), Formula (23), and Formula (24) i and P j satisfy Formula (13): In formula (13), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe segments containing gas compressor stations; x ij takes a value of 0 or 1. If the gas flow is from i to j, the value is 1; otherwise, it is 0, which is the gas flow direction variable. Formula (14) is as follows: P in Formula (22), Formula (22.1), Formula (23) and Formula (24) i and P j also satisfy Formula (14): In formula (14), P i is the pressure value of node i; P j is the pressure value of node j; CP is the set of pipe segments with compressor stations; x ij takes a value of 0 or 1. If the gas flow is from i to j, the value is 1; otherwise, it is 0, which is the gas flow direction variable; Formula (15) is as follows: P in Formula (22), Formula (22.1), Formula (23) and Formula (24) i and P j also satisfy Formula (15): In formula (15), P i is the pressure value of node i; P j is the pressure value of node i; P i ub is the upper limit of the pressure value of node i; y ij takes a value of 0 or 1. The bypass value of the compressor station takes a value of 1, otherwise 0; x ij takes a value of 0 or 1. If the gas flow is from i to j, it takes a value of 1, otherwise 0, which is the gas flow direction variable; Formula (16) is as follows: P in Formula (22), Formula (22.1), Formula (23), and Formula (24) i and P j also satisfy Formula (16): In formula (16), P i is the pressure value of node i; P j is the pressure value of node j; is the lower limit of the compression ratio of compressor station i; y ij takes a value of 0 or 1. When bypassing the compressor station, the value is 1; otherwise, it is 0. Equation (17) is as follows: P in Equation (22), Equation (22.1), Equation (23), and Equation (24) i satisfies Equation (17): In formula (17), N is the set of all network nodes; PI is the set of ordinary pipe segments (excluding compressor stations and regulating valves); P i is the pressure value of node i; P i lb is the lower limit of the pressure value of node i; P i ub is the upper limit of the pressure value of node i; The formula (18) is as follows: In formula (18), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 m³; VP is the set of pipe segments with regulating valves; PI is the set of ordinary pipe segments; is the upper limit of the transportation volume of the pipe segment from node i to node j; Formula (19) is as follows: In formula (19), the v ijm represents the satisfaction quantity d of the j-th requirement of node i ij indicating whether it is in the m-th ladder. If it is in the m-th ladder, it takes a value of 1; if it is not in the m-th ladder, it takes a value of 0, dimensionless Formula (20) is as follows: The said u ijm is a non - negative real - valued variable, representing the quantity that falls into the m - th step in d ij with the unit of 10,000 cubic meters; Formula (22) is as follows: In formula (22), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, independent of direction, with the unit of 10,000 m³; x ij takes the value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0. It is the air flow direction variable; f ij,0 is the flow value at the reference point, with the unit of 10,000 m³; P i is the pressure value of node i, with the unit of Pa; P j is the pressure value of node j, with the unit of Pa; where, λ is the friction coefficient between the gas and the inner wall of the pipe, dimensionless; Z is the gas compressibility factor, dimensionless; Δ * is the relative density of the gas, dimensionless; T is the average temperature of the pipe, in Kelvin; L is the length of the pipe, in meters; C0 takes the constant value of 0.03848; D is the inner diameter of the pipe, in meters m; where, the parameter where, g is the acceleration due to gravity; Z is the gas compressibility factor, dimensionless; T is the average temperature of the pipe, in Kelvin; R is the gas constant; h i is the elevation of the i-th section of the pipe, in meters; L i is the length of the i-th section of the pipe, in meters; β = a·Δh, where, a is the same as a in calculating θ, and Δh is the difference in elevation between the starting point and the ending point of the pipe section, in meters; Formula (22.1) is as follows: The parameter meanings in formula (22.1) are exactly the same as those in formula (22); Formula (23) is as follows: In formula (23), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 cubic meters; VP is the set of pipe segments with regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1, otherwise it is 0, which is the air flow direction variable; f ij,0 is the flow value at the base point, with the unit of 10,000 cubic meters; P i is the pressure value of node i, with the unit of Pascal; P j is the pressure value of node j, with the unit of Pascal; γ, θ, and β in formula (23) are the same as those in formula (22); Formula (24) is as follows: In formula (24), VP is the set of pipe segments containing regulating valves; x ij takes a value of 0 or 1. If the air flow is from i to j, the value is 1; otherwise, it is 0, which is the air flow direction variable; f ij,0 is the flow value at the reference point, with the unit of 10,000 m³; P i is the pressure value of node i, with the unit of Pa; P j is the pressure value of node j, with the unit of Pa; γ, θ, and β in formula (24) are the same as those in formula (22); Formula (25) is as follows: In formula (1), f ij also satisfies formula (25): In formula (25), f ij is the absolute value of the transportation volume under standard conditions between node i and node j, regardless of direction, with the unit of 10,000 m³; VP is the set of pipe segments with regulating valves; f ij,0 is the flow value at the base point, with the unit of 10,000 m³; ε is set artificially according to the convergence situation during the calculation, representing the difference between the unsolved solution and the initial solution; f ij,0 is the flow value at the base point, with the unit of 10,000 m³.

3. The method according to claim 2, wherein, In step 3), the proximity means that the solution result is less than 10 from the optimal value in the previous iteration -4 .

4. The method according to claim 2, wherein, In step 3), the "sufficiently small" means that the solution result is less than 10 -4 .

5. The method according to claim 2, wherein, The value range of the said ε is from 0.005 to 0.0005.