Regulation method for risk resilience of water engineering group considering uncertainty

Abstract

The application discloses a kind of water engineering group risk resilience regulation methods considering uncertainty, comprising: based on the water engineering group condition and hydrological data obtained, establish double-track dynamic hierarchical structure, and generate the resilience performance boundary of the lower limit of system performance and the current resilience performance state considering uncertainty;By identifying performance gap and tracking cross-level failure chain, determine the regulation potential of resilience failure engineering, and divide it into two kinds of failure engineering sets, controllable and uncontrollable;For the two kinds of failure engineering sets, respectively using multi-objective reinforcement learning guided hybrid gradient search algorithm for self-regulation optimization, and using collaborative-competitive co-evolution algorithm for collaborative regulation optimization, to output the optimal scheduling scheme set.The application can dynamically recognize system structure, accurately trace risk transmission, and significantly improve the risk active defense capability and system resilience of water engineering group under uncertainty disturbance.

Description

TECHNICAL FIELD

[0001] The present application relates to water conservancy dispatching technology, and in particular to a water project group risk resilience regulation method considering uncertainty. BACKGROUND

[0002] As a complex system composed of reservoirs, dikes, sluice stations and other water conservancy facilities, the water project group is the core physical carrier for flood control and disaster reduction, optimal allocation of water resources and ecological environment protection in a river basin. Therefore, how to improve the regulation and control capability of the water project group in the face of super-standard floods, forecast uncertainty and other strong disturbances has become a frontier scientific problem to be solved in the field of water conservancy, which has great theoretical and practical significance for ensuring the safety of the river basin and promoting sustainable development.

[0003] Currently, the dispatching optimization research of the water project group has made certain progress. Most of the existing technical solutions are based on a determined flood forecast process, and use conventional optimization algorithms such as linear programming, dynamic programming or genetic algorithm to optimize the dispatching rules of the project group. In terms of system description, a hydraulic model is usually constructed based on fixed and physical connection relationships (such as upstream and downstream river connection) between projects to simulate the flood evolution process. In terms of dispatching objectives, these methods often focus on a single and deterministic flood control index, such as minimizing the flood peak flow or the highest submerged water level of a certain key section downstream as the main optimization objective. In actual operation, many water project groups still rely on relatively fixed dispatching procedures or flood limit water level rules based on historical experience, and operate in stages and according to rules during floods to ensure basic flood control safety. These methods can play an effective role in dealing with conventional flood events within the standard and with high prediction accuracy, and are an important technical basis for the safe operation of the water project group.

[0004] However, with the deepening of the understanding of the complexity of the water project group system and the uncertainty of the operating environment, the existing technical solutions have exposed a series of deep-seated technical problems in dealing with future challenges. Mainly reflected in the lack of understanding of the dynamic characteristics of the system and the unclear description of the risk transmission mechanism, which limits the accuracy and foresight of the regulation and control decision.

[0005] Specifically, the current method relies too much on fixed physical topology to describe the system and fails to reveal the dynamic water power correlation between projects under extreme hydrological conditions, which causes the decision basis to be disconnected from the real state of the system. When a risk occurs, the existing technology can only identify isolated failure points, but cannot track the root cause of the failure and its complete chain of cross-level and multi-dimensional propagation in the system. The lack of these deep mechanisms makes it difficult to accurately implement regulation and control measures and also makes it impossible to effectively predict the potential cascading collapse risk. SUMMARY

[0006] The present application provides a water engineering group risk resilience regulation method considering uncertainty.

[0007] The technical scheme of the present application is a water engineering group risk resilience regulation method considering uncertainty, comprising:

[0008] Based on the obtained water engineering group working conditions and hydrological data, a double-track dynamic hierarchical structure is established, and a resilience performance boundary and a current resilience performance state are generated;

[0009] According to the double-track dynamic hierarchical structure, the resilience performance boundary and the current resilience performance state, the regulation potential of resilience failure engineering is determined, and it is divided into two types of failure engineering sets, i.e., controllable and uncontrollable;

[0010] For the two types of failure engineering sets, self-regulation optimization and collaborative regulation optimization are respectively performed, and an optimal scheduling scheme set is output.

[0011] Optionally, the double-track dynamic hierarchical structure is established, comprising:

[0012] The permanent engineering connection relationship is read from the water engineering group working condition data to construct a steady-state physical skeleton hierarchy;

[0013] In response to the real-time hydraulic state presented by the hydrological data, a hydraulic coupling tensor is established to quantify the dynamic correlation strength between the engineering, and a time-varying dynamic response hierarchy is generated according to the hydraulic coupling tensor;

[0014] And based on the hydraulic coupling tensor, the physical skeleton hierarchy and the dynamic response hierarchy are weighted and fused to obtain the double-track dynamic hierarchical structure.

[0015] Optionally, the hydraulic coupling tensor is constructed, comprising:

[0016] According to the real-time hydraulic state, the water level propagation relationship representing the mutual influence of water levels between engineering, the flow propagation relationship representing the mutual influence of flow, and the water level-flow coupling relationship representing the influence of water level on downstream flow are calculated respectively;

[0017] And the three relationships are quantified as a water level propagation matrix, a flow propagation matrix and a water level-flow coupling matrix respectively;

[0018] The water level propagation matrix, the flow propagation matrix and the water level-flow coupling matrix are combined by weighting to synthesize the hydraulic coupling tensor.

[0019] Optionally, the weighted fusion of the physical skeleton hierarchy and the dynamic response hierarchy based on the hydraulic coupling tensor comprises:

[0020] The spatial variance of the hydrodynamic coupling tensor is calculated, and a time-varying hierarchical fusion weight coefficient is determined according to the spatial variance.

[0021] The hierarchical fusion weight coefficient is used to perform weighted combination on the physical skeleton hierarchy and the dynamic response hierarchy, so as to generate a dual-track dynamic hierarchical structure.

[0022] Optionally, before determining the regulation potential of the toughness failure project, the method further comprises:

[0023] The current toughness performance state is compared with the toughness performance boundary, and a performance gap is calculated.

[0024] The failure risk probability of the performance gap being lower than a preset value under the uncertainty disturbance is counted.

[0025] When the failure risk probability and the absolute mean of the performance gap both exceed the corresponding probability threshold and performance gap threshold, respectively, the corresponding project is determined as a toughness failure node, and a toughness failure node set is formed.

[0026] Optionally, after obtaining the toughness failure node set, the method further comprises:

[0027] Based on the current toughness performance state and the toughness performance boundary, a node with a performance surplus is screened out to form a regulation source set.

[0028] For any toughness failure node in the toughness failure node set, a closed-loop collaborative path connecting the failure node and at least one regulation source is searched in the dual-track dynamic hierarchical structure.

[0029] At least one feasibility score of the closed-loop collaborative path is calculated to represent the regulation potential of the toughness failure node.

[0030] Optionally, calculating the feasibility score of the closed-loop collaborative path comprises:

[0031] At least two performance indicators of the closed-loop collaborative path are comprehensively evaluated, including path length, response total time, regulation efficiency, and net performance improvement.

[0032] Based on the evaluation results of the at least two performance indicators, the feasibility score is generated.

[0033] Optionally, the toughness failure projects are divided into two categories of regulatable and non-regulatable failure project sets, comprising:

[0034] The feasibility score is compared with a preset feasibility score threshold.

[0035] If the feasibility score is not less than the feasibility score threshold, the corresponding toughness failure node is divided into the regulatable failure project set.

[0036] Otherwise, the ductile failure nodes are divided into the uncontrollable failure engineering set.

[0037] Optionally, the self-regulation optimization and the collaborative regulation optimization are respectively performed on the two types of failure engineering sets, including:

[0038] For the controllable failure engineering set, a multi-objective reinforcement learning guided hybrid gradient search algorithm is used for optimization.

[0039] For the uncontrollable failure engineering set, a co-evolution algorithm of cooperation-competition is used for optimization, in combination with the matched auxiliary regulation engineering set.

[0040] Optionally, for the controllable failure engineering set, a multi-objective reinforcement learning guided hybrid gradient search algorithm is used for optimization, including:

[0041] Based on the physical constraints of the controllable failure engineering set, a safe exploration space for limiting the range of regulation actions is established.

[0042] In the safe exploration space, a multi-objective reinforcement learning agent is used to generate an initial regulation action.

[0043] And a gradient search algorithm is used to fine-tune the initial regulation action, while ensuring that the adjusted action does not exceed the boundary of the safe exploration space, to obtain an optimized regulation action sequence.

[0044] The beneficial effects can dynamically recognize the system structure, accurately trace the risk transmission, and significantly improve the risk active defense capability and system resilience of the water engineering group under uncertain disturbances. The related technical effects will be described later. BRIEF DESCRIPTION OF DRAWINGS

[0045] Fig. 1 is the flowchart of the present application.

[0046] Fig. 2 is the flowchart of the present application for establishing a dual-track dynamic hierarchical structure.

[0047] Fig. 3 is the flowchart of the present application for constructing a hydraulic coupling tensor.

[0048] Fig. 4 is the flowchart of the present application for weighting and fusing the physical skeleton hierarchy and the dynamic response hierarchy based on the hydraulic coupling tensor. DETAILED DESCRIPTION

[0049] In order to solve the above-mentioned problems existing in the prior art, the applicant has made in-depth analysis and found that:

