Optimal scheduling method of multi-source reservoir-canal system in irrigation district based on effective water level of canal tail
By embedding dynamic constraints and micro-simulation of the effective water level at the end of the canal into the water supply scheduling of the irrigation area, the water supply process of the canal system is optimized, the problems of ineffective water supply and tailwater loss are solved, and efficient and reliable water resource management is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- NANJING HYDRAULIC RES INST
- Filing Date
- 2026-05-14
- Publication Date
- 2026-06-12
Smart Images
Figure CN122198668A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the interdisciplinary technical fields of water resource optimization allocation, hydraulic digital twin and fluid dynamics simulation, and in particular relates to an optimization scheduling method for multi-source reservoir-canal system in irrigation districts based on the effective gravity flow water level at the end of the canal. Background Technology
[0002] In modern irrigation district operation and management, reservoirs and long-distance water conveyance canal systems have become core infrastructure for ensuring agricultural water security. With the increasing scarcity of water resources, irrigation district networks are gradually developing into multi-source reservoir-canal joint allocation systems, and the construction of smart water conservancy and digital twins is being widely promoted. In actual operation, irrigation district water supply heavily relies on gravity flow at the end of the canal system; whether the water level at the end of the canal can reach the effective conditions for gravity flow directly determines the feasibility of water supply.
[0003] To address the aforementioned issues, existing technologies typically employ macroscopic water balance methods to establish optimal scheduling models. This involves constructing system water balance equations based on the water supply process of reservoir groups, canal water conveyance losses, and user water demand processes. Macroscopic indicators such as total water supply and water shortage are used as objective functions, and linear programming and genetic algorithms are employed to solve for the optimal outflow process. This type of method generalizes canals as "pipeline" components with certain water conveyance losses, uses flow rate as the core decision variable, and calculates the degree of supply-demand matching through time integration.
[0004] However, the aforementioned existing technologies have obvious drawbacks: First, using the cumulative flow value as the assessment criterion makes it difficult to reflect the dynamic changes in the water level at the end of the canal. The optimized results show that the water volume is sufficient, but the actual water level is lower than the effective gravity flow water level, resulting in ineffective water supply. Second, the non-steady flow drainage process in the canal is ignored. After the gate is closed, the stagnant tailwater loses its gravity flow capacity as the water level drops, resulting in unusable tailwater loss. Third, the macroscopic model cannot couple the hydraulic response characteristics of the canal system. The optimized outflow process is difficult to execute in the actual system, and the scheduling scheme has poor feasibility. Summary of the Invention
[0005] Purpose of the Invention: The purpose of this invention is to provide an optimized scheduling method for multi-source reservoir-canal systems in irrigation districts based on the effective gravity-flow water level at the end of the canal. This method aims to solve the problems of ineffective water supply, tailwater loss, and poor feasibility of existing macro-water balance scheduling due to neglecting the dynamic changes of the effective gravity-flow water level at the end of the canal and the non-steady flow drainage process. By coupling the macro-scheduling of the reservoir group with the micro-simulation of the canal system, the dynamic gravity-flow water level at the end of the canal is dynamically embedded into the model for optimization, thereby achieving efficient utilization of water resources.
[0006] Technical solution: The optimized scheduling method for multi-source reservoir-canal system in irrigation districts based on the effective gravity-flow water level at the end of the canal, as described in this invention, includes the following steps:
[0007] S1. Obtain the topology and engineering parameters of the irrigation district reservoir group and backbone canal system. The backbone canal system includes several backbone channels. Combine remote sensing images, UAV images and elevation data to determine the minimum effective gravity flow water level at the end of each backbone channel to achieve gravity flow water discharge. Then, deduce the minimum water supply flow lower limit and the pre-storage flow threshold before formal water supply for each backbone channel under the condition of maintaining its effective gravity flow water level at the end, and form the scheduling bottom line constraint boundary.
[0008] S2. Divide the entire scheduling cycle into the channel pre-storage period and the reservoir water supply period, construct the initial channel storage capacity, and establish the channel end gravity flow discrimination function to determine whether the output water volume at the end of the backbone channel constitutes effective water supply, so as to obtain the initial water supply target required for macro-optimization, and at the same time give the target water demand of each backbone channel in the scheduling cycle.
[0009] S3. Using reservoir discharge flow, canal headwater flow, and gate opening / closing status as decision variables, establish a mixed integer linear programming model that includes minimizing discharge, minimizing canal distribution time, and minimizing gate opening time. Solve the model under the constraints of reservoir water balance, system supply and demand balance, canal flow balance, reservoir and canal boundary constraints, and water level and reservoir capacity variation constraints. Output the reservoir discharge process and the main canal headwater flow process for each time period.
[0010] S4. The discharge process of the reservoir group and the flow process at the head of the canal are used as the upstream boundary conditions of the one-dimensional unsteady flow dynamics model. The Saint-Venant equations with the control gate opening and closing are constructed. The boundary conditions are switched between the pre-storage period and the water supply period according to whether the water depth at the end of the canal reaches the minimum effective gravity flow level. The implicit difference scheme is used to solve the problem to obtain the actual water depth at the end of each main canal and the output water volume at the end. The output water volume at the end is integrated by the gravity flow discrimination function at the end of the canal to calculate the net effective water supply of each main canal.
[0011] S5. Calculate the deviation between the net effective water supply of each backbone channel and the target water demand, and use whether the maximum deviation among all backbone channels is less than the preset engineering tolerance threshold as the convergence criterion; when the convergence condition is not met and the maximum number of iterations is not reached, the initial water supply target is corrected according to the damping learning rate, and the process returns to step S3 to re-optimize until the convergence condition is met or the maximum number of iterations is reached.
[0012] The proposed optimization scheduling method for multi-source reservoir-canal systems in irrigation districts based on the effective gravity-flow water level at the end of the canal is as follows: Step S1 converts the minimum effective gravity-flow water level at the end of the canal, derived from remote sensing and elevation data, into a lower limit for water supply flow and a pre-storage flow threshold, thus constructing a bottom-line constraint for scheduling from the source. This effectively alleviates the problem of tailwater not reaching the ground and ineffective water supply caused by neglecting gravity-flow conditions in traditional macro-scheduling. Step S2, by dividing the pre-storage period and the water supply period and establishing a gravity-flow discrimination function at the end of the canal, links the dynamic water depth at the end with the water supply target, mitigating the problem of false water volume reporting caused by unsteady flow recession. Step S3 constructs a decision-making mechanism based on reservoir discharge, headwater diversion, and gate opening and closing. A mixed-integer linear programming model with multiple variables outputs the macroscopic scheduling process under various equilibrium constraints, significantly shortening the time for water distribution and gate operation, and improving system operating efficiency. Step S4 utilizes a one-dimensional unsteady flow dynamics model and the Saint-Venant equations, based on the boundary conditions of the canal end water depth switching between the pre-storage period and the water supply period, to accurately simulate the actual hydraulic response and net effective water supply of the canal system, overcoming the shortcomings of traditional water balance methods in characterizing dynamic changes in water level and drainage losses. Step S5 iteratively corrects the initial water supply target through deviation feedback and damping learning rate until the supply level difference of all channels meets the engineering tolerance, realizing the closed-loop coupling of macroscopic scheduling and microscopic hydraulic simulation. Overall, this method deeply embeds the dynamic changes of effective water level at the end of the canal and the unsteady flow process into the optimization framework, significantly suppressing ineffective water supply and tailrace losses, and improving the feasibility of multi-source reservoir-canal system scheduling schemes and the level of efficient water resource utilization.
[0013] Preferably, the minimum water supply flow limit of each backbone channel under the condition of maintaining its effective gravity-flow water level at the end, as described in step S1, is calculated using the Manning formula optimized with extreme value compensation:
[0014]
[0015] in, The minimum water supply flow rate for the main channel j is expressed in cubic meters per second. Let J be the Manning roughness of the backbone channel j, which is dimensionless; The width of the bottom of the main channel j is in meters. The minimum effective gravity-flow water level of the main channel j, in meters; The bottom slope of the main channel j is dimensionless.
[0016] This optimization step, by introducing an extreme value compensation coefficient for drainage loss, optimizes the Manning formula, enabling more accurate calculation of the minimum water supply flow required for each main channel to maintain an effective gravity-flow water level at the end. This effectively compensates for hydraulic losses during actual drainage and enhances the engineering applicability and reliability of the scheduling bottom line constraint.
[0017] Preferably, the formula for calculating the pre-stored flow threshold in step S1 is:
[0018]
[0019] in, The minimum flow threshold required for pre-storage of the main water supply channel before formal water supply, expressed in cubic meters per second; , is the pre-storage ratio coefficient for the main channel j, dimensionless, and its value range is determined according to the channel's scour resistance requirements and the need for stable water lifting.
[0020] This preferred step defines the pre-storage flow threshold before formal water supply by multiplying the pre-storage ratio coefficient by the minimum water supply flow lower limit. Under the premise of ensuring the safety of the channel against scour and the smooth water lifting, the initial storage capacity of the canal system can be established in advance, which can effectively alleviate the drastic fluctuations in water level and the lag in the end outlet caused by the sudden increase in flow in the early stage of water supply, and create conditions for the rapid formation of the effective gravity flow water level at the end of the canal.
[0021] Preferably, the initial storage capacity of the canal system in step S2 is determined by the geometric parameters of each main canal and the minimum effective gravity-flow water level. The initial water supply target is expressed as follows in the first iteration:
[0022]
[0023] in, The initial macroscopic water supply target for the first iteration is expressed in cubic meters. The target water demand for the irrigation area corresponding to the main channel j within the scheduling cycle is expressed in cubic meters. The width of the bottom of the main channel j is in meters. Here is the length of the main channel j, in meters; This refers to the total number of core channels; The first term represents the minimum effective gravity-flow water level of the main channel j, in meters; the second term indicates that the initial storage capacity of the canal system is included as part of the water supply compensation in the macro water supply target.
[0024] This optimization step incorporates the initial storage capacity of the canal system, determined by the geometric parameters of the backbone channels and the minimum effective gravity flow water level, as a water conveyance compensation term into the macroscopic initial water supply target of the first iteration. This achieves quantitative reservation of the canal system's water storage demand from the source, alleviating the problem of insufficient end-point water level or low water supply target caused by ignoring the storage process in traditional scheduling. It also provides an initial benchmark that is more consistent with the hydraulic physical process for the closed-loop iteration between subsequent macroscopic optimization and microscopic simulation.
[0025] Preferably, the formula for calculating the channel end gravity flow discrimination function in step S2 is:
[0026]
[0027] in, is the gravity flow discrimination index of the main channel j at time t. A value of 1 indicates that the water output at the end of the channel is an effective water supply, and a value of 0 indicates an ineffective water supply. The water depth at point L, the end of the main channel j, at time t is expressed in meters. The minimum effective gravity-flow water level of the main channel j is expressed in meters.
[0028] This optimization step constructs a binary gravity flow discrimination index based on the difference between the actual water depth at the end of the canal and the lowest effective gravity flow level. This enables real-time and quantitative determination of whether the water output at the end of the main canal constitutes effective water supply. In this way, it can accurately distinguish between effective water supply and ineffective water discharge in macro-optimization, alleviating the problem of miscounting non-gravity flow discharge as water supply in traditional water balance scheduling due to the inability to identify the end water level conditions. This lays the basis for the accurate calculation of net effective water supply.
[0029] Preferably, the objective function of the mixed-integer linear programming model described in step S3 is:
[0030]
[0031] in, Let be the discharge flow rate of reservoir i at time t, expressed in cubic meters per second. The headwater flow rate of the main channel j at time t is expressed in cubic meters per second. Let be the gate opening / closing state variable of the backbone channel j at time t; , , This is a weight vector; the first term minimizes the total system discharge; the second term introduces a time-series penalty factor t, which is related to the weights. The co-driving optimizer allocates water rights in the early stages of scheduling to minimize the channel water distribution time; the third term is based on the gate opening state variable. Apply penalty weights This drives the concentrated release of flow to reduce the total time the gate is open.
[0032] This optimization process constructs a three-term objective function with reservoir discharge flow, headworks diversion flow, and gate opening / closing status as decision variables. The first term minimizes the total system discharge to conserve water resources. The second term introduces a time-series penalty factor to drive water rights allocation to earlier stages to shorten the water distribution time. The third term reduces the frequency of operations by penalizing the gate opening status variable. The synergistic optimization of these three terms can significantly improve scheduling efficiency and reduce water transfer losses and operating costs while ensuring water supply demand.
[0033] Preferably, the constraints of the mixed-integer linear programming model in step S3 include at least one of the following:
[0034] Reservoir water balance constraints: ,in, Let be the reservoir capacity of reservoir i at time t, in cubic meters; Let be the reservoir capacity of reservoir i at time t-1, in cubic meters; The natural inflow rate of reservoir i at time t is expressed in cubic meters per second. Let be the discharge flow rate of reservoir i at time t, expressed in cubic meters per second. The discharge flow rate is expressed in cubic meters per second; Δt is the duration of each time period, expressed in seconds.
[0035] System water supply and demand constraints: That is, the total discharged water volume shall not be less than the sum of the target water requirements of each irrigation area; among which, The target water demand for the irrigation area corresponding to the main channel j within the scheduling cycle is expressed in cubic meters.
[0036] The water distribution volume in the channel meets the constraints: ,in, For channel j at any time The water flow rate is expressed in cubic meters per second. The water demand for irrigation sections is controlled by channel j, in cubic meters;
[0037] Reservoir discharge flow variation constraints: ,in, Let be the discharge flow rate of reservoir i at time t-1, in cubic meters per second; The maximum allowable flow variation for reservoir i is expressed in cubic meters per second.
[0038] Reservoir water level fluctuation constraints: ,in, Let be the water level of reservoir i at time t, in meters. Let be the water level of reservoir i during time period t-1, in meters. The maximum permissible water level fluctuation for reservoir i, in meters;
[0039] Reservoir water level constraints: ,in, , These are the lowest and highest operating water levels of reservoir i, respectively.
[0040] Reservoir capacity constraints: ,in, Let be the reservoir capacity of reservoir i at time t, in cubic meters; , These are the dead storage capacity and the beneficial storage capacity of reservoir i, respectively.
[0041] Channel traffic constraints: ,in, The headwater flow rate of the main channel j at time t is expressed in cubic meters per second. and These are the maximum and minimum design flow rates of the backbone channel j, respectively.
[0042] Channel water level constraints: ,in, The water level of the main channel j during time period t, in meters. and These are the maximum and minimum design water levels of the main channel j, respectively.
[0043] This optimization process establishes a multi-dimensional constraint system, including reservoir water balance, system supply and demand balance, channel flow balance, flow variation of reservoirs and channels, water level variation, and reservoir capacity boundaries. It fully embeds actual physical conditions such as water flow evolution time delay, engineering safety limits, and design capabilities into the optimization model, effectively ensuring the safety and feasibility of the macro-dispatch scheme in engineering operation and mitigating the risk of unexecutable discharge processes or channel overload caused by ignoring time delay effects or variation limits.
[0044] Preferably, the one-dimensional unsteady flow dynamics model described in step S4 is described by the Saint-Venant equations:
[0045]
[0046]
[0047] Where A is the cross-sectional area of the water passage, in square meters; Q is the flow rate, in cubic meters per second; Z is the water level, in meters; x is the spatial step size, in meters; t is the time, in seconds; and g is the gravitational acceleration, in meters per square second. For frictional slopes, dimensionless;
[0048] The Saint-Venant equations are discretized and solved using the Preissmann implicit difference scheme, and state variables are introduced to simulate the opening and closing behavior of the channel control gate: when the water depth at the end of the channel... When the time is determined to be the pre-storage period, the lower boundary condition is set to gate closure, causing upstream water to be converted into dynamic reservoir capacity to fill the channel; when When the water supply period begins, the gate opening boundary conditions are triggered to release the effective water volume; among which, The water depth at point L, the end of the main channel j, at time t is expressed in meters. The minimum effective gravity-flow water level of the main channel j is expressed in meters.
[0049] This optimization step employs the one-dimensional unsteady flow Saint-Venant equations and the Preissmann implicit difference scheme to accurately simulate the hydraulic response process of the irrigation canal system. By introducing a state variable based on whether the water depth at the end of the canal reaches the minimum effective gravity flow level, the automatic switching of boundary conditions between the pre-storage period and the water supply period is achieved: during the pre-storage period, the gate closure converts upstream water into dynamic reservoir capacity to fill the canal channel, and during the water supply period, the gate opening is triggered to release effective water volume. This dynamically depicts the complete hydraulic evolution process of the canal from water storage to gravity flow, effectively overcoming the limitation of the steady flow assumption in simulating water recession and water level changes, and providing a high-fidelity physical simulation basis for the accurate calculation of net effective water supply.
[0050] Preferably, the formula for calculating the net effective water supply of each backbone channel in step S4 is as follows:
[0051]
[0052] in, This represents the net effective water supply of the backbone channel j in the k-th iteration, in cubic meters. The output flow rate at the end of the backbone channel j at time t, expressed in cubic meters per second; The total duration of the scheduling cycle is in seconds; The gravity flow discrimination index of the main channel j at time t is used to isolate and integrate the effective water volume.
[0053] This preferred step achieves quantitative calculation of net effective water supply by integrating the product of the output flow at the end of the backbone channel and the gravity flow discrimination index within the scheduling cycle. The gravity flow discrimination index, as a switching function, contributes an integral value only when the water depth at the end of the channel meets the minimum gravity flow effective water level. This effectively isolates the non-gravity flow ineffective water volume during the pre-storage period and the water discharge process from a physical mechanism, providing real, reliable, and effective water supply feedback for the subsequent deviation correction of macro water supply targets.
[0054] Preferably, step S5, which involves feedback correction of the initial water supply target according to the damping learning rate, includes:
[0055] Calculate the water supply deviation of each main channel in the current k-th iteration:
[0056]
[0057] in, The target water demand for the irrigation area corresponding to the main channel j within the scheduling cycle is expressed in cubic meters. Let be the net effective water supply of backbone channel j in the k-th iteration, in cubic meters; calculate the maximum absolute deviation among all backbone channels:
[0058]
[0059] in, This refers to the total number of core channels;
[0060] like If the current iteration number k is greater than or equal to the preset engineering tolerance threshold ε, and the current iteration number k is less than the maximum iteration number K, then the macroscopic water supply target for the next iteration is updated using the following damping error feedback formula:
[0061]
[0062] in, Let α be the macroscopic water supply target for the k-th iteration, in cubic meters; α be the system damping learning rate, dimensionless, ranging from (0,1]; and the updated... As a new supply and demand balance constraint right-hand term in the mixed-integer linear programming model in step S3, it drives iterative optimization until convergence or the maximum number of iterations is reached.
[0063] This optimization step calculates the deviation between the net effective water supply of each backbone channel and the target water demand, uses the largest absolute deviation among all channels as the convergence criterion, and introduces a damping learning rate to provide feedback correction to the macro water supply target. The micro hydraulic simulation results are effectively embedded into the macro optimization iteration process. The setting of the damping learning rate can control the correction step size, avoid target quantity oscillation during the iteration process, and ensure the stability and convergence of the optimization process. Through closed-loop feedback, the optimized scheduling scheme that meets the gravity flow water supply demand of all channels is gradually approximated, and finally a high degree of synergy between macro water allocation and micro hydraulic response is achieved.
[0064] Beneficial Effects: Compared with existing technologies, this invention has the following significant advantages: 1. By coupling macroscopic scheduling of reservoir groups with microscopic simulation of canal systems, this invention incorporates the dynamic effective water level at the end of the canal as a core constraint into the model optimization, overcoming the shortcomings of traditional macroscopic water balance scheduling that ignores dynamic changes in water level and the non-steady flow recession process. This alleviates ineffective water supply and tailwater loss, significantly improving water resource utilization efficiency and the feasibility of scheduling schemes. 2. By reverse-deriving the minimum water supply flow lower limit and pre-storage flow threshold of each backbone canal under the condition of maintaining the effective water level at the end of the canal, a bottom-line constraint boundary for scheduling is formed. Combined with the reasonable division of the canal pre-storage period and water supply period, this effectively alleviates water supply failure caused by insufficient initial storage or excessively rapid recession, ensuring stable water supply at the end of the canal system. 1. Gravity-driven water flow is achieved; 2. A mixed-integer linear programming model and a one-dimensional unsteady flow dynamics model are coupled and iteratively solved, and the end-of-channel gravity flow discriminant function is used as the feedback correction basis. This enables the reservoir discharge and headwater water distribution plan to dynamically respond to hydraulic changes within the canal system. While minimizing the system discharge volume, water distribution time, and gate operation costs, the actual output water volume at the end of the channel meets the irrigation area's water demand requirements; 3. The macro water supply target is adaptively corrected through damping learning rate, and the maximum water supply deviation in all backbone channels is less than the engineering tolerance threshold is used as the convergence criterion. This ensures that the deviation between the net effective water supply and the actual water demand in each channel is controlled within the allowable range. The iterative process is stable and reliable, and the final output scheduling scheme has high precision and strong engineering practicality. Attached Figure Description
[0065] Figure 1 This is a schematic diagram of the method flow of the present invention;
[0066] Figure 2 This is a schematic diagram of the topology of the multi-source reservoir-canal system in the irrigation area according to the present invention;
[0067] Figure 3 This is a schematic diagram of the reservoir group storage process of the present invention;
[0068] Figure 4 This is a schematic diagram of the reservoir group discharge process according to the present invention;
[0069] Figure 5 This is a schematic diagram of the canal system flow process according to the present invention;
[0070] Figure 6 This is a schematic diagram of the cumulative water volume of the canal system according to the present invention. Detailed Implementation
[0071] The technical solution of the present invention will be further described below with reference to the accompanying drawings.
[0072] This invention provides an optimized scheduling method for a multi-source reservoir-canal system in an irrigation district based on the effective gravity-flow water level at the end of the canal. Figure 1As shown, it includes the following steps:
[0073] Step S1: Construction of physical boundaries and back calculation of gravity flow in the multi-source reservoir-canal system of the irrigation district
[0074] (1) Obtain the topological structure of the reservoir group and the backbone canal system in the irrigation area, analyze the relationship between the reservoir and canal system in the irrigation area, and clarify the water source situation, irrigation area division, water distribution direction, etc. in the irrigation area. Collect engineering data of the irrigation area, including the characteristic water level, storage capacity, water storage status, and historical storage and discharge data of the reservoirs; bottom width of the backbone canals. ,length Roughness bottom slope Physical parameters and historical water distribution data; water demand of each irrigation area controlled by the canals, planting structure, etc., to form an irrigation district database.
[0075] (2) Based on satellite and UAV remote sensing images and elevation (DEM) data of the irrigation area, the elevation relationship between the irrigation area and the control section at the end of the canal is mapped, and combined with the irrigation area management experience, the minimum effective gravity flow water level required for gravity flow at the end of each canal system is comprehensively determined. .
[0076] (3) The lower limit of gravity flow and the pre-stored flow are derived by reverse deduction of Manning formula to generate the bottom line constraint boundary of the optimization model in step 3.
[0077] ①Lower limit of gravity flow Channel J maintains the minimum effective gravity-flow water level. The required lower limit of physical flow rate serves as the flow rate condition that must be met for multi-source reservoir and canal water distribution.
[0078]
[0079] ②Pre-stored flow Based on the channel's erosion resistance requirements and the need for stable water supply, a small flow threshold is set for pre-storage in the channel before formal water supply.
[0080]
[0081] Let be the pre-storage ratio coefficient for channel j.
[0082] Step S2: Construction of Supply and Demand Optimization Scheduling Objectives and Phased Strategies for the Reservoir-Canal System
[0083] (1) Divide the entire scheduling cycle T into channel pre-storage period on the time axis. and reservoir water supply period Two non-overlapping subsets satisfying .
[0084] (2) Construct the initial storage capacity of the canal system as a compensation for water conveyance that cannot achieve gravity flow in the reservoir-canal system. Set the initial macroscopic water supply target for the first iteration (k=1) as:
[0085] In the formula Let J be the net water requirement (m³) for irrigation area j, where J is the number of irrigation areas, corresponding one-to-one with the control channel.
[0086] (3) Construct the self-flow discrimination function at the end of the channel ,when When the water distribution volume of channel j meets the gravity irrigation demand, it is the effective water volume (m³); when At that time, it is considered invalid water volume.
[0087]
[0088] In the formula Let L be the water depth (m) at the end of the channel at time t.
[0089] Step S3: Macroscopic Multi-Source Reservoir-Channel System Water Distribution Flow Optimization Decision
[0090] (1) Construct a system based on the reservoir discharge flow rate Headwater distribution flow and the opening and closing status of the gate For the mixed integer linear programming (MILP) model of the decision space, a minimization objective function is constructed that includes minimizing the discharge, minimizing the channel water distribution time, and minimizing the channel gate opening time.
[0091]
[0092] In the formula Let be the discharge flow rate (m³ / s) of reservoir i at time t. Let be the headwater flow rate (m³ / s) of channel j at time t. Let j be the gate opening state at time t; , , Let be the weight vector, where: the first term minimizes system leakage; the second term introduces a timing penalty factor t, where t represents time and increases linearly with the scheduling time, thus affecting the weights. The driver optimizer automatically allocates water rights to the early stages of scheduling to enable the channel to complete the water distribution task as quickly as possible; the third step involves monitoring the gate opening status. punishment This drives the centralized release of traffic, reducing the duration of channel gate opening.
[0093] (2) Construct a multidimensional constraint set: including reservoir water balance constraints, system water supply and demand constraints, channel water distribution constraints, reservoir discharge flow variation, reservoir water level variation, reservoir water level, reservoir capacity constraints, channel flow, water level constraints, etc. as constraints of the model.
[0094] Reservoir water balance constraints:
[0095] System water supply and demand constraints:
[0096] The water distribution volume in the channel meets the constraints:
[0097] Reservoir discharge flow variation constraints:
[0098] Reservoir water level fluctuation constraints:
[0099] Reservoir water level constraints:
[0100] Reservoir capacity constraints:
[0101] Channel traffic constraints:
[0102] Channel water level constraints:
[0103] in, Let be the reservoir capacity (m³) of reservoir i at time t; Natural inflow rate (m³ / s); The discharge flow rate is (m³ / s); Δt is the duration of each time period, in seconds. For channel j at any time The water flow rate is expressed in cubic meters per second. Control the water demand (m³) for irrigation area j in channel j; Let i be the maximum allowable flow rate variation of reservoir i (m³ / s). Let be the water level (m) of reservoir i during time period t. This refers to the fluctuation range of the reservoir water level; , These are the upper and lower limits of the reservoir water level (m); Let be the reservoir capacity (m³) of reservoir i during time period t. , These are the upper and lower limits of the reservoir capacity (i.e., the beneficial capacity and the dead capacity) (m³). , These are the maximum and minimum design flow rates (m³ / s) for channel j, respectively. The water level of channel j during time period t; , These represent the maximum and minimum design water levels (m) for channel j.
[0104] (3) Solve the above MILP model in multidimensional space using the branch and bound method, and output the decision tensor containing the discharge process of the reservoir group and the headwaters of the main channel. .
[0105] Step S4: Microscopic Channel System Hydrodynamic Simulation and Effective Water Volume Integration
[0106] (1) The decision tensor output in step S3 As upstream boundary conditions We constructed a set of one-dimensional unsteady flow partial differential equations for Saint-Venant hydrodynamics.
[0107]
[0108]
[0109] Where A is the cross-sectional area of the water passage; Q is the flow rate; Z is the water level; x is the spatial step size; t is the time; and g is the acceleration due to gravity. It is a frictional slope.
[0110] (2) When solving the partial differential equations, state variables are introduced to simulate the opening and closing behavior of the channel control gate: when the water depth at the end is... At that time, it was determined to be the pre-accumulation period of channel j. The lower boundary condition is set to gate closure, causing upstream water to be converted into dynamic reservoir capacity to fill the channel; when At that time, it was determined to be the water supply period. This triggers the gate opening boundary condition, releasing the effective water volume.
[0111] (3) The above equations with gate control boundary are discretized and solved using the Preissmann implicit difference scheme to obtain the array of actual physical water depths at the end of the canal. and terminal output water volume Substitute the channel end self-flow discrimination function into step S2(4). By using integral calculations to determine the net effective water supply at the end of each channel, the effective water supply of channel j in the k-th iteration can be expressed as: :
[0112]
[0113] Step S5: Closed-loop damping error feedback and iterative convergence
[0114] (1) Calculate the physical water supply deviation of each channel in the current k-th iteration: .
[0115] (2) Calculate the engineering tolerance threshold : Determine whether the maximum absolute deviation of each channel is less than the preset engineering tolerance threshold: If the conditions are met, the optimal scheduling of the multi-source reservoir-canal system is considered converged, and the current scheduling scheme is output as the final optimal scheduling scheme; if the conditions are not met, and the maximum number of iterations K has not been reached, the system damping learning rate is activated. This triggers the damping error feedback mechanism, updating the macro water supply target for the next round:
[0116]
[0117] (3) If If the maximum number of iterations is reached, termination is triggered, and the suboptimal approximation scheme with the smallest absolute value of deviation in the historical iteration record is output.
[0118] Taking an irrigation district with a complex series-parallel reservoir-canal system as an example, the implementation method is described in detail.
[0119] Step S1: Construction of physical boundaries and back calculation of gravity flow in the multi-source reservoir-canal system of the irrigation district
[0120] (1) Obtain the topological structure of the irrigation district's reservoir group and backbone canal system, and analyze the correlation between the reservoir and canal systems in the irrigation district, such as... Figure 2 As shown, in this embodiment, the system includes main reservoirs (reservoirs 1, 2, and 3) connected in series and side reservoirs (reservoirs 4 and 5) connected in parallel. Channel 1 is supplied by reservoir 1, channel 2 by reservoirs 1 and 2, channel 3 by reservoirs 1, 2, 3, and 4, and channel 4 by all five reservoirs. The maximum safe controlled flow rate of the main channel is limited to 115 m³ / s. Taking channel 2 as an example, its physical parameters are extracted: bottom width B = 20 m, channel length L = 90 km, bottom slope S0 = 0.0003, and Manning roughness n = 0.022. The scheduling cycle is set to T = 20 days. Through remote sensing and elevation determination, the elevation requirement of the terminal outlet water receiving area of channel 2 requires that its minimum effective gravity flow water level must reach [a certain value]. =0.7m, its actual farmland rigid water demand index Demand is =20 million m³. Table 1 shows the basic parameters of irrigation canals and the minimum effective gravity flow water level. Table 2 shows the basic parameters of reservoirs in the irrigation area, assuming that the dead storage capacity of each reservoir is 0 and the initial state of the canals has no channel water storage.
[0121] (2) Inversion flow limit: The system substitutes its physical characteristics into the Manning inversion formula optimized by extreme value compensation, and automatically calculates the lower limit of gravity flow required for Channel 2 to maintain a water level of 0.7m during the main water supply period. Simultaneously, by multiplying by the anti-scouring coefficient β=0.2, its initial pre-storage bed flow rate is set. =4.0m³ / s.
[0122] Table 1 Basic parameters of irrigation canals and minimum effective gravity-flow water level
[0123]
[0124] Table 2 Basic parameters of reservoirs in the irrigation area
[0125]
[0126] Step S2: Construction of Supply and Demand Optimization Scheduling Objectives and Phased Strategies for the Reservoir-Canal System
[0127] (1) The scheduling starts when the canal system is in a non-water storage state. The entire scheduling cycle T is divided into the canal pre-storage period on the time axis. and reservoir water supply period Two non-overlapping subsets satisfying The initial storage capacity of the canal system is constructed as a compensation portion for water conveyance where gravity flow is not possible in the reservoir-canal system. Taking Canal 2 as an example, it is 90 kilometers long. In the initial stage of sluice gate opening, a large amount of water will fill the dry canal bed. The system first calculates the static bottom tank capacity of Canal 2: (B×L× ) / 10000=(20×90000×0.7) / 10000=1.26 million m³. Therefore, in the first optimization iteration (k=1), the system spontaneously adjusts the macroscopic water demand target, setting the initial optimization target for channel 2 as: .
[0128] (2) Construct the self-flow discrimination function at the end of the channel ,when When the water distribution volume of channel j meets the gravity irrigation demand, it is the effective water volume (m³); when At that time, it is considered invalid water volume.
[0129]
[0130] In the formula Let L be the water depth (m) at the end of the channel at time t.
[0131] Step S3: Macroscopic Multi-Source Reservoir-Channel System Water Distribution Flow Optimization Decision
[0132] (1) Construct a system based on the reservoir discharge flow rate Headwater distribution flow and the opening and closing status of the gate For the mixed integer linear programming (MILP) model of the decision space, a minimization objective function is constructed that includes minimizing the discharge, minimizing the channel water distribution time, and minimizing the channel gate opening time.
[0133]
[0134] In the formula Let be the discharge flow rate (m³ / s) of reservoir i at time t. Let be the headwater flow rate (m³ / s) of channel j at time t. Let j be the gate opening state at time t; , , Let be the weight vector, and based on trial calculations, we take values of 100, 5, and 50 respectively. The first term minimizes system leakage; the second term introduces a timing penalty factor t, where t represents time and increases linearly with the scheduling time, thus affecting the weights. The driver optimizer spontaneously allocates water rights in the early stages of scheduling. For the same headwater flow rate, the longer the scheduling time t, the greater the penalty for this term, ensuring that water allocation is concentrated in the early stages to complete the allocation task as quickly as possible. The third term controls the gate opening status. punishment This drives the centralized release of traffic, reducing the duration of channel gate opening.
[0135] (2) Construct a multidimensional constraint set: including reservoir water balance constraints, system water supply and demand constraints, channel water distribution constraints, reservoir discharge flow variation, reservoir water level variation, reservoir water level, reservoir capacity constraints, channel flow, water level constraints, etc. as constraints of the model.
[0136] Reservoir water balance constraints:
[0137] System water supply and demand constraints:
[0138] The water distribution volume in the channel meets the constraints:
[0139] Reservoir discharge flow variation constraints:
[0140] Reservoir water level fluctuation constraints:
[0141] Reservoir water level constraints:
[0142] Reservoir capacity constraints:
[0143] Channel traffic constraints:
[0144] Channel water level constraints:
[0145] in, Let be the reservoir capacity (m³) of reservoir i at time t; Natural inflow rate (m³ / s); The discharge flow rate is (m³ / s). The duration of each time period, in seconds; For channel j at any time The water flow rate is expressed in cubic meters per second. Control the water demand (m³) for irrigation area j in channel j; Let i be the maximum allowable flow rate variation of reservoir i (m³ / s). Let be the water level (m) of reservoir i during time period t. This refers to the fluctuation range of the reservoir water level; , These are the upper and lower limits of the reservoir water level (m); Let be the reservoir capacity (m³) of reservoir i during time period t. , These are the upper and lower limits (m³) of the reservoir capacity, with the upper limit being the beneficial capacity of each reservoir. , These are the maximum and minimum design flow rates (m³ / s) for channel j, respectively. The water level of channel j during time period t; , These represent the maximum and minimum design water levels (m) for channel j. The maximum flow rate of the main canal is set at 115 m³ / s.
[0146] (3) The system inputs the macroscopic objective into a mixed integer linear programming (MILP) algorithm for optimization scheduling. It uses the branch and bound method to solve the MILP model in a multidimensional space and outputs a decision tensor that includes the discharge from the reservoir group and the flow process at the head of the main channel. The reservoir storage and release process is shown in Table 3 and... Figure 3 , Figure 4 When the main water source reservoir 1 needs to release a massive flow, due to the flow constraint of 115 m³ / s on the main stream, the system automatically triggers and quantifies the allocation of spatial commands, scheduling the downstream reservoirs 2 and 3 to reduce their capacity in advance and perform peak shaving and diversion when the flood peak arrives; at the same time, it directs the parallel reservoirs 4 and 5 to provide staggered lateral replenishment. This generates a flow decision result that meets the daily variation constraint.
[0147] Step S4: Microscopic Channel System Hydrodynamic Simulation and Effective Water Volume Integration
[0148] (1) The decision tensor output in step S3 The channel water diversion flow process is used as the upstream boundary condition for each channel. We constructed a set of one-dimensional unsteady flow partial differential equations for Saint-Venant hydrodynamics.
[0149]
[0150]
[0151] Where A is the cross-sectional area of the water passage; Q is the flow rate; Z is the water level; x is the spatial step size; t is the time; and g is the acceleration due to gravity. It is a frictional slope.
[0152] (2) When solving the partial differential equations, state variables are introduced to simulate the opening and closing behavior of the channel control gate: when the water depth at the end is... At that time, it was determined to be the pre-accumulation period of channel j. The lower boundary condition is set to gate closure, causing upstream water to be converted into dynamic reservoir capacity to fill the channel; when At that time, it was determined to be the water supply period. This triggers the gate opening boundary condition, releasing effective water volume. Taking Canal 2 as an example, when the equation reaches the end of Canal 2, since the initial water depth is less than 0.7m, the system determines the GFEI index to be 0, and the lower boundary condition is locked as gate closure. At this point, the evolving physical water flow is forcibly intercepted, and the channel water level rises rapidly; when the water depth at the end exceeds the 0.7m threshold, the system determines that the GFEI index jumps to 1, triggering the gate opening boundary condition and releasing the net agricultural water demand.
[0153] (3) The above equations with gate control boundary are discretized and solved using the Preissmann implicit difference scheme to obtain the array of actual physical water depths at the end of the canal. and terminal output water volume Substitute the channel end self-flow discrimination function into step S2(4). By using integral calculations to determine the net effective water supply at the end of each channel, the effective water supply of channel j in the k-th iteration can be expressed as: :
[0154]
[0155] At the end of the water supply cycle, the actual effective water volume released by the gate at the end of the channel during the kth iteration is calculated using the GFEI integral formula.
[0156] Step S5: Closed-loop damping error feedback and iterative convergence
[0157] (1) Calculate the physical water supply deviation of each channel in the current k-th iteration: .
[0158] (2) Calculate the engineering tolerance threshold : Determine whether the maximum absolute deviation of each channel is less than the preset engineering tolerance threshold: . The calculated value is set at 150,000 m³. If this condition is met, the multi-source reservoir-canal system optimization scheduling is considered converged, and the current scheduling scheme is output as the final optimal scheduling scheme. If this condition is not met, and the maximum number of iterations K has not been reached (K is set to 20), then the system damping learning rate is activated. , After trial calculations, the value was set to 0.75, triggering the damping error feedback mechanism and updating the macro water supply target for the next round. The new target is implicitly transmitted back to step S3, and the MILP central hub readjusts the release allocation and timing relay of the reservoir group.
[0159]
[0160] (3) If If the maximum number of iterations is reached, termination is triggered, and the suboptimal approximation scheme with the smallest absolute value of deviation in the historical iteration record is output.
[0161] Table 3. Storage and release processes of the reservoir group
[0162]
[0163] Table 4. Headwater Flow Process
[0164]
[0165] From Table 3, Figure 3 , Figure 6 Analysis shows that the total inflow of the system is 12 million m³, the reservoir group consumes 105.12 million m³ of storage capacity, and the effective water demand above the gravity flow level is 115.2 million m³, with only a small amount of water (approximately 1.92 million m³) remaining in the channels. This means that the total gross water volume released by the reservoir group is only slightly more than the demand, exceeding the amount of water stored in the channels. The effective water volume output at the end (the integral of GFEI=1) is sufficient to meet the demand indicators. As shown in Table 4, channel water intake is relatively concentrated in the early stages of the scheduling process. Corresponding to the reservoir release, the total water demand for Channel 3 is relatively small, with irrigation concentrated in the first two days. Channel 4's water intake on the last day is for the next scheduling cycle, and the actual water received is not included in this round of irrigation.
[0166] Depend on Figure 4 , Figure 5 It can be seen that by stimulating the potential of spatiotemporal relay scheduling of a series and parallel reservoir group, in complex water diversion projects where tributaries converge and flow carrying capacity is limited, peak staggering in both time and space effectively avoids the main canal and canal system from exceeding the flow limit.
Claims
1. A method for optimizing the scheduling of a multi-source reservoir-canal system in an irrigation district based on the effective gravity-flow water level at the end of the canal, characterized in that, Includes the following steps: S1. Obtain the topology and engineering parameters of the irrigation district reservoir group and backbone canal system. The backbone canal system includes several backbone channels. Combine remote sensing images, UAV images and elevation data to determine the minimum effective gravity flow water level at the end of each backbone channel to achieve gravity flow water discharge. Then, deduce the minimum water supply flow lower limit and the pre-storage flow threshold before formal water supply for each backbone channel under the condition of maintaining its effective gravity flow water level at the end, and form the scheduling bottom line constraint boundary. S2. Divide the entire scheduling cycle into the channel pre-storage period and the reservoir water supply period, construct the initial channel storage capacity, and establish the channel end gravity flow discrimination function to determine whether the output water volume at the end of the backbone channel constitutes effective water supply, so as to obtain the initial water supply target required for macro-optimization, and at the same time give the target water demand of each backbone channel in the scheduling cycle. S3. Using reservoir discharge flow, canal headwater flow, and gate opening / closing status as decision variables, establish a mixed integer linear programming model that includes minimizing discharge, minimizing canal distribution time, and minimizing gate opening time. Solve the model under the constraints of reservoir water balance, system supply and demand balance, canal flow balance, reservoir and canal boundary constraints, and water level and reservoir capacity variation constraints. Output the reservoir discharge process and the main canal headwater flow process for each time period. S4. The discharge process of the reservoir group and the flow process at the head of the canal are used as the upstream boundary conditions of the one-dimensional unsteady flow dynamics model. The Saint-Venant equations with the control gate opening and closing are constructed. The boundary conditions are switched between the pre-storage period and the water supply period according to whether the water depth at the end of the canal reaches the minimum effective gravity flow level. The implicit difference scheme is used to solve the problem to obtain the actual water depth at the end of each main canal and the output water volume at the end. The output water volume at the end is integrated by the gravity flow discrimination function at the end of the canal to calculate the net effective water supply of each main canal. S5. Calculate the deviation between the net effective water supply of each backbone channel and the target water demand, and use whether the maximum deviation among all backbone channels is less than the preset engineering tolerance threshold as the convergence criterion; when the convergence condition is not met and the maximum number of iterations is not reached, the initial water supply target is corrected according to the damping learning rate, and the process returns to step S3 to re-optimize until the convergence condition is met or the maximum number of iterations is reached.
2. The method according to claim 1, characterized in that, The reverse derivation of the minimum water supply flow limit for each backbone channel under the condition of maintaining the effective gravity flow water level at its end, as described in step S1, is performed using the Manning formula optimized with extreme value compensation: ;in, The minimum water supply flow rate for the main channel j is expressed in cubic meters per second. Let J be the Manning roughness of the backbone channel j, which is dimensionless; The width of the bottom of the main channel j is in meters. The minimum effective gravity-flow water level of the main channel j, in meters; The bottom slope of the main channel j is dimensionless.
3. The method according to claim 2, characterized in that, The formula for calculating the pre-stored flow threshold mentioned in step S1 is as follows: ;in, The threshold flow rate required for the main channel to store water before formal water supply, expressed in cubic meters per second; , is the pre-storage ratio coefficient for the main channel j, dimensionless, and its value range is determined according to the channel's scour resistance requirements and the need for stable water lifting.
4. The method according to claim 1, characterized in that, The initial storage capacity of the canal system mentioned in step S2 is determined by the geometric parameters of each main canal and the minimum effective gravity-flow water level. The initial water supply target is expressed as follows in the first iteration: ;in, The initial macroscopic water supply target for the first iteration is expressed in cubic meters. The target water demand for the irrigation area corresponding to the main channel j within the scheduling cycle is expressed in cubic meters. The width of the bottom of the main channel j is in meters. Here is the length of the main channel j, in meters; This refers to the total number of core channels; The first term represents the minimum effective gravity-flow water level of the main channel j, in meters; the second term indicates that the initial storage capacity of the canal system is included as part of the water supply compensation in the macro water supply target.
5. The method according to claim 1, characterized in that, The formula for calculating the channel end gravity flow discrimination function mentioned in step S2 is as follows: ;in, is the gravity flow discrimination index of the main channel j at time t. A value of 1 indicates that the water output at the end of the channel is an effective water supply, and a value of 0 indicates an ineffective water supply. The water depth at point L, the end of the main channel j, at time t is expressed in meters. The minimum effective gravity-flow water level of the main channel j is expressed in meters.
6. The method according to claim 1, characterized in that, The objective function of the mixed-integer linear programming model described in step S3 is: ;in, Let be the outflow rate of reservoir i at time t, expressed in cubic meters per second. The headwater flow rate of the main channel j at time t is expressed in cubic meters per second. Let be the gate opening / closing state variable of the backbone channel j at time t; , , This is a weight vector; the first term minimizes the total system discharge; the second term introduces a time-series penalty factor t, which is related to the weights. The co-driving optimizer allocates water rights in the early stages of scheduling to minimize the channel water distribution time; the third term is based on the gate opening state variable. Apply penalty weights This drives the concentrated release of flow to reduce the total time the gate is open.
7. The method according to claim 6, characterized in that, The constraints of the mixed-integer linear programming model described in step S3 include at least one of the following: Reservoir water balance constraints: ,in, Let be the reservoir capacity of reservoir i at time t, in cubic meters; Let be the reservoir capacity of reservoir i at time t-1, in cubic meters; The natural inflow rate of reservoir i at time t is expressed in cubic meters per second. Let be the outflow rate of reservoir i at time t, expressed in cubic meters per second. The discharge flow rate is expressed in cubic meters per second; Δt is the duration of each time period, expressed in seconds. System water supply and demand constraints: That is, the total discharged water volume shall not be less than the sum of the target water requirements of each irrigation area; among which, Δt represents the target water demand of the irrigation area corresponding to the main canal j within the scheduling cycle, in cubic meters; Δt is the duration of each time period, in seconds. The water distribution volume in the channel meets the constraints: ,in, For channel j at any time The water flow rate is expressed in cubic meters per second. The water demand for irrigation sections is controlled by channel j, in cubic meters; Reservoir discharge flow variation constraints: ,in, Let be the discharge flow rate of reservoir i at time t-1, in cubic meters per second; The maximum allowable flow variation for reservoir i is expressed in cubic meters per second. Reservoir water level fluctuation constraints: ,in, Let be the water level of reservoir i at time t, in meters. Let be the water level of reservoir i during time period t-1, in meters. The maximum permissible water level fluctuation for reservoir i, in meters; Reservoir water level constraints: ,in, , These are the lowest and highest operating water levels of reservoir i, respectively. Reservoir capacity constraints: ,in, Let be the reservoir capacity of reservoir i at time t, in cubic meters; , These are the dead storage capacity and the beneficial storage capacity of reservoir i, respectively. Channel traffic constraints: ,in, The headwater flow rate of the main channel j at time t is expressed in cubic meters per second. and These are the maximum and minimum design flow rates of the backbone channel j, respectively. Channel water level constraints: ,in, The water level of the main channel j during time period t, in meters. and These are the maximum and minimum design water levels of the main channel j, respectively.
8. The method according to claim 1, characterized in that, The one-dimensional unsteady flow dynamics model described in step S4 is described by the Saint-Venant equations: ; Where A is the cross-sectional area of the water passage, in square meters; Q is the flow rate, in cubic meters per second; Z is the water level, in meters; x is the spatial step size, in meters; t is the time, in seconds; and g is the gravitational acceleration, in meters per square second. For frictional slopes, dimensionless; The Saint-Venant equations are discretized and solved using the Preissmann implicit difference scheme, and state variables are introduced to simulate the opening and closing behavior of the channel control gate: when the water depth at the end of the channel... When the time is determined to be the pre-storage period, the lower boundary condition is set to gate closure, causing upstream water to be converted into dynamic reservoir capacity to fill the channel; when When the water supply period begins, the gate opening boundary conditions are triggered to release the effective water volume; among which, The water depth at point L, the end of the main channel j, at time t is expressed in meters. The minimum effective gravity-flow water level of the main channel j is expressed in meters.
9. The method according to claim 1, characterized in that, The formula for calculating the net effective water supply of each backbone channel mentioned in step S4 is as follows: ;in, This represents the net effective water supply of the backbone channel j in the k-th iteration, in cubic meters. The output flow rate at the end of the backbone channel j at time t, expressed in cubic meters per second; The total duration of the scheduling cycle is in seconds; The gravity flow discrimination index of the main channel j at time t is used to isolate and integrate the effective water volume.
10. The method according to claim 1, characterized in that, Step S5, which involves feedback correction of the initial water supply target according to the damping learning rate, includes: Calculate the water supply deviation of each main channel in the current k-th iteration: ;in, The target water demand for the irrigation area corresponding to the main channel j within the scheduling cycle is expressed in cubic meters. Let be the net effective water supply of backbone channel j in the k-th iteration, in cubic meters; calculate the maximum absolute deviation among all backbone channels: ;in, This refers to the total number of core channels; like If the current iteration number k is greater than or equal to the preset engineering tolerance threshold ε, and the current iteration number k is less than the maximum iteration number K, then the macroscopic water supply target for the next iteration is updated using the following damping error feedback formula: ;in, Let α be the macroscopic water supply target for the k-th iteration, in cubic meters; α be the system damping learning rate, dimensionless, ranging from (0,1]; and the updated... As a new supply and demand balance constraint right-hand term in the mixed-integer linear programming model in step S3, it drives iterative optimization until convergence or the maximum number of iterations is reached.