[0050] Existing methods generally rely on static system topology and cannot accurately capture the real-time changes in dynamic hydraulic coupling between projects during floods. The system structure of water engineering groups is not immutable. Under extreme flood conditions, due to factors such as river channel waterlogging, tributary backwater or inter-zone flood surge, temporary and significant coupling relationships may be formed between projects that have weak hydraulic connections under normal conditions, and the original key hydraulic transmission path may also change. Existing methods are based on fixed physical connections to build models and make decisions. This static understanding ignores the dynamic response characteristics of the system structure, which leads to the inability to identify the real dominant influence path and key node during floods, and thus may make suboptimal or even incorrect control judgments based on outdated or inaccurate system states.

[0051] Existing risk assessment and control response lack the ability to trace the deep root cause and structured transmission path of resilience failure. Traditional risk identification is usually point-like and single-dimensional, that is, to judge whether a single project's indicator (such as water level) is over-limit. This method cannot answer a deeper question: what is the root cause of a system-level resilience failure (such as a decrease in the recovery ability of the entire watershed), which performance (such as resistance, absorption) of which project node first fails, and how this failure is transmitted and amplified between different levels (from a single project to a subsystem to a large system) through the interaction between projects. Due to the lack of tracking and diagnosis capabilities of cross-level failure chains, existing methods are often lagging and lack of targetedness when facing systemic risks, making it difficult to achieve precise attack on the root cause of the risk and early intervention on potential cascading failures.

[0052] To this end, the following solutions are provided, and in combination Figs. 1 to 4 Various embodiments of the present application are described.

[0053] Embodiment one, describes the overall process of a water engineering group risk resilience control method considering uncertainty.

[0054] In this embodiment, a water engineering group risk resilience control method considering uncertainty includes:

[0055] Step 1, based on the obtained water engineering group working condition and hydrological data, a double-track dynamic hierarchical structure is established, and a resilience performance boundary and a current resilience performance state are generated.

[0056] This step corresponds to a system initialization and state evaluation phase. Specifically, the system first reads the basic data of the water engineering group, such as the geographical location of the project, the physical connection relationship, the reservoir capacity, the discharge capacity, and other working condition data, as well as real-time hydrological data such as water level and flow. Then, based on these data, a dual-track hierarchical structure is constructed that can reflect both the fixed physical topology between projects and the real-time hydraulic dynamic correlation. At the same time, through the simulation of historical extreme flood scenarios, the resilience performance boundary representing the performance bottom line of the system under the most unfavorable conditions is calculated and generated; and under the consideration of the composite disturbance of the current flood forecast uncertainty, the system resilience performance state corresponding to the current control scheme is evaluated.

[0057] In the second step, based on the hierarchical structure, performance boundary, and performance state, the control potential of resilience failure projects is determined, and they are divided into two categories: controllable and uncontrollable failure project sets.

[0058] This step corresponds to a risk identification and potential determination phase. The system first compares the current resilience performance state with the resilience performance boundary, quantifies the performance gap, and identifies the resilience failure project nodes by combining their probability of exceeding the safety threshold under uncertainty disturbance. For these failure nodes, a self-top-down tracking mechanism is used to locate and construct failure transmission chains across different levels using the dual-track hierarchical structure. Based on this, the system further analyzes whether the failure projects can restore their performance through internal or external collaborative control, i.e., determines their control potential, and finally divides all failure projects into two sets: controllable and uncontrollable, the latter of which will also be matched with corresponding auxiliary control resources.

[0059] In the third step, for the two sets of failure projects, self-control optimization and collaborative control optimization are performed respectively to output an optimal scheduling scheme set.

[0060] This step corresponds to a multi-objective scheduling optimization phase. For the controllable failure project set, the system uses an optimization algorithm that focuses on exploiting its own adjustment capability, such as the multi-objective reinforcement learning guided hybrid gradient search algorithm (MORL-HG), for adaptive control. For the uncontrollable failure project set, the system jointly matches the auxiliary control project set for it, and uses an optimization algorithm that focuses on multi-agent collaboration and resource competition, such as the co-evolution algorithm (CCEA-LF), for collaborative and synergistic control. The schemes generated by the two types of optimization algorithms are integrated and screened to output a set of Pareto optimal scheduling schemes for decision-makers to use.

[0061] In this embodiment, by constructing a dynamic hydraulic coupling tensor based on real-time hydrological data calculation, and establishing a double-track dynamic hierarchical structure accordingly, the technical defects of the existing method relying on static topology and unable to reflect the real hydraulic connection are solved. This scheme makes the cognition of system structure change from static to dynamic real-time mapping, which can capture the temporary strong coupling path formed by backwater, top support and other reasons during extreme floods. All subsequent control decisions are based on the most real system state, which improves the accuracy and pertinence of decision-making. By introducing the tracking and diagnosis mechanism of cross-level resilience failure chain, the problem that the existing risk assessment can only identify isolated failure points and cannot trace the risk deeply is solved. When systemic risk occurs, this method can trace and build a complete failure transmission path from top to bottom, accurately positioning the root cause of the problem and its specific failure dimension. The control response changes from blind response to precise attack on the source of risk, which can effectively prevent the occurrence of cascading failure. In summary, by dynamically cognizing the system structure and accurately tracing the risk transmission, the invention improves the risk proactive defense capability and overall resilience of the water engineering group under uncertain disturbance, making the scheduling and control more forward-looking, scientific and effective.

[0062] Embodiment two, this embodiment mainly describes the construction process of the double-track dynamic hierarchical structure, especially the calculation method of the hydraulic coupling tensor.

[0063] In this embodiment, the process of establishing a double-track dynamic hierarchical structure includes the following steps:

[0064] Step one, analyze the permanent engineering connection relationship included in the working condition data of the water engineering group, and construct a stable physical skeleton hierarchy.

[0065] The specific implementation is to read the engineering number, geographical location, upstream and downstream relationship, and permanent connection information defined by the hydraulic structures (such as dams, channels, tunnels) in the basic data of the water engineering group. Based on this information, an engineering topology adjacency matrix A is constructed physical , where if the matrix element A ij =1, it means that there is a direct and fixed physical connection between engineering i and engineering j, otherwise it is 0. On this basis, according to the natural river network structure of the basin, the single engineering node is clustered to form a subsystem (for example, the reservoir group on the same branch), and then each subsystem is combined to form a large system covering the whole basin, thereby forming a relatively stable physical skeleton hierarchy L physical with a single engineering-subsystem-large system three-layer structure.

[0066] Step two, in response to the real-time hydraulic state presented by the hydrological data, a hydraulic coupling tensor is established to quantify the dynamic correlation strength between engineering, and a time-varying dynamic response hierarchy is generated according to the hydraulic coupling tensor.

[0067] In the embodiment, the process of constructing a hydraulic coupling tensor is as follows:

[0068] According to the real-time hydraulic state, the water level propagation relationship representing the mutual influence of water levels between projects, the flow propagation relationship representing the mutual influence of flows, and the water level-flow coupling relationship representing the influence of water levels on downstream flows are calculated respectively.

[0069] The three relationships are quantified into a water level propagation matrix, a flow propagation matrix, and a water level-flow coupling matrix respectively.

[0070] The specific implementation is to read real-time water level data and flow data, calculate the hydraulic propagation time τ ij (t)=L ij / (g×h avg ) 0.5 , wherein L ij is the river distance from project i to project j, g is the acceleration of gravity, and h avg is the average water depth of the interval.

[0071] The water level propagation matrix H HH is calculated as follows: H HH,ij =(dH j / dH i )×exp(-τ ij (t) / T0); wherein H HH,ij represents the decayed intensity of the water level change of project i on the water level of project j; dH j / dH i is the partial derivative of the water level influence, which can be approximately calculated by a hydraulic model or a finite difference method; exp(-τ ij (t) / T0) is an exponential decay term representing the weakening of the influence with the propagation time; T0 is a characteristic time constant, which can be determined according to the flood propagation characteristics of the basin, for example, it can be 50% of the time required for a typical flood wave to propagate to the downstream control section, such as 6 hours.

[0072] The flow propagation matrix H QQ is calculated as follows: H QQ,ij =(dQ j / dQ i )×exp(-τ ij (t) / T0); wherein H QQ,ij represents the decayed intensity of the flow change of project i on the flow of project j; dQ j / dQ i is the partial derivative of the flow influence.

[0073] The water level-flow coupling matrix H HQ is calculated as follows: H HQ,ij =(dQ j / dHi ) x exp(-τ ij (t) / T0); wherein H HQ,ij represents the decayed intensity of the water level change of project i on the flow of project j; dQ j / dH i is the partial derivative of the water level-flow influence.

[0074] The water level propagation matrix, the flow propagation matrix, and the water level-flow coupling matrix are combined by weighting to synthesize the hydraulic coupling tensor.

[0075] The integrated hydraulic coupling intensity is calculated as: H ij (t) = w1 x |H HH,ij | + w2 x |H QQ,ij | + w3 x |H HQ,ij |; wherein H ij (t) is the integrated hydraulic coupling intensity between project i and j at time t; w1, w2, w3 are weight coefficients, which can be determined based on expert experience or hydraulic propagation characteristic analysis, for example, for a system mainly controlled by water level, w1 = 0.4, w2 = 0.3, w3 = 0.3 can be set; |*| represents taking the absolute value to ensure that the coupling intensity is non-negative. The N x N matrix obtained by this calculation is the hydraulic coupling tensor matrix H(t). On this basis, a coupling intensity threshold H threshold is set (for example, the 95th percentile of historical data statistics can be taken). When H ij (t) > H threshold , it is considered that there is temporary strong hydraulic coupling between project i and j, and a temporary edge is established in the dynamic graph. Subsequently, a community discovery algorithm such as spectral clustering is used to cluster this dynamically generated graph, thereby obtaining the time-varying dynamic response hierarchy L dynamic (t).

[0076] Step three, based on the hydraulic coupling tensor, the physical skeleton hierarchy and the dynamic response hierarchy are weighted and fused to jointly build a dual-track dynamic hierarchy structure.

[0077] In this step, the spatial variance of the hydraulic coupling tensor is calculated, and a time-varying hierarchy fusion weight coefficient is determined according to the spatial variance. Specifically, the spatial variance V ariance (H ij ) of all element values in the hydraulic coupling tensor matrix H(t) is calculated. This variance can reflect the uniformity or heterogeneity of the current overall system hydraulic connection. When the variance is small, it indicates that the system hydraulic state is stable, and at this time, more emphasis should be placed on the stable physical connection; on the contrary, when the variance is large, it indicates that there may be strong coupling caused by local floods or abnormal scheduling, and at this time, more attention should be paid to the dynamic response relationship. The hierarchy fusion weight coefficient is calculated as:

[0078] a(t) = 1 / (1 + K x V ariance (H ij )); wherein a(t) is the hierarchical fusion weight coefficient at time t, whose value range is (0, 1]; K is the sensitivity coefficient, used to adjust the sensitivity of variance to weight, which can be calibrated according to historical data, or take an empirical value, for example, K = 0.1.

[0079] Correspondingly, the hierarchical fusion weight coefficient is used to weight and combine the physical skeleton hierarchy and the dynamic response hierarchy to generate a dual-track dynamic hierarchy. The final hierarchy L final Through weighted fusion decision determination, for example, for a specific decision (such as judging whether two projects belong to the same subsystem), the decision result can be represented as Decision = a(t) x D physical + (1 - a(t)) x D dynamic , wherein D physical and D dynamic are the decision results based on the physical hierarchy and the dynamic hierarchy respectively. In this way, a dual-track hierarchical structure considering both static structure and dynamic response is generated.

[0080] Embodiment Three, this embodiment describes an index system for quantifying the resilience performance boundary and the current resilience performance state, especially the multi-hierarchy aggregation method and the weight determination method.

[0081] In this embodiment, the resilience performance boundary and the current resilience performance state are quantified based on a dimensionless performance index system including the following four dimensions: resistance performance; absorption performance; recovery performance; and adaptability performance.

[0082] First, a series of normalized structure characteristic factors need to be calculated, such as reservoir capacity utilization rate C u , safety elevation margin H s , response lag time of control action T d , historical failure frequency θ, etc., which are the basis for calculating four types of performance indicators.

[0083] Specifically, at the single project level, the calculation method of each dimension performance index is as follows:

[0084] Resistance performance Pr(single) = ω r1 x C u,norm + ω r2 x H s,norm + ω r3 x exp(-λ x T d,norm ); this index represents the initial ability of the project to resist disturbances, wherein C u,norm , H s,norm , T d,normare the normalized reservoir utilization rate, the safety elevation margin and the response lag time of the control action, respectively; ωa1, ωa2, ωa3 are the respective weights. r1 ,ω r2 ,ω r3 are the respective weights.

[0085] absorptive performance Pa (single) = (1- θ norm ) ωa1 × (1- α norm ) ωa2 × (1- η norm ) ωa3 ; the index represents the ability of the project to absorb and withstand damage when subjected to impact, wherein θ norm , α norm , η norm are the normalized historical failure frequency, disturbance sensitivity and functional coupling degree, respectively; ωa1, ωa2, ωa3 are the weights.

[0086] restoring force performance Ph (single) = κ norm ωh1 / (1+ Δt norm / T ref ) ωh2 ; the index represents the speed and degree of recovery of the project from the damaged state, wherein κ norm , Δt norm are the normalized recovery elastic coefficient and recovery time, respectively; T ref is the reference time constant; ωh1, ωh2 are the weights.

[0087] Adaptive force performance: adaptive force is more reflected at the system level, and a single project layer can not be calculated alone, or a comprehensive learning ability index can be used instead.

[0088] Further, based on the performance index system for quantification, further including the step of aggregating the indexes at the subsystem level.

[0089] In this embodiment, this step aims to reflect the short board effect and synergistic effect of system engineering, and the specific aggregation method includes but is not limited to:

[0090] subsystem resistance performance Pr (sub) =∑wi×P r,i (单) ; wherein wi is the importance weight of the i-th project in the subsystem, and the weighted average method is used for calculation.

[0091] The minimum value principle is adopted to determine the minimum value of the absorptive performance index of each single project in the subsystem as the absorptive performance of the subsystem, so as to reflect the short board effect. Namely, Pa (子) =min{P a,i (单)This method ensures that the absorption capacity of the subsystem is determined by its weakest link.

[0092] The harmonic mean of the recovery performance indicators of each single project in the subsystem is calculated using the principle of harmonic mean, and is determined as the recovery performance of the subsystem to reflect the synergistic recovery characteristics.

[0093] That is, Ph (子) =n / ∑(1 / P h,i (单) ), where n is the number of projects in the subsystem. Harmonic mean pays more attention to individuals with smaller values, and can reflect that the overall recovery efficiency is significantly affected by the slowest recovery project.

[0094] The subsystem adaptability performance Ps (子) =ωs1×ln(1+C c,norm ×D l,norm ); where C c,norm and D l,norm are normalized factors representing the command and coordination ability and scheduling flexibility of the subsystem, respectively.

[0095] Optionally, at the level of the large system, the subsystem performance indicators can be further aggregated, for example, the large system resistance can consider the spatial distribution of unevenness, represented as Pr (大) =μ(Pr (子) )-σ(Pr (子) ), that is, the mean minus the standard deviation of the subsystem resistance, so as to reduce the weight of the subsystem with large resistance dispersion.

[0096] Using the double short board effect, not only the weakest subsystem is considered, but also the aggregation degree of the vulnerable subsystem is reflected through the tail mean, and the large system absorption is represented as Pa (大) =min{P a,j (子)}×∑P a,j (子,弱) / k, that is, the shortest board absorption of the subsystem multiplied by its tail mean, and the default k=0.3n.

[0097] The calculation of the large system recovery force mainly uses the harmonic mean, while the coefficient of variation is used to punish the dispersion of the recovery speed, that is, Ph (大) =n / ∑(1 / P h,j (大) )×(1-σ(Pr (子) ) / μ(Pr (子) )).

[0098] The large system adaptability performance Ps (大) =α / μ(Ps (子) )+(1-α)×min{P s,j(子) The overall level and the bottleneck are balanced and controlled, and a is recommended to be 0.65.

[0099] In one preferred embodiment of the present embodiment, the determination method of all the weight coefficients (such as ωr1, ωa1, wi, etc.) is as follows: a combination of objective weighting method and subjective weighting method is adopted.

[0100] Specifically, a plurality of sets of actual response data of each project in historical flood events are collected, and a set of objective weights w j e is calculated based on the data by using an entropy weight method. Meanwhile, field experts are organized to compare and score the importance of each index, a judgment matrix is constructed by using an analytic hierarchy process (AHP), and another set of subjective weights w j a is calculated. The final combined weight calculation formula is: ω j = β × w j e + (1-β) × w j a ; wherein β is a preference coefficient, for example, β = 0.6 indicates that more emphasis is placed on the law reflected by objective data.

[0101] In the present embodiment, by using multi-level index aggregation and dynamic coupling, a resilience quantification system covering single projects, subsystems, and large systems is constructed, and the problem of difficult unified quantification of resilience characteristics of complex water engineering groups under multi-scale and multi-coupling relationships is solved. Based on the four-dimensional indexes of resistance, absorption, recovery, and adaptability, each level is dynamically aggregated to realize accurate characterization of resilience performance from the local to the global. The scheme can be dynamically adjusted with changes in engineering structure and evolution of running state, and is more suitable for application scenarios simulating system evolution and resilience strategy adjustment under extreme working conditions, and provides quantitative basis for risk prevention and control and resilience improvement of water engineering groups.

[0102] Embodiment Four, the present embodiment is used to describe the generation process of the resilience performance boundary, to illustrate how the performance lower limit as the risk assessment benchmark is established.

[0103] In the present embodiment, the generation process of the resilience performance boundary specifically includes the following steps:

[0104] Step one, constructing an extreme flood scenario set S extremeThe system reads a local historical flood database (e.g., containing all recorded flood process data in the past 50 years) and design standard flood data (e.g., 100-year, 1000-year design flood hydrograph published by the water department) in the region. By screening and combining these data, an extreme flood scenario set is constructed, which can represent the most severe challenges that the water project group may encounter. For example, the set can include: the flood event with the largest flood peak flow actually occurred in history, the flood event with the largest total flood volume in history, and the possible maximum flood (PMF) hydrograph derived from PMP (possible maximum precipitation), etc.

[0105] Step two, calculate the disturbance effect of extreme scenarios on normalized characteristic factors. Specifically, the extreme flood scenario set S extreme is constructed in step one as an external disturbance input. The system calculates the impact of these extreme scenarios on each normalized characteristic factor (such as reservoir capacity utilization rate C extreme , safety elevation margin H u,norm , etc.) in the system through a disturbance response function Φ(l)(S s,norm ). For example, for the safety elevation margin, an extreme flood scenario will reduce its margin by raising the reservoir water level. The updated characteristic factor value is calculated as x m ′=x m +Δx m,extreme ; where x m is the initial value of the characteristic factor; Δx m,extreme is the change of the characteristic factor under the extreme scenario; x m ′ is the disturbed characteristic factor value.

[0106] Step three, calculate the resilience performance value under extreme scenarios. Substitute all the characteristic factor values x m ′ disturbed by extreme scenarios obtained in step two into the multi-dimensional resilience performance index system constructed in example three. According to the calculation method defined in the system, the performance index values of the system under each extreme scenario are calculated at the single project-subsystem-system three levels, and the resistance-absorption-recovery-adaptation four dimensions.

[0107] Step four, determine the final resilience performance boundary R boundary . In this embodiment, in order to adopt a risk-averse conservative strategy, the resilience performance boundary is defined as the least favorable performance that the system can achieve under all possible extreme scenarios. The specific implementation is to take the minimum value of all performance values calculated in step three under each performance dimension. That is, R (l) boundary,k =mins∈S extreme {Pk,s(l)}; wherein R (l) boundary,k is the resilience performance boundary value of the kth performance dimension (e.g. absorption) at the lth level (e.g. subsystem level); S extreme is the set of extreme flood scenarios; P k,s(l) is the corresponding performance value calculated under the extreme scenario s. In this way, a set of performance baselines covering all levels and all dimensions is generated, which constitutes the resilience performance boundary. It provides a clear, physically meaningful reference benchmark for subsequent risk assessment.

[0108] Embodiment Five, this embodiment is used to describe how to construct a composite perturbation and assess the current resilience performance state in the presence of uncertainty.

[0109] In this embodiment, the assessment of a current resilience performance state is carried out under a composite perturbation, where the generation of the composite perturbation includes the following steps:

[0110] Step One, based on the flood forecast data and its statistical characteristics of historical forecast errors, a set of perturbation samples representing the forecast uncertainty is generated by random sampling. The specific implementation is that the system first reads the current latest flood forecast data, for example, the incoming flow forecast hydrograph Q forecast for the next 72 hours. At the same time, the system analyzes the historical forecast data and the measured data to extract the statistical characteristics of the forecast error. For example, it is found that the historical forecast error approximately obeys a normal distribution with mean μ err = -0.05 and standard deviation σ err = 0.15, i.e. ε Q ~ N(-0.05, 0.152), which indicates that there is a systematic tendency to underestimate the forecast. In order to efficiently generate perturbation samples that reflect this uncertainty, this embodiment preferably uses an improved Latin hypercube sampling method. This method divides the probability distribution space of the error into M (e.g. M = 100) equal probability intervals, and then independently randomly samples an error sample value ε i in each interval. Applying these error samples to the forecast values, M sets of perturbed flow sequences are generated: Q perturbed,i = Q forecast × (1 + εi), where i = 1...M. This set constitutes the perturbation samples representing the forecast uncertainty.

[0111] Step Two, from the historical and design flood data, the dynamic time warping algorithm is used to match the largest similar historical flood scenario most similar to the current forecast flood process.

[0112] Since a single forecast uncertainty often cannot fully capture the complex spatiotemporal evolution characteristics of a flood (such as the steepness of the flood peak, the multi-peak characteristics, etc.), the introduction of a similar real-occurred flood scenario can be an effective supplement. The specific implementation is as follows: first, for the current forecast flood process and each flood process in the historical database, a three-dimensional feature vector V1=(Q peak ,W total,γ ) is extracted respectively; wherein Q peak is the flood peak flow; W total is the total flood volume; and γ is a coefficient representing the shape of the process line (such as the ratio of the peak occurrence time to the flood duration). Then, the dynamic time warping (DTW) algorithm is used to calculate the similarity D DTW between the current forecast flood process line and each historical flood process line, which is especially suitable for comparing time series of different lengths or with time offset. The comprehensive similarity is calculated by a weighted formula: S total =w p ×exp(-∣Q p,h -Q p,c ∣ / Q p,c )+w w ×exp(-∣W t,h -W t,c ∣ / W t,c )+w d ×exp(-D DTW / D0);wherein the subscripts h and c represent history and current respectively; w p ,w w ,w d are weights; and D0 is a normalization constant. The historical flood with the maximum S total value is selected as the maximum similar historical flood scenario.

[0113] Step three, combine the perturbation samples and the maximum similar historical flood scenario to construct a composite perturbation. Specifically, the system combines the uncertainty perturbation samples generated in step one (representing the uncertainty in the forecast value) and the maximum similar historical flood scenario matched in step two (representing the uncertainty in the flood morphology). One preferred combination method is to superimpose the N groups of perturbation samples on the process line of the maximum similar historical flood scenario, thereby generating N groups of composite perturbation scenarios. Each group of composite perturbation scenarios contains both the randomness of the forecast error and the evolution morphology of the real flood.

[0114] On this basis, the N sets of complex disturbance scenarios are taken as inputs, and are applied to the water engineering group model one by one. For each set of scenarios, the resilience performance index of the system at each level and dimension is calculated according to the method of Example Three, so as to obtain N sets of current resilience performance state values. The set composed of the N sets of values is the current resilience performance state obtained through the final evaluation, which fully considers the uncertainty.

[0115] Example Six, this embodiment describes how to identify the resilience failure node, and how to track and build the failure transmission chain across levels.

[0116] In this embodiment, before determining the regulation potential of the resilience failure project, the following steps are further included:

[0117] Step One, compare the current resilience performance state with the resilience performance boundary, and calculate a performance gap. Specifically, the system compares the N sets of current resilience performance state values P current,i (where i = 1...N) obtained through the evaluation in Example Five with the resilience performance boundary R boundary generated in Example Four. The performance gap is defined as ΔP i = P current,i - R boundary . When ΔP i ≤ 0, it indicates that the performance of the project or system has dropped below the safety bottom line under the i th disturbance sample, and there is a failure tendency.

[0118] Step Two, count the failure risk probability that the performance gap is lower than the preset value under the uncertainty disturbance. Specifically, count the number of times that the performance gap ΔP i ≤ 0 in N simulations, and record it as N fail . The failure risk probability is calculated as P fail = N fail / N. The probability intuitively reflects the possibility of the occurrence of resilience failure of the project under the consideration of uncertainty.

[0119] Step Three, when the failure risk probability and the absolute mean value of the performance gap both exceed the corresponding probability threshold and performance gap threshold, respectively, the corresponding project is determined to be a resilience failure node, so as to form a resilience failure node set. The design of the double-threshold criterion is to avoid false positives due to accidental and slight performance fluctuations.

[0120] Specifically, a probability threshold θ p and a performance gap threshold δ p are set. The value of θ p can be determined according to the risk tolerance, for example, for a critical project, θ p = 5% can be taken, indicating that the failure probability exceeding 5% needs to be vigilant. The value of δ p can be determined according to the physical meaning of the performance index, for example, it is set to 10% of the performance boundary value, i.e. δ p = 0.1 × R boundaryAt the same time, the absolute mean value of performance gap ΔP avg =∑i∣ΔPi≤0∣ΔPi∣ / N fail is calculated in all failed samples, which reflects the average severity of failure. Only when P fail >θp and ΔP avg >δ p are met at the same time, the system will determine the project or subsystem as a resilience failure node. All determined nodes together constitute a resilience failure node set.

[0121] Further, after forming the resilience failure node set, the method further comprises: based on the dual-track dynamic hierarchical structure, using a top-down tracking mechanism, starting from the large system level, layer by layer locating the failed subsystems and single projects in the resilience failure node set, to construct a cross-level failure chain representing the failure transmission relationship. The implementation of this step is as follows:

[0122] Identify the large system layer failure trigger point: First check whether the comprehensive performance index of the large system level is determined to be failed. If so, mark the failed dimension (such as large system absorption failure), and record it as the starting point of the failure chain.

[0123] Locate the failure subsystem layer by layer: For each failure dimension of the large system layer, traverse all the subsystems under it. Using the dual-track hierarchical structure constructed in Example Two, check which subsystems are also determined to be failed in the same dimension. Calculate the contribution of each subsystem to the failure of the large system, and select the subsystems with contribution exceeding the threshold as the intermediate links of the failure chain.

[0124] Precisely locate the single project failure node: Within each located failure subsystem, further traverse all single projects it contains. Find out the single project nodes that are failed in the corresponding dimension. At the same time, use the comprehensive hydraulic coupling strength H ij (t) calculated in Example Two as a measure of failure transmission strength, identify the project inside the subsystem that is the main failure source (i.e. its transmission strength to other projects is the highest).

[0125] Construct a structured failure chain mapping: Integrate the failure node information of the above three levels to construct a cross-level failure chain like large system absorption failure-tributary A1 subsystem absorption failure-A1-2 reservoir absorption failure. At the same time, the type of failure chain (such as series, parallel, cascade) can be marked, and its key degree can be calculated to provide priority ranking for subsequent control decisions. In this way, the abstract system failure problem is transformed into a clear and traceable structured transmission path.

[0126] Embodiment Seven: This embodiment describes how a hierarchical resilience regulation response graph is constructed as a prerequisite for path search.

[0127] In this embodiment, in order to systematically search and evaluate regulation paths, a hierarchical resilience regulation response graph needs to be constructed first. The construction process of this graph is based on the failure engineering identified in the cross-level failure chain in Embodiment Six, and is constructed according to the characteristics of different levels. Specifically, it includes the following steps:

[0128] Step One: Construct the regulation response graph of the single engineering layer. In this layer, the nodes in the graph are defined as the three core performance dimension states of the single engineering, namely the resistance state, the absorption state, and the recovery state. The edges in the graph represent the coupling paths between these performance dimensions for regulation and conversion. For example, an action that represents improving the resistance performance through pre-discharge reservoir (a regulation action) can be modeled as an edge from the current resistance state node to a higher resistance state node. This edge will be assigned corresponding attributes, such as required time, resource consumption, and possible negative impact on other performance dimensions (such as absorption).

[0129] Step Two: Construct the regulation response graph of the subsystem layer. In this layer, the nodes in the graph are defined as the comprehensive performance states of each single engineering within the subsystem. The edges in the graph represent the collaborative regulation links between the engineering. The establishment of an edge needs to meet two basic conditions: spatial connectivity and timeliness. Specifically, only when two engineering have effective hydraulic connection in the double-track hierarchical structure constructed in Embodiment Two, and the regulation action of one engineering can have effective influence on the other within the required time (for example, within the critical time window of flood evolution), a collaborative regulation edge can be established between the two engineering nodes. For example, the water release action of upstream reservoir A can reach midstream water gate B within 6 hours and effectively raise its water level, so a edge can be established between A and B.

[0130] Step Three: Construct the regulation response graph of the large system layer. In this layer, the nodes in the graph usually represent the performance states of all key engineering in the whole basin. The edges in the graph mainly represent the long-distance regulation trigger links across regions and subsystems. Such links usually do not rely on direct hydraulic connection, but are achieved through global scheduling rules and instructions. For example, when it is monitored that the water level of a key section of the main stream is about to exceed the warning level, the dispatching center can issue an instruction to require multiple reservoir groups located in different tributaries to reduce the discharge at the same time. This scheduling decision based on global information forms a cross-regional regulation trigger edge between these scheduled reservoir group nodes.

[0131] Step four, assign regulation cost attributes to all edges. In order to be able to quantitatively evaluate in the subsequent path search, it is necessary to assign a set of regulation cost attributes to each edge in the above-mentioned all hierarchical response graph. In this embodiment, these attributes at least include: response time (the time from the start of the action to the expected effect), regulation efficiency (the amount of performance improvement that can be brought by unit resource consumption), and path complexity (the complexity of regulation instruction transmission and execution). These attribute values can be obtained by hydraulic model calculation, historical data statistics or expert evaluation, etc. Through the above steps, finally output a layered resilience regulation response graph containing multiple levels, which can comprehensively reflect various regulation possibilities and their corresponding costs in the system. The graph is the data basis for subsequent regulation potential determination.

[0132] Embodiment eight, this embodiment describes how to determine the regulation potential of resilience failure engineering, and how to classify it.

[0133] In this embodiment, after obtaining the set of resilience failure nodes, the process of determining the regulation potential of resilience failure engineering further includes:

[0134] Step one, based on the current resilience performance state and the resilience performance boundary, filter out the nodes with performance surplus to form a regulation source set.

[0135] Specifically, the system traverses all engineering nodes that have not been determined to be failed. Calculate its performance surplus degree ΔP surplus =P current -R boundary . Filter out those nodes whose performance surplus degree is greater than a preset surplus threshold (for example, θ surplus =0.2, indicating that the performance is at least 20% higher than the boundary), as potential regulation sources that can support other failed nodes. These regulation sources together constitute the regulation source candidate set S source .

[0136] Step two, for any resilience failure node in the set of resilience failure nodes, search for a closed-loop cooperative path connecting the failure node and at least one regulation source in the dual-track dynamic hierarchical structure.

[0137] In this embodiment, a closed-loop synergistic path refers to a control sequence starting from a failed node, passing through one or more control sources, and finally enabling the performance state of the failed node to return to the effective state interval above the resilience performance boundary. The specific implementation is to traverse the paths in the layered resilience control response graph constructed in Embodiment Seven from each resilience failed node using an improved depth-first search algorithm (DFS). The search process sets constraints such as the path length cannot exceed 5 hops (i.e., the number of engineering involved cannot exceed 5), the total response time of the path cannot exceed 6 hours, and the path must contain at least one node from the control source candidate set S source All paths that meet these constraints are recorded to form a candidate path set.

[0138] Step three, and calculate a feasibility score of the closed-loop synergistic path, which represents the control potential of the resilience failed node. In this embodiment, the method for calculating the feasibility score of the closed-loop synergistic path is to comprehensively evaluate at least two of the following performance indicators of the closed-loop synergistic path: path length, total response time, control efficiency, and net performance improvement; and generate the feasibility score based on the evaluation results of the at least two performance indicators.

[0139] A preferred feasibility score function is:

[0140] F(P i )=λ1×(1 / ∣P i ∣)+λ2×(1 / T total )+λ3×C efficiency (P i )+λ4×ΔΦ net (P i );

[0141] Where F(P i ) is the feasibility score of path P i ; ∣P i ∣ is the path length (number of nodes), and its reciprocal represents the path length penalty term; T total is the sum of the response times of all edges in the path, and its reciprocal represents the timeliness score; C efficiency (Pi) is the control efficiency, representing the ratio of the total performance transfer to the total capacity consumption of the control source; ΔΦ net (Pi) is the net performance improvement, representing the performance improvement obtained by the failed node minus the transmission loss existing in the path; λ1, λ2, λ3, λ4 are weight coefficients of each indicator, which can be set to λ1=0.2, λ2=0.3, λ3=0.3, λ4=0.2 according to the decision preference, and the sum is 1.

[0142] Step four, divide it into two types of failure engineering set, controllable and uncontrollable. The specific implementation is to compare the feasibility score with a preset feasibility score threshold; if the feasibility score is not less than the feasibility score threshold, the corresponding resilience failure node is divided into the controllable failure engineering set; otherwise, the resilience failure node is divided into the uncontrollable failure engineering set.

[0143] First, for all candidate paths of a failure node, arrange them in descending order according to their feasibility scores F(P i ), and select the path with the highest score as the optimal closed loop path. Then, compare the score F optimal of the optimal path with a preset feasibility score threshold θ loop . The threshold can be determined according to historical experience or simulation analysis, for example, θ loop =0.6.

[0144] If F optimal ≥θ loop , it is considered that there is a high-quality controllable path to restore the performance of the node, so the resilience failure node is divided into the controllable failure engineering set. If there is no any feasible closed loop path, or the highest score of all paths is also less than θ loop , it is considered that the performance cannot be effectively restored by internal coordination under the current system state, so it is divided into the uncontrollable failure engineering set. Alternatively, multiple thresholds can be set, for example, when the score is in a certain intermediate interval (such as 0.4≤F optimal <0.6), it can be divided into a weak controllable failure engineering set, indicating that it needs external stronger auxiliary support.

[0145] In this embodiment, by analyzing the multi-dimensional resilience capability state of the engineering node and identifying whether there is a closed loop self-help path in the control response graph, it is determined whether it has self-organization control potential, and then the failure engineering is divided into two types of controllable and uncontrollable. Unlike the traditional method which relies on fixed threshold and single path judgment, this scheme is based on the dynamic response characteristics of network structure, and identifies the closed loop by integrating multiple source control paths, which is more suitable for the multi-dimensional coupling characteristics in actual resilience recovery, and is suitable for multi-level, multi-region cross-dimensional control potential identification.

[0146] Embodiment nine provides a clear and explicit mathematical optimization model, that is, to define the objective function and constraint conditions that the algorithm needs to solve.

[0147] In this embodiment, before performing a specific control optimization algorithm, a multi-objective resilience balance scheduling model needs to be constructed first. The construction of this model is to read the controllable failure engineering set, the uncontrollable failure engineering set and the matching auxiliary control engineering set divided in embodiment eight, and then construct the optimization objective function hierarchically.

[0148] Single engineering layer optimization objective function: The objective of this level focuses on risk control and efficiency of a single engineering itself.

[0149] As a single engineering as a basic unit of resilience, its structure is single and its function is direct, and it is easy to be affected by lag response and boundary events, the key is to regulate the response speed and prevent single-point failure. The focus of optimization is the local regulation performance and critical state prevention and control capability. Therefore, the single engineering level aims to minimize failure risk, minimize response lag, and maximize resilience surplus, prioritizes dynamic stability and timeliness of single-point engineering regulation, and its objective function can be constructed as:

[0150] F 单 =min{a1×P fail +a2×T d -a3×ΔP res};

[0151] Where F 单 is the comprehensive optimization objective of the single engineering layer, aiming to minimize this value; P fail is the expected failure risk probability of the engineering after regulation; T d is the response lag time of the regulation action; ΔP res is the resilience surplus or redundancy of the engineering after regulation, the goal is to maximize this item, so the front is negative; a1, a2, a3 are the weight coefficients of each sub-objective.

[0152] Subsystem layer optimization objective function: The objective of this level focuses on the overall coordination and cost-effectiveness of the subsystem.

[0153] A subsystem is a collection of several single engineering units with certain autonomous regulation capabilities, but its regulation effect depends on the resilience coordination of internal single engineering units. It is necessary to focus on optimizing the internal regulation balance and synergy benefit of this level to prevent chain failure spread. Therefore, the subsystem level aims to minimize resilience imbalance, minimize redundant regulation cost, and maximize synergy efficiency to improve overall response consistency and resource allocation efficiency, and its objective function can be constructed as:

[0154] F 子 =min{b1×σ(P 子 )+b2×C cost -b3×E coop};

[0155] Where F 子 is the comprehensive optimization objective of the subsystem layer; σ(P 子 ) is the standard deviation of the resilience performance of each engineering in the subsystem, used to measure resilience imbalance, the goal is to minimize it; C costTotal cost paid for performing the coordinated regulation (e.g. water loss, power generation loss, etc.);E coop Synergistic effect, i.e. the additional performance gain brought by the coordinated regulation, the goal is to maximize this item; b1, b2, b3 are weight coefficients.

[0156] Optimization objective function at the large system level: the goal of this level focuses on the global risk exposure and the overall system resilience.

[0157] The large system is a whole network of cross-regional and cross-level linkage, with network linkage and cross-regional coupling attributes, but its structure is complex, the risk chain is long, and the regulation conflict is prone to occur. The coordination ability between subsystems and the overall comprehensive resilience level should be focused on to avoid local failure leading to global imbalance. Therefore, the large system level takes the minimization of system risk exposure, the minimization of regulation conflict, and the maximization of system resilience improvement as the goal, and emphasizes the improvement of overall situation and global coordination ability. Its objective function can be constructed as: 大 F exposure = min{c1×R conflict +c2×C 大 -c3×ΔP 大};

[0158] Where F exposure is the comprehensive optimization objective of the large system level; R conflict is the comprehensive risk exposure of the whole system after regulation; C 大 is the conflict degree of regulation actions between different subsystems or projects, the goal is to minimize the conflict; ΔP boundary is the total improvement of the comprehensive resilience performance of the whole large system, the goal is to maximize this item; c1, c2, c3 are weight coefficients.

[0159] In addition, the whole optimization model also needs to be solved under a series of constraint conditions. These constraint conditions at least include:

[0160] Water balance constraint: the inflow, outflow, and inter- interval inflow of each project and river reach must satisfy the water balance equation.

[0161] Regulation boundary constraint: all regulation actions (such as discharge flow, reservoir water level, etc.) must be within the safe operation boundary of the project design. For example, the reservoir water level cannot be higher than the high water level for flood control, nor can it be lower than the dead water level; the discharge flow cannot exceed the maximum discharge capacity of the flood discharge facility.

[0162] Resilience effectiveness constraint: the regulation scheme must ensure that the resilience performance value of the failed node is at least restored to above the resilience performance boundary R boundary .

[0163] By constructing the above complete mathematical model containing hierarchical objective function and various constraint conditions, clear and calculable optimization target and feasible region are provided for subsequent solution by using multi-objective reinforcement learning, collaborative co-evolution and other specific optimization algorithms.

[0164] Embodiment Ten, Multi-objective Self-regulation Optimization of Regulatable Failure Engineering.

[0165] This embodiment elaborates in detail the multi-objective reinforcement learning guided hybrid gradient search algorithm (MORL-HG) adopted for the regulatable failure engineering set.

[0166] In this embodiment, the multi-objective reinforcement learning guided hybrid gradient search algorithm is used to optimize the regulatable failure engineering set. The process specifically includes:

[0167] Step One, based on the physical constraints of the regulatable failure engineering set, a safe exploration space for limiting the range of regulation actions is established.

[0168] The specific implementation is that the system first reads the current running state of the regulatable failure engineering (such as water level Z, discharge flow Q, reservoir capacity V) and its own physical properties. Then, based on hydraulics and engineering safety criteria, a series of inequality constraints are established to define the safety boundary of regulation actions.

[0169] For example, these constraints can include:

[0170] Discharge capacity constraint: Q out ≤Kv×L×H 1.5 ; where Q out is the discharge flow action; Kv is the flow coefficient; L is the spillway width; H is the weir head. This constraint ensures that the scheduling instruction does not exceed the maximum discharge capacity of the engineering.

[0171] Water level change rate constraint: |dZ / dt|≤dZ max ; where dZ / dt is the water level change rate caused by regulation; dZ max is the maximum allowable water level change rate (e.g. 0.5 meters per hour) to prevent bank slope instability caused by rapid water level fluctuations.

[0172] Reservoir capacity constraint: V dead ≤V(t+1)≤V flood ; where V(t+1) is the expected reservoir capacity after executing the regulation action; V dead and V flood are the dead storage and flood control limited storage respectively.

[0173] The set of all continuous regulation actions (such as specific values of discharge flow) that satisfy the above physical constraints collectively constitutes the safe exploration space A safe (ss ), wherein s s represents the current state. The principle of this step is to fundamentally avoid the optimization algorithm from generating dangerous instructions that cause damage to the engineering entity in the exploration process by pre-defining an absolutely safe action space.

[0174] Step two, an initial control action is generated by a multi-objective reinforcement learning agent within the safe exploration space.

[0175] The specific implementation is to first initialize a multi-objective reinforcement learning (MORL) agent. The state space S of the agent can be a multi-dimensional feature vector containing the current water level, inflow, outflow, and four-dimensional resilience performance indicators, etc. The action space A a is set as a continuous control action.

[0176] Then, a multi-objective reward function is designed for the agent to guide its learning direction. The reward function can be designed as: r total = r risk + r performance + r robust ; wherein, r risk = -w1 x max(0, P fail - P threshold ), which is a penalty term that gives a negative reward when the expected failure risk exceeds the threshold;

[0177] r performance = w2 x (P current - P previous ), which is a performance reward term that encourages the agent to take actions that can improve resilience performance;

[0178] r robust = w3 x (1 - sigma(P scenarios ) / mu(P scenarios )), which is a robustness reward term that encourages the agent to find actions that perform stably under different uncertainty scenarios; w1, w2, w3 are weights, which can be set to 0.4, 0.4, 0.2, for example.

[0179] During training and decision-making, the agent can be pre-trained using historical control data to obtain an initial strategy pi0. At the decision-making moment, the agent outputs an initial control action a with good comprehensive performance based on the current state and its learned strategy.

[0180] Step three, a gradient search algorithm is used to fine-tune the initial control action while ensuring that the adjusted action does not exceed the boundary of the safe exploration space, to obtain an optimized control action sequence.

[0181] Since reinforcement learning is good at finding a good direction at the macro-strategy level, but may not be efficient in fine-tuning specific numerical values, gradient search is introduced to make local fine-tuning. The specific implementation is that after the agent generates an initial control action a, a gradient search is performed in a neighborhood of the action a to find a better action value. For example, by calculating the gradient of the immediate performance function J with respect to the action a, iterative updates are performed:

[0182] a' = a + a' x VJ(a); where a is the fine-tuned action; a' is the learning step size, for example, 0.01.

[0183] After each iterative update, it must be checked whether the adjusted action a' is still within the safe exploration space A safe (s) constructed in step one. If a' exceeds the boundary, it needs to be projected back to the nearest point on the safe boundary, so that the final output action is always absolutely safe.

[0184] By combining the global exploration ability of MORL and the local optimization ability of gradient search, a series of optimized control actions can be efficiently found under the premise of safety, forming the final self-regulation optimization scheme.

[0185] Embodiment eleven, this embodiment describes a collaborative control process for the non-controllable failure engineering set, including the matching mechanism of auxiliary engineering and the specific implementation of collaborative-competitive co-evolution algorithm (CCEA-LF).

[0186] In this embodiment, for the non-controllable failure engineering set, a matching auxiliary control engineering set is combined, and the collaborative-competitive co-evolution algorithm is used for optimization. The process specifically includes:

[0187] Step one, match an auxiliary control engineering set for the non-controllable failure engineering set.

[0188] First, a resilience demand vector is established for the failure engineering in the non-controllable failure engineering set, and a resilience supply vector is established for the candidate engineering in the system. Specifically, for a non-controllable failure engineering, its resilience demand vector V need can be composed of its gap values in each failure performance dimension, for example, V need =[ΔPr, ΔPa, ΔPh]. For other engineering in the system that has surplus capacity (i.e. potential auxiliary engineering), its resilience supply vector V supply can be composed of the performance surplus it can provide in each dimension.

[0189] Second, the matching degree between the resilience demand vector and the resilience supply vector is calculated, where the matching degree is determined based on the cosine similarity between the vectors and the hydraulic propagation distance between the failure engineering and the candidate engineering. A preferred matching degree calculation formula is: M ij = cos(V need,i , V supply,j ) x exp(-d ij / d0); where M ij is the matching degree of the candidate engineering j to the failure engineering i; cos(*) calculates the directional similarity of the demand vector and the supply vector, and the more consistent the direction is, the more suitable the supply is; the exp(*) term is a decay term based on the hydraulic propagation distance d ij , which indicates that the farther the distance is, the worse the support effect is, and d0 is a characteristic distance constant.

[0190] Finally, the candidate engineering set with the highest matching degree is selected to form the auxiliary control engineering set. The system will calculate the matching degree of all candidate engineering to each failure engineering, and select the Top-K engineering with the highest matching degree to form the exclusive auxiliary control engineering set.

[0191] Step two, for the set of uncontrolled failure engineering, a matching auxiliary control engineering set is combined to optimize using the co-evolution algorithm.

[0192] In this step, an independent evolution population is established for each auxiliary engineering in the auxiliary control engineering set, where each individual in the population encodes a support decision. For example, the size of a population can be set to 50, and each individual can be a string of codes that defines the support action of the engineering (such as the increase in the outflow, the scheduling time, etc.).

[0193] Then, in the iteration process, individuals are combined from the evolution populations to form a coordinated control scheme, and the evolution populations are evaluated and updated according to the comprehensive influence of the coordinated control scheme on the set of uncontrolled failure engineering to obtain the final coordinated control optimization scheme.

[0194] The specific implementation is that in each round of iteration:

[0195] Combination and evaluation: select the optimal individual (i.e. the optimal support decision) from the population of each auxiliary engineering, combine them to form a complete coordinated control scheme. Then evaluate the comprehensive performance improvement effect of the coordinated scheme on the target failure engineering through model simulation.

[0196] Fitness update: Based on the evaluation results, update the fitness of each participating individual. The fitness function here takes into account multi-level feedback, such as considering both the individual's direct contribution to failure engineering (local contribution) and the impact of the collaborative scheme on the overall subsystem or system resilience (system contribution), and dynamically adjusts the contribution weight.

[0197] Evolutionary operation: According to the updated individual fitness, each population independently performs genetic algorithm operations such as selection (e.g. tournament selection), crossover (e.g. arithmetic crossover) and mutation (e.g. Gaussian mutation) to generate a new generation of population.

[0198] Competition and alliance: In the case of limited resources (e.g. limited total available support water), introduce a competition mechanism, and individuals or collaborative combinations with higher fitness will have priority to obtain regulatory resources. At the same time, the system will dynamically adjust the alliance according to the collaborative effect, and move the engineering with lower contribution out of the auxiliary set.

[0199] This iterative process continues until the maximum number of iterations or the scheme converges. Finally, the algorithm outputs a set of collaborative regulation optimization schemes that include specific support actions, execution timing and collaborative paths for each auxiliary engineering.

[0200] In this embodiment, after obtaining the self-regulation optimization scheme set for the controllable engineering and the collaborative regulation optimization scheme set for the uncontrollable engineering respectively, the system also needs to perform the following steps to output the final optimal scheduling scheme set:

[0201] Step one, spatiotemporal coordination of scheme set. Specifically, the system places all the regulation actions in the two types of scheme sets in a unified spatiotemporal coordinate for checking to identify and eliminate potential conflicts.

[0202] For example, a collaborative regulation scheme may require upstream reservoir A to release water at time T to support downstream failure engineering C d , while a self-regulation scheme may require midstream water gate B between A and C d to reduce discharge at time T to improve its own resilience. These two actions may conflict in terms of water power. At this time, the system needs to adjust the execution time or action amplitude of one of the schemes according to the preset priority rules (e.g. the rule priority of ensuring key protection objects is higher) or through re-simulation evaluation to ensure the global feasibility of the integrated scheme.

[0203] Step two, non-inferior solution screening based on Pareto dominance. After coordination, all feasible solutions constitute a large pool of candidate solutions. Since the optimization objectives of the present application are multi-dimensional (e.g. risk, cost, efficiency in Example Nine), these objectives often conflict with each other, and there is no perfect solution that is optimal in all objectives. Therefore, the present embodiment adopts the Pareto optimality theory to screen the solutions.

[0204] Specifically, the system compares all solutions in the pool of candidate solutions two by two. If solution X is better than or equal to solution Y in all optimization objectives, and strictly better than solution Y in at least one objective, then solution X is said to dominate solution Y. After traversing all solutions and removing all solutions dominated by other solutions, the remaining set of solutions is the non-inferior solution set, also known as the Pareto optimal frontier. Each solution on this frontier is an optimal solution, because to improve any one of these solutions in any objective will necessarily come at the expense of at least another objective.

[0205] Step three, representative solution selection and solution set output. Finally, the system presents the Pareto optimal frontier obtained in step two to the decision maker. Optionally, the system can automatically select several representative solutions from the frontier according to pre-set decision preferences or balancing principles.

[0206] For example, a solution that focuses on risk aversion (i.e. the point on the frontier that minimizes risk), a solution that focuses on cost-effectiveness (i.e. the point on the frontier that minimizes cost), and a balanced solution (i.e. the inflection point on the frontier that is closest to the origin) can be provided. The final output of the optimal scheduling solution set is a structured data set that clearly lists the specific control parameters (e.g. the discharge hydrograph of each project) contained in each representative solution, the execution timing (the start and end times of each action), and the expected resilience improvement effect and corresponding cost.

[0207] Example Thirteen, assuming there is a simplified water project system consisting of an upstream reservoir A, a midstream sluice B, and a downstream control section C distributed along a river channel in sequence. This case will calculate the hydraulic coupling strength between each project in the system at a certain time t, and determine the weight coefficients for fusing the physical level and dynamic level.

[0208] Real-time hydrological data (t time)

[0209] Reservoir water level H A = 150.0 meters, H B = 85.0 meters.

[0210] Reservoir flow rate Q A = 1200 m 3 / s, QB = 1250 m 3 / s (interval inflow).

[0211] Physical and model parameters, including:

[0212] Inter-engine river distance LAB=30000 meters. Gravity acceleration g=9.8 m / s 2 . Interval average water depth havg=10.0 meters. Characteristic time constant T0=10800 seconds (3 hours). Sensitivity coefficient K=0.1. Tensor combination weight w1=0.4, w2=0.3, w3=0.3.

[0213] Step one: Calculate the hydraulic propagation time τ ij (t), according to the formula τ ij (t)=L ij / (g×h avg ) 0.5 :

[0214] The propagation time τ AB from the upstream reservoir A to the midstream sluice B is 30000 / (9.8x10.0) 0.5 =30000 / 9.899≈3030.5 seconds.

[0215] Step two: Calculate the partial derivative of each relationship matrix

[0216] In this case, the finite difference method is used to approximate the partial derivative. This requires a perturbation simulation through the hydraulic model. Suppose we apply a small perturbation to the working condition of the upstream reservoir A, and then observe the response of the midstream sluice B.

[0217] Suppose at the next moment, the discharge of the upstream reservoir A increases by a small amount ΔQ A =10 m 3 / s, causing a small change in its water level ΔH A =0.02 meters.

[0218] Through the hydraulic model (such as HEC-RAS or similar one-dimensional model) simulation, it is observed that after τ AB time, the disturbance propagates to the midstream sluice B, causing an increase in B's discharge ΔQ B =9.5 m 3 / s, and an increase in water level ΔH B =0.15 meters.

[0219] Based on this, the partial derivative approximation is calculated: dH B / dH A ≈ΔH B / ΔH A =0.15 / 0.02=7.5. dQB / dQ A ≈ΔQ B / ΔQ A =9.5 / 10=0.95。dQ B / dH A ≈ΔQ B / ΔH A =9.5 / 0.02=475。

[0220] Step three: Calculate the three types of coupling relationship matrix

[0221] According to the formula and the results of step one and two, calculate the influence element value of A to B in each relationship matrix:

[0222] Water level propagation matrix element H HH,AB :

[0223] H HH,AB =(dH B / dH A )×exp(-τ AB / T0)=7.5×exp(-3030.5 / 10800)=7.5×exp(-0.2806)=7.5×0.7553≈5.665。

[0224] Flow propagation matrix element H QQ,AB :

[0225] H QQ,AB =(dQ B / dQ A )×exp(-τ AB / T0)=0.95×0.7553≈0.718。

[0226] Water level-flow coupling matrix element H HQ,AB :

[0227] H HQ,AB =(dQ B / dH A )×exp(-τ AB / T0)=475×0.7553≈358.768。

[0228] Step four: Synthesize the comprehensive hydraulic coupling strength H ij (t)

[0229] According to the formula H ij (t)=w1×∣H HH,ij ∣+w2×∣H QQ,ij ∣+w3×∣H HQ,ij ∣:

[0230] Calculate the comprehensive hydraulic coupling strength HAB (t): H AB (t) = 0.4 x |5.665| + 0.3 x |0.718| + 0.3 x |358.768| = 2.266 + 0.2154 + 107.6304 = 110.11.

[0231] Assuming other elements are calculated by similar methods (e.g. B has less impact on A, self impact on self is not counted or 0, A's impact on C needs to be calculated in segments, etc.), a simplified 2x2 hydraulic coupling tensor matrix (only considering A and B) is finally obtained: H(t) = (0, 110.11; H BA , 0).

[0232] Where H BA is the reverse impact of B on A, usually small, in this example, set to 5.0.

[0233] Step five: calculate the hierarchical fusion weight coefficient a(t).

[0234] Calculate the spatial variance V ariance (H ij ) of all element values in the hydraulic coupling tensor matrix:

[0235] The elements in the matrix are {0, 110.11, 5.0, 0}. The mean H mean = (0 + 110.11 + 5.0 + 0) / 4 = 28.7775.

[0236] Where the variance V ariance (H ij ) = Σ (H ij - H mean ) 2 / (N2-1) = [(0-28.78) 2 + (110.11-28.78) 2 + (5.0-28.78) 2 + (0-28.78) 2 ] / (4-1) = 2945.57.

[0237] Calculate the hierarchical fusion weight coefficient a(t):

[0238] a(t) = 1 / (1 + K x V ariance (H ij )) = 1 / (1 + 0.1 x 2945.57) = 1 / (1 + 294.557) = 1 / 295.557 = 0.00338.

[0239] At a specific time t in this case, the calculated comprehensive hydraulic coupling strength H AB(t) ≈ 110.11, which is a higher value, indicating that at this time the operation state of the upstream reservoir A has a very strong hydraulic influence on the midstream sluice B.

[0240] The calculated hierarchical fusion weight coefficient a(t) ≈ 0.00338, which is a value very close to 0. This means that the weight of the physical skeleton hierarchy is only about 0.34%, while the weight of the dynamic response hierarchy (1-a(t)) is as high as about 99.66%.

[0241] This calculation result shows that at this moment, due to the existence of severe hydraulic dynamic correlation between projects (reflected in the maximum variance of coupling strength), the system decision should almost completely depend on the dynamic response hierarchy generated according to the real-time hydraulic state, and basically ignore the fixed physical connection relationship. This fully embodies the dynamic adaptive ability of the dual-track hierarchical structure of the present application.

[0242] According to one aspect of the present application, by introducing a hydraulic coupling tensor dynamically calculated based on real-time hydrological data, and using the spatial variance of the tensor to weight and fuse the static physical skeleton hierarchy and the time-varying dynamic response hierarchy, a dual-track dynamic hierarchical structure is constructed. The beneficial effects produced by this technical solution are that it changes the cognition of the water project group system structure from static description to dynamic real-time mapping, greatly improving the accuracy and timeliness of system state cognition. Specifically, under extreme flood scenarios, by inputting real-time water level and flow data, the method can quantify and identify temporary strong coupling paths that are not significant in physical connection but play a dominant role under the current hydraulic conditions due to factors such as river backwater and tributary backwater. This directly overcomes the technical defect in the background art that the fixed topological structure cannot capture the dynamic response of the system, resulting in a disconnection between the decision basis and the real system state. Ultimately, this accurate system state cognition provides a high-fidelity battle map for subsequent risk identification, path tracking and control decision-making, ensuring that the control instructions can be formulated based on the most real system hydraulic connection, thereby significantly improving the pertinence and effectiveness of the entire control scheme.

[0243] According to an aspect of the present application, by adopting a top-down tracking mechanism on the basis of a double-track dynamic hierarchical structure, a cross-hierarchical and structured conduction chain of the resilience failure is realized. The technical solution has the beneficial effect that it deepens the granularity of risk diagnosis from the traditional and isolated point failure judgment (such as single reservoir water level overrun) to the complete tracing of the systematic risk conduction chain. Specifically, when a system-level resilience index (such as the large system absorption capacity) fails, the method can input the failure signal and use the hierarchical structure and the hydraulic coupling relationship to locate the specific performance dimension (such as the insufficient resistance of a certain project) that causes the system failure at the subsystem or even single engineering level layer by layer in reverse. This directly solves the technical problems in the background art that the risk cannot be traced in depth and the failure conduction mechanism cannot be understood. By obtaining the structured failure chain of the large system-subsystem-single engineering, the decision maker not only knows where the problem is, but also deeply understands why and how the problem occurs, so as to formulate a control strategy that directly hits the root of the risk with surgical precision, effectively prevents the further conduction and amplification of the risk, and avoids the occurrence of cascading collapse.

[0244] According to an aspect of the present application, by searching for a closed-loop coordination path connecting the failure node and the performance surplus node on the hierarchical control response graph and calculating a quantitative feasibility score of the path that comprehensively considers the path length, response time, control efficiency and net performance improvement, the precise determination of the control potential of the failed project is realized. The technical solution has the beneficial effect that it provides a scientific and quantitative decision basis for subsequent scheduling and optimization of resource allocation, thereby improving the efficiency and rationality of the overall control. Specifically, in the specific scene of water engineering group flood control scheduling, the method can convert the abstract question of whether it can be saved into a calculable score. By comparing the score with a preset threshold, the system can automatically divide the failed projects into two categories: internal coordination control and external strong support. This avoids wasting valuable scheduling resources (such as the precious flood control storage capacity of the upstream reservoir) on failed projects with low internal coordination potential and low success rate of rescue. This data-driven classification enables the subsequent use of different optimization algorithms (such as self-optimization for controllable projects and coordination optimization for uncontrollable projects), realizes the individualized control strategy, and significantly improves the resource utilization efficiency and overall control performance of the water engineering group under complex flood conditions.

[0245] The above describes the preferred embodiments of the present application, but the present application is not limited to the specific details in the above embodiments. Within the technical concept of the present application, various equivalent transformations of the technical solutions of the present application can be made, and these equivalent transformations all belong to the protection scope of the present application.

Claims

1. A water project group risk resilience regulation method considering uncertainty, characterized in that, The method comprises the following steps: Based on the obtained water engineering group working conditions and hydrological data, a double-track dynamic hierarchical structure is established, and a resilience performance boundary and a current resilience performance state are generated; According to the double-track dynamic hierarchical structure, the resilience performance boundary and the current resilience performance state, the regulation potential of the resilience failure project is determined, and the resilience failure project is divided into two categories of regulatable and unregulatable failure project sets; For the two types of failure project sets, self-regulation optimization and collaborative regulation optimization are respectively performed, and an optimal scheduling scheme set is output; The double-track dynamic hierarchical structure is established, comprising: Read the permanent engineering connection relationship from the water engineering group working condition data to construct a steady-state physical skeleton hierarchy; In response to the real-time hydraulic state presented by the hydrological data, a hydraulic coupling tensor is established to quantify the dynamic correlation strength between projects, and a time-varying dynamic response hierarchy is generated based on the hydraulic coupling tensor; And based on the hydraulic coupling tensor, the physical skeleton hierarchy and the dynamic response hierarchy are weighted and fused to obtain the double-track dynamic hierarchical structure; The construction of the hydraulic coupling tensor comprises: According to the real-time hydraulic state, the water level propagation relationship representing the mutual influence of water levels between projects, the flow propagation relationship representing the mutual influence of flows, and the water level-flow coupling relationship representing the influence of water level on downstream flow are calculated respectively; And the three relationships are quantified as a water level propagation matrix, a flow propagation matrix and a water level-flow coupling matrix respectively; The water level propagation matrix, the flow propagation matrix and the water level-flow coupling matrix are combined to synthesize the hydraulic coupling tensor; Wherein, based on the hydraulic coupling tensor, the physical skeleton hierarchy and the dynamic response hierarchy are weighted and fused, comprising: Calculate the spatial variance of the hydraulic coupling tensor, and determine the time-varying hierarchical fusion weight coefficient according to the spatial variance; The hierarchical fusion weight coefficient is used to weight and combine the physical skeleton hierarchy and the dynamic response hierarchy to generate the double-track dynamic hierarchical structure.

2. The method according to claim 1, characterized in that, Before determining the regulation potential of the resilience failure project, it further comprises: Compare the current resilience performance state with the resilience performance boundary to calculate the performance gap; Statistical failure risk probability of performance gap below the preset value under uncertain disturbance; When the absolute mean values of the failure risk probability and the performance gap exceed the corresponding probability threshold and performance gap threshold respectively, the corresponding project is determined as a resilience failure node to form a resilience failure node set.

3. The method according to claim 2, characterized in that, After obtaining the resilience failure node set, the regulation potential of the resilience failure project is determined, further comprising: Based on the current resilience performance state and the resilience performance boundary, the performance surplus nodes are screened to form a regulation source set; For any resilience failure node in the resilience failure node set, search for a closed-loop collaborative path connecting the failure node and at least one regulation source in the double-track dynamic hierarchical structure; Calculate at least one feasibility score of the closed-loop collaborative path to represent the regulation potential of the resilience failure node.

4. The method according to claim 3, characterized in that, The calculation of the feasibility score of the closed-loop collaborative path comprises: Comprehensively evaluate at least two of the following performance indicators of the closed-loop collaborative path: path length, response total time, regulation efficiency, and net performance improvement; Based on the evaluation results of the at least two performance indicators, the feasibility score is generated.

5. The method according to claim 3, characterized in that, The resilience failure project is divided into two categories of regulatable and unregulatable failure project sets, comprising: The feasibility score is compared with a preset feasibility score threshold value; If the feasibility score is not less than the feasibility score threshold value, the corresponding toughness failure node is divided into a controllable failure engineering set; Otherwise, the toughness failure node is divided into an uncontrollable failure engineering set.

6. The method of claim 1, characterized in that, For the two types of failure engineering sets, self-regulation optimization and collaborative regulation optimization are respectively performed, including: For the controllable failure engineering set, a multi-objective reinforcement learning guided hybrid gradient search algorithm is used for optimization; For the uncontrollable failure engineering set, a collaborative-competitive co-evolution algorithm is used for optimization, combined with the matched auxiliary regulation engineering set.

7. The method according to claim 6, characterized in that, For the controllable failure engineering set, a multi-objective reinforcement learning guided hybrid gradient search algorithm is used for optimization, including: Based on the physical constraints of the controllable failure engineering set, a safe exploration space is established to limit the range of regulation actions; Within the safe exploration space, a multi-objective reinforcement learning agent is used to generate an initial regulation action; And a gradient search algorithm is used to fine-tune the initial regulation action, while ensuring that the adjusted action does not exceed the boundary of the safe exploration space, obtaining an optimized regulation action sequence.

Images